clear
load pvv
nm=max(hl_n(:,2:3));nn=max(nm);%網絡總節點數
bnn=size(hl_n);bn=bnn(1,1);%支路數，含對地支路
lnn=size(hl_l);ln=lnn(1,1);%負荷總節點數
gnn=size(hl_g);gn=gnn(1,1);%發電機節點數
pvn=size(hl_v);pv=pvn(1,1);pq=nn-pv-1;
%----------%求取導納矩陣
yy=zeros(nn+1,nn+1);YY=zeros(nn,nn);s=zeros(nn,nn);
j=sqrt(-1);
for m=1:bn
yy(hl_n(m,2)+1,abs(hl_n(m,3))+1)=1/(hl_n(m,4)+j*(hl_n(m,5)));
yy(abs(hl_n(m,3))+1,hl_n(m,2)+1)=yy(hl_n(m,2)+1,abs(hl_n(m,3))+1);
if hl_n(m,6)~=0
    yy(hl_n(m,2)+1,1)=yy(hl_n(m,2)+1,1)+(1-hl_n(m,6))*yy(hl_n(m,2)+1,abs(hl_n(m,3))+1)/hl_n(m,6)^2;
    yy(1,hl_n(m,2)+1)=yy(hl_n(m,2)+1,1);
    yy(abs(hl_n(m,3))+1,1)=yy(abs(hl_n(m,3))+1,1)+(-1+hl_n(m,6))*yy(hl_n(m,2)+1,abs(hl_n(m,3)))/hl_n(m,6);
    yy(1,abs(hl_n(m,3))+1)= yy(abs(hl_n(m,3))+1,1);
    yy(hl_n(m,2)+1,abs(hl_n(m,3))+1)=yy(hl_n(m,2)+1,abs(hl_n(m,3))+1)/hl_n(m,6);
end
end
for m=1:nn
    s=0;
    for n=1:nn
        YY(m,n)=-1*yy(m+1,n+1);
        s=s+yy(m+1,n);
    end
    YY(m,m)=s+yy(m+1,nn+1);
end
G=real(YY);B=imag(YY);
%-------------節點判斷與定值給定
v=ones(nn,1);o=zeros(nn,1);%電壓幅值/相角
e=ones(nn,1);f=zeros(nn,1);e(1)=1.06;%電壓實部/虛部
pg=zeros(nn,1);qg=zeros(nn,1);%發電機有功/無功
pl=zeros(nn,1);ql=zeros(nn,1);%負荷有功/無功
pn=zeros(nn,1);qn=zeros(nn,1);%網絡有功/無功
1*1==1;%1節點固定為平衡節點,考慮PV節點,
for l=1:pv
    s=hl_v(l,1);%pv節點號
    v(s)=hl_v(l,2);
end
for l=2:gn;
    s=hl_g(l,1);%PQ節點號，沒有考慮PV節點
    pg(s)=hl_g(l,2);%有功給定
    qg(s)=hl_g(l,3);%無功給定
end
 for l=1:ln;
     s=hl_l(l,1);%負荷節點號，
    pl(s)=hl_l(l,2);%有功給定
    ql(s)=hl_l(l,3);%無功給定
 end
%---------變量誤差
opq=ones(2*nn-2,1);
%------------------------------------------------------------------------------雅可比方程建立與迭代求解
h=0;
while(max(abs(opq))>0.0001);
    h=h+1;
    if h>100;
        disp('迭代次數大于100，退出');
        break;
    end
    %-------節點網絡功率
    for m=2:nn
        sump=0;
        sumq=0;
        a(m)=0;
        b(m)=0;
        for n=1:nn;
            sump=sump+(e(m)*(G(m,n)*e(n)-B(m,n)*f(n))+f(m)*(G(m,n)*f(n)+B(m,n)*e(n)));
            sumq=sumq+(f(m)*(G(m,n)*e(n)-B(m,n)*f(n))-e(m)*(G(m,n)*f(n)+B(m,n)*e(n)));
            a(m)=a(m)+(G(m,n)*e(n)-B(m,n)*f(n));
            b(m)=b(m)+(G(m,n)*f(n)+B(m,n)*e(n));
        end;
        pn(m)=sump;
        qn(m)=sumq;
    end;
  qn(pq+2:nn)=0;
   
  %------------------2(n-1)階雅可比矩陣
    for m=2:nn
        for n=2:nn;
            
            if m~=n
                Npv(m,n)=(G(m,n)*e(m)+B(m,n)*f(m));
                Lqv(m,n)=(G(m,n)*f(m)-B(m,n)*e(m));
                Hpo(m,n)=Lqv(m,n);
                Mqo(m,n)=-Npv(m,n);
                R(m,n)=0;
                S(m,n)=0;
            else
                Npv(m,m)= G(m,m)*e(m)+B(m,m)*f(m)+a(m);
                Lqv(m,m)=-B(m,m)*e(m)+G(m,m)*f(m)-b(m);
                Hpo(m,m)=-B(m,m)*e(m)+G(m,m)*f(m)+b(m);
                Mqo(m,m)=-G(m,m)*e(m)-B(m,m)*f(m)+a(m);
                R(m,n)=2*f(m);
                S(m,n)=2*e(m);
            end;
        end;
    end
        Hpo(1,:)=[];
        Hpo(:,1)=[];
        Npv(1,:)=[];
        Npv(:,1)=[];
        Mqo(1,:)=[];%Mqo已經刪去一行
        Mqo((pq+1):nn-1,:)=[];
        Mqo(:,1)=[];
        Lqv(1,:)=[];%Lqv已經刪去一行
        Lqv((pq+1):nn-1,:)=[];
        Lqv(:,1)=[];
        R(1:pq+1,:)=[];
        R(:,1)=[];
        S(1:pq+1,:)=[];
        S(:,1)=[];
  A=[ Hpo,Npv;Mqo, Lqv;R,S];%雅可比矩陣
        op=pg-pl-pn;
        oq=qg-ql-qn;
        ovv=v.*v-(e.*e+f.*f);
        op(1)=[];%刪除向量中與平衡節點對應的行
        oq(1)=[];
        oq(pq+1:nn-1)=[];%刪除向量中與pv節點對應的行
        ovv(1:pq+1)=[];%刪除向量中與pq節點,平衡節點對應的行
        opq=[op;oq;ovv];
        detx=inv(A)*opq;
        detf=detx(1:nn-1);
        dete=detx(nn:(2*nn-2));
        for i=2:nn
            f(i)=f(i)+detf(i-1);
            e(i)=e(i)+dete(i-1);
        end
end
    


       %----------------------------------------------------------------------------- %平衡節點功率
        sum=0;
        for m=1:nn
            u(m,1)=e(m)+j*f(m);
            sum=sum+YY(1,m)*u(m);
            v(m)=sqrt(e(m)^2+f(m)^2);
        end;
     s1=u(1)*conj(sum);
     pg(1)=real(s1);
     qg(1)=imag(s1);
        %支路功率，網損
        s=zeros(bn,bn);
  for m=1:nn
    for n=1:nn;
            s(m,n)=u(m)*conj(u(m)*yy(m+1,1)+[((u(m)-u(n))*yy(m+1,n+1))]);
     end;
           s(m,m)=0;
  end
  sh=0;th=0;lose=0;%'網絡總損耗'
  for i=1:nn
      sh=sh+pg(i)+j*qg(i);
      th=th+pl(i)+j*ql(i);
  end
  lose=sh-th     %'網絡總損耗'
  x=0;p=0;
  for i=1:nn
      x=x+pl(i);
  end
  for i=1:nn
      p=p+pg(i);
  end
 
       %----------------------輸出
        disp('-------------------------------計算結果---------------------');
        disp('1-節點/2-電壓值/3-負荷功率/4-發電機功率/');
        j=sqrt(-1);
        out=[[1:nn]',e+j*f,pl+j*ql,pg+j*qg]
        disp('平衡節點功率/');
        out1=s1
        disp('線路功率損耗');
        out=s
        disp('網絡總損耗');
        out=lose
         disp('輸電效率');
         out=x/p