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

1、/* 利用牛顿-拉夫逊迭代法(极坐标形式)，计算复杂电力系统潮流，具有收敛性好，收敛速度快等优点。所有参数应归算至标幺值下。/* 可计算最大节点数为100 ，可计算PQ,PV, 平衡节点/* 可计算非标准变比和平行支路*/#include#include#include#define M 100#defineNl 100int i,j,k,a,b,c;int t,l;double P,Q,H,J;int n,m,pq,pv;double eps;double aaM,bbM,ccM,ddM,max, rr,tt;double mo,c1,d1,c2,d2;double GMM,BMM,YMM;

2、double ykbMM,DM,dM,dUM;struct jdint num,ty;double p,q,S,U,zkj,dp,dq,du,dj;*/* 最大矩阵阶数*/* 迭代次数 */* 循环控制变量*/* 中间变量 */* 节点数 */* 支路数 */*PQ 节点数 */*PV 节点数 */* 迭代精度 */* 中间变量 */* 复数运算函数的返回值*/* 节点导纳矩阵中的实部、虚部及其模方值*/* 雅克比矩阵、不平衡量矩阵*/* 节点结构体*/* num 为节点号， ty 为节点类型*/* 节点有功、无功功率，功率模值，电压模值,阻抗角牛顿 - 拉夫逊中功率不平衡量、电压修正量 */

3、 jdM;struct zlint numb;/* 支路结构体 */*numb 为支路号*/int p1,p2;/* 支路的两个节点*/double kx;/* 非标准变比*/double r,x;/* 支路的电阻与电抗*/ zlM; FILE *fp1,*fp2;void data()/*读取数据*/int h,number;fp1=fopen(input.txt,r);fscanf(fp1,%d,%d,%d,%d,%lfn,&n,&m,&pq,&pv,&eps);/* 输入节点数,支路数,PQ节点数,PV节点数和迭代精度*/for(i=1;i=n;i+)/* 输入节点编号、类型、输入功率和

4、电压初值*/fscanf(fp1,%d,%d,&number,&h);if(h=1)/* 类型 h=1 是 PQ 节点 */fscanf(fp1,%lf,%lf,%lf,%lfn,&jdi.p,&jdi.q,&jdi.U,&jdi.zkj);jdi.num=number;jdi.ty=h;if(h=2)/* 类型 h=2 是 pv 节点 */fscanf(fp1,%lf,%lf,%lfn,&jdi.p,&jdi.U,&jdi.zkj);jdi.num=number;jdi.ty=h;jdi.q=-1.567;if(h=3)/* 类型 h=3 是平衡节点 */fscanf(fp1,%lf,%lf

5、n,&jdi.U,&jdi.zkj);jdi.num=number;jdi.ty=h;for(i=1;i=m;i+) /* 输入支路阻抗 */ fscanf(fp1,%d,%lf,%d,%d,%lf,%lfn,&zli.numb,&zli.kx,&zli.p1,&zli.p2,&zli.r,&zli.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);fpr

6、intf(fp2,节点数:%d支路数:%dPQ 节点数 :%dPV节点数:%d精度 :%fn,n,m,pq,pv,eps);fprintf(fp2,-n);for(i=1;i=pq;i+)fprintf(fp2,PQ 节点 :节点 %dP%d=%fQ%d=%fn,jdi.num,jdi.num,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,jdi.num,jdi.num,jdi.p,jdi.num,jdi.U,jdi.num,jdi.q);fprintf(fp2,平

7、衡节点:节点 %de%d=%ff%d=%fn,jdn.num,jdn.num,jdn.U,jdn.num,jdn.zkj);fprintf(fp2,-n);for(i=1;i=m;i+)fprintf(fp2,支路 %d:相关节点:%d,%d非标准变比:kx=%f R=%fX=%fn,i,zli.p1,zli.p2,zli.kx,zli.r,zli.x);fprintf(fp2,=n);/*-复数运算函数double mozhi(double a0,double b0)-*/* 复数求模值函数*/mo=sqrt(a0*a0+b0*b0);return mo;double ji(double a

8、1,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); retu

9、rn c2;return d2;/*-计算节点导纳矩阵-*/void Form_Y()for(i=1;i=n;i+)/* 节点导纳矩阵元素附初值*/for(j=1;j=n;j+)Gij=Bij=0;for(i=1;i=n;i+)/* 节点导纳矩阵的主对角线上的元素,非对地导纳加入相应的主对角线元素中(考虑非标准变比）*/for(j=1;j=m;j+)if(zlj.p1=i)if(zlj.kx=1)mozhi(zlj.r,zlj.x);if(mo=0)continue;shang(1,0,zlj.r,zlj.x);Gii+=c2;Bii+=d2;elsemozhi(zlj.r,zlj.x);if

10、(mo=0)continue;shang(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=1)mozhi(zlj.r,zlj.x);if(mo=0)continue;shang(1,0,zlj.r,zlj.x);Gii+=c2;Bii+=d2;elsemozhi(zlj.r,zlj.x);if(mo=0)continue;shang(1,0,zlj.r,zlj.x);Gii+=

11、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;mozhi(zlk.r,zlk.x);if(mo=0)continue;shang(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;shang(1,0,

12、zlk.r,zlk.x);Gij-=c2/zlk.kx;Bij-=d2/zlk.kx;Gji=Gij;Bji=Bij;/*-输出节点导纳矩阵-*/fprintf(fp2,nn*潮流计算过程*n);fprintf(fp2,n=n);fprintf(fp2,n节点导纳矩阵为:);for(i=1;i=n;i+)fprintf(fp2,n);for(k=1;k=0)fprintf(fp2,+j);fprintf(fp2,%f,Bik);elseBik=-Bik;fprintf(fp2,-j);fprintf(fp2,%fBik=-Bik;,Bik);fprintf(fp2,n-n);/*-牛顿拉夫逊-

13、*/void Calculate_Unbalanced_Para()for(i=1;i=n;i+)if(jdi.ty=1)/* 计算 PQ 节点不平衡量*/t=jdi.num;cct=ddt=0;for(j=1;j=n;j+)for(a=1;a=n;a+)if(jda.num=j)break;P=Q=0;P=jda.U*(Gtj*cos(jdi.zkj-jda.zkj)+Btj*sin(jdi.zkj-jda.zkj);Q=jda.U*(Gtj*sin(jdi.zkj-jda.zkj)-Btj*cos(jdi.zkj-jda.zkj);cct+=P;ddt+=Q;jdi.dp=jdi.p-jd

14、i.U*cct;jdi.dq=jdi.q-jdi.U*ddt;if(jdi.ty=2)/* 计算 PV 节点不平衡量*/t=jdi.num;cct=ddt=0;for(j=1;j=n;j+)for(a=1;a=n;a+)if(jda.num=j)break;P=Q=0;P=jda.U*(Gtj*cos(jdi.zkj-jda.zkj)+Btj*sin(jdi.zkj-jda.zkj);Q=jda.U*(Gtj*sin(jdi.zkj-jda.zkj)-Btj*cos(jdi.zkj-jda.zkj);cct+=P;ddt+=Q;jdi.dp=jdi.p-jdi.U*cct;jdi.q=jdi.

15、U*ddt;for(i=1;i=pq;i+)/* 形成不平衡量矩阵DM*/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;fprintf(fp2,nfprintf(fp2,ndp%d=%f,i,D2*t-1);dq%d=%f,i,D2*t);for(i=pq+1;i=n-1;i+)t=jdi.num;fprintf(fp2,ndp%d=%f,i,Dpq+t);void Form_Jacobi_Matr

16、ic()/* 形成雅克比矩阵*/for(i=1;i=pq;i+)/* 形成pq 节点子阵*/for(j=1;jn;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(j=pq)/* 求j=pq时的H、N、J、L */if(i!=j)/* 求i!=j时的H、N、J、L*/ykb2*i-12*j-1=Ui*Uj*(Gi1j1*sin(zi-zj)-Bi1j1*cos(zi-zj);/*H*/ykb2*i-12*j=Ui*Uj*(Gi1

17、j1*cos(zi-zj)+Bi1j1*sin(zi-zj);/*N*/ykb2*i2*j-1=-ykb2*i-12*j;ykb2*i2*j=ykb2*i-12*j-1;/*J*/*L*/else/* 求 i=j 时的 H 、 N 、 J、 L*/aai=0;bbi=0;for(k=1;kpq 时的 H 、 J */ykb2*i-1pq+j=Ui*Uj*(Gi1j1*sin(zi-zj)-Bi1j1*cos(zi-zj); /* H */ykb2*ipq+j=-Ui*Uj*(Gi1j1*cos(zi-zj)+Bi1j1*sin(zi-zj); /* J */for(i=pq+1;i=n-1;i

18、+)/* 形成pv节点子阵*/for(j=1;jn;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(j=pq)/* 求 jpq时的H*/if(i!=j)/* 求i!=j时的H*/ykbpq+ipq+j=Ui*Uj*(Gi1j1*sin(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*(Gi

19、1k1*sin(jdi.zkj-jdk.zkj)-Bi1k1*cos(jdi.zkj-jdk.zkj); aai=aai+H;ykbpq+ipq+j=-Ui*(aai-Ui*(Gi1i1*sin(jdi.zkj-jdi.zkj)-Bi1i1*cos(jdi.zkj-jdi.zkj);/*H*/*-输出雅克比矩阵-*/fprintf(fp2,nn雅克比矩阵为: );for(i=1;i=(2*pq+pv);i+)fprintf(fp2,n);for(j=1;j=2*pq+pv;j+)fprintf(fp2,%f,ykbij);void Solve_Equations()double lNlNl=0

20、; /定义 L 矩阵double uNlNl=0; /定义 U 矩阵double yNl=0; /定义数组Ydouble xNl=0; /定义数组Xdouble aNlNl=0; /定义系数矩阵double bNl=0; /定义右端项double sum=0;int i,j,k,s;int n;n=2*pq+pv;for(i=0; in; i+)/*求解修正方程组(LU分解法)*/for(j=0; jn; j+)aij=ykbi+1j+1;for(i=0; in; i+)bi=Di+1;for(i=0; in; i+)for(j=0; jn; j+)if(i=j) lij = 1;/* 初始化

21、矩阵l*/for(i=0;in;i+)u0i=(float)(a0i);for(k=0;kn-1;k+)for(i=k+1;in;i+)/* 开始 LU 分解 */* 第一步：对矩阵U 的首行进行计算/* 第二步：逐步进行LU 分解 */* 对 L 的第 k 列进行计算 */*/ for(s=0,sum=0;sn;s+) if(s!=k)sum+=lis*usk;lik=(float)(aik-sum)/ukk);for(j=k+1;jn;j+)/* 对U 的第k+1行进行计算*/for(s=0,sum=0;sn;s+)if(s!=k+1)sum+=lk+1s*usj;uk+1j=(float

22、)(ak+1j-sum);y0=b0 ;/* 回代法计算数组Y*/for(i=1;in;i+) for(j=0,sum=0;j=0;i-)for(j=n-1,sum=0;ji;j-)sum+=xj*uij;xi=(float)(yi-sum)/uii);for(i=1; i=n; i+)di=xi-1;max=fabs(d1);/* 选出最大的修正量的值*/for(i=1;imax)max=fabs(di);void Niudun_Lafuxun()int z=1;fprintf(fp2,n- 极坐标形式的牛顿do-拉夫逊迭代法结果:-n);max=1; if(z=eps)fprintf(fp

23、2,n迭代次数 :%dn,z);/* 开始迭代计算*/Calculate_Unbalanced_Para();Form_Jacobi_Matric();Solve_Equations();for(i=1;i=pq;i+)jdi.zkj+=d2*i-1;jdi.U+=d2*i*jdi.U;for(i=pq+1;i=n-1;i+)jdi.zkj+=dpq+i;fprintf(fp2,nn 输出 d ,dU: ); /* 输出修正量的值 */ for(c=1;c=n;c+)for(a=1;a=n;a+)if(jda.num=c)break;fprintf(fp2,n);if(jda.ty=1)fpr

24、intf(fp2,节点为%2dd =%8.5fdU=%8.5f,c,d2*a-1,d2*a);if(jda.ty=2)fprintf(fp2,节点为%2dd =%8.5f,c,dpq+a);fprintf(fp2,nn输出迭代过程中的电压值: );for(c=1;c=n;c+)for(a=1;a=n;a+)if(jda.num=c)break;fprintf(fp2,nU%d=%f,c,jda.U);fprintf(fp2, %f,jda.zkj);fprintf(fp2,nn-);z+; while(z=eps);/* 判断是否达到精度要求*/void Powerflow_Result()i

25、nt n1=jdn.num;fprintf(fp2,nn=nn);fprintf(fp2,*潮流计算结果*);fprintf(fp2,nn=nn);fprintf(fp2,n各节点电压值: );/* 输出各节点的电压值*/for(c=1;c=n;c+)for(a=1;a=n;a+)if(jda.num=c)break;fprintf(fp2,nU%d=%f,c,jda.U);fprintf(fp2, %f,jda.zkj);rr=tt=0;/* 计算节点的注入功率*/for(i=1;i=n;i+)int i1=jdi.num;ji(jdi.U,-jdi.zkj,Gn1i1,-Bn1i1);rr

26、+=c1;tt+=d1;ji(jdn.U,jdn.zkj,rr,tt);fprintf(fp2,nn各节点注入功率:n);for(i=1;i=0)fprintf(fp2,+j%fn,jdi.q);elsefprintf(fp2,-j%fn,-jdi.q);for(i=pq+1;i=0)fprintf(fp2,+j%fn,jdi.q);elsefprintf(fp2,-j%fn,-jdi.q);fprintf(fp2,平衡节点 :节点 %d,jdn.num,jdn.num);/* 平衡节点注入功率*/fprintf(fp2,S%d=%f,n1,c1);if(d1=0)fprintf(fp2,+j%f,d1);elsefprintf(fp2,-j%f,-d1);fprintf(fp2,nn线路功率 :n);rr=tt=0;for(i=1;i=m;i+)int i1=zli.p1;/* 计算线路功率*/int j1=zli.p2;aai=bbi=P=Q=0;for(a=1;a=n;a+)if(jda.num=i1)break;for(b=1;b=0)fprintf(fp2,+j%fn,d1);elsefprintf(fp2,-j%fn,-d1);aai+=c1;bbi+=d1;P=Q=0;ji(jdb.U,-jd