zhuchengxu.m

计算机潮流计算 应用的是牛拉法

clear all;
clc;
[dfile,pathname]=uigetfile('*.m','Select Data File');
if pathname == 0
  error(' you must select a valid data file')
else
  lfile =length(dfile);
  % strip off .m
  eval(dfile(1:lfile-2));
end                                  %读入数据


[nb,mb]=size(bus);         %重新编号
   [nl,ml]=size(line);
  nSW = 0;                 % number of swing bus counter
  nPV = 0;                 % number of PV bus counter
  nPQ = 0;                 % number of PQ bus counter
  for e = 1:nb,              % nb为总节点数
    type= bus(e,6);
    if type == 3,
       nSW = nSW + 1;     % increment swing bus counter
       SW(nSW,:)=bus(e,:);
    elseif type == 2,
        nPV = nPV +1;      % increment PV bus counter
       PV(nPV,:)=bus(e,:);
    else
        nPQ = nPQ + 1;     % increment PQ bus counter
       PQ(nPQ,:)=bus(e,:);
    end
  end; 
bus=[PQ;PV;SW];
  newbus=[1:nb]';
  nodenum=[newbus bus(:,1)];
  bus(:,1)=newbus; 
  for e=1:nl
      for l=1:2
          for k=1:nb
              if line(e,l)==nodenum(k,2)
                 line(e,l)=nodenum(k,1);
                 break
              end
          end
      end
  end 

  
  


[n,l]=size(bus);
m=0;
for k=1:n
    if bus(k,l)==1
        m=m+1;
    end
end
myf=fopen('D:\Downloads\matlab\work\output.txt','w');
fprintf(myf,'--------------节点导纳矩阵----------\n');
Y = y(bus,line);
[u,v]=size(Y);


for k=1:u,  
    for l=1:v,
        fprintf(myf, '%11.6f', real(Y(k,l)));
        fprintf(myf,'+j%11.6f  ',imag(Y(k,l)));
    end
    fprintf(myf,'\n');
end


r=1;
for e=1:10
    x=dPQ(Y,bus,n,m);
    dP=x(1:(n-1),1);
    dQ=x(n:(n-1+m),1);
    J=form_jac(bus,Y,n,m);
    fprintf(myf,'--------------第%d次迭代的结果----------\n',r);
    fprintf(myf,'-------------第%d次迭代的雅可比矩阵J----------\n',r);
    [u,v]=size(J);
for k=1:u
    for l=1:v
        fprintf(myf,'%11.6f',J(k,l));
    end
    fprintf(myf,'\n');
end
    fprintf(myf,'-------------第%d次迭代的功率偏差dP和dQ----------\n',r);
    for k=1:(n-1)
        fprintf(myf,'dP%d  %.6e\n',k,dP(k,1));
    end
    for k=1:m
        fprintf(myf,'dQ%d  %.6e\n',k,dQ(k,1));
    end
    U=zeros(m);
    for k=1:m
        U(k,k)=bus(k,2);
    end
    dang=zeros((n-1),1);
    dU=zeros(m,1);
    W=zeros((n-1+m),1);
    W=inv(J)*[dP;dQ];
    dang(1:(n-1),1)=W(1:(n-1),1);
    dU(1:m,1)=W(n:(n+m-1),1);
    dU=U*dU;
    fprintf(myf,'-------------第%d次迭代的节点相角和电压的偏差dx----------\n',r);
    for k=1:(n-1)
        fprintf(myf,'dang%d  %.6e\n',k,dang(k,1));
    end
    for k=1:m
        fprintf(myf,'dU%d  %.6e\n',k,dU(k,1));
    end
    bus(1:m,2)=bus(1:m,2)-dU;
    bus(1:(n-1),3)=bus(1:(n-1),3)-dang;
    fprintf(myf,'-------------第%d次迭代的节点相角delta(弧度为单位）和电压U----------\n',r);
    for k=1:(n-1)
        fprintf(myf,'ang%d  %f\n',k,bus(k,3));
    end
    for k=1:m
        fprintf(myf,'U%d  %f\n',k,bus(k,2));
    end
    r=r+1;
    if (max(abs(dU))<(10^(-10)))&(max(abs(dang))<(10^(-10)))
        break
    end
end




for k=(m+1):n
    for e=1:n
    bus(k,5)=bus(k,5)+bus(k,2)*bus(e,2)*(real(Y(k,e))*sin(bus(k,3)-bus(e,3))-imag(Y(k,e))*cos(bus(k,3)-bus(e,3)));
    end
end
for e=1:n
    bus(n,4)=bus(n,4)+bus(k,2)*bus(e,2)*(real(Y(k,e))*cos(bus(k,3)-bus(e,3))+imag(Y(k,e))*sin(bus(k,3)-bus(e,3)));
end

bus(:,1)=nodenum(:,2);                                  %改回编号。。。。。。。。
nbus=zeros(nb,mb);
for k=1:nb
    nbus(bus(k,1),:)=bus(k,:);
end
bus=nbus;

for e=1:nl
      for t=1:2
          for k=1:nb
              if line(e,t)==nodenum(k,1)
                 line(e,t)=nodenum(k,2);
                 break;
              end
          end
      end
end 
  

fprintf(myf,'---------------牛顿－拉夫逊法潮流计算结果----------\n');
fprintf(myf,'节点计算结果：\n');
fprintf(myf,'节点   节点电压    节点相角（角度）    节点注入功率\n');

for k=1:n
    fprintf(myf,'%2.d   %11.6f%11.6f%11.6f+j%11.6f\n',bus(k,1),bus(k,2),bus(k,3)*180/pi,bus(k,4),bus(k,5));
end


fprintf(myf,'线路计算结果：\n');
    fprintf(myf,'节点I    节点J          线路功率S(I,J)         线路功率S(J,I)         线路损耗dS(I,J)\n');
    [nl,ml]=size(line);
    for k=1:nl
        I=line(k,1);                       %读入线路参数
        J=line(k,2);
        Zt=line(k,3)+j*line(k,4);
        if Zt==0
           Yt=10^10;
        else
            Yt=1/Zt;
        end
        Ym=line(k,5)+j*line(k,6);
        K=line(k,7);
        if (K==0)&(J~=0)     % 普通线路: K=0;
            S(I,J)=(bus(I,2)^2)*(conj(Yt)+conj(Ym))-(bus(I,2)*cos(bus(I,3))+j*bus(I,2)*sin(bus(I,3)))*conj(bus(J,2)*cos(bus(J,3))+j*bus(J,2)*sin(bus(J,3)))*conj(Yt);
            S(J,I)=(bus(J,2)^2)*(conj(Yt)+conj(Ym))-(bus(J,2)*cos(bus(J,3))+j*bus(J,2)*sin(bus(J,3)))*conj(bus(I,2)*cos(bus(I,3))+j*bus(I,2)*sin(bus(I,3)))*conj(Yt);
            dS(I,J)=S(I,J)+S(J,I);
            fprintf(myf,'%2.d       %2.d       %f+j%f       %f+j%f       %f+j%f\n',I,J,real(S(I,J)),imag(S(I,J)),real(S(J,I)),imag(S(J,I)),real(dS(I,J)),imag(dS(I,J)));
        end

        if (K==0)&(J==0)               % 对地支路: K=0,J=0,R=X=0;
            S(I,I)=(bus(I,2)^2)*conj(Ym);
        fprintf(myf,'%2.d       0.000000       %f+j%f           0.000000                      %f+j%f\n',I,real(S(I,I)),imag(S(I,I)),real(S(I,I)),imag(S(I,I)));
    end

    if K<0                      % 变压器线路: Zt和Ym为折算到K侧的值,K在i侧
        S(I,J)=bus(I,2)^2*(conj(-K*Yt)+conj(Ym+(1+K)*Yt))-(bus(I,2)*cos(bus(I,3))+j*bus(I,2)*sin(bus(I,3)))*conj(bus(J,2)*cos(bus(J,3))+j*bus(J,2)*sin(bus(J,3)))*conj(Yt*(-K));
        S(J,I)=bus(J,2)^2*(conj(-K*Yt)+conj(Ym+(K^2+K)*Yt))-(bus(J,2)*cos(bus(J,3))+j*bus(J,2)*sin(bus(J,3)))*conj(bus(I,2)*cos(bus(I,3))+j*bus(I,2)*sin(bus(I,3)))*conj(Yt*(-K));
        dS(I,J)=S(I,J)+S(J,I);
        fprintf(myf,'%2.d       %2.d       %f+j%f       %f+j%f       %f+j%f\n',I,J,real(S(I,J)),imag(S(I,J)),real(S(J,I)),imag(S(J,I)),real(dS(I,J)),imag(dS(I,J)));
    end
    if K>0                 % 变压器线路: Zt和Ym为折算到i侧的值,K在j侧
        S(I,J)=bus(I,2)^2*(conj(Yt/K)+conj(Ym+Yt*(K-1)/K))-(bus(I,2)*cos(bus(I,3))+j*bus(I,2)*sin(bus(I,3)))*conj(bus(J,2)*cos(bus(J,3))+j*bus(J,2)*sin(bus(J,3)))*conj(Yt/K);
        S(J,I)=bus(J,2)^2*(conj(Yt/K)+conj(Ym+Yt*(1-K)/K/K))-(bus(J,2)*cos(bus(J,3))+j*bus(J,2)*sin(bus(J,3)))*conj(bus(I,2)*cos(bus(I,3))+j*bus(I,2)*sin(bus(I,3)))*conj(Yt/K);
         dS(I,J)=S(I,J)+S(J,I);
        fprintf(myf,'%2.d       %2.d       %f+j%f       %f+j%f       %f+j%f\n',I,J,real(S(I,J)),imag(S(I,J)),real(S(J,I)),imag(S(J,I)),real(dS(I,J)),imag(dS(I,J)));
    end 
    end
    fclose(myf);