电力系统潮流计算C程序

1、电力系统潮流计算注：这是一个基于N- R法的潮流计算通用程序，仅提供了子程序，需要做些处 理才能成为一个可运行的计算程序！此程序非我原创，仅与大家共享！ ！/*这里提供的是电力系统潮流计算机解法的五个子程序,采用的方法是*Newt on _Raphs on法.*程序中所用的变量说明如下：*N:网络节点总数.M:网络的PQ节点数.*L:网络的支路总数.N0:雅可比矩阵的行数.*N1:N0+1K:打印开关.K=1,则打印;否则,不打印.* K1:子程序PLSC中判断输入电压的形式.K仁1,则为极 座标形式.否则* 为直角坐标形式.* D:有功及无功功率误差的最大值.'* G(l,J):Yb

2、us 的电导元素(实 部).*B(I,J):Ybus 的电纳元素(虚部).S1(I):第I节点的起始节点号E1(I):*G1(I):第I支路的串联电导.B1(I):第I支路的串联电纳.*纳.C1(I):第I支路的pie型对称接地电*纳.C(I,J):第I节点J支路不对称接地电*CO(I):第I节点的接地电第I节点的终止节点号.* P(I):第I节点的注入有功功率Q(I):第I节点的注入无功功率.* P0(I):第I节点有功功率误差.Q0(I):第I节点无功功率误差.*差)V0(I):第I节点(PV节点)的电压误差(平方误*i差).*V(I):第I节点的电压误差幅值*E(I):第I节点的电压的实

3、部F(I):第I节点的电压的虚部.*素JM(l,J):Jacoby矩阵的第1行J列元*糸.*A(I,J):修正方程的增广矩阵，三角化矩阵的第1行J列元素,运算结*束后A矩阵的最后一列存放修正的解.*P1(I):第I支路由S1(I)节点注入的有功功率.*Q1(I):第I支路由S1(I)节点注入的无功功率.*P2(I):第I支路由E1(I)节点注入的有功功率.*Q2(I):第I支路由E1(I)节点注入的无功功率.*P3(I):第I支路的有功功率损耗.*Q3(I):第I支路的无功功率损耗.*ANGLE(I):第I节点电压的角度.*/#in elude <math.h>#i nclude

4、<stdio.h>#defi ne f1(i) (i-1)/*把习惯的一阶矩阵的下标转化为C语言数组下标*/#define f2(i,j,n) (i-1)*(n)+j-1)/*把习惯的二阶矩阵的下标转化为C语言数组下标*/*本子程序根据所给的支路导纳及有关信息，形成结点*导纳矩阵,如打印参数K=1,则输出电导矩阵G和电纳矩B*/void ybus(i nt n,i nt l,i nt m,float *g,float *b,float *g1,float *b1,float *c1,float *c,float *co,i nt k,i nt *s1,i nt *e1) exter

5、n FILE *file4;FILE *fp;int i,j,io,iO;int pos1,pos2;int st,e n;if(file4=NULL)fp=stdout;elsefp=file4; /* 输出到文件*/ /*初始化矩阵G,B */for(i=1;i<=n ;i+)for(j=1;j<=n ;j+)pos2=f2(i,j, n); gpos2=0;bpos2=0;/*计算支路导纳*/for(i=1;i<=l;i+)/* 计算对角元*/pos1=f1(i);st=s1pos1;e n=e1pos1; pos2=f2(st,st, n);gpos2+=g1pos1

6、; bpos2+=b1pos1+c1pos1; pos2=f2(e n,en,n);gpos2+=g1pos1;bpos2+=b1pos1+c1pos1;/*计算非对角元*/ pos2=f2(st,e n,n);gpos2-=g1pos1;bpos2-=b1pos1;gf2(e n,st, n) =gf2(st,e n,n); bf2(e n, st, n) =bf2(st,e n,n);/*计算接地支路导纳*/for(i=1;i<=n ;i+)/* 对称部分 */ bf2(i,i, n)+=cof1(i);/*非对称部分 */for(j=1;j<=l;j+)bf2(i,i, n)

7、+=cf2(i,j,l);if(k!=1)return; /* 如果K不为1,则返回;否则,打印导纳矩阵*/fprin tf(fp,"nBUS ADMITTANCE MATRIXY(BUS):");fprin tf(fp,"n*ARRAYG*");for(io=1;io <=n ;io+=5)i0=(io+4)> n?n :(io+4); fprin tf(fp,"n");for(j=io;j<=i0;j+)fprin tf(fp,"%13d",j);for(i=1;i<=n ;i+)fpr

8、in tf(fp,"n%2d",i);for(j=io;j<=i0;j+)fprin tf(fp,"%13.6f",gf2(i,j, n);fprin tf(fp,"n");fprin tf(fp,"n*ARRAYBfor(io=1;io <=n ;io+=5) i0=(io+4)> n?n :(io+4); fprin tf(fp,"n");for(j=io;jv=i0;j+)fprin tf(fp,"%13d",j);for(i=1;i<=n ;i+)fpr

9、in tf(fp,"n%2d",i);for(j=io;j<=i0;j+)fprin tf(fp,"%13.6f",bf2(i,j, n); fprin tf(fp,"n");fprin tf(fp,"n/*本子程序根据所给的功率及电压等数据*求出功率及电压误差量，并返回最大有功功率*以用于与给定误差比较.如打印参数K=1,则输*出P0,Q0（对PC结点）,V0（对PV结点）.*/ void dpqc(float *p,float *q,float *p0,float *q0,float *v,float *v0,i

10、ntm,*g,floatintn float*e,float *f,int k,float*b,float *dd)extern FILE *file4;FILE *fp;int i,j,l;int pos1,pos2;float a1,a2,d1,d;if(file4=NULL)fp=stdout; /* 输出到屏幕*/elsefp=file4;/* 输出到文件 */l=n-1;if(k=1)fprin tf(fp,"nCHANGE OFP0,V*2,P0(l),Q0(l),V0(l)");fprin tf(fp,"nIP0(I)Q0(I)");for

11、(i=1;i<=l;i+)a1=0;a2=0;pos1=f1(i);for(j=1;j<=n ;j+)/* a1,a2对应课本p185式(4-86)中括号内的式子*/pos2=f2(i,j, n);a1+=gpos2*ef1(j)-bpos2*ff1(j); a2+=gpos2*ff1(j)+bpos2*ef1(j);/* 计算式(4-86)(4-87)中的 deltaPi */ p0pos1=ppos1-epos1*a1-fpos1*a2; if(i <= m) /* 计算PQ结点中的deltaQi */q0pos1=qpos1-fpos1*a1+epos1*a2;else

12、 /* 计算PV结点中的deltaVi平方 */v0pos1=vpos1*vpos1-epos1*epos1-fpos1*fpos1; /*输出结果*/if(k=1)if(i<m)%2d%15.6e %15.6e",i,p0pos1,q0pos1);fprin tf(fp,"nelse if(i=m)fprin tf(fp,"n fprin tf(fp,"n%2d%15.6e %15.6e",i,p0pos1,q0pos1);IP0(l)VO(I)");elsefprintf(fp,"n%2d%15.6e %15.6e

13、",i,p0pos1,v0pos1);/*找到deltaP和deltaQ中的最小者，作为收敛指标，存在dd中*/d=0;for(i=1;i<=l;i+)pos1=f1(i);d1= p0pos1>0 ? pOpos1 : -pOpos1;if(d<d1)d=d1;if(i<=m)d1=q0pos1>0?q0pos1:-q0pos1;if(d<d1)d=d1;(*dd)=d;/* 本子程序根据节点导纳及电压求Jacoby矩阵,用于求* 电压修正量，如打印参数K=1,则输出Jacoby矩阵.* 对应于课本P186式(4-89)(4-90)*/void

14、jmcc(i nt m,i nt n ,i nt n0 float *e,float *f,float *g,float *b,float*jm,i nt k)extern FILE *file4;FILE *fp;int i,j,i1,io,i0,ns;int pos1,pos2;if(file4=NULL)fp=stdout;elsefp=file4;/*初始化矩阵jm */ for(i=1;i<=n 0;i+)for(j=1;j<=n 0;j+) jmf2(i,j, n0)=0;ns=n-1; /*去掉一个平衡结点*/* 计算式(4-89)(4-90) */for(i=1;i

15、<=n s;i+)/* 计算式(4-90) */for(i1=1;i1<=n;i1+)/* posl 是式(4-90)中的 j */ pos1=f1(i1);/* pos2 是式(4-90)中的 ij */ pos2=f2(i,i1, n);if(i<=m) /* i 是 PC结点 */*计算式(4-90)中的Jii等式右侧第一部分*/jmf2(2*i-1,2*i-1, n0)+=gpos2*fpos1+bpos2*epos1 ;/*计算式(4-90)中的Lii等式右侧第一部分*/ jmf2(2*i-1,2*i, n0)+=-gpos2*epos1+bpos2*fpos1;/

16、*计算式(4-90)中的Hii等式右侧第一部分*/ jmf2(2*i,2*i-1, n0)+=-gpos2*epos1+bpos2*fpos1;/*计算式(4-90)中的Nii等式右侧第一部分*/ jmf2(2*i,2*i, n0)+=-gpos2*fpos1-bpos2*epos1;/* pos2 是式(4-90)中的 ii */ pos2=f2(i,i, n);/* pos1 是式(4-90)中的 i */ pos1=f1(i);if(i<=m) /* i 是 PC结点 */* 计算式(4-90)中的Jii */jmf2(2*i-1,2*i-1, n0)+=-gpos2*fpos1+

17、bpos2*epos1 /* 计算式(4-90)中的Lii */jmf2(2*i-1,2*i, n0)+=gpos2*epos1+bpos2*fpoS1;/* 计算式(4-90)中的Hii */ jmf2(2*i,2*i-1, n0)+=-gpos2*epos1-bpos2*fpos1;/* 计算式(4-90)中的Jii */ jmf2(2*i,2*i, n0)+=-gpos2*fpos1+bpos2*epos1;if(i>m) /* PV 结点 */* 计算式(4-90)中的Rii */ jmf2(2*i-1,2*i-1, n0)=-2*epos1;/* 计算式(4-90)中的Sii

18、*/ jmf2(2*i-1,2*i, n0)=-2*fpos1;/* 计算式(4-89) */for(j=1;j<=n s;j+)if(j!=i)/* posl 是式(4-89)中的 i */ pos1=f1(i);/* pos2 是式(4-89)中的 ij */ pos2=f2(i,j, n);/* 计算式(4-89)中的Nij */ jmf2(2*i,2*j, n0)=bpos2*epos1-gpos2*fpos1;/* 计算式(4-89)中的Hij */ jmf2(2*i,2*j-1, n0)=-gpos2*epos1-bpos2*fpos1;if(i<=m) /* i 是

19、PC结点 */* 计算式(4-89)中的 Lij (=-Hij) */ jmf2(2*i-1,2*j, n0)=-jmf2(2*i,2*j-1, n0);/* 计算式(4-89)中的 Jij (=Nij) */ jmf2(2*i-1,2*j-1, n0)=jmf2(2*i,2*j, n0);else /* i 是 PV结点 */* 计算式(4-89)中的 Rij (=0) */ jmf2(2*i-1,2*j-1, n0)=0;/* 计算式(4-89)中的 Sij (=0) */ jmf2(2*i-1,2*j, n0)=0;if(k!=1)return;/*输出Jacoby矩阵*/fprin t

20、f(fp,"n JMATRIX(C)");for(io=1;io <=n 0;io+=5)i仁(io+4)> n0? n0:(io+4);fprin tf(fp,"n");for(j=io;j<=i1;j+)fprin tf(fp,"%10d",j);for(i=1;i<=n 0;i+)fprin tf(fp,"n%2d",i);for(j=io;jv=i1;j+)fprin tf(fp,"%12.6f",jmf2(i,j, nO);fprin tf(fp,"n

21、");/*本子程序用选列主元素的高斯消元法求解组*性方程组求各结点电压修正量,如打印参数K=1,则*输出增广矩阵变换中的上三角及电压修正量如果*无唯一解，则给出信息，并停止程序运行.*/void sevc ( float a, int nO, int k, int n1)extern FILE *file4;FILE *fp;int i,j,l, n2, n3, n4,i0,io,j1,i1;float t0,t,c;if(file4=NULL) fp=stdout;else fp=file4;for(i=1;i<=n 0;i+)l=i;for(j=i;j<=n 0;j+

22、)if( fabs(af2(j,i,n1) > fabs(af2(l,i,n1)l=j; /* 找到这行中的最大元*/if(l!=i) /* 行列交换*/for (j=i;j<=n 1;j+)t=af2(i,j, n1);af2(i,j, n1)=af2(l,j, n1);af2(l,j, n1)=t;if (fabs(af2(i,i, n1)-0)<1e-10) /*对角元近似于0,无解*/prin tf("nNo Solutio nn");exit (1);t0=af2(i,i, n1);for(j=i;j<=n 1;j+)/*除对角元*/af2

23、(i,j, n1)/=t0;if(i=n0)/*最后一行，不用消元*/con ti nue;/* 消元*/j1=i+1;for(i1=j1;i1<=n0;i1+)c=af2(i1,i, n1);for(j=i;j<=n 1;j+)af2(i1,j, n1) -= af2(i,j, n1) *c;if(k=1) /* 输出上三角矩阵*/fprin tf(fp,"nTria nglar An gme ntex Matrix ");for(io=1;io <=n 1;io+=5)i0=(io+4)> n1? n1:(io+4);fprin tf(fp,&q

24、uot;n");fprin tf(fp,"");for(i=io;i<=i0;i+)fprin tf(fp,"%12d",i);for(i=1;i<=n 0;i+)fprin tf(fp,"n");fprin tf(fp,"%2d",i);for(j=io;j<=i0;j+)fprin tf(fp,"%15.6f", af2(i,j, n1);/*回代求方程解*/n2=n1-2;for(i=1;i<=n 2;i+)n 3=n1-i;for(i 1=n 3;i1&

25、lt;=n0;i1+)n4=n0-i;af2( n4, n1, n1) -= af2(i1, n1, n1)*af2( n4,i1, n1); if(k!=1)return;/*输出电压修正值*/fprin tf(fp,"nVoltage correction E(i), F(i):");for(io=1;io <=n 0;io+=4)i1= (io+1)/2;i0=(io+3)/2)>( n0 /2)?( n0/2):(io+3)/2);fprin tf(fp,"n");for(j=i1;j<=i0;j+)fprin tf(fp,&q

26、uot;%16d%16d",j,j);i1 = 2*i0;fprin tf(fp,"n");for(i=io;i<=i1;i+)fprintf(fp,"%15.6f", af2(i,n1,n1);#define Pi 3.1415927/180void plsc(int n,int l,int m,float g,float b,float e,float f, int e1,int s1,float g1,float b1,float c1,float c,float co,float p1,float q1,float p2,floa

27、t q2,float p3,float q3,float p,float q,float v,float angle,int k1)extern FILE *file4;FILE *fp;float t1,t2,st,e n,cm,x,y,z,x1,x2,y1,y2;int i,i1,j,m1, ns,pos1,pos2,km;ns=n-1;if(file4=NULL)fp=stdout;elsefp=file4;fprin tf(fp,"nTHE RESULT ARE:");if(k1=1)for(i=0;i <n ;i+)an glei*=Pi;ei=vi*cos

28、(a nglei); fi=vi*si n(an glei);t1=0.0;t2=0.0;for(i=1;i<=n ;i+)pos 1= f1(i);pos2=f2( n,i, n); t1+=gpos2*epos1-bpos2*fpos1; t2+=gpos2*fpos1+bpos2*epos1;pos1=f1( n);ppos1=t1*epos1;qpos1=-t2*epos1;m1=m+1;for(i1=m1;i1<=n s;i1+)t1=0;t2=0;for(i=1;i<=n ;i+)pos 1= f1(i);pos2=f2(i1,i, n);t1+=gpos2*epos1-bpos2*fpos1;t2+=gpos2*fpos1+bpos2*epos1;pos1=f1(i1);qpos1=fpos1*t1-epos1*t2;for(i=0;i <n; i+)cm=coi;if(cm!=0)qi-=(ei*ei+fi*fi)*cm;fprin tf(fp,"nBUS DATA");fprin tf(fp,&q