基于二參數(shù)尋優(yōu)設(shè)計(jì)點(diǎn)的混合結(jié)構(gòu)可靠性分析算法

邱濤， 張建國(guó)， 邱繼偉， 魏娟， 游令非

(1.北京航空航天大學(xué) 可靠性與系統(tǒng)工程學(xué)院, 北京 100191;2.北京航空航天大學(xué) 可靠性與環(huán)境工程技術(shù)國(guó)防科技重點(diǎn)實(shí)驗(yàn)室, 北京 100191;3.中國(guó)兵器工業(yè)標(biāo)準(zhǔn)化研究所, 北京 100089)

0 引言

在機(jī)械產(chǎn)品可靠性分析過(guò)程中，受加工尺寸、材料性能等條件的限制，往往需要考慮大量的不確定性。對(duì)于部分統(tǒng)計(jì)信息不足的不確定變量，難以得到對(duì)應(yīng)的精確分布，可以采用區(qū)間變量對(duì)其進(jìn)行表述。因此，研究概率-區(qū)間共存的可靠性問(wèn)題具有一定的工程價(jià)值[1]。

針對(duì)概率-區(qū)間混合的可靠性問(wèn)題，Kang等[2]采用單層優(yōu)化模型將區(qū)間變量嵌入尋優(yōu)設(shè)計(jì)點(diǎn)(簡(jiǎn)稱MPP)的優(yōu)化問(wèn)題中。王騫等[3]通過(guò)構(gòu)建響應(yīng)面的方法求解最大失效概率。Xie等[4]提出一種單循環(huán)算法，將區(qū)間優(yōu)化的Karush-Kuhn-Tucker(KKT)最優(yōu)條件作為MPP搜索的約束條件。Du[5]提出了基于一次二階矩的統(tǒng)一不確定分析算法。姜潮等[6]考慮極限狀態(tài)帶的兩種不同構(gòu)形，采用了基于響應(yīng)面的混合可靠性算法。Elishakoff等[7]研究了不確定性的概率模型和凸模型混合問(wèn)題。Xie等[8]用高維模型表示功能函數(shù)，將混合可靠性問(wèn)題轉(zhuǎn)化為單步可靠性計(jì)算模型。宋向華等[9]利用非概率理論解決了混合可靠性問(wèn)題。劉瞻等[10]基于Kriging模型和重要抽樣法進(jìn)行了混合模型的可靠性分析。Du[11]提出了一種序列迭代算法，通過(guò)區(qū)間分析和概率分析的反復(fù)迭代求解，使隨機(jī)變量和區(qū)間變量同時(shí)收斂至最優(yōu)解。Zhang等[12]利用積分方法求解了概率-區(qū)間混合的失效概率。

目前，針對(duì)極限狀態(tài)函數(shù)非線性程度較高的可靠性問(wèn)題，計(jì)算結(jié)果的收斂穩(wěn)定性與效率是該問(wèn)題的關(guān)鍵。通過(guò)解耦雙層優(yōu)化模型，對(duì)外層進(jìn)行概率分析，對(duì)內(nèi)層進(jìn)行區(qū)間分析實(shí)現(xiàn)高效序列迭代。雖然一些概率分析算法(如Hasofer-Lind和Rackwitz-Fiessler(HLRF)算法等[13-16])能夠直接進(jìn)行可靠度計(jì)算，但存在結(jié)果不收斂的問(wèn)題。在概率分析時(shí)，通過(guò)引入兩個(gè)調(diào)整參數(shù)η和λ，分別調(diào)整搜索步長(zhǎng)和搜尋方向，改進(jìn)了傳統(tǒng)的尋優(yōu)MPP算法；在區(qū)間分析時(shí)，將功能函數(shù)進(jìn)行2次泰勒近似，并轉(zhuǎn)化為易于求解的2次規(guī)劃問(wèn)題，為解決概率-區(qū)間混合的可靠性問(wèn)題提供了一種新的算法。最后，通過(guò)兩個(gè)案例驗(yàn)證了該算法的精度、效率及穩(wěn)定性。

1 概率- 區(qū)間混合可靠性雙層優(yōu)化模型

機(jī)械結(jié)構(gòu)中既存在隨機(jī)變量又存在區(qū)間變量，則功能函數(shù)可以表示為

Z=g(X,Y)，

(1)

式中：X=[X1,X2,…,Xn]T為n維隨機(jī)向量；Y=[Y1,Y2,…,Ym]T為m維區(qū)間向量。當(dāng)構(gòu)件發(fā)生故障，即g(X,Y)≤0時(shí)，失效概率Pf定義為

Pf=Pr{g(X,Y)≤0}.

(2)

如圖1所示，X1、X2為隨機(jī)變量，功能函數(shù)因含區(qū)間變量Y，g(X,Y)=0在隨機(jī)空間中不再是一個(gè)曲面，而是由兩個(gè)邊界面構(gòu)成的極限狀態(tài)帶?？梢酝ㄟ^(guò)(3)式和(4)式來(lái)定義失效概率的范圍：

(3)

(4)

圖1 極限狀態(tài)區(qū)域Fig.1 Limit state area

通過(guò)將隨機(jī)變量轉(zhuǎn)換為獨(dú)立的標(biāo)準(zhǔn)正態(tài)隨機(jī)變量，上述問(wèn)題可以轉(zhuǎn)換為求解如下兩個(gè)優(yōu)化問(wèn)題，最大和最小可靠度指標(biāo)分別為

(5)

(6)

式中：U為隨機(jī)變量轉(zhuǎn)換到標(biāo)準(zhǔn)正態(tài)空間內(nèi)的值；G(U,Y)為隨機(jī)變量轉(zhuǎn)換到標(biāo)準(zhǔn)正態(tài)空間下功能函數(shù)的值；βU為可靠度指標(biāo)最大值，βL為可靠度指標(biāo)最小值，對(duì)應(yīng)的最大和最小失效概率分別為

Pfmax=Φ(-βL)，

(7)

Pfmin=Φ(-βU).

(8)

失效概率的范圍為Pfmin≤Pf≤Pfmax，但在實(shí)際工程問(wèn)題中，人們通常關(guān)心最大失效概率。因此本文所提算法將以求解最大失效概率Pfmax為例進(jìn)行混合可靠性分析。

2 混合可靠性分析

針對(duì)第1節(jié)(6)式中的概率-區(qū)間混合可靠性模型，提出一種基于二參數(shù)尋優(yōu)MPP的概率-區(qū)間混合可靠性分析算法。如圖2所示，該算法分為概率分析和區(qū)間優(yōu)化兩部分，二者的計(jì)算效率共同決定著混合可靠性的計(jì)算成本。在隨機(jī)變量迭代中，采用基于二參數(shù)尋優(yōu)MPP的概率分析迭代算法；在區(qū)間變量迭代中，首先將區(qū)間優(yōu)化問(wèn)題轉(zhuǎn)換為2次規(guī)劃問(wèn)題，再結(jié)合梯度投影法進(jìn)行區(qū)間優(yōu)化分析，求解區(qū)間最優(yōu)點(diǎn)。

圖2 序列迭代計(jì)算流程Fig.2 Calculation process of sequence iteration

為提高雙層優(yōu)化模型的計(jì)算效率，將雙層優(yōu)化模型解耦為概率分析和區(qū)間分析的序列迭代可靠性模型。一旦滿足收斂條件，則得到最大失效概率點(diǎn)U*，最大失效概率為

(9)

2.1 概率分析

2.1.1 一次二階矩算法和改進(jìn)的一次二階矩算法

Hasofer等[13]和Rackwitz等[14]提出了HLRF一次二階矩算法，其尋優(yōu)設(shè)計(jì)點(diǎn)為

Uk+1=Uk+dk，

(10)

式中：Uk為第k次迭代的MPP；Uk+1為第k+1次迭代的MPP；dk為搜索方向，

(11)

當(dāng)經(jīng)過(guò)有限次迭代結(jié)果收斂后，則收斂點(diǎn)即為MPP. 雖然上述HLRF算法迭代效率較高，但當(dāng)功能函數(shù)非線性程度較高時(shí)，此算法不能保證收斂。因此，基于不確定一維搜索準(zhǔn)則的改進(jìn)的一次二階矩(iHLRF)算法[17]被提出來(lái)，其迭代搜索設(shè)計(jì)點(diǎn)為

Uk+1=Uk+λkdk.

(12)

根據(jù)Armijo準(zhǔn)則：

(13)

(14)

ck為懲罰參數(shù)，且ck>0，mk(U)梯度可表示為

(15)

2.1.2 基于二參數(shù)尋優(yōu)MPP的概率分析

設(shè)當(dāng)前迭代步驟為第k+1次，在隨機(jī)變量迭代過(guò)程中，區(qū)間變量保持不變，故在以下尋優(yōu)MPP算法中省略區(qū)間變量。

(16)

圖3 功能函數(shù)及其在Uk點(diǎn)處的切平面Fig.3 Performance function and its tangent plane at point Uk

(17)

圖4 在平面G=0上尋優(yōu)MPPFig.4 Optimizing MPP on plane G=0

(18)

圖5 選定η確定平面G=ηG(Uk)Fig.5 Selecting η to determine the plane G=ηG(Uk)

圖6 在G=ηG(Uk)平面上尋優(yōu)MPPFig.6 Optimizing MPP on the plane G=ηG(Uk)

(19)

在平面G=ηG(Uk)上，原點(diǎn)連接到輔助點(diǎn)Hk+1的方向即為Uk+1的搜索方向，Uk+1的方向余弦矢量αk+1為

(20)

Uk+1=βk+1αk+1，

(21)

式中：βk+1為第k+1次迭代的可靠度指標(biāo)。

將(21)式代入(18)式，得

(22)

求解(22)式，得可靠度系數(shù)βk+1為

(23)

根據(jù)得到的可靠度系數(shù)βk+1，可由(21)式得第k+1次迭代的MPPUk+1.

通過(guò)引入兩個(gè)調(diào)整參數(shù)分別控制搜索步長(zhǎng)和搜索方向，通過(guò)參數(shù)η調(diào)整搜索步長(zhǎng)，使功能函數(shù)為高度非線性時(shí)，尋優(yōu)MPP收斂成為可能；通過(guò)參數(shù)λ調(diào)整搜索方向，來(lái)提高尋優(yōu)MPP的效率。在選擇λ時(shí)，可以將其設(shè)定為較大的值，一般取λ最大值為50，在具體的迭代中，可以設(shè)定λ=λ/c以自適應(yīng)更新迭代步長(zhǎng)，c為步長(zhǎng)的調(diào)整系數(shù)，是一個(gè)大于1的常數(shù)，根據(jù)數(shù)值分析，c最好取在1.1～1.3之間。如果λ的初始值較小，則不需要在整個(gè)迭代過(guò)程中改變?chǔ)?

2.2 區(qū)間分析

在第k+1次迭代中，利用概率分析求得Uk+1. 根據(jù)(6)式表示的雙層優(yōu)化模型，固定隨機(jī)變量Uk+1，區(qū)間分析可表示為(24)式的優(yōu)化問(wèn)題：

(24)

式中：Yk+1為第k+1次混合迭代的區(qū)間最優(yōu)點(diǎn)；Yk為第k次混合迭代的區(qū)間最優(yōu)點(diǎn)；YU和YL分別代表區(qū)間的上、下界矢量。

為提高計(jì)算效率，在進(jìn)行區(qū)間優(yōu)化之前，首先檢查優(yōu)化問(wèn)題的KKT條件，KKT優(yōu)化條件見(jiàn)(25)式和(26)式。若滿足KKT條件，則Yk+1=Yk；否則進(jìn)行區(qū)間優(yōu)化分析[18]。

對(duì)于(24)式中的優(yōu)化問(wèn)題，KKT條件為

(25)

(26)

式中：j=1,2,…,m表示m維區(qū)間變量的第j維。

在區(qū)間變量循環(huán)迭代中，固定隨機(jī)變量Uk+1，故在以下區(qū)間分析中省略隨機(jī)變量。設(shè)當(dāng)前迭代循環(huán)次數(shù)為s，當(dāng)前迭代點(diǎn)為Ys，功能函數(shù)進(jìn)行2次泰勒展開(kāi)，(24)式轉(zhuǎn)化為二次規(guī)劃問(wèn)題：

(27)

式中：Hs為海森矩陣。

由于(27)式為區(qū)間約束的二次規(guī)劃問(wèn)題，當(dāng)約束條件為區(qū)間約束時(shí)，將迭代方向投影到可行域上進(jìn)行計(jì)算就相當(dāng)方便，因此可以采用梯度投影法[19]高效地求得功能函數(shù)的極限值。首先將迭代方向投影于可行域，再在投影后的最速下降方向上搜索功能函數(shù)的第1個(gè)局部最優(yōu)點(diǎn)，最后利用Schur直接法求解最優(yōu)點(diǎn)。設(shè)最優(yōu)點(diǎn)為Ys(*)，則區(qū)間變量的迭代過(guò)程可表示為

Ys+1=Ys+γsYs(*)，

(28)

式中：γs為迭代步長(zhǎng)，為保證全局收斂，重復(fù)γs←γsρ，直至新迭代點(diǎn)Ys+1滿足Arimijo條件，即

(29)

θ=1×10-4，ρ=0.8.

若新迭代點(diǎn)Ys+1滿足KKT條件，則區(qū)間迭代停止，Yk+1=Ys+1；否則繼續(xù)區(qū)間迭代。

按照以上算法進(jìn)行迭代，如果滿足‖Uk+1-Uk‖≤ε1且|G(Uk+1,Yk+1)|≤ε2(ε1和ε2是足夠小的非負(fù)實(shí)數(shù))，則βL=βk+1，計(jì)算最壞情況下的失效概率Pfmax：

Pfmax=Φ(-βL).

(30)

綜上所述，本文提出的基于二參數(shù)尋優(yōu)MPP的混合可靠性分析算法流程如下：

1)將隨機(jī)變量X轉(zhuǎn)化為標(biāo)準(zhǔn)正態(tài)隨機(jī)變量U，得G(U,Y)=0；

2)初始迭代k=0，輸入初始值U0和Y0；

3)固定Yk，選擇η和λ的值(0≤η<1,λ>0)，在G=ηG(Uk)平面上確定Uk+1的搜索方向αk+1，并計(jì)算βk+1和Uk+1；

4)檢查KKT條件，如果滿足則Yk+1=Yk，轉(zhuǎn)步驟6，否則轉(zhuǎn)步驟5；

5)利用梯度投影法與Arimijo條件，計(jì)算Yk+1=Ys+1；

6)檢查收斂，如果滿足‖Uk+1-Uk‖≤ε1且|G(Uk+1,Yk+1)|≤ε2(ε1和ε2是足夠小的非負(fù)實(shí)數(shù))，則βL=βk+1，轉(zhuǎn)步驟7，否則k=k+1并轉(zhuǎn)步驟3；

7)計(jì)算最壞情況下的失效概率Pfmax=Φ(-βL)。

3 案例分析

為驗(yàn)證所提算法的精度、效率及穩(wěn)定性，本節(jié)通過(guò)兩個(gè)案例進(jìn)行分析。其中，收斂精度均取ε1=ε2=1×10-6，取λ>2，設(shè)定λ=λ/c，c=1.2以自適應(yīng)更新迭代步長(zhǎng)；N1和N2分別為混合可靠性分析循環(huán)序列迭代的迭代次數(shù)和函數(shù)調(diào)用次數(shù)。

3.1 懸臂梁失效概率計(jì)算

如圖7和圖8所示，在懸臂梁末端施加水平力Fx和垂直力Fy，懸臂梁長(zhǎng)度為l、橫截面寬度為w、高度為h，l(mm)、w(mm)和h(mm)服從正態(tài)分布。懸臂梁的最大應(yīng)力不超過(guò)屈服強(qiáng)度極限值S=380 MPa，F(xiàn)x和Fy為區(qū)間變量，其分布參數(shù)如表1所示。

圖7 懸臂梁Fig.7 Cantilever beam

圖8 懸臂梁橫截面Fig.8 Cross section of cantilever beam

表1 懸臂梁的不確定性特征Tab.1 Uncertainty characteristics of cantilever beam

圖9為本文算法與不同混合可靠性分析算法的收斂結(jié)果對(duì)比。由圖9可知，當(dāng)功能函數(shù)非線性程度較低時(shí)，相比于基于HLRF/iHLRF混合可靠性分析算法，本文所提算法同樣具有較高的收斂效率。表2為η、λ取不同值時(shí)的分析結(jié)果。由表2可知，隨著η的減小、λ的增大，本文所提算法精度基本不變，但收斂效率會(huì)逐漸提高。當(dāng)η在0～1之間取任意值時(shí)，λ可取大于0的任意值，且隨著λ的增大，收斂效率會(huì)提高；當(dāng)λ取大于0任意值并固定時(shí)，η可取0～1之間的任意值，且隨著η的減小，收斂效率會(huì)提高；當(dāng)η=0、λ→∞時(shí)，本文算法與基于HLRF的混合可靠性分析算法具有相同的計(jì)算效率，迭代次數(shù)均為4，函數(shù)調(diào)用次數(shù)為53.

圖9 本文算法與基于HLRF/iHLRF混合可靠性分析算法的收斂過(guò)程對(duì)比(η=0.1，λ=0.9)Fig.9 Comparison of the convergence process of the algorithm proposed in this paper and the mixed reliability analysis based on HLRF/iHLRF(η=0.1，λ=0.9)

為驗(yàn)證本文算法的精度，采取蒙特卡洛法(MCS)計(jì)算失效概率。首先將每個(gè)區(qū)間等分為10份，在區(qū)間變量的組合值下對(duì)隨機(jī)變量抽取1×108次，則功能函數(shù)調(diào)用1×1010次。表3為懸臂梁最壞情況下的失效概率。由表3可知，本文算法相對(duì)誤差為3.9%，與基于HLRF/iHLRF[19]混合可靠性分析算法的函數(shù)調(diào)用次數(shù)為同一數(shù)量級(jí)，收斂效率基本一致。綜上所述，本文算法通過(guò)對(duì)兩個(gè)參數(shù)的調(diào)整，既具有較高的計(jì)算精度，又具有較高的收斂效率。

表3 懸臂梁最壞情況下的失效概率Tab.3 Failure probabilities cantilever beam in the worst case

3.2 簡(jiǎn)支梁失效概率計(jì)算

如圖10所示，對(duì)簡(jiǎn)支梁施加均布載荷q和集中載荷F，梁的長(zhǎng)度為l，楊氏模量為E. 簡(jiǎn)支梁在C點(diǎn)的撓度不超過(guò)極限值δ=6 mm，q和F為區(qū)間變量，E、l、w和h服從正態(tài)分布，其分布參數(shù)特征如表4所示。本例為概率-區(qū)間混合的結(jié)構(gòu)可靠性問(wèn)題，求解簡(jiǎn)支梁在C點(diǎn)的撓度不超過(guò)極限值的最大可能失效概率。

圖10 簡(jiǎn)支梁Fig.10 Simple beam表4 簡(jiǎn)支梁的不確定性特征Tab.4 Uncertainty characteristics of simple beam

參數(shù)分布均值標(biāo)準(zhǔn)差下界上界l/mm正態(tài)81612w/mm正態(tài)1107h/mm正態(tài)904E/GPa區(qū)間200220q/N區(qū)間840895F/N區(qū)間16801825

表5為不同η、λ值時(shí)的失效概率分析結(jié)果。由表5可知，當(dāng)功能函數(shù)非線性程度較高時(shí)，隨著η取值減小，λ取值增大，分析結(jié)果存在不收斂的情況；當(dāng)η取值趨于0時(shí)，本文算法和基于HLRF的混合可靠性分析算法類似，會(huì)存在計(jì)算收斂不穩(wěn)定的問(wèn)題；當(dāng)η取值趨于1時(shí)，即收斂迭代步長(zhǎng)減小，能解決計(jì)算結(jié)果收斂不穩(wěn)定的問(wèn)題，使計(jì)算結(jié)果收斂；當(dāng)λ取值增大時(shí)，隨著沿梯度方向的步長(zhǎng)增加，收斂速度變快，但會(huì)引起收斂不穩(wěn)定的問(wèn)題，此時(shí)可以增加步長(zhǎng)調(diào)整系數(shù)c，以自適應(yīng)調(diào)整迭代步長(zhǎng)，解決計(jì)算收斂不穩(wěn)定的問(wèn)題。

表5 簡(jiǎn)支梁不同η、λ取值的分析結(jié)果Tab.5 Analyzed results of different η and λ values for simple beam

為驗(yàn)證本文算法的精度，將每個(gè)區(qū)間等分為10份，采用MCS在區(qū)間變量的組合值下對(duì)隨機(jī)變量抽取1×108次，則功能函數(shù)調(diào)用1×1011次。根據(jù)圖11的不同分析方法收斂結(jié)果對(duì)比圖及表6分析結(jié)果誤差與效率對(duì)比可知，與MCS相比，本文所提算法計(jì)算結(jié)果相對(duì)誤差為4.3%. 基于HLRF的混合可靠性分析算法，結(jié)果不收斂；基于iHLRF的混合可靠性分析算法，在迭代循環(huán)47次、函數(shù)調(diào)用526次后達(dá)到收斂；本文所提算法可在迭代循環(huán)50次、函數(shù)調(diào)用638次后達(dá)到收斂。以上對(duì)比結(jié)果表明，當(dāng)功能函數(shù)非線性程度較高時(shí)，通過(guò)選定合適的η和λ的值，本文算法能解決收斂不穩(wěn)定的問(wèn)題，同時(shí)具有較高的計(jì)算精度和效率。

圖11 非線性程度較高時(shí)本文算法(η=0.6，λ=10，c=1.2)與基于HLRF/iHLRF混合可靠性分析的收斂過(guò)程對(duì)比Fig.11 Comparison of convergence processes of the proposed algorithm (η=0.6，λ=10，c=1.2) and the mixed reliability analysis based on HLRF/iHLRF when the performance function having a high degree of nonlinearity 表6 簡(jiǎn)支梁最壞情況下的失效概率Tab.6 Failure probabilities of simple beam in the worst case

分析算法Pfmax相對(duì)誤差/%N1N2MCS0.011710111011基于HLRF混合可靠性分析算法基于iHLRF混合可靠性分析算法0.01224.345526本文算法0.01224.349638

4 結(jié)論

本文針對(duì)機(jī)械結(jié)構(gòu)中既有隨機(jī)變量又含區(qū)間變量的混合可靠性問(wèn)題，提出了基于二參數(shù)尋優(yōu)MPP的混合可靠性分析算法。所得結(jié)論如下：

1) 針對(duì)概率-區(qū)間混合的可靠性問(wèn)題，可將其轉(zhuǎn)化為雙層優(yōu)化模型，實(shí)現(xiàn)概率分析和區(qū)間分析共同作用的高效序列迭代算法，使計(jì)算結(jié)果具有較高的計(jì)算精度和效率。

2) 在概率分析時(shí)，通過(guò)引入兩個(gè)分別控制搜索方向和搜索步長(zhǎng)的調(diào)整參數(shù)，當(dāng)功能函數(shù)非線性程度較高時(shí)，同樣能使計(jì)算結(jié)果收斂，且具有較高的計(jì)算效率。