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

1、牛顿-拉夫逊迭代法极坐标潮流 计算C语言程序/*利用牛顿-拉夫逊迭代法（极坐标形式），计算复 杂电力系统潮流，具有收敛性好，收敛速度快等 优点。所有参数应归算至标幺值下。/*可计算最大节点数为100，可计算PQ,PV,平衡 节点*/*可计算非标准变比和平行支路*/#in clude<stdio.h>#in clude<math.h>#in clude<stdlib.h>#defi ne100/*最大矩阵阶数*/#defi neNl/*迭代次数*/100int/*循环控制变量*/int t,l;doublei,j,k,a,b,c;PQH,J;/*中间变量*/i

2、nt/*节点数*/n,m,/*支路数*/pq,/*PQ节点数*/pv；/*PV节点数*/doubleeps;/*迭代精度*/double aaM,bbM,ccM,ddM,max, rr,tt;/*中间变量*/doublemo,c1,d1,c2,d2;/*复数运算函数的返回值*/doubleGMM,BMM, YMM;/*节点导纳矩阵中的实部、虚部及其模方值*/doubleykbMM,DM,dM,dUM;/*雅克比矩阵、不平衡量矩阵*/structjd/*节点结构体*/intnu m,ty;/* num为节点号，ty为节点类型*/doublep,q,S,U,zkj,dp,dq,du,dj;/*节点

3、有功、无功功率，功率模值，电压模值,阻抗角牛顿-拉夫逊中功率不平衡量、电压修正量*/structzl/*支路结构体*/intnumb;/*numb为支路号*/intp1,p2；/*支路的两个节点*/doublekx;/*非标准变比*/doubler,x;/*支路的电阻与电抗*/ zlM;FILE *fp1,*fp2;voiddata()/*读取数据*/int h,n umber;fp仁 fope n("in put.txt","r");fsca nf(fp1,"%d,%d,%d,%d,%lfn",&n,&m,&

4、pq,&pv,&eps); /*输入节点数，支路数,PQ节点 数,PV节点数和迭代精度*/for(i=1;i<=n ;i+)/*输入节点编号、类型、输入功率和电压初值*/fscan f(fp1,"%d,%d",&number，&h);if(h=1)/*类型h=1是PQ节点*/fscan f(fp1,"%lf,%lf,%lf,%lfn",&jdi.p,&jdi.q,&jdi.U,&jdi.zkj);jdi. num=nu mber;jdi.ty=h;if(h=2)/*类型h=2是pv节点

5、*/fsca nf(fp1,",%lf,%lf,%lfn",&jdi.p,&jdi.U，&jdi.zkj);jdi. num=n umber;jdi.ty=h;jdi.q=-1.567;if(h=3)/*类型h=3是平衡节点*/fsca nf(fp1,",%lf,%lfn",&jdi.U,&jdi.zkj);jdi. num=n umber;jdi.ty=h;for(i=1;i<=m;i+)/*输入支路阻抗*/fsca nf(fp1,"%d,%lf,%d,%d,%lf,%lfn", &am

6、p;zli. numb, &zli.kx, &zli.p1, &zli.p2, &zli.r, &zli .x);fclose(fp1); if(fp2=fope n("output.txt","w")=NULL) prin tf(" can not ope n file!' n");exit(0);fprin tf（fp2,"电力系统潮流计算n "）;fprin tf(fp2,"*原始数据 *n");fprin tf(fp2,"= =n

7、");fprintf(fp2,"节点数:d 支路数:dPQ节点数:%d PV 节点数:%d精度:%fn",n,m,pq,pv,eps);fprin tf(fp2,"n");for(i=1;i<=pq;i+)fprintf(fp2,"PQ 节点：节点 %dP%d=%fQ%d=%fn",jdi.nu m,jdi. nu m,jdi.p,jdi. num,jdi.q);for(i=pq+1;i<=pq+pv;i+)fprintf(fp2,"PV 节点：节点 %dP%d=%fU%d=%f初值Q%d=%fn&qu

8、ot;,jdi.n um,jdi. nu m,jdi.p,jdi. num,jdi.U,jdi.n um,jdi.q);fprintf(fp2,"平衡节点： 节点%de%d=%ff%d=%fn",jdn . num,jd n. num,jd n.U,jd n. num,jd n.zkj);fprin tf(fp2,"n");for(i=1;i<=m;i+)fprintf(fp2,"支路 %d:相关节点:d,%d非标准变比：kx=%f R=%f X=%fn",i,zli.p1,zli.p2,zli.kx,zli.r,zli.x);

9、fprin tf(fp2,"=n")；/* 复数运算函数*/double mozhi(double aO,double b0) /*复数求模值函数*/ mo=sqrt(aO*aO+bO*bO);return mo;double ji(double a1,double b1,double a2,double b2)/*复数求积函数 a1为电压模值，a2为阻抗角，a3为导纳实部，a4为导纳虚部*/a1=a1*cos(b1);b1= a1*ta n( b1);c仁 a1*a2-b1*b2;d1=a1*b2+a2*b1;return c1;return d1;double sha

10、ng(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+)Gij=Bij=O；for(i=1;i<=n ;i+)/*节点导纳矩阵的主对角线上的元素 ，非对地导 纳加入相应的主对角线元素中(考虑非标准变 比)*/for(

11、j=1;j<=m;j+)if(zlj.p1=i)if(zlj.kx=1)mozhi(zlj.r,zlj.x); if(mo=0) continue;sha ng(1,0,zlj.r,zlj.x);Gii+=c2;Bii+=d2;elsemozhi(zlj.r,zlj.x);if(mo=0) continue;sha ng(1,0,zlj.r,zlj.x);Gii+=c2/zlj.kx+c2*(1-zlj.kx)/(zlj.kx*zlj. kx);Bii+=d2/zlj.kx+d2*(1-zlj.kx)/(zlj.kx*zlj .kx);else if(zlj.p2=i)if(zlj.kx

12、=1)mozhi(zlj.r,zlj.x);if(mo=0)con ti nue;sha ng(1,0,zlj.r,zlj.x);Gii+=c2;Bii+=d2;else mozhi(zlj.r,zlj.x); if(mo=0) continue;sha ng(1,0,zlj.r,zlj.x);Gii+=c2/zlj.kx+c2*(zlj.kx-1)/zlj.kx;Bii+=d2/zlj.kx+d2*(zlj.kx-1)/zlj.kx;for(k=1;k<=m;k+)/*节点导纳矩阵非主对角线上(考虑非标准变 比)的元素*/if(zlk.kx=1)i=zlk.p1;j=zlk.p2;mo

13、zhi(zlk.r,zlk.x);if(mo=0) continue; sha ng(1,0,zlk.r,zlk.x);Gij-=c2;Bij-=d2;Gji=Gij;Bji=Bij;elsei=zlk.p1;j=zlk.p2;mozhi(zlk.r,zlk.x); if(mo=0) continue; sha ng(1,0,zlk.r,zlk.x);Gij-=c2/zlk.kx;Bij-=d2/zlk.kx;Gji=Gij;Bji=Bij;/*输出节点导纳矩阵*/fprin tf(fp2,"nn*潮流计算过程*n");fprin tf(fp2,"n=n"

14、;)；fprintf(fp2,"n节点导纳矩阵为:")；for(i=1;i<=n ;i+)fprin tf(fp2,"n");for(k=1;k<=n ;k+)fprin tf(fp2,"%f",Gik);if(Bik>=0)fpri ntf(fp2,"+j");fprintf(fp2,"%f ",Bik);elseBik=-Bik;fprin tf(fp2,"-j"); fprintf(fp2,"%f Bik=-Bik;fprin tf(fp2,

15、"nn");/* 牛顿-*/void Calculatenbala nced_Para() for(i=1;i<=n ;i+)if(jdi.ty=1)节点不平衡量*/",Bik);拉夫逊/*计算PQt=jdi. num; cct=ddt=0;(21 二二 PN宀Eppn.曰 PTbspHbpsprEoyn.曰 PTdspHdpspf宀 oH + mpp CLH+S。 -(nzsphmzspdso。*巨&(nzsphmz-e:pdu0曰口0)*nspHosmzsphmzspdufm&+(NZsPHMZ-E:PDS8sE0)*nspHdoHOHd

16、三 eaiq (FEnu.巨 PE (+ecHV2L&04) (+ruHvr lhd04算PV节点不平衡量*/t=jdi. num;cct=ddt=O;for(j=1;j<=n ;j+)for(a=1;a<=n; a+)if(jda. num=j)break;P=Q=O;P=jda.U*(Gtj*cos(jdi.zkj-jda.zkj)+Btj *si n(jdi.zkj-jda.zkj);Q=jda.U*(Gtj*si n(jdi.zkj-jda.zkj)-Btj*cos(jd【i.zkj-jd【a.zkj);cct+=P;ddt+=Q;jd【i.dp=jd【i.p-jd

17、【i.U*cc【t;jdi.q=jdi.U*ddt;for(i=1;i<=pq;i+)/*形成不平衡量矩阵D【M*/D2*i-1=jdi.dp;D2*i=jdi.dq;for(i=pq+1;i<=n-1;i+)Dpq+i=jdi.dp;fprintf(fp2,"n不平衡量为:");/*输出不平衡量*/for(i=1;i<=pq;i+)t=jdi. num;fprin tf(fp2,"ndp%d=%f",i,D2*t-1);fprin tf(fp2,"ndq%d=%f",i,D2*t);for(i=pq+1;i<

18、=n-1;i+)t=jdi. num;fprin tf(fp2,"ndp%d=%f",i,Dpq+t);/*形成/*形成pqvoid Form_Jacobi_Matric()雅克比矩阵*/for(i=1;i<=pq;i+)节点子阵*/for(j=1;j< n;j+) int i1=jdi. num;int j1=jdj. num;double Ui=jdi.U;double Uj=jdj.U;double zi=jdi.zkj;double zj=jdj.zkj;if(jv=pq)/* 求jv=pq 时的 h、N、J、L */if(i!=j)/* 求 i!=j

19、时的 H、N、J、L*/ykb2*i-12*j-1=Ui*Uj*(Gi1j1*si n( zi-zj)-Bi1j1*cos(zi-zj);/* H */ykb2*i-12*j=Ui*Uj*(Gi1j1*cos(zi-zj)+Bi1j1*si n( zi-zj);/* N */ykb2*i2*j-1=-ykb2*i-12*j;/* J */ykb2*i2*j=ykb2*i-12*j-1;/* L */else/* 求i=j 时的 H、N、J、L*/aai=O;bbi=O;for(k=1;k<=n ;k+)int k1=jdk. num;H=J=O;H=jdk.U*(Gi1k1*si n(j

20、di.zkj-jdk.zkj)-Bi 1k1*cos(jdi.zkj-jd【k.zkj);J=jdk.U*(Gi1k1*cos(jdi.zkj-jdk.zkj)+Bi1k1*si n(jdi.zkj-jdk.zkj);aai=aai+H; bbi=bbi+J;ykb2*i-12*j-1=-Ui*(aai-Ui*(Gi1i1*si n(jd i.zkj-jdi.zkj)-Bi1i1*cos(jdi.zkj-jdi.zkj) );/*H*/ykb2*i2*j-1=Ui*(bbi-Ui*(Gi1i1*cos(jdi .zkj-jdi.zkj)+Bi1i1*s in (jdi.zkj-jdi.zkj)

21、;/*J*/ ykb2*i-12*j=ykb2*i2*j-1+2*Ui*Ui*Gi1i1;/*N*/ykb2*i2*j=-ykb2*i-12*j-1-2*Ui*Ui*Bi1i1;/*L*/else/* 求j>pq 时的 H、J */ykb2*i-1pq+j=Ui*Uj*(Gi1j1*si n( zi-zj)-Bi1j1*cos(zi-zj); /* H */ykb2*ipq+j=-Ui*Uj*(Gi1j1*cos(zi-zj)+Bi1j1*si n( zi-zj); /* J */for(i=pq+1;i<=n-1;i+)/*形成pv节点子阵*/for(j=1;j< n;j+

22、)int i1=jdi. num;int j1=jdj. num;double Ui=jdi.U;double Uj=jdj.U;double zi=jdi.zkj;double zj=jdj.zkj;if(j<=pq)/*求 j<=pq 时的 H、N */ykbpq+i2*j-1=Ui*Uj*(Gi1j1*si n( zi-zj)-Bi1j1*cos(zi-zj);/* H */ykbpq+i2*j=Ui*Uj*(Gi1j1*cos(zi-zj)+Bi1j1*si n( zi-zj);/* N */else/*求 j>pq 时的 H*/if(i!=j)/* 求 i!=j 时

23、的 H*/ykbpq+ipq+j=Ui*Uj*(Gi1j1*si n( zi-zj)-Bi1j1*cos(zi-zj); /* H */else/* 求 i=j 时的 H*/aai=0;for(k=1;k<=n ;k+)int k1=jdk. num;H=0;H=jdk.U*(Gi1k1*si n(jdi.zkj-jdk.zkj)-Bi1k1*cos(jdi.zkj-jdk.zkj); aai=aai+H;ykbpq+ipq+j=-Ui*(aai-Ui*(Gi1i1*s in (jd i.zkj-jdi.zkj)-Bi1i1*cos(jdi.zkj-jdi.zkj); /*H*/* 输出

24、雅克比矩阵*/fprintf(fp2,"nn雅克比矩阵为:");for(i=1;i<=(2*pq+pv);i+)fprin tf(fp2,"n");for(j=1;j<=2*pq+pv;j+)fprintf(fp2,"%f",ykbij);voidSolve_Equations()/*求解修正方程组(LU分解法)*/double lNINI=0; / 定义 L 矩阵double uNlNl=0; / 定义 U 矩阵double yNl=0; / 定义数组 Ydouble xNl=0; / 定义数组 Xdouble aNIN

25、I=0; / 定义系数矩阵 double bNI=0; II 定义右端项 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+)aij=ykbi+1j+1;for(i=0; i<n; i+)bi=Di+1;for(i=0; i<n; i+)/*初始化矩阵I*/ for(j=0; j< n; j+)if(i=j) lij = 1；for(i=0;i <n ;i+)*/u0i=(float)(a0i);/* 第的首行进行计算*/for(k=0;k< n-1;k

26、+)步进行LU分解*/ for(i=k+1;i <n ;i+)/*开始LU分解一步：对矩阵U/*第二步：逐/*对L的第k列进行计算*/ for(s=0,sum=0;s <n; s+) if(s!=k)sum+=lis*usk;lik=(float)(aik-sum)/ukk);/*对U的第for(j=k+1;j <n ;j+)k+1行进行计算*/ for(s=0,sum=0;s <n; s+)if(s!=k+1)sum+=lk+1s*usj;uk+1j=(float)(ak+1j-sum);yO=bO组Y*/*回代法计算数for(i=1;i <n ;i+) for

27、(j=0,sum=0;j<i;j+) sum+=yj*lij;yi=(float)(bi-sum);/*回代法计/*选出最大的修xn-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+=xj*uij;x【i=(float)(y【i-sum)/ui【i);for(i=1; i <=n; i+) di=xi-1; max=fabs(d1);正量的值*/for(i=1;i<=n ;i+) if(fabs(di)>max) max=fabs(di);void N

28、iudu n_Lafux un()int z=1;fprin tf(fp2,"n 极坐标形式的牛顿-拉夫逊迭代法结果：n");do max=1;if(z<Nl)&&(max>=eps)fprin tf(fp2,"n迭代次数：%dn",z);/*开始迭代计算*/Calculate_U nbala nced_Para();Form_Jacobi_Matric();Solve_Equati on s();for(i=1;i<=pq;i+)jdi.zkj+=d2*i-1;jdi.U+=d2*i*jdi.U;for(i=pq+1;

29、i<=n-1;i+)jdi.zkj+=dpq+i;fprintf(fp2,"nn输出 d 8 ,dU:");/*输出修正量的值*/for(c=1;c<=n; c+)for(a=1;a<=n; a+)if(jda. num=c)break;fprin tf(fp2,"n");if(jda.ty=1)fprintf(fp2,"节点为 %2dd 8 =%8.5fdU=%8.5f",c,d2*a-1,d2*a);if(jda.ty=2)fprintf(fp2,"节点为 %2dd 8 =%8.5f",c,d

30、pq+a);fprin tf(fp2,"nn输出迭代过程中的电压值:");for(c=1;c<=n; C+) for(a=1;a<=n; a+)if(jda. num=c) break;fprin tf(fp2,"nU%d=%f",c,jda.U);fprintf(fp2," / %f",jda.zkj);fprin tf(fp2,"nn");z+; while(z<NI)&&(max>=eps);/* 判断是否达到精度要求*/void Powerflow_Result()i

31、nt n 1=jd n. num;fprin tf(fp2,"nn=nn");fprin tf(fp2,"*潮流计算结果*");fprin tf(fp2,"nn=nn");fprintf(fp2,"n各节点电压值：")；/*输出各节点的电压值*/for(c=1;c<=n; C+) for(a=1;a<=n ;a+)if(jda. num=c)break;fprin tf(fp2,"nU%d=%f",c,jda.U);fprintf(fp2," / %f",jda.

32、zkj);rr=tt=0;/* 计算节点的注入功率*/for(i=1;i<=n ;i+)int i1=jdi. num;ji(jdi.U,-jdi.zkj,Gn1i1,-B n1i1);rr+=c1;tt+=d1;ji(jd【n.U,jd【n.zkj,rr,tt);fprintf(fp2,"nn各节点注入功率:n");for(i=1;i<=pq;i+)int i1=jdi. num;fprintf(fp2,"PQ 节点： 节点 %dS%d=%f",i1,i1,jdi.p); /*PQ 节点注入功if(jdi.q>=0)fprin tf(

33、fp2,"+j%fn",jdi.q);elsefprin tf(fp2,"-j%fn",-jdi.q);for(i=pq+1;i<=n-1;i+)int i1=jdi. num;fprintf(fp2,"PV 节点： 节点 %dS%d=%f",i1,i1,jdi.p); /*PV 节点注入功率*/if(jdi.q>=0)fprin tf(fp2,"+j%fn",jdi.q);elsefprin tf(fp2,"-j%fn",-jdi.q);fprintf(fp2,"平衡节点

34、： 节点 %d",jd n.nu m,jd n. num);/* 平衡节点注入功率*/fprin tf(fp2,"S%d=%f", n1,c1);if(d1>=0)fprin tf(fp2,"+j%f",d1);elsefprin tf(fp2,"-j%f",-d1);fprintf(fp2,"nn线路功率:n");rr=tt=0;for(i=1;i<=m;i+)int i1=zli.p1;/*计算线路功率*/in t j1=zli.p2;aai=bbi=P=Q=O;for(a=1;a<

35、=n; a+)if(jda. num=i1)break;for(b=1;b<=n; b+)if(jdb. num=j1)break;mozhi(zli.r,zli.x);if(mo=0)con ti nue;sha ng(1,0,zli.r,zli.x);ji(jda.U/zli.kx/zli.kx,-jda.zkj,c2,-d2);/*考虑非标准变比*/P+=c1;Q+=d1;ji(jdb.U/zli.kx,-jdb.zkj,c2,-d2);P-=c1;Q-=d1;ji(jda.U,jda.zkj,P,Q);fprintf(fp2,"线路 %d: %d-%d 的功率为：f",i,i1,j1,c1);if(d1>=0)fprin tf(fp2,"+j%fn",d1);elsefprin tf(fp2,"-j%fn",-d1)