CN106229988B - 一种基于Matlab的极坐标牛顿法潮流计算方法 - Google Patents

Abstract

本发明公开了一种基于Matlab的极坐标牛顿法潮流计算方法，采用矩阵运算和复数运算，并按照Matlab的编程特点形成雅可比矩阵及计算节点功率，包括以下步骤：计算雅可比初始计算矩阵；计算节点复功率；计算初始的雅可比矩阵分块子矩阵；用节点复功率对雅可比矩阵分块子矩阵对角元进行修正；由雅可比矩阵分块子矩阵形成雅可比矩阵；对雅可比矩阵进行调整。本发明在Matlab平台实现，便于科研人员使用Matlab提供的各种工具和函数对计算结果进行测试和分析。实践证明，本发明既方便了科研人员对程序进行编写、修改和调试，同时计算速度也基本接近了在C语言平台上实现的速度，为科研人员的科研工作提供了一个优秀的分析工具。

Description

一种基于Matlab的极坐标牛顿法潮流计算方法

技术领域

本发明涉及一种电力系统牛顿法潮流计算方法，特别是一种适合研究目的使用的极坐标牛顿法潮流计算方法。

背景技术

电力系统潮流计算是研究电力系统稳态运行的一项基本计算，它根据给定的运行条件和网络结构确定整个网络的运行状态。潮流计算也是电力系统其他分析的基础，如安全分析、暂态稳定分析等都要用到潮流计算。极坐标牛顿法潮流计算方法是一种最常用的潮流计算方法，科研人员经常以极坐标牛顿法潮流计算为基础进行进一步地研究。实用的商业软件采用稀疏矩阵技术和节点优化编号等高级技术。这些技术虽然能大幅度提高潮流计算的速度、降低内存占用量，但编程非常麻烦且难以修改和维护，不易增加新的功能，因而不适合科研人员用于研究目的使用。

Matlab软件以矩阵为最基本的数据单位，可以方便地处理各种矩阵和向量运算，也可以很方便自然地处理复数类型，其指令表达式与数学中常用的形式很接近，还有大量常见实用的函数，给编程带来很大便利。Matlab软件简单易用、代码短小易操作，易于编程和调试，计算功能强大，同时还具有非常强大的可视化图形处理和交互式功能，为科学研究以及工程应用提供了一种高效的编程工具，目前已经成为许多科学领域的基本工具和首选平台，在各种科学和工程计算领域得到了广泛的应用。为了适应越来越多的科研人员需要在Matlab平台上以极坐标牛顿法潮流计算为基础进行进一步地研究的需求，迫切需要一种基于Matlab软件的易于编程、修改和调试的极坐标牛顿法潮流计算方法。

根据电力系统节点的特点，潮流计算把电力系统节点分成3类：节点有功功率和无功功率已知、节点电压幅值和电压相角未知的节点称为PQ节点；节点有功功率和电压幅值已知、节点无功功率和电压相角未知的节点称为PV节点；节点电压幅值和电压相角已知，节点有功功率和无功功率未知的节点称为平衡节点。

牛顿法潮流计算分为两类：牛顿法潮流计算中节点电压采用极坐标表示时的计算方法，称为极坐标牛顿法潮流计算方法；牛顿法潮流计算中节点电压采用直角坐标表示时的计算方法，称为直角坐标牛顿法潮流计算方法。极坐标牛顿法潮流计算主要方程如下：

节点导纳矩阵为：

式中，Y_ik为节点导纳矩阵元素，当下标i≠k时，为节点i和节点k的互导纳，当下标i＝k时，为节点i的自导纳；n为节点数。

节点功率方程为：

式中，P_i、Q_i分别为节点i的节点有功功率和无功功率；U_i、U_k分别为节点i和节点k的节点电压幅值；θ_ik＝θ_i-θ_k，θ_i和θ_k分别为节点i和节点k的节点电压相角；G_ik、B_ik分别为节点导纳矩阵元素Y_ik的实部和虚部。

功率偏差方程为：

式中，ΔP_i、ΔQ_i分别为节点i的节点有功功率偏差和无功功率偏差；P_iS、Q_iS分别为节点i给定的节点注入有功功率和注入无功功率；m为PQ节点数。

修正方程组为：

式中，J为雅可比矩阵，H、N、M、L为雅可比矩阵的分块子矩阵。

雅可比矩阵各元素计算公式为：

当j≠i时

当j＝i时

或用下列公式计算：

式中，P_i、Q_i分别为节点i的有功功率和无功功率，按式(2)计算。

如图1-2所示，现有极坐标牛顿法潮流计算方法，主要包括以下步骤：

A、原始数据输入和电压初始化；

原始数据包括线路和变压器支路数据、节点注入有功功率和无功功率、节点电压幅值、节点无功补偿数据，以及收敛精度、最大迭代次数。

电压初始化采用平启动，即PV节点和平衡节点的节点电压幅值取给定值，PQ节点的节点电压幅值取1.0；所有节点电压的相角都取0.0。这里单位采用标幺值。

B、形成节点导纳矩阵；

根据输入的线路和变压器支路数据形成如式(1)所示的节点导纳矩阵。

C、形成雅可比矩阵；

按式(5)～式(16)计算雅可比矩阵的各元素。

D、计算节点功率及功率偏差；

按式(2)计算节点功率，按式(3)计算节点功率偏差。

E、解方程及修正节点电压幅值U和相角θ；

解修正方程组(4)，求出电压幅值修正量列向量ΔU及电压相角修正量列向量Δθ。

电压修正公式为：

式中，上标(t)表示第t次迭代的值；ΔU_i和Δθ_i分别为节点i的节点电压幅值修正量和节点电压相角修正量。

F、判断功率最大不平衡量|ΔP|_max和|ΔQ|_max是否都小于收敛精度ε；如果都小于收敛精度ε，进行步骤G，否则返回步骤C进行下一次迭代；

G、计算平衡节点的有功功率和无功功率及PV节点的无功功率，计算各支路有功功率和无功功率，结束。

直接采用上述原理实现的极坐标牛顿法潮流计算软件计算速度较慢，商业使用的极坐标牛顿法潮流计算软件采用稀疏矩阵技术和节点优化编号技术，比较复杂，不适合科研人员以此为基础进一步进行科学研究。因此，中国专利ZL201010509556.5提出了一种适合研究目的使用的牛顿法潮流计算方法，为以极坐标牛顿法潮流计算为基础进行进一步研究的科研人员提供一个易于修改和维护的牛顿法潮流计算方法，其特点如下：

(1)不采用稀疏矩阵技术和节点优化编号，大大降低了算法的实现难度；

(2)通过简单逻辑判断来避免不必要运算，提高了潮流计算的计算速度。

中国专利ZL201010509556.5所提出方法为以极坐标牛顿法潮流计算为基础进行进一步研究的科研人员提供了一个易于修改和维护的极坐标牛顿法潮流计算方法。该方法采用C语言等编译型编程语言实现时速度很快，但使用Matlab这类解释型编程语言实现时计算速度则很慢，同时该专利也没有充分利用Matlab擅长矩阵运算和复数运算的特点。因此需要一个充分利用Matlab的特点且计算快速的潮流计算方法供在Matlab平台上进行科学研究的科研人员使用。

发明内容

为解决现有技术存在的上述问题，本发明要提出一种基于Matlab的极坐标牛顿法潮流计算方法，可以充分利用Matlab特有的擅长矩阵运算和复数运算的特点，同时又有较快计算速度的潮流计算方法。

为了实现上述目的，本发明的技术方案如下：一种基于Matlab的极坐标牛顿法潮流计算方法，采用矩阵运算和复数运算。包括以下步骤：

A、原始数据输入和电压初始化；

B、形成节点导纳矩阵；

C、形成雅可比矩阵及计算节点功率；

Matlab擅长矩阵运算和复数运算，因此采用Matlab编程，应该推导出基于矩阵运算和复数运算的雅可比矩阵计算方法。

雅可比矩阵元素与节点类型有关，常规形成雅可比矩阵时要判断节点类型，根据节点类型决定哪些节点需要形成雅可比矩阵元素。这样处理对于通过循环来实现的算法，容易处理，但不适合通过矩阵整体运算形成雅可比矩阵的方法。因此，本发明形成雅可比矩阵时，不判断节点类型，对所有节点都形成雅可比矩阵元素，之后再去掉不需要的行和列。

下面推导采用矩阵运算计算雅可比矩阵元素和节点功率的公式。

对式(5)～式(8)的分析，可得

M_ij＝-N_ij (19)

L_ij＝H_ij (20)

因此，先求H_ij和N_ij，求出H_ij和N_ij后，自然就能得到M_ij和L_ij。

推导由矩阵运算计算雅可比矩阵的公式之前，先看看雅可比矩阵各元素如何用复数或相量表示。

从式(5)和式(6)可见，两式中都包含U_iU_j和θ_ij项，说明雅可比矩阵各元素中应该包含乘以的共轭的项，即包含项；观察式(5)和式(6)括号内的表达式，可知雅可比矩阵各元素中还应该包含Y_ij的共轭项。因此H_ij、N_ij、M_ij、L_ij可以由下式生成：

式中，上标(^)表示复数的共轭；为节点电压列向量。

式(21)可以看成由下式与节点导纳矩阵对应位置元素相乘得到的矩阵J⁰的第i行第j列的元素：

式(22)则可由节点电压列向量与其共轭转置相乘得到，因此雅可比矩阵分块子矩阵可由下式得到：

式中，J⁰为雅可比初始计算矩阵；上标H表示矩阵的共轭转置；.*表示两矩阵对应行列的元素相乘。

由J⁰得到初始雅可比矩阵分块子矩阵为：

H⁰＝-Im(J⁰) (24)

N⁰＝-Re(J⁰) (25)

M⁰＝Re(J⁰) (26)

L⁰＝-Im(J⁰) (27)

式中，H⁰、N⁰、M⁰、L⁰为初始雅可比矩阵的分块子矩阵；Re表示取矩阵元素的实部；Im表示取矩阵元素的虚部。

由式(24)～式(27)得到的初始雅可比矩阵分块子矩阵的非对角元已经是雅可比矩阵元素了，对角元还需要修正。

由式(13)～式(16)，并考虑到sinθ_ii＝0和cosθ_ii＝1，得：

H_ii＝-U_iU_i(G_iisinθ_ii-B_iicosθ_ii)+Q_i (28)

N_ii＝-U_iU_i(G_ii cosθ_ii+B_ii sinθ_ii)-P_i (29)

M_ii＝U_iU_i(G_ii cosθ_ii+B_ii sinθ_ii)-P_i (30)

L_ii＝-U_iU_i(G_iisinθ_ii-B_iicosθ_ii)-Q_i (31)

式(28)～式(31)中等式右侧的第1项就是H⁰、N⁰、M⁰、L⁰的对角元，因此只需对得到的H⁰、N⁰、M⁰、L⁰用P_i和Q_i修正，就能得到雅可比矩阵分块子矩阵对角元。

观察式(2)和式(6)，以及N⁰，可得

式中，为矩阵N⁰第i行第k列的元素。

观察式(2)和式(5)，以及H⁰，可得

由式(24)、式(25)、式(32)、式(33)得到节点复功率为：

式中，为节点i的复功率。

由各节点复功率组成的节点复功率列向量可以用Matlab的一个矩阵求和函数实现：

式中，为节点复功率列向量；sum为Matlab的矩阵求和函数；2表示对矩阵每一行的元素求和。

用节点复功率对雅可比矩阵分块子矩阵对角元进行修正如下：

形成雅可比矩阵及计算节点功率，包括以下步骤：

C1、计算雅可比初始计算矩阵J⁰；

C2、计算节点复功率；

C3、由J⁰计算初始的雅可比矩阵分块子矩阵；

C4、用节点复功率对初始雅可比矩阵分块子矩阵对角元进行修正；

C5、由修正后的雅可比矩阵分块子矩阵形成雅可比矩阵；

C6、对雅可比矩阵进行调整，去掉PV节点无功功率偏差对应的行以及平衡节点有功功率偏差和无功功率偏差对应的行；去掉PV节点电压幅值修正量对应的列以及平衡节点电压幅值修正量和电压相角修正量对应的列，结束。

D、计算节点功率偏差；

式(3)计算节点功率偏差公式写成矩阵运算的形成为：

式中，ΔP、ΔQ分别为节点有功功率偏差列向量和无功功率偏差列向量；P_S、Q_S分别为给定的节点注入有功功率列向量和注入无功功率列向量。

计算得到的节点功率偏差列向量ΔP和ΔQ中去掉PV节点无功功率偏差及平衡节点有功功率偏差和无功功率偏差。

E、解方程及修正电压幅值U和相角θ；

直接调用Matlab软件的解线性方程组算法解修正方程组(4)，求出电压幅值修正量列向量ΔU及电压相角修正量列向量Δθ。

对电压进行修正的式(17)和式(18)改写成矩阵形式为：

U^(t+1)＝U^(t)-ΔU^(t) (42)

θ^(t+1)＝θ^(t)-Δθ^(t) (43)

式中，上标(t)表示第t次迭代的值；ΔU和Δθ分别为节点电压幅值修正量列向量和节点电压相角修正量列向量。

F、判断功率最大不平衡量|ΔP|_max和|ΔQ|_max是否都小于收敛精度ε；如果都小于收敛精度ε，进行步骤G，否则返回步骤C进行下一次迭代。

G、计算平衡节点的有功功率和无功功率及PV节点的无功功率，计算各支路有功功率和无功功率，结束。

与现有技术相比，本发明具有以下有益效果：

1、本发明提出的方法在Matlab平台实现，便于科研人员使用Matlab提供的各种工具和函数对计算结果进行测试和分析。

2、本发明提出的方法采用矩阵运算和复数运算，减少了程序代码，简化了编程，使得程序更加清晰，便于科研人员修改程序、对程序进行调试和改进、添加新功能。

3、由于Matlab对矩阵运算实行优化，采用矩阵运算要比按矩阵元素循环运算编程快得多，同时直接调用Matlab的方程求解算法，也大大提高了计算速度。实践证明，本发明的方法既方便了科研人员对程序进行编写、修改和调试，同时计算速度也基本接近了在C语言平台上实现的速度，为科研人员的科研工作提供了一个优秀的分析工具。

附图说明

本发明共有附图4张。其中：

图1是现有极坐标牛顿法潮流计算的流程图。

图2是现有极坐标牛顿法形成雅可比矩阵的流程图。

图3是本发明极坐标牛顿法潮流计算的流程图。

图4是本发明形成雅可比矩阵和计算节点功率的流程图。

具体实施方式

下面结合附图对本发明进行进一步地说明，按照图3-4所示流程对一个修改后的445节点实际系统算例进行了计算。

445节点实际大型电力系统有445个节点，544条支路，含有大量的小阻抗支路。为了对各种方法进行比较，把这些小阻抗支路改为正常阻抗支路以满足各种方法的要求。

采用本发明和几种对比方法对445节点实际系统算例进行了计算，计算时收敛精度为0.00001。几种潮流计算算法分别为：

方法1：中国专利ZL201010509556.5方法，采用Matlab语言实现。

方法2：中国专利ZL201010509556.5方法，采用Matlab语言实现，但解方程直接调用Matlab的解方程算法。

方法3：中国专利ZL201010509556.5方法，采用Matlab语言实现，但解方程直接调用Matlab的解方程算法，同时对矩阵变量预定义维数。

方法4：本发明方法。

几种方法的计算时间见表1，计算时间不包括数据读入和输出及支路功率计算的时间。

表1几种极坐标牛顿法潮流计算计算时间比较

潮流计算方法	计算时间(s)
		方法1	18.303
方法2	3.095
		方法3	0.927
方法4	0.519

从表1可见，直接采用Matlab实现中国专利ZL201010509556.5方法，计算时间很长；采用Matlab实现专利ZL201010509556.5方法时，如果直接调用Matlab的解方程方法，计算速度就能大大提高，Matlab的解方程方法比较成熟稳定，有利于算法稳定性，Matlab的解方程方法的调用也非常简单，简化了编程，使程序更加清晰；如果同时对矩阵变量预定义维数，避免程序在执行过程中不断扩充矩阵的大小，可以大大提高计算速度。本发明的计算结果表明形成雅可比矩阵时采用矩阵运算技术可以进一步较大幅度地提高计算速度，程序也进一步简化。

本发明可以在任何版本的MATLAB编程语言实现，但建议使用较新版本的MATLAB语言。

本发明不局限于本实施例，任何在本发明披露的技术范围内的等同构思或者改变，均列为本发明的保护范围。

Claims (1)

1.一种基于Matlab的极坐标牛顿法潮流计算方法，包括以下步骤：

A、原始数据输入和电压初始化；

原始数据包括线路和变压器支路数据、节点注入有功功率和无功功率、节点电压幅值、节点无功补偿数据，以及收敛精度和最大迭代次数；

电压初始化采用平启动，即PV节点和平衡节点的节点电压幅值取给定值，PQ节点的节点电压幅值取1.0；所有节点电压的相角都取0.0；这里单位采用标幺值；

所述的PQ节点为有功功率和无功功率已知、电压幅值和电压相角未知的节点，所述的PV节点为有功功率和电压幅值已知、无功功率和电压相角未知的节点，所述的平衡节点为电压幅值和电压相角已知，节点的有功功率和无功功率未知的节点；

B、形成节点导纳矩阵；

根据输入的线路和变压器支路数据形成如式(1)所示的节点导纳矩阵；

式中，Y_ik为节点导纳矩阵元素，当下标i≠k时，为节点i和节点k的互导纳，当下标i＝k时，为节点i的自导纳；n为节点数；

其特征在于：

C、形成雅可比矩阵及计算节点功率；

C1、计算雅可比初始计算矩阵J⁰；

式中，J⁰为雅可比初始计算矩阵；上标^表示复数的共轭；为节点电压列向量；上标H表示矩阵的共轭转置；.*表示两矩阵对应行列的元素相乘；

C2、计算节点复功率；

由各节点复功率组成的节点复功率列向量用Matlab的一个矩阵求和函数实现：

式中，为节点复功率列向量；sum为Matlab的矩阵求和函数；2表示对矩阵每一行的元素求和；

C3、由J⁰计算初始的雅可比矩阵分块子矩阵；

由J⁰得到初始雅可比矩阵分块子矩阵为：

H⁰＝-Im(J⁰) (4)

N⁰＝-Re(J⁰) (5)

M⁰＝Re(J⁰) (6)

L⁰＝-Im(J⁰) (7)

式中，H⁰、N⁰、M⁰、L⁰为初始雅可比矩阵的分块子矩阵；Re表示取矩阵元素的实部；Im表示取矩阵元素的虚部；

C4、用节点复功率对初始雅可比矩阵分块子矩阵对角元进行修正；

用节点复功率对初始雅可比矩阵分块子矩阵对角元进行修正如下：

C5、由修正后的雅可比矩阵分块子矩阵形成雅可比矩阵；

式中，J为雅可比矩阵，H、N、M、L为修正后的雅可比矩阵的分块子矩阵；

C6、对雅可比矩阵进行调整，去掉PV节点无功功率偏差对应的行以及平衡节点有功功率偏差和无功功率偏差对应的行；去掉PV节点电压幅值修正量对应的列以及平衡节点电压幅值修正量和电压相角修正量对应的列，结束；

D、计算节点功率偏差；

按下式计算节点功率偏差：

式中，ΔP、ΔQ分别为节点有功功率偏差列向量和无功功率偏差列向量；P_S、Q_S分别为给定的节点注入有功功率列向量和注入无功功率列向量；

计算得到的节点功率偏差向量ΔP和ΔQ中去掉PV节点无功功率偏差及平衡节点有功功率偏差和无功功率偏差；

E、解方程及修正电压幅值U和相角θ；

由步骤C得到雅可比矩阵J和步骤D得到节点功率偏差向量ΔP和ΔQ构造成潮流计算的修正方程如下：

直接调用Matlab软件的解线性方程组算法解修正方程组(14)，求出电压幅值修正量向量ΔU及电压相角修正量向量Δθ；

按下式对节点电压幅值和相角进行修正：

U^(t+1)＝U^(t)-ΔU^(t) (15)

θ^(t+1)＝θ^(t)-Δθ^(t) (16)

式中，上标(t)表示第t次迭代的值；

F、判断功率最大不平衡量|ΔP|_max和|ΔQ|_max是否都小于收敛精度ε；如果都小于收敛精度ε，进行步骤G，否则返回步骤C进行下一次迭代；

G、计算平衡节点的有功功率和无功功率及PV节点的无功功率，计算各支路有功功率和无功功率，结束。