niudunchaoliu).txt

能进行牛顿法的电力电力系统的潮流计算而且它的收敛性很好已经进行了测试

%牛顿——拉夫逊法进行潮流计算 
%n=input('请输入节点数：n='); 
%nl=input('请输入支路数：nl='); 
%isb=input('请输入平衡母线节点号：isb='); 
%pr=input('请输入误差精度：pr='); 
%iternu=input('请输入迭代次数:iternu='); 
%B1=input('请输入由支路参数形成的矩阵： B1='); 
%B2=input('请输入各节点参数形成的矩阵： B2='); 
%以下为从数据文件读入数据 
clc %清屏 
clear %从内存中清除变量和函数 
fid = fopen('flowre1.dat','r'); 
n = fscanf(fid,'%d',[1,1]) %节点数 
n1= fscanf(fid,'%d',[1,1]) %支路数 
isb=fscanf(fid,'%d',[1,1]) %平衡母线节点号 
pr=fscanf(fid,'%g',[1,1])  %误差精度 
iternu=fscanf(fid,'%d',[1,1])%最大迭代次数 
B=fscanf(fid,'%g',[7,n]);  %支路数据 
BP=fscanf(fid,'%g',[8,n1]);  %节点数据 
B=B' 
BP=BP' 
B1=zeros(n1,7);B2=zeros(n,8);%B1支路参数形成的矩阵 B2各节点参数形成的矩阵 
for i=1:n1 
    k=1; 
    for z=1:6; 
      if (z==3)      
          B1(i,z)=complex(B(i,z),B(i,z+1));k=z+1; 
      else  
          if (z==4)  
              B1(i,z)=complex(0,B(i,k));  
          else 
              B1(i,z)=B(i,k); 
          end      
      end  
      k=k+1;     
  end   
end  
for i=1:n 
    k=1; 
    for z=1:6; 
      if (z==1)   
          B2(i,z)=complex(BP(i,z),BP(i,z+1));k=z+2; 
           
      else  
          if (z==2) 
              B2(i,z)=complex(BP(i,k),BP(i,k+1));k=k+2; 
          else 
              B2(i,z)=BP(i,k);   k=k+1;    
          end 
      end  
        
  end   
end  
fclose(fid);%关闭数据原文件 
Y=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);%初始化数组 
delta=zeros(1,n);S1=zeros(n1); 
%求导纳矩阵 
for i=1:n1 
    if B1(i,6)==0 
       p=B1(i,1);q=B1(i,2); 
    else 
       p=B1(i,2);q=B1(i,1); 
    end 
    Y(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5)); 
    Y(q,p)=Y(p,q); 
    Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2; 
    Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2; 
end 
disp('导纳矩阵 Y='); 
disp(Y) 
G=real(Y);B=imag(Y); 
fid=fopen('result.dat','a'); 
fprintf(fid,'\n导纳矩阵G 实部') 
fprintf(fid,'\n%g %g %g %g %g',G) 
fprintf(fid,'\n导纳矩阵B 虚部') 
fprintf(fid,'\n%g %g %g %g %g',B) 
 
for i=1:n 
    e(i)=real(B2(i,3)); 
    f(i)=imag(B2(i,3)); 
    V(i)=B2(i,4); 
end 
for i=1:n 
    S(i)=B2(i,1)-B2(i,2);%节点注入功率＝发电机注入节点的功率减去节点负荷功率 
    B(i,i)=B(i,i)+B2(i,5); 
end 
P=real(S);Q=imag(S); 
ICT1=0;IT2=1;N0=2*n;N=N0+1;a=0;overflag=0; 
while IT2~=0 
      IT2=0;a=a+1; 
       
      for i=1:n 
          if i~=isb%除平衡节点外 
             C(i)=0;D(i)=0; 
             for j1=1:n 
                 C(i)=C(i)+G(i,j1)*e(j1)-B(i,j1)*f(j1); 
                 D(i)=D(i)+G(i,j1)*f(j1)+B(i,j1)*e(j1); 
             end 
             P1=C(i)*e(i)+f(i)*D(i); 
             Q1=C(i)*f(i)-e(i)*D(i); 
%求'P,Q'     
             V2=e(i)^2+f(i)^2; 
             if B2(i,6)~=3 
                DP=P(i)-P1;%有功功率的误差 
                DQ=Q(i)-Q1;%无功功率的误差 
                for j1=1:n 
                    if j1~=isb&j1~=i 
                       X1=-G(i,j1)*e(i)-B(i,j1)*f(i); 
                       X2=B(i,j1)*e(i)-G(i,j1)*f(i); 
                       X3=X2; 
                       X4=-X1; 
                       p=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ; 
                       m=p+1; 
                       J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X4; 
                       J(m,q)=X2; 
                    elseif j1==i&j1~=isb 
                       X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i); 
                       X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i); 
                       X3=D(i)+B(i,i)*e(i)-G(i,i)*f(i); 
                       X4=-C(i)+G(i,i)*e(i)+B(i,i)*f(i); 
                       p=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ; 
                       m=p+1; 
                       J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X4; 
                       J(m,q)=X2; 
                    end 
                end 
             else 
                DP=P(i)-P1; 
                DV=V(i)^2-V2; 
                for j1=1:n 
                    if j1~=isb&j1~=i 
                       X1=-G(i,j1)*e(i)-B(i,j1)*f(i); 
                       X2=B(i,j1)*e(i)-G(i,j1)*f(i); 
                       X5=0;X6=0; 
                       p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV; 
                       m=p+1; 
                       J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6; 
                       J(m,q)=X2; 
                    elseif j1==i&j1~=isb 
                       X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i); 
                       X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i); 
                       X5=-2*e(i); 
                       X6=-2*f(i); 
                       p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV; 
                       m=p+1; 
                       J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6; 
                       J(m,q)=X2; 
                    end 
                end 
             end 
          end 
      end 
%求雅可比矩阵 
      for k=3:N0 
          k1=k+1;N1=N; 
          for k2=k1:N1 
              J(k,k2)=J(k,k2)./J(k,k); 
          end 
          J(k,k)=1; 
          if k~=3 
             k4=k-1; 
             for k3=3:k4 
                 for k2=k1:N1 
                  J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2); 
                 end 
                 J(k3,k)=0; 
             end 
             if k==N0  
                 break; 
             end 
                for k3=k1:N0 
                    for k2=k1:N1 
                        J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2); 
                    end 
                J(k3,k)=0; 
             end 
          else 
             for k3=k1:N0 
                 for k2=k1:N1 
                     J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2); 
                 end 
                 J(k3,k)=0; 
             end 
          end 
      end 
  
      for k=3:2:N0-1 
          L=(k+1)./2; 
          e(L)=e(L)-J(k,N); 
          k1=k+1; 
          f(L)=f(L)-J(k1,N); 
      end 
      for k=3:N0 
          DET=abs(J(k,N)); 
          if DET>=pr 
             IT2=IT2+1; 
          end 
      end 
      ICT2(a)=IT2; 
      ICT1=ICT1+1; 
      if ICT1>iternu  
          overflag=1; 
          break;     
      end  
end 
 
%用高斯消去法解"w=-J*V" 
if overflag==0 
disp('迭代次数：'); 
disp(ICT1); 
disp('没有达到精度要求的个数：'); 
disp(ICT2); 
for k=1:n 
    V(k)=sqrt(e(k)^2+f(k)^2); 
    delta(k)=atan(f(k)./e(k))*180./pi; 
    E(k)=e(k)+f(k)*j; 
end 
disp('各节点的实际电压标幺值E为(节点号从小到大排列)：'); 
disp(E); 
disp('各节点的电压大小V为(节点号从小到大排列)：'); 
disp(V); 
fprintf(fid,'\n各节点的电压大小V为(节点号从小到大排列):'); 
fprintf(fid,'\n%g %g %g %g %g ',V); 
 
disp('各节点的电压相角delta为(节点号从小到大排列)：'); 
disp(delta); 
fprintf(fid,'\n各节点的电压相角delta为(节点号从小到大排列)：'); 
fprintf(fid,'\n%g %g %g %g %g ',delta); 
 
for p=1:n 
    C(p)=0; 
    for q=1:n 
        C(p)=C(p)+conj(Y(p,q))*conj(E(q)); 
    end 
    S(p)=E(p)*C(p); 
end 
disp('各节点的功率S为(节点号从小到大排列)：'); 
disp(S); 
disp('各条支路的首端功率Si为(顺序同您输入B1时一致)：'); 
for i=1:n1 
    if B1(i,6)==0 
       p=B1(i,1);q=B1(i,2); 
    else 
       p=B1(i,2);q=B1(i,1); 
    end 
Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)*B1(i,5))-conj(E(q)))*conj(1./(B1(i,3)*B1(i,5)))); 
    disp(Si(p,q)); 
end 
disp('各条支路的末端功率Sj为(顺序同您输入B1时一致)：'); 
for i=1:n1 
    if B1(i,6)==0 
       p=B1(i,1);q=B1(i,2); 
    else 
       p=B1(i,2);q=B1(i,1); 
    end 
Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)./B1(i,5))-conj(E(p)))*conj(1./(B1(i,3)*B1(i,5)))); 
    disp(Sj(q,p)); 
end 
disp('各条支路的功率损耗DS为(顺序同您输入B1时一致)：'); 
for i=1:n1 
    if B1(i,6)==0 
       p=B1(i,1);q=B1(i,2); 
    else 
       p=B1(i,2);q=B1(i,1); 
    end 
    DS(i)=Si(p,q)+Sj(q,p); 
    disp(DS(i)); 
end 
fclose(fid); 
else  
   disp('在最大迭代次数范围内不收敛');  
   fclose(fid); 
end  
disp(J)