一种电力系统潮流计算的导纳矩阵计算方法与流程

本发明涉及一种电力系统潮流计算的导纳矩阵计算方法，特别是一种适合研究目的使用的基于matlab的潮流计算的导纳矩阵计算方法。

背景技术：

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

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

目前常用的牛顿法和快速分解法潮流计算方法都是以节点电压法为基础的，都需要求取节点导纳矩阵。

节点导纳矩阵为：

式中，yik为节点导纳矩阵元素，当下标i≠k时，为节点i和节点k之间的互导纳，当下标i＝k时，为节点i的自导纳；n为节点数。

可以采用追加支路法计算导纳矩阵元素，即依次扫描所有支路，每扫描一条支路在原有导纳矩阵元素基础上增加一个导纳增量。由于输电线路和变压器支路都属于支路，通常把输电线路和变压器支路都作为支路数据统一输入。作为区分，变压器支路的非标准变比侧的节点号加个负号。

输电线路采用如图2所示的π形等值电路，设第m条支路的串联支路导纳为

ym＝1/zm＝1/(rm+jxm)(2)

式中，rm，xm，zm分别为输电线路等值电路的电阻、电抗和阻抗。

追加第m条输电线路时，导纳矩阵元素的计算公式如下：

式中，下标i、j分别表示支路的首节点号im和末节点号jm去掉负号后的节点号，bm为输电线路等值电路的对地电纳，符号“←”表示右端计算结果赋值给左端变量。

变压器支路采用如图3所示的理想变压器串联一等值阻抗的等值电路表示，根据变比及等值阻抗处于位置不同，分为4种情况，图3(a)和图3(b)为等值阻抗位于标准变比侧(即1侧)，图3(c)和图3(d)为等值阻抗位于非标准变比侧(即km侧)。为了降低程序设计的复杂程度，通常把图3(c)和图3(d)的变比转换成(1/km)：1形式，从而把图3(c)和图3(d)所示的等值电路变成图3(a)和图3(b)所示的等值电路。图4(a)和图4(b)分别为图3(a)和图3(b)所示电路的π形等值电路，其中ym＝1/(rm+jxm)。

如图3(a)所示，当变比km位于首节点i侧时，根据其π形等值电路图4(a)得到第m条变压器支路的导纳矩阵元素的计算公式如下：

式中，km为变压器支路的变比。

如图3(b)所示，当变比km位于末节点j侧时，根据其π形等值电路图4(b)得到第m条变压器支路的导纳矩阵元素的计算公式如下：

如图1-5所示，现有潮流计算方法，主要包括以下步骤：

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

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

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

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

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

b、形成节点导纳矩阵；

形成节点导纳矩阵的步骤如下：

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

b2、取支路m的首节点号im、末节点号jm，并令i＝|im|、j＝|jm|；

b3、取支路m的电阻rm、电抗xm，并令ym＝1/(rm+jxm)；

b4、判断支路m的首节点号im、末节点号jm是否都大于0，如果不满足转至步骤b7；

b5、取支路m的对地电纳bm；

b6、按式(2)和式(3)计算输电线路两端节点对应的节点导纳矩阵元素；

b7、判断支路m的首节点号im是否小于0，如果不满足转至步骤b10；

b8、取变压器支路m的变比km；

b9、按式(2)和式(4)计算变压器支路两端节点对应的节点导纳矩阵元素；

b10、判断支路m的末节点号jm是否小于0，如果不满足转至步骤b13；

b11、取变压器支路m的变比km；

b12、按式(2)和式(5)计算变压器支路两端节点对应的节点导纳矩阵元素；

b13、令m＝m+1。

b14、判断m是否大于支路数l，如果m不大于l，则返回到步骤b2；否则，结束。

c、潮流计算迭代主程序；

根据采用的潮流计算的方法不同，可以采用极坐标牛顿法、直角坐标牛顿法、快速分解法进行潮流计算。

d、计算平衡节点的有功功率和无功功率及pv节点的无功功率；

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

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

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

直接采用上述原理实现的潮流计算软件计算速度较慢，商业使用的潮流计算软件采用稀疏矩阵技术和节点优化编号技术，比较复杂，不适合科研人员以此为基础进一步进行科学研究。因此，中国专利zl201710557623.2、zl201710557642.5和zl201710557622.8分别提出基于matlab的极坐标牛顿法潮流计算方法、直角坐标牛顿法潮流计算方法和快速分解法潮流计算方法，可以充分利用matlab特有的擅长矩阵运算和复数运算的特点，并采用matlab的稀疏矩阵技术和方程求解算法，设计出了简洁又有较快计算速度的潮流计算方法，为以潮流计算为基础进行进一步研究的科研人员提供一个易于修改和维护的潮流计算方法，其特点如下：

1、在matlab平台实现，便于科研人员使用matlab提供的各种工具和函数对计算结果进行测试和分析；

2、大部分函数采用矩阵运算和复数运算，减少了程序代码，简化了编程，使得程序更加清晰，便于科研人员修改程序、对程序进行调试和改进、添加新功能；

3、采用矩阵运算、matlab的稀疏矩阵技术并直接调用matlab的方程求解算法，大大提高了计算速度。

上述几项专利提出的基于matlab的潮流计算方法，为从事电力系统研究的科研人员提供了三种基于matlab平台的易于修改和维护且计算较为快速的潮流计算方法。这些方法采用matlab实现，充分利用matlab擅长矩阵运算和复数运算的特点，并使用matlab提供的稀疏矩阵技术和方程求解算法，大大简化了编程且提高了计算速度。但这几种潮流计算方法计算导纳矩阵时未实现矩阵运算，导纳矩阵计算速度相对较慢，仍有待进一步提高计算速度。

技术实现要素：

为解决现有技术存在的上述问题，本发明要提出基于matlab矩阵运算的导纳矩阵计算方法，充分利用matlab特有的擅长矩阵运算的特点，实现提高潮流计算的计算速度的目的。

为了实现上述目的，本发明的技术方案如下：一种电力系统潮流计算的导纳矩阵计算方法，采用矩阵运算。下面推导采用矩阵运算计算导纳矩阵的公式。

一个电力网络的导纳矩阵可以由节点支路关联矩阵和支路导纳求得，为了求节点支路关联矩阵，先规定电力网络中支路的方向，输电线路的方向为首节点指向末节点，变压器支路的方向为变压器非标准变比k侧的节点指向标准变比1侧的节点，图3(a)所示的变压器支路的方向为首节点指向末节点，图3(b)所示的变压器支路的方向为末节点指向首节点。

电力网络的节点支路关联矩阵为

式中，行号对应节点号，列号对应支路号，aik表示节点和支路的关联关系，aik＝1表示支路k从节点i连出，aik＝-1表示支路k向节点i连入，aik＝0表示支路k与节点i不直接相连，n为节点数，l为支路数。

电力网络中不含对地支路时，导纳矩阵为

y1＝aybda^t(7)

式中，ybd为串联支路导纳数组yb形成的对角矩阵，上标t表示矩阵的转置。

第m条输电线路等值电路的串联支路导纳为

ym＝1/zm＝1/(rm+jxm)(8)

式中，rm，xm，zm分别为输电线路等值电路的电阻、电抗和阻抗。

由图4(a)和图4(b)所示的变压器π形等值电路可知，第m条变压器支路π形等值电路的串联支路导纳为

y′m＝1/zm/km＝1/(rm+jxm)/km(9)

式中，rm，xm，zm，km分别为变压器支路等值电路的电阻、电抗、阻抗和变比。

输电线路可以看作变比km＝1的变压器支路，则输电线路和变压器支路π形等值电路的串联支路导纳可以统一写成式(9)的形式。

用矩阵表示支路π形等值电路的串联支路导纳数组为

yb＝1·/(r+jx)·/k(10)

式中，r为支路电阻数组，x为支路电抗数组，k为支路变比数组，输电线路的变比为1，“·/”表示两数组对应元素相除；

由于输电线路和变压器的π形等值电路中都包含对地支路，还应该把这些对地支路的导纳加在导纳矩阵的自导纳元素上。

为了计算各节点对地支路导纳，形成下列矩阵(11)和矩阵(12)。

式中，as表示有向支路首节点与支路的关联矩阵，行号对应节点号，列号对应支路号，aik表示节点和支路的关联关系，支路k从节点i连出，则aik＝1，否则，aik＝0。

式中，ae表示有向支路末节点与支路的关联矩阵，行号对应节点号，列号对应支路号，aik表示节点和支路的关联关系，支路k向节点i连入，则aik＝1，否则，aik＝0。

根据式(6)、式(11)和式(12)的定义，节点支路关联矩阵为

a＝as-ae(13)

由图2可知，输电线路等值电路两侧的对地导纳相同，为

ys＝ye＝jbm/2(14)

式中，ys为有向支路首节点对地电纳，ye为有向支路末节点对地电纳。

由图4(a)和4(b)可知，变压器支路π形等值电路非标准变比k侧节点的对地导纳为

变压器支路π形等值电路标准变比1侧节点的对地导纳为

综合式(14)～式(16)，并写成矩阵形式得到有向支路首、末节点的对地导纳数组为

ys＝yb·*(1-k)·/k+jb/2(17)

ye＝yb·*(k-1)+jb/2(18)

式中，b为对地电纳数组，变压器支路b值为0。

电力网络中各节点对地支路导纳数组y0为

y0＝as·ys+ae·ye(19)

电力网络支路的π形等值电路串联支路的节点导纳矩阵对角元素加上各节点对地支路导纳得到完整的节点导纳矩阵，计算公式如下：

yii＝y1ii+y0ii＝1,…,n(20)

本发明采用矩阵运算形成导纳矩阵，包括以下步骤：

b1、读支路首节点号数组i、末节点号数组j、电阻数组r、电抗数组x、对地电纳数组b、变压器变比数组k；

所述的首节点号数组i、末节点号数组j、电阻数组r、电抗数组x、对地电纳数组b、变压器变比数组k分别按顺序存放所有支路的首节点号im、末节点号jm、电阻rm、电抗xm、对地电纳bm、变压器变比km，变压器支路对地电纳bm为0，变压器非标准变比km侧的节点号加个负号，输电线路变比km设为1，其中下标m为支路序号；

b2、定义as和ae为n×l阶稀疏矩阵；

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

b4、判断支路m的首节点号im是否小于0，如果不满足转至步骤b6；

b5、令i＝|im|、j＝|jm|，转至步骤b7；

b6、令i＝|jm|、j＝|im|；

b7、令asim＝1、aejm＝1；

b8、令m＝m+1；

b9、判断m是否大于支路数l，如果m不大于l，则返回到步骤b4；

b10、计算支路π形等值电路串联支路导纳数组；

用式(10)计算支路π形等值电路串联支路导纳数组yb；

b11、计算节点支路关联矩阵；

用式(13)计算节点支路关联矩阵a；

b12、由串联支路形成初始导纳矩阵；

用式(7)计算由支路π形等值电路串联支路形成的初始导纳矩阵y1；

b13、计算有向支路首节点和末节点的对地导纳数组；

分别用式(17)和式(18)计算有向支路首节点和末节点的对地导纳数组ys和ye；

b14、形成各节点对地导纳数组；

用式(19)计算各节点对地导纳数组y0；

b15、设置节点计数m＝1；

b16、令ymm＝y1mm+y0m；

b17、令m＝m+1；

b18、判断m是否大于节点数n，如果m不大于n，则返回到步骤b16；否则，结束。

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

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

2、本发明提出的导纳矩阵计算采用矩阵运算和复数运算，减少了程序代码，简化了编程，使得程序更加清晰；使用矩阵运算也大大提高了计算速度。

附图说明

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

图1是现有潮流计算的流程图。

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

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

图4是变压器支路的π形等值电路图。

图5是现有潮流计算形成导纳矩阵的流程图。

图6是本发明形成导纳矩阵的流程图。

具体实施方式

下面结合附图对本发明进行进一步地说明，按照图1和图6所示流程对一个10428节点实际系统算例进行了计算，该算例有10428个节点、10436条支路。

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

方法1：中国专利zl201710557623.2方法，导纳矩阵计算采用循环结构及复数运算；

方法2：本发明方法，导纳矩阵计算采用复数运算和矩阵运算。

两种方法的潮流计算和导纳矩阵计算的计算时间见表1，潮流计算的计算时间不包括数据读入和输出的时间。

表1两种极坐标牛顿法潮流计算计算时间比较

从表1可见，中国专利zl201710557623.2方法计算导纳矩阵时间较长，占潮流计算计算时间的79.1％，占潮流计算大部分计算时间；本发明采用矩阵运算技术计算导纳矩阵明显提高了导纳矩阵的计算速度，导纳矩阵计算时间占潮流计算计算时间的42.4％。本发明导纳矩阵计算时间为专利201710557623.2方法的1/5。

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

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