基于梯度搜索法的非線性結構可靠性分析方法

陳鵬霏， 和鵬， 于泰龍， 劉巧伶

(1.長春工業大學 機電工程學院, 吉林 長春 130012; 2.吉林大學 機械與航空航天工程學院, 吉林 長春 130025)

0 引言

工程結構普遍具有非線性特征強和失效概率小的特點。并且，對于復雜的工程結構，其模型通常還具有隱式特征。如何高效、精確地計算隱式非線性結構的可靠度及可靠性靈敏度，一直是結構可靠性分析領域的研究熱點之一[1-2]。

目前，常見的結構可靠性分析方法，主要包括解析法、隨機模擬法和抽樣擬合法等。其中，一次二階矩法、二次二階矩法和高次高階矩法屬于傳統解析法[3-5]。一般需要根據結構功能函數的Taylor級數展開式，通過保留1階項，或者高階項來近似計算功能函數的均值和標準差，進而計算可靠度。因此，采用解析法的前提是需要已知功能函數的顯式表達式，故不適用于隱式結構可靠度的求解。隨機模擬法主要是以Monte Carlo法為代表，并以此為基礎，發展了重要抽樣(IS)法和自適應重要抽樣法等[6-8]。隨機模擬法對于顯式或者隱式工程結構的可靠度求解，均具有普遍的適用性，而且在解決強非線性問題時也很有效。然而，隨機模擬法的計算精度，是以大量的抽樣數據為基礎，因此在計算大型復雜結構的可靠度時，其計算效率或成本往往難以接受。針對上述問題，為了高效、精確地計算大型復雜工程結構的可靠度，抽樣擬合法[9-13]成為目前可靠性定量分析領域的研究熱點。該方法基本思想是，基于響應面法(RSM)、神經網絡、Kriging法等，利用有限個樣本點建立起輸入、響應之間的顯式函數關系，來替代原結構的隱式功能函數，再通過傳統解析法或隨機模擬法計算其可靠度。朱麗莎等[11]將有限元法和人工神經網絡技術相結合，研究了隨機轉子系統中工作參數對轉子系統振動可靠性的影響，并進行了排序。佟操等[12]采用Krigring模型提出了一種主動學習的可靠度計算方法，相比同類方法具有所需樣本數量少和計算精度高的特點。同時，目前可靠性分析領域采用的抽樣擬合法，基于非線性“黑箱”模型(二次多項式、Krigring模型、神經網絡模型等)，國內外諸多學者提出了各種選點策略來提高可靠性分析的精度。尤其是近年來，基于主動學習的非線性“黑箱”模型，通過在極限狀態方程或極限狀態表面附近，選擇新訓練樣本來擬合極限狀態方程，而非在整個狀態空間內擬合功能函數。因此，其計算精度和計算效率均得到很大的提高[2]。

本文以梯度搜索法為基礎，提出一種新的抽樣擬合法，來進行結構可靠性靈敏度分析。首先，通過高次梯度搜索法，反復迭代尋找極限狀態表面附近的訓練樣本點；之后，采用多項式函數或響應面函數擬合出結構的極限狀態方程，進行可靠度和可靠性靈敏度分析。同時，在數值和工程算例中，將本文所提方法與以往常用方法的分析結果，以及采用Monte Carlo法獲得的真值結果進行對比，證明了本文方法具有較高的計算精度和計算效率。

1 一次可靠性分析法

1.1 標準化、正態化狀態空間

假設向量x=[x1,x2,…,xn]中包含有n個服從任意分布且相互獨立的隨機變量，其功能函數為

zX=gx(x1,x2,…,xn)，

(1)

式中：zX表示狀態空間X內結構的功能函數；gx(·)為與向量x對應的函數方程。當極限狀態zX=0時，則極限狀態方程表示以x1,…,xn為坐標軸構成的狀態空間X中的一個超曲面。

在工程實際中，對于隨機變量xi(i=1,…,n)不服從正態分布的情況，可以采用映射變換的方式來實現隨機變量的標準化和正態化[3]。

根據映射變換法作映射變換，則有

Fxi(xi)=Φ(yi)，

(2)

(3)

(3)式代入(1)式，可得

(4)

式中：zY表示標準正態空間Y內結構的功能函數；gy(·)為與隨機向量y對應的函數方程。

所以，任意非正態分布隨機變量均可轉化為標準正態分布隨機變量。故此，本文僅討論相互獨立的正態分布隨機變量情形。

1.2 基于經典驗算點法的可靠性分析

圖1 極限狀態方程和驗算點示意圖Fig.1 Sketch map of limit state equation and checking point

若已知結構的功能函數，則失效概率[15]為

(5)

根據可靠度指標β的定義，可知驗算點對應的β*為

(6)

依據Taylor展開法，將(4)式所示的功能函數在驗算點y*處展開為線性化方程：

(7)

式中：?gy/?yi=?gx/?xi·?xi/?yi，?gx/?xi在隨機向量x的驗算點x*處計算，?xi/?yi在隨機向量y的驗算點y*處計算，根據(3)式可知：

(8)

φ(·)為標準正態概率密度函數，fxi(·)為變量xi的概率密度函數。

(7)式中右邊各項均為隨機向量y的線性組合，因此zYL也應服從正態分布，其均值μzYL和標準差σzYL分別為

(9)

式中：σyi為變量yi的標準差。

故結構的失效概率可近似表示為

Pf≈P{zYL≤0}=Φ(-βzYL)，

(10)

式中：βzYL為線性化極限狀態函數zYL的可靠度指標[16]，且

βzYL=μzYL/σzYL.

(11)

根據(7)式、(9)式和(11)式，可推導出標準正態空間Y內驗算點y*的坐標值：

(12)

式中：cosθyi稱為驗算點重要度系數，

(13)

最后，根據(3)式，可得基本隨機向量x的驗算點x*. 由推導過程可知，當采用一次可靠性分析方法計算失效概率Pf時，將極限狀態方程線性化是造成誤差的主要原因。

2 梯度搜索法計算樣本點

2.1 一次梯度搜索

由(12)式可知，一次可靠性分析中的驗算點可通過重要度系數逐步迭代獲得。因此，根據經典驗算點法(FORM)，關于非線性結構的一次梯度搜索過程如下：

1) 依據工程經驗找到失效域中一點，或直接將均值點x(1,0)確定為初始點，并將其轉換為標準正態空間內點y(1,0)；

3) 通過以上過程反復迭代nS1次，最終確定出貼近極限狀態曲面，且滿足精度要求的抽樣點y(1,nS1)作為一次梯度搜索的輸出響應點y*(1)，且可由(3)式將其轉化為狀態空間X中的點x*(1).

2.2 高次梯度搜索

通過一次梯度搜索獲得的響應點x*(1)，就是FORM中所說的驗算點[14]。根據點x*(1)建立的近似極限狀態方程屬于線性化方程，這是一次可靠性分析方法得到的失效概率Pf存在計算誤差的根本原因。因此，還需要在此基礎上獲得貼近極限狀態曲面的其他樣本點。具體過程如下：

1) 根據(7)式過響應點y*(1)建立結構極限狀態曲面的1階線性平面方程：

(14)

(14)式為標準正態空間Y中的一個超平面。

2) 以點y*(1)為中心，選取合適的r作為半徑，構建標準正態空間Y中的超球面方程：

(15)

(16)

式中：j=1,…,k-1,k+1,…,s-1,s+1,…,n.

圖2 基于梯度搜索法獲取極限狀態曲面附近樣本點流程框圖Fig.2 Flow chart of obtaining the sampling points near limit state surface by gradient search method

圖3 基于梯度搜索法的抽樣點分布與擬合函數示意圖Fig.3 Sketch map of sampling point distribution and fitting function based on gradient search method

超球面半徑r的選取，與參數空間內隨機變量的標準差有關，且影響結構失效概率計算過程中的效率和收斂性。當初次搜索時，可取步長系數α=1.0，之后逐次翻倍。經過若干次迭代搜索之后，若半徑r≥3rR時仍不能獲得令人滿意的結構失效概率，可將α值減半，即令α=0.5，重新進行取點搜索，并令α每次增長奇數倍；當半徑r再次超過3rR時，令α繼續減半，且每次增長奇數倍，直到獲得滿意的結構失效概率值為止。

2.3 隱式功能函數偏導數的計算

在獲取梯度信息和可靠性靈敏度分析過程中，偏導數的計算都至關重要。對于顯式功能函數，可直接通過求導來計算偏導數。而對于隱式功能函數，則一般通過差分法來計算函數gx(x)的偏導數[17]。這里，常用的差分法包括中心差分法和向前差分法兩種。中心差分法：

(17)

式中：Δxi表示向量x中變量xi的微小變化量；功能函數gxi+Δxi=g(x1,…,xi+Δxi,…,xn)；功能函數gxi-Δxi=g(x1,…,xi-Δxi,…,xn)。

向前差分法：

(18)

式中：功能函數gxi=g(x1,…,xn)。

因此，在計算隱式功能函數的偏導數時，兩種差分法所需的抽樣點數N*略有不同。采用中心差分法，N*=n0(2n+1)+1，n0為迭代次數；采用向前差分法，N*=n0(n+1)+1. 顯然，中心差分法所需抽樣點數較多，但其計算精度要優于向前差分法。

同時，在采用差分法計算偏導數時，差分步長的選取對于計算精度起到重要作用。步長選取過大，則誤差變大。但是，若步長選取過小，在實際工程中又會受到測試儀器精度等級的限制，影響輸出值的讀取精度，同樣會造成誤差。因此，在實際工程中，步長的選取需要參照測試儀器的最高測試精度等級，即在保證測量精度的前提下，步長應越小越好。本文在采用差分法計算偏導數時，考慮到各隨機參數的分布變化特征，以參數標準差的1/15作為求取偏導數的差分步長，來進行計算的。

3 可靠度及可靠性靈敏度計算

3.1 非線性極限狀態方程的擬合

(19)

(20)

因此，當抽樣的階次m無限大時，擬合函數方程所表達的超曲面，與結構真實的極限狀態方程之間，誤差將趨于0，即

(21)

3.2 計算可靠度及可靠性靈敏度

根據結構可靠性分析理論，失效概率Pf為基本隨機變量y=[y1,y2,…,yn]的聯合概率密度函數在失效域F內的積分[8]，即

(22)

式中：失效域F={y:y∈gy(y)≤0}。

(23)

(24)

式中：σt為第t個隨機變量的標準差。

對于復雜隱式極限狀態函數zY=gy(y)，失效域邊界gy(y)=0無法解析表達。此時，(22)式和(23)式的積分只能采用數字模擬的方法來進行計算，計算量非常大。但是，依據(19)式，可獲得結構失效域邊界函數的近似解析表達式，即y(y)=0. 因此，可直接采用數值積分的方法來進行求解，將大大減小計算量，提高計算效率。

4 算例

本節將通過數值和工程算例，來驗證所提方法的計算精度和成本效率。

由圖3可知，將均值點作為一次梯度搜索的初始點，可獲得對應的驗算點。根據該驗算點，即可獲得1階線性方程。選取適當步長，在1階線性方程上可繼續取2個2次抽樣初始點，繼續利用梯度搜索法可獲得另外2個對應的2次驗算點。根據這2個2次驗算點和上述驗算點(此時，共3個點)，可獲得拋物線方程，即2階曲線方程。將步長加倍后，在1階線性方程上繼續取2個3次抽樣初始點，于是又可獲得2個3次驗算點。根據3次驗算點、2次驗算點和驗算點(共5個點)，可擬合得4階曲線方程。至此，該4階曲線方程已與極限狀態方程非常近似。若滿足精度要求，即可跳出程序循環，結束計算。否則，還可以繼續搜索獲得更高階的擬合曲線方程。

采用本文所提方法算得的失效概率結果Pf、隨機變量靈敏度結果，以及所需要的抽樣次數N*如表1所示。為便于比較，表1中同時列出了采用Monte Carlo可靠性抽樣(MC_RS)法、IS法、RSM，以及FORM的計算結果。并且，將MC_RS大樣本模擬法獲得的計算結果作為精確解(即真值)，與其他方法相比較，可得到各種分析方法的絕對誤差和相對誤差(相對誤差=絕對誤差/真值)。

表1 算例1可靠性靈敏度分析結果和抽樣點數一覽表

由表1數據可知，在各種可靠性靈敏度計算方法中，與MC_RS法算得的真值比較，IS法計算結果的誤差較小，因此精度較高。但是，需要的抽樣次數也較大，因此并不適用于計算量較大的大型復雜結構可靠性靈敏度分析。而RSM采用含有交叉項的多項式函數來擬合出響應函數，盡管計算效率很高(對于算例1，根據中心復合法抽樣，在兩個變量情況下僅需5個樣本點)，但是RSM抽取樣本點的目的是為了構建結構的狀態方程，而非極限狀態方程，因此計算精度較差。由表1可知，對于強非線性結構，RSM的計算結果幾乎失真。FORM通過迭代計算，獲得極限狀態曲面上聯合概率密度值最大的點，并將其作為驗算點，把通過驗算點的1階線性平面近似看做結構的極限狀態曲面。對于弱非線性結構可靠度計算而言，FORM的計算精度和效率均尚可。但是，對于如算例1所示的強非線性結構，FORM的計算精度并不令人滿意。

為了更加清晰地比較RSM、FORM與本文所提方法在計算精度和效率(即抽樣點數)方面的差別，在采用本文方法分析算例1時將其所經歷的抽樣點數和對應結果列于表2中。

表2 算例1抽樣點數和對應結果一覽表

綜上所述，本文所提方法是在FORM的基礎上，通過梯度搜索不斷獲得強非線性結構極限狀態曲面上的樣本點，最后采用合適的擬合函數擬合極限狀態方程，計算可靠性靈敏度。表1中數據顯示，無論是失效概率，還是各參數的可靠性靈敏度，本文所提方法的計算精度，均與MC_RS法計算的真值非常貼近，因此能夠對工程實際產生有效的指導作用；同時與解決強非線性結構可靠性靈敏度分析的其他方法(如IS法)比較，本文所提方法的抽樣點數很少。所以，更加適用于計算工作量較大的大型復雜結構的可靠性靈敏度分析。

(25)

圖4 液壓缸運動的力學模型Fig.4 Mechanics model of hydraulic cylinder

當液壓缸低速運動，即油液的平均流速v較小時，由(25)式可知，其運動呈現出顯著的非線性特征。隨著v的變化，其運動狀態一般可分為邊界潤滑、混合潤滑和動壓潤滑3種情況[19]，對應的速度時域圖像，分別如圖5(a)、圖5(b)、圖5(c)所示。

圖5 不同潤滑區域的活塞瞬時速度時域圖像Fig.5 Time-domain images of piston instantaneous speed in different lubrication areas

g(XH)=tslim-ts，

(26)

采用本文所提方法，對(26)式進行可靠性靈敏度分析，計算結果如表3所示。同時，將MC_RS大樣本模擬法、IS法、RSM和FORM計算出的可靠性靈敏度結果值，一同列入表3中。在本算例的數值計算過程中，涉及偏導數的計算，采用的是中心差分法。雖然中心差分法的抽樣次數較向前差分法多(本例若采用向前差分法，則抽樣點數為27)，但計算精度較高。

表3 液壓缸活塞運動可靠性靈敏度計算結果及抽樣次數一覽表

通過對比表3中數據可知，采用本文所提方法計算的失效概率和可靠性靈敏度，與MC_RS大樣本模擬法獲得的真值相比，均比較接近，能夠正確反映液壓缸的運動可靠性，及各參數對可靠性的影響情況。同時，從抽樣次數來看，本文所提方法無疑具有更高的計算效率，尤其適用于計算工作量較大的大型非線性結構可靠性靈敏度分析。

此外，如表3所示，通過對比其他結構可靠性分析方法可知，由于IS法的計算精度和效率嚴重依賴于重要抽樣密度函數的選取，而重要抽樣密度函數一般需憑經驗選取。所以，在不明確失效域的情況下，重要抽樣法的計算精度和效率會大打折扣。因此，在本算例中IS法的計算精度和效率均不如本文方法。RSM雖然在計算成本(或抽樣點數)上少于本文方法，但是RSM法采用有限個樣本點擬合出的響應面函數，來對結構整個狀態空間內的功能函數進行替代，很多情況下其計算精度，尤其是各隨機參數靈敏度的計算精度難以滿足工程需要。在本算例中，RSM的失效概率計算精度最低；而其隨機參數靈敏度的計算結果，與真值結果(采用MC_RS法的計算結果)相比，由于數量級相差過大，故已失真。這里，FORM的抽樣點數少于本文方法。但是，如前所述，FORM屬于一次可靠性分析方法，由于其在計算過程中采用線性平面直接替代結構真實的極限狀態表面。因此，當結構的非線性特征越強時，其計算精度會越差。由于本算例屬于非線性工程結構算例，故FORM在計算精度上不如本文方法。

5 結論

本文針對強非線性結構，提出一種基于梯度搜索法的可靠性靈敏度分析方法。通過逐次迭代搜索極限狀態表面上的樣本點，擬合結構極限狀態方程計算可靠度及可靠性靈敏度。通過算例分析，得出主要結論如下：

1)與MC_RS和IS等隨機模擬法相比，本文所提方法具有計算成本小、效率高的特點，尤其適用于大型復雜工程結構的可靠性靈敏度分析。

2)與通過狀態空間中樣本點擬合結構功能函數或極限狀態方程的抽樣擬合法相比，本文所提方法具有更高的計算精度。