牛拉潮流程序

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);%将节点排序重新存储节点信息%%%%定义为节点数的方阵formatshortb=1;%排序标志位pq=0;%PQ节点标志位pv=0;%PV节点标志位ph=0;%平衡节点标志位%---------------------按照PQ,PV,平衡节点的次序排序各种节点%-------------------------------------统计PQ节点数0代表是pq节点fora=1:Total_of_BusifNODE(a,2)==0Z(b,:)=NODE(a,:);b=b+1;pq=pq+1;endend%-------------------------------------统计PV节点数2代表pv节点fora=1:Total_of_BusifNODE(a,2)==2Z(b,:)=NODE(a,:);b=b+1;pv=pv+1;endend%-------------------------------------统计平衡节点数3代表平衡节点fora=1:Total_of_BusifNODE(a,2)==3Z(b,:)=NODE(a,:);b=b+1;ph=ph+1;endendZZ2=Z;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%将节点进行重新排序%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%mm=zeros(n,1);fori=1:nmm(i,1)=i;endZ1(:,1)=mm(:,1);Branch1=zeros(nbr,2);fori=1:nifZ(i,1)~=Z1(i)forj=1:nbrifBranch(j,1)==Z(i,1)Branch1(j,1)=Z1(i);endifBranch(j,2)==Z(i,1)Branch1(j,2)=Z1(i);endendelseforj=1:nbrifBranch(j,1)==Z(i,1)Branch1(j,1)=Z(i,1);endifBranch(j,2)==Z(i,1)Branch1(j,2)=Z(i,1);endendendendBranch(:,1)=Branch1(:,1);Branch(:,2)=Branch1(:,2);Z(:,1)=Z1(:,1);j=sqrt(-1);%--------矩阵已经完成按照PQ，PV，平衡节点的顺序排列起来------YSNODE=Z;%保存排序后的原始节点数据%%%%=======================================================================Y=zeros(n,n);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%求互导纳%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%fori=1:nfort=1:nbrif(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));elseif(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));endendendend%%%%%%%%%%%%%%%%%%%%%%%%%求自导纳%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%fori=1:nfort=1:nbrif(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);elseifBranch(t,1)==i&&Branch(t,6)~=0%%%%%%变压器支路——————————且i为首节点%%%%%%%%Y(i,i)=Y(i,i)+1/(j*Branch(t,4));elseifBranch(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));endendendendend%%%%%%%%%若有并联电容器组，则自导纳要加上并联电容器的导纳%%%%%%%%%%%%%%%%%%%%%%fori=1:nifNODE(i,13)~=0Y(i,i)=Y(i,i)+j*NODE(i,13)endendYn=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);fori=1:n-1forj=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));endendfori=1:n-1Delta_P(i)=Delta_P(i)-(P(i)-PG(i));endDelta_P%%%%%%%%%%%%%%%%%%计算Delta_Q初始功率量%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%m=0;fori=1:n;ifType(i)==2;%计算PV节点的个数m=m+1;%m代表pv节点个数endendDelta_Q=zeros(1,n-m-1);fori=1:n-m-1forj=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));endendfori=1:n-m-1Delta_Q(i)=Delta_Q(i)-(Q(i)-QG(i));endDelta_QDelta_V=zeros(1,m);fori=1:mforj=1:nifType(j)==2Delta_V(i)=U0(j)^2-(e(j)^2+f(i)^2);endendendDelta_Vnum=0;disp(['第',num2str(num),'次时的Delta总的失配量为：'])%--------------------------进入循环体判断是否满足条件--------------------------------%%--------------------先算出最大值，作为判断是否收敛的依据----------------------------%DEL=[Delta_PDelta_Q];%Delta_PDelta_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_theDelta_u________________%%disp(['第一次最大失配量误差：',num2str(MAX)])%-------------循环判断-----------%ifMAX1e-004%判断依据disp('-----------------------下面开始下一次迭代过程！-----------------------')endwhileMAX1e-004num=num+1;%%%%%%%%%%%%%%%%%%%形成雅克比矩阵%%%%%%%%%%%%%%%%%%%%%%%%%%%%--------------------------先求非对角元素--(H)-------------%Hik=zeros(n-1,n-1);fori=1:n-1fork=1:n-1ifi~=ktheik=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));endendend%----------------------再求对角元素---（H）---------------------%fori=1:n-1fork=1:nifi~=ktheik=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));endendHik(i,i)=U_first(i)*Hik(i,i);endHik%-------------------------先求非对角元素---N-------------------%Nik=zeros(n-1,n-m-1);fori=1:n-1fork=1:n-m-1ifi~=ktheik=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));endendend%-------------------------再求对