nf2.m

电力系统可靠性中的潮流计算 使用牛啦法

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%main function nf
%to calculate power flow in normal state
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [databus,dataline,tp,tp1,Z0,M]=nf(NB,NL,databus,dataline)
[YB,B]=form_y(NB,NL,databus,dataline);
[PQnum,U,T,tp,tp1,YB,Z0,M]=nbus(NB,NL,databus,dataline,YB,B);
[U,T,Niterate]=iterate(NB,databus,tp,U,T,YB,PQnum);
if(Niterate==50)
  error('潮流不收敛');
else
  [databus,dataline,sl]=out(NL,PQnum,NB,YB,U,T,databus,dataline,tp,tp1);
end

%***********************************************************************
%subfunction form_y 
%to form the bus acceptance matrix YB
%***********************************************************************
function [YB,B]=form_y(NB,NL,databus,dataline)
% 用MATLAB 语言来实现潮流计算以及静态稳定和暂态稳定计算
% 节点个数NB 支路个数NL
% 节点数据矩阵 databus
% (节点类型 节点号 发电机有功 发电机无功 负荷有功 负荷无功 电压 相角）
% （节点类型 1 PQ节点 0 PV节点）
% 支路数据矩阵 dataline
% (起始节点号 终止节点号 支路电阻 支路电抗 支路变比 对地导纳）

% 形成支路导纳向量YL=(实部（YL）+虚部（YL）*i)
for k=1:NL
  temp1=dataline(k,3);
  temp2=dataline(k,4);
  BL(k)=-1/temp2;
  temp3=temp1/(temp1*temp1+temp2*temp2);
  temp4=temp2/(temp1*temp1+temp2*temp2);
  YL(k)=temp3-temp4*i;
end

% 形成节点导纳矩阵YB(=实部+虚部*i)
for k1=1:NB
  for k2=1:NB
     YB(k1,k2)=0.0+0.0*i;
     B(k1,k2)=0.0;
  end
end
for k=1:NL   
  n1=dataline(k,1);n2=dataline(k,2);n3=dataline(k,6);
  YB(n1,n1)=YB(n1,n1)+YL(k)+n3*i;
  B(n1,n1)=B(n1,n1)+BL(k);
  YB(n2,n2)=YB(n2,n2)+YL(k)+n3*i;
  B(n2,n2)=B(n2,n2)+BL(k);
  YB(n1,n2)=YB(n1,n2)-YL(k);
  YB(n2,n1)=YB(n2,n1)-YL(k);
  B(n1,n2)=B(n1,n2)-BL(k);
  B(n2,n1)=B(n2,n1)-BL(k);
end
for k=1:NL
 if(dataline(k,5)~=1)
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 n1=dataline(k,1);n2=dataline(k,2);n3=dataline(k,5);
 YB(n1,n1)=YB(n1,n1)+(1-n3)/(n3*n3)*YL(k);
 YB(n2,n2)=YB(n2,n2)+(n3-1)/n3*YL(k);
 YB(n1,n2)=0.0-1/n3*YL(k);
 YB(n2,n1)=0.0-1/n3*YL(k);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%the following is original code   
%%   n1=dataline(k,1);n2=dataline(k,2);n3=dataline(k,5);
%%   YB(n1,n1)=YB(n1,n1)+(1/(n3*n3)-1)*YL(k);
%%   YB(n1,n2)=0.0-1/n3*YL(k);
%%   YB(n2,n1)=0.0-1/n3*YL(k);
 end
end 

 
 
%*************************************************************************
%subfunction nbus
%to rearrange the order of bus 
%the first is PQ bus followed by PV bus, the last is swing bus. 
%the tp is incidence vector for old number to new; tp1 is for new to old 
%to form the bus-line related matrix M
%to form bus impendance matrix Z0 under the normal state
%*************************************************************************
function [PQnum,U,T,tp,tp1,YB,Z0,M]=nbus(NB,NL,databus,dataline,YB,B)
% 节点个数NB 支路个数NL
% 节点数据矩阵 databus
% (节点类型 节点号 发电机有功 发电机无功 负荷有功 负荷无功 电压 相角 对地导纳）
% （节点类型 1 PQ节点 0 PV节点）
% 支路数据矩阵 dataline
% (起始节点号 终止节点号 支路电阻 支路电抗 支路变比）

%PQ节点在前(1--PQnum)，PV节点其后，最后为平衡节点.
% tp 新编号对应原始编号 tp1 原始编号对应新编号
k1=1;k2=NB-1;
for k=1:NB
  if(databus(k,1)==1)
    tp(k1)=databus(k,2);
    tp1(databus(k,2))=k1;
    k1=k1+1;
  elseif(databus(k,1)==0)
    tp(k2)=databus(k,2);
    tp1(databus(k,2))=k2;
    k2=k2-1;
  else
    tp(NB)=databus(k,2);
    tp1(databus(k,2))=NB;
  end
end
PQnum=k1-1;
for k1=1:NB
  for k2=1:NB 
     Y(k1,k2)=YB(k1,k2);
     Bt(k1,k2)=B(k1,k2);
  end
end
for k1=1:NB
  U(k1)=databus(tp(k1),7);
  T(k1)=databus(tp(k1),8);
    E(k1)=U(k1)*cos(T(k1));
  F(k1)=U(k1)*sin(T(k1));
  for k2=1:NB
     YB(k1,k2)=Y(tp(k1),tp(k2));
     B(k1,k2)=Bt(tp(k1),tp(k2));
  end
end
%to form the bus-line related matrix M
%to form bus impendance matrix Z0 under the normal state
Bb=eye(NB-1,NB)*B*eye(NB-1,NB)';
Zt=zeros(1,NB-1);
Z0=[inv(Bb) Zt';Zt 0];
M=zeros(NB,NL);
for k1=1:NL
   d1=dataline(k1,1);
   d2=dataline(k1,2);
   M(tp1(d1),k1)=1;
   M(tp1(d2),k1)=-1;
end

   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%subfuction iterate
%to do iteration process
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [U,T,Niterate]=iterate(NB,databus,tp,U,T,YB,PQnum)
% 节点个数NB 支路个数NL
% 节点数据矩阵 databus
% (节点类型 节点号 发电机有功 发电机无功 负荷有功 负荷无功 电压 相角 对地导纳）
% （节点类型 1 PQ节点 0 PV节点）
% 支路数据矩阵 dataline
% (起始节点号 终止节点号 支路电阻 支路电抗 支路变比）
% tp 新编号对应原始编号 tp1 原始编号对应新编号
% min 收敛精度要求

min=0.0001;
for k1=1:50
  max=0.0;
  diatW=diatva(NB,U,T,YB,tp,databus,PQnum);
  J=jmatrix(NB,U,T,YB,PQnum);
%*******************************************************************
%解修正方程
%*******************************************************************
  diatU=0.0-J/diatW;
  for k2=1:2*(NB-1)
    if(abs(diatW(k2))>max)
      max=abs(diatW(k2));
    end
  end
  for k2=1:(NB-1)
    if(k2/(1+1)==0)
      F(k2)=F(k2)+diatU(k2);
    else
      E(k2)=E(k2)+diatU(k2);
    end
  end
  for k2=1:(NB-1)
     U(k2)=sqrt(E(k2)^2+F(k2)^2);
     T(k2)=atan(E(k2)/F(k2));
  end
  if(max<min),break;end
end
Niterate=k1;
 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%subfunction
%to calculate the power mismatch
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function diatW=diatva(NB,U,T,YB,tp,databus,PQnum)
% PQnum PQ节点个数
% U 电压（对应节点新编号）
% T 相角（对应节点新编号）
% diatP 有关偏差
% diatQ 无功偏差
% diatW 功率偏差
% YB 对应新编号的节点导纳矩阵
% databus 节点数据
% tp(k) 新编号k对应的原节点号tp(k)
for k2=1:NB-1 
  %t1(k2)=0.0;t2(k2)=0.0;t3(k2)=0.0;
  sump(k2)=0;
end
for k2=1:PQnum 
  sumq(k2)=0;
end
for k2=1:NB-1-PQnum
    sumv(k2)=0;
end

for k2=1:NB-1  
  for k3=1:NB %%%%%%%%%%%
    t1(k2)=E(k2)*(real(YB(k2,k3))*E(k3)-imag(YB(k2,k3))*F(k3))+F(k2)*(real(YB(k2,k3))*F(k3)+imag(YB(k2,k3))*E(k3));   %%%%%%%%%%
  sump(k2)=sump(k2)+t1(k2);

  end
end
 %全部PQ,PV点的Pi
for k2=1:PQnum 
  for k3=1:NB
      t2(k2)=F(k2)*(real(YB(k2,k3))*E(k3)-imag(YB(k2,k3))*F(k3))-E(k2)*(real(YB(k2,k3))*F(k3)+imag(YB(k2,k3))*E(k3));    %%%%%%%%%%%%%%%%
     sumq(k2)=sumq(k2)+t2(k2);
     
  end
end
%全部PQ点的Qi
for k2=1:NB-PQnum-1
  for k3=PQnum+1:NB-1
       t3(k3)=sqrt(E(k3)^2+F(k3)^2);
       sumv(k2)=sumv(k2)+t3(k3);
      
    end   
end
  
for k2=1:NB-1
  diatP(k2)=databus(tp(k2),3)-databus(tp(k2),5)-sump(k2)
end

for k2=1:PQnum    
  diatQ(k2)=databus(tp(k2),4)-databus(tp(k2),6)-sumq(k2)
end
for k2=1:NB-PQnum-1
   for k3=PQnum+1:NB-1
       diatV(k3)=databus(tp(k3),7)-sumv(k2)
   end
end
for k=1:2*(NB-1)
    for t=1:(NB-PQnum-1)
  if(k<=PQnum)
    diatW(k)=diatP(k);
    diatW(k+NB-1)=diatQ(k);
  end
  if((k>PQnum)&(k<=NB-1))
    diatW(k)=diatP(k)/databus(tp(k),7);end                %PV点在YB阵的编号是K吗？
  if((k>NB+PQnum-1)&(k<=(2*(NB-1))))
    diatW(k)=diatV(k-NB+1);end
 end 
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%subfuction jmatrix
%to form Jacobian matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function J=jmatrix(NB,U,T,YB,PQnum)
% J 雅可比矩阵 (NB+PQnum-1)*(NB+PQnum-1)维
% NB 节点个数
% U(k) 节点电压 k为新节点编号
% T(k) 节点相角
% YB 节点导纳矩阵 NB*NB维
for x=1:(NB-1)
  for j=1:2*(NB-1)
    if(rem(j,2)==0) 
        if(x==j/2)
            J(x,j)=-imag(YB(x,j/2))-databus(tp(j/2),4)-databus(tp(j/2),6)   %i点的无功
           
        else
            J(x,j)=-imag(YB(x,j/2));
        end
    else
         if(x==j/2)
           J(x,j)=real(YB(x,(j+1)/2))+databus(tp((j+1)/2),3)+databus(tp((j+1)/2),5)   %i点的有功 
         else
           J(x,j)=real(YB(x,(j+1)/2));
        end
    end
  end
end
for x=NB:(NB+PQnum-1)
  for j=1:2*(NB-1)
    if(rem(j,2)==0)  
        if((x-NB+1)==j/2)
            J(x,j)=real(YB(x-NB+1,j/2))-databus(tp(j/2),3)-databus(tp(j/2),5)%i点的有功 
        else
            J(x,j)=real(YB(x-NB+1,j/2))
        end        
    else
         if((x-NB+1)==j/2)
           J(x,j)=imag(YB(x-NB+1,(j+1)/2))-databus(tp((j+1)/2),4)-databus(tp((j+1)/2),6)  %i点的无功
        else
            J(x,j)=imag(YB(x-NB+1,(j+1)/2))
        end
    end
  end
end
for x=(NB+PQnum):2*(NB-1)
  for j=1:(NB-1)
    for k=1:(NB-PQnum-1)
      if(2*j-1==k)
          J(x,2*j+1)=-1.05
      else
        
          J(x,2*j)=0
      end
    end
  end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%subfunction out
%to calculate line flows and losses
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [databus,dataline,sl]=out(NL,PQnum,NB,YB,U,T,databus,dataline,tp,tp1)
% 平衡节点的有功和无功计算
databus(tp(NB),3)=databus(tp(NB),5);
databus(tp(NB),4)=databus(tp(NB),6);
for k=1:NB
  databus(tp(NB),3)=databus(tp(NB),3)+U(NB)*U(k)*real(YB(NB,k))*cos(T(k));
  databus(tp(NB),3)=databus(tp(NB),3)+U(NB)*U(k)*imag(YB(NB,k))*sin(-T(k));
  databus(tp(NB),4)=databus(tp(NB),4)+U(NB)*U(k)*real(YB(NB,k))*sin(-T(k));
  databus(tp(NB),4)=databus(tp(NB),4)-U(NB)*U(k)*imag(YB(NB,k))*cos(T(k));
  databus(tp(k),7)=U(k);
  databus(tp(k),8)=T(k);
end
%PV节点的无功计算
for k=PQnum:(NB-1)
  Q=0.0;
  for k1=1:NB
    Q=Q+U(k)*U(k1)*real(YB(k,k1))*sin(T(k)-T(k1));
    Q=Q-U(k)*U(k1)*imag(YB(k,k1))*cos(T(k)-T(k1));
  end
  databus(tp(k),4)=Q+databus(tp(k),6);
end
% 系统网损计算
pl=0.0;ql=0.0;
for k=1:NB
   pl=pl+databus(k,3)-databus(k,5);
   ql=ql+databus(k,4)-databus(k,6);
end
sl=pl+ql*i;
% 线路功率
for k=1:NL
  t1=tp1(dataline(k,1));
  t2=tp1(dataline(k,2));
  t3=dataline(k,3);
  t4=dataline(k,4);
  t5=cos(T(t1))+sin(T(t1))*i;
  t51=cos(T(t1))-sin(T(t1))*i;
  t6=cos(T(t2))+sin(T(t2))*i;
  t61=cos(T(t2))-sin(T(t2))*i;
  t7=-dataline(k,6);
  se=U(t1)*U(t1)*t7*i+U(t1)*t5*(U(t1)*t51-U(t2)*t61)/(t3-t4*i);
  es=U(t2)*U(t2)*t7*i+U(t2)*t6*(U(t2)*t61-U(t1)*t51)/(t3-t4*i);
  dataline(k,7)=real(se);  dataline(k,9)=real(es);
  dataline(k,8)=imag(se);  dataline(k,10)=imag(es);
end
save myout.mat databus dataline sl;
disp('The file myout.mat save the result for power flow under normal state');