niudun_lafuxun.m

电力系统计算节点数不多的潮流程序,数度比较快

[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 and convert to lower case
  dfile = lower(dfile(1:lfile-2));
  eval(dfile);
end
[nb,a]=size(bus);
[nl,c]=size(line);
[bus,line,nodenum]=newform(bus,line,nb,nl);            
Y=ybus(line,nb,nl);                                       
myf=fopen('test2.dat','w');
fprintf(myf,'－－－－－－－节点导纳矩阵Y：－－－－－－－－－\n')
for i=1:nb,
    for j=1:nb,
       fprintf(myf,' %11.6f ',Y(i,j))  ;
    end
    fprintf(myf,'\n') ; 
end
 

format long e
max1=100;
eps1=1.0e-10;
eps2=1.0e-10;
V=bus(:,2);
ang=bus(:,3);
ang1=[];
V1=[];

m=0;
for i=1:nb,
    if bus(i,6)==1
        m=m+1;
    end
end
n=nb-1-m; 

for i=1:nb-1,
    ang1(i,1)=bus(i,3);
    if bus(i,6)==1
        V1(i,1)=bus(i,2);
    end
end
angV=[ang1;V1]; 

for i=1:max1,
    A=form_jac(V,ang,Y,nb,m);                   
     b=form_ds(V,ang,Y,bus,nb,m);                   
    VdAV=gaussj(A,b);                           
    dAV=VdAV;
    for j=nb:nb-1+m,
        dAV(j,1)=VdAV(j,1)*V1(j-nb+1,1);
    end
    angV=angV-dAV;
    for x=1:nb-1,
        ang(x,1)=angV(x,1);
    end
    for y=nb:nb-1+m,
        V(y-nb+1,1)=angV(y,1);
    end
    V1=V;
  
  
  fprintf(myf,'第%d次迭代结果：\n',i)
    
  fprintf(myf,'－－－－－－－第%d次迭代的雅可比矩阵J－－－－－－－－－\n',i)
    
  for o=1:nb-1+m,
     for l=1:nb-1+m,
        fprintf(myf,' %11.6f ',A(o,l))  ;
     end
     fprintf(myf,'\n') ; 
  end

  fprintf(myf,'第%d次迭代的功率偏差dP和dQ\n',i);
    	for z=1:nb-1,
	           fprintf(myf, ' dP%d   %u\n  ',i,b(z));
        end
        for z=1:m,   
            fprintf(myf, ' dQ%d   %u\n   ',i,b(nb-1+z));
        end
     
            
            
  fprintf(myf,'第%d次迭代的节点相角和电压的偏差dA和dV\n',i);
    	for z=1:nb-1,
	           fprintf(myf, ' dA%d   %u\n    ',i,dAV(z));
           end
        for z=1:m,   
            fprintf(myf, ' dV%d   %u\n    ',i,dAV(nb-1+z));
        end
            fprintf(myf, '\n');    
            
   fprintf(myf,'第%d次迭代的节点相角和电压A和V\n',i);
    	for z=1:nb-1,
	           fprintf(myf, ' dA%d   %11.6f\n   ',i,angV(z));
        end
        for z=1:m,   
            fprintf(myf, ' dV%d   %11.6f\n   ',i,angV(nb-1+z));
        end
            fprintf(myf, '\n');           
    if (abs(dAV)<eps1&abs(b)<eps2)
        break
    end
end

P=zeros(nb,1);
Q=zeros(nb,1);
[P,Q]=form_s(V,ang,Y,nb,m,bus);
bus(:,2)=V;
bus(:,3)=ang;
bus(:,4)=P;
bus(:,5)=Q;
for i=1:nb,
    for j=1:nb,
       if nodenum(j,2)==bus(i,1),
           newbus(i,:)=bus(j,:);
           break
       end
   end
end
V=newbus(:,2);
ang=newbus(:,3);
P=newbus(:,4);
Q=newbus(:,5);
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

S=zeros(nb);
SS=zeros(nb);
delS=zeros(nb);
[S,SS,delS]=node(V,line,ang,nb); 
fprintf(myf, '--------------牛顿拉夫逊法潮流计算结果:---------------\n');
   fprintf(myf,'--------------节点计算结果-----------\n');
   fprintf(myf,'节点   节点电压    节点相角（角度）    节点注入功率\n');
    	for i=1:nb,
	        fprintf(myf, ' %d   %11.6f     %11.6f ',i,V(i),ang(i)*180/pi);
            fprintf(myf, '%11.6f+j%11.6f',P(i),Q(i));
            fprintf(myf, '\n');
        end
 	fprintf(myf, '--------------线路计算结果:-------------\n', b(i));
    fprintf(myf,'节点I   节点J        线路功率S(I,J)          线路功率S(J,I)            线路损耗dS(I,J)\n');
        for i=1:nl
            if line(i,2)==0
                line(i,2)=line(i,1);
            end
            fprintf(myf, '%d       %d      ',line(i,1),line(i,2));
            fprintf(myf, '%11.6f+j%11.6f    ',real(S(line(i,1),line(i,2))),imag(S(line(i,1),line(i,2))));
            fprintf(myf, '%11.6f+j%11.6f    ',real(SS(line(i,2),line(i,1))),imag(SS(line(i,2),line(i,1))));
            fprintf(myf, '%11.6f+j%11.6f    ',real(delS(line(i,1),line(i,2))),imag(delS(line(i,1),line(i,2))));
            fprintf(myf, '\n');
        end
    fclose(myf);