基于蒙特卡罗方法的机床主轴可靠性预测*

刘智键，刘宝林，胡远彪

(1.中国地质大学(北京) 工程技术学院， 北京 100083；2.国土资源部深部地质钻探技术重点实验室，北京 100083)

0 引言

随着智能制造技术的不断发展，机械加工技术和加工精度要求也越来越高，特别是现代智能化数控机床加工，由于实现了智能化技术，对数控机床主轴的可靠性提高了要求。

机床主轴可靠性预测是数控机床整体系统可靠性预测的重要环节之一，主轴可靠性与否对机床整体加工影响特别大，对提高生产效率、降低生产成本具有重要的意义[1]。因此，在维修机床故障过程中，机床主轴可靠性预测方法的研究显得尤为重要。首先，吉林大学机床实验室等对机床主轴进行寿命研究，以加速机床主轴寿命试验方法为基础，为主轴可靠性研究提供数据参考[2]；西安理工大学机床研究课题组提出了一种基于性能可靠性研究的方法，且应用到机床主轴可靠性能的研究[3-4]。以上两者对机床主轴可靠性分析存在两点不足之处：①加速机床寿命试验方法时间周期长，计算精度不高；②实验采集的数据样本量小，不足以验证方法的准确性。

针对以上不足之处，本文提出了一种基于蒙特卡罗模拟仿真的机床主轴可靠度求解方法。该方法借助于计算机技术建立主轴可靠性分析模型，将采集的大量样本数据传输计算机软件进行计算分析，实现了大样本数据简单、快速计算，节省了繁复的数学推理和演算过程，计算结果更加准确、可靠、通俗易懂，从而达到一种快速有效的检测效果。

1 蒙特卡罗理论

蒙特卡罗方法理论(简称MC)在很多方面得到了应用，例如数学、工程、以及经济数等等[5]，主要解决无法采用确定性模型解决实际问题的一种途径。它的主要思想是大量生成符合样本总体分布的随机数，是一种借助于计算机的统计学的模拟方法。

蒙特卡罗的理论是建立在统计学知识之上，假设P(T)是随机事件T发生的概率，对于该事件在N次实验过程中发生次数M次，此时事件T出现的频率是W(A)=M/N，则存在任意实数ε>0满足：

(1)

从公式(1)中可以看出，当N足够大时，事件T出现的概率为P(T)。此外，假设独立同分布变量是r，随机变量的期望值为E(rk)=μ，其方差为：

D(rk)=σ2(σ2≠0,k=1,2,3…)，将变量标准化，从而得到新的随机变量如下：

(2)

使其满足正态分布的性质，利用正态分布[6]的特点进行求解，新的随机变量的分布函数满足：

(3)

(1) 从而得到抽样样本r的分布函数近似服从正态分布:

(4)

(2) 抽样样本r的期望和方差随机变量公式的概率近似公式:

(5)

大数定律和中心极限定理[6]从理论上保证了蒙特卡罗方法理论的准确性，由大数定律可以看出蒙特卡罗随机产生的随机变量具有稳定性和收敛性，由中心极限定律可以分析蒙特卡罗方法产生误差等问题。

2 基于蒙特卡罗方法的逻辑思路

蒙特卡罗方法主要应用在工程、数学、经济学以及物理学上，随着计算机技术的不断发展，蒙特卡罗方法应用更加广泛。蒙特卡罗的求解问题思路是根据求解问题构造一个数学模型(概率统计模型)，将问题转化为求解模型的参数和特征向量，利用蒙特卡罗模拟求解模型的参数和特征值得出估计值[7]。图1是蒙特卡罗求解问题的主要思路图。

(1)首先，确定好需要解决分析的未知问题或者系统；

(2)建立系统分布的数学模型；

(3)利用蒙特卡罗(MC)方法模拟随机数；

(4)根据相关的参数条件来确定和选取相关的随机数；

(5)根据系统模型所服从的分布函数，求解出估计值。

图1 利用MC方法解决问题的逻辑思路图

3 主轴可靠性预测流程分析

主轴是数控机床主要的零部件，主轴的可靠性决定了机床的寿命，因为主轴是不容易破损，且寿命较长，主轴寿命检测实验会消耗较大的财力、人力、以及物力。因此采用蒙特卡罗方法仿真模主轴实验可以有效地解决此问题。图2表示利用蒙特卡罗方法分析机床主轴可靠性分析流程，集体步骤如下：

(1)实验过程中，通过采集机床主轴输出端电机的电压和电流信号，来计算出输出功率。为了保证实验数据的随机性，数据需要进行随机抽样。

(2)假设bi为某种概率分布的随机变量，利用数学模型找到主轴输出端电机功率信号与之对应的分布函数。

(3)结合表1并利用MATLAB[8]工具生成与变量bi对应的分布函数随机数，带入目标函数。

(4)设置主轴可靠性的判断条件，求解出满足条件的随机变量的个数，最后依据计算结果判断主轴的可靠性。

图2 基于蒙特卡罗主轴可靠性分析流程图

4 基于蒙特卡罗方法主轴可靠性预测实验分析

4.1 机床技术参数

本次实验分析过程中，采用学校机电实习中心的数控机床可靠性实验平台，由于主轴位置于机床内部部件，无法直接测量主轴振动信号，采用液压方式通过卡盘将加载棒与主轴相连接，利用电压和电路传感器采集输出端电机的电压和电流信号。机床详细技术参数如表2所示。

表2 高档数控机床可靠性能测试平台技术参数

4.2 主轴输出端电机功率时域图及置信区间

根据机床的工作原理，如果主轴的振动比较明显，则测得主轴输出端电机的运行功率(电压和电流乘积)的时域图振动变化比较明显。表2中有高档数控机床可靠性性能测试平台技术参数，实验过程中，假设采集的电流信号为I，电压信号为U，功率为P。根据试验台使用说明，次试验台的主轴输出端电机有A、B、C三相电压和三相电流，实验分析A相电压和电流，传感器测得的电压、电流信号需要经过公式(6)转换，从而得到电压、电流信号yv,yi：

(6)

(7)

(8)

(9)

在总体样本均值已知，方差未知情况下，结合统计学如公式(10)～公式(12)求出主轴输出端电机功率均值的区间估计。

(10)

(11)

(12)

其中，设置置信度α=0.05，可以求出数控机床主轴输出端电机功率的置信区间[9-10]为：

根据上述公式可以得到图3，试验台主轴输出端电机功率时域图中Y1和Y2求解公式如下：

已知采样频率为1024Hz，以1024的采样点为样本值，带入上述式子求得Y1=5814.90，Y2=6455.80。

图3 试验台主轴输出端电机功率时域图

图3中横坐标表示采样频率，纵坐标表示主轴电机输出功率值。

4.3 实验仿真及数据对比

利用蒙特卡罗方法生成符合抽样样本分布函数，设置不同的仿真次数，仿真结果如表3所示，待仿真实验结果趋于稳定时，根据实验结果判断株洲当前运行状态，从图4中可以看出试验台主轴输出端电机的功率最终稳定在6285.20，与图3电机功率时域图置信区间范围可以进行对比，试验台主轴正常工作下的功率处于置信区间内，得出当前状态下主轴的可靠性正常。

图4 基于蒙特卡罗仿真主轴输出端电机功率分布图

仿真次数n数控机床主轴输出端电机功率W60006286．370006283．580006285．290006285．2

如表3主轴输出端电机功率统计表所示，它是根据图4电机功率分布表统计出来的结果，从表3中可以看出，当我们实验仿真次数量增大时，主轴输出端电机功率不断减小，直到仿真次数大于8000之后，输出功率最终趋于6285.2W无偏估计值，此时我们还可以通过计算机强大计算能力，无限增大仿真次数，而结果依然处于稳定值不变。并且该无偏估计值正好处于样本均值的置信区间Y1、Y2之间，验证了计算机仿真结果数据的正确性。因此，本实验利用蒙特卡罗仿真产生符合抽样样本的总体分布的随机数，通过置信区间和功能函数计算出机床主轴可靠性的无偏估计，判断数控机床运行状态下主轴的可靠性。

5 结论

本文提出了一种基于蒙特卡罗模拟仿真的机床主轴可靠度求解方法。该方法借助于计算机技术计算出机床主轴输出端电机的功率并画出功率的时域图，通过置信区间和功能函数计算出机床主轴可靠性的无偏估计值，该无偏估计值正好处于样本均值的置信区间Y1、Y2之间(如图3所示)，验证了计算机仿真结果数据的正确性，实现了大样本数据简单、快速计算，节省了繁复的数学推理和演算过程，计算结果更加准确、可靠、通俗易懂，从而达到一种快速有效的检测效果。

[1] 何正嘉．机械设备运行可靠性评估的发展与思考[J]．机械工程学报，2016,50(2) :171-186.

[2] 刘瀚文．加工中心主轴可靠性试验研究[D]．长春: 吉林大学，2014.

[3] 訾佼佼,刘宏昭．基于退化量分布的电主轴可靠性评估[J]．中国机械工程，2016,25(6):807-812.

[4] 蒋喜,刘宏昭．基于伪寿命分布的电主轴极小子样可靠性研究[J]．振动与冲击，2015,32(19) :80-85.

[5] 肖刚,李天柁. 系统可靠性分析中的蒙特卡罗方法[M]. 北京:科学出版社,2013.

[6] 王梓坤.概率论基础及其应用[M].北京:北京师范大学出版社, 2006.

[7] 赵舜仁.关于随机变量和的特征函数[J].青岛建筑工程学院学报,2010 , 21(2):76-79.

[8] 陈桂明,戚红雨,潘伟.MATELAB数理统计[M].北京:科学出版社,2002.

[9] Stevens WL. Fiducial limits of the parameter of a discontinous distribution[J] . Biometrika,2016, 37(1-2):117-129.

[10] 籍永建，王红军. 基于EMD的主轴振动信号去噪方法研究[J]. 组合机床与自动化加工技术,2015(5):35-37.