1至m为PQ节点，m+1至n-1为PV节点，n为平衡节点

%节点排序：1至m为PQ节点，m+1至n-1为PV节点，n为平衡节点
clear all;
clc;
data;
Y=zeros(n);
for i=1:n
    if X(i,2)~=0;
        p=X(i,1);
        Y(p,p)=1./X(i,2);
    end
end
for i=1:nl
    if B1(i,6)==0
        p=B1(i,1);q=B1(i,2);
    else
        p=B1(i,2);q=B1(i,1);
    end
    Y(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5));
    Y(q,p)=Y(p,q);
    Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2;
    Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2;
end
disp('导纳矩阵Y=');
disp(Y);
%--------------------------------------------------------------------------
G=real(Y);B=imag(Y);
for i=1:n
    e(i)=real(B2(i,3));
    f(i)=imag(B2(i,3));
    V(i)=B2(i,4);%PV节点U的给定值
end
for i=1:n
    S(i)=B2(i,1)-B2(i,2);%注入功率
    B(i,i)=B(i,i)+B2(i,5);%无功补偿
end
Ps=real(S);Qs=imag(S);
num=0;control=1;J=zeros(2*n-2,2*n-2);b=zeros(2*n-2,1);
while control~=0
    for i=1:n-1
        C=0;
        D=0;
        for j=1:n
            C=C+G(i,j)*e(j)-B(i,j)*f(j);
            D=D+G(i,j)*f(j)+B(i,j)*e(j);
        end
            P(i)=e(i)*C+f(i)*D;
            Q(i)=f(i)*C-e(i)*D;
            V2=e(i)^2+f(i)^2;
            if i<=m
                DP=Ps(i)-P(i);
                DQ=Qs(i)-Q(i);%求PQ节点
                for j=1:n-1
                    if  i==j
                        X1=-C-G(i,i)*e(i)-B(i,i)*f(i);
                        X2=-D+B(i,i)*e(i)-G(i,i)*f(i);
                        X3=D+B(i,i)*e(i)-G(i,i)*f(i);
                        X4=-C+G(i,i)*e(i)+B(i,i)*f(i);
                        p=2*i-1;q=2*j-1;
                        u=p+1;v=q+1;
                        J(p,q)=X1;J(p,v)=X2;J(u,q)=X3;J(u,v)=X4;b(p,1)=DP;b(u,1)=DQ;
                    elseif i~=j
                        X1=-G(i,j)*e(i)-B(i,j)*f(i);
                        X2=B(i,j)*e(i)-G(i,j)*f(i);
                        X3=X2;
                        X4=-X1;
                        p=2*i-1;q=2*j-1;
                        u=p+1;v=q+1;
                        J(p,q)=X1;J(p,v)=X2;J(u,q)=X3;J(u,v)=X4;b(p,1)=DP;b(u,1)=DQ;              
                    end
                end
            else
                DP=Ps(i)-P(i);
                DV=V(i)^2-V2;%求PV节点
                for j=1:n-1
                    if  i==j
                        X1=-C-G(i,i)*e(i)-B(i,i)*f(i);
                        X2=-D+B(i,i)*e(i)-G(i,i)*f(i);
                        X5=-2*e(i);
                        X6=-2*f(i);
                        p=2*i-1;q=2*j-1;
                        u=p+1;v=q+1;
                        J(p,q)=X1;J(p,v)=X2;J(u,q)=X5;J(u,v)=X6;b(p,1)=DP;b(u,1)=DV;
                    elseif i~=j
                        X1=-G(i,j)*e(i)-B(i,j)*f(i);
                        X2=B(i,j)*e(i)-G(i,j)*f(i);
                        X5=0;
                        X6=0;
                        p=2*i-1;q=2*j-1;
                        u=p+1;v=q+1;
                        J(p,q)=X1;J(p,v)=X2;J(u,q)=X5;J(u,v)=X6;b(p,1)=DP;b(u,1)=DV;              
                    end                
                end
            end
        end
%--------------------------------------------------------------------------------
    J=-J;
    x=J\b;
%--------------------------------------------------------------------------------
    for i=1:n-1
    e(i)=e(i)+x(2*i-1);
    f(i)=f(i)+x(2*i);
    end
    ABS=abs(b);
    MAX=max(ABS);
    if MAX<=pr
        control=0;
    else
        num=num+1;
    end
    for i=1:n
        Dy(num,i)=sqrt(e(i)^2+f(i)^2);
    end
end
%--------------------------------------------------------------------------------
for i=1:n
    V(i)=sqrt(e(i)^2+f(i)^2);
    O(i)=atan(f(i)./e(i))*180./pi;
end
E=e+f*j;
disp('迭代次数=');
disp(num);
disp('各节点的实际电压标幺值E为（节点号从小到大排列）：');
disp(E);
disp('各节点的电压大小V为（节点号从小到大排列）：');
disp(V);
disp('各节点的电压相角O为（节点号从小到大排列）：');
disp(O);
for p=1:n
    C(p)=0;
    for q=1:n
        C(p)=C(p)+(Y(p,q))*(E(q));
    end
    S(p)=E(p)*conj(C(p));
end
disp('各节点的功率S为（节点号从小到大排列）：');
disp(S);
disp('各条支路的首端功率Si为（顺序同您输入B1时一样）：');
for i=1:nl
    if B1(i,6)==0
        p=B1(i,1);q=B1(i,2);
    else p=B1(i,2);q=B1(i,1);
    end
    Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4))+(conj(E(p)*B1(i,5))-conj(E(q)))*conj(1./(B1(i,3)*B1(i,5))));
    disp(Si(p,q));
end
disp('各条支路的末端功率Sj为（顺序同您输入B1时一样）：')
for i=1:nl
    if B1(i,6)==0
        p=B1(i,1);q=B1(i,2);
    else p=B1(i,2);q=B1(i,1);
    end
    Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4))+(conj(E(q)./B1(i,5))-conj(E(p)))*conj(1./(B1(i,3)*B1(i,5))));
    disp(Sj(q,p));
end
disp('各条支路的功率损耗DS为（顺序同您输入B1时一样）：');
for i=1:nl
    if B1(i,6)==0
        p=B1(i,1);q=B1(i,2);
    else p=B1(i,2);q=B1(i,1);
    end
    DS(i)=Si(p,q)+Sj(q,p);
    disp(DS(i));
end
for i=1:num
    K(i)=i;
end
disp('以下是每次迭代后个节点的电压值（如图所示）');
plot(K,Dy),xlabel('迭代次数'),ylabel('电压'),title('电压迭代次数曲线');

 

D137