基于Matlab稀疏矩阵的快速分解法潮流计算方法与流程

本发明涉及一种电力系统快速分解法潮流计算方法，特别是一种适合研究目的使用的快速分解法潮流计算方法。
背景技术：
：电力系统潮流计算是研究电力系统稳态运行的一项基本计算，它根据给定的运行条件和网络结构确定整个网络的运行状态。潮流计算也是电力系统其他分析的基础，如安全分析、暂态稳定分析等都要用到潮流计算。由于具有收敛可靠、计算速度快及内存需求少的优点，快速分解法成为当前潮流计算的主流方法之一，科研人员经常以快速分解法潮流计算为基础进行进一步地研究。实用的商业软件采用c语言等高级编程语言编写且使用稀疏矩阵技术和节点优化编号等高级技术。这些技术虽然能大幅度提高潮流计算的速度、降低内存占用量，但编程非常麻烦且难以修改和维护，不易增加新的功能，因而不适合科研人员用于研究目的使用。matlab软件以矩阵为最基本的数据单位，可以方便地处理各种矩阵和向量运算，也可以很方便自然地处理复数类型，其指令表达式与数学中常用的形式很接近，还有大量常见且实用的函数，给编程带来很大便利。matlab软件简单易用、代码短小易操作，易于编程和调试，计算功能强大，同时还具有非常强大的可视化图形处理和交互式功能，为科学研究以及工程应用提供了一种高效的编程工具，目前已经成为许多科学领域的基本工具和首选平台，在各种科学和工程计算领域得到了广泛的应用。为了适应越来越多的科研人员需要在matlab平台上以快速分解法潮流计算为基础进行进一步地研究的需求，迫切需要一种基于matlab软件的易于编程、修改和调试的快速分解法潮流计算方法。如图1-2所示，现有快速分解法潮流计算方法，主要包括以下步骤：a、输入原始数据和初始化电压；根据电力系统节点的特点，潮流计算把电力系统节点分成3类：节点有功功率和无功功率已知、节点电压幅值和电压相角未知的节点称为pq节点；节点有功功率和电压幅值已知、节点无功功率和电压相角未知的节点称为pv节点；节点电压幅值和电压相角已知，节点有功功率和无功功率未知的节点称为平衡节点。电压初始化采用平启动，即pv节点和平衡节点的电压幅值取给定值，pq节点的电压幅值取1.0；所有电压的相角都取0.0。这里相角单位为弧度，其他量单位采用标幺值。b、形成节点导纳矩阵；c、形成修正方程的系数矩阵b′和b″并进行因子表分解；潮流计算的基本方程是非线性方程组，通常采用逐次线性化方法迭代求解。线性化得到的方程称为修正方程，用来求电压幅值和相角的修正量。快速分解法修正方程是在极坐标牛顿法潮流计算修正方程基础上解耦并改进得到的。快速分解法修正方程为：b′δθ＝δp/v(1)b″δv＝δq/v(2)式中，δp/v和δq/v分别为有功功率和无功功率不平衡量除以电压幅值后的列向量；δv和δθ分别为电压幅值和电压相角修正量列向量；b′为导纳矩阵的虚部，但计算时不计及支路电阻、对地导纳和非标准变比，导纳矩阵中包含pq节点和pv节点相关的行和列；b″为导纳矩阵的虚部，仅包括与pq节点有关的行和列。d、设置迭代计数t＝0，设置收敛标志kp＝0，kq＝0；e、计算有功功率不平衡量δp；pq节点和pv节点的有功功率不平衡量为：式中，pis为节点i的注入有功功率；vi为节点i的电压幅值；θij＝θi-θj，θi、θj分别为节点i和节点j的电压相角；gij和bij分别为导纳矩阵元素的电导部分和电纳部分；n为节点数。求各节点中有功功率不平衡量绝对值最大的值，称为有功功率最大不平衡量，记为δpmax。f、判断有功功率最大不平衡量绝对值|δpmax|是否小于收敛精度ε；如果小于收敛精度ε，令kp＝1，转到步骤g；否则，解修正方程b'δθ＝δp/v，修正电压相角，令kp＝0，转到步骤h；求解修正方程b′δθ＝δp/v，得到δθ，按下式修正电压相角：θ(t+1)＝θ(t)-δθ(t)(4)式中，上标(t)表示第t次迭代；θ为节点电压相角列向量。g、判断kq是否等于1；如果kq＝1，转到步骤l；h、计算无功功率不平衡量δq；pq节点的无功功率不平衡量为：式中，qis为节点i的注入无功功率；m为pq节点数。求各节点中无功功率不平衡量绝对值最大的值，称为无功功率最大不平衡量，记为δqmax。i、判断无功功率最大不平衡量绝对值|δqmax|是否小于收敛精度ε；如果小于收敛精度ε，令kq＝1，转到步骤j；否则，解修正方程b"δv＝δq/v，修正电压幅值，令kq＝0，转到步骤k；求解修正方程b″δv＝δq/v，得到δv，按下式修正电压幅值列向量v：v(t+1)＝v(t)-δv(t)(6)j、判断kp是否等于1；如果kp＝1，转到步骤l；k、令t＝t+1，返回步骤e进行下一次迭代；l、计算平衡节点功率及pv节点的无功功率，计算支路功率，结束。步骤e和步骤f为p～θ迭代，即通过δp求δθ进而修正θ；步骤h和步骤i为q～v迭代，即通过δq求δv进而修正v。主流快速分解法都是按上述步骤设计方法，即先进行p～θ迭代，后进行q～v迭代。也有文献采用先进行q～v迭代，后进行p～θ迭代的方法。直接采用上述原理实现的快速分解法潮流计算软件计算速度较慢，商业使用的快速分解法潮流计算软件采用稀疏矩阵技术和节点优化编号技术，比较复杂，不适合科研人员以此为基础进一步进行科学研究。因此，中国专利zl201710056328.9提出了一种基于matlab的快速分解法潮流计算方法，可以充分利用matlab特有的擅长矩阵运算和复数运算的特点，设计出了简洁又有较快计算速度的潮流计算方法，其特点如下：(1)在matlab平台实现，便于科研人员使用matlab提供的各种工具和函数对计算结果进行测试和分析；(2)采用矩阵运算和复数运算，减少了程序代码，简化了编程，使得程序更加清晰，便于科研人员修改程序、对程序进行调试和改进、添加新功能；(3)采用矩阵运算并直接调用matlab的三角分解法方程求解算法，大大提高了计算速度。中国专利zl201710056328.9所提出方法以快速分解法潮流计算为基础，为进一步研究的科研人员提供了一个易于修改和维护的快速分解法潮流计算方法。该方法采用matlab实现，充分利用matlab擅长矩阵运算和复数运算的特点，但未使用稀疏矩阵技术，计算速度相对较慢，仍有待进一步提高计算速度。技术实现要素：为解决现有技术存在的上述问题，本发明要提出一种基于matlab稀疏矩阵技术的快速分解法潮流计算方法，在充分利用matlab特有的擅长矩阵运算和复数运算特点的基础上，采用matlab的稀疏矩阵技术，设计出具有较快计算速度的潮流计算方法。为了实现上述目的，本发明的技术方案如下：基于matlab稀疏矩阵的快速分解法潮流计算方法，采用matlab的稀疏矩阵技术，包括以下步骤：a、输入原始数据和初始化电压；电压初始化采用平启动，形成节点电压相量列向量同时形成节点电压幅值列向量初值v(0)和节点电压相角单位相量列向量初值b、记录相关节点类型的节点号；快速分解法修正方程组的方程个数及变量个数与电力系统的节点类型有关，p～θ迭代方程组中没有平衡节点有功功率不平衡量对应的方程和平衡节点相角变量；q～v迭代方程组中仅有pq节点无功功率不平衡量对应的方程和pq节点电压幅值变量。为了提高计算速度，形成方程组系数矩阵及方程右端向量时先不考虑节点类型，形成系数矩阵及方程右端向量后，再去掉无关的行和列。为此，设置两个数组记录有关节点类型的节点号，其中数组bt1记录pq节点和pv节点的节点号，数组bt2记录pq节点的节点号。记录相关节点类型的节点号使用matlab的find函数实现：bt1＝find(bus_type～＝vθ)(7)bt2＝find(bus_type＝＝pq)(8)式中，bus_type为节点类型列向量；～＝为不等于关系运算符；＝＝为等于关系运算符；vθ为平衡节点类型；pq为pq节点类型。c、形成节点导纳矩阵，并转化为稀疏的导纳矩阵y，按bt2提取导纳矩阵y的相应各行，形成仅包含pq节点对应行的稀疏导纳矩阵子矩阵ypq；形成节点导纳矩阵的步骤如下：c1、预定义导纳矩阵y的维数为n×n；c2、根据线路参数和变压器支路参数形成导纳矩阵y的元素；c3、根据无功补偿设备参数修正导纳矩阵y部分对角元素；c4、把导纳矩阵y转化为稀疏矩阵；d、形成修正方程的稀疏系数矩阵b′和b″并进行因子表分解；为了提高计算速度和简化程序，形成系数矩阵b′时不考虑节点类型，形成n阶方阵，然后再按数组bt1记录的节点号提取矩阵元素，去掉平衡节点对应的行和列，形成新的系数矩阵b′；按数组bt2记录的节点号提取矩阵ypq中pq节点对应的列，取矩阵元素的虚部，形成系数矩阵b″。直接调用matlab软件的lu分解法对系数矩阵b′进行三角分解形成下三角矩阵l1和上三角矩阵u1；对系数矩阵b″进行三角分解形成下三角矩阵l2和上三角矩阵u2。分解后得到的矩阵l1、u1、l2和u2都不包含无关的行和列，在迭代过程解方程时不用再提取矩阵元素。e、形成节点注入有功功率和无功功率向量；潮流计算迭代过程中，计算节点有功功率不平衡量向量和节点无功功率不平衡量向量时，要用到节点注入有功功率列向量ps和节点注入无功功率列向量qs，为了提高计算速度，先形成节点注入有功功率向量和节点注入无功功率向量。节点注入有功功率列向量为：ps＝pg-pl(9)式中，ps为节点注入有功功率列向量；pg为节点发电有功功率列向量；pl为节点负荷有功功率列向量。节点注入无功功率列向量为：qs＝qg-ql(10)式中，qs为节点注入无功功率列向量；qg为节点发电无功功率列向量；ql为节点负荷无功功率列向量。形成向量ps和qs时不考虑节点类型，然后再按数组bt1和bt2记录的节点号提取向量元素，去掉多余的元素。按数组bt1记录的节点号提取向量ps需要的元素，去掉平衡节点对应的元素，形成新的向量ps；按数组bt2记录的节点号提取向量qs需要的元素，去掉平衡节点和pv节点对应的元素，形成新的向量qs。f、设置迭代计数t＝0，设置收敛标志kp＝0，kq＝0；g、计算有功功率不平衡量δp，并求有功功率最大不平衡量δpmax；采用matlab矩阵运算和复数运算编程，需要推导出基于矩阵运算和复数运算的功率计算方法。定义节点i的复功率公式为：式中，为节点i的复功率；pi和qi分别为节点i的有功功率和无功功率；为节点电压相量；为节点电流相量的共轭，上标(^)表示复数的共轭。式(11)写成向量相乘的形式为：式中，为节点复功率列向量；为节点电压相量列向量；为节点电流相量的共轭值列向量；.*表示两向量对应的元素相乘。式(12)中，节点电流相量为：将式(13)代入式(12)，得：式中，为节点电压相量列向量；为节点电压相量的共轭值列向量；y为稀疏导纳矩阵；上标(^)表示复数的共轭。节点有功功率p为：式中，p为节点有功功率列向量；re表示取矩阵元素的实部。计算节点有功功率不平衡量的矩阵运算的形成为：式中，δp为节点有功功率不平衡量列向量；ps为节点注入有功功率列向量。求各节点中有功功率不平衡量绝对值最大的值，称为有功功率最大不平衡量，记为δpmax。h、判断有功功率最大不平衡量绝对值|δpmax|是否小于收敛精度ε；如果小于收敛精度ε，令kp＝1，转到步骤i；否则，解式(17)所示的修正方程，然后按式(18)修正电压相角，计算电压相量列向量b'δθ＝δp/v(17)θ(t+1)＝θ(t)-δθ(t)(18)式中，上标(t)表示第t次迭代的值；δθ为节点电压相角修正量列向量。使用步骤d形成的下三角矩阵l1和上三角矩阵u1直接调用matlab软件的解线性方程组算法解修正方程组(17)。计算电压相角后，按下式计算电压相量由于三角函数计算工作量较大，在q～v迭代时，角度没有更新，没有必要重新计算角度的正弦和余弦，因此把式(19)进行简化，用下式计算电压相量式中，为节点i的电压相角单位相量，在q～v迭代时不变，无须重新计算，为式(21)和式(20)写成矩阵运算形式分别为：式中，为节点电压相量列向量，v为节点电压幅值列向量，为节点电压相角单位相量列向量。令kp＝0，转到步骤j。i、判断kq是否等于1；如果kq＝1，转到步骤n；j、计算无功功率不平衡量δq，并求无功功率最大不平衡量δqmax；按式(14)计算节点复功率但q～v迭代时平衡节点和pv节点的无功功率不需要计算，因此把式(14)修改为：式中，为pq节点复功率列向量；为pq节点电压相量列向量；为节点电压相量的共轭值列向量；ypq为仅包含pq节点对应行的稀疏导纳矩阵子矩阵。按下式计算pq节点无功功率qpq：式中，qpq为pq节点无功功率列向量；im表示取矩阵元素的虚部。计算pq节点无功功率不平衡量矩阵运算的形成为：式中，δqpq为pq节点无功功率不平衡量列向量；qs为pq节点注入无功功率列向量。求各节点中无功功率不平衡量绝对值最大的值，称为无功功率最大不平衡量，记为δqmax。k、判断无功功率最大不平衡量绝对值|δqmax|是否小于收敛精度ε；如果小于收敛精度ε，令kq＝1，转到步骤l；否则，解式(27)所示的修正方程，然后按式(28)修正电压幅值，计算电压相量列向量b"δv＝δq/v(27)v(t+1)＝v(t)-δv(t)(28)式中，上标(t)表示第t次迭代的值；δv为节点电压幅值修正量列向量。使用步骤d形成的三角矩阵l2和u2直接调用matlab软件的解线性方程组算法解修正方程组(27)。计算电压幅值后，按式(23)计算电压相量列向量令kq＝0，转到步骤m。l、判断kp是否等于1；如果kp＝1，转到步骤n；m、令t＝t+1，返回步骤g进行下一次迭代；n、计算平衡节点的有功功率和无功功率及pv节点的无功功率，计算各支路有功功率和无功功率，结束。与现有技术相比，本发明具有以下有益效果：1、本发明提出的方法在matlab平台实现，便于科研人员使用matlab提供的各种工具和函数对计算结果进行测试和分析。2、本发明提出的方法采用矩阵运算和复数运算，减少了程序代码，简化了编程，使得程序更加清晰，便于科研人员修改程序、对程序进行调试和改进、添加新功能；使用矩阵运算也大大提高了计算速度。3、本发明采用matlab的稀疏矩阵技术，较大幅度地提高了计算速度，同时matlab的稀疏矩阵使用非常方便，可以像全矩阵一样用行列号直接使用稀疏矩阵的元素，也不需要设计稀疏存储结构。4、本发明修改了导纳矩阵的计算过程，进一步提高了计算速度。实践证明，本发明的方法既方便了科研人员对程序进行编写、修改和调试，同时计算速度也基本接近了在c语言平台上实现的速度，为科研人员的科研工作提供了一个优秀的分析工具。附图说明本发明共有附图4张。其中：图1是现有快速分解法潮流计算的流程图。图2是现有快速分解法潮流计算形成导纳矩阵的流程图。图3是本发明快速分解法潮流计算的流程图。图4是本发明快速分解法潮流计算形成导纳矩阵的流程图。具体实施方式下面结合附图及实施例对本发明进行进一步地说明。按照图3-4所示流程对一个修改后的445节点实际系统算例进行了计算。445节点实际大型电力系统有445个节点，544条支路，含有大量的小阻抗支路。为了对各种方法进行比较，把这些小阻抗支路改为正常阻抗支路以满足各种方法的要求。采用本发明和已有专利方法对445节点实际系统算例进行了计算，计算时收敛精度为0.00001。两种潮流计算方法分别为：方法1：中国专利zl201710056328.9方法，未采用稀疏矩阵技术；方法2：本发明方法。两种方法的潮流计算的计算时间见表1，计算时间不包括数据读入和输出及支路功率计算的时间。表1两种快速分解法潮流计算计算时间比较潮流计算方法计算时间(s)方法10.1487方法20.0575从表1可见，中国专利zl201710056328.9方法没有采用稀疏矩阵技术，计算时间较长；本发明采用稀疏矩阵技术并对某些计算过程根据稀疏矩阵的特点进行改进可以明显提高计算速度，同时采用matlab的稀疏矩阵技术比较简单方便。为了进一步优化算法，本发明在计算导纳矩阵和系数矩阵b'时，采用以下不同处理方法：方法1：导纳矩阵和系数矩阵b'都初始化为n×n阶稀疏矩阵，所有元素初值为零；方法2：导纳矩阵和系数矩阵b'都初始化为n×n阶稀疏矩阵，所有元素初值为零，预留n+2nbr个存储空间，其中nbr为支路数；方法3：导纳矩阵和系数矩阵b'都初始化为n×n阶全矩阵，所有元素初值为零，先计算补偿设备导纳，再计算支路导纳，形成导纳矩阵和系数矩阵b'后再转化为稀疏矩阵；方法4：导纳矩阵和系数矩阵b'都初始化为n×n阶全矩阵，所有元素初值为零，先计算支路导纳，再计算补偿设备导纳，形成导纳矩阵和系数矩阵b'后再转化为稀疏矩阵。几种不同处理导纳矩阵和系数矩阵b'方法的潮流计算的计算时间见表2，计算时间不包括数据读入和输出及支路功率计算的时间。表2几种快速分解法潮流计算计算时间比较潮流计算方法计算时间(s)方法10.0714方法20.0693方法30.0616方法40.0575从表2可见，直接把矩阵初始化为稀疏矩阵虽然比采用全矩阵速度快，但并不是最好的方法，先按全矩阵形成矩阵元素，再转化为稀疏矩阵是比较好的处理方法。在计算导纳矩阵时，先计算支路导纳，再用无功补偿设备参数修正导纳矩阵对角元比较快。本发明可以在任何版本的matlab编程语言实现，但建议使用较新版本的matlab语言。本发明不局限于本实施例，任何在本发明披露的技术范围内的等同构思或者改变，均列为本发明的保护范围。当前第1页12