main.m

利用该MATLAB程序可实现对一般电力系统的节点、支路及其损耗的潮流计算。

[dfile,pathname]=uigetfile('test11.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
global n;
global m;
[nb,mb]=size(bus);   %节点重新编号
[nl,ml]=size(line);
  nSW=0;                 % 平衡节点数目
  nPV=0;                 % PV节点数目
  nPQ=0;                 % PQ节点数目
  for i=1:nb,            % nb为总节点数
    type=bus(i,6);
    if type==3,
       nSW=nSW+1;     % 统计平衡节点数目
       SW(nSW,:)=bus(i,:);
    elseif type==2,
        nPV=nPV+1;    % 统计PV节点数目
       PV(nPV,:)=bus(i,:);
    else
        nPQ=nPQ+1;    % 统计PQ节点数目
       PQ(nPQ,:)=bus(i,:);
    end
  end;
  bus=[PQ;PV;SW];
  newbus=[1:nb]';
  f=bus(:,1);
  nodenum=[newbus bus(:,1)]; 
  bus(:,1)=newbus;
  for i=1:nl
      for j=1:2
          for k=1:nb
              if line(i,j)==nodenum(k,2)
                 line(i,j)=nodenum(k,1);
                 break
              end
          end
      end
  end
Y=y(bus,line); %形成节点导纳矩阵
K=0;     %迭代次数初值
Kmax=20;  %最大迭代次数
eps1=1.0e-10;
eps2=1.0e-10;
m=nPQ;n=nb;
Um=eye(m,m);
myf=fopen('output1.dat','w');
for K=1:Kmax
    for i=1:m
        for j=1:m
          if i==j
             Um(i,j)=bus(i,2);
          end
        end
    end
     b=dPQ(Y,bus);
     C=jac(bus,Y);
     dX=C\b';
     dx=dX';
     [nx,mx]=size(dx);
     for i=1:n-1     %计算相角
         bus(i,3)=bus(i,3)-dX(i,1); 
     end
     B=dx(nx,n:mx)*Um;  %计算电压差
     bus(1:m,2)=bus(1:m,2)-B';  %计算电压值
     dx(nx,n:mx)=B;     
        fprintf(myf,'--第%d次迭代时雅可比矩阵--',K)
         fprintf(myf, '\n');
            for i=1:(n+m-1)
                for j=1:(n+m-1)
                    fprintf(myf,'%8.6f', C(i,j));
                    fprintf(myf, '  ');
                end
                fprintf(myf, '\n');
            end    
        fprintf(myf,'--第%d次迭代时dPQ的误差--',K)
         fprintf(myf, '\n');
            for i=1:(n+m-1)
                fprintf(myf,'%8.6e', b(1,i));
                 fprintf(myf, '\n');
            end
             fprintf(myf, '\n');
           fprintf(myf,'--第%d次迭代时dx(误差)--',K) 
            fprintf(myf, '\n');
            for i=1:(n+m-1)
               fprintf(myf,'%8.6e', dX(i,1));       
                 fprintf(myf, '\n');
            end 
            fprintf(myf, '\n');
            fprintf(myf,'第%d次迭代后节点电压(仅PQ节点)',K)
            fprintf(myf, '\n');
             for i=1:m
               fprintf(myf,'%8.6f', bus(i,2));
                fprintf(myf, '\n');
             end
            fprintf(myf, '\n');
            fprintf(myf,'第%d次迭代后相角(角度)',K)
            fprintf(myf, '\n');
             for i=1:n
               fprintf(myf,'%8.6f', bus(i,3)*180/pi);
                fprintf(myf, '\n');
             end
            fprintf(myf, '\n');              
      if (max(abs(dx))<eps1)&(max(abs(b))<eps2)   %判断是否达到计算精度          
          break;
      end
end
%计算功率
P1=0;T=0;    %计算平衡节点的有功
for j=1:n
     T=bus(n,2)*bus(j,2)*(real(Y(n,j))*cos(bus(n,3)-bus(j,3))+imag(Y(n,j))*sin(bus(n,3)-bus(j,3)));
     P1=P1+T;
end
   bus(n,4)=P1;
for k=m+1:n   %计算PV节点、平衡节点的无功
    Q1=0;T=0;
    for j=1:n
        T=bus(k,2)*bus(j,2)*(real(Y(k,j))*sin(bus(k,3)-bus(j,3))-imag(Y(k,j))*cos(bus(k,3)-bus(j,3)));
        Q1=Q1+T;
    end
    bus(k,5)=Q1;
end
bus(:,1)=f;  %换回各节点、支路的初始编号
r=zeros(1,mb);
for t=1:n
    for l=t+1:n
        if bus(t,1)>bus(l,1)
           r=bus(t,:);bus(t,:)=bus(l,:);bus(l,:)=r;
        end
    end
end
for i=1:nl
      for j=1:2
          for k=1:nb
              if line(i,j)==nodenum(k,1)
                 line(i,j)=nodenum(k,2);
                 break
              end
          end
      end
end
Pf=loss(bus,line);  %计算支路潮流及损耗
    %将节点导纳矩阵、节点潮流计算结果写入文件output3
     fprintf(myf, '--节点导纳矩阵--\n');
     for k=1:n
       for j=1:n
       fprintf(myf,'%8.6f', real(Y(k,j)));
       fprintf(myf, '+i*');
       fprintf(myf,'%8.6f', imag(Y(k,j)));
       fprintf(myf, '  ');
     end
     fprintf(myf, '\n');
     end
  fprintf(myf, '--节点计算结果--\n');
  fprintf(myf, '--   bus(#) (volt) (ang（角度）)  (P)  (Q)   (bus type)--\n');
   for l=1:nb
     for j=1:mb
         if j==1|j==6
        fprintf(myf, '%8.1f ', bus(l,j));
         elseif j==3
        fprintf(myf, '%8.6f ', bus(l,j)*180/pi);
         else
        fprintf(myf, '%8.6f ', bus(l,j));    
         end
     end
     fprintf(myf, '\n');
   end 
   fprintf(myf, '--支路计算结果--\n');
   fprintf(myf, '--     (I)    (J)     Real(S(I,J))Imag(S(I,J)) Real(S(J,I)) Imag(S(J,I)) Real(delS(I,J))Imag(delS(I,J))--\n');
   for k=1:nl
     for j=1:5
         if j<=2
       fprintf(myf,'%8.1f', Pf(k,j));
       fprintf(myf, '   ');
         else
       fprintf(myf,'%8.6f', real(Pf(k,j)));
       fprintf(myf, '+i*');
       fprintf(myf,'%8.6f', imag(Pf(k,j)));
       fprintf(myf, '  ');
         end
     end
     fprintf(myf, '\n');
   end
   fclose(myf);