基于改进的雅可比矩阵的直角坐标牛顿法潮流计算方法与流程

本发明涉及一种电力系统的直角坐标牛顿法潮流计算方法，特别适合含小阻抗支路系统的直角坐标牛顿法潮流计算。
背景技术：
：电力系统潮流计算是研究电力系统稳态运行的一项基本计算，它根据电力系统给定的运行条件和网络结构确定整个电力系统的运行状态。潮流计算也是电力系统其他分析的基础，如安全分析、暂态稳定分析等都要用到潮流计算。由于具有收敛可靠、计算速度较快及内存需求适中的优点，牛顿法成为当前潮流计算的主流算法。牛顿法分为极坐标和直角坐标两种形式，其中直角坐标牛顿法潮流计算不需要三角函数计算，计算量相对小一些。在直角坐标牛顿法潮流计算中，节点i的电压采用直角坐标表示为：对正常电力网络，直角坐标牛顿法潮流计算具有良好的收敛性，但遇到含有小阻抗支路的病态网络时，直角坐标牛顿法潮流计算就可能发散。电力系统小阻抗支路可分为小阻抗线路和小阻抗变压器支路，在数学模型上线路可以看作变比为1:1的变压器，因此下面分析时仅以小阻抗变压器支路为例分析。小阻抗变压器模型见图1，变压器的非标准变比k位于节点i侧，阻抗位于标准变比侧。变压器阻抗zij＝rij+jxij很小，导纳为式中，yij、gjj、bjj分别为节点i和节点j之间小阻抗支路的导纳、电导和电纳；rij、xjj分别为节点i和节点j之间小阻抗支路的电阻和电抗。由于小阻抗支路lij的阻抗很小，支路的电压降也很小，因此变压器两端节点的电压应满足：如图2所示，现有直角坐标牛顿法潮流计算方法，主要包括以下步骤：a、输入原始数据和初始化电压根据电力系统节点的特点，潮流计算把电力系统节点分成3类：节点有功功率和无功功率已知、节点电压幅值和电压相角未知的节点称为pq节点；节点有功功率和电压幅值已知、节点无功功率和电压相角未知的节点称为pv节点；节点电压幅值和电压相角已知，节点有功功率和无功功率未知的节点称为平衡节点。电压初始化采用平启动，即pv节点和平衡节点的电压实部取给定值，pq节点的电压实部取1.0；所有电压的虚部都取0.0。这里单位采用标幺值。b、形成节点导纳矩阵设节点i和节点j原来的自电导与自电纳分别为gi0、bi0、gj0、bj0，在它们之间增加一条小阻抗支路后的自导纳和互导纳分别为：式中，yii、yjj分别为节点i和节点j的自导纳；yij为节点i和节点j之间的互导纳；rij、xjj分别为节点i和节点j之间小阻抗支路的电阻和电抗；k为节点i和节点j之间小阻抗支路的变比(如果是输电线支路，变比为1)；c、计算功率及电压偏差pq节点的功率偏差计算公式为：式中，n为电力系统的节点数，ei、fi分别为节点i的电压相量的实部和虚部，gim、bim分别为节点导纳矩阵元素yim的实部和虚部，pis、qis分别为节点i给定的注入有功功率和无功功率，为式中，pig为节点i的电源有功功率，pil为节点i的负荷有功功率，qig为节点i的电源无功功率，qil为节点i的负荷无功功率。pv节点的有功功率及电压偏差计算公式为：式中，vis为节点i给定的电压幅值。平衡节点不参与迭代计算，不需要计算功率偏差或电压偏差。平衡节点的电源有功功率和无功功率以及pv节点的电源无功功率采用式(9)计算。计算各节点功率或电压偏差的最大值，称为最大不平衡量，如果最大不平衡量的绝对值小于给定的收敛精度，转步骤f，否则执行步骤d。d、形成雅可比矩阵j雅可比矩阵j的元素(i≠j时)计算公式如下：pq节点按式(10)-(13)计算雅可比矩阵元素；pv节点按式(10)、(11)、(14)、(15)计算雅可比矩阵元素；平衡节点不计算雅可比矩阵元素。雅可比矩阵j的元素(i＝j时)计算公式如下：式中，ai、ci分别为节点i的计算注入电流相量的实部和虚部，为pq节点按式(16)-(19)计算雅可比矩阵元素；pv节点按式(16)、(17)、(20)、(21)计算雅可比矩阵元素；平衡节点不计算雅可比矩阵元素。e、解修正方程及修正电压实部e、虚部f修正方程为：式中，δw为修正方程不平衡量列向量，δp为有功功率不平衡量列向量，δq为无功功率不平衡量列向量，δv2为电压平方不平衡量列向量，j为雅可比矩阵，δe为电压相量实部列向量，δf为电压相量虚部列向量。电压修正公式为：式中，上标(t)表示第t次迭代。f、输出节点及支路数据。对正常电力网络，牛顿法潮流计算具有良好的收敛性，但遇到含有小阻抗支路的病态网络时，牛顿法潮流计算就可能发散。而电力系统中小阻抗支路普遍存在，收敛性是电力系统潮流计算这类非线性问题的最重要指标，计算不收敛就无法得到方程的解。因此改善直角坐标牛顿法潮流计算针对含有小阻抗支路电力系统的收敛性具有非常重要的意义。中国专利zl201410299531.5披露了一种修改常规直角坐标牛顿法潮流计算雅可比矩阵的方法，该方法依据是：潮流计算收敛时，式(6)中δpi、δqi都趋近于0，则ai和ci等于由节点注入功率给定值pis和qis计算出的ais和cis，即中国专利zl201410299531.5方法计算雅可比元素时采用由给定值pis和qis计算的ai和ci值，改善了潮流计算的收敛性，有效解决了含有电阻为0的小阻抗支路电力系统潮流计算的发散问题。但当小阻抗支路的电阻不为0时，该方法迭代次数增加，收敛性变差，甚至不收敛。中国专利zl201410315785.1提出了一种雅可比矩阵改变的直角坐标牛顿法潮流计算方法，该方法首次迭代和后续各次迭代采用不同的雅可比矩阵计算方法，首次迭代计算雅可比元素时采用由给定值pis和qis计算的ai和ci值，后续各次迭代计算雅可比元素时仍采用传统方法，有效解决了含电阻不为0的小阻抗支路电力系统潮流计算的发散问题，但该方法在电阻不为0的小阻抗支路较多时，迭代次数增加，收敛性变差。技术实现要素：为解决现有技术存在的上述问题，本发明要提出基于改进的雅可比矩阵的直角坐标牛顿法潮流计算方法，该方法可以改善直角坐标牛顿法潮流计算分析含有电阻不为0的小阻抗支路电力系统的收敛性。为了实现上述目的，本发明根据pv节点的电源无功功率不是给定值的特点，提出了基于改进的雅可比矩阵的直角坐标牛顿法潮流计算方法来改善潮流计算收敛性。本发明的首次迭代时pq节点采用由给定值pis和qis计算得到的ai和ci值形成雅可比矩阵元素，首次迭代时pv节点则采用给定值pis和qil计算得到的ai和ci值形成雅可比矩阵元素；后续各次迭代时所有节点都采用传统方法计算雅可比矩阵元素。本发明的技术方案如下：基于改进的雅可比矩阵的直角坐标牛顿法潮流计算方法，包括以下步骤：a、输入原始数据和初始化电压；b、形成节点导纳矩阵；c、设置迭代计数t＝0；d、计算功率及电压偏差，求最大不平衡量δwmax；e、判断最大不平衡量|δwmax|是否小于收敛精度ε；如果小于收敛精度ε，执行步骤i；否则，执行步骤f；f、形成雅可比矩阵；除首次迭代外，雅可比矩阵计算方法仍采用传统方法。首次迭代的雅可比矩阵计算方法根据节点类型采用不同方法。对于pq节点，因潮流计算收敛时，计算功率pi和qi与给定值pis和qis相等，故采用式(25)计算注入电流相量的实部和虚部效果较好；由于pv节点的电源无功功率qig不是给定的，原始数据输入时该值输入的是任意值，因此qis与潮流计算收敛时的计算无功功率qi相差很大，且误差不可预期，采用式(25)计算注入电流相量的实部和虚部效果可能较差，可采用式(26)计算。式中，pis为节点i给定的注入有功功率，qil为节点i的负荷无功功率。形成雅可比矩阵元素的具体步骤如下：f1、按式(10)-(15)计算i≠j时的雅可比矩阵元素；f2、令i＝1；f3、判断是否满足t＝0，如果不满足转步骤f7；f4、判断是否满足节点i是pq节点的条件，如果不满足转步骤f5；如果满足，则按式(25)计算节点i的注入电流相量的实部ai和虚部ci，然后转步骤f8；f5、判断是否满足节点i是pv节点的条件，如果不满足转步骤f9；f6、按式(26)计算节点i的注入电流相量的实部ai和虚部ci，转步骤f8；f7、按式(22)计算节点i的注入电流相量的实部ai和虚部ci；f8、按式(16)-(21)计算i＝j时的雅可比矩阵元素；f9、令i＝i+1；f10、判断i是否大于节点数n，如果i不大于n转步骤f3；否则转步骤g；g、解修正方程及修正各节点电压实部和虚部；h、令t＝t+1，返回步骤d进行下一次迭代；i、输出节点及支路数据。下面分两种情况分析首次迭代的情况。(1)小阻抗支路两端节点都是pq节点首次迭代时，与小阻抗支路有关的修正方程为：式中，ai、aj、bi、bj为与δek、δfk相关的项(k＝1,…,n且k≠i,j)；pi0、pj0、qi0、qj0为除小阻抗支路lij外节点的计算功率。式(27)-(30)中考虑到首次迭代时，电压为电压初值，即电压初值实部为1.0，虚部为0.0。得：|gij|和|bij|较大，与之相比，其它项较小，式(31)-(34)忽略较小量，得-(gij/k2)δei+(gij/k)δej+(bij/k2)δfi-(bij/k)δfj≈-gij/k2+gij/k(35)-gijδej+(gij/k)δei+bijδfj-(bij/k)δfi≈-gij+gij/k(36)(bij/k2)δei-(bij/k)δej+(gij/k2)δfi-(gij/k)δfj≈bij/k2-bij/k(37)bijδej-(bij/k)δei+gijδfj-(gij/k)δfi≈bij-bij/k(38)式(35)乘以bij与式(37)乘以gij相加，得式(39)中由于得δfi≈kδfj(40)由于初值则电压虚部修正后满足公式(2)。式(37)乘以bij，再与式(35)乘以gij相减，得式(41)中由于得δei/k2-δej/k≈1/k2-1/k(42)式(42)整理，得(1-δei)≈k(1-δej)(43)式(43)中，考虑电压实部初值首次迭代后电压实部为式(44)满足公式(2)。式(31)乘以k再加式(32)，得式(33)乘以k再加式(34)，得这样式(31)-(34)经过变换得到式(40)、(44)、(45)、(46)，而式(40)、(44)、(45)、(46)已经不存在小阻抗了，且满足小阻抗支路两端电压关系式(2)。由于小阻抗的影响已经不存在了，因此首次迭代时端点都是pq节点的小阻抗不会引起潮流计算的不收敛。(2)小阻抗支路两端节点分别为pv节点和pq节点设支路首端节点i为pv节点，末端节点j为pq节点。首次迭代时，与小阻抗支路有关的修正方程为式(47)-(50)中考虑到首次迭代时，电压为电压初值，即节点i电压初值实部为vis，虚部为0.0，节点j电压初值实部为1.0，虚部为0.0，得：-2visδei＝0(53)式(53)代入式(51)、式(52)、式(54)，得式(56)-(57)忽略较小量，得-gijδej+bijδfj-(bij/k)δfi≈-gij+gijvis/k(58)bijδej+gijδfj-(gij/k)δfi≈bij-bijvis/k(59)式(58)乘以bij与式(59)乘以gij相加，得式(60)中由于得δfi≈kδfj(61)由于初值则电压虚部修正后满足公式(2)。式(59)乘以bij与式(58)乘以gij相减，得式(62)中由于得1-δej≈vis/k(63)由式(63)，再考虑电压实部初值以及修正量首次迭代后电压实部为式(64)满足公式(2)。式(56)乘以(vis/k)再加式(55)，得式(61)和式(63)代入式(57)，得(cjs+bj0)δej+(-ajs+gj0)δfj+bj≈qjs+bj0-qj0(66)这样式(53)、式(55)～(57)经过变换得到式(61)、(64)、(65)、(66)，而式(61)、(64)、(65)、(66)已经不存在小阻抗了，且满足小阻抗支路两端电压关系式(2)。由于小阻抗的影响已经不存在了，因此首次迭代时端点分别为pv节点和pq节点的小阻抗不会引起潮流计算的不收敛。下面分两种情况分析第2次迭代的情况。(1)小阻抗支路两端节点都是pq节点第2次迭代时，与小阻抗支路有关的修正方程为：把式(22)代入到式(67)-(70)，得：考虑到首次迭代以后，小阻抗支路两端节点电压已经满足把此电压关系代入式(71)-(74)，得：式(75)-(78)忽略较小量，得-(gijej+bijfj)δei/k+(gijej+bijfj)δej+(bijej-gijfj)δfi/k+(gijfj-bijej)δfj≈0(79)-(gijej+bijfj)δej+(gijej+bijfj)δei/k+(bijej-gijfj)δfj+(gijfj-bijej)δfi/k≈0(80)(bijej-gijfj)δei/k+(gijfj-bijej)δej+(gijej+bijfj)δfi/k-(gijej+bijfj)δfj≈0(81)(bijej-gijfj)δej+(gijfj-bijej)δei/k+(gijej+bijfj)δfj-(gijej+bijfj)δfi/k≈0(82)式(80)乘以bij与式(82)乘以gij相加，得式(83)中由于得-fjδej+fjδei/k+ejδfj-ejδfi/k≈0(84)式(82)乘以bij，再与式(80)乘以gij相减，得式(85)中由于得ejδej-ejδei/k+fjδfj-fjδfi/k≈0(86)式(84)乘以ej与式(86)乘以fj相加，得式(87)中由于得δfi≈kδfj(88)由于首次迭代后有则修正后满足公式(2)。式(88)代入式(84)，得δei≈kδej(89)由于首次迭代后有则修正后满足公式(2)。式(75)加式(76)，得式(77)加式(78)，得这样式(75)～(78)经过变换得到式(88)～(91)，而式(88)～(91)已经不存在小阻抗了，且满足小阻抗支路两端电压关系式(2)。由于小阻抗的影响已经不存在了，因此第2次迭代时端点都是pq节点的小阻抗不会引起潮流计算的不收敛。(2)小阻抗支路两端节点分别为pv节点和pq节点设支路首端节点i为pv节点，末端节点j为pq节点。第2次迭代时，与小阻抗支路有关的修正方程除第3个方程改为式(92)外，其它3个方程仍为式(75)、式(76)和式(78)。除了如何得出关于无功偏差的方程式(91)外，其它推导过程同前。式(88)和式(89)代入式(78)，得这样式(75)、式(76)、式(92)、式(78)经过变换得到式(88)、(89)、(92)、(93)，而式(88)、(89)、(92)、(93)已经不存在小阻抗了，且满足小阻抗支路两端电压关系式(2)。由于小阻抗的影响已经不存在了，因此第2次迭代时端点分别为pv节点和pq节点的小阻抗不会引起潮流计算的不收敛。同理可证第2次后各次迭代时小阻抗不会引起潮流计算的不收敛。综上所述本发明方法能够使得潮流计算收敛，但对pv节点的不同处理可能影响迭代次数。由此可见，本发明解决了直角坐标牛顿法潮流计算在求解含有小阻抗支路系统时的收敛性问题。采用现有常规直角坐标牛顿法潮流计算不收敛时，本算法能够可靠收敛，且比现有专利技术迭代次数少。与现有技术相比，本发明具有以下有益效果：1、本发明通过首次迭代采用与以后各次迭代过程不同的雅可比矩阵计算方法，解决了直角坐标牛顿法潮流计算在求解含有小阻抗支路电力系统时的收敛性问题。在首次迭代过程中，pq节点采用节点注入功率给定值计算出的电流相量的实部和虚部计算雅可比元素，而pv节点仅采用节点注入有功功率给定值和负荷无功功率给定值计算出的电流相量的实部和虚部计算雅可比元素，不使用电源无功功率给定值，避免了pv节点的电源无功功率给定值不确定给潮流计算收敛性的不利影响。采用常规直角坐标牛顿法潮流计算不收敛时，本发明能够可靠收敛，且比现有专利技术迭代次数少。2、由于本发明不仅能有效解决常规直角坐标牛顿法潮流计算分析含有小阻抗支路电力系统的收敛性问题，同时也能对正常电力系统进行潮流计算，没有不良影响。附图说明本发明共有附图4张。其中：图1是电力系统小阻抗变压器模型示意图。图2是直角坐标牛顿法潮流计算的流程图。图3是现有发明直角坐标牛顿法潮流计算的流程图。图4是本发明直角坐标牛顿法潮流计算雅可比矩阵形成方法的流程图。具体实施方式下面结合附图对本发明进行进一步地说明。根据图1所示的小阻抗变压器模型，采用图2所示的直角坐标牛顿法潮流计算的流程图和图4所示的直角坐标牛顿法潮流计算雅可比矩阵形成方法的流程图，对一个实际大型电网进行了潮流计算。该实际大型电网有445个节点，含有大量的小阻抗支路。其中，x≤0.001的小阻抗支路有49条，x≤0.0001的小阻抗支路有41条，x≤0.00001的小阻抗支路有22条。其中阻抗值最小的是节点118和节点125之间的小阻抗支路l118-125为x＝0.00000001，变比k＝0.9565，k位于节点118侧。潮流计算的收敛精度为0.00001。为了验证本发明计算含电阻不为0的小阻抗支路电力系统的收敛性，把小阻抗支路l118-125、l60-122及l287-310的电阻改为r＝0.0001。作为对比，同时采用常规直角坐标牛顿法潮流算法及现有专利算法(专利号为zl201410315785.1，如图3所示)对该实际大型电网进行了潮流计算，迭代次数见表1。表1不同潮流方法的迭代结果方法常规算法zl201410315785.1算法本发明算法迭代结果不收敛7次收敛4次收敛由表1可见，对于修改后的445节点实际电力系统算例，常规直角坐标牛顿法潮流算法不收敛，本发明算法和专利zl201410315785.1算法都能够收敛，但本发明算法的迭代次数少3次。不同潮流计算方法各次迭代最大不平衡量见表2。单位为标幺值。表2不同潮流方法各次迭代最大不平衡量迭代序号常规算法zl201410315785.1算法本发明算法0-4754.570367135-4754.570367135-4754.5703671351-3451593.823720038-11.138394991-0.5071799002-886651.468310079-6.163450054-0.0685163943-222023.112200678-1.441071252-0.0013250884-55754.415245002-0.106199006-0.0000007925-13972.568194423-0.0063534556-6386.835620506-0.0001418637-6585.38761914-0.0000000628-378994.7769073519-98508.02584122610-37917.863557986由表2可知，3种方法首次迭代前最大不平衡量相同且很大。首次迭代后，现有专利方法最大不平衡量明显减少，迭代7次收敛；本发明方法最大不平衡量减少速度更快，迭代4次收敛；而常规方法的最大不平衡量则变大，最终发散。本发明可以采用任何一种编程语言和编程环境实现，如c语言、c++、fortran、delphi等。开发环境可以采用visualc++、borlandc++builder、visualfortran等。本发明不局限于本实施例，任何在本发明披露的技术范围内的等同构思或者改变，均列为本发明的保护范围。当前第1页12