返回课程总览

实验 11 · 理论工具

稳定性分析的数值工具

检查正定性、Lyapunov 差分、MDADT 子区间检查和 Schur 补，并认识 LMI 工作流。

完整讲义

Lab 11：稳定性分析的数值工具

前几个实验关注“模型是否按公式运行”。本实验直接读取 Lab 10 保存的异步故障轨迹，转向 稳定性分析中常见的数学对象：正定矩阵、Lyapunov 函数、切换计数、驻留时间、Schur 补和 线性矩阵不等式（LMI）。目标是读懂并检查这些对象，不要求完成控制器设计中的复杂 LMI 求解。

开始前复用 Lab 4 的 Git 检查点流程：用 git status 确认工作区干净；正反例、MDADT 子区间检查和输出表都通过后，按相同方法提交 Lab 11 检查点。

学习目标

  • 用特征值检查对称矩阵的正定性。
  • 沿数值轨迹计算 \(V(k)=\boldsymbol{x}(k)^{\mathsf T}P\boldsymbol{x}(k)\) 与 \(\Delta V(k)\)。
  • 从模式序列提取驻留区间，并检查基于模式的平均驻留时间（MDADT）。
  • 在低维例子中核对 Schur 补等价关系。
  • 理解 LMI 从理论条件到控制器增益的大致工作流。

目录结构

labs/lab10/results/
  analysis-input.mat
labs/lab11/
  stability_analysis.m
  countModeDwells.m
  checkMdadt.m
  results/
    lyapunov-history.csv
    dwell-summary.csv
    mdadt-check.csv
  figures/
    lyapunov-history.png
  notes.md
小知识：数值轨迹不能代替稳定性证明

一条轨迹上的 \(V(k)\) 下降，只说明这个初值和这段切换序列的数值表现。稳定性定理通常要求结论对一类初值、不确定性和允许切换信号都成立。仿真适合发现错误、展示现象和支持算例，但不能替代理论中的量词与不等式推导。

理论问题

正定矩阵、Lyapunov 函数与下降条件

给定对称矩阵：

\[ P=\begin{bmatrix}2.0&0.3\\0.3&1.0\end{bmatrix}. \]

若 \(P\) 的全部特征值为正，则 \(P\) 正定。二次型 \(V(\boldsymbol{x})=\boldsymbol{x}^{\mathsf T}P\boldsymbol{x}\) 可以看成一种加权状态大小： \(P\) 不仅为不同状态分量分配权重，也能通过非对角项表达分量之间的耦合。它满足：

\[ V(k)=\boldsymbol{x}(k)^{\mathsf T}P\boldsymbol{x}(k)>0 \quad\text{for }\boldsymbol{x}(k)\ne\boldsymbol{0},\qquad \Delta V(k)=V(k+1)-V(k). \]

沿某条轨迹计算 \(\Delta V(k)\) 是数值检查。理论证明则要从系统矩阵和允许条件出发，说明所需不等式对所有相关状态成立。

对一次固定的普通更新 \(\boldsymbol{x}(k+1)=A_q\boldsymbol{x}(k)\)，把更新式代入 \(V\) 可得：

\[ \Delta V=\boldsymbol{x}^{\mathsf T}(A_q^{\mathsf T}PA_q-P)\boldsymbol{x}. \]

只有 \(P\succ0\) 还不够；若 \(A_q^{\mathsf T}PA_q-P\prec0\)，也就是括号中的矩阵负定，才能保证这一次固定模型对所有非零状态 都有 \(\Delta V<0\)。这说明“\(V\) 是正的”和“\(V\) 沿系统下降”是两个不同条件。

本实验使用同一个 \(P\) 计算整条轨迹，相当于检查一个公共 Lyapunov 矩阵的候选结果。更一般 的切换系统理论可以为每个模式选择模式相关矩阵 \(P_i\)，并使用 \(V_i(\boldsymbol{x})=\boldsymbol{x}^{\mathsf T}P_i\boldsymbol{x}\)。这种选择可能降低单个矩阵 同时适合所有模式的限制，但切换时 Lyapunov 函数本身也会改变，必须额外约束模式之间的增长。

脉冲映射、控制器异步和模式切换都可能让 \(V\) 在个别时刻上升。稳定性分析不能直接删掉这些 上升点，而要把每段流动中的下降与事件处的增长放在一起估计。驻留时间条件限制切换发生的 频率或各模式累计停留时间，从而控制累计切换效应；这正是下一小节引入 MDADT 的原因。

基于模式的平均驻留时间

本实验采用基于模式的平均驻留时间（mode-dependent average dwell time, MDADT），分别记录 每个模式的进入次数和运行时间，而不是把所有切换合并成一个计数。控制理论文献中的连续时间 定义常用闭区间记号描述各模式的激活次数和总运行时间；本实验采用以下离散半开区间约定， 把长度为 \(N\) 的模式序列视为 \([0,N)\) 上的 \(N\) 个更新时刻。对 \([k_0,k_1)\) 和模式 \(i\)，定义：

  • \(T_i(k_0,k_1)\)：序列在该区间内处于模式 \(i\) 的更新步数；
  • \(N_i(k_0,k_1)\)：区间内部进入模式 \(i\) 的次数。

“区间内部进入”表示只计数满足 \(k_0<k<k_1\)、\(\sigma(k)=i\) 且 \(\sigma(k-1)\ne i\) 的位置。即使 \(\sigma(k_0)=i\)，区间起点的初始模式也不计为一次进入。 MDADT 条件逐模式写为：

\[ N_i(k_0,k_1)\le N_{0i}+\frac{T_i(k_0,k_1)}{\tau_i}. \]

\(N_{0i}\) 是模式 \(i\) 的抖振界，\(\tau_i\) 是该模式的平均驻留时间参数。定理中的 “对任意区间”不能被一次全区间检查替代，因此程序必须枚举模式序列的所有半开子区间 \([k_0,k_1)\)，其中 \(0\le k_0<k_1\le N\)。

Lab 10 的模式序列是长度为 15 的交替模式块。使用两组参数：

参数集 \(N_{0i}\) \(\tau_i\) 预期
expected-pass [1,1] [15,15] 每个模式的所有子区间都通过
full-only-trap [1,1] [20,20] 全区间通过，但至少一个子区间失败

第二组故意揭示“只检查全区间”的缺陷。例如区间 \([14,46)\) 内，模式 2 的进入次数为 2， 驻留步数为 16，因此 \(2>1+16/20\)；但完整 60 步区间对两个模式都通过。

Schur 补与 LMI

稳定性推导中常出现逆矩阵或二次矩阵项，它们不便直接交给 LMI 求解器。Schur 补可以在满足 正定前提时，把这类条件改写成等价的分块矩阵条件。若分块中的元素对待求矩阵变量是仿射的， 这个分块条件就是线性矩阵不等式（LMI），可以由数值优化软件检查可行性。

对对称分块矩阵：

\[ S=\begin{bmatrix}Q&R\\R^{\mathsf T}&W\end{bmatrix}. \]

当 \(W\) 正定时，\(S\) 正定等价于：

\[ W\succ0,\qquad Q-RW^{-1}R^{\mathsf T}\succ0. \]

实际代码不应显式计算逆矩阵；使用右除或线性方程求解更合适。本实验只给定数值矩阵并检查 等价关系，帮助理解如何从一个矩阵不等式得到分块 LMI；它不求解控制器，也不代替引理证明。

公式到代码

先读取 Lab 10 的固定输出，不把变量直接散落到当前工作区：

inputPath = fullfile(projectRoot, "..", "lab10", "results", "analysis-input.mat");
loaded = load(inputPath, "analysisInput");
x = loaded.analysisInput.x;
modeSequence = loaded.analysisInput.plantModes;

对应的课程路径是 labs/lab10/results/analysis-input.mat。如果文件不存在，应先完成并运行 Lab 10；本实验不存在替代示例轨迹，因为使用另一组数据会切断两次实验之间的验证关系。

正定性和 Lyapunov 历史：

assert(norm(P - P.', "fro") < 1e-12);
eigenvalues = eig(P);
assert(all(eigenvalues > 0));

V = sum(x .* (P * x), 1);
deltaV = diff(V);

Schur 补检查：

schurComplement = Q - R / W * R.';
assert(all(eig(W) > 0));
assert(all(eig(S) > 0));
assert(all(eig(schurComplement) > 0));

QBad = diag([0.01, 1]);
SBad = [QBad, R; R.', W];
schurComplementBad = QBad - R / W * R.';
assert(~all(eig(SBad) > 0));
assert(~all(eig(schurComplementBad) > 0));

countModeDwells.m 接收模式序列，返回每个连续区间的模式、开始索引、结束索引和长度。让数据结构表达数学定义，比在主脚本里散落多个计数变量更容易审查。

checkMdadt.m 的接口固定为：

resultTable = checkMdadt(modeSequence, chatterBounds, dwellTimes, parameterSetName);

返回表每行对应一个模式，列为 parameterSet、mode、chatterBound、 averageDwellTime、fullIntervalPass、allIntervalsPass、minimumMargin、 firstViolationStart 和 firstViolationEnd。违反量使用 \(N_{0i}+T_i/\tau_i-N_i\) 的最小值；无违反时，首次违反区间的两个端点记为 NaN。 违反区间端点使用公式中的零起始半开区间索引；例如输出 14 和 46 表示 \([14,46)\)。

非正定测试不应故意让主脚本中止。先单独构造反例，再断言它被正确拒绝：

PInvalid = P;
PInvalid(2, 2) = 0;
assert(~all(eig(PInvalid) > 0));

Schur 补同时测试正例和反例。正例保留给定的 Q、R、W；反例使用相同 R 与 W=1.2>0，只把 QBad=diag([0.01,1])。代码要明确断言正定一侧的两个判定都为真， 反例一侧的两个判定都为假，而不是只比较两个布尔值是否碰巧相等。

MATLAB 实验

  1. 运行 Lab 10，确认 analysis-input.mat 已生成；再由 Starter 中的固定路径加载它。
  2. 检查 P 的对称性与最小特征值。
  3. 计算 \(V(k)\)、\(\Delta V(k)\)，保存 lyapunov-history.csv。
  4. 用 countModeDwells 生成 dwell-summary.csv，人工核对前两个驻留区间。
  5. 调用 checkMdadt 检查 expected-pass 和 full-only-trap，把两张结果表合并后写入 mdadt-check.csv。
  6. 断言第一组所有子区间通过；再断言第二组全区间通过、但至少一个模式的所有子区间判定失败。
  7. 构造一个满足和一个不满足正定条件的 Schur 补例子，明确核对两侧真值。
  8. 绘制 \(V(k)\)，用标记指出模式切换位置。

若轨迹中的 \(\Delta V\) 偶尔为正，不要删除数据。切换系统的多 Lyapunov 分析可能允许切换时增长，但需要额外条件控制累计效应；本实验只记录现象并准确限定结论。

数学验证

完成以下对应检查：

数学对象 数值检查 不能由此推出
\(P\succ0\) 对称且最小特征值为正 \(P\) 自动适合任意系统
\(\Delta V(k)\) 对已生成轨迹逐点计算 对所有初值和切换信号下降
MDADT 给定模式序列的所有半开子区间 其他可能序列也满足约束
Schur 补 两侧特征值判定一致 已证明一般形式的引理

再用 PInvalid(2,2)=0 构造非正定矩阵，并用否定断言确认它被识别。这样验证反例而不会故意 中止主脚本。Schur 正例中，分块矩阵和 Schur 补的正定判定都应为真；QBad 反例中，两者 都应为假。

结果与记录

notes.md 建立一个“符号到程序变量”表，例如 \(P\) 对应 P、\(V(k)\) 对应 V(k+1) 数组元素、\(N_i\) 对应 entryCount、\(T_i\) 对应 timeInMode、\(\tau_i\) 对应 averageDwellTime。然后用一段文字区分三类材料：

  • 理论条件：一般矩阵不等式和适用假设。
  • 数值检查：对指定矩阵与轨迹的特征值、残差和计数。
  • 算例图表：帮助读者理解某组参数下的表现。

控制理论中的 LMI 工作流通常是：选择 Lyapunov 变量和辅助变量，按定理构造矩阵不等式， 交给求解器寻找可行解，再从变量恢复控制器增益 K。有些变换会先求 K L 一类变量， 再由非奇异矩阵 L 恢复 K。范数有界不确定性也常通过辅助标量和矩阵不等式处理。 本课程只解释这条链路，不安装优化工具箱，也不完成完整的 LMI 求解。

OpenCode 审计

复用 Lab 4 的工作流。输入 /models，确认当前模型是 DeepSeek V4 Pro，再按 Tab 进入 Plan 模式并输入：

请检查 labs/lab11/stability_analysis.m。

请把每个数值检查对应到正定性、Lyapunov 差分、MDADT 或 Schur 补的定义。
重点检查：
1. P 是否先检查对称性再判断特征值；
2. V 和 deltaV 的时间索引是否一致；
3. countModeDwells.m 是否漏掉最后一个驻留区间；
4. checkMdadt.m 是否排除区间起点的初始模式，并检查所有半开子区间；
5. expected-pass 是否全部通过，full-only-trap 是否只在子区间检查中暴露失败；
6. Schur 补代码是否避免显式 inv，并分别断言正例为真、反例为假；
7. 状态和模式是否确实来自 Lab 10 的 analysisInput；
8. notes.md 是否把单条数值轨迹误写成稳定性证明。

不要添加 Robust Control Toolbox 或 Optimization Toolbox 依赖。

若 OpenCode 建议通过放宽断言来让错误输入继续运行，应拒绝。这里的断言用于维护数学前提， 不是需要隐藏的异常。接受计划后按 Tab 进入 Build 模式。修改完成后，在第二个 PowerShell 终端运行 git status 和 git diff -- labs\lab11，并重新执行正例与反例。

完成标准

  • 正定和非正定测试都产生预期判定。
  • Lyapunov 表、驻留区间表与 MDADT 检查表能回溯到原始状态和模式序列。
  • expected-pass 对所有子区间通过，full-only-trap 能证明只检查全区间是不充分的。
  • Lab 11 不包含内置轨迹，删除 analysis-input.mat 后会明确要求先运行 Lab 10。
  • 能口头说明数值检查、理论条件和正式证明的区别。