廖子涵, 李賓賓*

(1.浙江大學 伊利諾伊大學厄巴納香檳校區(qū)聯(lián)合學院,海寧 314400;2.浙江大學 建筑工程學院,杭州 310058)

1 引言

實際工程中,結(jié)構(gòu)面臨幾何尺寸、材料屬性、邊界條件和荷載等不確定性因素。因此,可靠性分析對于結(jié)構(gòu)設(shè)計優(yōu)化和結(jié)構(gòu)安全性評估等工作具有重要意義。其中結(jié)構(gòu)失效概率的計算通常通過以下三種方式求解,(1) 近似估計法,如基于功能函數(shù)泰勒展開的可靠度分析方法。該方法對于高度非線性問題與多失效點問題不適用。(2) 蒙特卡洛積分法,包括直接蒙特卡洛模擬及馬爾科夫鏈蒙特卡洛模擬MCMC(Markov Chain Monte Carlo)法。此類方法具有良好的通用性與準確性,但計算代價較高。(3) 代理模型法。常用的代理模型包括響應面法[1]、高斯過程模型[2]、Kriging模型[3]和支持向量機[4,5]等。當參數(shù)維度較大時,該方法需要提取特征進行降維分析,這必然導致信息丟失且無法確保失效概率估計的無偏性。

為解決高維參數(shù)空間下失效區(qū)域概率估計問題,Au等[6]提出子集模擬算法,利用條件概率引入中間事件,自適應地將一個小的失效概率分解為一系列相對較大的條件概率的乘積。通過MCMC采樣,將樣本經(jīng)由中間事件進行一系列篩選與生成,滾動地不斷接近參數(shù)空間中的失效區(qū)域。其作為一種較為高效的數(shù)值計算方法,廣泛應用于各類結(jié)構(gòu)可靠度的計算。此外,由于可靠度問題(小概率事件)、優(yōu)化問題(極值事件)和貝葉斯更新問題(極限狀態(tài)方程可以看作一類特殊的似然函數(shù))之間存在聯(lián)系,子集模擬近年也通過改造逐漸應用在多個領(lǐng)域。Suo等[7]將非支配排序引入樣本排序的規(guī)則,實現(xiàn)多目標優(yōu)化。廣義子集模擬[8]將子集模擬法的應用推廣到多極限狀態(tài)方程與時變可靠性分析領(lǐng)域[9]。

子集模擬中,MCMC算法的使用在增強采樣效率的同時導致樣本之間存在相關(guān)性。這種相關(guān)性不僅會增大條件概率的計算方差,還會使得算法的遍歷性下降,樣本無法對失效區(qū)域進行有效覆蓋,進而影響單次模擬結(jié)果的可信度。目前降低樣本相關(guān)性主要體現(xiàn)在參數(shù)優(yōu)化和MCMC采樣算法優(yōu)化兩個方面。其中相關(guān)參數(shù)的優(yōu)化集中于條件概率的選取,條件概率越低,模擬收斂的速度越快,但馬氏鏈之間相關(guān)性較高,遍歷性下降。相反的,條件概率越高,算法遍歷性越好,但層與層之間的相關(guān)性上升,模擬所需計算量增大,子集模擬算法的優(yōu)勢削弱。

采樣器的更新優(yōu)化是近年來子集模擬發(fā)展的主要方向。在最初的子集模擬中,Au等[6]采用的是改進MH(Metropolis Hastings)算法。為了減少馬爾科夫鏈內(nèi)的相關(guān)性,Santoso等[10]提出重復采樣版本的MH采樣法,通過重復采樣保證產(chǎn)生的待選樣本與之前樣本不同,然而該方法的人為選擇性,導致MH算法中的可逆性條件不能滿足,生成的樣本雖然相關(guān)性降低,但分布與目標分布之間存在誤差[11]。此外,通過對種子樣本進行拉伸采樣[13],也能從馬爾科夫鏈的起始點降低鏈與鏈之間的相關(guān)性。更復雜地,Miao等[13]提出再生自適應子集模擬算法,通過隨機再生建議分布以及延遲拒絕算法,克服burn in問題并且盡可能多地獲得獨立樣本。然而,在子集模擬中馬爾科夫鏈的種子分布已經(jīng)自動滿足目標分布,因此再生建議分布理論上對于解決burn in問題不起作用[14]。為了方便在高維空間中對多種分布組成的概率密度函數(shù)進行蒙特卡羅采樣,Papaioannou等[11]提出自適應條件采樣aCS(adaptive Conditional Sampling),將變量θ通過Nataf變換[15]轉(zhuǎn)化為獨立和標準正態(tài)分布的變量u,并在標準正態(tài)空間U中進行采樣與數(shù)值積分。該方法在簡便性與準確性方面都較為突出,獲得了廣泛的應用。為了提升樣本探索空間的效率,將Hamiltonian采樣嵌入子集模擬法中[16,17],該方法基于功能函數(shù)的梯度矩陣,提升樣本向函數(shù)值較高的區(qū)域移動的速度,但同時會增加計算梯度矩陣的負擔。連續(xù)空間轉(zhuǎn)換[18]通過控制變量方法先給出一個粗略的失效概率估計值,再通過修正項進行細化修正,通過擴大樣本分布的方差以增加樣本的遍歷性。

2 子集模擬法

2.1 子集模擬法

定義n維結(jié)構(gòu)參數(shù)θ∈n及其概率密度函數(shù)fn(θ),對于代表結(jié)構(gòu)失效模式的極限狀態(tài)方程g(θ),結(jié)構(gòu)失效概率可表達為

(1)

式中IF(θ)為失效事件F={θ∈n∶g(θ)≤0}的指示函數(shù),即

(2)

式(1)可通過直接蒙特卡洛采樣積分方式近似求解,然而結(jié)構(gòu)失效一般為小概率事件,失效域僅占參數(shù)空間的極小部分,大量樣本點落在安全域,造成計算資源的嚴重浪費,特別是在參數(shù)θ的維度n很大的情況下。

子集模擬算法將結(jié)構(gòu)失效事件F表示為一系列逐級包含的失效事件的交集

(3)

F1?F2?…?FM=F

(4)

式中Fj={θ∈n∶g(θ)≤bj},b1>b2>…>bM=0。利用條件概率的性質(zhì),原結(jié)構(gòu)失效概率轉(zhuǎn)化為一系列條件失效概率的乘積

(5)

式中F0={θ∈n∶g(θ)≤∞}為確定事件,即Pr(F0)=1,Pr(F1|F0)=Pr(F1)。每一層的條件失效概率計算可參照式(1)進行,

(6)

對應的條件概率密度函數(shù)為

(7)

由于條件概率密度函數(shù)fn(θ|Fj-1)不是標準的概率分布,如何實現(xiàn)其快速有效采樣成為實現(xiàn)失效概率估計的關(guān)鍵。基于這一觀察,子集模擬可看作是一種自適應重要性采樣算法,從初始概率密度函數(shù)fn(θ)出發(fā),通過給定條件失效概率p0=Pr(F1|F0)(一般取為0.1或0.2)確定分界點b1和對應失效域F1,然后以失效域F1內(nèi)采樣點作為種子,基于MCMC原理實現(xiàn)fn(θ|F1)的采樣并進而確定新的分界點b2和失效域F2,以此類推,直至達到目標失效域FM=F。從上述過程可以看出,子集模擬選用MCMC采樣,每一層樣本基于上一層傳遞下來的種子樣本,從而實現(xiàn)條件概率密度函數(shù)fn(θ|Fj-1)的自適應采樣,且種子本身滿足分布fn(θ|Fj-1),大大提升了候選樣本的接受率。

圖1 樣本運動Fig.1 Sample movement

(8)

式中Tn(·|·)為轉(zhuǎn)移概率密度。通過隨機數(shù)與接受率比較后篩選得到候選樣本ξ。此時,候選樣本的總體分布介于fn(θ|Fj-1)與fn(θ)之間,近似于從fn(θ|Fj-1)擴散到fn(θ)的過渡階段。第二步,基于指示函數(shù)IFj-1(θ)對候選樣本ξ進行二次篩選,將超出邊界的樣本點收縮回初始點,得到更新后的樣本θk+1,即

(9)

(10)

即滿足細致平衡條件。式(10)利用IFj-1(θk)=1和函數(shù)min{·}的性質(zhì)。需要指出的是,直接應用MCMC采樣并不適用于高維隨機采樣,也就是存在高維災難問題,本文將簡要介紹兩種常用的高維采樣算法,從而使得子集模擬能夠解決高維可靠度估計問題。

2.2 高維隨機采樣

2.2.1 改進MH采樣

(11)

2.2.2 自適應條件采樣

(12)

3 混合采樣算法

3.1 橢圓切片采樣

橢圓切片采樣ES(Elliptical slice sampling)[19]為切片采樣在正態(tài)空間的一種特殊形式。與條件采樣和哈密頓采樣類似,首先對現(xiàn)有樣本u進行擾動

v=ucosβ+dusinβ

(13)

式中β～U[0,2π],du～Nn(0,1)。設(shè)定下一步角度選取的上下界[βmin,βmax]=[β1-2π,β1]。使用公式對侯選樣本進行篩選,若拒絕,則將β1設(shè)定為分布的上界或下界

βmin=β(β<0)
βmax=β(β≥0)

(14)

基于更新后的上下界均勻隨機取值β并重復以上步驟,直到樣本接受。該方法本質(zhì)上為一種收縮條件采樣(或正態(tài)分布下的軌跡收縮哈密頓采樣)。將邊界條件的篩選步驟與新樣本的生成結(jié)合在一起,保證新產(chǎn)生的樣本在滿足邊界條件的前提下互相不重復,增加了等效獨立樣本的數(shù)量,并且每次對角度的搜索從全分布出發(fā),增強了樣本的遍歷性。相對的,橢圓切片采樣收縮過程中拒絕步驟的嵌入也會帶來較大的計算量。如圖2所示,橢圓切片采樣在單次產(chǎn)生步驟的樣本范圍大于自適應條件采樣,帶來較好遍歷性的同時降低收縮步的接受概率。

圖2 產(chǎn)生步樣本分布Fig.2 Sample distribution in Generation step VS aCS sampler

通過分析可知,在可靠度分析領(lǐng)域,橢圓切片采樣范圍縮減次數(shù)的上限與問題維度和失效概率相關(guān)。每一次范圍的縮減,由于均勻取值的特性,每個維度縮減比例的期望為0.5。若該問題失效區(qū)域在每個維度上占比相等,子集模擬最后一層內(nèi)(所有種子樣本皆位于失效區(qū)域,所需范圍縮減次數(shù)最大)采樣范圍縮減次數(shù)的期望可以表示為

log0.5nPf=log0.5Pf/n

(15)

式中n為問題的維度?？梢钥闯?橢圓切片采樣的效率隨著問題的維度上升而升高(滿足上述假定的前提下),但若存在極端情況,即無論問題維度,失效區(qū)域的占比始終集中于一維,其余維度采樣范圍不受目標函數(shù)影響,則橢圓切片采樣的效率將不隨維度變化。

為保證算法效率,在橢圓切片采樣中設(shè)置范圍縮減次數(shù)的上限Nmax,以避免目標函數(shù)的過多調(diào)用?？紤]到上述極端情況的存在與結(jié)構(gòu)失效概率的量級,本文設(shè)置Nmax=20,其在極端情況下仍能采樣到失效區(qū)域(Pf≥0.520≈9.54×10-7)內(nèi)的樣本。

3.2 混合采樣子集模擬

本文基于上述采樣方法提出了一種搭配混合采樣器的子集模擬方法。與傳統(tǒng)子集模擬始終調(diào)用一種采樣器不同,該方法在一次模擬內(nèi)先后調(diào)用橢圓切片采樣與自適應條件采樣兩種采樣器。如圖3所示的示意性算例中,條件失效區(qū)域分別標記為A,B和C,隨著模擬的三個區(qū)域不斷縮小

圖3 自適應條件采樣器和混合采樣器Fig.3 aCS sampler VS Mixed sampler

(g(θ)

為避免使用橢圓切片采樣器的層數(shù)過多導致計算資源浪費,需要設(shè)置自適應的采樣器切換邊界以提升新算法的通用性。本文中當橢圓切片采樣器函數(shù)調(diào)用次數(shù)大于樣本數(shù)的4倍時,判斷橢圓切片采樣器效率過低,切換為自適應條件采樣器。

4 算例分析

援用4個算例(3個可自行控制維度的高維算例以及一個高度非線性算例),從不同角度檢驗新算法的特性并對比與其他算法的優(yōu)缺點。

4.1 n維線性超平面[20]

極限狀態(tài)方程為

(16)

式中Xi(i=1∶n)為獨立的標準正態(tài)分布變量。該算例中β值控制失效概率的大小,且失效概率不隨維度變化,又因為該算例為線性方程,使用FORM可以找到失效概率的解析值,給算法的測試提供了較為便利的環(huán)境,可較方便地用來研究算法在不同概率水平和不同維度下的性能。本文取β=6,對應的失效概率Pf=9.87×10-10[20]。

4.2 n維非線性超曲面1[21]

極限狀態(tài)方程的非線性集中于前兩維

(17)

式中Xi(i=1∶n)為獨立的標準正態(tài)分布變量。k為控制函數(shù)非線性的曲率參數(shù),k越大,失效概率越小,對算法要求越大。本文取k=-10,β0=4,對應的失效概率Pf=4.29×10-6[21]。

4.3 n維非線性超曲面2[22]

第二個非線性算例的極限狀態(tài)方程為以下形式,除第n維外,其他維度皆存在非線性。

(18)

式中Xi(i=1∶n-1)為獨立的標準正態(tài)分布變量,Xn滿足均值為1,標準偏差為0.2的正態(tài)分布。與4.1節(jié)線性超平面算例不同的是,該算例的失效概率與維度相關(guān)且失效概率的值無法通過計算直接得到。因此本文將多次蒙特卡羅積分計算結(jié)果的平均值作為失效概率的準確值[22]。

4.4 高度非線性問題[23]

最后一個算例為白噪聲激勵下具有不確定阻尼振子的雙自由度系統(tǒng),如圖4所示,極限狀態(tài)方程為

圖4 白噪聲激勵下的雙自由度系統(tǒng)Fig.4 Two degrees-of-freedom system under white noise excitation

(19)

表1 算例4涉及參數(shù)分布Tab.1 Distribution of parameters in example 4

由表2可知,所有算例在各類采樣器下運行5000次,保證失效概率浮動基本不影響偏差判斷,結(jié)果列入表2。對于算例1,除了改進MH采樣外,其余算法基本不存在偏差,橢圓切片器由于其采樣效率較低,結(jié)果的變異系數(shù)較大?；旌喜蓸悠鞯淖儺愊禂?shù)位于自適應條件采樣與橢圓切片采樣之間。在該算例中,所有維度的采樣范圍都隨著模擬層數(shù)的增加不斷縮減,因此即使該算例對應的失效概率較低,橢圓切片采樣仍能取得較好的效果。對于自適應條件采樣,該算例的極限狀態(tài)方程為線性,對于遍歷性要求較低,因此也能取得較好的效果。

表2 各算法數(shù)值積分結(jié)果Tab.1 Results of different MCMC algorithms

算例2中,自適應條件采樣的偏差與變異系數(shù)都明顯好于橢圓切片采樣,混合采樣偏差與自適應條件采樣接近,變異系數(shù)位于兩者之間,明顯優(yōu)于傳統(tǒng)的改進MH采樣。在該算例中,極限狀態(tài)方程的非線性與失效區(qū)域的收縮集中在前兩個維度,對于自適應條件采樣較為有利,而對于橢圓切片采

樣的效率不利。在模擬的最后幾層中,橢圓切片采樣在前兩維的采樣范圍縮減嚴重,造成樣本的大量重復,也導致其模擬結(jié)果較差。

對于第三個算例,橢圓切片采樣的偏差相較于自適應條件采樣更優(yōu),變異系數(shù)則較大?；旌喜蓸悠钆c橢圓切片采樣接近,變異系數(shù)位于兩者之間。在該算例中,極限狀態(tài)方程的非線性與失效區(qū)域的收分散于所有維度,對于采樣算法的遍歷性要求較高,混合采樣在保證遍歷性的前提下較橢圓切片采樣效率更高,因此數(shù)值積分結(jié)果的變異系數(shù)更低。

第四個算例為低維高非線性低失效概率算例,所有采樣器除橢圓切片采樣外偏差接近,其中自適應條件采樣的變異系數(shù)最小。該算例說明自適應采樣對于低維非線性問題依舊具有良好的遍歷性。與算例2類似,橢圓切片采樣對此類問題存在樣本點難以移動的問題,導致最終模擬結(jié)果不理想?；旌喜蓸拥男Ч麅?yōu)于橢圓切片采樣與傳統(tǒng)的改進MH采樣。

5 結(jié) 論

為解決單一采樣器帶來的不同問題下通用性較差的缺點,本文提出了一種混合采樣的子集模擬法。經(jīng)各類算例驗證,混合采樣器相比于其他單一類型的采樣器具有以下優(yōu)勢。

(1) 同時具有橢圓切片采樣器遍歷性較好與自適應條件采樣器采樣效率高的優(yōu)點。

(2) 通過調(diào)用次序不同,可有效避開兩種采樣器的短板,提升算法的魯棒性。

(3) 通過自適應的切換機制,自動選擇各階段適合的采樣器,對于不同問題具有良好的通用性。

(4) 繼承了兩種采樣器對于高維問題的適用性,對于算例的各類高維問題都有較優(yōu)秀的表現(xiàn)。

然而該算法依舊存在以下有待發(fā)展的方面。

(1) 橢圓切片采樣部分依舊不夠高效,導致同等函數(shù)調(diào)用次數(shù)下,積分結(jié)果的方差依舊大于單一自適應條件采樣,未來可考慮用替代模型減少函數(shù)調(diào)用次數(shù),增強算法效率。

(2) 自適應的采樣器切換邊界只基于函數(shù)調(diào)用次數(shù)與樣本數(shù)的比值,較粗糙。此外,層內(nèi)采樣過程也可考慮調(diào)用不同采樣器以達到更佳的效果。

目前算法只包含橢圓切片采樣與自適應條件采樣兩種采樣器,其他方法如哈密頓采樣和拉伸采樣等具有各自不同特點的采樣器也可以加入采樣器候選庫中,使得在子集模擬過程中,可針對模型各階段的特性選擇最適合的采樣方法。