通用電氣的貝克休斯公司 Baker Hughes（BHGE）是世界領先的全流程石油和天然氣公司，致力于尋求更佳的方式向世界輸送能源。BHGE 的數字團隊開發了企業級別的，依靠人工智能驅動的 SaaS 解決方案，以提高效率并減少石油和天然氣行業的非生產時間。要考慮關鍵任務問題，諸如預測燃氣輪機故障或優化石化廠等大型的系統；這些問題需要大規模構建和維護復雜的分析。為此，我們開發了一種分析驅動型戰略計劃，為客戶實現全公司數字化改造。

多年來一直幫助工業界解決最棘手問題，使我們掌握了最優雅持久的解決方案：

領域知識

傳統機器學習

概率技術

深度學習

之前難以處理的問題類別現在可以通過將一些看似毫無關系的技術并將它們與現代可擴展的軟硬件組合在一起部署來解決。我們與 TensorFlowProbability（TFP）團隊和 Google 的 Cloud ML 團隊的合作加速了我們大規模開發和部署這些技術的過程。

我們想將這些節點發生的創新通過這一系列的文章進行展示，希望以此激發在概率深度學習技術的工業應用中實現爆炸性增長。在這里你們可以看到我們在 Google Next 2018 展示的這些應用程序示例集。

需要概率深度學習

數十年來，基于物理的（即基于領域的）分析已經被成功地應用于設計和操作航空航天、汽車、石油和天然氣等行業的系統。它們提供了一個行之有效的方法，能夠從一些觀察中概括復雜的行為。牛頓的運動方程可以用來精確地預測原子和星系的運動。然而，對現實世界問題的創新解決方案需要將這些律法與解決未知因素的有效方法相結合。要捍衛科學，明確地追蹤我們信念的確定性至關重要。

我們舉一個采用射彈軌跡的 “簡單” 例子。為了預測拋射物落地的位置，進行準確預測所需的唯一觀察（數據）是投射拋射物的速度（速度和角度）。然而，在初始角度和風速中增加一個很小的不確定性（~5％）會給射彈落地的位置帶來巨大的不確定性（~200％），如下圖所示。

即使是簡單的系統，也可以通過不確定性進行準確的預測；想象一下優化大型非線性系統得有多么復雜！

Digital Twin 是 BHGE Analytics 團隊在解決實際工業問題方面經過多年磨練的結構。我們將其定義為物理系統的數字表示，不斷調整用以表示當前條件并預測未來狀態。如下圖所示，構建數字孿生所需的三大支柱是領域知識，概率推理和深度學習：

領域專業知識（即物理學）與傳統 ML 相結合，使我們能夠精確地解決已知問題，換句話說，已知的知識。通過傳統的 ML，我們指的是諸如多項式回歸，核密度方法和狀態空間估計方法（例如卡爾曼濾波器）之類的技術。在各種熱機械載荷下預測零件中裂紋的長度是領域專業知識的一個例子：它需要徹底了解問題的物理特性以及溫度和壓力等邊界條件的知識。即使很好地理解了這種現象，也需要幾種傳統的 ML 技術來計算特定于問題的系數（例如材料特性）。

概率推理提供了一種量化不確定性的系統方法，即已知的未知數。在上面的例子中，除了理解這些現象外，我們還需要考慮測量的不確定性和制造變異性等非模型因素。眾所周知，系統中存在會影響裂紋擴展的可變性（不確定性）。預測裂縫長度然后成為不確定性量化練習。更具體地說，當我們在測量中增加不確定性時，材料特性校準成為概率推理問題。

深度學習和現代 ML 有潛力識別和預測未知的模式和行為，也就是未知的未知。在裂縫傳播實例中，即使概率推理與領域模型相結合，我們也可以預測裂縫僅在具有已知不確定性的特定負載條件下傳播。考慮基于自動編碼器的異常檢測模型，該自動編碼器監視負載條件以及設備上的各種其他條件。這種深度學習模型可以捕獲基于物理的模型無法獲取的異常情況。我們將此預測表示為未知的未知，原因是因為我們沒有關于異常的任何數據來訓練模型。雖然深度學習模型訓練正常操作的數據（以及一些驗證異常），但它可以捕捉到之前未觀察到的正常行為的任何偏差。這是深度學習模型派上用場的眾多例子之一。可以說，在某些情況下，基于傳統技術（如 PCA 重建誤差）的簡單異常檢測模型也可以解決問題。根據我們的經驗，當已知的過程特征與這些方法的簡化假設一致時，更簡單的技術可以提供與深度學習模型類似甚至更好的性能，從而導致已知的異常。當異常事件突如其來毫無征兆時 - 當你之前從沒有看到特定的故障模式時，深度學習模型的表現真的出類拔萃。

在本文中，我們將重點關注在預測已知未知數的領域模型的概率推斷。我們將演示貝葉斯校準的能力，其中裂縫傳播問題被公式化為基于物理的概率推理模型。

概率深度學習使我們能夠在 “自學” 包中利用上面強調的所有功能。我們將在隨后的博客中討論使用 TFP 進行概率深度學習，其中包括模型差異，異常檢測，缺失數據估計和時間序列預測。

預測組件的壽命：概率性裂縫傳播示例

預測易于開裂的部件的壽命是一個老生常談的問題，斷裂力學界已經研究過。裂縫傳播模型是工程系統的預測和健康管理（PHM）解決方案的核心。題為《預測與工程系統健康管理：簡介》一書提供了一個很好的例子，說明現實世界數據如何用于校準工程模型。通過下面的示例，我們希望激發使用結合概率學習技術和工程領域模型的 “混合模型”。

疲勞裂紋擴展現象可以用巴黎定律建模。巴黎定律通過以下等式將裂紋擴展速率（da / dN）與應力強度因子（ΔK=Δσ√）聯系起來：

其中 a 是裂紋長度，N 是加載循環次數，σ 是應力，（C，m）是材料屬性。

將巴黎定律與特定幾何和加載配置相結合，我們得出裂縫大小的分析公式，作為加載循環的函數，如下所示：

其中 a? 是初始裂縫長度。通過了解給定應用中的 a? 和因子 Δσ√，可以使用等式 2 來計算未來裂縫的大小，假設應用將隨時間累積 N 個周期。可以將該預測的裂縫長度與安全操作的閾值進行比較。例如，如果預測的裂縫長度超過幾何極限（例如，部件厚度的一半），那就說明該修理了。

在給定裂縫長度 a 與加載循環 N 數據的情況下，需要針對每個物理組件校準參數 C 和 m。換句話說，考慮到現場測量的裂縫，我們希望推斷出 C 和 m，這樣就可以估算出在裂縫長度變大到一定程度，在部件失效的風險之前，可以安全地運行多少個循環。

在本例中，我們將使用由 Kim，An 和 Choi 編寫的 PHM 書中的樣本數據集，演示使用 TFP 進行 C 和 m 的概率校準。在 BHGE Digital，我們利用 Depend-on-Docker 項目實現分析開發的自動化。此處提供了此類自動化的示例，其中包含以下示例的完整代碼。

如下圖所示，我們使用了 Kim，An 和 Choi 的 PHM 書中的表 4.2 中提供的數據集。對于大多數裂縫傳播數據集而言，從觀察到的數據點來看，潛在趨勢并不十分明顯。

對于貝葉斯校準，我們需要定義校準變量的先驗分布。在實際應用中，這些先驗可以由主題專家通知。對于這個例子，我們假設 C 和 m 都是高斯和獨立的：

在 TFP 中，我們可以對此信息進行如下編碼：

1prio_par_logC = [-23., 1.1] # [location, scale] for Normal Prior

2prio_par_m = [4., 0.2] # [location, scale] for Normal Prior

3rv_logC = tfd.Normal(loc=0., scale=1., name='logC_norm')

4rv_m = tfd.Normal(loc=0., scale=1., name='m_norm')

我們為兩個變量定義了外部參數和標準正態分布，只是為了從標準化空間中進行采樣。因此，在計算裂縫模型時，我們需要對兩個隨機變量進行去標準化。

現在我們定義被校準的隨機變量的聯合對數概率以及由等式 2 定義的相關裂縫模型：

1def joint_log_prob(cycles, observations, y0, logC_norm, m_norm):

2# Joint logProbability function for both random variables and observations.

3# Some constants

4dsig = 75.

5B = tf.constant(dsig * np.pi**0.5, tf.float32)

6# Computing m and logC on original space

7logC = logC_norm * prio_par_logC[1]**0.5+ prio_par_logC[0] #

8m = m_norm * prio_par_m[1]**0.5 + prio_par_m[0]

9

10# Crack Propagation model

11crack_model =(cycles * tf.exp(logC) * (1 - m / 2.) * B**m + y0**(1-m / 2.))**(2. / (2. - m))

12y_model = observations - crack_model

13

14

15# Defining child model random variable

16rv_model = tfd.Independent(

17 tfd.Normal(loc=tf.zeros(observations.shape), scale=0.001),

18reinterpreted_batch_ndims=1, name = 'model')

19# Sum of logProbabilities

20return rv_logC.log_prob(logC_norm) + rv_m.log_prob(m_norm) + rv_model.log_prob(y_model)

最后，是時候設置采樣器并運行 TensorFlow 會話了：

1# This cell can take 12 minutes to run in Graph mode

2# Number of samples and burnin for the MCMC sampler

3samples = 10000

4burnin = 10000

5

6# Initial state for the HMC

7initial_state = [0., 0.]

8# Converting the data into tensors

9cycles = tf.convert_to_tensor(t_,tf.float32)

10observations = tf.convert_to_tensor(y_,tf.float32)

11y0 = tf.convert_to_tensor(y_[0], tf.float32)

12# Setting up a target posterior for our joint logprobability

13unormalized_target_posterior= lambda *args: joint_log_prob(cycles, observations, y0, *args)

14# And finally setting up the mcmc sampler

15[logC_samples, m_samples], kernel_results = tfp.mcmc.sample_chain(

16num_results= samples,

17num_burnin_steps= burnin,

18current_state=initial_state,

19kernel= tfp.mcmc.HamiltonianMonteCarlo(

20target_log_prob_fn=unormalized_target_posterior,

21step_size = 0.045,

22num_leapfrog_steps=6))

23

24

25# Tracking the acceptance rate for the sampled chain

26acceptance_rate = tf.reduce_mean(tf.to_float(kernel_results.is_accepted))

27

28# Actually running the sampler

29# The evaluate() function, defined at the top of this notebook, runs `sess.run()

30# in graph mode and allows code to be executed eagerly when Eager mode is enabled

31[logC_samples_, m_samples_, acceptance_rate_] = evaluate([

32logC_samples, m_samples, acceptance_rate])

33

34# Some initial results

35print('acceptance_rate:', acceptance_rate_)

值得注意的是，盡管我們開始使用 C 和 m 的兩個獨立高斯分布，但后驗分布是高度相關的。這種相關性的產生是因為解空間規定對于 m 的高值，唯一有物理意義的結果是 Cand 的小值，反之亦然。如果我們使用任意數量的確定性優化技術來找到適合此數據集的 C 和 m 的 “最佳擬合”，那么根據我們的起點和約束，我們最終會得到一些位于直線上的值。執行概率優化（也就是貝葉斯校準）為我們提供了可以解釋數據集的所有可能解決方案的全局視圖。

后驗抽樣后進行預測

現在到了最后一步，我們將后驗函數定義為對數正態分布，期望由物理模型定義：

已知損壞（在這種情況下為裂縫長度）是傾斜的，因此對數正常或 Gumbel 分布通常用于模擬損傷。在 TFP 中表達模型如下所示：

1def posterior(logC_samples, m_samples, time):

2n_s = len(logC_samples)

3n_inputs = len(time)

4

5# Some Constants

6dsig = 75.

7B = tf.constant(dsig * np.pi**0.5, tf.float32)

8

9# Crack Propagation model - compute in the log space

10

11y_model =(

12time[:,None] *

13tf.exp(logC_samples[None,:])*

14(1-m_samples[None,:]/2.0) *B**m_samples[None,:] +y0** (1-m_samples[None,:]/2.0))**(2. / (2. - m_samples[None,:]))

15noise = tfd.Normal(loc=0., scale=0.001)

16samples = y_model + noise.sample(n_s)[tf.newaxis,:]

17# The evaluate() function, defined at the top of this notebook, runs `sess.run()`

18# in graph mode and allows code to be executed eagerly when Eager mode is enabled

19samples_ = evaluate(samples)

20return samples_

21

22# Predict for a range of cycles

23time = np.arange(0, 3000, 100)

24y_samples = posterior(logC_samples_scale, m_samples_scale, time)

25print(y_samples.shape)

如下所示，使用混合物理概率模型的裂縫長度的 95％ 不確定界限的預測平均值。很顯然，該模型不僅捕獲平均行為，還捕獲每個時間點模型預測的不確定性估算。當我們偏離觀察時，模型預測中的不確定性呈指數增長趨勢。

下一步

我們選擇這個例子十分小心，來對概率模型做一個介紹。將概率模型應用于實際應用程序存在一些挑戰。即使在像這里所展示的簡單問題中，如果我們放松高斯先驗的假設，那么解決方案也會變得棘手。如果讀者感到好奇，可以嘗試使用統一先驗來觀察一下模型預測能力急劇下降的情況。

另一個微妙點是定義可能性的方式。我們選擇直接在預測誤差上定義它，而不是預測值。請注意，下面在 “基于預測誤差的可能性” 中描述的各個藍線是模型預測誤差的樣本。直接將公式更改為模型預測 - 而不是預測誤差 - 會產生非單調的非物理結果，如 “基于實際預測的可能性 ” 所示。我們稱這些結果是非物理的，因為裂縫長度不會隨著時間的推移而減小。如果我們只查看上面顯示的 “模糊” 百分位數視圖，我們可能不會注意到模型公式中可能導致較大預測誤差的細微差別。我們將在下一篇文章中解決其中一些更實際的挑戰及其緩解措施。

這是系列文章中的第一篇，旨在通過 TFP 擴展概率和深度學習技術在工業應用中的應用。我們（@sarunkarthi）很想聽聽關于您的應用程序的說明，并期待看到這些方法應用到更多獨特的方方面面。請繼續關注此文章供稿，了解有關異常檢測，缺失數據估算和變異推斷預測的更多更新的示例。