基于近似貝葉斯計算的包裝件模型選擇和參數估計

朱大鵬, 曹興瀟

(1. 蘭州交通大學 交通運輸學院, 蘭州 730070; 2. 蘭州交通大學 機電工程學院, 蘭州 730070)

包裝件主要由產品和緩沖材料構成,在運輸過程中,在外界振動載荷作用下,包裝件可能出現疲勞、首次穿越等方式的損壞,影響包裝件的運輸安全和可靠性。準確構建包裝件模型,是評價包裝件運輸安全性和可靠性的重要基礎,也是優化包裝設計,實現合理包裝的重要依據。

不考慮包裝件中產品具體動態特性,將產品看作剛性質量塊,僅考慮緩沖材料彈性特性和阻尼特性,包裝件可近似建模為單自由度支座激勵系統[1-2]。考慮到緩沖材料的黏彈性,文獻[3]采用復指數函數模型,文獻[4]采用分數階微分模型,分別建立了緩沖材料模型。以上模型能夠更加準確分析包裝件動態、準靜態特性。葛正浩等[5]采用有限元法對包裝件中緩沖材料建模,并識別出材料的模型參數。為準確分析包裝件的響應,需要研究緩沖材料靜態和動態力-變形特性,近年來,研究者從試驗分析、材料結構等角度出發,構建了緩沖材料的本構模型,這些模型解釋了緩沖材料在準靜態和動態壓縮載荷作用下緩沖材料產生的力學現象,為建立包裝件模型和進行響應分析奠定了重要基礎[6-8]。在文獻[9]中,作者對目前主要的運輸包裝系統的動力學模型進行了綜述。

對包裝件模型識別時,首先需要確定模型的類型,即給定一個包裝件,可能有多個類型的模型(如線性模型[10]、非線性彈性模型[11]、非線性阻尼模型、包含干摩擦的模型[12]、包含有遲滯元素的模型等)均可以在一定程度上表征其振動特性,需要在這些模型中選擇一種最準確的模型對包裝件建模和分析。結構模型選擇是一個重要工程問題,是實現準確構建結構模型的重要基礎,也是實現模型參數估計的重要前提。Fuentes等[13-14]將振動系統中可能的彈性特性和阻尼特性函數放入備選庫中,在稀疏性約束下采用Lasso回歸法優選出表征振動系統響應的基函數。De等[15]結合貝葉斯推斷和模型證偽,彌補了單一貝葉斯推斷和模型證偽法的缺點,提出了一種優選振動模型的方法。Safari等[16]在非線性振動系統共振頻率附近實施諧波激勵,根據系統自由相應,提出一種基于最優化方法的非線性振動系統模型選擇和參數識別的方法。Civera等[17-18]應用機電系統相似性原則,構建松弛向量擬合法識別系統特性識別的方法,并用于結構模態結構特性識別。

與其它振動結構相似,對包裝件建模時,由于以下因素,導致包裝件模型存在著一些誤差:

(1) 由于包裝件中緩沖材料本構模型未完全建立,緩沖材料在載荷作用下彈性變形和能耗機理仍需進一步探索,故只能采用各種常見模型近似表征包裝件振動特性,因此選擇的一定類型包裝件模型與真實包裝件特性之間必然存在著誤差[19]。

(2) 振動試驗時,由于只能采集有限長度的振動數據識別包裝件模型,造成模型識別的誤差。

(3) 采集的振動試驗數據中不可避免地存在著噪聲,在識別模型時,這些噪聲數據會影響建模準確性。

考慮以上因素,如果采用確定性模型對包裝件建模,則模型分析結果和真實響應之間必然存在較大誤差,為準確表征包裝件模型的不確定性,需將模型參數或模型結構的不確定性引入振動模型中。近年來,在不確定框架下,采用貝葉斯推斷構建振動系統模型在結構特性識別、故障診斷等得到了廣泛應用。Beck等[20-21]考慮了結構不確定性,采用貝葉斯推斷法建立了結構模型識別方法。Worden等[22]采用馬爾可夫鏈蒙特卡洛法識別振動系統的參數,采用偏差信息準則從備選模型中選擇出最佳模型。傳統的貝葉斯推斷框架的參數估計和模型選擇的主要難點在于需要分析似然函數,通常假定似然函數為高斯函數,采用Metropolis-Hasting算法迭代分析模型參數的馬爾可夫鏈[23-24]。但該算法存在一些缺點:首先,給定一個模型,其響應的似然函數不一定是高斯函數,對于非線性模型,似然函數類型未知,通常是非高斯函數;其次,采用馬爾可夫鏈蒙特卡洛法估計模型參數時,參數的馬爾可夫鏈收斂性對于狀態轉換函數的類型和參數非常敏感,雖然增加迭代次數可以在一定程度上改善收斂性,但參數的馬爾可夫鏈仍經常出現不能有效收斂的情況。為改善傳統的基于貝葉斯推斷的模型選擇和參數估計,Green等[25]總結了一些改進的貝葉斯推斷方法,Marjoram等[26]提出了近似貝葉斯計算法,該方法無需計算似然函數實現參數估計和模型選擇的方法。

本文結合近似貝葉斯計算法和序貫蒙特卡洛方法[27-28],提出一種包裝件模型選擇和參數估計方法,應用該方法避免了未知類型的似然函數的計算,有效避免了參數估計時初始參數選擇不當、狀態轉移函數參數選擇不當造成的參數選擇效率低下的問題,本文提出的方法可同時實現模型參數估計和模型選擇,具有較高的計算效率和穩健性。

1 貝葉斯推斷

考慮到包裝件模型類型和模型參數的不確定性,在給定的模型類型M和響應數據D的條件下,采用貝葉斯推斷識別包裝件模型參數

(1)

式中:θ為模型的參數向量;p(θ|D,M)為識別出的模型參數的后驗分布函數;p(D|θ,M)為模型的似然函數;p(θ|M)為模型參數的先驗分布;p(D|M)為歸一函數,以確保p(θ|D,M)為概率密度函數。故模型參數的后驗分布由下式確定

p(θ|D,M)∝p(D|θ,M)p(θ|M)

(2)

在貝葉斯推斷框架下,模型的不確定性由不確定的模型參數后驗分布描述,因此,模型在輸入數據x條件下的輸出y是不確定的,模型響應的分布p(y|x,M)可由下式確定

(3)

式中,p(y|x,θ,M)可在給定的模型類型和參數條件下通過數值模擬獲得。

利用式(1)或式(2)確定模型參數時,通常式中的似然函數類型未知,為計算方便,通常假定似然函數為一個高斯函數,但這與實際情況通常不相符,尤其對于非線性模型,其似然函數通常是非高斯的,因此這個近似假定經常造成計算誤差,甚至會導致錯誤結果,近似貝葉斯計算是解決該問題的有效方法,其基本思路是對比模型真實響應和模型預測響應之間的誤差,在此基礎上估計模型參數,可用下式表示

p(θ|D,M)≈pε(θ|D,M)∝p(θ|M)p(Δ(yreal,ysim)<ε|θ)

(4)

式中:ε為閾值;yreal和ysim分別為真實系統響應和模擬的系統響應;函數Δ為定義的系統真實響應和模型模擬響應之間的差異度量函數,該函數可根據實際需要定義,可表征時域響應、時域統計特征參數、頻域特征參數等的差異。在該基本思路的基礎上,人們提出了近似貝葉斯計算拒絕算法[29]和近似貝葉斯計算的馬爾可夫鏈蒙特卡洛算法,但這些算法仍存在著效率低、計算量大、收斂性不佳等問題。

2 基于序貫蒙特卡洛的近似貝葉斯計算

在式(4)中,如果ε足夠小,則p(Δ(yreal,ysim)<ε|θ)可近似認為是模型參數θ的后驗分布,簡寫為p(θ),由于p(θ)未知,無法直接采樣分析,本文采用重要性采樣分析θ的后驗分布p(θ)。我們假定一個建議分布η(θ),對建議分布進行蒙特卡洛采樣

(5)

式中:^為對真實分布的近似估計;N為采樣個數;Θ(i)為服從η分布的采樣點;δ為狄拉克δ函數,定義為

(6)

故后驗分布可由下式確定

(7)

式中,w為重要性權重,由下式定義

(8)

根據以上分析,對于未知的一個分布p(θ),如果不能直接采樣,可對建議分布η(θ)采樣,并分析重要性權重w,根據式(7)可實現對p(θ)的采樣。

由于模型參數θ的后驗分布p(θ)未知,假定的建議分布和真實的后驗分布通常差別較大,直接采用重要性采樣分析θ的后驗分布無法實現。因此,本文結合序貫蒙特卡洛法和重要性采樣求模型參數θ的后驗分布。該方法基本思路為:首先選取一系列數值逐步減小的閾值,εt(t=1,2,…,T),ε1>ε2>…>εT,定義中間分布

pt(θ)=p(θ|M)1(Δ(yreal,ysim,i(θ))≤εt)

(9)

式中:函數1(·)為指示函數,如果函數的變量為真,該函數值為1,否者函數值為0;ysim,i(θ)為選擇第i個參數時模擬的系統響應。根據近似貝葉斯計算,如果εT足夠小,則根據式(9)求得的pT(θ)即為參數的θ后驗分布。因此,利用基于序貫蒙特卡洛法的近似貝葉斯估計分析參數后驗分布的基本思路是通過逐步減小近似貝葉斯估計中的閾值,估計一系列參數的中間分布,當閾值逐步變小時,一系列的中間分布也逐步收斂至參數的后驗分布。本文根據式(7),采用重要性采樣法求參數θ的中間分布,為此,需定義一系列建議分布ηt(θ),由于隨著εt逐步變小,中間分布pt(θ)逐漸收斂至后驗分布pT(θ),根據序貫蒙特卡洛采樣原則,將ηt(θ)定義為對上一次中間分布pt-1(θ)數據擾動后得到的數據的分布

(10)

式中:wt-1由式(8)確定,wt-1(θt-1)=pt-1(θt-1)/ηt-1(θt-1);K為隨機擾動函數,通常采用高斯隨機擾動、均勻隨機擾動、Metropolis-Hasting擾動等。式(10)通常采用蒙特卡洛法近似實現

(11)

(12)

根據以上的推導,為求出參數的后驗分布,本文在參數的先驗分布和后驗分布之間設置T個中間分布,采用近似貝葉斯計算,通過逐漸減小閾值,使參數的先驗分布逐步地收斂至后驗分布。由于中間分布不能直接求出,設置了T個建議分布,將建議分布定義為上一個中間分布的擾動,根據式(12)求得的重要性權重,采用重要性采樣法求出各中間分布。因此,包裝件模型參數的近似貝葉斯計算法可總結如下:

算法1：

步驟1初始化中間分布的閾值ε1,ε2,…,εT,初始化中間分布指示數t=0;

步驟2初始化參數采樣序號i=1,設定參數總采樣數N;

步驟4將擾動后參數θ**(i)代入參數得先驗分布p(θ|M)中,如果p(θ**(i)|M)=0,返回步驟2重新采樣;

步驟5根據參數θ**(i)模擬包裝系統響應ysim,i,如果Δ(yreal,ysim,i)>εt,返回步驟2重新采樣;

步驟6令θt(i)=θ**(i),根據式(10)和式(12)計算重要性權重

(13)

步驟7如果i

以上算法中,θ*表示從上一次中間分布中采樣的參數,θ**表示經擾動后得到的待選參數采樣點。

3 考慮模型選擇的近似貝葉斯計算算法

對于給定的包裝系統,可能有n個類型的模型{M1,M2,…,Mn}可表征其振動響應特性,將m設定為模型編號,m=1,2,…,n,編號為m的模型包含km個參數,表示為:θ(m)=[θ(m)(1),θ(m)(2),…,θ(m)(km)]在這些待選模型中,需分析出能表征包裝件振動特性的最優模型。

為實現模型選擇,我們對算法1改進,首先根據包裝系統振動特性選出n個可能類型的模型,假定這些模型的先驗分布p(m)(通常假定為均勻分布)。在算法1中,我們添加模型指示參數m,從p(m)中選擇m,根據m確定模型類型,迭代分析模型參數后驗分布。具體算法如下:

算法2：

步驟1初始化中間分布的閾值ε1,ε2,…,εT,令t=0;

步驟2令i=1,設定參數總采樣數N;

步驟3從先驗分布p(m)中選擇模型指示參數m*;

步驟4如果t=0,直接從先驗分布p(θ|M(m*))中采樣,直接求得θ**(i),如果t≠0,則從t-1次采樣樣本{θ(m*)t-1}中采樣,采樣值θ*(i)乘以權重w(m*)t-1(i)后根據擾動函數K對采樣值擾動,得擾動后參數θ**(i);

步驟5將擾動后參數θ**(i)代入先驗分布p(θ|M(m*))中,如果p(θ**(i)|M(m*))=0,返回步驟2重新采樣;

步驟6根據參數θ**(i)模擬包裝系統響應ysim,i,如果Δ(yreal,ysim,i)>εt,返回步驟2重新采樣;

步驟7令mt(i)=m*,將θ**(i)添加至集合{θ(m*)t}中,計算重要性權重

(14)

步驟8如果i

給定待選模型Mm(m=1,2,…,n)和振動響應數據yreal,采用算法2可分析出模型指示參數m的后驗分布p(m|yreal),以及模型參數的后驗分布p(θi|yreal,m),i=1,2…,km,因此算法2可同時實現包裝系統的模型選擇和參數識別。在應用算法2時,如果模型mi不能良好表征包裝系統振動響應特性,則滿足Δ(yreal,ysim,i)≤εt的參數個數較少,即p(mi|yreal)較小,隨著εt逐漸減小,p(mi|yreal)變為0,說明模型mi可以剔除。t增至T時,如果在分布p(m|yreal)中僅剩一個類型的模型,該模型即為包裝系統的優選模型。

4 實例分析

圖1 試驗裝置示意圖

由于在包裝件垂向振動過程中,緩沖材料變形量很小,壓縮變形速度變化很有限,故本文不考慮包裝件彈性和阻尼的非線性,采用線性模型、包含有干摩擦的振動模型、Bouc-Wen模型(n=1,2,3,4)作為備選模型,這些模型的運動方程式如下。

線性模型M1:

(15)

式中:m為質量塊的質量;z為振動過程中緩沖材料的變形量,z=x-y;c和k分別為緩沖材料的阻尼系數和彈性系數。

包含有干摩擦的振動模型M2:

(16)

Bouc-Wen模型[30]M3～M6:

(17)

其中

(18)

式中:參數α、β和A為表征緩沖材料遲滯特性的參數,由于在式(17)和式(18)中;k和A均為緩沖材料的彈性系數,故這兩個參數可合并統一考慮。故在式(17)和式(18)所表示的包裝件滯回模型中,包含有線性彈性力、線性阻尼力、非線性滯回力。在試驗中,由于m為已知,故備選模型中的待識別參數可進行以下設定:對于M1,θ=[k,c];M2,θ=[k,c,F];M3～M6,θ=[A,c,α,β]。

本文在進行模型選擇和參數估計時,為確保算法準確實施,應對算法中的參數進行合理設置和選擇:

(1) 誤差函數Δ的選擇和計算。函數Δ有很多種計算方法,如真實數據和模擬數據時域統計矩誤差、頻域參數誤差、加權混合誤差等。本文選擇時域真實加速響應數據和模型模擬的加速響應數據之間的方差和作為誤差函數,采用歸一化均方誤差函數[31]作為誤差函數,定義為

(19)

(20)

采用均勻擾動雖然降低了模型選擇和參數識別的效率,但采用擾動方法有效提高了算法的可靠性。

(3) 中間分布閾值εi(i=1,2,…,T)的選擇。εi可采用人工設定、自適應設定等方法確定,本文首先確定一個較大的ε0值,以確保迭代分析能正常實施,設置系數e=εi/εi-1,為確保在迭代分析過程中能夠逐步分析出模型類型和模型參數,應選擇一個接近1的系數,本文選擇e=0.9。

在試驗過程中隨機采集10 s的振動激勵數據和響應數據,應用算法2進行模型選擇,算法中模型參數范圍和一些初始參數的設定為:k=[2 000,20 000];c=[0,200];F=[0,200];A=[2 000,20 000];α=[0,200];β=[-200,200];ε0=100;N=1 000;T=40。隨著T的增加,各模型后驗分布的變化如圖2所示。迭代過程中,滿足條件的各模型參數點在所有被接受的參數點中所占的比例變化如圖3所示。其中,Pr(Mm)用下式表示

圖2 模型后驗分布變化

圖3 模型概率變化圖

(21)

從圖2和3可知,當迭代次數達到27時,單自由度線性模型被剔除,迭代次數達到29時,包含有干摩擦的線性模型被剔除,迭代次數達到36(此時ε=2.5),在備選模型中,只有Bouc-wen(n=2)模型的誤差滿足近似Bayesian計算的要求,故本文中,根據試驗數據優選出的模型為Bouc-wen(n=2)模型。根據獲得的該模型參數的后驗分布,重新設定模型參數范圍,設置ε0=10,εT=1,N=1 000,采用算法1,經過22次迭代計算,得模型各參數的概率分布如圖4所示。各參數均值和方差如表1所示。

表1 識別出的模型參數

圖4 模型參數直方圖

識別出的模型參數可用于預測包裝件響應,根據圖4和表1中的參數分布,根據式(3),采用蒙特卡洛法模擬包裝件的加速度響應,模擬次數為1 000次,包裝件加速度響應模擬結果和真實記錄的包裝件加速度響應對比如圖5所示。從圖5可知,根據識別出的模型類型和模型參數可準確預測包裝件加速度響應。

圖5 包裝件加速度響應預測和真實加速度響應

5 結 論

(1) 基于重要性采樣和序貫蒙特卡洛,采用近似貝葉斯計算替代傳統的貝葉斯推斷,實現模型選擇和參數識別,根據本文提出的算法,采用試驗數據可同時實現模型類型選擇和參數識別。包裝件響應模擬結果表明,采用本文分析出的模型參數可準確預測包裝件加速度響應。

(2) 本文提出的算法避免了傳統貝葉斯計算中復雜的似然函數的計算,避免了對響應誤差呈高斯分布的假定造成的誤差,提高了分析準確性和計算效率,具有良好的通用性。

(3) 近似貝葉斯計算結果表明,與傳統的線性振動模型和包含有干摩擦的振動模型相比,采用Bouc-Wen模型能夠準確預測包裝件的響應,建議在包裝件模型中采用Bouc-Wen模型。

(4) 本文的模型選擇和參數估計算法是針對單自由度振動系統提出的,從算法中可以看出,該算法不受模型類型的限制,同樣適用于非線性彈性和非線性阻尼振動系統的模型選擇。對于多自由度包裝件,需要構建多自由度運動方程式,同時對包裝件中多個位置進行模型選擇和參數識別,但本文的算法適用每個位置的模型選擇。因此本文算法經改進,適用于非線性和多自由度包裝件模型選擇和參數識別。