基于性能退化预测的数控机床预防维修方法

代愽超 张英杰 李阳帆 陈 波

西安交通大学机械工程学院，西安，710049

0 引言

随着制造业整体实力的提高，通过寿命试验或加速寿命试验获取机床故障数据日趋困难且成本昂贵，为了保证数控机床加工质量的持续稳定性，采用合理的预防性维修方法已成为提高设备利用率和实现资产效率最大化的有效途径[1]。

生产控制方面，统计过程控制(statistical process control, SPC)作为一种集测定、记录、分析、控制等于一体的实用技术[2-4]，广泛应用于机床等设备的可靠性监测。其中普通控制图[2,5-6]、预控图[7]和选控图[8-9]直接关注评测目标的状态参数变化，而在状态参数检测设备不完善的生产线上，则可利用不合格品率控制图进行生产过程的控制[10-11]。

维护维修方面，大多数研究者将 “维修如新”作为假设，其原理简单、操作方便[12-14]。考虑到维修程度不同所带来的效果差异，“不完全维修”方面的研究也取得了一定的进展[15-16]。但“不完全维修”的维修程度难以掌控，因此一些学者认为将数控机床的维修视为“维修如旧”更为实用[17-18]。

目前在性能退化分析和预防维修决策方面的研究大多为单一方法研究，针对从状态预测到维修决策的全过程的方法依然较少，难以满足企业的实际需求。为此，本文提出了一种基于系统性能退化预测的数控机床预防维修方法。利用不合格品率控制图和Markov状态转移矩阵预测发生异常的时间、位置和原因，建立了Wiener退化模型来预测设备剩余寿命，计算并比较了维修和换新的后期收益效率，以实现对维修策略的指导。

1 基于SPC-Markov模型的机床状态分析

1.1 不合格品率控制图相关参数的计算

不合格品率控制图是建立在二项分布基础上的一种统计过程控制技术[10]。设样本的不合格品率为p，采用统计方法对其进行估计，表达式如下：

(1)

考虑到各样本中的产品数量可能不同，故对式(1)进行修正，得到

(2)

根据六西格玛(6σ)质量管理方法，结合二项分布的性质可计算得到中心线、上控制限和下控制限，其表达式分别如下：

(3)

1.2 基于Markov过程的异常概率预测

考虑到当前维修部位和维修水平会对下一次异常的发生造成影响，采用Markov状态转移矩阵进行各种异常模式的概率预测。

Markov状态转移矩阵如下：

(4)

其中，z为同一时刻最多可能的异常种类数；PIJ为t时刻处于I状态，到t+1时刻变为J状态的概率，且满足：

(5)

t时刻的状态空间向量

Et=(e1,e2,…,ez)

(6)

其中，e1,e2，…，ez的取值分别表示每一种异常是否发生，当取值为1时，表示该种异常发生；当取值为0时，该种异常不会发生。

则可求得t+1时刻各种异常情况的发生概率向量：

Gt+1=EtM

(7)

2 性能退化模型的建立及维修策略分析

2.1 基于Wiener过程的性能退化模型

Wiener退化过程模型[19-21]定义为

y(t)=a+μt+σBB(t)

(8)

式中，y(t)为t时刻的退化量测量值；a为初始退化量；μ为漂移系数；σB为扩散系数；B(t)为标准布朗运动值。

设i(i=1,2,…,N)为退化试验样本序号，N为退化试验样本总数；j(j=1,2,…,mi)为测量次数序号，mi为各样本的测量次数。进而可得到如下试验数据模型：

yij=ai+μitij+σBiBij

(9)

式中，tij、Bij、yij分别为第i个样本第j次的测量时间、布朗运动值和退化量测量值；ai、μi、σBi分别为第i个样本的初始退化量、漂移系数和扩散系数。

采用各样本分离的方法计算参数值，基本流程见图1。

图1 Wiener退化计算流程Fig.1 Evaluation process of Wiener degradation model

利用Wiener过程的高斯独立增量性可简化计算。令

(10)

j=2,3,…,mi

则

(11)

式中，Δtij、Δyij分别为样本i第j-1次和第j次测量之间的时间间隔和退化量增量。

建立似然函数：

(12)

先求得对数似然函数L(ai,μi,σBi)，再根据如下边界条件：

(13)

最终解得各样本参数的估计值分别为

(14)

(15)

(16)

将失效阈值记为Df，将y(t)首次达到Df的时间定义为机床寿命T，即

T=inf(t|y(t)=Df)t≥0

(17)

式中,inf(·)表示求集合中最大下界的操作。

退化初值a通常小于阈值Df，可靠度可描述为机床寿命T大于规定时间t的概率，即

(18)

实际求解时，选取参数a、μ、σB中变异系数较大的参数服从正态分布，以避免出现多重积分[20]。变异系数计算表达式如下：

(19)

式中,X为具体的变量(即本例中的a、μ和σB)；var(X)为X的方差；E(X)为X的期望。

2.2 基于剩余寿命和收益的维修策略分析

以“维修如旧”为前提，选取可靠度退化阈值Dr，认为系统平均寿命为退化曲线上的可靠度降低至Dr所需的总时间tA。则剩余寿命s(t)可表示为

s(t)=tA-t-t-

(20)

式中，t-为当前研究时段之前的总运行时间。

假设在机床正常工作期间，其日常运行的基础费用cm(单位时间)不变，则进行零部件维修和换新的后期收益效率分别为

(21)

(22)

式中，v为产品的出厂价格；w为单位时间产量；cr为每件产品所需原材料成本；Cf为维修费用；Cd为换新费用；tR为平均维修时间(MTTR)；tD为零部件更换平均用时(MTTD)。

通过比较式(21)、式(22)，将出现如下4种情况：

(1)ηf≥ηd且tR≤tD。ηf≥ηd表明在后续寿命内，相较于零部件换新，零部件维修产生的单位时间内的收益更高；又满足tR≤tD，表明维修造成的误工影响小。因此，宜进行关键零部件维修。

(2)ηf≥ηd且tR>tD。若当前生产任务紧迫则可先进行零部件更换，并在后续生产压力小的时候，对换下的零部件进行维修以备用。

(3)ηf<ηd且tR≤tD。若生产任务轻松，则直接更换零部件；若急需赶工且预期修一次能完成赶工任务，则建议维修，否则建议尽早更换零部件。

(4)ηf<ηd且tR>tD。此时更适合进行关键零部件的换新。

当tR和tD分别比s(t)和tA小很多时，利用式(21)、式(22)计算收益效率时可略去tR和tD以简化计算。此时若要满足ηf≥ηd，只需满足：

CftA≤Cds(t)

(23)

3 实例应用

3.1 生产线机床加工过程控制及状态预测

国内某活塞生产线的主要加工流程见图2。

图2 某活塞生产线加工流程Fig.2 Machining process of piston production line

本研究收集了该生产线30天(实际工作24天，第24天因检测出了大量返修品和报废品而被迫进行停机维修)的不合格品数据，选取前22天的数据，并计算得到不合格品率控制图参数，见表1，其中中心线值为0.049 6。

将各参数的计算结果绘制成控制图见图3，可以看出，由于每日产量的波动导致控制图的上下限不断变化，因此选择如下3个通用准则进行判异：①样本点距中心线大于3个标准差；②连续9点在中心线同一侧；③连续14点上下交错。

由图3可以看出，出现了4个异常点，其中点4、点21和点22均为数据超界点(情况①)，而点14则是因为连续超过9个点同侧(情况②)。由分析可知，点4为超下界点，而点6～14同样处于中心线的下侧，它们在统计规律上异常，但不合格品率偏低以及长期处于偏低态势都是对实际情况有利的现象，因此应查明引起这种情况的原因，并力图维持该状态以保证生产质量，但也要检测和分析是否是由检测偏差引起的。此外，点21为超上界点，表明系统可能出现了某种故障，点22的超界又加剧了这一可能性，因此预测系统将在第23个工作日发生异常。

表1 不合格品率控制图参数计算结果Tab.1 The calculation result of parameters for unqualified product rate control chart

图3 不合格品率控制图的预测结果Fig.3 Prediction results of unqualified product rate control chart

企业并未及时采取措施，且已将后续样本点添加到已有控制图中，结果见图4。

图4 不合格品率控制图的效果验证Fig.4 Effect validation of unqualified product rate control chart

图4是在图3的基础上绘制得到的，且计算中心线保持不变。由图4可以看出，点23回到了上下控制限范围内，但点22到点23是从超上界回落到中心线LC以下，跨度很大，生产过程极不稳定；点24则再次超出了上控制限。整个生产过程的不合格品率达到25%，报废量和返修量过大，企业无法接受。

利用Markov状态转移矩阵进行异常问题预测性诊断时，考虑到每台机床都有很多种潜在的异常类型，因数据量十分庞大，故无法详尽罗列完整的矩阵，本文仅分析下一次的异常最可能发生在哪一台机床。

收集得到该生产线的372条历史异常数据(实际状态转移371次)，见表2，其中，第Mi行Mj列的元素表示前一次异常发生在机床Mi，随后一次异常发生在机床Mj的情况出现的频数。

表2 各机床发生异常的转移关系Tab.2 Abnormity transmission relationship between machine tools

计算得到Markov状态转移矩阵：

M=

(24)

因收集不合格品数据前最近一次的异常发生在机床M7，故令该时间节点为t，则状态空间向量可表示为

Et=(0,0,0,0,0,0,1,0)

(25)

则t+1时刻各机床发生异常的概率向量可表示为

Gt+1=EtM=
(0.043,0.130,0.239,0.022,0.043,0.022,0.130,0.370)

(26)

依据式(26)，机床M8的异常概率最高，达到37%；机床M3的异常概率为23.9%以及机床M2、M7的异常概率均为13%，上述异常概率相对其他机床的异常概率也较高。由此可知，通过图3预测得到第23个工作日可能发生异常时，应优先对机床M8进行排查，若机床M8正常，则再检查机床M3、M2、M7，最后才检查其他机床。

在本生产实例中，正是M8机床的定位问题导致点22～24之间不合格品率的连续大幅度波动，这与模型预测的大概率事件相符合，表明该方法具有一定的准确性，可以引导异常排查，提高检修效率。

3.2 零部件剩余寿命预测及维修策略分析

根据预测第23个工作日发生异常，实际已工作22天(每天三班轮换，每班8 h)，即528 h。且预测最可能出现问题的是机床M8，若及时采取了相应检修措施，就会发现机床M8的定位元件出现了问题(但实际上是被迫停机后才检测出了定位问题)。

采用M8同类机床的定位精度数据进行分析。现有5台定位失效阈值Df=300 μm的样本机床，分别记作Si(i=1,2,3,4,5)，其定位精度检测结果见图5。

图5 定位精度的性能检测数据Fig.5 Performance testing data for positioning accuracy

利用式(14)～式(16)计算得到各样本初始退化量、漂移系数、扩散系数的估计值，见表3。表3中，样本5的初始退化量出现了负值，但考虑到有的机床需先使用一段时间后，定位精度才会提高到最佳状态，因此这里不做改动。

表3 各样本参数估计Tab.3 Parameter estimation of each sample

表4 参数估计值的变异性分析Tab.4 Variability analysis of the estimation values for parameters

(27)

(28)

式中，μa、σa分别为a的均值和标准差。

依据式(18)计算可靠度函数：

(29)

(30)

式中，φ(a)为a的概率密度函数。

取式(29)的积分区间为(-3σa,+3σa)，为补偿被忽略的区间，增加修正系数：

(31)

在MATLAB软件中绘制得到机床M8的可靠度退化曲线，见图6。

图6 机床M8的可靠度退化曲线Fig.6 Reliability degradation curve of M8

根据企业要求选取可靠度阈值Dr=0.6，由图6可知tA=1 300 h，又查得该定位元件在数据收集前实际已运行约21天，取t-=504 h，则有

s(t)=tA-t-t-=268 h

(32)

维修和换新所需时间均在3 h左右，可依据式(23)进行比较，得到优先维修的条件为

Cd/Cf≥4.85

(33)

对于引起本次故障的关键零部件，更换一次的平均费用约为2 360元，维修费用约为400元，满足式(33)。由2.2节中情况(1)可知：从经济性和时间效率两方面考虑，此时对关键定位元件进行维修比换新具有更好的综合效益。

求得tA后，根据式(20)，在当前成本条件下若要满足式(23)，只需满足：

t+t-≤1 079.66 h

(34)

这表明引起本次定位问题的关键零部件在总运行时间不超过1 079.66 h的情况下，优先选择维修是更加经济的选择。

4 结论

(1)提出了一种基于系统性能退化预测的数控机床预防维修方法，利用不合格品率控制图实现了生产线上机床设备的统计过程控制，降低了对专用检测工具的依赖性。

(2)基于Wiener退化模型，利用其增量的高斯独立性，采用各样本分离的方法求解模型参数，可较为准确地评估关键零部件的剩余寿命。

(3)基于“维修如旧”的假设，结合生产效益指导预防维修策略，且符合企业的实际利益需求。

(4)对国内某汽车活塞生产线数控机床的分析结果表明，所提方法的预测结果与实际情况基本一致，所构建的决策模型可为设备的预防性维修提供一定的参考。