CN108923421A - 基于低阶sod-irk算法的时滞电力系统关键特征值计算方法 - Google Patents

Abstract

本发明公开了基于低阶SOD‑IRK算法的时滞电力系统关键特征值计算方法，包括：将时滞电力系统模型状态变量分为时滞相关项和时滞无关项；将时滞电力系统模型的特征值问题转化成解算子的特征值问题；应用IRK方法对无穷维解算子进行离散化，得到有限维解算子离散化近似矩阵，通过计算解算子离散化近似矩阵的特征值，得到时滞电力系统降阶模型的关键特征值；计算解算子离散化近似矩阵的特征值时对解算子离散化近似矩阵进行旋转、放大预处理后采用稀疏特征值算法获得解算子离散化近似矩阵模值大于1的特征值；获得特定个数个模值最大的特征值；模值最大的特征值依次经过谱映射、反旋转、放大和牛顿校验，最终得到时滞电力系统模型的关键特征值。

Description

基于低阶SOD-IRK算法的时滞电力系统关键特征值计算方法

技术领域

本发明涉及电力系统技术领域，特别是涉及基于低阶SOD-IRK算法的时滞电力系统关键特征值计算方法，其中，SOD-IRK为“Solution Operator Discretization-ImplicitRunge Kutta”的英文缩写，对应的中文含义为：解算子隐式龙格库塔离散化。

背景技术

广域测量系统(Wide-Area Measurement System，WAMS)为互联电力系统提供实时广域测量信息。和本地信号相比，广域测量信息对于区间低频振荡模式显示出良好的观测性，比如联络线中的有功功率和相对转子角速度。将广域测量信号作为遥测反馈控制信号，广域阻尼控制器(Wide Area Damping Controllers，WADCs)对区间低频振荡有很好的抑制作用。然而，WAMS处理和传输广域信息的过程中会引入不可忽略的通讯时滞。这些时滞通常在几毫秒到几百毫秒之间，它们对系统小信号稳定性以及WADC设计的影响深远。因此，对于时滞信息物理融合的电力系统来说，得到有效的稳定性分析方法以及设计鲁棒控制器时必须考虑WAMS中的时滞问题。

近年来，越来越多学者开始关注将谱离散化方法用于时滞电力系统的特征分析。谱离散化方法可分为无穷小生成元离散化(Infinitesimal Generator Discretization，IGD)和解算子离散化(Solution Operator Discretization，SOD)两大类。中国发明专利基于显示IGD(Explicit IGD,EIGD)的大规模时滞电力系统特征值计算方法.201510055743.3.中国，提出了一种基于EIGD的大规模时滞电力系统特征值计算，这种计算方法属于第一类谱离散化方法，即无穷小生成元离散化方法。对于IGD类计算方法，通过一次计算只能得到给定点附近的几个特征值。为了获得所需的虚轴附近的关键振荡模态，往往需要不止一次计算。为了避免计算特征值时进行扫描，SOD类特征值计算方法可以解决这个问题。中国发明专利基于SOD-PS-R(Solution Operator Discretization-PesudoSpectrum and Rotation,解算子伪谱离散化和旋转)算法的时滞电力系统机电振荡模式计算方法.201510829129.8.[P].提出了一种基于SOD-PS-R的时滞电力系统关键特征值计算方法。

由于离散化方法的多样性，SOD类计算方法也可衍生出多种类型。本发明基于SOD-IRK算法，进一步运用降阶技术，改进了原来的算法，降低了运算时占用的计算机内存。

发明内容

为了解决现有技术的不足，本发明提供了基于低阶SOD-IRK算法的时滞电力系统关键特征值计算方法，用以得到大规模时滞电力系统的机电振荡模式。与原来的SOD-IRK算法相比，获得相同个数特征值时所需要的计算成本大大降低。

基于低阶SOD-IRK算法的时滞电力系统关键特征值计算方法，包括:

将时滞电力系统模型状态变量分为时滞相关项和时滞无关项，得到低阶时滞电力系统模型；

根据谱映射定理，低阶时滞电力系统模型的特征值与解算子特征值之间存在一一对应的关系，将时滞电力系统模型的特征值问题转化成解算子的特征值问题；

应用IRK方法对无穷维解算子进行离散化，得到有限维解算子离散化近似矩阵，通过计算解算子离散化近似矩阵的特征值，得到时滞电力系统降阶模型的关键特征值；

计算解算子离散化近似矩阵的特征值时对解算子离散化近似矩阵进行旋转、放大预处理后采用稀疏特征值算法获得解算子离散化近似矩阵模值大于1的特征值；

采用隐式重启动Arnoldi算法从解算子离散化近似矩阵的模值大于1的特征值中获得特定个数个模值最大的特征值；

模值最大的特征值依次经过谱映射、反旋转、放大和牛顿校验，最终得到时滞电力系统模型的关键特征值。

进一步优选的技术方案，降阶后的时滞电力系统模型为

式(4)中和分别为与时滞无关状态变量和与时滞相关状态变量的初始条件，并用和表示。其中，n₁为时滞无关状态变量的维数，n₂为时滞相关状态变量的维数。τ_i>0为时滞常数，i＝1,2,…,m。m为电力系统中时滞的个数。假设 为时滞常数，其中τ_max表示最大的时滞，为时滞无关项，为时滞相关项，为时滞无关项Δx₁(t)的导数，为时滞相关项Δx₁(t)的导数。其中，n₁+n₂＝n，n为时滞电力系统的维数，t为时间变量。分别为降阶后时滞电力系统模型的系统状态矩阵，为降阶后时滞电力系统模型的时滞状态矩阵。

进一步优选的技术方案，降阶后的时滞电力系统模型表示的线性化系统的特征方程为

式中：λ为特征值，和分别为系统降阶后特征值对应的与时滞无关和与时滞相关的右特征向量。

进一步优选的技术方案，所述解算子T(h):X→X是一种线性算子，将空间X上的初值状态转移到h时刻之后的系统状态；其中，h为转移步长，0≤h≤τ_max；

式(6)中包括常微分方程和转移方程两部分，其中，β为积分变量，ρ为变量，和分别为0和h+ρ时刻时滞电力系统的状态。

进一步优选的技术方案，所述降阶后的时滞电力系统模型的特征值μ与解算子特征值λ之间的关系如下：

式中：σ(T(h))表示算子的谱，\表示集合差运算。

进一步优选的技术方案，与解算子T(h)对应的离散化矩阵T_Ns表示如下：

式(8)中，N为IRK法的离散点数，s为IRK法的级数，和为零矩阵，为单位阵。

式(9)中：°表示Tracy–Singh积，A为IRK法的系数矩阵，为系统时滞状态矩阵；

式(10)中，为由拉格朗日插值系数决定的矩阵，γ为矩阵L_γ的序号。特别地，L₀和L_m分别表示如下：

式(11)中，为常系数矩阵，T表示对矩阵求转置。 式(10)中，E₂的表达式如下：

进一步优选的技术方案，当系统仅含有单个时滞即m＝1时，解算子的隐式龙格-库塔离散化过程中不再涉及拉格朗日插值，于是，Σ_Ns显式地表示如下：

式(16)中，为时滞电力系统的第m个时滞状态矩阵。

进一步优选的技术方案，所述旋转、放大预处理包括坐标轴旋转和特征值放大两个部分，坐标轴按照逆时针方向旋转θ角度后，则时滞电力系统的阻尼比小于给定常数ζ的关键特征值λ沿顺时针方向旋转θ弧度，对应着T_Ns的特征值由模值小于1变为模值大于1。

进一步优选的技术方案，在旋转-放大变换中，τ_i被变换原来的1/α，α表示放大倍数。从而，降阶前的系统矩阵被变换为h保持不变，

系统降阶以后，系统状态矩阵和时滞状态矩阵分别表示如下，i＝1,...,m：

式(18)中，

相应地，区间[-τ_max/α,0]重新划分为长度等于h的N'个子区间，N'＝τ_max/(αh)＝N/α，L_l(l＝0,1,…,m+1)被重新形成为最终，解算子的隐式龙格-库塔离散化矩阵T_Ns变为T_N′s：

式中：

进一步优选的技术方案，在Arnoldi算法中，设第j个Krylov向量表示为则第j+1个向量q_j+1可由矩阵T_N′s与向量q_j的乘积运算得到；

q_j+1＝T_N′sq_j＝(R′)^-1·Σ_N′s·q_j (24)

由于T_N′s具有的特殊逻辑结构，可知q_j+1的第sn+1:sn₁+sN'n₂个分量等于q_j的第sn₁+1:sn₁+s(N'-1)n₂个分量，即q_j+1(sn+1:sn₁+sN'n₂)＝q_j(sn₁+1:sn₁+s(N'-1)n₂,1)，q_j+1的第1:sn个分量，q_j+1(1:sn,1)，可以进一步分解为2个矩阵-向量乘积运算：

z＝Σ_N′s·q_j (25)

q_j+1(1:sn,1)＝(R′)^-1·z (26)

式中：为中间向量。

进一步优选的技术方案，式(25)的稀疏实现：

首先，从列的方向上将向量q_j(1:sn₁,1)、q_j(sn₁+1:sn,1)和q_j(sn+1:sn₁+sN'n₂,1)分别压缩为矩阵和即q_j(1:sn₁,1)＝vec(Q₁)、q_j(sn₁+1:sn₂,1)＝vec(Q₂)和q_j(sn+1:sn₁+sN′n₂,1)＝vec(Q₃)；

然后，将式(23)代入式(25)中，进而利用克罗内克积的性质，得：

进一步优选的技术方案，式(26)的稀疏实现：

采用迭代方法求解式(26)直至收敛；

式中：表示第k次迭代后q_j+1(1:sn,1)的近似值；

从列的方向上将向量和分别压缩为矩阵和即和则

进一步优选的技术方案，在计算得到解算子离散化矩阵的特征值μ″之后，依次经过谱映射、旋转-放大的反变换和牛顿校验后得到时滞电力系统的精确特征值λ，在牛顿校验前，λ由下列公式计算得到：

式中，μ″为解算子离散化矩阵的近似特征值，λ″为时滞电力系统的近似特征值。

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

第一、本发明提出的低阶SOD-IRK算法用于实际电力系统进行小干扰稳定性分析时，充分考虑了信号传输过程中时滞带来的影响以及系统的规模，只需一次计算特定个数的特征值，就可得到对系统进行稳定性分析所需的关键特征值。

第二、本发明的方法是基于低阶SOD-IRK算法大规模时滞电力系统的特征值计算方法，运用谱离散化方法设计框架，结合旋转-放大技术，通过计算有限维离散化矩阵的特征值，得到无限维时滞电力系统的特征值。

第三、本发明提出的低阶SOD-IRK算法在SOD-IRK算法的基础上，通过降阶技术大大减小了系统矩阵的维数，进而减小了计算量，提高了算法的计算效率。

附图说明

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

图1为基于低阶SOD-IRK的时滞电力系统关键特征值计算方法流程图。

具体实施方式

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

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

本申请的整体技术构思为：时滞电力系统模型由无时滞电力系统、广域阻尼控制器、反馈时滞以及控制时滞等四部分组成。将此模型线性化后，可得到考虑时滞影响的实际大规模电力系统小干扰稳定性分析的数学模型。把状态变量分为时滞相关项和时滞无关项，将此线性化模型降阶，从而得到维数较低的线性化方程。

降阶后解算子和时滞电力系统之间特征值集合的对应关系和降阶前是一致的。时滞电力系统频率在[0,2]Hz范围内的关键特征值λ被映射为解算子T(h)的特征值，其主要分布在单位圆的附近。

然后，用隐式龙格库塔(Implicit Runge Kutta,IRK)法把无穷维的解算子离散化为一个有限维的近似矩阵。计算这个离散化矩阵的模值最大特征值，然后由谱映射关系可得到系统的关键特征值。

为了充分利用解算子离散化矩阵的稀疏性，将离散化矩阵的主要子矩阵转换为常系数矩阵与系统状态矩阵的Kronecker积形式；为了优先得到系统的机电振荡模式并加快计算的收敛速度，采用旋转-放大预处理技术；充分利用解算子离散化矩阵和系统状态矩阵的稀疏性并应用Kronecker积的性质，实现了矩阵向量乘积的稀疏计算。

在此算法中，时滞电力系统离散化矩阵的阶数为s(n₁+Nn₂)，当面对大规模电力系统时，为了降低计算机占用内存，提高计算速度，需要采用稀疏特征值算法，如隐式重启动Arnoldi算法。在计算得到解算子的特征值之后，首先根据谱映射关系得到经过旋转和放大处理后时滞电力系统的特征值。然后，通过反旋转-放大技术，得到时滞电力系统机电振荡模式对应的关键特征值。最后，采用牛顿迭代法对近似特征值进行修正，得到时滞电力系统的精确特征值。

本申请的一种典型的实施方式中，如图1所示，基于SOD-IRK算法的时滞电力系统机电振荡模式计算方法，步骤如下：

S1：建立降阶后的时滞电力系统模型；把时滞电力系统状态变量分为时滞相关项和时滞无关项，对原来的时滞电力系统模型进行降阶，得到低维时滞电力系统模型；

S2：得到降阶后解算子和时滞电力系统之间特征值集合的对应关系，将时滞电力系统模型的特征值问题转化成解算子的特征值问题；应用IRK方法对无穷维解算子进行离散化，得到有限维解算子离散化近似矩阵；

S3：计算过程中采用旋转-放大技术预处理技术，并利用Kronecker积的性质进行稀疏实现。坐标轴旋转以后，阻尼比小于给定值的部分特征值λ，变换为解算子离散化近似矩阵模值大于1的特征值μ″；特征值放大以后，特征值μ″之间的相对距离增大，加快计算收敛；

S4：采用隐式重启动Arnoldi算法，从步骤S3的特征值μ″中获得解算子离散化近似矩阵的特定个数个特征值μ″，其模值最大；

解算子离散化近似矩阵的模值大于1的特征值对应时滞电力系统复平面右半平面上的特征值，所以优先获得模值最大的特征值可以加快判断系统的小干扰稳定性；

S5：从步骤S4中计算得到特征值μ″之后，依次将μ″经过谱映射、旋转-放大反变换和牛顿校验，最终得到时滞电力系统模型的关键特征值λ。根据谱映射关系，将离散化解算子的近似特征值转化为时滞电力系统的近似特征值；采用牛顿迭代法对近似特征值进行修正，得到时滞电力系统的精确特征值。

所述步骤S1中，降阶前时滞电力系统模型如下：

式中：为系统状态。τ_i>0(i＝1,2,…,m)为时滞常数。不失一般性地，假设 为时滞常数，其中τ_max表示最大的时滞。和分别为系统状态矩阵和时滞状态矩阵，前者为稠密矩阵，后者为稀疏矩阵。Δx(0)为系统状态变量的初始条件，用表示。

将状态变量Δx(t)分为与时滞无关项和与时滞相关项则(1)的第一式可改写为：

式中，n₁+n₂＝n，分别为降阶后时滞电力系统模型的系统状态矩阵，和为降阶后时滞电力系统模型的时滞状态矩阵，如下所示：

降阶后的时滞电力系统模型可写作

式(4)中和分别为与时滞无关状态变量和与时滞相关状态变量的初始条件，并用和表示。

式(4)表示的线性化系统的特征方程为

式中：λ为特征值，和分别为系统降阶后特征值对应的与时滞无关和与时滞相关的右特征向量。

步骤S2中，解算子T(h):X→X是一种线性算子，将空间X上的初值状态转移到h时刻之后的系统状态；其中，h为转移步长，0≤h≤τ_max；

式(6)中包括常微分方程和转移方程两部分。其中，β为积分变量，ρ为变量，和分别为0和h+θ时刻时滞电力系统的状态。

根据谱映射定理，所述时滞电力系统低阶模型的特征值μ与解算子特征值λ之间的关系如下：

式中：σ(T(h))表示算子的谱，\表示集合差运算。

复平面虚轴左侧对应着z平面中单位圆内部，为降阶时滞电力系统的稳定域。通过计算解算子模值最大的部分特征值，就可对时滞电力系统进行稳定性分析。若解算子特征值中存在模值比1大的特征值，即时滞电力系统的特征值实部大于零，那么该时滞电力系统处于不稳定状态；相反，若解算子的特征值模值全部小于1，即时滞电力系统的特征值实部全部小于零，那么该时滞电力系统为渐进稳定。若解算子特征值中模值最大的特征值恰好模值为1，即系统特征值实部为0，那么系统处于临界稳定状态。

解算子T(h)是一个无穷维度的线性算子。为了得到其特征值，首先采用IRK法对无穷维解算子T(h)进行离散化，得到一个有限维近似矩阵。通过计算解算子离散化近似矩阵的特征值，可得到时滞电力系统降阶模型的关键特征值。

与解算子T(h)对应的离散化矩阵T_Ns表示如下：

式(8)中，N为IRK法的离散点数，s为IRK法的级数，

式(9)中：°表示Tracy–Singh积，A为IRK法的系数矩阵。式(10)中，为由拉格朗日插值系数决定的矩阵，γ为矩阵L_γ的序号。特别地，L₀和L_m分别表示如下：

式(11)中，为常系数矩阵，T表示对矩阵求转置。 式(10)中，E₂的表达式如下：

特别地，当系统仅含有单个时滞(m＝1)时，解算子的隐式龙格-库塔离散化过程中不再涉及拉格朗日插值。于是，Σ_Ns可显式地表示如下：

式(16)中，为时滞电力系统的第m个时滞状态矩阵。

所述步骤S3中，利用旋转-放大预处理技术，包括坐标轴旋转和特征值放大两个部分。首先，坐标轴按照逆时针方向旋转θ角度后，则时滞电力系统的阻尼比小于给定常数ζ(ζ＝sinθ)的关键特征值λ沿顺时针方向旋转θ弧度，对应着T_Ns的特征值由模值小于1变为模值大于1。

在旋转-放大变换中，τ_i(i＝1,…,m)被变换原来的1/α，α表示放大倍数。从而，降阶前的系统矩阵被变换为h保持不变，

系统降阶以后，系统状态矩阵和时滞状态矩阵分别表示如下：

式(18)中，

相应地，区间[-τ_max/α,0]重新划分为长度等于h的N'个子区间，N'＝τ_max/(αh)＝N/α，L_l(l＝0,1,…,m+1)被重新形成为最终，解算子的隐式龙格-库塔离散化矩阵T_Ns变为T_N′s。

式中：

所述步骤S4中，矩阵T_N′s的维数为s(n₁+N′n₂)。面对大规模时滞电力系统时，矩阵阶数巨大。因此，为了得到所需要的关键特征值，需要采用稀疏特征值算法，如序贯法或子空间法。

所述步骤S4的步骤如下：

在Arnoldi算法中，最关键的操作就是在形成Krylov向量过程中的矩阵向量乘积。设第j个Krylov向量表示为则第j+1个向量q_j+1可由矩阵T_N′s与向量q_j的乘积运算得到。

q_j+1＝T_N′sq_j＝(R′)^-1·Σ_N′s·q_j (24)

由于T_N′s具有的特殊逻辑结构，可知q_j+1的第sn+1:sn₁+sN'n₂个分量等于q_j的第sn₁+1:sn₁+s(N'-1)n₂个分量，即q_j+1(sn+1:sn₁+sN'n₂)＝q_j(sn₁+1:sn₁+s(N'-1)n₂,1)。q_j+1的第1:sn个分量，q_j+1(1:sn,1)，可以进一步分解为2个矩阵-向量乘积运算：

z＝Σ_N′s·q_j (25)

q_j+1(1:sn,1)＝(R′)^-1·z (26)

式中：为中间向量。

接下来，本节将重点分析式(25)和式(26)的稀疏实现方法。

式(25)的稀疏实现：

首先，从列的方向上将向量q_j(1:sn₁,1)、q_j(sn₁+1:sn,1)和q_j(sn+1:sn₁+sN'n₂,1)分别压缩为矩阵和即q_j(1:sn₁,1)＝vec(Q₁)、q_j(sn₁+1:sn₂,1)＝vec(Q₂)和q_j(sn+1:sn₁+sN′n₂,1)＝vec(Q₃)。然后，将式(23)代入式(25)中，进而利用克罗内克积的性质，得：

式(26)的稀疏实现：

采用迭代方法(如IDR(s)方法)求解式(26)直至收敛。

式中：表示第k次迭代后q_j+1(1:sn,1)的近似值。

从列的方向上将向量和分别压缩为矩阵和即和则

所述步骤S5中，在计算得到解算子离散化矩阵的特征值μ″之后，依次经过谱映射、反旋转-放大和牛顿校验后得到时滞电力系统的精确特征值λ，在牛顿校验前，λ由下列公式计算得到：

式中，μ″为解算子离散化矩阵的近似特征值，λ″为时滞电力系统的近似特征值。

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

Claims (10)

1.基于低阶SOD-IRK算法的时滞电力系统关键特征值计算方法，其特征是，包括:

将时滞电力系统模型状态变量分为时滞相关项和时滞无关项，得到低阶时滞电力系统模型；

根据谱映射定理，低阶时滞电力系统模型的特征值与解算子特征值之间存在一一对应的关系，将时滞电力系统模型的特征值问题转化成解算子的特征值问题；

应用IRK方法对无穷维解算子进行离散化，得到有限维解算子离散化近似矩阵，通过计算解算子离散化近似矩阵的特征值，得到时滞电力系统降阶模型的关键特征值；

计算解算子离散化近似矩阵的特征值时对解算子离散化近似矩阵进行旋转、放大预处理后采用稀疏特征值算法获得解算子离散化近似矩阵模值大于1的特征值；

采用隐式重启动Arnoldi算法从解算子离散化近似矩阵的模值大于1的特征值中获得特定个数个模值最大的特征值；

模值最大的特征值依次经过谱映射、反旋转、放大和牛顿校验，最终得到时滞电力系统模型的关键特征值。

2.如权利要求1所述的基于低阶SOD-IRK算法的时滞电力系统关键特征值计算方法，其特征是，降阶后的时滞电力系统模型为

式(4)中和分别为与时滞无关状态变量和与时滞相关状态变量的初始条件，并用和表示，其中，n₁为时滞无关状态变量的维数，n₂为时滞相关状态变量的维数，τ_i>0为时滞常数，i＝1,2,…,m，m为电力系统中时滞的个数，假设 为时滞常数，其中τ_max表示最大的时滞，为时滞无关项，为时滞相关项，为时滞无关项Δx₁(t)的导数，为时滞相关项Δx₁(t)的导数，其中，n₁+n₂＝n，n为时滞电力系统的维数，t为时间变量，分别为降阶后时滞电力系统模型的系统状态矩阵，为降阶后时滞电力系统模型的时滞状态矩阵；

降阶后的时滞电力系统模型表示的线性化系统的特征方程为

式中：λ为特征值，和分别为系统降阶后特征值对应的与时滞无关和与时滞相关的右特征向量。

3.如权利要求1所述的基于低阶SOD-IRK算法的时滞电力系统关键特征值计算方法，其特征是，所述解算子T(h):X→X是一种线性算子，将空间X上的初值状态转移到h时刻之后的系统状态；其中，h为转移步长，0≤h≤τ_max；

式(6)中包括常微分方程和转移方程两部分，其中，β为积分变量，ρ为变量，和分别为0和h+θ时刻时滞电力系统的状态；

所述降阶后的时滞电力系统模型的特征值μ与解算子特征值λ之间的关系如下：

式中：σ(T(h))表示算子的谱，\表示集合差运算。

4.如权利要求3所述的基于低阶SOD-IRK算法的时滞电力系统关键特征值计算方法，其特征是，与解算子T(h)对应的离散化矩阵T_Ns表示如下：

式(8)中，N为IRK法的离散点数，s为IRK法的级数，和为零矩阵，为单位阵；

式(9)中：°表示Tracy–Singh积，A为IRK法的系数矩阵，为系统时滞状态矩阵；

式(10)中，为由拉格朗日插值系数决定的矩阵，γ为矩阵L_γ的序号，特别地，L₀和L_m分别表示如下：

式(11)中，为常系数矩阵，T表示对矩阵求转置， 式(10)中，E₂的表达式如下：

当系统仅含有单个时滞即m＝1时，解算子的隐式龙格-库塔离散化过程中不再涉及拉格朗日插值，于是，Σ_Ns显式地表示如下：

式(16)中，为时滞电力系统的第m个时滞状态矩阵。

5.如权利要求1所述的基于低阶SOD-IRK算法的时滞电力系统关键特征值计算方法，其特征是，旋转、放大预处理包括坐标轴旋转和特征值放大两个部分，坐标轴按照逆时针方向旋转θ角度后，则时滞电力系统的阻尼比小于给定常数ζ的关键特征值λ沿顺时针方向旋转θ弧度，对应着T_Ns的特征值由模值小于1变为模值大于1。

6.如权利要求5所述的基于低阶SOD-IRK算法的时滞电力系统关键特征值计算方法，其特征是，在旋转、放大变换中，τ_i被变换原来的1/α，α表示放大倍数，从而，降阶前的系统矩阵被变换为h保持不变，

系统降阶以后，系统状态矩阵和时滞状态矩阵分别表示如下，i＝1,...,m：

式(18)中，

相应地，区间[-τ_max/α,0]重新划分为长度等于h的N'个子区间，N'＝τ_max/(αh)＝N/α，L_i(l＝0,1,…,m+1)被重新形成为最终，解算子的隐式龙格-库塔离散化矩阵T_Ns变为T_N′s：

式中：

7.如权利要求1所述的基于低阶SOD-IRK算法的时滞电力系统关键特征值计算方法，其特征是，在Arnoldi算法中，设第j个Krylov向量表示为则第j+1个向量q_j+1可由矩阵T_N′s与向量q_j的乘积运算得到；

q_j+1＝T_N′sq_j＝(R′)^-1·Σ_N′s·q_j (24)

由于T_N′s具有的特殊逻辑结构，可知q_j+1的第sn+1:sn₁+sN'n₂个分量等于q_j的第sn₁+1:sn₁+s(N'-1)n₂个分量，即q_j+1(sn+1:sn₁+sN'n₂)＝q_j(sn₁+1:sn₁+s(N'-1)n₂,1)，q_j+1的第1:sn个分量，q_j+1(1:sn,1)，可以进一步分解为2个矩阵-向量乘积运算：

z＝Σ_N′s·q_j (25)

q_j+1(1:sn,1)＝(R′)^-1·z (26)

式中：为中间向量。

8.如权利要求7所述的基于低阶SOD-IRK算法的时滞电力系统关键特征值计算方法，其特征是，式(25)的稀疏实现：

首先，从列的方向上将向量q_j(1:sn₁,1)、q_j(sn₁+1:sn,1)和q_j(sn+1:sn₁+sN'n₂,1)分别压缩为矩阵和即q_j(1:sn₁,1)＝vec(Q₁)、q_j(sn₁+1:sn₂,1)＝vec(Q₂)和q_j(sn+1:sn₁+sN′n₂,1)＝vec(Q₃)；

然后，将式(23)代入式(25)中，进而利用克罗内克积的性质，得：

9.如权利要求7所述的基于低阶SOD-IRK算法的时滞电力系统关键特征值计算方法，其特征是，式(26)的稀疏实现：

采用迭代方法求解式(26)直至收敛；

式中：表示第k次迭代后q_j+1(1:sn,1)的近似值；

从列的方向上将向量和分别压缩为矩阵和即和则

10.如权利要求1所述的基于低阶SOD-IRK算法的时滞电力系统关键特征值计算方法，其特征是，在计算得到解算子离散化矩阵的特征值μ″之后，依次经过谱映射、反旋转-放大和牛顿校验后得到时滞电力系统的精确特征值λ，在牛顿校验前，λ由下列公式计算得到：

式中，μ″为解算子离散化矩阵的近似特征值，λ″为时滞电力系统的近似特征值。