电力系统计算编程

牛顿-拉夫逊法进行潮流计算

如图所示的一个电力网络，求系统中的潮流分布。计算精度要求各节点电压修正量不大于-510。

G

1 3 4 1.250-j3.75 10.000-j30.000

1.667-j5.000 1.667-j5.000

1.250-j3.7505(0) 2.500-j7.500 2 5

G

已知该系统中，节点1为平衡节点，保持U11.06j0 为定值，其余四个节点都是PQ节点，给定的注入功率分别为：

~0.20j0.20,S~0.45j0.15,S~0.40j0.05,S~0.60j0.10S2345计算步骤：

1. 形成节点导纳矩阵YB

2. 计算各节点功率的不平衡量

运用如下公式： Pi(0)5.00-j15.00 jn[eij1jn(0)(Gijej(0)Bijfj(0))fi(Gijfj(0)Bijej(0))]

Qi(0)j1[fi(0)(Gijej(0)Bijfi(0))ei(0)(Gijfj(0)Bijej(0))]

计算各节点功率Pi于是 △Pi(0)(0)、Qi(0)(0)

(0)=Pi─Pi；△Qi=Qi─Qi(0)

3. 计算雅可比矩阵中的各元素

先计算各节点注入电流：

(0) Ii=

=aii（0）+jbii（0）

然后计算雅可比矩阵各元素：

（1）当i≠j时雅可比矩阵的各个元素分别为： HijBijeiGijfi ；

NijGijeiBijfi ； Jij-Gijei-Bijfi ； LijBijeiGijfi ；

（2）当i=j时雅可比矩阵的各个元素分别为： HiiBiieiGiifibii;

NiiGiieiBiifaii; JiiGiieiBiifaii; LiiBiieiGiifibii; 列出雅可比矩阵J(0)。

4. 解修正方程式求各节点电压

求雅可比矩阵J(0)的逆矩阵（J(0)），再由

-1

（

J(0)）-1

求得

，则可以根据ei(1)ei(0)+△ei(0)，fi(1)fi(0)+△fi(0)得到各节点电压

的新值，即修正后的值，然后就可以开始第二次迭代，以此循环往复，直到满足各节点电压修正量不大于10-5的要求。 5. 计算平衡节点功率

平衡节点功率

：

和线路功率

1j=线路功率

：

YUj）

ij=）YUi

ji=）YUj

6. 将各节点的电压以极坐标的形式表示

用C语言编程的程序如下：

#include

#include

double divRe(double b1,double b2,double b3,double b4) {

double a1r;

a1r=(b1*b3+b2*b4)/(b3*b3+b4*b4);

return(a1r); }

double divIm(double b1,double b2,double b3,double b4) {

double a1i;

a1i=( b1*b4-b2*b3)/(b3*b3+b4*b4);

return(a1i); }

double mulRe(double b1,double b2,double b3,double b4) {

double a2r; a2r=b1*b3-b2*b4; return(a2r);

}

double mulIm(double b1,double b2,double b3,double b4) {

double a2i;

a2i=b2*b3+b1*b4; return(a2i); }

double detu(double b1,double b2) {

double udet;

udet=b1+b2;

return(udet); }

double Max(double a[],int n) { int i;

double max; for(i=0;iif(a[i]>a[i+1])

{max=a[i];a[i]=a[i+1];a[i+1]=max;} return(max);

}//找寻△ei[i]和△fi[i]中的最大值 void main() {

int i,j,k,h,km;

double t;

double eps,sumpi1,sumpi2,sumqi1,sumqi2,max,sumir,sumii,sumid,sumiq,I1r,I1i,l1d,l1q,l1m,l1n; double pi0[5],qi0[5],detpi[5],detqi[5],aii0[5],bii0[5],J0[8][8],detsi[8],detui[8],J0ni[8][8] ,H[4][4],N[4][4],J[4][4],L[4][4], euij[5][5], fuij[5][5];

static double ybr[5][5]={{6.250,-5.000,-1.250,0,0},{-5.000,10.834,-1.667,-1.667,-2.500},

{-1.250,-1.667,12.917,-10.000,0},{0,-1.667,-10.000,12.917,-1.250}, {0,-2.500,0,-1.250,3.750}};

static double ybi[5][5]={{-18.750,15.000,3.750,0,0},{15.000,-32.500,5.000,5.000,7.500}, {3.750,5.000,-38.750,30.000,0},{0,5.000,30.000,-38.750,3.750}, {0,7.500,0,3.750,-11.250}}; double ei0[5]={1.06,1.0,1.0,1.0,1.0}; double fi0[5]={0,0,0,0,0};

double pi[5]={0,0.2,-0.45,-0.4,-0.6}; double qi[5]={0,0.2,-0.15,-0.05,-0.1}; k=0; km=6; eps=0.00001; for(i=1;i<5;i++)

{printf(\"ei0[%d]=%f \ printf(\"fi1[%d]=%f\\n\ } do { k+=1;

printf(\"Now start...\\n\");

printf(\"The %d times\\n\ sumpi2=0; sumqi2=0;

for(i=1;i<5;i++)//第一个节点为平衡节点，从第二个节点开始求pi0[i] {for(j=0;j<5;j++) {

sumpi1=(ei0[i]*(ybr[i][j]*ei0[j]-ybi[i][j]*fi0[j])+fi0[i]*(ybr[i][j]*fi0[j]+ybi[i][j]*ei0[j]));

sumpi2+=sumpi1; }

pi0[i]=sumpi2;

printf(\"pi0[%d]=%f \ sumpi2=0;

}//完成pi0[i]的求解

for(i=1;i<5;i++)//第一个节点为平衡节点，从第二个节点开始求qi0[i] {for(j=0;j<5;j++)

{

sumqi1=(fi0[i]*(ybr[i][j]*ei0[j]-ybi[i][j]*fi0[j])-ei0[i]*(ybr[i][j]*fi0[j]+ybi[i][j]*ei0[j])); sumqi2+=sumqi1;

}

qi0[i]=sumqi2;

printf(\"qi0[%d]=%f \ sumqi2=0;

}//完成qi0[i]的求解 printf(\"\\n\"); for(i=1;i<5;i++) {

detpi[i]=pi[i]-pi0[i];

detqi[i]=qi[i]-qi0[i];

printf(\"detpi[%d]=%f \ printf(\"detqi[%d]=%f\\n\ }//计算△pi0[i]，△qi0[i]

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

{aii0[i]=divRe(pi0[i],qi0[i],ei0[i],fi0[i]); bii0[i]=divIm(pi0[i],qi0[i],ei0[i],fi0[i]); printf(\"aii0[%d]=%f \ printf(\"bii0[%d]=%f\\n\ }//求解aii0[i]，bii0[i] for(i=0;i<4;i++) { for(j=0;j<4;j++) if(i==j) {

H[i][j]=-ybi[i+1][j+1]*ei0[i+1]+ybr[i+1][j+1]*fi0[i+1]+bii0[i+1]; N[i][j]=ybr[i+1][j+1]*ei0[i+1]+ybi[i+1][j+1]*fi0[i+1]+aii0[i+1]; J[i][j]=-ybr[i+1][j+1]*ei0[i+1]-ybi[i+1][i+1]*fi0[i+1]+aii0[i+1]; L[i][j]=-ybi[i+1][j+1]*ei0[i+1]+ybr[i+1][j+1]*fi0[i+1]-bii0[i+1]; }

else {

H[i][j]=ybr[i+1][j+1]*fi0[i+1]-ybi[i+1][j+1]*ei0[i+1]; N[i][j]=ybr[i+1][j+1]*ei0[i+1]+ybi[i+1][j+1]*fi0[i+1]; J[i][j]=-ybi[i+1][j+1]*fi0[i+1]-ybr[i+1][j+1]*ei0[i+1]; L[i][j]=ybr[i+1][j+1]*fi0[i+1]-ybi[i+1][j+1]*ei0[i+1];

}

} //求取雅可比矩阵中的各元素 for(i=0;i<8;i++) for(j=0;j<8;j++) {

if(i%2==0&&j%2==0) J0[i][j]=H[i/2][j/2];

else if(i%2==0&&j%2!=0) J0[i][j]=N[i/2][(j-1)/2]; else if(i%2!=0&&j%2==0) J0[i][j]=J[(i-1)/2][j/2]; else J0[i][j]=L[i/2][(j-1)/2];

}//形成雅可比矩阵J0[i][j] for(i=0;i<8;i++) {

for(j=0;j<8;j++)

printf(\"J0[%d][%d]=%3.3f \ printf(\"\\n\");

}//显示雅可比矩阵J0[i][j] for(i=0;i<8;i++) {if(i%2==0) detsi[i]=detpi[(i+2)/2]; else detsi[i]=detqi[(i+1)/2];

printf(\"detsi[%d]=%f\\n\}//形成△pi[i]，△qi[i]的组合矩阵

for(i=0;i<8;i++)//雅可比矩阵J0[i][j]求逆矩阵 for(j=0;j<8;j++) {

if(i!=j) J0ni[i][j]=0; else

J0ni[i][j]=1;

}//逆矩阵J0ni[i][j]初始化，单位矩阵

for(i=0;i<8;i++) //将雅可比矩阵J0[i][j]化简为对角阵,逆矩阵跟随变化 {

for(j=0;j<8;j++) {

if(i!=j)

{

t=J0[j][i]/J0[i][i]; for(h=0;h<8;h++) {

J0[j][h]-=J0[i][h]*t; J0ni[j][h]-=J0ni[i][h]*t; } } } }

for(i=0;i<8;i++)//雅可比矩阵J0[i][j]对角元素化为1 if(J0[i][i]!=1) {

t=J0[i][i];

for(j=0;j<8;j++) J0ni[i][j]=J0ni[i][j]/t; }//求得逆矩阵J0ni[i][j] for(i=0;i<8;i++) { for(j=0;j<8;j++)

{printf(\"J0ni[%d][%d]=%1.5f \ printf(\"\\n\"); }//显示逆矩阵J0ni[i][j] for(i=0;i<8;i++) {

detui[i]=0; for(j=0;j<8;j++)

detui[i]-=J0ni[i][j]*detsi[j]; detui[i]*=(-1); }

for(i=0;i<8;i++) { printf(\"detui[%d]=%f \求得△ui[i] printf(\"\\n\"); for(i=1;i<5;i++) {

fi0[i]+=detui[2*i-2];

ei0[i]+=detui[2*i-1];

}//求得每一次迭代后各节点的电压 for(i=1;i<5;i++)

{printf(\"ei0[%d]=%f \ printf(\"fi1[%d]=%f\\n\ }

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

{pi[i]=detpi[i]+pi0[i]; qi[i]=detqi[i]+qi0[i]; }

max=Max(detui,8); printf(\"max=%f\\n\

}while(max>eps&&kfor(i=0;i<5;i++) {

I1r=mulRe(ybr[0][i],-ybi[0][i],ei0[i],-fi0[i]);

I1i=mulIm(ybr[0][i],-ybi[0][i],ei0[i],-fi0[i]);//复数运算第一步，两个复数相乘 sumir+=I1r; sumii+=I1i; }

pi[0]=mulRe(ei0[0],fi0[0],sumir,sumii);

qi[0]=mulIm(ei0[0],fi0[0],sumir,sumii);//复数运算第二步 printf(\"S1=%f+j%f\\n\求平衡点的功率 sumid=0; sumiq=0; for(i=1;i<6;i++) for(j=1;j<6;j++) {

euij[i-1][j-1]=detu(ei0[i-1],-ei0[j-1]); fuij[i-1][j-1]=detu(-fi0[i-1],fi0[j-1]);

l1d= mulRe(euij[i-1][j-1], fuij[i-1][j-1],- ybr[i-1][j-1], ybi[i-1][j-1]); l1q=mulIm(euij[i-1][j-1], fuij[i-1][j-1],- ybr[i-1][j-1], ybi[i-1][j-1]); l1m=mulRe(ei0[i-1], l1q,l1d, fi0[i-1]); l1n=mulIm(ei0[i-1], fi0[i-1],l1d, l1q);

printf(\"S[%d][%d]=%f+j%f\\ }//求线路功率 for(i=0;i<5;i++)

printf(\"u%d=%f∠%f\\n\3.14159); }//各节点电压的极坐标表示