clc;
% % t0 = clock;
% % for k=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] = test30_1;
[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);
%迭代計算
iter=100;
while (iter)
    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向量的元素擴展對角矩陣
    
    R=R(pv,pqv);                                                    %取出PV節點對應R
    S=S(pv,pqv);                                                    %取出PV節點對應S
    
    K=diag(dIre_dQ(pv));                                            %將K擴展成對角陣
    M=diag(dIim_dQ(pv));                                            %將M擴展成對角陣
    
    
    H=G+H;
    N=-B+N;
    J=B+J;
    L=G+L;
    
    %組成雅克比矩陣
    Jac0=[H,N;J,L];
    Jac1=[R,S];
    Matzr1=zeros(Lpq,Lpv);
    Matzr2=zeros(Lpv);
    Jac2=[Matzr1;K;Matzr1;M;Matzr2];
    
    Jac=[[Jac0;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);
    iter=iter-1;
end                                                                %迭代過程結束
% % end
% % et = etime(clock, t0);
% % % a=abs(Ubus);
Ubus;                                                              %得到各節點電壓值
Ss=0;
slack=bus(find(bus(:,2)==3),1);
for k=1:Ltol                                                       %得到平衡節點得注入功率
    Ss=Ss+Ubus(slack)*conj(Ubus(k))*conj(Ybus(slack,k));
end
disp('平衡節點注入功率Ss=')
disp(Ss)
%計算線路首端和末端的有功功率和無功功率