pvv.asv

1节点固定为平衡节点

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;
      
        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)));
        
        end;
        pn(m)=sump;
        qn(m)=sumq;
    end;
  p(1,1)=[];
  q(1,1)=[];
  q(pq+1:nn-1)=[];
  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