如何用TensorFlow概率編程工具包開(kāi)發(fā)深度學(xué)習(xí)

通用電氣的貝克休斯公司 Baker Hughes（BHGE）是世界領(lǐng)先的全流程石油和天然氣公司，致力于尋求更佳的方式向世界輸送能源。BHGE 的數(shù)字團(tuán)隊(duì)開(kāi)發(fā)了企業(yè)級(jí)別的，依靠人工智能驅(qū)動(dòng)的 SaaS 解決方案，以提高效率并減少石油和天然氣行業(yè)的非生產(chǎn)時(shí)間。要考慮關(guān)鍵任務(wù)問(wèn)題，諸如預(yù)測(cè)燃?xì)廨啓C(jī)故障或優(yōu)化石化廠等大型的系統(tǒng)；這些問(wèn)題需要大規(guī)模構(gòu)建和維護(hù)復(fù)雜的分析。為此，我們開(kāi)發(fā)了一種分析驅(qū)動(dòng)型戰(zhàn)略計(jì)劃，為客戶實(shí)現(xiàn)全公司數(shù)字化改造。

多年來(lái)一直幫助工業(yè)界解決最棘手問(wèn)題，使我們掌握了最優(yōu)雅持久的解決方案：

領(lǐng)域知識(shí)

傳統(tǒng)機(jī)器學(xué)習(xí)

概率技術(shù)

深度學(xué)習(xí)

之前難以處理的問(wèn)題類別現(xiàn)在可以通過(guò)將一些看似毫無(wú)關(guān)系的技術(shù)并將它們與現(xiàn)代可擴(kuò)展的軟硬件組合在一起部署來(lái)解決。我們與 TensorFlowProbability（TFP）團(tuán)隊(duì)和 Google 的 Cloud ML 團(tuán)隊(duì)的合作加速了我們大規(guī)模開(kāi)發(fā)和部署這些技術(shù)的過(guò)程。

我們想將這些節(jié)點(diǎn)發(fā)生的創(chuàng)新通過(guò)這一系列的文章進(jìn)行展示，希望以此激發(fā)在概率深度學(xué)習(xí)技術(shù)的工業(yè)應(yīng)用中實(shí)現(xiàn)爆炸性增長(zhǎng)。在這里你們可以看到我們?cè)?Google Next 2018 展示的這些應(yīng)用程序示例集。

需要概率深度學(xué)習(xí)

數(shù)十年來(lái)，基于物理的（即基于領(lǐng)域的）分析已經(jīng)被成功地應(yīng)用于設(shè)計(jì)和操作航空航天、汽車(chē)、石油和天然氣等行業(yè)的系統(tǒng)。它們提供了一個(gè)行之有效的方法，能夠從一些觀察中概括復(fù)雜的行為。牛頓的運(yùn)動(dòng)方程可以用來(lái)精確地預(yù)測(cè)原子和星系的運(yùn)動(dòng)。然而，對(duì)現(xiàn)實(shí)世界問(wèn)題的創(chuàng)新解決方案需要將這些律法與解決未知因素的有效方法相結(jié)合。要捍衛(wèi)科學(xué)，明確地追蹤我們信念的確定性至關(guān)重要。

我們舉一個(gè)采用射彈軌跡的 “簡(jiǎn)單” 例子。為了預(yù)測(cè)拋射物落地的位置，進(jìn)行準(zhǔn)確預(yù)測(cè)所需的唯一觀察（數(shù)據(jù)）是投射拋射物的速度（速度和角度）。然而，在初始角度和風(fēng)速中增加一個(gè)很小的不確定性（~5％）會(huì)給射彈落地的位置帶來(lái)巨大的不確定性（~200％），如下圖所示。

即使是簡(jiǎn)單的系統(tǒng)，也可以通過(guò)不確定性進(jìn)行準(zhǔn)確的預(yù)測(cè)；想象一下優(yōu)化大型非線性系統(tǒng)得有多么復(fù)雜！

Digital Twin 是 BHGE Analytics 團(tuán)隊(duì)在解決實(shí)際工業(yè)問(wèn)題方面經(jīng)過(guò)多年磨練的結(jié)構(gòu)。我們將其定義為物理系統(tǒng)的數(shù)字表示，不斷調(diào)整用以表示當(dāng)前條件并預(yù)測(cè)未來(lái)狀態(tài)。如下圖所示，構(gòu)建數(shù)字孿生所需的三大支柱是領(lǐng)域知識(shí)，概率推理和深度學(xué)習(xí)：

領(lǐng)域?qū)I(yè)知識(shí)（即物理學(xué)）與傳統(tǒng) ML 相結(jié)合，使我們能夠精確地解決已知問(wèn)題，換句話說(shuō)，已知的知識(shí)。通過(guò)傳統(tǒng)的 ML，我們指的是諸如多項(xiàng)式回歸，核密度方法和狀態(tài)空間估計(jì)方法（例如卡爾曼濾波器）之類的技術(shù)。在各種熱機(jī)械載荷下預(yù)測(cè)零件中裂紋的長(zhǎng)度是領(lǐng)域?qū)I(yè)知識(shí)的一個(gè)例子：它需要徹底了解問(wèn)題的物理特性以及溫度和壓力等邊界條件的知識(shí)。即使很好地理解了這種現(xiàn)象，也需要幾種傳統(tǒng)的 ML 技術(shù)來(lái)計(jì)算特定于問(wèn)題的系數(shù)（例如材料特性）。

概率推理提供了一種量化不確定性的系統(tǒng)方法，即已知的未知數(shù)。在上面的例子中，除了理解這些現(xiàn)象外，我們還需要考慮測(cè)量的不確定性和制造變異性等非模型因素。眾所周知，系統(tǒng)中存在會(huì)影響裂紋擴(kuò)展的可變性（不確定性）。預(yù)測(cè)裂縫長(zhǎng)度然后成為不確定性量化練習(xí)。更具體地說(shuō)，當(dāng)我們?cè)跍y(cè)量中增加不確定性時(shí)，材料特性校準(zhǔn)成為概率推理問(wèn)題。

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

在本文中，我們將重點(diǎn)關(guān)注在預(yù)測(cè)已知未知數(shù)的領(lǐng)域模型的概率推斷。我們將演示貝葉斯校準(zhǔn)的能力，其中裂縫傳播問(wèn)題被公式化為基于物理的概率推理模型。

概率深度學(xué)習(xí)使我們能夠在 “自學(xué)” 包中利用上面強(qiáng)調(diào)的所有功能。我們將在隨后的博客中討論使用 TFP 進(jìn)行概率深度學(xué)習(xí)，其中包括模型差異，異常檢測(cè)，缺失數(shù)據(jù)估計(jì)和時(shí)間序列預(yù)測(cè)。

預(yù)測(cè)組件的壽命：概率性裂縫傳播示例

預(yù)測(cè)易于開(kāi)裂的部件的壽命是一個(gè)老生常談的問(wèn)題，斷裂力學(xué)界已經(jīng)研究過(guò)。裂縫傳播模型是工程系統(tǒng)的預(yù)測(cè)和健康管理（PHM）解決方案的核心。題為《預(yù)測(cè)與工程系統(tǒng)健康管理：簡(jiǎn)介》一書(shū)提供了一個(gè)很好的例子，說(shuō)明現(xiàn)實(shí)世界數(shù)據(jù)如何用于校準(zhǔn)工程模型。通過(guò)下面的示例，我們希望激發(fā)使用結(jié)合概率學(xué)習(xí)技術(shù)和工程領(lǐng)域模型的 “混合模型”。

疲勞裂紋擴(kuò)展現(xiàn)象可以用巴黎定律建模。巴黎定律通過(guò)以下等式將裂紋擴(kuò)展速率（da / dN）與應(yīng)力強(qiáng)度因子（ΔK=Δσ√）聯(lián)系起來(lái)：

其中 a 是裂紋長(zhǎng)度，N 是加載循環(huán)次數(shù)，σ 是應(yīng)力，（C，m）是材料屬性。

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

其中 a? 是初始裂縫長(zhǎng)度。通過(guò)了解給定應(yīng)用中的 a? 和因子 Δσ√，可以使用等式 2 來(lái)計(jì)算未來(lái)裂縫的大小，假設(shè)應(yīng)用將隨時(shí)間累積 N 個(gè)周期?？梢詫⒃擃A(yù)測(cè)的裂縫長(zhǎng)度與安全操作的閾值進(jìn)行比較。例如，如果預(yù)測(cè)的裂縫長(zhǎng)度超過(guò)幾何極限（例如，部件厚度的一半），那就說(shuō)明該修理了。

在給定裂縫長(zhǎng)度 a 與加載循環(huán) N 數(shù)據(jù)的情況下，需要針對(duì)每個(gè)物理組件校準(zhǔn)參數(shù) C 和 m。換句話說(shuō)，考慮到現(xiàn)場(chǎng)測(cè)量的裂縫，我們希望推斷出 C 和 m，這樣就可以估算出在裂縫長(zhǎng)度變大到一定程度，在部件失效的風(fēng)險(xiǎn)之前，可以安全地運(yùn)行多少個(gè)循環(huán)。

在本例中，我們將使用由 Kim，An 和 Choi 編寫(xiě)的 PHM 書(shū)中的樣本數(shù)據(jù)集，演示使用 TFP 進(jìn)行 C 和 m 的概率校準(zhǔn)。在 BHGE Digital，我們利用 Depend-on-Docker 項(xiàng)目實(shí)現(xiàn)分析開(kāi)發(fā)的自動(dòng)化。此處提供了此類自動(dòng)化的示例，其中包含以下示例的完整代碼。

如下圖所示，我們使用了 Kim，An 和 Choi 的 PHM 書(shū)中的表 4.2 中提供的數(shù)據(jù)集。對(duì)于大多數(shù)裂縫傳播數(shù)據(jù)集而言，從觀察到的數(shù)據(jù)點(diǎn)來(lái)看，潛在趨勢(shì)并不十分明顯。

對(duì)于貝葉斯校準(zhǔn)，我們需要定義校準(zhǔn)變量的先驗(yàn)分布。在實(shí)際應(yīng)用中，這些先驗(yàn)可以由主題專家通知。對(duì)于這個(gè)例子，我們假設(shè) C 和 m 都是高斯和獨(dú)立的：

在 TFP 中，我們可以對(duì)此信息進(jìn)行如下編碼：

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')

我們?yōu)閮蓚€(gè)變量定義了外部參數(shù)和標(biāo)準(zhǔn)正態(tài)分布，只是為了從標(biāo)準(zhǔn)化空間中進(jìn)行采樣。因此，在計(jì)算裂縫模型時(shí)，我們需要對(duì)兩個(gè)隨機(jī)變量進(jìn)行去標(biāo)準(zhǔn)化。

現(xiàn)在我們定義被校準(zhǔn)的隨機(jī)變量的聯(lián)合對(duì)數(shù)概率以及由等式 2 定義的相關(guān)裂縫模型：

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)

最后，是時(shí)候設(shè)置采樣器并運(yùn)行 TensorFlow 會(huì)話了：

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_)

值得注意的是，盡管我們開(kāi)始使用 C 和 m 的兩個(gè)獨(dú)立高斯分布，但后驗(yàn)分布是高度相關(guān)的。這種相關(guān)性的產(chǎn)生是因?yàn)榻饪臻g規(guī)定對(duì)于 m 的高值，唯一有物理意義的結(jié)果是 Cand 的小值，反之亦然。如果我們使用任意數(shù)量的確定性優(yōu)化技術(shù)來(lái)找到適合此數(shù)據(jù)集的 C 和 m 的 “最佳擬合”，那么根據(jù)我們的起點(diǎn)和約束，我們最終會(huì)得到一些位于直線上的值。執(zhí)行概率優(yōu)化（也就是貝葉斯校準(zhǔn)）為我們提供了可以解釋數(shù)據(jù)集的所有可能解決方案的全局視圖。

后驗(yàn)抽樣后進(jìn)行預(yù)測(cè)

現(xiàn)在到了最后一步，我們將后驗(yàn)函數(shù)定義為對(duì)數(shù)正態(tài)分布，期望由物理模型定義：

已知損壞（在這種情況下為裂縫長(zhǎng)度）是傾斜的，因此對(duì)數(shù)正?；?Gumbel 分布通常用于模擬損傷。在 TFP 中表達(dá)模型如下所示：

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)

如下所示，使用混合物理概率模型的裂縫長(zhǎng)度的 95％ 不確定界限的預(yù)測(cè)平均值。很顯然，該模型不僅捕獲平均行為，還捕獲每個(gè)時(shí)間點(diǎn)模型預(yù)測(cè)的不確定性估算。當(dāng)我們偏離觀察時(shí)，模型預(yù)測(cè)中的不確定性呈指數(shù)增長(zhǎng)趨勢(shì)。

下一步

我們選擇這個(gè)例子十分小心，來(lái)對(duì)概率模型做一個(gè)介紹。將概率模型應(yīng)用于實(shí)際應(yīng)用程序存在一些挑戰(zhàn)。即使在像這里所展示的簡(jiǎn)單問(wèn)題中，如果我們放松高斯先驗(yàn)的假設(shè)，那么解決方案也會(huì)變得棘手。如果讀者感到好奇，可以嘗試使用統(tǒng)一先驗(yàn)來(lái)觀察一下模型預(yù)測(cè)能力急劇下降的情況。

另一個(gè)微妙點(diǎn)是定義可能性的方式。我們選擇直接在預(yù)測(cè)誤差上定義它，而不是預(yù)測(cè)值。請(qǐng)注意，下面在 “基于預(yù)測(cè)誤差的可能性” 中描述的各個(gè)藍(lán)線是模型預(yù)測(cè)誤差的樣本。直接將公式更改為模型預(yù)測(cè) - 而不是預(yù)測(cè)誤差 - 會(huì)產(chǎn)生非單調(diào)的非物理結(jié)果，如 “基于實(shí)際預(yù)測(cè)的可能性 ” 所示。我們稱這些結(jié)果是非物理的，因?yàn)榱芽p長(zhǎng)度不會(huì)隨著時(shí)間的推移而減小。如果我們只查看上面顯示的 “模糊” 百分位數(shù)視圖，我們可能不會(huì)注意到模型公式中可能導(dǎo)致較大預(yù)測(cè)誤差的細(xì)微差別。我們將在下一篇文章中解決其中一些更實(shí)際的挑戰(zhàn)及其緩解措施。

這是系列文章中的第一篇，旨在通過(guò) TFP 擴(kuò)展概率和深度學(xué)習(xí)技術(shù)在工業(yè)應(yīng)用中的應(yīng)用。我們（@sarunkarthi）很想聽(tīng)聽(tīng)關(guān)于您的應(yīng)用程序的說(shuō)明，并期待看到這些方法應(yīng)用到更多獨(dú)特的方方面面。請(qǐng)繼續(xù)關(guān)注此文章供稿，了解有關(guān)異常檢測(cè)，缺失數(shù)據(jù)估算和變異推斷預(yù)測(cè)的更多更新的示例。