潮流计算软件设计

1、电力系统潮流软件设计原始数据的输入目前计算机的速度和计算方法已经使我们能很快得到计算结果，但是上机以前的准备 工作却非常耗费时间，而且也容易出错。因此在程序设计时，必须尽可能减轻上机前的准 备工作，尽可能利用计算机代替人工繁琐的工作。在这里，原始数据的填写格式是很关键 的一个环节，它与程序使用的方便性和灵活性有着直接的关系。原始数据输入格式的设计，主要是从使用的角度出发，原则是简单明了，便于修改。 输入格式的简单明了就可以减轻数据填写的工作量，并减少或避免程序使用者在填写数据 时发生错误。电力系统潮流计算往往需要进行多种运行方式的调整和比较，因此在数据格 式上考虑计算过程中修改数据的方便性就显

2、得非常重要。以下所介绍的本潮流程序中所用到的6个信息：N：系统节点的总数；%：系统中支路数，即输电线路条数、变压器数的总和；Ng：发电机节点总数；M：负荷节点总数；向：系统平均电压，在迭代过程中，以它作为电压的初值；eps :迭代收敛所要求的精确度。潮流计算程序所需要的原始数据，分别归纳为下几个结构体：支路数据结构体Branch_Type，有5个数据成员，对应于每条支路的5个数据。该结 构体数组定义为：struct Branch_Typeint i,j;double R, X, YK; Branch400;当支路为输电线路时，这5个数据成员分别表示：i：输电线路一端的节点号；j：输电线路另一端

3、的节点号；R：输电线路的电阻；X：输电线路的电抗；Y0：输电线路充电电容的容纳。十丫0/2Y0/2 一z1:K图1-1网络支路的等值电路如图1-1 (a)所示。当支路为变压器支路时，这5个数据成员分别表示：i:变压器一端的节点号；J:变压器另一端的节点号，这两个节点号有一个带有负号，作为变压器支路的标志；R变压器的电阻；X：变压器的电抗（RT和XT都是归算到变压器标准变比侧的数值）；网：变压器的非标准变比（设在节点号为负的一侧）；变压器的模拟电路如图1-1 （b）所示。发电机节点和负荷节点数据分别定义为结构体数组，结构体有4个数据成员，其内容 是相同的。发电机节点数据定义为：struct Ge

4、nerator_Typedouble P, Q;int i;double V; Generator50;负荷节点数据定义为：struct Load_Typedouble P, Q;int i;double V; Load300;对于发电机节点，P、Q填正号；对于负荷节点，P、Q填负号。对于发电机节点和负荷节点，若为PQ节点，这些数据成员分别表示：P：节点的有功功率；Q：节点的无功功率；i：节点的编号；V：该节点正常运行的电压。若节点为PV节点时数据成员分别表示：P：节点的有功功率；Q：节点无功功率的上限；i：节点的编号；V：节点需要维持的电压，负号是PV节点的标志。在节点数据输入计算机后，为了

5、提高计算效率，应统计PV节点的总数Npv，并形成 PV节点结构体数组。PV节点结构体定义为：struct PVNode_Typedouble V;int i; PVNode100;每个PV节点有两个数据，第一个数据为PV节点的给定电压V，第二个数据为相应的 节点号，。在形成PV结构体数组的同时，把发电机节点或负荷节点数组中V前面负号去掉。数据优化在实用潮流计算程序中，对数据的输入次序应不加以限制，这样便于数据的填写和修 改。输入以后，在计算机内对数据再进行排队和整理。2.1负荷节点的优化为了简化程序，需要对负荷节点的顺序进行优化，优化后各节点按照节点编号的顺序 进行排列，该部分程序的实现框图为

6、图2-1。图2-1负荷节点优化框图其中定义了一个中间工作结构体数组LoadN300，该结构体定义为:struct Load_Typedouble P, Q;int i;double V; LoadN300;用来存放优化后的节点参数。2.2发电机节点的优化由于发电机节点的优化与负荷节点优化相似，因此只给出程序框图（图2-2）。图2-2发电机节点优化框图2.3对支路数据的排队整理优化为了使支路数据的排列方式适合形成导纳矩阵的上三角部分，整理过后支路数据按以 下次序排列：（1）两端节点号应把小号排在前边，大号排在后边。（2）各支路按其小节点号的顺序排列。实现这部分的框图为图2-3。首先对框图中的符号

7、做一介绍：for n=1 to Nbfalseabs(Branchn.i)abs(Branchn.j)truet=Branchn.iBranchn.i=Branchn.jBranchn.j=tc=0BranchN400V. 一个中间工作结构体数组，该结构体定义为：struct BranchN_Typeint i,j;double R, X, YK; BranchN400;用来存放优化后的支路数据。m伙:用来存放小节点号为i支路的另一节点号。sk:用来存放该支路原先编号。该函数首先将原始数据中节点按照（1）的要求优化，然后按照（2）的要求对支路顺 序进行再次优化。2.4对支路特殊情况的优化当某些

8、支路的电阻R大于电抗X的2倍时，影响到尸己分解法的收敛性，这时程序将 自动把支路分成两条支路（增加一个节点），满足PQ分解法要求电阻小于电抗的条件。由图2-4构成的程序来实现这部分功能。对框图中变量的说明：a：用来存放原始数据中支路的个数。i：为计数变量。该函数的实现思想是：当检查出某个支路满足上述情况时，则在该支路中增加一个节 点，节点号为N+1；增加的支路号为Nb+1，原支路仅保留电阻，新支路参数为原支路的 电抗，而对地支路不变。对支路特殊情况的优化框图图2-42.5平衡节点的优化要求把平衡节点排在最后，即作为第N个节点，同时为了简化程序，去掉一些判断还要求这个节点为负荷节点。如果平衡节点

9、没有负荷，则该节点的负荷功率填零。这样保证了节点N既是发电机节点，又是PV节点，又是负荷节点。false Generatori：P=0 andGeneratori：Q=0for i=1 to Nb、fa、abs.(Branchi.i)=jsef faabs(Branchi：j)=jlsetrueti=2truea=i, j=Generatori.iGeneratori.i=jGeneratora.i=Nfor i=1 to NlfalseLoadi.i=NtrueLoadi.i=jfor i=1 to Nbabs(Branchi.i)=Nabs.(Branchi：j)=Ntruesi=2fal

10、sefalseti=1truefalseBranchi：j0trueBranchi.j=-NBranchi.j=N*alsetrue_ _z fn 1 qpBranchi：i0trueBranchi.i=-jBranchi.i=jalse - .si=2truel _falseBranchi：j0trueBranchi.j=-jBranchi.j=j图2-5用来实现这一功能。对框图中变量的说明：J：用来存放原始数据中平衡节点的节点号。ti:为一个标识数组。ti=1表示第i条支路的小节点号为j，ti=2表示第i条支路的 大节点号为j。si：也为一个标识数组。si=1表示第i条支路的小节点号为N，

11、si=2表示第i条支 路的大节点号为N。该函数的设计思想是：先把平衡节点找出（如图2-5 A部分），然后与最后一个节点 进行交换（如图2-5 B、C部分），同时保持原有支路数据形式不变（如图2-5 D部分）。稀疏导纳矩阵的形成3.1基本公式当电力系统中i、j两点间输电线路的阻抗为勺时，节点i、j之间互导纳为Y.j=-=-yij(3-1)ij式中：y.是阻抗z.的倒数，即输电线串联支路的导纳；Yj是导纳矩阵中i行j列的非对角 元素。由于导纳矩阵的对称性，一般Yij=Yj支路i、j对导纳矩阵中i、j两行对角元素的影响可表示为如下的增量：(3-2)丫/3力二_!二%ij这里导纳矩阵对角元素Y.和匕也

12、就是节点i、j的自导纳。当节点i连有导纳为Y.的接地支路时，它对导纳矩阵的影响仅仅使i行对角元素增加 01如下的分量：Yf.(3-3)zK:1 ij(a)(b)图3-1变压器的等值电路当i、j两节点间的支路是非标准变比的变压器时如图3-1 (a)所示，我们可用n型等 值电路来模拟见图3-1 (b)所示，因此i、j之间互导纳可按下式计算：v 11Y=-=- y.lj Kz K ljiji、j节点自导纳分别有如下的增量：(3-5)(3-6),1Y.=y.11 K 2 jY.=y.jj jj3.2稀疏导纳矩阵的处理电力系统的导纳矩阵不仅具有对称性，而且具有稀疏性。当i、j节点之间没有直接联 系时，导

13、纳矩阵中非对角元素Y.及Y.应为零。由于导纳矩阵的对称性，在计算机中可以只存储其上三角部分或下三角部分，在本程序中，储存导纳矩阵的上三角部分及对角元素，因此，其中每个非对角元素Y.的下标都满ij足 ij。对角元素按节点编号顺序存放在对角元素结构体数组中：struct Yii_Typedouble G, B; Yii300;其中G存放对角元素的实部，B存放对角元素的虚部，每个数组的元素个数与系统节点数 相等。由于上三角矩阵中非对角元素和系统中不接地支路一一对应，非对角元素的个数等于 网络中不接地支路数。为了节约内存和提高计算速度，在计算机内存中只储存非零元素， 把非零非对角元素“挤实”在一起。为

14、了识别非对角元素的行号和列号，我们在每个元素 后存放相应的列下标。非对角元素结构体数组定义为：struct Yij_Typedouble G, B;int j; Yj400;按照这样的排列，取一个互导纳，就可以同时把该元素的列号取出来，为了判断该元 素的行号，需要借助于数组NYseq300o数组NYseq300按导纳矩阵行号的顺序存放各行 非对角元素的首地址（事实上，存放的是各行第一个非对角元素在导纳矩阵非零非对角元 素中的顺序号）。本程序中还定义导纳矩阵中各行非对角元素的个数NYsum300。一般我们有NYsumi= NYseq i+1 -NYseq i（3-7）由于非对角元素是逐行向下排列

15、的，所以就很容易判断出各非对角元素的行号。按照 对支路原始数据处理的要求，可以得出支路排列的顺序和导纳矩阵非对角元素的排列顺序 完全一样。因此只要顺序取出支路数据，按照式（3-1）求出倒数取负号之后，连同该支路 的大节点号（即列下标）顺序送入数组，就形成了导纳矩阵的上三角部分。3.3导纳矩阵形成过程及框图在本程序中，适应P-Q分解法的需要，导纳矩阵分为两步形成。第一步只用不接地支路构成导纳矩阵，不考虑接地支路（包括变压器非标准变比）的 影响。这里同时形成两个导纳矩阵，即通常意义上的系数矩阵Y和不考虑输电线路电阻的 系统导纳矩阵Y1，以适应BX法的要求，这两个导纳矩阵实际上只是半成品。Y主要用来

16、 形成BX法所要求的第一个因子表，当该因子表形成后，就在半成品的基础上把接地支路 及变压器非标准变比的影响加进去，形成完整的系统导纳矩阵。其中Y1用来形成BX法第 二个因子表，而Y将在整个迭带求解过程及线路潮流计算过程中发挥作用。只考虑不接地支路构成导纳矩阵的程序框图如图3-2所示。整个形成过程需要把不接 地支路扫描一遍，对每条不接地支路作两方面的工作。首先把阻抗求倒数并取负号后连同 大节点号送到导纳矩阵非对角元素数组Yj（对应于Y）和Yj1（对应于Y）中形成非对角 元素，然后把阻抗的倒数累加到该支路两端节点的自导纳上去见式（3-2）。为了累加形成对角元素，在计算开始时应对数组Yii和Yii1

17、清零见图3-2中、框。 框中NYsum为临时工作数组，定义为NYsum300，在其中累计导纳矩阵各行非对角元 素的个数，因此也需要预先清零。由于两个导纳矩阵的结构是相同的，共用一个NYsum数 组。(10)(11)(12)(13)(14)(15)(16)*i=abs(Branchn.i)j二abs(Brgnchn.j)R二Branchn.R, X=Branchn.XYK二Branchn.YK, Zmag2二R*R+X*XGij=R/Zmag2, Bij=-X/Zmag2b ij=-1.0/X;Branchn.i0falseor. - - - -Branchn.j0I trueYijn.G=-G

18、ij/YK, Yijn.B=-Bij/YKYij1n.G=0, Yij1n.B=-b_ij/YKYijn.j=j, Yij1n.j=jYijn.G=-Gij, Yijn.B=-BijYij1n.G=0, Yij1n.B=-b_ij对不接地支路的扫描用一个循环语句来完成图中框。框把支路的有关数据送进中间工作单元，因为支路为变压器节点号可能为负，所以 在这里对节点号取绝对值。在框中，把阻抗的倒数即支路导纳放到中间工作单元Gj、Bij中，而把支路电抗的 倒数放到中间工作单元。_j中。框判断支路是否为变压器支路。若为变压器支路，则导纳需除以变压器非标准变比 后再取负号送到导纳矩阵非对角元素数组Yij和

19、Yij1中；否则直接取负号送到导纳矩阵非 对角元素数组中框。框向Yij和Yj1数组送列号。这样就把支路阻抗数据变成了导纳矩阵的上三角部分。在(12)框中根据式(3-2)累计有关节点的自导纳。在(13)框中统计小节点号的不接地支路数目，这也就是导纳矩阵上三角部分每行非对角 元素的个数。至此，完成了一条不接地支路的处理；当循环由1做到Nb时形成了只考虑不接地支 路的导纳矩阵。(14)(16)框是由NYsum数组根据式(3-7)形成NYseq数组。3.4追加接地支路的程序框图(10)(12)(13)(14)图3-3追加对地支路框图(11)追加接地支路包括两部分内容，即追加对地电容支路和考虑变压器非标

20、准变比的影 响，其程序框图如图3-3所示。整个计算过程需要对支路数据再进行依次扫描，扫描是由一个循环语句来控制图中 框。在框中把支路有关数据送入中间工作单元。框根据节点号i、j的符号判断所取的 支路是输电线路还是变压器支路。当i、j中任一个为负时，为变压器支路，否则为输电线 路。当所取的支路为输电线路是，转入(12)(14)框，向相应的节点累计自导纳部分。当所取支路为变压器支路时，转入(11)框。在框中判断非标准变比设在支路的哪 一侧。如1章中所述，当i0时，非标准变比就设在i侧，否则设在j侧。由图3-1 (b) 1.，一 及式(3-5)可知，在非标准变比侧自导纳应累计&y.，但在形成不接地支

21、路的导纳矩阵时，该点自导纳累计了*，因此需要再追加累计(K-1)K=(1-k)yijO图中框就是 完成这些运算的。在非标准变比侧自导纳应累计y.，在形成不接地支路的导纳矩阵时，该 点自导纳同样累计了孔，因此需要再追加累计(K-1)笔二(l-K)Y，在框和(11)框中完成KKij这些运算。这样，顺次把支路数据扫描、处理一遍，就形成了描述网络的完整的导纳矩阵。稀疏系数矩阵线性方程式的求解4.1修正方程式的解法及计算公式在P-Q分解法潮流计算的迭代过程中，需要反复求解修正方程式p/v= B ev(4-1)q/v= B Av(4-2)如前所述，这两个方程的系数矩阵(B和B)在迭代过程中保持不变，只要求

22、对不断变化 的常数项(即误差项Ap/v、Aq/v)求解出相应的修正量Aev及Av。在这种情况下，可 以先将系数矩阵进行三角分解，然后只要用分解出的三角矩阵(或因子表)对不同常数项 进行前代及回代的运算，即可得到要求的修正量。系数矩阵的三角分解可以利用递推公式求得，也可以利用高斯消去法求得，这两种方 法在运算量及内存量上都是等效的。本程序利用高斯消去法对系数矩阵进行三角分解并形 成因子表的计算方法。在P-Q分解法潮流程序中，为了在迭代过程中轮流求解式(4-1)及式(4-2)需要形成两 个因子表，为此，可以把式(4-1)及式(4-2)统一为如下的形式：BAX=Al(4-3)在p-e迭代时，式中b即

23、b，Ax为Aev，Ai为Ap/v；在q-v迭代时，式中b为 b ，Ax 为Av，Ai 为Aq/v。以下简单归纳以下形成因子表及常数项进行前代及回代运算的有关公式，对于式(4-3) 的系数矩阵进行B三角分解以后，可以得到以下形式的因子表：B(1)B(1)BB121314111B1B (2)B (2)21B(1)2324221BB(1)1B (3)3132B (2)3433B1,n-1B21(2)B (4-4)1B(n - 2)n-1,n-1B B (i) B (2) B n-1,1n-1,2n-1,3n-1,4因子表中上三角部分元素组成了上三角矩阵U：U=1U121U13U231U14U24U3

24、4.Un-1Un-1U3：n-1:1Un - 2, n - 11其中元素(4-5)U.=B.(i)(ij)因子表中对角元素组成了对角矩阵D：DD22D3311D二D1-1,-1 -1其中元素.(4-6)11 B (i-1)式(4-5)中B.(i)为因子表中上三角部分i行j列元素，其上标(i)表示此元素由原来系 数矩阵元素B.经过i次运算得来。这i次运算中包括i-1次消去运算及一次规格化运算：Bjk)=Bjk-i) -B.k(k-1)Bkj(k)=Bij(k-1) -Bk-DUkj TOC o 1-5 h z lJ (k=1, 2,，i-1; j=k+1, k+2, .：n-1) (4-7)B.

25、二& (j=i+1, i+2, ., n-1)(4-8)ljB(i-1)ii由于系数矩阵B为对称矩阵，在因子表中不需要保留其下三角部分。在形成因子表的 过程中需要用到下三角部分的元素B) (kNUsumi_above-falsecount二count+1, nu=nu+1false true- Un_u.j=i -Btemp=Un u.value/Di above_ trueco.untNusumi_abo.veU falsej=Un_u.j, Bj=Bj-Btemp*Un_u.value count二count+1, n u=n u+1要求置1或置2,作为参数传个子程序。5：为PV节点数组的

26、计数变量，在形成第二个因子表时，用它可以判断系数矩阵 中应该去掉哪些列和哪些行。i：为因子表正在形成的行号变量。j：为列下标变量。n_：是因子表上三角矩阵元素计数变量。i_above：为消去行号计数变量（是被消行号计数变量，因子表按行消去过程中，依次 取行上面的到行）。i_pv：为PV节点的节点号变量。Btemp:为临时变量。count:为临时计数变量。形成因子表所需要的原始数据可由以下数组取得，这些数组的定义及内容详见1章和 3章：Load :负荷功率数组。PVNode： PV节点数组。3j：导纳矩阵非对角元素数组。NYseq :导纳矩阵各行非对角元素的首地址数组。Yii：导纳矩阵对角元素数

27、组。形成的因子表将放在以下数组中：NUsum :存放因子表上三角矩阵各行非对角元素数。Q：存放因子表的对角元素。U：存放因子表上三角矩阵元素，定义为结构体数组：struct U_Typedouble value;int j; U1300,U2300;其中存放该元素的数值，存放该元素的列下标。最后，在框图中还有一个很重要的工作数组8，形成因子表的运算主要在这个数组 中进行。现在分别介绍框图中各个部分的作用。图4-1中A、B、C三部分的作用是“传递”，通过这三部分的工作把系数矩阵中待消 行（i行）的元素按其下标稀疏的排列在工作数组B中。D框的作用是把工作数组B中的 待消行元素按照式（4-7）进行消

28、去运算。最后通过E框的工作把数组B中的元素按式（4-8） 进行格式化，并把数组B中非零元素搜集起来，紧密的排列到U数组中去。这样，从第一行（i=1）做到第N-1行（i=N-1），就形成了完整的因子表。整个计算过 程是一行号为循环变量的。以下将详细讨论图中各细框的工作情况。首先介绍当flag=1时即形成第一因子表时的 工作情况。A部分把工作数组B从i+1到N-1单元全部充零，并把导纳矩阵对角元素的虚部送到 B数组的第i个单元。B部分把导纳矩阵非对角元素的虚部按其列下标送到B数组的相应 单元中去。在flag=1时，程序不执行C部分中的运算，直接转入到D部分。至此，系数矩阵的第i行元素已稀疏的按其列

29、下标排列在B数组中，以下将在D部分 中按行对工作数组（即待消行i）中的元素进行消去运算。一般的说，在形成因子表第i行各元素时，工作数组应该与i-1行以前已形成的各行因子表元素进行消去运算。因此，在D部分中安排了消去行号i_above=1到i_above=i1 的循环过程。该框开始时，对n_u赋值1，为顺序取用因子表上三角矩阵已形成各元素作 好准备。由于因子表及系数矩阵B都具有稀疏的特性，消去过程比较复杂，以下用图所示的例 子来说明。B 33-1B 36-1B 38-1图4-2中因子表的第一行及第二行已经形成。为了清楚起见，在图4-1中已将这两行 因子表元素展开排列，实际上在数组U中它们是密集排

30、列的（见图4-3）。因子表第一行UUU因子表第二行1315-16U U221待消的第三行图4-2形成因子表时的消去过程在图4-2所示的例子中，第一行有三个非对角元素，第二行有三个非对角元素。由式（4-5）可知，图中各元素的具体意义为U =B （i）, U =B （1）,U =B 1313 ,1515 ,2828NUsum1=3 NUsum2=3U 13-U 15-U 16-U2U 28-U 29- TOC o 1-5 h z n_u=123456图4-3因子表上三角矩阵的存放形式工作数组中B33为系数矩阵第三行的对角元素，B6及B8为非对角元素。首先讨论因子表第一行与工作数组的消去过程。根据式

31、（4-7），可以写出 B1）=B-B X B/） 3/3/ 311/如上所述，(4-13)(4-14)式中：下标应该由待消行号3开始，即/=3, 4。因为系数矩阵是对称矩阵， 在因子表中不必保留下三角部分，因此上式中B31还必须利用式（4-9）求出：B31=D- B13（1）=D- U131111得到后，就可以进行消去运算。先从，即对角元素起B33（1）=B33-B31 13=B33-B31 U13 当/=4时：B34（1）=B34-B31 XB14（1）=B34-B31 U14 但在图4-2所示的例子中是零元素，因此上式变为B （1）=B3434这样，式（4-14）所表示的运算对工作数组（即

32、待消行）中的元素没有任何影响。因此，在 第一行与第三行进行消去运算时，只需要从因子表第一行中列下标为3的元素开始，顺次 取以后的元素（见图4-3），按照其列下标与工作数组中相应的元素进行消去运算。在本例 中，除了按式（4-13）进行计算以外，还应进行以下两次运算：B35（1）=B35-B31 X U15a（4-15）B36=B36-B31X U16 J式（4-15）中，B35为零（见图4-2），即系数矩阵中B35为零元素，但与第一行进行消去运算 后变成了非零元素，因而在工作数组（即待消行）中出现了一个注入元素。现在讨论第二行对第三行的消去运算。根据式（4-7），消去过程的计算公式为式中:B3j

33、 =B3k 一B32 Wj (4-16)B32(1)=D- XB23(2)=D- U232222由图4-2可知，在该例中U23=0，因此B32也等于零，这样式(4-16)变为 b3.(2)=b3(i)也就是说，由于U23是零元素，因此第二行不必对第三行进行消去运算。在一般情况下， 当工作数组中待消行号为，而第i (i i的限制，也就是说，式中导纳矩阵中的元素不仅要 考虑上三角矩阵中的元素，还要考虑对角元素及下三角矩阵中的元素，即考虑以下3中情 况： TOC o 1-5 h z JiJi的非对 角元素，而不能直接取出Ji的元素。为了取得J for i=1 to NA for i=1 to Nj=

34、Yijn.jA=-Yiii.BA=Yiii.G图5-1计算节点功率框图当求无功功率时:T Bi=Gj(10)(11)(12)(13)现在我们来讨论计算节点功率的程序框图(见图5-1)。首先，对图中符号作简单说明。NodalPower：存放节点功率的数组，每个节点占两个单元，分别存放P.和Q.。该数 组定义为 NodalPower2003。/lag：是一个标识变量，当进行P-0迭代时，应置1，NodalPowerik即节点i的有 功功率。当进行Q-V迭代时，/lag置2, NodalPowerik即为节点i的无功功率。NodalVoltage：节点电压结构体数组，定义为：struct Nodal

35、Voltage_Typedouble V,theta; NodalVoltage300;V为节点电压幅值，theta为节点电压相角。Yii、Yij、NYseq存放导纳矩阵(详见3章)。图中有两个函数sin和cos，它们分别表示正弦函数和余弦函数。在程序中为了累计各节点功率，在正式计算以前应将功率数组NodalPower的相应单 元清零，图5-1中、两框执行这个运算。当进P-0行迭代时，/lag在主程序中置1， 因此框中即为NodalPoweri1，所以在框中将各节点存放有功功率的单元清零。同样 道理，当进行Q-V迭代时，框中将各节点存放无功功率的单元清零。由于计算节点功率的过程以扫描导纳矩阵为

36、主要途径，因此程序包括两重循环。i循 环的主要作用是控制行号并累计对角元素对该节点(与行号对应)功率的影响。图中、 框与式(5-8)、式(5-9)及式(5-7)相对应。n循环的作用是按行取导纳矩阵非对角元素， 并累计这些元素对有关节点功率的影响。图中(11)框与式(5-12)、式(5-13)相对应。(13)框与 式(5-10)、式(5-11)相对应。由于导纳矩阵的上三角部分最后一行没有非对角元素，因此i循环做到N时，在框图 5-1中框即可转出，不必做下面n循环中的运算。如果不在框中转出即程序中不设 框的判断，就可能使计算发生错误。如前所述，导纳矩阵非对角元素首地址的定义为 NYseqN，当i=

37、N时，在框中NYseqi+1即NYseqN+1将是某一个随机数，因此由- 3框的循环体内运算的循环次数就失去了正确的控制。6.迭代过程电力系统潮流计算的迭代过程包括从送电压初值到求出系统各节点电压的全部计算 过程。本节介绍P-Q解法潮流计算的迭代过程。迭代过程是以送电压初值开始的。在程序中送电压初值包括两个内容，即向系统各PQ 节点送系统的平均电压岭和向各PV节点分别送其应维持的电压值。所有各节点电压角度。的初值都取0。由P-Q分解法的计算步骤可以看出，当求出系统各节点注入功率后，应该与其给定的 功率相比较，以得到各节点的功率误差。然后检查其中最大的功率误差是否满足收敛条件。 当需要继续进行迭

38、代时，可进一步用各节点功率误差构成修正方程式的常数项P/V或 QV图6-1所示框图就是用来在迭代过程中求系统最大功率误差和修正方程式常数项的。在介绍这个框图以前，我们首先归纳一下有关的计算公式。设节点i注入功率的计算值为W.，给定的负荷功率为W”该节点发电机的出力为W.， 则节点i的功率误差可按下式计算： TOC o 1-5 h z HYPERLINK l bookmark263 o Current Document W = Wr+W -W.(6-1)当节点i无负荷时，W=0，式(6-1)变为即1liAWz=W .-W.(6-2)当节点i不是发电机节点时，*广，式(6-1)变为gZAw=Wr-

39、W.(6-3)如果节点i为联络节点，则式(6-1)变为W.=-W.(6-4)以下介绍框图6-1。该程序框图包括求各节点功率误差、保留最大功率误差和求修正 方程式常数项3部分内容。计算各节点功率误差由(13)框完成，(15)(加)框的作用是留最 大功率误差和计算修正方程式的常数项。功率误差是按节点顺序逐点计算的，对每个节点都应判断是否为负荷节点或发电机节 点图中框和框。图中框为计算和判断作好准备，框取出相应节点的电压，用于 求修正方程式的常数项见(18)框。根据式(6-1)(6-4)不难看出3框中的计算内容。需要指出的是在这些框中不但 计算了各节点的功率误差，而且也为迭代收敛是打印输出节点功率作

40、好了准备，其中框 及框的作用将在6节中介绍。图中(15)(19)框执行留最大功率误差和求修正方程式常数项的运算。图中变量MaxError 存放每次迭代的最大功率误差，ErrorNode存放相应的节点号。在迭代过程中可以根据需 要打印出每次迭代时的MaxError及ErrorNode，这样有助于帮助计算人员分析影响迭代收 敛的因素。如前所述，当进行Q-V迭代时，修正方程式中不包括与PV节点有关方程式，因此， PV节点不应参与留最大功率误差和求常数项的运算，图中5框、(19)框对此作了处理。由图中可以看出，第N个节点(即平衡节点)也不参与留最大误差的运算，以避免对 收敛性作出错误的判断。平衡节点参

41、与前半部分的运算，主要是为了求出该节点实际发出 的功率，以备打印输出。最后，对求解线性方程组及修正电压部分作一说明。为此，在图6-2中比较详细的画 出了这部分程序的框图。当不满足收敛条件，即MaxErrorNeps时，需要修正节点电压矢量。在flag二1时，如 前所述，表示正在进行P-0迭代，因此应该用第一因子表对DI数组进行前代和回代的运 算，最后求出VQ6的列矩阵。当flag=2时，表示正在进行Q-V迭代，这时就应该用第二(19)(10)(11)(12)(13)(14)(15)(16)(17)(18)(20)因子表对DI求解，计算出站列矩阵。当修正量VQ6或V求出以后，可以由1到N-1顺序

42、对每个节点的电压进行修正。 当flag二1时，即进行P-0迭代时，由线性方程组求解得到的修正量必须除以系统平均电压 V。，否则当采用有名值进行计算是就要发生错误。图6-2修正电压程序框图7.支路功率计算与输出7.1支路功率计算潮流计算程序的各种不同解法的区别主要表现在如何求出节点电压的过程。当节点电 压求出以后，计算支路功率的方法对各种潮流计算程序都是一致的。首先讨论支路为输电线路的情况。其等值电路如图所示。当两端电压匕、V.已知时， 流经的电流为7（7-1）因此输电线路两端的功率分别为5疽卜（7-2）pWjM 七由于输电线路两端有等值对地电容，所以从节点i、.流入线路的无功功率还应减去线 路

43、的充电功率：Qi=Q-F 匕之、卜（7-3）.mx 匕l当不接地支路为变压器时，其模拟电路如图3-1（a）所示。由于图中理想变压器没有 损耗，所以从节点i、j流入变压器的功率就是流经支路阻抗乙.的功率。因此只要把非标准 变比侧电压折算到标准变比侧：（7-4）（7-5）（7-6）（7-7）澄 =! x v就可以按照式（7-1）、（7-2）类似的公式计算变压器支路的功率：n V- VI = .j 乙ijPV I、jj V Ij J当支路两端功率求出以后，支路的功率损耗可由下式求得:Pij=Pij+Pji Qij=Qj+Qji7.2节点数据的输出首先定义了一个结构体Nodedata_Type，用来存

44、储节点数据，其定义为:Nodedatai.i, Nodedatai.VNodedatai.theta, NodalPoweri1NodalPoweri2, Nodedatai.PNodedatai.Q打印 Nodedatai.i,Nodedatai.V, Nodedatai.theta,NodalPoweri1NodalPoweri2L Nodedatai.P.Nodedatai.Q图7-1打印节点数据框图struct Nodedata_Typeint i;double V,theta,P,Q; Nodedata300;其中i：为节点的编号，V：为节点的电压,theta：为节点电压相角，P：为

45、发电机有功功 率，Q：为发电机无功功率。节点数据的程序框图如图7-1所示，由图7-1中可以看出节点数据的处理是按优化后 节点编号顺序进行的，因此在进行节点数据输出时，必须对节点号进行还原。节点数据输 出的形式如表7-1所示。表中i为节点号，V代表节点电压的模值，0表示节点电压的角度， 鸟、Qt分别代表该节点负荷的有功功率及无功功率，P、Q表示该节点发电机的有功功率 及无功功率。I & IIi I V |0| Pl | Ql | Pg | Qg表7-1节点计算结果形式在框图7-1中除了对每个节点打印表7-1所列的7个数据框以外，还顺次检查性 系统各节点的电压值以寻求系统电压最低值及其节点号图中、

46、框。当节点为发电机节点时，应输出打印发电机的有功功率（P&）及无功功率（Q&），这 些数据是从Genepower数组中取得的图7T中框，该数组定义为GenePower503 o 在迭代过程中发电机的实际功率都已送到了 GenePower数组中，见6节中图6-1的框。各节点负荷功率（联络节点也作为零功率的负荷节点）是从节点功率NodalPower数 组中取得的。NodalPower数组中的数据是图5-1所示程序计算的结果，当节点不是发电机 节点时，其中存放的就是相应节点的实际负荷功率。当节点既是负荷节点又是发电机节点 时，其中存放着该节点发电机功率与负荷功率之差。为了打印输出这种节点的负荷功率，

47、 在图6-1中、框执行了向NodalPower数组中送负荷功率的运算。这样，每次迭代结 束后，NodalPower数组中存放的全都是节点的实际负荷功率。在图7-1中，A部分用来找出节点号为N的节点和节点号为J（原始数据中平衡节点 的节点号）的节点，并将两者的数据进行交换，以实现节点号的还原。B、C两部分是用来实现节点数据打印的。如果要打印全部节点数据，则执行C部分。 C部分中首先查找是否有数据优化增加的节点，如果有就跳过该节点不打印。如果要部分 打印节点数据，则执行B部分。B部分中按照输入的节点数和节点号进行查找和打印。最后将系统最低电压及其节点号Vmin、VminNode打印出来见7. 3节

48、图7-2中（14）框， 为运行、计算人员分析运行方式提供参考数据。7.3支路数据的输出支路数据输出的程序框图如图7-2所示。由图中框可以看出，支路数据的计算和输 出是按优化后的支路数据逐条进行的，因此在支路数据输出时，也必须对支路数据进行还 原。支路数据的输出形式如表7-2所示。表7-2中字母i、j表示支路两端的节点号，Pij、Qij、Pji、Qji分别表示支路首末端的 有功功率及无功功率。illi I j I Pi | Qij | P | Qj表7-2支路数据输出形式现在介绍图7-2。框中PLoss和QLoss存放系统总网损的有功部分和无功部分，为了累计系统的网损， 在框中首先对Ploss和

49、QLoss清零。、框控制循环次数，顺次对原始支路数据扫描一遍。PLoss=0, QLoss=0(10)(11)(12)(13)Vj二NodalVoltagej.V, theta二NodalVoltagej.thetaEj二Vj*cos(theta), Fj二Vj*sin(theta)falseBranchn.i0 orBranchn.j0truefalseBranchn.i0v trueEi二Ei/YK, Fi=Fi/YKYK=0DE二Ei-Ej;DF二Fi-Fj, Zmag2二R*R+X*X;Ir=(DE*R+DF*X)/7mag2Ii=(DF*R-DE*X)/Zmag2Ej二Ej/YK,

50、Fj=Fj/YKBranchPowern.Pij=Ir*Ei+Ii*Fi, BranchPowern.Qij=Ir*Fi-Ii*EiBranchPowern.Pji=-Ir*Ej-Ii*Fj,BranchPo】wern.Qji=Tr*Fj+Ii*Ej;BranchPowern.Qij=BranchPowern.Qij-Vi*Vi*YK/2.0BranchPowern.Qji=BranchPowern.Qji-Vj*Vj*YK/2.0PLoss二PLoss+BranchPowern.Pij+BranchPowern.PjiQLoss二QLoss+BranchPowern.Qij+BranchPo

51、wern.Qji图7-2打印支路数据框图（1）框取支路数据，将支路数据送入中间工作单元。、框将支路两端节点电压取出并化为直角坐标系统。for i=1 to Nb-false-ti=1-trueBranchPoweri.i=N, false，i=2.trueBranchPoweri.j=Nri = 1falsetrueBranchPoweri.i=b. false-ri=2 _trueBranchPoweri.i=b框判断所取支路是否为变压器支路。当所取支路为变压器支路时，在框中进一步 判断非标准变比设在哪一侧，然后折算电压图中框。（10）、（11）框按式（7-1）、（7-2）计算支路功率，（12）框按式（7-3）计入输电线路的充电功率。 当所取支路为变压器时，在框中YK已经清零，因此2框对变压器支路没有影响。（13）框按照式（7-7）累计全系统的网损。图7-2打印支路数据框图（2）在图7-2中，A部分用来分别找出一端节点号等于N和b （原始数据中平衡节点的节 点号）的支路，并对记录下这两个节点号的位置，程序中通过对t数组和尸数组置1和置0 来实现。当ti=1表示该支路的小节点号等于久当ti=2表示该支路的大节点号等于b； 当ri=1表示该支路的小节点号等于N；当ri=2