基于改進的動態(tài)Kriging模型的結構可靠度算法

魏娟， 張建國,*， 邱濤

(1. 北京航空航天大學 可靠性與系統工程學院, 北京 100083；2. 北京航空航天大學 可靠性與環(huán)境工程技術重點實驗室, 北京 100083)

當現有的可靠性分析方法應用于復雜的工程結構時，往往面臨巨大的挑戰(zhàn)[1]，功能函數通常是高度非線性甚至是隱式的，而且需要借助有限元分析(Finite Element Analysis，FEA)進行評估，計算量大，計算時間長[2-3]。一階可靠度算法(First Order Reliability Methodology，FORM)、二階可靠度算法(Second Order Reliability Methodology,SORM)已被廣泛用于估算結構系統失效的概率[4]，但是只能適用于功能函數表達式已知的情況。蒙特卡羅模擬(Monte Carlo Simulation,MCS)方法是結構可靠性分析中最簡單也是最廣泛使用的方法，計算精度高，適用于結構可靠性的隱式問題。但由于工程結構的失效概率通常較小，MCS須生成大量樣本以確保準確性，所帶來的計算成本可能會變得非常高，從而極大地限制了MCS在工程中的應用。

為了降低可靠性分析的計算成本，代理模型越來越多地被用來實現結構可靠性。代理模型作為一個近似模型可以找到輸入和輸出之間的潛在關系。如果與功能函數足夠接近，則代理模型可以完全表示功能函數，以準確評估性能函數的值。常用的方法有響應面方法(Response Surface Methodology，RSM)[5]、支持向量機[6]、神經網絡[7]、Kriging 模型[8]。其中，Kriging模型是一種方差估計最小的無偏估計的隨機過程算法，由于其對非線性函數的良好近似能力和獨特的誤差估計功能[9]，越來越多地用于可靠度的求解過程。Zhang等[10]運用Kriging模型進行結構可靠度分析，Liu等[11]運用Kriging模型和MCS方法進行齒輪接觸可靠性的分析。劉瞻[12]和陳志英[13]等將人工智能優(yōu)化算法融入到Kriging模型的構造中，進一步提高了預測精度。為了降低復雜結構的可靠性分析的計算成本并實現快速可靠性分析，盡可能減少樣本數量，Ma等[14]在采樣過程中應用學習函數，從而構造Kriging模型，Li等[15]提出一種基于最可能失效點(Most Probable Point,MPP)的局部抽樣方法，提高了Kriging的精度。Echard等[16]提出動態(tài)Kriging聯合MCS的可靠度計算方法。Xia等[17]提出了一種基于自適應動態(tài)泰勒Kriging的高效可靠性優(yōu)化算法。通過二元粒子群算法優(yōu)化選擇基函數，并確定具有期望擬合誤差的最小數量的采樣數據。

由于初始訓練的高質量樣本總是難以獲得，本文在訓練小樣本的基礎上，采用修改代理模型擬合誤差的動態(tài)更新機制。針對復雜結構隱式功能函數的可靠性分析問題，本文將混合粒子群-模擬退火(Particle Swarm Optimization-Simulated Annealing,PSOSA)算法應用到動態(tài)Kriging模型中相關參數的尋優(yōu)過程，并將改進的Kriging模型應用于可靠度的分析計算問題。該方法運用PSOSA算法高效求解最優(yōu)相關參數，構造Kriging模型；同時基于動態(tài)更新機制，逐步加入最佳樣本點，從而不斷修正 Kriging 模型的精度，并結合一次二階矩法進行可靠度計算。

1 動態(tài)Kriging模型

1.1 Kriging模型

近年來,Kriging模型作為一種新型的響應面模型技術在航天等工程優(yōu)化領域中得到了廣泛應用。對于k個n維樣本點及其對應的響應值Y=[G(x1),G(x2),…,G(xk)]T，Kriging模型采用高斯隨機過程模型，其插值結果定義為已知樣本函數響應值的線性加權[9]，即

(1)

式中：x=[x1,x2,…,xk]T；f(x)=[f1(x),f2(x),…,fp(x)]T為回歸多項式基函數向量，p為回歸多項式的數量；β=[β1,β2,…,βp]T為多項式參數向量；z(x)為服從正態(tài)分布N(0,σ2)的隨機過程，σ為標準差，協方差方程為

cov(z(xi),z(xj))=σ2R(xi,xj)

i,j=1,2,…,k

(2)

式中：R(xi,xj)為樣本點中任意2個輸入向量xi和xj的空間相關方程，通常采用高斯相關方程，其形式為

(3)

(4)

式中：R為相關矩陣，其中元素Rij=R(xi,xj)(i,j=1,2,…,k)，相關參數直接控制著Kriging模型的輸入輸出特性。

基于給定的樣本點，多項式參數向量β與隨機過程方差σ2的估計值計算式分別為

(5)

(6)

式中：F為由回歸多項式函數值組成的回歸矩陣[18]。

(7)

u=FTR-1r0-f

(8)

(9)

式中：r0=[R(x0,x1),R(x0,x2),…,R(x0,xk)]T。

1.2 動態(tài)更新機制

為了提高可靠性分析的效率，動態(tài)代理模型應符合以下條件：①初始樣本數量小；②用于小樣本學習的替代模型擬合；③每迭代一次，對FEA進行一次調用，以修正代理模型。

本文選擇動態(tài)Kriging模型來減少調用FEA的次數。代理模型的計算精度主要取決于對極限狀態(tài)面的逼近程度。一般來說，使用一些預設的訓練樣本很難構建具有高擬合精度的響應曲面，而響應曲面的動態(tài)更新可用于校正擬合誤差。動態(tài)Kriging可靠度具體算法流程如下：

步驟1建立初始的樣本集。通過拉丁超立方抽樣生成初始樣本點，并調用FEA計算其響應值。

步驟2構造動態(tài)Kriging。首先，通過優(yōu)化技術尋找最優(yōu)相關參數；然后，設置最優(yōu)基函數，并構造代理模型來近似目標函數。

步驟3利用DACE工具箱可求出模型對各個變量的梯度，結合FORM法計算MPP[12]。判斷MPP是否收斂。若是，求解可靠度等指標。若否，進入步驟4。

步驟4上述MPP作為新的樣本點加入初始樣本集，繼續(xù)迭代。

上述過程實現了可靠度計算過程中，Kriging模型在實際極限狀態(tài)響應面附近擬合誤差的自適應修正。減少了FEA的調用次數，這使計算成本大大降低，有利于減少模型對樣本的依賴程度，加快了收斂速度。

2 改進的動態(tài)Kriging模型

2.1 優(yōu)化算法

2.1.1 PSO算法

(10)

(11)

本文中系數w由式(12)確定：

(12)

式中：C=c1+c2。

PSO的工作原理如下：首先，粒子群由搜索區(qū)域內隨機分布的粒子組成，然后，PSO通過更新粒子的速度和位置來迭代尋找最優(yōu)解。當滿足收斂標準時，循環(huán)終止。但是，運用公式Pb和Pg來搜索最優(yōu)解，收斂速度快。然而所有粒子有移向最優(yōu)經驗點的趨勢，并在局部區(qū)域內搜索，所以可能會導致收斂過早的問題。因此，盡管PSO比其他方法收斂速度快，但結果可能精度不高。

2.1.2 SA算法

模擬退火算法起源于固體退火原理，模擬了高溫金屬降溫的熱力學過程。在搜索過程中具有概率突跳的能力，能夠有效地避免搜索過程中陷入局部極小解。試驗點從xi到xi+1的接受概率為

(13)

式中：y為目標函數，Δy=y(xi+1)-y(xi)和T為控制參數，可類比為溫度，溫度的降低實際上控制了接受概率。

設φj為每個變量接受的移動次數與試驗次數的比值，則擾動可寫為

Xtrial=Xcurrent+R′(-1,1)V

(14)

式中：R′(-1,1)為-1～1之間的隨機數；V為步長向量。Vbi是變量上下界之差的一半，可表示為

(15)

式中：Vi為V中的元素。

理論上證明，當采用非常緩慢的溫度下降率時，找到全局最小值的概率可以接近1。

最常見的溫度下降規(guī)律為

Ti+1=RTTi

(16)

式中：T0為初始溫度；RT為0.8～0.999 9之間的常數。

2.1.3 混合PSOSA算法

在混合PSOSA算法中，搜索最優(yōu)解時，速度矢量控制粒子在搜索空間內的移動方式。速度矢量有3部分組成：

本文PSOSA算法的主要思想是通過更新群體的群體最佳位置來改善群體的社會行為。當PSO循環(huán)中最優(yōu)解沒有改進時，舊的群體最佳位置將被替換為使用SA算法計算的新位置。新的最佳位置實際上向群體的社會領導者(即前一個全局最優(yōu)解所屬的粒子)發(fā)出信號，用于更新其方向。PSOSA算法充分結合了PSO算法全局搜索能力和SA算法的局部搜索能力，能達到很好的收斂結果。算法示意圖如圖1所示,具體算法如下：

步驟1初始化。m為粒子群粒子的數；imax為允許最大迭代次數；w為權重因子；c1和c2為學習因子；vmax為速度閾值，T0為初始溫度；RT為溫度降低參數。

步驟2粒子的位置。可行域均勻分布m個粒子；確定每個粒子的適應度值。

步驟4更新每個粒子的速度和位置：

步驟6判斷最優(yōu)解是否有所改進。如果第i+1次迭代的Pg沒有優(yōu)于第i次迭代的Pg，則進入步驟7；否則，更新溫度Ti+1=TiRT，返回步驟4。

步驟8重復整個過程直到收斂。

該算法主要有以下幾個特點：

1) 通過改進速度矢量的組成促使群體向整體最優(yōu)位置所在的搜索空間移動，這加快了收斂過程。

2) PSO算法和SA算法的結合提高了整個尋優(yōu)機制的質量和穩(wěn)定性。

3)在SA迭代中PSO速度公式中的認知成分得到了保留。

圖1 PSOSA算法示意圖Fig.1 Sketch map of PSOSA algorithm

2.2 Kriging建模流程

目前Kriging 建模過程普遍采用的模式搜索法對初始點十分敏感，容易造成求解θ*無法收斂，Kriging 預測精度很低。主要原因是模式搜索所采用的單點序列搜索方式決定了它只能從一個指定的起點出發(fā)，并收斂于一個距起點較近的局部最優(yōu)解。為了尋找最優(yōu)θ*，只能將起點接近于θ*，然而這很難確定。不僅如此，模式搜索法搜索路徑單一，不能根據似然函數的變化而變化，因此即使根據經驗信息，使起點接近θ*，搜索路徑也有可能偏離導致無法收斂[13]。

本文一方面運用PSOSA算法進行搜索，種群中各粒子在搜索最優(yōu)解的過程中實現信息共享，每個粒子不斷根據個體極值與群體極值的變化情況動態(tài)調整著自己的搜索方向，另外，當PSO循環(huán)沒有改善，原來全局最優(yōu)位置被通過SA計算得到的新位置代替，從而在缺乏先驗信息的前提下，也能搜索到最優(yōu)解θ*。另一方面,通過自適應動態(tài)更新機制減少調用FEA的次數，從而提高計算效率。

其中，運用PSOSA算法搜索極大似然意義下的最優(yōu)相關參數時：

1) 似然函數作為適應度函數。

2) 為轉化為無約束優(yōu)化，粒子采用對數形式。

綜上所述，本文的總體計算流程如圖2所示。

圖2 改進的動態(tài)Kriging模型流程圖Fig.2 Flowchart of improved dynamic Kriging model

3 案例分析

3.1 數值算例1

假設某結構的極限狀態(tài)函數[12]為

G(x)=x1x2-1500

(19)

式中：x1和x2的分布類型為正態(tài)分布，x1～N(38,3.8)，x2～N(54,2.7)，并且x1和x2相互獨立，具體分布參數見表1。本案例中，c1=1.5，c2=1.5，RT=0.8，w=2。

表1中，βr為可靠度系數，βr越大，可靠度越大。

對比MCS、RSM、經典Kriging和PSO-Kriging[13]方法，由表1結果可知，MCS方法選取樣本點108，計算失效概率結果為6.3×10-3；RSM方法處理非線性問題有一定的局限性，計算結果誤差較大；經典Kriging方法、PSO-Kriging方法與本文算法選取同樣數量的樣本數量，失效概率與MCS方法接近，但是經典Kriging方法與MCS方法相比，誤差為4.76%，PSO-Kriging方法誤差為4.13%，本文算法誤差為1.58%。由此看出本文算法精度相比于其他2種方法有所提高。

將式(19)轉化為標準正態(tài)分布變量空間下的功能函數，即

G(u1,u2)=Gx(σx1u1+ux1,σx2u2+ux2)=

(σx1u1+μx1)(σx2u2+ux2)-1 500

(20)

式中：u1和u2相互獨立，均服從標準正態(tài)分布。轉換后的極限狀態(tài)函數G(u1,u2)=0為一曲線，G(u1,u2)>0一側為安全域，G(u1,u2)<0一側為失效域，如圖3所示。根據本文所提算法，首先在給定隨機變量空間內，生成初始設計點，建立Kriging模型，并通過后續(xù)的動態(tài)更新機制逐漸增加樣本點。初始設計點和擬合過程的樣本點如圖4所示。由圖可知，通過更新機制可以有效選擇功能函數附近的樣本點，可以充分利用少量樣本點的信息擬合待求功能函數。

表1 不同方法結果對比(算例1)Table 1 Comparison of results of different methods (Example 1)

圖3 極限狀態(tài)函數Fig.3 Limit state function

圖4 算法說明圖Fig.4 Illustration of algorithm

3.2 數值算例2

將本文計算結果與其他方法進行對比，結果如表2所示。表2說明RSM在處理高度非線性問題時無法準確收斂和獲得正確結果。與其他方法相比，本文算法的結果與MCS結果最接近。其中MCS計算時間長達60.68 s，本文算法需時12.18 s。

表2 不同方法結果對比(算例2)Table 2 Comparison of results of different methods (Example 2)

3.3 工程案例

渦輪盤是現代航空發(fā)動機的核心部件，它在發(fā)動機燃燒室內受到高溫燃氣的推動，將燃氣的熱能轉化為機械能，驅動發(fā)動機的運轉。其可靠性水平的高低直接影響發(fā)動機的技術水平。

本文選擇某航空發(fā)動機渦輪盤為例，盤上銷釘沿周向均勻分布。參數不確定性的主要來源主要包括載荷、材料的隨機性。本文取轉速均值為9 550 r/min，銷軸的材料為3Cr13，輪盤的材料為TC11，葉片的載荷施加在銷軸上，合力均值為24 925 N。具體參數信息見表3。

渦輪盤的可靠性水平的高低主要取決于結構危險點處的應力應變分布。由于結構形狀、載荷為完全對稱的，本文選取輪盤的1/37為研究模型。首先以各變量的均值為設置參數，借助ANSYS 18.1 進行仿真計算，得到結果如圖5所示，輪盤與銷軸交界處應力應變水平最高，是結構的危險點，其中最大應力為0.834 GPa，最大應變?yōu)?.010 17。

當危險點的最大應力超過允許值Rmax=0.97 GPa，則結構發(fā)生失效，相應的功能函數可表示為

G=Rmax-S(ω,E1,ε1,ρ1,E2,ε2,ρ2,P)

(21)

上述功能函數為隱式函數。

采用本文所提算法進行可靠度計算，首先建立Kriging模型，其中相關參數θ*的迭代過程見圖6。

由圖6可以看出，隨著迭代次數的增加，本文所提算法在第15次迭代時基本尋得最優(yōu)相關參數θ*。

利用本文算法與ANSYS 18.1、 Isight進行聯合仿真計算，最終求得可靠度。為了驗證本文算法的有效性和工程實用性，將其與MCS方法進行對比，結果如表4所示。

表3 渦輪盤參數Table 3 Parameters of turbine disk

圖5 應力和應變計算結果Fig.5 Calculation results of stress and strain

圖6 θ*的迭代過程Fig.6 Iteration process of θ*

表4 計算結果對比Table 4 Comparison of calculation results

方法函數調用次數樣本點失效概率/10-3MCS方法1053.300本文算法15403.352

由表4可知，本文算法的結果與MCS結果十分接近，誤差僅為1.576%，驗證了本文算法的有效性以及準確性，并適用于工程實踐。

4 結 論

1) 將PSOSA算法引入Kriging的構造過程，PSO算法的全局搜索能力及SA算法跳出局部極小值的能力相結合，克服了經典Kriging過程的局限性和對初始樣本點的依賴，保證了極大似然意義下的最優(yōu)相關參數的求解精度，使尋優(yōu)過程更加高效和精確，從而有效保證了Kriging模型的預測精度。

2) 本文同時通過動態(tài)更新機制，盡可能地減少樣本點和調用有限元等數值計算的次數，能夠較好地解決功能函數隱式和高度非線性的問題，通過數值案例分析，相比于其他可靠度算法，能夠有效提高可靠度計算的效率和精度。通過工程案例分析，本文所提算法可應用于工程實際問題，尤其是功能函數隱式且復雜的問題。有一定的工程實用價值。