clc;
t0 = clock;
for i=1:100
j=sqrt(-1);
%% define named indices into bus, gen, branch matrices
    [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
        VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
    [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
        TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
        ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
    [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
        MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
        QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;

    %% default arguments
    [baseMVA, bus, gen, branch, areas, gencost] =case300;
    [i2e, bus, gen, branch]=ext2int(bus, gen, branch);
    [ref, pv, pq]=bustypes(bus, gen);
    pqv=[pq;pv];
    Lpq=length(pq);
    Lpv=length(pv);
    Lpqv=Lpq+Lpv;
    Ltol=Lpqv+1;
    %形成節點導納矩陣
    [Ybus, Yf, Yt]=makeYbus(baseMVA, bus, branch);
    %形成節點注入功率
    Sbus=makeSbus(baseMVA, bus, gen);
    G=real(Ybus);B=imag(Ybus);
    Ubus=ones(Ltol,1)*(1+0*j);
    on=find(gen(:, GEN_STATUS) > 0);                                %運行的機組
    gbus=gen(on, GEN_BUS);
    V=bus(:,VM) ; 
    V(gbus)=gen(on, VG) ;
    Ubus(ref)=V(ref);
    e=real(Ubus);
    f=imag(Ubus);
    Qi=(G*e-B*f).*f-(B*e+G*f).*e;
    Sbus(pv)=Sbus(pv)+j*Qi(pv);
    Y=Ybus(pqv,pqv);
    G=real(Y);
    B=imag(Y);
    Matzr1=zeros(Lpq,Lpv);
    Matzr2=zeros(Lpv);
    %迭代計算
    for iter=1:100
        absV2=abs(Ubus).^2;                                             %迭代過程中的同一量
        absV4=absV2.^2;

        detI=Ybus*Ubus-conj(Sbus./Ubus);                                %計算所有的不平衡量
        detV=bus(:,VM).^2-absV2;
        detb=[real(detI(pqv));imag(detI(pqv));detV(pv)];                %按照PQ、PV節點順序取出不平衡量
 
        %計算雅可比矩陣
        ebus=real(Ubus);
        fbus=imag(Ubus);
        Pbus=real(Sbus);
        Qbus=imag(Sbus);

        %計算雅可比矩陣中變量
        dIre_de=-(Pbus.*(fbus.^2-ebus.^2)-2*ebus.*fbus.*Qbus)./absV4;   %H
        dIre_df=-(Qbus.*(ebus.^2-fbus.^2)-2*ebus.*fbus.*Pbus)./absV4;   %N
        dIim_de=dIre_df;                                                %J
        dIim_df=-(Pbus.*(ebus.^2-fbus.^2)+2*ebus.*fbus.*Qbus)./absV4;   %L
        dV_de=-2*ebus;                                                  %R
        dV_df=-2*fbus;                                                  %S
        dIre_dQ=-fbus./absV2;                                           %K
        dIim_dQ=ebus./absV2;                                            %M

        %形成雅可比矩陣中在迭代過程中發生變化的量(對角元素)
        dIre_de=dIre_de(pqv);                                           %形成H的對角元素組成的向量
        dIre_df=dIre_df(pqv);                                           %形成N的對角元素組成的向量
        dIim_de=dIim_de(pqv);                                           %形成J的對角元素組成的向量
        dIim_df=dIim_df(pqv);                                           %形成L的對角元素組成的向量

        H=spdiags(dIre_de,0,Lpqv,Lpqv);                                 %將形成的H向量的元素擴展對角矩陣
        N=spdiags(dIre_df,0,Lpqv,Lpqv);                                 %將形成的H向量的元素擴展對角矩陣
        J=spdiags(dIim_de,0,Lpqv,Lpqv);                                 %將形成的H向量的元素擴展對角矩陣
        L=spdiags(dIim_df,0,Lpqv,Lpqv);                                 %將形成的H向量的元素擴展對角矩陣
        R=spdiags(dV_de,0,Ltol,Ltol);                                   %將形成的H向量的元素擴展對角矩陣
        S=spdiags(dV_df,0,Ltol,Ltol);                                   %將形成的H向量的元素擴展對角矩陣

        K=diag(dIre_dQ(pv));                                            %將K擴展成對角陣
        M=diag(dIim_dQ(pv));                                            %將M擴展成對角陣

        %組成雅克比矩陣
% %         Jac0=[(G+H),(-B+N);(B+J),(G+L)];
% %         Jac1=[R(pv,pqv),S(pv,pqv)];
        Jac2=[Matzr1;K;Matzr1;M;Matzr2];
        Jac1=[(G+H)         (-B+N)
              (B+J)         (G+L)
              R(pv,pqv)     S(pv,pqv)];
        Jac=[Jac1,Jac2];
                                       
        %建立修正方程并求解
        dx=-(Jac\detb);

        %收斂性的判斷
        if max(abs(dx))<1e-4;
           break
        end
        %修正量的求解
        ebus(pqv)=ebus(pqv)+dx(1:Lpqv);
        fbus(pqv)=fbus(pqv)+dx(Ltol:2*Lpqv);
        Ubus=ebus+j*fbus;
        Sbus(pv)=Sbus(pv)+j*dx(2*Lpqv+1:end);
    end                                                               %迭代過程結束
    end
    et=etime(clock, t0);
    clear;