CN104615882B - 基于eigd的大规模时滞电力系统特征值计算方法 - Google Patents

Abstract

本发明公开了基于EIGD的大规模时滞电力系统特征值计算方法，步骤如下：建立时滞电力系统模型，进行时滞电力系统模型的公式转化，转换为抽象柯西问题，进而将计算时滞电力系统模型的特征值转化成计算转化后公式的无穷小生成元的特征值；对无穷小生成元进行离散化，得到无穷小生成元的近似矩阵；近似矩阵通过位移处理，得到位移处理后的近似矩阵，位移处理后的近似矩阵经过逆变换，得到逆矩阵；从而将要求计算的部分特征值变为模值最大的部分特征值；采用Arnoldi算法来计算得到的逆矩阵的模值最大的部分特征值；从而得到近似矩阵的特征值为；利用牛顿迭代法进行校验，通过计算得到时滞电力系统精确的特征值和特征向量。

Description

基于EIGD的大规模时滞电力系统特征值计算方法

技术领域

本发明涉及一种基于EIGD的大规模时滞电力系统特征值计算方法，EIGD的中文解释：显示无穷小生成元离散化，EIGD的英文全称：Explicit Infinitesimal GeneratorDiscretization。

背景技术

广域测量系统(Wide-Area Measurement System，WAMS)的出现给大规模互联电力系统稳定分析与控制的发展带来新的契机。基于WAMS提供的广域信息的互联电网低频振荡控制，通过引入有效反映区间振荡模式的广域反馈信号，能够获得较好的阻尼控制性能，其为解决互联电网中的区域间低频振荡问题，进而提高系统的输电能力提供了新的控制手段，具有良好而又广泛的应用前景。

广域信号在由不同通信介质(如光纤、电话线、数字微波、卫星等)组成的WAMS通信网络中传输和处理时，存在几十到几百毫秒间变化的通信延时。时滞是导致系统控制律失效、运行状况恶化和系统失稳的一种重要诱因[Wu H X,Tsakalis K S,Heydt GT.Evaluation of time delay effects to wide-area power system stabilizerdesign.IEEE Trans.Power Syst.,2004,19(4):1935-1941.]。因此，利用广域测量信息进行电力系统闭环控制时，必须计及时滞的影响。

目前，考虑广域通信时滞的电力系统小干扰稳定性分析方法，总体上可分为时域方法和频率方法两大类。

(1)时域方法

时域方法，主要基于时滞依赖的Lyapunov-Krasovskii稳定定理和Razumikhin定理，并利用线性矩阵不等式(Linear Matrix Inequality,LMI)技术来设计鲁棒控制器，在保证系统稳定性的同时，能够得到系统能够承受的时滞上限。时域方法的缺点主要有以下两点：(1)基于Lyapunov的时滞依赖的稳定性判据，仅为系统稳定的充分条件，存在固有的保守性；(2)当应用于大规模电力系统时，需要与模型降阶技术相结合以降低计算量。但是，系统模型的精确性对稳定性的判别存在一定的影响。

(2)频域方法

●函数变换方法

在频域范围内，当考虑时滞影响后，电力系统线性化微分-代数方程对应的特征方程中标量时滞被转化为指数项。因而，特征方程是一个超越方程。为此，在求取电力系统的时滞稳定裕度时，通常采用函数变换的方法(如Rekasius变换[Rifat S,Nejat O.A novelstability study on multiple time-delay systems(MTDS)using the root clusteringparadigm.Proc of the American Control Conference,Boston,MA,2004,5422-5427.]、Lambert-W函数[Yi S,Nelson P W,Ulsoy A G.Time-delay systems:analysis andcontrol using the Lambert W function.Singapore:World Scientific PublishingCompany,2010:]、Pade近似，等)对超越项进行变换，避免直接求解特征方程的困难。

函数变换方法存在如下不足：(1)Rekasius变换方法，可以将超越的特征方程改写为一个代数方程，从而能够计算得到虚轴上所有特征值。然而，该方法不能得到复平面上其他特征值的分布。(2)当系统存在对称时滞且系统矩阵能够同时被三角化时，时滞系统的特征谱可以利用Lambert-W函数显示表示。可见，该方法存在较强的前提和假设条件。(3)Pade近似是一种常用的时滞环节处理方法。通过有理多项式来逼近时滞环节，进而可以方便地利用经典和现代控制理论设计广域阻尼控制律，并通过时域仿真验证控制器的有效性。然而，Pade近似的准确性对阻尼控制器设计的影响却未被深入分析和研究。

●谱离散化方法

近年来，在数值分析和计算数学领域，谱离散化方法大量研究用于计算时滞系统的部分特征值。文献[Milano F,Anghel M.Impact of time delays on power systemstability.IEEE Trans.Circuit Syst.-I,2012,59(4):889-900.]利用无穷小生成元离散化方法计算励磁系统中发电机端电压量测环节存在单一时滞情况下，系统的部分特征值。文献[Liang H,Choi B J,Zhuang W H,等.Stability enhancement of decentralizedInverter control through wireless communications in microgrids.IEEETrans.Smart Grid,2013,4(1):321-331.]利用无穷小生成元离散化方法分析在廉价通信方案下，当微电网中逆变器控制环节存在时滞时，系统的小干扰稳定性。

谱离散化方法存在如下不足：(1)无穷小生成元的近似矩阵，没有清晰、解析的表达形式，难以利用稀疏特征值方法计算得到系统的部分特征值；(2)没有利用系统线性化增广矩阵的稀疏特性，计算量较大，不能用于大规模电力系统时滞特征值的计算。

发明内容

本发明的目的就是为了解决上述问题，提供一种基于EIGD的大规模时滞电力系统特征值计算方法，它具有计算量小，内存要求低，计算速度快，计算效率高的优点，能够用于考虑广域时滞影响的大规模电力系统的特征值计算和时滞稳定性判别。

为了实现上述目的，本发明采用如下技术方案：

基于EIGD的大规模时滞电力系统特征值计算方法，包括如下步骤：

步骤(1)：建立时滞电力系统模型，依据时滞电力系统模型的状态变量在设定区间段内是连续可微的，进行时滞电力系统模型的公式转化，转换为抽象柯西问题，进而将计算时滞电力系统模型的特征值转化成计算转化后公式的无穷小生成元的特征值；为了得到无穷小生成元的特征值，对无穷小生成元进行离散化，从而得到无穷小生成元的近似矩阵

步骤(2)：近似矩阵通过位移处理，得到位移处理后的近似矩阵位移处理后的近似矩阵经过逆变换，得到逆矩阵从而将复平面上距离某个设定点最近的、设定数量个数的特征值，变换为设定数量个数的模值最大特征值；

步骤(3)：采用隐式重启动Arnoldi算法来计算步骤(2)得到的逆矩阵的模值最大特征值λ"；从而得到近似矩阵的特征值为λ；

步骤(4)：利用牛顿迭代法对步骤(3)得到的特征值为λ进行校验，通过计算得到时滞电力系统精确的特征值和特征向量。

所述步骤(1)的时滞电力系统模型如下：

式中：t为当前时刻；Δx(t)∈R^n×1为增量形式的电力系统状态变量向量，n为系统状态变量总数；表示电力系统状态变量向量的导数；τ_i>0为第i个时滞环节的时滞常数，i＝1,2,…,m；其中，最大的时滞表示为τ_max；为系统状态矩阵，为一个稠密矩阵；(i＝1,…,m)为系统时滞状态矩阵，为稀疏矩阵；

稠密矩阵和稀疏矩阵具体表示为：

式中：Α_i∈R^n×n，B_i∈R^n×l，C_i∈R^l×n，D_i∈R^l×l为线性化系统矩阵，l为系统代数变量总数。

所述步骤(1)的抽象柯西问题：

式中：u(t)为一个连续可微函数，且u(t)＝x(t+θ)，θ∈[-τ_max,0]。

所述步骤(1)的无穷小生成元的近似矩阵的计算公式为：

式中：Π_N∈R^{(N+1)n×(N+1)n}为伴随矩阵，Σ_N∈R^{(N+1)n×(N+1)n}为块上三角矩阵；

伴随矩阵Π_N和块上三角矩阵Σ_N具体表示为：

式中：表示克罗内克积；R_j(j＝0,1,…,N)由系统状态矩阵经过代数运算得到，即：

所述步骤(2)的位移处理之前需要做的步骤如下：

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

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

式(2)的等价增广形式为：

式中：w∈R^l×1为中间和辅助向量；

设I_n为n阶单位矩阵，则A'(λ)和B'(λ)具体表示为：

所述步骤(2)的位移处理步骤如下：

给定位移点s，将式(2)中的λ用λ'+s代替，则得到位移后的特征方程：

式中：

λ′＝λ-s

所述步骤(2)中：

设表示经过位移处理后的近似矩阵；逆矩阵的计算过程为

式中：(Σ′_N)^-1显示地表示如下：

(Σ′_N)^-1＝[Σ′₀；[0_Nn,I_Nn]]

Σ′₀＝(R₀′)^-1[I_n,R₁′,…,R′_N]

式中：0_Nn∈R^Nn×n为零矩阵，I_Nn∈R^Nn×n为单位矩阵；[；]和[,]分别表示在列和行的方向上进行向量或矩阵连接；A'(s)和B'(s)将式(4)和(5)中的λ用s代替后得到。

所述步骤(3)的步骤如下：

设第k个Krylov向量表示为q_k，则矩阵向量乘积的流程如下：

步骤(3-1)：利用稀疏矩阵乘法，计算得到

步骤(3-2)：从列的方向上，将r压缩为矩阵其中

步骤(3-3)：赋值z＝v₀；

步骤(3-4)：赋值j＝1；

步骤(3-5)：计算

步骤(3-6)：计算

步骤(3-7)：求和z＝z+p+w；

步骤(3-8)：将j增加1，重复步骤(3-5)-(3-7)，直至j＝N；

步骤(3-9)：计算q_k+1(1:n)＝(R₀′)^-1z；

步骤(3-10)：q_k+1((n+1):(N+1)n)＝r((n+1):(N+1)n)。

在步骤(3-10)中，q_k+1(i₁:i₂)和r(i₁:i₂)分别表示抽取q_k+1和r的第i₁到i₂个元素形成的列向量。

设利用隐式重启动Arnoldi算法计算得到的特征值为λ"，则稀疏近似矩阵对应的特征值为λ＝1/λ"+s＝λ′+s。此外，与λ对应的Krylov向量的前n个分量是特征向量v的良好的估计和近似。

所述步骤(4)的步骤为：

将式(3)线性化，得到牛顿迭代的修正方程式如下：

式中：λ^(k)为第k次迭代时特征值，v^(k)为第k次迭代时特征向量，w^(k)为第k次迭代时中间辅助向量的估计值；Δλ^(k)为第k次迭代时特征值的修正量；Δv^(k)为第k次迭代时特征向量的修正量；Δw^(k)为第k次迭代时中间辅助向量的估计值的修正量；f^(k)表示第k次迭代时增广形式特征方程的不平衡量。

为了保证特征向量v^(k)的唯一性，假设其第一个分量并在牛顿迭代过程中保持不变，相应地，对v^(k)中的其他元素进行规格化，为了避免式(10)表示的修正方程降阶并保留其稀疏特性，将Δv^(k)的第一个分量用Δλ^(k)代替。相应地，用g^'(k)代替J^'(k)的第一列，则得到修正后的雅可比矩阵为J^"(k)＝[g^'(k),J^'(k)(1:(n+l),2:(n+l))]；

经过上述处理后，修正方程式(10)变为：

式中：J^(k)为将J^"(k)的第一列用单位向量e₁代替后得到的矩阵，其与J^'(k)具有相同的稀疏结构；g^(k)为g^'(k)与e₁之差；

利用Sherman-Morrisony公式，得式(11)的解：

牛顿迭代收敛以后，v′^(k)的第一个元素就是待求的时滞电力系统的精确特征值λ。

所述步骤(4)如果计算得到的全部特征值λ都具有负实部，则时滞电力系统在稳态运行点处是小扰动稳定的；反之，若至少存在一个具有正实部的特征值，则系统在该点处是小扰动不稳定的。

本发明的有益效果：

1由上述流程看出，步骤(3-5)和步骤(3-9)的计算量最大、耗时最长。从而可以推知，整个流程的计算量近似正比于集合Ω_N中的离散点数N。

计算化无穷小生成元的近似矩阵和它的逆矩阵所具有的显示、稀疏特性，使得可以充分利用系统线性化增广系数矩阵Α_i，B_i，C_i，D_i(i＝0,1,…,m)的稀疏特性，以稀疏的方式实现步骤(3-5)和(3-9)，从而大大减小计算量和内存要求，提高了计算效率，能够用于考虑广域时滞影响的大规模电力系统特征值计算和时滞稳定性判别。

2能够高效、准确地计算得到，用于判别电力系统时滞稳定性判别所需的、复平面上最右侧的部分特征值。此外，由于不需要计算出系统的所有特征值，就可以实现电力系统时滞稳定性的判别。总起来说，该方法具有很高的计算效率。

附图说明

图1为本发明的流程图；

图2为四机两区域算例系统；

图3为N＝50时，和的特征值；

图4为N＝50时，和的特征值。

具体实施方式

下面结合附图与实施例对本发明作进一步说明。

1.时滞电力系统模型

考虑广域通信时滞影响后，电力系统可以如下的一组时滞微分方程组来描述：

式中：x∈R^n×1为电力系统的状态变量向量，n为系统状态变量总数。t为当前时刻，τ_i>0为第i个时滞环节的时滞常数，i＝1,2,…,m。其中，最大的时滞表示为τ_max。为系统状态矩阵，为一个稠密矩阵；为系统时滞状态矩阵，为稀疏矩阵。它们可具体表示为：

式中：Α_i∈R^n×n，B_i∈R^n×l，C_i∈R^l×n，D_i∈R^l×l为线性化系统矩阵，l为系统代数变量总数。

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

式中：λ和v分别为特征值和相应的右特征向量。

式(2)的等价增广形式为：

式中：w∈R^l×1为中间和辅助向量。设I_n为n阶单位矩阵，则A'(λ)和B'(λ)可具体表示为：

2.基于显示无穷小生成元离散化的时滞特征值计算方法------基本理论

由于时滞电力系统的状态变量x(t)在区间[-τ_max,∞)上是连续可微的，可以证明式(1)可以转化为下列抽象柯西问题：

式中：u(t)为一个连续可微函数，且u(t)＝x(t+θ),θ∈[-τ_max,0]。

由谱映射定理可知，式(1)表示的时滞系统的特征值等于无穷小生成元的特征值，因此可以通过计算的特征值，来得到时滞系统的特征值。然而，的特征值本质上是一个无穷维问题。因此，需要对进行离散化，从而得到其近似矩阵

给定正整数N，区间[-τ_max,0]上N+1个不同离散点形成的集合Ω_N定义为Ω_N:＝{θ_N,j,j＝1,…,N+1}，且有-τ_max≤θ_N,j≤…≤θ_N,N+1＝0。集合Ω_N中的N个非零离散点，选取为经过位移和归一化处理后的第N阶第一类切比雪夫多项式的零点。在X_N上，式(6)表示的抽象柯西问题可以被离散化，进而可以得到无穷小生成元的近似矩阵

式中：Π_N∈R^{(N+1)n×(N+1)n}为伴随矩阵，Σ_N∈R^{(N+1)n×(N+1)n}为块上三角矩阵。Π_N和Σ_N可具体表示为：

式中：表示克罗内克积。R_j(j＝0,1,…,N)由系统状态矩阵经过代数运算得到，即：

3.基于显示无穷小生成元离散化的时滞特征值计算方法------稀疏实现

3.1位移逆变换

电力系统机电振荡模式对应的特征值，通常位于复平面虚轴附近，决定着电力系统的小干扰稳定性。与其他特征值相比，这部分特征值的模值较小。因此，需要利用位移逆变换将近似矩阵从而将要求计算的部分特征值变为模值最大的部分特征值。

给定位移点s，将式(2)中的λ用λ'+s代替，则可以得到位移后的特征方程：

式中：

λ′＝λ-s

设表示经过位移处理后的近似矩阵。其逆矩阵为：

式中：(Σ′_N)^-1可以显示地表示如下：

(Σ′_N)^-1＝[Σ′₀；[0_Nn,I_Nn]]

Σ′₀＝(R₀′)^-1[I_n,R₁′,…,R′_N]

式中：0_Nn∈R^Nn×n为零矩阵，I_Nn∈R^Nn×n为单位矩阵。[；]和[,]分别表示在列和行的方向上进行向量或矩阵连接。A'(s)和B'(s)可以将式(4)和(5)中的λ用s代替后得到。

3.2矩阵向量乘积的稀疏实现

采用隐式重启动Arnoldi算法来计算经过位移逆变换后近似矩阵模值最大的部分特征值。在Arnoldi算法中，耗时最长、计算量最大的操作就是在形成Krylov向量过程中的矩阵向量乘积运算。设第k个Krylov向量表示为q_k，则矩阵向量乘积的流程如下：

步骤1：利用稀疏矩阵乘法，计算得到

步骤2：从列的方向上，将r压缩为矩阵其中

步骤3：赋值z＝v₀；

步骤4：赋值j＝1；

步骤5：计算

步骤6：计算

步骤7：求和z＝z+p+w；

步骤8：将j增加1，然后重复步骤5-步骤7，直至j＝N；

步骤9：计算q_k+1(1:n)＝(R₀′)^-1z；

步骤10：q_k+1((n+1):(N+1)n)＝r((n+1):(N+1)n)。

在步骤10中，q_k+1(i₁:i₂)和r(i₁:i₂)分别表示抽取q_k+1和r的第i₁到i₂个元素形成的列向量。

由上述流程可以看出，步骤5和步骤9的计算量最大、耗时最长。从而可以推知，整个流程的计算量近似正比于集合Ω_N中的离散点数N。

利用系统线性化增广系数矩阵的稀疏特性，步骤5和9可以以稀疏的方式实现，从而大大减小计算量，提高计算效率。

设利用隐式重启动Arnoldi算法计算得到的特征值为λ"，则稀疏近似矩阵对应的特征值为λ＝1/λ"+s＝λ′+s。此外，与λ对应的Krylov向量的前n个分量是特征向量v的良好的估计和近似。

3.3牛顿校验

为了改善计算得到的特征值和特征向量的精度，可以利用牛顿法，计算得到时滞电力系统精确的特征值和相应的特征向量。

将式(3)线性化，得到牛顿迭代的修正方程式如下：

式中：λ^(k)、v^(k)和w^(k)分别为第k次迭代时特征值、特征向量和中间辅助向量的估计值。前缀Δ表示它们各自的修正量。f^(k)表示第k次迭代时增广形式特征方程的不平衡量。

为了保证特征向量v^(k)的唯一性，假设其第一个分量并在牛顿迭代过程中保持不变。相应地，对v^(k)中的其他元素进行规格化，其值可能大于1。为了避免式(10)表示的修正方程降阶并保留其稀疏特性，将Δv^(k)的第一个分量用Δλ^(k)代替。相应地，用g^'(k)代替J^'(k)的第一列，则得到修正后的雅可比矩阵为J^"(k)＝[g^'(k),J^'(k)(1:(n+l),2:(n+l))]。经过上述处理后，修正方程式(10)变为：

式中：J^(k)为将J^"(k)的第一列用单位向量e₁代替后得到的矩阵，其与J^'(k)具有相同的稀疏结构。g^(k)为g^'(k)与e₁之差。

利用Sherman-Morrisony公式，可得式(11)的解：

牛顿迭代收敛以后，v′^(k)的第一个元素就是待求的时滞电力系统的精确特征值λ。

如果计算得到的全部特征值λ都具有负实部，则时滞电力系统在稳态运行点处是小扰动稳定的；反之，若至少存在一个具有正实部的特征值，则系统在该点处是小扰动不稳定的。

4.算例验证

利用四机两区算例系统，来验证本文提出的基于显示无穷小生成元离散化的大规模时滞电力系统的稀疏特征值计算方法的有效性。所有的分析均在Matlab中和在InterCore i5 4×3.4GHz 8GB RAM个人计算机上进行。

四机两区算例系统如图2所示。所有发电机采用高增益晶闸管励磁系统，并装设以转速偏差为输入的电力系统稳定器(PSS)。在此基础上，在G₁上装设以G₁和G₃相对转速偏差Δω₁₃为反馈信号的广域PSS，以进一步提高对区域间低频振荡模式的阻尼。。假设广域PSS的反馈时滞和输入时滞分别为τ₁＝150ms和τ₂＝100ms。系统状态变量和代数变量的维数分别为n＝56和l＝22。

4.1稀疏近似矩阵的逼近能力分析

为了验证本发明提出的基于显示无穷小生成元离散化的大规模时滞电力系统的稀疏特征值计算方法的有效性，首先分析无穷小生成元离散化得到的稀疏近似矩阵逼近无穷小生成元的能力和准确性。在四机两区系统中，当N＝50时，的维数为(N+1)n＝2856。首先，利用QR算法计算得到的全部特征值，如图3所示。接着，将这些特征值作为牛顿法的初始值，进行校验。其中一部分特征值能够收敛为无穷小生成元的特征值，如图3所示。

下面对和的特征值进行比较和分析。由图3可知，在Re(λ)∈[-48,0]、Im(λ)∈[0,820]的区域内，以及点(-100,0)附近，和的特征值能较好地吻合。图3中，剩余特征值为的特征值，经过牛顿迭代后，它们不能收敛为的相应的特征值。但是，它们均位于已经收敛的那部分特征值的左侧。通过上述分析可知，在一定程度上能较准确地逼近但这对于电力系统小干扰稳定性分析而言已经足够。

4.2方法的准确性分析

利用隐式重启动Arnoldi算法计算稀疏近似矩阵的部分特征值。由于四机两区算例系统的规模较小，可以通过为算法选择合适的参数，从而一次性得到实部大于某个数值、系统最右侧的部分特征值。当N＝50时，位移点选择为s＝0.05，计算得到其附近的r＝200个特征值，如4所示。由图可知，Arnoldi算法能够准确计算得到频率近似矩阵最大的特征值非常接近无穷小生成元的第19个特征值λ₁₉。

下面对本发明所提出的方法和文献[Breda D,Maset S,VermiglioR.Pseudospectral differencing methods for characteristic roots of delaydifferential equations.SIAM Journal on Scientific Computing,2005,27(2):482-495.]提出的基于隐式无穷小生成元离散化的特征值计算方法(Implicit InfinitesimalGenerator Discretization,IIGD方法)进行对比。利用IIGD方法计算得到s＝0.05附近200个特征值，如图4所示。与本发明提出的基于显示无穷小生成元离散化的特征值计算方法(EIGD)的结果相比，差别在于IIGD方法能够准确的求得的频率最大的特征值为λ₂₀，比EIGD求得的特征值多一个。

1.1计算时间分析

四种情况下，即N＝20，40，50和60时，利用QR计算近似矩阵的全部特征值，以及利用非稀疏IGD和稀疏IGD方法分别计算近似矩阵和在位移点s＝0.05附近r＝200个特征值。计算表1所示。

表1 QR、IIGD和EIGD算法所用的计算时间比较

	N＝20	N＝40	N＝50	N＝60
					QR	1.97s	13.71s	27.06s	47.70s
IIGD	10.14s	25.84s	34.96s	45.81s
					EIGD	6.24s	10.61s	12.69s	15.21s

首先，当N＝20时，虽然QR算法计算了稀疏近似矩阵的全部特征值，但是由于其在系统规模较小的情况下非常鲁棒且收敛收敛速度快，其比EIGD算法计算200个特征值所用的时间少得多。当N增加时，QR算法的计算量和计算时间急剧增加。相比较而言，EIGD算法的计算时间增加比较平缓。这是因为，在EIGD方法中，矩阵向量乘积的计算量与N大体上呈线性关系。其次，当N＝40，50和60时，EIGD方法所需的计算时间比QR和IIGD方法少。

上述虽然结合附图对本发明的具体实施方式进行了描述，但并非对本发明保护范围的限制，所属领域技术人员应该明白，在本发明的技术方案的基础上，本领域技术人员不需要付出创造性劳动即可做出的各种修改或变形仍在本发明的保护范围以内。

Claims (10)

1.基于EIGD的大规模时滞电力系统特征值计算方法，其特征是，包括如下步骤：

步骤(1)：建立时滞电力系统模型，依据时滞电力系统模型的状态变量在设定区间段内是连续可微的，进行时滞电力系统模型的公式转化，转换为抽象柯西问题，进而将计算时滞电力系统模型的特征值转化成计算转化后公式的无穷小生成元的特征值；为了得到无穷小生成元的特征值，对无穷小生成元进行离散化，从而得到无穷小生成元的近似矩阵

步骤(2)：近似矩阵通过位移处理，得到位移处理后的近似矩阵位移处理后的近似矩阵经过逆变换，得到逆矩阵从而将复平面上距离某个设定点最近的、设定数量个数的特征值，变换为设定数量个数的模值最大特征值；

步骤(3)：采用隐式重启动Arnoldi算法来计算步骤(2)得到的逆矩阵的模值最大特征值λ"；从而得到近似矩阵的特征值为λ；

步骤(4)：利用牛顿迭代法对步骤(3)得到的特征值λ进行校验，通过计算得到时滞电力系统精确的特征值和特征向量。

2.如权利要求1所述的基于EIGD的大规模时滞电力系统特征值计算方法，其特征是，所述步骤(1)的时滞电力系统模型如下：

式中：t为当前时刻；Δx(t)∈R^n×1为增量形式的电力系统状态变量向量，n为系统状态变量总数；表示电力系统状态变量向量的导数；τ_i>0为第i个时滞环节的时滞常数，i＝1,2,…,m；其中，最大的时滞表示为τ_max；为系统状态矩阵，为一个稠密矩阵；i＝1,…,m为系统时滞状态矩阵，为稀疏矩阵；

稠密矩阵和稀疏矩阵具体表示为：

式中：Α_i∈R^n×n，B_i∈R^n×l，C_i∈R^l×n，D_i∈R^l×l为线性化系统矩阵，l为系统代数变量总数。

3.如权利要求2所述的基于EIGD的大规模时滞电力系统特征值计算方法，其特征是，所述步骤(1)的抽象柯西问题：

式中：u(t)为一个连续可微函数，且u(t)＝x(t+θ)，θ∈[-τ_max,0]。

4.如权利要求3所述的基于EIGD的大规模时滞电力系统特征值计算方法，其特征是，所述步骤(1)的无穷小生成元的近似矩阵的计算公式为：

式中：Π_N∈R^{(N+1)n×(N+1)n}为伴随矩阵，Σ_N∈R^{(N+1)n×(N+1)n}为块上三角矩阵；

伴随矩阵Π_N和块上三角矩阵Σ_N具体表示为：

式中：表示克罗内克积；R_j由系统状态矩阵和系统时滞状态矩阵经过代数运算得到，i＝1,…,m，j＝0,1,…,N；即：

5.如权利要求4所述的基于EIGD的大规模时滞电力系统特征值计算方法，其特征是，所述步骤(2)的位移处理之前需要做的步骤如下：

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

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

式(2)的等价增广形式为：

式中：w∈R^l×1为中间和辅助向量；

设I_n为n阶单位矩阵，则A'(λ)和B'(λ)具体表示为：

6.如权利要求5所述的基于EIGD的大规模时滞电力系统特征值计算方法，其特征是，所述步骤(2)的位移处理步骤如下：

给定位移点s，将式(2)中的λ用λ'+s代替，则得到位移后的特征方程：

式中：

λ′＝λ-s

7.如权利要求6所述的基于EIGD的大规模时滞电力系统特征值计算方法，其特征是，所述步骤(2)中：

设表示经过位移处理后的近似矩阵；逆矩阵的计算过程为

式中：(Σ′_N)^-1显示地表示如下：

(Σ′_N)^-1＝[Σ′₀；[0_Nn,I_Nn]]

Σ′₀＝(R′₀)^-1[I_n,R′₁,…,R′_N]

式中：0_Nn∈R^Nn×n为零矩阵，I_Nn∈R^Nn×n为单位矩阵；[；]和[,]分别表示在列和行的方向上进行向量或矩阵连接；A'(s)和B'(s)将式(4)和(5)中的λ用s代替后得到。

8.如权利要求7所述的基于EIGD的大规模时滞电力系统特征值计算方法，其特征是，所述步骤(3)的步骤如下：

设第k个Krylov向量表示为q_k，则矩阵向量乘积的流程如下：

步骤(3-1)：利用稀疏矩阵乘法，计算得到

步骤(3-2)：从列的方向上，将r压缩为矩阵其中

步骤(3-3)：赋值z＝v₀；

步骤(3-4)：赋值j＝1；

步骤(3-5)：计算

步骤(3-6)：计算

步骤(3-7)：求和z＝z+p+w；

步骤(3-8)：将j增加1，重复步骤(3-5)-(3-7)，直至j＝N；

步骤(3-9)：计算q_k+1(1:n)＝(R′₀)^-1z；

步骤(3-10)：q_k+1((n+1):(N+1)n)＝r((n+1):(N+1)n)；

在步骤(3-10)中，q_k+1(i₁:i₂)和r(i₁:i₂)分别表示抽取q_k+1和r的第i₁到i₂个元素形成的列向量；

设利用隐式重启动Arnoldi算法计算得到的特征值为λ"，则稀疏近似矩阵对应的特征值为λ＝1/λ"+s＝λ′+s；此外，与λ对应的Krylov向量的前n个分量是特征向量v的良好的估计和近似。

9.如权利要求8所述的基于EIGD的大规模时滞电力系统特征值计算方法，其特征是，所述步骤(4)的步骤为：

将式(3)线性化，得到牛顿迭代的修正方程式如下：

式中：λ^(k)为第k次迭代时特征值，v^(k)为第k次迭代时特征向量，w^(k)为第k次迭代时中间辅助向量的估计值；Δλ^(k)为第k次迭代时特征值的修正量；Δv^(k)为第k次迭代时特征向量的修正量；Δw^(k)为第k次迭代时中间辅助向量的估计值的修正量；f^(k)表示第k次迭代时增广形式特征方程的不平衡量；

为了保证特征向量v^(k)的唯一性，假设其第一个分量在牛顿迭代过程中保持不变，相应地，对v^(k)中的其他元素进行规格化，为了避免式(10)表示的修正方程降阶并保留其稀疏特性，将Δv^(k)的第一个分量用Δλ^(k)代替；相应地，用g'^(k)代替J'^(k)的第一列，则得到修正后的雅可比矩阵为J"^(k)＝[g'^(k),J'^(k)(1:(n+l),2:(n+l))]；

经过上述处理后，修正方程式(10)变为：

式中：J^(k)为将J"^(k)的第一列用单位向量e₁代替后得到的矩阵，其与J'^(k)具有相同的稀疏结构；g^(k)为g'^(k)与e₁之差；

利用Sherman-Morrisony公式，得式(11)的解：

牛顿迭代收敛以后，v′^(k)的第一个元素就是待求的时滞电力系统的精确特征值λ。

10.如权利要求1所述的基于EIGD的大规模时滞电力系统特征值计算方法，其特征是，所述步骤(4)如果计算得到的全部特征值λ都具有负实部，则时滞电力系统在稳态运行点处是小扰动稳定的；反之，若至少存在一个具有正实部的特征值，则系统在该点处是小扰动不稳定的。