IEEE33节点潮流计算matlab程序——改进牛顿法潮流计算

改进牛顿法的基本原理

参考文献：一种新的配电网潮流算法——改进牛顿法-拉夫逊法

牛顿法是改进牛顿法的基础，对牛顿法作科学的近似，即雅可比矩阵做一些更改，使得每次计算得出的修正量都有所改变，但由于收敛精度恒定，最终计算出的结果误差较小，计算其结果可取，有意义。对系统进行条件假设：

1. 不存在对地支路(并联电容器组)。
2. 支路两端电压差值小。

假设（1）是当系统并联电容器组、恒定阻抗负载或较大导纳值时，这是不合理的，但当以上并联支路通过电压转换为节点功率或注入电流后，假设（1）将可以被接受。由于配电线路较短，潮流不大，假设（2）通常将会被默认接受。

由于假设（1）与（2）成立，又，则公式(22)中雅可比矩阵元素可简约为：

                      （31）

因为每个节点一般只与3-5个节点形成支路，固N、H、L、J 与导纳矩阵 Y 都有稀疏性，同时又具有对称性。可把（31）简化为：

                          （32）

                                （33）

其中， DB和DG为对角阵，仅由支路参数决定，对角元素分别为和；An-1为除源节点的节点—支路关联矩阵，为上三角阵，主对角元素取 1，非零非对角元素取-1，在程序中无须真正形成，可通过系统拓扑结构获得。

因此，公式（22）作进一步修改为：

           （34）

对节点与支路进行编号，可得An-1。本文按离根节点的距离，对支路进行分层编号，从而可以形成An-1阵。在已编号的网络中，支路方向指向源节点，节点与支路都要编号，支路入端节点编号为支路编号。

若定义：

                       （35）                                

                         （36）                                         

                          （37）

则公式(34)可写成

                                                （38）                                               （39）

                         （40）

公式(39)对应回代计算，而公式(40)对应计算的前推过程。在前推计算时可以求出E，W阵的逆为对角阵，对角元可用线路阻抗表示

                            （41）

                     （42） 

                                           （43）

Rij、Xij分别为支路 i—j 之间的电阻和电抗。

当系统接有分散的发电机时，若发电机的接点表征PQ节点，则算法不变；当表征PV节点，则公式（34）中未知，而已知为零，得：

                           （44）

其中，Xb部分元素已知，bx部分元素未知。且Xb已知元素与bx未知元素属方程组同一行。进而公式（44）可改为：

                          （45）

                          （46）

公式（45）中b1包括所有已知元素，将未知元素置为零。公式（46）中b2只有未知元素，其余元素则为零。这样就可解出X1与X2中的元素。

当系统中由环路时，即对重要的用户采用两端供电，这样就形成了一个环，这个环只可能在负荷点。假设i点为为环中被选的解裂点，分成m点和n点，则存在如下的边界条件：

，       

，        

对U进行一次行变换，对UT进行一次列变换，UDUT*X=b可以化成公式（44）。这样，对点m和n的处理与上述处理PV节点的方法相同。

病态网络不收敛的原因，一般是初值选取不当，也可能是雅可比矩阵自身缺陷所引起的。近似处理时，雅可比矩阵为UDUT形式，该阵被用来决定搜索的方向，它的线性潮流方程被用作前推回代的基础，以计算状态变量的修正增量。其中，D为对角阵，有助于避免显式形成，进而避免病态。U为仅由系统拓扑决定的上三角阵。改进牛顿法具有前推回代法的收敛性，但与前推回代法还是有很大区别，后者根本不需要计算潮流方程的偏微分，是根据欧姆定律、KVL和KCL，对网络进行前推回代，可求出状态变量的修正增量。

综上所述可知，改进牛顿法的优点在于，UDUT形式的矩阵不需要显式形成，而是直接进行前推回代，可避免雅可比矩阵和LU分解因子相关的可能的病态。另外，它是牛顿法，可以用于状态估计。在牛顿法中，潮流方程的偏微分就是雅可比矩阵，用以决定搜索方向， 再用LU分解的因子进行前推回代，以计算状态变量的修正增量。

以IEEE33节点为例进行算例计算分析：

IEEE33节点系统结构如下：

 IEEE33节点系统结构参收如下：

 潮流计算程序流程如下：

①选取恰当的电压及功率基准值，并得出电压、功率与阻抗标幺值。

②给支路与节点标号。本文使用的标号方法为对支路进行分层标号。根节点标号为零，支路入端节点作为支路编号。

③根据系统拓扑求取矩阵An-1，其主对角线元素为1。-1元素的确定:第m条支路入端节点为n时，其元素An-1（m，n）= -1。

④求得导纳矩阵

⑤赋予节点电压与相位初值。

⑥计算功率偏差并判断是否收敛。是，则输出，否，则进行下一步计算。

⑦求出电压与相位的修正量。

⑧对电压与相位进行修正并从步骤⑥重新开始计算。

程序流程图：

 部分程序展示如下： 

clc;
close all
Ub=12.66;%电压基准值kv
Sb=10;%电压基准值MVA
Zb=Ub*Ub/Sb;%阻抗基准值
......
%% 求导纳矩阵
Y=zeros(33,33);
for m=1:33 
for n=1:33
Y(m,m)=sum(1./Z0(m,:),2);
if n==m
Y(m,n)=Y(m,n);
else 
Y(m,n)=-1/Z0(m,n);
end
end
end
G=real(Y);
B=imag(Y);
........
%% 初始值
U=ones(33,1);%电压初值
thelta=zeros(33,1);%相位初值
Req=zeros(32,1);
Xeq=zeros(32,1);
P=-importdata('p.dbf')/1000/Sb;%读取dbf文件表格参数(pq有功量),并标幺化
deltaP=zeros(32,1);
Q=-importdata('q.dbf')/1000/Sb;%读取dbf文件表格参数(pq无功量)
PQ=zeros(32,1);
deltaQ=zeros(32,1);
Circulation=0;%循环次数
.....
%% 计算迭代
while Precision>10^-5%收敛精度
    Circulation=Circulation+1;
      
   for m=2:33 
        Sp=0;%计算功率偏差
       Sq=0;
       for n=1:33           
             Sp=Sp+U(m)*U(n)*(G(m,n)*cos(thelta(m)-thelta(n))+B(m,n)*sin(thelta(m)-thelta(n)));%
             Sq=Sq+U(m)*U(n)*( G(m,n)*sin(thelta(m)-thelta(n))-B(m,n)*cos(thelta(m)-thelta(n)));%  
       end
        PQ(m-1)=Sp+Sq*1i;        
       deltaP(m-1)=P(m-1)-Sp;
       deltaQ(m-1)=Q(m-1)-Sq;      
   end   
  Precision=max(abs(deltaQ));%收敛条件
 for m=1:32
    S(m)=deltaP(m)+deltaQ(m)*1i; %定义S
 end 
   SL=A\S;      %SL   
   %求矩阵W的逆
   Req(1)=X(1,2)/(U(1)*U(2)*cos(thelta(1)-thelta(2)));
   Xeq(1)=R(1,2)/(U(1)*U(2)*cos(thelta(1)-thelta(2)));
   ......
%% 输出结果  
U
thelta
Circulation
%% 画图
figure
plot(U)
xlabel('节点序号')
ylabel('节点电压标幺值')
title(' 改进牛顿法潮流计算电压分布')
figure
plot(thelta)
xlabel('节点序号')
ylabel('节点相位')
title(' 改进牛顿法潮流计算相位分布')

• 计算结果如下：
• 潮流计算电压分布：
• 

潮流计算各节点相位分布：