附录_潮流计算程序.

1、附录附录 牛顿-拉夫逊法潮流计算程序本附录所介绍的牛顿-拉夫逊法潮流计算程序采用极坐标形式，其中所涉及的计算公式与第四章中的式(4-42)-(4-59)完全相同，计算程序框图与图4-6基本一致。程序采用C语言编制。为了便于初学者阅读，节点导纳矩阵和雅可比矩阵都用满阵存储而未采用稀疏技巧；但为了适当照顾使用的方便性，在输入数据中对节点编号次序不作任何要求。下面先介绍输入文件的格式和要求，然后列出程序，最后说明潮流结果的输出。建议读者先从原始数据的输入中了解和熟悉它们在程序中对应的变量、结构体数组及其成员的名称，然后对照图4-6和式(4-42)-(4-59)仔细和耐心地阅读程序，最好能在计算机上亲

2、自实现和进行调试。在调试时，可以用下面给出的对应于例4-3系统的输入数据，以及程序运行中得出的中间结果，逐步与例4-3中所给出的中间结果进行对比，从而查找错误所在并进行改正。一、原始数据的输入程序通过“输入数据.txt”文件输入以下5个数据段。1. 信息(共6个)(1) 总节点数(变量num_node)(2) 线路和并联电容器总数(变量num_line)(3) 变压器支路总数(变量num_tran)(4) 发电机节点总数(变量num_gene)(5) 负荷节点总数(变量num_load)(6) 节点功率不平衡量的容许误差(变量error)2. 线路和并联电容器数据(结构体数组line)：每一线

3、路或并联电容器包括5个数据线路并联电容器成员I侧节点编号所接节点编号iJ侧节点编号同上j型等值电路电阻电容器电阻a型等值电路电抗电容器电抗(负数)b型等值电路一端电纳0.0cI侧节点编号和J侧节点编号可以对换；线路和并联电容器之间的次序可以任意，而且允许多条线路或多个电容器并联。3. 变压器支路数据(结构体数组tran)：每一变压器支路包括5个数据变压器支路成员1侧节点编号i2侧节点编号j电阻a电抗b非标准变比c变压器电阻、电抗和非标准变比与两侧节点编号之间的关系服从图2-29,即电阻和电抗在1侧而非标准变比在2侧；三绕组变压器需按图2-26化成3个变压器支路，其中3侧变压器支路的非标准变比为

4、1。变压器支路之间的次序可以任意。4. 发电机节点数据(结构体数组gene)：每一发电机节点包括5个数据发电机节点成员所在节点编号i节点种类j发出有功功率a发出无功功率b电压c对于PQ节点，节点种类为1，电压可给任意值；对于PV节点，节点种类为-1，发出无功功率可给任意值；对于平衡节点，节点种类为0，发出有功功率和发出无功功率都可给任意值。发电机节点之间的次序可以任意。5. 负荷节点数据(结构体数组load)：每一负荷节点包括3个数据负荷节点成员所在节点编号i吸收有功功率a吸收无功功率b负荷机节点之间的次序可以任意。对于例4-3中的电力系统，输入数据如下：4 4 1 2 2 0.000014

5、3 0.260331 0.495868 0.0258641 4 0.173554 0.330579 0.0172432 2 0.000000 -20.000000 0.000000 3 1 0.130165 0.247934 0.0129321 2 0.000000 0.166667 1.1282054 0 0.0 0.0 1.053 -1 0.2 0.0 1.052 0.5 0.34 0.15 0.1二、计算程序 /* * * * * * 牛顿-拉夫逊法潮流计算程序 * * * * * */# include <stdio.h># include <math.h>#

6、include <stdlib.h># define pnt 1 / 1 - 输出中间结果void read_data(); / 输入始数据void admt_matrix(); / 形成导纳矩阵void form_Jacobian(); / 形成雅可比矩阵, 计算功率误差void solv_Eqn(); / 求解修正方程式void node_flow(); / 输出节点潮流void branch_flow(); / 输出支路潮流double *newSpaceDouble2(int,int);void deleteSpaceDouble2(double *,int);int nu

7、m_node,num_line,num_tran,num_gene,num_load,iter;struct data *line,*tran,*gene,*load;double *G,*B,*Jacob;double *Um,*Ua,*P,*Q;double error_max;FILE *fin,*fou,*chk;struct dataint i;int j;double a;double b;double c;void main()int i,j,conv;double a,error;fin=fopen("输入数据.txt","r");if(

8、fin=NULL)printf(" 注意! 没有“输入数据.txt”文件n");exit(0);fou=fopen("潮流输出.txt","w");if(pnt=1) chk=fopen("中间结果.txt","w");/ 输入原始数据和形成节点导纳矩阵read_data();G=newSpaceDouble2(num_node,num_node);B=newSpaceDouble2(num_node,num_node);for(i=1;i<=num_node;i+)for(j=1;j&l

9、t;=num_node;j+)Gij=Bij=0.0;admt_matrix();/ 给定电压有效值和相位初值Um=new doublenum_node+1;Ua=new doublenum_node+1;for(i=1;i<=num_node;i+)Umi=1.0;Uai=0.0;for(i=1;i<=num_gene;i+)if(genei.j<=0) Umgenei.i=genei.c;iter=0;/ 形成雅可比矩阵计算功率误差Jacob=newSpaceDouble2(2*num_node,2*num_node+1);P=new doublenum_node+1;Q

10、=new doublenum_node+1;R2:form_Jacobian();/ 收敛判断error=0.0;for(i=1;i<=2*num_node;i+)if(fabs(Jacobi2*num_node+1)>error)error=fabs(Jacobi2*num_node+1);fprintf(fou,"n 迭代次数: %2d 最大功率误差: %11.6f",iter+1,error);if(error<error_max) / 潮流计算收敛conv=1;goto R1;if(iter>10) | (error>1.0e4) /

11、潮流计算不收敛fprintf(fou,"nn 潮流不收敛");goto nd;/ 求解修正方程式并修正电压solv_Eqn();for(i=1;i<=num_node;i+)a=Jacobi2*num_node+1;Uai=Uai+(-1*a);a=Jacobnum_node+i2*num_node+1;Umi=Umi-(Umi*a);if(pnt=1) && (iter<2)fprintf(chk,"nn 电压相位和有效值新值, 迭代 %dnn",iter+1);for(i=1;i<=num_node;i+)fprin

12、tf(chk," %3d %10.5f %8.5f n",i,Uai,Umi);iter=iter+1;goto R2;/ 输出潮流结果R1:node_flow();branch_flow();nd:fclose(fin);fclose(fou);if(pnt=1) fclose(chk);free(line); free(tran); free(gene); free(load);deleteSpaceDouble2(G,num_node);deleteSpaceDouble2(B,num_node);deleteSpaceDouble2(Jacob,2*num_node

13、);free(Um); free(Ua); free(P); free(Q);void read_data() / 输入始数据int i;fscanf(fin,"%d %d %d %d %d %lf",&num_node,&num_line,&num_tran,&num_gene,&num_load,&error_max);line=(struct data *)calloc(num_line+1,sizeof(struct data);tran=(struct data *)calloc(num_tran+1,sizeof(

14、struct data);gene=(struct data *)calloc(num_gene+1,sizeof(struct data);load=(struct data *)calloc(num_load+1,sizeof(struct data);for(i=1;i<=num_line;i+) fscanf(fin,"%d %d %lf %lf %lf", &linei.i,&linei.j,&linei.a,&linei.b,&linei.c);for(i=1;i<=num_tran;i+) fscanf(fi

15、n,"%d %d %lf %lf %lf",&trani.i,&trani.j,&trani.a,&trani.b,&trani.c);for(i=1;i<=num_gene;i+) fscanf(fin,"%d %d %lf %lf %lf",&genei.i,&genei.j,&genei.a,&genei.b,&genei.c);for(i=1;i<=num_load;i+) fscanf(fin,"%d%lf %lf",&loa

16、di.i,&loadi.a,&loadi.b);void admt_matrix() / 形成导纳矩阵int i,j;double r,x,b,kt;struct data *p,*end;end=line+num_line;for(p=line+1;p<=end;p+) / 线路i=p->i;j=p->j;r=p->a;x=p->b;b=r*r+x*x;r=r/b;x=-x/b;if(i=j)Gii+=r;Bii+=x;continue;b=p->c;Gij=Gij-r;Bij=Bij-x;Gji=Gij;Bji=Bij;Gii=Gii+r

17、;Bii=Bii+x+b;Gjj=Gjj+r;Bjj=Bjj+x+b;end=tran+num_tran; / 变压器for(p=tran+1;p<=end;p+)i=p->i;j=p->j;r=p->a;x=p->b;b=r*r+x*x;r=r/b;x=-x/b;kt=p->c;Gii+=r;Bii+=x;Gij=Gij-r/kt;Bij=Bij-x/kt;Gji=Gij;Bji=Bij;r=r/kt/kt;x=x/kt/kt;Gjj+=r;Bjj+=x;if(pnt=1)fprintf(chk,"nn 导纳矩阵中的非零元素nn");f

18、or(i=1;i<=num_node;i+) for(j=i;j<=num_node;j+) if(Gij!=0.0) | (Bij!=0.0)fprintf(chk," %3d %4d %16.5f %16.5fn",i,j,Gij,Bij); void form_Jacobian() / 形成雅可比矩阵, 计算功率误差int i,j,nu,ii,k,n2,kk;double vi,di,dij,vj,dj,p,q,b,g,gp,gq,lp,lq;double Hij,Lij,Nij,Mij,Hii,Nii,Mii,Lii,dp,dq;nu=2*num_nod

19、e+1;n2=2*num_node;for(i=1;i<=n2;i+)for(j=1;j<=nu;j+)Jacobij=0.0;for(i=1;i<=num_node;i+)vi=Umi;di=Uai;dp=0.0;dq=0.0;for(j=1;j<=num_node;j+) / 非对角元素if(j=i) continue;g=Gij;b=Bij;vj=Umj;dj=Uaj;dij=di-dj;Hij=-Umi*Umj*(g*sin(dij)-b*cos(dij);Lij=Hij;Jacobij=Hij;Jacobi+num_nodej+num_node=Lij;Nij

20、=-Umi*Umj*(g*cos(dij)+b*sin(dij);Mij=-Nij;Jacobij+num_node=Nij;Jacobi+num_nodej=Mij;p=Umj*(g*cos(dij)+b*sin(dij);q=Umj*(g*sin(dij)-b*cos(dij);dp=dp+p;dq=dq+q;g=Gii;b=Bii;Hii=vi*dq;Nii=-vi*dp-2*vi*vi*g;Mii=-vi*dp;Lii=-vi*dq+2*vi*vi*b;Jacobii=Hii;Jacobii+num_node=Nii;Jacobi+num_nodei=Mii;Jacobi+num_no

21、dei+num_node=Lii;Jacobinu=-vi*(dp+vi*g);Jacobi+num_nodenu=-vi*(dq-vi*b);Pi=vi*(dp+vi*g); / 节点注入有功Qi=vi*(dq-vi*b); / 节点注入无功 for(i=1;i<=num_load;i+)kk=loadi.i;lp=loadi.a;lq=loadi.b;Jacobkknu=-lp+Jacobkknu;Jacobkk+num_nodenu=-lq+Jacobkk+num_nodenu;for(i=1;i<=num_gene;i+)kk=genei.i;gp=genei.a;gq=g

22、enei.b;Jacobkknu=gp+Jacobkknu;Jacobkk+num_nodenu=gq+Jacobkk+num_nodenu;for(k=1;k<=num_gene;k+) / 去掉pv及平衡节点ii=genek.i;kk=genek.j;if(kk=0) / 平衡节点for(j=1;j<=n2;j+)Jacobiij=0.0;Jacobnum_node+iij=0.0;Jacobjii=0.0;Jacobjnum_node+ii=0.0;Jacobiiii=1.0;Jacobnum_node+iinum_node+ii=1.0;Jacobiinu=0.0;Jaco

23、bnum_node+iinu=0.0;if(kk<0) / pv节点for(j=1;j<=n2;j+)Jacobnum_node+iij=0.0;Jacobjnum_node+ii=0.0;Jacobnum_node+iinum_node+ii=1.0;Jacobnum_node+iinu=0.0;if(pnt=1) && (iter<2)fprintf(chk,"nn 雅可比矩阵中的非零元素, 迭代 %dnn",iter+1);for(i=1;i<=2*num_node;i+) for(j=1;j<=2*num_node+1;

24、j+) if(Jacobij!=0.0) && (Jacobij!=1.0)fprintf(chk," %3d %3d %15.5f n",i,j,Jacobij); void solv_Eqn() / 求解修正方程式int i,j,n2,nu,i1,k;double d,e;n2=2*num_node;nu=n2+1;for(i=1;i<=n2;i+) / 消去i1=i+1;d=1.0/Jacobii;for(j=i1;j<=nu;j+)e=Jacobij;if(e=0.0) continue;Jacobij=e*d;if(i=n2) cont

25、inue;for(j=i1;j<=n2;j+)e=Jacobji;if(e=0.0) continue;for(k=i1;k<=nu;k+)Jacobjk=Jacobjk-Jacobik*e;for(k=2;k<=n2;k+) / 回代i=n2-k+1;i1=i+1;for(j=i1;j<=n2;j+)Jacobinu=Jacobinu-Jacobij*Jacobjnu;if(pnt=1) && (iter<2)fprintf(chk,"nn 电压相位和有效值修正量, 迭代 %dnn",iter+1);for(i=1;i<

26、=num_node;i+) fprintf(chk," %3d %10.5f %8.5f n",i,-Jacobi2*num_node+1,-Jacobi+num_node2*num_node+1);void node_flow() / 输出节点潮流int i,j,k,ii,kk;double b1,b2,c1,c2;fprintf(fou,"nnn *-*-*- 潮 流 计 算 结 果 *-*-*-");fprintf(fou,"nn 节 点 潮 流 ");fprintf(fou,"nn no.i Um Ua PG QG&

27、quot;);fprintf(fou," PL QLnn");for(i=1;i<=num_node;i+)b1=b2=c1=c2=0.0;for(j=1;j<=num_gene;j+)ii=genej.i;kk=genej.j;if(i=ii) && (kk=0) / 平衡节点b1=Pii;b2=Qii;for(k=1;k<=num_load;k+)ii=loadk.i;if(i=ii)c1=loadk.a;c2=loadk.b;b1=b1+c1;b2=b2+c2;break;if(i=ii) && (kk=-1) / p

28、v节点b1=genej.a;b2=Qii;for(k=1;k<=num_load;k+)ii=loadk.i;if(i=ii)c1=loadk.a;c2=loadk.b;b2=b2+c2;break;for(j=1;j<=num_load;j+)ii=loadj.i;if(i=ii)c1=loadj.a;c2=loadj.b;break;fprintf(fou," %6d %10.5f %10.5f %10.5f %10.5f %10.5f %10.5fn",i,Umi,Uai*180.0/3.141592654,b1,b2,c1,c2);void brabch

29、_flow() / 输出支路潮流int i,j;double r,x,t,b,dij,cd,sd,ri,rj,xi,xj;double vi,vj,vij,pij,qij,pji,qji,dpb,dqb,ph,qh;struct data *p,*end;fprintf(fou,"nn 支 路 潮 流 ");fprintf(fou,"nn i j Pij Qij Pji");fprintf(fou," Qji dP dQnn");ph=0.0;qh=0.0;for(end=line+num_line,p=line+1;p<=en

30、d;p+)i=p->i;j=p->j;r=p->a;x=p->b;b=r*r+x*x;if(i=j)vi=Umi;b=vi*vi/b;pij=r*b;qij=x*b;pji=0.0;qji=0.0;dpb=pij;ph=ph+dpb;dqb=qij;qh=qh+dqb;elser=r/b;x=-x/b;b=p->c;dij=Uai-Uaj;vi=Umi;vj=Umj;vij=vi*vj;vi=vi*vi;vj=vj*vj;cd=vij*cos(dij);sd=vij*sin(dij);pij=vi*r-r*cd-x*sd;pji=vj*r-r*cd+x*sd;dp

31、b=pij+pji;ph=ph+dpb;qij=-vi*(b+x)+x*cd-r*sd;qji=-vj*(b+x)+x*cd+r*sd;dqb=qij+qji;qh=qh+dqb;fprintf(fou," %3d %3d %10.5f %10.5f %10.5f %10.5f %10.5f %10.5fn",i,j,pij,qij,pji,qji,dpb,dqb);for(end=tran+num_tran,p=tran+1;p<=end;p+)i=p->i;j=p->j;r=p->a;x=p->b;t=p->c;b=t*(r*r+x*

32、x);r/=b;x/=-b;b=t-1.0;ri=r*b;xi=x*b;rj=-ri/t;xj=-xi/t;vi=Umi;vj=Umj;vij=vi*vj;vi*=vi;vj*=vj;dij=Uai-Uaj;cd=vij*cos(dij);sd=vij*sin(dij);pij=vi*(ri+r)-r*cd-x*sd;pji=vj*(rj+r)-r*cd+x*sd;dpb=pij+pji;ph+=dpb;qij=-vi*(xi+x)+x*cd-r*sd;qji=-vj*(xj+x)+x*cd+r*sd;dqb=qij+qji; qh+=dqb;fprintf(fou," %3d %3

33、d %10.5f %10.5f %10.5f %10.5f %10.5f %10.5fn",i,j,pij,qij,pji,qji,dpb,dqb);fprintf(fou,"nn 系统总损耗:- 有功功率:%8.5f 无功功率:%8.5f",ph,qh);/ 分配二维双精度型数组double *newSpaceDouble2(int n1,int n2)int i;double *a=new double*n1+1;for(i=0;i<=n1;i+)ai=new doublen2+1;return a;/ 释放二维双精度型数组void deleteSpaceDouble2(doubl