牛頓拉夫遜極坐標形式的潮流上機c源程序

1、/*本程序利用牛頓-拉夫遜迭代法(極坐標形式)，計算復雜電力系統(tǒng)潮流，具有收斂性好，收斂速度快等優(yōu)點。所有參數(shù)應歸算至標幺值下。/*可計算最大節(jié)點數(shù)為100，可計算PQ,PV,平衡節(jié)點*/*可計算非標準變和平行支路*/#include<stdio.h>#include<math.h>#include<stdlib.h>#define M 100 /*最大矩陣階數(shù)*/#defineNl 100 /*迭代次數(shù)*/int i,j,k,a,b,c; /*循環(huán)控制變量*/int t,l; double P,Q,H,J; /*中間變量*/ int n, /*節(jié)點數(shù)*/

2、m, /*支路數(shù)*/ pq, /*PQ節(jié)點數(shù)*/ pv; /*PV節(jié)點數(shù)*/double eps; /*迭代精度*/double aaM,bbM,ccM,ddM,max, rr,tt; /*中間變量*/ double mo,c1,d1,c2,d2; /*復數(shù)運算函數(shù)的返回值*/ double GMM,BMM,YMM; /*節(jié)點導納矩陣中的實部、虛部及其模方值*/double ykbMM,DM,dM,dUM; /*雅克比矩陣、不平衡量矩陣*/struct jd /*節(jié)點結(jié)構(gòu)體*/ int num,ty; /* num為節(jié)點號，ty為節(jié)點類型*/ double p,q,S,U,zkj,dp,dq,

3、du,dj; /*節(jié)點有功、無功功率，功率模值，電壓模值,阻抗角 牛頓-拉夫遜中功率不平衡量、電壓修正量*/ jdM;struct zl /*支路結(jié)構(gòu)體*/ int numb; /*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

4、,%d,%lfn",&n,&m,&pq,&pv,&eps); /*輸入節(jié)點數(shù),支路數(shù),PQ節(jié)點數(shù),PV節(jié)點數(shù)和迭代精度*/ for(i=1;i<=n;i+) /*輸入節(jié)點編號、類型、輸入功率和電壓初值*/ 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=

5、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,%lfn",&jdi.U,&jdi.zkj); jdi.num=number; jdi.ty=h; for(i=1;i<=m;i+) /*輸入支路阻抗*/ fscanf(f

6、p1,"%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 華北電力大學 電力1104 譚程凱 201105030119n"

7、); fprintf(fp2," * 原始數(shù)據(jù) *n"); fprintf(fp2,"=n"); fprintf(fp2," 節(jié)點數(shù):%d 支路數(shù):%d PQ節(jié)點數(shù):%d PV節(jié)點數(shù):%d 精度:%fn", n,m,pq,pv,eps); fprintf(fp2," -n"); for(i=1;i<=pq;i+) fprintf(fp2," PQ節(jié)點: 節(jié)點%d P%d=%f Q%d=%fn", jdi.num,jdi.num,jdi.p,jdi.num,jdi.q); for(i=pq+

8、1;i<=pq+pv;i+) fprintf(fp2," PV節(jié)點: 節(jié)點%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," 平衡節(jié)點: 節(jié)點%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: 相關(guān)

9、節(jié)點:%d,%d 非標準變比:kx=%f R=%f X=%f n", i,zli.p1,zli.p2,zli.kx,zli.r,zli.x); fprintf(fp2," =n"); /*-復數(shù)運算函數(shù)-*/ double mozhi(double a0,double b0) /*復數(shù)求模值函數(shù)*/ mo=sqrt(a0*a0+b0*b0); return mo; double ji(double a1,double b1,double a2,double b2) /*復數(shù)求積函數(shù) a1為電壓模值，a2為阻抗角，a3為導納實部，a4為導納虛部*/ a1=a1*co

10、s(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) /*復數(shù)除法求商函數(shù)*/ c2=(a3*a4+b3*b4)/(a4*a4+b4*b4); d2=(a4*b3-a3*b4)/(a4*a4+b4*b4); return c2; return d2; /*-計算節(jié)點導納矩陣-*/ void Form_Y() for(i=1;i<=n;i+) /*節(jié)點導納矩陣元素附初值*/ for(j=

11、1;j<=n;j+) Gij=Bij=0; for(i=1;i<=n;i+) /*節(jié)點導納矩陣的主對角線上的元素,非對地導納加入相應的主對角線元素中(考慮非標準變比）*/ 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.k

12、x+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+

13、=d2/zlj.kx+d2*(zlj.kx-1)/zlj.kx; for(k=1;k<=m;k+) /*節(jié)點導納矩陣非主對角線上（考慮非標準變比）的元素*/ 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)

14、; Gij-=c2/zlk.kx; Bij-=d2/zlk.kx; Gji=Gij; Bji=Bij; /*-輸出節(jié)點導納矩陣-*/ fprintf(fp2,"nn * 潮流計算過程 *n"); fprintf(fp2,"n =n"); fprintf(fp2,"n 節(jié)點導納矩陣為:"); for(i=1;i<=n;i+) fprintf(fp2,"n "); for(k=1;k<=n;k+) fprintf(fp2,"%f",Gik); if(Bik>=0) fprintf(

15、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節(jié)點不平衡量*/ t=jdi.num; cct=ddt=0; for(j=1;j&l

16、t;=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節(jié)點不平衡量*/ t=jdi.num; cct=ddt=0; for(j=1;j

17、<=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;

18、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_Matr

19、ic() /*形成雅克比矩陣*/ for(i=1;i<=pq;i+) /*形成pq節(jié)點子陣*/ 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(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

20、*/ 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;k<=n;k+) int k1=jdk.num;H=J=0; H=jdk.U*(Gi1k1*sin(jdi.zkj-jdk.zkj)-Bi1k1*cos(jdi.zkj-jdk.zkj);J=jdk.U*(Gi1k1*cos(jdi.zkj-

21、jdk.zkj)+Bi1k1*sin(jdi.zkj-jdk.zkj); aai=aai+H; bbi=bbi+J; ykb2*i-12*j-1=-Ui*(aai-Ui*(Gi1i1*sin(jdi.zkj-jdi.zkj)-Bi1i1*cos(jdi.zkj-jdi.zkj); /*H*/ ykb2*i2*j-1=Ui*(bbi-Ui*(Gi1i1*cos(jdi.zkj-jdi.zkj)+Bi1i1*sin(jdi.zkj-jdi.zkj); /*J*/ ykb2*i-12*j=ykb2*i2*j-1+2*Ui*Ui*Gi1i1; /*N*/ ykb2*i2*j=-ykb2*i-12*j-

22、1-2*Ui*Ui*Bi1i1; /*L*/ else /*求j>pq時的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節(jié)點子陣*/ for(j=1;j<n;j+) int i1=jdi.num; int j1=jdj.num; double Ui=jdi.U; double Uj=jdj.U; doubl

23、e zi=jdi.zkj; double zj=jdj.zkj; if(j<=pq) /*求j<=pq時的H、N */ ykbpq+i2*j-1=Ui*Uj*(Gi1j1*sin(zi-zj)-Bi1j1*cos(zi-zj); /* H */ ykbpq+i2*j=Ui*Uj*(Gi1j1*cos(zi-zj)+Bi1j1*sin(zi-zj); /* N */ else /*求j>pq時的H*/ if(i!=j) /*求i!=j時的H*/ ykbpq+ipq+j=Ui*Uj*(Gi1j1*sin(zi-zj)-Bi1j1*cos(zi-zj); /* H */ else

24、/*求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+)

25、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; /定義數(shù)組Y double xNl=0; /定義數(shù)組X double aNlNl=0; /定義系數(shù)矩陣 double bNl=0; /定義右端項 double sum=0; int i,j,k,s; int n;n

26、=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+) /*初始化矩陣l*/ for(j=0; j<n; j+) if(i=j) lij = 1; for(i=0;i<n;i+) /*開始LU分解*/ u0i=(float)(a0i); /*第一步：對矩陣U的首行進行計算*/for(k=0;k<n-1;k+) /*第二步：逐步進行LU分解*/ for(i=k+1;i<n;i+) /*對L的第k

27、列進行計算*/ for(s=0,sum=0;s<n;s+) if(s!=k) sum+=lis*usk; lik=(float)(aik-sum)/ukk);for(j=k+1;j<n;j+) /*對U的第k+1行進行計算*/ for(s=0,sum=0;s<n;s+) if(s!=k+1) sum+=lk+1s*usj;uk+1j=(float)(ak+1j-sum);y0=b0 ; /*回代法計算數(shù)組Y*/ for(i=1;i<n;i+) for(j=0,sum=0;j<i;j+) sum+=yj*lij; yi=(float)(bi-sum);xn-1=(f

28、loat)(yn-1/un-1n-1); /*回代法計算數(shù)組X*/ for(i=n-2;i>=0;i-) for(j=n-1,sum=0;j>i;j-) sum+=xj*uij;xi=(float)(yi-sum)/uii); 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 Niudun_Lafuxun() int z=1; fprintf(fp2,"n -極坐標形式的牛頓-拉夫遜迭代法結(jié)

29、果:-n"); do max=1; if(z<Nl)&&(max>=eps) fprintf(fp2,"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

30、,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," 節(jié)點為 %2d d=%8.5f dU=%8.5f",c,d2*a-1,d2*a);if(jda.ty=2) fprintf(fp2," 節(jié)點為 %2d d=%8.5f",c,dpq+a);fprintf(fp2,"nn 輸出迭代過程中的電壓值: "); for(

31、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<Nl)&&(max>=eps); /*判斷是否達到精度要求*/ void Powerflow_Result()int n1=jdn.num;fprintf(fp2,"nn =nn");fprintf(

32、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,"n U%d=%f",c,jda.U);fprintf(fp2,"%f",jda.zkj);rr=tt=0; /*計算節(jié)點的注入功率*/for(i=1;i<=n;i+)int i1=jdi.n

33、um;ji(jdi.U,-jdi.zkj,Gn1i1,-Bn1i1);rr+=c1;tt+=d1;ji(jdn.U,jdn.zkj,rr,tt);fprintf(fp2,"nn 各節(jié)點注入功率:n"); for(i=1;i<=pq;i+) int i1=jdi.num; fprintf(fp2," PQ節(jié)點: 節(jié)點%d S%d=%f", i1,i1,jdi.p); /*PQ節(jié)點注入功率*/ if(jdi.q>=0) fprintf(fp2,"+j%fn",jdi.q); else fprintf(fp2,"-j%

34、fn",-jdi.q); for(i=pq+1;i<=n-1;i+) int i1=jdi.num; fprintf(fp2," PV節(jié)點: 節(jié)點%d S%d=%f", i1,i1,jdi.p); /*PV節(jié)點注入功率*/ if(jdi.q>=0) fprintf(fp2,"+j%fn",jdi.q); else fprintf(fp2,"-j%fn",-jdi.q); fprintf(fp2," 平衡節(jié)點: 節(jié)點%d",jdn.num,jdn.num); /*平衡節(jié)點注入功率*/ fprin

35、tf(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<=n;b+) if(jdb.num=j1)break;mozhi(zli.r,zli.x);if(mo=0)continue;shang(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,&q