# Joint Optimization of Production and Maintenance Using Monte Carlo Method and Metaheuristic Algorithms

Authors: Xiao-Zhi Ma; Wen-Yuan Lv
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

In the competitive business environment, manufacturers are seeking strategies to improve the product quality and the system reliability while reducing the costs. This paper addresses the problem of finding the optimal production and maintenance schedules for a deteriorating manufacturing system with the objective of minimizing the expected cost per unit time. The system consists of one machine which deteriorates with time and it may shift from an in-control state to an out-of-control state with a larger proportion of imperfect products. In addition, the hedging point policy is applied as the production-inventory control policy. The predictive maintenance is performed based on process inspections, whose sampling intervals are variable. To deal with the proposed problem, we build a joint model that coordinates production, inventory, maintenance, and quality control with 16 scenarios. Then we propose a novel approach speeding up the Monte Carlo simulation to calculate the objective function. Thus it becomes feasible to optimize the objective function by metaheuristic algorithms. Then we use the genetic algorithm to illustrate its feasibility. Next, the advantage of the proposed approach is verified by comparing with the traditional integral method. Finally, a sensitivity analysis with an orthogonal experiment is conducted to help managers find the factors with the most significant effect on the cost.


---

## Body

## 1. Introduction

In the competitive business environment, manufacturers are seeking strategies to improve the product quality and the system reliability while reducing the costs. Performance of a system is highly related to the coordination and cooperation of the subsystems. Production-inventory systems (PIS), maintenance systems (MS), and quality systems (QS) are three main subsystems of manufacturing systems with different goals.For the PIS, an important issue is how to control the production rate to keep suitable inventory level. In order to solve this issue, researchers have proposed two famous models: the economic production quantity (EPQ) [1] and the hedging point policy (HPP) [2]. The HPP is more suitable for systems subject to random failures and repairs. Because the HPP entails the build-up and preservation of a final product safety stock while the machine is operational in order to hedge against future shortages caused by machine failures, the optimality of the HPP has been demonstrated for failure and repair times described by homogeneous Markov processes in the case of constant demand rate [3] and stochastic demand rate [4]. For general failure and repair time distributions, the optimal control policy cannot be solved analytically [5, 6], but the structure of the optimal inventory policy can be approximated by the HPP [7–9]. Later, Rezg, Dellagi, and Chelbi [10] added one decision variable (when to start building up a buffer stock) to the traditional HPP, which can further decrease the total cost. Rather than HPP or EPQ, Souheil, Dellagi, and Rezg [11] used subcontracting to solve the issue and to satisfy a required service level.For the MS, there are enormous maintenance policies in the literature for various situations. The most common three types are corrective maintenance (CM), preventive maintenance (PM), and predictive maintenance (PRDM). CM is performed after failure to make the equipment continue to run. It is often a minimal repair [12], which just makes the equipment as good as old. PM is performed before failure to prevent the loss caused by sudden failure. It is usually performed based on time or age of the equipment [13]. PRDM is maintenance carried out when necessary, based on process inspections, diagnostic tests, or other means of condition monitoring [14–16]. The most widely used tool of process inspections is the X-bar control chart, whose sampling interval can be fixed or variable. In addition, Zied, Dellagi, and Rezg [17] showed the effect of the production rate variation on the optimal maintenance strategy. Later, Lv and Liu [18] proposed a maintenance policy for multicomponent systems with economic dependence and structural dependence.For the QS, manufacturers can either rework the imperfect products before sale with additional rework cost [19, 20] or sell the whole production quantity with warranty, where additional warranty cost is incurred [21]. In this paper, the manufacturer chooses the latter policy.Despite the conventional researches dealing separately with these subsystems, the integrated strategy is an appropriate approach to handle the interactions between them. For the integration of the three systems, relevant papers are summarized in Table1. In this table, “Monte Carlo” means that the Monte Carlo method [22] is used to calculate the objective function in a paper. “Metaheuristic method” means that the metaheuristic method [23], which includes genetic algorithms, particle swarm algorithms, etc., is used to optimize the objective function in a paper. “Orthogonal experiment” means that the orthogonal experiment is used to conduct the sensitivity analysis in a paper. It is noticeable that Mifdal, Hajej, and Dellagi [24] considered the integration of the PIS and the MS with a machine producing several products and satisfying some random demands, which is very complicated but meaningful.Table 1
Summarized literature review.
Papers Integrated systems Production-inventory model Sampling interval Solution approach Orthogonal experiment EPQ HPP Fixed Variable Monte Carlo Meta-heuristic method [10] MS & PIS  √ [28] MS & PIS √ [17] MS & PIS [11] MS & PIS [24] MS & PIS [29] MS & QS    √  √ [27] MS & QS   √   √ √ [30] MS & QS [31] MS & QS [32] MS & QS & PIS √   √  √ [21] MS & QS & PIS √ [33] MS & QS & PIS √  √ [19] MS & QS & PIS √ [34] MS & QS & PIS √ [35] MS & QS & PIS √ [36] MS & QS & PIS  √ [37] MS & QS & PIS √ [38] MS & QS & PIS  √   √ [20] MS & QS & PIS √ This paper MS & QS & PIS  √  √ √ √ √The remainder of this paper is organized as follows. In the next section, we specify targeted contributions of this paper. Section3 describes the notations and the problem. In Section 4, a model with 16 scenarios is developed. In Section 5, the developed model is solved by the Monte Carlo method and the genetic algorithm. In Section 6, an illustrative numerical example is presented. In addition, comparative and sensitivity analyses are also conducted in this section. Finally, Section 7 concludes the paper.


## 2. Targeted Contributions

By integrating all the three systems, this paper aims at minimizing the expected cost per unit time by determining 5 decision variables concerning the three systems. Comparing with the existing literature, the main contributions of this paper are the following:( 1 ) The process inspection is often integrated with the economic production quantity model as the production-inventory policy in the previous integration models. However, the hedging point policy is also a very important production-inventory policy and it often interacts with the process inspection in practice. To solve the manufacturing problem with the process inspection and the hedging point policy, this paper builds an optimization model with 16 scenarios.( 2 ) The Monte Carlo method is directly used to calculate the objective function in the process of a metaheuristic algorithm. Although some researchers have managed to use the Monte Carlo method to calculate objective functions in the process of genetic algorithms [25, 26], the process of genetic algorithms has to be modified because of the low speed of the Monte Carlo simulation. However, this paper proposes an approach speeding up the Monte Carlo simulation dramatically. This makes the direct application of the Monte Carlo method in the process of metaheuristic algorithms practical. As an experiment, this paper uses the Monte Carlo method to calculate the objective function in the process of the genetic algorithm, which is one of the most important metaheuristic algorithms.( 3 ) To assess the accuracy of the Monte Carlo method, this paper also calculates the expected cost per unit time by the integral method. Because 16 scenarios may occur, this calculation is very difficult. However, this paper solves it successfully by calculating conditional expectations by conditioning. This calculation method can be used in other manufacturing problems involving expectations of complex random variables.( 4 ) In the majority of previous papers, the sensitivity analysis was conducted by single factor experiments. In this paper, an orthogonal experiment is used to conduct the sensitivity analysis, which can make the sensitivity analysis more comprehensive with fewer trials.


## 3. Problem Statement

### 3.1. Notation

The notations used in this paper are defined as follows:Decision Variables H: Hedging level (safety stock) h 1: Duration of the first sampling interval L: Coefficient of the control limit k: Maximum number of samples in a production cycle n: Sample sizeModel Parameters τ: Random variable denoting the duration of the in-control state in a production cycle f τ ( t ): Probability density function of τ λ τ ( t ): Hazard rate function of τ X: Random variable denoting the time to failure of the machine λ X ( t ): Hazard rate function of X t p: Random variable denoting the preventive maintenance duration f p ( t ): Probability density function of tp t r: Random variable denoting the predictive maintenance duration f r ( t ): Probability density function of tr y: Number of samples needed to detect the state shift after it occurs in a production cycle, considering the type II error K: Setup cost d: Demand rate U m a x: Maximum production rate θ 1: Proportion of imperfect products when the production process is in the in-control state θ 2: Proportion of imperfect products when the production process is in the out-of-control state w: Warranty period λ 1 ( t ): Hazard rate function of perfect products λ 2 ( t ): Hazard rate function of imperfect products C w: Cost per product at each minimum repair for the warranty j: Index of sampling intervals μ 0: Mean value of the quality characteristic in the in-control state σ: Standard deviation of the quality characteristic δ: Magnitude of the shift in the mean of quality characteristic when the assignable cause occurs C F: Fixed cost per sample C V: Variable cost per sample C c: Cost of one corrective maintenance action C p: Cost of performing preventive maintenance per unit time C r: Cost of performing predictive maintenance per unit time C y: Cost of inspecting an alert signal h: Holding cost of a product unit during a unit of time C s: Shortage cost of a product unit h j: Length of the jth sampling interval w j: Time at the end of the jth sampling interval α: Probability of type I error β: Probability of type II error t b: Duration of the build-up of the buffer stock t d: Duration of the decrease of the buffer stock N i: Number of imperfect products produced in a production cycle N p: Number of perfect products produced in a production cycle C: Corrective maintenance cost in a production cycle P: Preventive maintenance cost in a production cycle R: Predictive maintenance cost in a production cycle B: Sampling cost in a production cycle A: Cost of inspecting alert signals in a production cycle I: Inventory cost in a production cycle F: Cost for the free minimal repair warranty incurred in a production cycle S: Shortage cost in a production cycle K s: Setup cost in a production cycle C a: Total cost of a production cycle T c: Length of a production cycleOther notations will be introduced where they are needed.


### 3.2. Problem Description

In the following, the policies for the production-inventory system, the maintenance system, and the quality system are explained.


#### 3.2.1. Production-Inventory System

The so-called hedging point policy is used to instantly control the production rate of the machineu(.) as follows:(1)ut=Umaxifxt<Handαt=1difxt=Handαt=10otherwisewhere x(t) is the inventory level at time t and α(t) denotes the state of the machine at time t. α(t)=1 if the machine is available for production. α(t)=0 if the machine is under maintenance.In addition, if the inventory level is greater than zero when the maintenance is complete, the machine will not be started until the inventory level is reduced to zero. If the inventory is not available to adjust demand during the interruption due to maintenance, customer demand is lost at a relevant cost.


#### 3.2.2. Maintenance System

This paper investigates an imperfect manufacturing process, which includes two states: the in-control state and the out-of-control state. The manufacturing process starts operating in the in-control state and one assignable cause may occur and lead to shifting process to the out-of-control state. The quality characteristic of products follows the normal distribution with mean valueμ0 and standard deviation σ in the in-control state. When the process shifts from the in-control state to the out-of-control state, the mean value changes from μ0 to μ1=μ0±δσ, where δ is the magnitude of quality shift, and the standard deviation is assumed to remain the same.A Shewhart X-bar control chart is used to monitor the quality characteristic with an alert signal to inform operators when the process shifts to the out-of-control state. When the process mean falls inside the control limits [LCL,UCL], the process is considered to be in the in-control state. Once it falls outside the control limits, an alert signal is generated since the assignable cause probably has affected the process. Here, LCL=μ0-L(σ/n) and UCL=μ0+L(σ/n), where L is the coefficient of the control limit and n is the sample size.The alert signal of the control chart triggers an error-free inspection to ascertain its authenticity. If the assignable cause indeed happens, the predictive maintenance (PRDM) is immediately carried out to restore the process to the in-control state [27]. If the alert signal is false (type I error), the process will continue operating without any maintenance action. Once the equipment fails, the corrective maintenance (CM) has to be performed at once to restore the equipment and the process. The CM is minimal repair and the repair time is negligible. If no PRDM is conducted after k samplings, there may be two possible conditions: (i) the process is in the in-control state and (ii) the process is in the out-of-control state but is undetected (type II error). Therefore, it is necessary to stop the production process to implement the preventive maintenance (PM) at the end of the (k+1)th sampling interval (without sampling). If the process is in the out-of-control state, the assignable cause will be detected during the PM. Thus the PM has to be replaced by the PRDM. The PM and the PRDM are considered to restore the equipment to the good-as-new condition.The lengths of sampling intervals are determined such that the integrated hazard over them is constant. That is,(2)∫wjwj+1λτtdt=∫0h1λτtdtwhere wj=∑i=1jhi,j=1,2,…,w0=0. τ is assumed to follow the Weibull distribution. Thus, its hazard rate function can be assumed to be λτ(t)=λνtν-1,λ>1,ν≥1. Therefore, the wj can be obtained:(3)wj=j1/v×h1.


#### 3.2.3. Quality System

When the system is in the in-control state, it produces some imperfect products. In the out-of-control state, the system produces more imperfect products than in the in-control state. Both the perfect and imperfect products will be sold with warranty. The free minimal repair warranty (FRW) policy is adopted. That is, all products sold can obtain free minimal repair when failing and the repair time is negligible. In addition, the manufacturer will pay for the repair.


## 3.1. Notation

The notations used in this paper are defined as follows:Decision Variables H: Hedging level (safety stock) h 1: Duration of the first sampling interval L: Coefficient of the control limit k: Maximum number of samples in a production cycle n: Sample sizeModel Parameters τ: Random variable denoting the duration of the in-control state in a production cycle f τ ( t ): Probability density function of τ λ τ ( t ): Hazard rate function of τ X: Random variable denoting the time to failure of the machine λ X ( t ): Hazard rate function of X t p: Random variable denoting the preventive maintenance duration f p ( t ): Probability density function of tp t r: Random variable denoting the predictive maintenance duration f r ( t ): Probability density function of tr y: Number of samples needed to detect the state shift after it occurs in a production cycle, considering the type II error K: Setup cost d: Demand rate U m a x: Maximum production rate θ 1: Proportion of imperfect products when the production process is in the in-control state θ 2: Proportion of imperfect products when the production process is in the out-of-control state w: Warranty period λ 1 ( t ): Hazard rate function of perfect products λ 2 ( t ): Hazard rate function of imperfect products C w: Cost per product at each minimum repair for the warranty j: Index of sampling intervals μ 0: Mean value of the quality characteristic in the in-control state σ: Standard deviation of the quality characteristic δ: Magnitude of the shift in the mean of quality characteristic when the assignable cause occurs C F: Fixed cost per sample C V: Variable cost per sample C c: Cost of one corrective maintenance action C p: Cost of performing preventive maintenance per unit time C r: Cost of performing predictive maintenance per unit time C y: Cost of inspecting an alert signal h: Holding cost of a product unit during a unit of time C s: Shortage cost of a product unit h j: Length of the jth sampling interval w j: Time at the end of the jth sampling interval α: Probability of type I error β: Probability of type II error t b: Duration of the build-up of the buffer stock t d: Duration of the decrease of the buffer stock N i: Number of imperfect products produced in a production cycle N p: Number of perfect products produced in a production cycle C: Corrective maintenance cost in a production cycle P: Preventive maintenance cost in a production cycle R: Predictive maintenance cost in a production cycle B: Sampling cost in a production cycle A: Cost of inspecting alert signals in a production cycle I: Inventory cost in a production cycle F: Cost for the free minimal repair warranty incurred in a production cycle S: Shortage cost in a production cycle K s: Setup cost in a production cycle C a: Total cost of a production cycle T c: Length of a production cycleOther notations will be introduced where they are needed.


## 3.2. Problem Description

In the following, the policies for the production-inventory system, the maintenance system, and the quality system are explained.


### 3.2.1. Production-Inventory System

The so-called hedging point policy is used to instantly control the production rate of the machineu(.) as follows:(1)ut=Umaxifxt<Handαt=1difxt=Handαt=10otherwisewhere x(t) is the inventory level at time t and α(t) denotes the state of the machine at time t. α(t)=1 if the machine is available for production. α(t)=0 if the machine is under maintenance.In addition, if the inventory level is greater than zero when the maintenance is complete, the machine will not be started until the inventory level is reduced to zero. If the inventory is not available to adjust demand during the interruption due to maintenance, customer demand is lost at a relevant cost.


### 3.2.2. Maintenance System

This paper investigates an imperfect manufacturing process, which includes two states: the in-control state and the out-of-control state. The manufacturing process starts operating in the in-control state and one assignable cause may occur and lead to shifting process to the out-of-control state. The quality characteristic of products follows the normal distribution with mean valueμ0 and standard deviation σ in the in-control state. When the process shifts from the in-control state to the out-of-control state, the mean value changes from μ0 to μ1=μ0±δσ, where δ is the magnitude of quality shift, and the standard deviation is assumed to remain the same.A Shewhart X-bar control chart is used to monitor the quality characteristic with an alert signal to inform operators when the process shifts to the out-of-control state. When the process mean falls inside the control limits [LCL,UCL], the process is considered to be in the in-control state. Once it falls outside the control limits, an alert signal is generated since the assignable cause probably has affected the process. Here, LCL=μ0-L(σ/n) and UCL=μ0+L(σ/n), where L is the coefficient of the control limit and n is the sample size.The alert signal of the control chart triggers an error-free inspection to ascertain its authenticity. If the assignable cause indeed happens, the predictive maintenance (PRDM) is immediately carried out to restore the process to the in-control state [27]. If the alert signal is false (type I error), the process will continue operating without any maintenance action. Once the equipment fails, the corrective maintenance (CM) has to be performed at once to restore the equipment and the process. The CM is minimal repair and the repair time is negligible. If no PRDM is conducted after k samplings, there may be two possible conditions: (i) the process is in the in-control state and (ii) the process is in the out-of-control state but is undetected (type II error). Therefore, it is necessary to stop the production process to implement the preventive maintenance (PM) at the end of the (k+1)th sampling interval (without sampling). If the process is in the out-of-control state, the assignable cause will be detected during the PM. Thus the PM has to be replaced by the PRDM. The PM and the PRDM are considered to restore the equipment to the good-as-new condition.The lengths of sampling intervals are determined such that the integrated hazard over them is constant. That is,(2)∫wjwj+1λτtdt=∫0h1λτtdtwhere wj=∑i=1jhi,j=1,2,…,w0=0. τ is assumed to follow the Weibull distribution. Thus, its hazard rate function can be assumed to be λτ(t)=λνtν-1,λ>1,ν≥1. Therefore, the wj can be obtained:(3)wj=j1/v×h1.


### 3.2.3. Quality System

When the system is in the in-control state, it produces some imperfect products. In the out-of-control state, the system produces more imperfect products than in the in-control state. Both the perfect and imperfect products will be sold with warranty. The free minimal repair warranty (FRW) policy is adopted. That is, all products sold can obtain free minimal repair when failing and the repair time is negligible. In addition, the manufacturer will pay for the repair.


## 3.2.1. Production-Inventory System

The so-called hedging point policy is used to instantly control the production rate of the machineu(.) as follows:(1)ut=Umaxifxt<Handαt=1difxt=Handαt=10otherwisewhere x(t) is the inventory level at time t and α(t) denotes the state of the machine at time t. α(t)=1 if the machine is available for production. α(t)=0 if the machine is under maintenance.In addition, if the inventory level is greater than zero when the maintenance is complete, the machine will not be started until the inventory level is reduced to zero. If the inventory is not available to adjust demand during the interruption due to maintenance, customer demand is lost at a relevant cost.


## 3.2.2. Maintenance System

This paper investigates an imperfect manufacturing process, which includes two states: the in-control state and the out-of-control state. The manufacturing process starts operating in the in-control state and one assignable cause may occur and lead to shifting process to the out-of-control state. The quality characteristic of products follows the normal distribution with mean valueμ0 and standard deviation σ in the in-control state. When the process shifts from the in-control state to the out-of-control state, the mean value changes from μ0 to μ1=μ0±δσ, where δ is the magnitude of quality shift, and the standard deviation is assumed to remain the same.A Shewhart X-bar control chart is used to monitor the quality characteristic with an alert signal to inform operators when the process shifts to the out-of-control state. When the process mean falls inside the control limits [LCL,UCL], the process is considered to be in the in-control state. Once it falls outside the control limits, an alert signal is generated since the assignable cause probably has affected the process. Here, LCL=μ0-L(σ/n) and UCL=μ0+L(σ/n), where L is the coefficient of the control limit and n is the sample size.The alert signal of the control chart triggers an error-free inspection to ascertain its authenticity. If the assignable cause indeed happens, the predictive maintenance (PRDM) is immediately carried out to restore the process to the in-control state [27]. If the alert signal is false (type I error), the process will continue operating without any maintenance action. Once the equipment fails, the corrective maintenance (CM) has to be performed at once to restore the equipment and the process. The CM is minimal repair and the repair time is negligible. If no PRDM is conducted after k samplings, there may be two possible conditions: (i) the process is in the in-control state and (ii) the process is in the out-of-control state but is undetected (type II error). Therefore, it is necessary to stop the production process to implement the preventive maintenance (PM) at the end of the (k+1)th sampling interval (without sampling). If the process is in the out-of-control state, the assignable cause will be detected during the PM. Thus the PM has to be replaced by the PRDM. The PM and the PRDM are considered to restore the equipment to the good-as-new condition.The lengths of sampling intervals are determined such that the integrated hazard over them is constant. That is,(2)∫wjwj+1λτtdt=∫0h1λτtdtwhere wj=∑i=1jhi,j=1,2,…,w0=0. τ is assumed to follow the Weibull distribution. Thus, its hazard rate function can be assumed to be λτ(t)=λνtν-1,λ>1,ν≥1. Therefore, the wj can be obtained:(3)wj=j1/v×h1.


## 3.2.3. Quality System

When the system is in the in-control state, it produces some imperfect products. In the out-of-control state, the system produces more imperfect products than in the in-control state. Both the perfect and imperfect products will be sold with warranty. The free minimal repair warranty (FRW) policy is adopted. That is, all products sold can obtain free minimal repair when failing and the repair time is negligible. In addition, the manufacturer will pay for the repair.


## 4. Model Development

First of all some key variables need to be calculated. The duration of the build-up of the buffer stock is(4)tb=HUmax-d.Under the FRW, a perfect product fails according to a nonhomogeneous Poisson process (NHPP) with rate λ1(t), while an imperfect product fails according to a NHPP with rate λ2(t). Therefore, given a perfect product, the mean number of failures within the warranty period w is ∫0wλ1(t)dt. On the other hand, given an imperfect product, the mean number of failures within the warranty period w is ∫0wλ2(t)dt. Consequently, the mean warranty costs of a perfect product and an imperfect product are, respectively,(5)Cpw=Cw∫0wλ1tdt,(6)Ciw=Cw∫0wλ2tdt.The whole cost of each sampling is (7)s=CF+nCV.The duration of the decrease of the buffer stock is (8)td=Hd.The probability of type I error is (9)α=2Φ-Lwhere Φ(.) represents the cdf of standard normal distribution. The probability of type II error is (10)β=ΦL-δn-Φ-L-δn.y is a geometric random variable with parameter 1-β. Thus,(11)Py=y0=βy0-11-β.


### 4.1. The 16 Scenarios

The period between the two successive starts of the good-as-new equipment is described as a production cycle. Thus, according to the policies proposed in Section3.2, 16 scenarios may occur in a production cycle. This section contributes to describing these scenarios in detail and formulating the conditional expectations of Ca and Tc for different scenarios.


#### 4.1.1. Scenario 1:0<τ<tb, 1≤y≤k, 0<tr<td

This scenario can be rewritten asτ=τ0,y=y0,tr=tr0 and 0<τ0<tb, 1≤y0≤k, 0<tr0<td. Then the behavior of inventory in this scenario can be described as Figure 1. In this scenario, the process shifts to out-of-control state before the hedging safety stock H is built. Due to the type II error, the alert signal is not sent until the y0th sampling. Then the PRDM is conducted for tr0. Because 0<tr0<td, shortage does not occur. The conditional expectations of Ca and Tc given that τ=τ0,y=y0,tr=tr0, that is, E(Ca∣τ=τ0,y=y0,tr=tr0) and E(Tc∣τ=τ0,y=y0,tr=tr0), when 0<τ0<tb, 1≤y0≤k, 0<tr0<td will be formulated.Figure 1
Behavior of inventory (Scenario 1).BecauseCa=C+Ks+P+R+B+A+I+F+S, E(Ca∣τ=τ0,y=y0,tr=tr0) can be obtained through calculating the conditional expectations of C,Ks,P,R,B,A,I,F, and S. The CM is minimal repair and the repair time is negligible; hence the machine fails according to an NHPP with rate λX(t). Therefore, the mean number of failures within the production period wy0 is ∫0wy0λX(t)dt. Consequently, the conditional expectations of C are(12)EC∣τ=τ0,y=y0,tr=tr0=Cc∫0wy0λXtdt.The conditional expectations of Ks are(13)EKs∣τ=τ0,y=y0,tr=tr0=K∫0wy0λXtdt+1.The PM is not conducted; hence(14)EP∣τ=τ0,y=y0,tr=tr0=0.The conditional expectations of R are(15)ER∣τ=τ0,y=y0,tr=tr0=Crtr0.The conditional expectations of B are(16)EB∣τ=τ0,y=y0,tr=tr0=y0s.The conditional expectations of A are(17)EA∣τ=τ0,y=y0,tr=tr0=Cy.The conditional expectations of I are(18)EI∣τ=τ0,y=y0,tr=tr0=h12tbH+wy0-tbH+12Htd.The conditional expectations of Ni are(19)ENi∣τ=τ0,y=y0,tr=tr0=tb-τ0Umax+wy0-tbdθ2+Umaxτ0θ1.The conditional expectations of Np are(20)ENp∣τ=τ0,y=y0,tr=tr0=tb-τ0Umax+wy0-tbd1-θ2+Umaxτ01-θ1.Hence,(21)EF∣τ=τ0,y=y0,tr=tr0=tb-τ0Umax+wy0-tbdθ2+Umaxτ0θ1Ciw+tb-τ0Umax+wy0-tbd1-θ2+Umaxτ01-θ1Cpw.Shortage does not occur; hence(22)ES∣τ=τ0,y=y0,tr=tr0=0.Above all, (23)ECa∣τ=τ0,y=y0,tr=tr0=Cc∫0wy0λXtdt+K∫0wy0λXtdt+1+Crtr0+y0s+Cy+h12tbH+wy0-tbH+12Htd+tb-τ0Umax+wy0-tbdθ2+Umaxτ0θ1Ciw+tb-τ0Umax+wy0-tbd1-θ2+Umaxτ01-θ1Cpw=Ea1τ0,y0,tr0where Ea1(τ0,y0,tr0) denotes E(Ca∣τ=τ0,y=y0,tr=tr0) when 0<τ0<tb, 1≤y0≤k, 0<tr0<td.From Figure1, it can be obtained that(24)ETc∣τ=τ0,y=y0,tr=tr0=wy0+td=Ec1y0where Ec1(y0) denotes E(Tc∣τ=τ0,y=y0,tr=tr0) when 0<τ0<tb, 1≤y0≤k, 0<tr0<td.


#### 4.1.2. Scenario 2:0<τ<tb, 1≤y≤k, tr>td

This scenario can be rewritten asτ=τ0,y=y0,tr=tr0 and 0<τ0<tb, 1≤y0≤k, tr0>td. Thus the behavior of inventory in this scenario can be described as Figure 2. In contrast to Scenario 1, shortage occurs since tr0>td. Hence(25)ES∣τ=τ0,y=y0,tr=tr0=Cstr0-tdd.Other costs are the same as Scenario 1. Therefore(26)ECa∣τ=τ0,y=y0,tr=tr0=Ea1τ0,y0,tr0+Cstr0-tdd=Ea2τ0,y0,tr0where Ea2(τ0,y0,tr0) denotes E(Ca∣τ=τ0,y=y0,tr=tr0) when 0<τ0<tb, 1≤y0≤k, tr0>td. From Figure 2, it can be obtained that(27)ETc∣τ=τ0,y=y0,tr=tr0=wy0+tr0=Ec2y0,tr0where Ec2(y0,tr0) denotes E(Tc∣τ=τ0,y=y0,tr=tr0) when 0<τ0<tb, 1≤y0≤k, tr0>td.Figure 2
Behavior of inventory (Scenario 2).


#### 4.1.3. Scenario 3:0<τ<tb, y≥k+1, 0<tr<td

This scenario can be rewritten asτ=τ0,y=y0,tr=tr0 and 0<τ0<tb, y0≥k+1, 0<tr0<td. Thus the behavior of inventory in this scenario can be described as Figure 3. The alert signal is not sent in the first k samples. Hence the preventive maintenance (PM) is implemented at the end of the (k+1)th sampling interval (without sampling). Then the assignable cause is detected during the PM. Thus, the PM is replaced by the PRDM. Because 0<tr0<td, shortage does not occur.Figure 3
Behavior of inventory (Scenario 3).The conditional expectations of all kinds of costs given thatτ=τ0,y=y0,tr=tr0 can be calculated using the same analogy as Scenario 1. And they are displayed in Tables 1 and 2 of the Annex (Supplementary 1). Hence(28)ECa∣τ=τ0,y=y0,tr=tr0=Cc∫0wk+1λXtdt+K∫0wk+1λXtdt+1+Crtr0+ks+h12tbH+wk+1-tbH+12Htd+tb-τ0Umax+wk+1-tbdθ2+Umaxτ0θ1Ciw+tb-τ0Umax+wk+1-tbd1-θ2+Umaxτ01-θ1Cpw=Ea3τ0,tr0where Ea3(τ0,tr0) denotes E(Ca∣τ=τ0,y=y0,tr=tr0) when 0<τ0<tb, y0≥k+1, 0<tr0<td.From Figure3, it can be obtained that(29)ETc∣τ=τ0,y=y0,tr=tr0=wk+1+td=Ec3where Ec3 denotes E(Tc∣τ=τ0,y=y0,tr=tr0) when 0<τ0<tb, y0≥k+1, 0<tr0<td.


#### 4.1.4. Scenario 4:0<τ<tb, y≥k+1, tr>td

This scenario can be rewritten asτ=τ0,y=y0,tr=tr0 and 0<τ0<tb, y0≥k+1, tr0>td. Thus the behavior of inventory in this scenario can be described as Figure 4. In contrast to Scenario 3, shortage occurs since tr0>td. Hence(30)ES∣τ=τ0,y=y0,tr=tr0=Cstr0-tdd.Other costs are the same as Scenario 3. Therefore(31)Ea4τ0,tr0=Ea3τ0,tr0+ES∣τ=τ0,y=y0,tr=tr0=Ea3τ0,tr0+Cstr0-tddwhere Ea4(τ0,tr0) denotes E(Ca∣τ=τ0,y=y0,tr=tr0) when 0<τ0<tb, y0≥k+1, tr0>td. From Figure 4, it can be obtained that(32)ETc∣τ=τ0,y=y0,tr=tr0=wk+1+tr0=Ec4tr0where Ec4(tr0) denotes E(Tc∣τ=τ0,y=y0,tr=tr0) when 0<τ0<tb, y0≥k+1, tr0>td.Figure 4
Behavior of inventory (Scenario 4).


#### 4.1.5. Scenario 5:tb<τ<w1, 1≤y≤k, 0<tr<td

This scenario can be rewritten asτ=τ0,y=y0,tr=tr0 and tb<τ0<w1, 1≤y0≤k, 0<tr0<td. Thus the behavior of inventory can be described as Figure 5. In this scenario, the process shifts to out-of-control state after the hedging safety stock H is built and before the first sampling. Due to the type II error, the alert signal is not sent until the y0th sampling. Then the PRDM is conducted for tr0. Because 0<tr0<td, shortage does not occur.Figure 5
Behavior of inventory (Scenario 5).The conditional expectations of all kinds of costs given thatτ=τ0,y=y0,tr=tr0 can be calculated using the same analogy as Scenario 1. And they are displayed in Tables 1 and 2 of the Annex. Therefore (33)ECa∣τ=τ0,y=y0,tr=tr0=Cc∫0wy0λXtdt+K∫0wy0λXtdt+1+Crtr0+y0s+Cy+h12tbH+wy0-tbH+12Htd+wy0-τ0dθ2+Umaxtb+τ0-tbdθ1Ciw+wy0-τ0d1-θ2+Umaxtb+τ0-tbd1-θ1Cpw=Ea5τ0,y0,tr0where Ea5(τ0,y0,tr0) denotes E(Ca∣τ=τ0,y=y0,tr=tr0) when tb<τ0<w1, 1≤y0≤k, 0<tr0<td.From Figure5, it can be obtained that(34)ETc∣τ=τ0,y=y0,tr=tr0=wy0+td=Ec5y0where Ec5(y0) denotes E(Tc∣τ=τ0,y=y0,tr=tr0) when tb<τ0<w1, 1≤y0≤k, 0<tr0<td.


#### 4.1.6. Scenario 6:tb<τ<w1, 1≤y≤k, tr>td

This scenario can be rewritten asτ=τ0,y=y0,tr=tr0 and tb<τ0<w1, 1≤y0≤k, tr0>td. Thus the behavior of inventory can be described as Figure 6. In contrast to Scenario 5, shortage occurs since tr0>td. Hence(35)ES∣τ=τ0,y=y0,tr=tr0=Cstr0-tdd.Other costs are the same as Scenario 5. Therefore(36)Ea6τ0,y0,tr0=Ea5τ0,y0,tr0+ES∣τ=τ0,y=y0,tr=tr0=Ea5τ0,y0,tr0+Cstr0-tddwhere Ea6(τ0,y0,tr0) denotes E(Ca∣τ=τ0,y=y0,tr=tr0) when tb<τ0<w1, 1≤y0≤k, tr0>td. From Figure 6, it can be obtained that(37)ETc∣τ=τ0,y=y0,tr=tr0=wy0+tr0=Ec6y0,tr0where Ec6(y0,tr0) denotes ETc∣τ=τ0,y=y0,tr=tr0 when tb<τ0<w1, 1≤y0≤k, tr0>td.Figure 6
Behavior of inventory (Scenario 6).


#### 4.1.7. Scenario 7:tb<τ<w1, y≥k+1, 0<tr<td

This scenario can be rewritten asτ=τ0,y=y0,tr=tr0 and tb<τ0<w1, y0≥k+1, 0<tr0<td. Thus the behavior of inventory can be described as Figure 7. The alert signal is not sent in the first k samples. Hence the preventive maintenance (PM) is implemented at the end of the (k+1)th sampling interval (without sampling). Then the assignable cause is detected during the PM. Thus the PM is replaced by PRDM. Because 0<tr0<td, shortage does not occur.Figure 7
Behavior of inventory (Scenario 7).The conditional expectations of all kinds of costs given thatτ=τ0,y=y0,tr=tr0 can be calculated using the same analogy as Scenario 1. And they are displayed in Tables 1 and 2 of the Annex. Therefore(38)ECa∣τ=τ0,y=y0,tr=tr0=Cc∫0wk+1λXtdt+K∫0wk+1λXtdt+1+Crtr0+ks+h12tbH+wk+1-tbH+12Htd+wk+1-τ0dθ2+Umaxtb+τ0-tbdθ1Ciw+wk+1-τ0d1-θ2+Umaxtb+τ0-tbd1-θ1Cpw=Ea7τ0,tr0where Ea7(τ0,tr0) denotes E(Ca∣τ=τ0,y=y0,tr=tr0) when tb<τ0<w1, y0≥k+1, 0<tr0<td.From Figure7, it can be obtained that(39)ETc∣τ=τ0,y=y0,tr=tr0=wk+1+td=Ec7where Ec7 denotes E(Tc∣τ=τ0,y=y0,tr=tr0) when tb<τ0<w1, y0≥k+1, 0<tr0<td.


#### 4.1.8. Scenario 8:tb<τ<w1, y≥k+1, tr>td

This scenario can be rewritten asτ=τ0,y=y0,tr=tr0 and tb<τ0<w1, y0≥k+1, tr0>td. Thus the behavior of inventory can be described as Figure 8. In contrast to Scenario 7, shortage occurs since tr0>td. Hence(40)ES∣τ=τ0,y=y0,tr=tr0=Cstr0-tdd.Other costs are the same as Scenario 7. Therefore(41)Ea8τ0,tr0=Ea7τ0,tr0+ES∣τ=τ0,y=y0,tr=tr0=Ea7τ0,tr0+Cstr0-tddwhere Ea8(τ0,tr0) denotes E(Ca∣τ=τ0,y=y0,tr=tr0) when tb<τ0<w1, y0≥k+1, tr0>td. From Figure 8, it can be obtained that(42)ETc∣τ=τ0,y=y0,tr=tr0=wk+1+tr0=Ec8tr0where Ec8 denotes E(Tc∣τ=τ0,y=y0,tr=tr0) when tb<τ0<w1, y0≥k+1, tr0>td.Figure 8
Behavior of inventory (Scenario 8).


#### 4.1.9. Scenario 9:wj<τ<wj+1(1≤j≤k-1), 1≤y≤k-j, tr<td

This scenario can be rewritten asτ=τ0,y=y0,tr=tr0 and wj<τ0<wj+1(1≤j≤k-1),1≤y0≤k-j,tr0<td. Thus the behavior of inventory can be described as Figure 9. In this scenario, the process shifts to out-of-control state in the (j+1)th sampling interval. After the shift, due to the type II error, the alert signal is not sent until the (j+y0)th sampling, while the alert signal may be sent in the first j+1 samplings due to the type I error. After the (j+y0)th sampling, the PRDM is conducted for tr0. Because 0<tr0<td, shortage does not occur.Figure 9
Behavior of inventory (Scenario 9).The conditional expectations of all kinds of costs given thatτ=τ0,y=y0,tr=tr0 can be calculated using the same analogy as Scenario 1. And they are displayed in Tables 1 and 2 of the Annex. Therefore(43)ECa∣τ=τ0,y=y0,tr=tr0=Cc∫0wj+y0λXtdt+K∫0wj+y0λXtdt+1+Crtr0+j+y0s+Cyjα+1+h12tbH+wj+y0-tbH+12Htd+wj+y0-τ0dθ2+Umaxtb+τ0-tbdθ1Ciw+wj+y0-τ0d1-θ2+Umaxtb+τ0-tbd1-θ1Cpw=Ea9τ0,y0,tr0where Ea9(τ0,y0,tr0) denotes E(Ca∣τ=τ0,y=y0,tr=tr0) when wj<τ0<wj+1(1≤j≤k-1),1≤y0≤k-j,tr0<td.From Figure9, it can be obtained that(44)ETc∣τ=τ0,y=y0,tr=tr0=wj+y0+td=Ec9y0where Ec9(y0) denotes E(Tc∣τ=τ0,y=y0,tr=tr0) when wj<τ0<wj+1(1≤j≤k-1),1≤y0≤k-j,tr0<td.


#### 4.1.10. Scenario 10:wj<τ<wj+1(1≤j≤k-1), 1≤y≤k-j, tr>td

This scenario can be rewritten asτ=τ0,y=y0,tr=tr0 and wj<τ0<wj+1(1≤j≤k-1),1≤y0≤k-j,tr0>td. Thus the behavior of inventory can be described as Figure 10. In contrast to Scenario 9, shortage occurs since tr0>td. Hence(45)ES∣τ=τ0,y=y0,tr=tr0=Cstr0-tdd.Other costs are the same as Scenario 9. Therefore(46)Ea10τ0,y0,tr0=Ea9τ0,y0,tr0+ES∣τ=τ0,y=y0,tr=tr0=Ea9τ0,y0,tr0+Cstr0-tddwhere Ea10(τ0,y0,tr0) denotes E(Ca∣τ=τ0,y=y0,tr=tr0) when wj<τ0<wj+1(1≤j≤k-1),1≤y0≤k-j,tr0>td. From Figure 10, it can be obtained that(47)ETc∣τ=τ0,y=y0,tr=tr0=wj+y0+tr0=Ec10tr0where Ec10(tr0) denotes E(Tc∣τ=τ0,y=y0,tr=tr0) when wj<τ0<wj+1(1≤j≤k-1),1≤y0≤k-j,tr0>td.Figure 10
Behavior of inventory (Scenario 10).


#### 4.1.11. Scenario 11:wj<τ<wj+1(1≤j≤k-1), y>k-j, tr<td

This scenario can be rewritten asτ=τ0,y=y0,tr=tr0 and wj<τ0<wj+1(1≤j≤k-1),y0>k-j,tr0<td. Thus the behavior of inventory can be described as Figure 11. No PRDM is conducted after k samplings. Hence the preventive maintenance (PM) is implemented at the end of the (k+1)th sampling interval (without sampling). Then the assignable cause is detected during the PM. Thus the PM is replaced by the PRDM. Because 0<tr0<td, shortage does not occur.Figure 11
Behavior of inventory (Scenario 11).The conditional expectations of all kinds of costs given thatτ=τ0,y=y0,tr=tr0 can be calculated using the same analogy as Scenario 1. And they are displayed in Tables 1 and 2 of the Annex. Therefore (48)ECa∣τ=τ0,y=y0,tr=tr0=Cc∫0wk+1λXtdt+K∫0wk+1λXtdt+1+Crtr0+ks+Cyjα+h12tbH+wk+1-tbH+12Htd+wk+1-τ0dθ2+Umaxtb+τ0-tbdθ1Ciw+wk+1-τ0d1-θ2+Umaxtb+τ0-tbd1-θ1Cpw=Ea11τ0,tr0where Ea11(τ0,tr0) denotes E(Ca∣τ=τ0,y=y0,tr=tr0) when wj<τ0<wj+1(1≤j≤k-1),y0>k-j,tr0<td.From Figure11, it can be obtained that(49)ETc∣τ=τ0,y=y0,tr=tr0=wk+1+td=Ec11where Ec11 denotes E(Tc∣τ=τ0,y=y0,tr=tr0) when wj<τ0<wj+1(1≤j≤k-1),y0>k-j,tr0<td.


#### 4.1.12. Scenario 12:wj<τ<wj+1(1≤j≤k-1), y>k-j, tr>td

This scenario can be rewritten asτ=τ0,y=y0,tr=tr0 and wj<τ0<wj+1(1≤j≤k-1),y0>k-j,tr0>td. Thus the behavior of inventory can be described as Figure 12. In contrast to Scenario 11, shortage occurs since tr0>td. Hence(50)ES∣τ=τ0,y=y0,tr=tr0=Cstr0-tdd.Other costs are the same as Scenario 11. Therefore(51)Ea12τ0,tr0=Ea11τ0,tr0+ES∣τ=τ0,y=y0,tr=tr0=Ea11τ0,tr0+Cstr0-tddwhere Ea12(τ0,tr0) denotes E(Ca∣τ=τ0,y=y0,tr=tr0) when wj<τ0<wj+1(1≤j≤k-1),y0>k-j,tr0>td. From Figure 12, it can be obtained that(52)ETc∣τ=τ0,y=y0,tr=tr0=wk+1+tr0=Ec12tr0where Ec12(tr0) denotes E(Tc∣τ=τ0,y=y0,tr=tr0) when wj<τ0<wj+1(1≤j≤k-1),y0>k-j,tr0>td.Figure 12
Behavior of inventory (Scenario 12).


#### 4.1.13. Scenario 13:wk<τ<wk+1,tr<td

This scenario can be rewritten asτ=τ0,tr=tr0 and wk<τ0<wk+1,tr0<td. Thus the behavior of inventory can be described as Figure 13. In this scenario, the process shifts to out-of-control state in the (k+1)th sampling interval. No PRDM is conducted after k samplings. Hence the preventive maintenance (PM) is implemented at the end of the (k+1)th sampling interval (without sampling). Then the assignable cause is detected during the PM. Thus the PM is replaced by PRDM. Because 0<tr0<td, shortage does not occur.Figure 13
Behavior of inventory (Scenario 13).The conditional expectations of all kinds of costs given thatτ=τ0,tr=tr0 can be calculated using the same analogy as Scenario 1. And they are displayed in Tables 1 and 2 of the Annex. Therefore (53)ECa∣τ=τ0,tr=tr0=Cc∫0wk+1λXtdt+K∫0wk+1λXtdt+1+Crtr0+ks+Cykα+h12tbH+wk+1-tbH+12Htd+wk+1-τ0dθ2+Umaxtb+τ0-tbdθ1Ciw+wk+1-τ0d1-θ2+Umaxtb+τ0-tbd1-θ1Cpw=Ea13τ0,tr0where Ea13(τ0,tr0) denotes E(Ca∣τ=τ0,tr=tr0) when wk<τ0<wk+1,tr0<td.From Figure13, it can be obtained that(54)ETc∣τ=τ0,tr=tr0=wk+1+td=Ec13where Ec13 denotes E(Tc∣τ=τ0,tr=tr0) when wk<τ0<wk+1,tr0<td.


#### 4.1.14. Scenario 14:wk<τ<wk+1,tr>td

This scenario can be rewritten asτ=τ0,tr=tr0 and wk<τ0<wk+1,tr0>td. Thus the behavior of inventory can be described as Figure 14. In contrast to Scenario 13, shortage occurs since tr0>td. Hence(55)ES∣τ=τ0,tr=tr0=Cstr0-tdd.Other costs are the same as Scenario 13. Therefore(56)Ea14τ0,tr0=Ea13τ0,tr0+ES∣τ=τ0,tr=tr0=Ea13τ0,tr0+Cstr0-tddwhere Ea14(τ0,tr0) denotes E(Ca∣τ=τ0,tr=tr0) when wk<τ0<wk+1,tr0>td. From Figure 14, it can be obtained that(57)ETc∣τ=τ0,tr=tr0=wk+1+tr0=Ec14tr0where Ec14(tr0) denotes E(Tc∣τ=τ0,tr=tr0) when wk<τ0<wk+1,tr0>td.Figure 14
Behavior of inventory (Scenario 14).


#### 4.1.15. Scenario 15:τ>wk+1,tp<td

This scenario can be rewritten asτ=τ0,tp=tp0 and τ0>wk+1,tp0<td. Thus the behavior of inventory can be described as Figure 15. In this scenario, the process keeps in the in-control state during the (k+1) sampling intervals. Hence the preventive maintenance (PM) is implemented at the end of the (k+1)th sampling interval (without sampling). Because tp0<td, shortage does not occur.Figure 15
Behavior of inventory (Scenario 15).The conditional expectations of all kinds of costs given thatτ=τ0,tp=tp0 can be calculated using the same analogy as Scenario 1. And they are displayed in Tables 1 and 2 of the Annex. Therefore(58)ECa∣τ=τ0,tp=tp0=Cc∫0wk+1λXtdt+K∫0wk+1λXtdt+1+Cptp0+ks+Cykα+h12tbH+wk+1-tbH+12Htd+Umaxtb+wk+1-tbdθ1Ciw+Umaxtb+wk+1-tbd1-θ1Cpw=Ea15tp0where Ea15(tp0) denotes E(Ca∣τ=τ0,tp=tp0) when τ0>wk+1,tp0<td.From Figure15, it can be obtained that(59)ETc∣τ=τ0,tp=tp0=wk+1+td=Ec15where Ec15 denotes E(Tc∣τ=τ0,tp=tp0) when τ0>wk+1,tp0<td.


#### 4.1.16. Scenario 16:τ>wk+1,tp>td

This scenario can be rewritten asτ=τ0,tp=tp0 and τ0>wk+1,tp0>td. Thus the behavior of inventory can be described as Figure 16. In contrast to Scenario 15, shortage occurs since tp0>td. Hence(60)ES∣τ=τ0,tp=tp0=Cstp0-tdd.Other costs are the same as Scenario 15. Therefore(61)Ea16tp0=Ea15tp0+ES∣τ=τ0,y=y0,tr=tr0=Ea15tp0+Cstp0-tddwhere Ea16(tp0) denotes E(Ca∣τ=τ0,tp=tp0) when τ0>wk+1,tp0>td. From Figure 16, it can be obtained that(62)ETc∣τ=τ0,tp=tp0=wk+1+tp0=Ec16tp0where Ec16(tp0) denotes E(Tc∣τ=τ0,tp=tp0) when τ0>wk+1,tp0>td.Figure 16
Behavior of inventory (Scenario 16).


### 4.2. Objective Function and Constraints

By applying the renewal theory, the total expected cost per unit timeETC can be obtained as follows:(63)ETC=ETCH,k,h1,L,n=ECaETc.Thus the nonlinear optimization problem can be obtained:(64)minETCs.t:(65)h1>tb,(66)H≤Hm,H∈N+,(67)k≤km,k∈N+,(68)0<L<Lm,(69)1≤n≤nm.The constraints are explained as follows. Inequality (65): in order to avoid overmuch sampling cost, the first sampling should be conducted after the hedging safety stock is built up. Inequality (66): in practice, the available inventory is finite. Inequality (67): in order to avoid excessive samplings, the value of k must be limited. Inequality (68): in order to decrease the error of the control chart, the value of L must be limited. Inequality (69): in practice, because of economic reasons, the sample size must be less than a definite number.


## 4.1. The 16 Scenarios

The period between the two successive starts of the good-as-new equipment is described as a production cycle. Thus, according to the policies proposed in Section3.2, 16 scenarios may occur in a production cycle. This section contributes to describing these scenarios in detail and formulating the conditional expectations of Ca and Tc for different scenarios.


### 4.1.1. Scenario 1:0<τ<tb, 1≤y≤k, 0<tr<td

This scenario can be rewritten asτ=τ0,y=y0,tr=tr0 and 0<τ0<tb, 1≤y0≤k, 0<tr0<td. Then the behavior of inventory in this scenario can be described as Figure 1. In this scenario, the process shifts to out-of-control state before the hedging safety stock H is built. Due to the type II error, the alert signal is not sent until the y0th sampling. Then the PRDM is conducted for tr0. Because 0<tr0<td, shortage does not occur. The conditional expectations of Ca and Tc given that τ=τ0,y=y0,tr=tr0, that is, E(Ca∣τ=τ0,y=y0,tr=tr0) and E(Tc∣τ=τ0,y=y0,tr=tr0), when 0<τ0<tb, 1≤y0≤k, 0<tr0<td will be formulated.Figure 1
Behavior of inventory (Scenario 1).BecauseCa=C+Ks+P+R+B+A+I+F+S, E(Ca∣τ=τ0,y=y0,tr=tr0) can be obtained through calculating the conditional expectations of C,Ks,P,R,B,A,I,F, and S. The CM is minimal repair and the repair time is negligible; hence the machine fails according to an NHPP with rate λX(t). Therefore, the mean number of failures within the production period wy0 is ∫0wy0λX(t)dt. Consequently, the conditional expectations of C are(12)EC∣τ=τ0,y=y0,tr=tr0=Cc∫0wy0λXtdt.The conditional expectations of Ks are(13)EKs∣τ=τ0,y=y0,tr=tr0=K∫0wy0λXtdt+1.The PM is not conducted; hence(14)EP∣τ=τ0,y=y0,tr=tr0=0.The conditional expectations of R are(15)ER∣τ=τ0,y=y0,tr=tr0=Crtr0.The conditional expectations of B are(16)EB∣τ=τ0,y=y0,tr=tr0=y0s.The conditional expectations of A are(17)EA∣τ=τ0,y=y0,tr=tr0=Cy.The conditional expectations of I are(18)EI∣τ=τ0,y=y0,tr=tr0=h12tbH+wy0-tbH+12Htd.The conditional expectations of Ni are(19)ENi∣τ=τ0,y=y0,tr=tr0=tb-τ0Umax+wy0-tbdθ2+Umaxτ0θ1.The conditional expectations of Np are(20)ENp∣τ=τ0,y=y0,tr=tr0=tb-τ0Umax+wy0-tbd1-θ2+Umaxτ01-θ1.Hence,(21)EF∣τ=τ0,y=y0,tr=tr0=tb-τ0Umax+wy0-tbdθ2+Umaxτ0θ1Ciw+tb-τ0Umax+wy0-tbd1-θ2+Umaxτ01-θ1Cpw.Shortage does not occur; hence(22)ES∣τ=τ0,y=y0,tr=tr0=0.Above all, (23)ECa∣τ=τ0,y=y0,tr=tr0=Cc∫0wy0λXtdt+K∫0wy0λXtdt+1+Crtr0+y0s+Cy+h12tbH+wy0-tbH+12Htd+tb-τ0Umax+wy0-tbdθ2+Umaxτ0θ1Ciw+tb-τ0Umax+wy0-tbd1-θ2+Umaxτ01-θ1Cpw=Ea1τ0,y0,tr0where Ea1(τ0,y0,tr0) denotes E(Ca∣τ=τ0,y=y0,tr=tr0) when 0<τ0<tb, 1≤y0≤k, 0<tr0<td.From Figure1, it can be obtained that(24)ETc∣τ=τ0,y=y0,tr=tr0=wy0+td=Ec1y0where Ec1(y0) denotes E(Tc∣τ=τ0,y=y0,tr=tr0) when 0<τ0<tb, 1≤y0≤k, 0<tr0<td.


### 4.1.2. Scenario 2:0<τ<tb, 1≤y≤k, tr>td

This scenario can be rewritten asτ=τ0,y=y0,tr=tr0 and 0<τ0<tb, 1≤y0≤k, tr0>td. Thus the behavior of inventory in this scenario can be described as Figure 2. In contrast to Scenario 1, shortage occurs since tr0>td. Hence(25)ES∣τ=τ0,y=y0,tr=tr0=Cstr0-tdd.Other costs are the same as Scenario 1. Therefore(26)ECa∣τ=τ0,y=y0,tr=tr0=Ea1τ0,y0,tr0+Cstr0-tdd=Ea2τ0,y0,tr0where Ea2(τ0,y0,tr0) denotes E(Ca∣τ=τ0,y=y0,tr=tr0) when 0<τ0<tb, 1≤y0≤k, tr0>td. From Figure 2, it can be obtained that(27)ETc∣τ=τ0,y=y0,tr=tr0=wy0+tr0=Ec2y0,tr0where Ec2(y0,tr0) denotes E(Tc∣τ=τ0,y=y0,tr=tr0) when 0<τ0<tb, 1≤y0≤k, tr0>td.Figure 2
Behavior of inventory (Scenario 2).


### 4.1.3. Scenario 3:0<τ<tb, y≥k+1, 0<tr<td

This scenario can be rewritten asτ=τ0,y=y0,tr=tr0 and 0<τ0<tb, y0≥k+1, 0<tr0<td. Thus the behavior of inventory in this scenario can be described as Figure 3. The alert signal is not sent in the first k samples. Hence the preventive maintenance (PM) is implemented at the end of the (k+1)th sampling interval (without sampling). Then the assignable cause is detected during the PM. Thus, the PM is replaced by the PRDM. Because 0<tr0<td, shortage does not occur.Figure 3
Behavior of inventory (Scenario 3).The conditional expectations of all kinds of costs given thatτ=τ0,y=y0,tr=tr0 can be calculated using the same analogy as Scenario 1. And they are displayed in Tables 1 and 2 of the Annex (Supplementary 1). Hence(28)ECa∣τ=τ0,y=y0,tr=tr0=Cc∫0wk+1λXtdt+K∫0wk+1λXtdt+1+Crtr0+ks+h12tbH+wk+1-tbH+12Htd+tb-τ0Umax+wk+1-tbdθ2+Umaxτ0θ1Ciw+tb-τ0Umax+wk+1-tbd1-θ2+Umaxτ01-θ1Cpw=Ea3τ0,tr0where Ea3(τ0,tr0) denotes E(Ca∣τ=τ0,y=y0,tr=tr0) when 0<τ0<tb, y0≥k+1, 0<tr0<td.From Figure3, it can be obtained that(29)ETc∣τ=τ0,y=y0,tr=tr0=wk+1+td=Ec3where Ec3 denotes E(Tc∣τ=τ0,y=y0,tr=tr0) when 0<τ0<tb, y0≥k+1, 0<tr0<td.


### 4.1.4. Scenario 4:0<τ<tb, y≥k+1, tr>td

This scenario can be rewritten asτ=τ0,y=y0,tr=tr0 and 0<τ0<tb, y0≥k+1, tr0>td. Thus the behavior of inventory in this scenario can be described as Figure 4. In contrast to Scenario 3, shortage occurs since tr0>td. Hence(30)ES∣τ=τ0,y=y0,tr=tr0=Cstr0-tdd.Other costs are the same as Scenario 3. Therefore(31)Ea4τ0,tr0=Ea3τ0,tr0+ES∣τ=τ0,y=y0,tr=tr0=Ea3τ0,tr0+Cstr0-tddwhere Ea4(τ0,tr0) denotes E(Ca∣τ=τ0,y=y0,tr=tr0) when 0<τ0<tb, y0≥k+1, tr0>td. From Figure 4, it can be obtained that(32)ETc∣τ=τ0,y=y0,tr=tr0=wk+1+tr0=Ec4tr0where Ec4(tr0) denotes E(Tc∣τ=τ0,y=y0,tr=tr0) when 0<τ0<tb, y0≥k+1, tr0>td.Figure 4
Behavior of inventory (Scenario 4).


### 4.1.5. Scenario 5:tb<τ<w1, 1≤y≤k, 0<tr<td

This scenario can be rewritten asτ=τ0,y=y0,tr=tr0 and tb<τ0<w1, 1≤y0≤k, 0<tr0<td. Thus the behavior of inventory can be described as Figure 5. In this scenario, the process shifts to out-of-control state after the hedging safety stock H is built and before the first sampling. Due to the type II error, the alert signal is not sent until the y0th sampling. Then the PRDM is conducted for tr0. Because 0<tr0<td, shortage does not occur.Figure 5
Behavior of inventory (Scenario 5).The conditional expectations of all kinds of costs given thatτ=τ0,y=y0,tr=tr0 can be calculated using the same analogy as Scenario 1. And they are displayed in Tables 1 and 2 of the Annex. Therefore (33)ECa∣τ=τ0,y=y0,tr=tr0=Cc∫0wy0λXtdt+K∫0wy0λXtdt+1+Crtr0+y0s+Cy+h12tbH+wy0-tbH+12Htd+wy0-τ0dθ2+Umaxtb+τ0-tbdθ1Ciw+wy0-τ0d1-θ2+Umaxtb+τ0-tbd1-θ1Cpw=Ea5τ0,y0,tr0where Ea5(τ0,y0,tr0) denotes E(Ca∣τ=τ0,y=y0,tr=tr0) when tb<τ0<w1, 1≤y0≤k, 0<tr0<td.From Figure5, it can be obtained that(34)ETc∣τ=τ0,y=y0,tr=tr0=wy0+td=Ec5y0where Ec5(y0) denotes E(Tc∣τ=τ0,y=y0,tr=tr0) when tb<τ0<w1, 1≤y0≤k, 0<tr0<td.


### 4.1.6. Scenario 6:tb<τ<w1, 1≤y≤k, tr>td

This scenario can be rewritten asτ=τ0,y=y0,tr=tr0 and tb<τ0<w1, 1≤y0≤k, tr0>td. Thus the behavior of inventory can be described as Figure 6. In contrast to Scenario 5, shortage occurs since tr0>td. Hence(35)ES∣τ=τ0,y=y0,tr=tr0=Cstr0-tdd.Other costs are the same as Scenario 5. Therefore(36)Ea6τ0,y0,tr0=Ea5τ0,y0,tr0+ES∣τ=τ0,y=y0,tr=tr0=Ea5τ0,y0,tr0+Cstr0-tddwhere Ea6(τ0,y0,tr0) denotes E(Ca∣τ=τ0,y=y0,tr=tr0) when tb<τ0<w1, 1≤y0≤k, tr0>td. From Figure 6, it can be obtained that(37)ETc∣τ=τ0,y=y0,tr=tr0=wy0+tr0=Ec6y0,tr0where Ec6(y0,tr0) denotes ETc∣τ=τ0,y=y0,tr=tr0 when tb<τ0<w1, 1≤y0≤k, tr0>td.Figure 6
Behavior of inventory (Scenario 6).


### 4.1.7. Scenario 7:tb<τ<w1, y≥k+1, 0<tr<td

This scenario can be rewritten asτ=τ0,y=y0,tr=tr0 and tb<τ0<w1, y0≥k+1, 0<tr0<td. Thus the behavior of inventory can be described as Figure 7. The alert signal is not sent in the first k samples. Hence the preventive maintenance (PM) is implemented at the end of the (k+1)th sampling interval (without sampling). Then the assignable cause is detected during the PM. Thus the PM is replaced by PRDM. Because 0<tr0<td, shortage does not occur.Figure 7
Behavior of inventory (Scenario 7).The conditional expectations of all kinds of costs given thatτ=τ0,y=y0,tr=tr0 can be calculated using the same analogy as Scenario 1. And they are displayed in Tables 1 and 2 of the Annex. Therefore(38)ECa∣τ=τ0,y=y0,tr=tr0=Cc∫0wk+1λXtdt+K∫0wk+1λXtdt+1+Crtr0+ks+h12tbH+wk+1-tbH+12Htd+wk+1-τ0dθ2+Umaxtb+τ0-tbdθ1Ciw+wk+1-τ0d1-θ2+Umaxtb+τ0-tbd1-θ1Cpw=Ea7τ0,tr0where Ea7(τ0,tr0) denotes E(Ca∣τ=τ0,y=y0,tr=tr0) when tb<τ0<w1, y0≥k+1, 0<tr0<td.From Figure7, it can be obtained that(39)ETc∣τ=τ0,y=y0,tr=tr0=wk+1+td=Ec7where Ec7 denotes E(Tc∣τ=τ0,y=y0,tr=tr0) when tb<τ0<w1, y0≥k+1, 0<tr0<td.


### 4.1.8. Scenario 8:tb<τ<w1, y≥k+1, tr>td

This scenario can be rewritten asτ=τ0,y=y0,tr=tr0 and tb<τ0<w1, y0≥k+1, tr0>td. Thus the behavior of inventory can be described as Figure 8. In contrast to Scenario 7, shortage occurs since tr0>td. Hence(40)ES∣τ=τ0,y=y0,tr=tr0=Cstr0-tdd.Other costs are the same as Scenario 7. Therefore(41)Ea8τ0,tr0=Ea7τ0,tr0+ES∣τ=τ0,y=y0,tr=tr0=Ea7τ0,tr0+Cstr0-tddwhere Ea8(τ0,tr0) denotes E(Ca∣τ=τ0,y=y0,tr=tr0) when tb<τ0<w1, y0≥k+1, tr0>td. From Figure 8, it can be obtained that(42)ETc∣τ=τ0,y=y0,tr=tr0=wk+1+tr0=Ec8tr0where Ec8 denotes E(Tc∣τ=τ0,y=y0,tr=tr0) when tb<τ0<w1, y0≥k+1, tr0>td.Figure 8
Behavior of inventory (Scenario 8).


### 4.1.9. Scenario 9:wj<τ<wj+1(1≤j≤k-1), 1≤y≤k-j, tr<td

This scenario can be rewritten asτ=τ0,y=y0,tr=tr0 and wj<τ0<wj+1(1≤j≤k-1),1≤y0≤k-j,tr0<td. Thus the behavior of inventory can be described as Figure 9. In this scenario, the process shifts to out-of-control state in the (j+1)th sampling interval. After the shift, due to the type II error, the alert signal is not sent until the (j+y0)th sampling, while the alert signal may be sent in the first j+1 samplings due to the type I error. After the (j+y0)th sampling, the PRDM is conducted for tr0. Because 0<tr0<td, shortage does not occur.Figure 9
Behavior of inventory (Scenario 9).The conditional expectations of all kinds of costs given thatτ=τ0,y=y0,tr=tr0 can be calculated using the same analogy as Scenario 1. And they are displayed in Tables 1 and 2 of the Annex. Therefore(43)ECa∣τ=τ0,y=y0,tr=tr0=Cc∫0wj+y0λXtdt+K∫0wj+y0λXtdt+1+Crtr0+j+y0s+Cyjα+1+h12tbH+wj+y0-tbH+12Htd+wj+y0-τ0dθ2+Umaxtb+τ0-tbdθ1Ciw+wj+y0-τ0d1-θ2+Umaxtb+τ0-tbd1-θ1Cpw=Ea9τ0,y0,tr0where Ea9(τ0,y0,tr0) denotes E(Ca∣τ=τ0,y=y0,tr=tr0) when wj<τ0<wj+1(1≤j≤k-1),1≤y0≤k-j,tr0<td.From Figure9, it can be obtained that(44)ETc∣τ=τ0,y=y0,tr=tr0=wj+y0+td=Ec9y0where Ec9(y0) denotes E(Tc∣τ=τ0,y=y0,tr=tr0) when wj<τ0<wj+1(1≤j≤k-1),1≤y0≤k-j,tr0<td.


### 4.1.10. Scenario 10:wj<τ<wj+1(1≤j≤k-1), 1≤y≤k-j, tr>td

This scenario can be rewritten asτ=τ0,y=y0,tr=tr0 and wj<τ0<wj+1(1≤j≤k-1),1≤y0≤k-j,tr0>td. Thus the behavior of inventory can be described as Figure 10. In contrast to Scenario 9, shortage occurs since tr0>td. Hence(45)ES∣τ=τ0,y=y0,tr=tr0=Cstr0-tdd.Other costs are the same as Scenario 9. Therefore(46)Ea10τ0,y0,tr0=Ea9τ0,y0,tr0+ES∣τ=τ0,y=y0,tr=tr0=Ea9τ0,y0,tr0+Cstr0-tddwhere Ea10(τ0,y0,tr0) denotes E(Ca∣τ=τ0,y=y0,tr=tr0) when wj<τ0<wj+1(1≤j≤k-1),1≤y0≤k-j,tr0>td. From Figure 10, it can be obtained that(47)ETc∣τ=τ0,y=y0,tr=tr0=wj+y0+tr0=Ec10tr0where Ec10(tr0) denotes E(Tc∣τ=τ0,y=y0,tr=tr0) when wj<τ0<wj+1(1≤j≤k-1),1≤y0≤k-j,tr0>td.Figure 10
Behavior of inventory (Scenario 10).


### 4.1.11. Scenario 11:wj<τ<wj+1(1≤j≤k-1), y>k-j, tr<td

This scenario can be rewritten asτ=τ0,y=y0,tr=tr0 and wj<τ0<wj+1(1≤j≤k-1),y0>k-j,tr0<td. Thus the behavior of inventory can be described as Figure 11. No PRDM is conducted after k samplings. Hence the preventive maintenance (PM) is implemented at the end of the (k+1)th sampling interval (without sampling). Then the assignable cause is detected during the PM. Thus the PM is replaced by the PRDM. Because 0<tr0<td, shortage does not occur.Figure 11
Behavior of inventory (Scenario 11).The conditional expectations of all kinds of costs given thatτ=τ0,y=y0,tr=tr0 can be calculated using the same analogy as Scenario 1. And they are displayed in Tables 1 and 2 of the Annex. Therefore (48)ECa∣τ=τ0,y=y0,tr=tr0=Cc∫0wk+1λXtdt+K∫0wk+1λXtdt+1+Crtr0+ks+Cyjα+h12tbH+wk+1-tbH+12Htd+wk+1-τ0dθ2+Umaxtb+τ0-tbdθ1Ciw+wk+1-τ0d1-θ2+Umaxtb+τ0-tbd1-θ1Cpw=Ea11τ0,tr0where Ea11(τ0,tr0) denotes E(Ca∣τ=τ0,y=y0,tr=tr0) when wj<τ0<wj+1(1≤j≤k-1),y0>k-j,tr0<td.From Figure11, it can be obtained that(49)ETc∣τ=τ0,y=y0,tr=tr0=wk+1+td=Ec11where Ec11 denotes E(Tc∣τ=τ0,y=y0,tr=tr0) when wj<τ0<wj+1(1≤j≤k-1),y0>k-j,tr0<td.


### 4.1.12. Scenario 12:wj<τ<wj+1(1≤j≤k-1), y>k-j, tr>td

This scenario can be rewritten asτ=τ0,y=y0,tr=tr0 and wj<τ0<wj+1(1≤j≤k-1),y0>k-j,tr0>td. Thus the behavior of inventory can be described as Figure 12. In contrast to Scenario 11, shortage occurs since tr0>td. Hence(50)ES∣τ=τ0,y=y0,tr=tr0=Cstr0-tdd.Other costs are the same as Scenario 11. Therefore(51)Ea12τ0,tr0=Ea11τ0,tr0+ES∣τ=τ0,y=y0,tr=tr0=Ea11τ0,tr0+Cstr0-tddwhere Ea12(τ0,tr0) denotes E(Ca∣τ=τ0,y=y0,tr=tr0) when wj<τ0<wj+1(1≤j≤k-1),y0>k-j,tr0>td. From Figure 12, it can be obtained that(52)ETc∣τ=τ0,y=y0,tr=tr0=wk+1+tr0=Ec12tr0where Ec12(tr0) denotes E(Tc∣τ=τ0,y=y0,tr=tr0) when wj<τ0<wj+1(1≤j≤k-1),y0>k-j,tr0>td.Figure 12
Behavior of inventory (Scenario 12).


### 4.1.13. Scenario 13:wk<τ<wk+1,tr<td

This scenario can be rewritten asτ=τ0,tr=tr0 and wk<τ0<wk+1,tr0<td. Thus the behavior of inventory can be described as Figure 13. In this scenario, the process shifts to out-of-control state in the (k+1)th sampling interval. No PRDM is conducted after k samplings. Hence the preventive maintenance (PM) is implemented at the end of the (k+1)th sampling interval (without sampling). Then the assignable cause is detected during the PM. Thus the PM is replaced by PRDM. Because 0<tr0<td, shortage does not occur.Figure 13
Behavior of inventory (Scenario 13).The conditional expectations of all kinds of costs given thatτ=τ0,tr=tr0 can be calculated using the same analogy as Scenario 1. And they are displayed in Tables 1 and 2 of the Annex. Therefore (53)ECa∣τ=τ0,tr=tr0=Cc∫0wk+1λXtdt+K∫0wk+1λXtdt+1+Crtr0+ks+Cykα+h12tbH+wk+1-tbH+12Htd+wk+1-τ0dθ2+Umaxtb+τ0-tbdθ1Ciw+wk+1-τ0d1-θ2+Umaxtb+τ0-tbd1-θ1Cpw=Ea13τ0,tr0where Ea13(τ0,tr0) denotes E(Ca∣τ=τ0,tr=tr0) when wk<τ0<wk+1,tr0<td.From Figure13, it can be obtained that(54)ETc∣τ=τ0,tr=tr0=wk+1+td=Ec13where Ec13 denotes E(Tc∣τ=τ0,tr=tr0) when wk<τ0<wk+1,tr0<td.


### 4.1.14. Scenario 14:wk<τ<wk+1,tr>td

This scenario can be rewritten asτ=τ0,tr=tr0 and wk<τ0<wk+1,tr0>td. Thus the behavior of inventory can be described as Figure 14. In contrast to Scenario 13, shortage occurs since tr0>td. Hence(55)ES∣τ=τ0,tr=tr0=Cstr0-tdd.Other costs are the same as Scenario 13. Therefore(56)Ea14τ0,tr0=Ea13τ0,tr0+ES∣τ=τ0,tr=tr0=Ea13τ0,tr0+Cstr0-tddwhere Ea14(τ0,tr0) denotes E(Ca∣τ=τ0,tr=tr0) when wk<τ0<wk+1,tr0>td. From Figure 14, it can be obtained that(57)ETc∣τ=τ0,tr=tr0=wk+1+tr0=Ec14tr0where Ec14(tr0) denotes E(Tc∣τ=τ0,tr=tr0) when wk<τ0<wk+1,tr0>td.Figure 14
Behavior of inventory (Scenario 14).


### 4.1.15. Scenario 15:τ>wk+1,tp<td

This scenario can be rewritten asτ=τ0,tp=tp0 and τ0>wk+1,tp0<td. Thus the behavior of inventory can be described as Figure 15. In this scenario, the process keeps in the in-control state during the (k+1) sampling intervals. Hence the preventive maintenance (PM) is implemented at the end of the (k+1)th sampling interval (without sampling). Because tp0<td, shortage does not occur.Figure 15
Behavior of inventory (Scenario 15).The conditional expectations of all kinds of costs given thatτ=τ0,tp=tp0 can be calculated using the same analogy as Scenario 1. And they are displayed in Tables 1 and 2 of the Annex. Therefore(58)ECa∣τ=τ0,tp=tp0=Cc∫0wk+1λXtdt+K∫0wk+1λXtdt+1+Cptp0+ks+Cykα+h12tbH+wk+1-tbH+12Htd+Umaxtb+wk+1-tbdθ1Ciw+Umaxtb+wk+1-tbd1-θ1Cpw=Ea15tp0where Ea15(tp0) denotes E(Ca∣τ=τ0,tp=tp0) when τ0>wk+1,tp0<td.From Figure15, it can be obtained that(59)ETc∣τ=τ0,tp=tp0=wk+1+td=Ec15where Ec15 denotes E(Tc∣τ=τ0,tp=tp0) when τ0>wk+1,tp0<td.


### 4.1.16. Scenario 16:τ>wk+1,tp>td

This scenario can be rewritten asτ=τ0,tp=tp0 and τ0>wk+1,tp0>td. Thus the behavior of inventory can be described as Figure 16. In contrast to Scenario 15, shortage occurs since tp0>td. Hence(60)ES∣τ=τ0,tp=tp0=Cstp0-tdd.Other costs are the same as Scenario 15. Therefore(61)Ea16tp0=Ea15tp0+ES∣τ=τ0,y=y0,tr=tr0=Ea15tp0+Cstp0-tddwhere Ea16(tp0) denotes E(Ca∣τ=τ0,tp=tp0) when τ0>wk+1,tp0>td. From Figure 16, it can be obtained that(62)ETc∣τ=τ0,tp=tp0=wk+1+tp0=Ec16tp0where Ec16(tp0) denotes E(Tc∣τ=τ0,tp=tp0) when τ0>wk+1,tp0>td.Figure 16
Behavior of inventory (Scenario 16).


## 4.1.1. Scenario 1:0<τ<tb, 1≤y≤k, 0<tr<td

This scenario can be rewritten asτ=τ0,y=y0,tr=tr0 and 0<τ0<tb, 1≤y0≤k, 0<tr0<td. Then the behavior of inventory in this scenario can be described as Figure 1. In this scenario, the process shifts to out-of-control state before the hedging safety stock H is built. Due to the type II error, the alert signal is not sent until the y0th sampling. Then the PRDM is conducted for tr0. Because 0<tr0<td, shortage does not occur. The conditional expectations of Ca and Tc given that τ=τ0,y=y0,tr=tr0, that is, E(Ca∣τ=τ0,y=y0,tr=tr0) and E(Tc∣τ=τ0,y=y0,tr=tr0), when 0<τ0<tb, 1≤y0≤k, 0<tr0<td will be formulated.Figure 1
Behavior of inventory (Scenario 1).BecauseCa=C+Ks+P+R+B+A+I+F+S, E(Ca∣τ=τ0,y=y0,tr=tr0) can be obtained through calculating the conditional expectations of C,Ks,P,R,B,A,I,F, and S. The CM is minimal repair and the repair time is negligible; hence the machine fails according to an NHPP with rate λX(t). Therefore, the mean number of failures within the production period wy0 is ∫0wy0λX(t)dt. Consequently, the conditional expectations of C are(12)EC∣τ=τ0,y=y0,tr=tr0=Cc∫0wy0λXtdt.The conditional expectations of Ks are(13)EKs∣τ=τ0,y=y0,tr=tr0=K∫0wy0λXtdt+1.The PM is not conducted; hence(14)EP∣τ=τ0,y=y0,tr=tr0=0.The conditional expectations of R are(15)ER∣τ=τ0,y=y0,tr=tr0=Crtr0.The conditional expectations of B are(16)EB∣τ=τ0,y=y0,tr=tr0=y0s.The conditional expectations of A are(17)EA∣τ=τ0,y=y0,tr=tr0=Cy.The conditional expectations of I are(18)EI∣τ=τ0,y=y0,tr=tr0=h12tbH+wy0-tbH+12Htd.The conditional expectations of Ni are(19)ENi∣τ=τ0,y=y0,tr=tr0=tb-τ0Umax+wy0-tbdθ2+Umaxτ0θ1.The conditional expectations of Np are(20)ENp∣τ=τ0,y=y0,tr=tr0=tb-τ0Umax+wy0-tbd1-θ2+Umaxτ01-θ1.Hence,(21)EF∣τ=τ0,y=y0,tr=tr0=tb-τ0Umax+wy0-tbdθ2+Umaxτ0θ1Ciw+tb-τ0Umax+wy0-tbd1-θ2+Umaxτ01-θ1Cpw.Shortage does not occur; hence(22)ES∣τ=τ0,y=y0,tr=tr0=0.Above all, (23)ECa∣τ=τ0,y=y0,tr=tr0=Cc∫0wy0λXtdt+K∫0wy0λXtdt+1+Crtr0+y0s+Cy+h12tbH+wy0-tbH+12Htd+tb-τ0Umax+wy0-tbdθ2+Umaxτ0θ1Ciw+tb-τ0Umax+wy0-tbd1-θ2+Umaxτ01-θ1Cpw=Ea1τ0,y0,tr0where Ea1(τ0,y0,tr0) denotes E(Ca∣τ=τ0,y=y0,tr=tr0) when 0<τ0<tb, 1≤y0≤k, 0<tr0<td.From Figure1, it can be obtained that(24)ETc∣τ=τ0,y=y0,tr=tr0=wy0+td=Ec1y0where Ec1(y0) denotes E(Tc∣τ=τ0,y=y0,tr=tr0) when 0<τ0<tb, 1≤y0≤k, 0<tr0<td.


## 4.1.2. Scenario 2:0<τ<tb, 1≤y≤k, tr>td

This scenario can be rewritten asτ=τ0,y=y0,tr=tr0 and 0<τ0<tb, 1≤y0≤k, tr0>td. Thus the behavior of inventory in this scenario can be described as Figure 2. In contrast to Scenario 1, shortage occurs since tr0>td. Hence(25)ES∣τ=τ0,y=y0,tr=tr0=Cstr0-tdd.Other costs are the same as Scenario 1. Therefore(26)ECa∣τ=τ0,y=y0,tr=tr0=Ea1τ0,y0,tr0+Cstr0-tdd=Ea2τ0,y0,tr0where Ea2(τ0,y0,tr0) denotes E(Ca∣τ=τ0,y=y0,tr=tr0) when 0<τ0<tb, 1≤y0≤k, tr0>td. From Figure 2, it can be obtained that(27)ETc∣τ=τ0,y=y0,tr=tr0=wy0+tr0=Ec2y0,tr0where Ec2(y0,tr0) denotes E(Tc∣τ=τ0,y=y0,tr=tr0) when 0<τ0<tb, 1≤y0≤k, tr0>td.Figure 2
Behavior of inventory (Scenario 2).


## 4.1.3. Scenario 3:0<τ<tb, y≥k+1, 0<tr<td

This scenario can be rewritten asτ=τ0,y=y0,tr=tr0 and 0<τ0<tb, y0≥k+1, 0<tr0<td. Thus the behavior of inventory in this scenario can be described as Figure 3. The alert signal is not sent in the first k samples. Hence the preventive maintenance (PM) is implemented at the end of the (k+1)th sampling interval (without sampling). Then the assignable cause is detected during the PM. Thus, the PM is replaced by the PRDM. Because 0<tr0<td, shortage does not occur.Figure 3
Behavior of inventory (Scenario 3).The conditional expectations of all kinds of costs given thatτ=τ0,y=y0,tr=tr0 can be calculated using the same analogy as Scenario 1. And they are displayed in Tables 1 and 2 of the Annex (Supplementary 1). Hence(28)ECa∣τ=τ0,y=y0,tr=tr0=Cc∫0wk+1λXtdt+K∫0wk+1λXtdt+1+Crtr0+ks+h12tbH+wk+1-tbH+12Htd+tb-τ0Umax+wk+1-tbdθ2+Umaxτ0θ1Ciw+tb-τ0Umax+wk+1-tbd1-θ2+Umaxτ01-θ1Cpw=Ea3τ0,tr0where Ea3(τ0,tr0) denotes E(Ca∣τ=τ0,y=y0,tr=tr0) when 0<τ0<tb, y0≥k+1, 0<tr0<td.From Figure3, it can be obtained that(29)ETc∣τ=τ0,y=y0,tr=tr0=wk+1+td=Ec3where Ec3 denotes E(Tc∣τ=τ0,y=y0,tr=tr0) when 0<τ0<tb, y0≥k+1, 0<tr0<td.


## 4.1.4. Scenario 4:0<τ<tb, y≥k+1, tr>td

This scenario can be rewritten asτ=τ0,y=y0,tr=tr0 and 0<τ0<tb, y0≥k+1, tr0>td. Thus the behavior of inventory in this scenario can be described as Figure 4. In contrast to Scenario 3, shortage occurs since tr0>td. Hence(30)ES∣τ=τ0,y=y0,tr=tr0=Cstr0-tdd.Other costs are the same as Scenario 3. Therefore(31)Ea4τ0,tr0=Ea3τ0,tr0+ES∣τ=τ0,y=y0,tr=tr0=Ea3τ0,tr0+Cstr0-tddwhere Ea4(τ0,tr0) denotes E(Ca∣τ=τ0,y=y0,tr=tr0) when 0<τ0<tb, y0≥k+1, tr0>td. From Figure 4, it can be obtained that(32)ETc∣τ=τ0,y=y0,tr=tr0=wk+1+tr0=Ec4tr0where Ec4(tr0) denotes E(Tc∣τ=τ0,y=y0,tr=tr0) when 0<τ0<tb, y0≥k+1, tr0>td.Figure 4
Behavior of inventory (Scenario 4).


## 4.1.5. Scenario 5:tb<τ<w1, 1≤y≤k, 0<tr<td

This scenario can be rewritten asτ=τ0,y=y0,tr=tr0 and tb<τ0<w1, 1≤y0≤k, 0<tr0<td. Thus the behavior of inventory can be described as Figure 5. In this scenario, the process shifts to out-of-control state after the hedging safety stock H is built and before the first sampling. Due to the type II error, the alert signal is not sent until the y0th sampling. Then the PRDM is conducted for tr0. Because 0<tr0<td, shortage does not occur.Figure 5
Behavior of inventory (Scenario 5).The conditional expectations of all kinds of costs given thatτ=τ0,y=y0,tr=tr0 can be calculated using the same analogy as Scenario 1. And they are displayed in Tables 1 and 2 of the Annex. Therefore (33)ECa∣τ=τ0,y=y0,tr=tr0=Cc∫0wy0λXtdt+K∫0wy0λXtdt+1+Crtr0+y0s+Cy+h12tbH+wy0-tbH+12Htd+wy0-τ0dθ2+Umaxtb+τ0-tbdθ1Ciw+wy0-τ0d1-θ2+Umaxtb+τ0-tbd1-θ1Cpw=Ea5τ0,y0,tr0where Ea5(τ0,y0,tr0) denotes E(Ca∣τ=τ0,y=y0,tr=tr0) when tb<τ0<w1, 1≤y0≤k, 0<tr0<td.From Figure5, it can be obtained that(34)ETc∣τ=τ0,y=y0,tr=tr0=wy0+td=Ec5y0where Ec5(y0) denotes E(Tc∣τ=τ0,y=y0,tr=tr0) when tb<τ0<w1, 1≤y0≤k, 0<tr0<td.


## 4.1.6. Scenario 6:tb<τ<w1, 1≤y≤k, tr>td

This scenario can be rewritten asτ=τ0,y=y0,tr=tr0 and tb<τ0<w1, 1≤y0≤k, tr0>td. Thus the behavior of inventory can be described as Figure 6. In contrast to Scenario 5, shortage occurs since tr0>td. Hence(35)ES∣τ=τ0,y=y0,tr=tr0=Cstr0-tdd.Other costs are the same as Scenario 5. Therefore(36)Ea6τ0,y0,tr0=Ea5τ0,y0,tr0+ES∣τ=τ0,y=y0,tr=tr0=Ea5τ0,y0,tr0+Cstr0-tddwhere Ea6(τ0,y0,tr0) denotes E(Ca∣τ=τ0,y=y0,tr=tr0) when tb<τ0<w1, 1≤y0≤k, tr0>td. From Figure 6, it can be obtained that(37)ETc∣τ=τ0,y=y0,tr=tr0=wy0+tr0=Ec6y0,tr0where Ec6(y0,tr0) denotes ETc∣τ=τ0,y=y0,tr=tr0 when tb<τ0<w1, 1≤y0≤k, tr0>td.Figure 6
Behavior of inventory (Scenario 6).


## 4.1.7. Scenario 7:tb<τ<w1, y≥k+1, 0<tr<td

This scenario can be rewritten asτ=τ0,y=y0,tr=tr0 and tb<τ0<w1, y0≥k+1, 0<tr0<td. Thus the behavior of inventory can be described as Figure 7. The alert signal is not sent in the first k samples. Hence the preventive maintenance (PM) is implemented at the end of the (k+1)th sampling interval (without sampling). Then the assignable cause is detected during the PM. Thus the PM is replaced by PRDM. Because 0<tr0<td, shortage does not occur.Figure 7
Behavior of inventory (Scenario 7).The conditional expectations of all kinds of costs given thatτ=τ0,y=y0,tr=tr0 can be calculated using the same analogy as Scenario 1. And they are displayed in Tables 1 and 2 of the Annex. Therefore(38)ECa∣τ=τ0,y=y0,tr=tr0=Cc∫0wk+1λXtdt+K∫0wk+1λXtdt+1+Crtr0+ks+h12tbH+wk+1-tbH+12Htd+wk+1-τ0dθ2+Umaxtb+τ0-tbdθ1Ciw+wk+1-τ0d1-θ2+Umaxtb+τ0-tbd1-θ1Cpw=Ea7τ0,tr0where Ea7(τ0,tr0) denotes E(Ca∣τ=τ0,y=y0,tr=tr0) when tb<τ0<w1, y0≥k+1, 0<tr0<td.From Figure7, it can be obtained that(39)ETc∣τ=τ0,y=y0,tr=tr0=wk+1+td=Ec7where Ec7 denotes E(Tc∣τ=τ0,y=y0,tr=tr0) when tb<τ0<w1, y0≥k+1, 0<tr0<td.


## 4.1.8. Scenario 8:tb<τ<w1, y≥k+1, tr>td

This scenario can be rewritten asτ=τ0,y=y0,tr=tr0 and tb<τ0<w1, y0≥k+1, tr0>td. Thus the behavior of inventory can be described as Figure 8. In contrast to Scenario 7, shortage occurs since tr0>td. Hence(40)ES∣τ=τ0,y=y0,tr=tr0=Cstr0-tdd.Other costs are the same as Scenario 7. Therefore(41)Ea8τ0,tr0=Ea7τ0,tr0+ES∣τ=τ0,y=y0,tr=tr0=Ea7τ0,tr0+Cstr0-tddwhere Ea8(τ0,tr0) denotes E(Ca∣τ=τ0,y=y0,tr=tr0) when tb<τ0<w1, y0≥k+1, tr0>td. From Figure 8, it can be obtained that(42)ETc∣τ=τ0,y=y0,tr=tr0=wk+1+tr0=Ec8tr0where Ec8 denotes E(Tc∣τ=τ0,y=y0,tr=tr0) when tb<τ0<w1, y0≥k+1, tr0>td.Figure 8
Behavior of inventory (Scenario 8).


## 4.1.9. Scenario 9:wj<τ<wj+1(1≤j≤k-1), 1≤y≤k-j, tr<td

This scenario can be rewritten asτ=τ0,y=y0,tr=tr0 and wj<τ0<wj+1(1≤j≤k-1),1≤y0≤k-j,tr0<td. Thus the behavior of inventory can be described as Figure 9. In this scenario, the process shifts to out-of-control state in the (j+1)th sampling interval. After the shift, due to the type II error, the alert signal is not sent until the (j+y0)th sampling, while the alert signal may be sent in the first j+1 samplings due to the type I error. After the (j+y0)th sampling, the PRDM is conducted for tr0. Because 0<tr0<td, shortage does not occur.Figure 9
Behavior of inventory (Scenario 9).The conditional expectations of all kinds of costs given thatτ=τ0,y=y0,tr=tr0 can be calculated using the same analogy as Scenario 1. And they are displayed in Tables 1 and 2 of the Annex. Therefore(43)ECa∣τ=τ0,y=y0,tr=tr0=Cc∫0wj+y0λXtdt+K∫0wj+y0λXtdt+1+Crtr0+j+y0s+Cyjα+1+h12tbH+wj+y0-tbH+12Htd+wj+y0-τ0dθ2+Umaxtb+τ0-tbdθ1Ciw+wj+y0-τ0d1-θ2+Umaxtb+τ0-tbd1-θ1Cpw=Ea9τ0,y0,tr0where Ea9(τ0,y0,tr0) denotes E(Ca∣τ=τ0,y=y0,tr=tr0) when wj<τ0<wj+1(1≤j≤k-1),1≤y0≤k-j,tr0<td.From Figure9, it can be obtained that(44)ETc∣τ=τ0,y=y0,tr=tr0=wj+y0+td=Ec9y0where Ec9(y0) denotes E(Tc∣τ=τ0,y=y0,tr=tr0) when wj<τ0<wj+1(1≤j≤k-1),1≤y0≤k-j,tr0<td.


## 4.1.10. Scenario 10:wj<τ<wj+1(1≤j≤k-1), 1≤y≤k-j, tr>td

This scenario can be rewritten asτ=τ0,y=y0,tr=tr0 and wj<τ0<wj+1(1≤j≤k-1),1≤y0≤k-j,tr0>td. Thus the behavior of inventory can be described as Figure 10. In contrast to Scenario 9, shortage occurs since tr0>td. Hence(45)ES∣τ=τ0,y=y0,tr=tr0=Cstr0-tdd.Other costs are the same as Scenario 9. Therefore(46)Ea10τ0,y0,tr0=Ea9τ0,y0,tr0+ES∣τ=τ0,y=y0,tr=tr0=Ea9τ0,y0,tr0+Cstr0-tddwhere Ea10(τ0,y0,tr0) denotes E(Ca∣τ=τ0,y=y0,tr=tr0) when wj<τ0<wj+1(1≤j≤k-1),1≤y0≤k-j,tr0>td. From Figure 10, it can be obtained that(47)ETc∣τ=τ0,y=y0,tr=tr0=wj+y0+tr0=Ec10tr0where Ec10(tr0) denotes E(Tc∣τ=τ0,y=y0,tr=tr0) when wj<τ0<wj+1(1≤j≤k-1),1≤y0≤k-j,tr0>td.Figure 10
Behavior of inventory (Scenario 10).


## 4.1.11. Scenario 11:wj<τ<wj+1(1≤j≤k-1), y>k-j, tr<td

This scenario can be rewritten asτ=τ0,y=y0,tr=tr0 and wj<τ0<wj+1(1≤j≤k-1),y0>k-j,tr0<td. Thus the behavior of inventory can be described as Figure 11. No PRDM is conducted after k samplings. Hence the preventive maintenance (PM) is implemented at the end of the (k+1)th sampling interval (without sampling). Then the assignable cause is detected during the PM. Thus the PM is replaced by the PRDM. Because 0<tr0<td, shortage does not occur.Figure 11
Behavior of inventory (Scenario 11).The conditional expectations of all kinds of costs given thatτ=τ0,y=y0,tr=tr0 can be calculated using the same analogy as Scenario 1. And they are displayed in Tables 1 and 2 of the Annex. Therefore (48)ECa∣τ=τ0,y=y0,tr=tr0=Cc∫0wk+1λXtdt+K∫0wk+1λXtdt+1+Crtr0+ks+Cyjα+h12tbH+wk+1-tbH+12Htd+wk+1-τ0dθ2+Umaxtb+τ0-tbdθ1Ciw+wk+1-τ0d1-θ2+Umaxtb+τ0-tbd1-θ1Cpw=Ea11τ0,tr0where Ea11(τ0,tr0) denotes E(Ca∣τ=τ0,y=y0,tr=tr0) when wj<τ0<wj+1(1≤j≤k-1),y0>k-j,tr0<td.From Figure11, it can be obtained that(49)ETc∣τ=τ0,y=y0,tr=tr0=wk+1+td=Ec11where Ec11 denotes E(Tc∣τ=τ0,y=y0,tr=tr0) when wj<τ0<wj+1(1≤j≤k-1),y0>k-j,tr0<td.


## 4.1.12. Scenario 12:wj<τ<wj+1(1≤j≤k-1), y>k-j, tr>td

This scenario can be rewritten asτ=τ0,y=y0,tr=tr0 and wj<τ0<wj+1(1≤j≤k-1),y0>k-j,tr0>td. Thus the behavior of inventory can be described as Figure 12. In contrast to Scenario 11, shortage occurs since tr0>td. Hence(50)ES∣τ=τ0,y=y0,tr=tr0=Cstr0-tdd.Other costs are the same as Scenario 11. Therefore(51)Ea12τ0,tr0=Ea11τ0,tr0+ES∣τ=τ0,y=y0,tr=tr0=Ea11τ0,tr0+Cstr0-tddwhere Ea12(τ0,tr0) denotes E(Ca∣τ=τ0,y=y0,tr=tr0) when wj<τ0<wj+1(1≤j≤k-1),y0>k-j,tr0>td. From Figure 12, it can be obtained that(52)ETc∣τ=τ0,y=y0,tr=tr0=wk+1+tr0=Ec12tr0where Ec12(tr0) denotes E(Tc∣τ=τ0,y=y0,tr=tr0) when wj<τ0<wj+1(1≤j≤k-1),y0>k-j,tr0>td.Figure 12
Behavior of inventory (Scenario 12).


## 4.1.13. Scenario 13:wk<τ<wk+1,tr<td

This scenario can be rewritten asτ=τ0,tr=tr0 and wk<τ0<wk+1,tr0<td. Thus the behavior of inventory can be described as Figure 13. In this scenario, the process shifts to out-of-control state in the (k+1)th sampling interval. No PRDM is conducted after k samplings. Hence the preventive maintenance (PM) is implemented at the end of the (k+1)th sampling interval (without sampling). Then the assignable cause is detected during the PM. Thus the PM is replaced by PRDM. Because 0<tr0<td, shortage does not occur.Figure 13
Behavior of inventory (Scenario 13).The conditional expectations of all kinds of costs given thatτ=τ0,tr=tr0 can be calculated using the same analogy as Scenario 1. And they are displayed in Tables 1 and 2 of the Annex. Therefore (53)ECa∣τ=τ0,tr=tr0=Cc∫0wk+1λXtdt+K∫0wk+1λXtdt+1+Crtr0+ks+Cykα+h12tbH+wk+1-tbH+12Htd+wk+1-τ0dθ2+Umaxtb+τ0-tbdθ1Ciw+wk+1-τ0d1-θ2+Umaxtb+τ0-tbd1-θ1Cpw=Ea13τ0,tr0where Ea13(τ0,tr0) denotes E(Ca∣τ=τ0,tr=tr0) when wk<τ0<wk+1,tr0<td.From Figure13, it can be obtained that(54)ETc∣τ=τ0,tr=tr0=wk+1+td=Ec13where Ec13 denotes E(Tc∣τ=τ0,tr=tr0) when wk<τ0<wk+1,tr0<td.


## 4.1.14. Scenario 14:wk<τ<wk+1,tr>td

This scenario can be rewritten asτ=τ0,tr=tr0 and wk<τ0<wk+1,tr0>td. Thus the behavior of inventory can be described as Figure 14. In contrast to Scenario 13, shortage occurs since tr0>td. Hence(55)ES∣τ=τ0,tr=tr0=Cstr0-tdd.Other costs are the same as Scenario 13. Therefore(56)Ea14τ0,tr0=Ea13τ0,tr0+ES∣τ=τ0,tr=tr0=Ea13τ0,tr0+Cstr0-tddwhere Ea14(τ0,tr0) denotes E(Ca∣τ=τ0,tr=tr0) when wk<τ0<wk+1,tr0>td. From Figure 14, it can be obtained that(57)ETc∣τ=τ0,tr=tr0=wk+1+tr0=Ec14tr0where Ec14(tr0) denotes E(Tc∣τ=τ0,tr=tr0) when wk<τ0<wk+1,tr0>td.Figure 14
Behavior of inventory (Scenario 14).


## 4.1.15. Scenario 15:τ>wk+1,tp<td

This scenario can be rewritten asτ=τ0,tp=tp0 and τ0>wk+1,tp0<td. Thus the behavior of inventory can be described as Figure 15. In this scenario, the process keeps in the in-control state during the (k+1) sampling intervals. Hence the preventive maintenance (PM) is implemented at the end of the (k+1)th sampling interval (without sampling). Because tp0<td, shortage does not occur.Figure 15
Behavior of inventory (Scenario 15).The conditional expectations of all kinds of costs given thatτ=τ0,tp=tp0 can be calculated using the same analogy as Scenario 1. And they are displayed in Tables 1 and 2 of the Annex. Therefore(58)ECa∣τ=τ0,tp=tp0=Cc∫0wk+1λXtdt+K∫0wk+1λXtdt+1+Cptp0+ks+Cykα+h12tbH+wk+1-tbH+12Htd+Umaxtb+wk+1-tbdθ1Ciw+Umaxtb+wk+1-tbd1-θ1Cpw=Ea15tp0where Ea15(tp0) denotes E(Ca∣τ=τ0,tp=tp0) when τ0>wk+1,tp0<td.From Figure15, it can be obtained that(59)ETc∣τ=τ0,tp=tp0=wk+1+td=Ec15where Ec15 denotes E(Tc∣τ=τ0,tp=tp0) when τ0>wk+1,tp0<td.


## 4.1.16. Scenario 16:τ>wk+1,tp>td

This scenario can be rewritten asτ=τ0,tp=tp0 and τ0>wk+1,tp0>td. Thus the behavior of inventory can be described as Figure 16. In contrast to Scenario 15, shortage occurs since tp0>td. Hence(60)ES∣τ=τ0,tp=tp0=Cstp0-tdd.Other costs are the same as Scenario 15. Therefore(61)Ea16tp0=Ea15tp0+ES∣τ=τ0,y=y0,tr=tr0=Ea15tp0+Cstp0-tddwhere Ea16(tp0) denotes E(Ca∣τ=τ0,tp=tp0) when τ0>wk+1,tp0>td. From Figure 16, it can be obtained that(62)ETc∣τ=τ0,tp=tp0=wk+1+tp0=Ec16tp0where Ec16(tp0) denotes E(Tc∣τ=τ0,tp=tp0) when τ0>wk+1,tp0>td.Figure 16
Behavior of inventory (Scenario 16).


## 4.2. Objective Function and Constraints

By applying the renewal theory, the total expected cost per unit timeETC can be obtained as follows:(63)ETC=ETCH,k,h1,L,n=ECaETc.Thus the nonlinear optimization problem can be obtained:(64)minETCs.t:(65)h1>tb,(66)H≤Hm,H∈N+,(67)k≤km,k∈N+,(68)0<L<Lm,(69)1≤n≤nm.The constraints are explained as follows. Inequality (65): in order to avoid overmuch sampling cost, the first sampling should be conducted after the hedging safety stock is built up. Inequality (66): in practice, the available inventory is finite. Inequality (67): in order to avoid excessive samplings, the value of k must be limited. Inequality (68): in order to decrease the error of the control chart, the value of L must be limited. Inequality (69): in practice, because of economic reasons, the sample size must be less than a definite number.


## 5. Solution Approach

The proposed model is a nonlinear program with a complicated objective function. Metaheuristic search algorithms are suitable for this kind of optimization problem. In the process of metaheuristic search algorithms, the objective function values of different decision variables must be calculated for many times. Thus, the method of calculating the objective function value is very critical. If the calculating speed is too low, the metaheuristic search algorithms will be unfeasible.In order to calculate the objective function,E(Ca) and E(Tc) must be calculated first. In the majority of papers, E(Ca) and E(Tc) are calculated by the integral method. However, because the proposed model has 16 scenarios, the integral calculation will be too complex and take too much time (see Section 6.2). In contrast, the Monte Carlo method can avoid the complex integral calculation. However, the Monte Carlo method needs enough repeat operations to guarantee its accuracy. In general, the repeat operations are realized by a loop in calculation software. The loop with enough repetition will also take too much time. In order to solve this problem, this paper proposes an approach of realizing the Monte Carlo simulation without a loop. As a result, the speed of the Monte Carlo simulation is increased enough to be used to calculate the objective function in the process of metaheuristic search algorithms. In this paper, the objective function is calculated by the Monte Carlo method and optimized by the genetic algorithm, which is one of the most important metaheuristic search algorithms.


### 5.1. Monte Carlo Method

To be specific, the Monte Carlo simulation withN trials is coded ingeniously in Matlab. Because Matlab is good at matrix manipulation, the N trials correspond to N-by-1 matrixes rather than N iterations of a loop. The operation process in Matlab is described briefly as follows.Step 1 (assign the model parameters).
Then the conditional expectations ofCa and Tc in 16 scenarios, Eai and Eci(1≤i≤16), can be calculated according to Section 4.1. Afterwards, Eai and Eci are all vectorized.Step 2.
The conditions of 16 scenarios are formulated according to Section4.1. For example, the condition of Scenario 1 can be formulated as follows.(70)C1=τ0<tb&y0≤k&tr0<td.The conditions of Scenarios 2-16 Ci(2≤i≤16) can also be formulated like C1. Note that Ci(1≤i≤16) are all logical value (0 or 1).Step 3.
N random numbers are generated for the distribution of τ. These random numbers are stored in a N-by-1 matrix called τ-. Similarly, the matrixes called X-, tp¯, tr¯, and y- can also be obtained. The random numbers generated for the distributions of X, tp, tr, and y are stored in X-, tp¯, tr¯, and y-, respectively.Step 4.
Substituteτ-, X-, tp¯, tr¯, and y- into the formulation of Eai, Eci, and Ci(1≤i≤16). Thus Eai, Eci, and Ci(1≤i≤16) become N-by-1 matrixes Eai¯, Eci¯, and Ci¯, respectively. Then the matrixes Ca¯ and Tc¯ can be obtained as follows.(71)Ca¯=∑i=116Eai¯.∗Ci¯.(72)Tc¯=∑i=116Eci¯.∗Ci¯.Here, each element of Ca¯ is the value of Ca in one trial and each element of Tc¯ is the value of Tc in one trial.Step 5.
The mean of elements inCa¯ is regarded as the approximate value of E(Ca). The mean of elements in Tc¯ is regarded as the approximate value of E(Tc). As long as the value of N is high enough, the precision of the two approximate values is enough.


### 5.2. Genetic Algorithms

The genetic algorithms (GA) originated from the model of [39] and have been a widely used metaheuristic search algorithm to solve the optimization problem in combinatorial optimization, machine learning, signal processing, etc. Based on the phenomenon in species evolution, including inheritance, mutation, natural selection, and hybridization, the GA uses the concept of chromosome to describe the solution of a problem. The initial population of chromosomes is generated by a coding process. Through the crossover and the mutation process, the chromosome with better fitness function (objective function) value wins the competition generation by generation. Finally, the optimal solution of the problem is found. In this paper, the GA is conducted by the Optimization App in matlab (R2016a).


## 5.1. Monte Carlo Method

To be specific, the Monte Carlo simulation withN trials is coded ingeniously in Matlab. Because Matlab is good at matrix manipulation, the N trials correspond to N-by-1 matrixes rather than N iterations of a loop. The operation process in Matlab is described briefly as follows.Step 1 (assign the model parameters).
Then the conditional expectations ofCa and Tc in 16 scenarios, Eai and Eci(1≤i≤16), can be calculated according to Section 4.1. Afterwards, Eai and Eci are all vectorized.Step 2.
The conditions of 16 scenarios are formulated according to Section4.1. For example, the condition of Scenario 1 can be formulated as follows.(70)C1=τ0<tb&y0≤k&tr0<td.The conditions of Scenarios 2-16 Ci(2≤i≤16) can also be formulated like C1. Note that Ci(1≤i≤16) are all logical value (0 or 1).Step 3.
N random numbers are generated for the distribution of τ. These random numbers are stored in a N-by-1 matrix called τ-. Similarly, the matrixes called X-, tp¯, tr¯, and y- can also be obtained. The random numbers generated for the distributions of X, tp, tr, and y are stored in X-, tp¯, tr¯, and y-, respectively.Step 4.
Substituteτ-, X-, tp¯, tr¯, and y- into the formulation of Eai, Eci, and Ci(1≤i≤16). Thus Eai, Eci, and Ci(1≤i≤16) become N-by-1 matrixes Eai¯, Eci¯, and Ci¯, respectively. Then the matrixes Ca¯ and Tc¯ can be obtained as follows.(71)Ca¯=∑i=116Eai¯.∗Ci¯.(72)Tc¯=∑i=116Eci¯.∗Ci¯.Here, each element of Ca¯ is the value of Ca in one trial and each element of Tc¯ is the value of Tc in one trial.Step 5.
The mean of elements inCa¯ is regarded as the approximate value of E(Ca). The mean of elements in Tc¯ is regarded as the approximate value of E(Tc). As long as the value of N is high enough, the precision of the two approximate values is enough.


## 5.2. Genetic Algorithms

The genetic algorithms (GA) originated from the model of [39] and have been a widely used metaheuristic search algorithm to solve the optimization problem in combinatorial optimization, machine learning, signal processing, etc. Based on the phenomenon in species evolution, including inheritance, mutation, natural selection, and hybridization, the GA uses the concept of chromosome to describe the solution of a problem. The initial population of chromosomes is generated by a coding process. Through the crossover and the mutation process, the chromosome with better fitness function (objective function) value wins the competition generation by generation. Finally, the optimal solution of the problem is found. In this paper, the GA is conducted by the Optimization App in matlab (R2016a).


## 6. Experimental Results

### 6.1. Numerical Example

τ, X and the times to failure of perfect and imperfect products are assumed to follow the Weibull distribution. tp and tr are assumed to follow the exponential distribution. Therefore, their probability density functions and hazard rate functions can be assumed as follows.(73)fτt=λvtv-1e-λtv.(74)λτt=λvtv-1.(75)λXt=λ0v0tv0-1.(76)λ1t=λ1v1tv1-1.(77)λ2t=λ2v2tv2-1.(78)fpt=1μpe-t/μp.(79)frt=1μre-t/μr.The values of the model parameters are recorded in Table 2, where the unit for time is a day and the unit for money is ¥.Table 2
Parameter values.
Parameter U m a x θ 1 θ 2 w μ p μ r C w δ C F C V C c C p C r C y h Value 100 0.15 0.35 4 10 20 100 2 1 0.1 60 700 800 5 0.1  Parameter C s K d λ 1 v 1 λ 2 v 2 λ 0 v 0 λ v H m k m L m n m Value 15 2 50 (1/36)2 2 (1/12)2 2 0.025 2 0.05 2 1000 10 4 100This optimization problem is solved with the Monte Carlo method and genetic algorithms. The optimization result isH∗=229, h1∗=102.9237, k∗=3, L∗=0.3590, n∗=5, ETC∗=576.5120.


### 6.2. Calculation Comparison

In this section, the integral method of calculating the objective function is described. Then the integral method and Monte Carlo method are compared through an orthogonal experiment.


#### 6.2.1. Integral Method

The expectation ofCa can be computed by conditioning τ:(80)ECa=∫0∞ECa∣τ=τ0fττ0dτ0=∫0tbECa∣τ=τ0fττ0dτ0+∫tbw1ECa∣τ=τ0fττ0dτ0+∑j=1k-1∫wjwj+1ECa∣τ=τ0fττ0dτ0+∫wkwk+1ECa∣τ=τ0fττ0dτ0+∫wk+1∞ECa∣τ=τ0fττ0dτ0.Let E0 denote the value of E(Ca∣τ=τ0) when 0<τ0<tb. Eb denotes E(Ca∣τ=τ0) when tb<τ0<w1. Ej denotes E(Ca∣τ=τ0) when wj<τ0<wj+1(1≤j≤k-1). Ek denotes E(Ca∣τ=τ0) when wk<τ0<wk+1. Ek+1 denote E(Ca∣τ=τ0) when τ0>wk+1. Thus(81)ECa=∫0tbE0fττ0dτ0+∫tbw1Ebfττ0dτ0+∑j=1k-1∫wjwj+1Ejfττ0dτ0+∫wkwk+1Ekfττ0dτ0+∫wk+1∞Ek+1fττ0dτ0.E0 can be computed by conditioning y:(82)E0=∑y0=1∞ECa∣y=y0,τ=τ0Py=y0∣τ=τ0=∑y0=1kECa∣y=y0,τ=τ0Py=y0∣τ=τ0+∑y0=k+1∞ECa∣y=y0,τ=τ0Py=y0∣τ=τ0.E01 denotes E(Ca∣y=y0,τ=τ0) when 0<τ0<tb,1≤y0≤k. E02 denotes E(Ca∣y=y0,τ=τ0) when 0<τ0<tb,y0≥k+1. Thus(83)E0=∑y0=1kE01Py=y0∣τ=τ0+∑y0=k+1∞E02Py=y0∣τ=τ0.E01 and E02 can be computed by conditioning tr:(84)E01=∫0∞ECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0=∫0H/dECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0+∫H/d∞ECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0=∫0H/dEa1τ0,y0,tr0frtr0dtr0+∫H/d∞Ea2τ0,y0,tr0frtr0dtr0,(85)E02=∫0∞ECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0=∫0H/dECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0+∫H/d∞ECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0=∫0H/dEa3τ0,tr0frtr0dtr0+∫H/d∞Ea4τ0,tr0frtr0dtr0.Eb can be computed by conditioning y:(86)Eb=∑y0=1∞ECa∣y=y0,τ=τ0Py=y0∣τ=τ0=∑y0=1kECa∣y=y0,τ=τ0Py=y0∣τ=τ0+∑y0=k+1∞ECa∣y=y0,τ=τ0Py=y0∣τ=τ0.Eb1 denotes E(Ca∣y=y0,τ=τ0) when tb<τ0<w1,1≤y0≤k. Eb2 denotes E(Ca∣y=y0,τ=τ0) when tb<τ0<w1,y0≥k+1. Thus(87)Eb=∑y0=1kEb1Py=y0∣τ=τ0+∑y0=k+1∞Eb2Py=y0∣τ=τ0.Eb1 and Eb2 can be computed by conditioning tr:(88)Eb1=∫0∞ECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0=∫0H/dECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0+∫H/d∞ECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0=∫0H/dEa5τ0,y0,tr0frtr0dtr0+∫H/d∞Ea6τ0,y0,tr0frtr0dtr0,(89)Eb2=∫0∞ECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0=∫0H/dECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0+∫H/d∞ECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0=∫0H/dEa7τ0,tr0frtr0dtr0+∫H/d∞Ea8τ0,tr0frtr0dtr0.Ej can be computed by conditioning y:(90)Ej=∑y0=1∞ECa∣y=y0,τ=τ0Py=y0∣τ=τ0=∑y0=1k-jECa∣y=y0,τ=τ0Py=y0∣τ=τ0+∑y0=k-j+1∞ECa∣y=y0,τ=τ0Py=y0∣τ=τ0.Ej1 denotes E(Ca∣y=y0,τ=τ0) when wj<τ0<wj+1,1≤y0≤k-j. Ej2 denotes E(Ca∣y=y0,τ=τ0) when wj<τ0<wj+1,y0≥k-j+1. Thus(91)Ej=∑y0=1k-jEj1Py=y0∣τ=τ0+∑y0=k-j+1∞Ej2Py=y0∣τ=τ0.Ej1 and Ej2 can be computed by conditioning tr:(92)Ej1=∫0∞ECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0=∫0H/dECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0+∫H/d∞ECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0=∫0H/dEa9τ0,y0,tr0frtr0dtr0+∫H/d∞Ea10τ0,y0,tr0frtr0dtr0,(93)Ej2=∫0∞ECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0=∫0H/dECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0+∫H/d∞ECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0=∫0H/dEa11τ0,tr0frtr0dtr0+∫H/d∞Ea12τ0,tr0frtr0dtr0.Ek can be computed by conditioning tr:(94)Ek=∫0∞ECa∣tr=tr0,τ=τ0frtr0dtr0=∫0tdECa∣tr=tr0,τ=τ0frtr0dtr0+∫td∞ECa∣tr=tr0,τ=τ0frtr0dtr0=∫0tdEa13τ0,tr0frtr0dtr0+∫td∞Ea14τ0,tr0frtr0dtr0.Ek+1 can be computed by conditioning tp:(95)Ek+1=∫0∞ECa∣tp=tp0,τ=τ0fptp0dtp0=∫0H/dECa∣tp=tp0,τ=τ0fptp0dtp0+∫H/d∞ECa∣tp=tp0,τ=τ0fptp0dtp0=∫0H/dEa15τ0,tp0fptp0dtp0+∫H/d∞Ea16τ0,tp0fptp0dtp0.Combing Equations (81)-(95), it can be obtained that(96)ECa=∫0tbE0fττ0dτ0+∫tbw1Ebfττ0dτ0+∑j=1k-1∫wjwj+1Ejfττ0dτ0+∫wkwk+1Ekfττ0dτ0+∫wk+1∞Ek+1fττ0dτ0where(97)E0=∑y0=1kE01Py=y0∣τ=τ0+∑y0=k+1∞E02Py=y0∣τ=τ0,(98)E01=∫0tdEa1τ0,y0,tr0frtr0dtr0+∫td∞Ea2τ0,y0,tr0frtr0dtr0,(99)E02=∫0tdEa3τ0,tr0frtr0dtr0+∫td∞Ea4τ0,tr0frtr0dtr0,(100)Eb=∑y0=1kEb1Py=y0∣τ=τ0+∑y0=k+1∞Eb2Py=y0∣τ=τ0,(101)Eb1=∫0tdEa5τ0,y0,tr0frtr0dtr0+∫td∞Ea6τ0,y0,tr0frtr0dtr0,(102)Eb2=∫0tdEa7τ0,tr0frtr0dtr0+∫td∞Ea8τ0,tr0frtr0dtr0,(103)Ej=∑y0=1k-jEj1Py=y0∣τ=τ0+∑y0=k-j+1∞Ej2Py=y0∣τ=τ0,(104)Ej1=∫0tdEa9τ0,y0,tr0frtr0dtr0+∫td∞Ea10τ0,y0,tr0frtr0dtr0,(105)Ej2=∫0tdEa11τ0,tr0frtr0dtr0+∫td∞Ea12τ0,tr0frtr0dtr0,(106)Ek=∫0tdEa13τ0,tr0frtr0dtr0+∫td∞Ea14τ0,tr0frtr0dtr0,(107)Ek+1=∫0tdEa15τ0,tp0fptp0dtp0+∫0tdEa16τ0,tp0fptp0dtp0.Similarly, E(Tc) can also be obtained:(108)ETc=∫0tbE0fττ0dτ0+∫tbw1Ebfττ0dτ0+∑j=1k-1∫wjwj+1Ejfττ0dτ0+∫wkwk+1Ekfττ0dτ0+∫wk+1∞Ek+1fττ0dτ0where(109)E0=∑y0=1kE01Py=y0∣τ=τ0+∑y0=k+1∞E02Py=y0∣τ=τ0,(110)E01=∫0tdEc1y0frtr0dtr0+∫td∞Ec2y0,tr0frtr0dtr0,(111)E02=∫0tdEc3frtr0dtr0+∫td∞Ec4tr0frtr0dtr0,(112)Eb=∑y0=1kEb1Py=y0∣τ=τ0+∑y0=k+1∞Eb2Py=y0∣τ=τ0,(113)Eb1=∫0tdEc5y0frtr0dtr0+∫td∞Ec6y0,tr0frtr0dtr0,(114)Eb2=∫0tdEc7frtr0dtr0+∫td∞Ec8tr0frtr0dtr0,(115)Ej=∑y0=1k-jEj1Py=y0∣τ=τ0+∑y0=k-j+1∞Ej2Py=y0∣τ=τ0,(116)Ej1=∫0tdEc9y0frtr0dtr0+∫td∞Ec10tr0frtr0dtr0,(117)Ej2=∫0tdEc11frtr0dtr0+∫td∞Ec12tr0frtr0dtr0,(118)Ek=∫0tdEc13τ0,tr0frtr0dtr0+∫td∞Ec14τ0,tr0frtr0dtr0,(119)Ek+1=∫0tdEc15fptp0dtp0+∫td∞Ec16tp0fptp0dtp0.Thus the values ofE(Ca) and E(Tc) are obtained and the objective function value can be calculated by Equation (63).


#### 6.2.2. Comparison of Calculation Methods

In this section, an orthogonal experiment is employed to compare the integral method and the Monte Carlo method. Table3 presents the high and low levels of the 5 input parameters in the orthogonal experiment, which are the decision variables of the objective function. Other parameter values are the same as those in Section 6.1. The L8 orthogonal array is used to carry out the experiment as shown in Table 4. The objective function (ETC) value is calculated by the two methods in every trial. The experiment result is also shown in Table 4. The relative error forETC value of the Monte Carlo method is between 0.0487% and 0.0713%. In addition, the computation speed of the Monte Carlo method is much higher than that of the integral method. So this paper uses the Monte Carlo method to calculate the objective function value in the process of the genetic algorithm.Table 3
Input parameters for orthogonal experiment.
Parameter H h 1 k L n Low 100 50 5 0.1 3 High 200 100 10 0.2 6Table 4
Comparison of the integral and Monte Carlo methods.
No. H h 1 k L n ETC value Time Integral method Monte Carlo method R e l a t i v e   e r r o r Integral method Monte Carlo method 1 100 50 5 0.1 3 650.1860 649.7387 0.0688% 7.538693s 0.061884s 2 100 50 5 0.2 6 650.2032 649.7520 0.0694% 2.378531s 0.094823s 3 100 100 10 0.1 3 578.1220 577.8375 0.0492% 7.272525s 0.117649s 4 100 100 10 0.2 6 578.1217 577.8374 0.0492% 4.626247s 0.060064s 5 200 50 10 0.1 6 641.2674 640.8161 0.0704% 6.499943s 0.070176s 6 200 50 10 0.2 3 641.2382 640.7808 0.0713% 4.793808s 0.058871s 7 200 100 5 0.1 6 577.2176 576.9364 0.0487% 4.002955s 0.054454s 8 200 100 5 0.2 3 577.2226 576.9414 0.0487% 2.219563s 0.052333s


### 6.3. Sensitivity Analysis

In this section, a Design of Experiments (DOE) study is conducted to investigate the sensitivity of the model parameters with the optimal expected cost per unit time (ETC∗).Table5 presents the high and low levels of the 8 input parameters in the DOE, involving production, inventory, maintenance, and quality. Other model parameters are the same as those in Section 6.1. Employing a fractional factorial design, the L16 orthogonal array is used to carry out the orthogonal experiment as shown in Table 6. The GA is used in every trial to obtain the optimal solution of the model. ETC∗ of every trial is also shown in Table 6. In order to analyze the effect of model parameters on ETC∗, the normal probability plot and main effect plot of the model parameters are depicted in Figures 17 and 18 by using Minitab.18. Figures 17 and 18 show that Cw and Cr have significant positive effect on ETC∗.Table 5
Parameter value for DOE.
Parameter U m a x C w C F C V C p C r C y h Low 100 100 1 0.1 700 800 5 0.1 High 200 200 2 0.2 1400 1600 10 0.2Table 6
Parameter values andETC∗ of 16 trials.
No. U m a x C w C F C V C p C r C y h E T C ∗ 1 100 100 1 0.1 700 800 5 0.1 576.5120 2 100 100 1 0.2 700 1600 10 0.2 685.7579 3 100 100 2 0.1 1400 800 10 0.2 578.1663 4 100 100 2 0.2 1400 1600 5 0.1 686.5931 5 100 200 1 0.1 1400 1600 5 0.2 889.3876 6 100 200 1 0.2 1400 800 10 0.1 771.2377 7 100 200 2 0.1 700 1600 10 0.1 889.7835 8 100 200 2 0.2 700 800 5 0.2 771.8040 9 200 100 1 0.1 1400 1600 10 0.1 685.8556 10 200 100 1 0.2 1400 800 5 0.2 577.7022 11 200 100 2 0.1 700 1600 5 0.2 687.0148 12 200 100 2 0.2 700 800 10 0.1 575.9532 13 200 200 1 0.1 700 800 10 0.2 771.9917 14 200 200 1 0.2 700 1600 5 0.1 888.6614 15 200 200 2 0.1 1400 800 5 0.1 771.6710 16 200 200 2 0.2 1400 1600 10 0.2 889.2610Figure 17
Normal probability plot of the effects.Figure 18
Main effects plot.According to the above analysis result, some suggestions for managers are proposed. In order to decrease the expected cost per unit time, the managers could focus on the process of predictive maintenance to decrease its cost per unit time. In addition, the managers could outsource the warranty to decrease the repair cost for warranty.


## 6.1. Numerical Example

τ, X and the times to failure of perfect and imperfect products are assumed to follow the Weibull distribution. tp and tr are assumed to follow the exponential distribution. Therefore, their probability density functions and hazard rate functions can be assumed as follows.(73)fτt=λvtv-1e-λtv.(74)λτt=λvtv-1.(75)λXt=λ0v0tv0-1.(76)λ1t=λ1v1tv1-1.(77)λ2t=λ2v2tv2-1.(78)fpt=1μpe-t/μp.(79)frt=1μre-t/μr.The values of the model parameters are recorded in Table 2, where the unit for time is a day and the unit for money is ¥.Table 2
Parameter values.
Parameter U m a x θ 1 θ 2 w μ p μ r C w δ C F C V C c C p C r C y h Value 100 0.15 0.35 4 10 20 100 2 1 0.1 60 700 800 5 0.1  Parameter C s K d λ 1 v 1 λ 2 v 2 λ 0 v 0 λ v H m k m L m n m Value 15 2 50 (1/36)2 2 (1/12)2 2 0.025 2 0.05 2 1000 10 4 100This optimization problem is solved with the Monte Carlo method and genetic algorithms. The optimization result isH∗=229, h1∗=102.9237, k∗=3, L∗=0.3590, n∗=5, ETC∗=576.5120.


## 6.2. Calculation Comparison

In this section, the integral method of calculating the objective function is described. Then the integral method and Monte Carlo method are compared through an orthogonal experiment.


### 6.2.1. Integral Method

The expectation ofCa can be computed by conditioning τ:(80)ECa=∫0∞ECa∣τ=τ0fττ0dτ0=∫0tbECa∣τ=τ0fττ0dτ0+∫tbw1ECa∣τ=τ0fττ0dτ0+∑j=1k-1∫wjwj+1ECa∣τ=τ0fττ0dτ0+∫wkwk+1ECa∣τ=τ0fττ0dτ0+∫wk+1∞ECa∣τ=τ0fττ0dτ0.Let E0 denote the value of E(Ca∣τ=τ0) when 0<τ0<tb. Eb denotes E(Ca∣τ=τ0) when tb<τ0<w1. Ej denotes E(Ca∣τ=τ0) when wj<τ0<wj+1(1≤j≤k-1). Ek denotes E(Ca∣τ=τ0) when wk<τ0<wk+1. Ek+1 denote E(Ca∣τ=τ0) when τ0>wk+1. Thus(81)ECa=∫0tbE0fττ0dτ0+∫tbw1Ebfττ0dτ0+∑j=1k-1∫wjwj+1Ejfττ0dτ0+∫wkwk+1Ekfττ0dτ0+∫wk+1∞Ek+1fττ0dτ0.E0 can be computed by conditioning y:(82)E0=∑y0=1∞ECa∣y=y0,τ=τ0Py=y0∣τ=τ0=∑y0=1kECa∣y=y0,τ=τ0Py=y0∣τ=τ0+∑y0=k+1∞ECa∣y=y0,τ=τ0Py=y0∣τ=τ0.E01 denotes E(Ca∣y=y0,τ=τ0) when 0<τ0<tb,1≤y0≤k. E02 denotes E(Ca∣y=y0,τ=τ0) when 0<τ0<tb,y0≥k+1. Thus(83)E0=∑y0=1kE01Py=y0∣τ=τ0+∑y0=k+1∞E02Py=y0∣τ=τ0.E01 and E02 can be computed by conditioning tr:(84)E01=∫0∞ECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0=∫0H/dECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0+∫H/d∞ECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0=∫0H/dEa1τ0,y0,tr0frtr0dtr0+∫H/d∞Ea2τ0,y0,tr0frtr0dtr0,(85)E02=∫0∞ECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0=∫0H/dECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0+∫H/d∞ECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0=∫0H/dEa3τ0,tr0frtr0dtr0+∫H/d∞Ea4τ0,tr0frtr0dtr0.Eb can be computed by conditioning y:(86)Eb=∑y0=1∞ECa∣y=y0,τ=τ0Py=y0∣τ=τ0=∑y0=1kECa∣y=y0,τ=τ0Py=y0∣τ=τ0+∑y0=k+1∞ECa∣y=y0,τ=τ0Py=y0∣τ=τ0.Eb1 denotes E(Ca∣y=y0,τ=τ0) when tb<τ0<w1,1≤y0≤k. Eb2 denotes E(Ca∣y=y0,τ=τ0) when tb<τ0<w1,y0≥k+1. Thus(87)Eb=∑y0=1kEb1Py=y0∣τ=τ0+∑y0=k+1∞Eb2Py=y0∣τ=τ0.Eb1 and Eb2 can be computed by conditioning tr:(88)Eb1=∫0∞ECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0=∫0H/dECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0+∫H/d∞ECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0=∫0H/dEa5τ0,y0,tr0frtr0dtr0+∫H/d∞Ea6τ0,y0,tr0frtr0dtr0,(89)Eb2=∫0∞ECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0=∫0H/dECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0+∫H/d∞ECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0=∫0H/dEa7τ0,tr0frtr0dtr0+∫H/d∞Ea8τ0,tr0frtr0dtr0.Ej can be computed by conditioning y:(90)Ej=∑y0=1∞ECa∣y=y0,τ=τ0Py=y0∣τ=τ0=∑y0=1k-jECa∣y=y0,τ=τ0Py=y0∣τ=τ0+∑y0=k-j+1∞ECa∣y=y0,τ=τ0Py=y0∣τ=τ0.Ej1 denotes E(Ca∣y=y0,τ=τ0) when wj<τ0<wj+1,1≤y0≤k-j. Ej2 denotes E(Ca∣y=y0,τ=τ0) when wj<τ0<wj+1,y0≥k-j+1. Thus(91)Ej=∑y0=1k-jEj1Py=y0∣τ=τ0+∑y0=k-j+1∞Ej2Py=y0∣τ=τ0.Ej1 and Ej2 can be computed by conditioning tr:(92)Ej1=∫0∞ECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0=∫0H/dECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0+∫H/d∞ECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0=∫0H/dEa9τ0,y0,tr0frtr0dtr0+∫H/d∞Ea10τ0,y0,tr0frtr0dtr0,(93)Ej2=∫0∞ECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0=∫0H/dECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0+∫H/d∞ECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0=∫0H/dEa11τ0,tr0frtr0dtr0+∫H/d∞Ea12τ0,tr0frtr0dtr0.Ek can be computed by conditioning tr:(94)Ek=∫0∞ECa∣tr=tr0,τ=τ0frtr0dtr0=∫0tdECa∣tr=tr0,τ=τ0frtr0dtr0+∫td∞ECa∣tr=tr0,τ=τ0frtr0dtr0=∫0tdEa13τ0,tr0frtr0dtr0+∫td∞Ea14τ0,tr0frtr0dtr0.Ek+1 can be computed by conditioning tp:(95)Ek+1=∫0∞ECa∣tp=tp0,τ=τ0fptp0dtp0=∫0H/dECa∣tp=tp0,τ=τ0fptp0dtp0+∫H/d∞ECa∣tp=tp0,τ=τ0fptp0dtp0=∫0H/dEa15τ0,tp0fptp0dtp0+∫H/d∞Ea16τ0,tp0fptp0dtp0.Combing Equations (81)-(95), it can be obtained that(96)ECa=∫0tbE0fττ0dτ0+∫tbw1Ebfττ0dτ0+∑j=1k-1∫wjwj+1Ejfττ0dτ0+∫wkwk+1Ekfττ0dτ0+∫wk+1∞Ek+1fττ0dτ0where(97)E0=∑y0=1kE01Py=y0∣τ=τ0+∑y0=k+1∞E02Py=y0∣τ=τ0,(98)E01=∫0tdEa1τ0,y0,tr0frtr0dtr0+∫td∞Ea2τ0,y0,tr0frtr0dtr0,(99)E02=∫0tdEa3τ0,tr0frtr0dtr0+∫td∞Ea4τ0,tr0frtr0dtr0,(100)Eb=∑y0=1kEb1Py=y0∣τ=τ0+∑y0=k+1∞Eb2Py=y0∣τ=τ0,(101)Eb1=∫0tdEa5τ0,y0,tr0frtr0dtr0+∫td∞Ea6τ0,y0,tr0frtr0dtr0,(102)Eb2=∫0tdEa7τ0,tr0frtr0dtr0+∫td∞Ea8τ0,tr0frtr0dtr0,(103)Ej=∑y0=1k-jEj1Py=y0∣τ=τ0+∑y0=k-j+1∞Ej2Py=y0∣τ=τ0,(104)Ej1=∫0tdEa9τ0,y0,tr0frtr0dtr0+∫td∞Ea10τ0,y0,tr0frtr0dtr0,(105)Ej2=∫0tdEa11τ0,tr0frtr0dtr0+∫td∞Ea12τ0,tr0frtr0dtr0,(106)Ek=∫0tdEa13τ0,tr0frtr0dtr0+∫td∞Ea14τ0,tr0frtr0dtr0,(107)Ek+1=∫0tdEa15τ0,tp0fptp0dtp0+∫0tdEa16τ0,tp0fptp0dtp0.Similarly, E(Tc) can also be obtained:(108)ETc=∫0tbE0fττ0dτ0+∫tbw1Ebfττ0dτ0+∑j=1k-1∫wjwj+1Ejfττ0dτ0+∫wkwk+1Ekfττ0dτ0+∫wk+1∞Ek+1fττ0dτ0where(109)E0=∑y0=1kE01Py=y0∣τ=τ0+∑y0=k+1∞E02Py=y0∣τ=τ0,(110)E01=∫0tdEc1y0frtr0dtr0+∫td∞Ec2y0,tr0frtr0dtr0,(111)E02=∫0tdEc3frtr0dtr0+∫td∞Ec4tr0frtr0dtr0,(112)Eb=∑y0=1kEb1Py=y0∣τ=τ0+∑y0=k+1∞Eb2Py=y0∣τ=τ0,(113)Eb1=∫0tdEc5y0frtr0dtr0+∫td∞Ec6y0,tr0frtr0dtr0,(114)Eb2=∫0tdEc7frtr0dtr0+∫td∞Ec8tr0frtr0dtr0,(115)Ej=∑y0=1k-jEj1Py=y0∣τ=τ0+∑y0=k-j+1∞Ej2Py=y0∣τ=τ0,(116)Ej1=∫0tdEc9y0frtr0dtr0+∫td∞Ec10tr0frtr0dtr0,(117)Ej2=∫0tdEc11frtr0dtr0+∫td∞Ec12tr0frtr0dtr0,(118)Ek=∫0tdEc13τ0,tr0frtr0dtr0+∫td∞Ec14τ0,tr0frtr0dtr0,(119)Ek+1=∫0tdEc15fptp0dtp0+∫td∞Ec16tp0fptp0dtp0.Thus the values ofE(Ca) and E(Tc) are obtained and the objective function value can be calculated by Equation (63).


### 6.2.2. Comparison of Calculation Methods

In this section, an orthogonal experiment is employed to compare the integral method and the Monte Carlo method. Table3 presents the high and low levels of the 5 input parameters in the orthogonal experiment, which are the decision variables of the objective function. Other parameter values are the same as those in Section 6.1. The L8 orthogonal array is used to carry out the experiment as shown in Table 4. The objective function (ETC) value is calculated by the two methods in every trial. The experiment result is also shown in Table 4. The relative error forETC value of the Monte Carlo method is between 0.0487% and 0.0713%. In addition, the computation speed of the Monte Carlo method is much higher than that of the integral method. So this paper uses the Monte Carlo method to calculate the objective function value in the process of the genetic algorithm.Table 3
Input parameters for orthogonal experiment.
Parameter H h 1 k L n Low 100 50 5 0.1 3 High 200 100 10 0.2 6Table 4
Comparison of the integral and Monte Carlo methods.
No. H h 1 k L n ETC value Time Integral method Monte Carlo method R e l a t i v e   e r r o r Integral method Monte Carlo method 1 100 50 5 0.1 3 650.1860 649.7387 0.0688% 7.538693s 0.061884s 2 100 50 5 0.2 6 650.2032 649.7520 0.0694% 2.378531s 0.094823s 3 100 100 10 0.1 3 578.1220 577.8375 0.0492% 7.272525s 0.117649s 4 100 100 10 0.2 6 578.1217 577.8374 0.0492% 4.626247s 0.060064s 5 200 50 10 0.1 6 641.2674 640.8161 0.0704% 6.499943s 0.070176s 6 200 50 10 0.2 3 641.2382 640.7808 0.0713% 4.793808s 0.058871s 7 200 100 5 0.1 6 577.2176 576.9364 0.0487% 4.002955s 0.054454s 8 200 100 5 0.2 3 577.2226 576.9414 0.0487% 2.219563s 0.052333s


## 6.2.1. Integral Method

The expectation ofCa can be computed by conditioning τ:(80)ECa=∫0∞ECa∣τ=τ0fττ0dτ0=∫0tbECa∣τ=τ0fττ0dτ0+∫tbw1ECa∣τ=τ0fττ0dτ0+∑j=1k-1∫wjwj+1ECa∣τ=τ0fττ0dτ0+∫wkwk+1ECa∣τ=τ0fττ0dτ0+∫wk+1∞ECa∣τ=τ0fττ0dτ0.Let E0 denote the value of E(Ca∣τ=τ0) when 0<τ0<tb. Eb denotes E(Ca∣τ=τ0) when tb<τ0<w1. Ej denotes E(Ca∣τ=τ0) when wj<τ0<wj+1(1≤j≤k-1). Ek denotes E(Ca∣τ=τ0) when wk<τ0<wk+1. Ek+1 denote E(Ca∣τ=τ0) when τ0>wk+1. Thus(81)ECa=∫0tbE0fττ0dτ0+∫tbw1Ebfττ0dτ0+∑j=1k-1∫wjwj+1Ejfττ0dτ0+∫wkwk+1Ekfττ0dτ0+∫wk+1∞Ek+1fττ0dτ0.E0 can be computed by conditioning y:(82)E0=∑y0=1∞ECa∣y=y0,τ=τ0Py=y0∣τ=τ0=∑y0=1kECa∣y=y0,τ=τ0Py=y0∣τ=τ0+∑y0=k+1∞ECa∣y=y0,τ=τ0Py=y0∣τ=τ0.E01 denotes E(Ca∣y=y0,τ=τ0) when 0<τ0<tb,1≤y0≤k. E02 denotes E(Ca∣y=y0,τ=τ0) when 0<τ0<tb,y0≥k+1. Thus(83)E0=∑y0=1kE01Py=y0∣τ=τ0+∑y0=k+1∞E02Py=y0∣τ=τ0.E01 and E02 can be computed by conditioning tr:(84)E01=∫0∞ECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0=∫0H/dECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0+∫H/d∞ECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0=∫0H/dEa1τ0,y0,tr0frtr0dtr0+∫H/d∞Ea2τ0,y0,tr0frtr0dtr0,(85)E02=∫0∞ECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0=∫0H/dECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0+∫H/d∞ECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0=∫0H/dEa3τ0,tr0frtr0dtr0+∫H/d∞Ea4τ0,tr0frtr0dtr0.Eb can be computed by conditioning y:(86)Eb=∑y0=1∞ECa∣y=y0,τ=τ0Py=y0∣τ=τ0=∑y0=1kECa∣y=y0,τ=τ0Py=y0∣τ=τ0+∑y0=k+1∞ECa∣y=y0,τ=τ0Py=y0∣τ=τ0.Eb1 denotes E(Ca∣y=y0,τ=τ0) when tb<τ0<w1,1≤y0≤k. Eb2 denotes E(Ca∣y=y0,τ=τ0) when tb<τ0<w1,y0≥k+1. Thus(87)Eb=∑y0=1kEb1Py=y0∣τ=τ0+∑y0=k+1∞Eb2Py=y0∣τ=τ0.Eb1 and Eb2 can be computed by conditioning tr:(88)Eb1=∫0∞ECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0=∫0H/dECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0+∫H/d∞ECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0=∫0H/dEa5τ0,y0,tr0frtr0dtr0+∫H/d∞Ea6τ0,y0,tr0frtr0dtr0,(89)Eb2=∫0∞ECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0=∫0H/dECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0+∫H/d∞ECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0=∫0H/dEa7τ0,tr0frtr0dtr0+∫H/d∞Ea8τ0,tr0frtr0dtr0.Ej can be computed by conditioning y:(90)Ej=∑y0=1∞ECa∣y=y0,τ=τ0Py=y0∣τ=τ0=∑y0=1k-jECa∣y=y0,τ=τ0Py=y0∣τ=τ0+∑y0=k-j+1∞ECa∣y=y0,τ=τ0Py=y0∣τ=τ0.Ej1 denotes E(Ca∣y=y0,τ=τ0) when wj<τ0<wj+1,1≤y0≤k-j. Ej2 denotes E(Ca∣y=y0,τ=τ0) when wj<τ0<wj+1,y0≥k-j+1. Thus(91)Ej=∑y0=1k-jEj1Py=y0∣τ=τ0+∑y0=k-j+1∞Ej2Py=y0∣τ=τ0.Ej1 and Ej2 can be computed by conditioning tr:(92)Ej1=∫0∞ECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0=∫0H/dECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0+∫H/d∞ECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0=∫0H/dEa9τ0,y0,tr0frtr0dtr0+∫H/d∞Ea10τ0,y0,tr0frtr0dtr0,(93)Ej2=∫0∞ECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0=∫0H/dECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0+∫H/d∞ECa∣tr=tr0,y=y0,τ=τ0frtr0dtr0=∫0H/dEa11τ0,tr0frtr0dtr0+∫H/d∞Ea12τ0,tr0frtr0dtr0.Ek can be computed by conditioning tr:(94)Ek=∫0∞ECa∣tr=tr0,τ=τ0frtr0dtr0=∫0tdECa∣tr=tr0,τ=τ0frtr0dtr0+∫td∞ECa∣tr=tr0,τ=τ0frtr0dtr0=∫0tdEa13τ0,tr0frtr0dtr0+∫td∞Ea14τ0,tr0frtr0dtr0.Ek+1 can be computed by conditioning tp:(95)Ek+1=∫0∞ECa∣tp=tp0,τ=τ0fptp0dtp0=∫0H/dECa∣tp=tp0,τ=τ0fptp0dtp0+∫H/d∞ECa∣tp=tp0,τ=τ0fptp0dtp0=∫0H/dEa15τ0,tp0fptp0dtp0+∫H/d∞Ea16τ0,tp0fptp0dtp0.Combing Equations (81)-(95), it can be obtained that(96)ECa=∫0tbE0fττ0dτ0+∫tbw1Ebfττ0dτ0+∑j=1k-1∫wjwj+1Ejfττ0dτ0+∫wkwk+1Ekfττ0dτ0+∫wk+1∞Ek+1fττ0dτ0where(97)E0=∑y0=1kE01Py=y0∣τ=τ0+∑y0=k+1∞E02Py=y0∣τ=τ0,(98)E01=∫0tdEa1τ0,y0,tr0frtr0dtr0+∫td∞Ea2τ0,y0,tr0frtr0dtr0,(99)E02=∫0tdEa3τ0,tr0frtr0dtr0+∫td∞Ea4τ0,tr0frtr0dtr0,(100)Eb=∑y0=1kEb1Py=y0∣τ=τ0+∑y0=k+1∞Eb2Py=y0∣τ=τ0,(101)Eb1=∫0tdEa5τ0,y0,tr0frtr0dtr0+∫td∞Ea6τ0,y0,tr0frtr0dtr0,(102)Eb2=∫0tdEa7τ0,tr0frtr0dtr0+∫td∞Ea8τ0,tr0frtr0dtr0,(103)Ej=∑y0=1k-jEj1Py=y0∣τ=τ0+∑y0=k-j+1∞Ej2Py=y0∣τ=τ0,(104)Ej1=∫0tdEa9τ0,y0,tr0frtr0dtr0+∫td∞Ea10τ0,y0,tr0frtr0dtr0,(105)Ej2=∫0tdEa11τ0,tr0frtr0dtr0+∫td∞Ea12τ0,tr0frtr0dtr0,(106)Ek=∫0tdEa13τ0,tr0frtr0dtr0+∫td∞Ea14τ0,tr0frtr0dtr0,(107)Ek+1=∫0tdEa15τ0,tp0fptp0dtp0+∫0tdEa16τ0,tp0fptp0dtp0.Similarly, E(Tc) can also be obtained:(108)ETc=∫0tbE0fττ0dτ0+∫tbw1Ebfττ0dτ0+∑j=1k-1∫wjwj+1Ejfττ0dτ0+∫wkwk+1Ekfττ0dτ0+∫wk+1∞Ek+1fττ0dτ0where(109)E0=∑y0=1kE01Py=y0∣τ=τ0+∑y0=k+1∞E02Py=y0∣τ=τ0,(110)E01=∫0tdEc1y0frtr0dtr0+∫td∞Ec2y0,tr0frtr0dtr0,(111)E02=∫0tdEc3frtr0dtr0+∫td∞Ec4tr0frtr0dtr0,(112)Eb=∑y0=1kEb1Py=y0∣τ=τ0+∑y0=k+1∞Eb2Py=y0∣τ=τ0,(113)Eb1=∫0tdEc5y0frtr0dtr0+∫td∞Ec6y0,tr0frtr0dtr0,(114)Eb2=∫0tdEc7frtr0dtr0+∫td∞Ec8tr0frtr0dtr0,(115)Ej=∑y0=1k-jEj1Py=y0∣τ=τ0+∑y0=k-j+1∞Ej2Py=y0∣τ=τ0,(116)Ej1=∫0tdEc9y0frtr0dtr0+∫td∞Ec10tr0frtr0dtr0,(117)Ej2=∫0tdEc11frtr0dtr0+∫td∞Ec12tr0frtr0dtr0,(118)Ek=∫0tdEc13τ0,tr0frtr0dtr0+∫td∞Ec14τ0,tr0frtr0dtr0,(119)Ek+1=∫0tdEc15fptp0dtp0+∫td∞Ec16tp0fptp0dtp0.Thus the values ofE(Ca) and E(Tc) are obtained and the objective function value can be calculated by Equation (63).


## 6.2.2. Comparison of Calculation Methods

In this section, an orthogonal experiment is employed to compare the integral method and the Monte Carlo method. Table3 presents the high and low levels of the 5 input parameters in the orthogonal experiment, which are the decision variables of the objective function. Other parameter values are the same as those in Section 6.1. The L8 orthogonal array is used to carry out the experiment as shown in Table 4. The objective function (ETC) value is calculated by the two methods in every trial. The experiment result is also shown in Table 4. The relative error forETC value of the Monte Carlo method is between 0.0487% and 0.0713%. In addition, the computation speed of the Monte Carlo method is much higher than that of the integral method. So this paper uses the Monte Carlo method to calculate the objective function value in the process of the genetic algorithm.Table 3
Input parameters for orthogonal experiment.
Parameter H h 1 k L n Low 100 50 5 0.1 3 High 200 100 10 0.2 6Table 4
Comparison of the integral and Monte Carlo methods.
No. H h 1 k L n ETC value Time Integral method Monte Carlo method R e l a t i v e   e r r o r Integral method Monte Carlo method 1 100 50 5 0.1 3 650.1860 649.7387 0.0688% 7.538693s 0.061884s 2 100 50 5 0.2 6 650.2032 649.7520 0.0694% 2.378531s 0.094823s 3 100 100 10 0.1 3 578.1220 577.8375 0.0492% 7.272525s 0.117649s 4 100 100 10 0.2 6 578.1217 577.8374 0.0492% 4.626247s 0.060064s 5 200 50 10 0.1 6 641.2674 640.8161 0.0704% 6.499943s 0.070176s 6 200 50 10 0.2 3 641.2382 640.7808 0.0713% 4.793808s 0.058871s 7 200 100 5 0.1 6 577.2176 576.9364 0.0487% 4.002955s 0.054454s 8 200 100 5 0.2 3 577.2226 576.9414 0.0487% 2.219563s 0.052333s


## 6.3. Sensitivity Analysis

In this section, a Design of Experiments (DOE) study is conducted to investigate the sensitivity of the model parameters with the optimal expected cost per unit time (ETC∗).Table5 presents the high and low levels of the 8 input parameters in the DOE, involving production, inventory, maintenance, and quality. Other model parameters are the same as those in Section 6.1. Employing a fractional factorial design, the L16 orthogonal array is used to carry out the orthogonal experiment as shown in Table 6. The GA is used in every trial to obtain the optimal solution of the model. ETC∗ of every trial is also shown in Table 6. In order to analyze the effect of model parameters on ETC∗, the normal probability plot and main effect plot of the model parameters are depicted in Figures 17 and 18 by using Minitab.18. Figures 17 and 18 show that Cw and Cr have significant positive effect on ETC∗.Table 5
Parameter value for DOE.
Parameter U m a x C w C F C V C p C r C y h Low 100 100 1 0.1 700 800 5 0.1 High 200 200 2 0.2 1400 1600 10 0.2Table 6
Parameter values andETC∗ of 16 trials.
No. U m a x C w C F C V C p C r C y h E T C ∗ 1 100 100 1 0.1 700 800 5 0.1 576.5120 2 100 100 1 0.2 700 1600 10 0.2 685.7579 3 100 100 2 0.1 1400 800 10 0.2 578.1663 4 100 100 2 0.2 1400 1600 5 0.1 686.5931 5 100 200 1 0.1 1400 1600 5 0.2 889.3876 6 100 200 1 0.2 1400 800 10 0.1 771.2377 7 100 200 2 0.1 700 1600 10 0.1 889.7835 8 100 200 2 0.2 700 800 5 0.2 771.8040 9 200 100 1 0.1 1400 1600 10 0.1 685.8556 10 200 100 1 0.2 1400 800 5 0.2 577.7022 11 200 100 2 0.1 700 1600 5 0.2 687.0148 12 200 100 2 0.2 700 800 10 0.1 575.9532 13 200 200 1 0.1 700 800 10 0.2 771.9917 14 200 200 1 0.2 700 1600 5 0.1 888.6614 15 200 200 2 0.1 1400 800 5 0.1 771.6710 16 200 200 2 0.2 1400 1600 10 0.2 889.2610Figure 17
Normal probability plot of the effects.Figure 18
Main effects plot.According to the above analysis result, some suggestions for managers are proposed. In order to decrease the expected cost per unit time, the managers could focus on the process of predictive maintenance to decrease its cost per unit time. In addition, the managers could outsource the warranty to decrease the repair cost for warranty.


## 7. Conclusion

The main contribution of this paper is the approach of realizing the Monte Carlo simulation. It speeds up the Monte Carlo simulation dramatically. Thus the direct application of the Monte Carlo method in the process of metaheuristic algorithms becomes practical. As a result, many optimization problems involving complex random variables can be solved expediently. As an application, this paper solves an optimization problem in the manufacturing industry. The problem integrates planning of the HPP and process inspections in a deteriorating manufacturing system, which is practical but has not been studied. Variable sampling intervals are used in the X-bar control chart and the three types of maintenance (CM, PM, and PRDM) are inspired from the industrial needs. The problem is modelled as a 16-scenario stochastic model and the objective function is calculated by the Monte Carlo method with the approach. Subsequently, the optimal solution is found by the genetic algorithm, which is one of the most important metaheuristic algorithms. In order to show the advantage of Monte Carlo method, the method is compared with the traditional integral method. Finally, a sensitivity analysis with an orthogonal experiment is conducted to help managers find the factors with the most significant effect on the expected cost per unit time.Incorporation of other factors such as multiple machines and products can be considered as future extensions of this work. Also, some more realistic assumptions, including several PM levels and limitation of PM manpower, can be incorporated. In addition, the Monte Carlo method with the approach can be used to calculate objective functions in the process of many other metaheuristic algorithms such as the particle swarm algorithm, the tabu search, and the simulated annealing.



---

*Source: 3670495-2019-05-30.xml*