Modeling condition-based maintenance and replacement strategies for an imperfect production-inventory system

Abstract

In this article, we consider an imperfect production-inventory system which produces a single type of product to meet the constant demand. The system deteriorates stochastically with usage and the deterioration process is modeled by a non-stationary gamma process. The production process is imperfect which means that the system produces some non-conforming items and the product quality depends on the degradation level of the production system. To prevent the system from deteriorating worse and improve the product quality, preventive maintenance is performed when the level of the system degradation reaches a certain threshold. However, the preventive maintenance is imperfect which cannot restore the system as good as new. Hence, the aging system will be replaced by a new one after some production cycles. The preventive maintenance cost, the replacement cost, the production cost, the inventory holding cost and the penalty cost of lost sales are considered in this article. The objective is to minimize the total cost per unit item which depends on two decision variables: the preventive maintenance threshold and the time at which the system is replaced. We derive the explicit expression of the total cost per unit item and the optimal joint policy can be obtained numerically. An illustrative example and sensitivity analysis are given to demonstrate the proposed model.

Introduction

For production-inventory systems, the economic production quantity (EPQ) model is applied to reduce the production- and inventory-related costs. In the earliest EPQ models, a common limited assumption is that machines never fail during a production run. However, most production systems suffer increasing wear or deteriorate with usage or aging process. Eventually, they fail from this deterioration without preventive maintenance (PM). The machine breakdowns and corrective maintenance have a great influence on production planning and inventory management. Therefore, an appropriate PM strategy should be selected to improve the system reliability with minor sacrifices during production run.

The production models with PM have been intensively studied by some researchers. Meller and Kim ¹ considered a production system and assumed that periodic PM can decrease the production operation’s failure rate in the long run. They used an embedded absorbing Markov chain analysis to determine the optimal inventory level which triggers the PM. Abboud ² considered a single machine production-inventory system and assumed that the times-to-failure and the repair times are geometrically distributed. A Markov chain model was developed and an efficient algorithm was proposed to find the economic manufacturing quantities. Chung et al. ³ extended the above model by considering the deteriorating items with stochastic machine unavailable time and shortage. Widyadana and Wee ⁴ developed two EPQ deteriorating item models of uniform distribution and exponential distribution of corrective and maintenance times. Wee and Widyadana ⁵ developed an EPQ model with stochastic PM time and rework process. Recoverable items (defined as items that do not meet quality standards but can be reprocessed) are reprocessed by the last in first out (LIFO) rule. Then, Wee and Widyadana ⁶ revisited the above model using the first in first out (FIFO) rule for rework process. With the development of sensor and information technologies, the degradation level of systems can be monitored to estimate the systems’ state and facilitate the prediction of failures.^7,8 Hence, condition-based maintenance (CBM) has been widely used in maintenance area due to its effectiveness and efficiency. However, few researches have been done to investigate the joint model for EPQ and CBM. Jafari and Makis ⁹ studied joint optimization of EPQ and CBM. The problems are formulated and solved in the semi-Markov decision process framework. Peng and Van Houtum ¹⁰ proposed a model to determine the production of lot-sizing and CBM policy by assuming a continuous degradation process of the system. Renewal theory is employed to compute the average long-run cost rate.

Note that the PM in the above literature is assumed to be perfect, that is, it can restore the machine to the state as good as new. However, in practice, it is not always the case. For more realistic, some researchers considered imperfect PM in production-inventory system models. Sheu and Chen ¹¹ developed an integrated model for the joint optimization of EPQ and level of PM. They assumed that after performing PM, the aging of the system is reduced in proportion to the PM level. Wang et al. ¹² extended the model by considering inspections which are used to see whether the process is in control. El-Ferik ¹³ assumed that the age-based PM is imperfect and the machine will be replaced at the Nth maintenance. The objective was to obtain the optimal PM interval and replacement policy such that the long-run average cost rate is minimized. Ouaret et al. ¹⁴ simultaneously determined the optimal production rate and replacement policy for a failure-prone manufacturing system by considering the imperfect maintenance. Semi-Markov processes and numerical methods were used to model the dynamics of the system and obtain the optimal control policies. Cheng et al. ¹⁵ studied a production system which consists of two machines and one buffer. The machine tends to deteriorate faster as the number of maintenance increases. An integrated model was developed for the joint determination of buffer size and PM interval. More recently, Cheng et al. ¹⁶ employed geometric process to describe the effects of imperfect maintenance. The aging machine will be replaced after N failures. An algorithm is proposed to obtain the optimal buffer size and replacement policy.

It is well known that the deterioration process influences not only the reliability of the production system but also the quality of products. Indeed, the machine state plays an important role in controlling the quality of produced items. ¹⁷ Ben-Daya ¹⁸ developed an EPQ model linking quality and PM. The production process may shift from “in-control” state to “out-of-control” state after random period of time and then produces some non-conforming items. The author jointly determined the EPQ and PM level and proved that performing PM can reduce the quality-related costs. Peter Chiu et al. ¹⁹ concerned with the algebraic derivation for the production lot-sizing problem considering random defective rate, rework process and backlogging. Radhoui et al. ¹⁷ studied a joint quality control and PM policy for an imperfect production system. The PM level was defined as the proportion of non-conforming items produced. Simulation and experimental design were used to determine the optimal buffer size and PM level simultaneously. Yang and Zeng ²⁰ proposed an equipment maintenance policy by integrating periodic equipment inspection and product quality control. Maintenance action is chosen based on the inspection result of the product quality shift and equipment degradation state. The optimal maintenance policy is obtained by genetic algorithm. Sana ²¹ considered a production-inventory model in an imperfect production process with time-varying demand. The author assumed that the unit production cost is a function of product reliability parameter and production rate which are two decision variables of the model. The optimal solutions are obtained using the Euler–Lagrange method. Pal et al. ²² assumed that non-conforming items are produced in the both “in-control” and “out-of-control” states. The probability of non-conforming items produced is smaller in the “in-control” state than in the “out-of-control” state. After regular production run time, the non-conforming items are reworked and the PM is performed. The optimal buffer inventory and production run time are obtained to minimize the unit production cost. Other production-maintenance models considering product quality can be seen in Chouikhi et al., ²³ Liao, ²⁴ and Tsao ²⁵ and Jeang. ²⁶

In the literature mentioned above, the states of production systems are generally divided into two classes: “in-control” state and “out-of-control” state. When the system stays in the “in-control” state, no or less non-conforming items are produced. While in the “out-of-control” state, the system produces more non-conforming items. This assumption reflects the fact that the degradation of systems influences the product quality. However, in practice, only two states may not be able to describe the degradation level of the system accurately. For example, the tool wear is a gradual deteriorating process. According to ISO 3685, ²⁷ based on the size of flank wear monitored (denoted by $V_B$ ), five states are generally recommended for tool wear: sharp ( $V_B\le 0.1\mathrm{mm}$ ), normal wear ( $0.1\mathrm{mm}<V_B\le 0.2\mathrm{mm}$ ), micro fracture ( $0.2\mathrm{mm}<V_B\le 0.3\mathrm{mm}$ ), macro wear ( $0.3\mathrm{mm}<V_B\le 0.4\mathrm{mm}$ ) and breakage ( $V_B>0.4\mathrm{mm}$ ). Since the tool state has a great impact on product quality, the proportion of non-conforming items produced increases as the tool wear level increases. To model the quality loss properly, a more accurate mathematical model needs to be built when addressing the combination of PM strategy and quality control. From the author’s literature search, no research on EPQ models considering quality loss with multi-state production process has been done, although it may be more realistic for some production systems in practice. In addition, age-based maintenance (ABM) policy is used in most of the above researches. However, when investigating production-inventory models with multi-state, CBM policy is more appropriate, although it may make the optimization model be more complicated.

In this article, CBM and replacement strategies are proposed for a production-inventory system taking into account degradation-level-dependent quality. The production system deteriorates gradually, and its deterioration process is modeled by gamma process which is widely used in maintenance models for deteriorating systems. ²⁸ The proportion of non-conforming items depends on the system’s degradation level monitored, that is, the more non-conforming items are produced when the system deteriorates to a higher level. PM actions are performed to control the system degradation and improve the product quality when the deterioration level reaches the threshold which is a decision variable of the model. Further assumed that the PM is imperfect, that is, it cannot bring the system to its initial state. Consequently, the aging system will be replaced after N maintenance actions, where N is another decision variable of the model. Our aim is to minimize the total cost per unit item produced consisting of PM, replacement, production, inventory holding and lost sales costs. The maintenance optimization model is formally developed and the optimal joint strategy is obtained numerically. The rest of this article is organized as follows. Section “Assumptions and notations” presents model assumptions and notations used. Section “Model description” provides a detailed model description. Section “Formulation of the model” is dedicated to the derivation of the explicit expression of the total cost per unit item. Section “An illustrative example and discussions” gives an illustrative example to obtain the optimal joint strategy numerically and carry out some sensitivity analysis. Finally, conclusions are drawn in the last section.

Assumptions and notations

The following assumptions and notations are used to develop the model.

Assumptions

• A1. At the beginning, a new production system is installed which produce one type of product to meet the constant demand.

• A2. The production system deteriorates with use and age and its degradation behavior is modeled by a stochastic process. The system failure occurs when the degradation level exceeds the predetermined threshold.

• A3. The production process is imperfect, that is, the system produces some defective or non-conforming items which cannot be reworked and are scrapped. In addition, the proportion of non-conforming items depends on the degradation level of the system, that is, the system produces non-conforming items at a higher rate as the degradation level increases.

• A4. The system is monitored continuously and perfectly, which reveals the true degradation level in real time. The PM is performed when the degradation level reaches the PM threshold. It is assumed that the PM is imperfect which cannot restore the system to the state as good as new. More details are available in the following section.

• A5. Since the PM is imperfect, the production system will be replaced by an identical new one after some production cycles.

• A6. The production rate of conforming items is greater than the demand rate. The initial inventory level is zero. In this article, the no-resumption (NR) policy ²⁹ is adopted which means that after the PM or replacement, a new production cycle starts when the on-hand inventory is depleted.

• A7. When demand is greater than the stock during PM time, shortage occurs. We assume shortage items are lost sales.

• A8. The model is studied over infinite planning horizon.

Notations

• P production rate of the system;

• D constant demand rate;

• $\alpha$  shape parameter of gamma process;

• $\beta$  scale parameter of gamma process;

• a ratio of the geometric process;

• $X(t)$  degradation level of the production system at time t;

• $X_i^{*}$  random variable denoting degradation state immediately after the ith PM;

• ${\xi }_i(1\le i\le M)$  thresholds of deteriorating process, where ${\xi }_M$ is the failure threshold;

• $b_i(1\le i\le M)$  production rate of non-conforming items when the degradation level $X(t)\in ({\xi }_{i-1},{\xi }_i]$ ;

• ${\tau }_x(x>0)$  time in which, for the first time, the system’s degradation exceeds the level x;

• $F_{{\tau }_x}$  probability distribution function of ${\tau }_x$ ;

• $f_{{\tau }_x}$  probability density function of ${\tau }_x$ ;

• $g_{X^{*}}^{(k)}$  probability density function of $X_i^{*}$ if PM threshold is ${\xi }_k$ ;

• $p\mathrm{and}q$  two parameters of Beta distribution;

• $T_i$  length of the ith production cycle;

• $t_{\mathrm{ij}}$  time interval between ${\tau }_{{\xi }_{j-1}}$ and ${\tau }_{{\xi }_j}$ in the ith production cycle;

• $T_{D_i}$  time to deplete all on-hand inventory after the ith production cycle;

• $T_{\mathrm{pm}}$  random variable denoting the duration of PM;

• $T_r$  random variable denoting the duration of replacement;

• $\phi$  probability density function of $T_{\mathrm{pm}}$ ;

• $\varphi$  probability density function of $T_r$ ;

• $\lambda$  parameter of exponential distribution function $\phi$ ;

• $\mu$  parameter of exponential distribution function $\varphi$ ;

• $c_{\mathrm{pm}}$  cost of PM per unit time;

• $c_r$  cost of replacement;

• $c_p$  cost of production per product including cost of raw materials and labors and so on;

• $c_h$  inventory holding cost for conforming items per unit per unit time;

• $c_l$  penalty cost of lost sales per item;

• $k(1\le k<M)$  decision variable which means the PM will be performed when $X(t)$ exceeds level ${\xi }_k$ ;

• N decision variable which means the system will be replaced after N production cycle;

• $C_{\mathrm{PM}}$  PM cost incurred in a renewal cycle;

• $C_R$  replacement cost incurred in a renewal cycle;

• $C_P$  production cost incurred in a renewal cycle;

• $C_{\mathrm{IH}}$  inventory holding cost incurred in a renewal cycle;

• $C_{\mathrm{LS}}$  penalty cost of lost sales incurred in a renewal cycle;

• $C(k,N)$  long-run average cost per conforming item (objective function).

Model description

Deterioration modeling

Consider a production system which deteriorates due to usage and age. The degradation evolution is modeled by a stochastic process. Let $X(t)$ denote the deterioration level of the system at time t. In this article, we assume that $X(t)$ is a gamma process which is widely used for modeling different degradation phenomena such as corrosion, erosion, crack growth, wear of structural components and so on. ²⁸ According to the definition of gamma process, $X(t)$ has a gamma density function with shape parameter $\alpha t$ and scale parameter $\beta$ which is given by

$$f_{\alpha t,\beta }(x)=\frac{{\beta }^{\alpha t}{x}^{\alpha t-1}}{\Gamma (\alpha t)}{e}^{-\beta x},x\ge 0$$

where $\Gamma (·)$ is the gamma function

$$\Gamma (a)=\underset{0}{\overset{\infty }{\int }}{u}^{a-1}{e}^{-u}\mathrm{du},a>0$$

The average deterioration rate is $m=\alpha /\beta$ and its variance is ${\sigma }^2=\alpha /{\beta }^2$ . Let ${\tau }_x$ denote the time in which, for the first time, the system’s degradation exceeds the level x, and then its distribution function $F_{{\tau }_x}$ has the following expression

$$F_{{\tau }_x}(t)=P({\tau }_x<t)=P(X(t)>x)=\frac{\Gamma (\alpha t,x\beta )}{\Gamma (\alpha t)},t>0$$

where $\Gamma (a,b)={\int }_b^{\infty }{u}^{a-1}{e}^{-u}\mathrm{du},a>0,b\ge 0$ .

The density function $f_{{\tau }_x}$ is given by

$$f_{{\tau }_x}(t)=\frac{\partial F_{{\tau }_x}(t)}{\partial t}=\frac{\alpha }{\Gamma (\alpha t)}{\displaystyle \underset{x\beta }{\overset{\infty }{\int }}\{\log (z)-\psi (\alpha t)\}{z}^{\alpha t-1}{e}^{-z}dz},t>0$$

where $\psi (a)=(\Gamma \prime (a)/\Gamma (a))=((\partial \log \Gamma (a))/\partial a)$ is the digamma function.

Let $y(<x)$ denote another degradation level and ${\tau }_y$ be the time in which the degradation level exceeds y. In the following model formulation, the distribution function or survival function of the random variable ${\tau }_x-{\tau }_y$ will be used. Note that ${\tau }_x-{\tau }_y\ne {\tau }_{x-y}$ due to the “overshoot behavior” of the gamma process. ³⁰ As given in Bérenguer et al., ³¹ the survival function of ${\tau }_x-{\tau }_y$ is as follows

$${\overline{F}}_{{\tau }_x-{\tau }_y}(t)=\{\begin{array}{c}-\int \underset{y<u<x,0<u+v<y,0<v}{\int }(\underset{0}{\overset{\infty }{\int }}{f}_{\alpha w,\beta }(u)\mathrm{dw})\frac{\partial f_{\alpha t,\beta }(v)}{\partial t}\mathrm{dudv},v>u>0;\hfill \\ \multicolumn{1}{c}{\underset{t}{\overset{\infty }{\int }}\underset{0}{\overset{v-t}{\int }}\alpha \underset{0}{\overset{x}{\int }}{f}_{\alpha u,\beta }(w)\mathrm{dw}\left(\underset{x-v}{\overset{\infty }{\int }}\frac{e^{-\beta z}}{z}\mathrm{dz}\right)\mathrm{dudv},y=x>0}\end{array}$$

(1)

Since expression (1) is too complicated for numerical computations, Bérenguer et al. ³¹ gave an approximation of the distribution function of ${\tau }_x-{\tau }_y$ as the distribution function of ${\tau }_{x-y-(1/2)\beta }$ when $x-y-(1/2)\beta >0$ . Therefore, $F_{{\tau }_x-{\tau }_y}$ and $f_{{\tau }_x-{\tau }_y}$ can be approximated as follows

$$F_{{\tau }_x-{\tau }_y}(t)\approx F_{{\tau }_{x-y-\frac{1}{2}\beta }}(t)=\frac{\Gamma (\alpha t,\beta (x-y)-(1/2))}{\Gamma (\alpha t)}$$

(2)

$$f_{{\tau }_x-{\tau }_y}(t)\approx f_{{\tau }_{x-y-\frac{1}{2}\beta }}(t)=\frac{\alpha }{\Gamma (\alpha t)}{\displaystyle \underset{\beta (x-y)-\frac{1}{2}}{\overset{\infty }{\int }}\{\log (z)-\psi (\alpha t)\}{z}^{\alpha t-1}{e}^{-z}dz}$$

(3)

Imperfect production process

In this article, the production process is assumed to be imperfect, that is, some produced items may be non-conforming. Furthermore, it is reasonable to assume that the rate of non-conforming items depends on the degradation level of the system, that is, the rate increases as the degradation level increases. Assume that there are M thresholds of the deterioration process, namely, ${\xi }_1<{\xi }_2<\cdots <{\xi }_M$ , where ${\xi }_M$ is the failure threshold of the system. When the degradation level of the production system $X(t)\in ({\xi }_{j-1},{\xi }_j]$ , $j=1,2,\dots ,M$ , the production rate of non-conforming items is assumed to be $b_j$ which is non-deceasing with respect to j.

Imperfect PM

According to assumption A4, a CBM is applied to this production system based on the information obtained from continuous monitoring. If the degradation level $X(t)$ exceeds the PM threshold ${\xi }_k$ , PM action starts immediately, where $k(1\le k<M)$ is a decision parameter. The PM cannot bring the system to the initial state. The meaning of imperfect PM is twofold: (1) there may be residual damage after PM and (2) the system tends to deteriorate faster as the number of maintenance increases.

Here, we assume that the distribution function of the degradation state immediately after the ith PM $X_i^{*}$ is a Beta distribution $g_{X^{*}}^{(k)}$ constrained to the interval $[0,{\xi }_k]$ and $g_{X^{*}}^{(k)}$ is given by Liao et al. ³²

$$g_{X^{*}}^{(k)}(x)=\frac{1}{{\xi }_k}\frac{\Gamma (p+q)}{\Gamma (p)\Gamma (q)}{\left(\frac{x}{{\xi }_k}\right)}^{p-1}{(1-\frac{x}{{\xi }_k})}^{q-1}·1_{\{0\le x\le {\xi }_k\}}$$

(4)

where the parameters $p>0$ and $q>0$ can be estimated using the Maximum Likelihood Estimators (MLE).

Remark 1

From the definition of Beta distribution, if the PM threshold is ${\xi }_k$ the degradation state $X_i^{*}$ after the ith PM is distributed within the interval $[0,{\xi }_k]$ . It can be seen that the effect of PM depends on the state at which PM is performed, that is, the maintenance effect would be better if the PM is performed at a better state. Actually, ${\mathrm{EX}}_i^{*}=(p{\xi }_k/p+q)$ and $\mathrm{Var}(X_i^{*})=((\mathrm{pq}{\xi }_k^2)/(p+q{)}^2(p+q+1))$ .

In addition, the system is assumed to tend to deteriorate faster as the number of PM increases. This phenomenon is common in practice. For example, the amplitude of the hydraulic coupler can reflect its deterioration states and the growth of the amplitude is accelerated after each repair (see the data in Tan et al. ³³ ). Here, we employ the accelerated deteriorating gamma process (ADGP) to model this stochastic process which is first introduced by Cheng et al. ¹⁵

The ith production cycle is referred to as the time duration from the beginning of the production after the (i − 1)th PM until the beginning of the ith PM action or replacement. Let ${\tau }_x^{(i)}$ denote the time in which, for the first time, the degradation exceeds the level x during the ith production cycle and ${\tau }_x^{(1)}={\tau }_x$ . We assume that ${\tau }_x^{(i)}$ is stochastically decreasing with respect to i which means that the deterioration process is accelerated after each PM. Furthermore, $\{{\tau }_x^{(i)},i=1,2,\dots \}$ forms a decreasing geometric process.^34,35 Let $F_{{\tau }_x^{(i)}}^{(i)}$ be the distribution function of ${\tau }_x^{(i)}$ , where $F_{{\tau }_x^{(1)}}^{(1)}=F_{{\tau }_x}$ , according to the definition of geometric process, we have the following

$$F_{{\tau }_x^{(i)}}^{(i)}(t)=F_{{\tau }_x^{(1)}}^{(1)}(a^{i-1}t)=F_{{\tau }_x}(a^{i-1}t)=\frac{\Gamma (a^{i-1}\alpha t,x\beta )}{\Gamma (a^{i-1}\alpha t)}$$

(5)

where a is the ratio of geometric process and $a\ge 1$ .

Consequently, the density function $f_{{\tau }_x^{(i)}}^{(i)}$ can be given by

$$f_{{\tau }_x^{(i)}}^{(i)}(t)=\frac{\partial F_{{\tau }_x^{(i)}}^{(i)}(t)}{\partial t}=\frac{a^{i-1}\alpha }{\Gamma (a^{i-1}\alpha t)}{\displaystyle \underset{x\beta }{\overset{\infty }{\int }}\{\log (z)-\psi (a^{i-1}\alpha t)\}{z}^{a^{i-1}\alpha t-1}{e}^{-z}dz,}t>0$$

(6)

Similar to expressions (2) and (3), $F_{{\tau }_x^{(i)}-{\tau }_y^{(i)}}^{(i)}$ and $f_{{\tau }_x^{(i)}-{\tau }_y^{(i)}}^{(i)}$ can be approximated as follows

$$F_{{\tau }_x^{(i)}-{\tau }_y^{(i)}}^{(i)}(t)\approx F_{{\tau }_{x-y--\frac{1}{2}\beta }^{(i)}}^{(1)}(a^{i-1}t)=\frac{\Gamma (a^{i-1}\alpha t,\beta (x-y)-(1/2))}{\Gamma (a^{i-1}\alpha t)}$$

(7)

$$f_{{\tau }_x^{(i)}-{\tau }_y^{(i)}}^{(i)}(t)\approx f_{{\tau }_{x-y-\frac{1}{2}\beta }^{(i)}}^{(1)}(t)=\frac{a^{i-1}\alpha }{\Gamma (a^{i-1}\alpha t)}{\displaystyle \underset{\beta (x-y)-\frac{1}{2}}{\overset{\infty }{\int }}\{\log (z)-\psi (a^{i-1}\alpha t)\}{z}^{a^{i-1}\alpha t-1}{e}^{-z}dz}$$

(8)

Remark 2

According to expression (5), the shape parameter and the scale parameter of the gamma process after the ith maintenance are $a^i$ and $\beta$ , respectively. Therefore, the average deterioration rate in the ith production cycle is $m_i=a^i/\beta$ and the variance is ${\sigma }_i^2=a^i/{\beta }^2$ . When the ratio $a>1$ , $m^i$ increases with respect to i which implies that the deterioration process is accelerated.

By assuming that the PM threshold is ${\xi }_4$ , Figure 1(a) shows the deterioration process of the system during N production cycles; Figure 1(b) illustrates the inventory level of conforming items. As shown in Figure 1, when the degradation level $X(t)\in ({\xi }_{j-1},{\xi }_j](j=1,2,3,4)$ , the accumulating rate of the inventory of conforming items is $(1-b_j)P-D$ . It can be seen from Figure 1(b) that the duration of the first PM $T_{\mathrm{pm}}$ is shorter than the depleting time of on-hand inventory $T_{D_1}$ ; thus, the second production cycle starts until the inventory is exhausted. The duration of the second PM is longer than the depleting time of inventory $T_{D_2}$ and then shortage occurs. After N production cycles, the aging system is replaced by an identical new one.

Figure 1.

The system degradation level and inventory level versus time

Formulation of the model

According to renewal reward theorem, the total expected quantity of conforming items produced can be derived by dividing the total expected cost per renewal cycle to the total expected quantity produced in a renewal cycle as the following

$$C(k,N)=\frac{E(\mathrm{total\; cost\; per\; renewal\; cycle})}{E(\mathrm{total\; quantity\; of\; conforming\; items\; produced\; per\; renewal\; cycle})}=\frac{E{C}_{Total}}{EQ}$$

(9)

Expected total quantity of conforming items produced in a renewal cycle

Let $Q_i(1\le i\le N)$ denote the quantity of conforming items produced in the ith production cycle within a renewal cycle, thus

$$\mathrm{EQ}=\sum _{i=1}^{N}E{Q}_i$$

(10)

According to assumption A3, the production rate of non-conforming items depends on the system’s degradation level, that is, the rate is $b_j(1\le j\le k)$ when the degradation level satisfies ${\xi }_{j-1}<X(t)\le {\xi }_j$ (here, we set ${\xi }_0=0$ ). Let $t_{\mathrm{ij}}={\tau }_{{\xi }_j}^{(i)}-{\tau }_{{\xi }_{j-1}}^{(i)}$ , where ${\tau }_{{\xi }_j}^{(i)}$ is the time at which, for the first time, the system’s degradation exceeds the level ${\xi }_j$ during the ith production cycle. The random variable $t_{\mathrm{ij}}$ denotes the production time in which the production rate of non-conforming items is $b_j$ . The distribution function of $t_{\mathrm{ij}}$ is $F_{{\tau }_{{\xi }_j}^{(i)}-{\tau }_{{\xi }_{j-1}}^{(i)}}^{(i)}$ and the expectation of it can be given by

$$E{t}_{\mathrm{ij}}=E[{\tau }_{{\xi }_j}^{(i)}-{\tau }_{{\xi }_{j-1}}^{(i)}]=\underset{0}{\overset{\infty }{\int }}{\overline{F}}_{{\tau }_{{\xi }_j}^{(i)}-{\tau }_{{\xi }_{j-1}}^{(i)}}^{(i)}(t)\mathrm{dt}$$

where ${\overline{F}}_{{\tau }_{{\xi }_j}^{(i)}-{\tau }_{{\xi }_{j-1}}^{(i)}}^{(i)}=1-F_{{\tau }_{{\xi }_j}^{(i)}-{\tau }_{{\xi }_{j-1}}^{(i)}}^{(i)}$ can be approximated by equation (7).

Based on the above discussion, we have the following

$$E{Q}_1=\sum _{j=1}^k(1-b_j)P·E{t}_{1j}$$

(11)

When $i\ge 2$ , the calculation of $E{Q}_i$ is slightly different from the calculation of $E{Q}_1$ due to the imperfect PM which cannot bring the system’s degradation state to zero but to a better state than the PM threshold. Let the random variable $X_{i-1}^{*}$ denote the degradation state immediately after the (i − 1)th PM whose distribution is constrained to interval $(0,{\xi }_k)$ . Conditioning on $X_{i-1}^{*}$ and applying the law of total probability, $E{Q}_i(i\ge 2)$ can be formulated as follows

$$E{Q}_i={\displaystyle \sum _{n=1}^{k}E(Q_i|{\xi }_{n-1}<X_{i-1}^{\ast }\le {\xi }_n)·P({\xi }_{n-1}<X_{i-1}^{\ast }\le {\xi }_n)}={\displaystyle \sum _{n=1}^k{\displaystyle \sum _{j=n}^k{\displaystyle \underset{{\xi }_{n-1}}{\overset{{\xi }_n}{\int }}(1-b_j)P·E{t}_{ij}·g_{X^{\ast }}^{(k)}(x)dx}}}$$

(12)

$$E{Q}_i={\displaystyle \sum _{n=1}^k{\displaystyle \underset{{\xi }_{n-1}}{\overset{{\xi }_n}{\int }}\left[(1-b_n)P·E{t}_{in}+{\displaystyle \sum _{j=n+1}^k(1-b_j)P·E{t}_{ij}}\right]g_{X^{\ast }}^{(k)}(x)dx}}$$

(12′)

where $E{t}_{\mathrm{in}}={\int }_0^{\infty }{\overline{F}}_{{\tau }_{{\xi }_n}^{(i)}-{\tau }_x^{(i)}}^{(i)}(t)\mathrm{dt}$ given that $x<X_{i-1}^{*}<x+\mathrm{dx}$ . Note that in the above calculations, we rewrite expression (12) as (12′). This is due to the fact that when ${\xi }_{n-1}<X_{i-1}^{*}\le {\xi }_n$ and $n\le j\le k$ , there is no uniform expression of $E{t}_{\mathrm{ij}}$ since $E{t}_{\mathrm{in}}={\int }_0^{\infty }{\overline{F}}_{{\tau }_{{\xi }_n}^{(i)}-{\tau }_x^{(i)}}^{(i)}(t)\mathrm{dt}$ while $E{t}_{\mathrm{ij}}={\int }_0^{\infty }{\overline{F}}_{{\tau }_{{\xi }_j}^{(i)}-{\tau }_{{\xi }_{j-1}}^{(i)}}^{(i)}(t)\mathrm{dt}(n+1\le j\le k)$ . Therefore, in expression (12′), we list $E{t}_{\mathrm{in}}$ separately.

Substituting equations (11) and (12′) into equation (10), we can obtain the expression of EQ as follows

$$\begin{array}{c}\mathrm{EQ}=\sum _{j=1}^k(1-b_j)P·E{t}_{1j}+\sum _{i=2}^N\sum _{n=1}^k\hfill \\ \underset{{\xi }_{n-1}}{\overset{{\xi }_n}{\int }}[(1-b_n)P·E{t}_{\mathrm{in}}+\sum _{j=n+1}^k(1-b_j)P·E{t}_{\mathrm{ij}}]g_{X^{*}}^{(k)}(x)\mathrm{dx}\hfill \end{array}$$

(13)

Expected total cost incurred in a renewal cycle

The total cost $C_{\mathrm{Total}}$ incurred in a renewal cycle includes the PM cost $C_{\mathrm{PM}}$ , replacement cost $C_R$ , production cost $C_P$ (e.g. raw materials cost and labor cost), inventory holding cost $C_{\mathrm{IH}}$ and penalty cost of lost sales $C_{\mathrm{LS}}$ .

The expected total PM cost incurred in a renewal cycle can be formulated as follows

$$E{C}_{\mathrm{PM}}=(N-1)c_{\mathrm{pm}}·E(T_{\mathrm{pm}})=(N-1)c_{\mathrm{pm}}\underset{0}{\overset{\infty }{\int }}t\phi (t)\mathrm{dt}$$

(14)

where $\phi (·)$ is the probability density function of the random variable $T_{\mathrm{pm}}$ .

The replacement cost incurred in a renewal cycle is as follows

$$E{C}_R=c_r$$

(15)

Denoting by $T_i(1\le i\le N)$ the length of the ith production cycle (i.e. ${\tau }_{{\xi }_k^{(i)}}^{(i)}$ ) and by ${\overline{F}}_{{\tau }_{{\xi }_k^{(i)}}}^{(i)}$ its distribution function, the expected total production cost incurred in a renewal cycle can be given by

$$E{C}_P=c_p·P·\sum _{i=1}^{N}E{T}_i=c_p·P·(E{T}_1+\sum _{i=2}^{N}E{T}_i)$$

(16)

where $E{T}_1={\int }_0^{\infty }{\overline{F}}_{{\tau }_{{\xi }_k}}(t)\mathrm{dt}$ and $E{T}_i={\int }_0^{{\xi }_k}{\int }_0^{\infty }{\overline{F}}_{{\tau }_{{\xi }_k}^{(i)}-{\tau }_x^{(i)}}^{(i)}(t)g_{X^{*}}^{(k)}(x)\mathrm{dtdx},2\le i\le N$ .

Let $S_i(1\le i\le N)$ denote the area under the curve of the inventory of conforming items in the ith production cycle (Figure 1); then the expected total inventory holding cost incurred in a renewal cycle can be modeled as follows

$$E{C}_{\mathrm{IH}}=\sum _{i=1}^N{c}_h·E{S}_i$$

(17)

According to Figure 1, the area $S_1$ can be divided into $(k+1)$ parts, that is, $S_{11},S_{12},\dots ,S_{1k},S_{1,k+1}$ . Therefore, we have the following

$$\begin{array}{cc}\multicolumn{1}{c}{E{S}_1}& =\sum _{j=1}^{k+1}E{S}_{1j}=\sum _{j=1}^{k}E{S}_{1j}+E{S}_{1,k+1}\hfill \\ \multicolumn{1}{c}{}& =\sum _{j=1}^{k}E\left\{\frac{1}{2}{t}_{1j}\{2\sum _{l=1}^{j-1}[(1-b_l)P-D]t_{1l}+[(1-b_j)P-D]t_{1j}\}\right\}+E\left\{\frac{{\{\sum _{l=1}^k[(1-b_l)P-D]t_{1l}\}}^2}{2D}\right\}\hfill \\ \multicolumn{1}{c}{}& =\sum _{j=1}^k\{\sum _{l=1}^{j-1}{r}_{l}E{t}_{1j}E{t}_{1l}+\frac{1}{2}{r}_{j}E(t_{1j}^2)\}+E\left\{\frac{{\left(\sum _{l=1}^k{r}_l{t}_{1l}\right)}^2}{2D}\right\}(\mathrm{here},\mathrm{we}\mathrm{set}{r}_l=[(1-b_l)P-D])\hfill \\ \multicolumn{1}{c}{}& =\sum _{j=1}^k\{\sum _{l=1}^{j-1}{r}_{l}E{t}_{1j}E{t}_{1l}+\frac{1}{2}{r}_{j}E(t_{1j}^2)\}+\left\{\left[\frac{\sum _{l=1}^k{r}_l^{2}E(t_{1l}^2)+2\sum _{l=1}^{k-1}\sum _{m=l+1}^k{r}_l{r}_{m}E{t}_{1l}E{t}_{1m}}{2D}\right]\right\}\hfill \end{array}$$

(18)

For $i\ge 2$ , the calculation of $E{S}_i$ is slightly different from the calculation of $E{S}_1$ since the initial degradation level of the production system in the ith production cycle is different from that of the first production cycle. More exactly, as mentioned above, $S_1$ consists of $(k+1)$ parts, while $S_i$ consists of $(n+1)$ parts given that the initial degradation level in the ith production cycle $X_{i-1}^{*}\in ({\xi }_{n-1},{\xi }_n],(1\le n\le k)$ . Therefore, conditioning on $X_{i-1}^{*}$ and applying the law of total probability, $E{S}_i(i\ge 2)$ can be formulated as follows

$$\begin{array}{cc}\multicolumn{1}{c}{E{S}_i}& =\sum _{n=1}^{k}E(S_i|{\xi }_{n-1}<X_{i-1}^{*}\le {\xi }_n)·P({\xi }_{n-1}<X_{i-1}^{*}\le {\xi }_n)=\sum _{n=1}^k\sum _{j=n}^{k+1}\underset{{\xi }_{n-1}}{\overset{{\xi }_n}{\int }}E{S}_{\mathrm{ij}}·g_{X^{*}}^{(k)}(x)\mathrm{dx}\hfill \\ \multicolumn{1}{c}{}& =\sum _{n=1}^k\{\sum _{j=n}^k\underset{{\xi }_{n-1}}{\overset{{\xi }_n}{\int }}E{S}_{\mathrm{ij}}·g_{X^{*}}^{(k)}(x)\mathrm{dx}+\underset{{\xi }_{n-1}}{\overset{{\xi }_n}{\int }}E{S}_{i,k+1}·g_{X^{*}}^{(k)}(x)\mathrm{dx}\}\hfill \\ \multicolumn{1}{c}{}& =\sum _{n=1}^k\{\sum _{j=n}^k\underset{{\xi }_{n-1}}{\overset{{\xi }_n}{\int }}E\left\{\frac{1}{2}{t}_{\mathrm{ij}}[2\sum _{l=n}^{j-1}{r}_l{t}_{\mathrm{il}}+r_j{t}_{\mathrm{ij}}]\right\}{g}_{X^{*}}^{(k)}(x)\mathrm{dx}+\underset{{\xi }_{n-1}}{\overset{{\xi }_n}{\int }}E\left\{\left(\frac{{\left(\sum _{l=n}^k{r}_l{t}_{\mathrm{il}}\right)}^2}{2D}\right)\right\}\{g_{X^{*}}^{(k)}(x)\mathrm{dx}\}\}\hfill \end{array}$$

(19)

$$\begin{array}{c}E{S}_i=\sum _{n=1}^k\{\underset{{\xi }_{n-1}}{\overset{{\xi }_n}{\int }}\{\frac{1}{2}{r}_{n}E(t_{\mathrm{in}}^2)+\sum _{j=n+1}^k[r_{n}E{t}_{\mathrm{ij}}E{t}_{\mathrm{in}}+\sum _{l=n+1}^{j-1}{r}_{l}E{t}_{\mathrm{ij}}E{t}_{\mathrm{il}}+\frac{1}{2}{r}_{j}E(t_{\mathrm{ij}}^2)]\}{g}_{X^{*}}^{(k)}(x)\mathrm{dx}+\hfill \\ \underset{{\xi }_{n-1}}{\overset{{\xi }_n}{\int }}\left\{\frac{[r_n^{2}E(t_{\mathrm{in}}^2)+\sum _{l=n+1}^k{r}_l^{2}E{(t_{\mathrm{il}})}^2+\sum _{m=n+1}^k{r}_n{r}_{m}E{t}_{\mathrm{in}}E{t}_{\mathrm{im}}+\sum _{l=n+1}^{k-1}\sum _{m=l+1}^k{r}_l{r}_m\mathrm{mE}{t}_{\mathrm{il}}E{t}_{\mathrm{im}}]}{2D}\right\}{g}_{X^{*}}^{(k)}(x)\mathrm{dx}\}\hfill \end{array}$$

(19′)

In the above calculations, we rewrite expression (19) as (19′) based on the similar consideration in the derivation of equation (12′). Substituting equations (18) and (19′) into equation (17), we can obtain the expression of $E{C}_{\mathrm{IH}}$ as follows

$$\begin{array}{c}E{C}_{IH}=c_h\cdot \{{\displaystyle \sum _{j=1}^k\left\{\left\{{\displaystyle \sum _{l=1}^{j-1}{r}_{l}E{t}_{1j}E}{t}_{1l}+\frac{1}{2}{r}_{j}E(t_{1j}^2)\right\}\right\}}+\left\{\frac{\left[{\displaystyle {\sum }_{l=1}^k{r}_l^2}E(t_{1l}^2)+2{\displaystyle {\sum }_{l=1}^{k-1}{\displaystyle {\sum }_{m=l+1}^k{r}_l{r}_{m}E{t}_{1l}E{t}_{1m}}}\right]}{2D}\right\}+\\ {\displaystyle \sum _{i=2}^N{\displaystyle \sum _{n=1}^k\left\{{\displaystyle \underset{{\xi }_{n-1}}{\overset{{\xi }_n}{\int }}\left\{\frac{1}{2}{r}_{n}E(t_{in}^2)+{\displaystyle \sum _{j=n+1}^k\left[r_{n}E{t}_{ij}E{t}_{in}+{\displaystyle \sum _{l=n+1}^{j-1}{r}_{l}E{t}_{ij}E{t}_{il}}+\frac{1}{2}{r}_{j}E(t_{ij}^2)\right]}\right\}{g}_{X^{*}}^{(k)}(x)dx}+{\displaystyle \underset{{\xi }_{n-1}}{\overset{{\xi }_n}{\int }}\frac{\left[r_n^{2}E(t_{in}^2)+{\displaystyle {\sum }_{l=n+1}^k{r}_l^{2}E{(t_{il})}^2}+{\displaystyle {\sum }_{m=n+1}^k{r}_n{r}_{m}E{t}_{in}E{t}_{im}}+{\displaystyle {\sum }_{l=n+1}^{k-1}{\displaystyle {\sum }_{m=l+1}^k{r}_l{r}_{m}mE{t}_{il}E{t}_{im}}}\right]}{2D}+}{g}_{X^{*}}^{(k)}(x)dx\right\}}}\end{array}$$

(20)

Next, we will derive the expression of the expected total penalty cost of lost sales $E{C}_{\mathrm{LS}}$ . Note that shortage occurs only when the duration of PM $T_{\mathrm{pm}}$ or the duration of replacement $T_r$ is greater than the time to deplete all on-hand inventory $T_{D_i}(1\le i\le N)$ . Thus, $E{C}_{\mathrm{LS}}$ can be modeled as follows

$$E{C}_{LS}=c_{l}D·\left\{{\displaystyle \sum _{i=1}^{N-1}E(T_{pm}-T_{D_i})·1_{\{T_{pm}>T_{D_i}\}}}+E(T_r-T_{D_N})·1_{\{T_r>T_{D_N}\}}\right\}$$

(21)

where $1_{\{A\}}$ is an indicator function, that is, $1_{\{A\}}=1$ if event A occurs, otherwise $1_{\{A\}}=0$ . To derive the explicit expression of $E{C}_{\mathrm{LS}}$ , we need to calculate the probability density function of random variables $T_{D_i}(1\le i\le N)$ . To this end, we first calculate the on-hand inventory of conforming items at the time just before the ith PM is performed. Let $I_i$ denote this inventory level, according to Figure 1, we have

$$I_1=\sum _{j=1}^k[(1-b_j)P-D]t_{1j}=\sum _{j=1}^k{r}_j{t}_{1j}$$

Consequently

$$T_{D_1}=\frac{I_1}{D}={\displaystyle \sum _{j=1}^k\frac{r_j}{D}{t}_{1j}}={\displaystyle \sum _{j=1}^k{u}_j{t}_{1j}}\left(\mathrm{here\; we\; set}{u}_j\mathrm{=}\frac{{\mathrm{r}}_{\mathrm{j}}}{\mathrm{D}}\right)$$

It can be seen that $T_{D_1}$ is a linear combination of random variables $t_{1j}(j=1,2,\dots ,k)$ . Let $h_1(t)$ be the density function of $T_{D_1}$ . Then, $h_1(t)$ is a k-fold convolution of the density functions of $u_j{t}_{1j}$ and we have

$$h_1(t)=\left({\displaystyle \prod _{j=1}^k{u}_j^{-1}}\right)·f_{{\tau }_{{\xi }_1}^{(1)}}^{(1)}\left(\frac{t}{u_1}\right)\ast f_{{\tau }_{{\xi }_2}^{(1)}-{\tau }_{{\xi }_1}^{(1)}}^{(1)}\left(\frac{t}{u_2}\right)\ast \cdots \ast f_{{\tau }_{{\xi }_k}^{(1)}-{\tau }_{{\xi }_{k-1}}^{(1)}}^{(1)}\left(\frac{t}{u_k}\right)$$

where $f_{{\tau }_{{\xi }_j}^{(1)}-{\tau }_{{\xi }_{j-1}}^{(1)}}^{(1)}$ can be approximated by equation (8).

When $i\ge 2$ , the calculation of the density function $h_i(t)$ of $T_{D_i}$ is slightly different from that of $T_{D_1}$ due to the non-zero initial degradation level before the ith production cycle. Conditioning on $X_{i-1}^{*}$ and applying the law of total probability, $h_i(t)$ can be obtained similarly which is given by

$$h_i(t);={\displaystyle \sum _{n=1}^k{\displaystyle \underset{{\xi }_{n-1}}{\overset{{\xi }_n}{\int }}\left({\displaystyle \prod _{j=n}^k{u}_j^{-1}}\right)}}·{\left[f_{{\tau }_{{\xi }_n}^{(i)}-{\tau }_x^{(i)}}^{(i)}\left(\frac{t}{u_n}\right)\ast f_{{\tau }_{{\xi }_{n+1}}^{(i)}-{\tau }_{{\xi }_n}^{(i)}}^{(i)}\left(\frac{t}{u_{n+1}}\right)\ast \cdots \ast f_{{\tau }_{{\xi }_k}^{(i)}-{\tau }_{{\xi }_{k-1}}^{(i)}}^{(i)}\left(\frac{t}{u_k}\right)\right]}^{g_{X^{\ast }}^{(k)}(x)dx}$$

Consequently, we can derive the expression of $E{C}_{\mathrm{LS}}$ as follows

$$\begin{array}{c}\multicolumn{1}{c}{E{C}_{\mathrm{LS}}=c_{l}D·\{\sum _{i=1}^{N-1}E(T_{\mathrm{pm}}-T_{D_i})·1_{\{T_{\mathrm{pm}}>T_{D_i}\}}+E(T_r-T_{D_N})·1_{\{T_r>T_{D_N}\}}\}}\\ \multicolumn{1}{c}{=c_{l}D·\{\sum _{i=1}^{N-1}\underset{0}{\overset{\infty }{\int }}\underset{x}{\overset{\infty }{\int }}(y-x)\phi (y)h_i(x)\mathrm{dydx}+\underset{0}{\overset{\infty }{\int }}\underset{x}{\overset{\infty }{\int }}(z-x)\phi (z)h_N(x)\mathrm{dydx}\}}\end{array}$$

(22)

where $\phi (·)$ and $\varphi (·)$ are the density functions of $T_{\mathrm{pm}}$ and $T_r$ , respectively.

Now, the objective function $C(k,N)$ can be obtained as follows

$$C(k,N)=\frac{E{C}_{\mathrm{Total}}}{EQ}=\frac{E{C}_{PM}+E{C}_R+E{C}_P+E{C}_{IH}+E{C}_{LS}}{EQ}$$

(23)

where $E{C}_{\mathrm{PM}}$ is given by equation (14), $E{C}_R$ is given by equation (15), $E{C}_P$ is given by equation (16), $E{C}_{\mathrm{IH}}$ is given by equation (20), $E{C}_{\mathrm{LS}}$ is given by equation (22) and EQ is given by equation (13).

The objective is to determine the optimal bivariate policy $(k^{*},N^{*})$ such that the average cost rate $C(k,N)$ is minimized. Due to the high complex nature of the above expressions, a closed-form solution cannot be derived. However, the optimization of $C(k,N)$ for a data set can be performed using numerical methods. When k is fixed, $C(k,N)$ is a function of N. When $1\le k\le M$ , we rewrite $C(k,N)$ as $C_k(N)$ . For a fixed k, we can find $N_k^{*}$ numerically such that $C_k(N_k^{*})$ is minimized. Consequently, we can find $N_1^{*},N_2^{*},\dots ,N_k^{*}$ , respectively, and the minimum average cost rate $C(k^{*},N^{*})$ is equal to $min\{C(1,N_1^{*}),C(2,N_2^{*}),\dots ,C(M,N_M^{*})\}$ .

An illustrative example and discussions

In this section, we provide some numerical results and sensitivity analysis to illustrate the model proposed. It is assumed that the distribution functions of the PM time $T_{\mathrm{pm}}$ and replacement time $T_r$ are exponential distributions $\phi (y)=1-e^{-\lambda y},y>0$ and $\varphi (z)=1-e^{-\mu z},z>0$ , respectively. The input data sets are given in Table 1.

Table 1.

Input data.

Parameter	Value	Parameter	Value	Parameter	Value
P	30 units/day	${\xi }_3$	3	q	1.2
D	26 units/day	${\xi }_4$	4	$\lambda$	0.5/day
$\alpha$	0.05	${\xi }_5$	5	$\mu$	0.25/day
$\beta$	1	$b_1$	3%	cpm	1000 US$/day
a	1.05	$b_2$	4%	cr	60,000 US$/time
M	4	$b_3$	5%	cp	50 US$/unit
${\xi }_1$	1	$b_4$	6%	ch	0.05 US$/unit/day
${\xi }_2$	2	p	0.4	cl	20 US$/unit

Using the computational software Mathematica 9.0, we can obtain the results presented in Table 2. For different PM threshold k (i.e. $k=1,2,3,4$ ), the optimal replacement policy $N_k^{*}$ , the expected total quantity of conforming items produced Q, the related costs and finally the average cost rate $C(k,N_k^{*})$ are reported in Table 2. The average cost rate as function of N for different k is illustrated in Figure 2. From Table 2 and Figure 2, we observe that when $k=4$ the optimal replacement policy $N_k^{*}=28$ with a minimum average cost rate C(4, 28) = 62.787 US$/unit. The solution process is implemented on DELL XPS with 2.2 GHz Core i7 processor, 8 GB of memory and Windows 7 operating system. We calculate 4 × 40 = 160 values of the objective function and it consumed 2012.8 s.

Table 2.

Computational results for different PM thresholds k.

k	$N_k^{*}$	Q (units)	CPM (US$/unit)	CR (US$/unit)	CP (US$/unit)	CIH (US$/unit)	CLS (US$/unit)	C(k, N*) (US$/unit)
1	26	10,407.2	4.804	5.765	52.694	0.187	0.030	63.480
2	25	15,543.4	3.088	3.860	56.325	0.337	0.022	63.632
3	27	21,381.1	2.432	2.806	57.102	0.549	0.021	62.910
4	28	26,772.8	2.017	2.241	57.785	0.716	0.018	62.787

Figure 2.

The average cost rate as function of N for different k.

It can be seen from Table 2 that the expected total quantity of conforming items produced Q increases with the increasing in PM threshold k. The PM, replacement and lost sales costs per conforming unit produced decrease as PM threshold increases. These can be explained by the fact that as PM threshold is increased, the production system works longer before PM is performed and more conforming items are produced in a renewal cycle which results in lower PM, replacement and lost sales costs per conforming unit produced. However, both the production and inventory holding costs per conforming unit produced increase as PM threshold increases. The former is due to the fact that as PM threshold is increased, the condition of the system becomes worse and worse. According to assumption A3, that is, the system produces non-conforming items at a higher rate as the degradation level increases; hence, more non-conforming items are produced and the proportion of defectives is increased which leads to a higher production cost per conforming unit produced. It can be interpreted that as PM threshold is increased, more items are produced and put into inventory. Therefore, the inventory holding cost per conforming unit produced is increased.

The sensitivity analysis is performed to investigate the influence of some model parameters on optimal bivariate policy as well as the expected cost rate per unit product. To this end, we vary one of the parameters at a time and keep other parameters unchanged. The results are reported in Table 3. From Table 3, we have the following observations:

• As $c_{\mathrm{pm}}$ increases, the optimal PM threshold $k^{*}$ increases, whereas the number of PM actions $N^{*}$ before replacement decreases. These phenomena can be illustrated as follows. When the PM cost per unit time increases, the optimal PM threshold should be increased to postpone the PM performance which can lower the total PM costs. Again, with the increase in $c_{\mathrm{pm}}$ , the optimal replacement policy $N^{*}$ should be reduced to decrease the number of PM actions before replacement which can also lower the total PM costs. In short, the former aims to reduce the frequency of PM actions and the latter aims to reduce the number of PM actions.

• Due to higher replacement cost $c_r$ , the optimal PM threshold and number of PM actions before replacement are greater, resulting in longer working time of the production system before its replacement. As longer time the production system works, more non-conforming items are produced and more items are put into inventory. Thus, the production and inventory holding cost per conforming unit produced are greater.

• When the production cost per unit $c_p$ is increased, the optimal PM threshold and number of PM actions before replacement decrease. This is due to the fact that as the production cost per unit increases, we need to decrease the proportion of defectives to lower the total production cost. To this end, the PM threshold should be reduced and the aging system should be replaced earlier to maintain the production system in good condition. However, the PM and replacement costs per unit produced increase as shown in Table 3.

• A higher inventory holding cost per unit item per unit time $c_h$ decreases the optimal PM threshold (or production run time) and the number of PM actions before replacement, resulting in lower production cost to compensate for the higher PM, replacement, inventory holding and lost sales costs.

• Finally, the optimal PM threshold and replacement policy are insensitive to the lost sales cost per unit item per unit time $c_l$ . This may be due to the fact that the possibility of the shortage occurrence is quite low in our model.

Table 3.

Sensitivity analysis.

Parameter values	k*	N*	Q (units)	CPM (US$/unit)	CR (US$/unit)	CP (US$/unit)	CIH (US$/unit)	CLS(US$/unit)	C(k, N*)(US$/unit)
cpm
500	1	34	11630.3	2.837	5.159	52.692	0.217	0.029	60.934
750	4	31	27933.9	1.611	2.148	57.768	0.715	0.018	62.259
1250	4	25	25418.5	2.360	2.360	57.806	0.718	0.019	63.264
1500	4	23	24392.7	2.706	2.459	57.824	0.719	0.019	63.728
cr
20,000	1	16	7971.9	3.763	2.509	52.701	0.218	0.034	59.225
40,000	1	22	9579.6	4.384	4.176	53.696	0.217	0.031	62.505
80,000	4	30	27566.6	2.104	2.902	57.773	0.715	0.018	63.513
100,000	4	33	28614.5	2.237	3.495	58.759	0.714	0.018	65.224
cp
30	4	27	26344.3	1.974	2.278	34.675	0.717	0.019	39.661
40	4	27	26344.3	1.974	2.278	46.233	0.717	0.019	51.219
60	1	26	10407.2	4.804	5.765	63.233	0.217	0.030	74.049
70	1	26	10407.2	4.804	5.765	73.772	0.217	0.030	84.588
ch
0.02	4	28	26772.8	2.017	2.241	57.785	0.286	0.018	62.348
0.035	4	28	26772.8	2.017	2.241	57.785	0.501	0.018	62.563
0.1	3	27	21381.1	2.432	2.806	57.102	1.098	0.021	63.459
0.15	1	26	10407.2	4.804	5.765	52.694	0.651	0.030	63.945
cl
10	4	28	26772.8	2.017	2.241	57.785	0.286	0.016	62.775
15	4	28	26772.8	2.017	2.241	57.785	0.286	0.017	62.776
25	4	28	26772.8	2.017	2.241	57.785	0.286	0.019	62.778
30	4	28	26772.8	2.017	2.241	57.785	0.286	0.020	62.779

To evaluate the cost saving potential of the proposed model, we compare our results with the results obtained from the ABM model. In the ABM model, the PM is performed on production system after every T production time. Since the PM is imperfect, the system will be replaced by a new identical one after N PMs or at the system’s failure whichever occurs first. T and N are two decision variables of the ABM model. The detailed descriptions and model formulations are given in Appendix 1.

For the policy $(T,N)$ , the input data are the same as given in Table 1. By Monte Carlo simulation implemented on MATLAB R2014a, the minimum average cost rate per unit C(T*, N*) = 69.164 US$/unit is obtained with the optimal PM interval $T^{*}=43\mathrm{days}$ and the optimal replacement policy $N^{*}=26$ . That is, the production system is preventively maintained after every 43 days and it will be replaced after 26 production cycles or at the system failure whichever occurs first. We can see that the average cost rate of ABM model is 69.164 US$/unit which is greater than that of CBM model 62.787 US$/unit.

For comparing the performance of the two policies: CBM policy and ABM policy, we compute the relative cost savings as follows

$$\Delta =\frac{C(T^{\ast },N^{\ast })-C(k^{\ast },N^{\ast })}{C(T^{\ast },N^{\ast })}\times 100\%$$

More comparable results of the two models are given in Table 4. For the 25 instances in numerical experiment, the average cost saving of CBM model is 9.4%. As expected, the CBM is more efficient than ABM due to the fact that the PM is performed based on the real-time state estimation rather than fixed time interval. The Another reason of cost reduction of CBM model lies in that the product quality is dependent on the system’s deterioration level, the CBM can control the product quality well (i.e. the PM is performed once the ratio of non-conforming items produced is higher than the warning line). In addition, as we can find in Table 4, the frequency of system replacement ( $N^{*}$ ) of ABM model is higher than that of CBM, which means that the CBM can prolong the life time of the production system and hence, it is more effective than ABM.

Table 4.

Comparable results obtained for two policies.

Parameter values	k*	N*	C(k*, N*)(US$/unit)	T* (day)	N*	C(T*, N*)(US$/unit)	$\Delta$ (%)
cpm
500	1	34	60.934	32	31	66.368	8.2
750	4	31	62.259	38	28	67.972	8.4
1000	4	28	62.787	43	26	69.164	9.2
1250	4	25	63.264	47	23	70.757	10.6
1500	4	23	63.728	52	21	71.812	11.3
cr
20,000	1	16	59.225	49	15	66.579	11.0
40,000	1	22	62.505	46	21	68.357	8.6
60,000	4	28	62.787	43	26	69.164	9.2
80,000	4	30	63.513	39	29	70.203	9.5
100,000	4	33	65.224	35	31	71.416	8.7
cp
30	4	27	39.661	46	25	45.247	12.3
40	4	27	51.219	45	25	57.325	10.7
50	4	28	62.787	43	26	69.164	9.2
60	1	26	74.049	41	26	81.425	9.1
70	1	26	84.588	40	26	92.793	8.8
ch
0.02	4	28	62.348	48	27	67.937	8.2
0.035	4	28	62.563	46	26	68.475	8.6
0.05	4	28	62.787	43	26	69.164	9.2
0.1	3	27	63.459	41	25	70.246	9.7
0.15	1	26	63.945	38	25	70.984	9.9
cl
10	4	28	62.775	43	26	69.105	9.2
15	4	28	62.776	43	26	69.143	9.2
20	4	28	62.787	43	26	69.164	9.2
25	4	28	62.778	43	26	69.186	9.3
30	4	28	62.779	43	26	69.192	9.3

Conclusion

This article proposes an integrated strategy for the joint optimization of CBM and replacement policy of an imperfect production-inventory system. The production system deteriorates gradually and its degradation process is modeled by a non-stationary gamma process. For more realistic, we assume in this article a degradation-level-dependent product quality which implies that the proportion of non-conforming items being produced increases as the system deteriorates worse. To decrease the quality loss, PM is performed when the degradation level reaches the threshold, but it cannot restore the system as good as new. Hence, the aging system will be replaced by an identical new one after the Nth PM. We derived the explicit expression of the total cost per unit item $C(k,N)$ and obtained the optimal joint policy $(k^{*},N^{*})$ numerically such that $C(k,N)$ is minimized. A comparable model, i.e., ABM model, is given to evaluate cost saving potential of our model. The results of experiment show that the CBM model is more efficient and effective than the ABM model.

As a future work, we can relax some assumptions. A straightforward relaxation consists to consider recoverable items and rework process. This relaxation allows to deal with the scheduling problem of rework process, that is, whether the non-conforming items are reworked in the same cycle or after M cycles. Other future works can consider the random demand rate in the same context which makes the model be more realistic.