# Condition-Based Predictive Order Model for a Mechanical Component following Inverse Gaussian Degradation Process

Authors: Cheng Wang; Jianxin Xu; Hongjun Wang; Zhenming Zhang
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

An efficient condition-based predictive spare ordering approach is the key to guarantee safe operation, improve service quality, and reduce maintenance costs under a predefined lower availability threshold. In this paper, we propose a condition-based predictive order model (CBPO) for a mechanical component, whose degradation path is modeled as inverse Gaussian (IG) process with covariate effect. The CBPO is dependent on the remaining useful life (RUL), random lead-time, speed-up lead-time degree, and availability threshold. RUL estimation is obtained through the IG degradation process at each inspection time. Both regular lead-time and expedited lead-time considered in RUL-based spare ordering policy can be cost-effective and reduce losses caused by unexpected failure. Speed-up lead-time can meet the urgent needs for spare parts on site. The decision variable of CBPO is the spare ordering time. Based on the CBPO under the lower availability threshold constraint, the objective of this study is to determine the optimal spare ordering time such that the expected cost rate is minimized. Finally, a case study of the mechanical spindle is presented to illustrate the proposed model and sensitivity analysis on critical parameters is performed.


---

## Body

## 1. Introduction

With the development of science and technology, the safety and effective perception of some products, especially for some high-end products, become more and more important. Product degradation condition prediction and efficient spare ordering are the key to guarantee safe operation, improve service quality, and reduce maintenance costs. It is of great importance to avoid the failure of a highly reliable product, especially for aircraft engine, high-end machine tools, high-speed train, wind turbines, and so on. In order to prevent or reduce product unexpected failure, various preventive maintenance (PM) policies have been proposed over the past several decades [1]. In PM policies, there are mainly two categories: time-based maintenance (TBM) policies and condition-based maintenance (CBM) policies. For TBM, much attention has been paid to lifetime distribution models. For CBM, decision is made according to degradation level or health condition. More recently, rapid development of computer-based condition monitoring (CBCM) technologies (e.g., advanced sensors) and Internet of Things (IoT) make it possible to track the product health in a real-time manner and have further facilitated efficient potential failure detection and CBM practices. Comparing with TBM, CBM is a more promising maintenance policy since it emphasizes combining data-driven reliability models with condition monitoring data. Therefore, CBM has received considerable attention in both academia and industry [2–9].Most CBM policies studied in the literatures assume that at any time there is an unlimited supply of available spares for replacement. However, this assumption is generally unrealistic and unpractical when available spares are limited and/or delivery lead times are much longer. When spares are expensive, scarce, and with higher and random lead times, it is important to consider shortage cost and holding cost. Therefore, when to order spare is a key factor to make tradeoff between shortage cost and holding cost. In other words, reasonable spare ordering time can minimize the entire maintenance cost in engineering practice. Motivated by addressing this type of problem, some order policies have been extensively researched. Wanget al. [10] proposed a joint optimization of CBM and spare ordering management for a single-unit system. Chien and Chen [11] presented an age-based spare ordering policy with random lead-time. Chien [12] presented an optimal spares ordering policy based on the optimal number of minimal repairs with random lead-time. Louitet al. [13] presented an order policy based on remaining useful life (RUL) estimation through the assessment of the component age and condition indicators. Godoyet al. [14] presented a graphic technique that considered a rule for decision based on both condition-based reliability function and a random/fixed lead-time. Panagiotidou [15] proposed a joint optimization of spares ordering and maintenance policies for multiple identical times. Wanget al. [16] proposed a prognostics-based spare ordering and system replacement policy with random lead-time for a deteriorating system. Chenet al. [17] proposed a joint optimization of replacement and spare ordering for critical rotary component based on collected condition monitoring signals. Caiet al. [18] proposed an appointment policy of spares based on (s,S) policy. Among them, prognostics-based order policy is few.In reality, prognostic largely focuses on estimating the RUL of the component by using life distribution or the available degradation data. Therefore, RUL estimation [19–26] can offer adequate lead-time for the maintainer to implement the essential maintenance actions ahead of failure. The impact of the dynamic environment (e.g., temperature, stress, humidity, or dust) on the RUL is indicated by covariates [3, 13, 27, 28]. In CBM policies, prognostic is based on product’s degradation process through available degradation data. Stochastic models with continuous-state are widely used to characterize the degradation process. Notable among them, the Wiener process [26, 29], the Gamma process [30], and the inverse Gaussian (IG) process [28, 31–33], including their variants, attract significant attention because of their nice mathematical properties and clear physical interpretations. CBM policies based on the Wiener process [26, 34, 35], the Gamma process [36, 37], and the IG process [38] have been extensively investigated. The IG process is a very effective method to characterize the random effects and covariates and has been frequently used to model the monotone degradation process. For most mechanical components, the path of degradation increases monotonically. Based on the merit of IG process, the mechanical component’s degradation path is modeled as IG process with covariates.This paper aims to develop a condition-based predictive order model (CBPO) for a mechanical component subject to IG degradation process with covariate effect. Covariates provide an indication of the external state and affect the degradation process of the mechanical component. When the degradation level of the component attains a predetermined threshold, the component fails. The proposed model depends on the RUL, random lead-time, speed-up lead-time degree, and availability threshold. The decision variable of the model is the spare ordering time. RUL estimation is obtained through the IG degradation process at each inspection time. Random lead-time includes both regular lead-time and expedited lead-time to be cost-effective and reduce losses caused by unexpected failure. Speed-up lead-time is triggered when the component failure occurs before the ordered spare is delivered. Based on the proposed model under the lower availability threshold constraint, what we aim is to find the optimal spare ordering time minimizing the expected cost rate. This paper differs from the existing works in following aspects. (a) Expedited order is considered in RUL-based spare ordering policy for the mechanical component subject to IG degradation process with covariate effect. (b) Speed-up lead-time is triggered when the component failure occurs before the ordered spare is delivered. (c) On the basis of cost rate objective function, we add the availability threshold constraint.The remainder of this paper is organized as follows. In Section2, we describe some assumptions. Section 3 introduces IG process and presents prediction method of RUL based on IG process. Based on RUL, a novel spare ordering model for the mechanical component is developed in Section 4. In Section 5, a case study of the mechanical spindle is given to illustrate the proposed order model and sensitivity analysis on critical parameters is performed. Section 6 concludes this paper and offers possible research in future.


## 2. Assumptions

The spare ordering problem is considered under the following assumptions:(i)
When degradation level of the component exceeds the predetermined failure thresholdl, the component fails.(ii)
The mechanical component's health condition can be detected perfectly by CBCM.(iii)
The mechanical component begins operating at time 0.(iv)
When the mechanical component is failed, the expedited order is triggered immediately. That is to say, it is not to consider the downtime and the corresponding shortage cost between the failed time and spare ordering time.(v)
Shortage cost per unit time is bigger than holding cost per unit time due to the component shutdown affecting the production progress and custom service negatively, i.e.,ρs>ρh.(vi)
The lead-time is a generally distributed random variable. There are two types of lead-time, i.e., regular lead-time and expedited lead-time. In regular lead-time distribution, letwt be the probability density function (PDF) of regular lead-time, Wt be the cumulative distribution function (CDF) of regular lead-time, W¯t be the survival function, μr be the mean of regular lead-time, and σr be the standard deviation of regular lead-time. In expedited lead-time distribution, let ht be the PDF of expedited lead-time, Ht be the CDF of expedited lead-time, H¯t be the survival function, μe be the mean of expedited lead-time, and σe be the standard deviation of expedited lead-time. And in general, we assume that μr is bigger than μe, i.e., μr>μe.(vii)
Replacements are made perfectly and do not affect the component's characteristics.(viii)
The spare in stock does not degenerate or fail.


## 3. Model Statements

### 3.1. IG Process

In practice, the performance of numerous mechanical components degrades over time, which can be modeled by a stochastic process. In our research, we will employ the IG process to represent monotonic degradation process. LetX(t) be degradation of the mechanical component at time t, and then degradation process {X(t),t>0} is called an IG process if it satisfies the following properties:(a)X(0)=0 with probability one;(b)X(t) has independent increments on nonoverlapping intervals, i.e., X(t)-X(s) and X(u)-X(v) are independent for t>s>u>v≥0;(c)X(t)-X(s) is subject to the IG distribution IG(μΛ(t)-Λ(s), ηΛ(t)-Λ(s)2) for t>s≥0, where Λ(t) is a monotone increasing function with Λ(0)=0, μ>0 and η>0 are constants, and IGa,b, a>0, b>0 denotes the IG distribution with PDF(1)fy;a,b=b2πy31/2exp⁡-by-a22a2y,where y>0 and CDF(2)Fy;a,b=Φbyya-1+exp⁡2baΦbyya+1,where y>0.To characterize the influence of the covariates on the degradation process, we use a model similar to a proportional hazards model (PHM) [39, 40]. The covariate modifies the baseline rates Λt as follows:(3)Λt,Zt=ψZtΛthZt,where ψZt is nonnegative function such that ψZ0=hZ0=1. For the sake of simplification, we assume hZt=1 and ψ· to be an exponential function which regulates the impacts of the covariate on the degradation.More precisely, for a given value of the covariate at timet, Zt=zt, the shape function of the IG process is defined by(4)Λt,zt=Λtexp⁡∑iγiZit,where γi is regression parameters defining the effects of the covariate of the IG process and Zt=Z1t,Z2t,... is a vector of time-dependent covariates.Let the constantll>0 be the failure threshold of the process X(t),t>0. We assume that the degradation follows an IG process model with function Λ(t) and parameters μ and η. Since X(0)=0, X(t) yields IGμΛ(t),ηΛ(t)2 with mean μΛ(t) and variance μ3Λ(t)/η.We define the first passage time to the failure thresholdl as T=inf⁡t≥0∣Xt≥l. Based on (2), then T follows the IG distribution with CDF given by(5)Flt=Pr⁡T<t=Pr⁡Xt≥l=1-Fl;μΛt,ηΛt2=ΦηlΛt,zt-lμ-exp⁡2ηΛt,ztμΦ-ηlΛt,zt+lμ,where parameters η and μ are positive and Φ(·) is the standard normal distribution function.To describe the IG process clearly, some assumptions are made as follows.(a) We assume that the shape functionΛ(t) is a linear function, i.e., Λ(t)=t.(b) The external force can affect the degradation process of the mechanical component. Based on (4), Λ(t) under the external force can be expressed by(6)Λt,zt=Λtexp⁡γZt.Let Zt denote the state of operating environment at time t, and Zt,t≥0 take values in E=0,1,2,...,r.(c) Under the external force, the degradation pathXt can be determined by the current stress and the accumulated degradation. The external force not only has an effect on the degradation speed but also the degradation volatility.According to (5), the mean of the lifetime, i.e., mean-time-to-failure (MTTF), can be easily formulated by(7)Tl=ET=∫0∞1-Fltdt=∫0∞Fl¯tdt.In addition, according to the nature of IG process, we denoteθ=μ,η,γ for the IG model parameters, and we obtain parameters estimation θ^=μ^,η^,γ^ by maximum likelihood estimation (MLE) at each time tk based on the CBCM data X0:k.


### 3.2. RUL Estimation

LetX0:k=x0,x1,…,xk denote the degradation observation at 0=t0<t1<⋯<tk, which can be irregularly spaced, and xi=Xti represents the degradation measurement of the mechanical component at time ti. Therefore, using the first passage time of the degradation process, we define the RUL, abbreviated as Rk that corresponds to Rk=inf⁡rk>0:Xtk+rk≥l∣X0:k.Based on the definition of the RUL and condition on the degradation dataX0:k, the distribution function of the RUL can be obtained by the following:(8)FRkrk;X0:k=Pr⁡Xtk+rk≥l∣X0:k=Pr⁡Xtk+rk-xk≥l-xk=Φηl-xkΛrk,zrk-l-xkμ-exp⁡2ηΛrk,zrkμ·Φ-ηl-xkΛrk,zrk+l-xkμ.For (8), the mean of Rk can be easily obtained as follows:(9)ERk=∫0∞1-FRkrk;X0:kdt=∫0∞F-Rkrk;X0:kdt,where F-Rk(rk;X0:k) denotes reliability function of the RUL.


## 3.1. IG Process

In practice, the performance of numerous mechanical components degrades over time, which can be modeled by a stochastic process. In our research, we will employ the IG process to represent monotonic degradation process. LetX(t) be degradation of the mechanical component at time t, and then degradation process {X(t),t>0} is called an IG process if it satisfies the following properties:(a)X(0)=0 with probability one;(b)X(t) has independent increments on nonoverlapping intervals, i.e., X(t)-X(s) and X(u)-X(v) are independent for t>s>u>v≥0;(c)X(t)-X(s) is subject to the IG distribution IG(μΛ(t)-Λ(s), ηΛ(t)-Λ(s)2) for t>s≥0, where Λ(t) is a monotone increasing function with Λ(0)=0, μ>0 and η>0 are constants, and IGa,b, a>0, b>0 denotes the IG distribution with PDF(1)fy;a,b=b2πy31/2exp⁡-by-a22a2y,where y>0 and CDF(2)Fy;a,b=Φbyya-1+exp⁡2baΦbyya+1,where y>0.To characterize the influence of the covariates on the degradation process, we use a model similar to a proportional hazards model (PHM) [39, 40]. The covariate modifies the baseline rates Λt as follows:(3)Λt,Zt=ψZtΛthZt,where ψZt is nonnegative function such that ψZ0=hZ0=1. For the sake of simplification, we assume hZt=1 and ψ· to be an exponential function which regulates the impacts of the covariate on the degradation.More precisely, for a given value of the covariate at timet, Zt=zt, the shape function of the IG process is defined by(4)Λt,zt=Λtexp⁡∑iγiZit,where γi is regression parameters defining the effects of the covariate of the IG process and Zt=Z1t,Z2t,... is a vector of time-dependent covariates.Let the constantll>0 be the failure threshold of the process X(t),t>0. We assume that the degradation follows an IG process model with function Λ(t) and parameters μ and η. Since X(0)=0, X(t) yields IGμΛ(t),ηΛ(t)2 with mean μΛ(t) and variance μ3Λ(t)/η.We define the first passage time to the failure thresholdl as T=inf⁡t≥0∣Xt≥l. Based on (2), then T follows the IG distribution with CDF given by(5)Flt=Pr⁡T<t=Pr⁡Xt≥l=1-Fl;μΛt,ηΛt2=ΦηlΛt,zt-lμ-exp⁡2ηΛt,ztμΦ-ηlΛt,zt+lμ,where parameters η and μ are positive and Φ(·) is the standard normal distribution function.To describe the IG process clearly, some assumptions are made as follows.(a) We assume that the shape functionΛ(t) is a linear function, i.e., Λ(t)=t.(b) The external force can affect the degradation process of the mechanical component. Based on (4), Λ(t) under the external force can be expressed by(6)Λt,zt=Λtexp⁡γZt.Let Zt denote the state of operating environment at time t, and Zt,t≥0 take values in E=0,1,2,...,r.(c) Under the external force, the degradation pathXt can be determined by the current stress and the accumulated degradation. The external force not only has an effect on the degradation speed but also the degradation volatility.According to (5), the mean of the lifetime, i.e., mean-time-to-failure (MTTF), can be easily formulated by(7)Tl=ET=∫0∞1-Fltdt=∫0∞Fl¯tdt.In addition, according to the nature of IG process, we denoteθ=μ,η,γ for the IG model parameters, and we obtain parameters estimation θ^=μ^,η^,γ^ by maximum likelihood estimation (MLE) at each time tk based on the CBCM data X0:k.


## 3.2. RUL Estimation

LetX0:k=x0,x1,…,xk denote the degradation observation at 0=t0<t1<⋯<tk, which can be irregularly spaced, and xi=Xti represents the degradation measurement of the mechanical component at time ti. Therefore, using the first passage time of the degradation process, we define the RUL, abbreviated as Rk that corresponds to Rk=inf⁡rk>0:Xtk+rk≥l∣X0:k.Based on the definition of the RUL and condition on the degradation dataX0:k, the distribution function of the RUL can be obtained by the following:(8)FRkrk;X0:k=Pr⁡Xtk+rk≥l∣X0:k=Pr⁡Xtk+rk-xk≥l-xk=Φηl-xkΛrk,zrk-l-xkμ-exp⁡2ηΛrk,zrkμ·Φ-ηl-xkΛrk,zrk+l-xkμ.For (8), the mean of Rk can be easily obtained as follows:(9)ERk=∫0∞1-FRkrk;X0:kdt=∫0∞F-Rkrk;X0:kdt,where F-Rk(rk;X0:k) denotes reliability function of the RUL.


## 4. Order Model

### 4.1. Order Process

To facilitate the implementation process of order model, the main procedures of the model at each time are summarized as follows.Step 1.
Parameters initialization: set the initial value of inspection timesk by choosing randomly or by the prior knowledge.Step 2.
Parameters estimation: estimate the degradation model parametersθ^=μ^,η^,γ^ by MLE according to the nature of IG process. The available estimation samples are obtained from the observed degradation data of the mechanical component till time tk, i.e., X0:k.Step 3.
RUL estimation: based on the real-time estimated parametersθ^=μ^,η^,γ^ in Step 2, we can derive the RUL distribution of the mechanical component at time tk immediately.Step 4.
Applying order model (i.e., CBPO): based on Step3, random lead-time distribution, and speed-up lead-time degree, we can seek the optimal spare ordering time tk0∗ that minimizes the expected cost rate under lower availability threshold by (20).Step 5.
Decision-making: according to the result of Step4, if the optimal spare ordering time tk0∗ is not equal to the next inspection time tk+1, then return to Step 1 and wait for the next inspection time for new decision. Otherwise, stop the procedure and place an order simultaneously. Figure 1 shows the flow chart of the implementation process of order model at each time tk.Figure 1
Flow chart of the implementation process of order model at each timetk.


### 4.2. Order Policy

The order policy is made up of regular order policy and expedited order policy. When the scheduled spare ordering time occurs before the component failure time, regular order policy is triggered. When the scheduled spare ordering time occurs after the component failure time, expedited order policy is triggered. The detailed order policy for the component is depicted in Figure2, where rk is a remaining time, T is the component failure time, Lr is regular lead-time, Le is expedited lead-time, tk is the current inspection time, tko is the scheduled spare ordering time, H is the holding time, and S is the shortage time.Figure 2
Possible order-replacement states of one cycle.State 1.
If the failed component occurs aftertko and the spare is delivered before the component failure, the failed component will not be replaced until the residual lifetime is equal to zero. Let Ps1 denote the probability of State 1 and corresponding CDF can be shown as(10)Ps1=Prtko+Lr<tk+rkand0<tko-tk<rk=∫tko-tk∞∫0tk+u-tkodWtdFRku,where t,u denote integral variable, respectively.State 2.
The failed component occurs aftertko and the spare is not delivered when the component failure occurs. To reduce the wait time after the component failure, we speed up the lead-time. The failed component is replaced by the spare as soon as the spare is delivered. Let Ps2 denote the probability of State 2 and corresponding CDF can be shown as(11)Ps2=Pr⁡tko+Lr>tk+rkand0<tko-tk<rk=∫tko-tk∞∫tk+u-tko∞dWtdFRku,where t,u denote integral variable, respectively.State 3.
The mechanical component failure occurs beforetko; an expedited order is made and the failed component is replaced by the spare as soon as the spare is delivered. Let Ps3 denote the probability of State 3 and corresponding CDF can be shown as(12)Ps3=Pr⁡tko>tk+rk=∫0∞∫0tko-tkdFRktdHu,where t, u denote integral variable, respectively.


### 4.3. Model of Objective Function

The main goal of this section is to build the expected cost rate model under lower availability threshold constraint.


#### 4.3.1. Cost Rate Model

(  1) The Expected Replacement Cycle Length. As shown in Figure 2 and formulation of the order-replacement states, a replacement cycle can be formulated as(13)U=tk+rk,iftko+Lr<tk+rkand0<tko-tk<rk,tko+Lr,iftko+Lr>tk+rkand0<tko-tk<rk,tko+Le,iftko>tk+rk,where U is a replacement cycle length.Since holding time only occurs in State 1, the expected holding timeEH during the replacement cycle is(14)EH=Etk+rk-tko-Lr=∫tko-tk∞∫0tk+u-tkotk+u-tko-tdWtdFRku,where t,u denote integral variable, respectively.Furthermore, since the shortage time may occur in State 2 and State 3, the expected shortage timeES during the replacement cycle is(15)ES=Etko+Lr-tk-rk1+κ+ELe=∫tko-tk∞∫tk+u-tko∞tko+t-tk-udWtdFRku1+κ+∫0∞uFRktko-tkhudu,0≤κ<∞where t,u denote integral variable, respectively, and κ denotes speed-up lead-time degree.Based on the above analyses and using expression (13), the expected replacement cycle length can be formulated as(16)EUtko=tk+∫0∞F-Rkrk;X0:kdt+ES=Tf+ES.In other words, the expected replacement cycle length can be denoted by the expected failure lifetime and the expected shortage time.(  2) The Expected Replacement Cycle Cost. As shown in Figure 2, the expected cost per replacement cycle can be derived, which includes inspection cost kρi, replacement cost C, regular ordering cost Cr, regular delivery cost Crd, expedited ordering cost Ce, expedited delivery cost Ced, speed-up lead-time cost κCed, holding cost ρhEH, shortage cost ρsES, and salvage value Csv, where Ce>Cr,Ced>Crd. Therefore, by (10), (11), (12), (14), and (15), the expected cost per replacement cycle EVtko can be expressed as(17)EVtko=k·ρi+C+Cr·Ps1+Ps2+Ce·Ps3+Crd·Ps1+Crd+κ·Ced·Ps2+Ced·Ps3+ρh·EH+ρs·ES-Csv.(  3) Cost Rate Model. Therefore, based on the renewal-reward theory in Ross [41], the expected cost rate is obtained as(18)CRtko=EVtkoEUtko,where tk<tko, EUtko and EVtko are, respectively, given by (16) and (17).


#### 4.3.2. Availability

Availability is an important performance index of repairable components and hereafter denoted byA. Based on key renewal theorem in Ross [41], A can be defined to be such that(19)Atko=Erun  timeErun  time+Edown  time=TlTl+ES=11+ES/Tlwhere Erun  time and Edown  time correspond to the expected run time and down time. Erun  time is equal to Tl and Edown  time is equal to ES. Based on (19), we can find that A decreases as ES increases.


#### 4.3.3. Cost Rate Model Subject to Availability

Consider the availability as a constraint attached to the cost rate model. Modifying (18), the cost rate model subject to availability is as follows:(20)CRtko=EVtkoEUtko,s.t.A≥A∗where A∗ is a predefined lower availability threshold. Based on (20), the optimal spare ordering time tko∗ can be obtained at each inspection time tk simultaneously. What we aim is to seek tko∗ by minimizing the expected cost rate under lower availability threshold.


## 4.1. Order Process

To facilitate the implementation process of order model, the main procedures of the model at each time are summarized as follows.Step 1.
Parameters initialization: set the initial value of inspection timesk by choosing randomly or by the prior knowledge.Step 2.
Parameters estimation: estimate the degradation model parametersθ^=μ^,η^,γ^ by MLE according to the nature of IG process. The available estimation samples are obtained from the observed degradation data of the mechanical component till time tk, i.e., X0:k.Step 3.
RUL estimation: based on the real-time estimated parametersθ^=μ^,η^,γ^ in Step 2, we can derive the RUL distribution of the mechanical component at time tk immediately.Step 4.
Applying order model (i.e., CBPO): based on Step3, random lead-time distribution, and speed-up lead-time degree, we can seek the optimal spare ordering time tk0∗ that minimizes the expected cost rate under lower availability threshold by (20).Step 5.
Decision-making: according to the result of Step4, if the optimal spare ordering time tk0∗ is not equal to the next inspection time tk+1, then return to Step 1 and wait for the next inspection time for new decision. Otherwise, stop the procedure and place an order simultaneously. Figure 1 shows the flow chart of the implementation process of order model at each time tk.Figure 1
Flow chart of the implementation process of order model at each timetk.


## 4.2. Order Policy

The order policy is made up of regular order policy and expedited order policy. When the scheduled spare ordering time occurs before the component failure time, regular order policy is triggered. When the scheduled spare ordering time occurs after the component failure time, expedited order policy is triggered. The detailed order policy for the component is depicted in Figure2, where rk is a remaining time, T is the component failure time, Lr is regular lead-time, Le is expedited lead-time, tk is the current inspection time, tko is the scheduled spare ordering time, H is the holding time, and S is the shortage time.Figure 2
Possible order-replacement states of one cycle.State 1.
If the failed component occurs aftertko and the spare is delivered before the component failure, the failed component will not be replaced until the residual lifetime is equal to zero. Let Ps1 denote the probability of State 1 and corresponding CDF can be shown as(10)Ps1=Prtko+Lr<tk+rkand0<tko-tk<rk=∫tko-tk∞∫0tk+u-tkodWtdFRku,where t,u denote integral variable, respectively.State 2.
The failed component occurs aftertko and the spare is not delivered when the component failure occurs. To reduce the wait time after the component failure, we speed up the lead-time. The failed component is replaced by the spare as soon as the spare is delivered. Let Ps2 denote the probability of State 2 and corresponding CDF can be shown as(11)Ps2=Pr⁡tko+Lr>tk+rkand0<tko-tk<rk=∫tko-tk∞∫tk+u-tko∞dWtdFRku,where t,u denote integral variable, respectively.State 3.
The mechanical component failure occurs beforetko; an expedited order is made and the failed component is replaced by the spare as soon as the spare is delivered. Let Ps3 denote the probability of State 3 and corresponding CDF can be shown as(12)Ps3=Pr⁡tko>tk+rk=∫0∞∫0tko-tkdFRktdHu,where t, u denote integral variable, respectively.


## 4.3. Model of Objective Function

The main goal of this section is to build the expected cost rate model under lower availability threshold constraint.


### 4.3.1. Cost Rate Model

(  1) The Expected Replacement Cycle Length. As shown in Figure 2 and formulation of the order-replacement states, a replacement cycle can be formulated as(13)U=tk+rk,iftko+Lr<tk+rkand0<tko-tk<rk,tko+Lr,iftko+Lr>tk+rkand0<tko-tk<rk,tko+Le,iftko>tk+rk,where U is a replacement cycle length.Since holding time only occurs in State 1, the expected holding timeEH during the replacement cycle is(14)EH=Etk+rk-tko-Lr=∫tko-tk∞∫0tk+u-tkotk+u-tko-tdWtdFRku,where t,u denote integral variable, respectively.Furthermore, since the shortage time may occur in State 2 and State 3, the expected shortage timeES during the replacement cycle is(15)ES=Etko+Lr-tk-rk1+κ+ELe=∫tko-tk∞∫tk+u-tko∞tko+t-tk-udWtdFRku1+κ+∫0∞uFRktko-tkhudu,0≤κ<∞where t,u denote integral variable, respectively, and κ denotes speed-up lead-time degree.Based on the above analyses and using expression (13), the expected replacement cycle length can be formulated as(16)EUtko=tk+∫0∞F-Rkrk;X0:kdt+ES=Tf+ES.In other words, the expected replacement cycle length can be denoted by the expected failure lifetime and the expected shortage time.(  2) The Expected Replacement Cycle Cost. As shown in Figure 2, the expected cost per replacement cycle can be derived, which includes inspection cost kρi, replacement cost C, regular ordering cost Cr, regular delivery cost Crd, expedited ordering cost Ce, expedited delivery cost Ced, speed-up lead-time cost κCed, holding cost ρhEH, shortage cost ρsES, and salvage value Csv, where Ce>Cr,Ced>Crd. Therefore, by (10), (11), (12), (14), and (15), the expected cost per replacement cycle EVtko can be expressed as(17)EVtko=k·ρi+C+Cr·Ps1+Ps2+Ce·Ps3+Crd·Ps1+Crd+κ·Ced·Ps2+Ced·Ps3+ρh·EH+ρs·ES-Csv.(  3) Cost Rate Model. Therefore, based on the renewal-reward theory in Ross [41], the expected cost rate is obtained as(18)CRtko=EVtkoEUtko,where tk<tko, EUtko and EVtko are, respectively, given by (16) and (17).


### 4.3.2. Availability

Availability is an important performance index of repairable components and hereafter denoted byA. Based on key renewal theorem in Ross [41], A can be defined to be such that(19)Atko=Erun  timeErun  time+Edown  time=TlTl+ES=11+ES/Tlwhere Erun  time and Edown  time correspond to the expected run time and down time. Erun  time is equal to Tl and Edown  time is equal to ES. Based on (19), we can find that A decreases as ES increases.


### 4.3.3. Cost Rate Model Subject to Availability

Consider the availability as a constraint attached to the cost rate model. Modifying (18), the cost rate model subject to availability is as follows:(20)CRtko=EVtkoEUtko,s.t.A≥A∗where A∗ is a predefined lower availability threshold. Based on (20), the optimal spare ordering time tko∗ can be obtained at each inspection time tk simultaneously. What we aim is to seek tko∗ by minimizing the expected cost rate under lower availability threshold.


## 4.3.1. Cost Rate Model

(  1) The Expected Replacement Cycle Length. As shown in Figure 2 and formulation of the order-replacement states, a replacement cycle can be formulated as(13)U=tk+rk,iftko+Lr<tk+rkand0<tko-tk<rk,tko+Lr,iftko+Lr>tk+rkand0<tko-tk<rk,tko+Le,iftko>tk+rk,where U is a replacement cycle length.Since holding time only occurs in State 1, the expected holding timeEH during the replacement cycle is(14)EH=Etk+rk-tko-Lr=∫tko-tk∞∫0tk+u-tkotk+u-tko-tdWtdFRku,where t,u denote integral variable, respectively.Furthermore, since the shortage time may occur in State 2 and State 3, the expected shortage timeES during the replacement cycle is(15)ES=Etko+Lr-tk-rk1+κ+ELe=∫tko-tk∞∫tk+u-tko∞tko+t-tk-udWtdFRku1+κ+∫0∞uFRktko-tkhudu,0≤κ<∞where t,u denote integral variable, respectively, and κ denotes speed-up lead-time degree.Based on the above analyses and using expression (13), the expected replacement cycle length can be formulated as(16)EUtko=tk+∫0∞F-Rkrk;X0:kdt+ES=Tf+ES.In other words, the expected replacement cycle length can be denoted by the expected failure lifetime and the expected shortage time.(  2) The Expected Replacement Cycle Cost. As shown in Figure 2, the expected cost per replacement cycle can be derived, which includes inspection cost kρi, replacement cost C, regular ordering cost Cr, regular delivery cost Crd, expedited ordering cost Ce, expedited delivery cost Ced, speed-up lead-time cost κCed, holding cost ρhEH, shortage cost ρsES, and salvage value Csv, where Ce>Cr,Ced>Crd. Therefore, by (10), (11), (12), (14), and (15), the expected cost per replacement cycle EVtko can be expressed as(17)EVtko=k·ρi+C+Cr·Ps1+Ps2+Ce·Ps3+Crd·Ps1+Crd+κ·Ced·Ps2+Ced·Ps3+ρh·EH+ρs·ES-Csv.(  3) Cost Rate Model. Therefore, based on the renewal-reward theory in Ross [41], the expected cost rate is obtained as(18)CRtko=EVtkoEUtko,where tk<tko, EUtko and EVtko are, respectively, given by (16) and (17).


## 4.3.2. Availability

Availability is an important performance index of repairable components and hereafter denoted byA. Based on key renewal theorem in Ross [41], A can be defined to be such that(19)Atko=Erun  timeErun  time+Edown  time=TlTl+ES=11+ES/Tlwhere Erun  time and Edown  time correspond to the expected run time and down time. Erun  time is equal to Tl and Edown  time is equal to ES. Based on (19), we can find that A decreases as ES increases.


## 4.3.3. Cost Rate Model Subject to Availability

Consider the availability as a constraint attached to the cost rate model. Modifying (18), the cost rate model subject to availability is as follows:(20)CRtko=EVtkoEUtko,s.t.A≥A∗where A∗ is a predefined lower availability threshold. Based on (20), the optimal spare ordering time tko∗ can be obtained at each inspection time tk simultaneously. What we aim is to seek tko∗ by minimizing the expected cost rate under lower availability threshold.


## 5. Case Study

### 5.1. An Application on CJK1630 Spindle

The spindle is the core component of a computer numerical control (CNC) machine tool, and its running state directly affects the overall quality of the machine tool. The spindle rotation accuracy is one of the important parameters and directly affects the cylindricity of machining parts, axial dimension precision, face shape precision, and surface roughness. The degradation level of the spindle rotation accuracy increases monotonically with the run time. Considering the good properties of the IG process, we describe the degradation process of spindle rotation accuracy by the IG process. The external environment, such as external force, can accelerate the degradation of the spindle rotation accuracy. Therefore, the degradation process of the spindle rotation accuracy can be reasonably described by IG process with covariate effect. We can predict the RUL of the spindle through the available degradation data, which is collected by monitoring the health condition of the spindle in real-time. The mechanical spindle, in practice, is commonly expensive, complex, and highly reliable. Meanwhile, its lead-time is usually not predictable. In order to ensure the completion of the production tasks on time, the availability of the spindle does not fall below the predefined threshold. Above all else, we choose the CJK1630, a CNC lathe, as the subject for the experiment and apply the CBPO to a case study of the mechanical spindle. In order to test the health condition of the mechanical spindle, a CBCM process was carried out by a set of dedicated test platforms each day. The test platform of the mechanical spindle is shown in Figure3 and the corresponding test parameters are shown in Table 1.Table 1
CJK1630 spindle parameters.
Speed (r/min) Radial force (N) Heat-balance temperature (C∘) 1400 140 55Figure 3
The test platform of the mechanical spindle.The inspection times are taken from zero to 90 with the measured frequencyΔt=1, and the failure threshold is chosen to be l=29.8 with the corresponding lifetime T=90. The degradation data of the spindle rotating accuracy with time variation is shown in Table 2 and the corresponding degradation path of the spindle is shown in Figure 4. In order to test the degradation of the spindle rotating accuracy under radial force, we consider that the covariate is Zti=0,1; i.e., Zti=0 indicates that at time ti the environment is “normal” and Zti=1 indicates that at time ti the environment is “external force”, where the classification criterion is according to the experience of the experts.Table 2
The degradation data with time variation.
CBCM time (d) Mean radial error (μm) CBCM time (d) Mean radial error (μm) 0 11.8 21 16.2 1 11.9 22 16.3 2 12.1 23 16.5 3 12.3 24 16.6 4 12.4 25 16.8 5 12.7 26 17.0 6 12.9 27 17.2 7 13.0 28 17.3 8 13.2 29 17.5 9 13.6 30 17.8 10 13.8 40 19.8 11 13.9 45 20.8 12 14.1 50 21.8 13 14.2 55 22.8 14 14.6 60 23.7 15 14.8 65 24.8 16 15.0 70 25.6 17 15.4 75 26.6 18 15.5 80 27.7 19 15.6 85 28.6 20 16.0 90 29.8Figure 4
The degradation path of the mechanical spindle.With the approach presented in Section3.1 and observed degradation data, the estimated degradation model parameters θ^=μ^,η^,γ^ can be obtained by MLE and partial results are selected to show in Table 3.Table 3
Partial parameter estimation results for CJK1630 spindle.
t k μ ^ k η ^ k γ ^ k 30 0.0157 0.0275 1.9820 35 0.0157 0.0276 1.9764 40 0.0156 0.0278 1.9722 45 0.0156 0.0279 1.9690 50 0.0158 0.0274 1.9876 55 0.0157 0.0276 1.9786For validating the effectiveness and applicability of the CBPO, the cost parameters are provided in Table4. Based on historical order practice of the spindle, we suggest the lognormal distribution be used as the distribution of the lead-time in this case study. The PDF of the expedited lead-time can be written as ht=exp⁡-ln⁡t-μh2/2σh2/tσh2π, where μh=2.7, σh=0.1 and the mean of the expedited lead-time is μe=15. Similarly, the PDF of the regular lead-time can be written as wt=exp⁡-ln⁡t-μw2/2σw2/tσw2π, where μw=3.4, σw=0.1 and the mean of the regular lead-time is μe=30. We plot ht and wt in Figure 5 for illustration, respectively. Note that ht and wt used here are obtained by consulting with supplier of the mechanical spindle.Table 4
Cost parameters in the case study.
C C r C e C r d C e d ρ i ρ h ρ s C s v 10 0.3 1 0.1 0.2 0.01 0.05 0.13 0.1Figure 5
The PDF of the lead-time.Next, we provide and discuss the decision-making results for the mechanical spindle in the case study. Specifically, decision-making results at each timetk mainly include the optimal spare ordering time tko∗ and corresponding expected cost rate CRtko∗. The decision-making results of the last 10 CBCM times by the CBPO are shown in Table 5 and the 45th data is selected to plot in Figure 6. From Table 5 and Figure 6, two observations can be made. First, the decision-making result can be real-time updated at each CBCM time. In other words, the CBPO is implemented dynamically according to the real-time health condition of the mechanical spindle. Second, the decision-making process is stopped at time tko=61 since the optimal spare ordering time is equal to the next inspection time tk=61. We will also find that the spare will be ordered at the optimal spare ordering time tko∗=61.Table 5
Decision-making results at the last 10 CBCM times under availability threshold (A∗=0.9700).
t k 45 46 47 48 49 50 51 52 53 54 t k o ∗ 61 61 61 61 61 61 60 60 61 61  C R t k o ∗ 0.1203 0.1204 0.1204 0.1206 0.1207 0.1208 0.1209 0.1210 0.1210 0.1212Figure 6
The decision-making result of the45th condition monitoring time by the CBPO.


### 5.2. Sensitivity Analysis on Critical Parameters

In this section, successive sensitivity analysis will be made forσr,σe,ρs,ρh,κ to explore their influences on the decision-making results. We take the 45th CBCM time as an example, i.e., tk=45. The results are summarized separately in Figures 7–9. Based on these experimental results, the following observation can be drawn.Figure 7
Sensitivity of optimal spare ordering time onσr and σe.Figure 8
Sensitivity of optimal spare ordering time onρs and ρh.Figure 9
Sensitivity of optimal spare ordering time onκ.( 1 ) In Figure 7, we can find that the optimal spare ordering time tko∗ deceases a little as the standard deviation σr increases. The reason for this may lie in that the probabilities that the ordered spare can be delivered earlier or later than its mean lead-time will be along with the variation of the standard deviation and these probabilities are highly correlated with the expected shortage time and holding time among the mechanical spindle's whole life. In the present case, the optimal spare ordering time tko∗decreases since the shortage cost is larger than holding cost. Meanwhile, the optimal spare ordering time tko∗ is almost unchanged as the standard deviation σe increases. The reason is simply that when the inspection time is equal to tk=45, the probability of State 3 is almost zero.( 2 ) In Figure 8, we can find that the optimal spare ordering time tko∗ decreases as the shortage cost per unit time ρs increases. The reason is that a spare unit for replacement should be ordered more early to avoid machine shut down when the mechanical component has a high shortage cost. However, the optimal spare ordering time tko∗ is almost unchanged as the holding cost per unit time ρh increases. This is also reasonable because tko increases under minimum cost rate as ρh increases while corresponding availability decreases. In order to make the availability not below the predefined threshold, one should place an order earlier. Therefore, the final result shows no change.( 3 ) In Figure 9, we can find that the optimal spare ordering time tko∗ first increases and then decreases as speed-up lead-time degree κ increases. Before the inflection point, as κ increases, the shortage time decreases and the availability increases. Since the speed-up cost is less than the cost of shortening the shortage time, tko∗needs to be moved backward. After the inflection point, with the increase of κ, the shortage time further decreases and the availability further increases. Because the speed-up cost at this time has less and less impact on shortening the shortage time, tko∗ needs to be moved forward.All the experimental results are intuitive and match our expectations.


## 5.1. An Application on CJK1630 Spindle

The spindle is the core component of a computer numerical control (CNC) machine tool, and its running state directly affects the overall quality of the machine tool. The spindle rotation accuracy is one of the important parameters and directly affects the cylindricity of machining parts, axial dimension precision, face shape precision, and surface roughness. The degradation level of the spindle rotation accuracy increases monotonically with the run time. Considering the good properties of the IG process, we describe the degradation process of spindle rotation accuracy by the IG process. The external environment, such as external force, can accelerate the degradation of the spindle rotation accuracy. Therefore, the degradation process of the spindle rotation accuracy can be reasonably described by IG process with covariate effect. We can predict the RUL of the spindle through the available degradation data, which is collected by monitoring the health condition of the spindle in real-time. The mechanical spindle, in practice, is commonly expensive, complex, and highly reliable. Meanwhile, its lead-time is usually not predictable. In order to ensure the completion of the production tasks on time, the availability of the spindle does not fall below the predefined threshold. Above all else, we choose the CJK1630, a CNC lathe, as the subject for the experiment and apply the CBPO to a case study of the mechanical spindle. In order to test the health condition of the mechanical spindle, a CBCM process was carried out by a set of dedicated test platforms each day. The test platform of the mechanical spindle is shown in Figure3 and the corresponding test parameters are shown in Table 1.Table 1
CJK1630 spindle parameters.
Speed (r/min) Radial force (N) Heat-balance temperature (C∘) 1400 140 55Figure 3
The test platform of the mechanical spindle.The inspection times are taken from zero to 90 with the measured frequencyΔt=1, and the failure threshold is chosen to be l=29.8 with the corresponding lifetime T=90. The degradation data of the spindle rotating accuracy with time variation is shown in Table 2 and the corresponding degradation path of the spindle is shown in Figure 4. In order to test the degradation of the spindle rotating accuracy under radial force, we consider that the covariate is Zti=0,1; i.e., Zti=0 indicates that at time ti the environment is “normal” and Zti=1 indicates that at time ti the environment is “external force”, where the classification criterion is according to the experience of the experts.Table 2
The degradation data with time variation.
CBCM time (d) Mean radial error (μm) CBCM time (d) Mean radial error (μm) 0 11.8 21 16.2 1 11.9 22 16.3 2 12.1 23 16.5 3 12.3 24 16.6 4 12.4 25 16.8 5 12.7 26 17.0 6 12.9 27 17.2 7 13.0 28 17.3 8 13.2 29 17.5 9 13.6 30 17.8 10 13.8 40 19.8 11 13.9 45 20.8 12 14.1 50 21.8 13 14.2 55 22.8 14 14.6 60 23.7 15 14.8 65 24.8 16 15.0 70 25.6 17 15.4 75 26.6 18 15.5 80 27.7 19 15.6 85 28.6 20 16.0 90 29.8Figure 4
The degradation path of the mechanical spindle.With the approach presented in Section3.1 and observed degradation data, the estimated degradation model parameters θ^=μ^,η^,γ^ can be obtained by MLE and partial results are selected to show in Table 3.Table 3
Partial parameter estimation results for CJK1630 spindle.
t k μ ^ k η ^ k γ ^ k 30 0.0157 0.0275 1.9820 35 0.0157 0.0276 1.9764 40 0.0156 0.0278 1.9722 45 0.0156 0.0279 1.9690 50 0.0158 0.0274 1.9876 55 0.0157 0.0276 1.9786For validating the effectiveness and applicability of the CBPO, the cost parameters are provided in Table4. Based on historical order practice of the spindle, we suggest the lognormal distribution be used as the distribution of the lead-time in this case study. The PDF of the expedited lead-time can be written as ht=exp⁡-ln⁡t-μh2/2σh2/tσh2π, where μh=2.7, σh=0.1 and the mean of the expedited lead-time is μe=15. Similarly, the PDF of the regular lead-time can be written as wt=exp⁡-ln⁡t-μw2/2σw2/tσw2π, where μw=3.4, σw=0.1 and the mean of the regular lead-time is μe=30. We plot ht and wt in Figure 5 for illustration, respectively. Note that ht and wt used here are obtained by consulting with supplier of the mechanical spindle.Table 4
Cost parameters in the case study.
C C r C e C r d C e d ρ i ρ h ρ s C s v 10 0.3 1 0.1 0.2 0.01 0.05 0.13 0.1Figure 5
The PDF of the lead-time.Next, we provide and discuss the decision-making results for the mechanical spindle in the case study. Specifically, decision-making results at each timetk mainly include the optimal spare ordering time tko∗ and corresponding expected cost rate CRtko∗. The decision-making results of the last 10 CBCM times by the CBPO are shown in Table 5 and the 45th data is selected to plot in Figure 6. From Table 5 and Figure 6, two observations can be made. First, the decision-making result can be real-time updated at each CBCM time. In other words, the CBPO is implemented dynamically according to the real-time health condition of the mechanical spindle. Second, the decision-making process is stopped at time tko=61 since the optimal spare ordering time is equal to the next inspection time tk=61. We will also find that the spare will be ordered at the optimal spare ordering time tko∗=61.Table 5
Decision-making results at the last 10 CBCM times under availability threshold (A∗=0.9700).
t k 45 46 47 48 49 50 51 52 53 54 t k o ∗ 61 61 61 61 61 61 60 60 61 61  C R t k o ∗ 0.1203 0.1204 0.1204 0.1206 0.1207 0.1208 0.1209 0.1210 0.1210 0.1212Figure 6
The decision-making result of the45th condition monitoring time by the CBPO.


## 5.2. Sensitivity Analysis on Critical Parameters

In this section, successive sensitivity analysis will be made forσr,σe,ρs,ρh,κ to explore their influences on the decision-making results. We take the 45th CBCM time as an example, i.e., tk=45. The results are summarized separately in Figures 7–9. Based on these experimental results, the following observation can be drawn.Figure 7
Sensitivity of optimal spare ordering time onσr and σe.Figure 8
Sensitivity of optimal spare ordering time onρs and ρh.Figure 9
Sensitivity of optimal spare ordering time onκ.( 1 ) In Figure 7, we can find that the optimal spare ordering time tko∗ deceases a little as the standard deviation σr increases. The reason for this may lie in that the probabilities that the ordered spare can be delivered earlier or later than its mean lead-time will be along with the variation of the standard deviation and these probabilities are highly correlated with the expected shortage time and holding time among the mechanical spindle's whole life. In the present case, the optimal spare ordering time tko∗decreases since the shortage cost is larger than holding cost. Meanwhile, the optimal spare ordering time tko∗ is almost unchanged as the standard deviation σe increases. The reason is simply that when the inspection time is equal to tk=45, the probability of State 3 is almost zero.( 2 ) In Figure 8, we can find that the optimal spare ordering time tko∗ decreases as the shortage cost per unit time ρs increases. The reason is that a spare unit for replacement should be ordered more early to avoid machine shut down when the mechanical component has a high shortage cost. However, the optimal spare ordering time tko∗ is almost unchanged as the holding cost per unit time ρh increases. This is also reasonable because tko increases under minimum cost rate as ρh increases while corresponding availability decreases. In order to make the availability not below the predefined threshold, one should place an order earlier. Therefore, the final result shows no change.( 3 ) In Figure 9, we can find that the optimal spare ordering time tko∗ first increases and then decreases as speed-up lead-time degree κ increases. Before the inflection point, as κ increases, the shortage time decreases and the availability increases. Since the speed-up cost is less than the cost of shortening the shortage time, tko∗needs to be moved backward. After the inflection point, with the increase of κ, the shortage time further decreases and the availability further increases. Because the speed-up cost at this time has less and less impact on shortening the shortage time, tko∗ needs to be moved forward.All the experimental results are intuitive and match our expectations.


## 6. Conclusions and Future Research

In this paper, CBPO is proposed for a mechanical component subject to IG degradation process with covariate effect. The CBPO, depending on RUL, random lead-time, speed-up lead-time degree, and availability threshold, can be real-time determined and make an accurate prediction for spare ordering time according to the actual health condition of the mechanical component in operation. The CBPO is illustrated as feasible and practical through a CNC lathe mechanical spindle case study. The sensitivity analysis results are intuitive and match our expectations. By using this model, the optimal spare ordering time with minimum expected cost rate under the lower availability threshold constraint can be calculated easily and effectively for enterprises.Inspection accuracy of the degradation characteristics for the mechanical components plays an important role in the CBPO. Engineers or managers should pay greater attention to improving the inspection accuracy. The higher inspection accuracy of the degradation characteristics for the mechanical components, the higher the prediction accuracy of spare ordering time. On the other hand, with reasonable speed-up lead-time degree, not only can the availability be improved, but maintenance costs can also be reduced. The approach proposed in this paper can provide guidance for some systems beyond mechanical components, which are subject to IG degradation process or other degradation processes. The approach can prompt reasonable arrangements for production, improve product quality and service quality, reduce maintenance costs, and enhance market competitiveness.Although the CBPO can provide satisfactory results, there are still some research points that need to be investigated further in the future and future research may focus on the following topics.(a) Multi-spare ordering model based on RUL for systems following stochastic degradation process.(b) The CBPO under variable working condition that is challenging and worthwhile.(c) Spare ordering model based on the joint effect of aging and cumulative degradation, which is more practical.



---

*Source: 9734189-2018-09-10.xml*