文章目录
文献学习笔记： Olgac N , Sipahi R . An exact method for the stability analysis of time-delayed linear time-invariant (LTI) systems[J]. IEEE Transactions on Automatic Control, 2002, 47(5):793-797.
1 研究的对象以及方法简述
我们研究的对象是一个线性时不变的时滞系统，分析它在不同时滞情况下的稳定性：
x
˙
=
A
x
+
B
x
(
t
−
τ
)
(1.1)
\dot x=Ax+Bx(t-\tau) \tag{1.1}
x˙=Ax+Bx(t−τ)(1.1)
对上述系统求拉普拉斯变换并求行列式等于0，得到系统的特征方程
∣
s
I
−
A
−
B
e
−
τ
s
∣
=
0
(1.2)
|sI-A-Be^{-\tau s}|=0 \tag{1.2}
∣sI−A−Be−τs∣=0(1.2)
系统能够收敛也就意味着特征根的实部都要位于复平面的左半平面，或者说特征根实部小于零。显然对于越高维的矩阵，时滞引入的指数项会大大提高求解难度，本文提出了一种有效的方式来简化求解过程，分析出使系统能够收敛的时滞区间。
方法简述： 时滞为0的情况下，求解系统的特征根是很简单的。系统通过定义root tendency来判断每一个时滞下特征根的移动，来判断每个时滞下是否出现有根向左/右半平面的移动，以此计算根的位置分布。
2 新定义的变量介绍
2.1 Crossing frequency
crossing frequency的叫法似乎并不是在这一篇论文中提出来的，我在读另一篇使用该方法的论文中，另一篇论文的作者称呼这个变量为crossing frequency，我感觉这个称呼非常贴切，故延用。
系统的特征方程(1.2)展开表示如下：
∑
k
=
0
n
a
k
(
s
)
⋅
e
−
k
τ
s
=
0
(2.1.1)
\sum_{k=0}^na_k(s)\cdot e^{-k\tau s}=0 \tag{2.1.1}
k=0∑nak(s)⋅e−kτs=0(2.1.1)
其中
a
k
(
s
)
a_k(s)
ak(s)是关于s的多项式，最高次数为n-k次。
可以证明式(2.1.1)对于任意的时滞
τ
\tau
τ，产生的虚根只会在有限多个虚根中取值，这里设虚根数量为m，虚根系数集合为
{
ω
c
}
=
{
ω
c
1
,
ω
c
2
,
…
,
ω
c
m
}
\{\omega_c\}=\{ \omega_{c1}, \omega_{c2}, \dots, \omega_{cm}\}
{ωc}={ωc1,ωc2,…,ωcm}，
其实严格来说应该是产生小于等于2m个特征根，因为对于一个非0的crossing frequency它所产生的一定是一对共轭虚根，但是为了方便表述，下面所谓的特征根我都特指共轭虚根中的一个
其下标c表示crossing，意为穿越虚轴，在后面的内容中我们会对这个穿越有更深刻且直观的体会。
这一系列的角频率被称为crossing frequency。
2.2 crossing frequency生成的时滞集合
时滞在特征方程中只位于指数项，在复平面上看，完全可以把这个指数当做一个单位向量，一旦i的系数经过一个
2
π
2\pi
2π，就会回到原来的位置。
可以想到，如果有一组参数
{
τ
k
,
ω
c
k
}
\{\tau_k,\omega_{ck}\}
{τk,ωck}使式(1.2)成立，那么必然还会有无限个等间距的时滞也能够使式(1.2)成立，显然这个间距为
2
π
/
ω
k
2\pi/\omega_k
2π/ωk。
所以每一个crossing frequency可以生成一套无限的时滞集合，其中任意的元素与生成它的crossing frequency能够满足式(1.2)
ω
c
k
⟹
{
τ
k
}
=
{
τ
k
1
,
τ
k
2
,
⋯
,
τ
k
∞
}
(2.2.1)
\omega_{ck}\Longrightarrow\{\tau_k\}=\{\tau_{k1}, \tau_{k2}, \cdots, \tau_{k\infty}\} \tag{2.2.1}
ωck⟹{τk}={τk1,τk2,⋯,τk∞}(2.2.1)
不过在使用时我们最关心的还是最小的大于零的
τ
\tau
τ，我们也通常取这个
τ
\tau
τ为上式中的
τ
k
1
\tau_{k1}
τk1，以此作为序列起点。
2.3 root sensitivities
至此我们已经知道了对于上述系统，任意一个时滞下会存在 m个特征根，他们分布在复平面的很多地方，由于我们已经知道了0时滞状态下的特征根分布，所以对于其他的时滞，我们不再关心特征根的新分布，而关心时滞由0到当下的时滞过程中，特征根有多少从左边移动到了右边，或者从右边移动到了左边。
也就是说，对于任意一个时滞，我们关心的其实是crossing frequency对应的m个虚根，它代表特征方程的m个特征根中，有小于等于m个正好处于虚轴上。
用我很喜欢的一个比喻，那就是m个特征根在复平面上到处周期性移动，留下m条闭合轨迹，轨迹会穿越虚轴，留下m个“穿越点”，当特征方程解出几个虚根的时候，就代表有几个特征根正好移动到虚轴上，那几个“穿越点”此时被“点亮”了。
我们非常关心这一个时滞
τ
n
o
w
\tau_{now}
τnow下，这m个虚根关于时滞
τ
\tau
τ在
τ
n
o
w
\tau_{now}
τnow下的导数值，这个导数值如果为正，那么紧接着这个根就会跑到右半平面，使系统更发散一点，这个导数值如果为负，那么紧接着这个根就会跑到左半平面，使系统更收敛一点。 这里的更发散和更收敛是一种形象化的表述，实际上是否收敛只取决于特征根是否全部都在左半平面。
root sensitivities就被定义为crossing frequency对应虚根对产生时的时滞的导数，即：
S
τ
s
∣
s
=
ω
c
k
i
,
τ
=
τ
k
l
=
d
s
d
τ
∣
s
=
ω
c
k
i
,
τ
=
τ
k
l
(2.3.1)
S^s_\tau|_{s=\omega_{ck}i,\ \tau=\tau_{kl}} = \frac{ds}{d\tau}|_{s=\omega_{ck}i,\ \tau=\tau_{kl}} \tag{2.3.1}
Sτs∣s=ωcki, τ=τkl=dτds∣s=ωcki, τ=τkl(2.3.1)
这里一定要理解一个很重要的概念！任意一个crossing frequency可以生成一组无限的时滞，而任意一个时滞也可以生成m个crossing frequency，这是理解上式导数值下标的关键。
在任意一个时滞下，产生了小于等于m个crossing frequency，同时，这一个时滞也成为了这几个crossing frequency的生成时滞集中的元素。因此，只有上式中下标的时滞是crossing frequency生成的时候，root sensitivities才有讨论的意义。这里的意义我想了很久，也正是想通的一瞬间完全弄明白了全文，它真的很关键，也真的很难用文字描述清楚
2.4 root tendency
root sensitivities本质上说明了位于虚轴临界点的特征根下一瞬间会向着哪里移动，它给出了一个二维方向，当然我们并不关心这个特征根会向虚轴正方向还是负方向移动，管他呢，我们只需要了解它下一时刻跑到了左半平面还是负半平面就好了，换句话说，我们只关心root sensitivities的实部是朝着正方向还是负方向。
定义root tendency：
R
T
∣
s
=
ω
c
k
i
=
s
g
n
[
R
e
(
S
τ
s
∣
s
=
ω
c
k
i
,
τ
=
τ
k
l
)
]
(2.4.1)
RT|_{s=\omega_{ck}i}=sgn[Re(S^s_\tau|_{s=\omega_{ck}i,\ \tau=\tau_{kl}} )] \tag{2.4.1}
RT∣s=ωcki=sgn[Re(Sτs∣s=ωcki, τ=τkl)](2.4.1)
它的意义相当明确，取root sensitivities的实部并取符号，显然我们得到以下结论：
{
R
T
∣
s
=
ω
c
k
i
=
+
1
,
ω
c
k
i
will move to right half plane
R
T
∣
s
=
ω
c
k
i
=
−
1
,
ω
c
k
i
will move to left half plane
(2.4.2)
\left\{ \begin{array}{c} RT|_{s=\omega_{ck}i}=+1, \omega_{ck}i\ \text{will move to right half plane} \\ RT|_{s=\omega_{ck}i}=-1, \omega_{ck}i\ \text{will move to left half plane} \end{array} \right. \tag{2.4.2}
{RT∣s=ωcki=+1,ωcki will move to right half planeRT∣s=ωcki=−1,ωcki will move to left half plane(2.4.2)
如果root tendency的性质仅仅止步于此的话，你一定会觉得它的作用仅仅只是标记一下特征根的移动方向罢了，但其实它最重要的性质会在3中介绍，这将直接决定了这种稳定性分析方法的优越性！
2.5 T（它仅仅是一个代换用的变量）
文献：Z. V. Rekasius, “A stability test for systems with delays,” in Proc. Joint
Automatic Control Conf., 1980, Paper No. TP9-A
上面这篇论文中介绍了一种代换方式，可以帮助解决特征方程中指数项难解的问题：
e
−
τ
s
=
1
−
T
s
1
+
T
s
,
τ
∈
R
+
,
T
∈
R
(2.5.1)
e^{-\tau s}=\frac{1-Ts}{1+Ts},\ \tau \in \mathbb{R^+},\ T\in \mathbb{R} \tag{2.5.1}
e−τs=1+Ts1−Ts, τ∈R+, T∈R(2.5.1)
(2.5.1)式在s为纯虚根的时候是左右严格相等的，且使得
τ
\tau
τ和
T
T
T有如下的代换关系：
τ
k
l
=
2
ω
c
k
[
t
a
n
−
1
(
ω
c
k
T
)
±
l
π
]
(2.5.2)
\tau_{kl}=\frac{2}{\omega_{ck}}[tan^{-1}(\omega_{ck}T)\pm l\pi] \tag{2.5.2}
τkl=ωck2[tan−1(ωckT)±lπ](2.5.2)
由于我们希望求解的是特征方程的纯虚根，所以可以将式(2.5.1)代入到(2.1.1)中，将指数项转化成多项式的形式，我们进行如下的化简：
∑
k
=
0
n
a
k
(
s
)
(
1
−
T
s
1
+
T
s
)
k
=
0
∑
k
=
0
n
a
k
(
s
)
(
1
+
T
s
)
n
−
k
(
1
−
T
s
)
k
=
0
∑
k
=
0
2
n
b
k
(
T
)
⋅
s
k
=
0
(2.5.3)
\sum_{k=0}^na_k(s)(\frac{1-Ts}{1+Ts})^k=0 \tag{2.5.3} \\ \sum_{k=0}^na_k(s)(1+Ts)^{n-k}(1-Ts)^k=0 \\ \sum_{k=0}^{2n}b_k(T)\cdot s^k=0
k=0∑nak(s)(1+Ts1−Ts)k=0k=0∑nak(s)(1+Ts)n−k(1−Ts)k=0k=0∑2nbk(T)⋅sk=0(2.5.3)
经过化简，关于s的多项式最高次数变成了2n次，且系数b是关于T的多项式，简化了运算。论文中给出结论，T与crossing frequency是一一对应的关系，具体证明可以参考论文原文。
3 关键的性质
由2中介绍的变量，我们已经清楚了在时滞变化的过程中，一共只有m个crossing frequency产生的特征根会发生穿越，并且发生穿越对应的时滞是可以通过累加一个固定的间隔
2
π
/
ω
c
k
2\pi /\omega_{ck}
2π/ωck预测出来的，这样我们就已经能够预测在哪些时滞下会发生根的穿越，导致稳定性的变化。
但是接下来要介绍一个更为关键的性质，那就是RT（root tendency）仅仅只跟crossing frequency有关，与时滞无关。
这个性质表明我们不光可以预测在哪些时滞下会发生根的穿越，还能预测到穿越的方向！
3.1 Root tendency与时滞无关
由式(2.1.1)两边微分化简可以得到：
d
s
d
τ
∣
s
=
ω
c
k
i
,
τ
=
τ
k
l
=
∑
j
=
0
n
−
j
s
a
j
e
−
j
τ
s
∑
j
=
0
n
[
d
a
j
d
s
−
j
a
j
τ
]
e
−
j
τ
s
∣
s
=
ω
c
k
i
,
τ
=
τ
k
l
(3.1.1)
\frac{ds}{d\tau}|_{s=\omega_{ck}i,\tau=\tau_{kl}}=\frac {\sum_{j=0}^n-jsa_je^{-j\tau s}} {\sum_{j=0}^n[\frac{da_j}{ds}-ja_j\tau]e^{-j\tau s}}|_{s=\omega_{ck}i,\tau=\tau_{kl}} \tag{3.1.1}
dτds∣s=ωcki,τ=τkl=∑j=0n[dsdaj−jajτ]e−jτs∑j=0n−jsaje−jτs∣s=ωcki,τ=τkl(3.1.1)
由于取s=
ω
c
k
i
\omega_{ck}i
ωcki，所以有：
e
−
j
τ
s
∣
s
=
ω
c
k
i
,
τ
=
τ
k
l
=
(
1
−
T
s
1
+
T
s
∣
s
=
ω
c
k
i
,
τ
=
τ
k
l
)
j
(3.1.2)
e^{-j\tau s}|_{s=\omega_{ck}i,\tau=\tau_{kl}} =(\frac{1-Ts}{1+Ts}|_{s=\omega_{ck}i,\tau=\tau_{kl}} )^j\tag{3.1.2}
e−jτs∣s=ωcki,τ=τkl=(1+Ts1−Ts∣s=ωcki,τ=τkl)j(3.1.2)
改写(3.1.1)为：
d
s
d
τ
∣
s
=
ω
c
k
i
,
τ
=
τ
k
l
=
s
∑
j
=
0
n
d
a
j
d
s
e
−
j
τ
s
∑
j
=
0
n
j
a
j
e
−
j
τ
s
−
τ
∣
s
=
ω
c
k
i
,
τ
=
τ
k
l
(3.1.3)
\frac{ds}{d\tau}|_{s=\omega_{ck}i,\tau=\tau_{kl}}=\frac {s} {\frac{\sum_{j=0}^n\frac{da_j}{ds}e^{-j\tau s}}{\sum_{j=0}^nja_je^{-j\tau s}} -\tau}|_{s=\omega_{ck}i,\tau=\tau_{kl}} \tag{3.1.3}
dτds∣s=ωcki,τ=τkl=∑j=0njaje−jτs∑j=0ndsdaje−jτs−τs∣s=ωcki,τ=τkl(3.1.3)
由此可以求出RT：
R
T
k
=
s
g
n
[
R
e
(
d
s
d
τ
)
]
=
s
g
n
[
I
m
(
∑
j
=
0
n
d
a
j
d
s
e
−
j
τ
s
/
∑
j
=
0
n
j
a
j
e
−
j
τ
s
)
]
(3.1.4)
RT_k=sgn[Re(\frac{ds}{d\tau})]\\ = sgn[Im(\sum_{j=0}^n\frac{da_j}{ds}e^{-j\tau s}/\sum_{j=0}^nja_je^{-j\tau s})] \tag{3.1.4}
RTk=sgn[Re(dτds)]=sgn[Im(j=0∑ndsdaje−jτs/j=0∑njaje−jτs)](3.1.4)
可见已经与时滞项无关，仅仅取决与crossing frequency，也就是说每一个特征根会朝哪个方向穿越虚轴在确定时滞之前就已经确定好了。
3.2 使用Routh-Hurwitz判据确定T， ω \omega ω, τ \tau τ
看式(2.5.3)，特征方程已经被转化成了关于s的多项式，且系数是关于T的多项式，根据Routh-Hurwitz判据，Routh表的第一列也全都是关于T的多项式表达式，我们知道第一列值发生多少次符号变化，就代表了特征方程有多少个根位于左半平面。我们记某一个T下，Routh表第一列发生的符号变化次数为NS（Number of Sign changes）。
论文中证明了存在一共m个T，会使得NS发生改变，所以我们需要从一个范围内扫过所有可能的T，看NS在哪些点发生了变化，由这m个T就可以代回到(2.5.3)计算出m个crossing frequency。我在看论文的时候在这里也非常懵逼，那我应该扫多大的T的区间？如果区间过大会使计算量大大提高！怎么办？论文示例中只把T从-0.6扫到了0.6。但是由于Routh表第一列全是关于T的多项式表达，对于常用的时滞模型不会产生过于高次的T，还是可以通过多项式的结构预先判断出T在那个范围外就不会发生NS变化了
所以说计算过程我们使用Routh表计算出m个T，再带回(2.5.3)计算出m个 ω c \omega_c ωc，再利用(2.5.2)计算出每个对应的时滞集合即可。
4 分析流程
- 建模，建立线性时不变的时滞系统(1.1)
- 计算特征方程(2.1.1)，并且代入(2.5.1)得到(2.5.3)
- 利用Routh表计算出 T c k T_{ck} Tck，注意，如果T=0时，NS变化量为n，那么这个T舍去不考虑
- T c k T_{ck} Tck代入(2.5.3)计算出 ω c k \omega_{ck} ωck
- T c k T_{ck} Tck与 ω c k \omega_{ck} ωck代入(2.5.2)计算出 { τ k } \{\tau_k\} {τk}，其实关键地，只需要计算出大于零的最小时滞 τ k 1 \tau_{k1} τk1与时滞间距 Δ τ k = 2 π / ω c k \Delta\tau_k = 2\pi /\omega_{ck} Δτk=2π/ωck
- 利用(3.1.4)计算出 R T k RT_k RTk
- 记系统位于右半平面的根的数量为NU(
τ
\tau
τ)（Number of Unstable roots），NU(0)只需要代入时滞为零就可以轻松解出，那么任意时刻
τ
\tau
τ的NU可以表示为
N U ( τ ) = N U ( 0 ) + ∑ k = 1 m ( τ − τ k 1 Δ τ k ) ⋅ U ( τ , τ k 1 ) ⋅ R T k (4.1) NU(\tau)=NU(0)+\sum_{k=1}^m(\frac{\tau-\tau_{k1}}{\Delta\tau_k})\cdot U(\tau,\tau_{k1})\cdot RT_k \tag{4.1} NU(τ)=NU(0)+k=1∑m(Δτkτ−τk1)⋅U(τ,τk1)⋅RTk(4.1)
其中 U ( τ , τ k 1 ) U(\tau,\tau_{k1}) U(τ,τk1)是一个用来判断特征根穿越虚轴的时候导致了多少根移动的函数，举个例子 ，如果crossing frequency=0，那么一次RT=1的穿越会导致NU增加一个，如果crossing frequency!=0，那么一次RT=1的穿越其实是两个共轭虚根的穿越，会导致NU增加两个。于是我们这样定义它：
U ( τ , τ k 1 ) = { 0 0 < τ < τ k 1 1 τ > τ k 1 , ω c k = 0 2 τ > τ k 1 , ω c k ≠ 0 (4.2) U(\tau,\tau_{k1})= \left\{\begin{array}{c} 0 & 0<\tau<\tau_{k1} \\ 1 & \tau>\tau_{k1},\ \omega_{ck}=0\\ 2 & \tau>\tau_{k1},\ \omega_{ck}\ne0 \end{array}\right. \tag{4.2} U(τ,τk1)=⎩⎨⎧0120<τ<τk1τ>τk1, ωck=0τ>τk1, ωck=0(4.2)
于是(4.1)式的意义就十分清晰了，从 τ = 0 \tau=0 τ=0开始计数，每过一个周期 Δ τ k \Delta\tau_{k} Δτk，会发生一次NU变化，统计从 τ = 0 \tau=0 τ=0到现在的所有变化，显然当NU=0时系统可以收敛，当NU ≠ \ne = 0时，系统不会收敛。
5 实例（均取自论文实例）
我们定义(1.1)中的矩阵AB分别为：
A
=
[
−
1
13.5
−
1
−
3
−
1
−
2
−
2
−
1
−
4
]
,
B
=
[
−
5.9
7.1
−
70.3
2
−
1
5
2
0
6
]
A=\left[ \begin{matrix} -1&13.5&-1 \\ -3&-1&-2 \\ -2&-1&-4\\ \end{matrix} \right] \ , B=\left[ \begin{matrix} -5.9&7.1&-70.3 \\ 2&-1&5 \\ 2&0&6\\ \end{matrix} \right]
A=⎣⎡−1−3−213.5−1−1−1−2−4⎦⎤ ,B=⎣⎡−5.9227.1−10−70.356⎦⎤
由此计算特征方程(2.1.1)，得到：
a
3
(
s
)
=
119.4
a
2
(
s
)
=
90.9
s
−
185.1
a
1
(
s
)
=
0.9
s
2
−
116.8
s
−
22.1
a
0
(
s
)
=
s
3
+
6
s
2
+
45.5
s
+
111.0
a_3(s)=119.4 \\ a_2(s)=90.9s - 185.1 \\ a_1(s)=0.9s^2-116.8s-22.1 \\ a_0(s)=s^3+6s^2+45.5s+111.0 \\
a3(s)=119.4a2(s)=90.9s−185.1a1(s)=0.9s2−116.8s−22.1a0(s)=s3+6s2+45.5s+111.0
代入Rekasius的代换公式，得到：
b
6
(
T
)
=
T
3
b
5
(
T
)
=
5.1
T
3
+
3
T
2
b
4
(
T
)
=
253.2
T
3
+
17.1
T
2
+
3
T
b
3
(
T
)
=
−
171.4
T
3
−
162.4
T
2
+
18.9
T
+
1
b
2
(
T
)
=
898.4
T
2
−
71.2
T
+
6.9
b
1
(
T
)
=
137.8
T
+
19.6
b
0
(
T
)
=
23.2
b_6(T)= T^3\\ b_5(T)= 5.1T^3+3T^2\\ b_4(T)= 253.2T^3+17.1T^2+3T\\ b_3(T)= -171.4T^3-162.4T^2+18.9T+1\\ b_2(T)= 898.4T^2-71.2T+6.9\\ b_1(T)= 137.8T+19.6\\ b_0(T)= 23.2\\
b6(T)=T3b5(T)=5.1T3+3T2b4(T)=253.2T3+17.1T2+3Tb3(T)=−171.4T3−162.4T2+18.9T+1b2(T)=898.4T2−71.2T+6.9b1(T)=137.8T+19.6b0(T)=23.2
改变T，从Routh表第一列计算NS变化可以得到下图：
一共有六处发生了NS的变化，但是注意，T=0时，NS变化了n=3，所以这个T我们舍去不考虑。
以此计算crossing frequency与时滞，RT和NU，我们可以得到如下表格：
| τ \tau τ | RT | ω \omega ω | T | NU |
|---|---|---|---|---|
| 0.1624 | 1 | 3.0347 | 0.0829 | 0 |
| 0.1859 | -1 | 2.9123 | 0.0953 | 2 |
| 0.2220 | 1 | 15.5032 | -0.4269 | 0 |
| 0.6273 | 1 | 15.5032 | -0.4269 | 2 |
| 0.8725 | 1 | 2.1109 | 0.6233 | 4 |
| 7.2070 | -1 | 0.8409 | -0.1332 | 6 |
| 40 |
注意，第i行的NU表示经过第i-1行之后NU的数量
那么由表格可以清晰看到
系统在：
0.1624
<
τ
<
0.1859
0.2220
<
τ
<
0.6273
0.1624<\tau<0.1859\\ 0.2220<\tau<0.6273
0.1624<τ<0.18590.2220<τ<0.6273
范围内是可以收敛的。
那么这个时候你可能会问了，看这个势头万一在时滞很大的时候又出现了NU=0怎么办？现在仅仅只讨论到时滞为7.2s的情况。其实很简单，如果时滞都达到了好几秒种，咱们现在就不应该讨论时滞系统了，而应该看看通讯设备是不是已经损坏，或者考虑换新的通讯方案了。
6 MATLAB实现代码
可以参考我的github，我的MATLAB代码写得并不好，只能说是可以完成代码功能而已。如有需要可以自取。
7 一些局限
可以看得出这个方法仍然是有些局限的，首先计算方面我们仍然相当于遍历一个区间内的T进行条件判断来筛选，这种数值计算会非常缓慢，我把精度设置到小数点后四位，需要计算几分钟才能出结果，并且总是存在精度限制，那我们就没办法确定在非常多小数点位数以后有没有可能存在一个解。
但是Olgac在2006年的一篇新的论文中提出了改进的办法，消除了这些局限。Olgac N , Sipahi R . An Improved Procedure in Detecting the Stability Robustness of Systems With Uncertain Delay[J]. IEEE Transactions on Automatic Control, 2006, 51(7):1164-1165.
留到下一次更新。