牛頓拉夫遜迭代法極坐標(biāo)潮流計(jì)算C語言程序

1、/*最大矩陣階數(shù)*/*最大矩陣階數(shù)*/*迭代次數(shù)*/*循環(huán)控制變量*/*中間變量*/*節(jié)點(diǎn)數(shù)*/*支路數(shù)*/*PQ節(jié)點(diǎn)數(shù)*/*PV節(jié)點(diǎn)數(shù)*/*迭代精度*/*中間變量*/*復(fù)數(shù)運(yùn)算函數(shù)的返回值*/*節(jié)點(diǎn)導(dǎo)納矩陣中的實(shí)部、虛部及其模方值/*雅克比矩陣、不平衡量矩陣*/*節(jié)點(diǎn)結(jié)構(gòu)體*/* num為節(jié)點(diǎn)號，ty為節(jié)點(diǎn)類型*/*節(jié)點(diǎn)有功、無功功率，功率模值，電壓模值，阻牛頓-拉夫遜中功率不平衡量、電壓/*支路結(jié)構(gòu)體*/*numb為支路號*/*支路的兩個節(jié)點(diǎn)*/*非標(biāo)準(zhǔn)變比*/*支路的電阻與電抗*/*讀取數(shù)據(jù)*/*利用牛頓-拉夫遜迭代法（極坐標(biāo)形式），計(jì)算復(fù)雜電力系統(tǒng)潮流，具有收斂性好，收斂速度快等優(yōu)點(diǎn)。

2、所 有 參 數(shù) 應(yīng) 歸 算 至 標(biāo) 幺 值 下/*可計(jì)算最大節(jié)點(diǎn)數(shù)為100，可計(jì)算PQ,PV，平衡節(jié)點(diǎn)*/ /*可計(jì)算非標(biāo)準(zhǔn)變比和平行支路*/#include#include#include#define M 100#define Nl 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;*/double ykbMM,DM,dM,dUM;struct jd( int

3、 num,ty;double p,q,S,U,zkj,dp,dq,du,dj;抗角修正量*/ jdM;struct zl( int numb;int p1,p2;double kx;double r,x; zlM;FILE *fp1,*fp2;void data()(int h,number;fp1=fopen(”,”r”);fscanf(fp1,”d,%d,%d,%d,%lfn”,&n,&m,&pq,&pv,&eps);/*輸入節(jié)點(diǎn)數(shù)，支路數(shù),PQ 節(jié)點(diǎn)數(shù),PV 節(jié)點(diǎn)數(shù)和迭代精度*/ for(i=1;i=n;i+)/*輸入節(jié)點(diǎn)編號、類型、輸入功率和電壓初值*/ ( fscanf(fp1,%d

4、,%d,&number,&h); if(h=1)/*類型 h=1 是 PQ 節(jié)點(diǎn)*/( 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é)點(diǎn)*/( fscanf(fp1,”,lf,%lf,%lfn”,&jdi.p,&jdi.U,&jdi.zkj); jdi.num=number; jdi.ty=h; jdi.q=; if(h=3)/*類型h=3是平衡節(jié)點(diǎn)*/( fscanf(fp1,”,lf,%lfn”,&jdi.U,&jdi.

5、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(”,w”)=NULL) ( printf( can not open file!n); exit(0); fprintf(fp2,”電力系統(tǒng)潮流計(jì)算n );fprintf(fp2,* 原始數(shù)據(jù) *n);fprintf(fp2,= =n);fprintf(fp2,節(jié)

6、點(diǎn)數(shù):d支路數(shù):d PQ節(jié)點(diǎn)數(shù):d PV節(jié)點(diǎn)數(shù):d精度:fn”,n,m,pq,pv,eps); fprintf(fp2, n);for(i=1;i=pq;i+)fprintf(fp2, PQ 節(jié)點(diǎn)： 節(jié)點(diǎn)d P%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é)點(diǎn)： 節(jié)點(diǎn)d P%d=%fU%d=%f 初值 Q%d=%fn,jdi.num,jdi.num,jdi.p,jdi.num,jdi.U,jdi.num,jdi.q);fprintf(fp2,平衡節(jié)點(diǎn)：節(jié)點(diǎn) d e%

7、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é)點(diǎn):d,%d 非標(biāo)準(zhǔn)變比:kx=%f R=%f X=%f n,i,zli.p1,zli.p2,zli.kx,zli.r,zli.x);fprintf(fp2,=n );/*復(fù)數(shù)運(yùn)算函數(shù)*/double mozhi(double a0,double b0)/* 復(fù)數(shù)求模值函數(shù) */(mo=sqrt(a0*a0+b0*b0);return mo;double ji(double a1,doub

8、le b1,double a2,double b2)/*復(fù)數(shù)求積函數(shù) a1 為電壓模值，a2為阻抗角，a3為導(dǎo)納實(shí)部，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);return c2;return d2;/*

9、計(jì)算節(jié)點(diǎn)導(dǎo)納矩陣*/void Form_Y() (for(i=1;i=n;i+)/*節(jié)點(diǎn)導(dǎo)納矩陣元素附初值*/for(j=1;j=n;j+)Gij=Bij=0;for(i=1;i=n;i+)/*節(jié)點(diǎn)導(dǎo)納矩陣的主對角線上的元素，非對地導(dǎo)納加入相應(yīng)的主對角線元素中(考慮非標(biāo)準(zhǔn)變比)*/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;

10、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+

11、c2*(zlj.kx-1)/zlj.kx;Bii+=d2/zlj.kx+d2*(zlj.kx-1)/zlj.kx;for(k=1;k=m;k+)/*節(jié)點(diǎn)導(dǎo)納矩陣非主對角線上(考慮非標(biāo)準(zhǔn)變比)的元素*/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,z

12、lk.x);Gij-=c2/zlk.kx;Bij-=d2/zlk.kx;Gji=Gij;Bji=Bij;/*輸出節(jié)點(diǎn)導(dǎo)納矩陣*/fprintf(fp2,nn* 潮流計(jì)算過程 *n);fprintf(fp2,n=n );fprintf(fp2,n節(jié)點(diǎn)導(dǎo)納矩陣為:)；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);/*牛

13、頓拉夫遜*/void Calculate_Unbalanced_Para() ( for(i=1;i=n;i+) ( if(jdi.ty=1)/*計(jì)算PQ節(jié)點(diǎn)不平衡量*/( 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;

14、 ddt+=Q; jdi.dp=jdi.p-jdi.U*cct; jdi.dq=jdi.q-jdi.U*ddt; if(jdi.ty=2)/*計(jì)算PV節(jié)點(diǎn)不平衡量*/( 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+

15、=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,ndp%d=%f,i,D2*t-1);fprintf(fp2,ndq%d=%f,i,D2*t);for(i=pq+1;i=n-1;i+)(t=jdi.num;fprintf(fp2,n dp%

16、d=%f,i,Dpq+t);void Form_Jacobi_Matric()(for(i=1;i=pq;i+) for(j=1;jn;j+) ( int i1=jdi.num;void Form_Jacobi_Matric()(for(i=1;i=pq;i+) 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)/*形成雅克比矩陣*/*形成pq節(jié)點(diǎn)子陣*/*求 j=pq 時的 H、N、J、L */

17、*求 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-z

18、j)-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é)點(diǎn)子陣 */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+i

19、pq+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(f

20、p2,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; kj+=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 5 ,dU: );/*輸出修正量的值*/for(c=1;c=n;c+)(for(a=1;a=n;a+)(i

21、f(jda.num=c)break;fprintf(fp2,n);if(jda.ty=1)fprintf(fp2, 節(jié)點(diǎn)為 %2d d 5 =%8.5fdU=%8.5f,c,d2*a-1,d2*a);if(jda.ty=2)fprintf(fp2, 節(jié)點(diǎn)為 %2d d 5 =%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, Z%f,jda.zkj);fprintf(fp2

22、,nn );z+; while(z=eps);/* 判斷是否達(dá)到精度要求 */void Powerflow_Result()(int n1=jdn.num;fprintf(fp2,nn=n n);fprintf(fp2,* 潮流計(jì)算結(jié)果 *);fprintf(fp2,nn =nn)；fprintf(fp2,n各節(jié)點(diǎn)電壓值：);/*輸出各節(jié)點(diǎ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, Z%f,jda.zkj);rr=tt=0;/*計(jì)算節(jié)點(diǎn)的注

23、入功率*/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各節(jié)點(diǎn)注入功率: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é)點(diǎn)： 節(jié)點(diǎn) d,jdn.num,

24、jdn.num); /* 平衡節(jié)點(diǎn)注入功率 */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;/* 計(jì)算線路功率 */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,-jdb.zkj,c2,-d2);/*考慮非標(biāo)準(zhǔn)變比*/P+=c1;Q+=d1;ji(jda.U/zli.kx,-