直角坐标系下的牛顿法潮流计算综述

直角坐标系下牛顿法潮流计算

1 计算原理

1.1 节点导纳矩阵的形成

在图1（a）的简单电力系统中，若略去变压器的励磁功率和线路电容，负荷用阻抗表示，便可以得到一个有5个节点（包括零电位点）和7条支路的等值网络，如图1（b）所示。将接于节点1和4的电势源和阻抗的串联组合变换成等值的电流源和导纳的并联组合，变得到图（c）的等值网络，其中I1y10E1和I4y40E4分别称为节点1和4的注 入电流源。

1(a)234y101y122y24y233y20(b)4y34y40Ė1y244y343y40Ė41y122y23İ1y´10y20İ4(c)图1 电力系统及其网络

以零电位点作为计算节点电压的参考点，根据基尔霍夫定律，可以写出4个独立节点的电流平衡方程如下

y12(U2U1)y20U2y23(U2U3)y24(U2U4)0y23(U3U2)y34(U3U4)0y24(U4U2)y34(U4U3)y40U4I4 （1） y10U1y12(U1U2)I1上述方程组经过整理可以写成

1

武汉理工大学《电力系统分析》课程设计说明书

Y11U1Y12U2I1Y21U1Y22U2Y23U3Y24U40Y32U2Y33U3Y34U40Y42U2Y43U3Y44U4I4 （2）

式中，Y11y10y12；Y22y20y23y24y12；Y33y23y34；Y44y40y24y34；

Y12Y21y12；Y23Y32y23；Y24Y42y24；Y34Y43y34。

一般的，对于有n个独立节点的网络，可以列写n个节点方程

Y11U1Y12U2Y21U1Y22U2Yn1U1Yn2U2也可以用矩阵写成

Y1nUnI1Y2nUnI2YnnUnIn （3）

Y11Y12Y1nU1I1YYY2nU2I22122Yn1Yn2YnnUnIn （4）

或缩写为

YUI （5）

矩阵Y称为节点导纳矩阵。它的对角线元素Yii称为节点i的自导纳，其值等于接于节点i的所有支路导纳之和。非对角线元素

Yij称为节点i、j 间的互导纳，它等于直接接于

Y0节点i、j间的支路导纳的负值。若节点i、j间不存在直接支路，则有ij。由此可知节点导纳矩阵是一个稀疏的对称矩阵。

1.2 直角坐标系下牛顿法潮流计算

采用直角坐标时，节点电压可表示为

ejf Viii导纳矩阵元素则表示为

2

武汉理工大学《电力系统分析》课程设计说明书

YijGijjBij

将上述表示式代入SiPijQiUiIiUinYUijjinj的右端，展开并分出实部和虚部，便得

n Piei(GijejBijfj)fi(GijfjBijej)j1j1

Qinn （6） fi(GijejBijfj)e(GijfjBijej) j1j1

假定系统中的第1，2，3···，m号节点为PQ节点，第i个节点的给定功率设为Pis和

Qis，对对该节点可列写方程

PiPisPiPisei(GijejBijfj)fi(GijfjBijej)0j1j1nnPiPisPiPisei(GijejBijfj)fi(GijfjBijej)0j1j1

nn

（i=1,2，···，m） （7） 假定系统中的第m+1，m+2，···，n-1号节点为PV节点，则对其中每一个节点可以列写方程

PiPisPiPisei(GijejBijfj)fi(GijfjBijej)0j1j1

Vi2Vis2Vi2Vis2(ei2fi2)0nn （i=m+1，m+2，···，n-1） （8） 第n号节点为平衡点，其电压Vnenjfn是给定的，故不参加迭代。

式（7）和式（8）总共包含了2（n-1）个方程，待求的变量有e1,f1,...,en1,fn1也是2（n-1）个。我们还可看到，方程（7）和式（8）已经具备了方程组的形式。因此，不难写出如下的修正方程式

WJV （9）

3

武汉理工大学《电力系统分析》课程设计说明书

式中

22WP1Q1...PmQmPm1Vm1...Pn1Vn1

TTVe1f1...emfmem1fm1...en1fn1

P1P1P1P1P1P1P1P1ef1emfmem1fm1en1fn11Q1Q1Q1Q1Q1Q1Q1Q1efefefef11mmm1m1n1n1PPPPPPPPmmmmmmmme1f1emfmem1fm1en1fn1QmQmQmQmQmQmQmQmef1emfmem1fm1en1fn11 JPm1Pm1Pm1Pm1Pm1Pm1Pm1Pm1ef1emfmem1fm1en1fn1122222222Vm1Vm1Vm1Vm1Vm1Vm1Vm1Vm1efefefef11mmm1m1n1n1PPPPPPPPn1n1n1n1n1n1n1n1e1f1emfmem1fm1en1fn122222222Vn1Vn1Vn1Vn1Vn1Vn1Vn1Vn1efefefef11mmm1m1n1n1上述方程中雅克比矩阵的各元素，可以对式（7）和式（8）求偏导数获得。当ij时

PiQi(GijeiBijfi)ejfjPiQi BijeiGijfi) （10）

fjejV2iV2i0ejfj当ji时

4

武汉理工大学《电力系统分析》课程设计说明书

nPi(GikekBikfk)GiieiBiifieik1nPi(GikfkBikek)BiieiGiififik1nQi(GikfkBikek)BiieiGiifieik1  （11） nQif(GikekBikfk)GiieiBiifiik1V2ie2eiiV2if2fii

修正方程式（11-48）还可以写成分块矩阵的形式

W1J11J12J1.n1 W2J21JV122J2.n1V2 W n1Jn1.1JJn1.2n1.n1Vn1式中，Wi和Vi都是二维列向量；Jij是22介方阵。

Veiif

i对于PQ节点

WPiiQ iPiPi JejfjijQ Qiiejfj对于PV节点

WPiiV2

i5

12） 13） （（武汉理工大学《电力系统分析》课程设计说明书

Piej JijV2iejPifj （14） V2ifj从表达式（1-7）~（1-11）可以看到，雅克比矩阵有以下特点：

（1）雅克比矩阵各元素都是节点电压的函数，它们的数值将在迭代过程中不断的改变。

（2）雅克比矩阵的子块Jij中的元素的表达式只用到导纳矩阵中的对应元素Yij。若

Yij0，则必有Jij0。因此，式（1-9）式中分块形式的雅克比矩阵同节点导纳矩阵一样稀疏，修正方程的求解同样可以用稀疏矩阵的求解技巧。

（3）无论在式（1-6）或式（1-9）中雅克比矩阵的元素或子块都不具有对称性。 用牛顿-拉夫逊法计算潮流的流程：首先要输入网络的原始数据以及各节点的给定值并形成节点导纳矩阵。输入节点电压初值ei(0)和fi(0)，置迭代计数k=0。然后开始进入牛顿法的迭代过程。在进行第k+1次迭代时，其计算步骤如下：

（1）按上一次迭代计算出的节点电压值e(k)和f(k)，利用式（7）和式（8）计算各类节点的不平衡量Pi(k)、Qi(k)和Vi2(k)。

（2）按条件校验收敛，即

maxPi(k)、Qi(k)、Vi2(k)<

如果收敛，迭代到此结束，转入计算各线路潮流和平衡节点的功率，并打印输出计算结果。不收敛则继续计算。

（3）利用式（10）和式（11）计算雅克比矩阵的各元素。 （4）解修正方程式（7）求节点电压的修正量ei(k)和fi(k)。 （5）修正各节点的电压

ei(k1)ei(k)ei(k),fi(k1)fi(k)fi(k)

（6）迭代计数加1，返回第一步继续迭代过程。

6

武汉理工大学《电力系统分析》课程设计说明书

2 计算过程

2．1 计算程序框图

图2 潮流计算程序框图

7

武汉理工大学《电力系统分析》课程设计说明书

2.2 节点导纳矩阵的形成

图3 节点导纳矩阵

根据第一章节点导纳矩阵的形成原理，设计题目可化为如图含有4节点的形式。 采用直角坐标时，节点电压可表示为

ejf Viii导纳矩阵元素则表示为

YijGijjBij

//****************计算导纳矩阵*******************

G[1][1]=z12r/(z12r*z12r+z12m*z12m)+k*k*z13r/(z13r*z13r+z13m*z13m)+z14r/(z14r*z14r+z14m*z14m); B[1][1]=-z12m/(z12r*z12r+z12m*z12m)-k*k*z13m/(z13r*z13r+z13m*z13m)-z14m/(z14r*z14r+z14m*z14m)+y140+y120;

G[2][2]=z12r/(z12r*z12r+z12m*z12m)+z24r/(z24r*z24r+z24m*z24m);

B[2][2]=-z12m/(z12r*z12r+z12m*z12m)-z24m/(z24r*z24r+z24m*z24m)+y240+y120; G[3][3]=z13r/(z13r*z13r+z13m*z13m); B[3][3]=-z13m/(z13r*z13r+z13m*z13m);

G[4][4]=z14r/(z14r*z14r+z14m*z14m)+z24r/(z24r*z24r+z24m*z24m);

8

武汉理工大学《电力系统分析》课程设计说明书

B[4][4]=-z14m/(z14r*z14r+z14m*z14m)-z24m/(z24r*z24r+z24m*z24m)+y240+y140; G[1][2]=G[2][1]=-z12r/(z12r*z12r+z12m*z12m); B[1][2]=B[2][1]=z12m/(z12r*z12r+z12m*z12m); G[1][3]=G[3][1]=-k*z13r/(z13r*z13r+z13m*z13m); B[1][3]=B[3][1]=k*z13m/(z13r*z13r+z13m*z13m); G[1][4]=G[4][1]=-z14r/(z14r*z14r+z14m*z14m); B[1][4]=B[4][1]=z14m/(z14r*z14r+z14m*z14m); G[2][3]=G[3][2]=0.0; B[2][3]=B[3][2]=0.0;

G[2][4]=G[4][2]=-z24r/(z24r*z24r+z24m*z24m); B[2][4]=B[4][2]=z24m/(z24r*z24r+z24m*z24m); G[3][4]=G[4][3]=0.0; B[3][4]=B[4][3]=0.0; for(i=1;i<5;i++) {for(j=1;j<5;j++)

{printf(\"%f+%fj\printf(\" \"); }

printf(\"\\n\");//形成节点导纳矩阵

将上述表示式代入SiPijQiUiIiUi便得

YUijjinj的右端，展开并分出实部和虚部，

Qifi(GijejBijfj)ei(GijfjBijej)j1j1 j1 Piei(GijejBijfj)fi(GijfjBijej)nj1nnn

2.3 计算各节点不平衡量

假定系统中的第1，2，3···，m号节点为PQ节点，第i个节点的给定功率设为Pis和

9

武汉理工大学《电力系统分析》课程设计说明书

Qis，对对该节点可列写方程

PiPisPiPisei(GijejBijfj)fi(GijfjBijej)0j1j1nnnnPiPisPiPisei(GijejBijfj)fi(GijfjBijej)0j1j1

（i=1,2，···，m）

假定系统中的第m+1，m+2，···，n-1号节点为PV节点，则对其中每一个节点可以列写方程

PiPisPiPisei(GijejBijfj)fi(GijfjBijej)0j1j1 Vi2Vis2Vi2Vis2(ei2fi2)0nn （i=m+1，m+2，···，n-1）

第n号节点为平衡点，其电压Vnenjfn是给定的，故不参加迭代，其计算程序如下：

//计算各节点不平衡量 loop1:

printf(\"迭代次数k1=%d\\n\for (i=1;i<5;i++) {float a=0,b=0; for(j=1;j<5;j++)

{a+=G[i][j]*e[j]-B[i][j]*f[j]; b+=G[i][j]*f[j]+B[i][j]*e[j]; }

P[i]=Ps[i]-(e[i]*a+f[i]*b);//计算有功功率的增量 Q[i]=Qs[i]-(f[i]*a-e[i]*b);//计算无功功率的增量

V32=V3s*V3s-e[3]*e[3]; }

printf(\"有功功率增量P[1]=%f\ printf(\" ,\");

printf(\"有功功率增量P[2]=%f\ printf(\" ,\");

printf(\"有功功率增量P[3]=%f\

printf(\"无功功率增量Q[1]=%f\ printf(\" ,\");

printf(\"无功功率增量Q[2]=%f\ printf(\" ,\");

printf(\"电压增量V32=%f\

printf(\"\\n\");

10

武汉理工大学《电力系统分析》课程设计说明书

2.4 雅克比矩阵计算

上述方程中雅克比矩阵的各元素，可以对计算各点不平衡量得公式中求偏导数获得。当ij时

PiQi(GijeiBijfi)ejfjPiQiBijeiGijfi) fjejV2iV2i0ejfj当ji时

nPi(GikekBikfk)GiieiBiifieik1nPi(GikfkBikek)BiieiGiififik1nQi(GikfkBikek)BiieiGiifieik1 nQi(GikekBikfk)GiieiBiififik1V2i2eiei2Vi2fifi

以下为程序，只列出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+=G[1][m]*e[m]-B[1][m]*f[m]; d+=G[1][m]*f[m]+B[1][m]*e[m];

11

武汉理工大学《电力系统分析》课程设计说明书

}

J[1*N-1][j*N-1]=-c-G[1][j]*e[1]-B[1][j]*f[1]; J[1*N-1][j*N]=-d+B[1][j]*e[1]-G[1][j]*f[1]; J[1*N][j*N-1]=d+B[1][j]*e[1]-G[1][j]*f[1]; J[1*N][j*N]=-c+G[1][j]*e[1]+B[1][j]*f[1]; } else

{J[1*N-1][j*N-1]=-G[1][j]*e[1]-B[1][j]*f[1]; J[1*N][j*N]=G[1][j]*e[1]-B[1][j]*f[1]; J[1*N-1][j*N]=B[1][j]*e[1]-G[1][j]*f[1]; J[1*N][j*N-1]=B[1][j]*e[1]-G[1][j]*f[1]; } }

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;ifor(i=1;iU[1][i]=J[1][i];

12

武汉理工大学《电力系统分析》课程设计说明书

L[i][1]=J[i][1]/U[1][1]; }

for(n=2;nfor(j=n;jfor(i=n;ib[1]=P[1];b[2]=Q[1];b[3]=P[2]; b[4]=Q[2];b[5]=P[3];b[6]=V32;

for(i=1;ifor(n=1;n<=i-1;n++)

sigma1+=L[i][n]*y[n]; y[i]=b[i]-sigma1; }

for(i=M-1;i>=1;i--) { sigma2=0;

for(n=i+1;nxe[1]=-x[1];xe[2]=-x[3];xe[3]=-x[5]; xf[1]=-x[2];xf[2]=-x[4];xf[3]=-x[6]; printf(\"节点电压:\\n\"); for(i=1;i<4;i++) {e[i]+=xe[i]; f[i]+=xf[i]; }

for(i=1;i<4;i++) {printf(\"e[%d]=\ printf(\"%f\ printf(\" ,\"); }

for(i=1;i<4;i++) {printf(\"f[%d]=\ printf(\"%f\ printf(\" ,\");

}

13

武汉理工大学《电力系统分析》课程设计说明书

2.6 计算网络中功率分布

最后要计算出平衡节点的功率和网络中的功率分布。输电线路功率的计算公式如下：

SijPijjQijViIijViyi0Vi(ViVj)yij

.*2*.***3 程序运行结果

经过C语言程序运行，可得以下结果。

图4 程序运行结果

14

武汉理工大学《电力系统分析》课程设计说明书

4 小结与心得体会

这次电力系统分析课程设计历时一个星期，在这一星期的日子里，不仅巩固了以前所学过的知识，而且学到了很多在书本上所没有学到过的知识。

在这次的课程设计中，对于C语言的各种功能终于有了一个比较全面和具体的认识，在亲自动手编写程序的过程中，发现了很多读程序时不能发现的漏洞。在编程之余，对于电力系统分析的计算方法也有了新的了解，诸如潮流计算，节点导纳法等。并且在设计过程中需要借组同学们的智慧，以及书籍，网络解决书本外的问题，在与同学讨论的过程中收获是非常大的，比如使用LU分解法求修正方程的方法就是经过讨论和查阅资料学习并使用的。所以光靠我自己的力量是很难完成任务的，也明白了三人行必有我师的道理。

通过这次课程设计使我更加体会到了理论与实际相结合的重要性，只有理论知识是远远不够的，在实践中可能会遇到各种各样的问题，不多经历就无法感受到这一点。要在实践中提高自己的动手能力和解决问题的能力，从而学以致用。最后，要感谢我的老师们的指导，他们不仅教会我专业必须掌握的知识技能，而且也使我懂得自主学习，持之以恒的道理，为将来的学习和工作夯实基础。

15

武汉理工大学《电力系统分析》课程设计说明书

参考文献

[1] 何仰赞等.电力系统分析上册[M]． 武汉：华中理工大学出版社. [2] 何仰赞等.电力系统分析下册[M]． 武汉：华中理工大学出版社. [3] 诸俊伟等.电力系统分析[M].

北京：中国电力出版社，1995.

[4] 周全仁等.电网计算与程序设计[M].长沙：湖南科学技术出版社,1983. [5]丁化成.单片机应用技术[A].

16

北京:北京航空航天大学出版社,2000.

武汉理工大学《电力系统分析》课程设计说明书

附录

C语言程序：

#include #include #include #define N 2 #define M 7 main()

{float z12r,z12m,y120,z13r,z13m,k,z14r,z14m,y140,z24r,z24m,y240,G[5][5],B[5][5],J[7][7]; float e[5]={0,1,1,1.1,1.05},f[5]={0},P[5],Q[5],Ps[5]={0,-0.30,-0.55,0.5},xe[4],xf[4]; float Qs[5]={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 L[M][M]={0},U[M][M]={0},sigma1,sigma2,b[M],y[M],x[M]; float p12,p13,p14,p21,p24,p31,p41,p42,q12,q13,q14,q21,q24,q31,q41,q42; printf(\"请输入z12的实部和虚部\\n\"); //输入电路中的阻抗 scanf(\"%f%f\printf(\"请输入z13的实部和虚部\\n\"); scanf(\"%f%f\printf(\"请输入z14的实部和虚部\\n\"); scanf(\"%f%f\printf(\"请输入z24的实部和虚部\\n\"); scanf(\"%f%f\printf(\"请输入y120的值\\n\"); scanf(\"%f\

printf(\"请输入y140的值\\n\"); scanf(\"%f\

printf(\"请输入y240的值\\n\");

17

武汉理工大学《电力系统分析》课程设计说明书

scanf(\"%f\

printf(\"请输入变比k的值\\n\"); scanf(\"%f\

//****************计算导纳矩阵*******************

G[1][1]=z12r/(z12r*z12r+z12m*z12m)+k*k*z13r/(z13r*z13r+z13m*z13m)+z14r/(z14r*z14r+z14m*z14m); B[1][1]=-z12m/(z12r*z12r+z12m*z12m)-k*k*z13m/(z13r*z13r+z13m*z13m)-z14m/(z14r*z14r+z14m*z14m)+y140+y120;

G[2][2]=z12r/(z12r*z12r+z12m*z12m)+z24r/(z24r*z24r+z24m*z24m);

B[2][2]=-z12m/(z12r*z12r+z12m*z12m)-z24m/(z24r*z24r+z24m*z24m)+y240+y120; G[3][3]=z13r/(z13r*z13r+z13m*z13m); B[3][3]=-z13m/(z13r*z13r+z13m*z13m);

G[4][4]=z14r/(z14r*z14r+z14m*z14m)+z24r/(z24r*z24r+z24m*z24m);

B[4][4]=-z14m/(z14r*z14r+z14m*z14m)-z24m/(z24r*z24r+z24m*z24m)+y240+y140; G[1][2]=G[2][1]=-z12r/(z12r*z12r+z12m*z12m); B[1][2]=B[2][1]=z12m/(z12r*z12r+z12m*z12m); G[1][3]=G[3][1]=-k*z13r/(z13r*z13r+z13m*z13m); B[1][3]=B[3][1]=k*z13m/(z13r*z13r+z13m*z13m); G[1][4]=G[4][1]=-z14r/(z14r*z14r+z14m*z14m); B[1][4]=B[4][1]=z14m/(z14r*z14r+z14m*z14m); G[2][3]=G[3][2]=0.0; B[2][3]=B[3][2]=0.0;

G[2][4]=G[4][2]=-z24r/(z24r*z24r+z24m*z24m); B[2][4]=B[4][2]=z24m/(z24r*z24r+z24m*z24m); G[3][4]=G[4][3]=0.0; B[3][4]=B[4][3]=0.0; for(i=1;i<5;i++) {for(j=1;j<5;j++)

{printf(\"%f+%fj\printf(\" \");

18

武汉理工大学《电力系统分析》课程设计说明书

}

printf(\"\\n\");//形成节点导纳矩阵

//******************************************* }

printf(\"\\n\");

//******************************************** //计算各节点不平衡量 loop1:

printf(\"迭代次数k1=%d\\n\for (i=1;i<5;i++) {float a=0,b=0; for(j=1;j<5;j++)

{a+=G[i][j]*e[j]-B[i][j]*f[j]; b+=G[i][j]*f[j]+B[i][j]*e[j]; }

P[i]=Ps[i]-(e[i]*a+f[i]*b);//计算有功功率的增量 Q[i]=Qs[i]-(f[i]*a-e[i]*b);//计算无功功率的增量

V32=V3s*V3s-e[3]*e[3]; }

printf(\"有功功率增量P[1]=%f\ printf(\" ,\");

printf(\"有功功率增量P[2]=%f\ printf(\" ,\");

printf(\"有功功率增量P[3]=%f\

printf(\"无功功率增量Q[1]=%f\ printf(\" ,\");

printf(\"无功功率增量Q[2]=%f\

19

武汉理工大学《电力系统分析》课程设计说明书

printf(\" ,\");

printf(\"电压增量V32=%f\

printf(\"\\n\");

//************筛选出最大值***********************

max=fabs(P[1])>fabs(P[2])?fabs(P[1]):fabs(P[2]); max=max>fabs(P[3])?max:fabs(P[3]); max=max>fabs(Q[1])?max:fabs(Q[1]); max=max>fabs(Q[2])?max:fabs(Q[2]); max=max>fabs(V32)?max:fabs(V32); printf(\"max=%f\\n\

//******************************************** 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+=G[1][m]*e[m]-B[1][m]*f[m]; d+=G[1][m]*f[m]+B[1][m]*e[m]; }

J[1*N-1][j*N-1]=-c-G[1][j]*e[1]-B[1][j]*f[1]; J[1*N-1][j*N]=-d+B[1][j]*e[1]-G[1][j]*f[1]; J[1*N][j*N-1]=d+B[1][j]*e[1]-G[1][j]*f[1]; J[1*N][j*N]=-c+G[1][j]*e[1]+B[1][j]*f[1];

20

武汉理工大学《电力系统分析》课程设计说明书

} else

{J[1*N-1][j*N-1]=-G[1][j]*e[1]-B[1][j]*f[1]; J[1*N][j*N]=G[1][j]*e[1]-B[1][j]*f[1]; J[1*N-1][j*N]=B[1][j]*e[1]-G[1][j]*f[1]; J[1*N][j*N-1]=B[1][j]*e[1]-G[1][j]*f[1]; } }

for(j=1;j<4;j++) {if(2==j) {float c=0,d=0; int m;

for(m=1;m<5;m++)

{c+=G[2][m]*e[m]-B[2][m]*f[m]; d+=G[2][m]*f[m]+B[2][m]*e[m]; }

J[2*N-1][j*N-1]=-c-G[2][j]*e[2]-B[2][j]*f[2]; J[2*N-1][j*N]=-d+B[2][j]*e[2]-G[2][j]*f[2]; J[2*N][j*N-1]=d+B[2][j]*e[2]-G[2][j]*f[2]; J[2*N][j*N]=-c+G[2][j]*e[2]+B[2][j]*f[2]; } else

{J[2*N-1][j*N-1]=-G[2][j]*e[2]-B[2][j]*f[2]; J[2*N][j*N]=G[2][j]*e[2]-B[2][j]*f[2]; J[2*N-1][j*N]=B[2][j]*e[2]-G[2][j]*f[2]; J[2*N][j*N-1]=B[2][j]*e[2]-G[2][j]*f[2]; } }

21

武汉理工大学《电力系统分析》课程设计说明书

for(j=1;j<4;j++) {if(3==j) {float c=0,d=0; int m;

for(m=1;m<5;m++)

{c+=G[3][m]*e[m]-B[3][m]*f[m]; d+=G[3][m]*f[m]+B[3][m]*e[m]; }

J[3*N-1][j*N-1]=-c-G[3][j]*e[3]-B[3][j]*f[3]; J[3*N-1][j*N]=-d+B[3][j]*e[3]-G[3][j]*f[3]; J[3*N][j*N-1]=-2*e[3]; J[3*N][j*N]=-2*f[3]; } else

{J[3*N-1][j*N-1]=-G[3][j]*e[3]-B[3][j]*f[3]; J[3*N-1][j*N]=B[3][j]*e[3]-G[3][j]*f[3]; J[3*N][j*N-1]=0; J[3*N][j*N]=0; } }

printf(\"雅克比矩阵是:\\n\"); for(i=1;i<7;i++) {for(j=1;j<7;j++) {printf(\"%f\printf(\" \"); }

printf(\"\\n\"); }

//********计算修正方程*************

22

武汉理工大学《电力系统分析》课程设计说明书

for(i=1;ifor(i=1;iU[1][i]=J[1][i]; L[i][1]=J[i][1]/U[1][1]; }

for(n=2;nfor(j=n;jfor(i=n;ib[1]=P[1];b[2]=Q[1];b[3]=P[2]; b[4]=Q[2];b[5]=P[3];b[6]=V32;

23

武汉理工大学《电力系统分析》课程设计说明书

for(i=1;ifor(n=1;n<=i-1;n++) sigma1+=L[i][n]*y[n]; y[i]=b[i]-sigma1; }

for(i=M-1;i>=1;i--) { sigma2=0;

for(n=i+1;nxe[1]=-x[1];xe[2]=-x[3];xe[3]=-x[5]; xf[1]=-x[2];xf[2]=-x[4];xf[3]=-x[6]; printf(\"节点电压:\\n\"); for(i=1;i<4;i++) {e[i]+=xe[i]; f[i]+=xf[i]; }

for(i=1;i<4;i++) {printf(\"e[%d]=\ printf(\"%f\ printf(\" ,\"); }

for(i=1;i<4;i++) {printf(\"f[%d]=\

24

武汉理工大学《电力系统分析》课程设计说明书

printf(\"%f\ printf(\" ,\"); }

printf(\"\\n\"); k1=k1+1; goto loop1; }

for(j=1;j<5;j++)

{a1+=G[4][j]*e[j]-B[4][j]*f[j]; b1+=G[4][j]*f[j]+B[4][j]*e[j]; }

P4=e[4]*a1+f[4]*b1; Q4=f[4]*a1-e[4]*b1;

printf(\"P4+Q4=%f+j%f\printf(\"\\n\");

p12=-2*e[1]*f[1]*y120-(e[1]*(e[1]-e[2])-f[1]*(f[1]-f[2]))*G[1][2]+(e[1]*(f[1]-f[2])+f[1]*(e[1]-e[2]))*B[1][2];

q12=-(e[1]*e[1]-f[1]*f[1])*y120+(e[1]*(e[1]-e[2])+f[1]*(f[1]-f[2]))*B[1][2]+(e[1]*(f[1]-f[2])+f[1]*(e[1]-e[2]))*G[1][2];

p13=(e[1]*e[1]-f[1]*f[1])*k*(k-1)*z13r/(z13r*z13r+z13m*z13m)+2*e[1]*f[1]*k*(k-1)*z13m/(z13r*z13r+z13m*z13m)-(e[1]*(e[1]-e[3])-f[1]*(f[1]-f[3]))*G[1][3]+(e[1]*(f[1]-f[3])+f[1]*(e[1]-e[3]))*B[1][3];

q13=(e[1]*e[1]-f[1]*f[1])*k*(k-1)*z13m/(z13r*z13r+z13m*z13m)-2*e[1]*f[1]*k*(k-1)*z13r/(z13r*z13r+z13m*z13m)+(e[1]*(e[1]-e[3])-f[1]*(f[1]-f[3]))*B[1][3]-(e[1]*(f[1]-f[3])+f[1]*(e[1]-e[3]))*G[1][3];

p14=-2*e[1]*f[1]*y140-(e[1]*(e[1]-e[4])-f[1]*(f[1]-f[4]))*G[1][4]+(e[1]*(f[1]-f[4])+f[1]*(e[1]-e[4]))*B[1][4];

q14=-(e[1]*e[1]-f[1]*f[1])*y140+(e[1]*(e[1]-e[4])-f[1]*(f[1]-f[4]))*B[1][4]+(e[1]*(f[1]-f[4])+f[1]*(e[1]-e[4]))*G[1][4];

25

武汉理工大学《电力系统分析》课程设计说明书

p21=-2*e[2]*f[2]*y120-(e[2]*(e[2]-e[1])-f[2]*(f[2]-f[1]))*G[2][1]+(e[2]*(f[2]-f[1])+f[2]*(e[2]-e[1]))*B[2][1];

q21=-(e[2]*e[2]-f[2]*f[2])*y120+(e[2]*(e[2]-e[1])+f[2]*(f[2]-f[1]))*B[2][1]+(e[2]*(f[2]-f[1])+f[2]*(e[2]-e[1]))*G[2][1];

p24=-2*e[2]*f[2]*y240-(e[2]*(e[2]-e[4])-f[2]*(f[2]-f[4]))*G[2][4]+(e[2]*(f[2]-f[4])+f[2]*(e[2]-e[4]))*B[2][4];

q24=-(e[2]*e[2]-f[2]*f[2])*y240+(e[2]*(e[2]-e[4])+f[2]*(f[2]-f[4]))*B[2][4]+(e[2]*(f[2]-f[4])+f[2]*(e[2]-e[4]))*G[2][4];

p31=(e[3]*e[3]-f[3]*f[3])*(1-k)*z13r/(z13r*z13r+z13m*z13m)+2*e[3]*f[3]*(1-k)*z13m/(z13r*z13r+z13m*z13m)-(e[3]*(e[3]-e[1])-f[3]*(f[3]-f[1]))*G[3][1]+(e[3]*(f[3]-f[1])+f[3]*(e[3]-e[1]))*B[3][1];

q31=-(e[3]*e[3]-f[3]*f[3])*(1-k)*z13m/(z13r*z13r+z13m*z13m)+2*e[3]*f[3]*(1-k)*z13r/(z13r*z13r+z13m*z13m)-(e[3]*(e[3]-e[1])-f[3]*(f[3]-f[1]))*B[3][1]+(e[3]*(f[3]-f[1])+f[3]*(e[3]-e[1]))*G[3][1];

p41=-2*e[4]*f[4]*y140-(e[4]*(e[4]-e[1])-f[4]*(f[4]-f[1]))*G[4][1]+(e[4]*(f[4]-f[1])+f[4]*(e[4]-e[1]))*B[4][1];

q41=-(e[4]*e[4]-f[4]*f[4])*y140+(e[4]*(e[4]-e[1])-f[4]*(f[4]-f[1]))*B[4][1]+(e[4]*(f[4]-f[1])+f[4]*(e[4]-e[1]))*G[4][1];

p42=-2*e[4]*f[4]*y240-(e[4]*(e[4]-e[2])-f[4]*(f[4]-f[2]))*G[4][2]+(e[4]*(f[4]-f[2])+f[4]*(e[4]-e[2]))*B[4][2];

q42=-(e[4]*e[4]-f[4]*f[4])*y240+(e[4]*(e[4]-e[2])+f[4]*(f[4]-f[2]))*B[4][2]+(e[4]*(f[4]-f[2])+f[4]*(e[4]-e[2]))*G[4][2];

printf(\"s12=%f+j%f\printf(\"s13=%f+j%f\printf(\"s14=%f+j%f\printf(\"s21=%f+j%f\printf(\"s24=%f+j%f\printf(\"s31=%f+j%f\printf(\"s41=%f+j%f\printf(\"s42=%f+j%f\

}

26