A QMM-MOP methodology for the maintenance scheduling of multistage manufacturing systems with a stream of deterioration

Abstract

In multistage manufacturing systems, a stream-of-deterioration (SOD) phenomenon poses two challenges for effective preventive-maintenance (PM) scheduling. First, the deterioration of each machine contributes to the deterioration of the final-product quality, and thus timely PM should be conducted to prevent excessive quality deterioration. Second, the deterioration of different machines leads to different degrees of deterioration in the final-product quality; thus, the PM of different machines will result in different degrees of improvement in the final-product quality. To address both challenges, a QMM-MOP methodology that adopts an interactive bi-level scheduling framework is proposed. In machine-level scheduling, a quality-integrated maintenance model (QMM) is developed by incorporating intermediate-product quality deterioration into the total cost to schedule timely PM for each individual machine. In system-level scheduling, maintenance-operation prioritization (MOP), based on a SOD-enabled quality-improvement factor, is proposed to select machines for PM. The case study shows that the proposed methodology can ensure a higher final-product quality with a lower total cost. The contribution of this paper is to develop a QMM-MOP methodology that integrates product-quality improvement into an interactive bi-level PM scheduling framework and enables MOP based on the-quality improvement factor to best improve the final-product quality.

Introduction

With the development of advanced manufacturing technologies, multistage manufacturing systems (MMSs), which produce products through sequential production stages, have been widely used in various industrial applications. ¹ Product quality and machine reliability are key performance indicators for MMSs and are significantly affected by the conditions of the machines. Preventive maintenance (PM) can improve the conditions of the machines and consequently improve the product quality and machine reliability.^2–5 However, PM activities consume resources and interrupt normal production, which may lead to considerable costs. Therefore, for MMSs, effective PM scheduling is required to ensure high product quality and machine reliability with the lowest possible PM cost.

The PM problem of MMSs belongs to the PM problems of multi-unit manufacturing systems, and many PM policies and models have been proposed for this type of problems.^6–13 However, existing PM policies and models for multi-unit manufacturing systems mostly focus on improving unit reliability and reducing maintenance cost, while they seldom consider the improvement of product quality. For MMSs, product quality is as important as machine reliability; thus, product-quality improvement should be considered in the PM scheduling of MMSs.

To develop PM policies for MMSs with considering the improvement of both product quality and machine reliability, a stream-of-deterioration (SOD) phenomenon^14,15 in the systems should be considered. The SOD phenomenon means that the deterioration of quality-related components (QRCs) of a machine will not only cause deviations in key product characteristics (KPCs) formed in the current stage, but also propagate to downstream stages and cause further deviations in the KPCs formed there. QRCs refer to machine components that have a significant impact on product quality, such as the cutting tools of a machine or the dies of a punching press. The SOD phenomenon poses two challenges for the PM scheduling of MMSs.

First, because of the SOD, the deviations of KPCs of the final products are the result of the combined effects of the QRC deterioration of machines in all stages. As each stage can introduce KPC deviations owing to the QRC deterioration, the final-product quality is subject to a high risk of deterioration. Therefore, timely PM should be conducted to decrease the product quality deterioration caused by the SOD.

To deal with this challenge, some researchers have proposed PM models where PM is performed to decrease the deterioration of product quality and machine reliability.^14–23 In these studies, mathematical models were developed to quantitatively describe the SOD phenomenon and further predict product-quality deterioration. The quality deterioration of intermediate and final products is converted into quality losses, which are further incorporated into the total cost of machine-level and system-level scheduling to drive PM for machines. The results of these studies show that by incorporating product quality loss into the total cost, more PM operations will be scheduled for machines to reduce the product-quality deterioration. More PM operations will increase PM costs, but will simultaneously bring about a greater reduction in product-quality loss.

Second, because of the SOD phenomenon, the deterioration of different machines will lead to different degrees of deterioration in the final-product quality. Therefore, the PM of different machines will result in different degrees of improvement in the final-product quality. This fact is ignored in existing quality-integrated PM models.^14–23 According to this fact, it is desirable to prioritize machines with greater abilities to improve the final-product quality through PM in the system-level scheduling. In addition, a quality-loss threshold (QLT) should be imposed on the final-product quality loss to trigger PM in the system-level scheduling. However, this challenge has not been addressed in the existing literature.

To deal with both challenges in the PM scheduling of MMSs, a QMM-MOP methodology adopting an interactive bi-level scheduling framework is proposed. In machine-level scheduling, a quality-integrated maintenance model (QMM) is developed to obtain the optimal PM interval for each machine, cycle by cycle, where the intermediate-product quality deterioration is converted into quality loss and then incorporated into the total cost to schedule timely PM. The optimal PM intervals of individual machines are fed into the system-level scheduling. In system-level scheduling, PM is triggered when the final-product quality loss reaches a predetermined QLT, or the scheduled PM time of a machine arrives, whichever comes first. If PM is triggered by the QLT, maintenance-operation prioritization (MOP) will be conducted, based on a quality improvement factor, to select machines for PM. This factor is developed, based on the SOD model, to evaluate the ability of each machine to improve final-product quality through PM. The decision results in system-level scheduling are fed back to the machine level for the next cycle scheduling.

The key differences between the proposed QMM-MOP methodology and existing methodologies cited previously are summarized as follows.

1. Conventional PM models^6–12 do not consider product-quality improvement in PM scheduling, while the QMM-MOP methodology considers product-quality improvement in PM scheduling at the machine and system levels.

2. Existing quality-integrated PM models^14–23 incorporate product-quality deterioration into the total cost to drive PM; however, they ignore the fact that the PM of different machines will result in different degrees of improvement in the final-product quality. The fact is considered in the QMM-MOP methodology. To this end, a quality-improvement factor is developed to evaluate the ability of each machine to improve the final-product quality through PM. This factor is further used to prioritize the PM operations of all machines in system-level scheduling. In the case-study section, the proposed methodology is compared with the above two types of models to demonstrate its superiority.

The reminder of this paper is organized as follows. Section 2 describes the SOD model. Section 3 presents the QMM-MOP methodology for the PM scheduling of MMSs, which includes the QMM for machine-level scheduling, MOP for system-level scheduling, and QMM-MOP for interactive bi-level scheduling. In Section 4, a case study is investigated using the proposed methodology, and results are discussed. In Section 5, conclusions are drawn and future work is presented.

Stream of deterioration (SOD) model

Consider a serial MMS applied to produce a kind of products during a finite planning horizon. In the serial MMS, the machines are subject to deterioration because of age and usage. Specifically, the deterioration of the QRCs of the machines will systematically lead to product quality deterioration and cause a SOD phenomenon. Lu and Zhou ¹⁴ developed a SOD model to quantitatively describe the SOD phenomenon. Under this model, the deviations of KPCs generated in each stage are considered as responses of three input factors: (i) the deviations of KPCs generated in upstream stages, (ii) the deterioration of QRCs of the machine in the current stage and (iii) the noise variables in the current stage. The noise variables refer to the factors that have significant impacts on product quality and vary randomly and cannot be controlled by maintenance. Examples of noise variables include the human factors, raw materials, environmental variations, etc. The noise variables are assumed to follow a multivariate normal distribution with zero means.

Figure 1 illustrates the SOD phenomenon in the serial MMS. In stage k, ${\mathbf{X}}_k(t)=[X_{k1}(t),\dots ,X_{k{n}_k}(t){]}^T$ denotes the deterioration states of the QRCs of machine k, where n_k denotes the number of variables measuring the deterioration states of QRCs of machine k; Z _k denotes the noise variables; and ${\mathbf{Y}}_k(t)=[Y_{1k}(t),\dots ,Y_{\mathrm{mk}}(t){]}^T$ denotes the deviations of the KPCs of the products output from stage k, where $Y_{\mathrm{sk}}\left(t\right)$ (s = 1, …, m) denotes the deviation of the sth KPC, and m denotes the number of KPCs of a final product. The noise variables follow a multivariate normal distribution with zero means, which can be represented as ${\mathbf{Z}}_k~N(0,{\mathrm{\Sigma }}_k^z)$ . From Figure 1, it can be seen that modeling the SOD is to mathematically describe the update of deviations of KPCs when the products flow along the stages. There are three situations to consider for modeling the update of ${\mathbf{Y}}_k(t)$ .

• Situation 1: The sth KPC is not formed in stage k and has not been formed in an upstream stage. In this situation, we have

$$Y_{\mathrm{sk}}\left(t\right)=0,k=1,\dots ,N$$

(1)

• Situation 2: The sth KPC is not formed in stage k but has been formed in an upstream stage. In this situation, we have

$$Y_{\mathrm{sk}}\left(t\right)={Y_s}_{(k-1)}\left(t\right),k=2,\dots ,N$$

(2)

• Situation 3: The sth KPC is formed in stage k, that is, $s\in O_k$ . In this situation, $Y_{\mathrm{sk}}\left(t\right)$ is the result of the combined effects of Y_k–1(t), X _k (t), and Z _k . The SOD model developed by Lu and Zhou ¹⁴ is adopted to mathematically describe this impact mechanism. The SOD model is given as

$$\begin{array}{c}\multicolumn{1}{c}{Y_{\mathrm{sk}}\left(t\right)={\phi }_{\mathrm{sk}}+\left[{{\mathbf{b}}_{\mathrm{sk}}}^T{{\mathbf{c}}_{\mathrm{sk}}}^T{{\mathbf{d}}_{\mathrm{sk}}}^T\right]\left[\begin{array}{c}{\mathbf{X}}_k\left(t\right)\hfill \\ \multicolumn{1}{c}{{\mathbf{Y}}_{k-1}\left(t\right)}\\ \multicolumn{1}{c}{{\mathbf{Z}}_k}\end{array}\right]}\\ \multicolumn{1}{c}{+\left[{{\mathbf{X}}_k}^T\left(t\right){{\mathbf{Y}}_{k-1}}^T\left(t\right){{\mathbf{Z}}_k}^T\right]\left[\begin{array}{ccc}0\hfill & {\mathbf{B}}_{\mathrm{sk}}\hfill & {\mathbf{C}}_{\mathrm{sk}}\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & {\mathbf{D}}_{\mathrm{sk}}\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill \end{array}\right]}\\ \multicolumn{1}{c}{\left[\begin{array}{c}{\mathbf{X}}_k\left(t\right)\hfill \\ \multicolumn{1}{c}{{\mathbf{Y}}_{k-1}\left(t\right)}\\ \multicolumn{1}{c}{{\mathbf{Z}}_k}\end{array}\right]+{\epsilon }_{\mathrm{sk}}}\\ s\in O_k;k=1,\dots ,N\hfill \end{array}$$

(3)

Figure 1.

Schematic illustration of SOD phenomenon in the serial MMS.

where ${\phi }_{\mathrm{sk}}$ is a constant, ${\mathbf{b}}_{\mathrm{sk}}$ , ${\mathbf{c}}_{\mathrm{sk}}$ , and ${\mathbf{d}}_{\mathrm{sk}}$ are vectors that define the linear effects of ${\mathbf{X}}_k\left(t\right)$ , ${\mathbf{Y}}_{k-1}\left(t\right)$ , and ${\mathbf{Z}}_k$ , respectively, ${\mathbf{B}}_{\mathrm{sk}}$ , ${\mathbf{C}}_{\mathrm{sk}}$ , and ${\mathbf{D}}_{\mathrm{sk}}$ are matrices that define the interactive effects between two of ${\mathbf{X}}_k\left(t\right)$ , ${\mathbf{Y}}_{k-1}\left(t\right)$ , and ${\mathbf{Z}}_k$ , and ${\epsilon }_{\mathrm{sk}}$ is the error term. ${\epsilon }_{\mathrm{sk}}$ is assumed to follow a normal distribution with zero mean, which is represented as ${\epsilon }_{\mathrm{sk}}~N(0,{{\sigma }_{\mathrm{sk}}}^2)$ . The parameters of the SOD model can be estimated based on data of products collected from the production process. The definition of data of products can be found in Section 3.2.1 of Lu and Zhou. ¹⁴

Based on the SOD model, a recursion formula that describes the transformation from ${{\mathbf{Y}}_k}_{-1}(t)$ to ${\mathbf{Y}}_k(t)$ can be derived as follows:

$${\mathbf{Y}}_k\left(t\right)={\mathbf{\Gamma }}_k\left(t\right){\mathbf{Y}}_{k-1}\left(t\right)+{\mathbf{L}}_k\left(t\right),k=1,\dots ,N$$

(4)

where ${\mathbf{\Gamma }}_k(t)=[{\mathbf{\Gamma }}_{1k}(t),\dots ,{\mathbf{\Gamma }}_{\mathrm{mk}}(t){]}^T$ with ${\mathbf{\Gamma }}_{\mathrm{sk}}(t)={{\mathbf{c}}_{\mathrm{sk}}}^T+{{\mathbf{X}}_k}^T(t){\mathbf{B}}_{\mathrm{sk}}+{{\mathbf{Z}}_k}^T{{\mathbf{D}}_{\mathrm{sk}}}^T$ , and ${\mathbf{L}}_k(t)=[L_{1k}(t),\dots ,L_{\mathrm{mk}}(t){]}^T$ with $L_{\mathrm{sk}}(t)={\phi }_{\mathrm{sk}}+{{\mathbf{b}}_{\mathrm{sk}}}^T{\mathbf{X}}_k(t)+{{\mathbf{d}}_{\mathrm{sk}}}^T{\mathbf{Z}}_k+{{\mathbf{X}}_k}^T(t){\mathbf{C}}_{\mathrm{sk}}{\mathbf{Z}}_k+{\epsilon }_{\mathrm{sk}}$ (s = 1, …, m). Y₀ represents the deviations of workblanks. It is assumed that ${\mathbf{Y}}_0\equiv 0$ for simplicity. Based on the recursion formula, it can be derived that

$$\begin{array}{c}\multicolumn{1}{c}{{\mathbf{Y}}_k(t)=\underset{i=1}{\overset{k}{\Pi }}{\mathbf{\Gamma }}_i(t){\mathbf{Y}}_0+\sum _{i=1}^{k-1}[\underset{l=i+1}{\overset{k}{\Pi }}{\mathbf{\Gamma }}_l(t)]{\mathbf{L}}_i(t)}\\ \multicolumn{1}{c}{+{\mathbf{L}}_k(t),k=1,\dots ,N}\end{array}$$

(5)

Based on equation (5), ${\mathbf{Y}}_k[t|X_1(t),...,X_k(t)]$ is the result of the combined effects of deterioration of the QRCs of machines in the previous k stages. ${\mathbf{Y}}_k[t|X_1(t),...,X_k(t)]$ will be converted into product-quality loss for PM optimization in Section 3.

Modeling the deterioration processes of QRCs

The deterioration processes of QRCs (i.e. the wear of cutting tool) are generally stochastic processes with independent and non-negative increments, which can be properly described by gamma process.^19,24 Assume that the deterioration processes of the QRCs of a machine are statistically independent. The deterioration states of the QRCs of machine k at t (t > 0) are represented as $X_{\mathrm{ki}}(t)(i=1,...,n_k)$ . As $\{X_{\mathrm{ki}}(t),t>0\}$ can be described by the gamma process, $X_{\mathrm{ki}}\left(t\right)$ follows a gamma distribution given by

$$\begin{array}{c}\multicolumn{1}{c}{f_{X_{\mathrm{ki}}\left(t\right)}\left[x_{\mathrm{ki}}|v_{\mathrm{ki}}\left(t\right),{\theta }_{\mathrm{ki}}\right]=\frac{{{\theta }_{\mathrm{ki}}}^{v_{\mathrm{ki}}\left(t\right)}{x_{\mathrm{ki}}}^{v_{\mathrm{ki}}\left(t\right)-1}\exp [-{\theta }_{\mathrm{ki}}{x}_{\mathrm{ki}}]}{\Gamma \left[v_{\mathrm{ki}}\left(t\right)\right]},}\\ \multicolumn{1}{c}{i=1,\dots ,n_k}\end{array}$$

(6)

where $v_{\mathrm{ki}}\left(t\right)$ and ${\theta }_{\mathrm{ki}}({\theta }_{\mathrm{ki}}>0)$ are the shape and scale parameters of the gamma distribution. The deterioration states of the QRCs are assumed to be 0 when the QRCs are in completely new states. That is $X_{\mathrm{ki}}\left(0\right)\equiv 0$ . Let $v_{\mathrm{ki}}\left(t\right)\equiv {\mu }_{\mathrm{ki}}t({\mu }_{\mathrm{ki}}>0)$ for $t>0$ . Then, the s-expectation and variance of $X_{\mathrm{ki}}\left(t\right)$ are given as

$$E\left[X_{\mathrm{ki}}\left(t\right)\right]=u_{\mathrm{ki}}t/{\theta }_{\mathrm{ki}},\mathrm{Var}\left[X_{\mathrm{ki}}\left(t\right)\right]={u_{\mathrm{ki}}t/{\theta }_{\mathrm{ki}}}^2$$

(7)

QMM-MOP methodology for the PM scheduling of a MMS

As stated previously, the deterioration of the QRCs of machines will systematically lead to product-quality deterioration. Thus, PM will be performed on machines to decrease the deterioration of product quality and machine reliability. In addition, whenever a machine fails, a minimal repair will be performed to restore it to the normal operating state. Frequent PM can ensure high product quality and machine reliability; however, it may cause excessively high PM costs. Thus, the QMM-MOP methodology is proposed to achieve the best trade-off between the PM cost and the cost caused by the deterioration of product quality and machine reliability.

The QMM-MOP methodology consists of interactive bi-level scheduling. In the machine-level scheduling, a QMM is developed to obtain the optimal PM intervals for each individual machine. In the QMM, the quality deterioration of intermediate products is converted into quality loss based on Taguchi’s loss function, and the quality loss, along with the minimal repair cost and PM cost, are integrated to form a total cost for PM time optimization.

In system-level scheduling, MOP is conducted based on the quality improvement factor to select machines for PM. The quality improvement factor evaluates the ability of each machine to improve the final-product quality through PM. Thus, the MOP enables the selection of the “best” machines for PM from the perspective of quality improvement. Note that the quality deterioration of the intermediate and final products can be predicted, based on the SOD model.

The following assumptions are made for the QMM-MOP methodology.

1. In a PM activity for a machine, all QRCs of the machine are repaired or replaced. Thus, the PM will restore the QRCs of the machine to completely new states. Meanwhile, those components of the machine that are judged to be in bad conditions are repaired or replaced. As a result, the PM will restore the machine condition to somewhere between “as good as new” and “as bad as old.”

2. Minimal repair does not change the machine condition or the deterioration states of QRCs. The duration of minimal repair is neglected.

3. The serial MMS is balanced. That is, the production rate of each machine is the same and equals the production rate of the system.

The QMM for machine-level scheduling

In this section, a QMM is developed to obtain the optimal PM interval for each machine, cycle by cycle, where intermediate-product quality deterioration is converted into quality loss and then incorporated into the total cost to schedule timely PM.

Evaluation of product quality loss

When developing the QMM for a machine, only the product-quality loss caused by the deterioration of the QRCs of the machine is considered. Hence, it is necessary to derive the deviations of the KPCs caused by the deterioration of the QRCs of each machine. To this end, the deviations of KPCs generated in the upstream stages are assumed to be 0. Then, the deviations of KPCs caused by the deterioration of the QRCs of machine k can be represented as $Y_{\mathrm{sk}}\left[t|{\mathbf{Y}}_{k-1}\left(t\right)=0\right]$ ( $s\in O_k$ ). $Y_{\mathrm{sk}}\left[t|{\mathbf{Y}}_{k-1}\left(t\right)=0\right]$ can be obtained by substituting ${\mathbf{Y}}_{k-1}\left(t\right)=0$ into equation (3):

$$\begin{array}{c}\multicolumn{1}{c}{Y_{\mathrm{sk}}\left[t|{\mathbf{Y}}_{k-1}\left(t\right)=0\right]={\phi }_{\mathrm{sk}}+{{\mathbf{b}}_{\mathrm{sk}}}^T{\mathbf{X}}_k\left(t\right)+{{\mathbf{d}}_{\mathrm{sk}}}^T{\mathbf{Z}}_k}\\ \multicolumn{1}{c}{+{{\mathbf{X}}_k}^T\left(t\right){\mathbf{C}}_{\mathrm{sk}}{\mathbf{Z}}_k+{\epsilon }_{\mathrm{sk}},s\in O_k}\end{array}$$

(8)

The product quality loss caused by the deviations of KPCs can be evaluated based on Taguchi’s loss function ²⁵ , which has the form $L(Y)=q(Y-\tau {)}^2$ . In the equation, $Y-\tau$ represents the deviation of a KPC and the constant parameter $q(q>0)$ is determined by financial considerations. Based on Taguchi’s loss function, the expected quality loss induced by $Y_{\mathrm{sk}}\left[t|{\mathbf{Y}}_{k-1}\left(t\right)=0\right]$ can be calculated as follows:

$$Q_k\left(t\right)=\sum _{s\in O_k}E\left\{q_s{Y_{\mathrm{sk}}}^2\left[t|{\mathbf{Y}}_{k-1}\left(t\right)=0\right]\right\}$$

(9)

where $q_s$ denotes the quality-loss coefficient of the sth KPC. By substituting equation (8) into equation (9), it can be derived that

$$\begin{array}{c}\multicolumn{1}{c}{Q_k\left(t\right)=\sum _{s=1}^m{q}_s}\\ \multicolumn{1}{c}{[{{\mathbf{v}}_k}^T{\mathrm{\Phi }}_{\mathrm{sk}}{\mathbf{v}}_k{t}^2+({{\mathbf{\gamma }}_{\mathrm{sk}}}^T{\mathbf{v}}_k+\sum _{i=1}^n{\varphi }_{\mathrm{ii}}^s{u}_{\mathrm{ki}}/{{\theta }_{\mathrm{ki}}}^2)t+{\omega }_{\mathrm{sk}}]}\end{array}$$

(10)

where ${\Phi }_{sk}={\left({\varphi }_{ij}^{sk}\right)}_{n_k\times n_k}\equiv C_{sk}{\Sigma }_k^z{C}_{sk}^T+b_{sk}{b}_{sk}^T$ , ${\mathbf{\gamma }}_{\mathrm{sk}}\equiv 2({\phi }_{\mathrm{sk}}{{\mathbf{b}}_{\mathrm{sk}}}^T+{{\mathbf{b}}_{\mathrm{sk}}}^T{\mathrm{\Sigma }}_k^z{{\mathbf{C}}_{\mathrm{sk}}}^T)$ , ${\omega }_{\mathrm{sk}}\equiv {{\sigma }_{\mathrm{sk}}}^2+{{\phi }_{\mathrm{sk}}}^2+{{\mathbf{d}}_{\mathrm{sk}}}^T{\mathrm{\Sigma }}_k^z{\mathbf{d}}_{\mathrm{sk}}$ , and ${\mathbf{v}}_k={[u_{k1}/u_{k1}{\theta }_{k1},...,u_{k{n}_k}/u_{k{n}_k}{\theta }_{k{n}_k}{\theta }_{k{n}_k}{\theta }_{k1},...,u_{k{n}_k}/u_{k{n}_k}{\theta }_{k{n}_k}{\theta }_{k{n}_k}]}^T$ . For the derivation method, please refer to Section 3.4.2 of Lu et al. ¹⁹

The QRCs of a machine are restored to completely new states after each PM. As a result, the quality-loss function will have the same form within each PM cycle. Note that a PM cycle is defined as the interval between consecutive PM operations.

Modeling of the machine-hazard rate

The deterioration of QRCs can directly influence their failure risk and can also affect the degradation processes of other components, owing to interactions between components. As a result, the deterioration states of QRCs may affect the failure risk of the entire machine. To evaluate the machine failure rate more accurately, it is desirable to integrate the deterioration states of the QRCs into the machine-hazard function and estimate their effects. The proportional hazard model (PHM) is commonly used to estimate the effects of covariates on the failure risk of a unit (e.g. a machine), where the covariates refer to the factors that may influence the failure risk of the unit. ²⁶ Clearly, the deterioration states of QRCs can be considered as covariates and further integrated into the machine-hazard function.

Given ${\mathbf{X}}_k\left(t\right)$ , that is, the deterioration states of QRCs, an integrated hazard function is developed based on the PHM for machine k:

$$h_k\left[t|{\mathbf{X}}_k\left(t\right)\right]={\lambda }_k\left(t\right)\mathrm{Exp}\left[{{\mathbf{\beta }}_k}^T{\mathbf{X}}_k\left(t\right)\right],k=1,\dots ,N$$

(11)

where ${\lambda }_k\left(t\right)$ represents the baseline hazard rate, depending on the operating age of the machine, and ${\mathbf{\beta }}_k$ is a vector of regression parameters defining the effects of QR deterioration on the machine-failure rate. Usually, the Weibull distribution is used to parameterize the baseline hazard rate, which is expressed as ${\lambda }_k\left(t\right)={\alpha }_k/{\alpha }_k{\eta }_k{\eta }_k{\left(t/t{\eta }_k{\eta }_k\right)}^{{\alpha }_k-1}$ .

PM will restore the machine condition to somewhere between “as good as new” and “as bad as old.” The virtual age method ²⁷ is adopted to evaluate the effect of PM on the machine-failure rate. Let $V_{k,j-1}$ denote the virtual age of machine k after the (j–1)th PM. Then, the virtual age of machine k after the jth PM can be calculated by

$$V_{\mathrm{kj}}=V_{k,j-1}+b_{\mathrm{kj}}{T}_{\mathrm{kj}},j=1,2,\dots$$

(12)

where $b_{\mathrm{kj}}(0<b_{\mathrm{kj}}<1)$ is the age reduction factor, and $V_{k0}\equiv 0$ .

Let $t$ denote the time since the last PM. Then within the jth PM cycle, the virtual age of machine k can be denoted as $V_{k,j-1}+t$ , and the deterioration states of the QRCs can be denoted as ${\mathbf{X}}_k\left(t\right)$ . Thus, the failure-rate function of machine k within the jth PM cycle is given as

$$\begin{array}{c}\multicolumn{1}{c}{h_{\mathrm{kj}}\left[t|{\mathbf{X}}_k\left(t\right)\right]={\lambda }_k(t+V_{k,j-1})\mathrm{Exp}\left[{{\mathbf{\beta }}_k}^T{\mathbf{X}}_k\left(t\right)\right]}\\ \multicolumn{1}{c}{=\frac{{\alpha }_k}{{{\eta }_k}^{{\alpha }_k}}{(t+V_{k,j-1})}^{{\alpha }_k-1}\mathrm{Exp}\left[{{\mathbf{\beta }}_k}^T{\mathbf{X}}_k\left(t\right)\right],}\\ \multicolumn{1}{c}{j=1,2,\dots }\end{array}$$

(13)

Considering the stochastic nature of ${\mathbf{X}}_k\left(t\right)$ , we can formulate $h_{\mathrm{kj}}\left(t\right)$ by the law of total probability as follows:

$$h_{\mathrm{kj}}\left(t\right)=\underset{{\mathbf{X}}_k(t)\ge 0}{\int }{h}_{\mathrm{kj}}\left[t|{\mathbf{X}}_k(t)\right]f[{\mathbf{X}}_k(t)]d{\mathbf{X}}_k(t)$$

(14)

where $f[{\mathbf{X}}_k(t)]$ is the joint probability-density function of $X_{k1}(t),\dots ,X_{k{n}_k}(t)$ . As $X_{k1}(t),\dots ,X_{k{n}_k}(t)$ are statistically independent of each other, $f[{\mathbf{X}}_k(t)]$ equals the product of the marginal probability-density functions of $X_{k1}(t),\dots ,X_{k{n}_k}(t)$ . In equation (14), by replacing $f[{\mathbf{X}}_k(t)]$ with the product of the marginal probability-density functions of $X_{k1}(t),\dots ,X_{k{n}_k}(t)$ , $h_{\mathrm{kj}}\left(t\right)$ can be derived as

$$h_{\mathrm{kj}}\left(t\right)=\frac{{\alpha }_k}{{{\eta }_k}^{{\alpha }_k}}{(t+V_{k,j-1})}^{{\alpha }_k-1}\underset{i=1}{\overset{n_k}{\Pi }}\left(\frac{{\theta }_{\mathrm{ki}}}{{\theta }_{\mathrm{ki}}-{\beta }_{\mathrm{ki}}}\right)$$

(15)

For the derivation method, please refer to Section 3.3 of Lu et al. ¹⁹ The failure-rate function will be used to calculate the minimal repair cost.

Quality-integrated maintenance modeling for individual machines

For each individual machine, the product-quality loss, minimal repair cost, and PM cost are integrated to form the total cost for PM time optimization. The PM cost includes the PM preparation cost (i.e. transportation of maintenance crew), production loss and PM operation cost (i.e. replacement of failed components). The aim of PM optimization is to find the optimal PM time that minimizes the cost per unit time over the current PM cycle, which is the ratio of the total cost to the time interval of the current cycle. The cost per unit time over a PM cycle is referred to as the cost-rate function.

The cost-rate function over the jth PM cycle for machine k is given as

$$\begin{array}{c}\multicolumn{1}{c}{C{r}_{\mathrm{kj}}\left(T_{\mathrm{kj}}\right)=\frac{\underset{0}{\overset{T_{\mathrm{kj}}}{\int }}r{Q}_k\left(t\right)\mathrm{dt}+c_k^{\mathrm{m}}\underset{0}{\overset{T_{\mathrm{kj}}}{\int }}{h}_{\mathrm{kj}}\left(t\right)\mathrm{dt}+(c^{\mathrm{s}}+c_k^{\mathrm{p}}+c^{\mathrm{e}}r{\tau }_{\mathrm{kj}})}{T_{\mathrm{kj}}+{\tau }_{\mathrm{kj}},}}\\ \multicolumn{1}{c}{j=1,2,\dots }\end{array}$$

(16)

In equation (16), $T_{\mathrm{kj}}$ denotes the time interval of the jth PM cycle for machine k. $r$ denotes the production rate of each machine, and ${\int }_0^{T_{\mathrm{kj}}}r{Q}_k\left(t\right)\mathrm{dt}$ represents the cumulative quality loss caused by the QRC deterioration of machine k over the jth PM cycle. $\underset{0}{\overset{T_{\mathrm{kj}}}{\int }}{h}_{\mathrm{kj}}\left(t\right)\mathrm{dt}$ represents the cumulative number of failures (or minimal repairs) of machine k over the jth PM cycle, and $c_k^{\mathrm{m}}\underset{0}{\overset{T_{\mathrm{kj}}}{\int }}{h}_{\mathrm{kj}}\left(t\right)\mathrm{dt}$ represents the cumulative minimal repair cost over the jth PM cycle. $(c^{\mathrm{s}}+c_k^{\mathrm{p}}+c^{\mathrm{e}}r{\tau }_{\mathrm{kj}})$ represents the PM-related cost for machine k, where $c^{\mathrm{s}}$ denotes the PM preparation cost, $c_k^{\mathrm{p}}$ denotes the PM operation cost, and $c^{\mathrm{e}}r{\tau }_{\mathrm{kj}}$ denotes the production loss during the jth PM operation. Note that for the serial MMS, the production rate of each machine is the same and equals that of the system.

From equation (13), it can be seen that the machine-failure rate continuously evolves with PM interventions. As a result, the cost-rate function continuously evolves with PM interventions. For the jth PM cycle, the optimal PM interval can be obtained by minimizing $C{r}_{\mathrm{kj}}\left(T_{\mathrm{kj}}\right)$ in equation (16). As a result, the PM scheduling for each individual machine is dynamic and progressive.

Maintenance-operation prioritization (MOP) for system-level scheduling

The optimal PM intervals of the individual machines are fed into the system-level scheduling. The principle of system-level scheduling is illustrated in Figure 2.

Figure 2.

System-level scheduling procedure.

In system-level scheduling, it is necessary to first determine when the PM will be initiated. PM is initiated when the scheduled PM time of a machine arrives, or the quality loss of the final products reaches the QLT, whichever comes first.

Let $t^{\mathrm{pm}}$ denote the time when the PM is initiated for the system, $t_k^s$ denote the scheduled PM time of machine k (which depends on the optimal PM intervals of machine k), and $t^{\mathrm{q}}$ denote the time when the final-product quality loss reaches the QLT. Then, we have

$$t^{\mathrm{pm}}=\min \{t^{\mathrm{q}},t_k^s:k=1,...,N\}$$

(17)

The time when the final-product quality loss reaches the QLT can be calculated, based on the final-product quality loss function and the QLT. Let $C^{\mathrm{q}}(C^{\mathrm{q}}>0)$ denote the QLT, and $Q_N\left(t\right)$ denote the final-product quality loss function. Then, we have

$$Q_N\left(t^{\mathrm{q}}\right)=C^{\mathrm{q}}$$

(18)

The final-product quality loss is caused by deviations of the KPCs of the final products. The deviations of the KPCs of the final products, denoted as ${\mathbf{Y}}_N[t|{\mathbf{X}}_1(t),...,{\mathbf{X}}_N(t)]$ , can be derived based on equation (5). Based on Taguchi’s loss function, the expected quality loss caused by the deterioration of the KPCs of the final products can be evaluated as follows:

$$\begin{array}{c}\multicolumn{1}{c}{Q_N(t)=E[{\mathbf{Y}}_N^T[t|{\mathrm{X}}_1(t),...,{\mathbf{X}}_N(t)]}\\ \multicolumn{1}{c}{\mathbf{q}{\mathbf{Y}}_N[t|{\mathbf{X}}_1(t),...,{\mathbf{X}}_N(t)]]}\\ \multicolumn{1}{c}{=E[\sum _{s=1}^m{Y_{\mathrm{sN}}}^2[t|{\mathbf{X}}_1(t),...,{\mathbf{X}}_N(t)]]}\end{array}$$

(19)

where $\mathbf{q}=\mathrm{diag}\{q_1,...,q_m\}$ is a diagonal matrix consisting of the quality-loss coefficients of all KPCs.

If the PM is triggered by the QLT, the quality loss of the final products must be reduced to below the QLT by performing PM on the machines. Based on the SOD effect, the PM of different machines will lead to different degrees of improvement in the final-product quality. When selecting machines for PM, it is desirable to give higher priority to the machines with greater quality-improvement abilities to best improve the final-product quality. Thus, a quality-improvement factor is proposed to evaluate the ability of each machine to improve the product quality though PM.

This factor is defined for a machine as the reduction of the final-product quality loss brought by the PM of this machine. Under this definition, the quality improvement factor for machine k is given by

$$I_q^k\left(t\right)=Q_N\left(t\right)-Q_N^k\left(t\right),k=1,\dots ,N$$

(20)

where $Q_N^k\left(t\right)$ represents the quality loss of the final products after the PM of machine k at time t. As the PM will restore the QRCs of machine k to completely new states, $Q_N^k\left(t\right)$ can be obtained using equation (19) by setting ${\mathbf{X}}_k(t)$ to 0.

The greater the quality-improvement factor for a machine, the greater the ability of the machine in quality improvement through PM. Therefore, the machine with the greatest quality-improvement factor is given the highest priority, and the machine with the least quality-improvement factor is given the lowest priority. Then, a priority queue $P_1,...,P_N$ can be obtained based on the quality improvement factor $I_q^k\left(t\right),k=1,..,N$ , where $P_1$ represents the machine with the highest priority; and $P_N$ represents the machine with the lowest priority. The first g machines in the priority queue, whose PM can reduce the quality loss of final products to below the QLT, will be selected for PM. Let $G^{\mathrm{q}}=\{P_1,...,P_g\}$ denote the first g machines in the priority queue. Then, $G^{\mathrm{q}}$ can be determined by

$$G^{\mathrm{q}}=\{\begin{array}{c}\{P_1\},Q_N^{P_1}\left(t\right)<C^{\mathrm{q}}\hfill \\ \multicolumn{1}{c}{\{P_1,P_2,\cdots ,P_g\},Q_N^{P_1,...,P_g}\left(t\right)<C^q<Q_N^{P_1,...,P_{g-1}}\left(t\right),g\ge 2}\end{array}$$

(21)

Each time PM is initiated for the system, it is desirable to select a group of machines for PM. This is because jointly maintaining multiple machines can save common PM preparation costs and system production losses compared with maintaining them separately. The maintenance time window (MTW) method ⁶ is adopted to select a group of machines for PM. When PM is initiated at $t^{\mathrm{pm}}$ , the machines whose scheduled PM times $t_k^{\mathrm{s}}$ fall within $[t^{\mathrm{pm}},t^{\mathrm{pm}}+T^{\mathrm{w}}]$ will be maintained together, where $T^{\mathrm{w}}$ is the MTW.

If PM is triggered by the QLT, the machines in $G^{\mathrm{q}}$ and the machines whose scheduled PM times fall within the MTW will be maintained simultaneously. The decision results at the system-level scheduling are fed back to the machine level for the next cycle scheduling. The optimal QLT and MTW can be obtained by minimizing the total cost of the system over the planning horizon.

QMM-MOP method for interactive bi-level scheduling

Combing the machine-level and system-level scheduling, a QMM-MOP methodology for interactive bi-level scheduling is formed. The flowchart of the interactive bi-level scheduling is illustrated by Figure 3.

• Step 1 Sequential PM pulling at the machine level. Start the maintenance scheduling for each machine from the first PM cycle $j_k=1$ . Pull the original PM intervals $T_{k{j}_k}^{*}(k=1,...,N)$ from the QMM for each machine.

• Step 2 Time point determination. From the first system-level PM cycle $J=1$ , determine the originally scheduled PM time for each machine in the system-level scheduling by

$$t_{\mathrm{kJ}}^s=T_{k{j}_k}^{*},k=1,\dots ,N$$

(22)

• Step 3 PM initiation determination. PM is initiated when the scheduled PM time of a machine arrives, or the quality loss of final products reaches the QLT. The quality loss of final products within the Jth system-level PM cycle can be evaluated by

$$\begin{array}{c}\multicolumn{1}{c}{Q_{\mathrm{NJ}}\left(t\right)=E}\\ \multicolumn{1}{c}{\left\{\sum _{s=1}^m{q}_s{Y}_{\mathrm{sN}}^2[t|{\mathbf{X}}_1(t-t_{1,j_1-1}^a),...,{\mathbf{X}}_N(t-t_{N,j_N-1}^a)]\right\},}\\ \multicolumn{1}{c}{t_{J-1}^{\mathrm{pm}}\le t\le t_J^{\mathrm{pm}}}\end{array}$$

(23)

Figure 3.

Flowchart of interactive bi-level scheduling for the MMS.

where $t_J^{\mathrm{pm}}$ denotes the time when the Jth PM is initiated, and $t_{k,j_k-1}^a(t_{k0}^a\equiv 0)$ denotes the actual PM time of machine k. $t^{\mathrm{q}}$ can be derived based on

$$Q_{\mathrm{NJ}}\left(t^{\mathrm{q}}\right)=C^{\mathrm{q}}$$

(24)

Determine the time when the Jth PM is initiated for the system by

$$t_J^{\mathrm{pm}}=\min \{t^{\mathrm{q}},t_{\mathrm{kJ}}^s:k=1,...,N\}$$

(25)

If $t_J^{\mathrm{pm}}=t^{\mathrm{q}}$ , the PM is triggered by the QLT, and go to step (4); otherwise, go to step (5).

• Step 4 Maintenance operation prioritization (MOP). Calculate the quality-improvement factor for all machines (denoted as $I_k^{\mathrm{q}}\left(t\right),k=1,...,N$ ) using equation (20), and then determine the first g machines in the priority queue (denoted as $G^q$ ) using equation (21). Go to step (5).

• Step 5 Maintenance operation combination and execution. Use MTW to combine PM operations, and then determine the group of machines for PM based on

$$G_J=\{\begin{array}{c}{G}^{\mathrm{q}}\cup \left\{k|t_J^{\mathrm{pm}}\le t_{\mathrm{kJ}}^{\mathrm{s}}\le t_J^{\mathrm{pm}}+T^{\mathrm{w}}:k=1,...,N\right\}\mathrm{if}{t}_J^{\mathrm{pm}}=t^{\mathrm{q}}\hfill \\ \multicolumn{1}{c}{\left\{k|t_J^{\mathrm{pm}}\le t_{\mathrm{kJ}}^{\mathrm{s}}\le t_J^{\mathrm{pm}}+T^{\mathrm{w}}:k=1,...,N\right\}\mathrm{if}{t}_J^{\mathrm{pm}}\ne t^{\mathrm{q}}}\end{array}$$

(26)

Execute PM operations on the machines in the group $G_J$ and update the statuses of the maintained machines $k\in G_J$ , that is.

$$\{\begin{array}{c}{j}_k=j_k+1\hfill \\ \multicolumn{1}{c}{t_{k,j_k-1}^a=t_J^{\mathrm{pm}}}\\ \multicolumn{1}{c}{V_{k,j_k-1}=V_{k,j_k-1}+b_{k,j_k-1}(t_{k,j_k-1}^a-t_{k,j_k-2}^a-{\tau }_{k,j_k-2})}\end{array},k\in G_J$$

(27)

Moreover, let ${\tau }_J({\tau }_0\equiv 0)$ denote the PM duration for the group $G_J$ . Assume that the PM operations of multiple machines can be implemented in parallel. Under this assumption, it can be known that the PM duration of a group equals to the maximum of PM durations of machines in the group. That is ${\tau }_J=\max \{{\tau }_{k{j}_k}:k\in G_J\}$ .

• Step 6 Planning horizon check. Check whether the scheduling time $T_s=t_J^{\mathrm{pm}}$ is greater than or equals to the planning horizon, $\mathrm{PH}$ . If yes, go to step (8) and terminate the PM scheduling. If no, go to step (7) and continue the PM scheduling.

• Step 7 Time point update. Fed the updated statuses of maintained machines back to machine-level scheduling and pull the original PM intervals $T_{k{j}_k}^{*}(k\in G_J)$ from the QMM for these machines. For the next system-level PM cycle, that is, $J=J+1$ , update the scheduled PM time of each machine by

$$t_{\mathrm{kJ}}^{\mathrm{s}}=\{\begin{array}{c}{t}_{J-1}^{\mathrm{pm}}+{\tau }_{J-1}+T_{k{j}_k}^{*},k\in G_{J-1}\hfill \\ \multicolumn{1}{c}{t_{k,J-1}^{\mathrm{s}}+{\tau }_{J-1},k\notin G_{J-1}}\end{array}$$

(28)

Then go back to step (3) to determine the time when the next system-level PM is initiated.

• Step 8 Maintenance scheduling termination. End the maintenance scheduling and output the PM schedules for the system, represented as $(t_J^{\mathrm{p}},G_J),J=1,...,W$ .

Based on the PM schedules of the system, the total cost of the system, including product-quality loss, minimal repair cost, and PM cost, over the planning horizon can be evaluated as follows:

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{TC}=\sum _{J=1}^W\underset{t_{J-1}^{\mathrm{pm}}+{\tau }_{J-1}}{\overset{t_J^{\mathrm{pm}}}{\int }}r{Q}_{\mathrm{NJ}}\left(t\right)\mathrm{dt}+\sum _{k=1}^N\sum _{j_k=1}^{w_k+1}}\\ \multicolumn{1}{c}{c_k^{\mathrm{m}}\underset{0}{\overset{t_{k{j}_k}^{\mathrm{a}}-t_{k,j_k-1}^{\mathrm{a}}}{\int }}{h}_{k{j}_k}\left(t\right)\mathrm{dt}+\sum _{J=1}^W(c^{\mathrm{s}}+c^{\mathrm{e}}r{\tau }_J+\sum _{k\in G_J}{c}_k^{\mathrm{p}})}\end{array}$$

(29)

where $w_k$ is the number of PM operations scheduled for machine k over the planning horizon, which satisfies $t_{k,w_{k+1}}^{\mathrm{a}}\equiv \mathrm{PH}$ ; and $W$ is the number of PM activities scheduled for the system over the planning horizon, which satisfies $t_W^{\mathrm{pm}}\equiv \mathrm{PH}$ . $\mathrm{TC}$ is a function of $C^{\mathrm{q}}$ and $T^{\mathrm{w}}$ . For each value of ( $C^{\mathrm{q}}$ , $T^{\mathrm{w}}$ ), the corresponding PM schedules for the system can be derived using the interactive bi-level scheduling procedure, and then $\mathrm{TC}$ can be evaluated with equation (29). The optimal $C^{\mathrm{q}}$ and $T^{\mathrm{w}}$ can be determined by minimizing $\mathrm{TC}$ .

From Figure 3, it can be obtained that the frequency count (i.e. execution times of statements) of the QMM-MOP method is T(N) = [(6L + 8)N + 4] m, where L denotes the number of candidates for the optimal PM time of each machine, and m denotes the number of PM cycles for the system. When the number of machines is considered as the variable, the time complexity of the proposed method is O(N). Apparently, the time complexity of the proposed method is linear with the number of machines. Therefore, the proposed method will not lead to large computational cost when applied to large scale MMSs. The proposed method adapts to large-scale systems.

Case study

Case overview

In this section, a case study is conducted to illustrate the effectiveness of the proposed QMM-MOP methodology. The case was cited from Lu and Zhou. ¹⁵ The case deals with a serial four-stage machining system employed to produce a kind of shaft sleeves. The shaft sleeve is illustrated in Figure 4. A shaft sleeve is completed through four stages of machining operations, as described in Table 1. The QRCs and the noise variables in each stage are presented in Table 2.

Figure 4.

Schematic illustration of the shaft sleeve.

Table 1.

Description of the machining operations.

Stage	Machine	Operation	KPCs
1	Milling machine	Mill Face B	1st KPC: distance between Face B and Face A
2	Lathe	Turn Outer cylinder D	2nd KPC: distance between Face C and Face B
			3rd KPC: diameter of Outer cylinder D
3	Boring machine	Bore Hole E	4th KPC: diameter of Hole E
4	Broaching machine	Broach Keyway F	5th KPC: distance between the top of Keyway F and bottom of Hole E;
			6th KPC: width of Keyway F

Table 2.

The QRCs and the noise variables in each stage.

Stage	Deterioration of QRCs of machines	Noise variables
1	X 11(t): the wear of the mill cutter	Z 11: the locating error
2	X 21(t): the wear on the major frank of the turning tool	Z 21: the locating error
	X 22(t): the wear on the minor frank of the turning tool	Z22: the tool setting error
3	X 31(t): the wear of the boring tool	Z 31: the tool setting error
		Z 32: the vibration of the boring tool
4	X 41(t): the wear along the height direction of the broach	Z 41: the locating error
	X 42(t): the wear along the width direction of the broach	Z 42: the expansion of the broach

The parameters of the SOD models are presented in equations (30)–(33). The parameters of the model errors are: σ₁₁ = 0.00082, σ₂₂ = 0.0008, σ₃₂ = 0.00085, σ₄₃ = 0.0009, σ₅₄ = 0.00086, σ₆₄ = 0.0008. The parameters for the noise variables in the four stages are: ${\mathrm{\Sigma }}_1^z=0.00171m{m}^2$ , ${\mathrm{\Sigma }}_2^z=\mathrm{diag}(0.00142,0.00126)m{m}^2$ , ${\mathrm{\Sigma }}_3^{\mathrm{z}}=0.0029m{m}^2$ , ${\mathrm{\Sigma }}_4^{\mathrm{z}}=\mathrm{diag}(0.0013,$ $0.00236)m{m}^2$ . Based on the SOD model, the deviations of KPCs of the final products can be derived by equation (5). Then, the final-product quality loss function can be derived by equation (19). The values for other parameters are shown in Table 3. In addition, q  = diag (2, 4, 16, 20, 5, 10) $/mm ² , r = 330 unit/day, and PH = 360 day. It is assumed that the age reduction factor of each PM activity for the same machine is the same. That is, $b_{k1}=b_{k2}=\cdots =b_{k{n}_k}$ . Note that the planning horizon in this case is different from that in the study. ¹⁵

Table 3.

Parameter values for the case.

k	${\mathbf{\mu }}_k$	${\mathbf{\theta }}_k$ (mm −1)	${\mathbf{\beta }}_k$	${\alpha }_k$	${\eta }_k$ (day)	$b_k$	${\tau }_k^{\mathrm{p}}$ (day)	$c_k^{\mathrm{p}}$ ($)	$c_k^{\mathrm{m}}$ ($)	$c^{\mathrm{s}}$ ($)	$c^{\mathrm{e}}$ ($/unit)
1	0.00142	3.5	0.015	2.3	120	0.08	0.5	3800	3700
2	(0.00163, 0.00166)	(3.3, 4.5)	(0.01, 0.02)	2.5	88	0.085	0.5	3530	3680	230	15
3	0.00163	4.2	0.019	2.8	126	0.083	0.5	3610	3550
4	(0.00149, 0.00318)	(4.4, 5.2)	(0.016, 0.011)	2.6	135	0.079	0.5	3840	3600

$$\begin{array}{c}\mathrm{Stage}1:\hfill \\ Y_{11}(t)=0.00036+0.987X_{11}(t)+0.935Z_{11}\hfill \\ \multicolumn{1}{c}{-0.0082X_{11}(t)Z_{11}+{\epsilon }_{11}}\\ Y_{21}(t)=0,Y_{31}(t)=0,Y_{41}(t)=0,Y_{51}(t)=0,\hfill \\ \multicolumn{1}{c}{Y_{61}(t)=0}\end{array}$$

(30)

$$\begin{array}{c}\mathrm{Stage}2:\hfill \\ Y_{22}(t)=0.00034-0.985X_{21}(t)+0.991Y_{11}(t)\hfill \\ \multicolumn{1}{c}{-0.926Z_{21}+0.0074X_{21}(t)Z_{21}+{\epsilon }_{22}}\\ Y_{32}(t)=-0.0004+2.012X_{22}(t)+1.923Z_{22}\hfill \\ \multicolumn{1}{c}{+0.0085X_{22}(t)Z_{22}+{\epsilon }_{32}}\\ Y_{12}(t)=Y_{11}(t),Y_{42}(t)=0,Y_{52}(t)=0,Y_{62}(t)=0\hfill \end{array}$$

(31)

$$\begin{array}{c}\mathrm{stage}3:\hfill \\ Y_{43}(t)=0.00032-1.984X_{31}(t)+0.843Z_{31}\hfill \\ \multicolumn{1}{c}{-0.0382X_{31}(t)Z_{31}+{\epsilon }_{43}}\\ Y_{13}(t)=Y_{12}(t),Y_{23}(t)=Y_{22}(t),Y_{33}(t)=Y_{32}(t),\hfill \\ \multicolumn{1}{c}{Y_{53}(t)=0,Y_{63}(t)=0}\end{array}$$

(32)

$$\begin{array}{c}\mathrm{stage}4:\hfill \\ Y_{54}(t)=-0.00052-0.988X_{41}(t)-0.707Y_{33}(t)\hfill \\ \multicolumn{1}{c}{+0.501Y_{43}(t)+0.945Z_{41}-0.0081X_{41}(t)Z_{41}+{\epsilon }_{54}}\\ Y_{64}(t)=0.0003-0.993X_{42}(t)+0.674Z_{42}\hfill \\ \multicolumn{1}{c}{+0.0085X_{42}(t)Z_{42}+{\epsilon }_{64}}\\ Y_{14}(t)=Y_{13}(t),Y_{24}(t)=Y_{23}(t),Y_{34}(t)=Y_{33}(t),\hfill \\ \multicolumn{1}{c}{Y_{44}(t)=Y_{43}(t)}\end{array}$$

(33)

$$\begin{array}{c}\multicolumn{1}{c}{{\mathbf{\mu }}_k=[{\mu }_{k1},...,{\mu }_{k{n}_k}],{\mathbf{\beta }}_k=[{\beta }_{k1},...,{\beta }_{k{n}_k}],}\\ \multicolumn{1}{c}{{\mathbf{\theta }}_k=[{\theta }_{k1},...,{\theta }_{k{n}_k}].}\end{array}$$

Obtainment of optimal PM schedule

The decision parameters of the proposed methodology are $C^{\mathrm{q}}$ and $T^{\mathrm{w}}$ . The optimal ( $C^{\mathrm{q}}$ , $T^{\mathrm{w}}$ ) and the corresponding optimal PM schedules can be obtained by minimizing the total cost TC. By varying the values of ( $C^{\mathrm{q}}$ , $T^{\mathrm{w}}$ ), the minimum value for TC can be obtained, based on the interactive bi-level scheduling procedure described in Section 3.3. It should be noted that when $C^{\mathrm{q}}$ is greater than or equals to 1.3, no PM is triggered by the QLT. When $T^{\mathrm{w}}$ is greater than or equals to 80, the PM operations of all machines are combined. Thus, the search ranges for $C^{\mathrm{q}}$ and $T^{\mathrm{w}}$ are set as $0<C^{\mathrm{q}}\le 1.3$ and $0<T^{\mathrm{w}}\le 80$ . The search-step lengths for $C^{\mathrm{q}}$ and $T^{\mathrm{w}}$ are 0.01 and 1, respectively. By a numerical search, the minimum value of the total cost is obtained as 200759. The corresponding optimal values for $C^{\mathrm{q}}$ and $T^{\mathrm{w}}$ are 1.13 and 28, respectively.

The optimal PM schedules are listed in Table 4. It can be found that seven PM activities are scheduled for the system. In each PM activity, multiple machines will be jointly maintained. This implies that the PM preparation cost and production loss will be reduced. For instance, in the second PM activity, four machines will be jointly maintained, which implies that three PM preparation costs and production losses will be avoided.

Table 4.

Optimal PM schedules for the entire system.

PM cycle ( $J$ )	1	2	3	4	5	6	7
PM time ( $t_J^{\mathrm{pm}}$ )	54.9	101.4	156.2	202.4	257.2	302.4	357.2
Machine 1		▪		▪		▪
Machine 2	▪	▪	▪	▪	▪	▪	▪
Machine 3	▪	▪	▪	▪	▪	▪	▪
Machine 4		▪		▪		▪

‘▪’ indicates that the machine in the row will be preventively maintained in the system PM cycle.

Table 5 shows the priority queue in the cases where the QLT triggers the PM. From Table 5, it can be observed that five PM activities are triggered by the QLT. This implies that the QLT plays an effective role in controlling the final-product quality. It can be also found that the priority queue is 3-2-4-1 in the five cases where the QLT triggers the PM. This indicates that the priority queue may largely depend on the physical properties of the serial MMSs. That is, the PM of some machines may inherently lead to greater improvement in the final-product quality in a serial MMS.

Table 5.

Priority queue in the cases where the QLT triggers PM.

PM cycle ( $J$ )	1	2	3	4	5	6	7
PM time ( $t_J^{\mathrm{pm}}$ )	54.9	101.4	156.2	202.4	257.2	302.4	357.2
P 1	3	3	3	Nan	3	Nan	3
P 2	2	2	2	Nan	2	Nan	2
P 3	4	4	4	Nan	4	Nan	4
P 4	1	1	1	Nan	1	Nan	1

Sensitivity analysis on the quality-loss coefficient

The quality-loss coefficient q is the most important parameter for the proposed methodology. Hence, a sensitivity analysis was carried out on this parameter. Table 6 shows the maintenance decision results for different values of q . In Table 6, V^q denotes the scaling ratio for q , M_s denotes the number of PM activities for the system, and M_k (k = 1, 2, 3, 4) denotes the number of PM operations for machine k.

Table 6.

Maintenance decision results under different q values.

V	Decision parameters		Number of PM activities/operations
	$C^{\mathrm{q}}$	$T^{\mathrm{w}}$	Ms	M 1	M 2	M 3	M 4
1	1.13	28	7	3	7	7	3
2	1.64	44	9	3	10	10	5
3	2.23	22	11	3	11	11	5
4	2.57	33	13	4	13	13	6
5	3.21	28	13	4	13	13	6

From Table 6, it can be observed that the optimal QLT increases with the increase of q ; however, the increasing rate of QLT is smaller than that of q . This means that with the increase of q , the quality control will become increasingly strict, and more PM activities will be triggered by the QLT. In fact, the results demonstrate that the number of PM activities for the system increases with q .

Table 6 also shows that, with an increase in q , the number of PM operations for the four machines all exhibit an increasing trend. This is because, with the increase in quality-loss coefficients, more PM operations are required to decrease the product-quality deterioration and thus prevent excessive quality loss.

Comparison with existing methodologies

In this section, the proposed QMM-MOP method is compared with two existing types of methods to demonstrate its superiority. As discussed in the introduction, the first type of method is the conventional PM models, and the second is the existing quality-integrated PM models.

For the conventional PM models, the MTW method is chosen for comparison with the proposed method. The only difference between them is whether product-quality improvement is considered in PM scheduling. In the MTW method, the MTW is used to group the PM operations in system-level scheduling. The PM schedules under the MTW method can be derived using the proposed QMM-MOP method by setting the quality-loss coefficients to 0. This is because, when the quality-loss coefficient is 0, product-quality deterioration will not impact PM scheduling.

For the existing quality-integrated PM models, the QMM-MTW method is chosen for comparison with the proposed method. The only difference between them is whether they consider the fact that the PM of different machines will result in different degrees of improvement in the final-product quality. In the QMM-MTW method, the QMM is adopted to obtain the optimal PM intervals for individual machines, and the MTW is used to group PM operations in system-level scheduling. The PM schedules under the QMM-MTW method can be derived using the proposed QMM-MOP method by setting QLT as +∞. When QLT equals +∞, the PM will never be triggered by the QLT, and thus MOP will never be conducted.

Table 7 shows a comparison of the three methods in terms of PM schedules. It can be found that the QMM-MOP method schedules more PM operations for all machines compared with the MTW method. This is because, the QMM-MOP method considers product-quality improvement in PM scheduling, and thus schedules more PM operations to reduce the product-quality deterioration. In addition, the QMM-MOP method schedules more PM operations for machine 3 than the QMM-MTW method. This is because, the QMM-MOP method conducts MOP based on the quality-improvement factor. As machine 3 is always given the highest priority, it is assigned more PM operations under the QMM-MOP method.

Table 7.

Comparison of PM schedules for the three methods.

Methods/PM schedule	Number of PM activities/operations
	Ms	M 1	M 2	M 3	M4
MTW method	6	3	4	3	3
QMM-MTW method	7	3	7	6	3
QMM-MOP method	7	3	7	7	3

Table 8 shows a comparison of the three methods in terms of cost performance. It can be found that the QMM-MOP method generates higher PM-related costs, but lower product-quality losses and minimal repair costs than the MTW method. The higher PM-related cost is because more PM operations are scheduled by the QMM-MOP method. More PM operations result in greater reduction in product-quality loss and minimal repair cost, which consequently results in a reduced total cost. Specifically, the proposed QMM-MOP method can achieve a 42.33% reduction in product-quality loss and a 19.57% reduction in total cost, compared with the MTW method.

Table 8.

Comparison of cost performance of the three methods.

Methods\costs	PM-related cost	Minimal repair cost	Quality loss	Total cost
MTW method	64,100	53,964	131,537	249,601
QMM-MTW method	88,225	33,752	80,320	202,297
QMM-MOP method	91,835	33,067	75,857	200,759

In addition, the QMM-MOP method generates higher PM-related costs than does the QMM-MTW method. This is because, the QMM-MOP method schedules more PM operations for machines with a greater ability to improve the final-product quality (through PM). More PM operations result in significant reduction in product-quality loss and consequently lead to a reduced total cost. Specifically, the proposed QMM-MOP method can achieve a 5.56% reduction in product-quality loss and a 0.76% reduction in total cost, compared with the QMM-MTW method.

Figure 5 shows the comparison of the three methods in terms of product-quality loss and the total cost under different values of q . The three methods generate the same product-quality loss and total cost in the special case where q equals 0. In general cases where q is greater than 0, the QMM-MOP method has a significantly lower product-quality loss and total cost than the MTW method. This demonstrates that it is beneficial to consider product-quality improvement in the PM scheduling of MMSs.

Figure 5.

Comparison of the three methods in terms of product-quality loss and total cost under different q values.

In addition, the QMM-MOP method has an obviously lower product-quality loss and a slightly lower total cost than the QMM-MTW method in general cases. This indicates that the QMM-MOP method can ensure a higher final-product quality with a lower total cost, compared with the QMM-MTW method. This suggests that conducting MOP, based on the quality improvement factor, can better improve the final-product quality.

Conclusions and future research

This paper proposes a QMM-MOP methodology to deal with the challenges posed by the SOD in PM scheduling of MMSs. This methodology adopts an interactive bi-level scheduling framework. In machine-level scheduling, a QMM is developed by integrating intermediate-product quality deterioration into the total cost to schedule timely PM for each individual machine. In system-level scheduling, MOP will be conducted, based on a SOD-enabled quality-improvement factor, to select machines for PM to best improve the final-product quality. The results of the machine-level scheduling and system-level scheduling interact dynamically.

A case study shows that the PM of some machines inherently leads to a greater improvement in the final-product quality, owing to the physical property of a serial MMS. In addition, an increase in the quality-loss coefficient drives more PM operations to be scheduled for all machines to decrease the product-quality deterioration. The QMM-MOP method is compared with two existing methods. The results demonstrate that compared with the two exiting methods, the QMM-MOP method can reduce the final-product quality loss and total cost. Two insights can be drawn for practitioners. First, it is beneficial to consider product-quality improvement in the PM scheduling of MMSs. Second, performing MOP, based on the quality improvement factor, can better improve the final-product quality.

This study assumes that the degradation processes of QRCs are stationary gamma processes. In practice, there are cases where the degradation processes of QRCs are non-stationary stochastic process. Thus, it is desirable to assume non-stationary gamma processes for the degradation processes of QRCs in future work. Under this assumption, it is necessary to dynamically update the gamma-process parameters. The parameters can be dynamically updated using the online degradation data of QRCs from sensors. In this paradigm, the maintenance model should be changed accordingly. In addition, the proposed maintenance model is for serial MMSs. In practice, more complex MMSs exist, such as serial-parallel MMSs. It is desirable to extend the proposed maintenance model to more complex MMSs in future work.