如何理解李雅普诺夫稳定性分析

J Pan

航空工程师

一、什么是平衡点
二、什么是李雅普诺夫稳定
三、李雅普诺夫第一法
四、李雅普诺夫第二法
五、MATLAB代码

一、什么是平衡点

一个控制系统就和一个社会一样，稳定性是首先要解决的重要问题，是其他一切工作的基础。稳定性问题的字面意思很好理解了，那就是系统在受到扰动后，能否能有能力在平衡态继续工作。大家都知道，历史上社会改革成本很高，且以失败者居多，从控制论的角度来看，就是对社会这个大系统的稳定性研究不够，导致扰动发生后，社会发散了。

稳定性是相对于平衡点而言的，那什么是平衡点呢？我们以发射火箭为例：

火箭的简化模型可以看成是一个倒立摆，如下图所示，在最低端施加控制力，来保持其在竖直方向的角度可控。

其状态方程如下：

\frac{d}{dt}\begin{bmatrix} \theta\\ \dot{\theta}\end{bmatrix}= \begin{bmatrix} \dot {\theta} \\ \frac{mgl}{J_t}sin\theta-\frac{\gamma}{J_t}\dot {\theta}+\frac{l}{J_t}cos\theta \cdot u \end{bmatrix}

其中 \gamma 为旋转摩擦系数， J_t=J+ml^2 , u 为施加的外力。现在我们要开始做实验了，假设对火箭不做任何控制，即 u=0 ，这时火箭的状态方程可进一步简化：

\frac{d}{dt}\begin{bmatrix} \theta\\ \dot{\theta}\end{bmatrix}= \begin{bmatrix} \dot {\theta} \\ \frac{mgl}{J_t}sin\theta-\frac{\gamma}{J_t}\dot {\theta} \end{bmatrix}

不失一般性，假设我们自己也造了一个火箭，其 mgl=J_t=1,\gamma=1 ，则上述方程变为：

\begin{bmatrix} \dot {x}_1(t)\\ \dot {x}_2(t)\end{bmatrix}= \begin{bmatrix} x_2 \\ sin(x_1)- x_2 \end{bmatrix}

注意，这个方程虽然看着简单，但确是一个非线性方程，解起来还是费点力气的，我们要借助MATLAB帮我们算一下，并画个图出来看一下，把 x_1 作为横坐标， x_2 作为纵坐标，然后随机选择一些初始点，看看向量 \begin{bmatrix} {x_1}& {x_2}\end{bmatrix}' 是怎么运动的，轨迹如下图所示：

可见，这个图形是比较有意思的，其中有三个比较明显的点，貌似是旋涡中心，无论初始条件是什么，最终都要平衡到这三个点上（实际很多，我们只计算了三个），这三个点的坐标目测是 \begin{bmatrix} {\pm n\pi}& {0}\end{bmatrix}' 。

我们从数学上再分析一下，什么是平衡点呢？——就是不再变化的点，那什么是不再变化呢？——导数为零呗，怎么样才能让导数为零呢？——状态方程的左边就是导数啊，让右边为零就可以了。

\begin{bmatrix} x_2 \\ sin(x_1)- x_2 \end{bmatrix}=\begin{bmatrix} 0\\ 0\end{bmatrix}

解这个方程还是比较容易的，它的解就是：

\begin{bmatrix} {x_1}\\ {x_2}\end{bmatrix}=\begin{bmatrix} {\pm n\pi}\\ {0}\end{bmatrix}

但是，如果你仔细看，还有一个点， \begin{bmatrix} {0}& {0}\end{bmatrix}' 肯定是数学解，但是似乎在图上并没有明显的显示出来，这是什么原因呢？——这代表着火箭竖直放置，且没有扰动，常识告诉我们这是一个极不稳定的点，就像你把铅笔立在桌子上，稍微风吹草动就倒了，而数值求解的时候，几乎寻找不到这个点。

那其他的平衡点又代表什么呢？——火箭水平躺着， \theta=\pm\pi ,而且不再变化， \dot {\theta}=0 ，这和我们的常识也是一致的。

可见平衡点就是系统状态不再发生变化的点，它可能不止一个，它也可能很脆弱，稍微有个扰动，就不稳定了。

二、什么是李雅普诺夫稳定

早在1892年，俄国有一个叫李雅普诺夫的学者发表了一篇著名的文章《运动稳定性一般》问题，建立了关于运动稳定的一般理论，光看这个文章的名字就不一般，也确实，在尔后百余年，这个理论在数学、力学和控制理论中全面开花，已经成为稳定性研究方向的基础性理论，俄罗斯人对于数学上和工程上的直觉确实令人赞叹。

李雅普诺夫稳定性理论研究的是在扰动下平衡点的稳定性问题。

简单来说，如果平衡状态 x_e 受到扰动后，仍然 停留在 x_e 附近，我们就称 x_e 在 李雅普诺夫意义下是稳定的 （Lyapunov stable）。

如果更进一步，如果平衡状态 x_e 受到扰动后， 最终都会收敛到 x_e ，我们就称 x_e 在李雅普诺夫意义下是 渐进稳定的 （Asymptotically stable）。

再进一步，如果平衡状态 x_e 受到 任何扰动 后， 最终都会收敛到 x_e ，我们就称 x_e 在李雅普诺夫意义下是 大范围内渐进稳定的 （Asymptotically stable in large）。

相反，如果平衡状态 x_e 受到 某种扰动 后， 状态开始偏离 x_e ，我们就称 x_e 在李雅普诺夫意义下是不 稳定的 （Unstable）。

示意图如下：

下面我们就分别具体看一下。

什么是李雅普诺夫意义下的稳定

请看状态方程：
\begin{bmatrix} \dot{x_1}\\ \dot{x_2}\end{bmatrix}= \begin{bmatrix} 0 &1\\ -1&0 \end{bmatrix}\begin{bmatrix} {x_1}\\ {x_2}\end{bmatrix}

很容易得到其时域的解为;

\bm{x}(t)=\bm{x}(0)\begin{bmatrix} {sin(t)}\\ {cos(t)}\end{bmatrix}

且平衡点为 \begin{bmatrix} {0}\\ {0}\end{bmatrix} ，我们随机取一些很小的扰动，比如 \bm{x}(0)=\begin{bmatrix} {0.01}\\ {0.01}\end{bmatrix} ， \bm{x}(0)=\begin{bmatrix} {0.02}\\ {0.02}\end{bmatrix} 以及 \bm{x}(0)=\begin{bmatrix} {0.03}\\ {0.03}\end{bmatrix} ，把 \begin{bmatrix} {x_1}\\ {x_2}\end{bmatrix} 在初始条件下的轨迹画出来，结果如下：

可以看出，如果初始扰动小，其状态轨迹的区间也小；相反，当初始扰动大的时候，状态轨迹的区间也变大；但是无论如何，状态轨迹的区间都是有限的，而且，如果想减小轨迹的区间，只要保证扰动在某范围内即可。比如对于本范例，要想 {\displaystyle \|x(t)-x_e \|< 0.01} ，只要保证初始扰动 {\displaystyle \|x_{0}-x_e \|< 0.01} 即可。

翻译成严谨的数学语言就是：对于任意的 \epsilon>0 ，存在 \delta>0 ，使得如果 {\displaystyle \|x(0)-x_e \|< \delta} ，则对于所有的 t>0 ，都有 {\displaystyle \|x(t)-x_e \|< \epsilon} 。

我们再举个例子看一下：

\begin{equation} \begin{aligned} \dot x_1(t)&=x_2+x_1(2-x_1^2-x_2^2)\\ \dot x_2(t)&=-x_1+x_2(2-x_1^2-x_2^2) \end{aligned} \end{equation}

求解稍微复杂一点，我们直接画出轨迹图如下：

显然 \begin{bmatrix} {0}\\ {0}\end{bmatrix} 是其平衡点。从轨迹图中可以看出，无论是小扰动（初始点在平衡零点附近），还是大扰动，状态轨迹最终都趋向一个固定的圆。我们套用一下李雅普诺夫稳定定义看一下，该系统是否稳定。我们不妨取 \epsilon=1 ，即在零点周围画个半径为1的圆，看看能否存在 \delta ，使得 {\displaystyle \|x(0)-x_e \|< \delta} ，则对于所有的 t>0 ，都有 {\displaystyle \|x(t)-x_e \|<1} 。显然， \delta 是不存在的，因此，该系统是不稳定的。

什么是渐进稳定

请看如下方程：
\begin{equation} \begin{aligned} \dot x_1(t)&=-x_1\\ \dot x_2(t)&=x_1+x_2-x_2^3 \end{aligned} \end{equation}

令 \begin{equation} \begin{aligned} \dot x_1(t)&=0\\ \dot x_2(t)&=0 \end{aligned} \end{equation} ，很容易计算该系统的平衡点为： \begin{bmatrix} {0}\\ {0}\end{bmatrix} ， \begin{bmatrix} {0}\\ {1}\end{bmatrix} ， \begin{bmatrix} {0}\\ {-1}\end{bmatrix} ，其轨迹图为：

可见，对于平衡点 \begin{bmatrix} {0}\\ {1}\end{bmatrix} ， \begin{bmatrix} {0}\\ {-1}\end{bmatrix} ，存在 \delta>0 ，比如图中的 \delta=0.1 ，在 {\displaystyle \|x(t)-x_e \|<0.1} 的区间内，状态轨迹最终都收敛到这两个平衡点上，因此是渐进稳定的。

而对于点 \begin{bmatrix} {0}\\ {0}\end{bmatrix} ，我们在其周围添加小扰动，发现无论扰动多么小，轨迹线都会偏离该平衡点，因此属于不稳定点。

什么是大范围渐进稳定

通过前面的例子，我们发现，同一个系统里面，可能有渐进稳定点，也可能有不稳定点，我们还不容易确定系统是否稳定。再来看一个简单的例子：

\begin{bmatrix} \dot{x_1}\\ \dot{x_2}\end{bmatrix}= \begin{bmatrix} 1 &-3\\ 5&-2\end{bmatrix}\begin{bmatrix} {x_1}\\ {x_2}\end{bmatrix}

其轨迹图为：

可见，对于任何扰动，最后都会收敛到一个平衡点，这就是大范围内渐进稳定。对于线性系统来说，很容易证明，如果平衡态是渐进稳定的，也必然是大范围渐进稳定的。

前面我们的文章说了，对于线性时不变系统，只要矩阵 A 的特征值具有负实部，那系统就是大范围渐进稳定的。

什么是不稳定

假设有一个系统，状态方程如下:
\begin{bmatrix} \dot{x_1}\\ \dot{x_2}\end{bmatrix}= \begin{bmatrix} 1 &3\\ -5&2 \end{bmatrix}\begin{bmatrix} {x_1}\\ {x_2}\end{bmatrix}

同样随机布置一些初始点，看看状态的轨迹如何：

很明显，所有初始状态的轨迹都呈现螺旋发散状，从数学上看，特征值 \lambda=1.5\pm 3.84i ，具有共轭根，但是实部是正的，因此发散。

再来看另外一个例子，状态方程为
\begin{bmatrix} \dot{x_1}\\ \dot{x_2}\end{bmatrix}= \begin{bmatrix} 4 &-2\\ 1&-3 \end{bmatrix}\begin{bmatrix} {x_1}\\ {x_2}\end{bmatrix}

状态量的轨迹为：

貌似也呈现发散状，但是和前面的例子有所不同，轨迹貌似呈现指数形式，计算状态方程特征值为 \lambda=\begin{bmatrix} 3.7\\-2.7\end{bmatrix} ，有一个不稳定的实根，导致系统发散。由以上两例可以看出，发散的轨迹可以有多种多样。

三、李雅普诺夫第一法

前面我们把平衡点分了几类，我们会发现，还是线性系统比较好计算，而且性能比较好，只要保证矩阵 A 具有负实部，就是大范围一致稳定的。因此，我们如果把方程的形式由

\dot {\bm{x}}(t)=f(\bm{x}(t))

改成

\dot {\bm{x}}(t)=A\bm{x}(t)

就会带来很多方便。这就需要将非线性系统在平衡态附近线性化，然后讨论线性化系统的特征值分布来研究原非线性系统的稳定性问题。这种方法，就是李雅普诺夫在他论文中提到的第一种方法，称之为第一法，也叫间接法。

我们再来分析一下前面所说的倒立摆的例子：

\begin{equation} \begin{aligned} \dot x_1(t)&=x_2\\ \dot x_2(t)&=sin(x_1)-x_2 \end{aligned} \end{equation}

这是一个典型的非线性方程，我们前面计算过了，其有多个平衡点 \begin{bmatrix} {\pm n\pi}& {0}\end{bmatrix}' ，我们不妨来研究一下 \begin{bmatrix} {\pi}& {0}\end{bmatrix}' 这个点，其附件的状态轨迹为：

现在将原来的非线性方程线性化，为方便起见，我们定义 z_1=x_1-\pi 以及 z_2=x_2 ，于是可以得到：

sin(\pi+z_1)=-sinz_1\approx -z_1

于是，原非线性方程就变为：

\begin{bmatrix} \dot{z_1}\\ \dot{z_2}\end{bmatrix}= \begin{bmatrix} 0 &1\\ -1&-1 \end{bmatrix}\begin{bmatrix} {z_1}\\ {z_2}\end{bmatrix}

在新的坐标系下，平衡点变为 \begin{bmatrix} {0}\\ {0}\end{bmatrix} ，其轨迹为：

可见，与原轨迹还是比较接近的。一般的书上，对于李雅普诺夫第一法都是一笔带过，其实在工程实践中，第一法应用非常多，比如复杂的飞机飞行控制，就是将飞机模型线性化成多个线性化模型进行设计，感兴趣的可参见 Design an LQR Servo Controller in Simulink 。

四、李雅普诺夫第二法

第二法就比较天才了，来源于一个朴素的想法：稳定的系统能量总是不断被耗散的，李雅普诺夫通过定义一个标量函数 V(\bm{x}) （通常能代表广义能量）来分析稳定性。这种方法的避免了直接求解方程，也没有进行近似线性化，所以也一般称之为直接法。如果标量函数 V(x) 满足：

V(\bm{x})=0 \ \text{if and only if } \bm{x}=0
V(\bm{x})>0 \ \text{if and only if } \bm{x}\ne0
\dot{V}(\bm{x})=\frac{d}{dt}V(\bm{x})=\sum_{i=1}^{n}\frac{\partial V}{\partial x_i}f_i(x)\leq 0\ \text{when}\ \bm{x} \ne 0

则称系统在李雅普诺夫意义下是稳定的，特别的，若 \bm{x} \ne 0 时，有 \dot{V}(\bm{x}) <0 ，则系统是渐进稳定的。举个例子：

\begin{equation} \begin{aligned} \dot {x}_1(t)&=x_2-x_1(x_1^2+x_2^2)\\ \dot {x}_2(t)&=-x_1-x_2(x_1^2+x_2^2) \end{aligned} \end{equation}

如果我们定义李雅普诺夫函数

V(\bm{x})=x_1^2+x_2^2

则有

\begin{equation} \begin{aligned} \dot {V}(\bm{x})&=2x_1 \dot x_1+2x_2 \dot x_2 \\ &=-2(x_1^2+x_2^2)^2 \end{aligned} \end{equation}
显然当 \bm{x} \ne 0 时，有 \dot {V}(\bm{x})<0 ，所以系统是渐进稳定的。

可见，如果能合理的选定李雅普诺夫函数，则非常容易的判断系统的稳定性。不过遗憾的是，对于复杂的系统，李雅普诺夫函数的选择可以称得上一门玄学，所以，对于工程师而言，笔者还是喜欢李雅普诺夫第一法。

五、MATLAB代码

鉴于很多知友对文章中插图的MATLAB代码感兴趣，先将部分代码附录如下，其余按格式更改即可。

首先是定义状态方程函数：

function d=dxdt(t,x)
d=[ x(2)+x(1)*(2-x(1)^2-x(2)^2); 
    -x(1)+x(2)*(2-x(1)^2-x(2)^2) ];

根据状态方程，画出变量轨迹：

figure('color','w');
hold on 
for theta=[0:20]*pi/10
    x0=3*[cos(theta);sin(theta)];%定义初始值数组
    [t,x]=ode45(@dxdt,[0:0.1:8],x0);
    plot(x(:,1),x(:,2),'linewidth',0.5)
    quiver(x(:,1),x(:,2),gradient(x(:,1)),gradient(x(:,2)),'linewidth',3.0);%增加轨迹方向箭头
for theta=[0:2:20]*pi/10
    x0=1e-5*[cos(theta);sin(theta)];
    [t,x]=ode45(@dxdt,[0:0.2:20],x0);