作者 Arun Subramaniyan，BHGE Digital 数据科学与分析副总裁

贝克休斯，通用电气公司 (BHGE)，是全球领先的全方位油气公司，其使命是寻找更好的方法为世界提供能源。BHGE Digital 团队开发了企业级、人工智能驱动的 SaaS 解决方案，以提高效率并减少石油和天然气行业的非生产性时间。考虑任务关键型问题，例如预测燃气轮机的故障或优化大型系统（如石化工厂）；这些问题需要构建和维护大规模的复杂分析。我们制定了以分析为驱动的战略，为我们的客户实现全公司范围的数字化转型。

通过多年来解决工业界面临的最棘手问题，我们了解到最优雅、最持久的解决方案存在于以下交叉点：

• 领域知识
• 传统机器学习 (ML)
• 概率技术
• 深度学习

以前无法解决的问题类别现在可以通过结合看似无关的技术并使用现代可扩展的软件和硬件来解决。我们与 TensorFlow Probability (TFP) 团队以及 Google 的 Cloud ML 团队的合作加快了我们在大规模开发和部署这些技术的旅程。

我们打算通过本系列博客展示这些交叉点上发生的创新，并希望激励工业应用概率深度学习技术的寒武纪大爆发。我们在 2018 年 Google Next 上展示的这些应用程序的示例集可以在这里查看 这里.

对概率深度学习的需求

基于物理的（即基于领域的）分析已成功地用于航空航天、汽车和石油和天然气等各种行业的系统设计和运营数十年。它们提供了一种强大的方法，可以从少量观察中概括复杂的行为。牛顿的运动方程可用于精确预测原子和星系的运动。但是，现实世界问题的创新解决方案需要将这些定律与有效的方法相结合，以解释未知因素。明确跟踪我们信念的确定性对于可辩护的科学至关重要。

例如，以“简单”的抛射体轨迹为例。为了预测抛射体着陆的位置，唯一需要进行准确预测的观察（“数据”）是抛射体的速度（速度和角度）。但是，在初始角度和风速中添加少量不确定性（~5%）会导致抛射体着陆位置产生很大的不确定性（~200%），如下图所示。
即使对于一个简单的系统，在存在不确定性的情况下进行准确预测也很困难；想象一下优化大型非线性系统的复杂性！

数字孪生是 BHGE Analytics 团队多年来解决实际工业问题而磨练出的一个概念。我们将其定义为物理系统的数字表示，不断调整以表示当前条件并预测未来状态。如下图所示，构建数字孪生所需的三个支柱是领域知识、概率推理和深度学习

• 领域专业知识（即物理学）与传统 ML 相结合，使我们能够精确地解决已知问题，即已知已知。我们指的是诸如多项式回归、核密度方法和状态空间估计方法（例如卡尔曼滤波器）之类的技术。预测在各种热机械载荷下组件裂纹的长度是领域专业知识的一个例子：它需要对问题物理学的透彻理解以及对边界条件（如温度和压力）的了解。即使现象已被充分理解，也需要几种传统 ML 技术来计算特定于问题（例如材料特性）的系数。
• 概率推理提供了一种量化不确定性的系统方法，即已知未知。在上面的例子中，除了了解现象之外，我们还需要考虑测量中的不确定性以及制造可变性等未建模因素。众所周知，系统中存在会影响裂纹生长的可变性（不确定性）。然后，预测裂纹长度成为不确定性量化练习。更具体地说，当我们将不确定性添加到测量值时，材料特性校准成为概率推理问题。
• 深度学习和现代 ML 有潜力识别和预测未知模式和行为，即未知未知。在裂纹扩展示例中，即使将概率推理与领域模型结合起来，我们也只能在具有已知不确定性的特定载荷条件下预测裂纹扩展。考虑一个基于自动编码器的异常检测模型，该模型监视载荷条件以及设备上的其他各种条件。这个深度学习模型可以捕捉基于物理的模型无法识别的异常条件。我们将此预测称为未知未知的原因是我们没有关于异常的任何数据来训练模型。虽然深度学习模型在正常操作数据（以及一些用于验证的异常）上进行训练，但它会捕获任何以前未观察到的偏离正常行为的行为。这是深度学习模型派上用场的众多示例之一。有人可能会争辩说，在某些情况下，基于传统技术（例如 PCA 重构误差）的简单异常检测模型也能奏效。根据我们的经验，当已知过程特征与这些方法的简化假设一致时，更简单的技术可以提供与深度学习模型类似甚至更好的性能，从而导致已知异常。当异常事先未知时——当您以前没有见过特定的故障模式时——深度学习模型真正发挥作用。

在本文中，我们将重点介绍带有领域模型的概率推理，以预测已知未知。我们将用一个公式为基于物理的概率推理模型的裂纹扩展问题证明贝叶斯校准的强大功能。

概率深度学习使我们能够在一个“自学习”包中利用上面突出显示的所有功能。我们在随后的博客中将使用 TFP 讨论概率深度学习，这些博客将涵盖模型差异、异常检测、缺失数据估计和时间序列预测。

预测组件的寿命：概率裂纹扩展示例

预测易于开裂的组件的寿命是一个由来已久的问题，一直被断裂力学界反复研究。裂纹扩展模型是工程系统预后与健康管理 (PHM) 解决方案的核心。恰如其分地名为工程系统预后与健康管理：导论一书提供了一个很好的例子，说明如何使用现实世界的数据来校准工程模型。通过以下示例，我们希望推动使用结合概率学习技术和工程领域模型的“混合模型”。

疲劳裂纹扩展现象可以用Paris 定律建模。Paris 定律通过以下方程将裂纹扩展速率 (da/dN) 与应力强度因子 (ΔK=Δσ√𝜋a) 联系起来

其中a是裂纹长度，N是载荷循环次数，σ是应力，(C, m) 是材料特性。

对于特定的几何形状和载荷配置，对 Paris 定律进行积分，我们就可以得到裂纹尺寸作为载荷循环函数的解析公式，如下所示

其中aₒ是初始裂纹长度。通过了解aₒ和给定应用程序上的因子Δσ√𝜋，可以使用方程 2 计算未来裂纹的尺寸，假设应用程序将随着时间推移累积N个循环。这种预测的裂纹长度可以与安全运行的阈值进行比较。例如，如果预测的裂纹长度超过几何极限（例如，组件厚度的二分之一），则需要进行维修。

参数C和m需要针对每个物理组件进行校准，给定裂纹长度 a 与载荷循环N数据。换句话说，给定现场测量的裂纹，我们希望推断C和m，以便我们能够估计在裂纹长度变得足够大以致造成组件失效之前，我们能够安全运行的循环次数。

在本例中，我们将使用Kim、An 和 Choi 撰写的 PHM 书籍中的样本数据集，演示使用 TFP 对C和m进行概率校准。在 BHGE Digital，我们利用Depend-on-Docker 项目来自动化分析开发。此处提供此类自动化的示例此处，以及以下示例的完整代码。

如下图所示，我们有Kim、An 和 Choi 撰写的 PHM 书籍第 4.2 表中提供的数据集。正如大多数裂纹扩展数据集的典型情况一样，从观察到的数据点中并不完全清楚潜在的趋势。

对于贝叶斯校准，我们需要定义校准变量的先验分布。在实际应用中，这些先验信息可以由主题专家提供。对于本示例，我们将假设C和m都是高斯分布且独立的

在 TFP 中，我们可以使用以下方法编码这些信息

prio_par_logC = [-23., 1.1] # [location, scale] for Normal Prior
prio_par_m = [4., 0.2] # [location, scale] for Normal Prior
rv_logC = tfd.Normal(loc=0., scale=1., name='logC_norm')
rv_m = tfd.Normal(loc=0., scale=1., name='m_norm')

为了从归一化空间采样，我们已经为两个变量定义了外部参数和标准正态分布。因此，在计算裂纹模型时，我们需要对两个随机变量进行反归一化。

现在我们定义了被校准的随机变量的联合对数概率以及由公式 2 定义的关联裂纹模型。

def joint_log_prob(cycles, observations, y0, logC_norm, m_norm):
    # Joint logProbability function for both random variables and observations.
    # Some constants
    dsig = 75.
    B = tf.constant(dsig * np.pi**0.5, tf.float32)
    # Computing m and logC on original space
    logC = logC_norm * prio_par_logC[1]**0.5+ prio_par_logC[0] # 
    m = m_norm * prio_par_m[1]**0.5 + prio_par_m[0]

    # Crack Propagation model
    crack_model =(cycles * tf.exp(logC) * (1 - m / 2.) * B**m + y0**(1-m / 2.))**(2. / (2. - m))
    y_model = observations - crack_model


    # Defining child model random variable
    rv_model = tfd.Independent(
        tfd.Normal(loc=tf.zeros(observations.shape), scale=0.001),
       reinterpreted_batch_ndims=1, name = 'model')
    # Sum of logProbabilities
    return rv_logC.log_prob(logC_norm) + rv_m.log_prob(m_norm) + rv_model.log_prob(y_model)

最后，是时候设置采样器并运行 TensorFlow 会话了。

# This cell can take 12 minutes to run in Graph mode
# Number of samples and burnin for the MCMC sampler
samples = 10000
burnin = 10000

# Initial state for the HMC
initial_state = [0., 0.]
# Converting the data into tensors
cycles = tf.convert_to_tensor(t_,tf.float32)
observations = tf.convert_to_tensor(y_,tf.float32)
y0 = tf.convert_to_tensor(y_[0], tf.float32)
# Setting up a target posterior for our joint logprobability
unormalized_target_posterior= lambda *args: joint_log_prob(cycles, observations, y0, *args)
# And finally setting up the mcmc sampler
[logC_samples, m_samples], kernel_results = tfp.mcmc.sample_chain(
    num_results= samples, 
    num_burnin_steps= burnin,
    current_state=initial_state,
    kernel= tfp.mcmc.HamiltonianMonteCarlo(
        target_log_prob_fn=unormalized_target_posterior,
        step_size = 0.045, 
        num_leapfrog_steps=6))


# Tracking the acceptance rate for the sampled chain
acceptance_rate = tf.reduce_mean(tf.to_float(kernel_results.is_accepted))

# Actually running the sampler
# The evaluate() function, defined at the top of this notebook, runs `sess.run()` 
# in graph mode and allows code to be executed eagerly when Eager mode is enabled
[logC_samples_, m_samples_, acceptance_rate_] = evaluate([
    logC_samples, m_samples, acceptance_rate])

# Some initial results
print('acceptance_rate:', acceptance_rate_)

有趣的是，尽管我们最初为C和m使用了两个独立的正态分布，但后验分布却高度相关。这种相关性是由于解空间要求对于m的高值，唯一物理上有意义的结果位于C的小值处，反之亦然。如果我们要使用任何数量的确定性优化技术来找到适合该数据集的C和m的“最佳拟合”，我们将最终得到一些取决于我们的起点和约束条件的直线上的值。执行概率优化（又称贝叶斯校准）为我们提供了可以解释数据集的所有可能解的全局视图。

对预测进行后验采样

现在是最后的行动，我们将定义一个后验函数作为对数正态分布，其期望由物理模型定义。

众所周知，损伤（在本例中为裂纹长度）是偏斜的，因此通常使用对数正态分布或 Gumbel 分布来模拟损伤。在 TFP 中表达模型如下所示。

def posterior(logC_samples, m_samples, time):
    n_s = len(logC_samples)
    n_inputs = len(time)

    # Some Constants
    dsig = 75.
    B = tf.constant(dsig * np.pi**0.5, tf.float32)

    # Crack Propagation model - compute in the log space

    y_model =(
        time[:,None]  *
        tf.exp(logC_samples[None,:])*
        (1-m_samples[None,:]/2.0) *B**m_samples[None,:] +y0** (1-m_samples[None,:]/2.0))**(2. / (2. - m_samples[None,:]))
    noise = tfd.Normal(loc=0., scale=0.001)
    samples = y_model + noise.sample(n_s)[tf.newaxis,:]
    # The evaluate() function, defined at the top of this notebook, runs `sess.run()` 
    # in graph mode and allows code to be executed eagerly when Eager mode is enabled
    samples_ = evaluate(samples)
    return samples_

# Predict for a range of cycles
time = np.arange(0, 3000, 100)
y_samples = posterior(logC_samples_scale, m_samples_scale, time)
print(y_samples.shape)

下面显示了使用混合物理概率模型预测的裂纹长度的 95% 不确定性边界内的平均值。显然，该模型不仅捕获了平均行为，而且还估计了每个时间点的模型预测的不确定性。随着我们远离观测值，模型预测的不确定性呈指数增长。

下一步

我们仔细选择此示例是为了提供对概率模型的温和介绍。将概率模型应用于现实世界应用程序存在一些挑战。即使在像这里显示的这样一个简单的问题中，如果我们放宽先验是高斯的假设，解就会变得难以处理。更具探究性的读者可以尝试使用均匀先验来查看模型预测能力的显着下降。

另一个微妙之处是似然函数的定义方式。我们选择在预测误差而不是直接预测值上定义它。请注意，下面“基于预测误差的似然”中描绘的各个蓝色线是模型预测误差的样本。将公式更改为直接对模型预测进行建模——而不是预测误差——会产生非单调的、非物理的结果，如“基于实际预测的似然”中所示。我们称这些结果为非物理结果，因为裂纹长度不能随时间减少。如果我们只查看上面显示的“模糊”百分位视图，我们可能不会注意到模型公式的细微差异，这会导致巨大的预测误差。我们将在接下来的博客中讨论一些更实际的挑战及其缓解措施。

这是旨在扩展使用 TFP 用于工业应用的概率和深度学习技术的博客系列中的第一篇。我们 (@sarunkarthi) 非常乐意了解您的应用，并期待看到这些方法以独特的方式使用。请继续关注此博客供稿，以获取有关使用变分推断进行异常检测、缺失数据估计和预测的更多更新和示例。

致谢

这篇博客是 Google 和 BHGE 多个团队辛勤工作和深入合作的结果。我们特别感谢 Josh Dillon、Mike Shwe、Scott Fitzharris、Alex Walker、Christian Hillaire 和 Fabio Nonato 的众多编辑、结果、代码片段以及最重要的是他们的热情支持。