pfc_nl.m

％ 该Matlab程序基于牛顿－拉夫逊算法

function [U,S,NetLoss,num_Iter]=PFC_NL(Y_B,PQ,PV,UA)
% 该程序基于牛顿－拉夫逊算法，用于计算已知导纳矩阵、PQ节点、PV节点、平衡节点(UA)的电力网络潮流
% U - 各节点母线电压    S - 各节点注入功率    S_net - 电力网络总损耗
% PQ_P - 实算PQ节点注入有功功率    PQ_Q - 实算PQ节点注入无功功率
% delt_PQ_P - 实算PQ节点有功功率修正值     delt_PQ_Q －实算PQ节点无功功率修正值
% delt_UA_P - 实平衡节点有功功率修正值     delt_U_2 - 实平衡节点电压平方修正值
% delt_PQV - 实算P Q U^2修正值   J - 雅可比矩阵
% e - 电压实部   f - 电压虚部    delt_ef - 电压实部与虚部修正值

% Y_B=input('请输入节点导纳矩阵:');
% PQ=input('请输入PQ节点信息矩阵(第一行:序号, 第二行:注入功率值:);
% PV=input('请输入PV节点信息矩阵(第一行:序号, 第二行:有功功率值, 第三行:节点电压大小);
% UA=input('请输入UA节点信息矩阵(第一行:序号, 第二行:节点电压):');

clc
%初始化参数
N=size(Y_B,2);%节点总数
L=size(PQ,2);%PQ节点数
M=size(PV,2);%PV节点数
U=1.000+0.000i*ones(1,N);
S=zeros(1,N);
S(PQ(1,:))=PQ(2,:);
U(UA(1,:))=UA(2,:);
if ~isempty(PV)
   U(PV(1,:))=PV(3,:);
   S(PV(1,:))=PV(2,:);
end
G=real(Y_B);
B=imag(Y_B);
e=real(U);
f=imag(U);
num_Iter=0;

%牛顿-拉夫逊迭代
while(1)
    %对PQ节点,计算delt_P,delt_Q
    delt_PQV=[];
    for k=1:L
        PQ_P(PQ(1,k))=sum(e(PQ(1,k))*(G(PQ(1,k),:).*e-B(PQ(1,k),:).*f)+f(PQ(1,k))*(G(PQ(1,k),:).*f+B(PQ(1,k),:).*e));
        PQ_Q(PQ(1,k))=sum(f(PQ(1,k))*(G(PQ(1,k),:).*e-B(PQ(1,k),:).*f)-e(PQ(1,k))*(G(PQ(1,k),:).*f+B(PQ(1,k),:).*e));
    end
    delt_PQ_P(PQ(1,:))=real(S(PQ(1,:)))-PQ_P(PQ(1,:));
    delt_PQ_Q(PQ(1,:))=imag(S(PQ(1,:)))-PQ_Q(PQ(1,:));
    
    %对PV节点,计算delt_UA_P,delt_U_2
    if ~isempty(PV)
       for k=1:M
           UA_P(k)=sum(e(PV(1,k))*(G(PV(1,k),:).*e-B(PV(1,k),:).*f)+f(PV(1,k))*(G(PV(1,k),:).*f+B(PV(1,k),:).*e));
           U_2(k)=e(PV(1,k))^2+f(PV(1,k))^2;
       end
       delt_UA_P=real(S(PV(1,:)))-UA_P;
       delt_U_2=U(PV(1,:)).^2-U_2;
    end
    %生成delt_PQV
    for k=1:L
        delt_PQV=[delt_PQV delt_PQ_P(PQ(1,k))];
        delt_PQV=[delt_PQV delt_PQ_Q(PQ(1,k))];
    end
    if ~isempty(PV)
       for k=1:M
           delt_PQV=[delt_PQV delt_UA_P(k)];
           delt_PQV=[delt_PQV delt_U_2(k)];
       end
    end
  
    %生成雅克比矩阵
    I(PQ(1,:))=(PQ_P(PQ(1,:))-PQ_Q(PQ(1,:))*i)./U(PQ(1,:));
    if ~isempty(PV)
       I(PV(1,:))=Y_B(PV(1,:),:)*U';
   end
    Ia=real(I);
    Ib=imag(I);
    for t=1:L+M % t行标
        %j=i
        if t<=L%PQ节点
            J(2*t-1,2*t-1)=-B(PQ(1,t),PQ(1,t))*e(PQ(1,t))+G(PQ(1,t),PQ(1,t))*f(PQ(1,t))+Ib(PQ(1,t));%H
            J(2*t-1,2*t)=G(PQ(1,t),PQ(1,t))*e(PQ(1,t))+B(PQ(1,t),PQ(1,t))*f(PQ(1,t))+Ia(PQ(1,t));%N
            J(2*t,2*t-1)=-G(PQ(1,t),PQ(1,t))*e(PQ(1,t))-B(PQ(1,t),PQ(1,t))*f(PQ(1,t))+Ia(PQ(1,t));%J
            J(2*t,2*t)=-B(PQ(1,t),PQ(1,t))*e(PQ(1,t))+G(PQ(1,t),PQ(1,t))*f(PQ(1,t))-Ib(PQ(1,t));%L
        elseif t>L%PV节点 ~isempty(PV)&
            J(2*t-1,2*t-1)=-B(PV(1,t-L),PV(1,t-L))*e(PV(1,t-L))+G(PV(1,t-L),PV(1,t-L))*f(PV(1,t-L))+Ib(PV(1,t-L));%H
            J(2*t-1,2*t)=G(PV(1,t-L),PV(1,t-L))*e(PV(1,t-L))+B(PV(1,t-L),PV(1,t-L))*f(PV(1,t-L))+Ia(PV(1,t-L));%N
            J(2*t,2*t-1)=2*f(PV(1,t-L));%R
            J(2*t,2*t)=2*e(PV(1,t-L));%S
        end
        %j~=i
        if t<=L%PQ节点   s列标
           for s=1:(L-t) %PQ--PQ 上三角阵  
               J(2*t-1,2*t-1+2*s)=-B(PQ(1,t),PQ(1,t+s))*e(PQ(1,t))+G(PQ(1,t),PQ(1,t+s))*f(PQ(1,t));%H
               J(2*t-1,2*t+2*s)=G(PQ(1,t),PQ(1,t+s))*e(PQ(1,t))+B(PQ(1,t),PQ(1,t+s))*f(PQ(1,t));%N
               J(2*t,2*t-1+2*s)=-J(2*t-1,2*t+2*s);%J
               J(2*t,2*t+2*s)=J(2*t-1,2*t+2*s-1);%L
           end
           for s=1:t-1%PQ--PQ 下三角阵
               J(2*t-1,2*t-2*s-1)=-B(PQ(1,t-s),PQ(1,t))*e(PQ(1,t))+G(PQ(1,t-s),PQ(1,t))*f(PQ(1,t));
               J(2*t-1,2*t-2*s)=G(PQ(1,t-s),PQ(1,t))*e(PQ(1,t))+B(PQ(1,t-s),PQ(1,t))*f(PQ(1,t));
               J(2*t,2*t-2*s-1)=-J(2*t-1,2*t-2*s);
               J(2*t,2*t-2*s)=J(2*t-1,2*t-2*s-1);
           end
           for s=1:M %PQ--PV 上三角  
               J(2*t-1,2*(L+s)-1)=-B(PQ(1,t),PV(1,s))*e(PQ(1,t))+G(PQ(1,t),PV(1,s))*f(PQ(1,t));%H
               J(2*t-1,2*(L+s))=G(PQ(1,t),PV(1,s))*e(PQ(1,t))+B(PQ(1,t),PV(1,s))*f(PQ(1,t));%N
               J(2*t,2*(L+s)-1)=-J(2*t-1,2*(L+s));%J
               J(2*t,2*(L+s))=J(2*t-1,2*(L+s)-1);%L
           end
        elseif t>L%PV节点 ~isempty(PV)
           for s=1:L %下三角 PV--PQ
               J(2*t-1,2*s-1)=-B(PV(1,t-L),PQ(1,s))*e(PV(1,t-L))+G(PV(1,t-L),PQ(1,s))*f(PV(1,t-L));%H
               J(2*t-1,2*s)=G(PV(1,t-L),PQ(1,s))*e(PV(1,t-L))+B(PV(1,t-L),PQ(1,s))*f(PV(1,t-L));%N
               J(2*t,2*s-1)=0;%R
               J(2*t,2*s)=0;%S
           end
           for s=1:M-1%上三角 PV--PV
               J(2*t-1,2*t-1+2*s)=-B(PV(1,t-L),PV(1,t-L+s))*e(PV(1,t-L))+G(PV(1,t-L),PV(1,t-L+s))*f(PV(1,t-L));%H
               J(2*t-1,2*t+2*s)=G(PV(1,t-L),PV(1,t-L+s))*e(PV(1,t-L))+B(PV(1,t-L),PV(1,t-L+s))*f(PV(1,t-L));%N
               J(2*t,2*t+2*s-1)=0;%R
               J(2*t,2*t+2*s)=0;%S
           end
           for s=1:t-L-1%下三角 PV--PV
               J(2*t-1,2*t-2*s-1)=-B(PV(1,t-L),PV(1,t-L-s))*e(PV(1,t-L))+G(PV(1,t-L),PV(1,t-L-s))*f(PV(1,t-L));%H
               J(2*t-1,2*t-2*s)=G(PV(1,t-L),PV(1,t-L-s))*e(PV(1,t-L))+B(PV(1,t-L),PV(1,t-L-s))*f(PV(1,t-L));%N
               J(2*t,2*t-2*s-1)=0;%R
               J(2*t,2*t-2*s)=0;%S
           end
        end 
    end
    %J
    delt_ef=J\delt_PQV';
    
    if find(abs(delt_ef)>1e-6)&num_Iter<10
       f(PQ(1,:))=f(PQ(1,:))+delt_ef(1:2:2*L-1)';
       e(PQ(1,:))=e(PQ(1,:))+delt_ef(2:2:2*L)';
       if ~isempty(PV)
          f(PV(1,:))=f(PV(1,:))+delt_ef(2*L+1:2:2*L+2*M-1)';
          e(PV(1,:))=e(PV(1,:))+delt_ef(2*L+2:2:2*L+2*M)';   
       end
   else
       break;
   end
   num_Iter=num_Iter+1;
end

%各节点母线电压
U=e+f*i;
%计算平衡节点功率
for k=1:size(UA,2)
    S(UA(1,k))=U(UA(1,k)).*sum((conj(Y_B(UA(1,k),:)).*conj(U)));
end
%计算PV节点的无功功率Q
if ~isempty(PV)
   for k=1:M
       S(PV(1,k))=S(PV(1,k))+(f(PV(1,k))*(G(PV(1,k),:)*e'-B(PV(1,k),:)*f')-e(PV(1,k))*(G(PV(1,k),:)*f'+B(PV(1,k),:)*e'))*i;
   end    
end
%计算线路功率
for w=1:N
    for v=1:N
        S_circuitrY_B(w,v)=U(w)*(conj(U(w))*conj(sum(-Y_B(w,:)))+(conj(U(w))-conj(U(v)))*conj(-Y_B(w,v)));
    end
end
%S_circuitrY_B
NetLoss=sum(sum(S_circuitrY_B));