牛拉法潮流計算(共10頁)

1、精選優(yōu)質(zhì)文檔-傾情為你奉上%本程序的功能是用牛拉法進行潮流計算%原理介紹詳見鞠平著電氣工程%默認數(shù)據(jù)為鞠平著電氣工程例8.4所示數(shù)據(jù)%B1是支路參數(shù)矩陣%第一列和第二列是節(jié)點編號。節(jié)點編號由小到大編寫%對于含有變壓器的支路，第一列為低壓側節(jié)點編號，第二列為高壓側節(jié)點編號%第三列為支路的串列阻抗參數(shù)，含變壓器支路此值為變壓器短路電抗%第四列為支路的對地導納參數(shù)，含變壓器支路此值不代入計算%第五烈為含變壓器支路的變壓器的變比，變壓器非標準電壓比%第六列為變壓器是否是否含有變壓器的參數(shù)，其中“1”為含有變壓器，“0”為不含有變壓器%B2為節(jié)點參數(shù)矩陣%第一列為節(jié)點注入發(fā)電功率參數(shù)%第二列為節(jié)點負荷功

2、率參數(shù)%第三列為節(jié)點電壓參數(shù)%第四列%第五列%第六列為節(jié)點類型參數(shù)，“1”為平衡節(jié)點，“2”為PQ節(jié)點，“3”為PV節(jié)點參數(shù)%X為節(jié)點號和對地參數(shù)矩陣%第一列為節(jié)點編號%第二列為節(jié)點對地參數(shù)clear;clc;num=input('是否采用默認數(shù)據(jù)？(1-默認數(shù)據(jù);2-手動輸入)');if num=1 n=4; n1=4; isb=4; pr=0.00001; B1=1 2 0.1667i 0 0.8864 1;1 3 0.1302+0.2479i 0.0258i 1 0;1 4 0.1736+0.3306i 0.0344i 1 0;3 4 0.2603+0.4959i 0.0

3、518i 1 0; B2=0 0 1 0 0 2;0 -0.5-0.3i 1 0 0 2;0.2 0 1.05 0 0 3;0 -0.15-0.1i 1.05 0 0 1; X=1 0;2 0.05i;3 0;4 0;else n=input('請輸入節(jié)點數(shù):n='); n1=input('請輸入支路數(shù):n1='); isb=input('請輸入平衡節(jié)點號:isb='); pr=input('請輸入誤差精度:pr='); B1=input('請輸入支路參數(shù):B1='); B2=input('請輸入節(jié)點參數(shù):

4、B2='); X=input('節(jié)點號和對地參數(shù):X=');endTimes=1; %迭代次數(shù)%創(chuàng)建節(jié)點導納矩陣Y=zeros(n);for i=1:n1 if B1(i,6)=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)-

5、B1(i,5)/B1(i,3); Y(q,p)=Y(p,q); Y(p,p)=Y(p,p)+B1(i,5)/B1(i,3)+(1-B1(i,5)/B1(i,3); Y(q,q)=Y(q,q)+B1(i,5)/B1(i,3)+(B1(i,5)*(B1(i,5)-1)/B1(i,3); endendfor i=1:n1 Y(i,i)=Y(i,i)+X(i,2); %計及補償電容電納enddisp('導納矩陣為：');disp(Y); %顯示導納矩陣 %初始化OrgS、DetaSOrgS=zeros(2*n-2,1);DetaS=zeros(2*n-2,1); %創(chuàng)建OrgS，用于存

6、儲初始功率參數(shù)h=0;j=0;for i=1:n %對PQ節(jié)點的處理 if i=isb&B2(i,6)=2 %不是平衡點&是PQ點 h=h+1; for j=1:n %公式8-74 %Pi=ei*(Gij*ej-Bij*fj)+fi*(Gij*fj+Bij*ej) %Qi=fi*(Gij*ej-Bij*fj)-ei*(Gij*fj+Bij*ej) 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)

7、*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 endendfor i=1:n %對PV節(jié)點的處理，注意這時不可再將h初始化為0 if i=isb&B2(i,6)=3 %不是平衡點&是PV點 h=h+1; for j=1:n %公式8-

8、75-a %Pi=ei*(Gij*ej-Bij*fj)+fi*(Gij*fj+Bij*ej) %Qi=fi*(Gij*ej-Bij*fj)-ei*(Gij*fj+Bij*ej) 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)+imag(B2(i,3)*(real(Y(i,j)*real

9、(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 endend%創(chuàng)建PVU 用于存儲PV節(jié)點的初始電壓PVU=zeros(n-h-1,1);t=0;for i=1:n if B2(i,6)=3 t=t+1; PVU(t,1)=B2(i,3); endend%創(chuàng)建DetaS，用于存儲有功功率、無功功率和電壓幅值的不平衡量h=0;for i=1:n %對PQ節(jié)點的處理 if i=isb&B2(i,6)=2 h=h+1; DetaS(2*

10、h-1,1)=real(B2(i,2)-OrgS(2*h-1,1); %delPi DetaS(2*h,1)=imag(B2(i,2)-OrgS(2*h,1); %delQi endendt=0;for i=1:n %對PV節(jié)點的處理，注意這時不可再將h初始化為0 if i=isb&B2(i,6)=3 h=h+1; t=t+1; DetaS(2*h-1,1)=real(B2(i,1)-OrgS(2*h-1,1); %delPi DetaS(2*h,1)=real(PVU(t,1)2+imag(PVU(t,1)2-real(B2(i,3)2-imag(B2(i,3)2; %delUi e

11、ndend% DetaS%創(chuàng)建I，用于存儲節(jié)點電流參數(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);%conj求共軛 endend%創(chuàng)建Jacbi(雅可比矩陣)Jacbi=zeros(2*n-2);h=0;k=0;for i=1:n %對PQ節(jié)點的處理 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)*r

12、eal(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)=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)

13、*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)的指針置于初始值，以確保雅可比矩陣換行 k=0; end end end endendk=0;for i=1:n %對PV節(jié)點的處理 if B2(i,6)=3 h=h+1; for j=1:n if j=isb k=k+1; if i=j %對角元

14、素的處理 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,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

15、,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 endenddisp('初始雅可比矩陣為：');disp(Jacbi); %求解修正方程，獲取節(jié)點電壓的不平衡量DetaU=zeros(2*n-2,1);DetaU=inv(Jacbi)*DetaS; %inv矩陣求逆

16、% DetaU%修正節(jié)點電壓j=0;for i=1:n %對PQ節(jié)點處理 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é)點的處理 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); endend% B2%開始循環(huán)*while abs(max(DetaU)>pr OrgS=zeros(2*n-2,1); h=0; j=0; for i=1:n if i=

17、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)+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

18、(i,j)*imag(B2(j,3)+imag(Y(i,j)*real(B2(j,3); end end end for 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)*imag(B2(j,3)+imag(Y(i,j)*real(B2(j,3); OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3)

19、*(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 end end% OrgS %創(chuàng)建DetaS h=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); DetaS(2*h,1)=imag(B2(i,2)-OrgS(2*h,1); end end t=0; for i=1:n if i=i

20、sb&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,1)=real(B2(i,1)-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; end end% DetaS %創(chuàng)建I 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)/

21、conj(B2(i,3); end end% I %創(chuàng)建Jacbi Jacbi=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,

22、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-

23、1); end if k=(n-1) k=0; end end end end end k=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(B

24、2(i,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 end end% Jacbi DetaU=zeros(2*n-2,1); DetaU=inv(

25、Jacbi)*DetaS;% DetaU %修正節(jié)點電壓 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); end end for 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); end end% B2 Times=Times+1; %迭代次數(shù)加1enddisp('迭代次數(shù)為：');disp(Times);disp('收斂時電壓修正量

26、為：');disp(DetaU);for k=1:n E(k)=B2(k,3); e(k)=real(E(k); f(k)=imag(E(k); V(k)=sqrt(e(k)2+f(k)2); sida(k)=atan(f(k)./e(k)*180./pi;end%= 計算各輸出量 =disp('各節(jié)點的實際電壓標幺值E為(節(jié)點號從小到大排列)：');disp(E); %顯示各節(jié)點的實際電壓標幺值E用復數(shù)表示disp('-');disp('各節(jié)點的電壓大小V為(節(jié)點號從小到大排列)：');disp(V); %顯示各節(jié)點的電壓大小V的模值d

27、isp('-');disp('各節(jié)點的電壓相角sida為(節(jié)點號從小到大排列)：');disp(sida); %顯示各節(jié)點的電壓相角for p=1:n C(p)=0; for q=1:n C(p)=C(p)+conj(Y(p,q)*conj(E(q); %計算各節(jié)點的注入電流的共軛值 end S(p)=E(p)*C(p); %計算各節(jié)點的功率 S = 電壓 X 注入電流的共軛值enddisp('各節(jié)點的功率S為(節(jié)點號從小到大排列)：');disp(S); %顯示各節(jié)點的注入功率Sline=zeros(n1,5);disp('-'

28、);disp('各條支路的首端功率Si為(順序同您輸入B1時一致)：');for i=1:n1 p=B1(i,1); q=B1(i,2); Sline(i,1)=B1(i,1); Sline(i,2)=B1(i,2); if B1(i,6)=0 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); Siz(i)=Si(p,q); else Si(p,q)=E(p)*(conj(E(p)*(1-B1(i,5)/B1(i,3)+(conj(E(p)-conj(E(q)*(B1(i,5)/B1(i,3); Siz