企业官网建设流程全解析

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_{k+1} = 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)) end

LMI 方法的威力在于它能处理不确定性和时变参数。比如系统矩阵 $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或ode23s：

opts = 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*(x'x) 是 SOS，且 -Vdot 是 SOS

SOS 方法在近十年发展很快，对于中小规模的多项式系统，它往往能找到手工方法找不到的李雅普诺夫函数。但它的局限也很明显：只适用于多项式系统，且随着状态维数增加，计算复杂度呈指数增长。

7. 我个人的一点使用心得

李雅普诺夫稳定性分析这套东西，理论优美，但落地的时候处处是细节。我自己的习惯是：先数值后理论，先仿真后证明。拿到一个系统，先用ode45或者ode15s跑几条轨迹看看大致行为，心里有个底；然后找平衡点、做线性化、用lyap或 LMI 求解；最后再回到仿真验证理论结论。这个顺序比反过来效率高得多，因为仿真能帮你快速排除掉那些明显不稳定的情况，省得在理论上白费功夫。

另外，MATLAB 的符号计算和数值计算要配合着用。符号计算帮你推导公式、验证手算结果，数值计算帮你处理大规模问题和实际数据。两者不是替代关系，而是互补关系。我见过有人非要用符号计算解一个十维系统的 Lyapunov 方程，结果内存直接爆掉；也见过有人对着一个简单的二阶系统硬写数值迭代，明明solve一行就能出结果。

最后提醒一句：任何稳定性结论都要有验证环节。不管是理论推导还是数值求解，最后都要回到仿真或者实验上确认。理论告诉你"应该稳定"，仿真告诉你"实际稳定"，两者对上，你才能放心。对不上，那就是有假设被违反了，得回去查。这个查的过程，往往比结论本身更有价值。