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^TPPA-Q ) 如果得到的 ( P ) 正定则系统渐近稳定。MATLAB里用lyap函数直接求解语法是P lyap(A, Q)。注意这里的符号约定MATLAB的lyap求解的是 ( APPA^TQ0 ) 和教科书上的写法可能差一个转置用之前最好用简单例子验证一下方向。我通常取 ( QI ) 单位阵因为这样计算最省事而且不影响判据的结论——只要存在一个正定的 ( 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(k1)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 ) 在 ( b0 ) 时负半定。负半定只能说明稳定不能说明渐近稳定因为 ( \dot{V}0 ) 可能出现在 ( \omega0 ) 但 ( \theta\neq 0 ) 的地方。要证明渐近稳定还需要用拉萨尔不变性原理进一步分析。这个例子说明物理能量函数是个好的起点但不一定一步到位。4.2 二次型加修正项的试凑方法对于没有明显物理能量的系统二次型 ( Vx^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在状态空间上的分布); % 检查是否有dV0的区域 if any(dV(:) 0) disp(存在dV0的区域需进一步分析); else disp(所有网格点上dV0); end这段代码的关键是网格范围的选择。范围太小可能漏掉远处的发散区域范围太大计算量上去而且可能包含平衡点附近以外的无效区域。我的经验是先根据系统的物理意义估计状态变量的合理范围然后在这个范围外扩20%左右。网格步长取范围的1%到2%太粗会漏细节太细计算慢。5. 完整案例一个二阶非线性系统的分析全流程5.1 系统描述与平衡点求解考虑系统 ( \dot{x}_1-x_1x_2^2 ) ( \dot{x}_2-x_2-x_1x_2 ) 。先找平衡点令两个方程为零解得 ( x_10, x_20 ) 是唯一平衡点。在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]);确认平衡点后下一步是构造李雅普诺夫函数。这个系统没有明显的物理能量先试二次型 ( Vx_1^2x_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^22x_1x_2^2 ) 。这个表达式在原点附近是负定的因为二次项主导三次项是高阶小量。但远离原点时三次项可能让 ( \dot{V} ) 变正。需要确定一个区域使得在这个区域内 ( \dot{V}0 ) 。5.3 确定吸引域范围要找到 ( \dot{V}0 ) 的区域可以令 ( \dot{V}0 ) 解出边界。对于这个例子( \dot{V}-2x_1^2-2x_2^22x_1x_2^20 ) 整理得 ( x_1^2x_2^2x_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; % 找dV0的等高线 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_10.1, x_20.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} ) 的表达式保留下来后续调参时随时可以拿出来重新评估这才是正确的使用姿势。