基于改進(jìn)的動(dòng)態(tài)Kriging模型的結(jié)構(gòu)可靠度算法

魏娟， 張建國,*， 邱濤

(1. 北京航空航天大學(xué) 可靠性與系統(tǒng)工程學(xué)院, 北京 100083；2. 北京航空航天大學(xué) 可靠性與環(huán)境工程技術(shù)重點(diǎn)實(shí)驗(yàn)室, 北京 100083)

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

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

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

1 動(dòng)態(tài)Kriging模型

1.1 Kriging模型

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

(1)

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

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

i,j=1,2,…,k

(2)

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

(3)

(4)

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

基于給定的樣本點(diǎn)，多項(xiàng)式參數(shù)向量β與隨機(jī)過程方差σ2的估計(jì)值計(jì)算式分別為

(5)

(6)

式中：F為由回歸多項(xiàng)式函數(shù)值組成的回歸矩陣[18]。

(7)

u=FTR-1r0-f

(8)

(9)

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

1.2 動(dòng)態(tài)更新機(jī)制

為了提高可靠性分析的效率，動(dòng)態(tài)代理模型應(yīng)符合以下條件：①初始樣本數(shù)量?。虎谟糜谛颖緦W(xué)習(xí)的替代模型擬合；③每迭代一次，對FEA進(jìn)行一次調(diào)用，以修正代理模型。

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

步驟1建立初始的樣本集。通過拉丁超立方抽樣生成初始樣本點(diǎn)，并調(diào)用FEA計(jì)算其響應(yīng)值。

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

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

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

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

2 改進(jìn)的動(dòng)態(tài)Kriging模型

2.1 優(yōu)化算法

2.1.1 PSO算法

(10)

(11)

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

(12)

式中：C=c1+c2。

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

2.1.2 SA算法

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

(13)

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

設(shè)φj為每個(gè)變量接受的移動(dòng)次數(shù)與試驗(yàn)次數(shù)的比值，則擾動(dòng)可寫為

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

(14)

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

(15)

式中：Vi為V中的元素。

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

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

Ti+1=RTTi

(16)

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

2.1.3 混合PSOSA算法

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

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

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

步驟2粒子的位置?？尚杏蚓鶆蚍植糾個(gè)粒子；確定每個(gè)粒子的適應(yīng)度值。

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

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

步驟8重復(fù)整個(gè)過程直到收斂。

該算法主要有以下幾個(gè)特點(diǎn)：

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

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

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

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

2.2 Kriging建模流程

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

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

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

1) 似然函數(shù)作為適應(yīng)度函數(shù)。

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

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

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

3 案例分析

3.1 數(shù)值算例1

假設(shè)某結(jié)構(gòu)的極限狀態(tài)函數(shù)[12]為

G(x)=x1x2-1500

(19)

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

表1中，βr為可靠度系數(shù)，βr越大，可靠度越大。

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

將式(19)轉(zhuǎn)化為標(biāo)準(zhǔn)正態(tài)分布變量空間下的功能函數(shù)，即

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

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

(20)

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

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

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

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

3.2 數(shù)值算例2

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

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

3.3 工程案例

渦輪盤是現(xiàn)代航空發(fā)動(dòng)機(jī)的核心部件，它在發(fā)動(dòng)機(jī)燃燒室內(nèi)受到高溫燃?xì)獾耐苿?dòng)，將燃?xì)獾臒崮苻D(zhuǎn)化為機(jī)械能，驅(qū)動(dòng)發(fā)動(dòng)機(jī)的運(yùn)轉(zhuǎn)。其可靠性水平的高低直接影響發(fā)動(dòng)機(jī)的技術(shù)水平。

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

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

當(dāng)危險(xiǎn)點(diǎn)的最大應(yīng)力超過允許值Rmax=0.97 GPa，則結(jié)構(gòu)發(fā)生失效，相應(yīng)的功能函數(shù)可表示為

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

(21)

上述功能函數(shù)為隱式函數(shù)。

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

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

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

表3 渦輪盤參數(shù)Table 3 Parameters of turbine disk

圖5 應(yīng)力和應(yīng)變計(jì)算結(jié)果Fig.5 Calculation results of stress and strain

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

表4 計(jì)算結(jié)果對比Table 4 Comparison of calculation results

方法函數(shù)調(diào)用次數(shù)樣本點(diǎn)失效概率/10-3MCS方法1053.300本文算法15403.352

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

4 結(jié) 論

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

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