CN109033022B - 基于低阶sod-lms算法的时滞电力系统特征值计算方法 - Google Patents

Abstract

本发明公开了基于低阶SOD‑LMS算法的时滞电力系统特征值计算方法，包括：针对时滞电力系统模型进行降阶处理得到低阶时滞电力系统模型；通过线性多步法方法对解算子进行离散化，得到解算子的离散化矩阵；利用坐标旋转预处理方法，将解算子的离散化矩阵旋转并进行近似得到旋转后的解算子离散化近似矩阵，将低阶时滞电力系统模型的阻尼比小于给定阻尼比ζ的关键特征值λ，变换为对应的旋转后的解算子离散化近似矩阵模值大于1的特征值μ″；从旋转后的解算子离散化近似矩阵中计算得到解算子模值最大的部分的近似特征值μ″；得到特征值μ″之后，对μ″进行反变换和牛顿校验，最终得到时滞电力系统模型的特征值λ。

Description

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

技术领域

本发明涉及电力系统技术领域，特别是涉及基于低阶SOD-LMS算法的时滞电力系统特征值计算方法，SOD-LMS为“Solution Operator with Linear Multi-Step”的英文缩写，中文含义：线性多步离散化。

背景技术

互联同步发电机之间的机电振荡模态的稳定性是系统安全运行所必需的。电力系统中出现的机电振荡按振荡范围及振荡频率的大小可大致分为两种类型：局部振荡模态和区间振荡模态。传统的电力系统稳定器(Power System Stabilizer,PSS)使用本地量测信息形成反馈控制，使系统可有效处理本地低频振荡。然而，由于传统PSS仅对本地信息进行量测，因此其不能很好地抑制区间低频振荡问题。

随着广域量测系统(Wide-Area Measurement System,WAMS)的发展，采用广域测量信号反馈的广域阻尼控制器(Wide-Area Damping Controller,WADC)可以有效提高区间低频振荡模式的阻尼水平。然而，相对于采用本地量测信号的阻尼控制，在广域阻尼控制中，时滞明显增加。在实际应用中，系统时滞从几十到几百毫秒不等。在系统存在时滞的情况下，不考虑时滞的WADC性能会恶化，甚至不能正常工作。

在电力系统中，机电振荡问题通常利用小干扰稳定性进行分析。时滞电力系统的小干扰稳定分析方法大体上可分为时域法和频域法两大类。时域法主要基于Lyapunov-Krasovskii稳定性定理和Razumikhin定理构建时滞依赖稳定性判据(Lyapunov泛函)，进而判定系统的时滞稳定性。而频域法，旨在通过计算系统线性化状态矩阵的特征值，分析系统在运行点附近的渐进稳定性。根据对时滞环节处理方式的不同，频域法可进一步分为预测法、函数变换法和特征分析法。常用的函数变换法包括Rekasius变换(或双线性变换、特征根聚类)、Lambert-W函数、Padé近似等。预测补偿法对受控对象的动态特性进行估计，用预估模型进行补偿，从而消除时滞对系统的影响。预测法包括Smith预估和模型预测法。特征分析法通过求解系统的特征值，进而沿用经典的特征值分析方法来分析系统的小干扰稳定性，进而可以对WADC进行优化。特征分析法既不需要对时滞在特征方程中引入的指数项进行特殊处理，也不需要对系统模型进行降阶处理。因此，相较于其他几种方法，特征分析法更为直接、准确，是时滞电力系统稳定性分析的理想工具。

中国发明专利基于Padé近似的时滞电力系统特征值计算与稳定性判别方法.201210271783.8:[P].利用Pade近似多项式逼近时滞环节，进而计算系统最右侧的特征值，进而判断时滞系统的稳定性。中国发明专利基于EIGD的大规模时滞电力系统特征值计算方法.201510055743.3:[P].提出了一种基于显式IGD(Explicit IGD,EIGD)方法的大规模时滞电力系统特征值计算方法。利用计算得到的系统最右侧的特征值，可以判断系统在固定时滞下的稳定性。

存在的问题是：已有的技术均未对矩阵进行降阶处理，在算法中对高阶矩阵进行处理所需的计算量显著高于对低阶矩阵进行相同的处理。

发明内容

为了解决现有技术的不足，本发明提供了基于低阶SOD-LMS算法的时滞电力系统特征值计算方法，用以得到大规模时滞电力系统的机电振荡模式。SOD-LMS算法只需计算解算子离散化近似矩阵中设定个数个模值最大的特征值，通过一次计算，就能够得到大规模时滞电力系统的机电振荡模式。

基于低阶SOD-LMS算法的时滞电力系统特征值计算方法，包括：

针对时滞电力系统模型进行降阶处理得到低阶时滞电力系统模型，并根据低阶时滞电力系统模型，将其转化为对应的解算子，低阶时滞电力系统模型的特征值与低阶时滞电力系统模型的解算子特征值为一一对应关系，从而将计算时滞电力系统特征值的问题转化为计算对应解算子特征值的问题；

通过线性多步法方法对解算子进行离散化，得到解算子的离散化矩阵；

利用坐标旋转预处理方法，将解算子的离散化矩阵旋转并进行近似得到旋转后的解算子离散化近似矩阵，将低阶时滞电力系统模型的阻尼比小于给定阻尼比ζ的关键特征值λ，变换为对应的旋转后的解算子离散化近似矩阵模值大于1的特征值μ″；

得到特征值μ″之后，对μ″进行反变换和牛顿校验，最终得到时滞电力系统模型的特征值λ。

进一步优选的技术方案，将时滞电力系统模型的时滞微分代数方程变换为时滞微分方程，根据时滞微分方程将状态变量分为与时滞无关项及与时滞相关项，从而得到低阶时滞电力系统模型。

进一步优选的技术方案，所述低阶时滞电力系统模型如下：

式中：

和

分别为与时滞无关和与时滞相关的电力系统状态变量向量，n₁和n₂分别为与时滞相关和与时滞无关的系统状态变量个数，t为当前时刻，τ_i>0(i＝1,2,…,m)为时滞常数，并假设0<τ₁<τ₂<…<τ_i…<τ_m＝τ_max为时滞环节的时滞常数，其中m为系统中时滞的个数，τ_max为最大的时滞，

和

为系统状态矩阵

的不同分块，均为稠密矩阵，

(i＝1,…,m)为系统时滞状态矩阵，为稀疏矩阵，Δx₁(t)和Δx₂(t)分别为t时刻系统状态变量x₁和x₂的增量，Δx₂(t-τ_i)为t-τ_i时刻x₂的增量，

和

分别为t时刻系统状态变量x₁和x₂导数的增量，Δx_0,1和Δx_0,2分别为系统状态变量x₁和x₂的初始值，可分别记为

和

进一步优选的技术方案，所述低阶时滞电力系统模型对应的线性化系统的特征方程为：

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

为系统状态矩阵，

为与不同时滞相对应的系统时滞状态矩阵。

进一步优选的技术方案，根据低阶时滞电力系统模型，将其转化为对应的解算子，具体为：

设

为由区间[-τ_max,0]到n维实数空间

映射的连续函数构成的巴拿赫空间，解算子

定义为将空间X上的初值条件

转移到h时刻之后时滞电力系统解分段的线性算子；其中，h为转移步长，0≤h≤τ_m；

其中，s为积分变量，θ为变量，

和

分别为0和h+θ时刻时滞电力系统的状态。

进一步优选的技术方案，解算子

的特征值μ与时滞系统的特征值λ之间具有如下关系：

式中：

表示解算子的谱，\表示集合差运算。

进一步优选的技术方案，与解算子

对应的离散化矩阵T_N表示如下：

式中：

和

分别为对应阶数的零矩阵，I_(k-1)n和

分别为对应阶数的对角单位矩阵。

式(10)和(11)中，N为离散化维数，N＝Q+k+s-，其中k为LMS法的系数，s-为Nordsieck插值的系数，而

h为k步LMS法的步长，α_k和β_k为k步LMS法的系数，

其中，

中的元素完全由LMS法的系数α_k决定，

中的元素由LMS法的系数β_k和步长h共同决定，而

中的元素由拉格朗日插值系数和LMS法的系数β_k和步长h共同决定，

和

分别为对应阶数的零矩阵，

和I_n分别为对应阶数的对角单位矩阵，

为Kronecker积运算。

进一步优选的技术方案，利用坐标旋转预处理方法，将解算子的离散化矩阵旋转并进行近似得到旋转后的解算子离散化近似矩阵，具体为：

将坐标轴逆时针旋转θ角度，将低阶时滞电力系统模型对应的线性化系统的特征方程中的λ用λ′e^-jθ代替，可以得到坐标轴旋转后的特征方程：

式中：

λ′＝λe^-jθ (13)

τ′_i＝τ_ie^jθ,i＝1,...,m (14)

对坐标旋转后的τ′_i(i＝1,…,m)进行近似，因此，式(14)变为：

τ′_i＝τ_ie^jθ≈τ_i,i＝1,...,m (16)

式(12)变为：

式中：

为和

分别为λ′和与其对应的右特征向量v的近似值。

进一步优选的技术方案，坐标轴旋转预处理后，λ'与其

的特征值μ'之间的映射关系为：

进一步优选的技术方案，为了改善算法的收敛性，对特征值μ'进行非线性放大，从而增大不同特征值之间的相对距离，假设对μ'进行α次乘方，则由式(18)得：

式中：

由式(19)可知，保持

和τ_i(i＝1,…,m)不变，将h增大α倍即可实现旋转-放大预处理，这表明，对μ'的放大处理，等价于将转移步长h增大α倍；

由于转移步长h被变换原来的α倍，区间[-τ_max,0]被重新划分为长度等于或小于αh的Q'个子区间，

N变为N'＝Q'+k+s-，l_i(i＝1,…,m)被重新形成为

T_N变为T_N′：

式中：

T_N′的谱与μ"之间有如下对应关系：

μ″∈σ(T′_N)\{0} (26)

进一步优选的技术方案，采用迭代特征值算法计算T_N′的设定个数个模值最大的特征值在应用解算子对应的离散化矩阵T_N′求解大规模电力系统的时滞特征值时步骤为：

设第j个Krylov向量表示为

则第j+1个向量q_j+1可通过矩阵T_N′与向量q_j的乘积计算得到：

由于T_N′具有特殊的逻辑结构，q_j+1的第n+1:(kn₁+N'n₂)个分量与q_j的第1:(k-1)n个和第kn-n₂+1:(kn₁+(N'-1)n₂)个分量之间存在一一对应关系，即：

而q_j+1的第1:n个分量，即q_j+1(1:n,1)，可以进一步分解为2个矩阵-向量乘积运算：

z＝Σ_N′·q_j (29)

式中：

为中间向量。

进一步优选的技术方案，利用稀疏实现方法计算式(29)，首先需要从列的方向上将向量q_k压缩为矩阵

和

即q_k＝[vec(Q₁)^T,vec(Q₂)^T]^T；

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

式(31)中，z的计算量主要由稠密矩阵-向量乘积

决定，为了减少复数运算，将式(15)代入，则其可以改写为

其中

可将其计算步骤分解为如下两个步骤：

D₀w＝-C₀q_k+1 (32)

z₀＝A′₀q_k+1+B′₀w (33)

将其改写为矩阵形式，可得：

进一步优选的技术方案，利用稀疏实现方法计算式(30)，首先需要将式(15)和式(1)代入式(24)中，得：

式中：A′₀,B′₀与A₀,B₀具有完全相同的稀疏特性，

A′₀＝α_kI_n-hβ_kαA₀e^-jθ (36)

B′₀＝-hβ_kαB₀e^-jθ (37)

然后利用矩阵之和的求逆公式计算

可得：

于是，式(30)中

的稀疏实现可以分解为如下两个步骤：

[D₀-C₀(A′₀)^-1B′₀]w＝-C₀(A′₀)^-1z (39)

q_k+1＝(A′₀)^-1(z-B′₀w) (40)

将其改写为矩阵形式，可得：

进一步优选的技术方案，在计算得到μ″之后，依次经过谱映射、坐标反旋转和牛顿校验后得到时滞电力系统的特征值λ，在牛顿校验前，估计特征值

的计算公式为：

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

第一、本发明提出的低阶SOD-LMS算法给出了低阶情况下解算子离散化矩阵及其变换形式。

第二、本发明提出的低阶SOD-LMS算法用于计算实际系统机电振荡模式对应的关键特征值时，能够大幅降低解算子离散化矩阵的维数，提高计算效率。第三、本发明提出的低阶SOD-LMS算法给出了低阶情况下对矩阵求解过程进行稀疏实现的具体形式。

附图说明

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

图1为基于低阶SOD-LMS电力系统机电振荡模式计算方法的流程图。

具体实施方式

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

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

本申请的一种典型的实施方式中，如图1所示，提供了基于低阶SOD-LMS算法的时滞电力系统特征值计算方法，包括如下步骤：

S1：建立低阶时滞系统模型；利用电力系统中受时滞影响的状态变量维数远小于总状态变量维数这一性质，将原有的时滞电力系统模型转化为低阶时滞电力系统模型。并根据低阶时滞电力系统模型，将其转化为对应的解算子，从而将计算电力系统特征值的问题转化为计算对应解算子特征值的问题；

S2：通过线性多步法方法对解算子进行离散化，得到解算子的离散化矩阵；

S3：利用坐标旋转预处理方法，将解算子的离散化矩阵旋转并进行近似得到旋转后的解算子离散化近似矩阵，将时滞电力系统模型的阻尼比小于给定阻尼比ζ的关键特征值λ，变换为对应的旋转后的解算子离散化近似矩阵模值大于1的特征值μ″；

S4：采用隐式重启动Arnoldi算法，从步骤S3中旋转后的解算子离散化近似矩阵中计算得到解算子模值最大的部分的近似特征值μ″；

S5：从步骤S4中得到特征值μ″之后，对μ″进行反变换和牛顿校验，最终得到时滞电力系统模型的特征值λ，即为时滞电力系统的机电振荡模式。

所述步骤S1中，首先给出DCPPS的时滞微分代数方程：

式中：x、y分别为系统状态变量和代数变量，τ_i>0(i＝1,2,…,m)为时滞常数，并假设

其中τ_max表示最大的时滞。

将式(1)变换为时滞微分方程，并考虑初值，有：

然后，将状态变量分为与时滞无关项

和与时滞相关项

可将式(2)的第一式可改写为：

式中，n₁+n₂＝n，

所述步骤S1中的低阶时滞电力系统模型如下：

式中：

和

分别为与时滞相关和与时滞无关的电力系统状态变量向量，n₁和n₂分别为与时滞相关和与时滞无关的系统状态变量个数。t为当前时刻。0<τ₁<τ₂<…<τ_i…<τ_m＝τ_max为时滞环节的时滞常数。

和

为系统状态矩阵

的不同分块，均为稠密矩阵。

(i＝1,…,m)为系统时滞状态矩阵，为稀疏矩阵。Δx₁(t)和Δx₂(t)分别为t时刻系统状态变量x₁和x₂的增量，Δx₂(t-τ_i)为t-τ_i时刻x₂的增量，

和

分别为t时刻系统状态变量x₁和x₂导数的增量。Δx_0,1和Δx_0,2分别为系统状态变量x₁和x₂的初始值，可分别记为

和

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

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

所述步骤S2中，设

为由区间[-τ_max,0]到n维实数空间

映射的连续函数构成的巴拿赫空间。解算子

定义为将空间X上的初值条件

转移到h时刻之后时滞电力系统解分段的线性算子；其中，h为转移步长，0≤h≤τ_m；

其中，s为积分变量，θ为变量，

和

分别为0和h+θ时刻时滞电力系统的状态。

所述时滞电力系统模型的特征值与时滞电力系统模型的解算子特征值之间的关系：

由谱映射定理可知，解算子

的特征值μ与时滞系统的特征值λ之间具有如下关系：

式中：

表示解算子的谱，\表示集合差运算。

通过式(8)可知，系统正实部的特征值则被映射到单位圆之外，而系统的负实部特征值被映射到解算子单位圆之内。如果解算子至少存在一个单位圆之外的特征值，则可以原系统是不稳定的，如果解算子所有特征值的模值均位于单位圆之内，则原系统是稳定的。

解算子

是表征X→X映射的无穷维线性算子。为了计算其特征值，首先采用线性多步法(Linear MultiStep,LMS)对解算子

进行离散化，得到与

对应的有限维的近似矩阵T_N，进而计算T_N的特征值并得到原系统机电振荡模式对应的特征值。

所述步骤S2的步骤如下：

与解算子

对应的离散化矩阵T_N表示如下：

式中：

和

分别为对应阶数的零矩阵，I_(k-1)n和

分别为对应阶数的对角单位矩阵。

式(10)和(11)中，N为离散化维数，N＝Q+k+s_-，其中k为LMS法的系数，s_-为Nordsieck插值的系数，而

h为k步LMS法的步长。α_k和β_k为k步LMS法的系数。

其中，

中的元素完全由LMS法的系数α决定，

中的元素由LMS法的系数β和步长h共同决定，而

中的元素由拉格朗日插值系数和LMS法的系数β和步长h共同决定，

和

分别为对应阶数的零矩阵，

和I_n分别为对应阶数的对角单位矩阵，

为Kronecker积运算。

所述步骤S3中，首先将坐标轴逆时针旋转θ角度，将电力系统中阻尼比小于最大计算阻尼比ζ(ζ＝sinθ)的关键特征值λ，变换为T_N模值大于1的特征值μ。

将式(6)中的λ用λ′e^-jθ代替，可以得到坐标轴旋转后的特征方程：

式中：

λ′＝λe^-jθ (13)

τ′_i＝τ_ie^jθ,i＝1,...,m (14)

由式(14)可知，当θ≠0时，坐标旋转后τ′_i(i＝1,…,m)变为复数。然而，离散化过程中需要将区间[-τ_max,0]划分为N＝Q+k+s_-个长度为h的子区间，其中

由于h是正实数、Q为正整数，因此τ′_max＝τ′_m必须为实数。显然，实际上为复数的τ′_max不能满足这一要求。因此，需要对坐标旋转后的τ′_i(i＝1,…,m)进行近似。因此，式(14)变为：

τ′_i＝τ_ie^jθ≈τ_i,i＝1,...,m (16)

式(12)变为：

式中：

为和

分别为λ′和与其对应的右特征向量v的近似值。

坐标轴旋转预处理后，λ'与其

的特征值μ'之间的映射关系为：

然后，为了改善算法的收敛性，可对特征值μ'进行非线性放大，从而增大不同特征值之间的相对距离。假设对μ'进行α次乘方，则由式(18)得：

式中：

由式(19)可知，保持

和τ_i(i＝1,…,m)不变，将h增大α倍即可实现旋转-放大预处理。这表明，对μ'的放大处理，可以等价于将转移步长h增大α倍。

由于转移步长h被变换原来的α倍，区间[-τ_max,0]被重新划分为长度等于(或小于)αh的Q'个子区间，

N变为N'＝Q'+k+s_-。l_i(i＝0,1,…,m+1)被重新形成为

最终，T_N变为T_N′：

式中：

T_N′的谱与μ"之间有如下对应关系：

μ″∈σ(T′_N)\{0} (26)

所述步骤S4中，矩阵T_N′的阶数为N'n。对于大规模电力系统，这个阶数将非常大。因此，必须采用迭代特征值算法计算T_N′的设定个数个模值最大的特征值在应用解算子对应的离散化矩阵T_N′求解大规模电力系统的时滞特征值时，所述步骤S4的步骤如下：

在迭代特征值算法中，最关键的操作就是在形成Krylov向量过程中的矩阵向量乘积。设第j个Krylov向量表示为

则第j+1个向量q_j+1可通过矩阵T_N′与向量q_j的乘积计算得到：

由于T_N′具有特殊的逻辑结构，q_j+1的第n+1:(kn₁+N'n₂)个分量与q_j的第1:(k-1)n个和第kn-n₂+1:(kn₁+(N'-1)n₂)个分量之间存在一一对应关系，即：

而q_j+1的第1:n个分量，即q_j+1(1:n,1)，可以进一步分解为2个矩阵-向量乘积运算：

z＝Σ_N′·q_j (29)

式中：

为中间向量。

由于直接计算式(29)和式(30)的计算量较大，可以利用稀疏实现方法降低计算量。下面分别介绍利用稀疏实现方法计算两者的步骤。

利用稀疏实现方法计算式(29)，首先需要从列的方向上将向量q_k压缩为矩阵

和

即q_k＝[vec(Q₁)^T,vec(Q₂)^T]^T。然后，将式(24)代入式(29)中，进而利用克罗内克积的性质，得：

式(31)中，z的计算量主要由稠密矩阵-向量乘积

决定。为了减少复数运算，将式(15)代入，则其可以改写为

其中

可将其计算步骤分解为如下两个步骤：

D₀w＝-C₀q_k+1 (32)

z₀＝A′₀q_k+1+B′₀w (33)

将其改写为矩阵形式，可得：

利用稀疏实现方法计算式(30)，首先需要将式(15)和式(1)代入式(24)中，可得：

式中：A′₀,B′₀与A₀,B₀具有完全相同的稀疏特性，

A′₀＝α_kI_n-hβ_kαA₀e^-jθ (36)

B′₀＝-hβ_kαB₀e^-jθ (37)

然后利用矩阵之和的求逆公式计算

可得：

于是，式(30)中

的稀疏实现可以分解为如下两个步骤：

[D₀-C₀(A′₀)^-1B′₀]w＝-C₀(A′₀)^-1z (39)

q_k+1＝(A′₀)^-1(z-B′₀w) (40)

将其改写为矩阵形式，可得：

所述步骤S5中，在计算得到μ″之后，依次经过谱映射、坐标反旋转和牛顿校验后得到时滞电力系统的特征值λ，在牛顿校验前，估计特征值

的计算公式为：

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

Claims (9)

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

针对时滞电力系统模型进行降阶处理得到低阶时滞电力系统模型，并根据低阶时滞电力系统模型，将其转化为对应的解算子，低阶时滞电力系统模型的特征值与低阶时滞电力系统模型的解算子特征值为一一对应关系，从而将计算时滞电力系统特征值的问题转化为计算对应解算子特征值的问题；

通过线性多步法对解算子进行离散化，得到解算子的离散化矩阵；

利用坐标旋转预处理方法，将解算子的离散化矩阵旋转并进行近似得到旋转后的解算子离散化近似矩阵，将低阶时滞电力系统模型的阻尼比小于给定阻尼比ζ的关键特征值λ，变换为对应的旋转后的解算子离散化近似矩阵模值大于1的特征值μ″；

得到特征值μ″之后，对μ″进行反变换和牛顿校验，最终得到时滞电力系统模型的特征值λ；

将时滞电力系统模型的时滞微分代数方程变换为时滞微分方程，根据时滞微分方程将状态变量分为与时滞无关项及与时滞相关项，从而得到低阶时滞电力系统模型；

所述低阶时滞电力系统模型如下：

式中：

和

分别为与时滞无关和与时滞相关的电力系统状态变量向量，n₁和n₂分别为与时滞无关和与时滞相关的系统状态变量个数，n₁+n₂＝n，n为描述电力系统动态特性的微分代数方程组中状态变量的个数，t为当前时刻，τ_i>0，i＝1,…,m，为时滞常数，并假设0<τ₁<τ₂<…<τ_i…<τ_m＝τ_max为时滞环节的时滞常数，其中m为系统中时滞的个数，τ_max为最大的时滞，

和

为系统状态矩阵

的不同分块，均为稠密矩阵，

为稠密的系统状态矩阵，

为系统时滞状态矩阵，为稀疏矩阵，Δx₁(t)和Δx₂(t)分别为t时刻系统状态变量x₁和x₂的增量，Δx₂(t-τ_i)为t-τ_i时刻x₂的增量，

和

分别为t时刻系统状态变量x₁和x₂导数的增量，Δx_0,1和Δx_0,2分别为系统状态变量x₁和x₂的初始值，可分别记为

和

所述低阶时滞电力系统模型对应的线性化系统的特征方程为：

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

为系统状态矩阵，

为与不同时滞相对应的系统时滞状态矩阵。

2.如权利要求1所述的基于低阶SOD-LMS算法的时滞电力系统特征值计算方法，其特征是，根据低阶时滞电力系统模型，将其转化为对应的解算子，具体为：

设

为由区间[-τ_max,0]到n维实数空间

映射的连续函数构成的巴拿赫空间，解算子

定义为将空间X上的初值条件

转移到h时刻之后时滞电力系统解分段的线性算子；其中，h为转移步长，0≤h≤τ_max；

其中，s为积分变量，θ为变量，

表示对状态变量Δx求导，

和

分别为0和h+θ时刻时滞电力系统的状态；

解算子

的特征值μ与时滞系统的特征值λ之间具有如下关系：

λ＝1/(hlnμ),μ∈σ(T(h))\{0} (8)

式中：

表示解算子的谱，h为k步LMS法的步长，\表示集合差运算。

3.如权利要求2所述的基于低阶SOD-LMS算法的时滞电力系统特征值计算方法，其特征是，与解算子

对应的离散化矩阵T_N表示如下：

式中：

和

分别为对应阶数的零矩阵，I_(k-1)n和

分别为对应阶数的对角单位矩阵；

式(10)和(11)中，N为离散化维数，N＝Q+k+s_-，其中k为LMS法的系数，s_-为Nordsieck插值的系数，而

h为k步LMS法的步长，α_k和β_k为k步LMS法的系数，

其中，

中的元素完全由LMS法的系数α_k决定，

中的元素由LMS法的系数β_k和步长h共同决定，而

中的元素由拉格朗日插值系数和LMS法的系数β_k和步长h共同决定，

和

分别为对应阶数的零矩阵，

和I_n分别为对应阶数的对角单位矩阵，

为Kronecker积运算。

4.如权利要求3所述的基于低阶SOD-LMS算法的时滞电力系统特征值计算方法，其特征是，利用坐标旋转预处理方法，将解算子的离散化矩阵旋转并进行近似得到旋转后的解算子离散化近似矩阵，具体为：

将坐标轴逆时针旋转θ角度，将低阶时滞电力系统模型对应的线性化系统的特征方程中的λ用λ′e^jθ代替，可以得到坐标轴旋转后的特征方程：

式中：

λ′＝λe^-jθ (13)

τ_i′＝τ_ie^jθ,i＝1,...,m (14)

对坐标旋转后的τ′_i，i＝1,…,m，进行近似，因此，式(14)变为：

τ′_i＝τ_ie^jθ≈τ_i,i＝1,...,m (16)

式(12)变为：

式中：

为和

分别为λ′和与其对应的右特征向量v的近似值；

坐标轴旋转预处理后，λ'与其

的特征值μ'之间的映射关系为：

5.如权利要求4所述的基于低阶SOD-LMS算法的时滞电力系统特征值计算方法，其特征是，为了改善算法的收敛性，对特征值μ'进行非线性放大，从而增大不同特征值之间的相对距离，假设对μ'进行α次乘方，则由式(18)得：

式中：

由式(19)可知，保持

和τ_i，i＝1,…,m，不变，将h增大α倍即可实现旋转-放大预处理，这表明，对μ'的放大处理，等价于将转移步长h增大α倍；

由于转移步长h被变换原来的α倍，区间[-τ_max,0]被重新划分为长度等于或小于αh的Q'个子区间，

N变为N'＝Q'+k+s_-，

被重新形成为

T_N变为T_N′：

式中：

T_N′的谱与μ"之间有如下对应关系：

μ″∈σ(T′_N)\{0} (26)

其中，T_N为SOD-LMS方法所形成的解算子的离散化矩阵，T_N′为经过旋转放大预处理后的T_N矩阵，Γ_N′描述当前时刻t到时刻t+h的系统状态转移关系，为解算子离散化矩阵T_N′的子矩阵，A′₁₁，A′₁₂，A′₂₁，A′₂₂，

分别为经过旋转放大预处理后的相应矩阵，

分别为经过旋转放大预处理后的相应向量。

6.如权利要求5所述的基于低阶SOD-LMS算法的时滞电力系统特征值计算方法，其特征是，采用迭代特征值算法计算T_N′的设定个数个模值最大的特征值在应用解算子对应的离散化矩阵T_N′求解大规模电力系统的时滞特征值时步骤为：

设第j个Krylov向量表示为

则第j+1个向量q_j+1可通过矩阵T_N′与向量q_j的乘积计算得到：

由于T_N′具有特殊的逻辑结构，q_j+1的第n+1:(kn₁+N'n₂)个分量与q_j的第1:(k-1)n个和第kn-n₂+1:(kn₁+(N'-1)n₂)个分量之间存在一一对应关系，即：

而q_j+1的第1:n个分量，即q_j+1(1:n,1)，可以进一步分解为2个矩阵-向量乘积运算：

z＝Σ_N′·q_j (29)

式中：

为中间向量。

7.如权利要求6所述的基于低阶SOD-LMS算法的时滞电力系统特征值计算方法，其特征是，利用稀疏实现方法计算式(29)，首先需要从列的方向上将向量q_k压缩为矩阵

和

即q_k＝[vec(Q₁)^T,vec(Q₂)^T]^T；

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

式(31)中，z的计算量主要由稠密矩阵-向量乘积

决定，Q₁为与时滞无关的电力系统状态矩阵，Q₂为与时滞相关的电力系统状态矩阵，为了减少复数运算，将式(15)代入，则其可以改写为

其中

可将其计算步骤分解为如下两个步骤：

D₀w＝-C₀q_k+1 (32)

z₀＝A′₀q_k+1+B′₀w (33)

将其改写为矩阵形式，可得：

为稠密的系统状态矩阵，可由系统增广状态矩阵

和

计算得到：

状态变量划分后，系统的增广状态矩阵

可以写成分块形式：

式中：

分别为

的分块子矩阵，A₀、B₀、C₀和D₀为系统的增广状态矩阵。

8.如权利要求7所述的基于低阶SOD-LMS算法的时滞电力系统特征值计算方法，其特征是，利用稀疏实现方法计算式(30)，首先，考虑将

的计算式：

和式(15)代入式(24)中，得：

式中：A′₀,B′₀分别为A₀,B₀的变换矩阵，A′₀,B′₀与A₀,B₀具有完全相同的稀疏特性，

A′₀＝α_kI_n-hβ_kαA₀e^-jθ (36)

B′₀＝-hβ_kαB₀e^-jθ (37)

然后利用矩阵之和的求逆公式计算

可得：

于是，式(30)中

的稀疏实现可以分解为如下两个步骤：

[D₀-C₀(A′₀)^-1B′₀]w＝-C₀(A′₀)^-1z (39)

q_k+1＝(A′₀)^-1(z-B′₀w) (40)

将其改写为矩阵形式，可得：

9.如权利要求5所述的基于低阶SOD-LMS算法的时滞电力系统特征值计算方法，其特征是，在计算得到μ″之后，依次经过谱映射、坐标反旋转和牛顿校验后得到时滞电力系统的特征值λ，在牛顿校验前，估计特征值

的计算公式为：

Images