CN108808703A - 基于低阶igd-irk的时滞电力系统小干扰稳定性分析方法 - Google Patents

Abstract

本发明公开了基于低阶IGD‑IRK的时滞电力系统小干扰稳定性分析方法，包括：建立时滞电力系统数学模型；根据时滞电力系统状态变量是否与时滞有关，时滞电力系统的微分方程改写为与时滞有关的部分及与时滞无关的部分；获得与重组后的时滞电力系统状态方程相对应的无穷小生成元；针对无穷小生成元基于隐式龙格‑库塔法低阶离散化，得到无穷小生成元低阶离散化近似矩阵；针对无穷小生成元低阶离散化近似矩阵进行位移逆变化，得到无穷小生成元低阶离散化近似逆矩阵；针对上述逆矩阵进行稀疏特征值计算，之后针对稀疏特征值再进行反变换和牛顿校验，得到时滞电力系统的精确特征值。

Description

基于低阶IGD-IRK的时滞电力系统小干扰稳定性分析方法

技术领域

本发明涉及电力系统技术领域，特别是涉及基于低阶IGD-IRK的时滞电力系统小干扰稳定性分析方法，IGD-IRK为英文“Infinitesimal Generator Discretization WithImplicit Runge-Kutta”的缩写，中文含义为：无穷小生成元隐式龙格-库塔离散化。

背景技术

当能源与电力的需求持续增加，建立跨区、跨国的大型互联电网是当代电力发展趋势。然而在互联初期，复杂的电气结构和薄弱的输电环节使得电力系统更容易出现区域间低频振荡现象。区间振荡的范围较广，振荡的机群之间的联系较为复杂，是制约系统小干扰稳定性的关键。以本地信号为输入的电力系统稳定器(Power System Stabilizer，PSS)能很好地抑制局部振荡模式，但是难以抑制区间振荡，从而难以保证大规模互联电网的稳定性。

广域测量系统(Wide-Area Measurement System，WAMS)的出现和发展，能够实时远距离采集动态参数，使得采用全局的、远端的信息实现电力系统稳定控制成为可能。广域测量系统主要由三部分组成，相量测量单元(Phasor Measurement Unit，PMU)，通信网络，以及监控系统。PMU同步测量系统各枢纽点的状态量，通过通信系统通道传送到监控系统，实现系统的保护控制与检测。经过分析处理的量测数据，不仅能检测系统低频振荡信息，还可以提升系统的阻尼水平和远距离大容量输电能力。

但是，广域量测数据在广域量测系统由不同的通信介质(如光纤、数字微波、电力线等)组成的通信系统中传输和处理时，存在明显的时延现象。时滞在控制系统中引入了一个滞后的相位，可能会使得系统的控制律失效，恶化系统的性能，甚至导致系统从小扰动稳定变为不稳定。因此分析时滞对于电力系统小干扰稳定性的影响，对于设计能够有效抑制区域间低频振荡的广域阻尼控制器，提升电力系统的安全、稳定运行水平，具有重要的理论意义和应用价值。

在现代电力系统中，时滞电力系统小干扰稳定性分析的关键就是关注系统的机电振荡问题。常用的时滞电力系统稳定性分析方法有函数变换法，时域法，预测补偿法等等。中国发明专利基于Padé近似的时滞电力系统特征值计算与稳定性判别方法，申请号为201210271783.8:[P].利用Padé近似多项式逼近时滞环节，进而计算系统最右侧的关键特征值，并判断系统的时滞稳定性。而特征分析法是电力系统小干扰稳定性分析的基本而有效的方法，已经形成了比较成熟和完善的理论。计算得到时滞电力系统的特征值，就可以运用经典的特征分析的思路和理论框架来分析系统小干扰稳定性，进行广域阻尼控制器的优化设计。

近年来，基于谱离散化的时滞系统部分特征值计算方法，开始用于时滞电力系统的稳定性分析。中国发明专利基于EIGD的大规模时滞电力系统特征值计算方法.201510055743.3[P].提出了一种基于显示IGD(Explicit Infinitesimal GeneratorDiscretization,EIGD)的大规模时滞电力系统特征值计算。利用计算得到时滞电力系统的特征值，进而判定系统的小干扰稳定性。

然而，上述时滞电力系统谱离散化计算方法生成的离散化矩阵维数很大，尽管采用稀疏算法进行计算，使其适用于大规模系统的特征值计算，但是矩阵固有的维数问题没有得到解决。

发明内容

为了解决现有技术的不足，本发明提供了基于低阶IGD-IRK的时滞电力系统小干扰稳定性分析方法，用以进行大规模时滞电力系统的稳定性分析，本发明采用的低阶IGD-IRK算法，通过移去与时滞无关的不需要离散化的状态变量，实现对离散化矩阵的降阶，降低了离散化矩阵的维数，解决了矩阵维数较大对于特征值计算效率的影响，从而使算法能高效的计算出时滞电力系统的最右侧的关键特征值。

基于低阶IGD-IRK的时滞电力系统小干扰稳定性分析方法，包括：

建立时滞电力系统数学模型，线性化处理得到时滞电力系统的微分方程；

根据时滞电力系统状态变量是否与时滞有关，时滞电力系统的微分方程改写为与时滞有关的部分及与时滞无关的部分，继而基于与时滞的相关性重组状态变量，获得重组后的时滞电力系统状态方程；

获得与重组后的时滞电力系统状态方程相对应的无穷小生成元，依据无穷小生成元谱映射原理，将计算时滞电力系统机电振荡模式的问题转化为计算无穷小生成元的特征值问题；

针对无穷小生成元基于隐式龙格-库塔法低阶离散化，得到无穷小生成元低阶离散化近似矩阵；

针对无穷小生成元低阶离散化近似矩阵进行位移逆变化，得到无穷小生成元低阶离散化近似逆矩阵；

针对上述逆矩阵进行稀疏特征值计算，之后针对特征值再进行反变换和牛顿校验，得到时滞电力系统的精确特征值，特征值λ则对应时滞电力系统的机电振荡模式。

进一步优选的技术方案，经线性化后的时滞电力系统模型为：

式中：n为系统状态变量总数，Δx(t)为t时刻时系统状态变量的增量，为系统状态变量导数的增量，Δx(t-τ_i)为t-τ_i时刻时系统状态变量的增量，τ_i>0为时滞常数，i＝1,…,m，m为时滞的个数，且满足0<τ₁<τ₂<…<τ_i…<τ_max，其中τ_max为最大的时滞，Δx(0)为系统状态变量的初始值，简写为 为稠密的系统状态矩阵，是稀疏的系统时滞状态矩阵。

进一步优选的技术方案，根据时滞电力系统状态变量是否与时滞有关，将时滞电力系统的微分方程改写成两个部分，即分为与时滞无关部分和与时滞相关部分 n₁为与时滞无关部分的状态变量的维数，n₂为与时滞相关部分的状态变量的维数，且满足n₁+n₂＝n，则描述系统动态特性的时滞微分方程重写为：

式(2)中：和分别为与时滞无关部分状态变量导数的增量和与时滞相关部分状态变量导数的增量，Δx⁽¹⁾(t-τ_i)和Δx⁽²⁾(t-τ_i)分别为t-τ_i时刻时与时滞无关部分系统状态变量的增量和与时滞相关部分系统状态变量的增量，和分别是由状态矩阵和经过重写为与时滞有关和无关的部分得到的矩阵，i＝1,…,m：

式中：和是重写后矩阵中的分块矩阵；

与式(2)对应的系统的特征方程表示为：

式中：λ为时滞电力系统的特征值，v为特征值对应的右特征向量。

进一步优选的技术方案，利用无穷小生成元将式(2)在Banach空间映射为抽象柯西方程：

式中：θ∈[-τ_max,0]，为Δx_t的导数。进一步优选的技术方案，时滞电力系统模型的特征值与无穷小生成元特征值之间的关系为：

式中：σ(·)表示无穷小生成元的谱，λ为时滞电力系统的特征值，该式说明时滞电力系统的特征值与无穷小生成元的特征值是一一对应的。

进一步优选的技术方案，基于隐式龙格-库塔离散化基本理论，在区间[-τ_max,0]上，建立Ns+1个离散点的集合Ω_N，N为给定正整数，s为隐式龙格-库塔法的级数，然后在每个离散点处，离散式(5)，估计Δx_t的值，生成无穷小生成元离散化矩阵；

式中：N_i为在第i个子区间以h_i为离散步长得到的离散点θ_j,i的个数，为在时滞区间[-τ_i,-τ_i-1]的离散点θ_j,i的基础上，采用s级IRK法的横坐标对每个子区间作进一步划分，得到的N_is个离散点构成的集合，i＝1,…,m，c_l为s级龙格-库塔法的横坐标；

只离散与时滞有关的部分，由此得与无穷小生成元对应的低阶离散化矩阵矩阵的阶数为(n+Nsn₂)，表示如下：

式中：s为隐式龙格-库塔法的级数，N为给定正整数，是n₂阶单位矩阵，是Nsn₂×n₁的零矩阵，表示克罗内克积运算，是矩阵的第一块行；子矩阵是高度稀疏的矩阵即：

式中：h_i为每个离散化子区间的步长，i＝1,…,m，ω＝A^-11_s∈^s×1， A为隐式龙格-库塔法中RadauⅡA方法的系数，w_i表示矩阵W中的第i列元素，即

进一步优选的技术方案，为了充分利用克罗内克积的性质和状态矩阵的稀疏特性，从而进行稀疏特征值计算，∑_Ns表示为：

式中：为单位向量，为重组后的稠密的状态矩阵，为重组后时滞状态矩阵的后n₂列，分别可以表示如下：

进一步优选的技术方案，利用位移逆变换技术，将模值较小的这部分特征值转换为主特征值，优先计算出系统的关键模态，即，将λ'+s代入式(4)替代λ得位移之后的特征方程，表示如下：

将矩阵中的第一块行的用进行替代，i＝0,1,…,m，得到位移操作后的无穷小生成元的离散化矩阵的逆矩阵表示为：

式(16)中：

式(17)中：是位移变换后得到的矩阵，是位移变换后时滞矩阵的后n₂列；

式中：为位移变换后时滞矩阵的后n₂列、n₂行组成的矩阵。

进一步优选的技术方案，采用隐式重启动Arnoldi算法来求取中的部分特征值，在Arnoldi算法中，利用与向量乘积形成Krylov子空间，设第k个Krylov向量为则第k+1个Krylov向量q_k+1计算如下：

采用诱导降维迭代方法来计算q_k+1，于是，式(19)写为：

式中：是第l次迭代后q_k+1的近似值。进一步优选的技术方案，采用诱导降维迭代方法来计算q_k+1，具体步骤如下：

将中的后(Ns+1)n₂行的元素，按照列的方向重新进行排列，得矩阵即接着，取矩阵Q₁的后Ns列，得到矩阵即最后，式通过克罗内克积的性质计算得到且。

于是，式(20)转换为：

在式(21)中，利用幂法计算矩阵向量乘积

进一步优选的技术方案，经过稀疏特征值计算后得到的的特征值为λ"，则的近似特征值为：

因此通过牛顿迭代法，得到时滞电力系统的精确特征值λ和与特征值对应的特征向量v。

与现有技术相比，本发明的有益效果是：

第一、本发明提出的低阶IGD-IRK算法可以适用于计算实际大规模时滞电力系统的机电振荡模式对应的特征值，充分考虑了电力系统的规模和时滞对于系统的影响。

第二、本发明提出的低阶IGD-IRK算法的优势在于在原有的非低阶的IGD-IRK计算方法的基础之上，对离散化矩阵进行降阶，使其维数大幅减少。因为在降阶过程中不涉及到简化，因此计算出的特征值与降阶前基本相同，在保证准确性的同时，能高效的计算出系统的特征值，减少计算量与计算时间。

第三、本发明的方法是基于低阶IGD-IRK算法的时滞电力系统小干扰稳定性分析方法，核心创新点体现在将以下各种技术结合在一起，采用降阶的思想，将无限维的时滞电力系统转换为维数较低的离散化矩阵，然后采用位移逆变换预处理技术、IRA稀疏特征值算法、迭代求解矩阵逆向量的乘积三个核心技术，高效计算出系统近似特征值，最后经牛顿迭代计算出精确特征值。

附图说明

构成本申请的一部分的说明书附图用来提供对本申请的进一步理解，本申请的示意性实施例及其说明用于解释本申请，并不构成对本申请的不当限定。

图1为低阶IGD-IRK算法中的离散点集合Ω_N；

图2为基于低阶IGD-IRK算法的时滞电力系统小干扰稳定性的特征分析方法的流程图；

图3为16机68节点算例系统。

具体实施方式

应该指出，以下详细说明都是例示性的，旨在对本申请提供进一步的说明。除非另有指明，本文使用的所有技术和科学术语具有与本申请所属技术领域的普通技术人员通常理解的相同含义。

需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本申请的示例性实施方式。如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，还应当理解的是，当在本说明书中使用术语“包含”和/或“包括”时，其指明存在特征、步骤、操作、器件、组件和/或它们的组合。

本申请的整体技术构思为：

首先要建立时滞电力系统的小干扰稳定性分析的线性化模型，将模型在稳态运行点附近线性化，由此可得用于时滞电力系统的线性化模型。

低阶算法的原理是对系统状态变量分成两个部分：与时滞有关的部分和与时滞无关的部分。然后根据对于系统状态变量的划分，对时滞电力系统的时滞微分方程进行重组，微分方程也可以重组为与时滞有关和无关两个部分。

通过谱映射定理，即将时滞系统的全部特征值转换为无穷维Banach(巴拿赫)空间时滞系统算子的谱。也就是说，求解无穷小生成元的特征值即为求解时滞电力系统的关键特征值。

依据低阶算法的思想，因为与时滞无关的部分对之前的系统状态无任何影响，所以可以认为与时滞无关部分的离散化是没有必要的。基于隐式龙格-库塔法的离散化理论对进行离散化，只将与时滞相关的部分提取出来，得到低阶的近似矩阵的阶数为(n+Nsn₂)，其维数与无时滞电力系统的矩阵维数近似，所以时滞电力系统的关键特征值可以有效率地被计算出来。

为了充分利用离散化矩阵和系统增广状态矩阵的稀疏特性，需要对离散化矩阵的第一块行显示重组为克罗内克积的形式，以便于稀疏计算的时候减少计算时间与计算量。

想要得到的低频振荡模式一般位于复平面虚轴附近，通过位移逆变换预处理技术，这部分模值较小的特征值转换为主特征值。

对于所述时滞电力系统离散化近似矩阵的阶数为(n+Nsn₂)，对于大规模系统，通常采用隐式重启动Arnoldi算法等稀疏算法求解部分特征值。在稀疏计算时，充分利用克罗内克积的性质以及系统状态矩阵的稀疏特性。最后，利用牛顿校验进行修正，得到时滞电力系统的精确特征值。

本申请的一种典型的实施方式中，如图1所示，关于低阶IGD-IRK算法是一种分段离散化方案，需要建立离散点的集合Ω_N。首先，在每个时滞区间[-τ_i,-τ_i-1]划分N_i个子区间；然后，用s级龙格-库塔法的横坐标对每个子区间做进一步划分；最终，得到具有Ns+1个离散点的集合Ω_N。为隐式龙格-库塔离散化方案奠定基础。

图3为低阶IGD-IRK算法的时滞电力系统小干扰稳定性分析方法，包括如下步骤：

S1：建立时滞电力系统线性化模型，得到系统的微分方程。

S2：根据系统状态变量是否与时滞有关，将系统的微分方程改写重组为两个部分，从而系统状态矩阵也相应地进行改写与重组；然后基于无穷小生成元谱映射原理，将计算时滞电力系统机电振荡模式的问题转化为计算无穷小生成元的特征值问题。

S3：基于隐式龙格-库塔离散化方案的理论，建立Ns+1个离散点的集合(其中，N为给定正整数，s为隐式龙格-库塔法的级数)；只针对与时滞有关的部分，对无穷小生成元进行离散化，除去与时滞无关的部分，从而生成无穷小生成元低阶的离散化矩阵，并将矩阵的第一块行改写成克罗内克积的形式。

S4：利用三种核心技术改进算法使其适用于大规模时滞电力系统，计算出系统的低频振荡模式，三种核心技术改进算法分别为：位移逆预处理技术将阻尼较小的特征值转换为主特征值；IRA算法进行稀疏特征值计算以及IDR(s)方法计算矩阵逆向量的乘积运算。

S5：从步骤S4中得到特征值λ″之后，经过反变换和牛顿校验，得到时滞电力系统的精确特征值λ，特征值λ则对应时滞电力系统的机电振荡模式。

所述步骤S1中，经线性化后的时滞电力系统模型如下：

式中：n为系统状态变量总数。Δx(t)为t时刻时系统状态变量的增量，为系统状态变量导数的增量，Δx(t-τ_i)为t-τ_i时刻时系统状态变量的增量。τ_i>0(i＝1,…,m)为时滞常数，m为时滞的个数，且满足0<τ₁<τ₂<…<τ_i…<τ_max，其中τ_max为最大的时滞。Δx(0)为系统状态变量的初始值，简写为 为稠密的系统状态矩阵。是稀疏的系统时滞状态矩阵。

所述步骤S2中，根据系统状态变量是否与时滞有关，可以将系统的微分方程改写成两个部分，即分为与时滞无关部分和与时滞相关部分n₁为与时滞无关部分的状态变量的维数，n₂为与时滞相关部分的状态变量的维数，且满足n₁+n₂＝n。则描述系统动态特性的时滞微分方程可重写为：

式(2)中：和分别为与时滞无关部分状态变量导数的增量和与时滞相关部分状态变量导数的增量，Δx⁽¹⁾(t-τ_i)和Δx⁽²⁾(t-τ_i)分别为t-τ_i时刻时与时滞无关部分系统状态变量的增量与时滞相关部分系统状态变量的增量。和分别是由状态矩阵和经过重写为与时滞有关和无关的部分得到的矩阵：

式中：和是重写后的分块矩阵。

与式(2)对应的系统的特征方程可以表示为：

式中：λ为特征值，v为特征值对应的右特征向量。

利用无穷小生成元将式(2)在Banach空间映射为抽象柯西方程：

式中：θ∈[-τ_max,0]。为Δx_t的导数。

所述时滞电力系统模型的特征值与无穷小生成元特征值之间的关系为：

式中：σ(·)表示无穷小生成元的谱，λ为电力系统的特征值。该式说明时滞电力系统的特征值即为无穷小生成元的特征值。

因为对于的特征值的求解是无穷维问题，需要采用隐式龙格-库塔离散化方法对进行离散化，得到有限维近似矩阵计算其特征值。

所述步骤S3：基于隐式龙格-库塔离散化基本理论，在区间[-τ_max,0]上，建立Ns+1个离散点的集合Ω_N。然后在每个离散点处，离散式(5)，估计Δx_t的值，生成无穷小生成元离散化矩阵。矩阵的维数会因为离散化方案的不同，离散点集合的不同，具有不同的维数以及不同的结构。

式中：N_i为在第i个子区间以h_i为离散步长得到的离散点θ_j,i的个数，在时滞区间[-τ_i,-τ_i-1]的离散点θ_j,i的基础上，采用s级IRK法的横坐标对每个子区间作进一步划分，得到的N_is个离散点构成的集合，i＝1,…,m。c_l为s级龙格-库塔法的横坐标。

根据降阶的基本思想，只离散与时滞有关的部分，由此可以得与无穷小生成元对应的低阶离散化矩阵矩阵的阶数为(n+Nsn₂)，表示如下：

式中：s为隐式龙格-库塔法的级数，N为给定正整数，是n₂阶单位矩阵，是Nsn₂×n₁的零矩阵，表示克罗内克积运算。是矩阵的第一块行；子矩阵是高度稀疏的矩阵即：

式中：h_i(i＝1,…,m)为每个离散化子区间的步长，ω＝A^-11_s∈^s×1，1_s＝[1,…,1]^T∈^s×1，W＝-A^-1∈^s×s。A为隐式龙格-库塔法中RadauⅡA方法的系数。w_i表示矩阵W中的第i列元素，即W＝-A^-1＝[w₁,w₂,…,w_s]∈^s×s。

为了充分利用克罗内克积的性质和状态矩阵的稀疏特性，从而进行稀疏特征值计算，∑_Ns可以表示为：

式中：为单位向量，为重组后的稠密的状态矩阵，为重组后时滞状态矩阵的后n₂列，分别可以表示如下。

所述步骤S4中，首先利用位移逆变换技术，可以将模值较小的这部分特征值转换为主特征值，优先计算出系统的关键模态。即，将λ'+s代入式(6)替代λ可以得位移之后的特征方程，表示如下：

将矩阵中的第一块行的用进行替代，可以得到位移操作后的无穷小生成元的离散化矩阵的逆矩阵可以表示为：

式(16)中：

式(17)中：是位移变换后得到的矩阵，是位移变换后时滞矩阵的后n₂列；

式中：为位移变换后的时滞矩阵的后n₂列和n₂行。

然而，传统的特征值计算方法并不能适用于大规模时滞电力系统，因此需要采用隐式重启动Arnoldi算法(Implicitly Restarted Arnoldi algorithm,IRA)来求取中的部分特征值。在IRA算法中，计算量最大的操作是利用与向量乘积形成Krylov子空间。设第k个Krylov向量为则第k+1个Krylov向量q_k+1可计算如下：

为避免直接求解利用诱导降维(Induced Dimension Reduction,IDR(s))迭代方法来计算q_k+1。于是，式(19)可写为：

式中：是第l次迭代后q_k+1的近似值。

迭代求解的好处是保持了的稀疏特性，具体步骤如下。

将中的后(Ns+1)n₂行的元素，按照列的方向重新进行排列，可以得矩阵即接着，取矩阵Q₁的后Ns列，得到矩阵即最后，式可以通过克罗内克积的性质计算得到。于是，式(20)可转换为：

在式(21)中，可以利用幂法计算矩阵向量乘积

所述步骤S5中，经过稀疏特征值计算后得到的的特征值为λ"，则的近似特征值为：

因此通过牛顿迭代法，得到时滞电力系统的精确特征值λ和与特征值对应的特征向量v。

为了更好的说明本申请技术方案的有效性，下面给出了更为详细的算例分析

利用16机68节点算例系统，来验证本发明提出的基于低阶IGD-IRK算法的时滞电力系统小干扰稳定性的特征分析方法的有效性。所有的分析均在Matlab中和在Inter3.4GHz 16GB RAM台式计算机上进行。IRA、IDR(s)和牛顿法的收敛精度均为10^-6。IDR(s)中的参数s选择为4。

16机68节点算例系统的单线图如图3所示。特征值分析表明，系统存在两个弱阻尼区间低频振荡模式这两个模式分别表现为G₁-G₁₃相对于G₁₂-G₁₆之间的振荡、G₁-G₉相对于G₁₀-G₁₃之间的振荡。为增强模式阻尼，分别在发电机G₂和G₅上安装广域阻尼控制器。反馈信号分别取为G₂和G₁₅的相对转速，G₅和G₁₃的相对转速。控制器参数分别为：K_{s_1}＝20，T_{w_1}＝10s，T_{1_1}＝0.411s，T_{2_1}＝0.479s，T_{3_1}＝1.0s，T_{4_1}＝0.155s；K_{s_2}＝20，T_{w_2}＝10s，T_{1_2}＝0.01s，T_{2_2}＝0.54s，T_{3_2}＝0.707s，T_{4_2}＝0.081s。两个控制器的反馈和输出时滞分别为τ_f1＝150ms，τ_o1＝90ms，τ_f2＝70ms和τ_o2＝40ms。该系统的状态变量和代数变量分别为n＝200和l＝448。

选择不同的参数变量，在算例系统上分别对原始的IGD-IRK算法和低阶IGD-IRK算法进行验证，来证明本发明提出的低阶IGD-IRK算法的高效性。

为了计算出时滞电力系统的低频振荡模态，分别取N＝30，40和50，s＝2与N＝30，s＝3，计算出在位移点j10附近的r＝30个特征值，测量指标与计算时间如表1所示。

表1 IIGD算法和低阶IIGD算法的测度指标和计算时间

表中的N_IDR表示稀疏特征值计算过程中，利用迭代方法(IDR(s))求解矩阵逆与向量乘积的总迭代次数。由表1可以看出，低阶IGD-IRK算法与原始的IGD-IRK算法相比，维数至少减少了90％左右，说明低阶算法使得离散化矩阵的维数大幅降低。

接着，对比两种算法的计算量与计算时间，在相同的参数下，虽然两种算法的N_IDR差别不大，但是原始的IGD-IRK算法的计算时间是低阶IGD-IRK算法的4～6倍，由此可知，低阶IGD-IRK算法通过对离散化矩阵维数的降阶，减少了计算时间，大幅提高了算法的计算效率，使其能够高效适用于大规模时滞电力系统分析。

以上所述仅为本申请的优选实施例而已，并不用于限制本申请，对于本领域的技术人员来说，本申请可以有各种更改和变化。凡在本申请的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本申请的保护范围之内。

Claims (10)

1.基于低阶IGD-IRK的时滞电力系统小干扰稳定性分析方法，其特征是，包括：

建立时滞电力系统数学模型，线性化处理得到时滞电力系统的微分方程；

根据时滞电力系统状态变量是否与时滞有关，时滞电力系统的微分方程改写为与时滞有关的部分及与时滞无关的部分，继而基于与时滞的相关性重组状态变量，获得重组后的时滞电力系统状态方程；

获得与重组后的时滞电力系统状态方程相对应的无穷小生成元，依据无穷小生成元谱映射原理，将计算时滞电力系统机电振荡模式的问题转化为计算无穷小生成元的特征值问题；

针对无穷小生成元基于隐式龙格-库塔法低阶离散化，得到无穷小生成元低阶离散化近似矩阵；

针对无穷小生成元低阶离散化近似矩阵进行位移逆变化，得到无穷小生成元低阶离散化近似逆矩阵；

针对上述逆矩阵进行稀疏特征值计算，之后针对特征值再进行反变换和牛顿校验，得到时滞电力系统的精确特征值，特征值λ则对应时滞电力系统的机电振荡模式。

2.如权利要求1所述的基于低阶IGD-IRK的时滞电力系统小干扰稳定性分析方法，其特征是，经线性化后的时滞电力系统模型为：

式中：n为系统状态变量总数，Δx(t)为t时刻时系统状态变量的增量，为系统状态变量导数的增量，Δx(t-τ_i)为t-τ_i时刻时系统状态变量的增量，τ_i>0为时滞常数，i＝1,…,m，m为时滞的个数，且满足0<τ₁<τ₂<…<τ_i…<τ_max，其中τ_max为最大的时滞，Δx(0)为系统状态变量的初始值，简写为 为稠密的系统状态矩阵，是稀疏的系统时滞状态矩阵。

3.如权利要求1所述的基于低阶IGD-IRK的时滞电力系统小干扰稳定性分析方法，其特征是，根据时滞电力系统状态变量是否与时滞有关，将时滞电力系统的微分方程改写成两个部分，即分为与时滞无关部分和与时滞相关部分n₁为与时滞无关部分的状态变量的维数，n₂为与时滞相关部分的状态变量的维数，且满足n₁+n₂＝n，则描述系统动态特性的时滞微分方程重写为：

式(2)中：和分别为与时滞无关部分状态变量导数的增量与时滞相关部分状态变量导数的增量，Δx⁽¹⁾(t-τ_i)和Δx⁽²⁾(t-τ_i)分别为t-τ_i时刻时与时滞无关部分系统状态变量的增量和与时滞相关部分系统状态变量的增量，和分别是由状态矩阵和经过重写为与时滞有关和无关的部分得到的矩阵，i＝1,…,m：

式中：和是重写后矩阵中的分块矩阵；

与式(2)对应的系统的特征方程表示为：

式中：λ为时滞电力系统的特征值，v为特征值对应的右特征向量。

4.如权利要求3所述的基于低阶IGD-IRK的时滞电力系统小干扰稳定性分析方法，其特征是，利用无穷小生成元将式(2)在Banach空间映射为抽象柯西方程：

式中：θ∈[-τ_max,0]，为Δx_t的导数。

5.如权利要求1所述的基于低阶IGD-IRK的时滞电力系统小干扰稳定性分析方法，其特征是，时滞电力系统模型的特征值与无穷小生成元特征值之间的关系为：

式中：σ(·)表示无穷小生成元的谱，λ为时滞电力系统的特征值，该式说明时滞电力系统的特征值与无穷小生成元的特征值是一一对应的。

6.如权利要求4所述的基于低阶IGD-IRK的时滞电力系统小干扰稳定性分析方法，其特征是，基于隐式龙格-库塔离散化基本理论，在区间[-τ_max,0]上，建立Ns+1个离散点的集合Ω_N，N为给定正整数，s为隐式龙格-库塔法的级数，然后在每个离散点处，离散式(5)，估计Δx_t的值，生成无穷小生成元离散化矩阵；

式中：N_i为在第i个子区间以h_i为离散步长得到的离散点θ_j,i的个数，为在时滞区间[-τ_i,-τ_i-1]的离散点θ_j,i的基础上，采用s级IRK法的横坐标对每个子区间作进一步划分，得到的N_is个离散点构成的集合，i＝1,…,m，c_l为s级龙格-库塔法的横坐标；

只离散与时滞有关的部分，由此得与无穷小生成元对应的低阶离散化矩阵矩阵的阶数为(n+Nsn₂)，表示如下：

式中：s为隐式龙格-库塔法的级数，N为给定正整数，是n₂阶单位矩阵，是Nsn₂×n₁的零矩阵，表示克罗内克积运算，是矩阵的第一块行；子矩阵是高度稀疏的矩阵。

7.如权利要求6所述的基于低阶IGD-IRK的时滞电力系统小干扰稳定性分析方法，其特征是，为了充分利用克罗内克积的性质和状态矩阵的稀疏特性，从而进行稀疏特征值计算，∑_Ns表示为：

式中：为单位向量，为重组后的稠密的状态矩阵，为重组后时滞状态矩阵的后n₂列，分别可以表示如下：

。

8.如权利要求3所述的基于低阶IGD-IRK的时滞电力系统小干扰稳定性分析方法，其特征是，利用位移逆变换技术，将模值较小的这部分特征值转换为主特征值，优先计算出系统的关键模态，即，将λ'+s代入式(4)替代λ得位移之后的特征方程，表示如下：

将矩阵中的第一块行的用进行替代，i＝0,1,…,m，得到位移操作后的无穷小生成元的离散化矩阵 的逆矩阵表示为：

式(16)中：

式(17)中：是位移变换后得到的矩阵，是位移变换后时滞矩阵的后n₂列；

式中：为位移变换后的时滞矩阵的后n₂列和n₂行。

9.如权利要求8所述的基于低阶IGD-IRK的时滞电力系统小干扰稳定性分析方法，其特征是，采用隐式重启动Arnoldi算法来求取中的部分特征值，在Arnoldi算法中，利用与向量乘积形成Krylov子空间，设第k个Krylov向量为则第k+1个Krylov向量q_k+1计算如下：

采用诱导降维迭代方法来计算q_k+1，于是，式(19)写为：

式中：是第l次迭代后q_k+1的近似值。

10.如权利要求9所述的基于低阶IGD-IRK的时滞电力系统小干扰稳定性分析方法，其特征是，经过稀疏特征值计算后得到的的特征值为λ"，则的近似特征值为：

因此通过牛顿迭代法，得到时滞电力系统的精确特征值λ和与特征值对应的特征向量v。