求解电力系统潮流计算结果不收敛问题

Question

以下的代码是计算电力系统潮流分析的MATLAB实现，针对给出的数据计算结果却拉满了迭代次数，求各位帮忙找找原因
之前尝试使用其他方式的电力系统潮流分析的方法进行计算，使用的是同样的数据，也就是说这组数据是可以计算的，但是本代码最后却不收敛，也就是说代码应该是存在问题的，用别的方法算出来差不多是迭代6次

% (bus#)(volt)(ang)(p)(q)(bus type)
clear;
bus=[
   1   1.00  0.00   -1.60  -0.80  1;
   2   1.00  0.00   -2.00  -1.00  1;
   3   1.00  0.00   -3.70  -1.30  1;
   4   1.05  0.00    5.00    0.00  2;
   5   1.05  0.00    0.00    0.00  3];

% b#1 b#2  ( R )( X )( G )( B ) ( K )
line = [
   1   2 0.04 0.25  0.0 0.25    0;
   1   3 0.10 0.35  0.0 0.0     0;
   2   3 0.08 0.30  0.0 0.25    0;
   5   3 0.00 0.03  0.0 0.0  1.05;
   4   2 0.00 0.015 0.0 0.0  1.05];

function [bus,line,nPQ,nPV,nodenum]=Sort(bus,line)
[nb,mb]=size(bus);
[nl,ml]=size(line);
  global nSW;
  global nPV;
  global nPQ;

  nSW = 0;                 % number of swing bus counter
  nPV = 0;                 % number of PV bus counter
  nPQ = 0;                 % number of PQ bus counter
  for i = 1:nb             % nb为总节点数
    type= bus(i,6);
    if type == 3
       nSW = nSW + 1;     % increment swing bus counter
       SW(nSW,:)=bus(i,:);
    elseif type == 2
        nPV = nPV +1;      % increment PV bus counter
       PV(nPV,:)=bus(i,:);
    else
        nPQ = nPQ + 1;     % increment PQ bus counter
       PQ(nPQ,:)=bus(i,:);
    end
  end
  
  bus=[PQ;PV;SW];
  newbus=[1:nb]';
  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 

function Y= y(bus,line)
% Purpose:   build admittance matrix Y from the line data
%
% Input:     bus  - bus data
%            line - line data
%
% Output:    Y    - admittance matrix

[nb,mb]=size(bus);
[nl,ml]=size(line);

Y=zeros(nb,nb);         % 对导纳矩阵赋初值0 
for k=1:nl
    I=line(k,1);                       %读入线路参数
    J=line(k,2);
    Zt=line(k,3)+j*line(k,4);
    Yt=1/Zt;
    Ym=line(k,5)+j*line(k,6);
    K=line(k,7);
    
    if (K==0)&(J~=0)                 % 普通线路: K=0;
        Y(I,I)=Y(I,I)+Yt+Ym;
        Y(J,J)=Y(J,J)+Yt+Ym;
        Y(I,J)=Y(I,J)-Yt;
        Y(J,I)=Y(I,J);
    end
    if (K==0)&(J==0)               % 对地支路: K=0,J=0,R=X=0;
        Y(I,I)=Y(I,I)+Ym;
    end
    if K>0                        % 变压器线路: Zt和Ym为折算到i侧的值,K在j侧
        Y(I,I)=Y(I,I)+Yt+Ym;
        Y(J,J)=Y(J,J)+Yt/K/K;
        Y(I,J)=Y(I,J)-Yt/K;
        Y(J,I)=Y(I,J);
    end
    if K<0                      % 变压器线路: Zt和Ym为折算到K侧的值,K在i侧
        Y(I,I)=Y(I,I)+Yt+Ym;
        Y(J,J)=Y(J,J)+K*K*Yt;
        Y(I,J)=Y(I,J)+K*Yt;
        Y(J,I)=Y(I,J);
    end
end

function [delP,delQ] =dPQ(Y,bus,nPQ,nPV)
[n,m]=size(Y);
[nb,mb]=size(bus);
G=zeros(n,m);
B=zeros(n,m);
for i=1:n
    for j=1:m
        G(i,j)=real(Y(i,j));
        B(i,j)=imag(Y(i,j));
    end
end
temp_delP=0;
temp_delQ=0;
for i=1:nb-1
    for j=1:m
     temp_delP=temp_delP+bus(i,2)*bus(j,2)*(G(i,j)*cos(bus(i,3)-bus(j,3))+B(i,j)*sin(bus(i,3)-bus(j,3)));
    end
      delP(i,1)= bus(i,4)-temp_delP;
end

for i=1:nPQ
    for j=1:m
      temp_delQ=temp_delQ+bus(i,2)*bus(j,2)*(G(i,j)*sin(bus(i,3)-bus(j,3))-B(i,j)*cos(bus(i,3)-bus(j,3)));
    end
      delQ(i,1)=bus(i,5)-temp_delQ;
end

function J=form_jac(bus,Y,nPQ,nPV)
[n,m]=size(Y);
[nb,mb]=size(bus);
G=zeros(n,m);
B=zeros(n,m);
H=zeros(nb-1,nb-1);
N=zeros(nb-1,nPQ);
K=zeros(nPQ,nb-1);
L=zeros(nPQ,nPQ);

for i=1:n
    for j=1:m
        G(i,j)=real(Y(i,j));
        B(i,j)=imag(Y(i,j));
    end
end

for i=1:nb-1
    for j=1:nb-1
        if i~=j
            H(i,j)=-bus(i,2)*bus(j,2)*(G(i,j)*sin(bus(i,3)-bus(j,3))-B(i,j)*cos(bus(i,3)-bus(j,3)));
        elseif i==j
          % H(i,j)=bus(i,2).^2*B(i,i)+bus(i,5);
            for jj=1:nb
                if jj~=i
                   H(i,j)= H(i,j)+bus(jj,2) *(G(i,jj)*sin(bus(i,3)-bus(jj,3))-B(i,jj)*cos(bus(i,3)-bus(jj,3)));
                end
            end
            H(i,j)=bus(i,2)*H(i,j);
        end 
    end
end

for i=1:nb-1
    for j=1:nPQ
        if i~=j
             N(i,j)=-bus(i,2)*bus(j,2)*(G(i,j)*cos(bus(i,3)-bus(j,3))+B(i,j)*sin(bus(i,3)-bus(j,3)));
        elseif i==j
            % N(i,i)=-bus(i,2).^2*G(i,i)-bus(i,4);  
            for jj=1:nb
                if jj~=i
                   N(i,i)= N(i,i)+bus(jj,2) *(G(i,jj)*cos(bus(i,3)-bus(jj,3))+B(i,jj)*sin(bus(i,3)-bus(jj,3)));
                end
            end
            N(i,i)=-bus(i,2).*N(i,i)-2*bus(i,2).^2.*G(i,i);
        end      
    end
end

for i=1:nPQ
    for j=1:nb-1
        if i~=j
             K(i,j)=bus(i,2)*bus(j,2)*(G(i,j)*cos(bus(i,3)-bus(j,3))+B(i,j)*sin(bus(i,3)-bus(j,3)));
        elseif i==j
            % K(i,i)=bus(i,2).^2*G(i,i)-bus(i,4);  
            K(i,i)=N(i,i)+2*bus(i,2).^2.*G(i,i);
        end      
    end
end

for i=1:nPQ
    for j=1:nPQ
        if i~=j
             L(i,j)=-bus(i,2)*bus(j,2)*(G(i,j)*sin(bus(i,3)-bus(j,3))-B(i,j)*cos(bus(i,3)-bus(j,3)));
        elseif i==j
             %L(i,i)=bus(i,2).^2*B(i,i)-bus(i,5);  
            L(i,i)= -H(i,i)+2*bus(i,2).^2.*B(i,i);
        end      
    end
end

%合成雅可比矩阵
J=[H N;K L];

data1;
Max=100;
%[bus,line]=Openfile;
%tic;
[nb,mb]=size(bus);
[nl,ml]=size(line);
[bus,line,nPQ,nPV,nodenum]=Sort(bus,line);
Y= y(bus,line);
myf = fopen('Output.m','w');
fclose(myf);


format long
eps = 1.0e-10;
for t=1:Max
   [delP,delQ] =dPQ(Y,bus,nPQ,nPV);
   J = form_jac(bus,Y,nPQ,nPV);
   UD = zeros(nPQ,nPQ);
   
   for i=1:nPQ
       UD(i,i) = bus(i,2);
   end
   
   A = J;
   b = [delP;delQ];
   %x=Gauss(A,b);
   x=A\b;
   dx = x(1:(nb-1),1);
   dU = UD*(x(nb:(nb+nPQ-1),1));
   
   bus(1:nPQ,2) = bus(1:nPQ,2)-dU;
   bus(1:nb-1,3) = bus(1:nb-1,3)-dx;
   sprintf('第%d次的迭代结果：\n',t)
  % disp(J);
   %disp(dx);
   %disp(dU);
   disp(bus);
   if (max(abs(dU))<eps)&&(max(abs(dx))<eps) 
       break
   end
end

Answer 1

引自强大而免费的微信小程序：皆我百晓生。

问题分析

潮流计算是电力系统的一种重要分析方法，用于计算电网中各节点的电压幅值和相角。而潮流计算的迭代过程最终应该收敛于某个稳定的状态，然而有时会出现不收敛的情况。导致不收敛的原因有很多，可能是数据问题、初始值问题、迭代步长设置问题等。

根据给出的MATLAB代码，推测可能出现不收敛问题的原因有：

1. 初始值问题：初始电压相角和幅值的设置可能不合适，导致迭代的起点离收敛点较远。
2. 数据配置问题：节点类型、线路参数的配置可能不正确或者有误，导致计算过程出现不收敛的情况。
3. 迭代步长设置问题：迭代步长可能过大或者过小，影响收敛速度和稳定性。

接下来，在代码分析的基础上，逐步详细解释其中的问题，给出相应的解决方案。

代码分析

首先，代码的前部分定义了bus和line的数据，并且包含了一些辅助函数。

接下来，代码中的Sort函数用于对bus和line数据进行排序和编号。不收敛的问题与排序和编号无关，所以这部分代码可以暂时不关注。

然后，y函数用于构建导纳矩阵Y，根据给出的线路参数计算Y的元素。这部分代码看起来没有问题。

再然后，dPQ函数用于计算节点注入有功和无功功率与节点电压幅值和相角之间的误差。根据给出的方法计算出delP和delQ。这部分代码也看起来没有问题。

接下来，form_jac函数用于构建雅可比矩阵J，J矩阵是潮流计算中的重要参数，它描述了电压相角和幅值的变化与节点注入功率之间的关系。这部分代码的逻辑看起来也没有问题。

然后，代码进入主循环部分，使用迭代的方法计算电压相角和幅值直到收敛。主要的迭代过程有以下几个步骤：

1. 计算delP和delQ，即节点注入功率的计算误差。
2. 构建雅可比矩阵J。
3. 求解线性方程组J*x = b，得到更新量dx和dU。
4. 根据更新量dx和dU更新节点电压相角和幅值。
5. 判断迭代是否收敛，如果没有收敛则继续迭代，否则结束循环。

解决方案

根据代码分析和问题分析，可以考虑以下几个方面进行解决：

1. 初始值设置： 可以尝试调整初始电压相角和幅值的设置，使其更接近收敛点。可以根据电力系统的实际情况、历史数据或者其他方法来设置合适的初始值。

2. 数据配置检查： 需要仔细检查bus和line数据的配置是否正确、完整。可以通过打印输出或者断点调试的方式检查数据是否符合预期。可以特别注意节点类型的设置是否正确，线路参数的设置是否准确。

3. 迭代步长设置： 可以尝试调整迭代步长的设置，如果步长过大可能导致迭代振荡或者发散，步长过小可能导致迭代收敛速度过慢。可以根据实际情况进行调整。

4. 迭代次数限制： 可以对迭代次数进行限制，达到最大迭代次数后强制结束迭代。可以在循环中加入判断语句，当迭代次数超过一定阈值时结束迭代。

5. 输出调试信息： 可以在迭代过程中添加打印语句，输出每次迭代结果、J矩阵、dx、dU等信息，方便查看和分析迭代过程中的变化和误差。

综上所述，可以根据上述解决方案逐一进行尝试和调整，从而解决电力系统潮流计算不收敛的问题。根据实际情况不同，可能需要进行多次尝试和调整才能找到最合适的解决方案。