C语言计算潮流程序.docVIP

节点数：3 支路数：3 计算精度：0.00010支路1： 0.0300+j0.0900 12支路2： 0.0200+j0.0900 23支路3： 0.0300+j0.0900 31节点1：PQ节点，S（1）=-0.5000-j0.2000节点2：PQ节点，S（2）=-0.6000-j0.2500节点3：平衡节点，U（3）=1.00000.0000n=5;nl=5;isb=1;pr=0.00001; B1=1 2 0.03i 0 1.05 0;2 3 0.08+0.3i 0.5i 1 0;2 4 0.1+0.35i 0 1 0;3 4 0.04+0.25i 0.5i 1 0;3 5 0.015i 0 1.05 1;B2=0 0 1.05 1.05 0 1;0 3.7+1.3i 1.05 0 0 2;0 2+1i 1.05 0 0 2;0 1.6+0.8i 1.05 0 0 2;5 0 1.05 1.05 0 3;X=1 0;2 0;3 0;4 0;5 0;na=3;Y=zeros(n);YI=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);O=zeros(1,n);for i=1:n if X(i,2)=0; p=X(i,1); Y(p,p)=1./X(i,2); endendfor i=1:nl if B1(i,6)=0 p=B1(i,1);q=B1(i,2); else p=B1(i,2);q=B1(i,1); end Y(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5); YI(p,q)=YI(p,q)-1./B1(i,3); Y(q,p)=Y(p,q); YI(q,p)=YI(p,q); Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)2)+B1(i,4)./2; YI(q,q)=YI(q,q)+1./B1(i,3); Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2; YI(p,p)=YI(p,p)+1./B1(i,3);endG=real(Y);B=imag(YI);BI=imag(Y);for i=1:n S(i)=B2(i,1)-B2(i,2); BI(i,i)=BI(i,i)+B2(i,5);endP=real(S);Q=imag(S);for i=1:n e(i)=real(B2(i,3); f(i)=imag(B2(i,3); V(i)=B2(i,4);endfor i=1:n if B2(i,6)=2 V(i)=sqrt(e(i)2+f(i)2); O(i)=atan(f(i)./e(i); endendfor i=2:n if i=n B(i,i)=1./B(i,i); else IC1=i+1; for j1=IC1:n B(i,j1)=B(i,j1)./B(i,i); end B(i,i)=1./B(i,i); for k=i+1:n for j1=i+1:n B(k,j1)=B(k,j1)-B(k,i)*B(i,j1); end end endendp=0;q=0;for i=1:n if B2(i,6)=2 p=p+1;k=0; for j1=1:n if B2(j1,6)=2 k=k+1; A(p,k)=BI(i,j1); end end endendfor i=1:na if i=na A(i,i)=1./A(i,i); else k=i+1; for j1=k:na A(i,j1)=A(i,j1)./A(i,i); end A(i,i)=1./A(i,i); for k=i+1:na for j1=i+1:na A(k,j1)=A(k,j1)-A(k,i)*A(i,j1); end end endendICT2=1;ICT1=0;kp=1;kq=1;K=1;DET=0;ICT3=1;while ICT2=0|ICT3=0ICT2=0;ICT3=0;for i=1:n if i=isb C(i)=0; for k=1:n C(i)=C(i)+V(k)*(G(i,k)*cos(O(i)-O(k)+BI(i,k)*sin(O(i)-O(k); end DP1(i)=P(i)-V(i)*C(i); DP(i)=DP1(i)./V(i); DET=abs(DP1(i); if DET=pr ICT2=ICT2+1; end endendNp(K)=ICT2; if ICT2=0 for i=2:n DP(i)=B(i,i)*DP(i); if i=n IC1=i+1; for k=IC1:n DP(k)=DP(k)-B(k,i)*DP(i); end else for LZ=3:i L=i+3-LZ; IC4=L-1; for MZ=2:IC4 I=IC4+2-MZ; DP(I)=DP(I)-B(I,L)*DP(L); end end end end for i=2:n O(i)=O(i)-DP(i); end kq=1;L=0; for i=1:n if B2(i,6)=2 C(i)=0;L=L+1; for k=1:n C(i)=C(i)+V(k)*(G(i,k)*sin(O(i)-O(k)-BI(i,k)*cos(O(i)-O(k); end DQ1(i)=Q(i)-V(i)*C(i); DQ(L)=DQ1(i)./V(i); DET=abs(DQ1(i); if DET=pr ICT3=ICT3+1; end end end else kp=0; if kq=0; L=0; for i=1:n if B2(i,6)=2 C(i)=0;L=L+1; for k=1:n C(i)=C(i)+V(k)*(G(i,k)*sin(O(i)-O(k)-BI(i,k)*cos(O(i)-O(k); end DQ1(i)=Q(i)-V(i)*C(i); DQ(L)=DQ1(i)./V(i); DET=abs(DQ1(i); end end end end Nq(K)=ICT3; if ICT3=0 L=0; for i=1:na DQ(i)=A(i,i)*DQ(i); if i=na for LZ=2:i L=i+2-LZ; IC4=L-1; for MZ=1:IC4 I=IC4+1-MZ; DQ(I)=DQ(I)-A(I,L)*DQ(L); end end else IC1=i+1; for k=IC1:na DQ(k)=DQ(k)-A(k,i)*DQ(i); end end end L=0; for i=1:n if B2(i,6)=2 L=L+1; V(i)=V(i)-DQ(L); end end kp=1; K=K+1; else kq=0; if kp=0 K=K+1; end end for i=1:n Dy(K-1,i)=V(i); endenddisp(迭代次数)disp(K);disp(每次没有达到精度要求的有功功率个数为);disp(Np);disp(每次没有达到精度要求的无功功率个数为);disp(Nq);for k=1:n E(k)=V(k)*cos(O(k)+V(k)*sin(O(k)*j; O(k)=O(k)*180./pi;enddisp(各节点的电压标么值E为);disp(E);disp(各节点的电压V大小);disp(V);disp(各节点的电压相角O);disp(O);for p=1:n C(p)=0; for q=1:n C(p)=C(p)+conj(Y(p,q)*conj(E(q); end S(p)=E(p)*C(p);enddisp(各节点的功率为);disp(S);disp(各条支路的首端功率为);for i=1:nl if B1(i,6)=0 p=B1(i,1);q=B1(i,2); else p=B1(i,2);q=B1(i,1); end Si(p,q)=E(p)*(conj(E(p)*conj(B1(i,4)./2)+(conj(E(p)*B1(i,5)-conj(E(q)*conj(1./(B1(i,3)*B1(i,5); disp(Si(p,q);enddisp(各条支路的末端功率为);for i=1:nl if B1(i,6)=0 p=B1(i,1);q=B1(i,2); else p=B1(i,2);q=B1(i,1); end Sj(q,p)=E(q)*(conj(E(q)*conj(B1(i,4)./2)+(conj(E(q)./B1(i,5)-conj(E(p)*conj(1./(B1(i,3)*B1(i,5); disp(Sj(q,p);enddisp(各条支路的功率损耗为);for i=1:nl if B1(i,6)=0 p=B1(i,1);q=B1(i,2); else p=B1(i,2);q=B1(i,1); end DS(i)=Si(p,q)+Sj(q,p); disp(DS(i);endfor i=1:K Cs(i)=i; for j=1:n Dy(K,j)=Dy(K-1,j); endenddisp(每次迭代后各节点的电压值如图所示);plot(Cs,Dy)xlabel(迭代次数)ylabel(电压)title(电压迭代次数曲线);题号：2 节点数：4 支路数：4 计算精度：0.00010 支路 1： 0.0200+j0.0800 13支路 2： 0.0400+j0.1200 14支路 3： 0.0500+j0.1400 24 支路 4： 0.0400+j0.1200 34 节点1：PQ节点， S(1)=-0.6000-j0.2500节点2：PQ节点， S(2)=-0.8000-j0.3500节点3：PV节点， P(3)=0.4000 V(3)=0.9500节点4：平衡节点， U(4)=1.00000.0000#include stdio.h#include math.h#define M 20 /*节点数、支路数极限值*/#define N 5 /*迭代次数极限值*/int n,m,dd=0,pq=0,pv=0,tt,qy;float eps; /*节点数、支路数、对地支路数、pq节点数、pv节点数、控制打印、互联网数、精度*/struct jiedian /*节点数据结构体*/int s; /*节点类型：1-PQ节点；2-PV节点；3-平衡节点*/float p,q,e,f,v; /*节点的有功、无功、电压实部、虚部、电压辐值*/jiedianM; /*如引用节点1的无功，则为jiedian1.q */struct zhilu /*支路数据结构体*/int p1,p2,s; /*支路两端节点号，支路类型：1-普通支路；2-变压器支路；3-对地支路*/float r,x,b,kt; /*支路的电阻、电抗、导纳、变压器非标准变比*/zhiluM;struct hulianwang /*互联网结构体*/int num,pv; /*互联区域号,指定的pv节点号*/ float p,eps2; /*规定的有功功率及其允许误差*/ int count; /*每个互联区域包括的节点数*/ int aM; /*包括的节点*/ hulianwangM;static float GMM,BMM,G1MM,B1MM; /*节点导纳阵矩阵*/static float ykb2*M2*M; /*节点导纳阵矩阵*/static float yinzi2*M2*M; /*因子表*/static float PMM,QMM; /*潮流计算结果*/FILE *fp1,*fp2; /*文件输入、输出指针*/ void input() /*从文件in.txt 输入线路基本参数、节点数据、支路数据*/ int i,j,h; /*循环变量、节点（支路）类型*/ /*打开电网参数文件*/ fp1=fopen(in.txt,r); if (fp1=NULL) printf(Can not open file in.txt !n); exit(0); /*读取节点数,支路数,互联网数,PQ节点数PV节点数和精度*/ fscanf(fp1,%d,%d,%d,%fn,&n,&m,&qy,&eps); /*读取节点信息*/ for(i=1;i=n;i+) fscanf(fp1,%d,&h); if(h=1) /*类型h=1是PQ节点*/ pq+; fscanf(fp1,%f,%f,%f,%fn,&jiediani.p,&jiediani.q,&jiediani.e,&jiediani.f); jiediani.s=1; jiediani.v=sqrt(jiediani.e*jiediani.e+jiediani.f*jiediani.f); else if(h=2) /*类型h=2是pv节点*/ pv+; fscanf(fp1,%f,%f,%f,%fn,&jiediani.p,&jiediani.v,&jiediani.e,&jiediani.f); jiediani.s=2; jiediani.q=0; else /*类型h=3是平衡节点*/ jiediani.p=0; jiediani.q=0; jiediani.e=1; jiediani.f=0; jiediani.v=1; jiediani.s=3; /*读取支路信息*/ for(i=1;i=m;i+) fscanf(fp1,%d,&h); if(h=1) /*类型h=1是普通支路*/ fscanf(fp1,%d,%d,%f,%f,%fn,&zhilui.p1,&zhilui.p2,&zhilui.r,&zhilui.x,&zhilui.b); zhilui.kt=1; zhilui.s=1; if(h=2) /*类型h=2是变压器支路*/ fscanf(fp1,%d,%d,%f,%f,%fn,&zhilui.p1,&zhilui.p2,&zhilui.r,&zhilui.x,&zhilui.kt); zhilui.s=2; if(h=3) /*类型h=3是接地支路*/ fscanf(fp1,%d,%d,%f,%f,%fn,&zhilui.p1,&zhilui.p2,&zhilui.r,&zhilui.x,&zhilui.b); zhilui.kt=1; zhilui.s=3; dd+; /*读互联网信息*/ if(qy!=0) for(i=1;i=qy;i+) /*输入互联网状况*/ fscanf(fp1,%d,%d,%f,%f,%d,&hulianwangi.num,&hulianwangi.pv,&hulianwangi.p,&hulianwangi.eps2,&hulianwangi.count); for(j=1;j=hulianwangi.count;j+) fscanf(fp1,%d,&(hulianwangi.aj); fclose(fp1); /*打开输出结果文件*/ fp2=fopen(out.txt,w); if(fp2=NULL) printf(Can not open file!n); exit(0); fprintf(fp2,n* 原始数据如下 *nn ); fprintf(fp2, 节点数:%2d 支路数:%2d 对地支路数:%2d PQ节点数:%2d PV节点数:%2d 精度1:%f n,n,m,dd,pq,pv,eps); fprintf(fp2,n-nn ); for(i=1;i=pq;i+) fprintf(fp2, 节点%2d PQ节点 P%d=%f Q%d=%fn,i,i,jiediani.p,i,jiediani.q); for(i=pq+1;i=pq+pv;i+) fprintf(fp2, 节点%2d PV节点 P%d=%f V%d=%fn,i,i,jiediani.p,i,jiediani.v); fprintf(fp2, 节点%2d 平衡节点n,i); fprintf(fp2,n-nn ); for(i=1;i=m;i+) if(zhilui.s=1) fprintf(fp2, 支路%2d：普通支路 相关节点:%2d,%2d R=%f X=%f B=%fn,i,zhilui.p1,zhilui.p2,zhilui.r,zhilui.x,zhilui.b); else if(zhilui.s=2) fprintf(fp2, 支路%2d：变压器支路 相关节点:%2d,%2d R=%f X=%f Kt=%fn,i,zhilui.p1,zhilui.p2,zhilui.r,zhilui.x,zhilui.kt); else fprintf(fp2, 支路%2d： 对地支路 相关节点:%2d R=%f X=%f B=%fn,i,zhilui.p1,zhilui.r,zhilui.x,zhilui.b); for(i=1;i=qy;i+) fprintf(fp2, 互联网区域:%2d 指定pv节点:%2d 规定有功功率:%8.5f 允许误差:%8.5f 包括的节点: ,hulianwangi.num,hulianwangi.pv,hulianwangi.p,hulianwangi.eps2); for (j=1;j=hulianwangi.count;j+) fprintf(fp2,%2d ,hulianwangi.aj);fprintf(fp2,n); void youhua() /*利用节点数据、支路数据形成新编号*/ int i,k,j,jdN; /*jdN记录节点*/ struct jiedian tem; /*中间临时节点所连支路数*/ for(i=1;i=n;i+) /*对节点连接支路数赋0*/ jdi=0; for(i=1;i=n;i+) /*考虑每个节点*/ for(j=1;j=m;j+) /*考虑每条支路*/ if(zhiluj.p1=i|zhiluj.p2=i)&(zhiluj.s!=3) /*如果是这条支路的节点且不是接地支路*/ jdi+; for(i=1;ipq;i+) /*对pq节点按所连支路的个数进行重新编号*/ for(j=i+1;jjdj) tem=jiediani; jiediani=jiedianj; jiedianj=tem; for(k=1;k=m;k+) /*更新支路信息*/ if(zhiluk.p1=j) zhiluk.p1=i; else if(zhiluk.p2=j) zhiluk.p2=i; else if(zhiluk.p1=i) zhiluk.p1=j; else if(zhiluk.p2=i) zhiluk.p2=j; for(i=pq+1;ipq+pv;i+) /*对pv节点按所连支路的个数进行重新编号*/ for(j=i+1;jjdj) tem=jiediani; jiediani=jiedianj; jiedianj=tem; for(k=1;k=m;k+) /*更新支路信息*/ if(zhiluk.p1=j) zhiluk.p1=i; else if(zhiluk.p2=j) zhiluk.p2=i; else if(zhiluk.p1=i) zhiluk.p1=j; else if(zhiluk.p2=i) zhiluk.p2=j; fprintf(fp2,n* 节点优化结果如下 *nn ); for(i=1;i=n;i+) if(jiediani.s=1) fprintf(fp2, 节点%2d PQ节点 P%d=%f Q%d=%f n,i,i,jiediani.p,i,jiediani.q); if(jiediani.s=2) fprintf(fp2, 节点%2d PV节点 P%d=%f V%d=%f n,i,i,jiediani.p,i,jiediani.v); if(jiediani.s=3) fprintf(fp2, 节点%2d 平衡节点 n,i); void form_y() /*利用支路数据形成Y，注意对地支路*/ int i,j,k; float S; for(i=1;i=n;i+) for(j=1;j=n;j+) Gij=Bij=0; for(i=1;i=m;i+) /*节点导纳矩阵的主对角线上的导纳*/ for(j=1;j=n;j+) if(zhilui.p1=j|zhilui.p2=j) S=zhilui.r*zhilui.r+zhilui.x*zhilui.x; if(S=0) continue; Gjj+=zhilui.r/S; Bjj+=-zhilui.x/S; if(zhilui.s=1) /*如果是普通支路*/ Bjj+=zhilui.b/2; else if(zhilui.s=2) /*如果是普通变压器支路*/ if(zhilui.p1=j) Gjj+=(zhilui.r/S*(1-zhilui.kt)/(zhilui.kt*zhilui.kt); Bjj+=(-zhilui.x/S*(1-zhilui.kt)/(zhilui.kt*zhilui.kt); else if(zhilui.p2=j) Gjj+=(zhilui.r/S*(zhilui.kt-1)/zhilui.kt; Bjj+=(-zhilui.x/S*(zhilui.kt-1)/zhilui.kt; else if(zhilui.s=3) /*如果是对地支路*/ Bjj+=zhilui.b; for(k=1;k=m;k+) /*节点导纳矩阵非主对角线上的导纳*/ i=zhiluk.p1; j=zhiluk.p2; S=zhiluk.r*zhiluk.r+zhiluk.x*zhiluk.x; if(S=0) continue; Gij+=-zhiluk.r/S; Bij+=zhiluk.x/S; if(zhiluk.kt!=1.0) Gij/=zhiluk.kt; Bij/=zhiluk.kt; Gji=Gij; Bji=Bij; fprintf(fp2,nn* 节点导纳矩阵为 *n); for(i=1;i=n;i+) fprintf(fp2,n ); for(j=1;j=n;j+) fprintf(fp2,%8.5f+j%8.5f ,Gij,Bij); void form_j() /*利用节点数据和Y形成J*/ float ei,fi,a=0,b=0; int i1,j1,k1,i,j,k; for(i=1;i=2*(pq+pv)+1;i+) for(j=1;j=2*(pq+pv)+1;j+) ykbij=0; for(i=1;i=pq;i+) for(j=1;jn;j+) i1=i; j1=j; ei=jiediani.e; fi=jiediani.f; if(i!=j) /*求i!=j时的H、N、J、L*/ ykb2*i-12*j-1=Gi1j1*ei+Bi1j1*fi; /* H */ ykb2*i-12*j=-Bi1j1*ei+Gi1j1*fi; /* N */ ykb2*i2*j-1=-Bi1j1*ei+Gi1j1*fi; /* J */ ykb2*i2*j=-Gi1j1*ei-Bi1j1*fi; /* L */ else /*求i=j时的H、N、J、K*/ a=0;b=0; for(k=1;k=n;k+) if(k!=i) k1=k; a=a+Gi1k1*jiediank.e-Bi1k1*jiediank.f; b=b+Gi1k1*jiediank.f+Bi1k1*jiediank.e; ykb2*i-12*j-1=2*Gi1i1*jiediani.e+a; /*H*/ ykb2*i2*j=-2*Bi1i1*jiediani.f+a; /*L*/ ykb2*i-12*j= 2*Gi1i1*jiediani.f+b; /*N*/ ykb2*i2*j-1=-2*Bi1i1*jiedi