basep_qpfc.m

运用P-Q分解法自编的潮流计算程序,可计算简单电力系统的潮流

%请输入节点数：n=5
%请输入支路数：nl=5
%请输入平衡母线节点号：isb=1
%请输入误差精度：pr=0.0001
%请输入由支路参数形成的矩阵：B1=[1  2  0.03i      0    1.05 0
%                              2  3  0.08+0.3i  0.5i 1    0
%                              2  4  0.1+0.35i  0    1    0
%                              3  4  0.04+0.25i 0.5i 1    0
%                              3  5  0.015i     0    1.05 1]
%第一列和第二列是节点编号，编号由小到大排列。
%第三列为支路串联阻抗参数，第四列为支路对地导纳参数。
%第五列为含有变压器参数支路变压器变比，第六列为是否含有变压器的参数，其中1为含有变压器，0为不含变压器
%请输入各节点参数形成的矩阵：B2=[0  0         1.05   1.05  0  1
%                              0  3.7+1.3i  1.05   0     0  2
%                              0  2+1i      1      0     0  2
%                              0  1.6+0.8i  1.05   0     0  2
%                              5  0         1.05   1.05  0  3]
%第一列为节点注入发电功率参数，第二列为节点负荷功率参数
%第三列为节点电压参数，
%第六列为节点类型参数，1为平衡节点，2为PQ节点，3为PV节点
%请输入P-Q节点数：na=3


%%P-Q分解法潮流计算
n=input('请输入节点数：n=');
nl=input('请输入支路数：nl=');
isb=input('请输入平衡母线节点号：isb=');
pr=input('请输入误差精度：pr=');
B1=input('请输入由支路参数形成的矩阵：B1=');
B2=input('请输入各节点参数形成的矩阵：B2=');
na=input('请输入P-Q节点数：na=');
Y=zeros(n);YI=zeros(n);%YI是不加接地支路的导纳矩阵，Y是追加接地支路的导纳矩阵
e=zeros(1,n);f=zeros(1,n);%
V=zeros(1,n);O=zeros(1,n);%
for i=1:nl
      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));
      YI(p,q)=YI(p,q)-1./(B1(i,3));
      Y(q,p)=Y(p,q);
      YI(q,p)=YI(p,q);
      Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2;
      YI(q,q)=YI(q,q)+1./(B1(i,3));
      Y(p,p)=Y(p,p)+1./(B1(i,3))+B1(i,4)./2;
      YI(p,p)=YI(p,p)+1./(B1(i,3));
end
%求导纳矩阵
G=real(Y);B=imag(YI);BI=imag(Y);
for i=1:n   %计算功率
      S(i)=B2(i,1)-B2(i,2);
      BI(i,i)=BI(i,i)+B2(i,5);
end
P=real(S);Q=imag(S);
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          %计算PQ节点电压幅值和角度
      if B2(i,6)==2
      V(i)=sqrt(e(i)^2+f(i)^2);
      O(i)=atan(f(i)./e(i));
      end
end
for i=2:n                %形成导纳矩阵因子表
      if i==n
          B(i,i)=1./B(i,i);
      else
          IC1=i+1;
          for j1=IC1:n
              B(i,j1)=B(i,j1)./B(i,i);
          end
          B(i,i)=1./B(i,i);
          for k=i+1:n
              for j1=i+1:n
                  B(k,j1)=B(k,j1)-B(k,i)*B(i,j1);
              end
          end
      end
end
p=0;q=0;%p，q重新赋初值0
for i=1:n
      if B2(i,6)==2
          p=p+1;k=0;
          for j1=1:n
               if B2(j1,6)==2
                   k=k+1;
                   A(p,k)=BI(i,j1);
               end
          end
      end
end
for i=1:na
      if i==na
          A(i,i)=1./A(i,i);
      else
          k=i+1;
          for j1=k:na
              A(i,j1)=A(i,j1)./A(i,i);
          end
           A(i,i)=1./A(i,i);
           for k=i+1:na
               for j1=i+1:na
                A(k,j1)=A(k,j1)-A(k,i)*A(i,j1);
               end
           end
      end
end
ICT2=1;ICT1=0;kp=1;kq=1;K=1;DET=0;ICT3=1;  %        
while ICT2~=0|ICT3~=0                     %~=表示不等于，|表示或
      ICT2=0;ICT3=0;
      for i=1:n
          if i~=isb
              C(i)=0;
              for k=1:n
                  C(i)=C(i)+V(k)*(G(i,k)*cos(O(i)-O(k))+BI(i,k)*sin(O(i)-O(k)));%功率误差方程
              end
              DP1(i)=P(i)-V(i)*C(i);
              DP(i)=DP1(i)./V(i);
              DET=abs(DP1(i));
              if DET>=pr
                  ICT2=ICT2+1;
              end
          end
      end
      Np(k)=ICT2;
      if ICT2~=0        %
          for i=2:n
              DP(i)=B(i,i)*DP(i);
              if i~=n
                  IC1=i+1;
                  for k=IC1:n
                      DP(k)=DP(k)-B(k,i)*DP(i);
                  end
              else
                  for LZ=3:i
                      L=i+3-LZ;
                      IC4=L-1;
                      for MZ=2:IC4
                          I=IC4+2-MZ;
                          DP(I)=DP(I)-B(I,L)*DP(L);
                      end
                  end
              end
          end
          for i=2:n
              O(i)=O(i)-DP(i);
          end
          kq=1;L=0;
          for i=1:n
              if B2(i,6)==2
                  C(i)=0;L=L+1;
                  for k=1:n
                      C(i)=C(i)+V(k)*(G(i,k)*sin(O(i)-O(k))-BI(i,k)*cos(O(i)-O(k)));
                  end
                  DQ1(i)=Q(i)-V(i)*C(i);
                  DQ(L)=DQ1(i)./V(i);
                  DET=abs(DQ1(i));
              if DET>=pr
                  ICT3=ICT3+1;
              end
          end
      end
else
      kp=0;
      if kq~=0;
          L=0;
          for i=1:n
              if B2(i,6)==2
                  C(i)=0;L=L+1;
                  for k=1:n
                       C(i)=C(i)+V(k)*(G(i,k)*sin(O(i)-O(k))-BI(i,k)*cos(O(i)-O(k)));
                   end
                   DQ1(i)=Q(i)-V(i)*C(i);
                   DQ(L)=DQ1(i)./V(i);
                   DET=abs(DQ1(i));
               end
           end
       end
end
Nq(K)=ICT3;
if ICT3~=0
       L=0;
       for i=1:na
           DQ(i)=A(i,i)*DQ(i);
           if i==na
               for LZ=2:i
                   L=i+2-LZ;
                   IC4=L-1;
                   for MZ=1:IC4
                       I=IC4+1-MZ;
                       DQ(I)=DQ(I)-A(I,L)*DQ(L);
                   end
               end
           else
               IC1=i+1;
               for k=IC1:na
                   DQ(k)=DQ(k)-A(k,i)*DQ(i);
               end
           end
       end
       L=0;
       for i=1:n
           if B2(i,6)==2
               L=L+1;
               V(i)=V(i)-DQ(L);
           end
       end
       kp=1;
       K=K+1;
else
       kq=0;
       if kp~=0
           K=K+1;
       end
end
end
disp('迭代次数');
disp(K);
disp('每次没有达到精度要求的有功功率个数为');
disp(Np);
disp('每次没有达到精度要求的无功功率个数为');
disp(Nq);
for k=1:n
      E(k)=V(k)*cos(O(k))+V(k)*sin(O(k))*j;
      O(k)=O(k)*180./pi;
end
disp('各节点的电压标幺值E为(节点号从小到达排列)');
disp(E);
disp('各节点的电压大小V为(节点号从小到达排列)');
disp(V);
disp('各节点的电压相角O为(节点号从小到达排列)');
disp(O);
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:nl
      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:nl
      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:nl
      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