Integrated production planning and condition-based maintenance considering uncertain demand and random failures

Abstract

This article proposes a proactive approach for the integrated problem of production planning and condition-based maintenance under uncertain environment, that is, the exogenous uncertainty of demand and the endogenous uncertainty of failure. Considering that the production system has limited capacity and suffers degradation with usage, an integration model of production planning and condition-based maintenance policy is proposed, in which the maintenance threshold and production quantity are proactively decided simultaneously. The robust and stable solution of the complicated non-linear model is obtained by a simulation-based multi-population genetic algorithm. Numerical results not only show the impacts of different factors on the model, but also prove the superiority of the developed model and algorithm in uncertain environment.

Introduction

In fact, the manufacturing systems always show a natural interdependence between production and maintenance activities, which raises the problem of how to jointly optimize the production and maintenance to avoid conflicts and minimize total cost. In particular, it might be more complicated to optimize the joint problem in an uncertain environment, which is a problem that manufacturers must face in pursuit of competitiveness and efficiency under increasing market pressure. In the last decades, Aghezzaf et al. ¹ first developed an integrated model of capacitated lot-sizing problem (CLSP) and interval-based preventive maintenance (PM) plan for a single component, and the integrated model had been demonstrated to be cost saving. Subsequently, some extensions based on CLSP and PM have been studied. Fitouhi and Nourelfath ² and Lu et al. ³ modeled the integrated problem with run-based PM plan and obtained further cost saving. An iterative method was proposed in Zhao et al. ⁴ with consideration of the expected consumption of capacity by random failures. Recently, Wang et al. ⁵ explored the integration problem of CLSP and imperfect PM with consideration of nonconforming items in production. These studies generally employ the expected failure numbers to evaluate the impact of random failures, which, however, ignore the stability of the plan in reality.

Lots of uncertainties exist in the tactical production planning, including but not limited to uncertain demand and random failures. Both uncertain demand and random failures can have an effect on the service level. On the one hand, the random failures caused by system degradation may lead to low productivity and increased downtime. On the other hand, the uncertain demand would increase uncompleted orders and lower the service level. Appropriate maintenance actions and production plan could improve the system availability and keep competitiveness of the manufacturer. Thus, a robust production plan and corresponding maintenance policy are needed to reduce the impacts of the uncertainties and minimize the total cost.

To deal with these uncertainties, many production systems have adopted the proactive approaches, which develop an initial plan to guide production activities by incorporating the knowledge of the uncertainty at the decision stage, and adjust the initial plan through reactive procedures to maintain its feasibility and performance when unforeseen situations occur. Generally, a deterministic plan constructed through predictive approach always has an optimum expected performance, but it is not likely to be the optimum one once the uncertainties are realized. So, the objective of the reactive procedures is to generate a feasible plan to the present situation that deviates as little as possible from the initial baseline plan, which is also referred to as the stability objective.

Several studies have been reported to deal with the failure uncertainty in production problem proactively. Goren and Sabuncuoglu ⁶ generated robust and stable schedules in a single-machine environment subject to machine breakdowns. The authors developed two surrogate measures for robustness and stability, which considered both busy and repair time distributions. Al-Hinai and ElMekkawy ⁷ addressed the problem of finding robust and stable solutions for the flexible job shop scheduling problem with random machine breakdowns. Cui et al. ⁸ and Lu et al. ⁹ proactively decided the production scheduling and PM simultaneously to optimize the bi-objective of robustness and stability considering that the breakdowns would affect the stability of the machine. However, relatively little attention has been given to failure uncertainty in tactical level. On the other hand, the demand uncertainty in tactical production planning has received much attention and gained a lot of insights. Jing et al. ¹⁰ developed a fuzzy mixed-integer linear programming model, which considers the uncertainties of market demands, to address a capacitated dynamic lot-sizing problem with remanufacturing. Guan et al. ¹¹ proposed a branch and cut algorithm to solve stochastic uncapacitated lot-sizing problems under demand uncertainty. Aghezzaf et al. ¹² presented and discussed three different models to generate robust tactical production plans in a multi-stage production system considering the random demand. An extensive review of production planning under uncertainty can be found in Mula et al. ¹³

For a degradation system, its condition can be monitored continuously or discretely. Based on the condition, proper maintenance can be employed to restore the system, and this is called condition-based maintenance (CBM). The problems of modeling and optimizing CBM have been widely investigated in literature. Among them, the control limits and/or inspection intervals of CBM have been optimized in many models in order to minimize the long-run cost.^14–17 Cox’s proportional hazards model (PHM) is one of the most commonly used models to describe degrading system in the CBM framework, which considers both system age and condition variables. Some works have been done in CBM to apply the PHM to derive the optimal control-limit maintenance policy.^18–22 Compared with PM, CBM monitoring the actual system condition is more effective in reducing cost and uncertainty of the stochastic degrading system.^23,24

The purpose of this article is to jointly optimize the CBM policy and production planning for a degrading system considering uncertain demand and random failures. The motivation for this work is the need for better production planning and maintenance policy under the uncertain environment in real industrial applications. In fact, it is more beneficial to make real-time PM decisions based on the inspection result every period. The main contribution of this article is to develop a proactive approach for the integration problem of production planning and CBM to obtain a robust plan under uncertain environment. As far as we know, no research integrating production and maintenance has considered the fact that the uncertainty of failure and demand can affect the stability of the production plan. In addition, a multi-population genetic algorithm (MPGA) with discrete event simulation is developed to solve this stochastic programming problem.

The remainder of this article is organized as follows. “Problem statement” presents the statement of problem. The formulation of integration model is developed in “Model formulation.” The solution method is proposed in “Solution algorithm.”“Numerical results and discussion” presents a case study and some comparison experiments to demonstrate the performance of the model. Conclusions are drawn in “Conclusion.”

Problem statement

The uncertainty in tactical planning problems has been considered by many manufacturing plants. Suppose that a production system is required to produce a set of products during a finite time horizon. The demand of every product is assumed to be mutually independent, and the probability distribution of the demands can be evaluated according to their past experience. Furthermore, assuming that the system is subject to failure, PM is needed to reduce failures and improve availability. Therefore, managers should not only determine the production quantity of each product in each period, but also determine the corresponding maintenance actions to minimize production and maintenance costs under the uncertain environment.

The PHM is used to describe the degradation of the system and calculate the system hazard rate. The system is assumed to degrade with operating time. Thus, random failure in production depends on the production decisions, which means that it becomes endogenous uncertainty. ²⁵ The relationship between system degradation and production quantity will be discussed in detail in the next section. Whenever a random failure occurs, a corrective maintenance (CM) is conducted immediately to restore it without changing the system’s usage and condition. Continuous monitoring of the system is impossible due to economic and technological constraints. So, inspections and PM decisions are made at the beginning of every period, and the inspections are assumed to be instantaneous and perfect. Based on the inspection information, the hazard rate is calculated. If the hazard rate exceeds a predetermined threshold, a PM is carried out to restore the system to the “as good as new” state.

In practice, exact demand and failures are discovered only after the production plan is put into execution. Unsatisfied demand will occur when the system capacity is not enough to execute the initial plan or the real demand is increased. Meanwhile, redundant items will be produced when the real demand is reduced. Managers must react quickly to the changes by carrying additional storage or increasing capacity by outsourcing, which would result in additional costs.

Model formulation

Degradation model and maintenance model

Hazard rate is widely used to represent the condition of the system in reliability theory. A PHM model which combines a Weibull baseline hazard rate function and system state covariate is concerned in this batch production model. The hazard rate is a product of a Weibull probability distribution function of the system usage and an exponential function of the degradation process, which is shown as follows

$$h\left(t,D\left(t\right)\right)=\frac{\delta }{\lambda }{\left(\frac{t}{\lambda }\right)}^{\delta -1}exp\left[\theta D\left(t\right)\right]$$

(1)

where $D(t)$ is a stochastic variable that indicates the system degradation state and $t$ is the usage variable. Hazard rate function can be parametrized by specifying a form of the baseline hazard rate due to the adaptability of the Weibull distribution. The degradation of the system is assumed to follow gamma process. Then, the probability density distribution of $D(t)$ is

$$f_{\alpha ·t,\beta }\left(x\right)=\frac{1}{\mathrm{\Gamma }(\alpha ·t)}{\beta }^{\alpha ·t}{x}^{\alpha ·t-1}\exp (-\beta x)$$

(2)

where $\alpha ·t$ and $\beta$ are the shape function and scale parameter, respectively. As the aforementioned assumption, the degradation is operation-dependent. This means that $t$ in the distribution function refers to the operating time of the system, which depends on the production decision. The uncertainty of the degradation process is endogenous, where the parameter realizations of the probability distribution are influenced by the production decisions.

To avoid disruptions in production, the inspection is carried out at the beginning of every period. The hazard rate can be calculated by equation (1) using the usage time $t_k$ and degradation state $D(t_k)$ after inspection. The maintenance policy is to preventively replace the system when its hazard rate exceeds a predetermined threshold $Q$ . Under the failure limit policy, the proposed CBM policy can be defined with a binary variable $u_k$ as

$$u_k=\{\begin{array}{c}1,h(t_k,D\left(t_k\right))\ge Q\hfill \\ \multicolumn{1}{c}{0,h(t_k,D\left(t_k\right))<Q}\end{array};\forall k$$

(3)

Here, the variable $u_k$ equals 1 represents that a PM action is performed at the beginning of $k$ -th period, or 0 otherwise. After PM, the system is restored to the “as good as new” condition, which means that the degradation state $D(t_k)$ and the usage time $t_k$ are both reset to 0.

Although the CBM policy can decrease failure probability and increase the system reliability, random failures may happen in production since the hazard rate is not 0. CM is conducted as soon as a failure occurs to resume production without changing the system’s usage time and degradation state. After CM, the production continues till the batch is completed.

Considering the random failure numbers in every period, the total maintenance time in each period can be defined as

$$M{T}_k=T_p·u_k+T_c·F_k;\forall k$$

(4)

Considering the inspection cost, the total maintenance cost in each period can be written as

$$M{C}_k=c_I+c_p·u_k+c_c·F_k;\forall k$$

(5)

Production model

In this problem, a set of products $P$ need to be produced during a finite horizon of $T$ periods. The bold face ${\mathit{d}}_{i,k}$ is used to represent the random demands so that they can be distinguished from their particular realizations $d_{i,k}$ . The demand $d_{i,k}$ is to be satisfied at the end of period $k$ ( $k=1,2,\dots ,T$ ). Note that there may be no sufficient capacity in every period for the system to satisfy additional demands. So, improving the robustness of the plan through safety stock in every period is not realistic. A robust predictive plan which can smooth the unpredictable fluctuations through backorder and inventory is needed. In production model, the total production cost in each period is evaluated as

$$\begin{array}{c}P{C}_k=\hfill \\ \sum _{i=1}^P(c_{i,k}·x_{i,k}+h_{i,k}·I_{i,k}+s_{i,k}·K_{i,k}+b_{i,k}·B_{i,k});\forall k\hfill \end{array}$$

(6)

Other constraints of production model are shown in the integrated model in next sub-section.

Integrated model

Let $\mathrm{TC}$ be the total cost of the initial plan and $T{C}^r$ be the corresponding realized total cost. Since ${\mathit{d}}_{i,k}$ are random variables, so is the total cost of the realized plan. From the point of robustness, what really matters is the performance of the realized plan rather than the initial plan. Therefore, the expected realized total cost $E(T{C}^r)$ is used to be the robustness measure in this article. Considering stability, the realized plan should have minimal deviation cost (i.e. penalty cost caused by the changes in initial plan). Both random failures and uncertain demand could change the initial plan due to the changes in supply and demand. Thus, the deviation between the supply capability (i.e. inventory in previous period and production in current period) and the realized demand in every period is used to evaluate the variable cost. Considering the hedge of uncertainties over continuous periods, the cumulative production and demand variables are used to calculate. As a result, the expected sum of deviation cost $\sum _{k=1}^T\sum _{i=1}^P[{\vartheta }^{+}·{\mathrm{\Delta }}_{i,k}^{+}+{\vartheta }^{-}·{\mathrm{\Delta }}_{i,k}^{-}]$ is used to be the stability measure. The planned cumulative production quantity may deviate from the cumulative demand by amounts ${\mathrm{\Delta }}_{i,k}^{+}$ or ${\mathrm{\Delta }}_{i,k}^{-}$ . An additional storage cost ${\vartheta }^{+}$ per unit for the deviation ${\mathrm{\Delta }}_{i,k}^{+}$ and an outsourcing cost ${\vartheta }^{-}$ per unit for the deviation ${\mathrm{\Delta }}_{i,k}^{-}$ are needed to be paid. When considering robustness and stability at the same time, the sum of expected realized total cost and expected deviation cost is minimized. The joint model can be finally established as follows:

Minimize

$$Z=E\left(T{C}^r\right)+E\left(\xi \right)$$

(7)

Subject to equations (1)–(6)

$$T{C}^r=\sum _{k=1}^T(P{C}_k+M{C}_k)$$

(8)

$$\xi =\sum _{k=1}^T\sum _{i=1}^P[{\vartheta }^{+}·{\mathrm{\Delta }}_{i,k}^{+}+{\vartheta }^{-}·{\mathrm{\Delta }}_{i,k}^{-}];\forall i,k$$

(9)

$${\mathrm{\Delta }}_{i,k}^{+}=\max \{0,\sum _{v=1}^k{X}_{i,v}-\sum _{v=1}^k{d}_{i,v}\};\forall i,k$$

(10)

$${\mathrm{\Delta }}_{i,k}^{-}=\max \{0,\sum _{v=1}^k{d}_{i,v}-\sum _{v=1}^k{X}_{i,v}\};\forall i,k$$

(11)

$$X_{i,k}\le K_{i,k}·W;\forall i,k$$

(12)

$$\sum _{i=1}^P(r_{i,k}·X_{i,k})<G_k;\forall k$$

(13)

$$P{T}_k=\sum _{i=1}^P(r_{i,k}·x_{i,k});\forall k$$

(14)

$$P{T}_k+M{T}_k\le G_k;\forall k$$

(15)

$$x_{i,k}\le K_{i,k}·W;\forall i,k$$

(16)

$$x_{i,k}-I_{i,k}+I_{i,k-1}+B_{i,k}-B_{i,k-1}=d_{i,k};\forall i,k$$

(17)

$$u_k=\{\begin{array}{c}1,h(t_k,D\left(t_k\right))\ge Q\hfill \\ \multicolumn{1}{c}{0,h(t_k,D\left(t_k\right))<Q}\end{array};\forall k$$

(18)

$$t_{k+1}=(1-u_k)t_k+\sum _{i=1}^P(r_{i,k}·x_{i,k});\forall k$$

(19)

$$x_{i,k}=g(X_{i,k},M{T}_k);\forall i,k$$

(20)

$$X_{i,k},x_{i,k},I_{i,k},B_{i,k}\ge 0;K_{i,k},u_k\in \{0,1\};\forall i,k$$

(21)

$$B_{i,0}=0;I_{i,0}=0;\forall i$$

(22)

The objective of the joint model is to minimize the sum of expected total cost and expected deviation cost. The total cost is the sum of realized production cost and maintenance cost in every period in Constraint (8). The sum of deviation cost is calculated in Constraint (9). Constraints (10) and (11) are the formulas of cumulative deviation between the planned production quantity and realized demand. Constraints (12) and (16) are the setup constraints. The capacity constraints for the initial plan and the real plan are shown in Constraints (13) and (15), respectively. Constraint (14) is the formula of total production time in every period. Constraint (17) is the standard inventory balance equation of the real production plan. Constraint (18) represents the CBM policy. Constraint (19) is the iterative formula of the system usage time at the beginning of every period. Constraint (20) represents the mapping relationship between the initial plan and the realized plan based on the maintenance time and corresponding constraints, which is explained in detail in Algorithm 3. Constraints (21) and (22) represent non-negative, binary and initialization variables, respectively.

Solution algorithm

Preliminary analysis

The aim is to construct an initial production plan and a maintenance policy that minimize the cost (7) for a capacitated system with known demand distributions and degradation process. Obviously, the integrated problem of production and maintenance is very complicated even assuming the demand is deterministic. Although there are several algorithms in literature to solve the deterministic problem efficiently, it is difficult to find the optimal decisions to guarantee the robustness and stability of the plan at the same time.

Considering that demand variables are independent of each other, the stability of the plan is only decided by the planned production quantity. The predictive demands that minimize the deviation cost $\xi$ can be easily determined by using Algorithm 1.

Algorithm 1. Determine the optimal predictive demands
1	for $k=1,2,\dots ,T$ :
2	for $i=1,2,\dots ,P$ :
3	use the convolution of demand distribution to derive the probability distribution of $\sum _{v=1}^k{d}_{i,v}$
4	set $p_{i,k}=\frac{{\vartheta }^{-}}{{\vartheta }^{-}+{\vartheta }^{+}}$
5	determine the cumulative demand $d\in {\mathbb{R}}^{+}$ such that $P(\sum _{v=1}^k{d}_{i,v}\le d-1)<p_{i,k}\le P\left(\underset{v=1}{\overset{k}{\sum _{v=1}^k{d}_{i,v}\le d}}\right)$
6	if ( $k>1$ ) set ${\overline{d}}_{i,k}^{*}=d-\sum _{v=1}^{k-1}{\overline{d}}_{i,v}^{*}$
7	else ${\overline{d}}_{i,k}^{*}=d$
8	end
9	end
10	end

If the predictive demand can be perfectly satisfied by the planned production quantity, in other words, the predictive demand is used as deterministic demand to decide the initial production plan, then the minimal deviation cost is achieved. The value for ${\overline{d}}_{i,k}^{*}$ can be derived from the formula shown in step 5 of Algorithm 1, which uses $p_{i,k}$ that is calculated in step 4. The formula presented in step 5 is nothing else than the classical fractional formula encountered in the newsboy problem. ²⁶ And, ${\vartheta }^{-}$ and ${\vartheta }^{+}$ are equivalent to the cost coefficients that used in newsboy problem. Please note that the obtained predictive demand through Algorithm 1 only plays a role in minimizing the deviation cost; therefore, a search algorithm on the predictive demand to minimize the objective function (7) is proposed later.

After the predictive demand is determined, the only uncertainty in this problem when deriving an initial plan is the random failure. The hazard rate $h(t,D(t))$ depends on both usage time $t$ and the system degradation level $D(t)$ . The system failure can be described as a model of a stochastic point process. Since $D(t)$ is a positive and increasing random variable, the process of the failure is a doubly stochastic Poisson process with a stochastic intensity depending on the system degradation level^27,28. Hence the probability of the failure numbers in period $k$ is given by

$$\mathrm{Pr}\{F_k=n\}=E\left[\frac{1}{n!}{\left[\mathrm{\Lambda }\left(t_k\right)\right]}^n\exp (-\mathrm{\Lambda }\left(t_k\right))\right]$$

(23)

where

$$\mathrm{\Lambda }\left(t_k\right)=\underset{t_k}{\overset{t_{k+1}}{\int }}h(\tau ,D\left(\tau \right))d\tau$$

(24)

The expected failure numbers $E(F_k)$ can be written as

$$\begin{array}{c}\multicolumn{1}{c}{E\left(F_k\right)=\sum _{n=0}^{\infty }\{n·E\left[\frac{1}{n!}{\left[\mathrm{\Lambda }\left(t_k\right)\right]}^n\exp (-\mathrm{\Lambda }\left(t_k\right))\right]\}}\\ =E\left(\mathrm{\Lambda }\left(t_k\right)\right)=E_{\left\{D\left(\tau \right)\right\}}\left\{\underset{t_k}{\overset{t_{k+1}}{\int }}h(\tau ,D\left(\tau \right))d\tau \right\}\hfill \\ =\underset{t_k}{\overset{t_{k+1}}{\int }}{E}_{\left\{D\left(\tau \right)\right\}}\left\{h(\tau ,D\left(\tau \right))\right\}d\tau \hfill \end{array}$$

(25)

Since the degradation level $D\left(\tau \right)$ follows a gamma process, the expected value of hazard rate based on the initial degradation level $D\left(t_k\right)$ at time $t_k$ is given as

$$\begin{array}{c}{E}_{\left\{D\left(\tau \right)\right\}}\left\{h(\tau ,D\left(\tau \right))\right\}=E_{\left\{D\left(\tau \right)\right\}}\left\{\frac{\delta }{\lambda }{\left(\frac{\tau }{\lambda }\right)}^{\delta -1}\exp \left[\theta D\left(\tau \right)\right]\right\}\hfill \\ =\underset{0}{\overset{\infty }{\int }}\frac{\delta }{\lambda }{\left(\frac{\tau }{\lambda }\right)}^{\delta -1}\exp \left[\theta (D\left(t_k\right)+x)\right]f\left(x\right)\mathrm{dx}\hfill \end{array}$$

(26)

where $f(x)$ is the probability density function of the increment of degradation level $x$ . Then

$$\begin{array}{c}{E}_{\left\{D\left(\tau \right)\right\}}\left\{h(\tau ,D\left(\tau \right))\right\}\hfill \\ =\frac{\delta }{\lambda }{\left(\frac{\tau }{\lambda }\right)}^{\delta -1}\exp \left(\theta D\left(t_k\right)\right)\underset{0}{\overset{\infty }{\int }}\exp \left(\theta x\right)f\left(x\right)\mathrm{dx}\hfill \\ =\frac{\delta }{\lambda }{\left(\frac{\tau }{\lambda }\right)}^{\delta -1}\exp \left(\theta D\left(t_k\right)\right)\hfill \\ \underset{0}{\overset{\infty }{\int }}\frac{{\beta }^{\alpha \tau }{x}^{\alpha \tau -1}\exp \left[(\theta -\beta )x\right]}{\mathrm{\Gamma }\left(\alpha \tau \right)}\mathrm{dx}\hfill \\ =\exp \left(\theta D\left(t_k\right)\right)\frac{\delta }{\lambda }{\left(\frac{\tau }{\lambda }\right)}^{\delta -1}\frac{{\beta }^{\alpha \tau }}{{(\beta -\theta )}^{\alpha \tau }}\hfill \\ \underset{0}{\overset{\infty }{\int }}\frac{{(\beta -\theta )}^{\alpha \tau }{x}^{\alpha \tau -1}\exp [-(\beta -\theta )x]}{\mathrm{\Gamma }\left(\alpha \tau \right)}\mathrm{dx}\hfill \\ =\exp \left(\theta D\left(t_k\right)\right)\frac{\delta }{\lambda }{\left(\frac{\tau }{\lambda }\right)}^{\delta -1}\frac{{\beta }^{\alpha \tau }}{{(\beta -\theta )}^{\alpha \tau }}\hfill \end{array}$$

(27)

Based on the equations (25) and (27), the expected failure numbers, which are used to evaluate the maintenance cost and time to derive an initial plan on the basis of deterministic demand and known maintenance threshold, can be calculated. And, the system reliability function is written as

$$R\left(t\right)=E_{\{D\left(\tau \right),0\le \tau \le t\}}\left[\exp (-\underset{0}{\overset{t}{\int }}h(\tau ,D\left(\tau \right))d\tau )\right]$$

(28)

According to the property that has been proved by Yashin, ²⁹ the failure cumulative distribution function becomes

$$\begin{array}{c}\multicolumn{1}{c}{H\left(t\right)=1-R\left(t\right)}\\ =1-\exp (-\underset{0}{\overset{t}{\int }}{E}_{\left\{D\left(\tau \right)\right\}}\left\{h(\tau ,D\left(\tau \right))\right\}d\tau )\hfill \end{array}$$

(29)

Given the initial usage $t_0$ and state $D(t_0)$ , the cumulative distribution function is as follows

$$\begin{array}{c}H[t|(t_0,D\left(t_0\right))]=\hfill \\ 1-\exp \{-\exp \left(\theta D\left(t_0\right)\right)\frac{\delta }{\lambda }\underset{t_0}{\overset{t_0+t}{\int }}{\left(\frac{\tau }{\lambda }\right)}^{\delta -1}\frac{{\beta }^{\alpha \tau }}{{(\beta -\theta )}^{\alpha \tau }}d\tau \}\hfill \end{array}$$

(30)

Based on equation (30), the inversion sampling is used to sample the failure time to calculate the failure numbers in the simulation.

MPGA

As the formulated model is a complicated stochastic mixed-integer programming and involves endogenous uncertainty, it is impossible to be solved effectively by commercial solvers. Thus, a MPGA is developed in this article, which has the effectiveness when solving the combinatorial optimization problems.^30–32 The flowchart of the proposed MPGA is shown in Figure 1.

Figure 1.

Flowchart of the proposed MPGA.

Coding strategy

The computational complexity of the integrated model forces us to use simulation instead of an analytical method to evaluate the objective function (7). Before performing the simulation, an initial plan is needed to be a baseline, in which the planned production quantity and maintenance threshold are the needed decision variables. Based on the previous analysis, the initial plan, with the objective to minimize the expected total cost $\mathrm{TC}$ , can be derived from a deterministic integrated model by a heuristic algorithm. So, the decision variables in MPGA become the predictive demands and the maintenance threshold. Since the predictive demands obtained by Algorithm 1 can only minimize the deviation cost, a neighborhood search on ${\overline{d}}_{i,k}$ in MPGA is performed to accelerate the search process while ensuring its quality. The real number code is used to represent the predictive demands within limited range. The maintenance threshold is represented by a 12-bit binary code. Figure 2 shows an example of chromosome. The binary number “000000000000” and “111111111111” are used to represent 0 and 1, for example, the threshold “000000011001” can be calculated as $\left(000000011001\right)/\left(111111111111\right)=0.0063$ . Then, the two sub-chromosomes are combined to form a chromosome.

Figure 2.

An example of the chromosome representation.

Population initialization

In MPGA, three different subpopulations are initialized with different methods. Each subpopulation contains $Y(y=1,\dots ,Y)$ chromosomes. The real coding part of the chromosomes in first subpopulation are randomly generated in the range of $[0.7*{\overline{d}}_{i,k},{\overline{d}}_{i,k}]$ ; in second subpopulation, the real coding part are randomly generated in the range of $[0.9*{\overline{d}}_{i,k},1.1*{\overline{d}}_{i,k}]$ ; and in third subpopulation, the real coding part are randomly generated in the range of $[{\overline{d}}_{i,k},1.3*{\overline{d}}_{i,k}]$ . The binary coding part in all subpopulations is randomly generated.

Selection, crossover and mutation

The roulette wheel selection method is adopted to pick up the parent chromosomes to generate the next generation. The selecting probability of a chromosome $y$ is proportional to its fitness value ${\mathrm{\Theta }}_y$ , which is the reciprocal of the objective function. Then, a specially designed genetic operator is proposed for the proposed chromosome representation. First, a crossover operator takes action to execute the gene recombination. Specifically, different sub-chromosomes use different crossover operators: the real coding sub-chromosome adopts uniform crossover and the binary coding sub-chromosome adopts single-point crossover. A random crossover rate $p^c$ is adopted in crossover operator. Second, two mutation operators including simple mutation and inversion mutation are applied in the two sub-chromosomes. The adaptive mutation rate $p_y^m$ of chromosome $y$ is adopted as follows

$$p_y^m=p^m*(1-\frac{{\mathrm{\Theta }}_y}{\sum _y{\mathrm{\Theta }}_y})$$

(30)

in which $p^m$ represents the constant mutation rate. Since all genes generated by these operators are within their domain, there is no need for repair operations.

Upgrading and immigration

The best chromosomes from all subpopulations are recorded as the elite population, which is updated for every generation to improve the fitness value. This is the upgrading process. To make different subpopulations co-evolutionary, an immigration strategy is proposed to keep chromosome information exchange between individuals from different subpopulations. For each subpopulation, the best chromosome is chosen to substitute the worst one in the next subpopulation.

Stopping criteria

If the best fitness value in the elite population does not change for a given number of consecutive generations, or the total number of generations reaches a given upper bound, the algorithm terminates, and the latest best fitness value and its corresponding chromosome are considered as the outputs of the MPGA.

Decoding algorithm

Due to the minimum deviation cost in objective function that can be obtained on the premise of meeting the predictive demands, the main purpose in the decoding algorithm is to use the predictive demands and fixed maintenance threshold to generate an initial plan $\Pi (X_{i,k},Q)$ with the objective of minimizing the expected total cost, and then to calculate objective function (7) based on the initial plan through simulation. The general steps for generating $\Pi (X_{i,k},Q)$ are outlined in Algorithm 2.

Algorithm 2. Generate initial plan $\Pi (X_{i,k},Q)$
Input: predictive demands ${\overline{d}}_{i,k}$ and fixed maintenance threshold $Q$
Output: initial plan $\Pi (X_{i,k},Q)$
1	set $l=1$ and substitute the demand $d_{i,k}$ with ${\overline{d}}_{i,k}$ in equation (17)
2	relax capacity constraint (15) by substituting the maintenance time $\{{\mathrm{MT}}_k^0\}$ with $\{0\}$ and set the objective as the production cost (6) to obtain a CLSP model under Constraints (14)-(17), (21)-(22)
3	while $(l<L)$
4	solve the CLSP model by CPLEX to obtain a trial production plan $\{x_{i,k}^l\}$
5	calculate the expected maintenance time $\{{\mathrm{MT}}_k^l\}$ based on the known $Q$ and $\{x_{i,k}^l\}$ through  equations (4), (5), (19), (25) and (27)
6	calculate $\mathrm{Gap}\left(l\right)=\frac{\sum _k({\mathrm{MT}}_k^l-{\mathrm{MT}}_k^{l-1})}{\sum _k{\mathrm{MT}}_k^l}$
7	if ( $\mathrm{Gap}\left(l\right)<\epsilon$ ) break and set $Q=Q$ , $X_{i,k}=x_{i,k}^l$
8	else update the CLSP model by substituting the maintenance time with $\{{\mathrm{MT}}_k^l\}$
9	end
10	$l=l+1$
11	end

First, the original model is transformed to be an integrated production and maintenance model (M₁) under deterministic demand and random failures. The difficulty is that the expected maintenance time depends on the production time in every period, resulting in a chicken-and-egg conundrum between production quantities and maintenance time. With the objective to minimize the expected total cost, M₁ can be solved by using the iterative algorithm in Zhao et al. ⁴ to generate an initial plan. The procedure is described in Algorithm 2, in which the simplified model is solved by CPLEX. By this algorithm, a near-optimal plan can be found without constraint violation. The convergence of Algorithm 2 is illustrated in Zhao et al. ⁴ After determining the initial plan $\Pi (X_{i,k},Q)$ , the objective function is calculated through the simulation described in Algorithm 3.

Algorithm 3. Calculate the objective function through simulation
Input: initial plan $\Pi (X_{i,k},Q)$
Output: the value of objective function (7)
1	sample the realized demand $d_{i,k}$
2	for $k=1,2,\dots ,T$ :
3	inspect the system and determine the usage $t_k$ and state $D(t_k)$
4	calculate the hazard rate $h(t_k,D(t_k))$ and determine $u_k$ through equation (3)
5	set $j=1$ and $t_0=t_k,D(t_0)=D(t_k)$ in equation (30)
6	sample the failure time ${\tau }_j$ based on equation (30)
7	if $(\sum _j{\tau }_j<\sum _{i=1}^P(r_{i,k}·X_{i,k}))$ set $j=j+1$ and update $t_0$ with ${\tau }_j$ in equation (30), then go to step 6.
8	else set $F_k=j-1$
9	end
10	calculate the realized maintenance cost $M{C}_k$ and time $M{T}_k$ in equations (4) and (5)
11	if $(M{T}_k+\sum _{i=1}^P(r_{i,k}·X_{i,k})\le G_k)$ set $x_{i,k}=X_{i,k}$
12	else initial plan is uncompleted, set $\mathrm{\Delta }G=M{T}_k+\sum _{i=1}^P(r_{i,k}·X_{i,k})-G_k$ , calculate the realized production plan by $x_{i,k}=X_{i,k},i\ne P;x_{P,k}=X_{P,k}-\lceil \frac{\mathrm{\Delta }G}{r_{P,k}}\rceil$
13	end
14	calculate the realized inventory and backorder quantities based on $x_{i,k}$ in equation (17)
15	calculate the realized production cost $P{C}_k$ based on equation (6)
16	end
17	calculate the objective function (7)

Numerical results and discussion

Parameters selection

In this article, all experiments have been performed on an Intel Core i7-3.6 GHz with 12 GB RAM computer, with algorithms compiled by MATLAB. The investigated degradation system has the characteristics given in Table 1. The planning horizon is composed of $T=8$ periods. And, there are $P=3$ kinds of products needed to be produced in lots to satisfy the demands in every period. The distribution functions of the three products’ demands are uniform, normal and normal distributions. The additional storage cost ${\vartheta }^{+}$ and outsourcing cost ${\vartheta }^{-}$ are 2¥/unit/period and 8¥/unit/period, respectively. The related parameters in production are presented in Table 2. In MPGA, set $p^c\in [0.85,0.95]$ , $p^m=0.2$ . The system capacity is determined by the total needed processing time over the horizon and average capacity utilization $\rho$ with the below formula

$$G_k=\frac{\sum _k\sum _i{r}_{i,k}{\overline{d}}_{i,k}}{\rho T}$$

(30)

Table 1.

Parameters of the degradation system.

Parameters	Values
$\alpha$	0.14
$\beta$	1
$\delta$	2
$\lambda$	500
$\theta$	0.01
$c_I$	5¥
$c_p$	200¥
$c_c$	100¥
$T_p$	20 h
$T_c$	10 h

Table 2.

Parameters of the three products in production.

Product	1	2	3
$d_{i,k}$ (units)	U(40,60)	N(20,5)	N(60,10)
$r_{i,k}$ (h/unit)	2	3	1
$s_{i,k}$ (¥/unit)	5	8	6
$c_{i,k}$ (¥/unit)	1	2	3
$h_{i,k}$ (¥/unit)	3	3	3
$b_{i,k}$ (¥/unit)	10	10	10

Since sampling and simulation are embedded in the procedure of MPGA, the computational time will depend mainly on the sample size $S$ . To obtain high-quality solution within a reasonable time, the sample size is a key parameter that requires appropriate selection. On the one hand, too few samples may cause the inaccuracy of the statistical value of the simulation; on the other hand, too many samples will bring huge computational time. In this article, the confidence interval approach is adopted to analyze the minimum sample size to guarantee a certain level of confidence for a given deviation. The minimum sample size $S$ can be estimated using the following equation

$$S=\min \{e:\frac{t_{e-1,1-\alpha /2}\sqrt{s^2\left(e\right)/e}}{\left|\overline{Z}\left(e\right)\right|}\le ϵ\}$$

(30)

where $e$ is the number of replications from the start, $t_{e-1,1-\alpha /2}$ is the $(1-\alpha )$ percentile of the t-student distribution with $e-1$ degrees of freedom, $s^2(e)$ is the sample variance from $e$ replications, $\overline{Z}(e)$ is the mean value of the objective function for the current replication and $ϵ$ is the specified relative deviation.

Three test instances with the parameters given in Tables 1 and 2 are generated to determine the sample size. For each instance with different decision variables, let ${\overline{Z}}_q$ and ${\phi }_q$ denote, respectively, the cumulative mean value and the relative deviation of the $q$ -th instance ( $q=1,2,3$ ). The results of the instances under 95% confidence level and $ϵ=4\mathrm{\%}$ relative deviation are shown in Figure 3. Clearly, the mean value tends to be smooth and the deviation tends to decrease as the sample size increases for each instance. And, all instances have a lower relative deviation than 4% when the sample size exceeds 100, which indicates that the minimum sample size is 100.

Figure 3.

The cumulative mean value and relative deviation of the three instances.

Results and sensitivity analysis

In this proposed example, the capacity utilization $\rho$ is set to 0.9 and the maximum generations of program is set to 50, and the program will stop if the best fitness value is not improved by 10 consecutive generations. The minimum and maximum objective values of each generation are plotted in Figure 4. The convergence of the proposed MPGA is evident. The optimal maintenance threshold $Q$ is 0.0058 and the value of the objective function is 6009.34¥, including the expected realized total cost of 4643.86¥ and the expected deviation cost of 1365.48¥.

Figure 4.

The evolution process and convergence of MPGA.

When considering the influence of system parameters on the expected costs, the main concerns are the parameters of demand variability, degradation rate and capacity tightness. The value of $\sigma ({\mathit{d}}_{i,k})$ that equals to the standard deviation of the demand distribution is used to represent the demand variability. The value is changed by multiplying a factor $\omega$ to test different demand variabilities. The factor $\omega$ is set to 0.5, 0.75, 1, 1.25 and 1.5. Then, five scenarios are generated and the corresponding standard deviations of the demand distribution are $0.5*\sigma ({\mathit{d}}_{i,k})$ , $0.75*\sigma ({\mathit{d}}_{i,k})$ , $\sigma ({\mathit{d}}_{i,k})$ , $1.25*\sigma ({\mathit{d}}_{i,k})$ and $1.5*\sigma ({\mathit{d}}_{i,k})$ , respectively. The degradation rate is easy defined by the scale parameter $\lambda$ in the Weibull hazard rate function of the PHM. Set $\lambda =300,400,500,600,700$ to generate five new different scenarios under original demand variability to test the impact of degradation rate only. The capacity tightness is the last factor tested in the experiment. The capacity utilization factor $\rho$ is used to define the capacity tightness with different values of 0.85, 0.875, 0.9, 0.925 and 0.95. Correspondingly, five new different scenarios are generated under the original demand variability and degradation rate that are presented in “Parameters selection.”

Fifteen instances with different parameters are realized with the developed MPGA. The expected realized total cost, the expected deviation cost and the objective value are compared under different factors. In Figure 5, there is a visible trend that all costs are increasing with the increase of factor $\omega$ when other parameters are same. The reason behind it is that the demand variability has a direct effect on the deviation cost, and higher $\omega$ means more fluctuations in production which may lead to more production cost. This result shows that demand variability has a great impact on the robustness and stability of the solution, and further indicates that the proposed model may not be applicable for the situation where demand variability is huge. The impact of the degradation rate based on the parameter $\lambda$ is shown in Figure 6, which shows that the expected realized total cost and the objective function are decreasing with the increase of parameter $\lambda$ , while the expected deviation cost has no obvious change. It can be explained by the fact that $\lambda$ is the scale parameter in the Weibull hazard function, and bigger $\lambda$ means less failures in the same time interval, which may decrease the maintenance cost and the impact of random failures on production capacity. In the same way, Figure 7 shows that all costs have a slow increasing trend with the increase of capacity utilization factor $\rho$ . The reason is that the capacity tightness will affect the initial plan, and tighter capacity may lead to more unsatisfied predictive demand, resulting in increased production cost and deviation cost. In summary, when only considering robustness, the three factors, especially the demand variability, have a significant impact on the solutions, while the degradation rate is less sensitive than the other factors when only considering stability.

Figure 5.

The impact of different demand variability factor $\omega$ .

Figure 6.

The impact of different deterioration rate factor $\lambda$ .

Figure 7.

The impact of different capacity tightness factor $\rho$ .

Performance analysis on the algorithm

The performance of the proposed MPGA is compared with that of a standard genetic algorithm (GA), and a simulated annealing (SA) ³³ on several instances. According to the results from the previous experiments, it is confirmed that the robustness and stability of the solution are mainly affected by demand variability and capacity tightness, so the comparison tests are conducted under different demand variabilities and capacity tightness. Three different levels of $\omega$ and $\rho$ are chosen, that is, $\omega =0.5,1,1.5$ and $\rho =0.85,0.9,0.95.$ Therefore, a total of nine (= 3²) tests are generated. For each test, these three algorithms are implemented using the same coding and decoding strategies as described in “MPGA.” To control the computational time, the stopping criteria of these algorithms are set to a maximum running time of 1800 s.

In Table 3, the expected realized total cost, the expected deviation cost and the objective function of the three algorithms are calculated. From the results, MPGA is always better than the other algorithms in all comparative tests, and GA is better than SA, which indicates that the global search ability in GA has a great influence on the solution of this problem. And, these results validate the good performance of the proposed MPGA in solving this model. With the increase of factor $\omega$ , all tests in comparison with all algorithms show an increasing trend, and the same trend can be found when considering the increase of factor $\rho$ . Furthermore, to compare the relative improvement of MPGA over other algorithms, the gap ratio, calculated as $\mathrm{gap}\left(\mathrm{GA}\right)=\mathrm{GA}-\mathrm{MPGA}/\mathrm{GA}-\mathrm{MPGA}\mathrm{MPGA}\mathrm{MPGA}$ and $\mathrm{gap}\left(\mathrm{SA}\right)=$ $\mathrm{SA}-\mathrm{MPGA}/\mathrm{SA}-\mathrm{MPGA}\mathrm{MPGA}\mathrm{MPGA}$ , of all comparative tests is given in Table 4. The results show that MPGA obtains better robustness than other algorithms in all tests, and obtains better stability except only two tests compared to GA. Nevertheless, the average improvement of stability is higher than robustness in both comparisons, which means that MPGA is better at optimizing the stability. Comparing the improvement under different factors, there is no visible trend in the results, which can be explained by the near-optimal solutions. To sum up, the developed MPGA has a good performance for this model considering that the average improvement in objective function is 2.51% and 10.40% over GA and SA, respectively.

Table 3.

Results of the three algorithms under different $\omega$ and $\rho$ .

$\rho$	$\omega$	MPGA			GA			SA
		$E(T{C}^r)$ (¥)	$E(\xi )$ (¥)	$Z$ (¥)	$E(T{C}^r)$ (¥)	$E(\xi )$ (¥)	$Z$ (¥)	$E(T{C}^r)$ (¥)	$E(\xi )$ (¥)	$Z$ (¥)
0.85	0.5	3867.30	702.30	4569.59	3903.27	721.55	4624.82	4196.20	897.20	5093.40
	1	4588.98	1331.65	5920.63	4650.39	1391.39	6041.77	4899.05	1393.05	6292.10
	1.5	5339.12	2038.62	7377.73	5651.31	2174.81	7826.11	5984.40	2283.40	8267.80
0.9	0.5	3904.09	719.59	4623.68	4017.13	767.13	4784.26	4261.40	891.40	5152.80
	1	4643.86	1365.48	6009.34	4712.44	1406.44	6118.88	5020.95	1469.95	6490.90
	1.5	5521.02	2100.52	7621.54	5575.62	2140.62	7716.24	6104.80	2323.80	8428.60
0.95	0.5	3923.11	763.11	4686.21	4122.19	746.64	4868.83	4284.20	898.20	5182.40
	1	4777.09	1399.18	6176.27	4817.57	1429.45	6247.02	5213.80	1533.80	6747.60
	1.5	5466.52	2212.52	7679.04	5619.36	2189.36	7808.72	6255.33	2491.33	8746.66

MPGA: multi-population genetic algorithm; GA: genetic algorithm; SA: simulated annealing.

Table 4.

Improvement of MPGA over other algorithms in cost saving under different $\omega$ and $\rho$ .

$\rho$	$\omega$	$\mathrm{gap}(\mathrm{GA})$			$\mathrm{gap}(\mathrm{SA})$
		$E(T{C}^r)$	$E(\xi )$	$Z$	$E(T{C}^r)$	$E(\xi )$	$Z$
0.85	0.5	0.93%	2.74%	1.21%	7.50%	24.34%	11.46%
	1	1.34%	4.49%	2.05%	5.35%	0.12%	6.27%
	1.5	5.85%	6.68%	6.08%	5.89%	4.99%	12.06%
0.9	0.5	2.90%	6.61%	3.47%	6.08%	16.20%	11.44%
	1	1.48%	3.00%	1.82%	6.55%	4.52%	8.01%
	1.5	0.99%	1.91%	1.24%	9.49%	8.56%	10.59%
0.95	0.5	5.07%	-2.16%	3.90%	3.93%	20.30%	10.59%
	1	0.85%	2.16%	1.15%	8.22%	7.30%	9.25%
	1.5	2.80%	-1.05%	1.69%	11.32%	13.79%	13.90%
Average		2.47%	2.71%	2.51%	7.15%	11.12%	10.40%

Performance analysis on the model

This part of experiment is the numerical investigation on the proposed integrated model. It has two parts: the first part concerns the performance of the model in uncertain environment, and the second part compares the performance of the model with that of a separate decision model.

To investigate the performance under uncertainty, three different methods are compared. The first one is the proposed proactive approach (M₁) with the aim to obtain a robust and stable plan under uncertain environment. According to the mean values of demands, the second method (M₂) solves a deterministic model to obtain a baseline plan and maintenance threshold. And, the solutions are put into simulation to calculate the expected real costs and deviation costs. The third method (M₃) is a posterior method which provides a lower bound for the original problem. Decisions are determined based on the real demand and failure information in a particular instance. To compare with the previous two methods, each test solves 100 instances with different perfect information and calculates the expected real costs (no deviation costs) of these instances. The corresponding total costs of these three methods are denoted as Z₁, Z₂ and Z₃. In fact, the gap between Z₁ and Z₂ evaluates the value of stochastic solution, and the gap between Z₁ and Z₃ evaluates the value of perfect information. ³⁴

Table 5 illustrates the experiment results for nine different tests under different factors of $\omega$ and $\rho$ . Z₁ is always smaller than Z₂ in all tests, which validates that the proposed model can generate a better plan at handling the uncertainties than that obtained by deterministic method. The gap between Z₁ and Z₂ is positive and it increases significantly when $\omega$ increases to 1.5. Although the gap shows an increasing trend with the increase of $\rho$ , factor $\omega$ plays a more important role than $\rho$ in this comparison. These results reflect that the uncertainties in decision-making should not be ignored. The gap between Z₁ and Z₃ is negative and it shows an increasing trend with the increase of $\omega$ . In fact, Z₃ does not change much under different factors. The increased gap reflects that the impact of demand variability is more significant than capacity tightness. According to the analysis, the proposed proactive approach has better performance in dealing with small uncertainties, and it has more advantages than the deterministic method in dealing with large uncertainties.

Table 5.

Comparison results with different decision models under different $\omega$ and $\rho$ .

$\rho$	$\omega$	Z1(¥)	Z2(¥)	Z3(¥)	Z2–Z1	(Z2–Z1)/Z1	Z3–Z1	(Z3–Z1)/Z1
0.85	0.5	4569.59	7116.39	3179.33	2546.80	0.5573	–1390.26	–0.3042
	1	5920.63	8134.14	3181.04	2213.51	0.3739	–2739.59	–0.4627
	1.5	7377.73	12958.78	3262.79	5581.05	0.7565	–4114.94	–0.5578
0.9	0.5	4623.68	7252.10	3201.08	2628.42	0.5685	–1422.60	–0.3077
	1	6009.34	8552.74	3209.53	2543.40	0.4232	–2799.81	–0.4659
	1.5	7621.54	13816.52	3235.00	6194.98	0.8128	–4386.54	–0.5755
0.95	0.5	4686.21	7280.64	3251.71	2644.43	0.5704	–1384.50	–0.2986
	1	6176.27	9164.28	3263.67	2988.01	0.4838	–2912.60	–0.4716
	1.5	7679.04	16514.32	3317.92	8835.28	1.1506	–4361.12	–0.5679

The second part of the experiment compares the performance of the proposed integration model with a separate decision (M₄) method, in which the manufacturer decides the production plan based on a predetermined maintenance threshold. The maintenance threshold is individually optimized to minimize the average long-run maintenance cost with the objective function $G\left(T\right)=c_p+c_c·{\int }_0^{T}h(t,D\left(t\right))\mathrm{dt}/c_p+c_c·{\int }_0^{T}h(t,D\left(t\right))\mathrm{dt}TT,$ where $h(t,D(t))$ is calculated in equation (27). Taking the derivative of $G(T)$ to get the optimal time $T^{*}$ , the optimal hazard rate threshold is $h(T^{*})=0.0072$ . Based on the maintenance threshold, the production plan is optimized by the developed MPGA under the robustness and stability objective. The comparison results under different factors are shown in Table 6.

Table 6.

Comparison results between integration model and separate decision method under different $\omega$ and $\rho$ .

$\rho$	$\omega$	$E(T{C}^r)$ (¥)			$E(\xi )$ (¥)			$Z$ (¥)
		M 1	M 4	GAP	M 1	M 4	GAP	M 1	M 4	GAP
0.85	0.5	3867.30	4141.74	7.10%	702.30	832.74	18.57%	4569.59	4974.48	8.86%
	1	4588.98	4722.09	2.90%	1331.65	1577.09	18.43%	5920.63	6299.18	6.39%
	1.5	5339.12	5536.42	3.70%	2038.62	2423.42	18.88%	7377.73	7959.84	7.89%
0.9	0.5	3904.09	4261.46	9.15%	719.59	883.83	22.82%	4623.68	5145.29	11.28%
	1	4643.86	4820.00	3.79%	1365.48	1620.00	18.64%	6009.34	6440.00	7.17%
	1.5	5521.02	5630.75	1.99%	2100.52	2522.75	20.10%	7621.54	8153.50	6.98%
0.95	0.5	3923.11	4370.41	11.40%	763.11	1008.36	32.14%	4686.21	5378.77	14.78%
	1	4777.09	5097.45	6.71%	1399.18	1711.78	22.34%	6176.27	6809.23	10.25%
	1.5	5466.52	5741.19	5.02%	2212.52	2656.61	20.07%	7679.04	8397.80	9.36%

From Table 6, it can be found that both robustness and stability of the integration model are better than the separate decision method. When comparing the relative improvement of the joint model that is calculated as $\mathrm{GAP}=M_4-M_1/M_4-M_1{M}_1{M}_1$ , it can be seen that the improvements in stability are higher than the improvements in robustness in all tests although both improvements are remarkable. When $\omega =0.5$ and $\rho =0.95$ , the highest improvement of 14.78% is obtained in all tests. The gap of the objective function shows a decreasing trend with the increase of $\omega$ and an increasing trend with the increase of $\rho$ . In total, from the overview of Table 6, it is concluded that the integration model performs better under small demand variability and high capacity tightness situations.

Conclusion

In this article, a stochastic optimization framework is proposed to deal with the integration problem of batch production and CBM with uncertain demand and random failures. The uncertainty of the input parameters and the information about the reactive policy followed at execution time are considered proactively at the decision stage. The integration model needs to determine the production quantity for each product and maintenance policy simultaneously to form the complete initial plan. The deterministic CLSP is already non-deterministic polynomial-time (NP)-hard, thus the integration problem with CBM under uncertain environment is even more complex. A simulation-based MPGA is developed to optimize the problem for the purpose of high robustness and stability. The algorithm is implemented and tested with other algorithms under different experiments and has been proved efficient in getting good initial plans. The economic benefits in both robustness and stability under uncertain environment validate the necessity to consider the uncertainties in decision-making. This methodology is suitable for production and maintenance decisions in manufacturing system with multi-product, multi-period uncertain demand, which is common in semiconductor, automotive and electronic industries.

Other interesting directions for further research can be followed in integration problem of production and maintenance. For example, if the reactive policy is changed, for example, re-planning policy, the decision procedure will be more complex and the performance may be better.