A joint design of production run length, maintenance policy and control chart with multiple assignable causes

Journal of Manufacturing Systems 42 (2017) 44–56

Contents lists available at ScienceDirect

Journal of Manufacturing Systems journal homepage: www.elsevier.com/locate/jmansys

Technical Paper

A joint design of production run length, maintenance policy and control chart with multiple assignable causes Ali Salmasnia∗ , Behnam Abdzadeh, Mohammadreza Namdar Department of Industrial Engineering, Faculty of Technology and Engineering, University of Qom, Iran

a r t i c l e

i n f o

Article history: Received 19 March 2016 Received in revised form 25 October 2016 Accepted 16 November 2016 Keywords: Economic production quantity Statistical process monitoring Maintenance Multiple assignable causes Particle swarm optimization

a b s t r a c t Although economic production quantity, statistical process monitoring and maintenance are three major concepts in process optimization of industrial environments, they have been often investigated separately in literature. Furthermore, in studies that consider these three concepts simultaneously, it is assumed that there is only one assignable cause in the production process. This simplified assumption is unlikely to occur in real production processes due to the usual complexity of manufacturing systems, which may lead to a poor performance in both economic and statistical criteria if the assignable cause originating the shift is different from the one anticipated at the design of the chart. To overcome these drawbacks, this paper develops an integrated model ofeconomic production quantity, statistical process monitoring and maintenance in the presence ofmultiple assignable causes. The particle swarm optimization algorithm is used to minimize the expected total cost per production cycle, subject to statistical quality constraints. Also, a comparative study is performed to illustrate the effect of considering multiple assignable causes on model’s costs. Finally, a sensitivity analysis is conducted on the expected total cost per production cycle and process variable values to extend insights into the matter. © 2016 The Society of Manufacturing Engineers. Published by Elsevier Ltd. All rights reserved.

1. Introduction Rapidly changing markets and the extension of product variety have increased the requirement for more complex equipment. To get the best output of expensive manufacturing processes, and to meet the quality challenge, equipment must be maintained in a suitable operating condition. This circumstance demands more effective maintenance policies. In this situation, the role of the equipment condition in controlling quantity, quality and cost of production becomes more evident and important than ever. Therefore, economic production quantity (EPQ), product quality, and maintenance are three interrelated problems. As mentioned earlier, there is a need for developing approaches that could capture the interdependence between these three main aspects of the modern manufacturing processes. Nevertheless, most of the researchers separately investigated three concepts of inventory, quality, and maintenance. For instance, Chen and Yang [4] presented a model for economic design control chart that only considers statistical properties of the process. Similar to this study, several other papers including Lee et al. [18] and Nenes et al. [23] studied economic design of control charts. Also, there are other papers such as Cheng and Chou [5], De la Torre Gutierrez and Pham [8], and Mehmood

∗ Corresponding author. E-mail address: [email protected] (A. Salmasnia).

et al. [22] that applied control charts with considering special rules like Western Electric rules to detect out-of-control state. Some of the researchers concentrated only on the maintenance issue. Several models were developed for various types of maintenance strategies. In this regard, Lee and Cha [16] implemented periodic preventive maintenance policies for a deteriorating repairable system. They studied detailed properties of the optimal policies that minimize the long-run expected cost. More recent studies about maintenance in industrial environments can be observed in Zhou et al. [34], and Nguyen et al. [24]. The integration of statistical process monitoring (SPM) and maintenance have attracted attention in recent years. In this regard, Linderman et al. [19] constructed a generalized analytic model that incorporates SPM and maintenance policy to minimize total expected cost. Zhou and Zhu [33] developed models that integrated maintenance activities with the design of control chart. They indicated that integrated models perform better than the two stand-alone models. Xiang [31] developed an integrated model of SPM and preventive maintenance (PM) for a manufacturing process that deteriorates according to a discrete-time Markov chain. Moreover, several other researches such as Yin et al. [32], Liu et al. [20], Le and Tan [15], investigated maintenance policy and statistical process monitoring concurrently. Some other approaches integrated SPM and inventory issues. Most of these researches aimed to determine design parameters of control chart in a way that inventory and quality costsbe

http://dx.doi.org/10.1016/j.jmsy.2016.11.003 0278-6125/© 2016 The Society of Manufacturing Engineers. Published by Elsevier Ltd. All rights reserved.

A. Salmasnia et al. / Journal of Manufacturing Systems 42 (2017) 44–56

45

minimized. In this regard, Rosenblatt and Lee [28] presented a comparative study of continuous and periodic inspection policies in deteriorating production systems. They also investigated the effects of imperfect production on the optimal production cycle time. Makis and Fung [21] demonstrated the effect of machine failures on the optimal lot size and the number of samplings in a production cycle. The widespread papers optimized production run-length and maintenance activities, simultaneously. Lee and Rosenblatt [17] conducted study on the joint problem of production planning and maintenance schedule. In this research, it is assumed that the process restoration cost depends on the existence shortages and the delay in detecting them. In another study, Wen et al. [30] suggested an EPQ model for a deteriorating system integrated with predictive maintenance strategy. Gan et al. [11] and Jiang et al. [13] are other studies that have been conducted in this context. Although considering inventory concept in optimizing the production process has a significant role in reducing the manufacturing cost, very few papers in the literature of SPM and maintenance considered this issue. In this regard, Pan et al. [25] presented an integrated EPQ model based on a control chart for an imperfect production process. This model assumes that equipment goes to the out-of-control state because of only a single assignable cause. However, in the real situations, the process might go to an out-of-control state due to different assignable causes. Ben-Daya and Makhdoum [2] developed an integrated production and quality model under various preventive maintenance policies. This model similar to Pan et al. [25] assumes that only one type of assignable cause can occur during the production process. Moreover, this model designs control chart just based on economic criteria, ignoring statistical properties of the process. Rahim and Ben-Daya [27] is the other paper that has been proposed in this context. According to literature review Table 1 summarizes the characteristics of the existing researches in the literature. To fill research gaps and overcome the above-mentioned drawbacks, this study integrates the concepts of SPM, maintenance and inventory in a unified model. It aims to minimize the expected total cost (ETC) of production process subject to statistical constraints. The ETC contains the inventory holding, ordering, maintenance, sampling and quality control costs. Moreover, the proposed model in contrast to the other approaches in this field considers nonuniform sampling in a way that integrated failure rate over all sampling intervals would be the same value. Also to make the model more adapted to real manufacturing situation, the process under consideration can go to out-of-control state due to several types of assignable causes. The rest of this paper is organized as follows: In the next section, the problem will be defined in detail. The Section 3 represents the proposed model for the joint design of production run length, maintenance policy and control chart with multiple assignable causes. In the Section 4, solution approach is explained. The Section 5 which is called experimental results consists of three sub-sections: (1) numerical example, (2) comparative study and (3) sensitivity analysis. Finally, in Section 6 conclusions and further perspectives are given based on the obtained results in the numerical example.

state due to machine deterioration, fatigue and etc. Consequently, the quality loss cost that is imposed to the manufacturer extremely increases because of producing more non-conforming outputs. This paper investigates an imperfect manufacturing process, which includes two states: (1) in-control state and (2) out-ofcontrol state. A Shewhart X-bar control chart is used to monitor a quality characteristic with an alert signal to inform operators when the process shifts to the out-of-control state. The manufacturing process starts operating in the in-control state and one type of assignable cause may occur during the cycle time. An assignable cause happens due to the process nature, machine failure, operator errors, undesirable materials and so on that ultimately leads to shifting process to the out-of-control state. In contrast to most of existing approaches, the suggested model integrates economic production quantity, statistical process monitoring and maintenance in the presence of multiple assignable causes. Furthermore, the proposed model sets sampling intervals in a way that integrated failure rate over all sampling intervals would be the same value. The imperfect process that mentioned earlier consists of three possible scenarios. These scenarios are introduced based on situations that may occur for a production run. Scenario1 happens if production process always remains in an in-control state from the beginning of the production cycle until the end. When scenario1 occurs, aplanned preventive maintenance is implemented at the end of the production cycle. Scenario 2 is a situation in which process shifts during production cycle to an out-of-control state and control chart detects this shift before conducting the planned preventive maintenance. In scenario 2, a corrective maintenance (CM) is required to restore the process to as-good-as-new situation. Scenario 3 is similar to scenario 2 with this difference that control chart can’t detect the shift to the end of the production cycle. In many situations, the alarm wouldn’t be issued at the exact time of shift in the quality characteristic due to a delay in signaling the change by the control chart. In this scenario at the end of the production cycle, when PM activity is conducted, the shift is identified and then PM activities are replaced with CM ones. The graphical representation of three scenarios is presented in Fig. 1 briefly.

2. Problem definition

scale parameter 0 , where 0 =

2.1. Notations Before developing the proposed model mathematically, the notations used to formulate the problem are presented in Table 2. As demonstrated in Table 2, notations are divided into three parts: indices, decision variables and parameters. 2.2. Assumptions and definition The considered assumptions in the proposed model are as follows: 1. Each cycle of the process starts in the in-control state. 2. The time before occurrence of the ith assignable cause follows a Weibull distribution which the probability density function  (PDF) is given by fi (t) = i t −1 e−i t t, i > 0,  ≥ 1, i = 1, ..., s and −1 the hazard rate is ri (t) = i t . 3. The time before occurring first of any assignable cause follows a Weibull distribution with the shape parameter  and the s 

i .This equation is proved in

i=1

Traditional EPQ models focus on inventory issueand make decisions based on storage and ordering costs with this assumption that the manufacturing process is perfect, which means quality defect and machine deterioration never happen. However, a perfect manufacturing process rarely can be found in real industrial environments. A production system often begins its operation in an in-control state and after a period of time, shifts to an out-of-control

Appendix A. f0 (t) = 0 t −1 e−0 t



(1)

4. If there is no alert signal after the kth sampling interval, a preventive maintenance is performed on the systemat the end of the th (k + 1) interval. An alert signal in the jth interval (0 < j < k), indicates that the process has been shifted to the out-of-control state.

46

A. Salmasnia et al. / Journal of Manufacturing Systems 42 (2017) 44–56

Table 1 Summarized literature review. Papers

[18] [23] [8] [4] [24] [33] [34] [21] [28] [22] [30] [19] [11] [32] [17] [15] [16] [31] [5] [20] [13] [25] [27] [2] This paper

IntegratedConcepts

SPM SPM SPM SPM Maintenance SPM/Maintenance Maintenance SPM/EMQ SPM/EMQ SPM EPQ/Maintenance SPM/Maintenance EPQ/Maintenance SPM/Maintenance EMQ/Maintenance SPM/Maintenance Maintenance SPM/Maintenance SPM/Inventory SPM/Maintenance EPQ/Maintenance SPM/Maintenance/EPQ SPM/Maintenance/EPQ SPM/Maintenance/EPQ SPM/Maintenance/EPQ

Assignable Cause

Sampling interval

Single √

Fixed

Multiple √

√

Time to shift distribution

Variable √ √

√

√

√

√

√ √ √

√ √

√

√

√

√

√ √

√ √ √

√ √ √

√

Weibull

√

Weibull √

Exponential

√

Exponential

√ √

Exponential

√

√

√

Weibull Exponential

√ √ √

√

Weibull Weibull Weibull Weibull

√ √ √

Economic/statistical √

Weibull Weibull Weibull

√

√ √

Economic √

Exponential Erlang

√

Type of design

√ √ √

Fig. 1. Graphical representation of three scenarios.

In this case, activities are carried out to find the assignable cause and eventually corrective maintenance is implemented. 5. The production run ends either with a true alarm or after th (k + 1) sampling interval. 6. The production process is monitored by taking samples with size n, at times h1 , h1 + h2 , h1 + h2 + h3 , ..., where hj is the jth sampling interval. Let Wj be the time of drawing jth sample,Wj = j 

hj , j = 1, 2, ... , W0 = 0. Non-uniform sampling, is conducted

i=1

in a way that integrated failure rate over all sampling intervals

would be the same value. As a result, hj , j = 1, 2, ... must satisfy Eq. (2).



h1

Wj+1

ri (t)dt = Wj

ri (t)dt

: ∀j = 1, 2, ...

(2)

0

Therefore, the Wj and hj can be obtained in Eq. (3) that is derived in Appendix B. 1

1

1

Wj = j  × h1 , hj = (j  − (j − 1)  ) × h1

(3)

A. Salmasnia et al. / Journal of Manufacturing Systems 42 (2017) 44–56 Table 2 Notations. Notation

Description

47

According to Eq. (2), hj is a non-increasing function of h1 , therefore h1 ≥ h2 ≥ h3 , ...,. Also Wj tends to infinity as j tends to infinity ( lim Wj = ∞). j→∞

Indices i j z

Index of assignable causes Index of sampling intervals Index of scenarios

decision variables h1 k L n

The duration of first sampling interval The number of sampling intervals in a perfect cycle time Control chart limit The sample size

parameter Ai A ARL0 ARL1j ARL1 ARLl ARLu B CCMi CF CPM CV Cy d D E(CM ) E(CQ ) E(Csampling ) IHC f (t) F(t) IHC p Pij P(SZ ) qij Q Qin Qouti rin routi s t T1 Tin Touti u0 Wj ˛ ˇi ıi i 0   i

The ith type of assignable cause The ordering cost per unit The average run length when the process is in the in-control state The average run length when ith assignable cause occurs The expected value of ARL1j The lower bound of ARL0 The upper bound of ARL1 The inventory holding cost per unit per year The cost of corrective maintenance activity when ith assignable cause occurs The fixed cost of sampling The cost of preventive maintenance activity The variable cost of sampling The cost of each false alarm The daily demand rate The annual demand rate The expected maintenance cost per production cycle The expected quality control cost per production cycle The expected sampling cost per production cycle The expected total cost per production cycle The time to shift probability density function The time to shift cumulative distribution function The inventory holding cost per production cycle The Production rate The conditional Probability that the ith assignable cause will occur during the jth sampling interval The probability of happening z th scenario The unconditional Probability that the ith assignable cause will occur during the jth sampling interval The economic production quantity The quality loss cost per unit when the process is in the in-control state The quality loss cost per unit when ith assignable cause occurs The expected number of sampling points when the process is in the in-control state The expected number of sampling points when the process is out-of-control due to occurrence of Ai Number of assignable causes The time to shift The time to detection and validation of the assignable cause The time that the process is in the in-control state in cycle time The time that the process is in the out-of-control state in cycle when ith assignable cause occurs The upper bound of sample size The time of jth sampling (the end of jth interval) Pr (exceeding control limits | process is in-control) Pr (not exceeding control limits | process is out-of-control due to the ith assignable cause occurred) The magnitude of the shift in the mean of quality characteristic when ith assignable cause occurs Scale parameter of the weibull distribution when ith assignable cause occurs The process mean in-control state Shape parameter of the weibull distribution The standard deviation of quality characteristic The expected in-control time during each sampling interval in which the Ai occurs

7. The following three can be neglected: (1) time to implement preventive and corrective maintenance, (2) time to find a false alarm, and (3) the time of sampling. 8. The quality characteristic is normally distributed. The occurrence of ith assignable cause leads to a mean shift from 0 to either 0 + ıi  or 0 − ıi .Since the probability of detecting positive and negative shifts are the same, this study only considers positive ones. 9. It is assumed that occurrence of assignable cause does not have effect on variance of quality characteristic. In other words, remains unchanged. 10. Although several types of assignable causes can be occurred in different cycles, when an assignable cause occurs in a cycle, no further assignable causes will happen in the same cycle. 11. The assignable causes are independent of each other, and each one occurs randomly. 12. It is perceived that the corrective maintenance restores the process to the as-good-as-new condition. When the process is identified out-of-control, the PM is inadequate and therefore is replaced by corrective maintenance. To derive the expected total cost in our model, the following terms are needed: a) Defining Pij (i = 1, 2, . . . , s , j = 1, 2, . . . , k + 1) as a conditional probability that the ith assignable cause will occur during the sampling interval hj , given that cause Ai does not happen until time Wj−1 .

 Wj Pij = P(Wj−1 < t < Wj |t > Wj−1 ) = 

=

e−i Wj−1 − e e

−i Wj−1



−i W  j

=1−e

−0

W

 ∞j−1

fi (t)dt

f (t)dt Wj−1 i

−i (w −w ) j

j−1

=1−e

−i h

1

(4)

Eq. (4) indicates that Pij is just a function of i , h1 and . Consequently, Pij = Pi for i = 1, 2, . . . , s , j = 1, 2, . . . , k + 1. b) Defining qij as an unconditional probability that the cause Ai will occur during the sampling interval hj .

Wj fi (t)dt = (1 − Pi )j−1 × Pi i = 1, 2, ..., s, j = 1, 2, ..., k + 1 (5)

qij = Wj−1

c) Let  be the expected duration of the in-control period within ij sampling interval hj , given that cause Ai has occurred during the jth sampling interval.

 Wj Wj−1

ij =

(t − Wj−1 ) × fi (t)dt (6)

qij

Afterwards, in each sampling interval that Ai occurs, the expected in-control time (i ) is given by: i =

k 

ij × qij =

j=1

k   j=1

Wj

(t − Wj−1 ) × fi (t)dt

(7)

Wj−1

c)

Defining P0j as the conditional probability that one of assignable causes will occur in the jth sampling interval, given that the process is in-control at time Wj−1 . Similar to Eq. (4), P0j = − h

1 − e 0 1 is only a function of 0 , h1 and , therefore P0j = P0 for j = 1, 2, . . . , k + 1.

48

A. Salmasnia et al. / Journal of Manufacturing Systems 42 (2017) 44–56

Defining q0j as an unconditional probability that the one of assignable causes will happen during the sampling interval. f0 (t)dt = (1 − P0 )j−1 × P0 , j = 1, 2, ..., k + 1

3. The proposed model In contrast to the perfect production process that is often used in the classic EPQ models; this paper suggests an imperfect process with multiple assignable causes, which consists of three possible scenarios. In such a process, the sampling, quality and the maintenance costs are considered in addition to the inventory holding cost (IHC) and ordering cost. Thus, this section contributes to describing the three scenarios in detail and formulating the cycle time and related costs of each scenario with considering multiple assignable causes. Scenario 1: The production process always remains in the in-control state from the start of the production cycle to the end. In other words, total time of being out-of-control is zero. In this situation, the preventive maintenance is implemented at the end of th (k + 1) sampling interval to guarantee the reliability and sustainability of the manufacturing process as illustrated in Fig. 2. 1 

E(Tin |S1 ) = Wk+1 = (k + 1) × h1

(10)

Scenario 2: As shown in Fig. 3, the process starts in the in-control state and is monitored by taking samples with size n at times h1 , h1 + h2 , . . . to check the status of the mean process. The process shifts to the th out-of-control state at a time between the jth and (j + 1) sampling due to occurrence of ith assignable cause Ai (i = 1 , 2 , . . . , s),and consequently the mean value of the quality characteristic changes from 0 to 0 + ıi 0 . Unfortunately, the control chart cannot respond to the occurrence of assignable cause immediately at th the (j + 1) sampling time, hence the process continues until the th (j + y) sampling that alert signal is sent. At this moment, the operators search to detect the assignable cause and conduct the corrective maintenance on the system. Then, the process is restored to as-good-as-new situation. In this scenario, the cycle time includes both the in-control and the out-of-control states in which the in-control time is assumed to follow a truncated Weibull distribution: f0 (t) f0 (t) = p(0 < t < Wk+1 ) F0 (Wk+1 )

f0 (t) 1−e

(11)

−0 W 

Therefore, the in-control time can be calculated by integrating this function: Wk

E(Tin |S2 ) =

t × f0 (t|0 < t < Wk+1 )dt

(12)

0

The out-of-control time consists of the following periods: (1) the time between the occurrence of the assignable cause Ai (i = 1 , 2 , . . . , s) and the next inspection time; (2) the time to trigger an alarm signal and; (3) the time to detect and validate the assignable cause. Therefore, the out-of-control time after happening Ai can be written as:  E(Touti |S2 ) =

k+1−j k  

( j=1

y=1



× E(Touti |S2 )

(14)

Scenario 3: As shown in Fig. 4 and similar to scenario 2, the process starts in the in-control state and shifts to out-of-control state at a time th between the jth and (j + 1) sampling due to occurrence of ith assignable cause Ai (i = 1 , 2 , . . . , s) that leads to shifting the mean of quality characteristic from 0 to 0 + ıi 0 . However, because of the restrictions of the control chart capability, the process continth ues with no signal alert until the end of (k + 1) sampling interval that assignable cause is detected during maintenance implementation. Thus the planned preventive maintenance has to be replaced by corrective maintenance to deal with this much tougher situation. The system is restored after maintenance implementation to as-good-as-new situation. Difference between scenarios 2 and 3 causes changing the upper limit of the integral from Wk to Wk+1 for computing expected incontrol time:



Wk+1

E(Tin |S3 ) =

t × f0 (t|0 < t < Wk+1 )dt

(15)

0

Therefore, the out-of-control time is obtained by subtracting the in-control time from the production time: 1

E(Tout |S3 ) = Wk+1 − E(Tin |S3 ) = (k + 1)  × h1 − E(Tin |S3 )

(16)

According to the defined scenarios, the occurrence probability of each scenario can be obtained as follows: Scenario 1 occurs when process remains in-control in k + 1 inspection intervals. Therefore,P(S1 ) is defined as the probability of occurrence a shift after at least k + 1 inspection interval (Eq. (17)). P(S1 ) = 1 − F0 (Wk+1 ) = e

−0 (k+1)h

(17)

1

P(S2 ) is equal to the conditional probability of occurring a shift before kth sampling given that the shift is detected before (k + 1)th  sampling. It must be noticed that  i is the probability of occurring 0

ith assignable cause and 1 − ˇik−1 is the probability of triggering an alarm signal when the process shifts to out-of-control state because of ith assignable cause. S  

i

P(S2 ) = F0 (Wk ) × (

i=1

0

× (1 − ˇik−1 ))

(18)

P(S3 ) is the probability of happening a shift before kth sampling while it is not detected until the end of cycle. Consequently, the occurrence probability of scenario 3 can be obtained as: P(S3 ) = P(0 < t < Wk+1 ) − P(S2 ) = F0 (Wk+1 ) − F0 (Wk )

k+1



i=1

0

(9)

E(Tout |S1 ) = 0

f0 (t|0 < t < Wk+1 ) =



i

(8)

Wj−1

=

s  

E(Tout |S2 ) =

Wj q0j =

As a result, the expected out-of-control time in this scenario is:

y−1

qij × (Wy+j−1 − Wj−1 ) × ˇi

× (1 − ˇi )) − i

+ T1

(13)

S  

i

×(

i=1

0

× (1 − ˇik−1 ))

(19)

With consideration ofthe cycle time and the probability of each scenario given, the sampling cost, the quality loss cost, the maintenance cost, the inventory holding cost and the ordering cost are defined as below. Fig. 5 shows the cost structure of this model. 3.1. Quality loss cost The quality loss cost is imposed to the manufacturer in both in-control and out-of-control states. However, it is an obvious fact that the quality loss cost extremely increases when the process goes

A. Salmasnia et al. / Journal of Manufacturing Systems 42 (2017) 44–56

49

Fig. 2. Graphical representation of scenario1 in the production cycle time.

Fig. 3. Graphical representation of scenario 2 in a production cycle.

Fig. 4. Graphical representation of scenario 3 in process cycle.

to out-of-control state because the probability of producing nonconforming items increases. To calculate the quality loss cost when the process is in-control (out-of-control), the expected number of products that are manufactured in the in-control (out-of-control) state must be multiplied by the Qin (Qouti ). Eventually, the expected quality loss cost is formulated in Eq. (23).

1

E(CQ |S1 ) = Qin × p × E(Tin |S1 ) = Qin × p × (k + 1)  × h1

E(CQ |S2 ) = Qin × p × E(Tin |S2 ) + p ×

s  

i

i=1

0

(20)

× (Qouti × E(Touti |S2 )) (21)

50

A. Salmasnia et al. / Journal of Manufacturing Systems 42 (2017) 44–56

Fig. 5. The cost structure of this model.

E(CQ |S3 ) = Qin × p × E(Tin |S3 ) + p × E(Tout |S3 ) ×

s  

i

i=1

0

× (Qouti ) (22)

E(CQ ) =

3 

E(CQ |Sz ) × P(Sz )

(23)

z=1

to out-of-control state in the jth interval, the number of sampling points when the process is in-control is equal to j-1. Therefore rin is computed as Eq. (28). If process shifts to an out-of-control state in jth interval, the number of sampling points in the out-of-control state can be obtained by subtracting j from sampling number in which a true signal is issued. Therefore routi can be calculated by Eq. (29).

rin =

k 

(j − 1)q0j

(28)

j=1

3.2. Sampling cost To compute the sampling cost in each scenario, the summation of the fixed and variable costs of each sampling (CF + (n × CV )) must be multiplied by the average number of sampling points. The average number of sampling points in scenarios 1 and 3 is equal to k because in these scenarios the production cycle continues to end (Eqs. (24) & (26)), whereas in scenario 2 production cycle is stopped when a true signal is issued. Therefore, the average number of sampling points in a given cycle can be calculated by summation of the expected number of sampling points in both out-of-control and incontrol states (Eq. (25)). Finally the total expected sampling cost per production cycle is accounted as Eq. (27). E(CSampling |S1 ) = (CF + (n × CV )) × k

(24)

E(CSampling |S3 ) = k × (CF + (n × CV ))

(25)

E(CSampling |S2 ) = (rin +

s  

i

i=1

E(CSampling ) =

3 

0

× routi )(CF + (n × CV ))

E(Csampling |Sz ) × P(Sz )

(26)

(27)

z=1

As mentioned in Table 2, rin is the expected number of sampling points when the process is in the in-control state and routi is the expected number of sampling points when the process is out-ofcontrol due to the occurrence of Ai . In the scenario 2 rin and routi are computed by Eqs. (28) and (29). It is obvious when the process goes

k  

routi = (

j=1



k+1−j

wj

(

fi (t)dt))(

wj−1

r(1 − ˇi )ˇi

r−1

)

(28)

r=1

3.3. Maintenance cost Maintenance cost consists of the false alarm cost, preventive and corrective maintenance costs and also depends on the scenario that happens. Since in the scenario1, the process is in-control throughout the production cycle, the corrective maintenance cost is definitely zero. However in the scenarios 2 and 3, maintenance cost entails corrective maintenance cost because the production process shifts to the out-of-control state. Considering the above descriptions, the maintenance cost in scenario1 is calculated by summation of the preventive maintenance cost and false alarm cost (Eq. (30)). To compute false alarm cost, the expected number of false alarms must be multiplied by the cost of each false alarm (Cy ). The expected number of false alarms depends on the average numbers 1 of sampling points (k) and the probability of type I error (˛ = ARL ). 0 In scenarios 2 and 3, the maintenance cost is formulated similar to scenario 1 with this difference that k and the PM cost replaced with rin and the CM cost, respectively (Eqs. (31) & (32)). With respect to the mentioned explanations, the maintenance cost can be obtained as Eq. (33). E(CM |S1 ) =

k × Cy + CPM ARL0

(30)

A. Salmasnia et al. / Journal of Manufacturing Systems 42 (2017) 44–56

rin × Cy  i + × CCMi ARL0 0

(31)

rin × Cy  i + × CCMi ARL0 0

(32)

s

E(CM |S2 ) =

i=1 s

E(CM |S3 ) =

i=1

E(CM ) =

3 

E(CM |Sz ) × P(Sz )

(33)

z=1

It should be mentioned that rin in scenario 3 is a little bit different from scenario 2, because in scenario3 process can shift to out-ofcontrol in the (k + 1)th interval and as a result the upper limit of summation changes to k + 1 (Eq. (34)). rin =

k+1 

(j − 1)q0j

(38)

4.1. PSO description

According to the classic EPQ model, inventory holding cost depends on inventory holding cost per unit per time unit (B), the cycle time (T), production(P) and demand rate(d). Thus inventory holding cost is formulated as follows: (35)

3.5. Ordering cost If the ordering cost per unit and the annual demand rate are depicted by A’ and D, respectively. The total ordering cost is calculated according to Eq. (36), considering classic EPQ model. Ordering Cost =

D × A D × A = P×T P × Wk+1

(36)

Eventually, the expected total cost and economic production quantity (EPQ) can be computed as demonstrated in Eqs. (37) and (38): ETC =

B × (P − d) × Wk+1 D × A + + E(CQ ) + E(CSampling ) + E(CM ) P × Wk+1 2 1

EPQ = Wk+1 × P = (k + 1)  × h1 × P 3.6. The objective function and constraints

In this section, with respect to costs that were described earlier, the objective function and constraints related to the model are explained as follows: MinETC s.t :

(39)

1 ≤ n ≤ u0

(39.1)

ARL0 > ARLl

(39.2)

ARL1 < ARLu

(39.3)

k ≥ kl

(39.4) +

h > 0, L > 0, n ∈ N , k ∈ N

+

4. Solution approach

(37)

3.4. Inventory holding cost

B × Wk+1 × (P − d) B × T × (P − d) = 2 2

1 Typically, in practical conditions because of economic reasons, the sample size must be less than a definite number u0 as shown in Eq. (39.1). 2 In order to decrease the number of false alarms without affecting the performance of control chart, ARL0 must be bigger than a predetermined value of ARLl , as illustrated in Eq. (39.2). 3 The control chart should detect a shift in the process as soon as possible. This occurs when ARL1 is less than a pre-determined value of ARLu , as shown in Eq. (39.3). 4 In order to assure process continuity, the number of sampling intervals in a perfect cycle must be bigger than a pre-determined value of kl . Eq. (39.4) shows this constraint mathematically.

The model that is explained in the last section has several complexities that prevent the model to be solved with exact methods. Three of the most important complexities are: (1) The mathematical model (39) includes both continuous and discrete decision variables; (2) it is obvious that solution space is discontinuous and non-convex; and (3) in the objective function, some decision variables are in limits of an integral or in the cumulative density function (CDF) of a standard normal distribution. One of the best ways to solve complex problems is using metaheuristic algorithms to achieve favorable results in a reasonable time. In the existing literature, there are many papers that used meta-heuristic algorithms for solving models similar to the proposed model in this research. For instance, Chih et al. [7] used PSO for an economic and economic statistical design of X-bar control chart. Faraz and Saniga [10] applied Genetic Algorithm (GA) for the economic statistical design of a T2 control chart with double warning lines. Also, Bashiri et al. [1] utilized multi-objective GA for an economic statistical design of X-bar control chart. In this study, Particle Swarm Optimization (PSO) is employed for optimizing the presented mathematical model (39) because of its good performance in optimizing non-linear models, unique searching mechanism, simplicity in concept, computational efficiency and ease in implementation [29]. It is one of the most popular metaheuristic algorithms that recently has been widely applied by some researchers such as Perez and Behdinan [26], Cheng et al. [6], and Ghodratnama et al. [12].

(34)

j=1

IHC =

51

(39.5)

As mentioned in Eq. (39), the objective function is equal to ETC that represents the total expected cost per production cycle. In addition, the necessary constraints (39.1–39.5) should be imposed in the mathematical programming to make the model more adapted to real manufacturing situation. These constraints are as follows:

PSO is a population-based meta-heuristic search algorithm, inspired by the social behavior of bird flocking or fish schooling. This algorithm utilizes collaboration among a population called particles, to find optimal decision variable values in a solution space. In PSO, particles fly through the problem space by following the current optimum particles. Each particle has a fitness value that is evaluated by the objective function and also has a velocity vector that determines the direction of the particle. PSO incorporates local and global searches to attain high search efficiency. This algorithm is initialized with a swarm of particles that have random positions and velocities. Then the algorithm searches for the optimal solution by updating particles based on the force of inertia and the two “best” values. Of these two values, the first one is considered as the best solution observed so far which is called global best (gbest). Another “best” value is the best value experienced by the ith particle which is called personal best (pbesti ). In other words, in this iterative process, the behavior of a particle is a compromise among three possible alternatives: (1) following its current velocity (2) going towards its personal best (3) going towards the global best. In each iteration, after finding the two best values, the particle will update its velocity and position.

52

A. Salmasnia et al. / Journal of Manufacturing Systems 42 (2017) 44–56

To summarize the PSO algorithm, it is assumed that N is the number of particles in the swarm.xit =

[n, h, l, k],vti and pbestit−1 represent the position(the vector of

decision variables), velocity and personal best of the ith particle in the iteration t. blo and bup are the lower and upper boundaries for decision variables, respectively. gbest t is considered as the best solution founded till iteration t. The parameters c1 and c2 are cognition and social learning factors, where the value of (c1 + c2 ) according to Kennedy et al. [14] is usually considered equal to 4. Eventually, w represents inertia weight that its value varies in each iteration. The value of inertia weight in the first iteration is a pre-determined value between 0 and 1. In each iteration, w is reduced to w × wdamp , where wdamp is a fixed factor definitely less than 1. The mentioned formulation for w leads to an increase in intensification of algorithm. • For each particle i = 1, . . ., N do:  Generate the initial value of each particle position with a uniformly distributed random vector: xit ∼U(blo , bup ), where blo and bup are the lower and upper limits of the search space, respectively. • Initialize the pbest of each particle to its initial position: pbesti → xi . • If (ETC(pi ) ≤ ETC(g)) update the gbest as: gbest → pbesti . • Generate the initial value of each particle velocity: vi ∼U(−|bup − blo |, |bup − blo |). • Until a stop criterion is met (e.g. number of iterations performed, or a solution with adequate objective function value is found), repeat:  For each particle i = 1, . . ., N do:   For each dimension d = 1, . . ., n do:  Generate random numbers: randp , randg ∼U(0, 1) • Update the particle velocity: vt = w.vt−1 + c1 randp (pbest t−1 − i i i • • • •

xit−1 ) + c2 randg (gbest t−1 − xit−1 )

Update the particle position: xit = xit−1 + vti If (ETC(xi ) ≤ ETC(pbesti )): Update the pbest of each particle If (ETC(pbesti ) ≤ ETC(gbest)) gbest Now gbest holds the best found solution.

In the rest of this section, a representation of a searching point by PSO algorithm in a two-dimensional feasible space is illustrated in Fig. 6 and the computational procedure of the PSO algorithm is summarized in Fig. 7.

4.2. Application of PSO for solving the proposed model The solution representation is a key factor in developing the PSO algorithm which could be a string of both integers and real numbers. The solution representation for the proposed model includes four-dimensional particles that each dimension refers to a certain decision variable. In the presented model, the sample size (n) and the number of sampling intervals in a perfect cycle time (k) are an integer, while the duration of first sampling interval (h1 ) and control chart limit (L) are real numbers. As mentioned in pseudo code presented for PSO earlier, for producing the initial value of each continuous decision variable, a uniformly distributed random value is generated between lower and upper limits of the considered decision variable. Furthermore, in order to generate an initial value for discrete decision variables, a random value from a uniform distribution in the interval [0,1] is generated. The values of discrete variables, i.e. the sample size (n)

and the number of sampling intervals in a perfect cycle time (k), are obtained according to Eqs. (40) and (41), respectively. n = min((nmin + floor((nmax − nmin + 1) × R1 )), nmax )

(40)

k = min((kmin + floor((kmax − kmin + 1) × R2 )), kmax )

(40)

where, kmin , kmax , nmin and nmax are the lower and upper limits of k and n, respectively. Furthermore, R1 and R2 are two random number that follows uniform distribution with R1 , R2 ∼U(0, 1) 5. Experimental results The purpose of this study is to minimize the ETC per production cycle, subject to the economic and statistical constraints. In this section, a numerical example is presented to demonstrate the applicability of the proposed model. The remainder of this section is divided into three sub-sections as follows: Section 5.1 presents the mathematical programming and the assigned values to model parameters. Then, PSO algorithm is used to optimize the presented mathematical programming and the obtained results are recorded. Section 5.2 represents a comparative study conducted between single and multiple assignable cause models. Finally in Section 5.3, a sensitivity analysis is implemented on three parameters: 1) the shape parameter of time to shift distribution, 2) the false alarm cost; and 3) the quality loss cost per unit when process is in the in-control state. 5.1. Numerical example In this section, an industrial example modified from Pan et al. [25] and Chen and Yang [4] is presented to demonstrate the application of the suggested model: a company sells a certain food product to a wholesaler in packages marked with the specific weight. According to Bisgaard et al. [3] it is assumed that the quality characteristic of the process is the weight of packages and assignable causes can only change the mean weight of packages. The values of main parameters related to manufacturing system are recorded in Table 3. Moreover, to make the model more suitable for real situation production, the constraints should be considered in the mathematical model. Generally, with respect to economic and statistical considerations, factors such as the sample size, ARL and the number of sampling intervals should be limited. Thus the proposed mathematical programming can be rewritten with respect to u0 = 20,ARLl = 100,ARLu = 10 and kl = 40. Min ETC 1 ≤ n ≤ 20 ARL0 > 100 ARL1 < 10 k ≥ 40 h > 0, L > 0, n ∈ N + , k ∈ N + The values assigned to the problem parameters are shown in Tables 3 and 4. As illustrated in Table 4,cri ,i and Qouti are presumed to be functions of ıi . When the ith assignable cause occurs,0 shift to 0 + ıi .i is a non-increasing function of ıi for i = 1, . . . , s.Qouti is proportional to the shift value that occurs in the product quality characteristic (ˇ0i = 1 − (3 − ıi ) + (−3 − ıi ) ≈ 1 − (3 − ıi ) for i = 1, . . . , s). In a manufacturing process, the greater shift makes a greater value for cri . It is assumed that there are seven types of assignable causes that (Ai , i = 1, 2, ..., 7) lead to 1, 1.5, 1.8, 2, 2.2, 2.5 and 3 shift value, respectively. It is considered that the assignable causes occur randomly and independently in the production process cycle. PDi indicates a prior distribution of ıi (ı1 = 1, ı2 = 1.5, ı3 = 1.8, ı4 =

A. Salmasnia et al. / Journal of Manufacturing Systems 42 (2017) 44–56

53

Fig. 6. Representation of a searching point by PSO algorithm [26].

Fig. 7. Complete computational procedure of the PSO algorithm [9].

Table 3 The parameter values. Parameter Value

T1 1.25

CF 20

Qin 210

CV 5

CPM 2400

CY 200

parameter Value

 2

s 7

S 60

B 10

D 10000

d 80

P 100

Table 4 The set of cost and probability parameters. Ai

ıi

ˇ0i

PDi

cri

Qouti

i

1 2 3 4 5 6 7 Base value

1 1.5 1.8 2 2.2 2.5 3 ı4 = 2

0.5000 0.3085 0.2119 0.1587 0.1151 0.0668 0.0228 ˇ04 = 0.1587

3.6650 8.7923 6.7426 4.2652 2.9151 5.3230 4.3512 PD4 = 4.2652

5000 11994.95 9198.636 5818.827 3976.944 7261.937 5936.153 cr1 = 5000

71.7964 210.5524 362.6424 500 667.6587 972.3504 1575.7 Qout4 = 500

0.02 0.04798 0.036795 0.023275 0.015908 0.029048 0.023745 i = 0.02

54

A. Salmasnia et al. / Journal of Manufacturing Systems 42 (2017) 44–56

2, ı5 = 2.2, ı6 = 2.5 and ı7 = 3). In this study PDi follows normal distribution. Since ˇ0i is in the reverse proportion to ıi and Qouti , consequently Qouti is proportional to ıi . Therefore Qout4 can be set as the “base case”, and for other cases Qouti is calculated proportion to Qout4 by Qouti = (ˇ0i /ˇ04 ) × Qout4 . In a similar way,1 can be set as the base case, and hence i = (PDi /PD1 ) × 1 and cri = (PDi /PD4 ) × cr1 . The value of cri , i and Qouti for different ıi are listed in Table 4. According to the characteristics of the model and based on criteria and requirements mentioned earlier in Section 4, PSO is coded and implemented in MATLAB2015. The optimal solution is attained as the follow:



L∗ , h∗1 , n, k ∗

⇒ EPQ =



W∗

=

k+1





3.6163, 0.4, 5, 40 , ETC ∗ = 93161 1 1 × P =  × h∗1 × P = 2 × (0.4) × (100) = 256.1

It is suggested that, during the production cycle, the first sampling with a sample size of 5 should be performed after spending 0.4 h. After 40 samplings, the preventive maintenance should be implemented, if the process remains in the in-control state throughout the entire cycle. Also, the control limit should be set at 3.6163 (this is typically set as 3 in the industrial applications). In this model, the organization not only satisfies the average daily (annual) demand, but also gains remarkable profit with a related production cost of $93161, which is optimal.

50 45 40 35 30 value 25 20 15 10 5 0

n L 11

9

3.5956 0.5

0

3.4419

0.4 0.4 1.5 2 Shape Parameter

1

h1

8 3.1128 0.4 2.5

k

3

Fig. 8. Effect of the shape parameter on the design parameters L, h1 , k and n. 50 45 45 40

44

40

35 30

n

value 25 20 15 11

L 9

10

5 3.5321 0 0.4 0

h1

8

3.4419

k

3.1618

0.4 200

400

600

800

0.401 1000

1200

the false alarm cost Fig. 9. Effect of Cy on the designparameters L, h1 , k and n.

5.2. Comparative study

50 45

To validate the effectiveness of the presented model, it is compared with the single assignable cause model in which ı,, Qout and cr are set to be 2, 0.02, 500 and 5000 respectively. As it can be inferred, these parameters are the same in value for ıi = 2 in the multiple assignable cause model. The results of solving both single and multiple assignable cause models using PSO algorithm are compared in Table 5. Table 5 indicates that the proposed multiple assignable cause model not only has less expected total cost in comparison with the single assignable cause model, but also results in better values for the statistical properties of ARL0 and ARL1 . In other words, the obtained results support this claim that the presented model is more efficient than the single assignable cause model.

40

5.3. Sensitivity analysis

43

40

40

46 40

40

35 30

n

value 25 20 15

L

10 5 0

4 3.7521 0

100

0.4 200

150

k 3.2622

3.4419

0.4 50

h1

10

9

250

0.4 300

350

Qin

Fig. 10. Effect of Qin on the design parameters L, h1 , k and n. 475.6

2,40,000

500 450 400

184420

1,90,000

350 300

256.1

In this section, the effect of parameters Cy (the false alarm cost), Qin (the quality loss cost per unit when the process is in the in-control state) and (the shape parameter of Weibull distribution) on the ETC, EPQ and the decision variables of the model are investigated. The variations in the ETC, EPQ and process variables are analyzed by changing one of the mentioned parameters (Cy ,Qin and ) while the other parameters remain fixed. The obtained results from sensitivity analysis are shown in Table 6 and Figs. 8–13 . It must be noticed that the most important item in sensitivity analysis is ETC that directly affects the profit of a company and its sensitivity to the parameters can play a substantial role in making decisions by a manufacturer. 1 According to the results obtained in Table 6, ETC∗ and EPQ ∗ are deeply influenced as the shape parameter value varies. It can be observed in Fig. 11 that with an increase in , the values of ETC∗ and EPQ ∗ decrease dramatically. For example, when  varies from 1.5 to 2, the optimal value of ETC∗ reduces 115000 units. On the other hand it is obvious that the decision variables are almost insensitive to the shape parameter of Weibull distribu-

value 1,40,000

250

180.1 98196

90,000

69214

200

ETC

150

EPQ

100 50

40,000

0

0.5

1

1.5 Shape parameter

2

2.5

3

0

Fig. 11. Effect of the shape parameter on the EPQ* and ETC*. 120000

106640

110550

98196

100000 80000 value 60000

ETC

40000

EPQ

20000 0

271.3 0

269

256.1 200

400

600 800 Cost of false alarm

1000

Fig. 12. Effect of Cy on the EPQ* and ETC*.

1200

A. Salmasnia et al. / Journal of Manufacturing Systems 42 (2017) 44–56

55

Table 5 The optimal results of single and multiple assignable cause models. Model

Optimal Results

Single AC Multiple AC

L∗

h∗1

k∗

n∗

ARL0

ARL1

ETC ∗

2.5 3.6163

2 0.4

40 40

4 5

657.61 23952

2.731 2.713

140050 93161

Table 6 The parameter sensitivity analysis. Parameter

Value

L∗

h∗1

k∗

n∗

ETC ∗

EPQ ∗



1.5 2 2.5

3.5956 3.4419 3.1128

0.4 0.4 0.4

40 40 43

11 9 8

184420 98196 69214

475.6 256.1 180.1

CY

20 200 1000

3.5321 3.4419 3.1618

0.4 0.4 0.401

45 40 44

11 9 8

110550 98196 106640

271.3 256.1 269

Qin

50 210 300

3.7521 3.4419 3.2622

0.4 0.4 0.4

46 40 40

4 9 10

40204 98196 135800

274.2 256.1 256.1

160000

respect to an increase in Cy is quite similar to positive shift in . It means h1 is completely insensitive to shift in Cy while n and L have a low sensitivity to Cy .

135800

140000 120000

98196

100000 value 80000

ETC

60000

40204

EPQ

40000 20000 0

274.2 0

50

256.1 100

150

200

256.1 250

300

350

Qin

Fig. 13. Effect of Qin on the EPQ* and ETC*.

tion (). Fig. 8 illustrates as the parameter  increases from 1.5 to 2 and from 2 to 2.5, the optimal value of h1 remains completely unchanged, while the L and n reduce mildly. On the other hand, it is known that a decrease in the control chart limit leads to improvement of type II error probability. As a result, simultaneous reduction of sample size and control chart limit can enhance both economic and statistical criteria when the shape parameter increases. The design variable k is also insensitive to  and it seems that the observed small variation in k is only because of random nature of PSO algorithm. 2 As expected, when Qin increases, the expected total cost grows extremely. For example as illustrated in Fig. 13, an increase in Qin from 50 to 300 leads to approximately 95000 unit increment in ETC*. Furthermore since EPQ* is a function of k and h1 , it remains almost unchanged while Qin increases. The depicted diagram in Fig. 10 shows that h1 is insensitive to changes in Qin . Because this variable is dependent on factors that reflect the current state of the process, such as ıi . Also the results show that k has little sensitivity to Qin . Eventually, it is clear that n and L with respect to an increase in Qin follow increasing and decreasing trends, respectively. 3 As show in Table 6 and Fig. 12, it can be noted that the values of ETC ∗ and EPQ ∗ are little sensitive to Cy compared to Qin and . The reason behind this is when Cy increases, variables of the model are modified in a way that ETC∗ and EPQ ∗ approximately remain unchanged while range of changes in ETC∗ and EPQ ∗ is lower than Qin and . Furthermore, the obtained results shows that Cy in contrast to the other parameters has non-linear effects on ETC ∗ and EPQ ∗ . According to Fig. 9, trends of n, h1 and k with

6. Conclusions To fill the gap between the simplified assumptions in traditional perfect production model and the real production conditions, this study proposed an integrated model for imperfect production process by incorporating three issues of statistical process monitoring, maintenance, and inventory. Moreover, the proposed method considered multiple assignable causes due to the usual complexity of manufacturing systems and variable sampling intervals in a way that integrated failure rate remains a constant value over all sampling intervals. PSO algorithm was applied to obtain optimal values of process variables that minimize the expected total cost per production cycle. The attained results indicated that the proposed method outperforms single assignable cause model with respect to both economic and statistical criteria. Eventually a sensitivity analysis was conducted on three parameters of: 1) shape parameter of time to shift distribution, 2) the false alarm cost, and 3) the quality loss cost per unit when process is in the in-control state. The results demonstrated that the shape parameter of time to shift distribution unlike the false alarm cost has significant effects on both ETC∗ and EPQ ∗ while the quality loss cost only influences ETC∗. To make further improvement, the process can be monitored by other control charts such as adaptive control chart, CUSUM or EWMA control chart to accelerate detecting small and median shifts. Furthermore, the presented model can be extended to simultaneously monitor mean and variance of quality characteristic under consideration.

Appendix A. Proof of definition (2) In the proposed model there are s types of assignable causes that are independent of each other and each one occurs randomly. Considering Ati as the time that process remains in the in-control state before ith assignable cause happens, the time before occurring earliest assignable cause is attained as follows:

56

A. Salmasnia et al. / Journal of Manufacturing Systems 42 (2017) 44–56

P(T > t) = P(At1 > t, At2 > t, ..., Ats > t) = P(At1 > t)P(At2 > t)...P(Ats > t)

s 

−( 





P(T > t) = e−1 t e−2 t ...e−s t = e 

= e−0 t ⇒ f0 (t) = 0 t −1 e−0 t

i )t 

i=1



Appendix B. Derivation of equation (3) Regarding Eq. (2), Wj and hj can be obtained as follows:

h1



Wj+1

ri (t)dt = Wj



Wj+1

i t −1 dt =

ri (t)dt ⇒ Wj

0

h1

 w

i t −1 dt ⇒ i t ␯

j+1

wj

 h1

= i t ␯

0

⇒ Wj+1  − Wj  = h1 

0 1

if : j = 1 ⇒ W2  = h1  + W1  ⇒ W2  = 2h1  ⇒ W2 = 2  h1 1

if : j = 2 ⇒ W3  = h1  + W2  ⇒ W3  = h1  + 2h1  ⇒ W3 = 3  h1 .. . 1

1

1

⇒ Wj = j  h1 ⇒ hj = Wj − Wj−1 ⇒ hj = j  h1 − (j − 1)  h1 1

1

⇒ hj = (j  − (j − 1)  )h1

References [1] Bashiri M, Amiri A, Doroudyan MH, Asgari A. Multi-objective genetic algorithm for economic statistical design of control chart. ScientiaIranica 2013;20:909–18. [2] Ben-Daya M, Makhdoum M. Integrated production and quality model under various preventive maintenance policies. J Oper Res Soc 1998;49:840–53. [3] Bisgaard S, Hunter WG, Pallesen L. Economic selection of quality manufactured product. Technometrics 1984;26:9–18. [4] Chen YS, Yang YM. Economic design of X-bar control charts with Weibull in-control times when there are multiple assignable causes. Int J Prod Econ 2002;77:17–23. [5] Cheng JC, Chou CY. A real-time inventory decision system using Western Electric run rules and ARMA control chart. Expert Syst Appl 2008;35:755–61. [6] Cheng L, Tsou CS, Yang DY. Cost-service trade-off analysis of reorder-point-lot-size inventory models. J Manuf Syst 2015;37:217–26. [7] Chih M, Yeh LL, Li FC. Particle swarm optimization for the economic and economic statistical designs of the control chart. Appl Soft Comput 2011;11:5053–67.

[8] De la Torre Gutierrez H, Pham DT. Estimation and generation of training patterns for control chart pattern recognition. Comput Ind Eng 2016;95:72–82. [9] Dos Santos Coelho L. An efficient particle swarm approach for mixed-integer programming in reliability–redundancy optimization applications. Reliab Eng Syst Saf 2009;94:830–7. [10] Faraz A, Saniga E. Economic statistical design of a T2 control chart with double warning lines. Qual Reliab Eng Int 2011;27:125–39. [11] Gan S, Zhang Z, Zhou Y, Shi J. Joint optimization of maintenance, buffer, and spare parts for a production system. Appl Math Modell 2015;39:6032–42. [12] Ghodratnama A, Jolai F, Tavakkoli-Moghaddam R. Solving a new multi-objective multi-route flexible flow line problem by multi-objective particle swarm optimization and NSGA-II. J Manuf Syst 2015;36:189–202. [13] Jiang Y, Chen M, Zhou D. Joint optimization of preventive maintenance and inventory policies for multi-unit systems subject to deteriorating spare part inventory. J Manuf Syst 2015;35:191–205. [14] Kennedy J, Eberhart RC, Shi Y. Swarm intelligence. San Francisco: Morgan Kaufmann Publishers; 2001. [15] Le MD, Tan CM. Optimal maintenance strategy of deteriorating system under imperfect maintenance and inspection using mixed inspection scheduling. Reliab Eng Syst Saf 2013;113:21–9. [16] Lee H, Cha JH. New stochastic models for preventive maintenance and maintenance optimization. Eur J Oper Res 2016;255:80–90. [17] Lee HL, Rosenblatt MJ. A production and maintenance planning model with restoration cost dependent on detection delay. IIE Trans 1989;21:368–75. [18] Lee PH, Torng CC, Liao LF. An economic design of combined double sampling and variable sampling interval control chart. Int J Prod Econ 2012;138:102–6. [19] Linderman K, McKone-Sweet KE, Anderson JC. An integrated systems approach to process control and maintenance. Eur J Oper Res 2005;164:324–40. [20] Liu L, Yu M, Ma Y, Tu Y. Economic and economic-statistical designs of an control chart for two-unit series systems with condition-based maintenance. Eur J Oper Res 2013;226:491–9. [21] Makis V, Fung J. An EMQ model with inspections and random machine failures. J Oper Res Soc 1998;49:66–76. [22] Mehmood R, Riaz M, Does RJ. Control charts for location based on different sampling schemes. J Appl Stat 2013;40:483–94. [23] Nenes G, Tasias KA, Celano G. A general model for the economic-statistical design of adaptive control charts for processes subject to multiple assignable causes. Int J Prod Res 2015;53:2146–64. [24] Nguyen DT, Dijoux Y, Fouladirad M. Analytical properties of an imperfect repair model and application in preventive maintenance scheduling. Eur J Oper Res 2016;256:439–53. [25] Pan E, Jin Y, Wang Sh, Cang T. An integrated EPQ model based on a control chart for an imperfect production process. Int J Prod Res 2012;50:6999–7011. [26] Perez RE, Behdinan K. Particle swarm approach for structural design optimization. Comput Struct 2007;85:1579–88. [27] Rahim AM, Ben-Daya M. A generalized economic model for joint determination of production run: inspection schedule and control chart design. Int J Prod Res 1998;36:277–89. [28] Rosenblatt MJ, Lee HL. Economic production cycles with imperfect production processes. IIE Trans 1986;18:48–55. [29] Talbi EG. Metaheuristics from design to implementation. John Wiley & Sons, Inc.; 2009. [30] Wen D, Ershun P, Ying W, Wenzhu L. An economic production quantity model for a deteriorating system integrated with predictive maintenance strategy. J Intell Manuf 2014:1–11, http://dx.doi.org/10.1007/s10845-014-0954-z. [31] Xiang Y. Joint optimization of control chart and preventive maintenance policies: a discrete-time Markov chain approach. Eur J Oper Res 2013;229:382–90. [32] Yin H, Zhang G, Zhu H, Deng Y, He F. An integrated model of statistical process control and maintenance based on the delayed monitoring. Reliab Eng Syst Saf 2015;133:323–33. [33] Zhou WH, Zhu GL. Economic design of integrated model of control chart and maintenance management. Math Comput Model 2008;47:1389–95. [34] Zhou X, Wu C, Li Y, Xi L. A preventive maintenance model for leased equipment subject to internal degradation and external shock damage. Reliab Eng Syst Saf 2016;154:1–7.