1、首先聲明一下，這些是從網(wǎng)站上轉(zhuǎn)載的，不是本人上編寫的基于MATLAB的電力系統(tǒng)潮流計算%簡單潮流計算的小程序，相關(guān)的原始數(shù)據(jù)數(shù)據(jù)數(shù)據(jù)輸入格式如下：%B1是支路參數(shù)矩陣，第一列和第二列是節(jié)點(diǎn)編號。節(jié)點(diǎn)編號由小到大編寫%對于含有變壓器的支路，第一列為低壓側(cè)節(jié)點(diǎn)編號，第二列為高壓側(cè)節(jié)點(diǎn)%編號，將變壓器的串聯(lián)阻抗置于低壓側(cè)處理。%第三列為支路的串列阻抗參數(shù)。%第四列為支路的對地導(dǎo)納參數(shù)。%第五烈為含變壓器支路的變壓器的變比%第六列為變壓器是否是否含有變壓器的參數(shù)，其中“1”為含有變壓器，%“0”為不含有變壓器。%B2為節(jié)點(diǎn)參數(shù)矩陣，其中第一列為節(jié)點(diǎn)注入發(fā)電功率參數(shù)；第二列為節(jié)點(diǎn)%負(fù)荷功率參數(shù)；第三列為

2、節(jié)點(diǎn)電壓參數(shù)；第六列為節(jié)點(diǎn)類型參數(shù)，其中%“1”為平衡節(jié)點(diǎn)，“2”為PQ節(jié)點(diǎn)，“3”為PV節(jié)點(diǎn)參數(shù)。%X為節(jié)點(diǎn)號和對地參數(shù)矩陣。其中第一列為節(jié)點(diǎn)編號，第二列為節(jié)點(diǎn)對地%參數(shù)。n=input(請輸入節(jié)點(diǎn)數(shù):n=);n1=input(請輸入支路數(shù):n1=);isb=input(請輸入平衡節(jié)點(diǎn)號:isb=);pr=input(請輸入誤差精度:pr=);B1=input(請輸入支路參數(shù):B1=);B2=input(請輸入節(jié)點(diǎn)參數(shù):B2=);X=input(節(jié)點(diǎn)號和對地參數(shù):X=);Y=zeros(n);Times=1; %置迭代次數(shù)為初始值%創(chuàng)建節(jié)點(diǎn)導(dǎo)納矩陣for i=1:n1 if B1(i,6)=

3、0 %不含變壓器的支路 p=B1(i,1); q=B1(i,2); Y(p,q)=Y(p,q)-1/B1(i,3); Y(q,p)=Y(p,q); Y(p,p)=Y(p,p)+1/B1(i,3)+0.5*B1(i,4); Y(q,q)=Y(q,q)+1/B1(i,3)+0.5*B1(i,4); else %含有變壓器的支路 p=B1(i,1); q=B1(i,2); Y(p,q)=Y(p,q)-1/(B1(i,3)*B1(i,5); Y(q,p)=Y(p,q); Y(p,p)=Y(p,p)+1/B1(i,3); Y(q,q)=Y(q,q)+1/(B1(i,5)2*B1(i,3); endend

4、YOrgS=zeros(2*n-2,1);DetaS=zeros(2*n-2,1); %將OrgS、DetaS初始化%創(chuàng)建OrgS，用于存儲初始功率參數(shù)h=0;j=0;for i=1:n %對PQ節(jié)點(diǎn)的處理 if i=isb&B2(i,6)=2 h=h+1; for j=1:n OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)+imag(B2(i,3)*(real(Y(i,j)*imag(B2(j,3)+imag(Y(i,j)*real(B2(j,3); Org

5、S(2*h,1)=OrgS(2*h,1)+imag(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)-real(B2(i,3)*(real(Y(i,j)*imag(B2(j,3)+imag(Y(i,j)*real(B2(j,3); end endendfor i=1:n %對PV節(jié)點(diǎn)的處理，注意這時不可再將h初始化為0 if i=isb&B2(i,6)=3 h=h+1; for j=1:n OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-im

6、ag(Y(i,j)*imag(B2(j,3)+imag(B2(i,3)*(real(Y(i,j)*imag(B2(j,3)+imag(Y(i,j)*real(B2(j,3); OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)-real(B2(i,3)*(real(Y(i,j)*imag(B2(j,3)+imag(Y(i,j)*real(B2(j,3); end endendOrgS%創(chuàng)建PVU 用于存儲PV節(jié)點(diǎn)的初始電壓PVU=zeros(n-h-1,1);t=0;fo

7、r i=1:n if B2(i,6)=3 t=t+1; PVU(t,1)=B2(i,3); endendPVU%創(chuàng)建DetaS，用于存儲有功功率、無功功率和電壓幅值的不平衡量h=0;for i=1:n %對PQ節(jié)點(diǎn)的處理 if i=isb&B2(i,6)=2 h=h+1; DetaS(2*h-1,1)=real(B2(i,2)-OrgS(2*h-1,1); DetaS(2*h,1)=imag(B2(i,2)-OrgS(2*h,1); endendt=0;for i=1:n %對PV節(jié)點(diǎn)的處理，注意這時不可再將h初始化為0 if i=isb&B2(i,6)=3 h=h+1; t=t+1; Det

8、aS(2*h-1,1)=real(B2(i,2)-OrgS(2*h-1,1); DetaS(2*h,1)=real(PVU(t,1)2+imag(PVU(t,1)2-real(B2(i,3)2-imag(B2(i,3)2; endendDetaS%創(chuàng)建I，用于存儲節(jié)點(diǎn)電流參數(shù)i=zeros(n-1,1);h=0;for i=1:n if i=isb h=h+1; I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1)/conj(B2(i,3); endendI%創(chuàng)建Jacbi(雅可比矩陣)Jacbi=zeros(2*n-2);h=0;k=0;for i=1:n %

9、對PQ節(jié)點(diǎn)的處理 if B2(i,6)=2 h=h+1; for j=1:n if j=isb k=k+1; if i=j %對角元素的處理 Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+real(Y(i,j)*imag(B2(i,3)+imag(I(h,1); Jacbi(2*h-1,2*k)=real(Y(i,j)*real(B2(i,3)+imag(Y(i,j)*imag(B2(i,3)+real(I(h,1); Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1); Jacbi(2*h,2*k)=J

10、acbi(2*h-1,2*k-1)-2*imag(I(h,1); else %非對角元素的處理 Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+real(Y(i,j)*imag(B2(i,3); Jacbi(2*h-1,2*k)=real(Y(i,j)*real(B2(i,3)+imag(Y(i,j)*imag(B2(i,3); Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k); Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1); end if k=(n-1) %將用于內(nèi)循環(huán)的指針置于初始值，以確保雅可比矩陣換行

11、k=0; end end end endendk=0;for i=1:n %對PV節(jié)點(diǎn)的處理 if B2(i,6)=3 h=h+1; for j=1:n if j=isb k=k+1; if i=j %對角元素的處理 Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+real(Y(i,j)*imag(B2(i,3)+imag(I(h,1); Jacbi(2*h-1,2*k)=real(Y(i,j)*real(B2(i,3)+imag(Y(i,j)*imag(B2(i,3)+real(I(h,1); Jacbi(2*h,2*k-1)=2*imag(B2(i

12、,3); Jacbi(2*h,2*k)=2*real(B2(i,3); else %非對角元素的處理 Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+real(Y(i,j)*imag(B2(i,3); Jacbi(2*h-1,2*k)=real(Y(i,j)*real(B2(i,3)+imag(Y(i,j)*imag(B2(i,3); Jacbi(2*h,2*k-1)=0; Jacbi(2*h,2*k)=0; end if k=(n-1) %將用于內(nèi)循環(huán)的指針置于初始值，以確保雅可比矩陣換行 k=0; end end end endendJacbi%求

13、解修正方程，獲取節(jié)點(diǎn)電壓的不平衡量DetaU=zeros(2*n-2,1);DetaU=inv(Jacbi)*DetaS;DetaU%修正節(jié)點(diǎn)電壓j=0;for i=1:n %對PQ節(jié)點(diǎn)處理 if B2(i,6)=2 j=j+1; B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1); endendfor i=1:n %對PV節(jié)點(diǎn)的處理 if B2(i,6)=3 j=j+1; B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1); endendB2%開始循環(huán)*while abs(max(De

14、taU)prOrgS=zeros(2*n-2,1); %!初始功率參數(shù)在迭代過程中是不累加的，所以在這里必須將其初始化為零矩陣h=0;j=0;for i=1:n if i=isb&B2(i,6)=2 h=h+1; for j=1:n OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)+imag(B2(i,3)*(real(Y(i,j)*imag(B2(j,3)+imag(Y(i,j)*real(B2(j,3); OrgS(2*h,1)=OrgS(2*h,1)+im

15、ag(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)-real(B2(i,3)*(real(Y(i,j)*imag(B2(j,3)+imag(Y(i,j)*real(B2(j,3); end endendfor i=1:n if i=isb&B2(i,6)=3 h=h+1; for j=1:n OrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)+imag(B2(i,3)*(real(Y(i,j)

16、*imag(B2(j,3)+imag(Y(i,j)*real(B2(j,3); OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)-real(B2(i,3)*(real(Y(i,j)*imag(B2(j,3)+imag(Y(i,j)*real(B2(j,3); end endendOrgS%創(chuàng)建DetaSh=0;for i=1:n if i=isb&B2(i,6)=2 h=h+1; DetaS(2*h-1,1)=real(B2(i,2)-OrgS(2*h-1,1); De

17、taS(2*h,1)=imag(B2(i,2)-OrgS(2*h,1); endendt=0;for i=1:n if i=isb&B2(i,6)=3 h=h+1; t=t+1; DetaS(2*h-1,1)=real(B2(i,2)-OrgS(2*h-1,1); DetaS(2*h,1)=real(PVU(t,1)2+imag(PVU(t,1)2-real(B2(i,3)2-imag(B2(i,3)2; endendDetaS%創(chuàng)建Ii=zeros(n-1,1);h=0;for i=1:n if i=isb h=h+1; I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*s

18、qrt(-1)/conj(B2(i,3); endendI%創(chuàng)建JacbiJacbi=zeros(2*n-2);h=0;k=0;for i=1:n if B2(i,6)=2 h=h+1; for j=1:n if j=isb k=k+1; if i=j Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+real(Y(i,j)*imag(B2(i,3)+imag(I(h,1); Jacbi(2*h-1,2*k)=real(Y(i,j)*real(B2(i,3)+imag(Y(i,j)*imag(B2(i,3)+real(I(h,1); Jacbi(2*h,

19、2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1); Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1); else Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+real(Y(i,j)*imag(B2(i,3); Jacbi(2*h-1,2*k)=real(Y(i,j)*real(B2(i,3)+imag(Y(i,j)*imag(B2(i,3); Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k); Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-

20、1); end if k=(n-1) k=0; end end end endendk=0;for i=1:n if B2(i,6)=3 h=h+1; for j=1:n if j=isb k=k+1; if i=j Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+real(Y(i,j)*imag(B2(i,3)+imag(I(h,1); Jacbi(2*h-1,2*k)=real(Y(i,j)*real(B2(i,3)+imag(Y(i,j)*imag(B2(i,3)+real(I(h,1); Jacbi(2*h,2*k-1)=2*imag(B2(i

21、,3); Jacbi(2*h,2*k)=2*real(B2(i,3); else Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+real(Y(i,j)*imag(B2(i,3); Jacbi(2*h-1,2*k)=real(Y(i,j)*real(B2(i,3)+imag(Y(i,j)*imag(B2(i,3); Jacbi(2*h,2*k-1)=0; Jacbi(2*h,2*k)=0; end if k=(n-1) k=0; end end end endendJacbiDetaU=zeros(2*n-2,1);DetaU=inv(Jacbi)*D

22、etaS;DetaU%修正節(jié)點(diǎn)電壓j=0;for i=1:n if B2(i,6)=2 j=j+1; B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1); endendfor i=1:n if B2(i,6)=3 j=j+1; B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1); endendB2Times=Times+1; %迭代次數(shù)加1endTimes一個原始數(shù)據(jù)的例子節(jié)點(diǎn)數(shù) 5支路數(shù) 5平衡節(jié)點(diǎn)編號 5精度pr 0.000001B1(支路參數(shù)矩陣)1 2 0.04+0.25i 0.

23、5i 1 0;1 3 0.1+0.35i 0 1 0;2 3 0.08+0.30i 0.5i 1 0;4 2 0.015i 0 1.05 1;5 3 0.03i 0 1.05 1B2(節(jié)點(diǎn)參數(shù)矩陣)0 -1.6-0.8i 1 0 0 2;0 -2-1i 1 0 0 2;0 -3.7-1.3i 1 0 0 2;0 5+0i 1.05 1.05 0 3;0 0 1.05 1.05 0 1X(節(jié)點(diǎn)號和對地參數(shù))1 0;2 0;3 0;4 0;5 0 程序2：%本程序采用牛頓-拉夫訊法對一個五節(jié)點(diǎn)電力系統(tǒng)進(jìn)行潮流計算% n=input(請輸入系統(tǒng)節(jié)點(diǎn)數(shù)目:n=);% nl=input(請輸入系統(tǒng)支路數(shù)

24、目:nl=);% ph=input(請輸入平衡節(jié)點(diǎn)號:ph=);% jd=input(請輸入誤差精度：jd=);% B=input(請輸入由支路參數(shù)形成的矩陣:B=);% A=input(請輸入各節(jié)點(diǎn)參數(shù)形成的矩陣:A=);clear alln=5;nl=5;ph=5;jd=1E-6;%支路參數(shù)矩陣Br=1 2 0.04+0.25*j 0.25j 1; 1 3 0.1+0.35*j 0 1; 2 3 0.08+0.30*j 0.25*j 1; 2 4 0.015*j 0 1.05; 3 5 0.03*j 0 1.05;%節(jié)點(diǎn)參數(shù)矩陣:2-PQ節(jié)點(diǎn)，3-PV節(jié)點(diǎn)，4-平衡節(jié)點(diǎn)A=-1.6-0.8

25、*j 1 0 2; -2-1*j 1 0 2; -3.7-1.3*j 1 0 2; 5 1.05 1.05 3; 0 1.05 0 4;%雅克比矩陣形成Y=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);for i=1:nl p=Br(i,1); q=Br(i,2); Y(p,q)=Y(p,q)-1./(Br(i,3)*Br(i,5); Y(q,p)=Y(p,q); Y(q,q)=Y(q,q)+1./(Br(i,3)+Br(i,4); Y(p,p)=Y(p,p)+1./(Br(i,3)*Br(i,5)2)+Br(i,4);end disp(節(jié)點(diǎn)導(dǎo)

26、納矩陣為：Y=);disp(Y);G=real(Y);B=imag(Y);for i=1:n e(i)=real(A(i,2); f(i)=imag(A(i,2); V(i)=A(i,3);endfor i=1:n S(i)=A(i);endP=real(S);Q=imag(S);%雅克比矩陣求取Ci=0;a=1;NO=2*n;N=NO-1;while a=0 a=0; for i=1:n if i=ph C(i)=0; D(i)=0; for p=1:n C(i)=C(i)+G(i,p)*e(p)-B(i,p)*f(p); D(i)=D(i)+G(i,p)*f(p)+B(i,p)*e(p);

27、 end P1=e(i)*C(i)+f(i)*D(i); Q1=f(i)*C(i)-e(i)*D(i); V2=e(i)2+f(i)2; if A(i,4)=3 DP=P(i)-P1; DQ=Q(i)-Q1; for k=1:n%非平衡節(jié)點(diǎn)時，非對角線元素 if k=ph & k=i X1=-G(i,k)*f(i)+B(i,k)*e(i); X2=-G(i,k)*e(i)-B(i,k)*f(i); X3=-X2; X4=X1; p=2*i-1; q=2*k-1; J(p,q)=X1; J(p,N)=DP; m=p+1; J(m,q)=X3; J(m,N)=DQ; q=q+1; J(p,q)=X

28、2; J(m,q)=X4; %非平衡節(jié)點(diǎn)時對角線元素 elseif k=ph & k=i X1=-D(i)+B(i,i)*e(i)-G(i,i)*f(i); X2=-C(i)-G(i,i)*e(i)-B(i,i)*f(i); X3=-C(i)+G(i,i)*e(i)+B(i,i)*f(i); X4=D(i)+B(i,i)*e(i)-G(i,i)*f(i); p=2*i-1; q=2*k-1; J(p,q)=X1; J(p,N)=DP; m=p+1; J(m,q)=X3; J(m,N)=DQ; q=q+1; J(p,q)=X2; J(m,q)=X4; end end else DP=P(i)-P

29、1; DV=V(i)2-V2; for s=1:n if s=ph & s=i X1=-G(i,s)*f(i)+B(i,s)*e(i); X2=-G(i,s)*e(i)-B(i,s)*f(i); X5=0; X6=0; p=2*i-1; q=2*s-1; J(p,q)=X1; J(p,N)=DP; m=p+1; J(m,q)=X5; J(m,N)=DV; q=q+1; J(p,q)=X2; J(m,q)=X6; elseif s=ph & s=i X1=-D(i)+B(i,i)*e(i)-G(i,i)*f(i); X2=-C(i)-G(i,s)*e(i)-B(i,s)*f(i); X5=-2*

30、f(i); X6=-2*e(i); p=2*i-1; q=2*s-1; J(p,q)=X1; J(p,N)=DP; m=p+1; J(m,q)=X5; J(m,N)=DV; q=q+1; J(p,q)=X2; J(m,q)=X6; end end end end end disp(系統(tǒng)雅克比矩陣為：J=); disp(J);%對于高壓電力系統(tǒng)，最大元素出現(xiàn)在對角線位置，簡單計算如下 DJ=J(:,1:NO-2); DP=J(1:NO-2,N); DU=DJDP; disp(*); disp( ) disp(第幾次的修正值DU？); disp(Ci); disp(*); disp(DU); for i=1:NO-2 eps=abs(DP(i); if eps=jd; a=a+1; end end Ci=Ci+1; fo