nrflowcal.m

本程序用牛顿-拉夫逊法计算电力系统潮流

clc 
%本程序用来计算潮流,将平衡节点设置成最后一个节点。程序的例子是《电力系统分析》（下册 何仰赞 华中科技大学出版）P66
%例11-6。程序中的参数：n,nl,pr,max1,Z,A,都可视情况改变，或者用input函数在程序运行的时候输入。
n=4;%节点个数
nl=4;%支路条数
max1=100;%最大迭代次数
pr=0.00001;%精度
Z=[1 2 0.10+0.40j 0.01528j 0.01528j 0;
1 3 1.1 0.3j 0 2;
2 4 0.08+0.4j 0.01413j 0.01413j 0;
1 4 0.12+0.5j 0.01920j 0.01920j 0;];%由支路参数形成的矩阵Z=[节点号1 节点号2 支路阻抗 靠近节点1的对地导纳 靠近节点2的对地导纳 参数]'); 参数的详细资料请查阅本程序附带的文本。
A=[-0.30-0.18j 1 1;
-0.55-0.13j 1 1;
0.5 1.1 2;
0 1.05 3];%节点的功率，电压矩阵 A=[节点输入功率 节点电压 节点类型（PQ节点的节点类型是1；PV：2；平衡节点：3）]

for i=1:nl
    if Z(i,6)==1
        K=Z(i,3);Zt=Z(i,4);
        Z(i,3)=K*Zt;Z(i,4)=(1-K)/(K^2*Zt);Z(i,5)=(K-1)/(K*Zt);
    elseif Z(i,6)==2
         K=Z(i,3);Zt=Z(i,4);
         Z(i,3)=Zt/K;Z(i,4)=K*(K-1)/Zt;Z(i,5)=(1-K)/Zt;
    elseif Z(i,6)==3
         Zt=Z(i,3);K=Z(i,4);
         Z(i,3)=K*Zt;Z(i,4)=(K-1)/(K*Zt);Z(i,5)=(1-K)/(K^2*Zt);
    elseif Z(i,6)==4
         Zt=Z(i,3);K=Z(i,4);
         Z(i,3)=Zt/K;Z(i,4)=(1-K)/Zt;Z(i,5)=K*(K-1)/Zt;
    else
        
    end
end           %对含有变压器的支路进行处理，详见本程序附带的文本。
Y=zeros(n);
for i=1:nl
    p=Z(i,1);q=Z(i,2);
    Y(p,q)=-1/Z(i,3);
    Y(q,p)=Y(p,q);
    Y(p,p)=Y(p,p)+1/Z(i,3)+Z(i,4);
    Y(q,q)=Y(q,q)+1/Z(i,3)+Z(i,5);
end
disp('结点导纳矩阵Y为：');  %计算Y矩阵
disp(Y);
G=real(Y);B=imag(Y);
S=zeros(1,n-1);V=zeros(1,n);cta=zeros(1,n);
for i=1:n-1
     S(i)=A(i,1);
end
for i=1:n
    cta(i)=angle(A(i,2));
    V(i)=abs(A(i,2));
end
P=real(S);Q=imag(S);
for k=1:max1
    t1=1;t2=1;
    
    for i=1:n-1
        C(i)=0;D(i)=0;
        for j1=1:n
            C(i)=C(i)+V(i)*V(j1)*G(i,j1)*cos(cta(i)-cta(j1))+V(i)*V(j1)*B(i,j1)*sin(cta(i)-cta(j1));
            D(i)=D(i)+V(i)*V(j1)*G(i,j1)*sin(cta(i)-cta(j1))-V(i)*V(j1)*B(i,j1)*cos(cta(i)-cta(j1));
        end
        DP(t1)=P(i)-C(i);
        t1=t1+1;
        if A(i,3)==1
            DQ(t2)=Q(i)-D(i);
            t2=t2+1;
        end
    end
    t1=t1-1;t2=t2-1;
    DPQ=[DP';DQ'];                 %计算P,Q的偏差
    if norm(DPQ,inf)<pr            %如果DPQ小于规定的精度则停止迭代
        break;
    end
    H=zeros(t1,t1);N=zeros(t1,t2);K=zeros(t2,t1);L=zeros(t2,t2);
    for i=1:t1
        for j1=1:t1
            if j1~=i
                H(i,j1)=0-V(i)*V(j1)*(G(i,j1)*sin(cta(i)-cta(j1))-B(i,j1)*cos(cta(i)-cta(j1)));
            elseif j1==i
                H(i,j1)=V(i)^2*B(i,j1)+D(i);
            end
        end
    end
    for i=1:t1
        for j1=1:t2
            if j1~=i
                N(i,j1)=0-V(i)*V(j1)*(G(i,j1)*cos(cta(i)-cta(j1))+B(i,j1)*sin(cta(i)-cta(j1)));
            elseif j1==i
                N(i,j1)=0-V(i)^2*G(i,j1)+C(i);
            end
        end
    end
    for i=1:t2
        for j1=1:t1
            if j1~=i
                K(i,j1)=V(i)*V(j1)*(G(i,j1)*cos(cta(i)-cta(j1))+B(i,j1)*sin(cta(i)-cta(j1)));
            elseif j1==i
                K(i,j1)=V(i)^2*G(i,j1)-C(i);
            end
        end
    end
    for i=1:t2
        for j1=1:t2
            if j1~=i
               L(i,j1)=0-V(i)*V(j1)*(G(i,j1)*sin(cta(i)-cta(j1))-B(i,j1)*cos(cta(i)-cta(j1)));
            elseif j1==i
                L(i,j1)=V(i)^2*B(i,j1)-D(i);
            end
        end
    end
    J=[H,N;K,L];  %计算雅可比矩阵
    modify=-J\DPQ;
    Dcta=modify(1:t1,:);
    t3=V(:,1:t2);

    DV=diag(t3,0)*modify(t1+1:t1+t2,:);
    for i=1:t1
        cta(1,i)=cta(1,i)+Dcta(i,1);
    end
    t4=1;
    for i=1:t2
        if A(i,3)==1
            V(1,i)=V(1,i)+DV(t4,1);
            t4=t4+1;
        end
    end
end
disp('V=');
disp(V);
VV=zeros(1,n);
for i=1:n
    VV(i)=V(i)*cos(cta(i))+1j*V(i)*sin(cta(i));  %计算出复数表示的电压
end
for p=1:n
    c(p)=0;
    for q=1:n
        c(p)=c(p)+conj(Y(p,q))*conj(VV(q));
    end
    s(p)=VV(p)*c(p);                             %计算出各个节点的复功率
    Ploss(p)=real(s(p));                         %计算出各个节点的有功损耗
end
PlossSum=0;
for p=1:n
    PlossSum=PlossSum+Ploss(p);
end
disp('s=');
disp(s);