1、附錄附錄 牛頓-拉夫遜法潮流計(jì)算程序本附錄所介紹的牛頓-拉夫遜法潮流計(jì)算程序采用極坐標(biāo)形式，其中所涉及的計(jì)算公式與第四章中的式(4-42)-(4-59)完全相同，計(jì)算程序框圖與圖4-6基本一致。程序采用C語(yǔ)言編制。為了便于初學(xué)者閱讀，節(jié)點(diǎn)導(dǎo)納矩陣和雅可比矩陣都用滿陣存儲(chǔ)而未采用稀疏技巧；但為了適當(dāng)照顧使用的方便性，在輸入數(shù)據(jù)中對(duì)節(jié)點(diǎn)編號(hào)次序不作任何要求。下面先介紹輸入文件的格式和要求，然后列出程序，最后說(shuō)明潮流結(jié)果的輸出。建議讀者先從原始數(shù)據(jù)的輸入中了解和熟悉它們?cè)诔绦蛑袑?duì)應(yīng)的變量、結(jié)構(gòu)體數(shù)組及其成員的名稱，然后對(duì)照?qǐng)D4-6和式(4-42)-(4-59)仔細(xì)和耐心地閱讀程序，最好能在計(jì)算機(jī)上親

2、自實(shí)現(xiàn)和進(jìn)行調(diào)試。在調(diào)試時(shí)，可以用下面給出的對(duì)應(yīng)于例4-3系統(tǒng)的輸入數(shù)據(jù)，以及程序運(yùn)行中得出的中間結(jié)果，逐步與例4-3中所給出的中間結(jié)果進(jìn)行對(duì)比，從而查找錯(cuò)誤所在并進(jìn)行改正。一、原始數(shù)據(jù)的輸入程序通過(guò)“輸入數(shù)據(jù).txt”文件輸入以下5個(gè)數(shù)據(jù)段。1. 信息(共6個(gè))(1) 總節(jié)點(diǎn)數(shù)(變量num_node)(2) 線路和并聯(lián)電容器總數(shù)(變量num_line)(3) 變壓器支路總數(shù)(變量num_tran)(4) 發(fā)電機(jī)節(jié)點(diǎn)總數(shù)(變量num_gene)(5) 負(fù)荷節(jié)點(diǎn)總數(shù)(變量num_load)(6) 節(jié)點(diǎn)功率不平衡量的容許誤差(變量error)2. 線路和并聯(lián)電容器數(shù)據(jù)(結(jié)構(gòu)體數(shù)組line)：每一線

3、路或并聯(lián)電容器包括5個(gè)數(shù)據(jù)線路并聯(lián)電容器成員I側(cè)節(jié)點(diǎn)編號(hào)所接節(jié)點(diǎn)編號(hào)iJ側(cè)節(jié)點(diǎn)編號(hào)同上j型等值電路電阻電容器電阻a型等值電路電抗電容器電抗(負(fù)數(shù))b型等值電路一端電納0.0cI側(cè)節(jié)點(diǎn)編號(hào)和J側(cè)節(jié)點(diǎn)編號(hào)可以對(duì)換；線路和并聯(lián)電容器之間的次序可以任意，而且允許多條線路或多個(gè)電容器并聯(lián)。3. 變壓器支路數(shù)據(jù)(結(jié)構(gòu)體數(shù)組tran)：每一變壓器支路包括5個(gè)數(shù)據(jù)變壓器支路成員1側(cè)節(jié)點(diǎn)編號(hào)i2側(cè)節(jié)點(diǎn)編號(hào)j電阻a電抗b非標(biāo)準(zhǔn)變比c變壓器電阻、電抗和非標(biāo)準(zhǔn)變比與兩側(cè)節(jié)點(diǎn)編號(hào)之間的關(guān)系服從圖2-29,即電阻和電抗在1側(cè)而非標(biāo)準(zhǔn)變比在2側(cè)；三繞組變壓器需按圖2-26化成3個(gè)變壓器支路，其中3側(cè)變壓器支路的非標(biāo)準(zhǔn)變比為

4、1。變壓器支路之間的次序可以任意。4. 發(fā)電機(jī)節(jié)點(diǎn)數(shù)據(jù)(結(jié)構(gòu)體數(shù)組gene)：每一發(fā)電機(jī)節(jié)點(diǎn)包括5個(gè)數(shù)據(jù)發(fā)電機(jī)節(jié)點(diǎn)成員所在節(jié)點(diǎn)編號(hào)i節(jié)點(diǎn)種類(lèi)j發(fā)出有功功率a發(fā)出無(wú)功功率b電壓c對(duì)于PQ節(jié)點(diǎn)，節(jié)點(diǎn)種類(lèi)為1，電壓可給任意值；對(duì)于PV節(jié)點(diǎn)，節(jié)點(diǎn)種類(lèi)為-1，發(fā)出無(wú)功功率可給任意值；對(duì)于平衡節(jié)點(diǎn)，節(jié)點(diǎn)種類(lèi)為0，發(fā)出有功功率和發(fā)出無(wú)功功率都可給任意值。發(fā)電機(jī)節(jié)點(diǎn)之間的次序可以任意。5. 負(fù)荷節(jié)點(diǎn)數(shù)據(jù)(結(jié)構(gòu)體數(shù)組load)：每一負(fù)荷節(jié)點(diǎn)包括3個(gè)數(shù)據(jù)負(fù)荷節(jié)點(diǎn)成員所在節(jié)點(diǎn)編號(hào)i吸收有功功率a吸收無(wú)功功率b負(fù)荷機(jī)節(jié)點(diǎn)之間的次序可以任意。對(duì)于例4-3中的電力系統(tǒng)，輸入數(shù)據(jù)如下：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二、計(jì)算程序 /* * * * * * 牛頓-拉夫遜法潮流計(jì)算程序 * * * * * */# include <stdio.h># include <math.h>#

6、include <stdlib.h># define pnt 1 / 1 - 輸出中間結(jié)果void read_data(); / 輸入始數(shù)據(jù)void admt_matrix(); / 形成導(dǎo)納矩陣void form_Jacobian(); / 形成雅可比矩陣, 計(jì)算功率誤差void solv_Eqn(); / 求解修正方程式void node_flow(); / 輸出節(jié)點(diǎn)潮流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("輸入數(shù)據(jù).txt","r");if(

8、fin=NULL)printf(" 注意! 沒(méi)有“輸入數(shù)據(jù).txt”文件n");exit(0);fou=fopen("潮流輸出.txt","w");if(pnt=1) chk=fopen("中間結(jié)果.txt","w");/ 輸入原始數(shù)據(jù)和形成節(jié)點(diǎn)導(dǎo)納矩陣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;/ 形成雅可比矩陣計(jì)算功率誤差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 迭代次數(shù): %2d 最大功率誤差: %11.6f",iter+1,error);if(error<error_max) / 潮流計(jì)算收斂conv=1;goto R1;if(iter>10) | (error>1.0e4) /

11、潮流計(jì)算不收斂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;/ 輸出潮流結(jié)果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() / 輸入始數(shù)據(jù)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() / 形成導(dǎo)納矩陣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 導(dǎo)納矩陣中的非零元素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() / 形成雅可比矩陣, 計(jì)算功率誤差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+) / 非對(duì)角元素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); / 節(jié)點(diǎn)注入有功Qi=vi*(dq-vi*b); / 節(jié)點(diǎn)注入無(wú)功 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及平衡節(jié)點(diǎn)ii=genek.i;kk=genek.j;if(kk=0) / 平衡節(jié)點(diǎn)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節(jié)點(diǎn)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() / 輸出節(jié)點(diǎn)潮流int i,j,k,ii,kk;double b1,b2,c1,c2;fprintf(fou,"nnn *-*-*- 潮 流 計(jì) 算 結(jié) 果 *-*-*-");fprintf(fou,"nn 節(jié) 點(diǎn) 潮 流 ");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) / 平衡節(jié)點(diǎn)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節(jié)點(diǎn)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 系統(tǒng)總損耗:- 有功功率:%8.5f 無(wú)功功率:%8.5f",ph,qh);/ 分配二維雙精度型數(shù)組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;/ 釋放二維雙精度型數(shù)組void deleteSpaceDouble2(doubl