直角坐标系下的牛顿法潮流计算.doc
Four short words sum up what has lifted most successful individuals above the crowd: a little bit more.-author-date直角坐标系下的牛顿法潮流计算直角坐标系下的牛顿法潮流计算直角坐标系下牛顿法潮流计算1 计算原理1.1 节点导纳矩阵的形成在图1（a）的简单电力系统中，若略去变压器的励磁功率和线路电容，负荷用阻抗表示，便可以得到一个有5个节点（包括零电位点）和7条支路的等值网络，如图1（b）所示。将接于节点1和4的电势源和阻抗的串联组合变换成等值的电流源和导纳的并联组合，变得到图（c）的等值网络，其中和分别称为节点1和4的注 入电流源。 图1 电力系统及其网络以零电位点作为计算节点电压的参考点，根据基尔霍夫定律，可以写出4个独立节点的电流平衡方程如下 （1）上述方程组经过整理可以写成 （2）式中，；。一般的，对于有个独立节点的网络，可以列写个节点方程 （3）也可以用矩阵写成 （4）或缩写为 （5）矩阵称为节点导纳矩阵。它的对角线元素称为节点的自导纳，其值等于接于节点的所有支路导纳之和。非对角线元素称为节点、 间的互导纳，它等于直接接于节点、间的支路导纳的负值。若节点、间不存在直接支路，则有。由此可知节点导纳矩阵是一个稀疏的对称矩阵。1.2 直角坐标系下牛顿法潮流计算 采用直角坐标时，节点电压可表示为导纳矩阵元素则表示为将上述表示式代入的右端，展开并分出实部和虚部，便得 （6）假定系统中的第1，2，3···，m号节点为PQ节点，第i个节点的给定功率设为和，对对该节点可列写方程 （i=1,2，···，m） （7）假定系统中的第m+1，m+2，···，n-1号节点为PV节点，则对其中每一个节点可以列写方程 （i=m+1，m+2，···，n-1） （8）第n号节点为平衡点，其电压是给定的，故不参加迭代。式（7）和式（8）总共包含了2（n-1）个方程，待求的变量有也是2（n-1）个。我们还可看到，方程（7）和式（8）已经具备了方程组的形式。因此，不难写出如下的修正方程式 （9）式中上述方程中雅克比矩阵的各元素，可以对式（7）和式（8）求偏导数获得。当时 （10）当时 （11）修正方程式（11-48）还可以写成分块矩阵的形式 （12）式中，和都是二维列向量；是介方阵。对于PQ节点 （13）对于PV节点 （14）从表达式（1-7）（1-11）可以看到，雅克比矩阵有以下特点：（1） 雅克比矩阵各元素都是节点电压的函数，它们的数值将在迭代过程中不断的改变。（2） 雅克比矩阵的子块中的元素的表达式只用到导纳矩阵中的对应元素。若，则必有。因此，式（1-9）式中分块形式的雅克比矩阵同节点导纳矩阵一样稀疏，修正方程的求解同样可以用稀疏矩阵的求解技巧。（3） 无论在式（1-6）或式（1-9）中雅克比矩阵的元素或子块都不具有对称性。用牛顿-拉夫逊法计算潮流的流程：首先要输入网络的原始数据以及各节点的给定值并形成节点导纳矩阵。输入节点电压初值和，置迭代计数k=0。然后开始进入牛顿法的迭代过程。在进行第k+1次迭代时，其计算步骤如下：（1） 按上一次迭代计算出的节点电压值和，利用式（7）和式（8）计算各类节点的不平衡量、和。（2） 按条件校验收敛，即< 如果收敛，迭代到此结束，转入计算各线路潮流和平衡节点的功率，并打印输出计算结果。不收敛则继续计算。（3）利用式（10）和式（11）计算雅克比矩阵的各元素。（4）解修正方程式（7）求节点电压的修正量和。（5）修正各节点的电压（6）迭代计数加1，返回第一步继续迭代过程。2 计算过程21 计算程序框图图2 潮流计算程序框图2.2 节点导纳矩阵的形成图3 节点导纳矩阵根据第一章节点导纳矩阵的形成原理，设计题目可化为如图含有4节点的形式。采用直角坐标时，节点电压可表示为导纳矩阵元素则表示为/*计算导纳矩阵*G11=z12r/(z12r*z12r+z12m*z12m)+k*k*z13r/(z13r*z13r+z13m*z13m)+z14r/(z14r*z14r+z14m*z14m);B11=-z12m/(z12r*z12r+z12m*z12m)-k*k*z13m/(z13r*z13r+z13m*z13m)-z14m/(z14r*z14r+z14m*z14m)+y140+y120;G22=z12r/(z12r*z12r+z12m*z12m)+z24r/(z24r*z24r+z24m*z24m);B22=-z12m/(z12r*z12r+z12m*z12m)-z24m/(z24r*z24r+z24m*z24m)+y240+y120;G33=z13r/(z13r*z13r+z13m*z13m);B33=-z13m/(z13r*z13r+z13m*z13m);G44=z14r/(z14r*z14r+z14m*z14m)+z24r/(z24r*z24r+z24m*z24m);B44=-z14m/(z14r*z14r+z14m*z14m)-z24m/(z24r*z24r+z24m*z24m)+y240+y140;G12=G21=-z12r/(z12r*z12r+z12m*z12m);B12=B21=z12m/(z12r*z12r+z12m*z12m);G13=G31=-k*z13r/(z13r*z13r+z13m*z13m);B13=B31=k*z13m/(z13r*z13r+z13m*z13m);G14=G41=-z14r/(z14r*z14r+z14m*z14m);B14=B41=z14m/(z14r*z14r+z14m*z14m);G23=G32=0.0;B23=B32=0.0;G24=G42=-z24r/(z24r*z24r+z24m*z24m);B24=B42=z24m/(z24r*z24r+z24m*z24m);G34=G43=0.0;B34=B43=0.0;for(i=1;i<5;i+)for(j=1;j<5;j+)printf("%f+%fj",Gij,Bij);printf(" ");printf("n");/形成节点导纳矩阵将上述表示式代入的右端，展开并分出实部和虚部，便得2.3 计算各节点不平衡量假定系统中的第1，2，3···，m号节点为PQ节点，第i个节点的给定功率设为和，对对该节点可列写方程 （i=1,2，···，m） 假定系统中的第m+1，m+2，···，n-1号节点为PV节点，则对其中每一个节点可以列写方程 （i=m+1，m+2，···，n-1） 第n号节点为平衡点，其电压是给定的，故不参加迭代，其计算程序如下：/计算各节点不平衡量-loop1:printf("迭代次数k1=%dn",k1);for (i=1;i<5;i+)float a=0,b=0;for(j=1;j<5;j+)a+=Gij*ej-Bij*fj; b+=Gij*fj+Bij*ej;Pi=Psi-(ei*a+fi*b);/计算有功功率的增量Qi=Qsi-(fi*a-ei*b);/计算无功功率的增量V32=V3s*V3s-e3*e3;printf("有功功率增量P1=%f",P1); printf(" ,");printf("有功功率增量P2=%f",P2); printf(" ,");printf("有功功率增量P3=%f",P3);printf("无功功率增量Q1=%f",Q1); printf(" ,");printf("无功功率增量Q2=%f",Q2); printf(" ,");printf("电压增量V32=%f",V32);printf("n");2.4 雅克比矩阵计算上述方程中雅克比矩阵的各元素，可以对计算各点不平衡量得公式中求偏导数获得。当时当时以下为程序，只列出j=1的情况/*形成雅克比矩阵*for(j=1;j<4;j+)if(1=j)float c=0,d=0;int m;for(m=1;m<5;m+)c+=G1m*em-B1m*fm; d+=G1m*fm+B1m*em;J1*N-1j*N-1=-c-G1j*e1-B1j*f1;J1*N-1j*N=-d+B1j*e1-G1j*f1;J1*Nj*N-1=d+B1j*e1-G1j*f1;J1*Nj*N=-c+G1j*e1+B1j*f1;elseJ1*N-1j*N-1=-G1j*e1-B1j*f1; J1*Nj*N=G1j*e1-B1j*f1; J1*N-1j*N=B1j*e1-G1j*f1; J1*Nj*N-1=B1j*e1-G1j*f1;2.5 LU分解法求修正方程LU分解,又称Gauss消去法,可把任意方阵分解成下三角矩阵的基本变换形式(行交换)和上三角矩阵的乘积。其数学表达式为:A=LU。其中L为下三角矩阵的基本变换形式,U为上三角矩阵。若有矩阵Ax=b 把矩阵LU分解，求AX=b的问题就等价于求出A=LU后：因为Ly=b可求y，再因为Ux=y，可求出x。原始的求法x=A（-1）*b,某些情况下，如果矩阵A中的数非常小，我认为不是因为大数除以小数误差大么，1/A算出的误差会很大。但LU可以把A分解成两个都比A大的矩阵的乘积，1/L的误差比1/A小的多。求修正方程的程序如下：/*计算修正方程*for(i=1;i<M;i+)Lii=1;for(i=1;i<M;i+)U1i=J1i;Li1=Ji1/U11;for(n=2;n<M;n+)for(j=n;j<M;j+)sigma1=0;for(s=0;s<=n-1;s+)sigma1+=Lns*Usj;Unj=Jnj-sigma1;for(i=n;i<M;i+)sigma2=0;for(s=0;s<=n-1;s+)sigma2+=Lis*Usn;Lin=(Jin-sigma2)/Unn;b1=P1;b2=Q1;b3=P2;b4=Q2;b5=P3;b6=V32;for(i=1;i<M;i+)sigma1=0;for(n=1;n<=i-1;n+)sigma1+=Lin*yn;yi=bi-sigma1;for(i=M-1;i>=1;i-)sigma2=0;for(n=i+1;n<M;n+)sigma2+=Uin*xn;xi=(yi-sigma2)/Uii;xe1=-x1;xe2=-x3;xe3=-x5;xf1=-x2;xf2=-x4;xf3=-x6;printf("节点电压:n");for(i=1;i<4;i+)ei+=xei; fi+=xfi;for(i=1;i<4;i+)printf("e%d=",i); printf("%f",ei); printf(" ,");for(i=1;i<4;i+)printf("f%d=",i); printf("%f",fi); printf(" ,");2.6 计算网络中功率分布最后要计算出平衡节点的功率和网络中的功率分布。输电线路功率的计算公式如下：3 程序运行结果经过C语言程序运行，可得以下结果。图4 程序运行结果4 小结与心得体会这次电力系统分析课程设计历时一个星期，在这一星期的日子里，不仅巩固了以前所学过的知识，而且学到了很多在书本上所没有学到过的知识。在这次的课程设计中，对于C语言的各种功能终于有了一个比较全面和具体的认识，在亲自动手编写程序的过程中，发现了很多读程序时不能发现的漏洞。在编程之余，对于电力系统分析的计算方法也有了新的了解，诸如潮流计算，节点导纳法等。并且在设计过程中需要借组同学们的智慧，以及书籍，网络解决书本外的问题，在与同学讨论的过程中收获是非常大的，比如使用LU分解法求修正方程的方法就是经过讨论和查阅资料学习并使用的。所以光靠我自己的力量是很难完成任务的，也明白了三人行必有我师的道理。通过这次课程设计使我更加体会到了理论与实际相结合的重要性，只有理论知识是远远不够的，在实践中可能会遇到各种各样的问题，不多经历就无法感受到这一点。要在实践中提高自己的动手能力和解决问题的能力，从而学以致用。最后，要感谢我的老师们的指导，他们不仅教会我专业必须掌握的知识技能，而且也使我懂得自主学习，持之以恒的道理，为将来的学习和工作夯实基础。参考文献1 何仰赞等.电力系统分析上册M武汉：华中理工大学出版社.2 何仰赞等.电力系统分析下册M武汉：华中理工大学出版社.3 诸俊伟等.电力系统分析M.北京：中国电力出版社，1995.4 周全仁等.电网计算与程序设计M.长沙：湖南科学技术出版社,1983.5丁化成.单片机应用技术A.北京:北京航空航天大学出版社,2000.附录C语言程序：#include <stdio.h>#include <math.h>#include <stdlib.h>#define N 2#define M 7main()float z12r,z12m,y120,z13r,z13m,k,z14r,z14m,y140,z24r,z24m,y240,G55,B55,J77;float e5=0,1,1,1.1,1.05,f5=0,P5,Q5,Ps5=0,-0.30,-0.55,0.5,xe4,xf4;float Qs5=0,-0.18,-0.13,V3s=1.1,V4S=1.05;float V32,max,P4,Q4;float a1=0,b1=0;int i,j,n,s,k1=0;float LMM=0,UMM=0,sigma1,sigma2,bM,yM,xM;float p12,p13,p14,p21,p24,p31,p41,p42,q12,q13,q14,q21,q24,q31,q41,q42;printf("请输入z12的实部和虚部n"); /输入电路中的阻抗scanf("%f%f",&z12r,&z12m);printf("请输入z13的实部和虚部n");scanf("%f%f",&z13r,&z13m);printf("请输入z14的实部和虚部n");scanf("%f%f",&z14r,&z14m);printf("请输入z24的实部和虚部n");scanf("%f%f",&z24r,&z24m);printf("请输入y120的值n");scanf("%f",&y120);printf("请输入y140的值n");scanf("%f",&y140);printf("请输入y240的值n");scanf("%f",&y240);printf("请输入变比k的值n");scanf("%f",&k);/*计算导纳矩阵*G11=z12r/(z12r*z12r+z12m*z12m)+k*k*z13r/(z13r*z13r+z13m*z13m)+z14r/(z14r*z14r+z14m*z14m);B11=-z12m/(z12r*z12r+z12m*z12m)-k*k*z13m/(z13r*z13r+z13m*z13m)-z14m/(z14r*z14r+z14m*z14m)+y140+y120;G22=z12r/(z12r*z12r+z12m*z12m)+z24r/(z24r*z24r+z24m*z24m);B22=-z12m/(z12r*z12r+z12m*z12m)-z24m/(z24r*z24r+z24m*z24m)+y240+y120;G33=z13r/(z13r*z13r+z13m*z13m);B33=-z13m/(z13r*z13r+z13m*z13m);G44=z14r/(z14r*z14r+z14m*z14m)+z24r/(z24r*z24r+z24m*z24m);B44=-z14m/(z14r*z14r+z14m*z14m)-z24m/(z24r*z24r+z24m*z24m)+y240+y140;G12=G21=-z12r/(z12r*z12r+z12m*z12m);B12=B21=z12m/(z12r*z12r+z12m*z12m);G13=G31=-k*z13r/(z13r*z13r+z13m*z13m);B13=B31=k*z13m/(z13r*z13r+z13m*z13m);G14=G41=-z14r/(z14r*z14r+z14m*z14m);B14=B41=z14m/(z14r*z14r+z14m*z14m);G23=G32=0.0;B23=B32=0.0;G24=G42=-z24r/(z24r*z24r+z24m*z24m);B24=B42=z24m/(z24r*z24r+z24m*z24m);G34=G43=0.0;B34=B43=0.0;for(i=1;i<5;i+)for(j=1;j<5;j+)printf("%f+%fj",Gij,Bij);printf(" ");printf("n");/形成节点导纳矩阵/*printf("n");/*/计算各节点不平衡量loop1:printf("迭代次数k1=%dn",k1);for (i=1;i<5;i+)float a=0,b=0;for(j=1;j<5;j+)a+=Gij*ej-Bij*fj; b+=Gij*fj+Bij*ej;Pi=Psi-(ei*a+fi*b);/计算有功功率的增量Qi=Qsi-(fi*a-ei*b);/计算无功功率的增量V32=V3s*V3s-e3*e3;printf("有功功率增量P1=%f",P1); printf(" ,");printf("有功功率增量P2=%f",P2); printf(" ,");printf("有功功率增量P3=%f",P3);printf("无功功率增量Q1=%f",Q1); printf(" ,");printf("无功功率增量Q2=%f",Q2); printf(" ,");printf("电压增量V32=%f",V32);printf("n");/*筛选出最大值*max=fabs(P1)>fabs(P2)?fabs(P1):fabs(P2);max=max>fabs(P3)?max:fabs(P3);max=max>fabs(Q1)?max:fabs(Q1);max=max>fabs(Q2)?max:fabs(Q2);max=max>fabs(V32)?max:fabs(V32);printf("max=%fn",max);/*while (max>0.00001)/*形成雅克比矩阵*for(j=1;j<4;j+)if(1=j)float c=0,d=0;int m;for(m=1;m<5;m+)c+=G1m*em-B1m*fm; d+=G1m*fm+B1m*em;J1*N-1j*N-1=-c-G1j*e1-B1j*f1;J1*N-1j*N=-d+B1j*e1-G1j*f1;J1*Nj*N-1=d+B1j*e1-G1j*f1;J1*Nj*N=-c+G1j*e1+B1j*f1;elseJ1*N-1j*N-1=-G1j*e1-B1j*f1; J1*Nj*N=G1j*e1-B1j*f1; J1*N-1j*N=B1j*e1-G1j*f1; J1*Nj*N-1=B1j*e1-G1j*f1;for(j=1;j<4;j+)if(2=j)float c=0,d=0;int m;for(m=1;m<5;m+)c+=G2m*em-B2m*fm; d+=G2m*fm+B2m*em;J2*N-1j*N-1=-c-G2j*e2-B2j*f2;J2*N-1j*N=-d+B2j*e2-G2j*f2;J2*Nj*N-1=d+B2j*e2-G2j*f2;J2*Nj*N=-c+G2j*e2+B2j*f2;elseJ2*N-1j*N-1=-G2j*e2-B2j*f2;J2*Nj*N=G2j*e2-B2j*f2;J2*N-1j*N=B2j*e2-G2j*f2;J2*Nj*N-1=B2j*e2-G2j*f2;for(j=1;j<4;j+)if(3=j)float c=0,d=0;int m;for(m=1;m<5;m+)c+=G3m*em-B3m*fm; d+=G3m*fm+B3m*em;J3*N-1j*N-1=-c-G3j*e3-B3j*f3;J3*N-1j*N=-d+B3j*e3-G3j*f3;J3*Nj*N-1=-2*e3;J3*Nj*N=-2*f3;elseJ3*N-1j*N-1=-G3j*e3-B3j*f3; J3*N-1j*N=B3j*e3-G3j*f3; J3*Nj*N-1=0; J3*Nj*N=0;printf("雅克比矩阵是:n");for(i=1;i<7;i+)for(j=1;j<7;j+)printf("%f",Jij);printf(" ");printf("n");/*计算修正方程*for(i=1;i<M;i+)Lii=1;for(i=1;i<M;i+)U1i=J1i;Li1=Ji1/U11;for(n=2;n<M;n+)for(j=n;j<M;j+)sigma1=0;for(s=0;s<=n-1;s+)sigma1+=Lns*Usj;Unj=Jnj-sigma1;for(i=n;i<M;i+)sigma2=0;for(s=0;s<=n-1;s+)sigma2+=Lis*Usn;Lin=(Jin-sigma2)/Unn;b1=P1;b2=Q1;b3=P2;b4=Q2;b5=P3;b6=V32;for(i=1;i<M;i+)sigma1=0;for(n=1;n<=i-1;n+)sigma1+=Lin*yn;yi=bi-sigma1;for(i=M-1;i>=1;i-)sigma2=0;for(n=i+1;n<M;n+)sigma2+=Uin*xn;xi=(yi-sigma2)/Uii;xe1=-x1;xe2=-x3;xe3=-x5;xf1=-x2;xf2=-x4;xf3=-x6;printf("节点电压:n");for(i=1;i<4;i+)ei+=xei; fi+=xfi;for(i=1;i<4;i+)printf("e%d=",i); printf("%f",ei); printf(" ,");for(i=1;i<4;i+)printf("f%d=",i); printf("%f",fi); printf(" ,");printf("n");k1=k1+1;goto loop1;for(j=1;j<5;j+)a1+=G4j*ej-B4j*fj; b1+=G4j*fj+B4j*ej;P4=e4*a1+f4*b1;Q4=f4*a1-e4*b1;printf("P4+Q4=%f+j%f",P4,Q4);printf("n");p12=-2*e1*f1*y120-(e1*(e1-e2)-f1*(f1-f2)*G12+(e1*(f1-f2)+f1*(e1-e2)*B12;q12=-(e1*e1-f1*f1)*y120+(e1*(e1-e2)+f1*(f1-f2)*B12+(e1*(f1-f2)+f1*(e1-e2)*G12;p13=(e1*e1-f1*f1)*k*(k-1)*z13r/(z13r*z13r+z13m*z13m)+2*e1*f1*k*(k-1)*z13m/(z13r*z13r+z13m*z13m)-(e1*(e1-e3)-f1*(f1-f3)*G13+(e1*(f1-f3)+f1*(e1-e3)*B13;q13=(e1*e1-f1*f1)*k*(k-1)*z13m/(z13r*z13r+z13m*z13m)-2*e1*f1*k*(k-1)*z13r/(z13r*z13r+z13m*z13m)+(e1*(e1-e3)-f1*(f1-f3)*B13-(e1*(f1-f3)+f1*(e1-e3)*G13;p14=-2*e1*f1*y140-(e1*(e1-e4)-f1*(f1-f4)*G14+(e1*(f1-f4)+f1*(e1-e4)*B14;q14=-(e1*e1-f1*f1)*y140+(e1*(e1-e4)-f1*(f1-f4)*B14+(e1*(f1-f4)+f1*(e1-e4)*G14;p21=-2*e2*f2*y120-(e2*(e2-e1)-f2*(f2-f1)*G21+(e2*(f2-f1)+f2*(e2-e1)*B21;q21=-(e2*e2-f2*f2)*y120+(e2*(e2-e1)+f2*(f2-f1)*B21+(e2*(f2-f1)+f2*(e2-e1)*G21;p24=-2*e2*f2*y240-(e2*(e2-e4)-f2*(f2-f4)*G24+(e2*(f2-f4)+f2*(e2-e4)*B24;q24=-(e2*e2-f2*f2)*y240+(e2*(e2-e4)+f2*(f2-f4)*B24+(e2*(f2-f4)+f2*(e2-e4)*G24;p31=(e3*e3-f3*f3)*(1-k)*z13r/(z13r*z13r+z13m*z13m)+2*e3*f3*(1-k)*z13m/(z13r*z13r+z13m*z13m)-(e3*(e3-e1)-f3*(f3-f1)*G31+(e3*(f3-f1)+f3*(e3-e1)*B31;q31=-(e3*e3-f3*f3)*(1-k)*z13m/(z13r*z13r+z13m*z13m)+2*e3*f3*(1-k)*z13r/(z13r*z13r+z13m*z13m)-(e3*(e3-e1)-f3*(f3-f1)*B31+(e3*(f3-f1)+f3*(e3-e1)*G31;p41=-2*e4*f4*y140-(e4*(e4-e1)-f4*(f4-f1)*G41+(e4*(f4-f1)+f4*(e4-e1)*B41;q41=-(e4*e4-f4*f4)*y140+(e4*(e4-e1)-f4*(f4-f1)*B41+(
