牛拉潮流程序

1、clear;clc;%+輸+入帶有變壓器的支路矩陣中各節(jié)點(diǎn)對應(yīng)各變比 % function%打開矩陣 文件%= NODE,Branch=OpDF_;Node=NODE;%節(jié)點(diǎn)號%節(jié)點(diǎn)類型%各支路信息保存在BR中%支路變壓器變比， 0代表沒有變壓器%節(jié)點(diǎn)數(shù)%支路數(shù)%取節(jié)點(diǎn)矩陣的行和列%bu矩陣的行數(shù)即節(jié)點(diǎn)數(shù)%取%支路矩陣的行和列%支路 branch 矩陣行數(shù)即支路數(shù)N=Node(:,1);Type=Node(:,2);BR=BranchK=Branch(:,6); n=length(N); nbr=length(K);Total_of_Bus1=size(NODE);Total_of_Bus=T

2、otal_of_Bus1(1,1) Total_of_Branch1=size(Branch);Total_of_Branch=Total_of_Branch1(1,1);%將節(jié)點(diǎn)排序重新存儲節(jié)點(diǎn)信息%定%義為 節(jié)點(diǎn)數(shù)的方陣Z=zeros(Total_of_Bus1);format short%排序標(biāo)志位%P（節(jié)點(diǎn)標(biāo)志位%P節(jié)點(diǎn)標(biāo)志位%平衡節(jié)點(diǎn)標(biāo)志位按照PQ, PV,平衡節(jié)點(diǎn)的次序排序各種節(jié)點(diǎn)-統(tǒng)計(jì)PQ節(jié)點(diǎn)數(shù)0代表是pq節(jié)點(diǎn)b=1;pq=0;pv=0;ph=0;%for a=1:Total_of_Bus if NODE(a,2) = 0Z(b,:)=NODE(a,:);b=b+1;pq=pq+

3、1;endend統(tǒng)計(jì)PV節(jié)點(diǎn)數(shù)2代表pv節(jié)點(diǎn)統(tǒng)計(jì)平衡節(jié)點(diǎn)數(shù) 3 代表平衡節(jié)點(diǎn)%for a=1:Total_of_Bus if NODE(a,2) = 2Z(b,:)=NODE(a,:);b=b+1;pv=pv+1;endend%for a=1:Total_of_Bus if NODE(a,2) = 3 Z(b,:)=NODE(a,:); b=b+1; ph=ph+1; end end ZZ2=Z;%將%節(jié)%點(diǎn)%進(jìn)%行%重%新%排%序% % mm=zeros(n,1);for i=1:nmm(i,1)=i;endZ1(:,1)=mm(:,1);Branch1=zeros(nbr,2);for i

4、=1:nif Z(i,1)=Z1(i)for j=1:nbrif Branch(j,1)=Z(i,1)Branch1(j,1)=Z1(i);endif Branch(j,2)=Z(i,1)Branch1(j,2)=Z1(i);endendelsefor j=1:nbrif Branch(j,1)=Z(i,1)Branch1(j,1)=Z(i,1);endif Branch(j,2)=Z(i,1)Branch1(j,2)=Z(i,1);endendendendBranch(:,1)=Branch1(:,1);Branch(:,2)=Branch1(:,2);Z(:,1)=Z1(:,1);j=sq

5、rt(-1);% 矩陣已經(jīng)完成按照 PQ PV,平衡節(jié)點(diǎn)的順序排列起來YSNODE=Z; %保存排序后的原始節(jié)點(diǎn)數(shù)據(jù) %=Y=zeros(n,n);%求%互%導(dǎo)%納% for i=1:nfor t=1:nbr%非變壓器支路 %if (Branch(t,1)=i|Branch(t,2)=i)&& Branch(t,6)=0 Y(Branch(t,1),Branch(t,2)=-1/(Branch(t,3)+j*Branch(t,4); Y(Branch(t,2),Branch(t,1)=Y(Branch(t,1),Branch(t,2);%變壓器支路 %else if (Bran

6、ch(t,1)=i|Branch(t,2)=i)&&Branch(t,6)=0 Y(Branch(t,1),Branch(t,2)=(-1/(j*Branch(t,4)/Branch(t,6); Y(Branch(t,2),Branch(t,1)=Y(Branch(t,1),Branch(t,2);endendendend%求%自%導(dǎo)%納% % for i=1:nfor t=1:nbr%非變壓器支路 %if (Branch(t,1)=i|Branch(t,2)=i)&& Branch(t,6)=0 Y(i,i)=Y(i,i)+1/(Branch(t,3)+j*B

7、ranch(t,4)+(1/2)*j*Branch(t,5);else if Branch(t,1)=i&&Branch(t,6)=0%變%壓% 器支路為首節(jié)點(diǎn) %Y(i,i)=Y(i,i)+1/(j*Branch(t,4);%變%壓%器_支路且 ielse if Branch(t,2)=i&&Branch(t,6)=0為末節(jié)點(diǎn) % Y(i,i)=Y(i,i)+(1/(j*Branch(t,4)/(Branch(t,6)*Branch(t,6);endendendendend%若%有%并聯(lián)電容器組，則自導(dǎo)納要加上并聯(lián)電容器的導(dǎo)納i=1:nif NODE(i,13

8、)=0Y(i,i)=Y(i,i)+j*NODE(i,13)endfor end%實(shí)部，即電導(dǎo)%虛部，即電納Y n=length(N); G=real(Y); B=imag(Y);%給%定%初%始%的%電%壓%值%與相位值%U_first=Z(:,3);%初始電壓幅值phase_first=Z(:,4);%初始相位值e=U_first.*cos(phase_first);f=U_first.*sin(phase_first);%節(jié)點(diǎn)負(fù)荷有功分量 %節(jié)點(diǎn)負(fù)荷無功分量 %發(fā)電機(jī)發(fā)出的有功 %發(fā)電機(jī)發(fā)出的無功 %節(jié)點(diǎn)電壓都的初始值%算%l%P 初始功率量 % P=Z(:,5);Q=Z(:,6);PG=

9、Z(:,7);QG=Z(:,8);U0=Z(:,9);Delta_P=zeros(1,n-1);for i=1:n-1for j=1:nDelta_P(i)=Delta_P(i)-e(i)*(G(i,j)*e(j)-B(i,j)*f(j)+f(j)*(G(i,j)*f(j)+B(i,j)*e(j); endendfor i=1:n-1Delta_P(i)=Delta_P(i)-(P(i)-PG(i);endDelta_P%«%)初始功率量 % m=0;%十算PV節(jié)點(diǎn)的個(gè)數(shù) %n代表pv節(jié)點(diǎn)個(gè)數(shù)for i=1:n;if Type(i)=2;m=m+1;end end Delta_Q=z

10、eros(1,n-m-1); for i=1:n-m-1for j=1:nDelta_Q(i)=Delta_Q(i)-f(i)*(G(i,j)*e(j)-B(i,j)*f(j)+e(i)*(G(i,j)*f(j)+B(i,j)*e(j); end end for i=1:n-m-1Delta_Q(i)=Delta_Q(i)-(Q(i)-QG(i); end Delta_Q Delta_V=zeros(1,m); for i=1:mfor j=1:nif Type(j)=2Delta_V(i)=U0a)A2-(eQ)A2+f(i)A2);end end endDelta_V num=0;disp

11、( ' 第' ,num2str(num), ' 次時(shí)的 Delta 總的失配量為： % 進(jìn)入循環(huán)體判斷是否滿足條件' )% 先算出最大值，作為判斷是否收斂的依據(jù)DEL=Delta_P Delta_Q;%Delta_P Delta_Q_MAX =max(abs(DEL);MAXTheta_first=zeros(1,n);U_f=U_first'Delta_F_E1=Theta_first(1:n-1) U_f(1:n-m-1); Delta_F=Delta_F_E1'Delta_Cor=Delta_F_E1;Delta_u%disp( '

12、 第一次最大失配量誤差：%if MAX>1e-004 disp(Delta_the%' ,num2str(MAX)循環(huán)判斷 % 判斷依據(jù)面開始下一次迭代過程！endwhile MAX>1e-004num=num+1;%形%成%雅%克%比%矩%陣 %先求非對角元素 -(H)%Hik=zeros(n-1,n-1);for i=1:n-1for k=1:n-1 if i=k theik=Theta_first(i)-Theta_first(k); Hik(i,k)=-U_first(i)*U_first(k)*(G(i,k)*sin(theik)-B(i,k)*cos(theik

13、);endendend再求對角元素H)%for i=1:n-1for k=1:nif i=k theik=Theta_first(i)-Theta_first(k); Hik(i,i)=Hik(i,i)+U_first(k)*(G(i,k)*sin(theik)-B(i,k)*cos(theik);endend Hik(i,i)=U_first(i)*Hik(i,i);end先求非對角元素 -NHik %Nik=zeros(n-1,n-m-1);for i=1:n-1for k=1:n-m-1 if i=k theik=Theta_first(i)-Theta_first(k); Nik(i,

14、k)=-U_first(i)*U_first(k)*(G(i,k)*cos(theik)+B(i,k)*sin(theik);endendend再求對角元素%for i=1:n-m-1for k=1:nif i=k theik=Theta_first(i)-Theta_first(k); Nik(i,i)=Nik(i,i)+U_first(k)*(G(i,k)*cos(theik)+B(i,k)*sin(theik);endend Nik(i,i)=-U_first(i)*Nik(i,i)-2*U_first(i)*U_first(i)*G(i,i);end先求非對角元素M)Nik %Mik=

15、zeros(n-m-1,n-1);for i=1:n-m-1for k=1:n-1 if i=k theik=Theta_first(i)-Theta_first(k);Mik(i,k)=U_first(i)*U_first(k)*(G(i,k)*cos(theik)+B(i,k)*sin(theik);endendend再求對角元素%for i=1:n-m-1for k=1:nif i=k theik=Theta_first(i)-Theta_first(k); Mik(i,i)=Mik(i,i)+U_first(k)*(G(i,k)*cos(theik)+B(i,k)*sin(theik)

16、;endend Mik(i,i)=-U_first(i)*Mik(i,i);endMikL)% 先求非對角元素Lik=zeros(n-m-1,n-m-1);for i=1:n-m-1for k=1:n-m-1 if i=k theik=Theta_first(i)-Theta_first(k);Lik(i,k)=-U_first(i)*U_first(k)*(G(i,k)*sin(theik)-B(i,k)*cos(theik);endendend再求對角元素%for i=1:n-m-1for k=1:nif i=k theik=Theta_first(i)-Theta_first(k); L

17、ik(i,i)=Lik(i,i)+U_first(k)*(G(i,k)*sin(theik)-B(i,k)*cos(theik);endend Lik(i,i)=-U_first(i)*Lik(i,i)+2*U_first(i)*U_first(i)*B(i,i);endLik% 至此雅可比矩陣已經(jīng)形成 % 開始構(gòu)造 Delta_f ； Delta_ekacb=Hik Nik;Mik Lik; kacb%雅%克%比%矩%陣 %-修正各個(gè)量,包括 e,f,P,Q,UA2(重要！ )-%DEL=DEL'Delta_F_E=(-1*inv(kacb)*DEL;Delta_F=Delta_F_

18、E' Delta_Cor=Delta_F+Delta_Cor;Theta_first(1,1:n-1)=Delta_Cor(1,1:n-1);Theta_first(1,n)=0;%初%始%相%角%的%修%正%Theta_first=Theta_first'%修%正%后的角度值 %Theta_firstDelta_C=Delta_Cor'U_first(1:n-m-1,1)=Delta_C(n:2*n-m-2,1);U_first%修%正%后%的%電壓值 %e=U_first.*cos(Theta_first);f=U_first.*sin(Theta_first);計(jì)

19、算修正Delta_Pe f %- Delta_P=zeros(1,n-1); for i=1:n-1 for k=1:nDelta_P(i)=Delta_P(i)-e(i,1)*(G(i,k)*e(k,1)-B(i,k)*f(k,1)-f(i,1)*(G(i,k)*f(k,1)+B(i,k )*e(k,1);endendfor i=1:n-1Delta_P(i)=Delta_P(i)-(P(i,1)-PG(i,1);endDelta_P- 計(jì)算完成 % 計(jì)算 Delta_Q%Delta_P%Delta_Q=zeros(1,n-m-1);for i=1:n-m-1for k=1:nDelta_Q

20、(i)=Delta_Q(i)-f(i)*(G(i,k)*e(k)-B(i,k)*f(k)+e(i)*(G(i,k)*f(k)+B(i,k)*e(k); endendfor i=1:n-m-1Delta_Q(i)=Delta_Q(i)-(Q(i)-QG(i);endDelta_QDEL=Delta_P Delta_Q;disp( ' 第 ' ,num2str(num), ' 次時(shí)的 Delta 總的失配量為： ' )% DEL% 繼續(xù)判斷最大值MAX =max(abs(DEL);Theta_first=Theta_first'end%求%平%衡%節(jié)%點(diǎn)%的有功功率和無功功率%Ps0=0;i=n;for t=1:ntheij=Theta_first(i)-Theta_first(t);Ps0=Ps0+U_first(t)*(G(i,t)*cos(theij)+B(i,t)*sin(theij);endPs0=U_first(i)*Ps0;Z(i,7)=Ps0;Qs0=0;i=n;for t=1:ntheij=Theta_first(i)-Theta_first(t);Qs0=Qs0+U_first(t)*(G(i,t)*sin(theij)-B(i,t)*cos(thei