CN103838965B - 基于广义特征值的时滞稳定上限计算系统及其计算方法 - Google Patents

Abstract

本发明公开了电力系统稳定分析技术领域中的一种基于广义特征值的时滞稳定上限计算系统及其计算方法。系统包括顺序相连的数据采集模块、时滞系统处理模块、时滞上限求解模块和结果输出模块；方法包括：采集建立时滞系统状态方程所需的网络结构参数、发电机频率和发电机功角；建立时滞系统状态方程，并对时滞系统状态方程中的参数矩阵进行降阶处理，得到降阶后的时滞系统状态方程；生成基于改进自由权矩阵的时滞稳定判据；利用时滞稳定判据求解时滞稳定上限。本发明能够有效降低时滞稳定上限计算过程中的保守性，具有很好的正确性和有效性。

Description

基于广义特征值的时滞稳定上限计算系统及其计算方法

技术领域

本发明属于电力系统稳定分析技术领域，尤其涉及一种基于广义特征值的

时滞稳定上限计算系统及其计算方法。

背景技术

基于广域信息的反馈控制器的设计存在时滞问题，而这必然会造成控制器的控制效果降低，甚至出现负阻尼的情况。因此，迫切需要对系统的时滞稳定上限进行研究。

目前国内外学者在系统时滞稳定上限的研究方面取得了大量有益成果，主要可以分为3类：时域法、频域法和直接法。时域法可判定系统在特定场景下是否稳定，但在稳定程度、时滞稳定上限等信息的获取方面还需要进一步研究。频域法通过在实数空间中搜索时滞系统的关键特征根，能够在一定程度上揭示时滞系统变化规律，但计算量较大，求解速度有待提升。直接法借助Lyapunov理论和线性矩阵不等式(Linear Matrix Inequality,LMI)技术，可同时考虑时滞的随机变动、存在切换环节等情况，适用范围更为广泛，但该方法具有一定的保守性。

针对以上问题，本发明提出一种基于广义特征值的时滞稳定上限计算系统及其计算方法。本发明从电力系统实际出发，能够有效解决直接法求解时滞上限过程中保守性较高的问题。首先利用读入的广域信号数据建立时滞系统状态方程，并在有效保留系统中低频振荡成分的基础上降阶系统。其次，该发明形成基于改进自由权矩阵的时滞稳定判据，在此基础上，将时滞稳定判据等价变换，利用广义特征值法求解系统的时滞稳定上限。基于IEEE4机11节点系统和IEEE16机68节点系统仿真表明，本发明可以较好的降低传统方法的保守性，具有很好的有效性和正确性。

发明内容

本发明的目的在于，提供一种基于广义特征值的时滞稳定上限计算系统及其计算方法，用于解决直接法求解时滞上限过程中保守性较高的问题。

为了实现上述目的，本发明提出的技术方案是，一种基于广义特征值的时滞稳定上限计算系统，其特征是所述系统包括顺序相连的数据采集模块、时滞系统处理模块、时滞上限求解模块和结果输出模块；

所述数据采集模块用于采集建立时滞系统状态方程所需的网络结构参数、发电机频率和发电机功角，并将采集的数据发送至时滞系统处理模块；

所述时滞系统处理模块用于建立时滞系统状态方程，并对时滞系统状态方程中的参数矩阵进行降阶处理；

所述时滞上限求解模块用于生成基于改进自由权矩阵的时滞稳定判据，并利用时滞稳定判据求解时滞稳定上限；

所述结果输出模块用于输出时滞稳定上限结果。

一种基于广义特征值的时滞稳定上限计算方法，其特征是所述方法包括：

步骤1：采集建立时滞系统状态方程所需的网络结构参数、发电机频率和发电机功角；

步骤2：建立时滞系统状态方程，并对时滞系统状态方程中的参数矩阵进行降阶处理，得到降阶后的时滞系统状态方程；

步骤3：生成基于改进自由权矩阵的时滞稳定判据；

步骤4：利用时滞稳定判据求解时滞稳定上限。

所述降阶后的时滞系统状态方程为其中，x(t)为降阶后的电力系统状态向量；

Α为降阶后的电力系统状态矩阵；

Α_d为降阶后的电力系统时滞矩阵；

为降阶后的电力系统状态量对应的状态值；

d(t)为时滞，0≤d(t)≤h且

h为时滞稳定上限且h>0；

μ为时滞最大变化率。

所述基于改进自由权矩阵的时滞稳定判据为：

$\left[\begin{array}{ccccc}\mathrm{\&Phi;}& \mathrm{hN}& \mathrm{hS}& \mathrm{hM}& {\mathrm{hA}}_c^T(Z_1+Z_2)\\ {\mathrm{hN}}^T& {-\mathrm{hZ}}_1& 0& 0& 0\\ {\mathrm{hS}}^T& 0& {-\mathrm{hZ}}_1& 0& 0\\ {\mathrm{hM}}^T& 0& 0& {-\mathrm{hZ}}_2& 0\\ h(Z_1^T+Z_2^T)A_c& 0& 0& 0& -h(Z_1+Z_2)\end{array}\right]<0;$

其中， $\mathrm{\&Phi;}={\mathrm{\&Phi;}}_1{\mathrm{\&Phi;}}_2{\mathrm{\&Phi;}}_2^T;$

${\mathrm{\&Phi;}}_1=\left[\begin{array}{ccc}\mathrm{PA}+A^{T}P+Q+R& {\mathrm{PA}}_d& 0\\ A_d^{T}P& -(1-\mathrm{\&mu;})Q& 0\\ 0& 0& -R\end{array}\right];$

Φ₂=[N+M-N+S-M-S]；

A_c=[AA_d0]；

N、M和S为改进自由权矩阵；

改进自由权矩阵N满足 ${2\mathrm{\&zeta;}}_1^T\left(t\right)N[x\left(t\right)-x(t-d\left(t\right))-{\&Integral;}_{t-d\left(t\right)}^t\stackrel{\&CenterDot;}{x}\left(s\right)\mathrm{ds}]=0;$

改进自由权矩阵S满足 ${2\mathrm{\&zeta;}}_1^T\left(t\right)S[x(t-d\left(t\right))-x(t-h)-{\&Integral;}_{t-h}^{t-d\left(t\right)}\stackrel{\&CenterDot;}{x}\left(s\right)\mathrm{ds}]=0;$

改进自由权矩阵M满足 ${2\mathrm{\&zeta;}}_1^T\left(t\right)M[x\left(t\right)-x(t-h)-{\&Integral;}_{t-h}^t\stackrel{\&CenterDot;}{x}\left(s\right)\mathrm{ds}]=0;$

?₁(t)=[x^T(t)x^T(t-d(t))x^T(t-h)]^T；

P、Q、R、Z₁和Z₂为待定矩阵。

所述利用时滞稳定判据求解时滞稳定上限包括：

子步骤A1：将时滞稳定判据等价变换为

$\left[\begin{array}{ccccc}\mathrm{\&Phi;}& N& S& M& A_c^T(Z_1+Z_2)\\ N^T& {-Y}_1& 0& 0& 0\\ S^T& 0& {-Y}_1& 0& 0\\ M^T& 0& 0& {-Y}_2& 0\\ (Z_1^T+Z_2^T)A_c& 0& 0& 0& -(Y_1+Y_2)\end{array}\right]<0;$

其中，Y₁和Y₂为附加矩阵，并且Y₁=Y₁ ^T≥0，Y₂=Y₂ ^T≥0，

$\left[\begin{array}{cc}{Y}_1& 0\\ 0& Y_2\end{array}\right]<v\left[\begin{array}{cc}{Z}_1& 0\\ 0& Z_2\end{array}\right],v=1/h;$

子步骤A2：以v最小为目标，以时滞稳定判据等价变换后的矩阵不等式 $\left[\begin{array}{ccccc}\mathrm{\&Phi;}& N& S& M& A_c^T(Z_1+Z_2)\\ N^T& {-Y}_1& 0& 0& 0\\ S^T& 0& {-Y}_1& 0& 0\\ M^T& 0& 0& {-Y}_2& 0\\ (Z_1^T+Z_2^T)A_c& 0& 0& 0& -(Y_1+Y_2)\end{array}\right]<0$ 以及 $\left[\begin{array}{cc}{Y}_1& 0\\ 0& Y_2\end{array}\right]<\left[\begin{array}{cc}{Z}_1& 0\\ 0& Z_2\end{array}\right]$ 为约束条件，计算待定矩阵P、Q、R、Z₁、Z₂、Y₁和Y₂，进而得到时滞稳定上限h。

本发明利用改进自由权矩阵建立时滞系统稳定判据，借助广义特征值法对系统时滞上限进行求解，其能够有效降低时滞稳定上限计算过程中的保守性，具有很好的正确性和有效性。

附图说明

图1是基于广义特征值的时滞稳定上限计算系统结构图；

图2是IEEE4机11节点系统结构图；

图3是IEEE4机11节点SMA降阶前后特征根对比图；其中，(a)是全阶开环状态矩阵和降阶开环状态矩阵的特征根对比图，(b)是全阶闭环状态矩阵和降阶闭环状态矩阵的特征根对比图；

图4是IEEE4机11节点不同延时情况下发电机1-4相对功角动态响应曲线图；

图5是IEEE4机11节点不同延时情况下发电机2-3相对功角动态响应曲线图；

图6是IEEE4机11节点系统G1与G4功角差各时滞时间下的阻尼比结果表；

图7是IEEE4机11节点系统G2与G3功角差各时滞时间下的阻尼比结果表；

图8是IEEE16机68节点系统结构图；

图9是IEEE16机68节点Schur降阶前后频率响应对比图；其中，(a)是全阶开环状态矩阵和降阶开环状态矩阵的特征根对比图，(b)是全阶闭环状态矩阵和降阶闭环状态矩阵的特征根对比图；

图10是IEEE16机68节点不同延时情况下发电机1-16相对功角动态响应曲线图；

图11是IEEE16机68节点不同延时情况下发电机3-14相对功角动态响应曲线图；

图12是IEEE16机68节点不同延时情况下发电机10-15相对功角动态响应曲线图；

图13是16机系统G1与G16功角差各时滞时间下的阻尼比表；

图14是16机系统G3与G14功角差各时滞时间下的阻尼比表；

图15是16机系统G10与G15功角差各时滞时间下的阻尼比表。

具体实施方式

下面结合附图，对优选实施例作详细说明。应该强调的是，下述说明仅仅是示例性的，而不是为了限制本发明的范围及其应用。

实施例1

图1是本发明提供的基于广义特征值的时滞稳定上限计算系统结构图。如图1所示，本发明提供的基于广义特征值的时滞稳定上限计算系统结构图包括顺序相连的数据采集模块、时滞系统处理模块、时滞上限求解模块和结果输出模块。

数据采集模块用于采集网络结构参数、发电机频率和发电机功角，并将采集的数据发送至时滞系统处理模块。其中，网络结构参数是用于建立时滞系统状态方程所需的参数，而发电机频率和发电机功角则是时滞系统状态方程中的状态向量。

时滞系统处理模块用于建立时滞系统状态方程，并对时滞系统状态方程中的参数矩阵进行降阶处理，得到降阶后的时滞系统状态方程。

时滞上限求解模块用于生成基于改进自由权矩阵的时滞稳定判据；并对所述时滞稳定判据进行等价变换，利用广义特征值法求解系统的时滞稳定上限。

结果输出模块用于输出时滞稳定上限结果。

本发明提供的基于广义特征值的时滞稳定上限计算方法包括：

步骤1：采集建立时滞系统状态方程所需的网络结构参数、发电机频率和发电机功角。

网络结构参数包括时滞系统中线路的阻抗、导纳、发电机的内阻抗和负荷的等值阻抗。发电机频率用于计算时滞系统转速，发电机功角和时滞系统转速的变化量为时滞系统状态方程的状态向量。

步骤2：建立时滞系统状态方程，并对时滞系统状态方程中的参数矩阵进行降阶处理，得到降阶后的时滞系统状态方程。

多输入多输出电力系统的状态方程可表示为：

$\left\{\begin{array}{c}{\stackrel{\&CenterDot;}{x}}^{\&prime;}\left(t\right)=A^{\&prime;}{x}^{\&prime;}\left(t\right)+B^{\&prime;}{u}^{\&prime;}\left(t\right)\\ u^{\&prime;}\left(t\right)=K_1^{\&prime;}{x}^{\&prime;}\left(t\right)\end{array}\right.---\left(1\right)$

其中，x′(t)∈Rⁿ为电力系统状态向量，u′(t)∈R^m为电力系统控制输入向量，A′∈R^n×n为电力系统状态矩阵，B′∈R^n×m为电力系统控制矩阵。

通过状态反馈后得到相应的闭环系统为：

${\stackrel{\&CenterDot;}{x}}^{\&prime;}\left(t\right)=C^{\&prime;}{x}^{\&prime;}\left(t\right)---\left(2\right)$

其中，C′为闭环状态矩阵，当系统经过状态反馈时，闭环状态矩阵为C′=A′+B′K′₁，其中，K′₁∈R^m×n为各附加控制器的综合状态反馈矩阵。

在实际电力系统中，控制输入向量通过SCADA/WAMS系统向各控制器传达，信号传递过程必然存在一定时滞，则相应闭环系统可描述为：

${\stackrel{\&CenterDot;}{x}}^{\&prime;}\left(t\right)=A^{\&prime;}{x}^{\&prime;}\left(t\right)+B^{\&prime;}{K}_1^{\&prime;}{x}^{\&prime;}(t-d\left(t\right))---\left(3\right)$

由式(3)可知，电力系统时滞矩阵A′_d=B′K′₁。对于含时滞环节的系统，其状态方程有如下形式：

其中，矩阵A′和A′_d分别为电力系统状态矩阵和电力系统时滞矩阵，h为时滞稳定上限。公式(4)中时滞d(t)满足条件：

0≤d(t)≤h (5)

$\stackrel{\&CenterDot;}{d}\left(t\right)\&le;\mathrm{\&mu;}---\left(6\right)$

公式(4)-式(6)即为时滞系统状态方程，μ为时滞最大变化率。

实际系统分析中，仅对低频振荡特征根相应的特征向量中与系统状态量Δω（发电机转速变化量）和Δδ（发电机功角变化量）相对应的元素感兴趣，以便了解在Δω_i中所含的该振荡模式分量的相对幅值及相位。常用的系统降阶方法包括SMA降阶方法和Schur平衡降阶方法。由于系统降阶方法已经是本领域常用的方法，因此本发明只以选择模式分析法SMA为例，对系统降阶做简单介绍。

对于系统状态方程将其按下式划分

$\left[\begin{array}{c}{\stackrel{\&CenterDot;}{X}}_1\\ {\stackrel{\&CenterDot;}{X}}_2\end{array}\right]=\left[\begin{array}{cc}{A}_{11}& A_{12}\\ A_{21}& A_{22}\end{array}\right]\left[\begin{array}{c}{X}_1\\ X_2\end{array}\right]---\left(7\right)$

其中，X₁=[Δω^T,Δδ^T]为保留变量，X₂为其他变量，待消去。

由式(7)可消去X₂，得：

${\stackrel{\&CenterDot;}{X}}_1=[A_{11}+A_{12}{(\mathrm{pI}-A_{22})}^{-1}{A}_{21}]X_1---\left(8\right)$

其中，I为单位阵，p为微分算子。

将上式改写为

${\stackrel{\&CenterDot;}{X}}_r=A_r\left(p\right)X_r---\left(9\right)$

其中，X_r=X₁为保留变量，A_r(p)为运算形式的“降阶”系统系数阵。

由式(7)-(9)可以得到两个重要性质：

(1)如果p=λ₁(i=1,2,…,N)为式(7)相应系统特征根，即|λ_iI-A|=0,则_p=λ_i也为式(8)或(9)形式上降阶的系统特征根，即亦有|λ_iI-A_r(λ_i)|=0，特征根不发生变化，系统模式不变。

(2)对于原系统，λ_i的特征向量u_i，有Au_i=λ_iu_i。设降阶系统λ_i相应的特征向量为u_ri，即A_r(λ_i)u_ri=λ_iu_ri，则u_ri和u_i中保留变量X_r相对应的元素相等，即特征向量的相应元素不变。因此，在X_r保留变量处去观察同一模式λ_i的振荡时，相对幅值即相位不变，或者说模态不变。这样，所关心频带输入输出特性被完整的保留下来。

通过降阶处理，可以得到降阶后的时滞系统状态方程为：

公式（10）中，x(t)是降阶后的电力系统状态向量，Α为降阶后的电力系统状态矩阵，Α_d为降阶后的电力系统时滞矩阵，为降阶后的电力系统状态量对应的状态值，d(t)、h和μ的含义同公式（4）且满足式（5）和式（6）。

步骤3：生成基于改进自由权矩阵的时滞稳定判据。

构造如下形式的Lyapunov-Krasovskii泛函（李雅普诺夫-克拉索夫斯基泛函）：

$\begin{array}{c}V\left(x\right)=x^T\left(t\right)\mathrm{Px}\left(t\right)+{\&Integral;}_{t-d\left(t\right)}^t{x}^T\left(s\right)\mathrm{Qx}\left(s\right)\mathrm{ds}\\ +{\&Integral;}_{t-h}^t{x}^T\left(s\right)\mathrm{Rx}\left(s\right)\mathrm{ds}+{\&Integral;}_{-h}^0{\&Integral;}_{t+\mathrm{\&theta;}}^t{\stackrel{\&CenterDot;}{x}}^T\left(s\right)(Z_1+Z_2)\stackrel{\&CenterDot;}{x}\left(x\right)\mathrm{dsd\&theta;}\end{array}---\left(11\right)$

公式（11）中，P=P^T>0，Q=Q^T≥0，R=R^T≥0和Z_i=Z_i ^T≥0(i=1,2)都是待定矩阵。

由Newton-Leibniz公式（牛顿-莱布尼兹公式）可知，对于改进自由权矩阵N、S和M，式(12)-式(14)成立：

${2\mathrm{\&zeta;}}_1^T\left(t\right)N[x\left(t\right)-x(t-d\left(t\right))-{\&Integral;}_{t-d\left(t\right)}^t\stackrel{\&CenterDot;}{x}\left(s\right)\mathrm{ds}]=0---\left(12\right)$

${2\mathrm{\&zeta;}}_1^T\left(t\right)S[x(t-d\left(t\right))-x(t-h)-{\&Integral;}_{t-h}^{t-d\left(t\right)}\stackrel{\&CenterDot;}{x}\left(s\right)\mathrm{ds}]=0---\left(13\right)$

${2\mathrm{\&zeta;}}_1^T\left(t\right)M[x\left(t\right)-x(t-h)-{\&Integral;}_{t-h}^t\stackrel{\&CenterDot;}{x}\left(s\right)\mathrm{ds}]=0---\left(14\right)$

其中，ζ₁(t)=[x^T(t)x^T(t-d(t))x^T(t-h)]^T。

计算式(11)中V(x)关于t的导数为：

$\begin{array}{c}\stackrel{\&CenterDot;}{V}\left(x\right)={2x}^T\left(t\right)P\stackrel{\&CenterDot;}{x}\left(t\right)+x^T\left(s\right)\mathrm{Qx}\left(s\right)-(1-\stackrel{\&CenterDot;}{d}\left(t\right))x^T(t-d\left(t\right))\mathrm{Qx}(t-d\left(t\right))\\ +x^T\left(t\right)\mathrm{Rx}\left(t\right)-x^T(t-h)\mathrm{Rx}(t-h)+h{\stackrel{\&CenterDot;}{x}}^T\left(t\right)(Z_1+Z_2)\stackrel{\&CenterDot;}{x}\left(t\right)\\ -{\&Integral;}_{t-h}^t{\stackrel{\&CenterDot;}{x}}^T\left(s\right)(Z_1+Z_2)\stackrel{\&CenterDot;}{x}\left(s\right)\mathrm{ds}\end{array}---\left(15\right)$

将式(10)中第一个式子及式(12)-式(14)代入式(15)，并加入必要的松散项可得：

$\begin{array}{c}\stackrel{\&CenterDot;}{V}\left(x\right)\&le;x^T\left(t\right)(\mathrm{PA}+A^{T}P)x\left(t\right)+x^T\left(s\right)(Q+R)x\left(s\right)\\ -(1-\mathrm{\&mu;})x^T(t-d\left(t\right))\mathrm{Qx}(t-d\left(t\right))-x^T(t-h)\mathrm{Rx}(t-h)\\ +h{\stackrel{\&CenterDot;}{x}}^T\left(t\right)(Z_1+Z_2)\stackrel{\&CenterDot;}{x}\left(t\right)-{\&Integral;}_{t-h}^t{\stackrel{\&CenterDot;}{x}}^T\left(s\right)Z_1{\stackrel{\&CenterDot;}{x}}^T\left(s\right)\mathrm{ds}\\ -{\&Integral;}_{t-h}^{t-d\left(t\right)}{\stackrel{\&CenterDot;}{x}}^T\left(s\right)Z_1{\stackrel{\&CenterDot;}{x}}^T\left(s\right)\mathrm{ds}-{\&Integral;}_{t-h}^t{\stackrel{\&CenterDot;}{x}}^T\left(s\right)Z_2{\stackrel{\&CenterDot;}{x}}^T\left(s\right)\mathrm{ds}\\ +{2\mathrm{\&zeta;}}_1^T\left(t\right)N[x\left(t\right)-x(t-d\left(t\right))-{\&Integral;}_{t-d\left(t\right)}^t\stackrel{\&CenterDot;}{x}\left(s\right)\mathrm{ds}]\\ +{2\mathrm{\&zeta;}}_1^T\left(t\right)S[x(t-d\left(t\right))-x(t-h)-{\&Integral;}_{t-d\left(t\right)}^t\stackrel{\&CenterDot;}{x}\left(s\right)\mathrm{ds}]\\ +{2\mathrm{\&zeta;}}_1^T\left(t\right)M[x\left(t\right)-x(t-h)-{\&Integral;}_{t-h}^t\stackrel{\&CenterDot;}{x}\left(s\right)\mathrm{ds}]\\ \&le;{\mathrm{\&zeta;}}_1^T\left(t\right)[\mathrm{\&Phi;}+{\mathrm{\&Phi;}}_s]{\mathrm{\&zeta;}}_1\left(t\right)-{\&Integral;}_{t-d\left(t\right)}^t{\mathrm{\&theta;}}_1\mathrm{ds}-{\&Integral;}_{t-h}^{t-d\left(t\right)}{\mathrm{\&theta;}}_2\mathrm{ds}-{\&Integral;}_{t-h}^t{\mathrm{\&theta;}}_3\mathrm{ds}\end{array}----\left(16\right)$

${\mathrm{\&Phi;}}_s={\mathrm{hA}}_c^T(Z_1+Z_2)A_c+{\mathrm{hNZ}}_1^{-1}{N}^T+{\mathrm{hSZ}}_1^{-1}{S}^T+{\mathrm{hMZ}}_1^{-1}{M}^T---\left(17\right)$

${\mathrm{\&theta;}}_1=[{\mathrm{\&zeta;}}_1^T\left(t\right)N+\stackrel{\&CenterDot;}{x}\left(t\right)Z_1]Z_1^{-1}[N^T{\mathrm{\&zeta;}}_1\left(t\right)+Z_1\stackrel{\&CenterDot;}{x}\left(t\right)]---\left(18\right)$

${\mathrm{\&theta;}}_2=[{\mathrm{\&zeta;}}_1^T\left(t\right)S+\stackrel{\&CenterDot;}{x}\left(t\right)Z_1]Z_1^{-1}[S^T{\mathrm{\&zeta;}}_1\left(t\right)+Z_1\stackrel{\&CenterDot;}{x}\left(t\right)]---\left(19\right)$

${\mathrm{\&theta;}}_3=[{\mathrm{\&zeta;}}_1^T\left(t\right)M+\stackrel{\&CenterDot;}{x}\left(t\right)Z_2]Z_2^{-1}[M^T{\mathrm{\&zeta;}}_1\left(t\right)+Z_1\stackrel{\&CenterDot;}{x}\left(t\right)]---\left(20\right)$

上述公式中， $\mathrm{\&Phi;}={\mathrm{\&Phi;}}_1+{\mathrm{\&Phi;}}_2+{\mathrm{\&Phi;}}_2^T.$

${\mathrm{\&Phi;}}_1=\left[\begin{array}{ccc}\mathrm{PA}+A^{T}P+Q+R& {\mathrm{PA}}_d& 0\\ A_d^{T}P& -(1-\mathrm{\&mu;})Q& 0\\ 0& 0& -R\end{array}\right],$ Φ₂=[N+M-N+S-M-S]。

A_c=[AA_d0]。

$N^T=\left[\begin{array}{ccc}{N}_1^T& N_2^T& N_3^T\end{array}\right],$ $S^T=\left[\begin{array}{ccc}{S}_1^T& S_2^T& S_3^T\end{array}\right],$ $M^T=\left[\begin{array}{ccc}{M}_1^T& M_2^T& M_3^T\end{array}\right].$

N₁、N₂和N₃为矩阵N的分块矩阵且N₁的维度与PA+A^TP+Q+R的维度相等，N₂的维度与的维度相等，N₃的维度与R的维度相等。

S₁、S₂和S₃为矩阵S的分块矩阵且S₁的维度与PA+A^TP+Q+R的维度相等，S₂的维度与的维度相等，S₃的维度与R的维度相等。

M₁、M₂和M₃为矩阵M的分块矩阵且M₁的维度与PA+A^TP+Q+R的维度相等，M₂的维度与的维度相等，M₃的维度与R的维度相等。

考虑到式(11)中Z_i=Z_i ^T≥0，(i=1,2)，因此公式(18)-式(20)中θ_i=θ_i ^T≥0，(i=1,2,3)。若Φ+Φ_s≤0，则式(16)中对于Φ+Φ_s，利用Schur补后可得式(10)表征的时滞系统稳定判据如下：

对于给定标量h>0和μ，若存在P=P^T>0，Q=Q^T≥0，R=R^T≥0和Z_i=Z_i ^T≥0(i=1,2)， $N^T=\left[\begin{array}{ccc}{N}_1^T& N_2^T& N_3^T\end{array}\right],$ $S^T=\left[\begin{array}{ccc}{S}_1^T& S_2^T& S_3^T\end{array}\right]$ 和 $M^T=\left[\begin{array}{ccc}{M}_1^T& M_2^T& M_3^T\end{array}\right],$ 使得如下线性矩阵不等式成立：

$\left[\begin{array}{ccccc}\mathrm{\&Phi;}& \mathrm{hN}& \mathrm{hS}& \mathrm{hM}& {\mathrm{hA}}_c^T(Z_1+Z_2)\\ {\mathrm{hN}}^T& {-\mathrm{hZ}}_1& 0& 0& 0\\ {\mathrm{hS}}^T& 0& {-\mathrm{hZ}}_1& 0& 0\\ {\mathrm{hM}}^T& 0& 0& {-\mathrm{hZ}}_2& 0\\ h(Z_1^T+Z_2^T)A_c& 0& 0& 0& -h(Z_1+Z_2)\end{array}\right]<0---\left(21\right)$

则对于同时满足条件式(5)和式(6)的时滞系统(10)渐进稳定。

步骤4：利用时滞稳定判据求解时滞稳定上限。

式(21)表征的线性矩阵不等式仅能判定系统是否稳定，而无法获取系统时滞稳定上限等信息。考虑到广义特征值法能够求解优化问题的全局最小值，因此，本发明提出利用广义特征值法计算系统的时滞稳定上限。由于式(21)不是标准的广义特征值形式，无法直接利用广义特征值法进行求解，因此，做如下变换，将式(21)同时左乘与右乘式(22)，如下：

$\left[\begin{array}{ccccc}I& 0& 0& 0& 0\\ 0& I/h& 0& 0& 0\\ 0& 0& I/h& 0& 0\\ 0& 0& 0& I/h& 0\\ 0& 0& 0& 0& I/h\end{array}\right]>0---\left(22\right)$

可得：

$\left[\begin{array}{ccccc}\mathrm{\&Phi;}& N& S& M& A_c^T(Z_1+Z_2)\\ N^T& {-\mathrm{vZ}}_1& 0& 0& 0\\ S^T& 0& {-\mathrm{vZ}}_1& 0& 0\\ M^T& 0& 0& {-\mathrm{vZ}}_2& 0\\ (Z_1^T+Z_2^T)A_c& 0& 0& 0& -v(Z_1+Z_2)\end{array}\right]<0---\left(23\right)$

其中，v=1h，I为单位阵。

为了求解最大时滞稳定上限h，即最小的v，本发明引入附加矩阵Y_i=Y_i ^T≥0(i=1,2)，该矩阵需满足式(24)成立：

$\left[\begin{array}{cc}{Y}_1& 0\\ 0& Y_2\end{array}\right]<v\left[\begin{array}{cc}{Z}_1& 0\\ 0& Z_2\end{array}\right]---\left(24\right)$

再将式(24)代入式(23)可得：

$\left[\begin{array}{ccccc}\mathrm{\&Phi;}& N& S& M& A_c^T(Z_1+Z_2)\\ N^T& {-Y}_1& 0& 0& 0\\ S^T& 0& {-Y}_1& 0& 0\\ M^T& 0& 0& {-Y}_2& 0\\ (Z_1^T+Z_2^T)A_c& 0& 0& 0& -(Y_1+Y_2)\end{array}\right]<0---\left(25\right)$

由此，时滞稳定上限问题已经转化为以下优化问题：

以v最小为目标，以 $\left[\begin{array}{ccccc}\mathrm{\&Phi;}& N& S& M& A_c^T(Z_1+Z_2)\\ N^T& {-Y}_1& 0& 0& 0\\ S^T& 0& {-Y}_1& 0& 0\\ M^T& 0& 0& {-Y}_2& 0\\ (Z_1^T+Z_2^T)A_c& 0& 0& 0& -(Y_1+Y_2)\end{array}\right]<0$ 和 $\left[\begin{array}{cc}{Y}_1& 0\\ 0& Y_2\end{array}\right]<v\left[\begin{array}{cc}{Z}_1& 0\\ 0& Z_2\end{array}\right]$ 为约束条件，计算待定矩阵P、Q、R、Z₁、Z₂、Y₁和Y₂，进而得到时滞稳定上限h。公式表示如下：

min v

P，Q，R，Z_i，Y_i

s.t.

$\begin{array}{c}\left[\begin{array}{cc}{Y}_1& 0\\ 0& Y_2\end{array}\right]<v\left[\begin{array}{cc}{Z}_1& 0\\ 0& Z_2\end{array}\right]\\ \left[\begin{array}{ccccc}\mathrm{\&Phi;}& N& S& M& A_c^T(Z_1+Z_2)\\ N^T& {-Y}_1& 0& 0& 0\\ S^T& 0& {-Y}_1& 0& 0\\ M^T& 0& 0& {-Y}_2& 0\\ (Z_1^T+Z_2^T)A_c& 0& 0& 0& -(Y_1+Y_2)\end{array}\right]<0\end{array}---\left(26\right)$

通过求解式(26)，可以计算出以矩阵P、Q、R、Z₁、Z₂、Y₁和Y₂为变量，以式(25)和式(26)为约束的最小v。最终，利用h=1/v可以求出时滞稳定上限。

实施例2

基于MATLAB仿真软件搭建的图2所示的IEEE4机11节点系统，发电机采用6阶详细模型，励磁系统采用快速励磁，基准模型下的负荷采用50%恒阻抗和50%恒电流模型。首先，通过模态分析法得到四机系统的状态矩阵，并利用SMA方法分别对开、闭环状态矩阵进行降阶，如图3所示。

图3(a)是全阶开环状态矩阵和降阶开环状态矩阵的特征根对比图，(b)是全阶闭环状态矩阵和降阶闭环状态矩阵的特征根对比图。由图3可以看出，开、闭环状态矩阵在利用SMA进行降阶后，均保留了系统中低频振荡成分。因此，利用降阶后系统矩阵能够求解系统允许的最大时滞。将降阶后状态矩阵Α和时滞矩阵Α_d代入式(26)，求得最大时滞边界h=281.88ms。其中，降阶后状态矩阵Α如下：

$A=\left[\begin{array}{cccccc}0& 0& 0& 376.9& 0& 0\\ 0& 0& 0& 0& 376.9& 0\\ 0& 0& 0& 0& 0& 376.9\\ -0.073& 0.065& 0.004& -0.730& 0.272& 0.076\\ 0.058& -0.087& 0.009& 1.160& -0.343& -0.134\\ 0.008& 0.011& -0.085& -0.020& 0.047& -0.554\end{array}\right].$

降阶后时滞矩阵Α_d如下：

$A_d=\left[\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& -0.234& -0.839& 0.010\\ 0& -0.0011& 0.0010& -0.348& -1.362& -0.138\\ 0& 0.0010& 0& 0.049& -0.290& -0.638\end{array}\right].$

时滞分别设置为h₁=50ms，h₂=100ms，h₃=288.88ms和h₄=350ms。所观察的物理量为发电机G1和G4之间的功角差，以及G2和G3之间的功角差，分别如图4和图5所示。

由图4和图5均可看出，在未加广域阻尼控制时，发电机之间的功角发生了严重的低频振荡。加入广域阻尼控制后，在不考虑时滞的情况下，低频振荡可得到有效抑制；然而随着时滞的增加，阻尼效果随之减弱，在最大时滞288ms的情况下，虽然仍存在一定的阻尼，但其阻尼比已经降到10%以下，说明此时控制器已经不满足控制要求。

利用prony算法对各时滞下的功角差曲线进行阻尼比分析，结果如图6和图7所示。

实施例3

基于MATLAB仿真软件搭建的图8所示的IEEE16机68节点系统进一步考查基于广义特征值法的最大时滞求解方法的有效性和通用性。该系统重要联络线为区域4和区域5的联络线1-2，1-27和8-9。发电机采用6阶详细模型，励磁采用IEEE-DC1型励磁，负荷模型15%恒有功功率，25%的恒有功电流和15%的恒无功功率，25%的恒无功功率和60%恒阻抗。首先利用Schur平衡降阶方法，在保证所关心频带的输入输出特性保持不变的前提下，对系统进行降阶，如图9所示。

图9(a)中实线为开环全阶系统的频率响应，虚线为开环降阶系统的频率响应，可以看出，降阶系统和全阶系统的输入输出特性相同。图9(b)中实线为闭环全阶系统的频率响应，虚线为闭环降阶系统的频率响应，可以看出，闭环降阶系统同样保留了全阶系统的特性。因此，利用降阶后的系统矩阵求解全阶系统的时滞稳定上限具有有效性和可行性。

降阶后的状态矩阵Α_re与降阶后的时滞矩阵Α_dre代入式(26)，可得最大的时滞上限为h=93ms。时滞分别设置为h₁=50ms和h₂=93ms。以发电机G1和G16之间的功角差，G3和G14之间的功角差与G10和G15之间的功角差，为所观察的物理量，分别如图10，图11和图12所示。其中，降阶后的状态矩阵Α_re如下：

$A_{\mathrm{re}}=\left[\begin{array}{cccccccccc}-16.1234& 191.6215& -34.7234& -93.6034& -46.9813& -40.7245& 148.6905& -9.4037& 0.5458& -0.1871\\ -8.5822& 73.8511& 38.3796& 28.2807& -31.1118& -3.5328& 3.0915& -0.8165& 0.4670& 0.3871\\ 8.6901& -56.6390& -119.0312& -154.2303& 89.4859& 8.3981& 7.6545& 47.8142& -16.0572& 20.4179\\ -16.3933& 128.3466& 100.7138& 91.8190& -63.4074& 1.0582& -30.2333& -5.4016& 1.8643& -2.0558\\ -21.5589& 183.0718& 41.8903& -25.8071& -3.0065& 15.5680& -58.5954& 54.2491& -17.9821& 24.6753\\ 8.9164& -85.9699& 41.1758& 105.8473& -47.0510& -8.1517& 54.5981& 35.2028& -2.4424& 8.9442\\ -3.3296& 27.0818& 6.4952& -4.4047& 2.3259& 5.3976& -8.7728& 33.8530& -9.1379& 13.1160\\ 0.6362& -4.3077& -5.9066& -6.7725& 3.3248& -1.2021& -0.6273& -11.9533& 5.6377& -0.8597\\ 0.3450& -3.4204& 2.0537& 4.4665& -1.9128& -0.9482& 0.6125& -6.2861& 4.3435& 6.7691\\ 0.4673& -4.1903& 0.9437& 3.9651& -2.1654& -0.7490& 2.0926& -2.3038& -5.0688& -4.4607\end{array}\right].$

降阶后的时滞矩阵Α_dre如下：

$A_{\mathrm{dre}}=\left[\begin{array}{cccccccccc}-0.0168& 0.1387& 0.0401& 0.0041& -0.0139& 0.0056& -0.0470& -0.0035& 0& 0.0023\\ 0.1610& -1.3994& -0.2093& 0.2968& 0.0753& -0.0122& 0.3261& 0.1445& -0.0264& 0.0095\\ 0.0414& -0.2064& -0.5163& -0.6159& 0.2242& -0.0837& 0.2789& -0.2764& 0.0783& -0.0886\\ -0.0247& 0.3531& -0.7123& -1.1881& 0.4669& -0.0484& -0.1044& -0.5840& 0.1524& -0.1677\\ 0& 0.0947& 0.4236& 0.6927& -0.5043& -0.0954& 0.6054& 0.4081& -0.1020& 0.1251\\ 0.1195& -0.6120& 0.4210& 1.2198& -0.9955& -0.3345& 1.9667& 0.7972& -0.1851& 0.2245\\ -0.1383& -0.0010& 0.5268& -0.0149& 1.2620& 0.8615& -4.0937& -0.4567& 0.0720& -0.1352\\ -0.5098& 3.6470& -2.0410& -5.4194& 3.2488& 0.7903& -5.5637& -3.1691& 0.7561& -0.8615\\ 0.5506& -3.8102& 1.9174& 5.4325& -3.4277& -0.9281& 6.2442& 3.2459& -0.7673& 0.8833\\ 1.0523& -7.1184& 4.5809& 11.8740& -7.5930& -1.9539& 13.1229& 7.0687& -1.6838& 1.9525\end{array}\right].$

随着时滞的增加，阻尼效果随之减弱，在最大时滞93ms的情况下，虽然存在一定的阻尼效果，但是其阻尼比已经降到10%以下，此时控制器同样也不能满足控制的要求。利用prony方法对各时滞下的功角差曲线进行阻尼比分析，结果如图13、图14和图15所示。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求的保护范围为准。

Claims (3)

1.一种基于广义特征值的时滞稳定上限计算方法，其特征是所述方法包括：

步骤1：采集建立时滞系统状态方程所需的网络结构参数、发电机频率和发电机功角；

步骤2：建立时滞系统状态方程，并对时滞系统状态方程中的参数矩阵进行降阶处理，得到降阶后的时滞系统状态方程；

步骤3：生成基于改进自由权矩阵的时滞稳定判据；

步骤4：利用时滞稳定判据求解时滞稳定上限；

所述降阶后的时滞系统状态方程为

其中，x(t)为降阶后的电力系统状态向量；

Α为降阶后的电力系统状态矩阵；

Α_d为降阶后的电力系统时滞矩阵；

为降阶后的电力系统状态量对应的状态值；

d(t)为时滞，0≤d(t)≤h且

h为时滞稳定上限且h>0；

μ为时滞最大变化率；

所述基于改进自由权矩阵的时滞稳定判据为：

$\left[\begin{array}{ccccc}\mathrm{\&Phi;}& hN& hS& hM& {\mathrm{hA}}_c^T(Z_1+Z_2)\\ {\mathrm{hN}}^T& -{\mathrm{hZ}}_1& 0& 0& 0\\ {\mathrm{hS}}^T& 0& -{\mathrm{hZ}}_1& 0& 0\\ {\mathrm{hM}}^T& 0& 0& -{\mathrm{hZ}}_2& 0\\ h(Z_1^T+Z_2^T)A_c& 0& 0& 0& -h(Z_1+Z_2)\end{array}\right]<0;$

其中，

${\mathrm{\&Phi;}}_1=\left[\begin{array}{ccc}PA+A^{T}P+Q+R& {\mathrm{PA}}_d& 0\\ A_d^{T}P& -(1-\mathrm{\&mu;})Q& 0\\ 0& 0& -R\end{array}\right];$

Φ₂＝[N+M -N+S -M-S]；

A_c＝[A A_d 0]；

N、M和S为改进自由权矩阵；

改进自由权矩阵N满足

改进自由权矩阵S满足

改进自由权矩阵M满足

P、Q、R、Z₁和Z₂为待定矩阵。

2.根据权利要求1所述的计算方法，其特征是所述利用时滞稳定判据求解时滞稳定上限包括：

子步骤A1：将时滞稳定判据等价变换为

其中，Y₁和Y₂为附加矩阵，并且Y₁＝Y₁ ^T≥0，Y₂＝Y₂ ^T≥0，v＝1/h；

子步骤A2：以v最小为目标，以时滞稳定判据等价变换后的矩阵不等式以及为约束条件，计算待定矩阵P、Q、R、Z₁、Z₂、Y₁和Y₂，进而得到时滞稳定上限h。

3.一种用于执行权利要求1所述方法的基于广义特征值的时滞稳定上限计算系统，其特征是所述系统包括顺序相连的数据采集模块、时滞系统处理模块、时滞上限求解模块和结果输出模块；

所述数据采集模块用于采集建立时滞系统状态方程所需的网络结构参数、发电机频率和发电机功角，并将采集的数据发送至时滞系统处理模块；

所述时滞系统处理模块用于建立时滞系统状态方程，并对时滞系统状态方程中的参数矩阵进行降阶处理；

所述时滞上限求解模块用于生成基于改进自由权矩阵的时滞稳定判据，并利用时滞稳定判据求解时滞稳定上限；

所述结果输出模块用于输出时滞稳定上限结果。