牛拉潮流程序.docx

编号：时间：2021年x月x日书山有路勤为径，学海无涯苦作舟页码：第9页 共9页clear;clc;%+输入带有变压器的支路矩阵中各节点对应各变比% function%=NODE,Branch=OpDF_; %打开矩阵 (test.m)文件Node=NODE;N=Node(:,1); %节点号Type=Node(:,2); %节点类型BR=Branch %将支路信息保存在BR中%K=Branch(:,6); %支路变压器变比，0代表没有变压器n=length(N); %节点数nbr=length(K); %支路数Total_of_Bus1=size(NODE);%取节点矩阵的行和列Total_of_Bus=Total_of_Bus1(1,1)%bus矩阵的行数即 节点数Total_of_Branch1=size(Branch);%取支路矩阵的行和列Total_of_Branch=Total_of_Branch1(1,1);% 支路branch矩阵行数即 支路数Z=zeros(Total_of_Bus1); %将节点排序重新存储节点信息%定义为 节点数的方阵format shortb=1; %排序标志位pq=0; %PQ节点标志位pv=0; %PV节点标志位ph=0; %平衡节点标志位%-按照PQ,PV,平衡节点的次序排序各种节点%-统计PQ节点数 0代表是pq节点for a=1:Total_of_Bus if NODE(a,2) = 0 Z(b,:)=NODE(a,:); b=b+1; pq=pq+1; endend%-统计PV节点数2代表pv节点for a=1:Total_of_Bus if NODE(a,2) = 2 Z(b,:)=NODE(a,:); b=b+1; pv=pv+1; endend%-统计平衡节点数 3代表平衡节点for a=1:Total_of_Bus if NODE(a,2) = 3 Z(b,:)=NODE(a,:); b=b+1; ph=ph+1; endendZZ2=Z;%将节点进行重新排序%mm=zeros(n,1);for i=1:n mm(i,1)=i;endZ1(:,1)=mm(:,1);Branch1=zeros(nbr,2);for i=1:n if Z(i,1)=Z1(i) for j=1:nbr if Branch(j,1)=Z(i,1) Branch1(j,1)=Z1(i); end if Branch(j,2)=Z(i,1) Branch1(j,2)=Z1(i); end end else for j=1:nbr if Branch(j,1)=Z(i,1) Branch1(j,1)=Z(i,1); end if Branch(j,2)=Z(i,1) Branch1(j,2)=Z(i,1); end end endendBranch(:,1)=Branch1(:,1);Branch(:,2)=Branch1(:,2);Z(:,1)=Z1(:,1); j=sqrt(-1); %-矩阵已经完成按照PQ，PV，平衡节点的顺序排列起来-YSNODE=Z; %保存排序后的原始节点数据%= Y=zeros(n,n);%求互导纳%for i=1:n for 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 (Branch(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); end end endend%求自导纳% for i=1:n for 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*Branch(t,4)+(1/2)*j*Branch(t,5); else if Branch(t,1)=i&Branch(t,6)=0 %变压器支路且i为首节点% Y(i,i)=Y(i,i)+1/(j*Branch(t,4); else if Branch(t,2)=i&Branch(t,6)=0%_变压器支路且i为末节点% Y(i,i)=Y(i,i)+(1/(j*Branch(t,4)/(Branch(t,6)*Branch(t,6); end end end endend %若有并联电容器组，则自导纳要加上并联电容器的导纳%for i=1:n if NODE(i,13)=0 Y(i,i)=Y(i,i)+j*NODE(i,13) endendY 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);%计算Delta_P初始功率量%P=Z(:,5); %节点负荷有功分量Q=Z(:,6); %节点负荷无功分量PG=Z(:,7); %发电机发出的有功QG=Z(:,8); %发电机发出的无功U0=Z(:,9); %节点电压都的初始值Delta_P=zeros(1,n-1);for i=1:n-1 for j=1:n Delta_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-1 Delta_P(i)=Delta_P(i)-(P(i)-PG(i);endDelta_P%计算Delta_Q初始功率量%m=0;for i=1:n; if Type(i)=2; %计算PV节点的个数 m=m+1; %m代表pv节点个数 endendDelta_Q=zeros(1,n-m-1);for i=1:n-m-1 for j=1:n Delta_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); endendfor i=1:n-m-1 Delta_Q(i)=Delta_Q(i)-(Q(i)-QG(i);endDelta_QDelta_V=zeros(1,m);for i=1:m for j=1:n if Type(j)=2 Delta_V(i)=U0(j)2-(e(j)2+f(i)2); end endendDelta_Vnum=0;disp(第,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_the Delta_u_%disp(第一次最大失配量误差：,num2str(MAX)%-循环判断-% if MAX1e-004 % 判断依据 disp(-下面开始下一次迭代过程！-)endwhile MAX1e-004 num=num+1;%形成雅克比矩阵%-先求非对角元素-(H)-%Hik=zeros(n-1,n-1);for i=1:n-1 for 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); end endend%-再求对角元素-（H）-%for i=1:n-1 for k=1:n if 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); end end Hik(i,i)=U_first(i)*Hik(i,i);endHik%-先求非对角元素-N-%Nik=zeros(n-1,n-m-1);for i=1:n-1 for k=1:n-m-1 if i=k theik=Theta_first(i)-Theta_first(k); Nik(i,k)=-U_first(i)*U_first(k)*(G(i,k)*cos(theik)+B(i,k)*sin(theik); end endend%-再求对角元素-%for i=1:n-m-1 for k=1:n if 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); end end Nik(i,i)=-U_first(i)*Nik(i,i)-2*U_first(i)*U_first(i)*G(i,i);endNik%-先求非对角元素-（M）-%Mik=zeros(n-m-1,n-1);for i=1:n-m-1 for 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); end endend %-再求对角元素-% for i=1:n-m-1 for k=1:n if 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); end end Mik(i,i)=-U_first(i)*Mik(i,i);endMik%-先求非对角元素-（L）-%Lik=zeros(n-m-1,n-m-1);for i=1:n-m-1 for 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); end endend%-再求对角元素-%for i=1:n-m-1 for k=1:n if i=k theik=Theta_first(i)-Theta_first(k); Lik(i,i)=Lik(i,i)+U_first(k)*(G(i,k)*sin(theik)-B(i,k)*cos(theik); end end 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,U2(重要！)-% DEL=DEL;Delta_F_E=(-1*inv(kacb)*DEL;Delta_F=Delta_F_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_first %修正后的角度值%Delta_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);e f%-计算修正Delta_P%Delta_P=zeros(1,n-1);for i=1:n-1 for k=1:n Delta_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-1 Delta_P(i)=Delta_P(i)-(P(i,1)-PG(i,1);endDelta_P%-Delta_P-计算完成-%-计算Delta_Q-%Delta_Q=zeros(1,n-m-1);for i=1:n-m-1 for k=1:n Delta_Q(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-1 Delta_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:n theij=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:n theij=Theta_first(i)-Theta_first(t); Qs0=Qs0+U_first(t)*(G(i,t)*sin(theij)-B(i,t)*cos(theij);endQs0=U_first(i)*Qs0;Z(i,8)=Qs0;%计算PV节点的无功功率%Qv0=zeros(1,n);for i=n-m:n-1 for t=1:n theij=Theta_first(i)-Theta_first(t); Qv0(i)=Qv0(i)+U_first(t)*(G(i,t)*sin(theij)-B(i,t)*cos(theij); endQv0(i)=U_first(i)*Qv0(i);endfor i=n-m:n-1 Z(i,8)=Qv0(i);endZj=sqrt(-1);for i=1:n Angle(i)=Theta_first(i)*180/pi; P_in(i)=Z(i,7)-Z(i,5)+j*(Z(i,8)-Z(i,6);endfor i=1:n Flow(i,1)=Z(i,1); Flow(i,2)=U_first(i); Flow(i,3)=Angle(i); Flow(i,4)=P_in(i);endFlow(:,1)=Z2(:,1);disp(#节点电压和注入功率#)disp( 节点号 节点电压幅值 节点电压相位 节点注入功率)Flowdisp(=完*成=)test.m的例子5节点%Bus 信息 %1编号 2类型 3U *4 5 P(L) 6Q(L) 7P(G) 8Q(G) 9Uo 10Qmax 11Qmin 12(串并联电容电抗器)NODE=1 0 1.0000 0 0.8055 0.5320 0 0 1 0 0 0 0; 2 0 1.0000 0 0.1800 0.1200 0 0 1 0 0 0 0; 3 0 1.0000 0 0 0 0 0 1 0 0 0 0; 4 2 1.0000 0 0 0 0.5 0 1.0000 0 0 0 0; 5 3 1.0000 0 0 0 0 0 1.0000 0 0 0 0;%-% Branch 信息 %i j R X B K(1)Branch=1 2 0.0250 0.0800 0.1400 0.0000; 1 3 0.0300 0.1000 0.1800 0.0000; 2 3 0.0200 0.0600 0.1000 0.0000; 4 2 0.0000 0.1905 0.0000 1.0522; 5 3 0.0000 0.1905 0.0000 1.0522; %-第 9 页 共 9 页