本发明涉及一种电力系统潮流计算的导纳矩阵计算方法,特别是一种适合研究目的使用的基于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。C、潮流计算迭代主程序;根据采用的潮流计算的方法不同,可以采用极坐标牛顿法、直角坐标牛顿法、快速分解法进行潮流计算。D、计算平衡节点的有功功率和无功功率及PV节点的无功功率;平衡节点的有功功率和无功功率及PV节点的无功功率未知,需要计算求出。E、计算各支路有功功率和无功功率;F、输出计算结果,结束。直接采用上述原理实现的潮流计算软件计算速度较慢,商业使用的潮流计算软件采用稀疏矩阵技术和节点优化编号技术,比较复杂,不适合科研人员以此为基础进一步进行科学研究。因此,中国专利CN201710557623.2、CN201710557642.5和CN201710557622.8分别提出基于Matlab的极坐标牛顿法潮流计算方法、直角坐标牛顿法潮流计算方法和快速分解法潮流计算方法,可以充分利用Matlab特有的擅长矩阵运算和复数运算的特点,并采用Matlab的稀疏矩阵技术和方程求解算法,设计出了简洁又有较快计算速度的潮流计算方法,为以潮流计算为基础进行进一步研究的科研人员提供了3种易于修改和维护的潮流计算方法;中国专利CN201710942153.1在上述专利的基础上提出了基于关联矩阵和矩阵运算的导纳矩阵计算方法,进一步提高了潮流计算的速度。但中国专利CN201710942153.1方法计算导纳矩阵时所使用的关联矩阵未实现矩阵运算,导纳矩阵计算速度相对较慢,仍有待进一步提高计算速度。技术实现要素:为解决现有技术存在的上述问题,本发明要提出基于Matlab矩阵运算的导纳矩阵计算方法,充分利用Matlab特有的擅长矩阵运算的特点,实现提高潮流计算的计算速度的目的。为了实现上述目的,本发明的技术方案如下:基于Matlab矩阵运算的潮流计算导纳矩阵计算方法,完全采用矩阵运算计算导纳矩阵。下面推导完全采用矩阵运算计算导纳矩阵的公式。变压器支路的π形等值电路不是对称的,追加支路法计算导纳矩阵元素时,支路的首末节点的自导纳增量不同。为了计算方便,采用有向支路概念,先规定电力网络中支路的方向,输电线路的方向为首节点指向末节点,变压器支路的方向为变压器非标准变比km侧的节点指向标准变比1侧的节点,图3(a)所示的变压器支路的方向为首节点指向末节点,图3(b)所示的变压器支路的方向为末节点指向首节点。根据定义,有向支路的首末节点号数组为:I0=(I<0)·*abs(I)+(I>0)·*abs(J)(6)J0=(I<0)·*J+(I>0)·*I(7)式中,I、J分别为支路原来的首、末节点号数组,其中变压器支路的非标准变比km侧的节点号为负数,I<0表示I数组元素分别与0比较,成立为1,不成立为0,结果仍为一个数组,I0、J0分别为有向支路首、末节点号数组,abs为Matlab求数组元素绝对值函数,“.*”表示两数组对应元素相乘。第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为YB=1·/(R+jX)·/K(10)式中,R为支路电阻数组,X为支路电抗数组,K为支路变比数组,输电线路的变比为1,“·/”表示两数组对应元素相除。电力网络中的导纳矩阵的部分互导纳矩阵为Y1=sparse(I0,J0,-YB,n,n)(11)式中,sparse为Matlab形成稀疏矩阵函数,其参数分别为矩阵的行号数组、列号数组、元素值数组、行数、列数,I0、J0分别为有向支路首节点号数组、末节点号数组,n为节点数。式(11)形成的部分互导纳矩阵是由一半互导纳元素构成的矩阵,另外一半互导纳元素可以通过矩阵转置得到。利用sparse函数形成稀疏矩阵时,如果两个元素的行列号相同,这两个元素会合并。电力网络中导纳矩阵的自导纳元素由输电线路和变压器的π形等值电路串联支路和对地支路计算。由图2可知,输电线路等值电路两侧节点的对地导纳相同,两侧节点自导纳增量为ys=ye=y′m+jbm/2(12)式中,ys为有向支路首节点对地电纳,ye为有向支路末节点对地电纳。由图4(a)和4(b)可知,变压器支路π形等值电路非标准变比km侧节点的自导纳增量为变压器支路π形等值电路标准变比1侧节点的自导纳增量为ye=ym=y′mkm(14)综合式(12)~式(14),并写成矩阵形式得到有向支路首节点的自导纳增量数组YS、末节点的自导纳增量数组YE为YS=YB·/K+jB/2(15)YE=YB·*K+jB/2(16)式中,B为对地电纳数组,变压器支路B值为0。利用Matlab的形成稀疏矩阵函数sparse由YS和YE生成导纳矩阵的自导纳元素构成的对角矩阵为:Y0D=sparse(I0,I0,YS,n,n)+sparse(J0,J0,YE,n,n)(17)式中,Y0D为YS和YE生成的导纳矩阵的对角矩阵。使用式(17)生成矩阵时,连在相同节点的自导纳增量会合并,得到该节点总的自导纳。电力网络支路的π形等值电路串联支路的形成节点互导纳矩阵加上各节点自导纳形成的对角矩阵得到由支路数据计算出的节点导纳矩阵,计算公式如下:Y=Y1+Y1T+Y0D(18)式中,上标T表示矩阵的转置。本发明完全采用矩阵运算形成导纳矩阵,包括以下步骤: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、根据有向支路方向的定义由式(6)和式(7)形成有向支路的首末节点号数组I0和J0;B3、用式(10)计算支路π形等值电路串联支路导纳数组YB;B4、由串联支路形成部分互导纳矩阵;用式(11)计算由支路π形等值电路串联支路形成的部分互导纳矩阵Y1;B5、计算有向支路首节点和末节点的自导纳增量数组;分别用式(15)和式(16)计算有向支路首节点和末节点的自导纳增量数组YS和YE;B6、用式(17)由自导纳增量数组YS和YE形成自导纳对角矩阵Y0D;B7、形成节点导纳矩阵Y;由支路的π形等值电路串联支路形成的部分互导纳矩阵Y1和自导纳对角矩阵Y0D按式(18)计算节点导纳矩阵Y。B8、节点导纳矩阵Y的自导纳追加无功补偿导纳值。与现有技术相比,本发明具有以下有益效果:1、本发明提出的方法在Matlab平台实现,便于科研人员使用Matlab提供的各种工具和函数对计算结果进行测试和分析。2、本发明提出的导纳矩阵计算完全采用矩阵运算和复数运算,减少了程序代码,简化了编程,使得程序更加清晰;使用矩阵运算也大大提高了计算速度。附图说明本发明共有附图6张。其中:图1是现有潮流计算的流程图。图2是输电线路的等值电路图。图3是变压器支路的等值电路图。图4是变压器支路的π形等值电路图。图5是现有潮流计算形成导纳矩阵的流程图。图6是本发明形成导纳矩阵的流程图。具体实施方式下面结合附图对本发明进行进一步地说明,按照图1和图6所示流程对一个10428节点实际系统算例进行了计算,该算例有10428个节点、10436条支路。采用本发明方法和两项已有专利方法对10428节点实际系统算例进行了计算,潮流计算采用极坐标牛顿法,计算时相角单位为弧度,其他量采用标幺值,收敛精度为0.00001。3种潮流计算方法分别为:方法1:中国专利CN201710557623.2方法,导纳矩阵计算采用循环结构及复数运算;方法2:中国专利CN201710942153.1方法,导纳矩阵计算采用复数运算和矩阵运算,其中关联矩阵采用循环结构;方法3:本发明方法,完全采用矩阵运算和复数运算计算导纳矩阵。3种方法的潮流计算和导纳矩阵计算的计算时间见表1,潮流计算的计算时间不包括数据读入和输出的时间。表13种极坐标牛顿法潮流计算计算时间比较潮流计算方法潮流计算计算时间(s)导纳矩阵计算时方法12.6085间(s)2.0626方法20.94130.3987方法30.55060.0091从表1可见,中国专利CN201710557623.2方法计算导纳矩阵时间较长,占潮流计算计算时间的79.1%,占潮流计算大部分计算时间;中国专利CN201710942153.1方法采用关联矩阵和矩阵运算技术计算导纳矩阵明显提高了导纳矩阵的计算速度,导纳矩阵计算时间占潮流计算计算时间的42.4%;本发明完全采用矩阵运算计算导纳矩阵进一步提高了导纳矩阵的计算速度,导纳矩阵计算时间仅占潮流计算计算时间的1.65%。本发明导纳矩阵计算时间为专利CN201710942153.1方法的1/44。本发明可以在任何版本的MATLAB编程语言实现,但建议使用较新版本的MATLAB语言。本发明不局限于本实施例,任何在本发明披露的技术范围内的等同构思或者改变,均列为本发明的保护范围。当前第1页1 2 3