时滞电力系统 Lyapunov 稳定性分析方法

摘要

本发明属于电力系统技术领域，为解决原有时滞电力系统稳定分析计算效率低下的问题，有效降低时滞微分方程的维数，具有更高计算效率，为此，本发明采取的技术方案是，时滞电力系统Lyapunov稳定性分析方法，包括如下步骤：建立电力系统的带约束时滞微分方程(CTODE)模型：将系统状态按不考虑时滞影响的状态在前，考虑时滞影响的状态在后的方式重新排列整理，就得到原时滞系统所对应的带约束时滞微分方程(CTODE)模型；基于带约束时滞微分方程模型的新稳定判据。本发明主要应用于电力系统。

说明书

技术领域

本发明属于电力系统技术领域，涉及一种时滞电力系统Lyapunov稳定性分析方法。

背景技术

在自然界中，系统未来的发展趋势既取决于当前状态，也与过去状态有关，这类现象称 为时滞[1，2]。时滞现象广泛存在于电力系统各个环节，但传统控制信号主要取自本地量测 装置，时滞很小，通常不予考虑。但在广域环境下，远方量测环节的时滞非常明显，因此研 究时滞对电力系统稳定性的影响具有十分重要的现实意义[3,4]。

已有针对电力系统的时滞稳定研究，多基于如下时滞微分方程(TODE)模型来开展：

$\Delta \stackrel{·}{z}=\tilde{A}·\mathrm{\Delta z}+\underset{i=1}{\overset{k}{\Sigma }}{\tilde{A}}_i·{\mathrm{\Delta z}}_{\mathrm{\tau i}}$   (1)

其中：z∈Rⁿ为系统状态变量，是含n个元素的实数向量，n为状态变量个数，R为实数(下 同)；z_τi＝[z₁(t-τ_i),...,z_n(t-τ_i)]∈Rⁿ,τ_i∈R,i＝1,2,...,k为时滞系数，为常数矩阵。如：文献[4,5]给出了上述模型的详细推导过程，并借助其研究了电力系统时滞 小扰动稳定域和扩展时滞小扰动稳定域。文献[6-7]利用该模型，通过优化算法求解时滞系统 的特征值轨迹，进而研究时滞电力系统的小扰动稳定性。文献[8-9]采用该模型和Rekasius变 换，给出了一种精确求解时滞系统稳定裕度的有效方法。文献[10]则利用该模型，给出了一 种基于复矩阵变换的时滞稳定裕度求解方法，所得结果与文献[8-9]完全一致，但具有更高计 算效率，并能提供更多有用信息。文献[11-12]采用该模型研究了在考虑控制回路参数随机变 动时的电力系统概率时滞稳定性。[13-14]则采用该模型研究了广域反馈控制方法，以提高电 力系统的运行稳定性。上述研究，均假设时滞为常数或可表示为简单函数，当这一条件不满 足时，则只能采取Lyapunov方法来分析系统的时滞稳定性，由此又形成两类主要分析方法： 一类是基于Lyapunov-Razumikhin定理[1-2]的方法，如[15-16]；另一类则是基于 Lyapunov-Krasovskii定理[1-2]的方法，如[17-18]，由于Lyapunov-Krasovskii类方法考虑了时 滞轨迹的影响，因此一般具有更小的保守性，因而成为近年来关注的热点[1-2,19]。

上述方法在采用式(1)所示模型开展研究时，存在如下问题：现代电力系统的规模极其庞 大，因此其动态方程的维数极高，即式(1)中向量z和矩阵的维数都很高。但我们知道， 在进行电力系统广域控制器设计时，一般只需采集少量远方数据，如南方电网公司在进行直 流系统广域协调控制器设计时，仅实时采集直流两端数个关键动态参数，它们的传输时滞需 要考虑，而其数目远小于整个南方电网公司电力系统动态模型的维数[20]；再如基于广域测 量系统信息进行电力系统稳定器协调控制器设计时，控制器远程输入量也仅为系统几个关键 点的测量信息(如远方节点的频率或输电断面潮流)，其数目也远小于系统动态模型的维数[3, 21-22]。也即式(1)中真正起作用的时滞变量的数目将远小于系统状态变量的个数，换言之， 式(1)中矩阵中的非零元素极少。但已有方法在分析计算时，均未考虑这一情况， 而将认为满秩，由此造成很多无谓计算，严重影响了相关方法的计算效率。

发明内容

本发明旨在克服现有技术的不足，解决原有时滞电力系统稳定分析计算效率低下的问题， 为此，本发明采取的技术方案是，基于带约束时滞微分方程模型的时滞电力系统Lyapunov稳 定性分析方法，包括如下步骤：

建立电力系统的带约束时滞微分方程（CTODE）模型

时滞电力系统

$\stackrel{·}{z}=F(z,z_{\tau })$

其中：z＝[z₁,z₂,...,z_n]∈Rⁿ为系统状态向量，向量中的元素个数为n，Rⁿ表示n维实数 向量；z_τ＝(z_τ1,...,z_τi,...,z_τk)，其中的z_τi＝[z₁(t-τ_i),...,z_n(t-τ_i)]∈Rⁿ，τ_i∈R,i＝1,2,...,k为时 滞系数；

将系统状态按不考虑时滞影响的状态在前，考虑时滞影响的状态在后的方 式重新排列整理，就得到原时滞系统所对应的带约束时滞微分方程(CTODE)模型：

${\stackrel{·}{z}}_1=F_1(z_1,z_2)$

${\stackrel{·}{z}}_2=F_2(z_1,z_2,z_{2,\tau })$

其中：z＝[z₁,z₂]，为不考虑时滞影响的系统状态向量，是含n₁个元素的实数向量，n₁为不考虑时滞影响的状态变量的数目；为考虑 时滞影响的系统状态向量，是含n₂个元素的实数向量，n₂为考虑时滞影响的状态变量的数目， n＝n₁+n₂为状态向量z的元素个数；z_2,τ＝(z_2,τ1,...,z_2,τi,...,z_2,τk)，其中的时滞状态向量 τ_i∈R,i＝1,2,...,k为时滞系数；

进一步，在系统平衡点处对其线性化，可得

$\Delta {\stackrel{·}{z}}_1=A_{11}\Delta z_1+A_{12}{\mathrm{\Delta z}}_2$

$\Delta {\stackrel{·}{z}}_2=A_{21}{\mathrm{\Delta z}}_1+A_{22}{\mathrm{\Delta z}}_2+\underset{i=1}{\overset{k}{\Sigma }}{A}_{d,i}{\mathrm{\Delta z}}_{2,\mathrm{\tau i}}$

其中：$A_{11}=\frac{{\partial F}_1}{{\partial z}_1},A_{12}\frac{{\partial F}_1}{{\partial z}_2},A_{21}=\frac{{\partial F}_2}{{\partial z}_1},A_{22}=\frac{{\partial F}_2}{{\partial z}_2},A_{d,i}=\frac{{\partial F}_2}{{\partial z}_{2,\mathrm{\tau i}}}$

即得到时滞系统CTODE模型的线性化形式。

基于带约束时滞微分方程模型的新稳定判据

对于时滞系统的CTODE线性化模型

$\Delta {\stackrel{·}{z}}_1=A_{11}\Delta z_1+A_{12}{\mathrm{\Delta z}}_2$

$\Delta {\stackrel{·}{z}}_2=A_{21}{\mathrm{\Delta z}}_1+A_{22}{\mathrm{\Delta z}}_2+\underset{i=1}{\overset{k}{\Sigma }}{A}_{d,i}{\mathrm{\Delta z}}_{2,\mathrm{\tau i}}$

当k=1时，如下定理给出了该系统稳定的条件：

定理：给定标量τ₁>0，若存在如下对称正定矩阵：分别称为对称第一矩阵、第二矩阵、...第六矩阵，和任 意合适维数的矩阵P₁₂,N₁,N₂,X₁₂，分别称为一般第一矩阵、一般第二矩阵...一般第四矩阵，使 得下式成立，则时滞系统在时滞为τ₁时渐进稳定：

$\bar{\Phi }=\left(\begin{array}{ccc}{\bar{\Phi }}_{11}& {\bar{\Phi }}_{12}& {\bar{\Phi }}_{13}\\ {\bar{\Phi }}_{12}^T& {\bar{\Phi }}_{22}& {\bar{\Phi }}_{23}\\ {\bar{\Phi }}_{13}^T& {\bar{\Phi }}_{23}^T& {\bar{\Phi }}_{33}\end{array}\right)<0,$$\bar{\Psi }=\left(\begin{array}{ccc}{X}_{11}& X_{12}& N_1\\ X_{12}^T& X_{22}& N_2\\ N_1^T& N_2^T& Z\end{array}\right)>0,$

其中，

${\bar{\Phi }}_{11}=P_{11}{A}_{11}+A_{11}^T{P}_{11}+P_{12}{A}_{21}+A_{21}^T{P}_{12}^T+{\tau }_1{A}_{21}^T{\mathrm{ZA}}_{21}$

${\bar{\Phi }}_{12}=P_{11}{A}_{12}+P_{12}{A}_{22}+A_{11}^T{P}_{12}+{A_{21}^T{P}_{22}+\tau }_1{A}_{21}^T{\mathrm{ZA}}_{22}$

${\bar{\Phi }}_{13}=P_{12}{A}_{d,1}+{\tau }_1{A}_{21}^T{\mathrm{ZA}}_{d,1}$

${\bar{\Phi }}_{22}=P_{12}^T{A}_{12}+A_{12}^T{P}_{12}+P_{22}{A}_{22}+A_{22}^T{P}_{22}+Q+N_1$

${+N}_1^T+{\tau }_1{X}_{11}+{\tau }_1{A}_{22}^T{\mathrm{ZA}}_{22}$

${\bar{\Phi }}_{23}=P_{22}{A}_{d,1}-N_1+N_2+{\tau }_1{X}_{12}+{\tau }_1{A}_{22}^T{\mathrm{ZA}}_{d,1}$

${\bar{\Phi }}_{33}=-Q-N_2-N_2^T+{\tau }_1{X}_{22}+{\tau }_1{A}_{d,1}^T{\mathrm{ZA}}_{d,1}$

P₁₁,P₂₂,Q,Z,X₁₁,X₂₂和P₁₂,N₁,N₂,X₁₂矩阵为线性矩阵不等式系统的算法条件。

建立电力系统带约束时滞微分方程模型进一步具体为：

1、建立时滞电力系统的带约束时滞微分代数方程（CTDAE）模型

令：和分别代表系统中不考虑时滞影响的状态向量和代数向量，m₁为不考 虑时滞影响的代数量个数，和分别代表考虑时滞影响的状态向量和代数向量， m₂为考虑时滞影响的代数量个数，按不考虑时滞的相关量在前，考虑时滞的相关量在后，建 立微分代数方程模型，形式如下

${\stackrel{·}{z}}_1=F_1(z_1,z_2,y_1,y_2)$

${\stackrel{·}{z}}_2=F_2(z_1,z_2,y_1,y_2,z_{2,\tau },y_{2,\tau })$

0＝G₁(z₁,z₂,y₁,y₂)

0＝G₂(z₂,y₂)

0＝G_2,i(z_2,τi,y_2,τi),i＝1,2,...,k

式中，z_2,τ＝(z_2,τ1,z_2,τ2,...,z_2,τi,...,z_2,τk)，$z_{2,\mathrm{\tau i}}=[z_{21}\left({t-\tau }_i\right),...,z_{{2n}_2}\left({t-\tau }_i\right)]$为系统 的时滞状态变量所构成的向量；y_2,τ＝(y_2,τ1,y_2,τ2,...,y_2,τi,...,y_2,τk)， 为系统的时滞代数变量所构成的向量；G₁(·),G₂(·)对应当前时 刻的代数约束，G_2,i(·)则对应τ_i时刻前的代数约束。

2、在系统平衡点处对步骤1中的微分代数方程线性化，可得：

$\Delta {\stackrel{·}{z}}_1={\tilde{A}}_{11}{\mathrm{\Delta z}}_1+{\tilde{A}}_{12}{\mathrm{\Delta z}}_2+{\tilde{B}}_{11}{\mathrm{\Delta y}}_1+{\tilde{B}}_{12}{\mathrm{\Delta y}}_2$

$\Delta {\stackrel{·}{z}}_2={\tilde{A}}_{21}{\mathrm{\Delta z}}_1+{\tilde{A}}_{22}{\mathrm{\Delta z}}_2+{\tilde{B}}_{21}{\mathrm{\Delta y}}_1+{\tilde{B}}_{21}{\mathrm{\Delta y}}_2$

$+\underset{i=1}{\overset{k}{\Sigma }}{\tilde{A}}_{d,i}{\mathrm{\Delta z}}_{2,\mathrm{\tau i}}+\underset{i=1}{\overset{k}{\Sigma }}{\tilde{B}}_{d,i}{\mathrm{\Delta y}}_{2,\mathrm{\tau i}}$

$0={\tilde{C}}_{11}{\mathrm{\Delta z}}_1+{\tilde{C}}_{12}{\mathrm{\Delta z}}_2+{\tilde{D}}_{11}{\mathrm{\Delta y}}_1+{\tilde{D}}_{12}{\mathrm{\Delta y}}_2$

$0={\tilde{C}}_{22}{\mathrm{\Delta z}}_2+{\tilde{D}}_{22}{\mathrm{\Delta y}}_2$

$0={\tilde{C}}_{22,i}{\mathrm{\Delta z}}_{2,\mathrm{\tau i}}+{\tilde{D}}_{22,i}{\mathrm{\Delta y}}_{2,\mathrm{\tau i}},i=\mathrm{1,2},...k$

其中：${\tilde{A}}_{i,j}=\frac{{\partial F}_i}{{\partial z}_j},{\tilde{B}}_{i,j}=\frac{{\partial F}_i}{{\partial y}_j},{\tilde{C}}_{i,j}=\frac{{\partial G}_i}{{\partial z}_j},{\tilde{D}}_{i,j}=\frac{{\partial G}_i}{{\partial y}_j};i,j=\mathrm{1,2};$${\tilde{A}}_{d,i}=\frac{{\partial F}_2}{{\partial z}_{2,\mathrm{\tau i}}},{\tilde{B}}_{d,i}=\frac{{\partial F}_2}{{\partial y}_{2,\mathrm{\tau i}}},{\tilde{C}}_{22,i}=\frac{{\partial G}_{2,i}}{{\partial z}_{2,\mathrm{\tau i}}},{\tilde{D}}_{22,i}=\frac{{\partial G}_{2,i}}{{\partial y}_{2,\mathrm{\tau i}}}.$由于在[-τ_max,0]时间段内，隐函数定理总成立，则矩阵可逆，其中：

$\tilde{D}=\left(\begin{array}{cc}{\tilde{D}}_{11}& {\tilde{D}}_{12}\\ 0& {\tilde{D}}_{22}\end{array}\right)$

并令：

${\Omega }_{11}=-{\tilde{D}}_{11}^{-1},{\Omega }_{22}=-{\tilde{D}}_{22}^{-1}$

${\Omega }_{22,i}=-{\tilde{D}}_{22,i}^{-1},i=\mathrm{1,2},...,k$

由此可得：

Δy₁＝K₁₁Δx₁+K₁₂Δx₂

Δy₂＝K₂₂Δx₂

Δy_2,τi＝K_22,iΔx_2,τi＝Ω_22,iC_22,iΔx_2,τi

其中：

$K_{11}={\Omega }_{11}{\tilde{C}}_{11}$

$K_{12}={\Omega }_{11}({\tilde{C}}_{12}+{\tilde{D}}_{12}{\Omega }_{22}{\tilde{C}}_{22})$

$K_{22}={\Omega }_{22}{\tilde{C}}_{22}$

3、整理得到系统的带约束时滞微分方程（CTODE）模型如下：

$\Delta {\stackrel{·}{z}}_1=A_{11}{\mathrm{\Delta z}}_1+A_{12}{\mathrm{\Delta z}}_2$

$\Delta {\stackrel{·}{z}}_2=A_{21}{\mathrm{\Delta z}}_1+A_{22}{\mathrm{\Delta z}}_2+\underset{i=1}{\overset{k}{\Sigma }}{A}_{d,i}{\mathrm{\Delta z}}_{2,\mathrm{\tau i}}$

其中：

$A_{11}={\tilde{A}}_{11}+{\tilde{B}}_{11}{K}_{11}$

$A_{12}={\tilde{A}}_{12}+{\tilde{B}}_{11}{K}_{12}+{\tilde{B}}_{12}{K}_{22}$

$A_{21}={\tilde{A}}_{21}+{\tilde{B}}_{21}{K}_{11}$

$A_{22}={\tilde{A}}_{22}+{\tilde{B}}_{21}{K}_{12}+{\tilde{B}}_{22}{K}_{22}$

$A_{d,i}={\tilde{A}}_{d,i}+{\tilde{B}}_{d,i}{K}_{22,i}$

基于带约束时滞微分方程模型的新稳定判据进一步包括：

4、建立用于Lyapunov判稳的线性矩阵不等式系统

4.1定义一个标量τ₁，定义第一对称矩阵变量到第六对称矩阵变量分别为Q＝Q^T，Z＝Z^T，以及一般第一矩阵变量到一般第四矩 阵变量分别为P₁₂、X₁₂、N₁、N₂，式中上标“T”均指矩阵转置；

4.2描述一个线性矩阵不等式系统

P₁₁＞0，P₂₂＞0，X₁₁＞0，X₂₂＞0，Q＞0，Z＞0，$\bar{\Phi }=\left(\begin{array}{ccc}{\bar{\Phi }}_{11}& {\bar{\Phi }}_{12}& {\bar{\Phi }}_{13}\\ {\bar{\Phi }}_{12}^T& {\bar{\Phi }}_{22}& {\bar{\Phi }}_{23}\\ {\bar{\Phi }}_{13}^T& {\bar{\Phi }}_{23}^T& {\bar{\Phi }}_{33}\end{array}\right)<0,$

$\bar{\Psi }=\left(\begin{array}{ccc}{X}_{11}& X_{12}& N_1\\ X_{12}^T& X_{22}& N_2\\ N_1^T& N_2^T& Z\end{array}\right)>0,$

其中：

${\bar{\Phi }}_{11}=P_{11}{A}_{11}+A_{11}^T{P}_{11}+P_{12}{A}_{21}+A_{21}^T{P}_{12}^T+{\tau }_1{A}_{21}^T{\mathrm{ZA}}_{21}$

${\bar{\Phi }}_{12}=P_{11}{A}_{12}+P_{12}{A}_{22}+A_{11}^T{P}_{12}+{A_{21}^T{P}_{22}+\tau }_1{A}_{21}^T{\mathrm{ZA}}_{22}$

${\bar{\Phi }}_{13}=P_{12}{A}_{d,1}+{\tau }_1{A}_{21}^T{\mathrm{ZA}}_{d,1}$

${\bar{\Phi }}_{22}=P_{12}^T{A}_{12}+A_{12}^T{P}_{12}+P_{22}{A}_{22}+A_{22}^T{P}_{22}+Q+N_1$${+N}_1^T+{\tau }_1{X}_{11}+{\tau }_1{A}_{22}^T{\mathrm{ZA}}_{22}$

${\bar{\Phi }}_{23}=P_{22}{A}_{d,1}-N_1+N_2+{\tau }_1{X}_{12}+{\tau }_1{A}_{22}^T{\mathrm{ZA}}_{d,1}$

${\bar{\Phi }}_{33}=-Q-N_2-N_2^T+{\tau }_1{X}_{22}+{\tau }_1{A}_{d,1}^T{\mathrm{ZA}}_{d,1}$

式中上标“T”均指矩阵转置。

5、对于标量τ₁，利用线性矩阵不等式（LMI）求解工具，判断步骤4.2所给线性不等式 的可行性，如果可行，则当时滞常数为τ₁时，系统在平衡点附近渐进稳定。

本发明的技术特点及效果：

基于CTODE模型的Lyapunov-Krasovskii稳定分析方法，有效降低了时滞微分方程的维 数，待求变量数大为减少，因而具有更高计算效率。

附图说明

图1新方法加速比随维数的变化情况；

Fig.1Speedup of new criterion。

图2两方法待求变量数随矩阵A维数变化(n2=2)；

Fig.2Unknown variables of two criteria changing with dimension of matrix A(n2=2)。

图3两种方法的计算误差变化情况；

Fig.3Errors of two criteria。

图4工作原理流程框图。

具体实施方式

为解决原有时滞电力系统稳定分析计算效率低下的问题，本专利采用带约束时滞微分方 程模型CTODE来描述时滞电力系统的动态，基于CTODE模型进行Lyapunov稳定性分析， 本方法可以提高含时滞电力系统稳定性的分析效率。

一、方法原理说明

1、时滞电力系统的CTODE模型

对于时滞电力系统

$\stackrel{.}{z}=F\left(z{,z}_{\tau }\right)$   (2)

其中：z＝[z₁,z₂,...,z_n]∈Rⁿ为系统状态向量，向量中的元素个数为n，Rⁿ表示n维实数 向量(下同)；z_τ＝[z_τ1,...,z_τi,...,z_τk]，其中的z_τi＝[z₁(t-τ_i),...,z_n(t-τi₎]∈Rⁿ，τ_i∈R,i＝1,2,...,k 为时滞系数。

将系统状态按不考虑时滞影响的状态在前，考虑时滞影响的状态在后的方 式重新排列整理，就得到原时滞系统所对应的带约束时滞微分方程(CTODE)模型：

${\stackrel{·}{z}}_1=F_1(z_1,z_2)$  (3a)

${\stackrel{·}{z}}_2=F_2(z_1,z_2,z_{2\tau })$  (3b)

其中：z＝[z₁,z₂]，为不考虑时滞影响的系统状态向量，是含n₁个元素的实数向量，n₁为不考虑时滞影响的状态变量的数目；为考虑 时滞影响的系统状态向量，是含n₂个元素的实数向量，n₂为考虑时滞影响的状态变量的数目， n＝n₁+n₂为状态向量z的元素个数；z_2,τ＝(z_2,τ1,...,z_2,τi,...,z_2,τk)，其中的时滞状态向量 τ_i∈R,i＝1,2,...,k为时滞系数。（符号说明下同。）

$\Delta {\stackrel{·}{z}}_1=A_{11}\Delta z_1+A_{12}{\mathrm{\Delta z}}_2$  (4a)

$\Delta {\stackrel{·}{z}}_2=A_{21}{\mathrm{\Delta z}}_1+A_{22}{\mathrm{\Delta z}}_2+\underset{i=1}{\overset{k}{\Sigma }}{A}_{d,i}{\mathrm{\Delta z}}_{2,\mathrm{\tau i}}$  (4b)

其中：

$A_{11}=\frac{{\partial F}_1}{{\partial z}_1},A_{12}\frac{{\partial F}_1}{{\partial z}_2},A_{21}=\frac{{\partial F}_2}{{\partial z}_1},A_{22}=\frac{{\partial F}_2}{{\partial z}_2},A_{d,i}=\frac{{\partial F}_2}{{\partial z}_{2,\mathrm{\tau i}}}$

观察式(3)可以看出，系统时滞环节的动态由式(3b)决定，同时它又被完全约束在式(3a)对 应的动态流形(Dynamic Manifold)上，因此，可将式(3a)视作时滞环节的一个动态约束，因此 称该模型为带约束时滞微分方程CTODE。对线性模型(4)可做类似分析，不再赘述。

从CTODE模型的构成可看出，它将原时滞模型(2)的状态空间，通过变换分解为不含时 滞子空间和含时滞子空间，前者的状态变量为，其动态仅由微分方程(3a)来描述；后 者对应的状态量为，其动态则由一个滞后型时滞微分方程(3b)来描述。由于引起时滞 系统分析计算困难的根本原因在于时滞微分方程^{[1-2，22-33]}，对比式(2)和式(3)不难看出，前者 时滞微分方程的维数为n，而后者为n₂，且一般有n₂远小于n。由于对微分方程进行分析计算 相对简单，而在CTODE模型中，时滞微分方程的维数大为降低，由此必然能提升算法分析 计算的效率。

2、电力系统CTODE模型的建立方法

含时滞的电力系统动态，往往需用如下时滞微分代数方程(TDAE)来加以描述。

$\stackrel{·}{z}=F(z,z_{\tau }y,y_{\tau })$  (5a)

0＝G(z,y)    (5b)

0＝G_i(z_τi,y_τi),i＝1,2,...,k    (5c)

式中：z,z_τ的含义同式(2)；y∈R^m为系统代数变量，m为代数变量个数；y_τ＝(y_τ1,...,y_τi,...,y_τk)， 其中的y_τi＝[y₁(t-τ_i),...,y_m(t-τ_i)]∈R^m,i＝1,2,...,k为系统的时滞代数变量构成的向量；G(·)对应当 前时刻的代数约束；G_i(·)则对应τ_i时刻前的代数约束。为描述方便，令τ＝[τ₁,τ₂,...,τ_k]^T为系统 时滞向量，τ_max＝max(τ)为系统最大时滞，并约定系统在[-τ_max,0]时间段内均存在正常解，即在 此时间段内隐函数定理均成立^[1-2，19]。由此，理论上可从(5b)和(5c)中直接解出y,z和y_τi,z_τi间 的关系如下：

y＝h(z)

y_τi＝h_i(z_τi),i＝1,2,...,k

将其代入式(5)，即可得到：

$\stackrel{·}{z}=F(z,z_{\tau },h\left(z\right),h_1\left(z_{\tau 1}\right),...,h_k\left(z_{\mathrm{\tau k}}\right))=f\left(z{,z}_{\tau }\right)$  (6)

对于复杂的时滞电力系统，往往难以通过简单观察就能实现其时滞环节与其他动态环节 的有效解耦，为此采用如下方法推导其对应的线性时滞模型。设式(6)在其平衡点附近的线性 化模型为

$\Delta \stackrel{·}{z}=\tilde{A}·\mathrm{\Delta z}+\underset{i=1}{\overset{k}{\Sigma }}{\tilde{A}}_i·{\mathrm{\Delta z}}_{\mathrm{\tau i}}$  (7)

考虑到其中的矩阵i＝1,2,...,k非常稀疏，因此引入变换矩阵T，并令：

x＝T·Δz  (8a)

x_τi＝T·Δz_τi,i＝1,2,...,k  (8b)

将上述变换代入式(7)，经整理可得：

$\stackrel{·}{x}=\mathrm{Ax}+\underset{i=1}{\overset{k}{\Sigma }}{A}_i{x}_{\mathrm{\tau i}}$  (9)

其中：$A=T·\tilde{A}·T^{-1}$  (10a)

$A_i=T·{\tilde{A}}_i·T^{-1},i=\mathrm{1,2},...,k$  (10b)

同时，通过优选变换矩阵T，可使上式中的矩阵A_i具有如下形式：

$A_i=\left(\begin{array}{cc}0& 0\\ 0& A_{d,i}\end{array}\right)$  (11)

其中则式(9)重写为如下形式：

${\stackrel{·}{x}}_1=A_{11}{x}_1+A_{12}{x}_2$  (12a)

${\stackrel{·}{x}}_2=A_{21}{x}_1+A_{22}{x}_2+\underset{i=1}{\overset{k}{\Sigma }}{A}_{d,i}{x}_{2,\mathrm{\tau i}}$  (12b)

其中：矩阵$A=\left(\begin{array}{cc}{A}_{11}& A_{12}\\ A_{21}& A_{22}\end{array}\right),$$x={[x_1^T,x_2^T]}^T,x_1\in R^{n_1},x_2\in R^{n_2},x_{2,\mathrm{\tau i}}\in R^{n_2},n=n_1+n_2.$

3、基于CTODE模型的新稳定判据

对于时滞系统(12)，当k=1时，如下定理给出了该系统稳定的条件。

定理：给定标量τ₁>0，若存在如下对称正定矩阵：分别称为对称第一矩阵、第二矩阵、…第六矩阵，和任 意合适维数的矩阵P₁₂,N₁,N₂,X₁₂，分别称为一般第一矩阵、一般第二矩阵...一般第四矩阵，使 得下式成立，则时滞系统在时滞为τ₁时渐进稳定：

$\bar{\Phi }=\left(\begin{array}{ccc}{\bar{\Phi }}_{11}& {\bar{\Phi }}_{12}& {\bar{\Phi }}_{13}\\ {\bar{\Phi }}_{12}^T& {\bar{\Phi }}_{22}& {\bar{\Phi }}_{23}\\ {\bar{\Phi }}_{13}^T& {\bar{\Phi }}_{23}^T& {\bar{\Phi }}_{33}\end{array}\right)<0$  (13a)

$\bar{\Psi }=\left(\begin{array}{ccc}{X}_{11}& X_{12}& N_1\\ X_{12}^T& X_{22}& N_2\\ N_1^T& N_2^T& Z\end{array}\right)>0$  (13b)

其中，

${\bar{\Phi }}_{11}=P_{11}{A}_{11}+A_{11}^T{P}_{11}+P_{12}{A}_{21}+A_{21}^T{P}_{12}^T+{\tau }_1{A}_{21}^T{\mathrm{ZA}}_{21}$

${\bar{\Phi }}_{12}=P_{11}{A}_{12}+P_{12}{A}_{22}+A_{11}^T{P}_{12}+{A_{21}^T{P}_{22}+\tau }_1{A}_{21}^T{\mathrm{ZA}}_{22}$

${\bar{\Phi }}_{13}=P_{12}{A}_{d,1}+{\tau }_1{A}_{21}^T{\mathrm{ZA}}_{d,1}$

${\bar{\Phi }}_{22}=P_{12}^T{A}_{12}+A_{12}^T{P}_{12}+P_{22}{A}_{22}+A_{22}^T{P}_{22}+Q+N_1$

${+N}_1^T+{\tau }_1{X}_{11}+{\tau }_1{A}_{22}^T{\mathrm{ZA}}_{22}$

${\bar{\Phi }}_{23}=P_{22}{A}_{d,1}-N_1+N_2+{\tau }_1{X}_{12}+{\tau }_1{A}_{22}^T{\mathrm{ZA}}_{d,1}$

${\bar{\Phi }}_{33}=-Q-N_2-N_2^T+{\tau }_1{X}_{22}+{\tau }_1{A}_{d,1}^T{\mathrm{ZA}}_{d,1}$

P₁₁,P₂₂,Q,Z,X₁₁,X₂₂和P₁₂,N₁,N₂,X₁₂矩阵为线性矩阵不等式系统的算法条件。

证明：

对于系统(12)，当k=1时，构造如下Lyapunov-Krasovskii泛函：

V(t,τ₁)＝V₁(t)+V₂(t,τ₁)+V₃(t,τ₁)(14)

其中：

$V_1\left(t\right)=[x_1^T\left(t\right),x_2^T\left(t\right)]\left(\begin{array}{cc}{P}_{11}& P_{12}\\ P_{12}^T& P_{22}\end{array}\right)\left(\begin{array}{c}{x}_1\left(t\right)\\ x_2\left(t\right)\end{array}\right)$  (15a)

$V_2(t,{\tau }_1)={\int }_{t{\tau }_1}^t{x}_2^T\left(s\right)Q{x}_2\left(s\right)·\mathrm{ds}$  (15b)

$V_3(t,{\tau }_1)={\int }_{-{\tau }_1}^0{\int }_{t+ϵ}^t{\stackrel{·}{x}}_2^T\left(s\right)Z{\stackrel{·}{x}}_2\left(s\right)·\mathrm{ds}·\mathrm{dϵ}$  (15c)

沿系统轨迹对V(t,τ₁)求导，在其导函数推导过程中，引入如下松弛项：

$x_2^T\left(t\right)N_1[x_2\left(t\right)-x_2(t-{\tau }_1)-{\int }_{{t-\tau }_1}^t{\stackrel{·}{x}}_2\left(s\right)\mathrm{ds}]=0$

$x_2^T(t-{\tau }_1)N_2[x_2\left(t\right)-x_2(t-{\tau }_1)-{\int }_{{t-\tau }_1}^t{\stackrel{·}{x}}_2\left(s\right)\mathrm{ds}]=0$

${\tau }_1{x}_2^T\left(t\right)[X_{11}-X_{11}]x_2\left(t\right)=0$

${\tau }_1{x}_2^T(t-{\tau }_1)[X_{22}-X_{22}]x_2(t-{\tau }_1)=0$

${\tau }_1{x}_2^T\left(t\right)[X_{12}-X_{12}]x_2(t-{\tau }_1)=0$

经仔细推导，最后得出如下结果：

$\stackrel{·}{V}(t,{\tau }_1)={\eta }^T\bar{\Phi }\eta -{\int }_{{t-\tau }_1}^t{\xi }^T\bar{\Psi }\xi ·\mathrm{ds}$

其中：η＝[x₁,x₂,x₂(t-τ₁)]^T,不难看出，若式(13a)和(13b)成立，则有进一步由Lyapunov稳定性理论可知，此时系统(12)渐进稳定，定理得证。

二、基于带约束时滞微分方程模型（CTODE）的电力系统Lyapunov稳定性分析步骤

1、建立时滞电力系统的微分代数方程模型（CTDAE）

令：分别代表系统中不考虑时滞影响的状态向量和代数向量（n₁为不考 虑时滞影响的状态量个数，m₁为不考虑时滞影响的代数量个数），分别代表 考虑时滞影响的状态向量和代数向量（n₂为考虑时滞影响的状态量个数，m₂为考虑时滞影响 的代数量个数）。按不考虑时滞的相关量在前，考虑时滞的相关量在后，建立微分代数方程模 型，形式如下

${\stackrel{·}{z}}_1=F_1(z_1,z_2,y_1,y_2)$  (16a)

${\stackrel{·}{z}}_2=F_2(z_1,z_2,y_1,y_2,z_{2,\tau },y_{2,\tau })$  (16b)

0＝G₁(z₁,z₂,y₁,y₂)  (16c)

0＝G₂(z₂,y₂)  (16d)

0＝G_2,i(z_2,τi,y_2,τi),i＝1,2,...,k  (16e)

式中，z_2,τ＝(z_2,τ1,z_2,τ2,...,z_2,τi,...,z_2,τk)，为系 统的时滞状态变量所构成的向量；y_2,τ＝(y_2,τ1,y2_,τ2,...,y_2,τi,...,y_2,τk)， 为系统的时滞代数变量所构成的向量；G₁(·),G₂(·)对应当前时 刻的代数约束，G_2,i(·)则对应τ_i时刻前的代数约束。

2、在系统平衡点处对式(16)线性化，可得：

$\Delta {\stackrel{·}{z}}_1={\tilde{A}}_{11}{\mathrm{\Delta z}}_1+{\tilde{A}}_{12}{\mathrm{\Delta z}}_2+{\tilde{B}}_{11}{\mathrm{\Delta y}}_1+{\tilde{B}}_{12}{\mathrm{\Delta y}}_2$  (17a)

$\Delta {\stackrel{·}{z}}_2={\tilde{A}}_{21}{\mathrm{\Delta z}}_1+{\tilde{A}}_{22}{\mathrm{\Delta z}}_2+{\tilde{B}}_{21}{\mathrm{\Delta y}}_1+{\tilde{B}}_{21}{\mathrm{\Delta y}}_2$  (17b)

$+\underset{i=1}{\overset{k}{\Sigma }}{\tilde{A}}_{d,i}{\mathrm{\Delta z}}_{2,\mathrm{\tau i}}+\underset{i=1}{\overset{k}{\Sigma }}{\tilde{B}}_{d,i}{\mathrm{\Delta y}}_{2,\mathrm{\tau i}}$

$0={\tilde{C}}_{11}{\mathrm{\Delta z}}_1+{\tilde{C}}_{12}{\mathrm{\Delta z}}_2+{\tilde{D}}_{11}{\mathrm{\Delta y}}_1+{\tilde{D}}_{12}{\mathrm{\Delta y}}_2$  (17c)

$0={\tilde{C}}_{22}{\mathrm{\Delta z}}_2+{\tilde{D}}_{22}{\mathrm{\Delta y}}_2$  (17d)

$0={\tilde{C}}_{22,i}{\mathrm{\Delta z}}_{2,\mathrm{\tau i}}+{\tilde{D}}_{22,i}{\mathrm{\Delta y}}_{2,\mathrm{\tau i}},i\mathrm{1,2},...k$  (17e)

由于在[-τ_max,0]时间段内，隐函数定理总成立，则矩阵可逆，其中：

$\tilde{D}=\left(\begin{array}{cc}{\tilde{D}}_{11}& {\tilde{D}}_{12}\\ 0& {\tilde{D}}_{22}\end{array}\right)$  (18)

并令：

${\Omega }_{11}=-{\tilde{D}}_{11}^{-1},{\Omega }_{22}=-{\tilde{D}}_{22}^{-1}$  (19)

${\Omega }_{22,i}=-{\tilde{D}}_{22,i}^{-1},i=\mathrm{1,2},...,k$  (20)

由此可得：

Δy₁＝K₁₁Δx₁+K₁₂Δx₂  (21a)

Δy₂＝K₂₂Δx₂  (21b)

Δy_2,τi＝K_22,iΔx_2,τi＝Ω_22,iC_22,iΔx_2,τi  (21c)

其中：

$K_{11}={\Omega }_{11}{\tilde{C}}_{11}$  (22a)

$K_{12}={\Omega }_{11}({\tilde{C}}_{12}+{\tilde{D}}_{12}{\Omega }_{22}{\tilde{C}}_{22})$  (22b)

$K_{22}={\Omega }_{22}{\tilde{C}}_{22}$  (22c)

3、整理得到系统的带约束时滞微分方程模型（CTODE）如下：

$\Delta {\stackrel{·}{z}}_1=A_{11}{\mathrm{\Delta z}}_1+A_{12}{\mathrm{\Delta z}}_2$(23a)

$\Delta {\stackrel{·}{z}}_2=A_{21}{\mathrm{\Delta z}}_1+A_{22}{\mathrm{\Delta z}}_2+\underset{i=1}{\overset{k}{\Sigma }}{A}_{d,i}{\mathrm{\Delta z}}_{2,\mathrm{\tau i}}$  (23b)

其中：

$A_{11}={\tilde{A}}_{11}+{\tilde{B}}_{11}{K}_{11}$

$A_{12}={\tilde{A}}_{12}+{\tilde{B}}_{11}{K}_{12}+{\tilde{B}}_{12}{K}_{22}$

$A_{21}={\tilde{A}}_{21}+{\tilde{B}}_{21}{K}_{11}$

$A_{22}={\tilde{A}}_{22}+{\tilde{B}}_{21}{K}_{12}+{\tilde{B}}_{22}{K}_{22}$

$A_{d,i}={\tilde{A}}_{d,i}+{\tilde{B}}_{d,i}{K}_{22,i}$

4、建立用于Lyapunov判稳的线性矩阵不等式系统

4.1定义一个标量τ₁，定义对称正定矩阵变量$P_{11}=P_{11}^T,P_{22}=P_{22}^T,Q=Q^T,Z=Z^T,X_{11}=X_{11}^T,$定义任意维数矩阵变量P₁₂、X₁₂、N₁、N₂，式中上标“T”均指矩阵转置；

4.2描述一个线性矩阵不等式系统 P₁₁＞0，P₂₂＞0，X₁₁＞0，X₂₂＞0，Q＞0，Z＞0，$\bar{\Phi }=\left(\begin{array}{ccc}{\bar{\Phi }}_{11}& {\bar{\Phi }}_{12}& {\bar{\Phi }}_{13}\\ {\bar{\Phi }}_{12}^T& {\bar{\Phi }}_{22}& {\bar{\Phi }}_{23}\\ {\bar{\Phi }}_{13}^T& {\bar{\Phi }}_{23}^T& {\bar{\Phi }}_{33}\end{array}\right)<0,$$\bar{\Psi }=\left(\begin{array}{ccc}{X}_{11}& X_{12}& N_1\\ X_{12}^T& X_{22}& N_2\\ N_1^T& N_2^T& Z\end{array}\right)>0$

                                                   (24)

其中：

${\bar{\Phi }}_{11}=P_{11}{A}_{11}+A_{11}^T{P}_{11}+P_{12}{A}_{21}+A_{21}^T{P}_{12}^T+{\tau }_1{A}_{21}^T{\mathrm{ZA}}_{21}$

${\bar{\Phi }}_{12}=P_{11}{A}_{12}+P_{12}{A}_{22}+A_{11}^T{P}_{12}+{A_{21}^T{P}_{22}+\tau }_1{A}_{21}^T{\mathrm{ZA}}_{22}$

${\bar{\Phi }}_{13}=P_{12}{A}_{d,1}+{\tau }_1{A}_{21}^T{\mathrm{ZA}}_{d,1}$

${\bar{\Phi }}_{22}=P_{12}^T{A}_{12}+A_{12}^T{P}_{12}+P_{22}{A}_{22}+A_{22}^T{P}_{22}+Q+N_1$

${+N}_1^T+{\tau }_1{X}_{11}+{\tau }_1{A}_{22}^T{\mathrm{ZA}}_{22}$

${\bar{\Phi }}_{23}=P_{22}{A}_{d,1}-N_1+N_2+{\tau }_1{X}_{12}+{\tau }_1{A}_{22}^T{\mathrm{ZA}}_{d,1}$

${\bar{\Phi }}_{33}=-Q-N_2-N_2^T+{\tau }_1{X}_{22}+{\tau }_1{A}_{d,1}^T{\mathrm{ZA}}_{d,1}$

式中上标“T”均指矩阵转置。

5、对于标量τ₁，利用线性矩阵不等式（LMI）求解工具，判断步骤4.2所给线性不等式 的可行性，如果可行，则当时滞常数为τ₁时，系统(23)在平衡点附近渐进稳定。

三、有益效果

这里比较原有方法（基于TODE模型）和本专利方法（基于CTODE模型）的计算效果。

首先给出原有分析方法要用到的判据，然后采用一通用时滞系统，对新稳定分析方法计 算效率及效果进行验证；进而采用单机无穷大系统作进一步示例和验证。

1、原有方法

原有的基于时滞微分方程模型（TODE）的方法在进行Lyapunov稳定性分析时，对于时 滞系统

$\stackrel{·}{x}=\mathrm{Ax}+\underset{i=1}{\overset{k}{\Sigma }}{A}_{d,i}{x}_{\mathrm{\tau i}}$  (25)

当k=1时，如下定理给出了该系统稳定的条件。

定理[22]：给定标量τ₁>0，若存在对称正定矩阵P＝P^T＞0,Q＝Q^T＞0,Z＝Z^T＞0, 和任意合适维数的矩阵N₁,N₂,X₁₂，且式(26)成立，则时滞系统(25) 在时滞为τ₁时渐进稳定。

$\Phi =\left(\begin{array}{cc}{\Phi }_{11}+{\tau }_1{A}^T\mathrm{ZA}& {\Phi }_{12}+{\tau }_1{A}^{T}Z{A}_1\\ {\Phi }_{12}^T+{\tau }_1{A}_1^T\mathrm{ZA}& {\Phi }_{22}+{\tau }_1{A}_1^{T}Z{A}_1\end{array}\right)<0$  (26a)

$\Psi =\left(\begin{array}{ccc}{X}_{11}& X_{12}& N_1\\ X_{12}^T& X_{22}& N_2\\ N_1^T& N_2^T& Z\end{array}\right)>0$  (26b)

其中：

${\Phi }_{11=}\mathrm{PA}+A^{T}P+N_1+N_1^T+Q+{\tau }_1{X}_{11}$

${\Phi }_{12}={\mathrm{PA}}_1-N_1+N_2^T+{\tau }_1{X}_{12}$

${\Phi }_{22}=-N_2-N_2^T-Q+{\tau }_1{X}_{22}$

2、实施例1——通用时滞系统

考虑某一单时滞系统，其模型如用式(1)所示，相应的系数矩阵由式(27)给出。采用式(28) 所示变换矩阵，并采用式(8)和式(10)变换方法，可得式(12)所示最终结果，其中矩阵A,A_d,1示 于式(29)，此时n＝10,n₁＝8,n₂＝2。

$\tilde{A}=\left(\begin{array}{cccccccccc}-20.00& 0.800& 0.000& -0.900& 0.000& 0.000& 0.000& 1.000& 0.000& -3.000\\ -0.800& -10.60& 0.100& 0.000& -1.000& -2.000& 0.000& 0.900& 0.100& 1.000\\ 0.100& 5.200& -6.000& 0.990& -0.300& 2.000& 1.000& 3.000& -0.300& 0.000\\ 1.250.& 1.230& -0.320& -8.100& -0.100& 0.100& 0.520& 0.200& 0.000& -0.100\\ 1.980& -1.100& -0.008& -0.010& -8.789& 0.000& 0.090& -0.800& 0.200& -1.810\\ -0.100& 0.180& 0.000& -1.100& 0.000& -5.000& 1.000& 0.500& -0.080& 0.000\\ -0.220& 0.320& 3.000& -3.100& -1.000& -1.000& -6.000& 1.900& -2.000& 0.500\\ 0.000& 3.300& 0.500& -3.990& -1.300& 0.000& 1.200& -5.000& 2.500& 2.100\\ -0.150& -0.300& -0.100& 0.800& -0.220& 0.800& 0.200& 0.000& -5.600& 0.100\\ 1.000& 3.200& 1.000& -1.000& -1.000& 0.000& -2.000& 0.000& -3.000& -7.000\end{array}\right)$   (27a)

${\tilde{A}}_1=\left(\begin{array}{cccccccccc}0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000\\ 0.000& 0.000& -6.000& 0.000& 0.000& 0.000& 0.000& 3.000& 0.000& 0.000\\ 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000\\ 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000\\ 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000\\ 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000\\ 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000\\ 0.000& 0.000& 0.500& 0.000& 0.000& 0.000& 0.000& -5.000& 0.000& 0.000\\ 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000\\ 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000& 0.000\end{array}\right)$   (27b)

${T=T}^{-1}=\left(\begin{array}{cccccccccc}1& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 1& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 1\\ 0& 0& 1& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 1& 0& 0\end{array}\right)$   (28)

$A=\left(\begin{array}{cccccccccc}-20.00& 0.800& 0.000& -0.900& 0.000& 0.000& 0.000& -3.000& 0.000& -1.000\\ -0.800& -10.60& 1.000& 0.000& -1.000& -2.000& 0.000& 1.000& 0.100& 0.900\\ -1.500& -0.300& -5.600& 0.800& -0.220& -0.800& 0.200& 0.100& -0.100& 0.000\\ 1.250& 1.230& 0.000& -8.100& -0.100& 0.100& 0.500& -0.100& -0.320& 0.200\\ 1.980& -1.100& 0.200& -0.010& -8.789& 0.000& 0.090& -1.810& -0.008& -0.800\\ -0.100& 0.180& -0.080& -1.100& 0.000& -5.000& 1.000& 0.000& 0.000& 0.500\\ -0.220& 0.320& -2.000& -3.100& -1.000& -1.000& -6.000& 0.500& 3.000& 1.900\\ 1.000& 3.200& 3.000& 1.000& 1.000& 0.000& -2.000& -7.000& 1.000& 0.000\\ -0.100& 5.200& -0.300& -0.990& -0.300& 2.000& 1.000& 0.000& 0.000& 0.000\\ 0.000& 3.300& 2.500& -3.990& -1.300& 0.000& 1.200& 2.100& 0.000& 0.000\end{array}\right)$   (29a)

$A_{d,1}=\left(\begin{array}{cc}-6.000& 3.000\\ 0.500& -5.000\end{array}\right)$   (29b)

（1）两种稳定分析方法的计算结果比较

首先利用[10]方法求解系统的实际时滞稳定裕度为229.911ms。再分别用式(26)和式(13) 两种方法进行判稳，所得结果示于表1。可以看到，采用新方法，可得到与原有方法相近的 时滞稳定裕度，彼此误差不足0.2%。但新方法的计算用时仅为原有方法的1.85%，具有更高 的计算效率。为解释新的方法为什么具有更高的计算效率，我们分析两种方法所需确定的待 求变量数目：当采用判据式(26)时，系统待求的变量矩阵分别为P,Q,Z,X₁₁,X₂₂,X₁₂,N₁,N₂，这 些矩阵均为n×n方阵，前5个为对称矩阵，后三个为非对称方阵。则对于本算例系统，待求 变量数为：

$N_{\mathrm{TODE}}=5\times \frac{(n+1)\times n}{2}+\times n^2=575$

而对于式(13)的新方法，其待求变量矩阵仍为：P,Q,Z,X₁₁,X₂₂,X₁₂,N₁,N₂。但其中除矩阵P与 原判据(26)中待求变量P完全相同外，这里：

$P=\left(\begin{array}{cc}{P}_{11}& P_{22}\\ P_{12}^T& P_{22}\end{array}\right)$

为n×n对称矩阵，其他矩阵均为n₂×n₂方阵。因此新方法对应的待求变量数变为：

$N_{\mathrm{CTODE}}=\frac{(n+1)\times n}{2}+4\times \frac{(n_2+1)\times n_2}{2}+3\times {\left(n_2\right)}^2=79$

很明显，由于新方法的待求变量数较原方法大为减少，其具有更高计算效率就不难理解了。

表1两种方法计算结果比较

Tab.1Comparison between two criteria

（2）两种稳定方法计算效率比较

为进一步研究两种方法的计算效率与计算效果随系统规模的变化规律，首先保持矩阵A_d,1不变，按如下方式对矩阵A加以变动，以形成不同维数的时滞系统：取3≤r≤10，并由矩阵A 的右下角r×r方阵形成矩阵A_r，进而构建如下时滞系统系列：

$\stackrel{·}{x}=A_{r}x+\underset{i=1}{\overset{k}{\Sigma }}{A}_{r,i}{x}_{\mathrm{\tau i}},r=3,...,10$  (30)

其中：$A_{r,i}=\left(\begin{array}{cc}0& 0\\ 0& A_{d,i}\end{array}\right)\in R^{r\times r}$  (31)

当r＝8时，A_r矩阵即对应于式(29a)中的阴影部分，此时有n＝r＝8,n₁＝r-2＝6,n₂＝2。对 由式(30)形成的时滞系统系列，分别采用两种方法进行判稳，计算结果示于表2。图1绘出了 新方法加速比随系统维数的变化曲线。可以看出，随着系统维数的增加，新方法的加速效果 更为明显。

表2两种方法计算结果比较(不同时滞系统)

Tab.2Comparison between two criteria with different dimension time-delay systems

进一步，假设A_d,1维数不变，图2绘出了当矩阵A的维数n从1增至100时，两种方法所 需求解的未知变量数目的变化情况。不难看出，随着维数n的增加，新方法与原有方法之间 所需求解变量数目的差距不断增大，由此决定新的方法在高维情况下较原方法具有更好的计 算效率。

（3）交叉项对新判据的影响

观察式(29)我们可以看出，对于该系统，实际是假设在广域控制器设计时，全部使用远 程控制信号(由A_d,1的四个元素决定，它们需要考虑信号传输的延时)进行协调控制。但在实际 应用中，协调控制器往往需要协调远程量测信息和本地量测信息，为此引入交叉因子α，并令：

A_α＝A+α·A_i

A_α,i＝(1-α)·A_i

由此形成新的时滞系统系列：

$\stackrel{·}{x}=A_{\alpha }x+\underset{i=1}{\overset{k}{\Sigma }}{A}_{\alpha ,i}{x}_{\mathrm{\tau i}}$

其中，矩阵A_i由式(11)给出。当α=0时，即为式(29)所给原系统，此时协调控制信号只取远方 信息；当α=1时，则控制器只取本地信息，将不再考虑控制回路的信号延时的影响。

我们令α在0～0.60之间变动，表3给出了α取不同数值时两种方法的计算结果，图3给出 了两种方法误差随α变动的情况。从这些结果可以看出如下规律：两种方法计算结果非常接近， 误差基本在5%以内，且多数误差在1%以下；两种方法计算误差随α变动的趋势一致，一种 方法误差随α增大而变大时，另一种误差也将变大，反之亦然。

表3计算结果随交叉系数α的变化情况

Tab.3Results of new criterion with differentαvalues

3、实施例2——单机无穷大系统

单机无穷大系统及其参数的取值与文献[4]相同，表4给出了基于原有TODE模型稳定判 据(原方法)和基于CTODE模型稳定判据(新方法)分别求解系统时滞稳定裕度的计算结果，表 中的加速比p定义为：

$p=\frac{t_{\mathrm{TODE}}-t_{\mathrm{CTODE}}}{t_{\mathrm{TODE}}}\times 100\%$

其中，t_TODE,t_CTODE为采用两种模型求解时滞稳定裕度时所用的时间。从中可以看出如下规律：

1)P_m在1.00-1.04之间变动时，原方法计算结果的保守性要好于新方法；而P_m在1.05-1.13 之间变动时，新方法计算结果的保守性要优于原方法。

2)采用新方法的计算用时更少，不同的计算场景可实现8.47%-32.6%的加速比。但与表 1计算结果相比，本算例计算速度提升很有限，其原因可用下式(32)-式(34)来加以解释：它们 给出Pm=1.0时，系统的矩阵和变换矩阵T。

$\tilde{A}=\left(\begin{array}{cccc}0.0000& 376.9911& 0.0000& 0.0000\\ -0.0963& -0.5000& -0.0801& 0.0000\\ -0.0480& 0.0000& -0.1667& 0.1000\\ 0.0000& 0.0000& 0.0000& -1.0000\end{array}\right)$  (32a)

${\tilde{A}}_1=\left(\begin{array}{cccc}0.0000& 0.0000& 0.0000& 0.0000\\ 0.0000& 0.0000& 0.0000& 0.0000\\ 0.0000& 0.0000& 0.0000& 0.0000\\ 38.0187& 0.0000& -95.2560& 0.0000\end{array}\right)$  (32b)

$A=\left(\begin{array}{cccc}-0.5000& -0.0963& 0.0801& 0.0000\\ 376.9911& 0.0000& 0.0000& 0.0000\\ 0.0000& -0.0480& -0.1667& 0.1000\\ 0.0000& 0.0000& 0.0000& -1.0000\end{array}\right)$  (33a)

  (33b)

$T=T^{-1}\left(\begin{array}{cccc}0& 1& 0& 0\\ 1& 0& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right)$   从式(33b)可看出，变换后得到的A_d,1矩阵为3维方阵，则n=4,n₁=1,n₂=3。因此，采用原方法 和新方法求解时待求变量数分别为：

$N_1=5\times \frac{(n+1)\times n}{2}+3\times n^2=98$

$N_2=\frac{(n+1)\times n}{2}+4\times \frac{(n_2+1)\times n_2}{2}+3\times {\left(n_2\right)}^2=61$

表4单机无穷大系统计算结果

从上述计算结果可以看出，基于CTODE模型的Lyapunov-Krasovskii稳定分析方法，与 原有方法具有相似的计算精度，却具有更高的计算效率。

在电力系统广域测量平台(WAMS)中，可借助本发明方法，对电力系统的时滞安全性进 行快速判别。

参考文献

1.K.Q Gu,V.L.Kharitonov,Chen Jie,Stability of time-delay systems[M],Boston,2003.

2.M.Lakshmanan,D.V.Senthilkumar,Dynamics of nonlinear time-delay systems[M],Springer,Berlin,2011.

3.H.X.Wu,K.S.Tsakalis and G.T.Heydt,Evaluation of time delay effects to wide-area power system stabilizer  design[J],IEEE Trans.on Power Systems,2004,19(4):1935-1941.

4.Jia Hongjie,Yu Xiaodan,Yu Yixin and Wang Chengshan,Power system small signal stability region with time  delay[J],Int.Journal of Electrical Power&Energy Systems,2008,30(1):16-22.

5.Jia Hongjie,Yu Xiaodan,Wang Chengshan,et al.Study on power system extended small signal stability  region(DE-SSSR)in time delay space,Proc.of15th IEEE MELECON,2010.4.25-4.28,Valletta,Malta,Spain, Vol.1:1569-1574.

6.A.Saffet,Computation of time delay margin for power system small-signal stability,European Trans.on  Electrical Power,2009,19(7):949-968.

7.贾宏杰，谢星星和余晓丹，考虑时滞影响的电力系统小扰动稳定域，电力系统自动化，2006,30(21): pp.1-5.

Jia Hongjie,Xie Xingxing,Yu Xiaodan,Power system small signal stability region with time delay considered, Automation of Electric Power Systems,2006,30(21):1-5.

8.Nejat O,Rifat S,An exact method for the stability analysis of time-delayed linear time-invariant(LTI)systems, IEEE Trans.on Automatic Control,2002,47(5):793-797.

9.Jia Hongjie,Cao Xiaodong,Yu Xiaodan,Zhang Pei,A simple approach to determine power system delay  margin,Proc.of IEEE General Meeting,2007.6.24-6.28,Tampa,Florida.Vol.1-10:799-805.

10.Jia Hongjie,Yu Xiaodan,A simple method for power system stability analysis with multiple time delays,Proc. of IEEE General Meeting,2008.7.20-7.24,Pittsburgh,PA USA,Vol.1-11:4004-4010.

11.A.Saffet,C.O.Nwankpa,Probabilistic evaluation of small-signal stability of power systems with time delays, Int.Review of Electrical Engineering,Part B,2010,5(1):205-216.

12.An Haiyun,Jia Hongjie,Yu Xiaodan,An improved delay-dependent robust stability criterion and application  to power system stability analysis with time delays,Proc of IEEE POWERCON,2010.10.24-10.28,Hang  Zhou,China,Paper No:5666129

13.A.M.Hasan,P.Minwon,Y.In-Keun,Enhancement of transient stability by fuzzy logic-controlled SMES  considering communication delay,Int.Journal of Electrical Power & Energy Systems,2009,31(7-8):402-408.

14.Lu Chao,Han Yingduo,Wu Xiaochen,Field experiments of wide area damping controllers for multiple  HVDC links,Proc of IEEE Asia Pacific Conference on Circuits and Systems,2008.11.30-12.3,Macao,China: 627-630.

15.M.Jankovic,Control Lyapunov-Razumikhin functions and robust stabilization of time delay systems,IEEE  Trans.on Automatic Control,2001,46(7):1048-1060.

16.S.Dashkovskiy,H.R.Karimi,M.Kosmykov,A Lyapunov-Razumikhin approach for stability analysis of  logistics networks with time-delays,International Journal of Systems Science,2012,43(5):845-853.

17.E.R Mahboobi,S.K.Y.Nikravesh,Stabilising predictive control of non-linear time-delay systems using  control Lyapunov-Krasovskii functionals,IET Control Theory&Applications,2009,3(10):1395-1400.

18.E.Fridman,Stability of systems with uncertain delays:a new complete Lyapunov-Krasovskii functional,IEEE Trans.on Automatic Control,2006,51(5):885-890.

19.Wu Min,He Yong and She Jinhua,Stability analysis and robust control of time-delay systems,Science Press, Beijing,2009.

20.李鹏，贺静波，石景海等．交直流并联大电网广域阻尼控制技术理论与实践[J]，南方电网技术，2008, 2(4):13-17.

21.C.W.Taylor,The future in on-line security assessment and wide-area stability control.Proc.of IEEE PES  Winter Meeting,2000.1.23-1.27,Singapore,Vol.1:78-83.

22.B.Chaudhuri,R.Majumder,B.Pal,Wide area measurement based stabilizing control of power system  considering signal transmission delay,Proc.of IEEE PES General Meeting,2005.6.12-6.16,San Francisco, CA,Vol.2:1447-1452.

23.J.K.Hale and S.M.V.Lunel,Introduction to functional differential equations[M].Springer-Verlag,New York, 1993.

24.Liu Fang,Yokoyama Ryuichi,Zhou Yicheng,et al,Free-weighting matrices-based robust wide-area FACTS  control design with considering signal time delay for stability enhancement of power systems,IEEJ Trans.on  Electrical & Electronic Engineering,2012,7(1):31-39.