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

朱大鵬, 曹興瀟

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

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

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

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

與其它振動結(jié)構(gòu)相似,對包裝件建模時,由于以下因素,導(dǎo)致包裝件模型存在著一些誤差:

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

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

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

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

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

1 貝葉斯推斷

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

(1)

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

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

(2)

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

(3)

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

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

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

(4)

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

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

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

(5)

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

(6)

故后驗分布可由下式確定

(7)

式中,w為重要性權(quán)重,由下式定義

(8)

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

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

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

(9)

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

(10)

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

(11)

(12)

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

算法1：

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

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

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

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

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

(13)

步驟7如果i

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

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

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

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

算法2：

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

步驟2令i=1,設(shè)定參數(shù)總采樣數(shù)N;

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

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

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

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

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

(14)

步驟8如果i

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

4 實例分析

圖1 試驗裝置示意圖

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

線性模型M1:

(15)

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

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

(16)

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

(17)

其中

(18)

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

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

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

(19)

(20)

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

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

在試驗過程中隨機采集10 s的振動激勵數(shù)據(jù)和響應(yīng)數(shù)據(jù),應(yīng)用算法2進行模型選擇,算法中模型參數(shù)范圍和一些初始參數(shù)的設(shè)定為: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所示。迭代過程中,滿足條件的各模型參數(shù)點在所有被接受的參數(shù)點中所占的比例變化如圖3所示。其中,Pr(Mm)用下式表示

圖2 模型后驗分布變化

圖3 模型概率變化圖

(21)

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

表1 識別出的模型參數(shù)

圖4 模型參數(shù)直方圖

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

圖5 包裝件加速度響應(yīng)預(yù)測和真實加速度響應(yīng)

5 結(jié) 論

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

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

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

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