main1.m

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

function[Upoint,S]=main1(Y,S)
if(~exist('Y'))%平衡节点下标为（1，1）
S=[0,0.2+0.2i,-0.45-0.15i,-0.4-0.05i,-0.6-0.1i];
Y=[6.25-18.75i	-5+15i	-1.25+3.75i	0 0;
-5+15i	10.834-32.5i -1.667+5i -1.667+5i -2.5+7.5i;
-1.25+3.75i	-1.667+5i 12.917-38.75i	-10+30i	0;
0 -1.667+5i	-10+30i	12.917-38.75i -1.25+3.75i;
0 -2.5+7.5i	0 -1.25+3.75i 3.75-11.25i;];%默认
mm=5;
U=ones(1,mm);% 定义电压值，包括平衡节点
U(1)=1.06;% 定义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);% 定义相位角
deltaU=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

   BP=zeros(nn,nn);BPP=zeros(mm,mm);   %定义Ｂ＇和Ｂ＇＇
 BP(2:nn,2:nn)=B(2:nn,2:nn);BPP(2:mm,2:mm)=B(2:mm,2:mm);
  
 
 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
   deltaron=(-pinv(BP)*(deltaP./U)')'./U;   %解算矩阵方程BP
   ro=ro+deltaron;%根据delta修正ro
    for a=2: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
   deltaU=(-pinv(BPP)*(deltaQ./U)')'; %解算矩阵方程BP
   U=U+deltaU;%根据delta修正U
   solution=[deltaU(1:mm),deltaron];
    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);%计算网络损耗
  SL=zeros(nn,nn);%计算线路潮流SL
   for a=1:nn
       for b=1:nn
   SL(a,b)=U(a).*(conj(U(a))-conj(U(b))).*conj(Y(a)-Y(b));
       end
   end
   save result.mat Upoint S SL;%输出给result.txt
   Upoint 
   S
   SL