工业 AI：使用 TensorFlow Probability 预测已知的未知问题

文 / Venkatesh Rajagopalan，数据科学与分析总监和 Arun Subramaniyan，BHGE Digital 数据科学与分析部副总裁

昨天，我们发布的此系列的第一篇文章中，我们介绍了我们的分析理念：将领域知识、概率方法、传统机器学习 (ML) 和深度学习技术结合起来，解决工业领域中某些最难以解决的问题。我们还用构建裂纹扩展混合模型的示例说明了此方法。与 Google 的 TensorFlow Probability (TFP) 团队和 Cloud ML 团队的持续协作加快了我们大规模开发和部署这些技术的速度。通过将物理知识与看似无关的技术相结合，并使用现代的可扩展软硬件部署它们，之前难以解决的问题现在都可以解决。

在本文中，我们将展示分析方法的另一个应用，即使用新增可用信息更新现有模型。我们会使用《工程系统预测与健康管理：入门》(Prognostics and Health Management of Engineering Systems: An Introduction) 一书中提供的裂纹扩展示例。

更新模型的需求

构建预测模型的目的在于预测系统未来的状态。这些预测模型可以是基于系统的已知物理知识完全由数据驱动的模型，也可以是数据与物理知识相结合的混合模型。在工业领域，这种模型用于估算组件（例如涡轮叶片）或系统（例如电动潜油泵）的剩余使用寿命。研究人员可以使用实际的现场数据（如适用），或者使用来自模拟或对照实验的数据开发这些模型。

当使用模拟或对照实验数据开发的模型建好时，研究人员需要使用现场数据对其进行验证和优化。但是，模型建好后，再单独使用现场数据重新构建模型通常并不可行，因为特定时间的现场数据数量往往非常有限。此外，研究人员通常希望保留模型 “内核”，仅修改其系数。因此，分析方法需要满足能够持续更新模型以提升模型准确度的要求。

我们会考虑以下情况：正如我们在第一篇文章中的操作，我们根据物理知识和已知材料属性构建了一个裂纹扩展初始模型。在使用该模型预测裂纹尺寸时，模型的初期预测合理准确，但随着时间的推移，我们发现预测和观测结果之间存在显著差异。

因此，利用新的可用裂纹测量数据更新现有模型，有助提升模型的准确度。这种情况在工业领域非常常见；如果材料属性发生改变，或在组件中引入了新的设计（例如燃气轮机的叶片），则不会存在可用的现场数据。本文将探讨如何采用一种基于无迹卡尔曼滤波 (UKF) 的方法更新裂纹扩展初始模型。

关于数据的简要讨论

《工程系统预测与健康管理：入门》中提供的这个数据集包含裂纹测量结果和周期， 如下图所示：

虽然此数据选自教材示例，但它展示的是现场数据中的某些常见特征，例如：

1. 噪声

与周期相关的裂纹尺寸总体趋于上升，但并非完全没有变化。某些实例的裂纹尺寸观测结果小于前一个实例的测量结果。出现这种现象的原因很可能是测量错误，这也是现场数据的共同特征。换言之，此数据非常嘈杂。

2. 样本量小

该数据集仅由 16 个数据点组成，不足以构建任何有意义的数据驱动型模型。因此，要构建有用的模型，充分利用对裂纹扩展物理知识的理解至关重要。

为了使情况进一步复杂化，在实际应用中，我们甚至不可能从同一个组件或系统中获得逐渐退化的数据。可用于建模的数据通常来源于处在不同退化阶段的不同资源。因此，混合建模方法对工业分析至关重要。

为了解释更新模型的方法，我们将数据集分成三个部分：

• 初始数据集 — 此数据集中的初始模型预测合理准确。它由前六个数据点组成 [周期 0 到 500]

• 更新数据集 — 此数据集中的初始模型预测不准确。从其名称可知，此数据集用于更新初始模型。它由接下来的五个数据点组成 [周期 600 到 1000]

• 测试数据集 — 此数据集用于验证更新后模型的预测结果。它由最后五个数据点组成 [周期 1100 到 1500]

使用最小二乘解获得的初始模型参数为 logC = -22.135 和 m =3.765，其中 C 和 m 为材料属性。初始数据集的模型预测如下图所示：

关于测量值的模型预测相当准确。我们发现初始数据集的均方误差 (MSE) 为：MSE(i)=3.93 10^-7。

模型在更新数据集对应时间段内的裂纹尺寸预测如下图所示。

显然，初始模型的预测不再准确。更新数据点上的 MSE 变为 MSE(u)=5.67 10^-6，几乎是初始误差 MSE(i) 的 14 倍！而且，随着周期延长，模型预测与更新数据集中观察到的裂纹尺寸之间存在明显差异。因此，我们需要优化模型，以使其更加准确。

更新方法

可以用于模型更新的技术之一是卡尔曼滤波 (KF) 及其变体，即扩展卡尔曼滤波 (EKF) 和无迹卡尔曼滤波 (UKF)。

采用 KF 时，我们做出以下假设：

• 过程模型为线性模型

• 测量模型为线性模型

• 初始状态为高斯分布

如果采用 EKF 和 UKF，可以放宽前两个假设。在卡尔曼滤波器中，执行的中央运算是通过系统动力学传播高斯随机变量 (GRV)。在 EKF 中，GRV 接近状态分布，然后通过非线性系统的一阶线性化分析传播。这种传播可能会为转化后的 GRV 后验均值和协方差带来很大误差，可能造成次优性能，有时还会出现滤波发散。

UKF 使用确定性采样方法来解决此问题。GRV 再次接近状态分布，但现在我们使用由谨慎选择的采样点（叫做 Sigma 点）构成的极小集来表征。当 Sigma 点通过非线性系统传播时，它们可以捕捉到任何非线性特征的 3 阶精度的后验均值和协方差（泰勒级数展开）。与此相反，EKF 仅达到一阶精度。值得注意的是，UKF 的计算复杂性与 EKF 的相同。

如前文所述，我们会利用 UKF 来更新裂纹扩展初始模型。为此，需要在 UKF 框架中抛出此问题。我们将需要更新的参数 C 和 m 视为 UKF 的 “状态”。按照惯例，状态用 x 表示，并建模为随机漫步过程。之前描述的裂纹扩展模型为测量方程。综上所述，UKF 的过程和测量模型可以表示为：

其中 a 为裂纹尺寸， u 为周期数， h 代表裂纹扩展模型， w 为过程噪声， v 为测量噪声。假定过程噪声和测量噪声是均值为零的高斯白噪声。

基于 UKF 的模型更新方法中的步骤总结如下：

• 初始化状态向量及其不确定性（协方差）

• 生成基于状态的 Sigma 点

• 传播 Sigma 点并计算状态的预测平均值和协方差

• 计算基于预测平均值和协方差的 Sigma 点

• 计算每个 Sigma 点的输出值

• 将预测输出值作为第 5 步单个输出值的加权和计算

• 计算预测输出值的不确定性（协方差）

• 计算预测状态和预测输出值之间的互协方差

• 将新息作为测得的输出值和预测输出值之间的差计算

• 使用互协方差和输出协方差计算滤波增益

• 更新状态向量及其协方差

• 针对更新模型使用的全部数据点重复第 2 到 11 步

有兴趣的读者可以通过

https://en.wikipedia.org/wiki/Kalman_filter#Unscented_Kalman_filter，了解 UKF 背后的数学知识。

我们已将基于 UKF 的裂纹扩展模型更新方法作为 Depend On Docker 项目实施，您可以在

https://github.com/bhgedigital/ModelUpdate 中找到该项目，也可以在 Google Colab 中查看。为清晰起见，下面只提供了算法关键步骤的代码片段。

生成 Sigma 点

可以使用下列代码在 TensorFlow 中生成 Sigma 点：

## Inputs: State Vector - x_prev, Covariance Matrix - cov_prev, # States - num_st, UKF Parameter - kappa
## Output: Sigma Points Matrix - sig

[eigval, eigvec] = tf.linalg.eigh(cov_prev)                                   # Eigen Decomposition of the covariance matrix        
S_tf = tf.diag(tf.sqrt(eigval))                                              
sqrt_cov = tf.matmul(tf.matmul(eigvec, S_tf), tf.matrix_inverse(eigvec))      # Square-root of the covariance matrix
eta = tf.sqrt((alpha ** 2) * (num_st + kappa))                                # UKF scaling factor
sqrt_sc = tf.scalar_mul(eta, sqrt_cov)                                        # Scaled square-root
sig = np.matlib.repmat(x_prev, 1, 2*num_st+1)
sig[:, 1:num_st+1] = sig[:, 1:num_st+1] + sqrt_sc
sig[:, num_st+1:2*num_st+1] = sig[:, num_st+1:2*num_st+1] - sqrt_sc           # Sigma points

预测平均值和协方差

在生成 Sigma 点后，将通过过程模型传播这些点。然后，经过传播的 Sigma 点将用于计算预测平均值和协方差。

## Inputs:  Sigma Point Matrix - sig, UKF Weights - w0_mean, w0_cov, wi, # States - num_st, Process covariance - Q
## Outputs: Predicted State Vector - x_minus, Predicted Covariance - cov_minus

x_minus = w0_mean * sig[:, 0] + w_i * tf.reduce_sum(sig[:, 1:2 * num_st + 1], axis=1)
x_minus = tf.reshape(x_minus, (num_st, 1))
temp1 = tf.reshape(sig[:, 0], (num_st, 1))
# Predict Covariance
cov_minus = w0_cov * (temp1 - x_minus) * tf.transpose((temp1 - x_minus))
for i in range(1, 2 * num_st + 1):
 temp1 = tf.reshape(sig[:, i], (num_st, 1))
 cov_minus = cov_minus + w_i * (temp1 - x_minus) * tf.transpose((temp1 - x_minus))
 cov_minus = cov_minus + Q

计算所有 Sigma 点的输出值

所有 Sigma 点输出值的计算方法如下：

## Inputs: Sigma Point Matrix - sig, Current Crack Size - init_crack, Stress Constant - dsig,
##         Cycles Interval - diff_cycles, # States - num_st
## Output: Output Matrix - gam

gam = np.zeros((1, 2*num_st+1))
for indx in range(0, 2*num_st+1):
 logC = sig[0, indx]
 m = sig[1, indx]
 gam[:, indx] = predict_crack(diff_cycles, logC, m, init_crack, dsig)

计算输出值、互协方差和输出协方差

预测输出值、互协方差和输出协方差的计算方法如下：

## Inputs:  Output Matrix - gam, Sigma Point Matrix - sig, State Vector - x_prev, UKF Parameters - w0_mean, w0_cov, wi,
##          # States - num_st, Measurement Covariance - R
## Outputs: Cross Covariance - cov_xy, Output Covariance - cov_yy, Output Vector - yhat

yhat = w0_mean * gam[:, 0] + w_i * tf.reduce_sum(gam[:, 1:2 * num_st + 1])            # Output vector
temp2 = tf.reshape(sig[:, 0], (num_st, 1))
cov_yy = w0_cov * (gam[:, 0] - yhat) * tf.transpose(gam[:, 0] - yhat)                 # Output covariance
cov_xy = w0_cov * (temp2 - x_prev) * tf.transpose(gam[:, 0] - yhat)                   # Cross covariance
for i in range(1, 2*num_st+1):
 temp2 = tf.reshape(sig[:, i], (num_st, 1))
 cov_yy = cov_yy + w_i * (gam[:, i] - yhat) * tf.transpose(gam[:, i] - yhat)
 cov_xy = cov_xy + w_i * (temp2 - x_prev) * tf.transpose(gam[:, i] - yhat)

cov_yy = (cov_yy + tf.transpose(cov_yy))/2                                            # Ensure symmetry
cov_yy = cov_yy + R

更新状态向量及其协方差

最后，更新状态向量及其协方差。

## Inputs : Predicted State Vector - x_minus, Predicted State Covariance - cov_minus, Predicted Output - yhat,
#           Measured Output - y, Output Covariance cov_yy, Cross Covariance - cov_xy
## Outputs: Updated State Vector - x_hat, Updated State Covariance - cov_hat

k_gain = cov_xy * np.linalg.inv(np.matrix(cov_yy))                            # calculate filter gain
innov = y - y_hat                                                             # innovation
x_hat = x_minus + k_gain * innov                                              # update state
cov_hat = cov_minus - k_gain * cov_yy * tf.transpose(k_gain)                  # update covariance
cov_hat = (cov_hat + tf.transpose(cov_hat)) / 2                               # ensure symmetry

结果和讨论

如前文所述，更新数据集用于优化使用 UKF 的初始模型。我们发现更新后的模型参数为： logC = -22.2037 和 m = 3.5762。模型更新前后的预测如下图所示：

我们可以观察到，与初始模型相比，更新后模型的预测更加准确。但这是一个非常明显的结果，因为我们预测时使用的数据集与更新模型时使用的数据集相同。仅当验证预测时使用的数据集与更新模型时使用的数据集不同时，才能真正衡量更新后模型的准确度。更重要的是，更新后的模型应该能比初始模型更加准确地预测裂纹日后的演进趋势。否则，即使是更新后的模型也不会特别有用。

因此，我们使用更新后的模型预测了测试数据集对应周期的裂纹尺寸。更新后的模型和初始模型预测的裂纹尺寸如下图所示：

基于 UKF 的模型更新方法不仅可以估算模型参数的平均值，还可以计算与估算值相关的协方差。换言之，它可以提供状态变量的联合分布。我们利用 TensorFlow Probability 的 MultivariateNormalFullCovariance 函数创建了示例，用于计算与预测输出值相关的不确定性，即本例中的裂纹尺寸。上图还绘出了预测裂纹尺寸的 95% 可信预测区间。

从图中不难看出，在预测裂纹日后的演进趋势方面，更新后的模型明显比初始模型更加准确。这表明 UKF 能够从更新数据集的测量值中提取相关信息，以修改模型参数，从而使模型更加准确有用。这里的另一个注意事项是，此方法对噪声测量的稳健性。在更新数据集中，与裂纹的其他测量值相比，周期 600 和 700 对应的测量值是异常值。但这些异常值对模型更新过程的影响非常微小，因为 UKF 不会尝试将这些数据点的预测结果与观测值之间的误差最小化。这正是我们期望的结果，要获得准确的更新后模型，这种对噪声测量的稳健性至关重要。

最后，我们使用初始模型和更新后模型预测到 2200 周期的裂纹尺寸。预测结果和更新后模型的不确定性界限如下图所示。

从上述图表中，我们可以明显看出，如果使用初始模型，裂纹预测会非常不准确。事实上，当找到大于修复阈值的裂纹的真实可能性明显很小 (<< 5 %) 时，初始模型会导致系统在 2100 周期生成错误警报。

总结

正如我们在此文章系列的第一篇文章中所述，我们介绍了如何解决 “已知的未知问题”。在第一篇文章中，我们介绍了如何将领域知识与传统 ML 相结合，计算关于材料属性的已知的已知问题。这里展示的是更现实环境中的相同示例，其中的所有数据都非预先可用。本文中讨论的基于 UKF 的预测模型更新方法在数据有限和嘈杂的实际应用中非常有效。该方法可以兼容模型中的非线性特征，而且计算效率很高。

后续步骤

我们在第一篇文章构建了这些概念，在本文将已知的已知问题扩展到了已知的未知问题。下一篇将会介绍未知的未知问题，从而将三篇整合到一起。在接下来的文章中，我们会解决几个关于贝叶斯模型选择和模型优化的开放式问题。

致谢

本文是 Google 和 BHGE 多个团队辛勤工作与深入协作的成果。我们要特别感谢 Josh Dillon、Mike Shwe、Scott Fitzharris、Alex Walker、Shyam Sivaramakrishnan 和 Fabio Nonato，感谢他们多次编辑并提供诸多研究成果和代码片段，最重要的是，感谢他们的热心支持。

更多 AI 相关阅读：