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

/*利用牛顿-拉夫逊迭代法(极坐标形式)，计算复杂电力系统潮流，具有收敛性好，收敛速度快等优点。所有参数应归算至标幺值下。/*可计算最大节点数为100，可计算PQ,PV,平衡节点*/*可计算非标准变比和平行支路*/#include#include#include#define M 100 /*最大矩阵阶数*/#defineNl 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 aaM,bbM,ccM,ddM,max, rr,tt; /*中间变量*/ double mo,c1,d1,c2,d2; /*复数运算函数的返回值*/ double GMM,BMM,YMM; /*节点导纳矩阵中的实部、虚部及其模方值*/double ykbMM,DM,dM,dUM; /*雅克比矩阵、不平衡量矩阵*/struct jd /*节点结构体*/ int num,ty; /* num为节点号，ty为节点类型*/ double p,q,S,U,zkj,dp,dq,du,dj; /*节点有功、无功功率，功率模值，电压模值,阻抗角 牛顿-拉夫逊中功率不平衡量、电压修正量*/ jdM;struct zl /*支路结构体*/ int 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+) /*输入节点编号、类型、输入功率和电压初值*/ 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,%lfn,&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); fprintf(fp2, 节点数:%d 支路数:%d PQ节点数:%d PV节点数:%d 精度:%fn, n,m,pq,pv,eps); fprintf(fp2, -n); for(i=1;i=pq;i+) fprintf(fp2, PQ节点: 节点%d P%d=%f Q%d=%fn, jdi.num,jdi.num,jdi.p,jdi.num,jdi.q); for(i=pq+1;i=pq+pv;i+) fprintf(fp2, PV节点: 节点%d P%d=%f U%d=%f 初值Q%d=%fn, jdi.num,jdi.num,jdi.p,jdi.num,jdi.U,jdi.num,jdi.q); fprintf(fp2, 平衡节点: 节点%d e%d=%f f%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=%f X=%f n, 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 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+) 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; else mozhi(zlj.r,zlj.x); if(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; else mozhi(zlj.r,zlj.x); if(mo=0) continue; shang(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; mozhi(zlk.r,zlk.x); if(mo=0) continue; shang(1,0,zlk.r,zlk.x); Gij-=c2; Bij-=d2; Gji=Gij; Bji=Bij; else i=zlk.p1; j=zlk.p2; mozhi(zlk.r,zlk.x); if(mo=0) continue; shang(1,0,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); else Bik=-Bik; fprintf(fp2,-j); fprintf(fp2,%f ,Bik); Bik=-Bik; fprintf(fp2,n -n); /*-牛顿拉夫逊-*/ 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-jdi.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.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,n dp%d=%f,i,D2*t-1); fprintf(fp2,n dq%d=%f,i,D2*t); for(i=pq+1;i=n-1;i+)t=jdi.num;fprintf(fp2,n dp%d=%f,i,Dpq+t); void Form_Jacobi_Matric() /*形成雅克比矩阵*/ 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*(Gi1j1*cos(zi-zj)+Bi1j1*sin(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=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+) /*形成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*(Gi1k1*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() /* 求解修正方程组 (LU分解法)*/ double lNlNl=0; /定义L矩阵 double uNlNl=0; /定义U矩阵 double yNl=0; /定义数组Y double xNl=0; /定义数组X double aNlNl=0; /定义系数矩阵 double bNl=0; /定义右端项 double sum=0; int i,j,k,s; int n;n=2*pq+pv; for(i=0; in; i+) for(j=0; jn; j+) aij=ykbi+1j+1; for(i=0; in; i+) bi=Di+1; for(i=0; in; i+) /*初始化矩阵l*/ for(j=0; jn; j+) if(i=j) lij = 1; for(i=0;in;i+) /*开始LU分解*/ u0i=(float)(a0i); /*第一步：对矩阵U的首行进行计算*/for(k=0;kn-1;k+) /*第二步：逐步进行LU分解*/ for(i=k+1;in;i+) /*对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)(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 -极坐标形式的牛顿-拉夫逊迭代法结果:-n); do max=1; if(z=eps) fprintf(fp2,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) fprintf(fp2, 节点为 %2d d=%8.5f dU=%8.5f,c,d2*a-1,d2*a);if(jda.ty=2) fprintf(fp2, 节点为 %2d d=%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,n U%d=%f,c,jda.U);fprintf(fp2,%f,jda.zkj);fprintf(fp2,nn -); z+; while(z=eps); /*判断是否达到精度要求*/ void Powerflow_Result()int 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,n U%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+=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); else fprintf(fp2,-j%fn,-jdi.q); for(i=pq+1;i=0) fprintf(fp2,+j%fn,jdi.q); else fprintf(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);else fprintf(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); else fprintf(fp2,-j%fn,-d1);aai+=c