%本程序的功能是用牛拉法進(jìn)行潮流計(jì)算%原理介紹詳見鞠平著電氣工程%默認(rèn)數(shù)據(jù)為鞠平著電氣工程例8.4所示數(shù)據(jù)%B1是支路參數(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”為不含有變壓器%B2為節(jié)點(diǎn)參數(shù)矩陣%第一列為節(jié)點(diǎn)注入發(fā)電功率參數(shù)%第二列為節(jié)點(diǎn)負(fù)荷功率參數(shù)%第三列為節(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)號(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 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(請(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: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)-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); %計(jì)及補(bǔ)償電容電納enddisp(導(dǎo)納矩陣為：);disp(Y); %顯示導(dǎo)納矩陣 %初始化OrgS、DetaSOrgS=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ì)PQ節(jié)點(diǎn)的處理 if i=isb&B2(i,6)=2 %不是平衡點(diǎn)&是PQ點(diǎn) 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)*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 %對(duì)PV節(jié)點(diǎn)的處理，注意這時(shí)不可再將h初始化為0 if i=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,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 endend%創(chuàng)建PVU 用于存儲(chǔ)PV節(jié)點(diǎn)的初始電壓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，用于存儲(chǔ)有功功率、無功功率和電壓幅值的不平衡量h=0;for i=1:n %對(duì)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); %delPi DetaS(2*h,1)=imag(B2(i,2)-OrgS(2*h,1); %delQi endendt=0;for i=1:n %對(duì)PV節(jié)點(diǎn)的處理，注意這時(shí)不可再將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 endend% DetaS%創(chuàng)建I，用于存儲(chǔ)節(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);%conj求共軛 endend%創(chuàng)建Jacbi(雅可比矩陣)Jacbi=zeros(2*n-2);h=0;k=0;for i=1:n %對(duì)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 %對(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(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 %非對(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)=-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 %對(duì)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 %對(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)=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; end if k=(n-1) %將用于內(nèi)循環(huán)的指針置于初始值，以確保雅可比矩陣換行 k=0; end end end endenddisp(初始雅可比矩陣為：);disp(Jacbi); %求解修正方程，獲取節(jié)點(diǎn)電壓的不平衡量DetaU=zeros(2*n-2,1);DetaU=inv(Jacbi)*DetaS; %inv矩陣求逆% DetaU%修正節(jié)點(diǎn)電壓j=0;for i=1:n %對(duì)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 %對(duì)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); 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=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(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)*(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=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,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)/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,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-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(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,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(Jacbi)*DetaS;% 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); 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(收斂時(shí)電壓修正量為：);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%= 計(jì)算各輸出量 =disp(各節(jié)點(diǎn)的實(shí)際電壓標(biāo)幺值E為(節(jié)點(diǎn)號(hào)從小到大排列)：);disp(E); %顯示各節(jié)點(diǎn)的實(shí)際電壓標(biāo)幺值E用復(fù)數(shù)表示disp(-);disp(各節(jié)點(diǎn)的電壓大小V為(節(jié)點(diǎn)號(hào)從小到大排列)：);disp(V); %顯示各節(jié)點(diǎn)的電壓大小V的模值disp(-);disp(各節(jié)點(diǎn)的電壓相角sida為(節(jié)點(diǎn)號(hào)從小到大排列)：);disp(sida); %顯示各節(jié)點(diǎn)的電壓相角for p=1:n C(p)=0; for q=1:n C(p)=C(p)+conj(Y(p,q)*conj(E(q); %計(jì)算各節(jié)點(diǎn)的注入電流的共軛值 end S(p)=E(p)*C(p); %計(jì)算各節(jié)點(diǎn)的功率 S = 電壓 X 注入電流的共軛值enddisp(各節(jié)點(diǎn)的功率S為(節(jié)點(diǎn)號(hào)從小到大排列)：);disp(S); %顯示各節(jié)點(diǎn)的注入功率Sline=zeros(n1,5);disp(-);disp(各條支路的首端功率Si為(順序同您輸入B1時(shí)一致)：);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