基于近似贝叶斯计算的包装件模型选择和参数估计

朱大鹏, 曹兴潇

(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) 本文的模型选择和参数估计算法是针对单自由度振动系统提出的,从算法中可以看出,该算法不受模型类型的限制,同样适用于非线性弹性和非线性阻尼振动系统的模型选择。对于多自由度包装件,需要构建多自由度运动方程式,同时对包装件中多个位置进行模型选择和参数识别,但本文的算法适用每个位置的模型选择。因此本文算法经改进,适用于非线性和多自由度包装件模型选择和参数识别。