牛拉潮流程序

1、clear;clc;%+输+入带有变压器的支路矩阵中各节点对应各变比 % function%打开矩阵 文件%= NODE,Branch=OpDF_;Node=NODE;%节点号%节点类型%各支路信息保存在BR中%支路变压器变比， 0代表没有变压器%节点数%支路数%取节点矩阵的行和列%bu矩阵的行数即节点数%取%支路矩阵的行和列%支路 branch 矩阵行数即支路数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);%将节点排序重新存储节点信息%定%义为 节点数的方阵Z=zeros(Total_of_Bus1);format short%排序标志位%P（节点标志位%P节点标志位%平衡节点标志位按照PQ, PV,平衡节点的次序排序各种节点-统计PQ节点数0代表是pq节点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统计PV节点数2代表pv节点统计平衡节点数 3 代表平衡节点%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;%将%节%点%进%行%重%新%排%序% % 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);% 矩阵已经完成按照 PQ PV,平衡节点的顺序排列起来YSNODE=Z; %保存排序后的原始节点数据 %=Y=zeros(n,n);%求%互%导%纳% 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%求%自%导%纳% % 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%变%压% 器支路为首节点 %Y(i,i)=Y(i,i)+1/(j*Branch(t,4);%变%压%器_支路且 ielse if Branch(t,2)=i&&Branch(t,6)=0为末节点 % Y(i,i)=Y(i,i)+(1/(j*Branch(t,4)/(Branch(t,6)*Branch(t,6);endendendendend%若%有%并联电容器组，则自导纳要加上并联电容器的导纳i=1:nif NODE(i,13

8、)=0Y(i,i)=Y(i,i)+j*NODE(i,13)endfor end%实部，即电导%虚部，即电纳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);%节点负荷有功分量 %节点负荷无功分量 %发电机发出的有功 %发电机发出的无功 %节点电压都的初始值%算%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节点的个数 %n代表pv节点个数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), ' 次时的 Delta 总的失配量为： % 进入循环体判断是否满足条件' )% 先算出最大值，作为判断是否收敛的依据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)循环判断 % 判断依据面开始下一次迭代过程！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% 至此雅可比矩阵已经形成 % 开始构造 Delta_f ； Delta_ekacb=Hik Nik;Mik Lik; kacb%雅%克%比%矩%阵 %-修正各个量,包括 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);计

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- 计算完成 % 计算 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), ' 次时的 Delta 总的失配量为： ' )% DEL% 继续判断最大值MAX =max(abs(DEL);Theta_first=Theta_first'end%求%平%衡%节%点%的有功功率和无功功率%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