潮流计算软件设计
电力系统潮流软件设计1. 原始数据的输入目前计算机的速度和计算方法已经使我们能很快得到计算结果，但是上机以前的准备 工作却非常耗费时间，而且也容易出错因此在程序设计时，必须尽可能减轻上机前的准 备工作，尽可能利用计算机代替人工繁琐的工作在这里，原始数据的填写格式是很关键 的一个环节，它与程序使用的方便性和灵活性有着直接的关系原始数据输入格式的设计，主要是从使用的角度出发，原则是简单明了，便于修改 输入格式的简单明了就可以减轻数据填写的工作量，并减少或避免程序使用者在填写数据 时发生错误电力系统潮流计算往往需要进行多种运行方式的调整和比较，因此在数据格 式上考虑计算过程中修改数据的方便性就显得非常重要以下所介绍的本潮流程序中所用到的 6 个信息：N：系统节点的总数；Nb：系统中支路数，即输电线路条数、变压器数的总和；Ng：发电机节点总数；Nl：负荷节点总数；V0：系统平均电压，在迭代过程中，以它作为电压的初值；eps ：迭代收敛所要求的精确度潮流计算程序所需要的原始数据，分别归纳为下几个结构体：支路数据结构体Branch_Type,有5个数据成员，对应于每条支路的5个数据该结 构体数组定义为：struct Branch_Type{int i, j;double R, X, YK;} Branch[400];当支路为输电线路时，这5个数据成员分别表示：i：输电线路一端的节点号；j：输电线路另一端的节点号；R：输电线路的电阻；X：输电线路的电抗；Y0：输电线路充电电容的容纳。
十 Y0/2z—ijY0/2"zj1:K(b)(a)图 1 - 1 网络支路的等值电路 如图1-1 (a)所示当支路为变压器支路时，这5个数据成员分别表示：i:变压器一端的节点号；j变压器另一端的节点号，这两个节点号有一个带有负号，作为变压器支路的标志； R：变压器的电阻；X：变压器的电抗（rt和XT都是归算到变压器标准变比侧的数值）；Y0：变压器的非标准变比（设在节点号为负的一侧）；变压器的模拟电路如图1-1 （b）所示发电机节点和负荷节点数据分别定义为结构体数组，结构体有4个数据成员，其内容 是相同的发电机节点数据定义为：struct Generator_Type{double P, Q;int i;double V;} Generator[50];负荷节点数据定义为：struct Load_Type{double P, Q;int i;double V;} Load[300];对于发电机节点， P、Q 填正号；对于负荷节点， P、Q 填负号 对于发电机节点和负荷节点，若为 PQ 节点，这些数据成员分别表示：P：节点的有功功率；Q：节点的无功功率；i：节点的编号；V该节点正常运行的电压。
若节点为 PV 节点时数据成员分别表示：P：节点的有功功率；Q：节点无功功率的上限；i：节点的编号；V：节点需要维持的电压，负号是PV节点的标志在节点数据输入计算机后，为了提高计算效率，应统计PV节点的总数Npv，并形成 PV节点结构体数组PV节点结构体定义为：struct PVNode_Type{double V;int i;} PVNode[100];每个PV节点有两个数据，第一个数据为PV节点的给定电压V，第二个数据为相应的 节点号i在形成PV结构体数组的同时，把发电机节点或负荷节点数组中V前面负号去掉2. 数据优化在实用潮流计算程序中，对数据的输入次序应不加以限制，这样便于数据的填写和修 改输入以后，在计算机内对数据再进行排队和整理2.1 负荷节点的优化为了简化程序，需要对负荷节点的顺序进行优化，优化后各节点按照节点编号的顺序 进行排列，该部分程序的实现框图为图 2-1图 2 - 1 负荷节点优化框图其中定义了一个中间工作结构体数组LoadN[300],该结构体定义为: struct Load_Type{double P, Q;int i;double V;} LoadN [300]; 用来存放优化后的节点参数。
2.2 发电机节点的优化由于发电机节点的优化与负荷节点优化相似，因此只给出程序框图（图 2-2）图 2 - 2 发电机节点优化框图2.3 对支路数据的排队整理优化为了使支路数据的排列方式适合形成导纳矩阵的上三角部分，整理过后支路数据按以 下次序排列：（1）两端节点号应把小号排在前边，大号排在后边2）各支路按其小节点号的顺序排列 实现这部分的框图为图 2-3首先对框图中的符号做一介绍：for n=l to Nbfalse———,abs(Rranch[n].i)〉abs(Rranch[n].j)truet 二Branch[n].iBranch[n].i二Branch[n].jRranch[n].j=tc=0图 2 - 3 支路5数据优化框图BranchN[400]:为一个中间工作结构体数组，该结构体定义为：struct BranchN_Type{int i, j;double R, X, YK;} BranchN[400]; 用来存放优化后的支路数据m[k]：用来存放小节点号为i支路的另一节点号s[k]:用来存放该支路原先编号该函数首先将原始数据中节点按照（1）的要求优化，然后按照（2）的要求对支路顺 序进行再次优化。
2.4 对支路特殊情况的优化当某些支路的电阻R大于电抗X的2倍时，影响到PQ分解法的收敛性，这时程序将 自动把支路分成两条支路（增加一个节点），满足PQ分解法要求电阻小于电抗的条件由图 2-4 构成的程序来实现这部分功能对框图中变量的说明：a：用来存放原始数据中支路的个数i：为计数变量该函数的实现思想是：当检查出某个支路满足上述情况时，则在该支路中增加一个节 点，节点号为N+1；增加的支路号为Nb+1，原支路仅保留电阻，新支路参数为原支路的 电抗，而对地支路不变图 2 - 4 对支路特殊情况的优化框图2.5 平衡节点的优化为了简化程序，去掉一些判断，要求把平衡节点排在最后，即作为第N个节点，同时 还要求这个节点为负荷节点如果平衡节点没有负荷，则该节点的负荷功率填零这样保证了节点 N 既是发电机节点，又是 PV 节点，又是负荷节点AB'alsot[i]=1DtrueLoad[i].i=NCt[i]=1truefalseseabs(Branch[i].i)二NtrueBranch[i].j二-Nlsefalses[i]==1abs(Branch[i].j)=Nabs(Branch[i].j)=jtruetrues[i]==2truetrueGenerator[i].i=Nfalse-' •-••3■a-bs(Branch[i].i)=jBranch[i].i〈O falsefalsfalseBranch[i].i<0trueBranch[i].i二Branch[i].j<0 ""q 1 •='" = ■false Generator[i].P=0 and''■Generator[i].Q=0'''falseBranch[i].j<0truetruefalse .j!=Nfalse1T trii pJ truea二i, j二Generator「i].itrueGenerator「i].i=jGenerator「a].i=NBranch「i].j二-jBranch「i].i二—NLoad「i].i=jfor i=l to Nltrues「i]=2for i=l to NbBranch「i].i=iBranch「i].i=ifor i=l to Ngt「i]=2for i=l to NgBranch「i].i=NBranch「i].i=N»for i=1 to Nbfor i=1 to Nb图 2 - 5 平衡节点优化框图图 2 - 5 用来实现这一功能。
对框图中变量的说明：j用来存放原始数据中平衡节点的节点号t[i]:为一个标识数组t[i]=1表示第i条支路的小节点号为j, t[i]=2表示第i条支路的 大节点号为 js[i]：也为一个标识数组s[i]=1表示第i条支路的小节点号为N s[i]=2表示第i条支 路的大节点号为 N该函数的设计思想是：先把平衡节点找出（如图 2-5 A 部分），然后与最后一个节点 进行交换（如图 2-5 B、C 部分），同时保持原有支路数据形式不变（如图 2-5 D 部分）1ijj3. 稀疏导纳矩阵的形成3.1 基本公式当电力系统中i、j两点间输电线路的阻抗为z..时，节点i、j之间互导纳为 ij(3-1)Y 1 Y..=-—=ij z ijij式中：y..是阻抗z..的倒数，即输电线串联支路的导纳；Y..是导纳矩阵中i行j列的非对角 ij ij ij元素由于导纳矩阵的对称性，一般Yij=Yji支路 i、 j 对导纳矩阵中 i、 j 两行对角元素的影响可表示为如下的增量(3-2)△Y 二△二丄二yii jj z ijij这里导纳矩阵对角元素Y..和Y..也就是节点i、j的自导纳ii ij0i当节点i连有导纳为Y0•的接地支路时，它对导纳矩阵的影响仅仅使i行对角元素增加 0i如下的分量：△Yii=Y0i(3-3)(a)(b) 变压器的等值电路当i、j两节点间的支路是非标准变比的变压器时［如图3-1 (a)所示］,我们可用n型等 值电路来模拟［见图3-1 (b)所示］,因此i、j之间互导纳可按下式计算：Y..二-丄=-—y..ij Kz K ijiji、 j 节点自导纳分别有如下的增量：(3-5)(3-6)1△Y^——y..ii K 2 ij△Yjj=yjjjj jj3.2 稀疏导纳矩阵的处理电力系统的导纳矩阵不仅具有对称性，而且具有稀疏性。
当 i、j 节点之间没有直接联 系时，导纳矩阵中非对角元素Y厂及Y..应为零由于导纳矩阵的对称性，在计算机中可以只存储其上三角部分或下三角部分，在本程序中，储存导纳矩阵的上三角部分及对角元素，因此，其中每个非对角元素Y.•的下标都满 ij 足 i < j 对角元素按节点编号顺序存放在对角元素结构体数组中：struct Yii_Type{double G, B;} Yii[300];其中G存放对角元素的实部，B存放对角元素的虚部，每个数组的元素个数与系统节点数 相等由于上三角矩阵中非对角元素和系统中不接地支路一一对应，非对角元素的个数等于 网络中不接地支路数Nb为了节约内存和提高计算速度,在计算机内存中只储存非零元素, 把非零非对角元素“挤实”在一起为了识别非对角元素的行号和列号，我们在每个元素 后存放相应的列下标非对角元素结构体数组定义为：struct Yij_Type{double G, B;int j;} Yij[400];按照这样的排列，取一个互导纳，就可以同时把该元素的列号取出来，为了判断该元 素的行号，需要借助于数组NYseq[300]数组NYseq[300]按导纳矩阵行号的顺序存放各行 非对角元素的首地址（事实上，存放的是各行第一个非对角元素在导纳矩阵非零非对角元 素中的顺序号）。
本程序中还定义导纳矩阵中各行非对角元素的个数NYsum[300]一般我们有NYsum[i]= NYseq[i+1]-NYseq[i] （3-7）由于非对角元素是逐行向下排列的，所以就很容易判断出各非对角元素的行号按照 对支路原始数据处理的要求，可以得出支路排列的顺序和导纳矩阵非对角元素的排列顺序 完全一样因此只要顺序取出支路数据，按照式 （3-1）求出倒数取负号之后，连同该支路 的大节点号（即列下标）顺序送入塔数组，就形成了导纳矩阵的上三角部分ij3.3 导纳矩阵形成过程及框图在本程序中，适应P-Q分解法的需要，导纳矩阵分为两步形成第一步只用不接地支路构成导纳矩阵，不考虑接地支路（包括变压器非标准变比）的 影响这里同时形成两个导纳矩阵，即通常意义上的系数矩阵Y和不考虑输电线路电阻的 系统导纳矩阵Y'，以适应BX法的要求，这两个导纳矩阵实际上只是半成品Y主要用来 形成 BX 法所要求的第一个因子表，当该因子表形成后，就在半成品的基础上把接地支路 及变压器非标准变比的影响加进去，形成完整的系统导纳矩阵其中Y'用来形成BX法第 二个因子表，而Y将在整个迭带求解过程及线路潮流计算过程中发挥作用。
只考虑不接地支路构成导纳矩阵的程序框图如图 3-2 所示整个形成过程需要把不接 地支路扫描一遍，对每条不接地支路作两方面的工作首先把阻抗求倒数并取负号后连同 大节点号送到导纳矩阵非对角元素数组Yj （对应于Y）和Yij1 （对应于Y'）中形成非对角 元素，然后把阻抗的倒数累加到该支路两端节点的自导纳上去[见式（3-2）]为了累加形成对角元素，在计算开始时应对数组Yii和Yii1清零[见图3-2中⑴、⑵框] ⑵框中NYsum为临时工作数组，定义为NYsum[300]，在其中累计导纳矩阵各行非对角元 素的个数，因此也需要预先清零由于两个导纳矩阵的结构是相同的，共用一个NYsum数 组1)(2)(3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15)图 3 - 2 形成不接地支路的导纳矩阵对不接地支路的扫描用一个循环语句来完成［图中⑶框］⑷框把支路的有关数据送进中间工作单元，因为支路为变压器节点号可能为负，所以 在这里对节点号取绝对值在⑸框中，把阻抗的倒数即支路导纳放到中间工作单元Gj、Bj中，而把支路电抗的 倒数放到中间工作单元b_j中⑺框判断支路是否为变压器支路。
若为变压器支路，则导纳需除以变压器非标准变比 后再取负号送到导纳矩阵非对角元素数组 Yij 和 Yij1 中；否则直接取负号送到导纳矩阵非 对角元素数组中［⑻框］⑼框向Yij和Yij1数组送列号这样就把支路阻抗数据变成了导纳矩阵的上三角部分 在⑽〜(12)框中根据式(3-2)累计有关节点的自导纳在(13)框中统计小节点号的不接地支路数目，这也就是导纳矩阵上三角部分每行非对角 元素的个数至此，完成了一条不接地支路的处理；当循环由 1 做到 Nb 时形成了只考虑不接地支 路的导纳矩阵14)〜(16)框是由NYsum数组根据式(3-7)形成NYseq数组3.4 追加接地支路的程序框图(1)(2)(3)(4)(5)(6)(7)(8)(9)(10)图 3 - 3 追加对地支路框图(12)(13)(14)追加接地支路包括两部分内容，即追加对地电容支路和考虑变压器非标准变比的影 响，其程序框图如图 3-3 所示整个计算过程需要对支路数据再进行依次扫描，扫描是由一个循环语句来控制［图中⑴ 框］在⑵框中把支路有关数据送入中间工作单元⑵框根据节点号i、j的符号判断所取的 支路是输电线路还是变压器支路当i、j中任一个为负时，为变压器支路，否则为输电线 路。
当所取的支路为输电线路是，转入（12）〜（14）框，向相应的节点累计自导纳部分当所取支路为变压器支路时，转入⑷〜（11）框在⑷框中判断非标准变比设在支路的哪 一侧如1章中所述，当i<0时，非标准变比就设在i侧，否则设在j侧由图3-1 （b）1及式（3-5）可知，在非标准变比侧自导纳应累计丄y，但在形成不接地支路的导纳矩阵时，K 2 ij该点自导纳累计了雪，因此需要再追加累计（丄-1）蓦=（1-丄）Y图中⑹〜⑼框就是K K K K ij完成这些运算的在非标准变比侧自导纳应累计y.,在形成不接地支路的导纳矩阵时，该 ij点自导纳同样累计了耳，因此需要再追加累计（K-1）县= （1-K）Y，在⑽框和（11）框中完成 K K .,这些运算这样，顺次把支路数据扫描、处理一遍，就形成了描述网络的完整的导纳矩阵4. 稀疏系数矩阵线性方程式的求解4.1 修正方程式的解法及计算公式在 P - Q 分解法潮流计算的迭代过程中，需要反复求解修正方程式△P/V= B' △OV (4-1)△Q/V= B " △V (4-2)如前所述，这两个方程的系数矩阵(B'和B")在迭代过程中保持不变，只要求对不断变化 的常数项(即误差项Ap/v、Aq/v)求解出相应的修正量AOV及在这种情况下，可 以先将系数矩阵进行三角分解，然后只要用分解出的三角矩阵(或因子表)对不同常数项 进行前代及回代的运算，即可得到要求的修正量。
系数矩阵的三角分解可以利用递推公式求得，也可以利用高斯消去法求得，这两种方 法在运算量及内存量上都是等效的本程序利用高斯消去法对系数矩阵进行三角分解并形 成因子表的计算方法在P-Q分解法潮流程序中，为了在迭代过程中轮流求解式(4-1)及式(4-2)需要形成两 个因子表,为此，可以把式(4-1)及式(4-2)统一为如下的形式：B^X二△/ (4-3)在P-O迭代时，式中B即B'，△ 为AOV, △/为AP/V；在Q-V迭代时，式中B为 B"，△为△#, △/ 为△Q/V以下简单归纳以下形成因子表及常数项进行前代及回代运算的有关公式，对于式(4-3) 的系数矩阵进行B三角分解以后，可以得到以下形式的因子表：B (1)12B (1)13B (1)14… B (1)1,n-11B (2)231B (2)24… B (2)2,n-1B(1)22B (1)32B (3)34...B ⑶3,n-1B (2)33B (1) n-1,2B (2) n-1,3B (3) n -1,41B (n-2)n -1, n -1111B21B31Bn-1,1(4-4)因子表中上三角部分元素组成了上三角矩阵U「1 U U U …U12 13141,n-11UU… U23242,n-11U… U343,n-11Un-2,n-11其中元素(4-5)Uij=Bij(i) (i 这 i 次运算中包括 i-1 次消去运算及一次规格化运算：ijB (k)=B (k-1) -B (k-1)B (k)=B (k-1) -B (k-1)Uij ij ik kj ij ik kj(k=1, 2, …, i-1; j=k+1, k+2, …, n-1) (4-7)B (i -1)B..(i)= Q (j=i+l, i+2,…，n-1) (4-8)ij B (i-1)ii由于系数矩阵B为对称矩阵，在因子表中不需要保留其下三角部分在形成因子表的 过程中需要用到下三角部分的元素Bik(i)伙＜i)时，可以按下式由上三角矩阵相应的元素求 出：Bik(k-1)=Bki(k)Bkk(k-1)=Uki1kk(4-9)利用式(4-6)〜(4-9)就可以由对称系数矩阵B有规律的计算出因子表的所有元素 对式(4-3)来说，当系数矩阵为对称矩阵时，其具体计算公式为：△l.(i)二△"-1)-U.•△l.(i-1)(j=2, 3, •…n-1; i=1, 2,…,j-1) (4-10)j j ij i△l(j)=D.Al(j-1) (4-11)显然，顺次取因子表各元素参加一次运算就完成了前代过程的计算。 回代过程计算公 式为：(4-12)△x.二△")-灯 U □ x (j=n-1, n-2,…，1) i i ij jj=i+14.2 因子表形成程序框图由于式(4-3)中的系数矩阵和导纳矩阵 B 具有相同的结构和稀疏性，因而在一般情况 下因子表的上三角矩阵U也是稀疏的但由于在消去过程中可能会增加一些新的非零元素, 所以U和B的结构不完全相同为了节约内存和提高计算速度，在因子表中将不保留零元素，而把上三角矩阵中非零 元素紧凑的排列在一起，对每个上三角矩阵U的非零元素都附一个列下标，并对其中每一 行都给出非零元素的数目以下框图(图 4-1)来说明形成因子表的具体过程 首先对框图中符号作一说明flag：是一个标识变量，当flag=1时，该框图所表示的程序形成第一个因子表(与B'相 对应)当flag二2时形成第二个因子表(与B"相对应泊flag在该子程序调用前根据主程序F⑴count=1(1)⑵⑵⑶⑶B「count]=0B[i]二Yii[i].Btrue flag=2 and i=i_pv -true'■ GQunt〉Nusum「i_abo.ve卜計 false■Fa 1 qqtrue■ oount>NUsum[i_above]-falsetrue"■ U[n_u].j=i _ "4 falsefalsetrueB[j]=0j二U「n_u].j, B[j]=B[j]-Btemp*U[n_u].value oount=oQunt+1, n u=n u+1U[n_u].value=B[j]*Btemp, U[n_u].j=j oount=oount+1, n u二n u+1Btemp=1.0/B「i], D「i]二Btempoount=0n_pv=n_pv+1 i_pv二PVNode「n_pv].i *(nusum+i)=O;D「i]=Ofor count二NYseq「i] to NYseq「i+1]Tj二Yij「oount].j. B「j]二Yij「oount].Bi二PVNode「oount].i, B「j]=On_pv=1i pv二PvNode「1].i)=U「n u].value/D「i above]for count二i+1 to N—1oount=oount+1. n u二n u+1for i=1 to N—1truefalseflag=2for oount=1 to Npvfor j二i+1 to NTn u=1. i above=1i above=i above+1Nusum「i]=ooun1图 4 — 1 形成因子表框图要求置 1或置 2，作为参数传个子程序。 n_pv：为PV节点数组的计数变量，在形成第二个因子表时，用它可以判断系数矩阵 中应该去掉哪些列和哪些行i：为因子表正在形成的行号变量j为列下标变量n_u：是因子表上三角矩阵元素计数变量i_above:为消去行号计数变量（是被消行号计数变量，因子表按行消去过程中，依次 取行上面的到行）i_pv：为PV节点的节点号变量Btemp:为临时变量count:为临时计数变量形成因子表所需要的原始数据可由以下数组取得，这些数组的定义及内容详见 1 章和 3 章：Load :负荷功率数组PVNode： PV节点数组Yj：导纳矩阵非对角元素数组NYseq:导纳矩阵各行非对角元素的首地址数组Yii：导纳矩阵对角元素数组形成的因子表将放在以下数组中：NUsum :存放因子表上三角矩阵各行非对角元素数D：存放因子表的对角元素U：存放因子表上三角矩阵元素，定义为结构体数组：struct U_Type{double value;int j;} U1[300],U2[300];其中存放该元素的数值，存放该元素的列下标最后，在框图中还有一个很重要的工作数组B[N],形成因子表的运算主要在这个数组 中进行现在分别介绍框图中各个部分的作用。 图 4-1 中 A、B、C 三部分的作用是“传递”，通过这三部分的工作把系数矩阵中待消 行（i行）的元素按其下标稀疏的排列在工作数组B中D框的作用是把工作数组B中的 待消行元素按照式（4-7）进行消去运算最后通过E框的工作把数组B中的元素按式（4-8） 进行格式化，并把数组B中非零元素搜集起来，紧密的排列到U数组中去这样，从第一行（i=l）做到第N-1行（i=N-l）,就形成了完整的因子表整个计算过 程是一行号为循环变量的以下将详细讨论图中各细框的工作情况首先介绍当flag=1时即形成第一因子表时的 工作情况A部分把工作数组B从i+1到N-1单元全部充零，并把导纳矩阵对角元素的虚部送到 B数组的第i个单元B部分把导纳矩阵非对角元素的虚部按其列下标送到B数组的相应 单元中去在flag=1时，程序不执行C部分中的运算，直接转入到D部分至此，系数矩阵的第i行元素已稀疏的按其列下标排列在B数组中，以下将在D部分 中按行对工作数组（即待消行i）中的元素进行消去运算一般的说，在形成因子表第 i 行各元素时，工作数组应该与 i-1 行以前已形成的各行因子表元素进行消去运算因此，在 D 部分中安排了消去行号 i_above=1 到 i_above=i-1 的循环过程。 该框开始时，对 n_u 赋值 1，为顺序取用因子表上三角矩阵已形成各元素作 好准备由于因子表及系数矩阵 B 都具有稀疏的特性，消去过程比较复杂，以下用图所示的例 子来说明BBB图 4 - 2 中因子表的第一行及第二行已经形成为了清楚起见，在图 4 - 1 中已将这两行 因子表元素展开排列，实际上在数组 U 中它们是密集排列的(见图 4-3)因子表第一行 ►UUU因子表第二行 ► 13-U―24 15 16U28U——29待消的第三行 ►图 4-2 形成因子表时的消去过程在图 4-2 所示的例子中，第一行有三个非对角元素，第二行有三个非对角元素由式 (4-5)可知，图中各元素的具体意义为U13=B13(1), U15=B15(1), „, U28=B28(2)NUsum[l]=3 NUsum[2]=3十1315162^28十29n_u=1 2 3 4 5 6图 4-3 因子表上三角矩阵的存放形式工作数组中b33为系数矩阵第三行的对角元素，b36及b38为非对角元素首先讨论因子表第一行与工作数组的消去过程根据式(4-7)，可以写出B3j(1)=B3j-B31 X Bij⑴式中：下标应该由待消行号3开始，即j=3, 4…。 因为系数矩阵是对称矩阵，如上所述, 在因子表中不必保留下三角部分，因此上式中B31还必须利用式(4-9)求出：B =丄 B (1)= — U(4-13)(4-14)31 D 13 D 1311 11得到后，就可以进行消去运算先从，即对角元素起B (1)=B -B XB (1)=B -B XU当 j=4 时：33 33 31 13 33 31 13B34(1)=B34-B31XB14(1)=B34-B31XU14但在图 4-2所示的例子中是零元素，因此上式变为B34(1)=B34这样，式(4-14)所表示的运算对工作数组(即待消行)中的元素没有任何影响因此，在 第一行与第三行进行消去运算时，只需要从因子表第一行中列下标为3的元素开始，顺次 取以后的元素(见图4-3)，按照其列下标与工作数组中相应的元素进行消去运算在本例中，除了按式(4-13)进行计算以外，还应进行以下两次运算：B (1)=B -B x u “35 35 31 15a (4-15)B (i)=B —B X U >36 36 31 16式(4-15)中，B35为零(见图4-2)，即系数矩阵中B35为零元素，但与第一行进行消去运算 后变成了非零元素，因而在工作数组(即待消行)中出现了一个注入元素。 现在讨论第二行对第三行的消去运算根据式(4-7)，消去过程的计算公式为式中：B3j(2)=B3k ⑴-B32⑴'Bj ⑵(4-16)B32(1)= D ' B23⑵讣'U2322 22由图4-2可知，在该例中U23=0,因此B32(i)也等于零，这样式(4-16)变为B3j(2)=B3j(1)也就是说，由于u23是零元素，因此第二行不必对第三行进行消去运算在一般情况下， 当工作数组中待消行号为i,而第i'(/y)行中没有下标为i的元素时，即可省去i'行对i行 的消去运算过程在图4-1的D部分中，为了判断i'行是否需要对工作数组中的待消行进行消去运算， 必须对因子表i'行中元素逐个检查其列下标是否等于待消行的行号i［D部分中⑴］当查到 列下标等于i的元素时，首先按照式(4-9)求出Bii' (i'-1)［D部分中⑵］,然后用该行剩余元素 (包括列下标为i的元素)按列下标对工作数组中相应的元素进行消去运算［D部分中⑶］如果在因子表的i'行中没有列下标为i的元素，那么D部分中⑴一直检查到底，直接转入 i'+1 行对 i 行的消去过程图4-1中E部分的作用是最终形成因子表的第i行首先形成因子表的对角元素［E部 分中⑴］ ，然后按式(4-8)对工作中的非零元素进行规格化计算，并把这些元素及其列下标 顺序排列到因子表上三角矩阵数组U中。 E部分中⑵、⑶累计并形成因子表第i行非零非 对角元素的个数至此因子表的第i行元素全部形成以下讨论形成第二个因子表的特点由 P-Q 分解法的基本原理可知，第二个因子表是在 Q-V 迭代过程中求解修正方程式 (4-2)用的当系统中有r个PV节点时，式(4-2)的系数矩阵就要去掉与PV节点有关的r 行和r列，这样系数矩阵就降为N-r阶因此第二个因子表的上三角矩阵实质上也就是N-r 阶在图4-1中，F部分的作用是在形成第二个因子表时跳过与PV节点有关的行，C框的 作用是去掉与PV节点有关的列，这样第二个因子表的上三角矩阵实际上就降为N-r阶关于形成因子表的过程就介绍到这里4.3 线性方程组的求解过程与框图对系数矩阵进行三角分解并形成因子表后，即可利用式(4-10)〜(4-12)对线性方程组 不同的常数项进行前代和回代运算，从而求得相应的解以下举例说明程序框图 4-5 是如何实现对线性方程组进行求解的，在这里应特别注意 因子表稀疏性的影响设有五阶线性方程组，其因子表结构如图 4-4所示，因子表上三角矩阵元素是按行存 放的，所以在应用时按行取用也比较方便在前代过程中先取因子表的第一行元素进行运一 D 一~UUU11121315DUU222325D=D，U=UU3334 35DU4445D551—图 4-4 因子表结构举例算，即令i=1，这样式(4-10)变为：△l.(1)=^I.-U1jA/1 j=2, 3, 4, 5)具体计算为:△⑴二△I4-UQI1二△△4⑴二△I5-U/I1显然，由于本例中U14=0 （见图4-4）,因此对厶行就可以省去以上的运算。 因此以上的运 算只有3次（即j=2, 3, 5）是有具体意义的在用△/］对其他常数项进行计算以后，即可按 式（4-11）进行规格化运算：△mu(1)(2)(3)(4)(5)(6)(7)(8)(9)(10)(11)图 4-5 修正方程式直接解程序框图然后，再取因子表上三角矩阵第二行元素进行运算，即令i=2，这时式（4-10）变为△l.（2尸△/.⑴―U2•△厶⑴ （j=3, 4, 5）j j 2j 2具体计算为△厶⑵二乂⑴-⑴△i4（2）=Ai4（i）-u24Ai2（1）=Ai4 ⑴△i5⑵二△厶⑴-空乂 ⑴由于本例中U24=0 （见图4-4）,因此对厶行又可以省去以上的运算因此，以上的运 算只有两次（即 j=3, 5）是有实际意义的进行以上运算以后，又可以对第二个常数项△/$）进行规格化运算： △I2（2）=D22^I2（1）由以上的讨论，我们可以对前代过程的运算归纳以下的规律：（1） 顺次取因子表上三角矩阵中的第i行（由1到N-1）元素参加运算，每行的运算 次数等于因子表中该行非零非对角元素的数目2） 参加运算的元素由所取因子表元素的行号（i）及列号（j）来确定，其计算内容 为△i=Ai-u..Ai.（3） 当因子表i行元素运算结束以后，即可对常数项第i个元素进行规格化运算：△i=d.Ai.i ii i根据以上规律，不难看出图中框⑴〜⑹的意义。 图中⑺〜（11）框执行回代运算，程序结构和⑴〜⑹框相似，但有两点不同：（1） 回代过程运算是从后向前做的，因此⑺框的循环是从N-1到1同时，因子表上 三角矩阵元素也是倒取的，它的地仍由 n_u 控制在前代过程中 U 矩阵中的元素是顺取的， 前代过程结束后， n_u 所记的地址恰好超过了最后一个元素的位置因此在回代过程中， 在取U矩阵元素时，只要每次向前移一个单元［见图4-5中⑽框］,就可以得到计算所要求 的元素2） 回代过程计算公式的内容和前代过程不同：△I.=AI.-U.AI.i i ij j 不难看出，前代过程实际上是按列消去过程，回代过程则是按行进行计算的回代结束后△/数组中的内容就是线性方程组的解5. 迭代过程中节点功率的计算5.1 基本公式分析在迭代过程中需要不断根据迭代出来的电压值计算各节点的功率，然后与给定的功率 比较，以判断是否收敛当不收敛时，还要以两者的差值求出修正方程式的常数项，以便 进行迭代计算节点功率计算量在整个迭代过程中占很大比重，因此，节点功率的计算速 度对整个 P-Q 分解法潮流程序的效率有较大的影响节点功率的计算公式在电压的极坐标系统中应为Pi=Vi 工V. (G cos0 + B sin0 )、i i j (/ (/ (/ ijj©卜 (5-1)Qi=Vi 工V. (G cos0 - B sin0 )丿i i / i/ i/ i/ i//©i式中：jg i表示工 后面的j点应与i节点相连，并包括j=i的情况。 必须注意，在上式中，并没有 />i 的限制，也就是说，式中导纳矩阵中的元素不仅要 考虑上三角矩阵中的元素，还要考虑对角元素及下三角矩阵中的元素，即考虑以下3中情 况：(1 )j>i(2) j 由式(5-2)可以看出，导纳矩阵对角元素G..、B..只对节点i的注入功率有影响:ii ii(5-3)△Pi' = V2 x Gii△Qi' = V2 x(-Bii 丿」式中：△P.'、△Q.'分别为对角元素G..、B„对节点i注入功率的影响因此，在求从节点 i i ii ii1 到节点 n 各节点的整个过程中，导纳矩阵对角元素只起一次作用但是，导纳矩阵的非对角元素G..及B..不仅在计算节点i的注入功率时起作用：ij ij△P." =V.V.(G.cos3..+B.sin3..)i i j ij ij ij ij(5-4)△Q." =V.V.(G..sin3-B.cos3..) Ji i j ij ij ij ij 而且，在求j点注入功率时也起作用事实上，由式(5-2)可知△P.” =VV.(G.cos0..+B.sin0..)[(5-5)由于导纳矩阵的对称性：△Qj" 評叫厂弩叫)Jj J l' j JI JI ji又因Gji=Gij,Bji=-Bijcose =cose ,sine =-sineji ij ji ij因此，式(5-5)可以改写为：△p." =V.V.(G.cos3-B.sin3..) >(5-6)i i j ij ij ij ij△Qi" zyQfinOjB 严e))等号式(5-6)与式(5-4)除了等号右边某些项的符号不同外，其结构完全一样。 通过以上的讨论，我们在程序中，不必按节点一个一个分别求其注入功率，而是把全 部节点注入功率的计算包括在一个统一的积累过程中为此，可以顺序取导纳矩阵中的非 对角元素，分别按照式(5-4)及式(5-6)同时积累i点与j点的注入功率，取对角元素按照 式(5-3)只积累相应点的节点功率这样，对导纳矩阵所有元素扫描一遍，就得到了个节 点的注入功率5.2 节点功率的计算过程及框图P-Q分解法的主要特点是P-e和Q-V分别轮流进行迭代，因此，求节点功率的程序 也要求轮流计算节点有功功率及节点无功功率由式(5-3)、式(5-4)、式(5-6)可以看出，有功功率的积累公式和无功功率的公式在结构上是一样的，为了简化程序可以把它们统一为以下的公式对角元素对节点功率的影响：当求△Pi'时，式中：△W' =V.2 x A.iii(5-7)当求△ Qi'时，式中：Ai=Gii[i]i ii(5 -8)Ai= -Bii[i](5-9)非对角元素对节点功率的影响：△W"=VV(Acose+Bsine)i i j i ij i ij(5-10)当求有功功率时：△W'" =VV(A cose -Bsine )j i j i ij i ij(5 -11)A =G ，B =B(5 -12)ii i ij(1)(2)(3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)for i=1 to NNodalPower[i][flag]=0for i=1 to NVi二NodalVoltage[i].VA二Yij[n].G, B=Yij[n].B、flag=ltrueVVA=Yii[i].GA=—Yii[i].B< VNodalPower[i][flag]=NodalPower[i][flag]+Vi*Vi*Aflag=1*j=Yij[n].jVV=Vi*NodalVoltage[j].Vtheta二NodalVoltage[i].theta-NodalVoltage[j].thetaA=A*VV*cos(theta)B二B*W*sin( theta)NodalPower[i][flag]二NodalPower[i][flag]+A+BNodalPower[j][flag]二NodalPower[j][flag]+A—B当求无功功率时：图 5 — 1 计算节点功率框图Ai=Bij， Bi=Giji ij i ij(5—13)现在我们来讨论计算节点功率的程序框图(见图 5-1)。 首先，对图中符号作简单说明NodalPower：存放节点功率的数组，每个节点占两个单元，分别存放P.和Q.该数 ii 组定义为 NodalPower[200][3]flag：是一个标识变量，当进行P-0迭代时，应置1，NodalPower[i][k]即节点i的有 功功率当进行Q-V迭代时，flag置2, NodalPower[i][k]即为节点i的无功功率NodalVoltage：节点电压结构体数组，定义为：struct NodalVoltage_Type{double V,theta;} NodalVoltage[300];V为节点电压幅值，theta为节点电压相角Yii、Yij、NYseq 存放导纳矩阵(详见3 章)图中有两个函数sin和cos,它们分别表示正弦函数和余弦函数在程序中为了累计各节点功率，在正式计算以前应将功率数组 NodalPower 的相应单 元清零，图5-1中⑴、⑵两框执行这个运算当进P-0行迭代时，flag在主程序中置1, 因此⑵框中即为NodalPower[i][1]，所以在⑵框中将各节点存放有功功率的单元清零同样 道理，当进行Q-V迭代时，⑵框中将各节点存放无功功率的单元清零。 由于计算节点功率的过程以扫描导纳矩阵为主要途径，因此程序包括两重循环 i 循 环的主要作用是控制行号并累计对角元素对该节点(与行号对应)功率的影响图中⑹、 ⑺框与式(5-8)、式(5-9)及式(5-7)相对应n循环的作用是按行取导纳矩阵非对角元素， 并累计这些元素对有关节点功率的影响图中(11)框与式(5-12)、式(5-13)相对应13)框与 式(5-10)、式(5-11)相对应由于导纳矩阵的上三角部分最后一行没有非对角元素，因此i循环做到N时，在框图 5-1中⑻框即可转出，不必做下面n循环中的运算如果不在⑻框中转出[即程序中不设⑻ 框]的判断]，就可能使计算发生错误如前所述，导纳矩阵非对角元素首地址的定义为 NYseq[N]，当i=N时，在⑼框中NYseq[i+1]即NYseq[N+1]将是某一个随机数，因此由⑽- 3框的循环体内运算的循环次数就失去了正确的控制6. 迭代过程电力系统潮流计算的迭代过程包括从送电压初值到求出系统各节点电压的全部计算 过程本节介绍P-Q解法潮流计算的迭代过程迭代过程是以送电压初值开始的在程序中送电压初值包括两个内容，即向系统各 PQ 节点送系统的平均电压V0和向各PV节点分别送其应维持的电压值。 所有各节点电压角度0的初值都取0°由 P-Q 分解法的计算步骤可以看出，当求出系统各节点注入功率后，应该与其给定的 功率相比较，以得到各节点的功率误差然后检查其中最大的功率误差是否满足收敛条件 当需要继续进行迭代时，可进一步用各节点功率误差构成修正方程式的常数项AP/V或厶 Q/V图 6-1 所示框图就是用来在迭代过程中求系统最大功率误差和修正方程式常数项的 在介绍这个框图以前，我们首先归纳一下有关的计算公式设节点i注入功率的计算值为W.,给定的负荷功率为W},该节点发电机的出力为W.,i li gi则节点 i 的功率误差可按下式计算：△W=W +W -W (6-1)i }i gi i当节点i无负荷时，W}=0,式(6-1)变为}i△Wi=Wgi-Wi (6-2)当节点i不是发电机节点时，Wg=0,式(6-1)变为'△W=W -W (6-3)i }i i如果节点i为联络节点，则式(6-1)变为△Wi=-Wi (6-4)以下介绍框图 6-1该程序框图包括求各节点功率误差、保留最大功率误差和求修正 方程式常数项3部分内容计算各节点功率误差由⑶〜(13)框完成，(15)〜⑶)框的作用是留最 大功率误差和计算修正方程式的常数项。 功率误差是按节点顺序逐点计算的，对每个节点都应判断是否为负荷节点或发电机节 点［图中⑶框和⑼框］图中⑴框为计算和判断作好准备，⑵框取出相应节点的电压，用于 求修正方程式的常数项［见(18)框］根据式(6T)〜(6-4)不难看出⑶〜3框中的计算内容需要指出的是在这些框中不但 计算了各节点的功率误差，而且也为迭代收敛是打印输出节点功率作好了准备，其中⑺框 及⑽框的作用将在6节中介绍图中(15)〜(19)框执行留最大功率误差和求修正方程式常数项的运算图中变量MaxError 存放每次迭代的最大功率误差， ErrorNode 存放相应的节点号在迭代过程中可以根据需 要打印出每次迭代时的MaxError及ErrorNode,这样有助于帮助计算人员分析影响迭代收 敛的因素如前所述，当进行Q-V迭代时，修正。
