基于Matlab矩阵运算的快速分解法系数矩阵计算法

摘要

本发明公开了一种基于Matlab矩阵运算的快速分解法系数矩阵计算法，完全采用矩阵运算计算修正方程系数矩阵B'。包括以下步骤：读支路首节点号数组I、末节点号数组J、电抗数组X；由支路电抗计算支路电纳数组B

说明书

技术领域

本发明涉及一种电力系统快速分解法潮流计算的修正方程系数矩阵计算方法，特别是一种适合研究目的使用的基于Matlab的快速分解法潮流计算的修正方程系数矩阵计算方法。

背景技术

电力系统潮流计算是研究电力系统稳态运行的一项基本计算，它根据给定的运行条件和网络结构确定整个网络的运行状态。潮流计算也是电力系统其他分析的基础，如安全分析、暂态稳定分析等都要用到潮流计算。由于具有收敛可靠、计算速度快及内存需求少的优点，快速分解法成为当前潮流计算的主流方法之一，科研人员经常以快速分解法潮流计算为基础进行进一步地研究。实用的商业软件采用稀疏矩阵技术和节点优化编号等高级技术。这些技术虽然能大幅度提高潮流计算的速度、降低内存占用量，但编程非常麻烦且难以修改和维护，不易增加新的功能，因而不适合科研人员用于研究目的使用。

Matlab软件以矩阵为最基本的数据单位，可以方便地处理各种矩阵和向量运算，也可以很方便自然地处理复数类型，其指令表达式与数学中常用的形式很接近，还有大量常见实用的函数，给编程带来很大便利。Matlab软件简单易用、代码短小易操作，易于编程和调试，计算功能强大，同时还具有非常强大的可视化图形处理和交互式功能，为科学研究以及工程应用提供了一种高效的编程工具，目前已经成为许多科学领域的基本工具和首选平台，在各种科学和工程计算领域得到了广泛的应用。为了适应越来越多的科研人员需要在Matlab平台上以快速分解法潮流计算为基础进行进一步地研究的需求，迫切需要一种基于Matlab软件的易于编程、修改和调试的快速分解法潮流计算方法。

快速分解法潮流计算节点功率方程为：

式中，P_i、Q_i分别为节点i的节点有功功率和无功功率；U_i、U_k分别为节点i和节点k的节点电压幅值；θ_ik＝θ_i-θ_k，θ_i和θ_k分别为节点i和节点k的节点电压相角；G_ik、B_ik分别为节点导纳矩阵元素Y_ik的实部和虚部。

功率不平衡量方程为：

式中，ΔP_i、ΔQ_i分别为节点i的节点有功功率不平衡量和无功功率不平衡量；P_is、Q_is分别为节点i给定的节点注入有功功率和注入无功功率；p为PQ节点数。

快速分解法潮流计算就是通过求解方程(2)，得到各节点电压的幅值和相角，进而计算出流过各支路的功率的过程。

如图1所示，现有快速分解法潮流计算方法，主要包括以下步骤：

A、原始数据输入和电压初始化；

原始数据包括线路和变压器支路数据、节点注入有功功率和无功功率、节点电压幅值、节点无功补偿数据，以及收敛精度、最大迭代次数；

输电线路和变压器支路都作为支路数据统一输入，作为区分，变压器支路的非标准变比侧的节点号加个负号；

输电线路采用如图2所示的π形等值电路，图中，r_m、x_m和b_m分别为输电线路等值电路的电阻、电抗和对地电纳。变压器支路采用如图3所示的理想变压器串联一等值阻抗的等值电路表示，根据变比及等值阻抗处于位置不同，分为4种情况，图3(a)和图3(b)为等值阻抗位于标准变比侧(即1侧)，图3(c)和图3(d)为等值阻抗位于非标准变比侧(即k_m侧)。为了降低程序设计的复杂程度，通常把图3(c)和图3(d)的变比转换成(1/k_m):1形式，从而把图3(c)和图3(d)所示的等值电路变成图3(a)和图3(b)所示的等值电路。

根据电力系统节点的特点，潮流计算把电力系统节点分成3类：节点有功功率和无功功率已知、节点电压幅值和电压相角未知的节点称为PQ节点；节点有功功率和电压幅值已知、节点无功功率和电压相角未知的节点称为PV节点；节点电压幅值和电压相角已知，节点有功功率和无功功率未知的节点称为平衡节点。

电压初始化采用平启动，即PV节点和平衡节点的电压幅值取给定值，PQ节点的电压幅值取1.0；所有电压的相角都取0.0。这里相角单位为弧度，其他量单位采用标幺值。

B、形成节点导纳矩阵；

C、形成修正方程的系数矩阵B'和B”并进行因子表分解；

潮流计算的基本方程是非线性方程组，通常采用逐次线性化方法迭代求解。线性化得到的方程称为修正方程，用来求电压幅值和相角的修正量。快速分解法修正方程是在极坐标牛顿法潮流计算修正方程基础上解耦并改进得到的。

快速分解法修正方程为：

B′Δθ＝ΔP/U (3)

B″ΔU＝ΔQ/U(4)

式中，ΔP/U和ΔQ/U分别为有功功率和无功功率不平衡量除以电压幅值后的列向量；ΔU和Δθ分别为电压幅值和电压相角修正量列向量；B′为导纳矩阵的虚部，但计算时不计及支路电阻、对地导纳和非标准变比，导纳矩阵中包含PQ节点和PV节点相关的行和列；B″为导纳矩阵的虚部，仅包括与PQ节点有关的行和列。

形成修正方程的系数矩阵B'的步骤如下：

C1、设置支路计数m＝1；

C2、取支路m的首节点号i_m、末节点号j_m、电抗x_m，并令i＝|i_m|、j＝|j_m|；

C3、计算修正方程系数矩阵B'的元素；

采用追加支路法计算修正方程系数矩阵B'元素，即依次扫描所有支路，每扫描一条支路在原有系数矩阵B'元素基础上增加一个增量。

计算系数矩阵B'元素的公式为：

式中，符号“←”表示右端计算结果赋值给左端变量。

C4、令m＝m+1。

C5、判断m是否大于支路数l，如果m不大于l，则返回到步骤C2；否则，转步骤D。

D、快速分解法潮流计算迭代主程序；

E、计算平衡节点的有功功率和无功功率及PV节点的无功功率；

平衡节点的有功功率和无功功率及PV节点的无功功率未知，需要计算求出。

F、计算各支路有功功率和无功功率；

G、输出计算结果，结束。

直接采用上述原理实现的潮流计算软件计算速度较慢，商业使用的潮流计算软件采用稀疏矩阵技术和节点优化编号技术，比较复杂，不适合科研人员以此为基础进一步进行科学研究。因此，中国专利CN201710557622.8提出基于Matlab的快速分解法潮流计算方法，可以充分利用Matlab特有的擅长矩阵运算和复数运算的特点，并采用Matlab的稀疏矩阵技术和方程求解算法，设计出了简洁又有较快计算速度的快速分解法潮流计算方法，为以快速分解法潮流计算为基础进行进一步研究的科研人员提供了一个易于修改和维护的快速分解法潮流计算方法。但该方法计算修正方程系数矩阵B'时未实现矩阵运算，系数矩阵B'的计算速度相对较慢。为此，中国专利CN201710942196.X提出了用关联矩阵和矩阵运算计算修正方程系数矩阵B'的方法，提高了潮流计算的计算速度，但计算系数矩阵B'时使用的关联矩阵仍采用循环结构，计算速度仍有待于进一步提高。

发明内容

为解决现有技术存在的上述问题，本发明要提出基于Matlab矩阵运算的快速分解法潮流计算修正方程系数矩阵B'的计算方法，充分利用Matlab特有的擅长矩阵运算的特点，实现提高潮流计算的计算速度的目的。

为了实现上述目的，本发明的技术方案如下：基于Matlab矩阵运算的快速分解法系数矩阵计算法，完全采用矩阵运算计算修正方程系数矩阵B'。下面推导完全采用矩阵运算计算修正方程系数矩阵B'的公式。

根据定义，快速分解法修正方程系数矩阵B'元素由支路电抗计算得到。

忽略支路电阻及对地电纳，并设变压器的变比为1，得支路电纳数组B_B为

B_B＝-1·/X>

式中，X为支路电抗数组，“·/”表示两数组对应元素相除；

系数矩阵B'是非常稀疏的矩阵，其部分非对角元素组成的矩阵B1可以使用Matlab的sparse函数由支路电纳数组B_B生成，为

B1＝sparse(abs(I),abs(J),-B_B,n,n)(7)

式中，sparse为Matlab形成稀疏矩阵函数，其参数分别为矩阵的行号数组、列号数组、元素值数组、行数、列数，I、J分别为支路首、末节点号数组，其中变压器支路的非标准变比k_m侧的节点号为负数，n为节点数。

式(7)形成了系数矩阵B'一半的非对角元素，另外一半的非对角元素可以通过矩阵转置得到。利用sparse函数形成稀疏矩阵时，如果两个元素的行列号相同，这两个元素会合并。

系数矩阵B'的非对角矩阵为：

B′_OD＝B1+B1^T>

式中，上标T表示矩阵的转置。

由于系数矩阵B'的非对角矩阵B′_OD的对角元素为0，因此B′_OD的每行元素相加，就得到了系数矩阵B'对应行对角元素的相反数，矩阵每行元素相加可以由下式实现：

BD＝-B′_OD×1_n×1>

式中，BD为系数矩阵B'的对角元素组成n维列向量，1_n×1为元素都为1的n维列向量。

由列向量BD生成系数矩阵B'的对角矩阵B′_D，可以采用Matlab的sparse函数实现如下：

B′_D＝sparse(1:n,1:n,BD,n,n)(10)

式中，1:n表示形成从1到n、步长为1的数组。

系数矩阵B'的非对角矩阵和对角矩阵相加，得到系数矩阵B'为

B′＝B′_OD+B′_D>

本发明采用矩阵运算方法形成修正方程系数矩阵B'，包括以下步骤：

C1、读支路首节点号数组I、末节点号数组J、电抗数组X；

所述的首节点号数组I、末节点号数组J、电抗数组X分别按顺序存放所有支路的首节点号i_m、末节点号j_m、电抗x_m，其中下标m为支路序号；

支路包括输电线路和变压器支路，为了区分二者，变压器支路非标准变比k_m侧的节点号加个负号。

C2、用式(6)由支路电抗计算支路电纳数组B_B；

C3、用式(7)由支路电纳数组B_B形成系数矩阵B'部分非对角元素组成的矩阵B1；

C4、用式(8)形成系数矩阵B'的非对角矩阵；

C5、用式(9)形成系数矩阵B'的对角元素组成的n维列向量BD；

C6、用式(10)生成系数矩阵B'的对角矩阵；

C7、用式(11)计算系数矩阵B'。

与现有技术相比，本发明具有以下有益效果：

1、本发明提出的方法在Matlab平台实现，便于科研人员使用Matlab提供的各种工具和函数对计算结果进行测试和分析。

2、本发明提出的快速分解法修正方程系数矩阵B'的计算完全采用矩阵运算，减少了程序代码，简化了编程，使得程序更加清晰；完全采用矩阵运算计算修正方程系数矩阵B'也大大提高了计算速度。

附图说明

本发明共有附图5张。其中：

图1是现有快速分解法潮流计算的流程图。

图2是输电线路的等值电路图。

图3是变压器支路的等值电路图。

图4是现有快速分解法修正方程系数矩阵B'计算的流程图。

图5是本发明快速分解法修正方程系数矩阵B'计算的流程图。

具体实施方式

下面结合附图对本发明进行进一步地说明，按照图1和图5所示流程对一个10428节点实际系统算例进行了计算，该算例有10428个节点、10436条支路。由于存在电阻大于电抗的小阻抗支路，传统的快速分解法潮流计算不收敛，采用串联补偿法对小阻抗支路进行处理，增加了24个支路和24个节点，此时系统规模为10452个节点、10460条支路。

采用本发明方法和两项已有专利方法对上述实际系统算例进行了计算，计算时相角单位为弧度，其他量采用标幺值，收敛精度为0.00001。3种潮流计算方法分别为：

方法1：中国专利CN201710557622.8方法，修正方程系数矩阵B'采用循环结构计算；

方法2：中国专利CN201710942196.X方法，修正方程系数矩阵B'采用关联矩阵和矩阵运算，其中关联矩阵采用循环结构；

方法3：本发明方法，修正方程系数矩阵B'完全采用矩阵运算。

3种方法的潮流计算和修正方程系数矩阵B'计算的计算时间见表1，潮流计算的计算时间不包括数据读入和输出的时间。

表1 3种快速分解法潮流计算计算时间比较

潮流计算方法潮流计算计算时间(s)系数矩阵B'计算时间(s)方法14.34161.3930方法23.17970.2322方法32.94220.0038

从表1可见，中国专利CN201710557622.8方法计算系数矩阵B'时间较长，占潮流计算计算时间的1/3，占潮流计算较大一部分计算时间；中国专利CN201710942196.X方法采用关联矩阵和矩阵运算技术计算系数矩阵B'明显提高了系数矩阵B'的计算速度，系数矩阵B'计算仅占潮流计算计算时间的1/14。本发明完全采用矩阵运算技术计算系数矩阵B'进一步提高了系数矩阵B'的计算速度，系数矩阵B'计算时间仅占潮流计算计算时间的1/774。本发明系数矩阵B'计算时间为专利CN201710942196.X方法的1/61。

本发明可以在任何版本的MATLAB编程语言实现，但建议使用较新版本的MATLAB语言。

本发明不局限于本实施例，任何在本发明披露的技术范围内的等同构思或者改变，均列为本发明的保护范围。