牛頓拉夫遜法潮流計(jì)算matlab程序.doc

%電力系統(tǒng)的潮流計(jì)算,以下程序參考文獻(xiàn) 電力系統(tǒng)畢業(yè)設(shè)計(jì)中國水利電力出版社%(該文獻(xiàn)用極坐標(biāo)下的牛頓拉夫遜方法實(shí)現(xiàn)，在此為了與課本一致做了修改)%為了計(jì)算方便將原來的下標(biāo)做以下修改: S2 S3 S4 S5 U2 U3 U4 U5 改為S1 S2 S3 S4 U1 U2%U3 U4 ,即原題的平衡點(diǎn)1就變?yōu)楝F(xiàn)在的平衡點(diǎn)5%1.形成節(jié)點(diǎn)導(dǎo)納矩陣,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.500+7.500j;yb25=-1.250+3.750j;yb21=-1.667+5.000j;yb22=12.917-38.750j;yb23=-10.000+30.000j;yb24=0.000-0.000j;yb35=0.000-0.000j;yb31=-1.667+5.000j;yb32=-10.000+30.000j;yb33=12.917-38.750j;yb34=-1.250+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;%計(jì)算各節(jié)點(diǎn)功率的不平衡量設(shè)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);%設(shè)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 %計(jì)算Pi,Qi,設(shè)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 %設(shè)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 %計(jì)算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 %計(jì)算Hii,Nii,Jii,Lii,由公式4-44b 左側(cè)公式實(shí)現(xiàn),sum(Ai)，sum(Bi)用于實(shí)現(xiàn)公式中的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 %設(shè)雅可比矩陣為JJ，以下語句用來實(shí)現(xiàn)雅可比矩陣中對角線上元素H N J L 的排列 for m=1:N1 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 %以下語句用于實(shí)現(xiàn)雅可比矩陣非對角線上元素的排列 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 %設(shè)由P,Q的改變量組成的81矩陣為PQ，由E,F的改變量組成的81矩陣為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,Uend%計(jì)算S(5),也就是題目中的S1，即平衡節(jié)點(diǎn)功率for m=1:N1+1 I(m)=YB(5,m)*U(m);endS(5)=U(5)*sum(conj(I)%設(shè)網(wǎng)絡(luò)總損耗為Ss,計(jì)算輸電效率 efficiencyfor m=1:N1+1 S0(m)=S(m); P(m)=real(S(m);endSs=sum(S0)efficiency=(a