1. 从一个被问烂了的问题说起为什么仿真曲线收敛了系统却不一定稳定做控制或者动态系统分析的人几乎都绕不开一个场景辛辛苦苦搭好模型跑出一条响应曲线看着它慢慢趋于平缓心里一块石头落地觉得稳了。结果换一组初始条件或者把某个参数往上调了百分之十曲线直接发散到天上去。这种看起来稳、实际上不稳的翻车我在带项目和帮人看模型的时候见过太多次。问题的根子在于时域仿真只能告诉你某一条轨迹在某个特定条件下的表现它给不出对所有可能的初始状态都成立的结论。而稳定性这个概念本质上是一个全局性的、关于系统内在结构的判断。李雅普诺夫稳定性分析就是干这件事的——它不依赖你跑多少条曲线而是通过构造一个类似能量函数的标量函数从数学结构上判断平衡点附近轨迹的收敛性质。这篇内容我打算把 MATLAB 环境下做李雅普诺夫稳定性分析的完整链路讲透。核心关键词包括李雅普诺夫函数、平衡点、渐近稳定、线性矩阵不等式LMI、Lyapunov 方程、数值求解。适合正在学现代控制理论的学生、做非线性系统分析的工程师以及需要用仿真佐证理论结论的研究人员。不管你是第一次接触这个概念还是已经会用lyap函数但说不清背后逻辑下面这些内容应该都能对上你的需求。我会从最容易被忽略的平衡点求解讲起一路走到线性系统的 Lyapunov 方程数值解、非线性系统的构造技巧、基于 LMI 的自动化搜索最后落到实际项目里怎么验证和避坑。全程用 MATLAB 代码说话能跑、能改、能复现。2. 平衡点都没找对后面全是白算很多人拿到系统方程第一反应就是打开 Simulink 连线或者直接写个 ODE 求解器跑数值解。但李雅普诺夫分析的第一步根本不是仿真而是确定平衡点。稳定性永远是相对于某个平衡点而言的脱离了平衡点谈稳定就像说这个人站得稳却不说他站在哪里一样没有意义。2.1 平衡点的数学定义与常见误区对于一个连续时间自治系统 $\dot{x} f(x)$平衡点 $x_e$ 满足 $f(x_e) 0$。注意这里的关键词是自治——如果系统显含时间 $t$即 $\dot{x} f(x, t)$那平衡点的定义会复杂得多通常需要讨论的是时变解而非固定点。初学者最容易犯的错就是对着一个非自治系统硬套自治系统的结论。另一个高频误区是默认原点就是平衡点。线性系统 $\dot{x} Ax$ 确实永远有原点这个平衡点但非线性系统未必。比如 $\dot{x} x^2 - 1$平衡点在 $x \pm 1$原点根本不是平衡点。如果你对着这个系统在原点附近做线性化得到的雅可比矩阵是 $2x|_0 0$线性化完全失效这时候还硬用线性方法判断结论必然是错的。在 MATLAB 里求平衡点符号计算是最直接的路子。用 Symbolic Math Toolbox 把方程写出来solve一把梭syms x1 x2 real f1 x2; f2 -sin(x1) - 0.5*x2; eq [f1 0, f2 0]; sol solve(eq, [x1, x2], ReturnConditions, true); disp(sol.x1) disp(sol.x2)这段代码对应的是一个带阻尼的单摆系统。解出来你会看到平衡点是 $(k\pi, 0)$ 这一族$k$ 取整数。物理上很好理解摆锤竖直向下和竖直向上都是平衡位置但显然一个稳一个不稳。这就是为什么必须先找全平衡点再逐个分析——你不能只盯着原点看。2.2 数值求解平衡点的实用技巧符号解不出来的情况太常见了尤其是系统稍微复杂一点solve直接返回一个空结构体或者一堆root占位符。这时候得转向数值方法。fsolve是 MATLAB 里最常用的非线性方程数值求解器但它有个脾气对初值极其敏感。我的经验是先用相图或者网格扫描大致定位平衡点所在的区域再把这个区域里的点作为fsolve的初值。下面这段代码演示了网格扫描加精修的流程f (x) [x(2); -sin(x(1)) - 0.5*x(2)]; [X1, X2] meshgrid(-2*pi:0.5:2*pi, -3:0.5:3); Fnorm arrayfun((a,b) norm(f([a;b])), X1, X2); % 找出函数范数接近零的网格点作为候选初值 candidates [X1(Fnorm 0.1), X2(Fnorm 0.1)]; equilibria []; for k 1:size(candidates, 1) [xe, ~, exitflag] fsolve(f, candidates(k,:), ... optimoptions(fsolve, Display, off, TolFun, 1e-10)); if exitflag 0 equilibria [equilibria; xe]; end end % 去重 equilibria uniquetol(equilibria, 1e-4, ByRows, true); disp(equilibria)这里有个细节值得说uniquetol的去重容差不能设得太小否则同一个平衡点因为数值误差会被当成好几个也不能太大否则相邻的平衡点会被合并。我一般取1e-4到1e-3之间具体看系统尺度。如果系统变量量级差异很大比如一个变量是角度、另一个是电流最好先做归一化再扫描。提示fsolve的exitflag大于零只代表收敛到一个解不代表这个解就是你要的平衡点。一定要回代验证 $|f(x_e)|$ 是否足够小我通常要求残差范数小于 $10^{-8}$ 才认账。3. 线性系统的李雅普诺夫方程从理论到 lyap 函数的落地线性时不变系统是李雅普诺夫分析里最友好的一类因为存在一个充要条件而且这个条件可以化成一个标准的矩阵方程来求解。这部分是整篇内容的地基地基打不牢后面非线性那套根本没法展开。3.1 Lyapunov 方程的来龙去脉考虑线性系统 $\dot{x} Ax$。我们想找一个二次型李雅普诺夫函数 $V(x) x^T P x$其中 $P$ 是正定对称矩阵。沿着系统轨迹对 $V$ 求时间导数$$\dot{V}(x) \dot{x}^T P x x^T P \dot{x} x^T (A^T P P A) x$$要让 $\dot{V}(x)$ 负定就需要 $A^T P P A -Q$其中 $Q$ 是任意正定对称矩阵。这个方程就是大名鼎鼎的Lyapunov 方程。理论结论很干净系统 $\dot{x} Ax$ 渐近稳定的充要条件是对任意给定的正定对称 $Q$Lyapunov 方程存在唯一正定对称解 $P$。实际用的时候为了省事通常直接取 $Q I$因为如果 $Q I$ 时解不出正定的 $P$换别的 $Q$ 也救不回来——这个结论的证明依赖于 Lyapunov 方程解对 $Q$ 的连续依赖性这里不展开记住结论就行。3.2 lyap 与 dlyap 的正确打开方式MATLAB 的 Control System Toolbox 提供了lyap函数直接求解连续时间 Lyapunov 方程。用法简单到令人发指A [0 1; -2 -3]; Q eye(2); P lyap(A, Q); % 注意lyap 求解的是 A*P P*A Q 0 disp(P) eig(P)这里有个极其容易踩的坑MATLAB 的lyap函数定义的方程形式是 $AP PA^T Q 0$而不是我们理论推导里的 $A^T P PA Q 0$。两者差了一个转置。所以上面代码里我写的是lyap(A, Q)把 $A$ 转置一下传进去得到的 $P$ 才对应我们想要的 $A^T P PA -Q$。如果你不转置直接写lyap(A, Q)解出来的 $P$ 是另一个方程的解虽然它可能也是正定的但对应的李雅普诺夫函数 $V x^T P x$ 对原系统未必成立。这个错误我在审别人代码的时候至少见过五次而且因为结果看起来正常特别难被发现。验证环节不能省。解出 $P$ 之后一定要做两件事一是检查 $P$ 是否对称正定用eig(P)看特征值是否全正二是回代验证残差residual A*P P*A Q; disp(norm(residual)) % 应该接近机器精度如果残差在 $10^{-10}$ 量级以下说明求解没问题。如果残差很大要么是矩阵接近奇异要么是系统本身不稳定导致方程病态。3.3 离散系统的 dlyap 与采样周期陷阱离散系统 $\dot{x}$ 换成 $x_{k1} A_d x_k$李雅普诺夫方程变成 $A_d^T P A_d - P -Q$对应的函数是dlyap。这里有个隐蔽的坑连续系统稳定离散化之后不一定稳定这取决于采样周期的选择。我做过一个测试一个连续系统极点实部是 $-0.5$用零阶保持器离散化采样周期从 $0.01$ 秒一路加到 $2$ 秒。当采样周期超过某个阈值时离散系统的极点会跑到单位圆外dlyap解出来的 $P$ 就不再正定了。这个阈值和系统的时间常数直接相关经验上采样频率至少要达到系统带宽的十倍以上才比较保险。A [0 1; -2 -3]; B [0; 1]; for Ts [0.01, 0.1, 0.5, 1.0, 2.0] Ad expm(A*Ts); try P dlyap(Ad, eye(2)); fprintf(Ts%.2f, min eig(P)%.4f\n, Ts, min(eig(P))); catch fprintf(Ts%.2f, dlyap 求解失败\n, Ts); end end跑一遍你会发现随着Ts增大min(eig(P))逐渐减小甚至变负。这个实验比任何理论推导都更能让人记住采样周期的重要性。4. 非线性系统的李雅普诺夫函数没有万能公式但有章法到了非线性这块事情就变得艺术起来了。线性系统有充要条件、有现成函数非线性系统没有通用的构造方法。李雅普诺夫第二法直接法只告诉你如果存在这样的函数那么稳定但怎么找这个函数它不管。这是很多人卡住的地方。4.1 克拉索夫斯基方法用雅可比矩阵碰运气克拉索夫斯基方法算是最机械化的一种构造思路。它的核心是如果雅可比矩阵 $J(x) \partial f / \partial x$ 在某个区域内满足 $J(x) J^T(x)$ 负定那么系统在该区域内渐近稳定而且李雅普诺夫函数可以直接取 $V(x) f^T(x) f(x)$。这个方法的好处是不需要猜雅可比矩阵是现成的。坏处是条件太强很多系统满足不了。但对于一些结构比较规整的系统它往往能快速给出结论。syms x1 x2 real f [x2; -x1 - x2^3]; J jacobian(f, [x1, x2]); S J J; % 检查 S 是否负定 eigS eig(S); disp(eigS)对于这个系统$S$ 的特征值一个是 $-2$另一个是 $-2 - 6x_2^2$在整个平面上都负所以全局渐近稳定$V f^T f$ 就是一个合法的李雅普诺夫函数。但如果你把阻尼项改成 $-x_2$$S$ 的特征值就会出现与 $x_2$ 相关的项在某些区域可能变正克拉索夫斯基方法就失效了。4.2 变量梯度法与待定系数手工构造的实战套路当克拉索夫斯基方法不灵时变量梯度法是另一个常用手段。它的思路是假设 $\dot{V}$ 的梯度形式然后通过积分还原出 $V$同时要求旋度条件满足保证积分与路径无关。这个方法在教科书里有详细推导但实际手算起来很繁琐。我在项目里更常用的是待定系数法根据系统的物理意义猜一个 $V$ 的形式比如机械系统用动能加势能电路系统用电容储能加电感储能然后代入 $\dot{V}$ 验证不满足就调整系数。这个过程在 MATLAB 里可以用符号计算加速syms x1 x2 a b c real V a*x1^2 b*x1*x2 c*x2^2; f1 x2; f2 -x1 - x2^3; Vdot diff(V, x1)*f1 diff(V, x2)*f2; Vdot simplify(Vdot); % 收集 x1^2, x1*x2, x2^2, x2^4 等项的系数 collect(Vdot, [x1, x2])把Vdot展开后你会得到关于 $a, b, c$ 的一组约束。目标是让所有项的系数都非正或者负定。对于这个例子取 $a 1, b 0, c 1$ 就能让 $\dot{V} -2x_2^4$虽然只是半负定但结合拉萨尔不变性原理仍然可以推出渐近稳定。注意$\dot{V}$ 半负定不等于渐近稳定必须配合拉萨尔不变集定理。很多人在这里直接下结论说渐近稳定是错的。拉萨尔定理要求不变集里除了原点没有别的轨迹这个条件要单独验证。4.3 基于 LMI 的自动化搜索让求解器帮你找 P对于可以写成线性矩阵不等式形式的稳定性问题MATLAB 的 Robust Control Toolbox 或者 YALMIP 工具箱能帮大忙。核心思想是把找正定矩阵 $P$ 使得某个矩阵不等式成立这件事交给凸优化求解器。以线性系统为例稳定性条件 $A^T P PA 0, P 0$ 本身就是一个 LMI 可行性问题。用 YALMIP 写出来是这样的% 需要安装 YALMIP 和任一 SDP 求解器如 SeDuMi、SDPT3 A [0 1; -2 -3]; P sdpvar(2, 2); Constraints [P eye(2)*1e-6, A*P P*A -eye(2)*1e-6]; optimize(Constraints, [], sdpsettings(solver, sedumi, verbose, 0)); if double(Constraints(1)) Pval double(P); disp(Pval) disp(eig(Pval)) endLMI 方法的威力在于它能处理不确定性和时变参数。比如系统矩阵 $A$ 在一个凸多面体里变化你可以对每个顶点写一个 LMI求解器会找一个对所有顶点都成立的公共 $P$。这是线性系统lyap函数做不到的。不过 LMI 也不是银弹。求解器的数值精度、问题的规模、约束的冗余程度都会影响结果。我遇到过规模稍大的问题状态维数超过 20求解器直接报数值错误这时候要么降维要么换用交替方向乘子法之类的分布式算法。5. 仿真验证把理论结论和数值实验对上理论推导给出稳定的结论之后必须用仿真验证。这不是走过场而是因为理论推导的假设在数值实现里可能被破坏比如线性化忽略了高阶项、离散化引入了额外动态、求解器精度不够等等。5.1 相图与向量场最直观的稳定性可视化对于二维系统相图是最有说服力的验证工具。MATLAB 里画相图可以用quiver画向量场再叠加几条从不同初始条件出发的轨迹f (t, x) [x(2); -x(1) - x(2)^3]; [X1, X2] meshgrid(-3:0.3:3, -3:0.3:3); U X2; V -X1 - X2.^3; figure; quiver(X1, X2, U, V, Color, [0.7 0.7 0.7]); hold on; for x0 [-2 -1; -1 2; 1 -2; 2 1; 0.5 0.5] [~, X] ode45(f, [0 20], x0); plot(X(:,1), X(:,2), LineWidth, 1.5); end xlabel(x_1); ylabel(x_2); title(相图与轨迹); grid on;看相图的时候重点观察三件事轨迹是否都朝原点汇聚、有没有闭合轨道极限环、有没有轨迹跑到无穷远。如果所有轨迹都收敛到原点那渐近稳定的结论就站得住。5.2 李雅普诺夫函数的等高线与导数场光看轨迹收敛还不够最好把李雅普诺夫函数本身也画出来看看它的等高线是不是包裹着轨迹以及 $\dot{V}$ 在轨迹上是不是单调递减。V (x1, x2) x1.^2 x2.^2; Vdot (x1, x2) -2*x2.^4; [X1, X2] meshgrid(-2:0.1:2, -2:0.1:2); contour(X1, X2, V(X1, X2), 20, LineWidth, 1); hold on; % 叠加一条轨迹 [~, X] ode45(f, [0 15], [1.5; -1]); plot(X(:,1), X(:,2), r, LineWidth, 2); % 检查 Vdot 沿轨迹的符号 Vdot_along Vdot(X(:,1), X(:,2)); fprintf(Vdot 最大值: %.4f\n, max(Vdot_along));如果Vdot沿轨迹的最大值小于零或者等于零但只在原点取到那渐近稳定的结论就得到了数值支持。这个检查比单纯看轨迹收敛更严格因为它直接验证了李雅普诺夫函数的单调性。5.3 数值精度与刚性系统的处理有个坑我必须单独拎出来说刚性系统用ode45跑出来的结果可能完全不可信。刚性系统的特征是时间常数差异极大ode45为了满足精度会疯狂缩小步长跑到天荒地老或者直接报步长过小的警告。判断系统是否刚性可以看雅可比矩阵特征值的实部比值。如果最大实部和最小实部的比值超过 $10^3$基本就是刚性系统了。这时候应该换用ode15s或ode23sopts odeset(RelTol, 1e-8, AbsTol, 1e-10); [t1, X1] ode45(f, [0 20], [1; 1], opts); [t2, X2] ode15s(f, [0 20], [1; 1], opts); fprintf(ode45 步数: %d, ode15s 步数: %d\n, length(t1), length(t2));对于刚性系统ode15s的步数往往比ode45少一两个数量级而且结果更可靠。我见过有人用ode45跑刚性系统曲线看起来收敛了实际上是因为数值耗散把真实的不稳定给抹平了换成ode15s立刻发散。这种假收敛最害人。6. 几个让我印象深刻的翻车现场与排查思路理论讲完了代码也给了最后这部分我想聊几个实际踩过的坑。这些坑在教科书里不会写但每一个都能让你白干好几天。6.1 线性化点选错导致的假稳定有个项目里系统在原点附近线性化后极点都在左半平面lyap解出来的 $P$ 正定一切看起来完美。但实际仿真时从稍远一点的初始条件出发系统直接发散。排查了半天才发现原点的稳定域非常小线性化只在原点附近的一个极小邻域内有效稍微远一点高阶项就主导了动态。这个教训是线性化分析给出的稳定性是局部的稳定域的大小必须单独估计。估计方法可以用反向轨迹法或者用李雅普诺夫函数的等高线找最大的不变集。在 MATLAB 里可以数值搜索使 $\dot{V} 0$ 成立的最大区域% 估计稳定域找 V(x) c 中最大的 c 使得该等高线内 Vdot 0 c_values linspace(0.1, 10, 100); max_c 0; for c c_values % 在 V(x) c 的等高线上采样 theta linspace(0, 2*pi, 200); % 这里以 V x1^2 x2^2 为例 x1 sqrt(c)*cos(theta); x2 sqrt(c)*sin(theta); Vdot_vals arrayfun((a,b) Vdot(a,b), x1, x2); if all(Vdot_vals 0) max_c c; else break; end end fprintf(估计的稳定域半径平方: %.4f\n, max_c);6.2 数值求解 Lyapunov 方程时的病态问题当系统矩阵 $A$ 的特征值实部有正有负或者非常接近虚轴时Lyapunov 方程会变得病态。lyap函数可能返回一个条件数极大的 $P$虽然特征值都是正的但数值上极不稳定。判断病态的方法是看 $P$ 的条件数P lyap(A, eye(size(A))); cond_P cond(P); fprintf(P 的条件数: %.2e\n, cond_P);如果条件数超过 $10^{10}$解出来的 $P$ 基本不可信。这时候可以尝试用lyapchol函数它基于 Cholesky 分解数值稳定性更好但要求系统必须是稳定的不稳定系统直接报错。另一个办法是手动缩放系统矩阵把状态变量归一化到相近的量级。6.3 离散化引入的伪稳定与伪不稳定前面提过采样周期的影响这里再补一个更隐蔽的情况用c2d离散化时离散化方法的选择会影响稳定性判断。零阶保持器和双线性变换Tustin给出的离散系统极点位置不同在某些采样周期下一个方法判断稳定另一个可能判断不稳定。我的做法是对同一个连续系统用多种离散化方法分别算一遍如果结论一致那基本可信如果出现分歧说明采样周期选得不合适需要减小采样周期重新评估。A [0 1; -10 -1]; for method {zoh, tustin, matched} for Ts [0.01, 0.05, 0.1] Ad c2d(ss(A, [], [], []), Ts, method{1}).A; stable all(abs(eig(Ad)) 1); fprintf(方法%s, Ts%.2f, 稳定%d\n, method{1}, Ts, stable); end end跑一遍这个代码你会看到不同方法在不同采样周期下的稳定性判断确实会有差异。这个实验能帮你建立起对采样周期和离散化方法的直觉。6.4 李雅普诺夫函数存在但找不到的尴尬最后说一个哲学层面的坑李雅普诺夫第二法是充分条件不是必要条件。系统稳定不代表你能找到一个简单的多项式形式的李雅普诺夫函数。有些系统的李雅普诺夫函数是非多项式的、甚至是不可微的。遇到这种情况不要死磕。可以退而求其次用数值方法验证稳定域或者用反例搜索来证伪。如果实在需要严格的稳定性证明可以考虑用平方和SOS优化它能搜索多项式形式的李雅普诺夫函数比手工构造强大得多但计算成本也高得多。% SOS 方法示意需要 SOSTOOLS 工具箱 % 这里只给出思路具体语法参考 SOSTOOLS 文档 % 目标找多项式 V(x) 使得 V - epsilon*(xx) 是 SOS且 -Vdot 是 SOSSOS 方法在近十年发展很快对于中小规模的多项式系统它往往能找到手工方法找不到的李雅普诺夫函数。但它的局限也很明显只适用于多项式系统且随着状态维数增加计算复杂度呈指数增长。7. 我个人的一点使用心得李雅普诺夫稳定性分析这套东西理论优美但落地的时候处处是细节。我自己的习惯是先数值后理论先仿真后证明。拿到一个系统先用ode45或者ode15s跑几条轨迹看看大致行为心里有个底然后找平衡点、做线性化、用lyap或 LMI 求解最后再回到仿真验证理论结论。这个顺序比反过来效率高得多因为仿真能帮你快速排除掉那些明显不稳定的情况省得在理论上白费功夫。另外MATLAB 的符号计算和数值计算要配合着用。符号计算帮你推导公式、验证手算结果数值计算帮你处理大规模问题和实际数据。两者不是替代关系而是互补关系。我见过有人非要用符号计算解一个十维系统的 Lyapunov 方程结果内存直接爆掉也见过有人对着一个简单的二阶系统硬写数值迭代明明solve一行就能出结果。最后提醒一句任何稳定性结论都要有验证环节。不管是理论推导还是数值求解最后都要回到仿真或者实验上确认。理论告诉你应该稳定仿真告诉你实际稳定两者对上你才能放心。对不上那就是有假设被违反了得回去查。这个查的过程往往比结论本身更有价值。