CN104933639B - 一种针对大规模电力系统的小干扰稳定性快速分析方法 - Google Patents

Abstract

本发明公开了一种针对大规模电力系统的小干扰稳定性快速分析方法，该方法首先获取电力系统网络及网络参数、动态元件模型及模型参数，形成系统的微分代数方程组DAE，在系统稳态点(x₀,y₀)处对DAE进行线性化，得到对应的线性动态系统，然后并行计算系统关键振荡模态特征值，最后基于得到的关键振荡模态特征值计算结果得到系统关键振荡模态的阻尼比。通过与阻尼比临界值ζ₀对比判断系统的稳定性；若存在非稳定振荡模态，可根据计算得到的特征值及其特征向量计算系统系统变量对该模态的参与因子，进而确定与非稳定振荡模态强相关的动态元件。本方法能够对大规模电力系统的小干扰稳定性进行快速分析，满足实时性要求。

Description

一种针对大规模电力系统的小干扰稳定性快速分析方法

技术领域

本发明属于电力系统的稳定分析与控制技术领域，尤其涉及一种大规模电力系统小干扰稳定性的快速分析方法。

背景技术

小干扰稳定是电力系统安全稳定运行的基础，小干扰稳定分析是电力系统运行分析与控制的基本和重要功能模块。随着智能电网的发展，新能源并网发电容量和新型负荷接入容量日趋增大，对电力系统的运行带来了新的特点，其中主要包括系统运行波动性和随机性的增强。这些特点要求在更短的时间尺度内，对系统的小干扰稳定性进行判定和分析，对可能发生的低频振荡现象进行及时预警，并采取预防校正措施。

同时，随着区域间互联电力系统规模不断扩大以及接入系统的元件类型日趋丰富，电力系统的动态行为建模及其稳定性分析和控制成为极具挑战性的难题。在电力系统的小干扰稳定分析领域，目前实际电力系统动态模型微分方程的数量多达10⁴及以上。对于如此庞大的系统，基于全特征值理论的小干扰稳定分析方法已难以适应于小干扰稳定分析对于实时性的要求。对于采用部分特征值算法的小干扰稳定分析方法，在大规模系统的实际应用中，由于计算复杂度高，同样难以满足稳定分析的实时性要求。另外，对于采用部分特征值算法的分析方法，在实施过程中，可能存在遗漏振荡模态、算法收敛性差等问题，导致产生不可靠的稳定性分析结果，难以应用于生产实践。

本发明基于已有的关于部分特征值算法及其在电力系统小干扰稳定分析中的应用研究，提出针对大规模电力系统的小干扰稳定性快速分析方法。该方法可灵活采用两种常用的谱变换预处理技术和隐式重启动Arnoldi、Krylov-Schur算法，同时利用变换点在复平面上位置的相互独立性采用并行计算技术加速关键振荡模态特征值的求解。以上特点使得本发明提出的方法在获得准确可靠振荡模态信息的同时，能够显著降低分析计算时间，满足大规模电力系统小干扰稳定在线分析的要求。

发明内容

本发明的目的在于针对现有技术的不足，提出一种针对大规模电力系统的小干扰稳定性快速分析方法。

本发明的目的是通过以下技术方案来实现的：一种针对大规模电力系统的小干扰稳定性快速分析方法，该方法包括以下步骤：

第一步：获取电力系统网络及网络参数、动态元件模型及模型参数，形成系统的微分代数方程组(DAE)，在系统稳态点(x₀,y₀)处对DAE进行线性化，得到对应的线性动态系统。

第二步：并行计算系统关键振荡模态特征值：首先确定系统临界阻尼比值ζ₀、振荡模态最大频率f₁、振荡模态最小频率f₂，通过以上3个参数确定待求解特征值所在的复平面区域。根据该复平面区域和计算线程数量，确定用于Cayley变换的位移点、Shift-Invert变换的位移点集合及每个位移点对应的特征值搜索区域。在Cayley变换的位移点进行Cayley变换，在Shift-Invert变换的位移点进行Shift-Invert变换后，启用计算线程使用部分特征值计算方法计算位移点搜索区域内的目标特征值及其左右特征向量。

第三步：基于第二步得到的关键振荡模态特征值计算结果，得到系统关键振荡模态的阻尼比。通过与阻尼比临界值ζ₀对比，若所有关键振荡模态阻尼比大于临界值，则系统小干扰稳定，所有关键振荡模态为稳定模态。若存在非稳定振荡模态，可根据计算得到的特征值及其特征向量计算系统系统变量对该模态的参与因子，进而确定与非稳定振荡模态强相关的动态元件。

进一步地，所述步骤1中，所述系统微分代数方程组(DAE)具有如下形式：

其中，x为系统状态变量，y为系统代数变量，f为系统微分方程组，g为系统代数方程组。

在系统平衡点(x₀,y₀)处对上述DAE方程组线性化，得到线性动态系统的状态空间模型具有如下形式：

其中Δx为系统状态变量增量，Δy为系统代数变量增量，该线性动态系统的稳定性质由状态矩阵A_S＝ A-BD^-1C的特征值表征。

进一步地，所述步骤2中，所述计算线程使用部分特征值计算方法为隐式重启动Arnoldi算法、Krylov-Schur算法或Inexact Jacobi-Davidson算法。

进一步地，所述步骤2中，所述ζ₀取值为3％-15％，f₁取值为0Hz-0.1Hz,f₂取值为2.5Hz-5Hz。

进一步地，所述步骤3中，所述的阻尼比由下式确定：

其中，σ为特征值实部，ω为特征值虚部。

进一步地，所述步骤3中，所述的参与因子第i个元素p_i由下式计算：

p_i＝u_iv_i，

式中，u_i为左特征向量第i个元素，v_i为右特征向量第i个元素，p_i越大表示第i个状态变量与该振荡模态相关性越强。

本发明的有益效果是，本发明提出了针对大规模电力系统的小干扰稳定性快速分析方法。与已有的技术相比，本发明提出的方法主要有以下改进：

1、该方法结合了两种成熟高效的谱变换处理技术，Cayley变换和Shift-Invert变换，有效加速了系统关键振荡模态特征值的收敛速度，使得小干扰稳定分析结果具有可靠性；

2、该方法提出的并行计算框架基于谱变换技术中所采用变换点的相互独立性，对每个变换点相应复平面区域内的特征值进行并行求解，从而显著降低大规模电力系统小干扰稳定分析所需计算时间，满足实时性要求。

附图说明

图1是大规模电力系统小干扰稳定分析中关键振荡模态特征值的并行计算方法的流程图；

图2是电力系统小干扰稳定分析中关键振荡模态特征值的分布区域；

图3是Shift-Invert变换点及每个变换点对应的特征值搜索区域；

图4是本方法采用的并行计算框架示意图；

图5是本方法所采用的并行计算框架中管理计算线程(主线程)的工作流程图；

图6是本方法应用于一个实际电力系统的小干扰稳定分析结果。

具体实施方式

针对大规模电力系统的小干扰稳定性快速分析方法包括如下步骤：

第一步：获取电力系统网络及网络参数、动态元件模型及模型参数，形成系统的微分代数方程组(DAE)，在系统稳态点(x₀,y₀)处对DAE进行线性化，得到对应的线性动态系统。

所述系统微分代数方程组(DAE)具有如下形式：

其中，x为系统状态变量，y为系统代数变量，f为系统微分方程组，g为系统代数方程组。

在系统平衡点(x₀,y₀)处对上述DAE方程组线性化，得到线性动态系统的状态空间模型具有如下形式：

其中Δx为系统状态变量增量，Δy为系统代数变量增量，该线性动态系统的稳定性质由状态矩阵A_S＝ A-BD^-1C的特征值表征。

第二步：并行计算系统关键振荡模态特征值。

在该步骤中，首先确定系统临界阻尼比值ζ₀、振荡模态最大频率f₁、振荡模态最小频率f₂，通过以上3个参数确定待求解特征值所在的复平面区域。根据该复平面区域和计算线程数量，确定用于Cayley变换的位移点、Shift-Invert变换的位移点集合及每个位移点对应的特征值搜索区域。在Cayley变换的位移点进行Cayley变换，在Shift-Invert变换的位移点进行Shift-Invert变换后，启用计算线程使用部分特征值计算方法(包括但不限于隐式重启动Arnoldi算法、Krylov-Schur算法、Inexact Jacobi-Davidson算法)计算位移点搜索区域内的目标特征值及其左右特征向量。

ζ₀取值为3％-15％，f₁取值为0Hz-0.1Hz,f₂取值为2.5Hz-5Hz，以上三个参数的具体取值由小干扰稳定分析需求确定。

所述的Cayley变换具有如下形式：

M_s＝(A_s-σ₁I)(A_s-σ₂I)^-1

其中，M_s为谱变换后的状态矩阵，A_s为系统状态矩阵，I为与A_s维度相同的单位对角矩阵，σ₁为反变换点，σ₂为变换点。σ₁、σ₂变换点对应的特征值搜索区域为阻尼比小于ζ₁的复平面区域，σ₁与σ₂的中位线与ζ₁确定的射线重合，系统临界阻尼ζ₀与ζ₁满足关系：

ζ₁＜ζ₀

所述的Shift-Invert变换具有如下形式：

M_s＝(A_s-σI)^-1

其中，M_s为谱变换后的状态矩阵，A_s为系统状态矩阵，I为与A_s维度相同的单位对角矩阵,σ为位移变换点，σ对应的特征值搜索区域为复平面上以σ圆心和r为半径的圆形区域。所有Shift-Invert变换点的搜索区域能够完整覆盖以ζ₀、ζ₁、f₁、f₂为边界的复平面区域。

第三步：基于第二步得到的关键振荡模态特征值计算结果，得到系统关键振荡模态的阻尼比。通过与阻尼比临界值ζ₀对比，若所有关键振荡模态阻尼比大于临界值，则系统小干扰稳定，所有关键振荡模态为稳定模态。若存在非稳定振荡模态，可根据计算得到的特征值及其特征向量计算系统系统变量对该模态的参与因子，进而确定与非稳定振荡模态强相关的动态元件(发电机组、感应电动机负荷等)。

所述的阻尼比由下式确定：

上式中，σ为特征值实部，ω为特征值虚部。

所述的参与因子第i个元素p_i由下式计算：

p_i＝u_iv_i

上式中，u_i为左特征向量第i个元素，v_i为右特征向量第i个元素，p_i越大表示第i个状态变量与该振荡模态相关性越强。

以下结合图表，对本发明的实施例作详细说明，该发明的流程图如图1所示。

实施例：

为了验证本发明所提出的针对大规模电力系统的小干扰稳定性快速分析方法的可用性和计算效率，发明人使用C编程语言和OpenMP并行编程模型开发实现了基于Cayley、Shift-Invert谱变换技术和Krylov-Schur算法的大规模电力系统关键振荡模态特征值分析程序，并使用一台装配有64核心Intel Xeon E7-8837和1TB内存的对称多处理机(Symmetric Multi-Processor)完成了本实施例的测试和验证。在本实施例中，第三方软件包KLU、Intel Math Kernel Library被用于实现矩阵分解等基础矩阵运算。表1所示的三个电力系统算例被用于本实施例中的测试与验证。

表1：测试算例的系统参数

本发明所提出的计算方法主要利用了大规模电力系统关键振荡模态特征值分析计算的下述特点：

关键振荡模态特征值的保守分布区域包括阻尼比小于ζ₀的左半复平面和完整的右半复平面，参见附图2的阴影区域A、B。ζ₁右侧的关键振荡模态特征值可通过Cayley变换映射为单位圆外的主导特征值，其余区域的特征值通过多变换点的Shift-Invert变换映射为相应谱平面的主导特征值。Krylov-Schur算法优先收敛于主导特征值的性质使得谱变换技术能够有效加速目标特征值的收敛。另外，由于与不同变换点对应的特征值求解具有良好的解耦性(参见附图3)，通过并行计算技术(本实施例中采用OpenMP并行编程模型)可有效加速问题求解。同时，本发明提出的并行计算框架易于应用于其它部分特征值算法，例如隐式重启动Arnoldi、非精确Jacobi-Davidson算法。

表2给出了本发明所提出方法应用于表1中3个测试系统的计算耗时结果。表2的数据说明本发明所提出的小干扰稳定性快速分析方法能够获得良好的并行计算加速比，具有良好的并行加速性能。该优势主要归因于以下3个性质：

性质一：该方法充分利用了Cayley、Shift-Invert预处理技术中变换点位置的独立性，在所提出的并行计算框架下可有效挖掘部分特征值算法的并行计算应用潜力；

性质二：应用Cayle、Shift-Invert预处理后，所有待求解的关键振荡模态特征值分区域映射为主导特征值，对子空间类部分特征值算法(例如，Krylov-Schur算法、隐式重启动Arnoldi)的收敛性有明显增强作用；

性质三：在具体程序实现中，充分利用系统状态矩阵的稀疏性质，可有效减少内存需求，提高计算效率。

表2：并行计算时间及加速比

与此同时，为了说明本发明提出的方法计算结果的准确性，我们针对三个测试算例的关键振荡模态特征值残差进行了检验，结果如表3所示。

表3：关键振荡模态的特征值残差

由表3和附图6中的计算结果可得，本发明所提出的针对大规模电力系统的小干扰稳定性快速分析方法的可靠性(关键振荡模态特征值的准确性)在启用不同数量线程的情况下保持一致，进一步表明该方法具有良好的并行扩展性能。附图6的计算结果说明，测试系统3在运行点处存在3个负阻尼振荡模态(对应特征值位于右半复平面)，系统处于小干扰不稳定状态，需要采取稳定控制措施。

综上所述，本发明提出的针对大规模电力系统的小干扰稳定性快速分析方法能够高效可靠地计算系统关键振荡模态信息，可以对大规模电力系统的小干扰稳定性进行快速分析，满足在线应用的实时性要求。

Claims (6)

1.一种针对大规模电力系统的小干扰稳定性快速分析方法，其特征在于，该方法包括以下步骤：

第一步：获取电力系统网络及网络参数、动态元件模型及模型参数，形成系统的微分代数方程组(DAE)，在系统稳态点(x₀,y₀)处对DAE进行线性化，得到对应的线性动态系统;

第二步：并行计算系统关键振荡模态特征值：首先确定系统临界阻尼比值ζ₀、振荡模态最大频率f₁、振荡模态最小频率f₂，通过以上3个参数确定待求解特征值所在的复平面区域;根据该复平面区域和计算线程数量，确定用于Cayley变换的位移点、Shift-Invert变换的位移点集合及每个位移点对应的特征值搜索区域; 在Cayley变换的位移点进行Cayley变换，在Shift-Invert变换的位移点进行Shift-Invert变换后，启用计算线程使用部分特征值计算方法计算位移点搜索区域内的目标特征值及其左右特征向量;

第三步：基于第二步得到的关键振荡模态特征值计算结果，得到系统关键振荡模态的阻尼比; 通过与阻尼比临界值ζ₀对比，若所有关键振荡模态阻尼比大于临界值，则系统小干扰稳定，所有关键振荡模态为稳定模态; 若存在非稳定振荡模态，可根据计算得到的特征值及其特征向量计算系统系统变量对该模态的参与因子，进而确定与非稳定振荡模态强相关的动态元件。

2.根据权利要求1所述的方法，其特征在于，所述步骤1中，所述系统微分代数方程组(DAE)具有如下形式：

其中，x为系统状态变量，y为系统代数变量，f为系统微分方程组，g为系统代数方程组;

在系统平衡点(x₀,y₀)处对上述DAE方程组线性化，得到线性动态系统的状态空间模型具有如下形式：

其中Δx为系统状态变量增量，Δy为系统代数变量增量， ; 该线性动态系统的稳定性质由状态矩阵A_S=A-BD^-1C的特征值表征。

3.根据权利要求1所述的方法，其特征在于，所述步骤2中，所述计算线程使用部分特征值计算方法为隐式重启动Arnoldi算法、Krylov-Schur算法或Inexact Jacobi-Davidson算法。

4.根据权利要求1所述的方法，其特征在于，所述步骤2中，所述ζ₀取值为3%-15%，f₁取值为0Hz-0.1Hz,f₂取值为2.5Hz-5Hz。

5.根据权利要求1所述的方法，其特征在于，所述步骤3中，所述的阻尼比由下式确定：

其中，σ为特征值实部，ω为特征值虚部。

6.根据权利要求1所述的方法，其特征在于，所述步骤3中，所述的参与因子第i个元素p_i由下式计算：

p_i＝u_iv_i，

式中，u_i为左特征向量第i个元素，v_i为右特征向量第i个元素，p_i越大表示第i个状态变量与该振荡模态相关性越强。