基于自適應Kriging代理模型的交叉熵重要抽樣法

史朝印，呂震宙，李璐祎，王燕萍

西北工業大學 航空學院，西安 710072

結構可靠性分析的主要任務在于計算結構失效概率[1]，迄今為止，研究者們已經提出了許多高效的可靠性計算方法。基于設計點的方法[2-4]是最為常見的一種可靠性分析方法，其中，一次二階矩方法[4]最具有代表性，它通過將極限狀態函數在設計點處進行泰勒展開，保留展開式一階項，從而實現對極限狀態函數的線性近似。但該方法需要計算極限狀態函數在設計點處的梯度，且此方法的準確度嚴重依賴于所研究問題的線性特征。對于具有復雜失效域的問題，其失效域往往是數量眾多的不連通域，且每個失效域相對于整個取值域占比較小，失效域的個數難以預先得到，設計點也難以預先求解。因此該方法無法解決復雜失效域、多設計點等問題。數字模擬方法中最基本的蒙特卡羅模擬(Monte Carlo Simulation, MCS)方法通過抽樣并判斷樣本所處狀態，計算結構失效概率。對于實際工程問題，失效概率通常較小(10-3量級及更小)，導致MCS方法往往需要大量的樣本(一般為102～104/Pf，Pf為結構真實失效概率)，才能獲得收斂解，計算效率低下。因此，研究者們提出了一系列改進數字模擬方法。其中，重要抽樣(Importance Sampling, IS)法[5-8]是一種應用廣泛的方差縮減方法，它通過構建重要抽樣密度函數，使所抽取的樣本盡可能多地落入失效域中，從而大幅縮小抽樣樣本池規模，提高失效概率的計算效率。

重要抽樣方法的關鍵在于重要抽樣密度函數的構建。Melcher[9]提出了基于設計點的方法，通過將原隨機變量概率密度函數的中心平移至設計點處構建重要抽樣密度函數。但該方法無法處理多設計點、復雜失效域的問題。吳斌等[10]提出了擴大方差的方法，這種方法難以確定方差擴大的倍數。Priebe和Marchette[11]提出了混合重要抽樣法，對每個失效模式單獨構造重要抽樣密度函數后加權混合，有效解決了多失效模式問題。Au和Beek[12]提出了基于核密度估計的方法，通過馬爾可夫鏈對失效域進行預抽樣，獲得一定量的失效樣本后，通過核密度估計得到重要抽樣密度函數；但核密度估計方法計算成本較大，且精度難以保證。Dubourg等[13]提出了元模型重要抽樣(Meta-IS)方法，使用Kriging模型近似理論上最優重要抽樣密度函數，同時，將失效概率轉化為失效概率估計值和修正系數乘積的形式；但該方法在Kriging構建時收斂標準的物理意義并不明確，且需要大量調用極限狀態函數以計算樣本的真實響應值。Cadini等[14]提出了改進的Meta-IS-AK2方法，但該方法通過單高斯模型獲得的重要抽樣密度函數針對多失效域問題效率較低。Kurtz和Song[15]提出了交叉熵重要抽樣(Cross Entropy Importance Sampling, CE-IS)方法，使用高斯混合模型作為重要抽樣函數，通過交叉熵原理，迭代更新高斯混合模型(Gaussian Mixture Model, GMM)的參數，逐步逼近最優重要抽樣密度函數。該方法較Meta-IS-AK2獲得的重要抽樣密度函數有更高的抽樣效率，但在GMM的參數計算過程中，需大量調用極限狀態函數計算樣本點處真實輸出響應值，使得最終的計算效率大幅下降。且Kurtz和Song[15]所提方法以落入失效域的重要樣本比率作為收斂準則，沒有考慮失效概率估計值的變異系數，對于部分問題，會存在解不穩定，迭代次數過多，甚至無收斂解的問題。

對于實際工程問題，極限狀態函數往往是基于有限元模型[16]的隱式函數，每次調用都會帶來大量的計算量，因此，應盡可能減少極限狀態函數的調用次數。代理模型[17-22]可以有效避免極限狀態函數多次調用問題，其通過一定的選點準則，篩選出對失效概率計算影響較大的點，添加至模型訓練集，構建代理模型預測樣本點處的響應值，計算結構失效概率。其中，Kriging模型由于其在可靠性分析領域存在特有的優勢受到了最為廣泛的關注。Kriging模型假設樣本點之間服從高斯隨機過程，通過極大似然估計模型參數值。由于高斯過程的特性，Kriging模型不僅可以預測樣本點處的響應值，還可以預測樣本點處的方差值。因此，Kriging模型有別于其他代理模型，對于樣本點的預測具有隨機屬性，可以提供樣本點的統計信息，與可靠性分析通過輸入變量統計信息求解輸出響應部分統計信息的思想相吻合。基于此，研究者們提出了一系列的學習函數，自適應地篩選出訓練樣本點，以提高Kriging模型的構建效率。常用的自適應學習函數有期望可行性函數[23]、U學習函數[24]、H學習函數[25]、期望最大學習函數[26]等。其中，U學習函數綜合考慮樣本點距極限狀態面的距離和樣本點的預測誤差，篩選出最易錯誤判斷失效狀態的點，簡單易行，本文選擇U函數作為學習函數。通過自適應Kriging(Adaptive Kriging, AK)代理模型方法，來用少量訓練樣本即可構建高精度的收斂的代理模型替代極限狀態函數，顯著提高失效概率的計算效率。常用的AK-MCS[24]依舊需要構建大規模的樣本池，對于小失效概率問題，計算效率依舊低下。

本文將CE-IS方法與AK方法相結合，通過少量訓練樣本，自適應地擬合極限狀態面，減少CE-IS的極限狀態函數調用次數，同時通過GMM，構建高效的重要抽樣函數，大幅縮小AK-MCS方法的樣本池規模。相較于現有方法，本文所提方法有以下創新性：

1) 改進CE-IS方法收斂條件，避免冗余迭代，同時增強CE-IS方法的適用性。

2) 通過AK模型擬合極限狀態函數，可顯著減少CE-IS方法的計算量。

3) 通過GMM近似最優重要抽樣密度函數，相較于現有的重要抽樣密度函數的構造方法，更為高效。

1 基于交叉熵的重要抽樣法

1.1 重要抽樣法

結構可靠性分析的關鍵是失效概率的計算，對于已知結構，其極限狀態函數為g(x)，其中，x為隨機輸入變量，其概率密度函數為fX(x)，結構的失效域為F={x|g(x)≤0}，失效概率Pf即為輸入隨機變量落入失效域的概率，可由積分求得，即

(1)

式(1)通常無解析解，故常采用MCS方法作為參照解。引入失效域指示函數IF(x):

(2)

則通過MCS方法估計結構失效概率的表達式為

(3)

式中：xi(i=1,2,…,N)為N個按fX(x)抽取的MCS樣本點;Rd為d維的實數域。

(4)

(5)

理論上可以推導最優的重要抽樣密度函數hopt(x)為

hopt(x)=IF(x)·fX(x)/Pf

(6)

以hopt(x)作為重要抽樣密度函數時，失效概率估計值的方差為零。通過對最優重要抽樣密度函數進行近似擬合，即可得到抽樣效率最高的重要抽樣密度函數。

1.2 高斯混合模型

假設隨機向量x服從K元高斯混合分布，其概率密度函數p(x)為

(7)

式中：K為高斯混合成分個數；f(x|μk,Σk)為第k個混合成分的混合高斯分布密度函數，μk和Σk分別為相應均值向量和協方差矩陣;πk為混為混合成分系數，由密度函數的歸一化條件知，πk應滿足:

(8)

0≤πk≤1

(9)

式中:混合系數πk表示隨機變量x來自于第k個混合成分的概率。

記二值隨機變量zk(k=1,2,…,K)表示在給定樣本x的情況下，其是否由第k個混合變量生成，zk的后驗概率為

(10)

對于一個未知K元混合高斯分布模型

(11)

其未知參數v共3×K組，每一元高斯分布的參數為{πk,μk,Σk}，則

v={πk,μk,Σk；k=1,2,…,K}

(12)

通過交叉熵準則，可以迭代求出K元混合高斯分布的參數，使其盡可能地逼近理論最優重要抽樣密度函數。

1.3 交叉熵準則

在信息論中，Kullback-Leibler交叉熵(K-L Cross-Entropy，KL-CE)D,簡稱交叉熵(Cross-Entropy，CE)常用來度量函數的差異，其定義為

(13)

由式(13)可知，交叉熵越小，函數f(x)與函數g(x)之間的差異就越小。為了使p(x;v)逼近最優重要抽樣密度函數hopt(x)，計算它們之間的交叉熵:

lnp(x;v)dx

(14)

將式(6)代入式(14)可得

(15)

arg maxEX[IF(x)·ln(p(x;v))]

(16)

式(16)為嵌套形式，通常無解析解，一般采用迭代算法求解，逐步更新GMM參數，直至滿足收斂要求。

為了提高計算效率，每一步迭代也采用重要抽樣法進行計算，即，使用前一步迭代所得GMMp(x;w)作為當前重要抽樣密度函數，w為其參數。式(16)可改寫為

p(x;w)dx=arg maxEX～p(x;w)[IF(x)·

ln(p(x;v))·W(x)]

(17)

(18)

(19)

將GMM的定義式(7)代入式(19)，得到關于模型參數v的方程：

(20)

分別對3×K組模型參數求導，以均值向量μk(k=1,2,…,K)為例，可建立方程為

(21)

解得

(22)

同理可得Σk和πk(k=1,2,…,K)的表達式為

(23)

(24)

式(22)～式(24)即為基于交叉熵原理得到的GMM參數更新準則，通過多次迭代，GMM即可較為準確地近似最優重要抽樣密度函數hopt(x)。但注意到，式(22)～式(24)中失效域指示函數IF(x)的計算需要調用結構極限狀態函數，而多次迭代會使得計算量顯著加大，因此，本文將AK方法與CE-IS方法相結合，從而提高算法的效率。

2 CE-IS-AK算法

2.1 Kriging模型

Y=f(x)Tβ+z(x)

(25)

式中：f(x)=[f1(x)f2(x) …fM(x)]T為M×1維基函數向量;β=[β1β2…βM]T為相應的系數向量；z(x)為零均值高斯隨機過程，協方差函數為

(26)

(27)

(28)

R=

(29)

由上述過程可知，Kriging模型的構建即為相關函數參數θ的確定，采取最大似然估計計算最優參數為

(30)

(31)

[FR-1r(x)-f(x)]T(FTR-1F)-1·

[FR-1r(x)-f(x)]}

(32)

2.2 自適應Kriging模型

為了高效地選擇模型訓練點，通常需采取自適應的方法，選擇恰當學習函數與收斂條件。本文選擇的U學習函數表達式為

(33)

U函數表征了樣本輸出響應y=g(x)的正負狀態被誤分的概率。U函數值越小，樣本點被錯誤分類的概率越大，故而篩選出U函數值較小的樣本點作為新的訓練點：

(34)

式中：S為候選樣本池。

U(x)=2時，樣本點被正確分類的概率為1-Φ(2)=97.7%，故可以認為當U(x)≥2時，所構建的Kriging模型對樣本點的正負號判斷有超過97.7%的概率預測正確。因此，確定Kriging模型構建的收斂條件為

(35)

2.3 本文所提算法

本文所提方法結合了CE-IS方法與AK方法，避免了CE-IS方法計算量過大的缺陷，并提出了新的收斂準則，避免了冗余迭代，增強了CE-IS算法的適用性。同時，較AK-MCS方法縮減了樣本池，使小失效概率的計算更為高效。以下給出本文所提方法的詳細計算步驟。

4) 用Kriging模型計算樣本池S(l)中樣本點U函數值。

(36)

7) 計算抽樣效率：

(37)

判斷η是否大于指定值η*，如滿足，則執行第8)步；否則，根據式(22)～式(24)更新GMM參數，增加迭代次數l=l+1，轉跳至第3)步。

8) 擴大樣本池規模，如：N=10N，轉跳至第3)步。

方法流程圖如圖1所示。

圖1 CE-IS-AK方法流程圖Fig.1 Flow chart of CE-IS-AK method

2.4 一些討論及說明

文獻[15]以其所提ρ分位數作為收斂準則，ρ分位數的概念與抽樣效率η概念相同，均表征重要樣本落入失效域中的比例；以其作為收斂準則，即以抽樣效率η作為收斂準則，存在以下3個主要問題：

1) 對于部分復雜失效域問題，由于理論最優重要抽樣密度函數的高度非線性以及不連續性，GMM難以十分準確地逼近理論最優重要抽樣密度函數，單純以抽樣效率作為收斂準會造成方法迭代次數過多，甚至不收斂。

2) 對于部分失效概率較大的問題，只需較少的樣本點即可得到穩定的解，以抽樣效率作為收斂指標會使算法進行多次無意義迭代。

3) 對于部分極小失效概率問題，即使使用重要抽樣的方法，也需要較大規模的樣本池，僅通過抽樣效率難以判斷是否得到了穩定的收斂解。

本文以失效概率估計值的變異系數作為方法收斂準則，可以準確衡量解的穩定性，以及樣本池規模是否恰當，同時避免方法出現冗余迭代。

本文構建AK模型用于擬合極限狀態函數，因此，每次更新重要樣本池時，無需重新構建模型，僅需在原有Kriging模型的基礎上，通過學習函數添加少量訓練樣本點，修正現有模型，即可完成收斂。

根據文獻[15]的研究，本文建議高斯混合模型中混合成分數目K的取值應滿足如下要求：

max{d，Ncomponent}≤K≤5%ηNNcomponent

(38)

式中：Ncomponent為結構失效模式數目，對于單模式問題Ncomponent=1。建議迭代樣本池規模N為103～105。

對于復雜失效問題，往往難以預先判斷其失效域個數，因此建議選擇d

3 算例分析

本節給出5個算例證明所提算法的高效與適用性，第1個算例為小失效概率問題；第2個算例為串聯系統問題，具有多失效域；第3個算例為復雜多失效域問題，不連續失效域較多，失效邊界復雜；第4個算例為一個動態非線性振子問題；第5個算例為工程有限元算例。本文將所提方法與AK-MCS和CE-IS等方法所得結果進行對比。

3.1 算例1

一個非線性單設計點問題的極限狀態函數為

g(x1，x2)=0.5(x1-2)2-1.5(x2-5)3-3

(39)

輸入隨機變量為標準正態變量，圖2給出了算例一通過CE-IS-AK方法最終所得的重要樣本點和極限狀態面以及真實極限狀態面。表1給出了算例一CE-IS-AK方法具體迭代結果，表2給出了算例一的各方法計算對比結果，表中，Ncall為極限狀態函數調用次數，Ncondidate為備選樣本點總數。

針對此算例，本文選擇混合成分的個數K=20，每次迭代的樣本池規模為N=103，收斂條件為ε*=5%，η*=80%。從算例1結果可看出，本文所提方法具有較高的效率，這是由于GMM模型本身的特性，可以對不連續的多失效域進行近似。通過3次迭代，方法即收斂。相較于CE-IS方法，本文所提方法在大幅減小極限狀態函數調用次數的情況下，依舊保持了與其相當的高抽樣效率；與AK-MCS方法相比，本文在極限狀態函數調用次數基本相當的情況下，大幅度縮減了樣本池規模，在實際計算中，大幅提升了計算效率。

圖2 算例1 CE-IS-AK方法所得重要樣本點及 極限狀態響應面Fig.2 Approximate LSF and important samples obtained by CE-IS-AK of Example 1

表1 算例1 CE-IS-AK方法迭代結果

Table 1Results of iteration of CE-IS-AK of Example 1

迭代次數NcallNcandidateη(l)/%P^(l)f/10-5Cov[P^(l)f]/%l=012+51 000320.4532.09l=161 000672.1713.75l=241 000882.852.86

表2 算例1的計算結果Table 2 Results of Example 1

3.2 算例2

串聯系統具有4個失效域，其極限狀態函數為

g(x1,x2)=

(40)

式中：輸入隨機變量均為標準正態變量。

圖3給出了算例2通過CE-IS-AK方法最終所得的重要樣本點、極限狀態面以及真實極限狀態面。表3給出了算例2 CE-IS-AK方法具體迭代結果，表4給出了各方法計算對比結果。

針對此算例，本文選擇混合成分的個數K=50，每次迭代的樣本池規模為N=103，ε*=5%，迭代門限值為η*=80%。

CE-IS方法只考慮抽樣效率，故出現了冗余迭代的情況。而由于本文提出了更加合理的收斂條件，在失效概率估計值的變異系數滿足要求后即收斂，雖然最終得到的抽樣效率略低于CE-IS方法，但求解精度并無顯著區別，且計算量有所減小，若再進行一次迭代即可得到與CE-IS方法相近的求解情況。與AK-MCS方法相比，本文所提方法雖然調用極限狀態函數次數有所增加，但樣本池規模遠遠小于AK-MCS方法，避免了對非常大一部分安全域內的樣本點進行無意義的Kriging模型的反復多次預測。

圖3 算例2 CE-IS-AK方法所得重要樣本點及 極限狀態響應面Fig.3 Approximate LSF and important samples obtained by CE-IS-AK of Example 2

表3 算例2 CE-IS-AK方法迭代結果

Table 3 Results of iteration of CE-IS-AK of Example 2

迭代次數NcallNcandidateη(l)/%^Pf/10-3Cov[^P(l)f]/%l=030+351 000541.7132.09l=1421 000672.226.75l=2371 000762.223.66l=3181 000882.222.94

表4 算例2的計算結果Table 4 Results of Example 2

注：(3)*表示迭代到l=2，(4)*表示迭代到l=3。

3.3 算例3

一經典的復雜高度非線性多失效域問題的極限狀態函數為

(41)

式中：輸入隨機變量均為標準正態變量。表5給出了算例3的各種方法的計算結果。圖4給出通過所提方法最終所得的重要樣本點、近似極限狀態面以及真實極限狀態面。

針對算例3，本文選擇混合成分的個數K=100，每次迭代的樣本池規模為N=103，ε*=5%，η*=80%。此算例為一高度非線性多失效域算例，失效域數量較多，面積均較小且較為分散，理論上難以獲得抽樣效率較高的重要抽樣密度函數，但此算例失效概率較大，因此并不需要特別大量的重要樣本點即可得到收斂解。CE-IS方法由于僅以抽樣效率作為收斂指標，故沒有得到收斂的解。本文所提方法雖然最終抽樣效率較之前算例有明顯下降，但依舊得到了較為準確的收斂解，同時，抽樣效率仍為所有對比方法中最高且樣本池規模最小。故而，相對于其他對比方法而言，仍然具有較大優勢。

表5 算例3的計算結果Table 5 Results of Example 3

圖4 算例3 CE-IS-AK方法所得重要樣本點及 極限狀態響應面Fig.4 Approximate LSF and important samples obtained by CE-IS-AK of Example 3

3.4 算例4

一單自由度無阻尼彈簧振子模型，系統結構示意圖見圖5，系統的極限狀態函數為

(42)

針對此算例，選擇混合成分的個數K=100，每次迭代的樣本池規模為N=103，ε*=5%，η*=80%。此算例的各方法對比結果證明了本文所提算法對較高維數的問題依然可以高效處理。

圖5 單自由度無阻尼彈簧振子模型Fig.5 Model of single degree of freedom undamped spring oscillator

表6 算例4輸入隨機變量分布參數

Table 6Distribution parameters of input random variables of Example 4

隨機變量分布類型均值標準差m正態分布10.05c1正態分布10.1c2正態分布0.10.01r正態分布0.50.05F1正態分布10.2t1正態分布10.2

表7 算例4的計算結果Table 7 Results of Example 4

3.5 算例5

以桁架中點位移u(X)為研究對象，結構極限狀態函數為：

g(X)=τ-u(X)

(43)

式中:τ=12 cm為位移u(X)的閾值。

相應地，結構失效域為F={X|g(X)≤0}。算例5無法通過結構力學知識解析得到u(X)的具體表達式，因此，采用有限元仿真對其進行計算。

圖6 桁架結構示意圖Fig.6 A sketch of truss

表8 算例5輸入隨機變量分布參數

Table 8Distribution parameters of input random variables of Example 5

隨機變量分布類型均值標準差E1,E2/Pa對數正態2.1×10112.1×1010A1/m2對數正態2×10-32×10-4A2/m2對數正態1×10-31×10-4P1,P2,…,P6/N耿貝爾5×1047.5×103

本文采用規模為1×106的MCS樣本池，雖然在計算過程中需大量調用有限元模型，但由于該結構較為簡單，有限元模型計算迅速，因此，仍可通過MCS方法得到參照解。表9中給出了該工程算例的計算結果。同時，各方法均由同一計算機進行計算驗證，其計算時間T一并列出。

表9 算例5的計算結果Table 9 Results of Example 5

針對此算例，選擇混合成分的個數K=100，每次迭代的樣本池規模為N=103，ε*=5%，η*=80%。

3.6 一些討論及說明

本文所提算法可高效處理復雜失效域與小失效概率耦合問題，即結構系統的失效域較復雜且同時失效概率也較小的問題。

對于小失效概率問題，常規的數字模擬方法，需要規模巨大的樣本池來保證計算的精度。所提方法使用基于高斯模型的重要抽樣法，大幅度縮減了樣本池規模。同時，引入自適應Kriging代理模型輔助計算，減少了功能函數的調用次數，進一步提高了計算效率。算例1所示結果也證明了本文所提方法能在精度范圍要求內高效處理小失效概率問題。

對于復雜失效域問題，往往是串并聯系統問題、隱式函數問題，甚至是有限元模型問題，針對此類問題，常規代理模型方法雖然能減少計算成本巨大的功能函數的調用次數，但由于樣本池規模所限，計算時間成本依舊較大。本文所提方法引入基于高斯混合模型的重要抽樣法，極大地縮減了樣本池的規模，使計算效率進一步提高。算例2、算例3、算例4和算例5的計算結果也驗證了本文所提方法的高效性。

4 結 論

1) 本文針對復雜失效域和小失效概率耦合的失效概率計算問題，提出了交叉熵重要抽樣結合自適應Kriging模型的CE-IS-AK方法。通過交叉熵準則指導高斯混合模型逐步逼近理論最優重要抽樣密度函數，參數更新中使用AK模型近似計算樣本響應值，從而大大減少了計算過程中的結構極限狀態函數調用次數，在保證方法計算精度的同時提高了CE-IS方法的計算效率。

2) 本文改進了CE-IS方法的收斂準則，以失效概率估計值的變異系數作為最終指標，抽樣效率作為參考指標，避免了冗余迭代，也克服了CE-IS方法對于部分復雜失效域問題無法得到收斂解的缺點，擴大了方法的適用范圍。

3) 相較于CE-IS方法，由于AK模型的引入，結構極限狀態函數的調用次數大幅縮減，計算效率顯著提高。相較于AK-MCS方法，由于采用了高斯混合模型代替重要抽樣密度函數，抽樣效率較MCS方法得到了顯著提高，樣本池規模大幅縮小，適用于小失效概率問題。同時，由于利用了高斯混合模型，方法對于多失效域和復雜失效域等問題也有良好的適用性。最后，算例分析結果也證明了本文所提方法的高效性。