基于典型变量分析的故障识别贡献图研究

1. 引言

1.1 故障识别背景与挑战

在现代工业过程中，故障被定义为过程操作期间发生的任何异常事件。随着制造设施日益集成化和大规模化，故障以非直观方式动态传播的可能性增加，可能对设备、生命和环境造成重大危害。大多数制造设施缺乏精确的第一性原理动态模型，因此工业中的过程监控系统通常基于历史数据库中收集的测量数据构建。当存在大量强相关的过程变量时，故障识别任务可能相当具有挑战性。

1.2 数据驱动方法在过程监控中的应用

统计过程监控（SPM）应用多变量数据驱动方法对过程数据进行故障检测和诊断，在过去二十年中在学术研究和工业实践中都很流行。PCA、PLS和其他改进方法等数据驱动方法用于描述在正常过程条件下收集的数据的特征。这些方法是降维技术，将高维过程数据投影到低得多的维空间中。

1.3 现有贡献图方法的局限性

贡献图是确定哪些变量与统计量不再处于正常操作条件（NOC）最强相关的最流行技术。过程变量的贡献越高，表明特定变量中与故障相关的偏差越大。然而，传统的贡献图只检查一个观测点（时间点）的贡献，对于时间序列数据，需要多个贡献图来说明多个观测点。

1.4 本文研究目标与创新点

本文旨在将贡献图应用于状态空间（CVA模型中保留的状态）和残差空间（模型中剩余的状态），并使用二维彩色图展示贡献。主要创新包括：

• 开发基于典型状态空间和残差空间的两种贡献图
• 将故障变量分类为SSFVs和RSFVs
• 采用二维热力图进行可视化，提供比传统一维图更丰富、更可靠的信息

2. 典型变量分析理论基础

2.1 CVA基本概念与数学原理

2.1.1 典型变量定义与优化目标

CVA是一种多变量降维技术，旨在从两组变量（如输入和输出）中提取具有最大相关性的线性组合。假设有两组变量：输入向量$\bm{x} \in \mathbb{R}^m$和输出向量$\bm{y} \in \mathbb{R}^n$，它们的协方差矩阵分别为$\bm{\Sigma}_{xx}$、$\bm{\Sigma}_{yy}$，互协方差矩阵为$\bm{\Sigma}_{xy}$。

CVA的目标是找到投影矩阵$\bm{J}$和$\bm{L}$，使得投影后的变量：

满足：

• $\bm{c}$和$\bm{d}$的协方差矩阵为单位矩阵（即变量间独立且标准化）
• $\bm{c}$和$\bm{d}$之间的相关性最大化

2.1.2 奇异值分解求解方法

通过奇异值分解（SVD）：

得到：

• 典型相关系数矩阵$\bm{D} = \bm{\Sigma}$
• 投影矩阵：$$\bm{J} = \bm{U}^\mathsf{T} \bm{\Sigma}_{xx}^{-1/2}, \quad \bm{L} = \bm{V}^\mathsf{T} \bm{\Sigma}_{yy}^{-1/2}$$

2.2 CVA状态空间模型构建

2.2.1 过去向量与未来向量定义

给定时间序列输出数据$\bm{y}(t) \in R^{m_y}$和输入数据$\bm{u}(t) \in R^{m_u}$，线性状态空间模型为：

其中$\bm{x}(t) \in R^d$是$d$维状态向量，$\bm{v}(t)$和$\bm{w}(t)$是独立白噪声过程。

过去向量定义为：

未来向量定义为：

2.2.2 状态向量估计方法

状态向量可以通过CVA状态估计为：

其中$\bm{J}_d = \bm{U}_d^\mathsf{T} \widehat{\bm{\Sigma}}_{pp}^{-1/2}$，$\bm{U}_d$包含公式(3)中$\bm{U}$的前$d$列。

2.2.3 滞后阶数与系统阶次选择

由于实际中只能获得有限数量的数据，向量$\bm{p}(t)$和$\bm{f}(t)$通常被截断。系统的真实状态阶次$d$以及最佳的滞后阶数$l$和$h$通常是未知的。解决方案是使用不同滞后阶数拟合多个ARX模型，并采用Akaike信息准则（AIC）等模型选择标准来评估不同模型的拟合优度与复杂度，选择使得AIC值最小的滞后阶数$l$和$h$。

2.3 CVA过程监控统计量

2.3.1 状态空间统计量(T_d²)

该统计量衡量状态空间内的变异，若超出控制限，表示系统状态异常。

2.3.2 残差空间统计量(T_e²)

其中：

3. 基于CVA的贡献图方法

3.1 贡献图基本原理

贡献图是一种广泛使用的故障识别工具。其基本逻辑是，一个变量对故障统计量的"贡献"越大，就越有可能是导致该统计量超限的故障变量。在CVA框架下，我们利用载荷矩阵$\bm{J}_d$和$\bm{J}_e$来"回溯"到原始变量。

3.2 CVA贡献度计算推导

3.2.1 状态空间贡献度(c^d)

其中$c^{d}_{p_{k}}(t) = \bm{x}_d^{\mathrm{T}}(t)(p_{k}(t)\bar{\bm{J}}_{d,k}^{\mathrm{T}})^{\mathrm{T}}$是数据向量$\bm{p}(t)$中第$k$个元素$p_k(t)$对$T_d^2$统计量的贡献。

3.2.2 残差空间贡献度(c^e)

过程变量$u_m$的贡献计算方式与$y_m$相同。

3.2.3 动态时序特性的考虑

由于$\bm{p}(t)$包含多个时间滞后的数据，对于一个特定的过程变量$y_m$，它的贡献是其所有滞后观测值$y_m(t-1), y_m(t-2), ..., y_m(t-l)$的贡献之和。这是CVA贡献图与PCA贡献图的一个重要区别，它明确地考虑了变量的动态时序特性。

3.3 故障变量分类方法

3.3.1 状态空间故障变量(SSFV)

由$c^d$识别。这类变量与模型中已有状态的显著变化有关。意味着过程的动态特性未发生根本改变，但这些状态偏离了正常工况。模型本身对这些变量仍然是有效的。

3.3.2 残差空间故障变量(RSFV)

由$c^e$识别。这类变量表明过程中出现了新的动态特性或状态，这些是原有CVA模型在正常工况下未曾捕捉到的。模型本身已无法充分描述这些变量背后的新动态。

3.3.3 组合贡献度计算方法

由于CVA基于贡献可能因$\sum_{pp}^{-1/2}$的求逆而过于敏感，可以通过组合两种贡献来减少这种敏感性，例如取平均值：

这种组合贡献度在实践中被证明更加鲁棒和可靠。

4. 二维贡献图可视化技术

4.1 传统一维贡献图的局限性

传统的贡献图是1D的条形图，只显示在某个时间点上各变量的贡献值。在动态系统中，由于故障传播效应，不同时间点的1D贡献图可能显示不同的最主要故障变量，导致诊断不可靠。此外，1D图无法提供故障传播路径和持续时间的信息。

4.2 二维贡献图构建方法

4.2.1 数据矩阵组织方式

将所有时间点的"贡献快照"按时间顺序堆叠起来，形成一个时间 × 变量的矩阵。矩阵的每一列对应一个时间点的1D贡献图，每一行对应一个变量在所有时间点的贡献序列。

4.2.2 热力图颜色映射

以变量编号为纵轴，时间为横轴，贡献度以颜色深浅表示。通常使用连续的颜色梯度，如从浅色（低贡献）到深色（高贡献），使高贡献变量在图中呈现为明亮的色带。

4.3 二维贡献图优势分析

4.3.1 时间维度全局视角

2D贡献图提供完整时间序列的全局视角，避免了"瞬时快照"可能带来的误导。操作人员可以一目了然地看到整个故障期间所有变量的贡献变化。

4.3.2 故障传播路径可视化

通过观察色带在时间和变量维度上的变化，可以清晰地可视化故障的传播路径，区分是短暂的扰动还是持续的故障。

4.3.3 自动诊断时间窗口识别

2D贡献图无需先验知识即可自动捕捉最佳诊断时间窗口，能够识别出那些在整个故障期间持续高贡献的关键变量。

5. 田纳西-伊斯曼过程案例研究

5.1 实验设置与数据预处理

5.1.1 过程变量说明

田纳西-伊斯曼过程（TEP）是一个著名的化工过程基准测试问题，包含41个测量变量和12个操纵变量。过程包括三个主要单元（反应器、分离器和汽提塔），从四种反应物（A、C、D、E）生产两种产品（G和H）。

5.1.2 故障场景设计

TEP提供了21种预编程故障，用于评估过程控制和监控方法。每个测试数据集包含960个观测值，每3分钟采样一次。每个数据集开始时没有故障，故障在第160个样本后发生。

5.1.3 数据标准化处理

为避免特定变量不适当地主导降维过程，所有数据都采用标准化的"z-scoring"方法进行处理，即每个变量减去其均值并除以其标准差。

5.2 故障1详细分析

5.2.1 故障机理说明

故障1是流4中A/C进料比的阶跃变化。具体来说，流4中C的组成从51 mol%增加到54 mol%，而A的组成从48.5 mol%减少到45.5 mol%。这些变化导致流5中A减少，进而通过控制回路调节增加了流1中的A进料。

5.2.2 贡献图结果对比

基于状态空间的CVA贡献图（$c^d$）和基于残差空间的CVA贡献图（$c^e$）都显示大量"噪声"，这是由于依赖单一统计量的CVA基于贡献的过度敏感性。组合贡献图（$c = (c^d + c^e)/2$）显示出更加清晰和可靠的结果。

5.2.3 与传统方法比较

CVA基于贡献正确识别了变量$x_{34}$作为故障变量，而Liu的优化PCA基于方法[20]未能识别该变量。变量$x_{34}$（流9中的组分F）最终稳定在高于其NOC值的稳态值。

5.3 故障4详细分析

5.3.1 故障特征分析

故障4是反应器冷却水入口温度的阶跃变化，导致反应器冷却水流速（$x_{51}$）的阶跃变化。当故障在第160个样本发生时，反应器温度（$x_9$）突然升高，随后通过控制回路得到补偿。

5.3.2 状态空间与残差空间对比

对于故障4，基于状态空间的贡献显示多个变量出现大的瞬态偏差，而基于残差空间的贡献仅持续显示$x_{51}$的贡献。控制系统能够补偿故障对过程操作的影响，因此在51个变量中唯一持续的变化是反应器冷却水流速$x_{51}$。

5.3.3 RSFV识别优势

故障4的故障变量$x_{51}$是一个RSFV，如残差空间贡献图所示。在状态空间中观察到来自正常操作条件的最小干扰，而故障4导致残差空间中一个变量的大偏差。

5.4 全部21种故障综合分析

5.4.1 故障分类统计

• 故障1、2、5、6、8、10、12、13、14、16、17、18和20改变了状态空间和残差空间的属性，导致贡献$c^d$和$c^e$都可用于识别故障变量。
• 故障3、9、15和21的属性未被状态空间或残差空间捕获。这些故障使用任何仅基于历史数据的方法都很难检测。
• 故障4、7、11和19可以被归类为第三类故障，其中只有残差空间受到显著影响。

5.4.2 SSFV与RSFV分布规律

除故障10、12、16和20外，残差空间中出现的故障变量数量多于状态空间中的数量。换句话说，异常事件往往创建RSFV而不是SSFV。故障变量更可能出现在残差空间中表明测量噪声的特性已经改变和/或过程中已经创建了新状态，并且用于模型的现有状态$x_d(t)$不再足以描述过程动态或输入-输出关系。

5.4.3 方法有效性验证

对于某些故障，一些变量，如故障1中的$x_1$和$x_{18}$，表现为SSFV或RSFV，但不同时表现为两者。其他一些变量，如故障2中的$x_{47}$，既是SSFV又是RSFV，这表明故障变量在CVA模型的状态和残差空间中都有显著变化。

6. 结论与展望

6.1 主要研究成果总结

本文提出了基于CVA的贡献图用于识别与故障密切相关的变量，其中包括基于状态和残差空间变化的贡献。故障变量被分类为SSFV和RSFV，这为了解每个故障的特性提供了一些见解。所提出的方法已在TEP的所有21个故障中得到验证，仿真结果表明，基于CVA的贡献图识别的故障变量可以影响状态空间、残差空间或两者的统计量；并且观察到异常事件更常与残差空间中的故障变量相关，而不是状态空间中的故障变量。

6.2 方法优势与局限性

基于CVA的贡献图能够有效处理序列相关性，考虑过程动态特性，并通过二维可视化提供更全面的故障信息。然而，CVA基于贡献可能因协方差矩阵求逆而过于敏感，需要组合贡献来缓解这一问题。

6.3 未来研究方向

6.3.1 权重优化策略

案例研究中使用的组合贡献对状态和残差空间使用了相等的权重（$c = (c^d + c^e)/2$）。未来的一个考虑是，如果对状态和残差空间使用不同的权重，CVA基于故障识别是否会得到改进。

6.3.2 先验知识融合

最优权重的定义可能需要假设关于故障的先验知识，例如它们的概率或可能受其影响的变量数量的上限。

6.3.3 工业应用拓展

将CVA基于贡献图和二维可视化技术应用于更广泛的工业过程，验证其在复杂工业环境中的有效性和鲁棒性。

附录

A. 田纳西-伊斯曼过程变量表

B. 全部故障识别结果汇总表

详细列出21种故障的SSFV和RSFV识别结果，展示CVA基于贡献图在各种故障场景下的性能。

C. 二维贡献图实现伪代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

% 预处理测试数据（自动缩放）
X = pretreat(X);
% 计算贡献
[m,n] = size(X);
CONT = zeros(m,n);
for k = 1:m % 对于每个新观测值k
    t = P * X(k,:);
    cont = zeros(a,n);
    for i = 1:a % 对于每个主成分
        for j = 1:n % 对于每个过程变量
            cont(i,j) = t(i) * X(k,j) * P(i,j) / σ_i^2;
            cont(i,j) = cont(i,j) * (cont(i,j) > 0);
        end
    end
    CONT(k,:) = sum(cont, 1);
end
% 绘制二维贡献图
imagesc(CONT);
xlabel('新观测值（小时）');
ylabel('过程变量');

Canonical variate analysis-based contributions for fault identification

文章主要内容

将贡献图应用于状态空间（CVA模型中保留的状态）和残差空间（模型中剩余的状态），并且使用二维彩色图展示贡献。

本文两种贡献图

• 基于典型状态空间
• 基于残差空间
  分别分类为
• SSFVs 状态空间故障变量
• RSFVs 残差空间故障变量
  仿真得到的基本结论是异常事件更常与残差空间中的故障变量相关

背景

1. 现有的工业过程监控系统通常基于历史测量数据构建：
   大多数工业系统没有精确的第一性原理动态模型
2. 存在大量强相关的变量时，故障识别任务可能具有挑战性
3. 降维技术可以提高识别故障能力：PCA、PLS、CVA
4. PCA、PLS适合i.i.d.的变量（不考虑时间延迟的情况下估计协方差，从而解释观测值的变化）
5. 序列相关的变量，零时间延迟协方差矩阵不能完全代表整个变化
6. 对于序列相关的变量，构建带时间延迟协方差矩阵，再运用PCA和PLS方法，提取动态模型。
7. 利用与零特征值相对应的协方差矩阵特征向量，建立多变量自回归外输入（ARX）模型。缺点ARX对于动态过程不灵活性。
8. CVA考虑了序列相关,状态变量在零时间延迟下是统计独立的
9. 贡献图用于确定哪些变量与统计量与离开正常状态最强相关。过程变量贡献越高特定变量与故障变化相关越强。
10. 改进：基于PCA（并不总能正确识别）、基于T^2和Q统计量（比单一统计量有效）、引入置信限（PCA残差的模糊效应可能误导故障变量的判定）、通过缺失数据法最大化组合指标降低（对非故障变量无模糊效应的贡献图）。
11. 之前的研究比较数据驱动方法和模型驱动方法，简单故障识别通过贡献图可以轻松实现，复杂故障判定需要额外信息。
12. 基于重构的方法，从异常时间子空间确定故障变量，使用识别指数（重构SPE/故障SPE）识别故障
13. 即使使用滞后数据基础PCA和PLS也无法产生最准确的动态模型。
14. 贡献图和状态空间结合使用。同时考虑CVA。

contribution map:贡献图

贡献图是一种广泛使用的故障识别工具。其基本逻辑是，一个变量对故障统计量的“贡献”越大，就越有可能是导致该统计量超限的故障变量。

CVA:典型成分分析

• 降维技术
• 系统动态过程
• 序列相关性
  但是在这个文章之前，
  没有关于它用于工业过程的故障识别的有效性的广泛研究
  识别与故障最相关的变量（根原因和故障溯源）

CVA 的统计方法

CVA 是一种多变量降维技术，旨在从两组变量（如输入和输出）中提取具有最大相关性的线性组合。

• 假设有两组变量：
  • 输入向量 $\bm{x} \in \mathbb{R}^m$
  • 输出向量 $\bm{y} \in \mathbb{R}^n$
• 它们的协方差矩阵分别为 $\bm{\Sigma}_{xx}$、$\bm{\Sigma}_{yy}$，互协方差矩阵为 $\bm{\Sigma}_{xy}$。

CVA 的目标是找到投影矩阵 $\bm{J}$ 和 $\bm{L}$，使得投影后的变量：

满足：

• $\bm{c}$ 和 $\bm{d}$ 的协方差矩阵为单位矩阵（即变量间独立且标准化）
• $\bm{c}$ 和 $\bm{d}$ 之间的相关性最大化

通过奇异值分解（SVD）：

得到：

• 典型相关系数矩阵 $\bm{D} = \bm{\Sigma}$
• 投影矩阵：$$\bm{J} = \bm{U}^\mathsf{T} \bm{\Sigma}_{xx}^{-1/2}, \quad \bm{L} = \bm{V}^\mathsf{T} \bm{\Sigma}_{yy}^{-1/2}$$

CVA 状态向量

CVA 也可用于构建动态系统的状态空间模型，适用于具有时序相关性的过程数据。

• 考虑状态空间模型：$$\begin{aligned} \bm{x}(t+1) &= \bm{A} \bm{x}(t) + \bm{B} \bm{u}(t) + \bm{v}(t) \\ \bm{y}(t) &= \bm{C} \bm{x}(t) + \bm{D} \bm{u}(t) + \bm{E} \bm{v}(t) + \bm{w}(t) \end{aligned}$$
• 定义：
  • 过去向量 $\bm{p}(t)$：包含过去的输入和输出数据
  $\boldsymbol{p}(t) = \left[ \boldsymbol{y}^{\mathrm{T}}(t-1), \boldsymbol{y}^{\mathrm{T}}(t-2), \ldots, \boldsymbol{u}^{\mathrm{T}}(t-1), \boldsymbol{u}^{\mathrm{T}}(t-2), \ldots \right]^{\mathrm{T}}$
  • 未来向量 $\bm{f}(t)$：包含当前和未来的输出数据
  $\boldsymbol{f}(t) = \left[ \boldsymbol{y}^{\mathrm{T}}(t), \boldsymbol{y}^{\mathrm{T}}(t+1), \ldots \right]^{\mathrm{T}}$

通过将 CVA 应用于 $\bm{p}(t)$ 和 $\bm{f}(t)$，可以估计状态向量：

其中 $\bm{U}_d$ 是 SVD 中前 $d$ 列组成的矩阵，$d$ 是系统阶次。

数据向量的实际截断

由于在实际应用中，我们只能获得有限长度的数据，因此必须将向量 $\bm{p}(t)$ 和 $\bm{f}(t)$ 进行截断：

其中：

• $l$：过去向量的**滞后阶数**，决定了历史窗口的大小。
• $h$：未来向量的**预测步长**，决定了未来窗口的大小。

系统阶次与滞后阶数的选择

系统的真实状态阶次 $d$ 以及最佳的滞后阶数 $l$ 和 $h$ 通常是未知的。论文指出，一个有效的解决方案是：

1. 使用不同滞后阶数拟合多个带外生变量的自回归模型（ARX模型）。
2. 采用赤池信息量准则（Akaike Information Criterion, AIC） 等模型选择标准来评估不同模型的拟合优度与复杂度。
3. 选择使得AIC值最小的滞后阶数 $l$ 和 $h$，从而在模型精度和复杂性之间取得最佳平衡。

通过这种截断和优化选择，状态向量得以在实际中计算。

过程监控统计量

基于 CVA 模型，可以构造两个关键统计量用于过程监控：

1. 状态空间统计量 $T_d^2(t)$

• 衡量状态空间内的变异
• 若超出控制限，表示系统状态异常

2. 残差空间统计量 $T_e^2(t)$

其中：

• $\bm{U}_e$ 是 SVD 中剩余的列
• 衡量残差空间内的变异
• 若超出控制限，表示噪声特性变化或出现新状态

基于CVA的贡献度

1. 贡献度的计算推导

利用CVA的载荷矩阵 $\bm{J}_d$ 和 $\bm{J}_e$ 来“回溯”到原始变量。

• 状态空间贡献度 $c^d$:
  我们从状态空间统计量 $T_d^2$ 出发：

  $$T_d^2(t) = \bm{x}_d^\mathsf{T}(t)\bm{x}_d(t) = \bm{x}_d^\mathsf{T}(t) (\bm{J}_d \bm{p}(t))$$

  这个公式可以解读为：统计量 $T_d^2$ 是由状态向量 $\bm{x}_d(t)$ 与通过载荷矩阵 $\bm{J}_d$ 投影过来的过去数据 $\bm{p}(t)$ 共同作用的结果。

  因此，我们可以将 $T_d^2$ 分解为向量 $\bm{p}(t)$ 中每一个元素 $p_k(t)$ 的贡献之和：

  $$c^{d}(t) = \sum_{k=1}^{(m_{y}+m_{u})l} \underbrace{\bm{x}_d^{\mathsf{T}}(t)(p_{k}(t) \bar{\bm{J}}_{d,k}^{\mathsf{T}})^{\mathsf{T}}}_{c^{d}_{p_k}}$$

  其中，$c^{d}_{p_k}$ 就是过去数据向量中第 $k$ 个元素 $p_k(t)$ 对 $T_d^2$ 统计量的贡献。

• 关键点：处理动态与滞后:
  由于 $\bm{p}(t)$ 包含了多个时间滞后的数据（见公式(9)），对于一个特定的过程变量 $y_m$，它的贡献是其所有滞后观测值 $y_m(t-1), y_m(t-2), ..., y_m(t-l)$ 的贡献之和：

  $$c^{d}_{y_{m}}(t) = \sum_{j=1}^{l} \bm{x}_d^{\mathsf{T}}(t)(y_{m}(t-j) \bar{\bm{J}}_{d,m_{j}})^{\mathsf{T}}$$

  这是CVA贡献图与PCA贡献图的一个重要区别，它明确地考虑了变量的动态时序特性。

• 残差空间贡献度 $c^e$:
  残差空间贡献度的计算方式与状态空间完全类似，只是将状态向量 $\bm{x}_d(t)$ 和载荷矩阵 $\bm{J}_d$ 替换为残差向量 $\bm{x}_e(t)$ 和载荷矩阵 $\bm{J}_e$：

  $$c^{e}_{y_{m}}(t) = \sum_{j=1}^{l} \bm{x}_e^{\mathsf{T}}(t)(y_{m}(t-j) \bar{\bm{J}}_{e,m_{j}})^{\mathsf{T}}$$

2. 贡献度的融合与故障变量分类

• 问题：单一贡献度的敏感性
  论文指出，由于在计算状态向量时需要求逆矩阵 $\widehat{\Sigma}_{pp}^{-1/2}$，这可能导致基于单一统计量（$T_d^2$ 或 $T_e^2$）的贡献图过于敏感，产生“噪声”并可能误判。

• 解决方案：组合贡献度
  为了克服这个问题，作者借鉴了PCA/PLS中的思想，提出将两种贡献度结合起来使用，例如取平均值：

  $$c = (c^d + c^e)/2$$

  这种组合贡献度在实践中被证明更加鲁棒和可靠。

• 创新的故障分类：SSFV 与 RSFV
  这是本文的一个核心贡献。通过分别考察 $c^d$ 和 $c^e$，作者将识别出的故障变量分为两类：

  • 状态空间故障变量（SSFV）: 由 $c^d$ 识别。这类变量与模型中已有状态的显著变化有关。意味着过程的动态特性未发生根本改变，但这些状态偏离了正常工况。模型本身对这些变量仍然是有效的。
  • 残差空间故障变量（RSFV）: 由 $c^e$ 识别。这类变量表明过程中出现了新的动态特性或状态，这些是原有CVA模型在正常工况下未曾捕捉到的。模型本身已无法充分描述这些变量背后的新动态。

  这种分类为操作人员提供了更深层次的故障洞察，有助于理解故障的本质。

3. 2D贡献图

• 传统的贡献图是1D的条形图，只显示在某个时间点上各变量的贡献值。
• 本文采用二维热力图，以变量编号为纵轴，时间为横轴，贡献度以颜色深浅表示。
• 优势：这种可视化方法可以清晰地展示故障的传播路径和持续时间，能够区分出是短暂的扰动还是持续的故障，从而提供比1D图更丰富、更可靠的信息。

Two-Dimensional Contribution Map for Fault Identification:二维贡献图

将所有时间点的“贡献快照”按时间顺序堆叠起来，形成一个 时间 × 变量 的矩阵，并用热力图进行可视化。

二维贡献图的核心优点

1. 提供时间维度全局视角，避免“瞬时快照”的误导
2. 无需先验知识，自动捕捉最佳诊断时间窗口
3. 直观揭示故障的传播路径与动态特性

对比

参考文献

1. Jiang, B., Huang, D., Zhu, X., Yang, F., & Braatz, R. D. (2015). Canonical variate analysis-based contributions for fault identification. Journal of Process Control, 26, 17-25.
2. Zhu, X., & Braatz, R. D. (2014). Two-dimensional contribution map for fault identification. IEEE Control Systems Magazine, 34(5), 72-77.
3. Chiang, L. H., Russell, E. L., & Braatz, R. D. (2001). Fault Detection and Diagnosis in Industrial Systems. Springer Verlag.
4. Qin, S. J. (2012). Survey on data-driven industrial process monitoring and diagnosis. Annual Reviews in Control, 36(2), 220-234.

李雅普诺夫稳定性理论完整教程

一、引言：稳定性的重要性

控制系统中，由于外部扰动和内部参数变化的存在，系统的状态会偏离预设轨迹。系统的稳定性是其正常工作的基本前提，即系统在受到扰动后能够自行恢复到平衡状态或保持在平衡点附近。

李雅普诺夫稳定性理论由俄国数学家亚历山大·李雅普诺夫于1892年提出，是现代控制理论的基石。该理论分为两大类：

• 间接法（第一方法）：通过线性化模型分析稳定性；
• 直接法（第二方法）：通过构造能量函数直接判断稳定性。

二、李雅普诺夫稳定性的定义

考虑非线性系统：

设 $x(t)$ 是系统在初始条件 $x(t_0) = x_0$ 下的解。若对任意 $\varepsilon > 0$，存在 $\delta(\varepsilon) > 0$，使得当初始状态 $\tilde{x}_0$ 满足 $\| \tilde{x}_0 - x_0 \| < \delta$ 时，系统的解 $\tilde{x}(t)$ 对所有 $t > t_0$ 都满足：

则称解 $x(t)$ 是稳定的。否则，称为不稳定的。

• 大范围稳定：若 $\delta$ 可以任意大，即从任意初始状态出发，系统都稳定；
• 渐近稳定：若系统稳定，且存在 $\delta$ 使得 $\tilde{x}(t) \to x(t)$；
• 全局渐近稳定：若系统在整个状态空间内都是渐近稳定的。

三、李雅普诺夫间接法（第一方法）

基本思想

通过系统在平衡点处的线性化模型判断其局部稳定性。

步骤

1. 求平衡点 $x_e$：解 $f(x_e) = 0$；
2. 计算雅可比矩阵：$$A = \left. \frac{\partial f}{\partial x^T} \right|_{x = x_e}$$
3. 分析线性化系统 $\dot{y} = A y$（其中 $y = x - x_e$）；
4. 根据 $A$ 的特征值判断稳定性：
   • 所有特征值实部为负 ⇒ 渐近稳定；
   • 存在特征值实部为正 ⇒ 不稳定；
   • 存在实部为零的特征值 ⇒ 无法判断（需进一步分析非线性项）。

例题分析

例题1： 判断系统在原点处的稳定性

解：

1. 验证原点是平衡点
2. 计算雅可比矩阵：

3. 特征值为 $-\frac{1}{2} \pm j\frac{\sqrt{3}}{2}$，实部均为负
4. 结论： 原点是渐近稳定平衡点

优缺点

• 计算简便，适合快速分析局部稳定性；
• 仅适用于局部稳定性分析，对临界情况无效。

四、李雅普诺夫直接法（第二方法）

基本思想

构造一个标量函数 $V(x)$（称为李雅普诺夫函数），通过分析 $V(x)$ 及其沿系统轨迹的导数 $\dot{V}(x)$ 来判断稳定性。

标量函数的定号性

• 正定：$V(0) = 0$，且 $V(x) > 0 \ \forall x \neq 0$；
• 半正定：$V(x) \geq 0$；
• 负定：$-V(x)$ 正定；
• 半负定：$-V(x)$ 半正定；
• 不定：无固定符号。

稳定性判据

1. 稳定：$V(x)$ 正定，$\dot{V}(x)$ 半负定；
2. 渐近稳定：$V(x)$ 正定，$\dot{V}(x)$ 负定；
3. 全局渐近稳定：满足渐近稳定，且 $V(x)$ 径向无界；
4. 不稳定：$V(x)$ 正定，$\dot{V}(x)$ 正定。

五、李雅普诺夫函数的构造方法与实例

1. 二次型函数法（线性与轻度非线性系统）

例题2： 线性系统分析

解法1： 取 $V(x) = x_1^2 + x_2^2 > 0$

半负定，需进一步验证：当 $\dot{V}(x) = 0$ 时，$x_2 = 0$，代入原方程得 $x_1 = 0$，故除原点外 $\dot{V}(x)$ 不恒为零 ⇒ 渐近稳定

解法2： 取 $V(x) = 1.5x_1^2 + x_2^2 + x_1x_2 = x^T \begin{bmatrix} 1.5 & 0.5 \\ 0.5 & 1 \end{bmatrix} x > 0$

负定 ⇒ 全局渐近稳定

结论： 同一系统可构造不同的李雅普诺夫函数，得到相同结论但证明难度不同。

例题3： 非线性系统分析

构造： 取 $V(x) = 3x_1^2 + 4x_2^2 > 0$

负定，且当 $\|x\| \to \infty$ 时 $V(x) \to \infty$ ⇒ 全局渐近稳定

技巧： 对于含高阶项的非线性系统，选择适当的二次型函数可消去交叉项。

2. 克拉索夫斯基方法

适用系统： $\dot{x} = f(x)$，其中 $f(0) = 0$

方法：

• 构造 $V(x) = f^T(x) P f(x)$，$P > 0$；
• 若 $F^T(x) + F(x)$ 负定（其中 $F(x) = \frac{\partial f}{\partial x^T}$），则原点渐近稳定；
• 若 $\|x\| \to \infty$ 时 $\|f(x)\| \to \infty$，则为全局渐近稳定。

例题4：

构造：

1. 计算雅可比矩阵：

2. 分析 $F^T(x) + F(x)$ 的定号性：

• 顺序主子式：$\Delta_1 = x_2 - 1$，$\Delta_2 = 1 - x_2 - x_1^2$
• 在区域 $x_2 < 1 - x_1^2$ 内，$\Delta_1 < 0$，$\Delta_2 > 0$ ⇒ 负定

结论： 原点在该区域内是渐近稳定的

3. 变量梯度法（Schultz-Gibson）

步骤：

• 假设 $\nabla V = [\nabla_1, \dots, \nabla_n]^T$；
• 通过旋度条件 $\frac{\partial \nabla_i}{\partial x_j} = \frac{\partial \nabla_j}{\partial x_i}$ 确定系数；
• 沿路径积分得 $V(x)$；
• 验证 $V(x)$ 正定，$\dot{V}(x)$ 负定。

例题5：

构造步骤：

1. 设梯度向量：

2. 由旋度条件得 $a_{12} = a_{21}$

3. 取 $a_{12} = a_{21} = 0$，$a_{11} = a_{22} = 1$

4. 计算：

在 $1 - 2x_1x_2 > 0$ 区域内负定

5. 积分得：

结论： 原点在局部区域内是渐近稳定的

4. 偶函数法

适用系统： 具有对称结构的系统

方法：

• 构造 $V(x) = \sum a_i(x_i)$，其中 $a_i(x_i)$ 是偶正定函数；
• 利用奇偶性简化 $\dot{V}(x)$ 的符号判断。

例题6：

构造：

1. 设 $V(x) = a(x_1) + b(x_2)$，其中 $a(x_1), b(x_2)$ 为偶函数

2. 取：

3. 积分得：

4. 最终：

结论： 原点渐近稳定（非全局）

六、线性系统的稳定性分析

线性定常系统：$\dot{x} = Ax$

• 渐近稳定 ⇔ $A$ 的所有特征值实部为负；
• 稳定 ⇔ 所有特征值实部非正，且虚轴上的特征值对应一阶约当块；
• 不稳定 ⇔ 存在实部为正的特征值。

Routh-Hurwitz判据

• 通过特征方程的系数构造Routh表；
• 第一列全正 ⇒ 系统稳定；
• 变号次数 = 右半平面根的数量。

例题7： 特征方程 $2s^6 + 5s^5 + 3s^4 + 4s^3 + 6s^2 + 14s + 7 = 0$

构造Routh表，第一列出现负号，说明有右半平面根，系统不稳定。

李雅普诺夫方程判据

• 系统渐近稳定 ⇔ 对任意 $Q > 0$，存在 $P > 0$ 满足：$$A^TP + PA = -Q$$
• 可通过求解李雅普诺夫方程判断稳定性。

例题8： 判定系统 $\dot{x} = \begin{bmatrix} 0 & 1 \\ -1 & -1 \end{bmatrix} x$ 的稳定性

取 $Q = I$，解李雅普诺夫方程得：

故系统渐近稳定。

七、参数稳定性与稳定域

单参数稳定域

例题9： 单位负反馈系统，开环传递函数：

特征方程：

根据Routh判据：

双参数稳定域

例题10：

特征方程：

稳定条件：

八、吸引域分析

对于非线性系统，确定吸引域是重要而困难的任务。

保守吸引域估计方法：
在 $\dot{V}(x) = 0$ 的边界上求 $V(x)$ 的最小值 $V_{\min}$，则

就是一个保守的吸引域。

九、方法总结与对比

叶下落着
摇着落下来
阳光下翻着金色
红的黄的绿的叶
来自风的方向，
时时刻刻飘来
似乎这已经是永远
但雨也会停的
泪也会干的

伤了膝盖，还是不太稳重。
站定之后，感觉缓过来没那么疼之后，
一瞬间感觉好累啊。

用多阶段决策的术语描述最短路问题

已知

• 状态集 ( S_k, k=1,2,\cdots,6 )
• 允许决策集 ( U_k(s_k), \forall s_k \in S_k, k=1,2,\cdots,5 )
• 状态转移方程 ( s_{k+1} = T_k(s_k, u_k), \forall s_k \in S_k, u_k(s_k) \in U_k(s_k), k=1,2,\cdots,5 )
• 阶段指标函数 ( d_k(s_k, u_k), \forall s_k \in S_k, u_k(s_k) \in U_k(s_k), k=1,2,\cdots,5 )

问题

求 ( p_{1,5} \in P_{1,5} ) 使下述过程指标函数达到最小
[ V_{1,5}(s_1, p_{1,5}) = \sum_{k=1}^{5} d_k(s_k, u_k) ]

一

泪似乎不能轻易流下
不然会失去分量

心确然顾不了这么多
匆然就随着起落

恍惚、眩晕
我在这里
但是除了壳
已经空无一物

二

那边叶子
轻薄地扇了扇
飞走了

三

一直以为
我从来都没有变。
变的太多了。

模拟-离心力-检测法（径向）

模拟-科里奥利力-MEMS陀螺仪测转速（横向）

MEMS:Micro-Electro-Mechanical System。微机电结构。简单说就是，微型化的机械电子元件。
本质上说就是利用了一个惯性传感器。
关键：位移和转速成正比。
K的标定。

数字-位置编码-光电码盘转速检测法

光电耦合器

• 透射型
• 反射型
  光电码盘
• 绝对码
  • 二进制码
  • 循环码（格雷码）
• 增量码
  • 单相 只能测转速（计数脉冲个数）
  • 三相 AB相位差可以知道转动方向，C相脉冲位置设为零点，可以读出转角。
    霍尔计数方式（霍尔编码器）
    将编码的器件换成永磁体，耦合器换成霍尔传感器

空间滤波器式检测法

检测旋转物体表面不规则光反射分布。
从频率计算。

测量特点对比

精度对比

标准

成本

出题

1.光电码盘计算转速
2.MEMS计算转速
3.

动机

多面体模型

• 多面体：由一组线性约束定义的可行域。
• 极点（“绿点”）：由若干起作用约束（不等式成为等式）确定的唯一解，即可行域的顶点。
• 最优解必为某一极点（讨论三种情况：内点；边界且目标梯度与边界正交；边界且不正交）。

标准型

• 标准形式：最小化或最大化线性目标，约束为等式并附非负性限制。
• 松弛变量与引入基变量用于转化不等式为等式。

单纯形法

• 迭代在极点间移动以改善目标值，直到达到最优极点。
• 特殊情况：退化（基变量为0）、多重最优解、无界解、无可行解。
• 初始可行解的获取：人工变量法、两阶段法等。
• 对偶单纯形法：当基解满足最优性但不满足可行性时使用。

基于逆矩阵的单纯形迭代

• 使用基矩阵逆更新基解与灵敏系数，减少每步计算量。
• 通过列替换和更新公式维持效率。

对偶问题与互补关系

• 每个线性规划都有对应的对偶问题（标准型与对偶型）。
• 互补松弛（互补性条件）：原问题与对偶问题的基解在最优时满足互补松弛条件。
• 弱对偶性：任一原问题可行解的目标值都是对偶问题最优目标的下界/上界（取决于最大/最小化）。
• 强对偶性：若一方存在最优解，则另一方也存在且两者最优目标值相等。

影子价格与灵敏度分析

• 影子价格（对偶变量）表示约束右端项微小变动对目标值的边际影响。
• 灵敏度分析包括系数变化、右端项变化和基结构稳定性分析。
• 参数化线性规划：分析目标或约束参数在区间内变化对最优解的影响。

整数与0-1线性规划

• 割平面法：通过加入有效割约束逐步切除非整点解。
• 分枝定界法：枚举式搜索结合界值剪枝以求整数最优解。
• 0-1问题常见模型：选择性约束（p-选q等形式）、背包、指派等。

动态规划

• 分阶段模型：将问题分解为若干阶段和状态转移。
• 最优性原理：子问题的最优解构成原问题最优解的部分。
• 主要算法：值迭代与策略迭代（用于确定最优策略与价值函数）。
• 适用场景：序列决策、资源分配、最短路径、背包的分阶段解法等。

动机

线性二分类问题，找到最优分类面

最优分类面

• 优点：容忍局部扰动
• 定义：离两类样本最小距离最大
  由定义容易得到，对称；而且是一个min-max问题

拓展

• 拓展至线性不可分问题
• 离群值的处理

最优超平面只由支持向量决定

超平面的求解

问题的重写-有约束的优化问题

原始问题与等价问题。可以归结到凸优化问题。

一般凸优化问题

最小化目标函数；约束不等式组；约束等式组。
约束构成了可行域。
归结到一般的凸优化问题。

使用拉格朗日乘数法求解

对于一个有约束的优化问题，使用拉格朗日乘数法转化到无约束优化问题。
对乘子求拉格朗日函数最大化，使得其满足约束。
再求最小值则转化为原问题。
此时是一个min-max优化问题，不易求解，它的对偶问题max-min

对偶问题

三个定理。
定理1：原始问题是对偶问题的上界。
定理2：在一定条件下，两个问题的最优解相等。
定理3：满足定理2的一定条件，$W^*$是原问题最优解，$\alpha^*$、$\beta^*$是对偶问题最优解，只要它们一起满足KKT条件。

对偶形式计算

非线性问题-核函数

如果原始空间样本特征是有限维的，那么一定存在一个高维特征空间使得样本变得可分。

• 寻找低维向高维映射函数（基函数）
• 利用核函数简化计算，核函数是一类广义内积，是这种映射的隐式表征

常见核函数

• 线性
• 多项式
• RBF
• sigmoid

非线性问题-软间隔