ABSTRACT

In order to overcome the shortcoming of traditional Lyapunov direct method which is conservative in determining the stability region of power system, a method for calculating the stability region of power system based on special points is studied. The main idea is to study and analyze the three stable forms of energy function V, and find their respective stable regions, and take the union of the stable regions corresponding to various stable forms as the stable regions of the system. In the process of solving, the boundary of the stable region is approximated by the special point method. The simulation results of a single machine infinite bus system show that the method is effective and improves the accuracy of stability region calculation.

Keywords:Power System, Lyapunov’s Direct Method, Stability Region, Union, Special points

基于并集及特殊点求解的稳定域边界分析方法

黄 罡

华北电力大学电气与电子工程学院，北京

收稿日期：2019年8月5日；录用日期：2019年8月20日；发布日期：2019年8月27日

摘 要

针对传统的李雅普诺夫直接法确定电力系统稳定域保守性强的缺点，研究了一种基于特殊点计算电力系统稳定域的方法。主要思路是：通过研究分析能量函数V的3种稳定形式，并求取它们各自对应的稳定域；将各种稳定形式对应的稳定域的并集作为系统的稳定域。在求解过程中，利用特殊点法对稳定域边界进行逼近。以单机无穷大系统作为案例进行仿真验证，结果表明了该方法的有效性，提高了稳定域计算的精度。

关键词 :电力系统，李雅普诺夫直接法，稳定域，并集，特殊点

Copyright © 2019 by author and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

http://creativecommons.org/licenses/by/4.0/

1. 引言

稳定性问题是电力系统需要解决的最核心问题之一，即系统在受到各种扰动后能否保持稳定运行。电力系统在受到大扰动后，由于发电机组的机械功率输入与电磁功率输出的不平衡，引起发电机转子间产生相对摇摆，若系统不能恢复到某一平衡状态，则系统将失去稳定。电力系统的稳定域表征了该系统的抗干扰能力，是量化系统稳定性的重要指标。只要能够将电力系统的工作状态保持在稳定域范围内，就不会有失稳的风险，因此稳定域的确定对于保障电力系安全运行是至关重要的 [1] [2] [3] [4] [5] 。

稳定域的计算方法很多，一种方法是基于Lyapunov能量函数对稳定域近似估计 [6] [7] ，但得到的稳定域保守性较强，另外，实际的电力系统并不一定存在直接法对应的Lyapunov型能量函数。文献 [8] 基于蒙特卡洛方法对稳定域进行近似计算，初始点取得越多，得到的稳定域精度越高；但初始点取得过多，会导致计算机运算时间太长。文献 [9] [10] [11] 提出基于可达集理论确定电力系统的稳定域，该方法首先在电力系统的状态空间找到一个较小的稳定区域并将其作为目标集，然后逆时间求解该目标集得到可达集，将可达集作为电力系统的稳定域。但该方法依靠最优控制求取Hamilton-Jacobi方程的最优解，而寻找最优解的过程较繁琐且耗时较长。

近年来，随着现代微分动力系统稳定边界理论的发展，学者们提出了主导不稳定平衡点(CUEP)处稳定边界的超平面近似 [12] 和二次近似 [13] 。文献 [14] 基于Taylor级数提出了一个对稳定流形的二次近似，但该近似是非精确的，需要动能校正。文献 [15] 提出利用递推方法近似计算稳定域边界各项系数的方法，但未有实际算例。

李雅普诺夫稳定性理论只描述了系统稳定的一种形式，因此只给出了系统稳定的充分条件，本文针对Lyapunov直接法的保守性，通过描述能量函数V的其他2种典型的稳定形式，并求取各种稳定形式对应的稳定域，取它们的并集作为稳定域的估计，从而使所求稳定域更加逼近真实稳定域。在求取其他稳定形式对应的稳定域时，通过求取稳定边界上的特殊点，从而实现对稳定域边界的近似。仿真算例验证了上述理论的正确性，结果显示本文采用的方法可以有效地估计稳定域边界，且近似程度较高。

2. 基本原理

李雅普诺夫直接法通过构建李雅普诺夫函数，李雅普诺夫函数满足图1(a)中曲线1所示的稳定形式，进而求得系统对应第一种稳定形式的稳定域 $R_1$ 。由于李雅普诺夫定理是系统稳定性的充分条件，因此求得的稳定域也是偏保守的。系统的实际稳定域R应该是所有稳定形式相对应的稳定域的并集，数学表达形式为 $R=R_1\cup R_2\cup R_3\cup \cdots \cup R_n$ ，如图1(b)所示。

图1. 典型的稳定形式及稳定域构成示意图

动力系统趋于稳定有多种表达方式，图1(a)给出了3种典型的稳定形式。曲线1对应函数值变化率始终为负值，即运动过程中函数值单调下降，那么系统的受扰运动会达到其平衡状态，而且函数值在运动过程中都是非负的。李雅普诺夫稳定定理主要对曲线1所示的稳定形式进行分析求解。曲线2的稳定形式与曲线1是相反的。曲线2对应函数值变化率始终为正值，即运动过程中函数值单调上升，那么系统的受扰运动会达到其平衡状态，而且函数值在运动过程中都是非正的。曲线3对应振荡衰减的形式，该种形式下系统最终也将达到其平衡状态。

本文通过对图1(a)中的曲线2和曲线3所示的稳定形式进行分析求解，使所求的稳定域更加逼近真实的稳定域。对于第2种稳定形式对应的稳定域 $R_2$ ，经证明，与第一种稳定形式对应的稳定域是一样的，即 $R_1=R_2$ ，如图1(b)所示。

3. 稳定域边界分析

本文以单机无穷大母线系统为例进行分析求解。不计阻尼因素时，采用初积分法可构造出不计阻尼的能量型李雅普诺夫函数，其构造过程可参考参考文献 [16] ，进而求得系统不计阻尼时的稳定域。

其中稳定域边界的表达式如下所示：

$\frac{1}{2}M{\omega }^2-P_M\left(\delta -{\delta }_s\right)-P_{em3}\left(\cos \delta -\cos {\delta }_s\right)=V_{cr}$ (1)

式中， $V_{cr}=V\left({\delta }_u,0\right)=V\left(\text{π}-{\delta }_s,0\right)=P_{em3}\left[2\cos {\delta }_s-\left(\text{π}-2{\delta }_s\right)\sin {\delta }_s\right]$ 。

当计及阻尼因素时，在 $\left(\delta ,\omega \right)$ 平面上，系统的真实稳定域表示如图2所示。稳定域边界经过不稳定平衡点 ${\delta }_u$ 。

真实的稳定域边界求解困难，目前还没有相关文献求出稳定域边界具体的数学表达形式。对于计及阻尼的稳定域构成，如图1(b)所示，真实稳定域的 $R_1$ 和 $R_2$ 部分可以利用李雅普诺夫直接法求得，现在考虑求取 $R_3$ 部分，将 $R_1$ 、 $R_2$ 和 $R_3$ 的并集作为系统的稳定域，使所求稳定域更加逼近真实稳定域。直接求取 $R_3$ 比较困难，考虑通过求取稳定域边界上的特殊点，进而实现稳定域边界的近似。

图2. 有阻尼的单机无穷大系统的稳定域

3.1. 稳定域R₃的分析

对于稳定域 $R_3$ ，得到一个数学形式的稳定性定理。若在原点的小邻域R存在一个标量函数 $H\left(x\right)$ ，具有如下条件：如果沿着轨线 $x\left(t\right)$ ，函数 $H\left(x\right)$ 的变化趋势呈现如图1(a)所示的振荡衰减形式，这表示随着时间的增加，轨线 $x\left(t\right)$ 沿着收缩的曲面最终收敛到原点。那么，原点 $x=0$ 在域内就是渐进稳定的，小邻域R为第3种稳定形式对应的稳定域 $R_3$ ，如图3所示。

图3. 稳定域R₃示意图

对于单机无穷大母线系统，在 $\left(\delta ,\omega \right)$ 平面上，构造满足第3种稳定形式的函数 $H\left(\delta ,\omega \right)$ ，其数学表达式如下：

$H\left(\delta ,\omega \right)=H_1\left(\delta ,\omega \right)H_2\left(\delta ,\omega \right)=\left[\frac{1}{2}M{\omega }^2-P_M\left(\delta -{\delta }_s\right)-P_{em3}\left(\cos \delta -\cos {\delta }_s\right)\right]\cdot \left(\delta -{\delta }_s\right)$ (2)

以下对函数 $H\left(\delta ,\omega \right)$ 进行解释说明，函数 $H\left(\delta ,\omega \right)$ 由2部分组成 $H_1\left(\delta ,\omega \right)$ 和 $H_2\left(\delta ,\omega \right)$ ，其中 $H_1\left(\delta ,\omega \right)$ 对时间的变化率为

$\frac{\text{d}{H}_1}{\text{d}t}=\frac{\partial H_1}{\partial \delta }\frac{\text{d}\delta }{\text{d}t}+\frac{\partial H_1}{\partial \omega }\frac{\text{d}\omega }{\text{d}t}=\left(-P_M+P_{em3}\sin \delta \right)\frac{\text{d}\delta }{\text{d}t}+M\frac{\text{d}\delta }{\text{d}t}\frac{\text{d}\omega }{\text{d}t}=\left(-P_M+P_{em3}\sin \delta +M\frac{\text{d}\omega }{\text{d}t}\right)\frac{\text{d}\delta }{\text{d}t}=-D{\left(\frac{\text{d}\delta }{\text{d}t}\right)}^2<0$

而 $H_2\left(\delta ,\omega \right)=\delta -{\delta }_s$ 在故障后稳定平衡点 $\left({\delta }_s,0\right)$ 的小邻域R内沿着轨线是呈振荡形式变化的。

因此 $H\left(\delta ,\omega \right)$ 在故障后稳定平衡点 $\left({\delta }_s,0\right)$ 的小邻域R内沿着轨线是呈现振荡衰减的形式变化的。这说明故障后稳定平衡点 $\left({\delta }_s,0\right)$ 在小邻域R内是渐进稳定的，小邻域R为第3种稳定形式对应的稳定域 $R_3$ 。直接求取 $R_3$ 比较困难，考虑通过求取稳定域边界上的特殊点，进而实现稳定域边界的近似。

3.2. 特殊点的求解

下面介绍在 $\left(\delta ,\omega \right)$ 平面确定稳定域 $R_3$ 边界上特殊点的过程，在稳定域 $R_3$ 内， $H_2\left(\delta ,\omega \right)$ 沿着轨线会满足振荡特性，轨线经过横轴 $\omega =0$ 时， $H_2\left(\delta ,\omega \right)$ 达到极大值或极小值，对应的函数 $H\left(\delta ,\omega \right)$ 会达到其波峰或波谷值，且经过横轴 $\omega =0$ 时，功角 $\delta$ 的最大值为 ${\delta }_u$ ，因为如果功角 $\delta$ 一旦超过了 ${\delta }_u$ ，系统将失去稳定，沿着轨线函数 $H\left(\delta ,\omega \right)$ 将不体现振荡特性，因此，在稳定域 $R_3$ 内， $\delta \le {\delta }_u$ ，且点 $C\left({\delta }_u,0\right)$ 位于稳定边界上，如图4(b)所示。稳定域 $R_3$ 内，在点C处函数 $H\left(\delta ,\omega \right)$ 达到其最大值，即图4(a)中的点c对应的纵坐标 $H_c$ ，沿着其他任何轨线函数 $$ 在振荡过程中的最大值都小于 $H_c$ 。而当轨线经过纵轴 $\delta ={\delta }_s$ 时，对应的函数 $H\left(\delta ,\omega \right)$ 会达到其过零点。沿着稳定边界，函数 $H\left(\delta ,\omega \right)$ 会呈现图4(a)所示的振荡特性，因此稳定边界上必定存在一点，其对应的函数 $H\left(\delta ,\omega \right)$ 值等于0，即图4(b)中的点A。当轨线从图4(b)中稳定边界上的点A经过一段时间运动到点C的位置，相应的函数值 $H\left(\delta ,\omega \right)$ 从0变化到 $H_c$ ，变化过程对应图4(a)中的从点a变化到点c的过程。

图4. 特殊点的确定

仅采用A，C两点无法准确地描述从点A到点C之间的稳定边界，因此考虑取A，C两点横坐标的中点 $\delta =\frac{{\delta }_s+{\delta }_u}{2}=\frac{\text{π}}{2}$ ，将A到C之间的部分稳定边界上横坐标为 $\delta =\frac{\text{π}}{2}$ 的点记作点B，如图4(b)所示，其函数值 $H\left(\delta ,\omega \right)$ 对应图4(a)中的点b的纵坐标。利用3个特殊点 $A,B,C$ ，对局部稳定边界逼近。当然，若要提高近似程度，可以在AC之间的稳定边界上求取更多的特殊点，本文不做详细分析。

下面介绍特殊点 $A,B,C$ 的求取过程。

点B到点C之间 $\omega$ 的变化过程可以看成是关于功角 $\delta$ 的函数，即 $\omega =\omega \left(\delta \right)$ 。函数 $\omega \left(\delta \right)$ 满足2个条件：(1) $\omega \left(\delta \right)$ 经过点C(2) $\omega \left(\delta \right)$ 满足转子运动方程。为了求解方便，将点C转化为原点，即令

$x=\delta -{\delta }_u$ (3)

则系统运动方程转变为下式所示。

$\{\begin{array}{l}\stackrel{\dot{}}{x}=\omega \\ M\stackrel{\dot{}}{\omega }=P_M-D\omega -P_{em}\sin \left(x+{\delta }_u\right)\end{array}$ (4)

在 $\left(x,\omega \right)$ 平面上， $\omega$ 的变化过程用数学表达式表示如下：

$\omega =\omega \left(x\right)$ (5)

$\omega$ 对时间的导数可以表示如下式所示。

$\stackrel{\dot{}}{\omega }=\frac{\partial \omega }{\partial x}\stackrel{\dot{}}{x}$ (6)

将式(7)代入式(5)，可得下式

$P_M-D\omega -P_{em}\sin \left(x+{\delta }_u\right)=M\frac{\partial \omega }{\partial x}\stackrel{\dot{}}{x}$ (7)

在 $\left(x,\omega \right)$ 平面的点 $C\left(0,0\right)$ 附近对 $\omega \left(x\right)$ 进行二阶泰勒展开：

$\omega \left(x\right)=a_{1}x+\frac{1}{2}{a}_2{x}^2$ (8)

将式(9)代入式(8)中求解，可以求得各次项系数 $a_1$ 、 $a_2$ ，于是可得 $\omega \left(x\right)$ 的二阶近似代数表达式。将点B的横坐标 $x=\frac{\text{π}}{2}-{\delta }_u$ 代入式(9)，可以求得点B的纵坐标 ${\omega }_B$ ，其数学表达式如下式所示。

${\omega }_B=a_1\left(\frac{\text{π}}{2}-{\delta }_u\right)+\frac{1}{2}{a}_2{\left(\frac{\text{π}}{2}-{\delta }_u\right)}^2$ (9)

于是，求得特殊点B的坐标为 $\left(\frac{\text{π}}{2},{\omega }_B\right)$ 。这里值得注意，若想得到精度更高的点B的坐标，可以对 $\omega \left(x\right)$ 在点 $C\left(0,0\right)$ 附近进行三阶及以上的泰勒展开，相应的求解过程将更加复杂，此处不再详细描述。

求得特殊点B后，对特殊点A进行求解。类似过程与上述类似，在此不再赘述。

特殊点 $A,B,C$ 都位于稳定域 $R_3$ 边界上，求得特殊点 $A,B,C$ 后，考虑稳定边界必定是平滑的，可以用平滑的曲线将特殊点 $A,B,C$ 连接起来，作为局部稳定边界的近似。

因此，系统稳定域的一种估计可给出如下：

$R_{est}=R_1\cup R_2\cup R_3$

相比于传统的李雅普诺夫方法利用 $R_1$ 部分直接估计稳定域，稳定域的估计范围得到了扩大，计算精度提高了。

4. 仿真验证

本文采用文献 [17] 中的单机—无穷大系统为例，验证本文所提出的稳定域边界分析方法。初始条件 ${\omega }_0=2\text{π}f,{\delta }_0=0.6027\text{rad}$ ； $t=0\text{s}$ 时双回输电线中的一回线的始端发生两相短路接地，故障期间 $x_2=2.8,P_{em2}=0.504$ ；故障切除后 $x_3=1.041,P_{em3}=1.35$ ；其他参数为 $T_J=8.18,P_m=1$ 。

在 $\left(\delta ,\omega \right)$ 平面，利用本文的方法求得故障后的局部稳定域边界如图5(a)所示。

图5(a)为利用本文方法求得的单机系统稳定域边界近似示意图。在不稳定平衡点附近对稳定域边界进行了局部近似，图中点A、B是利用二阶泰勒展开方法求得的特殊点，进而得到局部稳定边界，如图中实线1所示。点A’、B’是利用三阶泰勒展开方法求得的特殊点，进而得到局部稳定边界，如图中虚线2所示。另外，不计阻尼时的稳定域被包含在计及阻尼的稳定域近似边界内部，这与通常的认识也是一致的。

图5. 单机无穷大系统的稳定边界近似及其故障示意图

利用特殊点求得的局部稳定域边界，对故障的临界切除时间进行估计，并依据此指标，比较局部稳定边界1和局部稳定边界2对真实稳定边界的逼近效果。临界切除时间由持续故障轨迹与稳定域近似边界的交点确定，交点所对应的时间可以近似为该故障的临界切除时间。图5(b)给出了利用稳定域边界近似来计算临界切除时间的机理。

利用近似边界求得的临界故障切除时间结果如表1所示。

表1. 临界切除时间结果比较

由表1可知，利用时域仿真法所得的临界切除时间为0.2456s~0.2457 s，可作为真实稳定边界对应的故障临界切除时间。利用故障后局部稳定边界1所得的故障临界切除时间近似则为0.2546~0.2547 s，误差为3.66%，且偏乐观；而采用局部稳定边界2所得的故障临界切除时间近似则为0.2458s~0.2459 s，误差为0.08%，且偏乐观；而基于李雅普诺夫直接法求得的稳定域，求得的临界切除时间近似则为0.2194s~0.2195 s，该近似误差为−10.67%，且过于保守。

从表1可以看出：根据对故障临界切除时间的估计，本文所提方法相比于李雅普诺夫直接法，求得的稳定域近似边界对真实稳定域边界的逼近效果更好，且误差大大减小了。此外，在求解特殊点的过程中，采用的泰勒展开级数越高，求得的稳定域精度越高。

5. 结语

本文针对李雅普诺夫直接法求取电力系统稳定域的保守性，对能量函数V的其他2种稳定形式进行了研究分析，在此基础上，提出了一种基于并集及特殊点求解的稳定域边界分析方法，将各种稳定形式对应的稳定域的并集作为系统的稳定域，并且在求解过程中利用特殊点法对电力系统稳定域边界进行逼近，仿真算例验证了该方法能够正确估计系统的稳定域边界，且极大地改善了李雅普诺夫直接法估计稳定域的保守性。

文章引用

黄 罡. 基于并集及特殊点求解的稳定域边界分析方法
Stability Region Boundary Analysis Method Based on Union Set and Special Points Solution[J]. 智能电网, 2019, 09(04): 164-171. https://doi.org/10.12677/SG.2019.94018

参考文献