CN108321821A - 基于sod-irk的时滞电力系统稳定性判别方法 - Google Patents

Abstract

本发明公开了基于SOD‑IRK的时滞电力系统稳定性判别方法，包括：建立时滞电力系统模型；采用不同的隐式龙格‑库塔方法对解算子进行离散化，得到解算子的离散化矩阵；将无限维的特征值问题转化为有限维的特征值问题；采用隐式Arnoldi算法计算解算子的离散化矩阵模值最大的特征值的近似值；计算过程中采用旋转‑放大方法进行预处理，采用Kronecker积的性质进行稀疏实现；根据谱映射关系，将解算子离散化矩阵模值最大的特征值的近似值转化为时滞电力系统模型的近似特征值；采用牛顿迭代法对近似特征值进行修正，得到时滞电力系统的精确特征值；根据精确特征值阻尼比的大小来判断时滞电力系统的稳定性。

Description

基于SOD-IRK的时滞电力系统稳定性判别方法

技术领域

本发明涉及基于解算子龙格-库塔离散化(Solution Operator DiscretizationMethods with Implicit Runge-Kutta,SOD-IRK)的时滞电力系统稳定性判别方法。

背景技术

随着全球能源互联网的兴起，互联电力系统的规模逐渐增大，区间低频振荡问题更加显著。传统的解决方案是安装电力系统稳定器(Power System Stabilizer,PSS)，但是由于其反馈控制信号来源于当地，不能有效阻尼互联电力系统的区间振荡。

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

广域信号在由不同通信介质(如光纤、电话线、数字微波、卫星等)组成的WAMS通信网络中传输和处理时，存在几十到几百毫秒间变化的通信延时。时滞是导致系统控制律失效、运行状况恶化和系统失稳的一种重要诱因。因此，利用广域测量信息进行电力系统闭环控制时，必须计及时滞的影响。

在现代电力系统中，小干扰稳定主要关注的是机电振荡问题。以状态空间模型为基础的特征值分析法是研究机电振荡的强有力工具。目前，研究人员已经提出了许多计算大规模电力系统关键特征值子集的有效方法，主要包括基于降阶系统的选择模式分析法，AESOPS算法和S-矩阵法，幂法、牛顿法、Rayleigh商迭代法等序贯法以及同时迭代法、Arnoldi算法、重新分解的双重迭代和Jacobi-Davidson方法等子空间迭代法。这些方法在计算部分特征值时都利用了增广状态矩阵的稀疏性，多数方法都是通过对原系统进行谱变换从而改变特征谱的分布，然后求取系统的特征值，再通过反变换得到原系统的关键特征值。但是，以上提到的方法均未考虑时滞的影响。中国发明专利基于Padé近似的时滞电力系统特征值计算与稳定性判别方法.201210271783.8[P].利用Pade近似多项式逼近时滞环节，进而计算系统最右侧的关键特征值，并判断系统的时滞稳定性。中国发明专利基于EIGD的大规模时滞电力系统特征值计算方法.201510055743.3.[P].提出了一种基于显示IGD(Explicit IGD,EIGD)的大规模时滞电力系统特征值计算。利用计算得到的系统最右侧的关键特征值，可以判断系统在固定时滞下的稳定性。这些时滞稳定性判别方法，均需要通过多次扫描[0.1,2.5]Hz低频振荡频率范围内、靠近虚轴的关键特征值，才能判断系统的时滞稳定性。中国发明专利基于SOD-PS-R算法的时滞电力系统机电振荡模式计算方法.201510829129.8.[P].提出了一种基于解算子伪谱离散化(SOD-PS)的大规模时滞电力系统特征值计算方法，其离散化矩阵维数较大，计算效率较低。

发明内容

为了解决现有技术的不足，本发明提供了基于SOD-IRK的时滞电力系统稳定性判别方法；

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

基于SOD-IRK的时滞电力系统稳定性判别方法，包括：

步骤(1)：建立时滞电力系统模型；依据时滞电力系统模型的特征值与时滞电力系统模型的解算子特征值之间的谱映射关系，将计算时滞电力系统模型的特征值转化成计算解算子的特征值；从而将判断时滞电力系统稳定性的问题转化为计算解算子的模值最大的特征值问题；

步骤(2)：采用不同的隐式龙格-库塔方法对解算子进行离散化，得到解算子的离散化矩阵；将无限维的特征值问题转化为有限维的特征值问题；

步骤(3)：采用隐式Arnoldi算法计算步骤(2)得到的解算子的离散化矩阵模值最大的特征值的近似值；计算过程中采用旋转-放大方法进行预处理，采用Kronecker积的性质进行稀疏实现；

步骤(4)：根据谱映射关系，将解算子离散化矩阵模值最大的特征值的近似值转化为时滞电力系统模型的近似特征值；

步骤(5)：采用牛顿迭代法对近似特征值进行修正，得到时滞电力系统的精确特征值；

步骤(6)：根据精确特征值阻尼比的大小来判断时滞电力系统的稳定性。

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

式中：为电力系统的状态变量向量，n为系统状态变量总数。t为当前时刻。0<τ₁<τ₂<…<τ_i…<τ_m为时滞环节的时滞常数，其中最大的时滞为τ_m。为系统状态矩阵，为稠密矩阵。为系统时滞状态矩阵，为稀疏矩阵。Δx(t)为t时刻系统状态变量的增量，Δx(t-τ_i)为t-τ_i时刻系统状态变量的增量，为t时刻系统状态变量导数的增量。Δx(0)为系统状态变量的初始值(即初始条件)，并简写为

式中：为高度稀疏的系统增广状态矩阵。

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

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

式(1.3)所示时滞电力系统的解算子X→X定义为将X上t时刻初值条件映射为t+h时刻系统状态的线性算子，即：

式中：是由区间[-τ_max,0]到n维实数空间映射的连续函数构成的巴拿赫空间，并赋有上确界范数。h为转移步长，0≤h≤τ_max。

在频域，与时滞电力系统的谱之间存在如下指数关系

式中：μ为的特征值；σ(·)表示谱；\表示集合差运算。

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

首先，将时滞区间[-τ_max,0]分为长度为h＝τ_max/N的N个子区间，[θ_N,θ_N–1],…,[θ₁,θ₀]；然后，利用p阶s级IRK的s个横坐标c₁,…,c_s对每个子区间进行离散化；最终，得到具有Ns个元素的集合Ω_Ns：

用Ω_Ns代替[-τ_max,0]，从而连续空间X被转化为离散空间令表示t＝δh+θ_j+c_qh(δ＝0,1；j＝1,…,N；q＝1,…,s)处Δx(t)的近似值，则X_Ns上系统的状态变量可表示为Δx^(δ)，即：

式中：和分别表示定义在[-τ_max,0]和[h-τ_max,h]上的系统状态变量。

结合解算子定义(1.4)和时滞系统状态方程(1.1)，可知Δx_h(t)是一个分段函数，其离散化原理如下：

(2-1)当t∈[-τ_max,-h]，Δx_h(t)为初始状态利用解算子的转移特性(1.4)，可以得到t＝-jh+c_qh处Δx(t+h)的估计值为公式(1.8)：其中，q＝1,…,s；j＝2,…,N；

将式(1.8)写成向量形式，得：

(2-2)当t∈[-h,0]，t+h-τ_i<0，i＝1,…,m，故Δx_h(t)是常微分方程(1.10)的解。

SOD-IRK的基本思想，就是得到公式(1.9)和公式(1.10)后，利用不同的IRK方法推导得到不同的Δx⁽¹⁾与Δx⁽⁰⁾之间的显式表达式，得到解算子离散化矩阵T_Ns。所述IRK方法，包括：Radau-IIA方法、Gauss-Legendre/SDIRK方法、Radau-IA方法、Lobatto-ⅢA/C方法或Lobatto-ⅢB方法。

p阶s级IRK递推公式的系数(A,b,c)可用Butcher表示：

系数(A,b,c)有如下特点：

1)0≤c₁<…<c_s≤1；

2)如果方法是自洽的consistent，则

采用不同的IRK离散化方案，得到的解算子离散化矩阵T_Ns是不同的。接下来，以Radau-IIA方法为例，详细推导解算子离散化矩阵T_Ns的表达式；

Radau-IIA法递推公式的系数(A,b,c)具有如下特点：

1)A可逆；

2)0<c₁<…<c_s＝1；

3)b＝(a_s1,…,a_ss)。

利用IRK求解(1.10)必须首先计算得到Δx(c_kh-τ_i)，i＝1,…,m-1；k＝1,…,s；考虑到点集c_kh-τ_i通常不属于离散集合Ω_N或h+Ω_N，利用c_kh-τ_i附近的p+1个点处Δx的近似解进行中间拉格朗日插值，从而得到Δx(c_kh-τ_i)的近似值

式中：l_r(c_kh-τ_i)表示由点c_kh-τ_i附近p+1个点计算得到的第r个拉格朗日插值系数。

其中，表示大于或等于p/2的最小整数。γ_i为正整数且满足不等式：

令

式中：为由拉格朗日插值系数决定的常数矩阵，表示Kronecker积。

将(1.12)代入(1.10)，进而利用基于Radau-IIA离散化方案的IRK方法可得t＝–h+c_qh(q＝1,…,s)处Δx_h(t)的估计值

将(1.15)改写成向量形式，得：

式中：

将R_Ns和Σ_Ns写为常数矩阵与系统状态矩阵的Kronecker积之和的形式：

式中：为常系数矩阵，

综合(1.9)和(1.16)，可得Δx⁽¹⁾与Δx⁽⁰⁾之间的关系式：

Δx⁽¹⁾＝T_NsΔx⁽⁰⁾ (1.19)

式中：X_Ns→X_Ns为解算子的IRK离散化矩阵，

Γ_Ns＝[I_(N-1)sn 0_(N-1)sn×sn] (1.21)

至此，无穷维解算子的特征值问题已被转化为有限维离散化矩阵T_Ns的特征值问题。

其他IRK离散化方案

解算子的其他IRK离散化方案与Radau-IIA离散化方案的不同之处，对应两种处理方式，总结为以下两点：

第一种处理方式：当c_s≠1时，增加离散点c_s+1＝1以得到每个子区间最右端点处系统状态的估计值。

第二种处理方式：当c₁＝0时，每个子区间左端点与相邻子区间的右端点重合，即θ_j–1+c₁h＝θ_j+c_sh或θ_j–1+c₁h＝θ_j+c_s+1h，j＝2,…,N。因而，需要去掉重复的离散点。

Gauss-Legendre/SDIRK离散化方案:

Gauss-Legendre/SDIRK离散化的系数c满足：0<c₁<…<c_s<1。

在解算子的离散化过程中，需要采用上述第一种处理方式。

最终，T_Ns及其子矩阵R_Ns、Σ_Ns、Γ_Ns的维数分别为N(s+1)n×N(s+1)n、(s+1)n×(s+1)n、(s+1)n×N(s+1)n、(N–1)(s+1)n×N(s+1)n。

它们的表达式分别为：

Γ_Ns＝[I_(N-1)(s+1)n 0_{(N-1)(s+1)n×(s+1)n}] (1.24)

式中：

Radau-IA离散化方案:

Radau-ⅠA的系数c满足：0＝c₁<…<c_s<1。

在解算子的离散化过程中，需要采用上述第二种处理方式。最终，T_Ns及其子矩阵R_Ns、Σ_Ns、Γ_Ns的维数分别为(Ns+1)n×(Ns+1)n、(s+1)n×(s+1)n、(s+1)n×(Ns+1)n、(N–1)sn×(Ns+1)n。

Σ_Ns和Γ_Ns的表达式分别为：

式中：

Lobatto-ⅢA/C离散化方案:

Lobatto-ⅢA/C离散化的系数(A,b,c)有如下特点：

1)0＝c₁<…<c_s＝1；

2)b＝(a_s1,…,a_ss)^T。

在解算子的离散化过程中，需要同时采用上述两种处理方式。

最终，T_Ns及其子矩阵R_Ns、Σ_Ns、Γ_Ns的维数分别为(Ns–N+1)n×(Ns–N+1)n、sn×sn、sn×(Ns–N+1)n、(N–1)(s–1)n×(Ns–N+1)n。

Σ_Ns和Γ_Ns的表达式分别为：

式中：系数矩阵由与状态Δx(c_qh+h-τ_i),q＝1,…,s,对应的拉格朗日插值系数决定。和分别表示为：

Lobatto-ⅢB离散化方案:

Lobatto-ⅢB法的系数(A,b,c)有如下特点：

1)0＝c₁<…<c_s＝1；

2)b≠(a_s1,…,a_ss)^T。

应用Lobatto-ⅢB离散化方案时，将矩阵A的最后一行系数替换为b，剩下的公式推导与Lobatto-ⅢA/C方法一致。

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

采用旋转-放大方法进行预处理，一方面将时滞电力系统阻尼比小于ζ的特征值λ乘以旋转因子e^-jθ，ζ＝sin(θ)，从而将时滞系统阻尼比小于ζ的特征值λ变换到s平面虚轴右侧，进而利用(1.5)将这些特征值映射到z平面单位圆之外；另一方面，将z平面的特征值进行α次乘方，以增加这些特征值之间相对距离(即分散性)。

预处理后的特征值λ′与时滞电力系统原始特征值λ之间的关系为：

将(1.32)代入(1.3)，可得预处理后的特征方程：

式中：

为了保证预处理后时滞仍然为实数，式(1.35)对时滞做了近似处理。因此，实际计算得到的λ′是(1.33)解的近似值。

预处理后，T_Ns变为T_N′s：

式中：将R_Ns和Σ_Ns表达式中的替换为得到R_N′s和Σ_N′s，；

T_N′s的特征值μ′与λ′之间的谱映射关系为：

采用隐式重启动Arnoldi算法(implicitly restarted Arnoldi,IRA)来计算T_N′s模值最大的r个特征值μ′。在计算过程中，计算量最大的操作是迭代生成Krylov向量。

采用Radau-IIA离散化方案，阐述由第j个Krylov向量高效计算第j+1个Krylov向量q_j+1的方法。

q_j+1＝T_N′sq_j (1.38)

由于T_N′s的次对角为单位矩阵，q_j+1的后(N′–1)sn个元素通过直接将q_j的元素向后位移sn得到。前sn个元素q_j+1(1:sn)的计算分解为(1.39)和(1.40)。

z＝Σ_N′sq_j (1.39)

式中：为中间向量。

令q_j＝vec(Q)，其中vec(·)表示将矩阵按列重排为列向量的操作。

利用Kronecker积的性质，公式(1.41)给出了公式(1.39)的高效实现。

由于为高度稀疏矩阵，其计算量很小。

由于难以获得R_N′s的逆，式(1.40)只能通过迭代求解直至收敛。

其中，第l次迭代需要计算：

式中：其中

式(1.42)中涉及与的乘积运算是每次IRA迭代计算量最大之处，k＝1,…,s，通过利用公式(1.2)中系统增广状态矩阵的稀疏性(矩阵是稀疏的是指矩阵数值为0的元素数目远远多于非0元素的数目，并且非0元素分布没有规律)计算得到。

所述步骤(4)的步骤为：在频域，的特征值μ与时滞电力系统的特征值λ之间的关系可用式(1.43)表示。其中，将s平面虚轴右侧的λ映射到z平面单位圆之外。

式中：μ为的特征值；σ(·)表示谱；\表示集合差运算。

所述步骤(5)的步骤为：在计算得到T_N′s模值最大的r个特征值μ′后，通过式(1.37)的逆变换以及式(1.32)得到系统特征值的近似

特征向量的近似直接取为与对应Krylov向量的前n个元素。

以为初始值，利用稀疏实现的牛顿法求解公式(1.3)，得到时滞电力系统的准确特征值和特征向量{λ,v}。

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

若存在某个特征值λ的阻尼比小于0，则时滞电力系统处于小干扰不稳定的状态；

若所有特征值的最小阻尼比等于0，则时滞电力系统处于临界稳定的状态；

若所有特征值的阻尼比均大于0，则时滞电力系统处于渐进稳定的状态。

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

第一、本发明提出的SOD-IRK算法用于计算实际系统机电振荡模式对应的关键特征值时，充分考虑了实际系统的规模，以及通信时滞的影响。

第二、本发明提出的SOD-IRK算法只需计算解算子离散化近似矩阵的模值最大的设定个数特征值，就可以得到时滞电力系统机电振荡模式对应的特征值。

第三、本发明提出的SOD-IRK算法综合了多种隐式龙格-库塔离散化方案，拓展了大规模时滞电力系统关键特征值的计算方法，并且其得到的关键特征值有助于广域阻尼控制器的设计。

第四、本发明的方法是基于SOD-IRK降低了SOD-PS-R方法得到的离散化矩阵维数，大大提高了时滞系统关键特征值的计算效率，特别是规模较大的系统。

附图说明

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

图1为本发明的流程图。

图2为离散点集合Ω_Ns

具体实施方式

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

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

如图1所示，基于SOD-IRK的时滞电力系统稳定性判别方法，包括：

步骤(1)：建立时滞电力系统模型；依据时滞电力系统模型的特征值与时滞电力系统模型的解算子特征值之间的谱映射关系，将计算时滞电力系统模型的特征值转化成计算解算子的特征值；从而将判断时滞电力系统稳定性的问题转化为计算解算子的模值最大的特征值问题；

步骤(2)：采用不同的隐式龙格-库塔方法对解算子进行离散化，得到解算子的离散化矩阵；将无限维的特征值问题转化为有限维的特征值问题；

步骤(3)：采用隐式Arnoldi算法计算步骤(2)得到的解算子的离散化矩阵模值最大的特征值的近似值；计算过程中采用旋转-放大方法进行预处理，采用Kronecker积的性质进行稀疏实现；

步骤(4)：根据谱映射关系，将解算子离散化矩阵模值最大的特征值的近似值转化为时滞电力系统模型的近似特征值；

步骤(5)：采用牛顿迭代法对近似特征值进行修正，得到时滞电力系统的精确特征值；

步骤(6)：根据精确特征值阻尼比的大小来判断时滞电力系统的稳定性。

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

式中：为电力系统的状态变量向量，n为系统状态变量总数。t为当前时刻。0<τ₁<τ₂<…<τ_i…<τ_m为时滞环节的时滞常数，其中最大的时滞为τ_m。为系统状态矩阵，为稠密矩阵。为系统时滞状态矩阵，为稀疏矩阵。Δx(t)为t时刻系统状态变量的增量，Δx(t-τ_i)为t-τ_i时刻系统状态变量的增量，为t时刻系统状态变量导数的增量。Δx(0)为系统状态变量的初始值(即初始条件)，并简写为

式中：为高度稀疏的系统增广状态矩阵。

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

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

式(1.3)所示时滞电力系统的解算子X→X定义为将X上t时刻初值条件映射为t+h时刻系统状态的线性算子，即：

式中：是由区间[-τ_max,0]到n维实数空间映射的连续函数构成的巴拿赫空间，并赋有上确界范数。h为转移步长，0≤h≤τ_max。

在频域，与时滞电力系统的谱之间存在如下指数关系

式中：μ为的特征值；σ(·)表示谱；\表示集合差运算。

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

首先，将时滞区间[-τ_max,0]分为长度为h＝τ_max/N的N个子区间，[θ_N,θ_N–1],…,[θ₁,θ₀]；然后，利用p阶s级IRK的s个横坐标c₁,…,c_s对每个子区间进行离散化；最终，得到具有Ns个元素的集合Ω_Ns，如图2所示，

用Ω_Ns代替[-τ_max,0]，从而连续空间X被转化为离散空间令表示t＝δh+θ_j+c_qh(δ＝0,1；j＝1,…,N；q＝1,…,s)处Δx(t)的近似值，则X_Ns上系统的状态变量可表示为Δx^(δ)，即：

式中：和分别表示定义在[-τ_max,0]和[h-τ_max,h]上的系统状态变量。

结合解算子定义(1.4)和时滞系统状态方程(1.1)，可知Δx_h(t)是一个分段函数，其离散化原理如下：

(2-1)当t∈[-τ_max,-h]，Δx_h(t)为初始状态利用解算子的转移特性(1.4)，可以得到t＝-jh+c_qh处Δx(t+h)的估计值为公式(1.8)：其中，q＝1,…,s；j＝2,…,N；

将式(1.8)写成向量形式，得：

(2-2)当t∈[-h,0]，t+h-τ_i<0，i＝1,…,m，故Δx_h(t)是常微分方程(1.10)的解。

SOD-IRK的基本思想，就是得到公式(1.9)和公式(1.10)后，利用不同的IRK方法推导得到不同的Δx⁽¹⁾与Δx⁽⁰⁾之间的显式表达式，得到解算子离散化矩阵T_Ns。所述IRK方法，包括：Radau-IIA方法、Gauss-Legendre/SDIRK方法、Radau-IA方法、Lobatto-ⅢA/C方法或Lobatto-ⅢB方法。

p阶s级IRK递推公式的系数(A,b,c)可用Butcher表示：

系数(A,b,c)有如下特点：

1)0≤c₁<…<c_s≤1；

2)如果方法是自洽的consistent，则

采用不同的IRK离散化方案，得到的解算子离散化矩阵T_Ns是不同的。接下来，以Radau-IIA方法为例，详细推导解算子离散化矩阵T_Ns的表达式；

Radau-IIA法递推公式的系数(A,b,c)具有如下特点：

1)A可逆；

2)0<c₁<…<c_s＝1；

3)b＝(a_s1,…,a_ss)。

利用IRK求解(1.10)必须首先计算得到Δx(c_kh-τ_i)，i＝1,…,m-1；k＝1,…,s；考虑到点集c_kh-τ_i通常不属于离散集合Ω_N或h+Ω_N，利用c_kh-τ_i附近的p+1个点处Δx的近似解进行中间拉格朗日插值，从而得到Δx(c_kh-τ_i)的近似值

式中：l_r(c_kh-τ_i)表示由点c_kh-τ_i附近p+1个点计算得到的第r个拉格朗日插值系数。

其中，表示大于或等于p/2的最小整数。γ_i为正整数且满足不等式：

令

式中：为由拉格朗日插值系数决定的常数矩阵，表示Kronecker积。

将(1.12)代入(1.10)，进而利用基于Radau-IIA离散化方案的IRK方法可得t＝–h+c_qh(q＝1,…,s)处Δx_h(t)的估计值

将(1.15)改写成向量形式，得：

式中：

将R_Ns和Σ_Ns写为常数矩阵与系统状态矩阵的Kronecker积之和的形式：

式中：为常系数矩阵，

综合(1.9)和(1.16)，可得Δx⁽¹⁾与Δx⁽⁰⁾之间的关系式：

Δx⁽¹⁾＝T_NsΔx⁽⁰⁾ (1.19)

式中：X_Ns→X_Ns为解算子的IRK离散化矩阵，

Γ_Ns＝[I_(N-1)sn 0_(N-1)sn×sn] (1.21)

至此，无穷维解算子的特征值问题已被转化为有限维离散化矩阵T_Ns的特征值问题。

其他IRK离散化方案

解算子的其他IRK离散化方案与Radau-IIA离散化方案的不同之处，对应两种处理方式，总结为以下两点：

第一种处理方式：当c_s≠1时，增加离散点c_s+1＝1以得到每个子区间最右端点处系统状态的估计值。

第二种处理方式：当c₁＝0时，每个子区间左端点与相邻子区间的右端点重合，即θ_j–1+c₁h＝θ_j+c_sh或θ_j–1+c₁h＝θ_j+c_s+1h，j＝2,…,N。因而，需要去掉重复的离散点。

Gauss-Legendre/SDIRK离散化方案:

Gauss-Legendre/SDIRK离散化的系数c满足：0<c₁<…<c_s<1。

在解算子的离散化过程中，需要采用上述第一种处理方式。

最终，T_Ns及其子矩阵R_Ns、Σ_Ns、Γ_Ns的维数分别为N(s+1)n×N(s+1)n、(s+1)n×(s+1)n、(s+1)n×N(s+1)n、(N–1)(s+1)n×N(s+1)n。

它们的表达式分别为：

Γ_Ns＝[I_(N-1)(s+1)n 0_{(N-1)(s+1)n×(s+1)n}] (1.24)

式中：

Radau-IA离散化方案:

Radau-ⅠA的系数c满足：0＝c₁<…<c_s<1。

在解算子的离散化过程中，需要采用上述第二种处理方式。最终，T_Ns及其子矩阵R_Ns、Σ_Ns、Γ_Ns的维数分别为(Ns+1)n×(Ns+1)n、(s+1)n×(s+1)n、(s+1)n×(Ns+1)n、(N–1)sn×(Ns+1)n。

Σ_Ns和Γ_Ns的表达式分别为：

式中：

Lobatto-ⅢA/C离散化方案:

Lobatto-ⅢA/C离散化的系数(A,b,c)有如下特点：

1)0＝c₁<…<c_s＝1；

2)b＝(a_s1,…,a_ss)^T。

在解算子的离散化过程中，需要同时采用上述两种处理方式。

最终，T_Ns及其子矩阵R_Ns、Σ_Ns、Γ_Ns的维数分别为(Ns–N+1)n×(Ns–N+1)n、sn×sn、sn×(Ns–N+1)n、(N–1)(s–1)n×(Ns–N+1)n。

Σ_Ns和Γ_Ns的表达式分别为：

式中：系数矩阵由与状态Δx(c_qh+h-τ_i),q＝1,…,s,对应的拉格朗日插值系数决定。和分别表示为：

Lobatto-ⅢB离散化方案:

Lobatto-ⅢB法的系数(A,b,c)有如下特点：

1)0＝c₁<…<c_s＝1；

2)b≠(a_s1,…,a_ss)^T。

应用Lobatto-ⅢB离散化方案时，将矩阵A的最后一行系数替换为b，剩下的公式推导与Lobatto-ⅢA/C方法一致。

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

采用旋转-放大方法进行预处理，一方面将时滞电力系统阻尼比小于ζ的特征值λ乘以旋转因子e^-jθ，ζ＝sin(θ)，从而将时滞系统阻尼比小于ζ的特征值λ变换到s平面虚轴右侧，进而利用(1.5)将这些特征值映射到z平面单位圆之外；另一方面，将z平面的特征值进行α次乘方，以增加这些特征值之间相对距离(即分散性)。

预处理后的特征值λ′与时滞电力系统原始特征值λ之间的关系为：

将(1.32)代入(1.3)，可得预处理后的特征方程：

式中：

为了保证预处理后时滞仍然为实数，式(1.35)对时滞做了近似处理。因此，实际计算得到的λ′是(1.33)解的近似值。

预处理后，T_Ns变为T_N′s：

式中：将R_Ns和Σ_Ns表达式中的替换为得到R_N′s和Σ_N′s，；

T_N′s的特征值μ′与λ′之间的谱映射关系为：

采用隐式重启动Arnoldi算法(implicitly restarted Arnoldi,IRA)来计算T_N′s模值最大的r个特征值μ′。在计算过程中，计算量最大的操作是迭代生成Krylov向量。

采用Radau-IIA离散化方案，阐述由第j个Krylov向量高效计算第j+1个Krylov向量q_j+1的方法。

q_j+1＝T_N′sq_j (1.38)

由于T_N′s的次对角为单位矩阵，q_j+1的后(N′–1)sn个元素通过直接将q_j的元素向后位移sn得到。前sn个元素q_j+1(1:sn)的计算分解为(1.39)和(1.40)。

z＝Σ_N′sq_j (1.39)

式中：为中间向量。

令q_j＝vec(Q)，其中vec(·)表示将矩阵按列重排为列向量的操作。

利用Kronecker积的性质，公式(1.41)给出了公式(1.39)的高效实现。

由于为高度稀疏矩阵，其计算量很小。

由于难以获得R_N′s的逆，式(1.40)只能通过迭代求解直至收敛。

其中，第l次迭代需要计算：

式中：其中

式(1.42)中涉及与的乘积运算是每次IRA迭代计算量最大之处，k＝1,…,s，通过利用公式(1.2)中系统增广状态矩阵的稀疏性(矩阵是稀疏的是指矩阵数值为0的元素数目远远多于非0元素的数目，并且非0元素分布没有规律)计算得到。

所述步骤(4)的步骤为：在频域，的特征值μ与时滞电力系统的特征值λ之间的关系可用式(1.43)表示。其中，将s平面虚轴右侧的λ映射到z平面单位圆之外。

式中：μ为的特征值；σ(·)表示谱；\表示集合差运算。

所述步骤(5)的步骤为：在计算得到T_N′s模值最大的r个特征值μ′后，通过式(1.37)的逆变换以及式(1.32)得到系统特征值的近似

特征向量的近似直接取为与对应Krylov向量的前n个元素。

以为初始值，利用稀疏实现的牛顿法求解公式(1.3)，得到时滞电力系统的准确特征值和特征向量{λ,v}。

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

若存在某个特征值λ的阻尼比小于0，则时滞电力系统处于小干扰不稳定的状态；

若所有特征值的最小阻尼比等于0，则时滞电力系统处于临界稳定的状态；

若所有特征值的阻尼比均大于0，则时滞电力系统处于渐进稳定的状态。

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

Claims (10)

1.基于SOD-IRK的时滞电力系统稳定性判别方法，其特征是，包括：

步骤(1)：建立时滞电力系统模型；依据时滞电力系统模型的特征值与时滞电力系统模型的解算子特征值之间的谱映射关系，将计算时滞电力系统模型的特征值转化成计算解算子的特征值；从而将判断时滞电力系统稳定性的问题转化为计算解算子的模值最大的特征值问题；

步骤(2)：采用不同的隐式龙格-库塔方法对解算子进行离散化，得到解算子的离散化矩阵；将无限维的特征值问题转化为有限维的特征值问题；

步骤(3)：采用隐式Arnoldi算法计算步骤(2)得到的解算子的离散化矩阵模值最大的特征值的近似值；计算过程中采用旋转-放大方法进行预处理，采用Kronecker积的性质进行稀疏实现；

步骤(4)：根据谱映射关系，将解算子离散化矩阵模值最大的特征值的近似值转化为时滞电力系统模型的近似特征值；

步骤(5)：采用牛顿迭代法对近似特征值进行修正，得到时滞电力系统的精确特征值；

步骤(6)：根据精确特征值阻尼比的大小来判断时滞电力系统的稳定性。

2.如权利要求1所述的基于SOD-IRK的时滞电力系统稳定性判别方法，其特征是，所述步骤(1)的时滞电力系统模型如下：

式中：为电力系统的状态变量向量，n为系统状态变量总数；t为当前时刻；0<τ₁<τ₂<…<τ_i…<τ_m为时滞环节的时滞常数，其中最大的时滞为τ_m；为系统状态矩阵，为稠密矩阵；为系统时滞状态矩阵，为稀疏矩阵；Δx(t)为t时刻系统状态变量的增量，Δx(t-τ_i)为t-τ_i时刻系统状态变量的增量，为t时刻系统状态变量导数的增量；

Δx(0)为系统状态变量的初始值，并简写为

式中：为高度稀疏的系统增广状态矩阵；

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

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

式(1.3)所示时滞电力系统的解算子X→X定义为将X上t时刻初值条件映射为t+h时刻系统状态的线性算子，即：

式中：是由区间[-τ_max,0]到n维实数空间映射的连续函数构成的巴拿赫空间，并赋有上确界范数；h为转移步长，0≤h≤τ_max；

在频域，与时滞电力系统的谱之间存在如下指数关系

式中：μ为的特征值；σ(·)表示谱；\表示集合差运算。

3.如权利要求2所述的基于SOD-IRK的时滞电力系统稳定性判别方法，其特征是，所述步骤(2)的步骤为：

首先，将时滞区间[-τ_max,0]分为长度为h＝τ_max/N的N个子区间，[θ_N,θ_N–1],…,[θ₁,θ₀]；然后，利用p阶s级IRK的s个横坐标c₁,…,c_s对每个子区间进行离散化；最终，得到具有Ns个元素的集合Ω_Ns：

用Ω_Ns代替[-τ_max,0]，从而连续空间X被转化为离散空间令表示t＝δh+θ_j+c_qh(δ＝0,1；j＝1,…,N；q＝1,…,s)处Δx(t)的近似值，则X_Ns上系统的状态变量可表示为Δx(^δ)，即：

式中：和分别表示定义在[-τ_max,0]和[h-τ_max,h]上的系统状态变量；j＝1，…，N，δ＝0，1；

结合解算子定义(1.4)和时滞系统状态方程(1.1)，可知Δx_h(t)是一个分段函数，其离散化原理如下：

(2-1)当t∈[-τ_max,-h]，Δx_h(t)为初始状态利用解算子的转移特性(1.4)，得到t＝-jh+c_qh处Δx(t+h)的估计值为公式(1.8)：其中，q＝1,…,s；j＝2,…,N；

将式(1.8)写成向量形式，得：

(2-2)当t∈[-h,0]，t+h-τ_i<0，i＝1,…,m，故Δx_h(t)是常微分方程(1.10)的解；

得到公式(1.9)和公式(1.10)后，利用IRK方法推导得到解算子离散化矩阵T_Ns；

所述IRK方法，包括：Radau-IIA方法、Gauss-Legendre/SDIRK方法、Radau-IA方法、Lobatto-ⅢA/C方法或Lobatto-ⅢB方法。

4.如权利要求3所述的基于SOD-IRK的时滞电力系统稳定性判别方法，其特征是，

p阶s级IRK递推公式的系数(A,b,c)可用Butcher表示：

系数(A,b,c)有如下特点：

1)0≤c₁<…<c_s≤1；

2)如果方法是自洽的consistent，则i＝1,…,s；

假设IRK方法是Radau-IIA方法，解算子离散化矩阵T_Ns的表达式如下：

Radau-IIA法递推公式的系数(A,b,c)具有如下特点：

1)A可逆；

2)0<c₁<…<c_s＝1；

3)b＝(a_s1,…,a_ss)；

利用IRK求解(1.10)必须首先计算得到Δx(c_kh-τ_i)，i＝1,…,m-1；k＝1,…,s；考虑到点集c_kh-τ_i通常不属于离散集合Ω_N或h+Ω_N，利用c_kh-τ_i附近的p+1个点处Δx的近似解进行中间拉格朗日插值，从而得到Δx(c_kh-τ_i)的近似值

式中：l_r(c_kh-τ_i)表示由点c_kh-τ_i附近p+1个点计算得到的第r个拉格朗日插值系数；

其中，表示大于或等于p/2的最小整数；γ_i为正整数且满足不等式：

令i＝1,…,m-1：

式中：为由拉格朗日插值系数决定的常数矩阵，i＝1,…,m-1；表示Kronecker积；

将(1.12)代入(1.10)，进而利用基于Radau-IIA离散化方案的IRK方法可得t＝–h+c_qh处Δx_h(t)的估计值q＝1,…,s；

将(1.15)改写成向量形式，得：

式中：

将R_Ns和Σ_Ns写为常数矩阵与系统状态矩阵的Kronecker积之和的形式：

式中：为常系数矩阵，

综合(1.9)和(1.16)，得Δx⁽¹⁾与Δx⁽⁰⁾之间的关系式：

Δx⁽¹⁾＝T_NsΔx⁽⁰⁾ (1.19)

式中：X_Ns→X_Ns为解算子的IRK离散化矩阵，

Γ_Ns＝[I_(N-1)sn 0_(N-1)sn×sn] (1.21)

至此，无穷维解算子的特征值问题已被转化为有限维离散化矩阵T_Ns的特征值问题。

5.如权利要求3所述的基于SOD-IRK的时滞电力系统稳定性判别方法，其特征是，

解算子的其他IRK离散化方案与Radau-IIA离散化方案的不同之处，对应两种处理方式，总结为以下两点：

第一种处理方式：当c_s≠1时，增加离散点c_s+1＝1以得到每个子区间最右端点处系统状态的估计值；

第二种处理方式：当c₁＝0时，每个子区间左端点与相邻子区间的右端点重合，即θ_j–1+c₁h＝θ_j+c_sh或θ_j–1+c₁h＝θ_j+c_s+1h，j＝2,…,N；因而，需要去掉重复的离散点。

6.如权利要求5所述的基于SOD-IRK的时滞电力系统稳定性判别方法，其特征是，

Gauss-Legendre/SDIRK离散化方案:

Gauss-Legendre/SDIRK离散化的系数c满足：0<c₁<…<c_s<1；

在解算子的离散化过程中，需要采用上述第一种处理方式；

最终，T_Ns及其子矩阵R_Ns、Σ_Ns、Γ_Ns的维数分别为N(s+1)n×N(s+1)n、(s+1)n×(s+1)n、(s+1)n×N(s+1)n、(N–1)(s+1)n×N(s+1)n；

表达式分别为：

Γ_Ns＝[I_(N-1)(s+1)n 0_{(N-1)(s+1)n×(s+1)n}] (1.24)

式中：＝0,…,m；

Radau-IA离散化方案:

Radau-ⅠA的系数c满足：0＝c₁<…<c_s<1；

在解算子的离散化过程中，需要采用上述第二种处理方式；最终，T_Ns及其子矩阵R_Ns、Σ_Ns、Γ_Ns的维数分别为(Ns+1)n×(Ns+1)n、(s+1)n×(s+1)n、(s+1)n×(Ns+1)n、(N–1)sn×(Ns+1)n；

Σ_Ns和Γ_Ns的表达式分别为：

式中：i＝1,…,m；

Lobatto-ⅢA/C离散化方案:

Lobatto-ⅢA/C离散化的系数(A,b,c)有如下特点：

1)0＝c₁<…<c_s＝1；

2)b＝(a_s1,…,a_ss)^T；

在解算子的离散化过程中，需要同时采用上述两种处理方式；

最终，T_Ns及其子矩阵R_Ns、Σ_Ns、Γ_Ns的维数分别为(Ns–N+1)n×(Ns–N+1)n、sn×sn、sn×(Ns–N+1)n、(N–1)(s–1)n×(Ns–N+1)n；

Σ_Ns和Γ_Ns的表达式分别为：

式中：＝1,…,m；系数矩阵＝1,…,m,由与状态Δx(c_qh+h-τ_i),q＝1,…,s,对应的拉格朗日插值系数决定；和分别表示为：

Lobatto-ⅢB离散化方案:

Lobatto-ⅢB法的系数(A,b,c)有如下特点：

1)0＝c₁<…<c_s＝1；

2)b≠(a_s1,…,a_ss)^T；

应用Lobatto-ⅢB离散化方案时，将矩阵A的最后一行系数替换为b，剩下的公式推导与Lobatto-ⅢA/C方法一致。

7.如权利要求4所述的基于SOD-IRK的时滞电力系统稳定性判别方法，其特征是，

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

采用旋转-放大方法进行预处理，一方面将时滞电力系统阻尼比小于ζ的特征值λ乘以旋转因子e^-jθ，ζ＝sin(θ)，从而将时滞系统阻尼比小于ζ的特征值λ变换到s平面虚轴右侧，进而利用(1.5)将这些特征值映射到z平面单位圆之外；另一方面，将z平面的特征值进行α次乘方，以增加这些特征值之间相对距离；

预处理后的特征值λ′与时滞电力系统原始特征值λ之间的关系为：

将(1.32)代入(1.3)，可得预处理后的特征方程：

式中：

为了保证预处理后时滞仍然为实数，式(1.35)对时滞做了近似处理；因此，实际计算得到的λ′是(1.33)解的近似值；

预处理后，T_Ns变为T_N′s：

式中：将R_Ns和Σ_Ns表达式中的替换为i＝1,…,m，得到R_N′s和Σ_N′s，；

T_N′s的特征值μ′与λ′之间的谱映射关系为：

采用隐式重启动Arnoldi算法(implicitly restarted Arnoldi,IRA)来计算T_N′s模值最大的r个特征值μ′；在计算过程中，计算量最大的操作是迭代生成Krylov向量；

采用Radau-IIA离散化方案，阐述由第j个Krylov向量高效计算第j+1个Krylov向量q_j+1的方法；

q_j+1＝T_N′sq_j (1.38)

由于T_N′s的次对角为单位矩阵，q_j+1的后(N′–1)sn个元素通过直接将q_j的元素向后位移sn得到；前sn个元素q_j+1(1:sn)的计算分解为(1.39)和(1.40)；

z＝Σ_N′sq_j (1.39)

式中：为中间向量；

令q_j＝vec(Q)，其中vec(·)表示将矩阵按列重排为列向量的操作；

利用Kronecker积的性质，公式(1.41)给出了公式(1.39)的高效实现；

由于为高度稀疏矩阵，其计算量小；

由于难以获得R_N′s的逆，式(1.40)只能通过迭代求解直至收敛；

其中，第l次迭代需要计算：

式中：其中

式(142)中涉及与的乘积运算是每次IRA迭代计算量最大之处，k＝1,…,s，通过利用公式(1.2)中系统增广状态矩阵的稀疏性计算得到。

8.如权利要求7所述的基于SOD-IRK的时滞电力系统稳定性判别方法，其特征是，

所述步骤(4)的步骤为：在频域，的特征值μ与时滞电力系统的特征值λ之间的关系可用式(1.43)表示；其中，将s平面虚轴右侧的λ映射到z平面单位圆之外；

式中：μ为的特征值；σ(·)表示谱；\表示集合差运算。

9.如权利要求8所述的基于SOD-IRK的时滞电力系统稳定性判别方法，其特征是，

所述步骤(5)的步骤为：在计算得到T_N′s模值最大的r个特征值μ′后，通过式(1.37)的逆变换以及式(1.32)得到系统特征值的近似

特征向量的近似直接取为与对应Krylov向量的前n个元素；

以为初始值，利用稀疏实现的牛顿法求解公式(1.3)，得到时滞电力系统的准确特征值和特征向量{λ,v}。

10.如权利要求9所述的基于SOD-IRK的时滞电力系统稳定性判别方法，其特征是，

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

若存在某个特征值λ的阻尼比小于0，则时滞电力系统处于小干扰不稳定的状态；

若所有特征值的最小阻尼比等于0，则时滞电力系统处于临界稳定的状态；

若所有特征值的阻尼比均大于0，则时滞电力系统处于渐进稳定的状态。