Chapter9.3：线性系统稳定性分析及综合实例

该系列博客主要讲述Matlab软件在自动控制方面的应用，如无自动控制理论基础，请先学习自动控制系列博文，该系列博客不再详细讲解自动控制理论知识。
自动控制理论基础相关链接：https://blog.csdn.net/qq_39032096/category_10287468.html?spm=1001.2014.3001.5482
博客参考书籍：《MATLAB/Simulink与控制系统仿真》。

---

---

3.线性系统稳定性分析及综合实例

李雅普诺夫理论提出了判断系统稳定性的两种方法：一种是利用线性系统微分方程的解来判断系统稳定性，称之为李雅普诺夫第一法或间接法；另一种先利用经验和技巧来构造李雅普诺夫函数，进而利用李雅普诺夫函数来判断系统稳定性，称之为李雅普诺夫第二法或直接法；

3.1 李雅普诺夫意义下的稳定性

设系统方程为：
$$\dot{x}=f(x,t)$$
式中：$x$为$n$维状态向量，且显含时间变量$t$，$f(x,t)$为线性或非线性、定常或时变的$n$维向量函数，其展开式为：
$${\dot{x}}_i=f_i(x_1,x_2,\cdots ,x_n,t);i=1,2,\cdots ,n$$
假定方程的解为$x(t;x_0,t_0)$，式中$x_0、t_0$分别为初始状态向量和初始时刻，初始条件$x_0$必满足$x(t_0;x_0,t_0)=x_0$;

1. 平衡状态

   对于所有$t$，满足：
   $${\dot{x}}_e=f(x_e,t)=0$$
   的状态$x_e$称为平衡状态；

   线性定常系统$\dot{x}=Ax$，其平衡状态满足$A{x}_e=0$，当$A$为非奇异矩阵时，系统只有唯一的零解，即只存在一个位于状态空间原点的平衡状态，若$A$为奇异矩阵，则系统存在有无穷多个平衡状态，对于非线性系统，可能有一个或多个平衡状态；

2. 李雅普诺夫意义下的稳定性

   设系统初始状态位于以平衡状态$x_e$为球心、$\delta$为半径的闭球域$S(\delta )$内，即：
   $$||x_0-x_e||\le \delta ,t=t_0$$
   若能使系统方程的解$x(t;x_0,t_0)$在$t\to \infty$的过程中，都位于以$x_e$为球心、任意规定的半径为$\epsilon$的闭球域$S(\epsilon )$内，即：
   $$||x(t;x_0,t_0)-x_0||\le \epsilon ,t\ge t_0$$
   则称系统的平衡状态$x_e$在李雅普诺夫意义下是稳定的；

   式中：$||\cdot ||$为欧几里得范数，其几何意义是空间距离的尺度，如$||x_0-x_e||$表示状态空间$x_0$点至$x_e$点之间距离的尺度，数学表达式为：
   $$||x_0-x_e||=\sqrt{(x_{10}-x_{1e}{)}^2+\cdots +(x_{n0}-x_{ne}{)}^2}$$

3. 渐近稳定性

   若系统的平衡状态$x_e$不仅具有李雅普诺夫意义下的稳定性，且有：
   $$\underset{t\to \infty }{\lim }||x(t;x_0,t_0)-x_e||=0$$
   则称此平衡状态是渐近稳定的；

4. 大范围(全局)渐近稳定性

   当初始条件扩展至整个状态空间，且平衡状态均具有渐近稳定性时，称此平衡状态是大范围渐近稳定的，此时$\delta \to \infty ,S(\delta )\to \infty$；当$t\to \infty$时，由状态空间中任一点出发的轨迹都收敛至$x$；

   对于严格线性的系统，如果它是渐近稳定的，则必定是大范围渐近稳定的，对于非线性系统来说，其稳定性与初始条件密切相关，系统渐近稳定不一定是大范围渐近稳定；

5. 不稳定性

   如果对于某个实数$\epsilon >0$和任一个实数$\delta >0$，不管这两个实数有多小，在$S(\delta )$内总存在一个状态$x_0$，使得由这一状态出发的轨迹超出$S(\epsilon )$，则平衡状态$x_e$称为是不稳定的；

3.2 李雅普诺夫第一法(间接法)

李雅普诺夫第一法是利用状态方程解的特性来判断系统稳定性的方法，适用于线性定常、线性时变及非线性函数可线性化的情况；

定理9：对于线性定常系统$\dot{x}=Ax,x(0)=x_0,t\ge 0$，有：

• 系统的每一平衡状态是在李雅普诺夫意义下稳定的充分必要条件：$A$的所有特征值均具有非正(负或零)实部，且具有零实部的特征值为$A$的最小多项式的单根；
• 系统的唯一平衡状态$x_e=0$是渐近稳定的充分必要条件：$A$的所有特征值均具有负实部；

3.3 李雅普诺夫第二法(直接法)

李雅普诺夫提出，可虚构一个能量函数，被称为李雅普诺夫函数，一般它与$x_1,x_2,\cdots ,x_n$及$t$有关，记为$V(x,t)$;

若不显含$t$，则记为$V(x)$；它是一个标量函数，考虑到能量总是大于零，故为正定函数；能量衰减特性用$\dot{V}(x,t)$或$\dot{V}(x)$表示；李雅普诺夫第二法用$V$及$\dot{V}$的符号特征，直接对平衡状态稳定性做出判断，无需求出系统状态方程的解；

对于线性系统，通常用二次型函数$x^{T}Px$作为李雅普诺夫函数；

1. 标量函数定号性简单介绍

   正定性：标量函数$V(x)$对所有在域$S$中的非零状态$x$有$V(x)>0$且$V(0)=0$，则在域$S$内的标量函数$V(x)$称为是正定的；

   负定性：如果$-V(x)$是正定函数，则标量函数$V(x)$称为负定函数；

   正半定性：如果标量函数$V(x)$除了原点及某些状态处等于零外，在域$S$内的所有状态都是正定的，则$V(x)$称为正半定函数；

   负半定性：如果$-V(x)$是正半定函数，则标量函数$V(x)$称为负半定函数；

   不定性：如果在域$S$内，不论域$S$多么小，$V(x)$即可为正值，也可为负值，则标量函数$V(x)$称为不定函数；

2. 李雅普诺夫第二法主要定理

   定理10(定常系统大范围渐近稳定判别定理1)：对于定常系统：
   $$\dot{x}=f(x),t\ge 0$$
   其中：$f(0)=0$，如果存在一个具有连续一阶导数的标量函数$V(x)，V(0)=0$，且对于状态空间$X$中的一切非零点$x$满足如下条件：

   • $V(x)$为正定；
   • $\dot{V}(x)$为负定；
   • 当$||x||\to \infty$时$V(x)\to \infty$.

   则系统的原点平衡点状态是大范围渐近稳定的.

   实例分析：

   Example1： 设系统状态方程为：
   $$\begin{array}{rl}& {\dot{x}}_1=x_2-x_1(x_1^2+x_2^2)\\ & {\dot{x}}_2=-x_1-x_2(x_1^2+x_2^2)\end{array}$$
   确定系统的稳定性.

   解：

   原点$(x_1=0,x_2=0)$是该系统唯一的平衡状态；

   选定正定标量函数$V(x)$为：
   $$V(x)=x_1^2+x_2^2$$
   则沿任意轨迹$V(x)$对时间的导数：
   $$\dot{V}(x)=2x_1{\dot{x}}_1+2x_2{\dot{x}}_2=-2(x_1^2+x_2^2{)}^2$$
   是负定的.

   说明$V(x)$沿任意轨迹是连续减小的，因此$V(x)$是一个李雅普诺夫函数.

   由于当$||x||\to \infty$时$V(x)\to \infty$，所以系统在原点处的平衡状态是大范围渐近稳定的.

   定理11(定常系统大范围渐近稳定判别定理2)：对于定常系统，如果存在一个具有连续一阶导数的标量$V(x)，V(x)=0$，且对状态空间$X$中的一切非零点$x$满足如下条件：

   • $V(x)$为正定；
   • $\dot{V}(x)$为负半定；
   • 对任意$x\in X，\dot{V}(x(t;x_0,0))\not\equiv 0$；
   • 当$||x||\to \infty$时，$V(x)\to \infty$.

   则系统的原点平衡状态是大范围渐近稳定的.

   实例分析：

   Example2： 已知定常系统状态方程为：
   $$\begin{array}{rl}& {\dot{x}}_1=x_2\\ & {\dot{x}}_2=-x_1-(1+x_2{)}^2{x}_2\end{array}$$
   确定系统的稳定性.

   解：

   原点$(x_1=0,x_2=0)$为系统唯一的平衡状态.

   取$V(x)=x_1^2+x_2^2$，且有：

   $V(x)=x_1^2+x_2^2$为正定，

   $\dot{V}(x)=2x_1{\dot{x}}_1+2x_2{\dot{x}}_2=-2x_2^2(1+x_2{)}^2$.容易看出，除了①$x_1$任意，$x_2=0$；②$x_1$任意，$x_2=-1$时，$\dot{V}(x)=0$外，均有$\dot{V}(x)<0$，所以，$\dot{V}(x)$为负半定.

   检查是否$\dot{V}(x(t;x_0,0))\not\equiv 0$.考虑到使得$\dot{V}(x)=0$的可能性只有上述①和②两种情况，所有问题归结为判断这两种情况是否为系统的受扰运动解.

   考察情况①：设$\overline{x}(t;x_0,0)={\left[\begin{array}{cc}{x}_1(t)& 0\end{array}\right]}^T$，则由$x_2(t)\equiv 0$可导出${\dot{x}}_2=0$，将此代入系统状态方程可得：
   $$\begin{array}{rl}& {\dot{x}}_1(t)=x_2(t)=0\\ & 0={\dot{x}}_2(t)=-(1+x_2(t){)}^2{x}_2(t)-x_1(t)=-x_1(t)\end{array}$$
   表明，除了点$(x_1=0,x_2=0)$外，$\overline{x}(t;x_0,0)=[x_1(t)0{]}^T$不是系统的受扰运动解.

   考察情况②：设$\overline{x}(t;x_0,0)=[x_1(t)-1{]}^T$，则由$x_2(t)=-1$可导出${\dot{x}}_2(t)=0$，将此代入系统状态方程可得：
   $$\begin{array}{rl}& {\dot{x}}_1(t)=x_2(t)=-1\\ & 0={\dot{x}}_2(t)=-(1+x_2(t){)}^2{x}_2(t)-x_1(t)=-x_1(t)\end{array}$$
   这是一个矛盾的结果，表明$\overline{x}(t;x_0,0)=[x_1(t)-1{]}^T$也不是系统受扰运动解.

   综合以上分析，$\dot{V}(x(t;x_0,0))\not\equiv 0$.

   当$||x||\to \infty$时，显然有$V(x)=||x|{|}^2\to \infty$.

   因此，系统的原点平衡状态是大范围渐近稳定的.

   定理12(不稳定判别定理)：对于定常系统，如果存在一个具有连续一阶导数的标量函数$V(x)$($V(0)=0$)，和围绕原点的域$\Omega$，使得对于一切$x\in \Omega$和一切$t\ge t_0$满足如下条件：

   • $V(x)$为正定；
   • $\dot{V(x)}$为正定；

   则系统平衡状态为不稳定.

3.4 线性定常系统的李雅普诺夫稳定性分析

1. 线性定常连续系统渐近稳定性的判别

   设线性定常系统状态方程为：
   $$\dot{x}=Ax,x(0)=x_0,t\ge 0$$
   $A$为非奇异矩阵.故原点是唯一平衡状态.

   设取正定二次型函数$V(x)=x^{T}Px$作为可能的李雅普诺夫函数，考虑到系统状态方程，则有：
   $$\dot{V}(x)={\dot{x}}^{T}Px+x^{T}P\dot{x}=x^T(A^{T}P+PA)x$$
   令：
   $$A^{T}P+PA=-Q\text{\hspace{1em}(1)}$$
   有：
   $$\dot{V}(x)=-x^{T}Qx$$
   线性定常连续系统渐近稳定的充分必要条件：给定一正定矩阵$P$，存在着满足式$(1)$的正定矩阵$Q$，$x^{T}Px$是该系统的一个李雅普诺夫函数，式$(1)$称为李雅普诺夫矩阵代数方程。使用中常选取$Q$阵为单位阵或对角线阵。

   定理13：线性定常系统$\dot{x}=Ax,x(0)=x_0,t\ge 0$的原点平衡状态$x_e=0$为渐近稳定的充分必要条件：对于任意给定的一个正定对称矩阵$Q$，有唯一的正定对称矩阵$P$使式$(1)$成立.

   实例分析：

   Example3： 已知线性定常连续系统状态方程为：
   $${\dot{x}}_1=x_2,{\dot{x}}_2=2x_1-x_2$$
   用李雅普诺夫方程判定系统的渐近稳定性.

   解：

   先用特征值判据判断，系统状态方程为：
   $$\dot{x}=\left[\begin{array}{cc}0& 1\\ 2& -1\end{array}\right]x,A=\left[\begin{array}{cc}0& 1\\ 2& -1\end{array}\right]$$

   $$|\lambda I-A|=\left|\begin{array}{cc}\lambda & -1\\ -2& \lambda +1\end{array}\right|={\lambda }^2+\lambda -2=(\lambda -1)(\lambda +2)=0$$

   特征值为：$-2,1$，故系统不稳定.

   令：
   $$A^{T}P+PA=-Q=-I$$

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

   则有：
   $$\left[\begin{array}{cc}0& 2\\ 1& -1\end{array}\right]\left[\begin{array}{cc}{P}_{11}& P_{12}\\ P_{21}& P_{22}\end{array}\right]+\left[\begin{array}{cc}{P}_{11}& P_{12}\\ P_{21}& P_{22}\end{array}\right]\left[\begin{array}{cc}0& 1\\ 2& -1\end{array}\right]=\left[\begin{array}{cc}-1& 0\\ 0& -1\end{array}\right]$$
   展开：
   $$4P_{12}=-1,2P_{12}-2P_{22}=-1,P_{11}-P_{12}+2P_{22}=0$$
   解得：
   $$P=\left[\begin{array}{cc}{P}_{11}& P_{12}\\ P_{21}& P_{22}\end{array}\right]=\left[\begin{array}{cc}-\frac{3}{4}& -\frac{1}{4}\\ -\frac{1}{4}& \frac{1}{4}\end{array}\right]$$
   由于$P_{11}=-\frac{3}{4}，\det P=-\frac{1}{4}<0$，故$P$不定，可知系统非渐近稳定.由特征值判据知系统是不稳定的.

2. 线性定常离散系统渐近稳定性的判别

   设线性定常离散系统状态方程为：
   $$x(k+1)=\Phi x(k),x(0)=x_0;k=0,1,2,\cdots ,$$
   式中：$\Phi$阵为非奇异；原点是平衡状态.

   取正定二次型函数：
   $$V(x(k))=x^T(k)Px(k)$$
   以$\Delta V(x(k))$代替$\dot{V}(x)$，有：
   $$\Delta V(x(k))=V(x(k+1))-V(x(k))$$
   考虑到线性定常离散系统状态方程式有：
   $$\Delta V(x(k))=x^T(k)({\Phi }^{T}P\Phi -P)x(k)$$
   令：
   $${\Phi }^{T}P\Phi -P=-Q\text{\hspace{1em}(2)}$$
   有：
   $$\Delta V(x(k))=-x^T(k)Qx(k)$$
   定理14：线性定常离散系统渐近稳定的充分必要条件：给定任一正定对称矩阵$Q$，存在一个正定对称矩阵$P$使式$(2)$成立；

   $x^T(k)Px(k)$是系统的一个李雅普诺夫函数，式$(2)$称为离散的李雅普诺夫代数方程，通常可取$Q=I$.

   如果$\Delta V(x(k))$沿任一解的序列不恒为零，则$Q$可取为正半定矩阵.

3. 李雅普诺夫稳定性判据汇总表：

   

3.5 MATLAB/SIMULINK在李雅普诺夫稳定性分析中的应用

• $\mathrm{P}=\mathrm{lyap}(\mathrm{A},\mathrm{Q})$：$\mathrm{A}$是已知系统的状态矩阵，$\mathrm{Q}$是给定的正定对称矩阵，$\mathrm{P}$是李雅普诺夫方程$\mathrm{AP}+\mathrm{P}{\mathrm{A}}^{\mathrm{T}}+\mathrm{Q}=0$的解，即正定实对称矩阵$\mathrm{P}$；
• $\mathrm{P}=lyap2(\mathrm{A},\mathrm{Q})$：$\mathrm{A}$是已知系统的状态矩阵，$\mathrm{Q}$是给定的正定对称矩阵，$\mathrm{P}$是李雅普诺夫方程$\mathrm{AP}+\mathrm{P}{\mathrm{A}}^{\mathrm{T}}+\mathrm{Q}=0$的解，即正定实对称矩阵$\mathrm{P}$；$lyap2()$采用特征值分解法求解李雅普诺夫方程，运算速度比$\mathrm{lyap}()$快；
• $\mathrm{P}=\mathrm{dlyap}(\mathrm{A},\mathrm{Q})$：$\mathrm{A}$是已知系统的状态矩阵，$\mathrm{Q}$是给定的正定对称矩阵，$\mathrm{P}$是李雅普诺夫方程$\mathrm{AP}+\mathrm{P}{\mathrm{A}}^{\mathrm{T}}+\mathrm{Q}=0$的解，即正定实对称矩阵$\mathrm{P}$，用于离散系统求解；

【实战】：已知系统状态矩阵$A=\left[\begin{array}{cc}0& 1\\ -1& -1\end{array}\right]$，给定正定对称矩阵$Q=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]$，确定李雅普诺夫方程的解$\mathrm{P}$;

解：

% 实例Chapter9.3.5
clc;clear;A=[0,1;-1,-1];              % 状态矩阵
Q=[1,0;0,1];                % 给定的正定对称矩阵
P=lyap(A',Q)% 	求解得到的矩阵P
%	P =
%    	1.5000    0.5000
%    	0.5000    1.0000

3.6 综合实例及MATLAB/SIMULINK应用

实验要求：已知线性定常系统状态空间表达式如下：
$$\{\begin{array}{ll}& \left[\begin{array}{c}{x}_1^{\prime }\\ x_2^{\prime }\\ x_3^{\prime }\end{array}\right]=\left[\begin{array}{ccc}0& 1& 0\\ 0& -2& 1\\ -1& 0& -1\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right]+\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right]u\\ & \\ & y=\left[\begin{array}{ccc}1& 0& 0\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right]\end{array}$$
使用$\mathrm{MATLAB}$进行以下分析：

1. 分析系统的可控性、可观测性；
2. 对系统进行非奇异线性变换，使其状态转移矩阵对角化，进而分析变换后系统的可控性、可观性；
3. 分析系统的稳定性，绘制系统的阶跃响应曲线；
4. 使用$\mathrm{LTI}\mathrm{Viewer}$工具绘制上述系统的阶跃响应和冲激响应曲线；

解：

【$STEP1$】：建立控制系统的数学模型。

% 实例Chapter9.3.6
clc;clear;% 建立系统的数学模型
A=[0,1,0;0,-2,1;-1,0,-1];
B=[0;0;1];C=[1,0,0];D=0;
sys=ss(A,B,C,D)

% 控制系统数学模型结果
sys =A = x1  x2  x3x1   0   1   0x2   0  -2   1x3  -1   0  -1B = u1x1   0x2   0x3   1C = x1  x2  x3y1   1   0   0D = u1y1   0

【$STEP2$】：系统的可控性、可观性。

% 控制系统的可控性、可观性
control_matrix=ctrb(A,B);
rank_control=rank(control_matrix);if rank_control<3disp('系统不可控！');
elsedisp('系统完全可控！');
endobserve_matrix=obsv(A,C);
rank_observe=rank(observe_matrix);if rank_observe<3disp('系统不可观！');
elsedisp('系统完全可观！');
end

% 系统的可控性、可观性结果
系统完全可控！
系统完全可观！

【$STEP3$】：对系统进行线性对角化变换。

% 系统进行线性对角化变换
[V,S]=eig(A);               % 求取矩阵特征值和特征向量
if rank(V)==3               % 矩阵满秩，进行对角化变换AA=inv(V)*A*V;BB=inv(V)*B;CC=C*V;disp('系统对角化变换后的新系统模型：');sys_diag=ss(AA,BB,CC,D)
elsedisp('这个系统不可对角化！');
end

% 线性对角化变换后系统模型
系统对角化变换后的新系统模型：sys_diag =A = x1                   x2                   x3x1        -0.338+0.562i      -5e-16-2.5e-16i   2.78e-17+6.18e-16ix2  -3.89e-16+2.78e-16i        -0.338-0.562i   -8.33e-17-8.4e-16ix3   3.33e-16+1.95e-16i   3.33e-16-2.21e-16i      -2.32-4.99e-17iB = u1x1      0.593+0.398ix2      0.593-0.398ix3  -0.619-1.57e-18iC = x1            x2            x3y1  -0.459-0.39i  -0.459+0.39i        -0.379D = u1y1   0

【$STEP4$】：线性变换后系统可控性、可观性。

% 线性变换后控制系统的可控性、可观性
control_matrix=ctrb(AA,BB);
rank_control=rank(control_matrix);if rank_control<3disp('变换后系统不可控！');
elsedisp('变换后系统完全可控！');
endobserve_matrix=obsv(AA,CC);
rank_observe=rank(observe_matrix);if rank_observe<3disp('变换后系统不可观！');
elsedisp('变换后系统完全可观！');
end% 线性变换后可控性、可观性结果
% 变换后系统完全可控！
% 变换后系统完全可观！

【$STEP5$】：分析系统的稳定性。

% 利用李雅普诺夫方法分析系统的稳定性
P=lyap(A',eye(3))% 系统稳定性结果：P矩阵为正定矩阵，系统渐近稳定。
%P =
%    3.5000    1.5000    0.5000
%    1.5000    1.0000    0.5000
%    0.5000    0.5000    1.0000

【$STEP6$】：绘制系统的阶跃响应曲线。

% 绘制系统的阶跃响应曲线
step(sys)
set(findobj(get(gca,'Children'),'LineWidth',0.5),'LineWidth',1.5);
title('控制系统阶跃响应','FontSize',15);

【$STEP7$】：使用$\mathrm{LTI}\mathrm{Viewer}$工具进行分析。

ltiview;

【完整代码】

% 实例Chapter9.3.6
clc;clear;% STEP1：建立系统的数学模型
A=[0,1,0;0,-2,1;-1,0,-1];
B=[0;0;1];C=[1,0,0];D=0;
sys=ss(A,B,C,D)% STEP2：控制系统的可控性、可观性
control_matrix=ctrb(A,B);
rank_control=rank(control_matrix);if rank_control<3disp('系统不可控！');
elsedisp('系统完全可控！');
endobserve_matrix=obsv(A,C);
rank_observe=rank(observe_matrix);if rank_observe<3disp('系统不可观！');
elsedisp('系统完全可观！');
end% STEP3：系统进行线性对角化变换
[V,S]=eig(A);               % 求取矩阵特征值和特征向量
if rank(V)==3               % 矩阵满秩，进行对角化变换AA=inv(V)*A*V;BB=inv(V)*B;CC=C*V;disp('系统对角化变换后的新系统模型：');sys_diag=ss(AA,BB,CC,D)
elsedisp('这个系统不可对角化！');
end% STEP4：线性变换后控制系统的可控性、可观性
control_matrix=ctrb(AA,BB);
rank_control=rank(control_matrix);if rank_control<3disp('变换后系统不可控！');
elsedisp('变换后系统完全可控！');
endobserve_matrix=obsv(AA,CC);
rank_observe=rank(observe_matrix);if rank_observe<3disp('变换后系统不可观！');
elsedisp('变换后系统完全可观！');
end% STEP5：利用李雅普诺夫方法分析系统的稳定性
P=lyap(A',eye(3))% STEP6：绘制系统的阶跃响应曲线
step(sys)
set(findobj(get(gca,'Children'),'LineWidth',0.5),'LineWidth',1.5);
title('控制系统阶跃响应','FontSize',15);