牛拉法潮流計(jì)算

1、牛拉法潮流計(jì)算%本程序的功能是用牛拉法進(jìn)行潮流計(jì)算原理介紹詳見鞠平著電氣工程%默認(rèn)數(shù)據(jù)為鞠平著電氣工程例8.4所示數(shù)據(jù)土是支路參數(shù)矩陣%第一列和第二列是節(jié)點(diǎn)編號(hào)。節(jié)點(diǎn)編號(hào)由小到大編寫%對(duì)于含有變壓器的支路，第一列為低壓側(cè)節(jié)點(diǎn)編號(hào)，第二列為高壓側(cè)節(jié)點(diǎn)編號(hào)第 三列為支路的串列阻抗參數(shù)，含變壓器支路此值為變壓器短路電抗第四列為支路的對(duì)地 導(dǎo)納參數(shù)，含變壓器支路此值不代入計(jì)算第五烈為含變壓器支路的變壓器的變比，變壓 器非標(biāo)準(zhǔn)電壓比%第六列為變壓器是否是否含有變壓器的參數(shù)，其中“ 1”為含有變壓器，“0”為不 含有變壓器2為節(jié)點(diǎn)參數(shù)矩陣%第一列為節(jié)點(diǎn)注入發(fā)電功率參數(shù)第二列為節(jié)點(diǎn)負(fù)荷功率參數(shù)第三列為節(jié)點(diǎn)電

2、壓 參數(shù) %第四列 %第五列%第六列為節(jié)點(diǎn)類型參數(shù)，“1”為平衡節(jié)點(diǎn)，“2”為PQ節(jié)點(diǎn)，“3”為PV節(jié)點(diǎn)參數(shù)%乂為節(jié)點(diǎn)號(hào)和對(duì)地參數(shù)矩陣第一列為節(jié)點(diǎn)編號(hào)第二列為節(jié)點(diǎn)對(duì)地參數(shù)clear; clc;num=input（是否采用默認(rèn)數(shù)據(jù)？（1-默認(rèn)數(shù)據(jù);2-手動(dòng)輸入）； 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.0518i 1 0;B2=0 0 1 0 0

3、 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; elsen=input（請(qǐng)輸入節(jié)點(diǎn)數(shù):n=）;n1=input（請(qǐng)輸入支路數(shù):n1=）;isb=input（請(qǐng)輸入平衡節(jié)點(diǎn)號(hào):isb=）;pr=input（請(qǐng)輸入誤差精度:pr=）;B1=input（請(qǐng)輸入支路參數(shù):B1=）;B2=input（請(qǐng)輸入節(jié)點(diǎn)參數(shù):B2=）;X=input（節(jié)點(diǎn)號(hào)和對(duì)地參數(shù):X=）; endTimes=1; %迭代次數(shù)%創(chuàng)建節(jié)點(diǎn)導(dǎo)納矩陣Y=zeros(n); for i=1:n1if B1(i,

4、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)-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(

5、i,3) + (B1(i,5)*(B1(i,5)-1)/B1(i,3); end endfor i=1:n1Y(i,i)=Y(i,i)+X(i,2); %計(jì)及補(bǔ)償電容電納 enddisp(導(dǎo)納矩陣為：)；disp(Y);%顯示導(dǎo)納矩陣%初始化 OrgS、DetaS OrgS=zeros(2*n-2,1); DetaS=zeros(2*n-2,1);%創(chuàng)建OrgS，用于存儲(chǔ)初始功率參數(shù)h=0; j=0;for i=1:n %對(duì)？。節(jié)點(diǎn)的處理if i=isb&B2(i,6)=2%不是平衡點(diǎn)&是 PQ 點(diǎn)h=h+1;forj=1:n%公式8-74%Pi=ei*(Gij*ej-Bij*fj)+fi*(

6、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)=0rgS(2*h,1)+imag(B2(i,3)*(real(Y(i,j)*real(B2(j,3)- imag(Y(i,j)*imag(B2(j,3

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

8、,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(i,j)*imag(B2(j,3)+imag(Y(i,j)*real(B2(j,3); end end end%創(chuàng)建PVU用于存儲(chǔ)PV節(jié)點(diǎn)的初始電壓PVU=zer

9、os(n-h-1,1); t=0; for i=1:nif B2(i,6)=3t=t+1;PVU(t,1)=B2(i,3);end end%創(chuàng)建DetaS，用于存儲(chǔ)有功功率、無(wú)功功率和電壓幅值的不平衡量h=0;for i=1:n %對(duì) PQ 節(jié)點(diǎn)的處理 if i=isb&B2(i,6)=2h=h+1;DetaS(2*h-1,1)=real(B2(i,2)-OrgS(2*h-1,1); TlPiDetaS(2*h,1)=imag(B2(i,2)-OrgS(2*h,1); TlQi end end t=0;for i=1:n %對(duì)PV節(jié)點(diǎn)的處理，注意這時(shí)不可再將h初始化為0ifi=isb&B2(i

10、,6)=3h=h+1;t=t+1;DetaS(2*h-1,1)=real(B2(i,1)-OrgS(2*h-1,1); TlPiDetaS(2*h,1)=real(PVU(t,1)廣2+imag(PVU(t,1)2-real(B2(i,3)2-imag(B2(i,3)*2; TlUi end end % DetaS%創(chuàng)建I,用于存儲(chǔ)節(jié)點(diǎn)電流參數(shù)i=zeros(n-1,1); h=0; for i=1:n if i=isbh=h+1;I(h,1) = (OrgS(2*hT,1)-OrgS(2*h,1)*sqrt(T)/conj(B2(i,3);%conj 求共軛 end end%創(chuàng)建 Jacbi

11、(雅可比矩陣)Jacbi=zeros(2*n-2); h=0; k=0;for i=1:n %對(duì) PQ 節(jié)點(diǎn)的處理 if B2(i,6)=2h=h+1;forj=1:nif j=isbk=k+1;if i=j %對(duì)角元素的處理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(

12、2*h-1,2*k)+2*real(I(h,1);Jacbi(2*h,2*k)=Jacbi(2*hT,2*kT)-2*imag(I(h,1);else %非對(duì)角元素的處理Jacbi(2*hT,2*kT)=-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);endif k=(n-

13、1) %將用于內(nèi)循環(huán)的指針置于初始值，以確保雅可比矩陣換行 k=0;endendend end end k=0;for i=1:n %對(duì)PV節(jié)點(diǎn)的處理 j=1:nif j=isbif B2(i,6)=3 k=k+1;for i=1:n %對(duì)PV節(jié)點(diǎn)的處理 j=1:nif j=isbif B2(i,6)=3 k=k+1;h=h+1;forif i=j %對(duì)角元素的處理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(

14、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 %非對(duì)角元素的處理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;endif k=(n-1) %將用于內(nèi)循環(huán)的指針置于初始值，以確保雅可比矩陣換行 k=0;endendend end enddisp(初始雅可比矩陣為：)；disp(Jacbi);%求解修正方程，獲取節(jié)點(diǎn)電壓的不平衡量DetaU=zeros(2*n-2,1);DetaU=inv(Jacbi)*DetaS; %in v 矩陣求逆 % DetaU%修正節(jié)點(diǎn)電壓j=0;for i=1:n %MPQ 節(jié)點(diǎn)處理if B2(i,6)=2j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);end endfor i=1:n %對(duì)PV 節(jié)點(diǎn)的處理