代愽超 張英杰 李陽帆 陳 波

西安交通大學(xué)機(jī)械工程學(xué)院，西安，710049

0 引言

隨著制造業(yè)整體實(shí)力的提高，通過壽命試驗(yàn)或加速壽命試驗(yàn)獲取機(jī)床故障數(shù)據(jù)日趨困難且成本昂貴，為了保證數(shù)控機(jī)床加工質(zhì)量的持續(xù)穩(wěn)定性，采用合理的預(yù)防性維修方法已成為提高設(shè)備利用率和實(shí)現(xiàn)資產(chǎn)效率最大化的有效途徑[1]。

生產(chǎn)控制方面，統(tǒng)計(jì)過程控制(statistical process control, SPC)作為一種集測定、記錄、分析、控制等于一體的實(shí)用技術(shù)[2-4]，廣泛應(yīng)用于機(jī)床等設(shè)備的可靠性監(jiān)測。其中普通控制圖[2,5-6]、預(yù)控圖[7]和選控圖[8-9]直接關(guān)注評測目標(biāo)的狀態(tài)參數(shù)變化，而在狀態(tài)參數(shù)檢測設(shè)備不完善的生產(chǎn)線上，則可利用不合格品率控制圖進(jìn)行生產(chǎn)過程的控制[10-11]。

維護(hù)維修方面，大多數(shù)研究者將 “維修如新”作為假設(shè)，其原理簡單、操作方便[12-14]?？紤]到維修程度不同所帶來的效果差異，“不完全維修”方面的研究也取得了一定的進(jìn)展[15-16]。但“不完全維修”的維修程度難以掌控，因此一些學(xué)者認(rèn)為將數(shù)控機(jī)床的維修視為“維修如舊”更為實(shí)用[17-18]。

目前在性能退化分析和預(yù)防維修決策方面的研究大多為單一方法研究，針對從狀態(tài)預(yù)測到維修決策的全過程的方法依然較少，難以滿足企業(yè)的實(shí)際需求。為此，本文提出了一種基于系統(tǒng)性能退化預(yù)測的數(shù)控機(jī)床預(yù)防維修方法。利用不合格品率控制圖和Markov狀態(tài)轉(zhuǎn)移矩陣預(yù)測發(fā)生異常的時間、位置和原因，建立了Wiener退化模型來預(yù)測設(shè)備剩余壽命，計(jì)算并比較了維修和換新的后期收益效率，以實(shí)現(xiàn)對維修策略的指導(dǎo)。

1 基于SPC-Markov模型的機(jī)床狀態(tài)分析

1.1 不合格品率控制圖相關(guān)參數(shù)的計(jì)算

不合格品率控制圖是建立在二項(xiàng)分布基礎(chǔ)上的一種統(tǒng)計(jì)過程控制技術(shù)[10]。設(shè)樣本的不合格品率為p，采用統(tǒng)計(jì)方法對其進(jìn)行估計(jì)，表達(dá)式如下：

(1)

考慮到各樣本中的產(chǎn)品數(shù)量可能不同，故對式(1)進(jìn)行修正，得到

(2)

根據(jù)六西格瑪(6σ)質(zhì)量管理方法，結(jié)合二項(xiàng)分布的性質(zhì)可計(jì)算得到中心線、上控制限和下控制限，其表達(dá)式分別如下：

(3)

1.2 基于Markov過程的異常概率預(yù)測

考慮到當(dāng)前維修部位和維修水平會對下一次異常的發(fā)生造成影響，采用Markov狀態(tài)轉(zhuǎn)移矩陣進(jìn)行各種異常模式的概率預(yù)測。

Markov狀態(tài)轉(zhuǎn)移矩陣如下：

(4)

其中，z為同一時刻最多可能的異常種類數(shù)；PIJ為t時刻處于I狀態(tài)，到t+1時刻變?yōu)镴狀態(tài)的概率，且滿足：

(5)

t時刻的狀態(tài)空間向量

Et=(e1,e2,…,ez)

(6)

其中，e1,e2，…，ez的取值分別表示每一種異常是否發(fā)生，當(dāng)取值為1時，表示該種異常發(fā)生；當(dāng)取值為0時，該種異常不會發(fā)生。

則可求得t+1時刻各種異常情況的發(fā)生概率向量：

Gt+1=EtM

(7)

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

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

Wiener退化過程模型[19-21]定義為

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

(8)

式中，y(t)為t時刻的退化量測量值；a為初始退化量；μ為漂移系數(shù)；σB為擴(kuò)散系數(shù)；B(t)為標(biāo)準(zhǔn)布朗運(yùn)動值。

設(shè)i(i=1,2,…,N)為退化試驗(yàn)樣本序號，N為退化試驗(yàn)樣本總數(shù)；j(j=1,2,…,mi)為測量次數(shù)序號，mi為各樣本的測量次數(shù)。進(jìn)而可得到如下試驗(yàn)數(shù)據(jù)模型：

yij=ai+μitij+σBiBij

(9)

式中，tij、Bij、yij分別為第i個樣本第j次的測量時間、布朗運(yùn)動值和退化量測量值；ai、μi、σBi分別為第i個樣本的初始退化量、漂移系數(shù)和擴(kuò)散系數(shù)。

采用各樣本分離的方法計(jì)算參數(shù)值，基本流程見圖1。

圖1 Wiener退化計(jì)算流程Fig.1 Evaluation process of Wiener degradation model

利用Wiener過程的高斯獨(dú)立增量性可簡化計(jì)算。令

(10)

j=2,3,…,mi

則

(11)

式中，Δtij、Δyij分別為樣本i第j-1次和第j次測量之間的時間間隔和退化量增量。

建立似然函數(shù)：

(12)

先求得對數(shù)似然函數(shù)L(ai,μi,σBi)，再根據(jù)如下邊界條件：

(13)

最終解得各樣本參數(shù)的估計(jì)值分別為

(14)

(15)

(16)

將失效閾值記為Df，將y(t)首次達(dá)到Df的時間定義為機(jī)床壽命T，即

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

(17)

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

退化初值a通常小于閾值Df，可靠度可描述為機(jī)床壽命T大于規(guī)定時間t的概率，即

(18)

實(shí)際求解時，選取參數(shù)a、μ、σB中變異系數(shù)較大的參數(shù)服從正態(tài)分布，以避免出現(xiàn)多重積分[20]。變異系數(shù)計(jì)算表達(dá)式如下：

(19)

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

2.2 基于剩余壽命和收益的維修策略分析

以“維修如舊”為前提，選取可靠度退化閾值Dr，認(rèn)為系統(tǒng)平均壽命為退化曲線上的可靠度降低至Dr所需的總時間tA。則剩余壽命s(t)可表示為

s(t)=tA-t-t-

(20)

式中，t-為當(dāng)前研究時段之前的總運(yùn)行時間。

假設(shè)在機(jī)床正常工作期間，其日常運(yùn)行的基礎(chǔ)費(fèi)用cm(單位時間)不變，則進(jìn)行零部件維修和換新的后期收益效率分別為

(21)

(22)

式中，v為產(chǎn)品的出廠價(jià)格；w為單位時間產(chǎn)量；cr為每件產(chǎn)品所需原材料成本；Cf為維修費(fèi)用；Cd為換新費(fèi)用；tR為平均維修時間(MTTR)；tD為零部件更換平均用時(MTTD)。

通過比較式(21)、式(22)，將出現(xiàn)如下4種情況：

(1)ηf≥ηd且tR≤tD。ηf≥ηd表明在后續(xù)壽命內(nèi)，相較于零部件換新，零部件維修產(chǎn)生的單位時間內(nèi)的收益更高；又滿足tR≤tD，表明維修造成的誤工影響小。因此，宜進(jìn)行關(guān)鍵零部件維修。

(2)ηf≥ηd且tR>tD。若當(dāng)前生產(chǎn)任務(wù)緊迫則可先進(jìn)行零部件更換，并在后續(xù)生產(chǎn)壓力小的時候，對換下的零部件進(jìn)行維修以備用。

(3)ηf<ηd且tR≤tD。若生產(chǎn)任務(wù)輕松，則直接更換零部件；若急需趕工且預(yù)期修一次能完成趕工任務(wù)，則建議維修，否則建議盡早更換零部件。

(4)ηf<ηd且tR>tD。此時更適合進(jìn)行關(guān)鍵零部件的換新。

當(dāng)tR和tD分別比s(t)和tA小很多時，利用式(21)、式(22)計(jì)算收益效率時可略去tR和tD以簡化計(jì)算。此時若要滿足ηf≥ηd，只需滿足：

CftA≤Cds(t)

(23)

3 實(shí)例應(yīng)用

3.1 生產(chǎn)線機(jī)床加工過程控制及狀態(tài)預(yù)測

國內(nèi)某活塞生產(chǎn)線的主要加工流程見圖2。

圖2 某活塞生產(chǎn)線加工流程Fig.2 Machining process of piston production line

本研究收集了該生產(chǎn)線30天(實(shí)際工作24天，第24天因檢測出了大量返修品和報(bào)廢品而被迫進(jìn)行停機(jī)維修)的不合格品數(shù)據(jù)，選取前22天的數(shù)據(jù)，并計(jì)算得到不合格品率控制圖參數(shù)，見表1，其中中心線值為0.049 6。

將各參數(shù)的計(jì)算結(jié)果繪制成控制圖見圖3，可以看出，由于每日產(chǎn)量的波動導(dǎo)致控制圖的上下限不斷變化，因此選擇如下3個通用準(zhǔn)則進(jìn)行判異：①樣本點(diǎn)距中心線大于3個標(biāo)準(zhǔn)差；②連續(xù)9點(diǎn)在中心線同一側(cè)；③連續(xù)14點(diǎn)上下交錯。

由圖3可以看出，出現(xiàn)了4個異常點(diǎn)，其中點(diǎn)4、點(diǎn)21和點(diǎn)22均為數(shù)據(jù)超界點(diǎn)(情況①)，而點(diǎn)14則是因?yàn)檫B續(xù)超過9個點(diǎn)同側(cè)(情況②)。由分析可知，點(diǎn)4為超下界點(diǎn)，而點(diǎn)6～14同樣處于中心線的下側(cè)，它們在統(tǒng)計(jì)規(guī)律上異常，但不合格品率偏低以及長期處于偏低態(tài)勢都是對實(shí)際情況有利的現(xiàn)象，因此應(yīng)查明引起這種情況的原因，并力圖維持該狀態(tài)以保證生產(chǎn)質(zhì)量，但也要檢測和分析是否是由檢測偏差引起的。此外，點(diǎn)21為超上界點(diǎn)，表明系統(tǒng)可能出現(xiàn)了某種故障，點(diǎn)22的超界又加劇了這一可能性，因此預(yù)測系統(tǒng)將在第23個工作日發(fā)生異常。

表1 不合格品率控制圖參數(shù)計(jì)算結(jié)果Tab.1 The calculation result of parameters for unqualified product rate control chart

圖3 不合格品率控制圖的預(yù)測結(jié)果Fig.3 Prediction results of unqualified product rate control chart

企業(yè)并未及時采取措施，且已將后續(xù)樣本點(diǎn)添加到已有控制圖中，結(jié)果見圖4。

圖4 不合格品率控制圖的效果驗(yàn)證Fig.4 Effect validation of unqualified product rate control chart

圖4是在圖3的基礎(chǔ)上繪制得到的，且計(jì)算中心線保持不變。由圖4可以看出，點(diǎn)23回到了上下控制限范圍內(nèi)，但點(diǎn)22到點(diǎn)23是從超上界回落到中心線LC以下，跨度很大，生產(chǎn)過程極不穩(wěn)定；點(diǎn)24則再次超出了上控制限。整個生產(chǎn)過程的不合格品率達(dá)到25%，報(bào)廢量和返修量過大，企業(yè)無法接受。

利用Markov狀態(tài)轉(zhuǎn)移矩陣進(jìn)行異常問題預(yù)測性診斷時，考慮到每臺機(jī)床都有很多種潛在的異常類型，因數(shù)據(jù)量十分龐大，故無法詳盡羅列完整的矩陣，本文僅分析下一次的異常最可能發(fā)生在哪一臺機(jī)床。

收集得到該生產(chǎn)線的372條歷史異常數(shù)據(jù)(實(shí)際狀態(tài)轉(zhuǎn)移371次)，見表2，其中，第Mi行Mj列的元素表示前一次異常發(fā)生在機(jī)床Mi，隨后一次異常發(fā)生在機(jī)床Mj的情況出現(xiàn)的頻數(shù)。

表2 各機(jī)床發(fā)生異常的轉(zhuǎn)移關(guān)系Tab.2 Abnormity transmission relationship between machine tools

計(jì)算得到Markov狀態(tài)轉(zhuǎn)移矩陣：

M=

(24)

因收集不合格品數(shù)據(jù)前最近一次的異常發(fā)生在機(jī)床M7，故令該時間節(jié)點(diǎn)為t，則狀態(tài)空間向量可表示為

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

(25)

則t+1時刻各機(jī)床發(fā)生異常的概率向量可表示為

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

(26)

依據(jù)式(26)，機(jī)床M8的異常概率最高，達(dá)到37%；機(jī)床M3的異常概率為23.9%以及機(jī)床M2、M7的異常概率均為13%，上述異常概率相對其他機(jī)床的異常概率也較高。由此可知，通過圖3預(yù)測得到第23個工作日可能發(fā)生異常時，應(yīng)優(yōu)先對機(jī)床M8進(jìn)行排查，若機(jī)床M8正常，則再檢查機(jī)床M3、M2、M7，最后才檢查其他機(jī)床。

在本生產(chǎn)實(shí)例中，正是M8機(jī)床的定位問題導(dǎo)致點(diǎn)22～24之間不合格品率的連續(xù)大幅度波動，這與模型預(yù)測的大概率事件相符合，表明該方法具有一定的準(zhǔn)確性，可以引導(dǎo)異常排查，提高檢修效率。

3.2 零部件剩余壽命預(yù)測及維修策略分析

根據(jù)預(yù)測第23個工作日發(fā)生異常，實(shí)際已工作22天(每天三班輪換，每班8 h)，即528 h。且預(yù)測最可能出現(xiàn)問題的是機(jī)床M8，若及時采取了相應(yīng)檢修措施，就會發(fā)現(xiàn)機(jī)床M8的定位元件出現(xiàn)了問題(但實(shí)際上是被迫停機(jī)后才檢測出了定位問題)。

采用M8同類機(jī)床的定位精度數(shù)據(jù)進(jìn)行分析?，F(xiàn)有5臺定位失效閾值Df=300 μm的樣本機(jī)床，分別記作Si(i=1,2,3,4,5)，其定位精度檢測結(jié)果見圖5。

圖5 定位精度的性能檢測數(shù)據(jù)Fig.5 Performance testing data for positioning accuracy

利用式(14)～式(16)計(jì)算得到各樣本初始退化量、漂移系數(shù)、擴(kuò)散系數(shù)的估計(jì)值，見表3。表3中，樣本5的初始退化量出現(xiàn)了負(fù)值，但考慮到有的機(jī)床需先使用一段時間后，定位精度才會提高到最佳狀態(tài)，因此這里不做改動。

表3 各樣本參數(shù)估計(jì)Tab.3 Parameter estimation of each sample

表4 參數(shù)估計(jì)值的變異性分析Tab.4 Variability analysis of the estimation values for parameters

(27)

(28)

式中，μa、σa分別為a的均值和標(biāo)準(zhǔn)差。

依據(jù)式(18)計(jì)算可靠度函數(shù)：

(29)

(30)

式中，φ(a)為a的概率密度函數(shù)。

取式(29)的積分區(qū)間為(-3σa,+3σa)，為補(bǔ)償被忽略的區(qū)間，增加修正系數(shù)：

(31)

在MATLAB軟件中繪制得到機(jī)床M8的可靠度退化曲線，見圖6。

圖6 機(jī)床M8的可靠度退化曲線Fig.6 Reliability degradation curve of M8

根據(jù)企業(yè)要求選取可靠度閾值Dr=0.6，由圖6可知tA=1 300 h，又查得該定位元件在數(shù)據(jù)收集前實(shí)際已運(yùn)行約21天，取t-=504 h，則有

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

(32)

維修和換新所需時間均在3 h左右，可依據(jù)式(23)進(jìn)行比較，得到優(yōu)先維修的條件為

Cd/Cf≥4.85

(33)

對于引起本次故障的關(guān)鍵零部件，更換一次的平均費(fèi)用約為2 360元，維修費(fèi)用約為400元，滿足式(33)。由2.2節(jié)中情況(1)可知：從經(jīng)濟(jì)性和時間效率兩方面考慮，此時對關(guān)鍵定位元件進(jìn)行維修比換新具有更好的綜合效益。

求得tA后，根據(jù)式(20)，在當(dāng)前成本條件下若要滿足式(23)，只需滿足：

t+t-≤1 079.66 h

(34)

這表明引起本次定位問題的關(guān)鍵零部件在總運(yùn)行時間不超過1 079.66 h的情況下，優(yōu)先選擇維修是更加經(jīng)濟(jì)的選擇。

4 結(jié)論

(1)提出了一種基于系統(tǒng)性能退化預(yù)測的數(shù)控機(jī)床預(yù)防維修方法，利用不合格品率控制圖實(shí)現(xiàn)了生產(chǎn)線上機(jī)床設(shè)備的統(tǒng)計(jì)過程控制，降低了對專用檢測工具的依賴性。

(2)基于Wiener退化模型，利用其增量的高斯獨(dú)立性，采用各樣本分離的方法求解模型參數(shù)，可較為準(zhǔn)確地評估關(guān)鍵零部件的剩余壽命。

(3)基于“維修如舊”的假設(shè)，結(jié)合生產(chǎn)效益指導(dǎo)預(yù)防維修策略，且符合企業(yè)的實(shí)際利益需求。

(4)對國內(nèi)某汽車活塞生產(chǎn)線數(shù)控機(jī)床的分析結(jié)果表明，所提方法的預(yù)測結(jié)果與實(shí)際情況基本一致，所構(gòu)建的決策模型可為設(shè)備的預(yù)防性維修提供一定的參考。