曹汝男, 孙志礼, 张毅博, 王 健

(东北大学 机械工程与自动化学院, 辽宁 沈阳 110819)

近年来,实际工程中机械结构的可靠性问题备受关注,为此国内外学者提出了多种可靠性分析方法.一次一阶矩法(FORM)和二次二阶矩法(SORM)将非线性功能函数线性化,通过计算近似得到功能函数的失效概率[1-2];但FORM和SORM适用范围有限,只适用于显式功能函数,且对于高维非线性问题精度较差.Monte Carlo方法[3]将求解可靠性的多维积分问题转化为数学期望的方式,通过大量抽样计算失效概率,但因为计算量大、效率低,多用于检验新方法的准确性.20世纪90年代众多学者提出了代理模型法[4-5],如多项式响应面法[6]、神经网络法[7]、支持向量机[8],以及Kriging模型法[9].代理模型法通过插值或者回归方法构建隐式功能函数的替代模型,随后对该替代模型进行失效概率的计算.Kriging模型是一种高效的插值方法,以最小方差无偏估计保证插值精度,适用于求解非线性隐式函数问题.Echard等[10]将Kriging和Monte Carlo方法相结合,提出了一种具有主动学习功能的可靠性分析方法,即AK-MCS法.孙志礼等[11]考虑到联合概率密度函数和Kriging标准差对失效概率精度的共同影响,提出了新的学习函数VF以及相应的学习停止条件.

AK-MCS等传统方法在每次迭代中只能选取一个最佳样本点,即在每次更新Kriging模型时只能使用一台电脑进行一次仿真分析;这样不但效率低,也不能充分利用实验室的资源.为了在保证计算精度的前提下,进一步提高计算效率,本文基于自适应学习函数VF对AK-MCS法进行改进,结合k-means聚类分析方法,在每次迭代中将Kriging模型及候补样本点分组,选取多个样本点从多方位提高Kriging模型精度,减少学习过程的迭代次数;同时实现并行计算,降低仿真总时间,提高计算速度.将改进方法(AK-MCS-K方法)应用于某一类型火炮协调器中,对火炮协调器的运动精度可靠性进行分析.

1 理论基础

1.1 Monte Carlo法

假设极限状态方程将整个区域分为两个部分,即失效域F和安全域S:失效域F={x|G(x)≤0},反之系统处在安全域内S={x|G(x)>0};x是随机变量,G(x)是功能函数的响应值.失效概率Pf可以表示为

(1)

式中f(x)是随机变量的联合密度函数.

使用Monte Carlo法,首先要抽取NMC个随机样本点xj(j=1,2,…,NMC),然后将样本点代入功能函数中得到对应的响应值G(xj).统计其中G(xj)≤0的数量,记为Nf.Nf与NMC的比值即为失效概率近似值.

(2)

1.2 Kriging方法

Kriging由参数化模型和随机过程两部分组成.Kriging模型中极限状态函数为

(3)

式中:βh是回归系数;gh(x)表示变量x的多项式;z(x)服从高斯分布N(0,σ2),其协方差可表示为

Cov[z(xi),z(xj)]=σ2R(xi,xj;θ) .

(4)

式中:σ2是高斯过程的方差;R(xi,xj;θ)是xi和xj的相关函数;θ是相关函数的参数.

已知样本集SDoE=[x1,x2,…,xN]及对应的响应值Y=[y1,y2,…,yN]T,G(x)的最小方差无偏估计为

(5)

rT(x)R-1r(x)] .

(6)

式中:

r(x)=[R(x1,x;θ),…,R(xN,x;θ)]T;

u(x)=GTR-1r(x)-g(x);

R=(R(xi,xj;θ))N×N;

G=[g(x1),g(x2),…,g(xN)]T.

2 AK-MCS-K方法

2.1 自适应学习函数VF和学习停止条件

(7)

针对样本点的选取,AK-MCS方法提出了学习函数U：

U(x)=|μG(x)|/σG(x) .

(8)

选取min(U(x))的点更新Kriging模型可以有效地提高模型精度.

VF(x)=σG·f(x) .

(9)

每次迭代选取VF(x)最大的点,调用真实功能函数得到对应的响应值,用以更新Kriging模型,直到精度满足要求.学习函数VF对应的学习停止条件为

eε=Nu/Nf≤[e] .

(10)

式中:Nu为符号估计错误点数的期望;Nf为失效点的个数;eε为失效概率误差的期望值；[e]为eε的阈值,根据工程需要选取.

2.2 k-means算法

k-means算法[12]是由Mac Que提出的一种应用广泛的经典聚类分析方法.

k-means算法计算步骤如下:

①随机抽取k个数据点作为聚类中心;

②计算每个点到聚类中心的距离,按照距离最小原则进行分类;

③计算k个组的中心点所在位置;

④如果前后两次中心点位置相同,则输出结果,否则返回步骤②.

2.3 AK-MCS-K方法原理

为了进一步提升计算效率,本文对基于学习函数VF的主动学习方法进行改进.为了减少迭代次数,拟在每一次迭代中选取多个最佳样本点,同时加入多个样本点对Kriging模型进行拟合.如图1所示,以每次取4个样本点为例,虚线表示真实的功能函数,实线表示拟合曲线,叉号表示候选样本点集,六角星表示选择的最佳样本点.如果直接按照使VF取最大值的标准选点,图中选取的4个样本点中,有3个点距离过于接近,并不能充分有效地提高Kriging模型的精度.因此,引入k-means算法,先对候补样本点集进行分类,将其分为4部分,然后在每一部分中选取满足学习函数VF的最佳样本点,结果如图2所示.

应用k-means方法可以成功避免选取样本点过密的情况,并且可以在多个位置同时提高Kriging模型的精度.

2.4 AK-MCS-K方法计算流程

步骤1 应用拉丁超立方抽样选取N0个初始样本点xDoE=[x1,x2,…,xN0];

步骤2 计算相应的响应值y=[y1,y2,…,yN0],并使用DACE工具箱建立Kriging模型;

步骤3 使用Monte Carlo方法建立样本空间S=[x1,x2,…,xMC];

步骤5 使用聚类分析方法将候补样本点进行分组,共分为k个组;

步骤6 使用学习函数VF在每一组中选取VF(x)值最大的点作为最佳样本点;

步骤7 将选取的k个最佳样本点加入到样本空间中并调用真实功能函数计算其响应值,重新构建Kriging模型;

步骤8 根据学习停止条件式(10)判断Kriging模型精度是否满足要求,如果不满足则返回步骤4;

值得注意的是,在使用Monte Carlo 法计算失效概率时,应满足式(11),其中[δ]为阈值,可根据精度需要进行取值.

(11)

2.5 算例分析:多维非线性系统

为验证AK-MCS-K方法的可行性和优势所在,选用一个非线性无阻尼单自由度系统[13-14]为研究对象进行可靠性分析.如图3所示,k1,k2为弹簧刚度,m为小车质量.

其功能函数为

(12)

式中R为给定的阈值,6个变量相互独立,满足表1所示正态分布.

表1 随机变量分布Table 1 Distribution of random variables

图4和图5分别表示失效概率估计值和停止条件eε随迭代次数的变化曲线.图中4条曲线分别代表AK-MCS-EFF,AK-MCS-U,AK-MCS-VF和AK-MCS-K(k=2)4种方法.可以发现,AK-MCS-K方法得到的失效概率估计值可以很好地收敛在标准值附近,并且其收敛速度快于其他3种方法.

由表2可见,对于AK-MCS-K方法,随着k值的增加,其对应的迭代次数逐渐降低,计算效率得到提高,但并行计算使用的电脑也会增加,所以应根据实际情况选取适当k值.

表2 算例结果统计表Table 2 Statistical table of the results of an example

3 火炮协调器定位精度可靠性分析

3.1 火炮协调器简介

协调器的结构如图6所示,主要包括耳轴、减速箱、协调支臂、液压油缸和摆弹油缸,以及平衡机等.在工作时,协调器接住供弹仓推出的炮弹,协调到指定角度,然后在翻转油缸的作用下将托弹盘连同炮弹快速地翻转到输弹线,将炮弹输送到炮膛.判定协调器是否失效的一个重要依据就是能否精确地将炮弹送到炮膛.本文针对协调器对炮弹协调过程的定位精度进行可靠性分析.

3.2 协调器刚柔耦合模型可靠性分析

对于实际工程中的可靠性分析问题,为节约成本,方便操作,多使用虚拟样机技术代替物理样机进行研究.协调器刚柔耦合建模过程如图7所示.

使用ADAMS软件对协调器协调过程进行动力学仿真,得到炮弹质心位移曲线.如图8所示,协调器摆弹臂的整个协调运动过程共1.2 s,只需采集t=1.2 s时刻炮弹质心位移,并以此为标准值ξ.

建立定位精度可靠性分析的状态函数:

G(X)=γ-|ξ-ξ(r1,r2,r3,r4,lx,ly)| .

(13)

式中:γ为允许的误差范围,可以根据精度要求适当取值,此处取6.0 mm;r1～r4表示孔半径;lx和ly表示协调臂x和y方向上的制造误差;ξ为不考虑原始误差的摆弹臂质心位移；ξ(r1,r2,r3,r4,lx,ly)表示考虑原始误差后得到的摆弹臂质心位移,其结果通过仿真得到.

式(13)中G(X)与6个随机变量间的函数关系无法用公式明确地表达出来,只能通过仿真得到,因此状态函数属于多维隐式函数.各随机变量分布如表3所示.

将ADAMS软件建立的协调器刚柔耦合参数化模型,转换成.cmd文件格式,搭建Simulink-ADAMS联合仿真平台用于可靠性分析.采用AK-MCS-K方法进行计算,其中初始样本点数为15,k值取4,[e]=0.05.结果如图9和图10所示.

表3 随机变量分布类型及参数
Table 3Distribution types and parameters of randomvariables

mm

从图9和图10中的分析结果可以得知,在迭代次数约为55次,对应的样本点数约为15+55×4=235时,失效概率已经收敛,同时失效概率误差的期望值eε降到0.05以下,满足了工程中的精度要求,并最终得到失效概率为1.83×10-3.值得注意的是,调用ADAMS软件进行1次协调过程的仿真时间约为26 min,虽然其样本点数为235,但是仅迭代了55次,由于使用多台计算机并行运算,仿真所用时间也仅为55×26(min)=23.8(h).因此,将AK-MCS-K用于求解仿真耗时的可靠性问题可以有效地提高计算效率,降低时间成本.

4 结 论

1) AK-MCS-K方法在保证精度的条件下对原始AK-MCS方法每次迭代选取一个最佳样本点的不足进行改进,使得每次迭代可以选取多个最佳样本点,从而大大减少迭代次数.

2) 在选取样本点过程中,采用k-means方法将候补样本点分成k个组,同时将Kriging模型分为k个部分,每部分选取一个最佳样本点,成功避免了选取样本点过密的情况,可以有效地提高Kriging模型的精度.

3) 在最佳样本点的仿真过程中,采用并行运算,多台电脑同时进行分析,一次仿真的时间,可以完成k个最佳样本点对应的响应值的仿真,有效地提高了计算效率.