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所述方法的基于广义特征值的时滞稳定上限计算***，其特征是所述***包括顺序相连的数据采集模块、时滞***处理模块、时滞上限求解模块和结果输出模块；

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

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

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

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