[发明专利]基于Matlab稀疏矩阵的快速分解法潮流计算方法

摘要

本发明公开了一种基于Matlab稀疏矩阵的快速分解法潮流计算方法，采用Matlab的稀疏矩阵技术，在Matlab平台实现，便于科研人员使用Matlab提供的各种工具和函数对计算结果进行测试和分析。本发明采用矩阵运算和复数运算，减少了程序代码，简化了编程，使得程序更加清晰，便于科研人员修改程序、对程序进行调试和改进、添加新功能；使用矩阵运算也大大提高了计算速度。本发明采用Matlab的稀疏矩阵技术，较大幅度地提高了计算速度，同时Matlab的稀疏矩阵使用非常方便，可以像全矩阵一样用行列号直接使用稀疏矩阵的元素，也不需要设计稀疏存储结构。本发明修改了导纳矩阵的计算过程，进一步提高了计算速度。

主权项

1.基于Matlab稀疏矩阵的快速分解法潮流计算方法，其特征在于：包括以下步骤：A、输入原始数据和初始化电压；电压初始化采用平启动，形成节点电压相量列向量同时形成节点电压幅值列向量初值V⁽⁰⁾和节点电压相角单位相量列向量初值B、记录相关节点类型的节点号；快速分解法修正方程组的方程个数及变量个数与电力系统的节点类型有关，P～θ迭代方程组中没有平衡节点有功功率不平衡量对应的方程和平衡节点相角变量；Q～V迭代方程组中仅有PQ节点无功功率不平衡量对应的方程和PQ节点电压幅值变量；设置两个数组记录有关节点类型的节点号，其中数组bt₁记录PQ节点和PV节点的节点号，数组bt₂记录PQ节点的节点号；记录相关节点类型的节点号使用Matlab的find函数实现：bt₁＝find(bus_type～＝Vθ)                 (1)bt₂＝find(bus_type＝＝PQ)                (2)式中，bus_type为节点类型列向量；～＝为不等于关系运算符；＝＝为等于关系运算符；Vθ为平衡节点类型；PQ为PQ节点类型；C、形成节点导纳矩阵，并转化为稀疏的导纳矩阵Y，按bt₂提取导纳矩阵Y的相应各行，形成仅包含PQ节点对应行的稀疏导纳矩阵子矩阵Y_PQ；形成节点导纳矩阵的步骤如下：C1、预定义导纳矩阵Y的维数为n×n；C2、根据线路参数和变压器支路参数形成导纳矩阵Y的元素；C3、根据无功补偿设备参数修正导纳矩阵Y部分对角元素；C4、把导纳矩阵Y转化为稀疏矩阵；D、形成修正方程的稀疏系数矩阵B′和B″并进行因子表分解；形成系数矩阵B′时不考虑节点类型，形成n阶方阵，然后再按数组bt₁记录的节点号提取矩阵元素，去掉平衡节点对应的行和列，形成新的系数矩阵B′；按数组bt₂记录的节点号提取矩阵Y_PQ中PQ节点对应的列，取矩阵元素的虚部，形成系数矩阵B″；直接调用Matlab软件的LU分解法对系数矩阵B′进行三角分解形成下三角矩阵L₁和上三角矩阵U₁；对系数矩阵B″进行三角分解形成下三角矩阵L₂和上三角矩阵U₂；分解后得到的矩阵L₁、U₁、L₂和U₂都不包含无关的行和列，在迭代过程解方程时不用再提取矩阵元素；E、形成节点注入有功功率和无功功率向量；先形成节点注入有功功率向量和节点注入无功功率向量；节点注入有功功率列向量为：P_s＝P_G‑P_L                            (3)式中，P_s为节点注入有功功率列向量；P_G为节点发电有功功率列向量；P_L为节点负荷有功功率列向量；节点注入无功功率列向量为：Q_s＝Q_G‑Q_L                           (4)式中，Q_s为节点注入无功功率列向量；Q_G为节点发电无功功率列向量；Q_L为节点负荷无功功率列向量；形成向量P_s和Q_s时不考虑节点类型，然后再按数组bt₁和bt₂记录的节点号提取向量元素，去掉多余的元素；按数组bt₁记录的节点号提取向量P_s需要的元素，去掉平衡节点对应的元素，形成新的向量P_s；按数组bt₂记录的节点号提取向量Q_s需要的元素，去掉平衡节点和PV节点对应的元素，形成新的向量Q_s；F、设置迭代计数t＝0，设置收敛标志K_P＝0，K_Q＝0；G、计算有功功率不平衡量ΔP，并求有功功率最大不平衡量ΔP_max；采用Matlab矩阵运算和复数运算编程，需要推导出基于矩阵运算和复数运算的功率计算方法；定义节点i的复功率公式为：式中，为节点i的复功率；P_i和Q_i分别为节点i的有功功率和无功功率；为节点电压相量；为节点电流相量的共轭，上标(^)表示复数的共轭；式(5)写成向量相乘的形式为：式中，为节点复功率列向量；为节点电压相量列向量；为节点电流相量的共轭值列向量；.*表示两向量对应的元素相乘；式(6)中，节点电流相量为：将式(7)代入式(6)，得：式中，为节点电压相量列向量；为节点电压相量的共轭值列向量；Y为稀疏导纳矩阵；上标(^)表示复数的共轭；节点有功功率P为：式中，P为节点有功功率列向量；Re表示取矩阵元素的实部；计算节点有功功率不平衡量的矩阵运算的形成为：式中，ΔP为节点有功功率不平衡量列向量；P_s为节点注入有功功率列向量；求各节点中有功功率不平衡量绝对值最大的值，称为有功功率最大不平衡量，记为ΔP_max；H、判断有功功率最大不平衡量绝对值|ΔP_max|是否小于收敛精度ε；如果小于收敛精度ε，令K_P＝1，转到步骤I；否则，解式(11)所示的修正方程，然后按式(12)修正电压相角列向量θ，计算电压相量列向量B'Δθ＝ΔP/V                            (11)θ^(t+1)＝θ^(t)‑Δθ^(t)                         (12)式中，上标(t)表示第t次迭代的值；Δθ为节点电压相角修正量列向量；ΔP/V为有功功率不平衡量除以电压幅值后的列向量；使用步骤D形成的下三角矩阵L₁和上三角矩阵U₁直接调用Matlab软件的解线性方程组算法解修正方程(11)；计算电压相角后，按下式计算电压相量由于三角函数计算工作量较大，在Q～V迭代时，角度没有更新，没有必要重新计算角度的正弦和余弦，因此把式(13)进行简化，用下式计算电压相量式中，为节点i的电压相角单位相量，在Q～V迭代时不变，无须重新计算，为式(15)和式(14)写成矩阵运算形式分别为：式中，为节点电压相量列向量，V为节点电压幅值列向量，为节点电压相角单位相量列向量；θ为节点电压相角列向量；令K_P＝0，转到步骤J；I、判断K_Q是否等于1；如果K_Q＝1，转到步骤N；J、计算无功功率不平衡量ΔQ，并求无功功率最大不平衡量ΔQ_max；按式(8)计算节点复功率但Q～V迭代时平衡节点和PV节点的无功功率不需要计算，因此把式(8)修改为：式中，为PQ节点复功率列向量；为PQ节点电压相量列向量；为节点电压相量的共轭值列向量；Y_PQ为仅包含PQ节点对应行的稀疏导纳矩阵子矩阵；按下式计算PQ节点无功功率Q_PQ：式中，Q_PQ为PQ节点无功功率列向量；Im表示取矩阵元素的虚部；计算PQ节点无功功率不平衡量矩阵运算的形成为：式中，ΔQ_PQ为PQ节点无功功率不平衡量列向量；Q_s为PQ节点注入无功功率列向量；求各节点中无功功率不平衡量绝对值最大的值，称为无功功率最大不平衡量，记为ΔQ_max；K、判断无功功率最大不平衡量绝对值|ΔQ_max|是否小于收敛精度ε；如果小于收敛精度ε，令K_Q＝1，转到步骤L；否则，解式(21)所示的修正方程，然后按式(22)修正电压幅值，计算电压相量列向量B"ΔV＝ΔQ/V                           (21)V^(t+1)＝V^(t)‑ΔV^(t)                         (22)式中，上标(t)表示第t次迭代的值；ΔV为节点电压幅值修正量列向量；ΔQ/V为无功功率不平衡量除以电压幅值后的列向量；使用步骤D形成的三角矩阵L₂和U₂直接调用Matlab软件的解线性方程组算法解修正方程(21)；计算电压幅值后，按式(17)计算电压相量列向量令K_Q＝0，转到步骤M；L、判断K_P是否等于1；如果K_P＝1，转到步骤N；M、令t＝t+1，返回步骤G进行下一次迭代；N、计算平衡节点的有功功率和无功功率及PV节点的无功功率，计算各支路有功功率和无功功率，结束。