nrloop.m

程序实现用牛顿－拉夫逊迭代法进行潮流计算。不仅有潮流计算

%*******************************
%Filename:NRloop.m
%Author:Hweel_Zheng(郑奕辉)     
%First created:2008.08.23      
%Last mended:2008.08.25
%******************************
%Newton-Raphson迭代
[deltaPQ,Jac,iP,iQ]=DPQ_jac(bus,Y);
cta=bus(1:iP,3);
U=bus(1:iQ,2);
ctaU=[cta;U];
G=real(Y);B=imag(Y);
multiplierU=[ones(iP,1);U];
pre=1.0e-10;%精度
kmax=10;%迭代最大次数

outputY;%输出节点导纳矩阵

for iloop=1:kmax
    if(max(abs(deltaPQ))<pre)
        break
    end
    DctaU=Jac\deltaPQ;
    DctaU=DctaU.*multiplierU;   %因为模左除后得到的是ΔU/U，所以要乘上乘子U
    ctaU=ctaU-DctaU;
    cta=ctaU(1:iP);
    U=ctaU(iP+1:iP+iQ);
    multiplierU=[ones(iP,1);U]; %Δδ部分不变，所以乘上1
    bus(1:iP,3)=cta;
    bus(1:iQ,2)=U;
    output_iterate;%输出迭代过程中的相关数据
    [deltaPQ,Jac,iP,iQ]=DPQ_jac(bus,Y);
end
%----------------------------------------------------------
%下面求PV节点的无功注入和SW节点的有功和无功注入
cta=bus(:,3);
U=bus(:,2);
PVQ=0;SWP=0;SWQ=0;
for i=1:nb
    if(bus(i,6)==2)
        PVQ=0;
        for jl=1:nb
            PVQ=PVQ+U(i)*U(jl)*(G(i,jl)*sin(cta(i)-cta(jl))-...
                B(i,jl)*cos(cta(i)-cta(jl)));
        end
        bus(i,5)=PVQ;   %求出各PV节点无功功率
    end
end
for jl=1:nb
    SWP=SWP+U(nb)*U(jl)*(G(nb,jl)*cos(cta(nb)-cta(jl))+...
                B(nb,jl)*sin(cta(nb)-cta(jl)));
    SWQ=SWQ+U(nb)*U(jl)*(G(nb,jl)*sin(cta(nb)-cta(jl))-...
                B(nb,jl)*cos(cta(nb)-cta(jl)));
end
bus(nb,4)=SWP;bus(nb,5)=SWQ;    %求出SW节点无功和有功功率
%至此，求出了每个节点上的电压幅值，相角，节点注入有功和无功
%----------------------------------------------------------
%下面求线路潮流和损耗
V=U.*cos(cta)+j.*U.*sin(cta);
S=zeros(nl,3);Sij=S(:,1);Sji=S(:,2);DSij=S(:,3);
for i=1:nl
    li=line(i,1);lj=line(i,2);K=line(i,7);
    Zt=line(i,3)+j*line(i,4);
    if li~=0&lj~=0
        Yt=1/Zt;
    end
    Ym=line(i,5)+j*line(i,6);
    if li~=0&lj~=0&K==0    %普通线路
        Iij=V(li)*(Yt+Ym)-V(lj)*Yt;
        Iji=V(lj)*(Ym+Yt)-V(li)*Yt;
        Sij(i)=V(li)*conj(Iij);
        Sji(i)=V(lj)*conj(Iji);
        DSij(i)=Sij(i)+Sji(i);
    elseif li~=0&lj~=0&K>0 %变压器线路K在j侧
        Iij=V(li)*(Ym+Yt)-V(lj)*(Yt/K);
        Iji=V(lj)*(Yt/(K^2))-V(li)*(Yt/K);
        Sij(i)=V(li)*conj(Iij);
        Sji(i)=V(lj)*conj(Iji);
        DSij(i)=Sij(i)+Sji(i);
    elseif li~=0&lj~=0&K<0 %变压器线路K在i侧
        Iij=V(li)*(Ym+Yt)-V(lj)*(Yt*K);
        Iji=V(lj)*(Yt*K^2)-V(li)*(Yt*K);
        Sij(i)=V(li)*conj(Iij);
        Sji(i)=V(lj)*conj(Iji);
        DSij(i)=Sij(i)+Sji(i);
    else        %接地线路
        Iij=V(li)*Ym;
        Sij(i)=V(li)*conj(Iij);
        Sji(i)=0;
        DSij(i)=Sij(i)+Sji(i);
    end
end
S=[Sij,Sji,DSij];
lineS=[line,S];
%-----------------------------------------------------------
%下面把节点顺序还原为原来的状态
for i=1:nb
    for k=1:nb
        if nodenum(k,2)==i
            tempbus(i,:)=bus(k,:);
        end
    end
end
tempbus(:,3)=tempbus(:,3).*(180/pi);    %变换为角度
tempbus(:,1)=[1:nb]';

for i=1:nl
    for k=1:2
        if lineS(i,k)~=0
            lineS(i,k)=nodenum(lineS(i,k),2);
        end
    end
end
line=lineS(:,1:7);
flow=[lineS(:,1),lineS(:,2),lineS(:,8:10)];