1. 从一个“反直觉”的仿真现象说起
如果你做过控制系统仿真，大概率遇到过这种场景：控制器在时域响应里看着挺稳，超调不大、调节时间也能接受，可一旦把相轨迹画出来，或者把初始条件稍微改一改，系统就像脱缰的野马，状态变量直接发散到无穷大。更让人头疼的是，有些非线性系统在小扰动下收敛得漂漂亮亮，扰动稍微大一点就彻底失控。这类问题，靠看阶跃响应曲线是看不出来的，得换一套工具——李雅普诺夫稳定性分析。
李雅普诺夫稳定性分析在MATLAB里的实现，核心就两件事：一是构造一个标量函数（李雅普诺夫函数），二是判断它沿系统轨迹的导数是否满足符号条件。听起来简单，但真正动手做的时候，难点全在“怎么构造”和“怎么验证”上。线性系统还好，有现成的李雅普诺夫方程可以解；非线性系统就麻烦了，构造方法五花八门，没有一个万能公式。这篇文章面向的是已经学过自动控制原理、手上有点MATLAB基础、但一碰到非线性稳定性分析就发怵的读者。我会把线性系统和非线性系统的分析路径分开讲，重点放在实操层面——怎么在MATLAB里把李雅普诺夫函数构造出来、怎么用数值方法验证、遇到不收敛的情况怎么排查。全文的代码都可以直接复制运行，参数我会解释清楚为什么这么选。
2. 李雅普诺夫稳定性到底在判断什么
2.1 稳定性的三种层次：稳定、渐近稳定、指数稳定
很多人把“稳定”当成一个非黑即白的概念，实际上李雅普诺夫框架下至少分三个层次。稳定（Lyapunov stable）指的是：初始状态稍微偏一点，后续轨迹不会跑远，但也不保证一定回到原点，可能在一个等幅振荡的闭轨上转圈。渐近稳定（asymptotically stable）要求更强：不仅不跑远，最终还得收敛到平衡点。指数稳定（exponentially stable）则要求收敛速度有下界，形如 ( |x(t)| \leq k|x(0)|e^{-\lambda t} ) 。
为什么要在意这个区分？因为工程上很多场景只要求稳定就够了，比如某些振荡器设计，你就是要它维持等幅振荡；但大多数控制场景要求渐近稳定，否则稳态误差消不掉。在MATLAB里做数值验证时，这三种层次的判断标准不一样：稳定看轨迹是否有界，渐近稳定看轨迹是否趋于零，指数稳定还要拟合收敛速率。我见过不少人拿一个渐近稳定的系统去套指数稳定的判据，结果怎么调参数都不对，问题出在概念层次没对齐。
2.2 李雅普诺夫第二法的核心逻辑
李雅普诺夫第二法（直接法）的精髓在于：不求解微分方程，而是找一个“能量类”的标量函数 ( V(x) ) ，通过它的符号性质来判断稳定性。这个函数需要满足两个条件：( V(x) ) 在平衡点附近正定，( \dot{V}(x) ) 沿系统轨迹负定（或半负定）。物理直觉是：如果系统的“广义能量”一直在衰减，那状态就不可能跑到无穷远去。
这里有个容易踩的坑：( V(x) ) 不是唯一的。同一个系统可以构造出无穷多个满足条件的李雅普诺夫函数，有的保守、有的精确。对于线性系统，二次型函数 ( V(x)=x^TPx ) 是标准选择，因为可以转化为解李雅普诺夫方程；对于非线性系统，二次型不一定够用，可能需要引入高次项、三角函数项甚至分段函数。MATLAB在这个环节的角色是：帮你验证你构造的 ( V(x) ) 到底满不满足条件，而不是替你自动构造。这一点必须清醒，否则会陷入“换个函数试试”的盲目试错。
2.3 为什么MATLAB是验证李雅普诺夫条件的合适工具
手工验证 ( \dot{V}(x) ) 的符号性质，对于简单系统还能凑合，一旦维数上去或者非线性项复杂，符号计算量会爆炸。MATLAB的优势在于三点：第一，符号计算工具箱可以直接对 ( V(x) ) 求梯度、代入系统方程、化简 ( \dot{V}(x) ) 的表达式；第二，数值计算能力可以在状态空间网格上批量计算 ( V(x) ) 和 ( \dot{V}(x) ) 的值，用可视化手段直观判断符号；第三，控制系统工具箱提供了lyap和dlyap函数，线性系统的李雅普诺夫方程求解一行代码搞定。
不过要注意，MATLAB的符号计算在复杂表达式上可能化简不彻底，这时候需要结合数值验证交叉确认。我的习惯是：符号计算给出解析形式，数值计算在网格上扫一遍，两者结论一致才放心。如果只依赖符号计算，可能因为化简规则的问题漏掉某些区域；只依赖数值计算，又可能因为网格太粗错过局部信息。
3. 线性系统的李雅普诺夫方程实操
3.1 连续时间李雅普诺夫方程的标准解法
线性定常系统 ( \dot{x}=Ax ) 的稳定性判据是：给定任意正定对称矩阵 ( Q ) ，解李雅普诺夫方程 ( A^TP+PA=-Q ) ，如果得到的 ( P ) 正定，则系统渐近稳定。MATLAB里用lyap函数直接求解，语法是P = lyap(A, Q)。注意这里的符号约定：MATLAB的lyap求解的是 ( AP+PA^T+Q=0 ) ，和教科书上的写法可能差一个转置，用之前最好用简单例子验证一下方向。
我通常取 ( Q=I ) （单位阵），因为这样计算最省事，而且不影响判据的结论——只要存在一个正定的 ( P ) 就行，( Q ) 的具体形式不影响稳定性判断。下面这段代码演示了一个二阶系统的完整流程：
% 定义系统矩阵 A = [0 1; -2 -3]; % 取Q为单位阵 Q = eye(2); % 求解李雅普诺夫方程 P = lyap(A, Q); % 检查P的正定性 eig_P = eig(P); disp('P的特征值：'); disp(eig_P); if all(eig_P > 0) disp('P正定，系统渐近稳定'); else disp('P非正定，系统不稳定或临界稳定'); end这段代码跑出来，P的特征值都是正的，说明系统渐近稳定。你可以把A改成[0 1; 2 3]试试，会看到P出现负特征值，系统不稳定。这种对比测试是理解判据的好办法。
3.2 离散时间系统的处理差异
离散系统 ( x(k+1)=A_dx(k) ) 的李雅普诺夫方程是 ( A_d^TPA_d-P=-Q ) ，MATLAB用dlyap求解。这里有个细节：连续系统的稳定性要求特征值实部为负，离散系统要求特征值模小于1。在构造测试用例时，离散系统的A_d不能直接拿连续系统的A来用，需要先做离散化，比如用零阶保持器：Ad = expm(A*Ts)，其中Ts是采样周期。
采样周期的选择会影响离散系统的稳定性。Ts太大，原本稳定的连续系统离散化后可能变得不稳定。我做过一个测试：连续系统极点在-10附近，Ts取0.01秒时离散系统稳定，取0.5秒时直接发散。所以在做离散系统分析前，先确认采样周期是否合理，否则后面的李雅普诺夫分析全是白费功夫。
3.3 用特征值交叉验证结果
lyap函数给出的P是否正定，本质上等价于A的特征值是否都在左半平面。我习惯用eig(A)做一次交叉验证，确保没有因为数值精度问题导致误判。特别是当A的特征值靠近虚轴时，lyap的解可能对Q的选择很敏感，这时候换几个不同的Q试试，如果结论一致才可信。
还有一个实用技巧：如果A的特征值实部有正有负，系统本身就是不稳定的，没必要再解李雅普诺夫方程。先用eig快速筛一遍，能省不少时间。只有当特征值都在左半平面但你想量化“稳定程度”时，才需要进一步分析P的条件数——条件数越大，说明系统对扰动越敏感，鲁棒性越差。
4. 非线性系统：构造李雅普诺夫函数的实战路径
4.1 从物理能量出发的构造思路
非线性系统没有通用的李雅普诺夫函数构造方法，但有一类系统特别好处理：机械系统、电路系统这类有明确物理能量表达式的。比如一个单摆系统，能量函数就是动能加势能，直接拿来做李雅普诺夫函数候选。MATLAB在这里的作用是帮你验证这个候选函数是否满足 ( \dot{V}<0 ) 。
以单摆为例，状态方程是 ( \dot{\theta}=\omega ) ，( \dot{\omega}=-\frac{g}{L}\sin\theta-\frac{b}{m}\omega ) 。取 ( V=\frac{1}{2}\omega^2+\frac{g}{L}(1-\cos\theta) ) ，这是动能加势能的形式。在MATLAB里用符号计算求 ( \dot{V} ) ：
syms theta omega g L b m real % 系统方程 dtheta = omega; domega = -g/L*sin(theta) - b/m*omega; % 李雅普诺夫函数候选 V = 0.5*omega^2 + g/L*(1-cos(theta)); % 求V对时间的导数 dV = diff(V, theta)*dtheta + diff(V, omega)*domega; dV = simplify(dV); disp('dV/dt = '); disp(dV);跑出来结果是 ( \dot{V}=-\frac{b}{m}\omega^2 ) ，在 ( b>0 ) 时负半定。负半定只能说明稳定，不能说明渐近稳定，因为 ( \dot{V}=0 ) 可能出现在 ( \omega=0 ) 但 ( \theta\neq 0 ) 的地方。要证明渐近稳定，还需要用拉萨尔不变性原理进一步分析。这个例子说明：物理能量函数是个好的起点，但不一定一步到位。
4.2 二次型加修正项的试凑方法
对于没有明显物理能量的系统，二次型 ( V=x^TPx ) 是最常用的起点。但非线性系统的 ( \dot{V} ) 里会出现高次项，二次型不一定能压住它们。这时候需要在二次型基础上加修正项，比如三次项、四次项，甚至交叉项。试凑的过程在MATLAB里可以半自动化：先固定二次型部分，然后逐项添加修正项，观察 ( \dot{V} ) 的符号是否改善。
我通常的做法是：先把 ( \dot{V} ) 的表达式展开，按幂次分组，看哪一组的符号不确定。然后针对性地添加修正项去抵消那些“坏”项。这个过程可能需要反复迭代，MATLAB的符号计算能帮你快速看到每次修改后的效果。但要注意，修正项不能破坏 ( V ) 本身的正定性，否则整个判据就失效了。
4.3 数值网格验证：符号计算之外的保险
符号计算给出的 ( \dot{V} ) 表达式，有时候因为化简规则的限制，看起来符号不确定，但实际上在某个区域内是负定的。这时候数值网格验证就派上用场了。思路很简单：在平衡点附近取一个状态空间网格，逐点计算 ( V ) 和 ( \dot{V} ) 的值，用颜色图或者等高线图展示符号分布。
% 定义网格范围 [x1, x2] = meshgrid(-2:0.1:2, -2:0.1:2); % 计算V和dV（以某个具体系统为例） V = x1.^2 + x2.^2; % 这里替换成实际的V表达式 dV = -2*x1.^2 - 2*x2.^2 + 0.5*x1.^2.*x2; % 替换成实际的dV表达式 % 可视化dV的符号 figure; surf(x1, x2, dV); xlabel('x1'); ylabel('x2'); zlabel('dV/dt'); title('dV/dt在状态空间上的分布'); % 检查是否有dV>0的区域 if any(dV(:) > 0) disp('存在dV>0的区域，需进一步分析'); else disp('所有网格点上dV<=0'); end这段代码的关键是网格范围的选择。范围太小，可能漏掉远处的发散区域；范围太大，计算量上去而且可能包含平衡点附近以外的无效区域。我的经验是：先根据系统的物理意义估计状态变量的合理范围，然后在这个范围外扩20%左右。网格步长取范围的1%到2%，太粗会漏细节，太细计算慢。
5. 完整案例：一个二阶非线性系统的分析全流程
5.1 系统描述与平衡点求解
考虑系统 ( \dot{x}_1=-x_1+x_2^2 ) ，( \dot{x}_2=-x_2-x_1x_2 ) 。先找平衡点：令两个方程为零，解得 ( x_1=0, x_2=0 ) 是唯一平衡点。在MATLAB里用solve验证：
syms x1 x2 real eq1 = -x1 + x2^2 == 0; eq2 = -x2 - x1*x2 == 0; sol = solve([eq1, eq2], [x1, x2]); disp('平衡点：'); disp([sol.x1, sol.x2]);确认平衡点后，下一步是构造李雅普诺夫函数。这个系统没有明显的物理能量，先试二次型 ( V=x_1^2+x_2^2 ) 。
5.2 构造候选函数并计算导数
syms x1 x2 real % 系统方程 f1 = -x1 + x2^2; f2 = -x2 - x1*x2; % 候选李雅普诺夫函数 V = x1^2 + x2^2; % 计算dV/dt dV = diff(V, x1)*f1 + diff(V, x2)*f2; dV = simplify(dV); disp('dV/dt = '); disp(dV);跑出来 ( \dot{V}=-2x_1^2-2x_2^2+2x_1x_2^2 ) 。这个表达式在原点附近是负定的，因为二次项主导，三次项是高阶小量。但远离原点时，三次项可能让 ( \dot{V} ) 变正。需要确定一个区域，使得在这个区域内 ( \dot{V}<0 ) 。
5.3 确定吸引域范围
要找到 ( \dot{V}<0 ) 的区域，可以令 ( \dot{V}=0 ) 解出边界。对于这个例子，( \dot{V}=-2x_1^2-2x_2^2+2x_1x_2^2=0 ) ，整理得 ( x_1^2+x_2^2=x_1x_2^2 ) 。这个方程描述了一条闭合曲线，曲线内部就是吸引域的一个估计。在MATLAB里可以数值求解这条曲线：
% 在网格上计算dV [x1, x2] = meshgrid(-3:0.05:3, -3:0.05:3); dV = -2*x1.^2 - 2*x2.^2 + 2.*x1.*x2.^2; % 找dV=0的等高线 figure; contour(x1, x2, dV, [0 0], 'r', 'LineWidth', 2); hold on; contour(x1, x2, dV, [-1 -0.5 -0.1], 'b'); xlabel('x1'); ylabel('x2'); title('dV/dt的零水平集与吸引域估计'); grid on;红色曲线是 ( \dot{V}=0 ) 的边界，蓝色曲线是 ( \dot{V} ) 取负值的等高线。从图上可以读出，吸引域大致在 ( |x_1|<1.5, |x_2|<1.5 ) 的范围内。这个估计是保守的，实际吸引域可能更大，但保守估计在工程上更安全。
5.4 时域仿真验证
数值验证的最后一步是时域仿真：从吸引域内取几个初始条件，看轨迹是否收敛到原点。
% 定义系统方程 sys = @(t, x) [-x(1) + x(2)^2; -x(2) - x(1)*x(2)]; % 取几个初始条件 init_conditions = [0.5 0.5; -0.8 0.3; 1.0 -0.5; -1.2 -1.0]; figure; hold on; for i = 1:size(init_conditions, 1) [t, x] = ode45(sys, [0 20], init_conditions(i, :)); plot(x(:,1), x(:,2), 'LineWidth', 1.5); plot(init_conditions(i,1), init_conditions(i,2), 'o', 'MarkerSize', 8); end plot(0, 0, 'k*', 'MarkerSize', 12); xlabel('x1'); ylabel('x2'); title('相轨迹与吸引域验证'); grid on;如果所有轨迹都收敛到原点，说明吸引域估计是有效的。如果有轨迹发散，说明估计偏大，需要缩小范围重新分析。这个“构造-验证-修正”的循环，是非线性李雅普诺夫分析的标准工作流。
6. 常见问题与排查技巧实录
6.1 符号计算化简不彻底怎么办
MATLAB的simplify函数有时候不能把 ( \dot{V} ) 化简到最简形式，导致符号判断困难。我的应对策略是：先用expand展开，再用collect按变量幂次收集，最后用factor尝试因式分解。如果还是不行，就代入具体的数值点做抽样检查。比如取 ( x_1=0.1, x_2=0.1 ) 代入 ( \dot{V} ) 表达式，看数值符号是否与预期一致。多个抽样点结论一致，基本可以确认符号性质。
另一个技巧是用subs把符号表达式转换成数值函数句柄，然后用fmincon找 ( \dot{V} ) 在某个区域内的最大值。如果最大值小于零，说明在这个区域内 ( \dot{V} ) 负定。这个方法比纯符号计算更可靠，代价是计算量稍大。
6.2 李雅普诺夫方程解不出来的情况
lyap函数报错或者返回NaN，通常是因为A的特征值有零实部或者正实部。零实部意味着系统临界稳定，李雅普诺夫方程可能无解或者解不唯一。这时候需要检查系统模型是否有问题，比如是不是漏掉了阻尼项。正实部则说明系统本身不稳定，解李雅普诺夫方程没有意义。
还有一种情况是A的维数很高，lyap的计算量随维数立方增长，可能跑很久。这时候可以考虑用lyapchol函数，它利用 Cholesky 分解加速，适合大规模系统。但lyapchol要求A是稳定的，用之前先确认。
6.3 吸引域估计偏保守的改进方向
二次型李雅普诺夫函数给出的吸引域估计通常偏保守，因为二次型的等高线是椭圆，而实际吸引域可能是不规则形状。改进方向有两个：一是用更高次的李雅普诺夫函数，让等高线更贴合实际边界；二是用数值方法直接计算吸引域的边界，比如反向积分法或者水平集方法。后者在MATLAB里可以用ode45反向时间积分实现，但计算量较大，适合离线分析。
我在实际项目中遇到过一个案例：二次型估计的吸引域只有实际的三分之一，导致控制器参数设计过于保守，响应速度上不去。后来换成四次型李雅普诺夫函数，吸引域估计扩大了近一倍，控制器参数可以调得更激进，响应速度明显改善。这个经验说明：李雅普诺夫函数的保守性直接影响工程设计的裕度，值得花时间优化。
6.4 数值验证中的网格选择陷阱
网格范围太大，可能把平衡点附近以外的无效区域也算进去，导致误判；网格太小，可能漏掉远处的发散区域。我的做法是：先用较大的范围粗扫一遍，确定大致边界，然后在边界附近加密网格细扫。步长从0.1开始，逐步减小到0.01，观察结论是否稳定。如果步长减小后结论变了，说明之前的网格太粗，需要继续细化。
还有一个陷阱是：网格点上的 ( \dot{V} ) 值可能因为浮点误差出现极小的正值，这时候不要急着判定系统不稳定，先看这个正值是否在数值精度范围内（比如小于1e-10）。如果是，可以认为是零，不影响结论。
7. 工具选型与代码组织建议
7.1 符号计算与数值计算的配合策略
符号计算适合推导解析表达式，数值计算适合验证和可视化。我的工作流是：先用符号计算求出 ( V ) 和 ( \dot{V} ) 的表达式，然后用matlabFunction把符号表达式转换成函数句柄，供数值计算调用。这样既保留了符号计算的精确性，又利用了数值计算的速度。
% 符号计算 syms x1 x2 real V = x1^2 + x2^2; dV = -2*x1^2 - 2*x2^2 + 2*x1*x2^2; % 转换为函数句柄 V_func = matlabFunction(V, 'Vars', {[x1; x2]}); dV_func = matlabFunction(dV, 'Vars', {[x1; x2]}); % 数值计算时直接调用 x_test = [0.5; 0.5]; disp(['V = ', num2str(V_func(x_test))]); disp(['dV = ', num2str(dV_func(x_test))]);这种配合方式在复杂系统中特别有用，因为符号表达式可能很长，直接嵌入数值循环会拖慢速度。
7.2 代码模块化与复用
李雅普诺夫分析的代码可以分成三个模块：系统定义模块、李雅普诺夫函数构造模块、验证模块。系统定义模块负责写状态方程和求平衡点；构造模块负责试凑和符号计算；验证模块负责数值网格扫描和时域仿真。每个模块写成独立的函数文件，主脚本调用它们。这样换一个系统分析时，只需要改系统定义模块，其他模块直接复用。
我自己的代码库里有一个lyapunov_analysis文件夹，里面放了这几个模块的模板。每次新项目开始时，复制一份改改参数就能用，省去了重复写代码的时间。这个习惯坚持了几年，积累下来的模板覆盖了线性系统、多项式非线性系统、三角函数非线性系统等常见类型，效率提升很明显。
7.3 可视化技巧：让稳定性一目了然
稳定性分析的结果最终要给人看，可视化质量直接影响沟通效率。我常用的三种图：相轨迹图（展示收敛行为）、( \dot{V} ) 的等高线图（展示负定区域）、( V ) 的时间演化图（展示能量衰减）。相轨迹图用plot或quiver画，等高线图用contour或surf画，时间演化图用plot画。
配色上，我习惯用红色标注不稳定区域，蓝色标注稳定区域，黑色标注平衡点。图例和坐标轴标签要写清楚，避免读者猜。如果图要放到报告里，字号至少调到12，线宽至少1.5，否则打印出来看不清。
8. 从分析到设计的延伸思考
李雅普诺夫分析不只是“判断稳定性”的工具，它还可以反过来指导控制器设计。比如，如果你希望系统以某个速率收敛，可以在 ( \dot{V} ) 的条件里加入衰减率要求：( \dot{V}\leq-\alpha V ) ，这样解出来 ( V(t)\leq V(0)e^{-\alpha t} ) ，收敛速度就有了保证。这个思路在模型预测控制和自适应控制里用得很多。
另一个延伸方向是鲁棒性分析：如果系统有不确定性或者外部扰动，李雅普诺夫函数可以加上扰动项，分析在什么扰动范围内系统仍然稳定。这个在MATLAB里可以通过robustcontrol工具箱配合实现，但核心逻辑还是李雅普诺夫第二法。
我在实际项目里的体会是：李雅普诺夫分析的价值不在于给出一个“稳定/不稳定”的二值结论，而在于它提供了一个框架，让你能定量地讨论“有多稳定”“能抗多大扰动”“收敛多快”。这些定量信息才是控制器参数整定的依据。很多人做完分析就扔到一边，只记住一个结论，那就浪费了这个工具的大部分价值。把 ( V ) 和 ( \dot{V} ) 的表达式保留下来，后续调参时随时可以拿出来重新评估，这才是正确的使用姿势。
