main3.m

电力系统潮流计算 采取牛顿-拉弗逊法和PQ分解法 用matlab计算矩阵的功能

function[Upoint,S]=main(Y,S,U)
if(~exist('Y'))%平衡节点下标为（1，1）
S=[0 1.5 -1];
Y=[-20j 10j 10j;10j -20j 10j;10j 10j -20j;];%默认
mm=3;
U=ones(1,mm);% 定义电压值，包括平衡节点
U(1)=1;% 定义U的平衡节点
%其后面的n-mm个是当作PV节点的V使用
end
if(Y==1)
    load Y.mat;load S.mat;load U.mat;
end%读取Y，S，U
S(1)=0;
P=real(S);%定义P
Q=imag(S);%定义Q
   nn=length(P);
   pdQ=~~Q;mm=1+sum(pdQ);%判断有几个Q值
ro=zeros(1,nn);% 定义相位角
deltaUDU=zeros(1,mm);%其与Q的数目相同,但会留一个平衡节点的不用
   G=real(Y);
   B=imag(Y);%定义Y,G,B
   decide=ones(1,mm+nn);%定义最小值
   decide=decide.*10e-5;
  End=0;
 deltaP=zeros(1,nn);deltaQ=zeros(1,mm);%定义deltaP,deltaQ
   
   while(End==0)
   for a=1:nn%计算P的偏移量
        SumP=0;
       for b=1:nn
           SumP=SumP+U(b).*(G(a,b).*cos(ro(a)-ro(b))+B(a,b).*sin(ro(a)-ro(b)));
       end
       deltaP(a)=P(a)-U(a).*SumP;
   end
   for a=1:mm%计算Q的偏移量
        SumQ=0;
       for b=1:mm
           SumQ=SumQ+U(b).*(G(a,b).*sin(ro(a)-ro(b))-B(a,b).*cos(ro(a)-ro(b)));
       end
       deltaQ(a)=Q(a)-U(a).*SumQ;
   end

   Yak=zeros(mm+nn,mm+nn);%建立YKEBI矩阵
   for a=1:nn%H,L子阵和对角
       Hii=0;
       for b=1:nn
           deltaro=ro(a)-ro(b);
           Yak(a,b)=U(a).*U(b).*(G(a,b).*sin(deltaro)-B(a,b).*cos(deltaro));
           if((b<=mm)&&(a<=mm))
           Yak(a+nn,b+nn)=Yak(a,b);%若mm<nn则部分执行
           end
           Hii=Yak(a,b)+Hii;
       end
       Yak(a,a)=-(Hii-Yak(a,a));
       if(a<=mm)
        Yak(a+nn,a+nn)=-Yak(a,a)-2.*U(a).^2.*B(a,a);
       end   
  end
    for a=1:nn %N,J子阵和对角
      Nii=0;
       for b=1:mm
           if(a<=mm)
           deltaro=ro(a)-ro(b);
           Yak(a,b+nn)=U(a).*U(b).*(G(a,b).*cos(deltaro)+B(a,b).*sin(deltaro));
           Yak(a+nn,b)=-Yak(a,b+nn);
           Nii=Nii+Yak(a,b+nn);
           else
           deltaro=ro(a)-ro(b);    
           Yak(a,b+nn)=U(a).*U(b).*(G(a,b).*cos(deltaro)+B(a,b).*sin(deltaro));
           deltaro=-deltaro;
           Yak(b+nn,a)=U(a).*U(b).*(G(b,a).*cos(deltaro)+B(b,a).*sin(deltaro));
           end
       end
       if(a<=mm)
       Yak(a,a+nn)=Nii-Yak(a,a+nn);%还没有加上UG
       Yak(a+nn,a)=Yak(a,a+nn);
       Yak(a,a+nn)=Yak(a,a+nn)+2.*U(a).^2.*G(a,a);
       end
       end
   for a=1:nn+mm%赋零“父”字形
       Yak(1,a)=0;Yak(a,1)=0;
   end
   for a=1:nn+mm
       Yak(a,nn+1)=0;
       Yak(nn+1,a)=0;
   end
   deltaPQ=[deltaP,deltaQ];
   solution=zeros(1,mm+nn);
   solution=(pinv(Yak)*(deltaPQ'))';%结果是deltaro和deltaUDU的组合
   
   deltaron=zeros(1,nn);
   deltaU=zeros(1,mm);%deltaUDU已定义
   deltaron=solution(1,1:nn);        %取数
        deltaUDU=solution(1,(nn+1):(mm+nn));
      deltaU=deltaUDU.*U;
   U=U+deltaU;ro=ro+deltaron;  %delta值加上ro和U 
   if(solution<decide)
    End=1;
   end
   end
   Upoint=U.*exp(ro.*i);%节点电压
   S1=Upoint(1)*sum(conj(Y(1,:)).*conj(Upoint));%计算平衡节点功率
   S(1)=S1;deltaS=sum(S);%计算网络损耗
    save result.mat Upoint S;%输出给result.txt
    S