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

/*最大矩阵阶数*//*迭代次数*//*循环控制变量*//*最大矩阵阶数*//*迭代次数*//*循环控制变量*//*中间变量*//*节点数*//*支路数*//*PQ节点数*//*PV节点数*//*迭代精度*//*中间变量*//*复数运算函数的返回值*//*节点导纳矩阵中的实部、虚/*雅克比矩阵、不平衡量矩阵/*节点结构体*//*num为节点号，ty为节点类型/*节点有功、无功功率，功率模牛顿--拉夫逊中功率不平/*利用牛顿-拉夫逊迭代法（极坐标形式），计算复杂电力系统潮流，具有收敛性好，收敛速度快等优点。所有参数应归算至标幺值下。/*可计算最大节点数为100,可计算PQ，PV,平衡节点*//*可计算非标准变比和平行支路*/#include<stdio.h>#include<math.h>#include<stdlib.h>#defineM100#defineNl100inti,j,k,a,b,c;intt,l;doubleP,Q,H,J;intn,m,pq,pv;doubleeps;doubleaa[M],bb[M],cc[M],dd[M],max,rr,tt;doublemo,c1,d1,c2,d2;doubleG[M][M],B[M][M],Y[M][M];部及其模方值*/doubleykb[M][M],D[M],d[M],dU[M];*/structjd{intnum,ty;*/doublep,q,S,U,zkj,dp,dq,du,dj;值，电压模值,阻抗角衡量、电压修正量*/}jd[M];structzl{intnumb;intp1,p2;doublekx;doubler,x;}zl[M];FILE*fp1,*fp2;voiddata(){inth,number;/*支路结构体*//*numb为支路号*//*支路的两个节点*//*非标准变比*//*支路的电阻与电抗*//*读取数据*/fp1=fopen("input.txt","r");fscanf(fp1,"%d,%d,%d,%d,%lf\n",&n,&m,&pq,&pv,&eps); /*输入节点数,支路数,PQ节点数,PV节点数和迭代精度*/for(i=1;i<=n;i++) /*输入节点编号、类型、输入功率和电压初值*/{fscanf(fp1,"%d,%d",&number,&h);if(h==1) /*类型h=1是PQ节点*/{fscanf(fp1,"%lf,%lf,%lf,%lf\n",&jd[i].p,&jd[i].q,&jd[i].U,&jd[i].zkj);jd[i].num=number;jd[i].ty=h;}if(h==2) /*类型h=2是pv节点*/{fscanf(fp1,",%lf,%lf,%lf\n",&jd[i].p,&jd[i].U,&jd[i].zkj);jd[i].num=number;jd[i].ty=h;jd[i].q=-1.567;}if(h==3)/*类型h=3是平衡节点*/{fscanf(fp1,",%lf,%lf\n",&jd[i].U,&jd[i].zkj);jd[i].num=number;jd[i].ty=h;}}for(i=1;i<=m;i++) /*输入支路阻抗*/fscanf(fp1,"%d,%lf,%d,%d,%lf,%lf\n",&zl[i].numb,&zl[i].kx,&zl[i].p1,&zl[i].p2,&zl[i].r,&zl[i].x);fclose(fp1);if((fp2=fopen("output.txt","w"))==NULL){printf("cannotopenfile!\n");exit(0);}

fprintf(fp2,"fprintf(fp2,"原始数据*\n");电力系统潮流计算\nfprintf(fp2,"fprintf(fp2,"原始数据*\n");fprintf(fp2,"================================================================================\n");fprintf(fp2," 节点数：％d支路数:％dPQ节点数:％dPV节点数:％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=%fR=%fX=%f\n",i,zl[i].p1,zl[i].p2,zl[i].kx,zl[i].r,zl[i].x);fprintf(fp2,"=============\n");}/* 复数运算函数 */doublemozhi(doublea0,doubleb0) /*复数求模值函数*/{ mo=sqrt(a0*a0+b0*b0);returnmo;}doubleji(doubleo1,doubleb1,doubleo2,doubleb2) /*复数求积函数a1为电压模值，a2为阻抗角，a3为导纳实部，a4为导纳虚部*/{ a1=a1*cos(b1);b1=a1*tan(b1);c1=a1*a2-b1*b2;d1=a1*b2+a2*b1;returnc1;returnd1;}doubleshang(doublea3,doubleb3,doublea4,doubleb4)/*复数除法求商函数*/{ c2=(a3*a4+b3*b4)/(a4*a4+b4*b4);d2=(a4*b3-a3*b4)/(a4*a4+b4*b4);returnc2;returnd2;}/* 计算节点导纳矩阵 */voidForm_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);}}elseif(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",G[i][k]);if(B[i][k]>=0){fprintf(fp2,"+j");fprintf(fp2,"%f",B[i][k]);}else{B[i][k]=-B[i][k];fprintf(fp2,"-j");fprintf(fp2,"%f",B[i][k]);B[i][k]=-B[i][k];}}}fprintf(fp2,"\n \n");}/* 牛顿——拉夫逊 */voidCalculate_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",i,D[2*t-1]);fprintf(fp2,"\n dq[%d]=%f",i,D[2*t]);}for(i=pq+1;i<=n-1;i++){t=jd[i].num;fprintf(fp2,"\ndp[%d]=%f",i,D[pq+t]);}}voidForm_Jacobi_Matric(){for(i=1;i<=pq;i++)for(j=1;j<n;j++){voidForm_Jacobi_Matric(){for(i=1;i<=pq;i++)for(j=1;j<n;j++){inti1=jd[i].num;intj1=jd[j].num;doubleUi=jd[i].U;doubleUj=jd[j].U;doublezi=jd[i].zkj;doublezj=jd[j].zkj;if(j<=pq){/*形成雅克比矩阵*//*形成pq节点子阵*//*求j<=pq时的H、N、J、L*//*求/*求i!=j时的H、N、J、L*/{ykb[2*i-1][2*j-1]=Ui*Uj*(G[i1][j1]*sin(zi-zj)-B[i1][j1]*cos(zi-zj));/*H*/ykb[2*i-1][2*j]=Ui*Uj*(G[i1][j1]*cos(zi-zj)+B[i1][j1]*sin(zi-zj));/*N*//*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++){intk1=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;j<n;j++){inti1=jd[i].num;intj1=jd[j].num;doubleUi=jd[i].U;doubleUj=jd[j].U;doublezi=jd[i].zkj;doublezj=jd[j].zkj;if(jv=pq) /*求jv=pq时的H、N*/{ykb[pq+i][2*j-1]=Ui*Uj*(G[i1][j1]*sin(zi-zj)-B[i1][j1]*cos(zi-zj)); /*H*/ykb[pq+i][2*j]=Ui*Uj*(G[i1][j1]*cos(zi-zj)+B[i1][j1]*sin(zi-zj)); /*N*/}else /*求j>pq时的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++){intk1=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",ykb[i][j]);}}}voidSolve_Equations()/*求解修正方程组(LU分解法)*/{doublel[NI][NI]={0};〃定义L矩阵doubleu[NI][NI]={0};〃定义U矩阵doubley[NI]={0};〃定义数组丫doublex[NI]={0};〃定义数组XdoubIea[NI][NI]={0};//定义系数矩阵doubIeb[NI]={0};//定义右端项doubIesum=0;inti,j,k,s;

intn;n=2*pq+pv;for(i=0;i<n;i++){for(j=0;j<n;j++){a[i][j]=ykb[i+1][j+1];}}for(i=0;i<n;i++){b[i]=D[i+1];}for(i=0;i<n;i++){for(i=0;i<n;i++){for(j=0;j<n;j++){if(i==j)l[i][j]=1;}/*初始化矩阵I*/for(i=0;i<n;i++){u[0][i]=(fIoat)(a[0][i]);}for(i=0;i<n;i++){u[0][i]=(fIoat)(a[0][i]);}/*开始LU分解*//*第一步：对矩阵U的首行进行计算*/for(k=0;kvn-1;k++) /*第二步：逐步进行LU分解*/{for(i=k+1;ivn;i++) /*对L的第k列进行计算*/{for(s=0,sum=0;svn;s++){if(s!=k)sum+=I[i][s]*u[s][k];}I[i][k]=(fIoat)((a[i][k]-sum)/u[k][k]);}for(j=k+1;jvn;j++) /*对U的第k+1行进行计算*/{for(s=0,sum=0;svn;s++){if(s!=k+1)sum+=I[k+1][s]*u[s][j];}u[k+1][j]=(fIoat)((a[k+1][j]-sum));y[O]=b[O]; /*回代法计算数组Yfor(i=1;i<n;i++){for(j=O,sum=O;j<i;j++){sum+=y[j]*l[i][j];}y[i]=(float)(b[i]-sum);}x[n-1]=(floot)(y[n-1]/u[n-1][n-1]);/*回代法计算数组X*/for(i=n-2;i>=O;i--){for(j=n-1,sum=O;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]);}voidNiudun_Lafuxun(){intz=1;fprintf(fp2,"\n 极坐标形式的牛顿-拉夫逊迭代法结果: \n");do{max=1;if((z<Nl)&&(max>=eps)){fprintf(fp2,"\n迭代次数:%d\n",z); /*开始迭代计算*/}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输出d6,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 d6dU=%8.5f",c,d[2*a-1],d[2*a])；if(jd[a].ty==2)fprintf(fp2,"节点为%2d d6=%8.5f",c,d[pq+a])；}fprintf(fp2,"\n\n输出迭代过程中的电压值:")；for(c=1；c<=n；c++){for(a=1；a<=n；a++){if(jd[a].num==c)break；}fprintf(fp2,"\nU[%d]=%f",c,jd[a].U)；fprintf(fp2,"Z%f",jd[o].zkj);}fprintf(fp2,"\n\n z++;}while((z<Nl)&&(max>=eps)); /*判断是否达到精度要求*/=%8.5f");}voidPowerflow_Result()=%8.5f");{intn1=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,"\nU[%d]=%f",c,jd[a].U);fprintf(fp2,"Z%f",jd[o].zkj);}rr=tt=0; /*计算节点的注入功率*/for(i=1;i<=n;i++){inti1=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++){inti1=jd[i].num;fprintf(fp2,"PQ节点： 节点％dS[%d]=%f",i1,i1,jd[i].p);/*PQ节点注入功率*/if(jd[i].q>=0)fprintf(fp2,"+j%f\n",jd[i].q);elsefprintf(fp2,"-j%f\n",-jd[i].q);}for(i=pq+1;i<=n-1;i++){inti1=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",jd[i].q);elsefprintf(fp2,"-j%f\n",-jd[i].q);}fprintf(fp2," 平衡节点：节点％d",jd[n].num,jd[n].num);/*平衡节点注入功率*/fprintf(fp2," S[%d]=%f",n1,c1);if(d1>=0)fprintf(fp2,"+j%f",d1);elsefprintf(fp2,"-j%f",-d1);fprintf(fp2,"\n\n 线路功率:\n");rr=tt=0;for(i=1;i<=m;i++){inti1=zl[i].p1; /*计算线路功率*/intj1=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