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

所有参数应归算至标幺值下。

/*可计算最大节点数为100，可计算PQ,PV,平衡节点*//*可计算非标准变比和平行支路*/#include<stdio.h>#include<math.h>#include<stdlib.h>#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, /*PQ节点数*/pv; /*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; /* num为节点号，ty为节点类型*/ double p,q,S,U,zkj,dp,dq,du,dj; /*节点有功、无功功率，功率模值，电压模值,阻抗角牛顿--拉夫逊中功率不平衡量、电压修正量*/} jd[M];struct zl /*支路结构体*/{ int numb; /*numb为支路号*/int p1,p2; /*支路的两个节点*/double kx; /*非标准变比*/double r,x; /*支路的电阻与电抗*/ } zl[M];FILE *fp1,*fp2;void data() /* 读取数据*/{int h,number;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(" 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",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");}/*-------------------------------牛顿——拉夫逊-------------------------------*/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",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,"\n dp[%d]=%f",i,D[pq+t]);}}void Form_Jacobi_Matric() /*形成雅克比矩阵*/{for(i=1;i<=pq;i++) /*形成pq节点子阵*/for(j=1;j<n;j++){ int i1=jd[i].num;int j1=jd[j].num;double Ui=jd[i].U;double Uj=jd[j].U;double zi=jd[i].zkj;double zj=jd[j].zkj;if(j<=pq) /*求j<=pq时的H、N、J、L */{if(i!=j) /*求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 */ykb[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;j<n;j++){int i1=jd[i].num;int j1=jd[j].num;double Ui=jd[i].U;double Uj=jd[j].U;double zi=jd[i].zkj;double zj=jd[j].zkj;if(j<=pq) /*求j<=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++){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",ykb[i][j]);}}}void Solve_Equations() /* 求解修正方程组(LU分解法)*/ {double l[Nl][Nl]={0}; //定义L矩阵double u[Nl][Nl]={0}; //定义U矩阵double y[Nl]={0}; //定义数组Ydouble x[Nl]={0}; //定义数组Xdouble 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; 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++) /*初始化矩阵l*/{for(j=0; j<n; j++){if(i==j) l[i][j] = 1;}}for(i=0;i<n;i++) /*开始LU分解*/{ u[0][i]=(float)(a[0][i]);} /*第一步：对矩阵U的首行进行计算*/for(k=0;k<n-1;k++) /*第二步：逐步进行LU分解*/ { for(i=k+1;i<n;i++) /*对L的第k列进行计算*/{ for(s=0,sum=0;s<n;s++){ if(s!=k)sum+=l[i][s]*u[s][k];}l[i][k]=(float)((a[i][k]-sum)/u[k][k]);}for(j=k+1;j<n;j++) /*对U的第k+1行进行计算*/ { for(s=0,sum=0;s<n;s++){ if(s!=k+1)sum+=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<n;i++){ for(j=0,sum=0;j<i;j++){ sum+=y[j]*l[i][j];}y[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<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 输出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",c,d[2*a-1],d[2*a]);if(jd[a].ty==2)fprintf(fp2," 节点为%2d dδ=%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,"\n U[%d]=%f",c,jd[a].U);fprintf(fp2,"∠%f",jd[a].zkj);}fprintf(fp2,"\n\n ------------------------------------------------------------------------------");z++;} while((z<Nl)&&(max>=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",c,jd[a].U);fprintf(fp2,"∠%f",jd[a].zkj);}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",jd[i].q);elsefprintf(fp2,"-j%f\n",-jd[i].q);}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",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++){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",i,i1,j1,c1);if(d1>=0)fprintf(fp2,"+j%f\n",d1);elsefprintf(fp2,"-j%f\n",-d1);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",i,j1,i1,c1);if(d1>=0)fprintf(fp2,"+j%f\n",d1);elsefprintf(fp2,"-j%f\n",-d1);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",i,aa[i]);if(bb[i]>=0)fprintf(fp2,"+j%f\n",bb[i]);elsefprintf(fp2,"-j%f\n",-bb[i]);}fprintf(fp2,"\n\n 网络总损耗功率为: %f",rr); /*计算网络总损耗*/ if(tt>=0)fprintf(fp2,"+j%f\n",tt);elsefprintf(fp2,"-j%f\n",-tt);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<n;i++){jd[i].U=jd[i].U; jd[i].zkj=0;}Niudun_Lafuxun(); /*牛顿--拉夫逊迭代*/Powerflow_Result();printf("潮流计算结果存放于output.txt文件中!\n");}附录：输入文件按以下规则输入节点数，支路数，PQ节点数，PV节点数，精度节点编号，节点类型（1为PQ节点，2为PV节点，3为平衡节点）节点输入有功功率，无功功率(节点电压纵分量，横分量)支路编号，非标准变比（若没有则填1即可），相关节点，相关节点，支路电阻，支路电抗。