牛頓拉夫遜迭代法極坐標潮流計算C語言程序

1、/* 利用牛頓-拉夫遜迭代法(極坐標形式)，計算復(fù)雜電力系統(tǒng)潮流，具有收斂性好，收斂速度快等優(yōu)點。所有參數(shù)應(yīng)歸算至標幺值下。/* 可計算最大節(jié)點數(shù)為100 ，可計算PQ,PV, 平衡節(jié)點/* 可計算非標準變比和平行支路*/#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;*/* 最大矩陣階數(shù)*/* 迭代次數(shù) */* 循環(huán)控制變量*/* 中間變量 */* 節(jié)點數(shù) */* 支路數(shù) */*PQ 節(jié)點數(shù) */*PV 節(jié)點數(shù) */* 迭代精度 */* 中間變量 */* 復(fù)數(shù)運算函數(shù)的返回值*/* 節(jié)點導(dǎo)納矩陣中的實部、虛部及其模方值*/* 雅克比矩陣、不平衡量矩陣*/* 節(jié)點結(jié)構(gòu)體*/* num 為節(jié)點號， ty 為節(jié)點類型*/* 節(jié)點有功、無功功率，功率模值，電壓模值,阻抗角牛頓 - 拉夫遜中功率不平衡量、電壓修正量 */

3、 jdM;struct zlint numb;/* 支路結(jié)構(gòu)體 */*numb 為支路號*/int p1,p2;/* 支路的兩個節(jié)點*/double kx;/* 非標準變比*/double r,x;/* 支路的電阻與電抗*/ zlM; FILE *fp1,*fp2;void data()/*讀取數(shù)據(jù)*/int h,number;fp1=fopen(input.txt,r);fscanf(fp1,%d,%d,%d,%d,%lfn,&n,&m,&pq,&pv,&eps);/* 輸入節(jié)點數(shù),支路數(shù),PQ節(jié)點數(shù),PV節(jié)點數(shù)和迭代精度*/for(i=1;i=n;i+)/* 輸入節(jié)點編號、類型、輸入功率和

4、電壓初值*/fscanf(fp1,%d,%d,&number,&h);if(h=1)/* 類型 h=1 是 PQ 節(jié)點 */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 節(jié)點 */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 是平衡節(jié)點 */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,電力系統(tǒng)潮流計算n);fprintf(fp2,*原始數(shù)據(jù)*n);fprintf(fp2,= =n);fpr

6、intf(fp2,節(jié)點數(shù):%d支路數(shù):%dPQ 節(jié)點數(shù) :%dPV節(jié)點數(shù):%d精度 :%fn,n,m,pq,pv,eps);fprintf(fp2,-n);for(i=1;i=pq;i+)fprintf(fp2,PQ 節(jié)點 :節(jié)點 %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節(jié)點 :節(jié)點 %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、衡節(jié)點:節(jié)點 %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:相關(guān)節(jié)點:%d,%d非標準變比:kx=%f R=%fX=%fn,i,zli.p1,zli.p2,zli.kx,zli.r,zli.x);fprintf(fp2,=n);/*-復(fù)數(shù)運算函數(shù)double mozhi(double a0,double b0)-*/* 復(fù)數(shù)求模值函數(shù)*/mo=sqrt(a0*a0+b0*b0);return mo;double ji(double a

8、1,double b1,double a2,double b2) /* 復(fù)數(shù)求積函數(shù) a1 為電壓模值， a2 為阻抗角， a3 為導(dǎo)納實部， a4 為導(dǎo)納虛部 */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)/* 復(fù)數(shù)除法求商函數(shù)*/c2=(a3*a4+b3*b4)/(a4*a4+b4*b4); d2=(a4*b3-a3*b4)/(a4*a4+b4*b4); retu

9、rn c2;return d2;/*-計算節(jié)點導(dǎo)納矩陣-*/void Form_Y()for(i=1;i=n;i+)/* 節(jié)點導(dǎo)納矩陣元素附初值*/for(j=1;j=n;j+)Gij=Bij=0;for(i=1;i=n;i+)/* 節(jié)點導(dǎo)納矩陣的主對角線上的元素,非對地導(dǎo)納加入相應(yīng)的主對角線元素中(考慮非標準變比）*/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+)/* 節(jié)點導(dǎo)納矩陣非主對角線上（考慮非標準變比）的元素*/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;/*-輸出節(jié)點導(dǎo)納矩陣-*/fprintf(fp2,nn*潮流計算過程*n);fprintf(fp2,n=n);fprintf(fp2,n節(jié)點導(dǎo)納矩陣為:);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 節(jié)點不平衡量*/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 節(jié)點不平衡量*/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 節(jié)點子陣*/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節(jié)點子陣*/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; /定義數(shù)組Ydouble xNl=0; /定義數(shù)組Xdouble aNlNl=0; /定義系數(shù)矩陣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 ;/* 回代法計算數(shù)組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-拉夫遜迭代法結(jié)果:-n);max=1; if(z=eps)fprintf(fp

23、2,n迭代次數(shù) :%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,節(jié)點為%2dd =%8.5fdU=%8.5f,c,d2*a-1,d2*a);if(jda.ty=2)fprintf(fp2,節(jié)點為%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,*潮流計算結(jié)果*);fprintf(fp2,nn=nn);fprintf(fp2,n各節(jié)點電壓值: );/* 輸出各節(jié)點的電壓值*/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;/* 計算節(jié)點的注入功率*/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各節(jié)點注入功率: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,平衡節(jié)點 :節(jié)點 %d,jdn.num,jdn.num);/* 平衡節(jié)點注入功率*/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