%电力系统的潮流计算 ,以下程序参考文献 %(该文献用极坐标下的牛顿——拉夫逊方法实现，在此为了与课本一致做了修改 《电力系统毕业设计》中国水利电力出版社 ) %为了计算方便将原来的下标做以下修改 U2 %U3 U4 , 即原题的平衡点 1 就变为现在的平衡点 5 : S2 S3 S4 S5 U2 U3 U4 U5 改为 S1 S2 S3 S4 U1 %1.形成节点导纳矩阵 , yb55=6.250-18.750j;yb51=-5.000+15.000j;yb52=-1.250+3.750j;yb53=0.000-0.000j;yb54=0.000- 0.000j; yb15=-5.000+15.000j;yb11=10.834-32.500j;yb12=-1.667+5.000j;yb13=-1.667+5.000j;yb14=-2.5 00+7.500j; yb25=-1.250+3.750j;yb21=-1.667+5.000j;yb22=12.917-38.750j;yb23=-10.000+30.000j;yb24=0.0 00-0.000j; yb35=0.000-0.000j;yb31=-1.667+5.000j;yb32=-10.000+30.000j;yb33=12.917-38.750j;yb34=-1.2 50+3.750j; yb45=0.000-0.000j;yb41=-2.500+7.500j;yb42=0.000-0.000j;yb43=-1.250+3.750j;yb44=3.750-11. 250j; YB=[yb11 yb12 yb13 yb14 yb15; yb21 yb22 yb23 yb24 yb25 ;yb31 yb32 yb33 yb34 yb35; yb41 yb42 yb43 yb44 yb45 ;yb51 yb52 yb53 yb54 yb55]; %计算各节点功率的不平衡量设 U=E+jF ； Y=G+Bj ； E(1)=1.00;E(2)=1.00;E(3)=1.00;E(4)=1.00; F(1)=0;F(2)=0;F(3)=0;F(4)=0; G=real(YB);B=imag(YB); %设 S=P+Bj ； S(1)=0.20+0.20i;S(2)=-0.45-0.15i;S(3)=-0.40-0.05i;S(4)=-0.60-0.10i; P=real(S);Q=imag(S); k=0;precision=1; N1=4; while precision > 0.00000001 E(5)=1.06;F(5)=0; for m=1:N1 for n=1:N1+1 %计算 Pi,Qi,设 Pi=Pt;Qi=Qt 

Pt(n)=(E(m)*(G(m,n)*E(n)-B(m,n)*F(n))+F(m)*(G(m,n)*F(n)+B(m,n)*E(n))); Qt(n)=(F(m)*(G(m,n)*E(n)-B(m,n)*F(n))-E(m)*(G(m,n)*F(n)+B(m,n)*E(n))); end %设 P,Q 的改变量为 dP,dQ dP(m)=P(m)-sum(Pt); dQ(m)=Q(m)-sum(Qt); end for m=1:N1 for n=1:N1+1 %计算 Hij Nij Jij Lij H(m,n)=-B(m,n)*E(m)+G(m,n)*F(m); N(m,n)=G(m,n)*E(m)+B(m,n)*F(m); J(m,n)=-B(m,n)*F(m)-G(m,n)*E(m); L(m,n)=G(m,n)*F(m)-B(m,n)*E(m); end end for m=1:N1 for n=1:N1+1 Bi(n)=G(m,n)*F(n)+B(m,n)*E(n); Ai(n)=G(m,n)*E(n)-B(m,n)*F(n); end %计算 Hii,Nii,Jii,Lii, 由公式 4-44b 左侧公式实现 ,sum(Ai) ，sum(Bi) 用于实现公式中 的 sigerma 从 j 到 n 的求和 ; H(m,m)=sum(Bi)-(B(m,m)*E(m)+G(m,m)*F(m))+2*G(m,m)*F(m); N(m,m)=sum(Ai)-(G(m,m)*E(m)-B(m,m)*F(m))+2*G(m,m)*E(m); J(m,m)=-2*B(m,m)*F(m)+sum(Ai)-(G(m,m)*E(m)-B(m,m)*F(m)); L(m,m)=-2*B(m,m)*E(m)-(sum(Bi)-(B(m,m)*E(m)+G(m,m)*F(m))); end 

%设雅可比矩阵为 JJ，以下语句用来实现雅可比矩阵中对角线上元素 for m=1:N1 H N J L 的排列 JJ(2*m-1,2*m-1)=H(m,m); JJ(2*m-1,2*m)=N(m,m); JJ(2*m,2*m-1)=J(m,m); JJ(2*m,2*m)=L(m,m); end %以下语句用于实现雅可比矩阵非对角线上元素的排列 for m=1:N1 for n=1:N1 if m==n else H(m,n)=-B(m,n)*E(m)+G(m,n)*F(m); N(m,n)=G(m,n)*E(m)+B(m,n)*F(m); J(m,n)=-B(m,n)*F(m)-G(m,n)*E(m); L(m,n)=G(m,n)*F(m)-B(m,n)*E(m); JJ(2*m-1,2*n-1)=H(m,n); JJ(2*m-1,2*n)=N(m,n); JJ(2*m,2*n-1)=J(m,n); JJ(2*m,2*n)=L(m,n); end end end %设由 P,Q 的改变量组成的 8× 1 矩阵为 PQ，由 E,F 的改变量组成的 8× 1 矩阵为 dU for m=1:N1 PQ(2*m-1)=dP(m); PQ(2*m)=dQ(m); end dU=inv(JJ)*PQ'; precision=max(abs(dU)); 

for n=1:N1 F(n)=F(n)+dU(2*n-1); E(n)=E(n)+dU(2*n); end for n=1:N1+1 U(n)=E(n)+(F(n))*j; end k=k+1; k-1, dU=dU',PQ,U end %计算 S(5),也就是题目中的 S1，即平衡节点功率 for m=1:N1+1 I(m)=YB(5,m)*U(m); end S(5)=U(5)*sum(conj(I)) %设网络总损耗为 Ss,计算输电效率 efficiency for m=1:N1+1 S0(m)=S(m); P(m)=real(S(m)); end Ss=sum(S0) efficiency=(abs(P(3)+P(4)+P(2)))/(P(5)+(P(1)))*100 %计算线功率 S(m,n), 与课本中各元素的相对位置有所不同 for m=1:N1+1 for n=1:N1+1 S(m,n)=U(m)*(conj(U(m))-conj(U(n)))*conj(-YB(m,n)); end end 

S