基于貝葉斯推理的包裝件動力學模型優(yōu)化選擇研究

朱大鵬，余珍，曹興瀟

裝備防護

基于貝葉斯推理的包裝件動力學模型優(yōu)化選擇研究

朱大鵬a，余珍a，曹興瀟b

（蘭州交通大學 a.交通運輸學院 b.機電工程學院，蘭州 730070）

在多種類型的模型中挑選出最優(yōu)包裝件模型，并實現(xiàn)參數(shù)識別的方法。文中將包裝件模型構(gòu)建為參數(shù)不確定模型，在貝葉斯推理的框架下，采用馬爾可夫鏈蒙特卡洛法識別模型參數(shù)，采用偏差信息準則（DIC）計算各備選模型的DIC參數(shù)，選擇出最優(yōu)包裝件模型。在振動實驗臺用質(zhì)量塊–緩沖材料模擬包裝件并進行隨機振動測試，分析結(jié)果表明，Bouc–Wen（=2）模型為文中包裝系統(tǒng)的最佳模型。文中提出的基于貝葉斯推理的包裝件模型優(yōu)化選擇和參數(shù)識別方法考慮了模型不確定性，構(gòu)建的模型可準確預(yù)測包裝件在隨機振動下加速度響應(yīng)的時域信號。

模型選擇；貝葉斯推斷；馬爾可夫鏈蒙特卡洛；Metropolis–Hastings算法；偏差信息準則

包裝件在流通過程中，由于長時間受到振動載荷的作用，包裝件中的產(chǎn)品在振動載荷作用下可能產(chǎn)生損壞。為降低產(chǎn)品損壞概率，需構(gòu)建包裝件模型，分析包裝件在振動載荷下的響應(yīng)，找出影響包裝件響應(yīng)的因素，提出優(yōu)化包裝設(shè)計的方案。構(gòu)建準確的包裝件模型，是準確分析包裝件響應(yīng)、判斷包裝件可靠性和安全性的重要基礎(chǔ)。

在振動載荷下，如何準確構(gòu)建結(jié)構(gòu)模型并進行模型參數(shù)識別，在結(jié)構(gòu)響應(yīng)預(yù)測、結(jié)構(gòu)狀態(tài)監(jiān)控、故障診斷等領(lǐng)域，是一個關(guān)鍵問題，近年來該問題受到了國內(nèi)外研究者的廣泛關(guān)注。De等[1]考慮了結(jié)構(gòu)非線性和參數(shù)不確定，將結(jié)構(gòu)響應(yīng)劃分為兩部分：系統(tǒng)的名義線性部分和系統(tǒng)修正部分，其中修正部分的響應(yīng)用神經(jīng)網(wǎng)絡(luò)建模。Scheel等[2]提出了在諧波激勵和隨機激勵下非線性系統(tǒng)識別方法，該模型可準確模擬非線性系統(tǒng)在共振狀態(tài)的響應(yīng)。Esfahani等[3]提出了一種識別Bouc–Wen模型的黑箱模型，采用多項式解耦對模型進行了簡化。Chen等[4]和Jin等[5]構(gòu)建實驗系統(tǒng)，結(jié)合傳統(tǒng)線性系統(tǒng)測試數(shù)據(jù)和基于高速攝像系統(tǒng)的圖像處理結(jié)果，提出了非線性結(jié)構(gòu)識別方法。Civera等[6-7]將電路識別領(lǐng)域中的向量擬合法進行改進，提出用于結(jié)構(gòu)特性識別的快速松弛向量擬合法，該方法已成功用于工程結(jié)構(gòu)識別和建模中。朱大鵬等[8-9]將包裝系統(tǒng)簡化為單自由度振動系統(tǒng)，提出了隨機振動載荷和沖擊載荷條件下系統(tǒng)參數(shù)識別方法。由于可以采用多種模型描述一個給定的振動系統(tǒng)，在這些模型中，需要優(yōu)化選擇一種最佳模型，減小模型響應(yīng)誤差。Fuentes等[10]在單自由度系統(tǒng)中構(gòu)建各種可能的非線性彈性和阻尼函數(shù)庫，在稀疏約束條件下，采用標準Lasso回歸選擇單自由度模型中最優(yōu)的彈性和阻尼函數(shù)類型。Safari等[11]系統(tǒng)對比了目前常用的非線性系統(tǒng)識別方法，實例分析表明，混合式方法可準確識別出復(fù)雜的系統(tǒng)彈性和阻尼特性。

文中將包裝件模型簡化為單自由度系統(tǒng)，考慮到模型的不準確性、實驗數(shù)據(jù)受噪聲污染等因素，將模型構(gòu)建為參數(shù)不確定的模型，應(yīng)用貝葉斯推斷法推斷識別模型的不確定參數(shù)。包裝件模型可用多種模型類型表征，文中采用DIC準則（Deviance Information Criterion）對這些模型進行對比，比選出最優(yōu)包裝件模型。

1 包裝件模型的不確定性

將包裝件簡化為單自由度支座激勵系統(tǒng)，對該系統(tǒng)進行實驗測試，記錄輸入和輸出數(shù)據(jù)，并識別包裝件模型參數(shù)。由于以下原因，包裝件模型中存在著不確定性。

1）包裝件模型存在著誤差。雖然目前可采用各種線性和非線性模型對包裝件中的彈性特性和阻尼特性建模，但大部分的緩沖材料在載荷作用下變形和能耗機理尚不清晰，其本構(gòu)模型尚未完全建立，故目前采用的模型與真實的緩沖材料特性之間存在著誤差。

2）實驗過程中噪聲污染。在實驗數(shù)據(jù)記錄過程中不可避免會受到噪聲影響，這也會影響模型參數(shù)識別。

3）記錄的數(shù)據(jù)有限。實驗過程中只能記錄有限長度的數(shù)據(jù)，用這些數(shù)據(jù)表征包裝件特性時也會造成一定的誤差。

考慮到上述模型中存在的不確定性，文中將包裝件模型看作是具有不確定參數(shù)的振動模型，根據(jù)實驗數(shù)據(jù)識別出的模型參數(shù)不是確定的值，而是具有一定分布的不確定參數(shù)：

式中：為包裝件模型中的參數(shù)向量，對于線性彈性–阻尼模型，=[,]，對于Bouc–Wen模型，=[,,,]；為采集的實驗數(shù)據(jù)；為選擇的包裝件模型。由于包裝件模型的不確定性，應(yīng)用該模型預(yù)測包裝件響應(yīng)時，也具有不確定性：

式中：為包裝件受到的激勵。的分布可由式（3）確定。

2 包裝件模型的識別

對于不確定的包裝件模型，可采用Bayesian推斷法確定其模型參數(shù)：

式中：(|)稱為模型的似然函數(shù)；()為模型參數(shù)的先驗分布，分母()由式（5）確定。

式中，分母()是一個歸一化函數(shù)，確保式（4）獲得的函數(shù)是一個概率密度函數(shù)，故式（4）可簡化為式（6）。

可用式（6）估計包裝件模型參數(shù)分布，再進行歸一化處理。先驗分布函數(shù)()通常根據(jù)經(jīng)驗確定。對于似然函數(shù)，假定模型預(yù)測的包裝件響應(yīng)(t,)和真實記錄的包裝件響應(yīng)y之間的誤差為ε，即：

其中，為采樣點個數(shù)。假定式（7）中各采樣點誤差互相獨立且呈正態(tài)分布，則似然函數(shù)由式（8）確定。

其中，為模型響應(yīng)和真實響應(yīng)誤差過程的標準差。直接根據(jù)式（6）求出包裝件模型參數(shù)是很困難的，在Bayesian推斷框架下，可采用馬爾可夫鏈蒙特卡洛（Markov Chain Monte Carlo, MCMC）法[12]識別參數(shù)，該方法的基本原理是對參數(shù)構(gòu)建一個平穩(wěn)分布的馬爾可夫鏈，通過合理控制的變化，該馬爾可夫鏈的分布最終收斂于(|)，即應(yīng)用馬爾可夫鏈生成一系列參數(shù)[0,1,2,…]，隨著生成的參數(shù)數(shù)量的增加，可確保參數(shù)符合(|)分布。文中采用Metropolis–Hastings算法[13]設(shè)置參數(shù)馬爾可夫鏈，具體步驟如下：

1）設(shè)置馬爾可夫鏈總長度，令=1，估計初始參數(shù)值0，令(t)=0。

2）令=+1，根據(jù)狀態(tài)轉(zhuǎn)移函數(shù)(*|)生成*。

3）根據(jù)式（8）構(gòu)建的似然函數(shù)計算：

4）如果≥1，則接受*為新的參數(shù)，否則以概率接受*為新的參數(shù)；

5）如果*被接受為新的參數(shù)，則(t+1)=*，否則(t+1)=。

6）重復(fù)步驟2—5，直至=。

通過合理選擇和狀態(tài)轉(zhuǎn)移函數(shù)中的參數(shù)，參數(shù)馬爾可夫鏈可較快收斂至真實參數(shù)范圍，可用收斂后的馬爾可夫鏈近似估計真實的系統(tǒng)參數(shù)。

3 識別最優(yōu)模型

對于一個給定的包裝系統(tǒng)，可采用多種類型的模型分析該系統(tǒng)特性，在每個模型類型框架下，根據(jù)實驗數(shù)據(jù)，可識別出每種模型的參數(shù)。如，對于單自由度包裝件模型，可采用線性振動系統(tǒng)、非線性彈性振動系統(tǒng)、非線性阻尼振動系統(tǒng)、Bouc–Wen系統(tǒng)等對包裝件進行參數(shù)識別和響應(yīng)預(yù)測，因此，需要對比各類模型，優(yōu)選出最佳模型，用于對包裝件進行響應(yīng)預(yù)測。

為優(yōu)選模型，定義模型偏差參數(shù)()[14]：

()=?2log(|)+(10)

式中：為常數(shù)。似然函數(shù)(|)由式（8）確定，很顯然，模型偏差參數(shù)()越小，模型準確性越好，在對各模型進行對比時，常數(shù)可不予考慮。由于對模型參數(shù)識別時，為隨機變量，故()也具有隨機性，()的均值可通過蒙特卡洛法確定：

式中：為采樣個數(shù)；為采樣的參數(shù)值；[ ]為求期望值運算。定義模型參數(shù)有效數(shù)：

4 實例分析

對于線性振動模型，包裝系統(tǒng)運動方程式為：

式中：m，c和k分別為質(zhì)量、阻尼系數(shù)和彈性系數(shù)，由于實際中m值很容易獲得，故線性振動模型中，需識別的參數(shù)為θ=[k,c]。

對于Bouc–Wen模型[17]，包裝系統(tǒng)運動方程為：

其中，

在式（16）和式（17）中，參數(shù)和均代表材料彈性參數(shù)，在式（16）中可將和統(tǒng)一考慮，故在Bouc–Wen模型中需識別的參數(shù)為=[,,,]。

1）根據(jù)線性振動模型和Bouc–Wen模型基本特性，預(yù)估模型初始參數(shù)0。

2）根據(jù)0采用龍格庫塔法模擬系統(tǒng)加速度響應(yīng)(t,0)。

3）計算模擬響應(yīng)和真實響應(yīng)之間的誤差，根據(jù)式（8）計算似然函數(shù)。

4）在Metropolis–Hastings算法中重復(fù)步驟2和3，迭代分析求出參數(shù)的馬爾可夫鏈數(shù)據(jù)。

在以上算法步驟中，需注意以下事項：

1）在Metropolis–Hastings算法中，狀態(tài)轉(zhuǎn)移方程采用高斯函數(shù)，即馬爾可夫過程的下一個狀態(tài)根據(jù)一個高斯函數(shù)確定，其均值為當前狀態(tài)，即：(*|)～(,)，其中( )表示正態(tài)分布。文中確定各具體參數(shù)下一狀態(tài)時，需合理選擇正態(tài)分布中的值，該值選擇過大或過小都會造成馬爾可夫過程收斂速度慢的問題。

2）參數(shù)初始值0應(yīng)合理選擇，以提高參數(shù)馬爾可夫過程的收斂速度。

3）在計算能力允許的前提下，盡可能增加Metropolis–Hastings算法的迭代次數(shù)，即使參數(shù)初始值0的選擇和狀態(tài)轉(zhuǎn)移方程中的選擇存在誤差，隨著的增加，也可收斂至真實參數(shù)值。

根據(jù)以上算法步驟，設(shè)定迭代運算10 000次，分別識別線性振動模型和Bouc–Wen模型（=1,2,3,4）參數(shù)的馬爾可夫過程，去除燃燒期（為準確起見，文中選擇前4 000個數(shù)據(jù)作為燃燒期），各模型參數(shù)分布情況如表1所示。從表1結(jié)果可以看出，考慮到模型的不確定性，可將包裝件用參數(shù)不確定的模型表示，相應(yīng)地，在給定振動載荷下，包裝件的響應(yīng)也是不確定的值，是在一定范圍內(nèi)波動的。對于一個給定的包裝件，由于存在多個可描述其振動特性的模型，文中采用DIC準則優(yōu)選出最佳模型，根據(jù)識別出的模型參數(shù)計算各模型響應(yīng)，對比包裝件真實響應(yīng)分析各模型的響應(yīng)誤差，采用式（11）—（13）計算各模型的DIC值，如表1所示。從表1中可以看出，Bouc–Wen模型（=2）的DIC值最小，用該模型是這些模型中的最佳包裝件模型。文中識別出的Bouc–Wen模型（=2）參數(shù)分布如圖2所示。根據(jù)圖2的參數(shù)識別結(jié)果，可將各識別的參數(shù)結(jié)果近似看作正態(tài)分布，識別出各參數(shù)的均值和方差，如表1所示，這些參數(shù)識別結(jié)果可作為模型響應(yīng)預(yù)測的依據(jù)。

在得到包裝件的最優(yōu)模型及模型參數(shù)后，可用式（3）預(yù)測該不確定模型的響應(yīng)，首先對包裝件模型中的不確定參數(shù)離散化，將式（3）的積分處理為不確定離散參數(shù)下的模擬，模擬的包裝件響應(yīng)和實驗記錄的真實響應(yīng)如圖3所示，實線為實驗記錄的質(zhì)量塊的加速度響應(yīng)，虛線為對Bouc–Wen模型（=2）參數(shù)離散化處理后預(yù)測的包裝件響應(yīng)的包絡(luò)線。從圖3可以看出，利用該模型預(yù)測的包裝件響應(yīng)是一個不確定的振動響應(yīng)，該模型可準確預(yù)測包裝件的響應(yīng)，且考慮了包裝件模型中的一些不確定性。根據(jù)識別出的模型類型及參數(shù)，可對包裝件在運輸中的安全性和可靠性進行評價分析。

表1 各模型參數(shù)及DIC值

Tab.1 Parameters and DIC values of each model

圖2 Bouc–Wen模型（n=2）參數(shù)分布

圖3 模型預(yù)測的包裝件響應(yīng)和真實響應(yīng)對比

4 結(jié)語

文中考慮了包裝件模型中的不確定性，在貝葉斯推斷框架下構(gòu)建包裝件模型，采用馬爾可夫鏈蒙特卡洛法識別各備選模型參數(shù)，采用偏差信息準則優(yōu)選出包裝件最優(yōu)模型。文中得出主要結(jié)論如下：

1）考慮到模型的不確定性，在貝葉斯推斷框架下，沒有一個模型能夠精確無誤地描述包裝件動態(tài)特性，只能在一系列備選模型中選擇出最近似的模型。在文中的包裝系統(tǒng)中，Bouc–Wen模型能夠更好地描述包裝件的振動特性，該模型考慮了模型中的不確定因素，能夠準確預(yù)測包裝件的響應(yīng)。

2）文中提出的基于Metropolis–Hastings算法的馬爾可夫鏈蒙特卡洛參數(shù)識別方法需對參數(shù)初始值和狀態(tài)轉(zhuǎn)移函數(shù)中各參數(shù)的方差合理選擇，以提高參數(shù)收斂速度，減小計算量。

3）偏差信息準則綜合考慮了模型的準確性和復(fù)雜性，是優(yōu)選最優(yōu)模型的有效方法，該算法具有計算效率高、普適性好的特點，該方法不僅適用于包裝模型的選擇，在結(jié)構(gòu)、電氣、機械等領(lǐng)域也可廣泛應(yīng)用。

4）在備選模型中只選用了線性模型和Bouc–Wen模型，在今后研究中，可根據(jù)包裝件響應(yīng)的時域和頻域特性，在備選模型中增加其他類型的模型進行優(yōu)選，以提高模型選擇的準確性。

[1] DE S. Uncertainty Quantification of Locally Nonlinear Dynamical Systems Using Neural Networks[J]. Journal of Computing in Civil Engineering, 2021, 35(4): 04021009.

[2] SCHEEL M, KLEYMAN G, TATAR A, et al. Experimental Assessment of Polynomial Nonlinear State-Space and Nonlinear-Mode Models for Near-Resonant Vibrations[EB/OL]. 2020: arXiv: 2011.08520. https://arxiv. org/abs/2011.08520

[3] ESFAHANI A F, DREESEN P, TIELS K, et al. Parameter Reduction in Nonlinear State-space Identification of Hysteresis[J]. Mechanical Systems and Signal Processing, 2018, 104: 884-895.

[4] CHEN W, JANA D, SINGH A, et al. Measurement and Identification of the Nonlinear Dynamics of a Jointed Structure Using Full-Field Data, Part I: Measurement of Nonlinear Dynamics[J]. Mechanical Systems and Signal Processing, 2022, 166: 108401.

[5] JIN M, KOSOVA G, CENEDESE M, et al. Measurement and Identification of the Nonlinear Dynamics of a Jointed Structure Using Full-Field Data; Part Ⅱ-Nonlinear System Identification[J]. Mechanical Systems and Signal Processing, 2022, 166: 108402.

[6] CIVERA M, CALAMAI G, ZANOTTI FRAGONARA L. Experimental Modal Analysis of Structural Systems by Using the Fast Relaxed Vector Fitting Method[J]. Structural Control and Health Monitoring, 2021, 28(4): 2695.

[7] CIVERA M, CALAMAI G, FRAGONARA L Z. System Identification Via Fast Relaxed Vector Fitting for the Structural Health Monitoring of Masonry Bridges[J]. Structures, 2021, 30: 277-293.

[8] 朱大鵬. 隨機振動條件下包裝件動態(tài)特性識別方法[J]. 包裝工程, 2017, 38(17): 96-101.

ZHU Da-peng. Dynamic Property Identification Method for Packages under Random Vibration[J]. Packaging Engineering, 2017, 38(17): 96-101.

[9] 朱大鵬. 包裝系統(tǒng)非線性特性識別方法[J]. 包裝工程, 2016, 37(21): 46-51.

ZHU Da-peng. Identification Method of Nonlinear Properties in Packaging System[J]. Packaging Engineering, 2016, 37(21): 46-51.

[10] FUENTES R, NAYEK R, GARDNER P, et al. Equation Discovery for Nonlinear Dynamical Systems: a Bayesian Viewpoint[J]. Mechanical Systems and Signal Processing, 2021, 154: 107528.

[11] SAFARI S, MONSALVE J L. Benchmarking Optimisation Methods for Model Selection and Parameter Estimation of Nonlinear Systems[J]. Vibration, 2021, 4(3): 648-665.

[12] GREEN P L, WORDEN K. Bayesian and Markov Chain Monte Carlo Methods for Identifying Nonlinear Systems in the Presence of Uncertainty[J]. Philosophical Transactions Series A, Mathematical, Physical, and Engineering Sciences, 2015, 373(2051): 20140405.

[13] CHIB S, GREENBERG E. Understanding the Metropolis-Hastings Algorithm[J]. The American Statistician, 1995, 49(4): 327-335.

[14] SPIEGELHALTER D J, BEST N G, CARLIN B P, et al. Bayesian Measures of Model Complexity and Fit[J]. Journal of the Royal Statistical Society: Series b, 2002, 64(4): 583-639.

[15] SPIEGELHALTER D J, BEST N G, CARLIN B P, et al. The Deviance Information Criterion: 12 Years on[J]. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 2014, 76(3): 485-493.

[16] WORDEN K. Data Processing and Experiment Design for the Restoring Force Surface Method, Part I: Integration and Differentiation of Measured Time Data[J]. Mechanical Systems and Signal Processing, 1990, 4(4): 295-319.

[17] WEN Y K. Method for Random Vibration of Hysteretic Systems[J]. Journal of the Engineering Mechanics Division, 1976, 102(2): 249-263.

Kinetic Package Model Optimal Selection Using Bayesian Inference

ZHU Da-penga, YU Zhena, CAO Xing-xiaob

(a. School of Traffic and Transportation, b. School of Mechanical Engineering, Lanzhou Jiaotong University, Lanzhou 730070, China)

The work aims to select the optimal model from various models and achieve parameter identification. The package model was constructed as the parameter uncertainty model. The model parameters were identified with the Markov Chain Monte Carlo method in Bayesian inference framework. Deviation information criterion (DIC) was used to calculate the DIC parameters of each alternative model and select the optimal package model. The mass block-buffer material was used to simulate the package and the random vibration test was carried out on the vibration test bench. The analysis result showed that the Bouc-Wen model (=2) was the optimal model for the package. The model selection and parameter identification method proposed based on Bayesian inference takes the uncertainty of the model into account. The constructed model can accurately predict the time-domain signal of acceleration response of the package under random vibration.

model selection; Bayesian inference; Markov Chain Monte-Carlo (MCMC); Metropolis-Hastings algorithm; deviance information criterion (DIC)

TB485.3

A

1001-3563(2023)05-0238-06

10.19554/j.cnki.1001-3563.2023.05.030

2022–10–10

甘肅省自然科學基金（20JR5RA400）；蘭州交通大學–天津大學聯(lián)合創(chuàng)新基金項目（2022063）

朱大鵬（1977—），男，博士，教授，主要研究方向為運輸包裝、運輸安全。

責任編輯：曾鈺嬋