%*******************************
%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)];