牛顿-拉夫逊迭代法极坐标潮流计算C语言程序

/*利用牛顿-拉夫逊迭代法(极坐标形式)，计算复杂电力系统潮流，具有收敛性好，收敛速度快等优点。所有参数应归算至标幺值下。

/*可计算最大节点数为100，可计算PQ,PV,平衡节点*/ /*可计算非标准变比和平行支路*/ #include #include #include

#define M 100 /*最大矩阵阶数*/ #define Nl 100 /*迭代次数*/

int i,j,k,a,b,c; int t,l;

double P,Q,H,J; int n, m, pq, pv; double eps; double aa[M],bb[M],cc[M],dd[M],max, rr,tt; double mo,c1,d1,c2,d2; double G[M][M],B[M][M],Y[M][M]; 部及其模方值*/

double ykb[M][M],D[M],d[M],dU[M]; */

struct jd { int num,ty; */ double p,q,S,U,zkj,dp,dq,du,dj; 值，电压模值,阻抗角 衡量、电压修正量*/ } jd[M];

struct zl { int numb; int p1,p2; double kx; double r,x; } zl[M]; FILE *fp1,*fp2;

void data() { int h,number;

- /*循环控制变量*/ /*中间变量*/ /*节点数*/ /*支路数*/ /*PQ节点数*/ /*PV节点数*/ /*迭代精度*/

/*中间变量*/

/*复数运算函数的返回值*/ /*节点导纳矩阵中的实部、虚 /*雅克比矩阵、不平衡量矩阵 /*节点结构体*/ /* num为节点号，ty为节点类型 /*节点有功、无功功率，功率模 牛顿--拉夫逊中功率不平 /*支路结构体*/

/*numb为支路号*/ /*支路的两个节点*/ /*非标准变比*/ /*支路的电阻与电抗*/ /* 读取数据 */ .可修编-

. -

fp1=fopen(\"input.txt\

fscanf(fp1,\"%d,%d,%d,%d,%lf\\n\ /*输入节点数,支路数,PQ节点数,PV节点数和迭代精度*/

for(i=1;i<=n;i++) /*输入节点编号、类型、输入功率和电压初值*/ { fscanf(fp1,\"%d,%d\ if(h==1) /*类型h=1是PQ节点*/ { fscanf(fp1,\"%lf,%lf,%lf,%lf\\n\ jd[i].num=number; jd[i].ty=h; } if(h==2) /*类型h=2是pv节点*/ { fscanf(fp1,\ jd[i].num=number; jd[i].ty=h; jd[i].q=-1.567; } if(h==3) /*类型h=3是平衡节点*/ { fscanf(fp1,\ jd[i].num=number; jd[i].ty=h; } }

for(i=1;i<=m;i++) /*输入支路阻抗*/

fscanf(fp1,\"%d,%lf,%d,%d,%lf,%lf\\n\i].x);

fclose(fp1);

if((fp2=fopen(\"output.txt\ { printf(\" can not open file!\\n\"); exit(0); }

- .可修编-

. -

fprintf(fp2,\" 电力系统潮流计算\\n \"); fprintf(fp2,\" ********** 原始数据 *********\\n\");

fprintf(fp2,\"================================================================================\\n\");

fprintf(fp2,\" 节点数:%d 支路数:%d PQ节点数:%d PV节点数:%d 精度:%f\\n\ n,m,pq,pv,eps);

fprintf(fp2,\" ------------------------------------------------------------------------------\\n\"); for(i=1;i<=pq;i++)

fprintf(fp2,\" PQ节点: 节点%d P[%d]=%f Q[%d]=%f\\n\ jd[i].num,jd[i].num,jd[i].p,jd[i].num,jd[i].q); for(i=pq+1;i<=pq+pv;i++)

fprintf(fp2,\" PV节点: 节点%d P[%d]=%f U[%d]=%f 初值Q[%d]=%f\\n\ jd[i].num,jd[i].num,jd[i].p,jd[i].num,jd[i].U,jd[i].num,jd[i].q);

fprintf(fp2,\" 平衡节点: 节点%d e[%d]=%f f[%d]=%f\\n\ jd[n].num,jd[n].num,jd[n].U,jd[n].num,jd[n].zkj);

fprintf(fp2,\" -------------------------------------------------------------------------------\\n\"); for(i=1;i<=m;i++)

fprintf(fp2,\" 支路%d: 相关节点:%d,%d 非标准变比:kx=%f R=%f X=%f \\n\ i,zl[i].p1,zl[i].p2,zl[i].kx,zl[i].r,zl[i].x); fprintf(fp2,\"

==============================================================================\\n\"); }

/*------------------------------------复数运算函数--------------------------------------*/

double mozhi(double a0,double b0) /*复数求模值函数*/ { mo=sqrt(a0*a0+b0*b0); return mo; }

double ji(double a1,double b1,double a2,double b2) /*复数求积函数 a1为电压模值，a2为阻抗角，a3为导纳实部，a4为导纳虚部*/ { a1=a1*cos(b1); b1=a1*tan(b1); c1=a1*a2-b1*b2; d1=a1*b2+a2*b1; return c1; return d1; }

double shang(double a3,double b3,double a4,double b4) /*复数除法求商函数*/ { c2=(a3*a4+b3*b4)/(a4*a4+b4*b4);

- .可修编-

. -

d2=(a4*b3-a3*b4)/(a4*a4+b4*b4); return c2; return d2; }

/*--------------------------------计算节点导纳矩阵----------------------------------*/ void Form_Y() { for(i=1;i<=n;i++) /*节点导纳矩阵元素附初值*/ for(j=1;j<=n;j++) G[i][j]=B[i][j]=0; for(i=1;i<=n;i++) /*节点导纳矩阵的主对角线上的元素,非对地导纳加入相应的主对角线元素中(考虑非标准变比）*/ for(j=1;j<=m;j++) if(zl[j].p1==i) { if(zl[j].kx==1) { mozhi(zl[j].r,zl[j].x); if(mo==0) continue; shang(1,0,zl[j].r,zl[j].x); G[i][i]+=c2; B[i][i]+=d2; } else { mozhi(zl[j].r,zl[j].x); if(mo==0) continue; shang(1,0,zl[j].r,zl[j].x); G[i][i]+=c2/zl[j].kx+c2*(1-zl[j].kx)/(zl[j].kx*zl[j].kx); B[i][i]+=d2/zl[j].kx+d2*(1-zl[j].kx)/(zl[j].kx*zl[j].kx); } } else if(zl[j].p2==i) { if(zl[j].kx==1) { mozhi(zl[j].r,zl[j].x); if(mo==0) continue; shang(1,0,zl[j].r,zl[j].x);

- .可修编-

. -

G[i][i]+=c2; B[i][i]+=d2; } else { mozhi(zl[j].r,zl[j].x); if(mo==0) continue; shang(1,0,zl[j].r,zl[j].x); G[i][i]+=c2/zl[j].kx+c2*(zl[j].kx-1)/zl[j].kx; B[i][i]+=d2/zl[j].kx+d2*(zl[j].kx-1)/zl[j].kx; } }

for(k=1;k<=m;k++) /*节点导纳矩阵非主对角线上（考虑非标准变比）的元素*/ if(zl[k].kx==1) { i=zl[k].p1; j=zl[k].p2; mozhi(zl[k].r,zl[k].x); if(mo==0) continue; shang(1,0,zl[k].r,zl[k].x); G[i][j]-=c2; B[i][j]-=d2; G[j][i]=G[i][j]; B[j][i]=B[i][j]; } else { i=zl[k].p1; j=zl[k].p2; mozhi(zl[k].r,zl[k].x); if(mo==0) continue; shang(1,0,zl[k].r,zl[k].x); G[i][j]-=c2/zl[k].kx; B[i][j]-=d2/zl[k].kx; G[j][i]=G[i][j]; B[j][i]=B[i][j]; }

/*--------------------------输出节点导纳矩阵------------------------------*/

fprintf(fp2,\"\\n\\n ********* 潮流计算过程 *********\\n\");

- .可修编-

. -

fprintf(fp2,\"\\n

==============================================================================\\n\");

fprintf(fp2,\"\\n 节点导纳矩阵为:\"); for(i=1;i<=n;i++) { fprintf(fp2,\"\\n \"); for(k=1;k<=n;k++) { fprintf(fp2,\"%f\ if(B[i][k]>=0) { fprintf(fp2,\"+j\");

fprintf(fp2,\"%f \ } else { B[i][k]=-B[i][k]; fprintf(fp2,\"-j\");

fprintf(fp2,\"%f \ B[i][k]=-B[i][k]; } } }

fprintf(fp2,\"\\n ------------------------------------------------------------------------------\\n\"); }

/*-------------------------------牛顿——拉夫逊-------------------------------*/ void Calculate_Unbalanced_Para() { for(i=1;i<=n;i++) { if(jd[i].ty==1) /*计算PQ节点不平衡量*/ { t=jd[i].num; cc[t]=dd[t]=0; for(j=1;j<=n;j++) { for(a=1;a<=n;a++) { if(jd[a].num==j) break; }

- .可修编-

. -

P=Q=0;

P=jd[a].U*(G[t][j]*cos(jd[i].zkj-jd[a].zkj)+B[t][j]*sin(jd[i].zkj-jd[a].zkj));

Q=jd[a].U*(G[t][j]*sin(jd[i].zkj-jd[a].zkj)-B[t][j]*cos(jd[i].zkj-jd[a].zkj)); cc[t]+=P; dd[t]+=Q; } jd[i].dp=jd[i].p-jd[i].U*cc[t]; jd[i].dq=jd[i].q-jd[i].U*dd[t]; } if(jd[i].ty==2) /*计算PV节点不平衡量*/ { t=jd[i].num; cc[t]=dd[t]=0; for(j=1;j<=n;j++) { for(a=1;a<=n;a++) { if(jd[a].num==j) break; } P=Q=0;

P=jd[a].U*(G[t][j]*cos(jd[i].zkj-jd[a].zkj)+B[t][j]*sin(jd[i].zkj-jd[a].zkj));

Q=jd[a].U*(G[t][j]*sin(jd[i].zkj-jd[a].zkj)-B[t][j]*cos(jd[i].zkj-jd[a].zkj)); cc[t]+=P; dd[t]+=Q; } jd[i].dp=jd[i].p-jd[i].U*cc[t]; jd[i].q=jd[i].U*dd[t]; } }

for(i=1;i<=pq;i++) /*形成不平衡量矩阵D[M]*/ { D[2*i-1]=jd[i].dp; D[2*i]=jd[i].dq; }

for(i=pq+1;i<=n-1;i++) { D[pq+i]=jd[i].dp;

- .可修编-

. -

} fprintf(fp2,\"\\n 不平衡量为:\"); /*输出不平衡量*/ for(i=1;i<=pq;i++) { t=jd[i].num; fprintf(fp2,\"\\n dp[%d]=%f\ fprintf(fp2,\"\\n dq[%d]=%f\ } for(i=pq+1;i<=n-1;i++) { t=jd[i].num; fprintf(fp2,\"\\n dp[%d]=%f\ } }

void Form_Jacobi_Matric() /*形成雅克比矩阵*/ { for(i=1;i<=pq;i++) /*形成pq节点子阵*/ for(j=1;jykb[2*i][2*j-1]=-ykb[2*i-1][2*j]; /* J */ ykb[2*i][2*j]=ykb[2*i-1][2*j-1]; /* L */ } else /*求i=j时的H、N、J、L*/ { aa[i]=0;bb[i]=0; for(k=1;k<=n;k++) { int k1=jd[k].num; H=J=0;

- .可修编-

. -

H=jd[k].U*(G[i1][k1]*sin(jd[i].zkj-jd[k].zkj)-B[i1][k1]*cos(jd[i].zkj-jd[k].zkj)); J=jd[k].U*(G[i1][k1]*cos(jd[i].zkj-jd[k].zkj)+B[i1][k1]*sin(jd[i].zkj-jd[k].zkj)); aa[i]=aa[i]+H; bb[i]=bb[i]+J; }

ykb[2*i-1][2*j-1]=-Ui*(aa[i]-Ui*(G[i1][i1]*sin(jd[i].zkj-jd[i].zkj)-B[i1][i1]*cos(jd[i].zkj-jd[i].zkj))); /*H*/

ykb[2*i][2*j-1]=Ui*(bb[i]-Ui*(G[i1][i1]*cos(jd[i].zkj-jd[i].zkj)+B[i1][i1]*sin(jd[i].zkj-jd[i].zkj))); /*J*/ ykb[2*i-1][2*j]=ykb[2*i][2*j-1]+2*Ui*Ui*G[i1][i1]; /*N*/ ykb[2*i][2*j]=-ykb[2*i-1][2*j-1]-2*Ui*Ui*B[i1][i1]; /*L*/ } } else /*求j>pq时的H、J */ { ykb[2*i-1][pq+j]=Ui*Uj*(G[i1][j1]*sin(zi-zj)-B[i1][j1]*cos(zi-zj)); /* H */ ykb[2*i][pq+j]=-Ui*Uj*(G[i1][j1]*cos(zi-zj)+B[i1][j1]*sin(zi-zj)); /* J */ } }

for(i=pq+1;i<=n-1;i++) /*形成pv节点子阵*/

for(j=1;jpq时的H*/

- .可修编-

. -

{ if(i!=j) /*求i!=j时的H*/ ykb[pq+i][pq+j]=Ui*Uj*(G[i1][j1]*sin(zi-zj)-B[i1][j1]*cos(zi-zj)); /* H */ else /*求i=j时的H*/ { aa[i]=0; for(k=1;k<=n;k++) { int k1=jd[k].num; H=0;

H=jd[k].U*(G[i1][k1]*sin(jd[i].zkj-jd[k].zkj)-B[i1][k1]*cos(jd[i].zkj-jd[k].zkj)); aa[i]=aa[i]+H; }

ykb[pq+i][pq+j]=-Ui*(aa[i]-Ui*(G[i1][i1]*sin(jd[i].zkj-jd[i].zkj)-B[i1][i1]*cos(jd[i].zkj-jd[i].zkj))); /*H*/ } } } /*------------------------------输出雅克比矩阵--------------------------------*/ fprintf(fp2,\"\\n\\n 雅克比矩阵为: \"); for(i=1;i<=(2*pq+pv);i++) {

fprintf(fp2,\"\\n\"); for(j=1;j<=2*pq+pv;j++) { fprintf(fp2,\" %f\ } } }

void Solve_Equations() /* 求解修正方程组 (LU分解法)*/ {

double l[Nl][Nl]={0}; //定义L矩阵 double u[Nl][Nl]={0}; //定义U矩阵 double y[Nl]={0}; //定义数组Y double x[Nl]={0}; //定义数组X

double a[Nl][Nl]={0}; //定义系数矩阵 double b[Nl]={0}; //定义右端项 double sum=0; int i,j,k,s;

- .可修编-

. -

int n;

n=2*pq+pv;

for(i=0; ifor(j=0; ja[i][j]=ykb[i+1][j+1]; } }

for(i=0; ib[i]=D[i+1]; }

for(i=0; ifor(j=0; jif(i==j) l[i][j] = 1; } }

for(i=0;i{ u[0][i]=(float)(a[0][i]);} /*第一步：对矩阵U的首行进行计算*/

for(k=0;ksum+=l[i][s]*u[s][k]; }

l[i][k]=(float)((a[i][k]-sum)/u[k][k]); }

for(j=k+1;jsum+=l[k+1][s]*u[s][j]; }

u[k+1][j]=(float)((a[k+1][j]-sum));

- .可修编-

. -

} }

y[0]=b[0] ; /*回代法计算数组Y*/ for(i=1;i{ for(j=0,sum=0;jy[i]=(float)(b[i]-sum); }

x[n-1]=(float)(y[n-1]/u[n-1][n-1]); /*回代法计算数组X*/ for(i=n-2;i>=0;i--)

{ for(j=n-1,sum=0;j>i;j--) { sum+=x[j]*u[i][j];}

x[i]=(float)((y[i]-sum)/u[i][i]); }

for(i=1; i<=n; i++) {

d[i]=x[i-1];

}

max=fabs(d[1]); /*选出最大的修正量的值*/ for(i=1;i<=n;i++)

if((fabs(d[i]))>max) max=fabs(d[i]); }

void Niudun_Lafuxun() { int z=1; fprintf(fp2,\"\\n --------------极坐标形式的牛顿-拉夫逊迭代法结果:--------------\\n\"); do { max=1; if((z=eps)) { fprintf(fp2,\"\\n 迭代次数: %d\\n\ /*开始迭代计算*/ } Calculate_Unbalanced_Para(); Form_Jacobi_Matric();

- .可修编-

. -

Solve_Equations(); for(i=1;i<=pq;i++) { jd[i].zkj+=d[2*i-1]; jd[i].U+=d[2*i]*jd[i].U; }

for(i=pq+1;i<=n-1;i++) { jd[i].zkj+=d[pq+i]; }

fprintf(fp2,\"\\n\\n 输出 dδ,dU: \"); /*输出修正量的值*/ for(c=1;c<=n;c++) { for(a=1;a<=n;a++) { if(jd[a].num==c) break; } fprintf(fp2,\"\\n\"); if(jd[a].ty==1)

fprintf(fp2,\" 节点为 %2d dδ=%8.5f dU=%8.5f\ if(jd[a].ty==2)

fprintf(fp2,\" 节点为 %2d dδ=%8.5f\ } fprintf(fp2,\"\\n\\n 输出迭代过程中的电压值: \"); for(c=1;c<=n;c++) { for(a=1;a<=n;a++) { if(jd[a].num==c) break; } fprintf(fp2,\"\\n U[%d]=%f\ fprintf(fp2,\"∠%f\ } fprintf(fp2,\"\\n\\n ------------------------------------------------------------------------------\"); z++; } while((z=eps)); /*判断是否达到精度要求*/ }

void Powerflow_Result()

- .可修编-

. -

{ int n1=jd[n].num; fprintf(fp2,\"\\n\\n

==============================================================================\\n\\n\"); fprintf(fp2,\" ******潮流计算结果******\"); fprintf(fp2,\"\\n\\n

==============================================================================\\n\\n\"); fprintf(fp2,\"\\n 各节点电压值: \"); /*输出各节点的电压值*/ for(c=1;c<=n;c++) { for(a=1;a<=n;a++) { if(jd[a].num==c) break; } fprintf(fp2,\"\\n U[%d]=%f\ fprintf(fp2,\"∠%f\ } rr=tt=0; /*计算节点的注入功率*/ for(i=1;i<=n;i++) { int i1=jd[i].num; ji(jd[i].U,-jd[i].zkj,G[n1][i1],-B[n1][i1]); rr+=c1; tt+=d1; } ji(jd[n].U,jd[n].zkj,rr,tt); fprintf(fp2,\"\\n\\n 各节点注入功率:\\n\"); for(i=1;i<=pq;i++) { int i1=jd[i].num; fprintf(fp2,\" PQ节点: 节点%d S[%d]=%f\ i1,i1,jd[i].p); /*PQ节点注入功率*/ if(jd[i].q>=0) fprintf(fp2,\"+j%f\\n\ else fprintf(fp2,\"-j%f\\n\ }

for(i=pq+1;i<=n-1;i++) {

- .可修编-

. -

int i1=jd[i].num; fprintf(fp2,\" PV节点: 节点%d S[%d]=%f\ i1,i1,jd[i].p); /*PV节点注入功率*/ if(jd[i].q>=0) fprintf(fp2,\"+j%f\\n\ else fprintf(fp2,\"-j%f\\n\ }

fprintf(fp2,\" 平衡节点: 节点%d\ /*平衡节点注入功率*/ fprintf(fp2,\" S[%d]=%f\ if(d1>=0) fprintf(fp2,\"+j%f\ else fprintf(fp2,\"-j%f\ fprintf(fp2,\"\\n\\n 线路功率:\\n\"); rr=tt=0;

for(i=1;i<=m;i++) { int i1=zl[i].p1; /*计算线路功率*/ int j1=zl[i].p2; aa[i]=bb[i]=P=Q=0; for(a=1;a<=n;a++) { if(jd[a].num==i1) break; } for(b=1;b<=n;b++) { if(jd[b].num==j1) break; } mozhi(zl[i].r,zl[i].x); if(mo==0) continue; shang(1,0,zl[i].r,zl[i].x); ji(jd[a].U/zl[i].kx/zl[i].kx,-jd[a].zkj,c2,-d2); /*考虑非标准变比*/ P+=c1;Q+=d1; ji(jd[b].U/zl[i].kx,-jd[b].zkj,c2,-d2); P-=c1;Q-=d1; ji(jd[a].U,jd[a].zkj,P,Q); fprintf(fp2,\" 线路%d: %d--%d 的功率为: %f\ if(d1>=0)

- .可修编-

. -

fprintf(fp2,\"+j%f\\n\ else fprintf(fp2,\"-j%f\\n\ aa[i]+=c1; bb[i]+=d1; P=Q=0; ji(jd[b].U,-jd[b].zkj,c2,-d2); /*考虑非标准变比*/ P+=c1;Q+=d1; ji(jd[a].U/zl[i].kx,-jd[a].zkj,c2,-d2); P-=c1;Q-=d1; ji(jd[b].U,jd[b].zkj,P,Q); fprintf(fp2,\" 线路%d: %d--%d 的功率为: %f\ if(d1>=0) fprintf(fp2,\"+j%f\\n\ else fprintf(fp2,\"-j%f\\n\ aa[i]+=c1; bb[i]+=d1; rr+=aa[i]; tt+=bb[i]; }

fprintf(fp2,\"\\n\\n 线路损耗功率:\\n\"); /*计算线路功率损耗*/ for(i=1;i<=m;i++) { int i1=zl[i].p1; int j1=zl[i].p2; fprintf(fp2,\" 线路%d损耗的功率为: %f\ if(bb[i]>=0) fprintf(fp2,\"+j%f\\n\ else fprintf(fp2,\"-j%f\\n\ } fprintf(fp2,\"\\n\\n 网络总损耗功率为: %f\ /*计算网络总损耗*/ if(tt>=0) fprintf(fp2,\"+j%f\\n\ else fprintf(fp2,\"-j%f\\n\ fprintf(fp2,\"\\n

==============================================================================\\n\");

fprintf(fp2,\"\\n\\n ********* 潮流计算结束 *********\"); }

- .可修编-

. -

void main() { printf(\"请仔细阅读附录里的输入文件使用说明以便于你正确输入!\\n\"); data(); /*读取数据*/

Form_Y(); /*形成节点导纳矩阵*/ for(i=1;i<=pq;i++) /* U、zkj附初值*/ { jd[i].U=1; jd[i].zkj=0; }

for(i=pq+1;i附录：

输入文件按以下规则输入

节点数，支路数，PQ节点数，PV节点数，精度

节点编号，节点类型（1为PQ节点，2为PV节点，3为平衡节点） 节点输入有功功率，无功功率(节点电压纵分量，横分量)

支路编号，非标准变比（若没有则填1即可），相关节点，相关节点，支路电阻，支路电抗

- .可修编-