Microsoft Word - فهرست.doc DOI: 10.22094/JOIE.2018.242.1532 131 Three Hybrid Metaheuristic Algorithms for Stochastic Flexible Flow Shop Scheduling Problem with Preventive Maintenance and Budget Constraint Sadigh Raissi a,* , Ramtin Rooeinfar b ,Vahid Reza Ghezavati b a Islamic Azad University, South Tehran Branch, Tehran, Iran b School of Industrial Engineering, South Tehran Branch, Islamic Azad University, Tehran, Iran Received 20 2017; Revised 02 2018; Accepted 08 2018 Abstract Stochastic flexible flow shop scheduling problem (SFFSSP) is one the main focus of researchers due to the complexity arises from inherent uncertainties and also the difficulty of solving such NP-hard problems. Conventionally, in such problems each machine’s job process time may encounter uncertainty due to their relevant random behaviour. In order to examine such problems more realistically, fixed interval preventive maintenance (PM) and budget constraint are considered.PM activity is a crucial task to reduce the production efficiency. In the current research we focused on a scheduling problem which a job is processed at the upstream stage and all the downstream machines get busy or alternatively PM cost is significant, consequently the job waits inside the buffers and increases the associated holding cost. This paper proposes a new more realistic mathematical model which considers both the PM and holding cost of jobs inside the buffers in the stochastic flexible flow shop scheduling problem. The holding cost is controlled in the model via the budget constraint. In order to solve the proposedmodel, three hybrid metaheuristic algorithms are introduced. They include a couple of well-known metaheuristic algorithms which have efficient quality solutions in the literature. The two algorithms of them constructed byincorporationof the particle swarm optimization algorithm (PSO) and parallel simulated annealing (PSA) methods under different random generation policies. The third one enriched based on genetic algorithm (GA) with PSA. To evaluate the performance of the proposed algorithms, different numerical examples are presented. Computational experiments revealed that the proposed algorithms embedboth desirable accuracy and CPU time. Among them, the PSO-PSAП outperforms than other algorithms in terms of makespan and CPU time especially for large size problems. Keywords: Stochastic flexible flow shop; Budget constraint, Preventive maintenance; Genetic algorithm; Simulated annealing; Particle swarm optimization. 1. Introduction The flexible flow shop scheduling problem (FFSSP) consists of a flow manufacturing line with one or more parallel machine on some processing stages (or workstation) in series. Multiple products (or jobs) are produced in each stage. The objective function of FFSSP is on minimizing the total completions of all jobs (makespan or Cmax)(Hoogeveen et al., 1996; Gupta, 1998). This kind of issue often takes into account an NP-hard problem. This means that at a reasonable computational time, an optimal answer cannot be obtained. (Brucker and Kramer, 1995)have proved that a two-stage FFSSP remains NP-hard even there is only one machine on the first stage and two machines on the second stage.Most researchers have focused on the system which their relevant process time deploy a given deterministic value. However, in real circumstances, processing time of jobs at each stage is the main part of makespan. Most research focused on the system with deterministicprocessing time. However, in a real system, theprocessing time is arandom variable due to random behaviour of tool wearing, operator skill, material variability and so on (Koulamas and Kyparisis, 2000). According to Choi and Wang (2012) makespan estimationmight become invalid under differentcircumstances. Another important factor that effects on makespan is interruption caused by machine breakdowns (Fahmy and Sharif, 2009). One of the most important ways to deal with these types of random interruptions is on following preventive and scheduled tasks. Hence, unavailability of machines, due to preventive maintenance, could be incorporated as a set of constraints in the mathematical models. Taking into account such disruptions due to preventive maintenance and stochastic processing time makes the problem more real, but more complicated and has been kept on the research gap at yet.The significant contributions of this paper are to propose a more realistic scheduling model for such production and delivering an efficient solution methods.The proposed model integrates FFSSP with preventive maintenance (PM) and budget constraint under stochastic processing time. Hence,consideringqueues between consecutive stages is a respectable alternative. When a job is processed at upstream stage and all the downstream machines get to busy state or comes under *Corresponding author Email address: raissi@azad.ac.ir , Summer & Autumn 2019,Vol.12, Issue 2 Journal of Optimization in Industrial Engineering 131-147 November September September Sadigh Raissi et al. / Three Hybrid Metaheuristic Algorithms… 132 any PM activity, therefore the job should wait in the bufferswhichcause holding cost.Because the total budget in handhas a given level this leads to interrupt the operations and cause increasing of the makespan. Consequently, in the proposed model the jobs are scheduled in such a way that there is no interruption in their operations. By considering all of these affecting factors, the optimized makespan will be acquired by the model. Obviously, such FFSSP also ruins an NP-hard model. Therefore, three hybrid populations based metaheuristic algorithms are introduced as solution methods. The flexible flow shop scheduling problem has been studied extensively in the literature. Koulamas and Kyparisis(2000) developed a heuristic method for a two and a three-stage FFSSP with random processing time based on the makespan objective function. They proved the effect of the proposed heuristic by finding the solutions which were better than the available lower bounds. A Tabu search (TS) algorithm together with a procedure for a constructing a complete schedule to solve the FFSSP with limited buffers has been given by Wardono and Fathi(2004) that minimizes job completion time. Their algorithm acts based on the stage-oriented decomposition approach. Akrami et al. (2006) presented a mixed-integer linear mathematical programming for FFSSP. They assumed that there are limited buffers among stages. Two meta-heuristics, based on genetic algorithm (GA) and TS algorithm, were presented to optimality solve the model. A flow shop scheduling was investigated with limited intermediate buffers was studied by Wang and Tang (2009). They focused on the makespan minimizing as a performance criterion and they applied a TS algorithm to optimize the problem. In order to improve the diversity of the TS, a scatter search mechanism was applied. The computational results showed the high efficiency of their proposed hybrid metaheuristic. Al-Hinai and ElMekkawy(2011) proposed a flexible flow shop problem with random machine failures and a two stage hybrid GA was applied. The first stage was on minimizing the primary objective; the makespan, and the second stage was optimized the bi- objective function and integrates machines assignments and operations sequencing with the expected machine breakdown. Tran and Ng (2011)presented a water-flow algorithm to solve the FFSSP with intermediate buffers. They pooled the amount of precipitation and falling force to form flexible erosion capability. This work helped the erosion process of the algorithm to focus on exploiting promising regions strongly. They also utilized an improved procedure for constructing a complete schedule from a permutation that represents the sequence of jobs at the first stage of the scheduling problem. Kianfar et al. (2012) investigated a FFSSP with non-deterministic arrival of jobs and sequence dependent setup times. They used average tardiness of jobs as the objective function. To optimize the problem, they presented a novel dispatching rule and hybrid GA. A computer simulation model was also developed to evaluate the presented dispatching rule. The results showed that their proposed dispatching rule can lead to much better results in comparison with the traditional dispatching rules. Singh and Mahapatra (2012) proposed a novel PSO to solve FFSSP. An efficient mutation operator was embedded in PSO to prevent solutions from falling into the local optimums. The performance of the PSO was evaluated against GA by a set of test problems taken from the literature. According to the obtained results, the percentage deviation of the proposed PSO of the lower bound is equal to 2.961 and the same measure for GA is equal to3.559.In order to deal with uncertain job processing times in FFSSP’s, Choi and Wang (2012) presented a decomposition-based approach. The method combines two reactive approaches with a reactive- proactive approach. Lin and Ying (2013) proposed a hybrid algorithm based on the features of artificial immune systems and the annealing process of simulated annealing algorithms (SA) to optimize the FFSSP with limited buffer storage between stages. Almeder and Hartl(2013) studied the FFSSP with limited buffer storage in the metal–working industry. The partitioned the problem as a two stage problem. The first stage contains a single machine and its buffer. The semi-finished parts are stored in this buffer until a machine of the second stage is available. The second stage contains two parallel non- identical machines. They applied a hybrid approach based on discrete-event simulation and variable neighbourhood search to optimize the problem. The hybrid approach was led to an improvement between 3% and 10% compared with the current production plan of the company. By considering unrelated parallel machines, sequence- dependent setup times, probable reworks and different ready times, Rabiee et al. (2014) investigated the no-wait two-stage FFSSP. They proposed a novel hybrid algorithm based on imperialist competitive algorithm (ICA), SA, variable neighbourhood search (VNS) and GA. The performance of the proposed hybrid algorithm was evaluated against ICA, SA, VNS, GA and ant colony optimization (ACO) and high efficiency of their proposed hybrid algorithm was revealed. Arnaout (2014) tackled a rescheduling problem for the flow shop problem associated with stochastic processing and setup time. They proposed a new repair rule which and compared it with the existing algorithms. The results obtained from the computational experiments showed that the proposed repair rule can perform better those available algorithms in literature. Rahmani and Heydari (2014) studied FFSSP under uncertain processing times and unexpected arrivals of new jobs. They proposed a new approach to find robust schedules in this situation. Their approach was a proactive–reactive method which uses a two-step procedure. In order to consider stochastic processing time, Wang and Choi (2014) introduced a novel decomposition- based the Holonic approach (DBHA) to solve a FFSSP with stochastic processing time. Their proposed method was based on genetic algorithm control (GAC) and the shorting processing time contract net procontrol (SPT- CNP). Also, K-means clustering is utilized to divide machines into different clusters according to their Journal of Optimization in Industrial Engineering Vol.12, Issue 2, Summer & Autumn 2019, 131- 147 133 stochastic environments. The gained results showed that DBHA had better quality solutions against GAC and SPT- CNPA simulation based optimization approach for stochastic hybrid FFSSP in a real-world semiconductor back-end assembly facility was presented by Lin and Chen (2015). In their approach, the optimization strategy, based on genetic algorithm, was used to find the optimal assignment of the production line and machine type at each stage. The simulation model was used to evaluate the performance of solutions. They applied a real case study to prove the necessity of using simulation optimization approaches for practical applications. A novel algorithm was developed by Li and Pan (2015) to solve the hybrid flow shop with limited buffers. They combine two metaheuristic algorithms of artificial bee colony and TS and used makespan as the performance criterion. Their proposed algorithm was evaluated against several algorithms reported in the literature and the experimental results showed the high effectively and efficient performance of the proposed algorithm. Sangsawang et al. (2015) proposed two metaheuristic algorithms to optimize a two-stage re-entrant FFSSP. The first algorithm was a hybridization of GA and adaptive auto-tuning based on a fuzzy logic controller. The second one was a hybridization of particle swarm optimization (PSO) and Cauchy distribution. Experimental results revealed that both hybrid algorithms could present better solutions and are more powerful than the classical metaheuristics. The statistical analyses indicated that the hybrid PSO and hybrid GA can improve the best solutions in the literature by averages of 15.60% and 15.51%, respectively. Zabihzadeh and Rezaeian(2015) developed a mixed integer linear programming model for FFSSP. They assumed that there are some robots between stages for unloading, transferring and loading parts. Their proposed model has ability of determining the number of robots and jobs sequence. They applied two meta-heuristic algorithms; GA and ACO, to solve their model. Computational results showed that the GA is more efficient than ACO to optimize the model. Tang et al. (2016) presented a new model for dynamic FFSSP’s. By taking emergency maintenance, the proposed model minimized both objectives of energy consumption and makespan. Since they developed model was NP-hard, a multi-objective PSO algorithm to optimize the model. For finite capacity material requirement planning system in a flexible flow shop, Sukkerd and Wuttipornpun(2016) presented a hybrid GA and a Tabu search algorithm (HGATS). The results showed that the HGATS can outperform on comparing with the existing algorithms. Correspondingly, computational time of HGATS is acceptable when applied to real industrial systems. Rahmani and Ramezanian(2016) addressed a stochastic FFSSP in which new jobs arrive into the process as disruptions. They developed a mathematical model to minimize total weighted tardiness. A variable neighbourhood search algorithm was used to solve the model. The efficiency of the proposed algorithm evaluated through a set of test problems. González-Neira et al. (2016) presented a multi-criteria FFSSP in which criterion is quantitative and the other is qualitative. They assumed that job processing time deploys by a stochastic manner.The integral analysis method (IAM) was implemented to solve the relevant problem. In the IAM method, the problem first was introduced, then, cardinal analysis, ordinal analysis and integration analysis are done. Results showed that IAM method is able to select the alternatives with high efficiency in terms of both types of criteria. The most important criticism of scheduling problems is the gap between academic and practical problems. Even though the developed models have tried to be more realistic, but they fail to find the exact makespan in real life problems. This comes from the fact that the model cannot consider all the factors affected on makespan. As a result, there is a gap between obtained makespan by mathematical models and the makespan occurs in reality. In order to bridge this gap, this research presents an integrated mathematical model which is capable to consider all of the influential factors on makespan. The proposed model can simultaneously consider preventive maintenance, stochastic processing time and budget constraint. By integrating all of these subjects the model will reflect the real performance of a FFSSP. The rest of the paper is organized as follows: In section 2 the problem definition, mathematical model formulation and assumption of the model is presented. In section 3, the proposed hybrid metaheuristic algorithms in which three hybrid metaheuristic algorithms, including combining parallel SA with two types of PSO algorithms (PSO-PSAІ and PSO-PSAП), and a GA with parallel SA (GA-PSA) are specially explained. In section 4, computational results are presentedwhich actually compare the results of proposed metaheuristic algorithms. First, the small scale problem size is solved with CPLEX software. Then, the metaheuristic algorithms are used to find the near optimal solutions for the large scale of the problems. Analysis of the results and comparisons demonstrate the performance of the proposed solving methodologies on different problem sizes. Finally, conclusions and future researches are presented in the last section. 2. Mathematical Model Formulation In this section the problem under consideration is described. Consider a problem of J jobs and S working stages as shown in Figure 1. Each stage s has a number of 𝑁𝑠 identical parallel machines that operate in parallel. All jobs visit all stages from the first to the last stage. They are processed by one machine at each stage. Each job processed all the stages and every machine process maximum one job at a time. Thereis a buffer between two consecutive stages. The machines are under PM tasks. . Consequently, each machine could be available or unavailable at each time. All jobs are independent of each other and they are available at time zero. The processing time of jobs is considered to be stochastic. In order to Sadigh Raissi et al. / Three Hybrid Metaheuristic Algorithms… 134 model such randomly distributed processing time, the stochastic programing technique is applied in which each possible processing time is called as a scenario. The probability of occurrences of each scenario is shown by Pr(π). Therefore, the averaged processing time of each job at each stage is calculated according to its processing times in different scenarios and the probability of occurrences of each scenario. The objective is to set a sequence of jobs, allocating jobs to machines at each stage and determining the time between two consecutive maintenance activities in such a way that total completion time being minimized. Figure 2 shows an example of the Gantt chart of proposed stochastic flexible flow shop scheduling problem (SFFSSP) with 3 stages and 5 jobs. There is2, 3 and 2 machines in stages 1, 2 and 3, respectively. In machine 1 in stage 1, maintenance occurs after processing job 3. The other maintenance occurs in stages 2 and 3. The job 2 is processed in stage 1 from time 2 to 5 and is started in stage 2 from time 7. Therefore, job 2 is sent to the buffer of stage 1 and cause holding cost for each period of time, until machine 1 in stage 2 gets empty. Sometimes more than one job is placed in the buffers. For example, in the buffer of stage 2, from time 11 and 12, jobs 1 and 3 are placed in and increasing the holding cost for staying of each period of time in the buffer. Other constraints as well assumptions are listed as follows:  There is a set of jobs denoted J (j=1, 2,…, J) jobs which are available at time zero and no job may be cancelled before completion.  The set of L consecutive stages denoted by (s=1, 2,…, L),  Each stage is equipped with non-identical parallel machines denoted by (m=1, 2,…, 𝑁𝑠 ).  The set of R, denoted by (r=1,2,…, R) indicated the number of intermediate buffers.  All jobs have to process serially through all machines at each stage for only once time.  The interruption processing of each job on all machines at each stage is not acceptable. On the other hand, once a job is started, it must be processed to completion without any interruption either on or between machines.  All machines are continuously available at time zero, in another word; non-machines are failing at the starting time.  Each machine could process only one job at the same time, and each job has to visit each machine exactly once.  The preventive maintenance activities are performed on each machine at the fixed intervals (𝑃𝑚𝑠).  Once the preventive maintenance activity is carried out, there is no probability of a subsequent machine breakdown.  All jobsat each intermediate buffer havea same holding cost, but different at each stage. Fig.1.The framework of the flexible flow shop problem of this study 2.1. Notations To present the model using mathematical terms, consider the following notations. Indices s Indexofstages{s= 1,2,…, L} m Index of machines at s {m= 1, 2, …,𝑁𝑠} j Index of jobs{d, j= 1, 2,…, J} r Index of intermediate buffers{r= 1,2,…, R} h Index of job sequence{h, u=1 ,2,…,𝐾𝑚} n Index of maintenance activity which is done on machine j {n= 1,2,…,𝑉𝑚𝑠} 𝜋Index of probabilistic scenarios{𝜋 = 1, 2,…,𝛱} Parameters 𝑝𝑗𝑠(𝜋)The processing time of job j at stage s in scenario π 𝑃𝑟⁡(𝜋)The probability of occurrences scenario π ℎ𝑟𝑠𝑗 Holding cost of job jinintermediate buffer ofr at stage s 𝜇𝑚𝑠The repair rate of machine m at stage s 𝜆𝑚𝑠The failure rate of machine m at stage s 𝐷𝑚𝑠The duration time of maintenance activity of machinematstages �̅�𝑚𝑠The mean number of jobs that are performed on machine m at stage s 𝛽The minimum of availability of the system M An arbitrary large position number B Total budget Dependent decision variables 𝑆𝑗𝑠Starting time for the processing of job j at stage s 𝐶𝑗𝑠Completion time of job j at stage s 𝑃𝑚𝑠The time between two consecutive maintenance activities on machine m at stage s 𝑚𝑚𝑠The number of maintenance activities on machine m at stage s 𝑑𝑟𝑠𝑗 Waiting time of job j atintermediate bufferr at stage s 𝑍𝑗𝑟(𝑡) {1 ⁡if⁡job⁡𝑗⁡is⁡in⁡intermediate⁡buffer⁡𝑟⁡at⁡time⁡𝑡0 otherwise 𝑇𝑚𝑠The completion time of the last PM action on machine m at stage s 𝐴𝑚𝑠(𝑡)The availability of machine m at stage s at timet 𝐴𝑠(𝑡)The unavailability of stage s at time t 𝐴𝑠𝑦𝑠(𝑡)The unavailability of system at time t W 0 or 1 Journal of Optimization in Industrial Engineering Vol.12, Issue 2, Summer & Autumn 2019, 131- 147 135 Independent decision variables 𝑋𝑗ℎ𝑚𝑠 {1 ⁡if⁡job⁡𝑗⁡is⁡processed⁡in⁡sequence⁡ℎ⁡by⁡machine⁡𝑚⁡at⁡stage⁡𝑠0 ⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡otherwise 𝑌𝑛𝑚𝑠 {1 ⁡if⁡𝑛⁡th⁡maintenance⁡activity⁡on⁡machine⁡𝑚⁡at⁡stage⁡𝑠⁡is⁡executed0 otherwise 2.2 Mathematical model In our proposed model,availability of machine min stage s at time t under PM activities could be calculated according to Villemeur(1991) by Eq. 1. 𝐴𝑚𝑠(𝑡) = 𝜇𝑚𝑠𝜇𝑚𝑠+𝜆𝑚𝑠 + 𝜆𝑚𝑠𝜇𝑚𝑠+𝜆𝑚𝑠 𝑒𝑥𝑝[−(𝜇𝑚𝑠 + 𝜆𝑚𝑠)(𝑡 − 𝑇𝑚𝑠)]⁡⁡⁡⁡⁡⁡(1)re we considered system configuration as a parallel-series system in which machines ateach stage are parallel and stages are a series. A stage is said to be unavailable if all of its machines are unavailable. Therefore, the unavailability of stage s is calculated by Eq. 2 and unavailability of the total system by Eq. 3. 𝐴𝑠(𝑡) = ∏(1 − 𝐴𝑚𝑠(𝑡))𝑁𝑠𝑚=1 ⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡(2) 𝐴𝑠𝑦𝑠(𝑡) = 1 − ∏(1 − 𝐴𝑠(𝑡))⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡(3)𝐿𝑠=1 Therefore, the proposed model as an extension to basic model was taken fromGonzález-Neira et al. (2016) will be as follows. 𝑀𝑖𝑛⁡⁡{𝑚𝑎𝑥⁡(𝐶𝑗𝑠)}⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡(4) St: ∑ ∑ 𝑋𝑗ℎ𝑚𝑠𝐾𝑚ℎ=1 𝑁𝑠 𝑚=1 = 1⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡∀𝑗 ∈ 𝐽;𝑠 ∈ 𝐿⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡(5) ∑ ∑𝑋𝑗ℎ𝑚𝑠𝐽𝑗=1 𝑁𝑠 𝑚=1 = 1⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡∀𝑠 ∈ 𝐿;ℎ ∈ 𝐾𝑚⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡(6) 𝐶𝑗𝑠 ≥ 𝑆𝑗𝑠 + ∑𝑃𝑟⁡(𝜋)𝑝𝑗𝑠(𝜋)𝛱𝜋 ⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡(7) 𝐶𝑗𝑠−1 − 𝑀 (1 − 𝑍𝑗𝑟(𝑡)) ≤ 𝑡 ≤ 𝑆𝑗𝑠 + 𝑀 (1 − 𝑍𝑗𝑟(𝑡)) ∀𝑗 ∈ 𝐽;𝑠 ∈ 𝐿;∀𝑟 ∈ 𝑅|𝑠 = 𝑟⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡(8) 𝑆𝑗𝑠 + (1 − 𝑋𝑗ℎ𝑚𝑠)𝑀 ≥ 𝐶𝑑𝑠 − (1 − 𝑋𝑑𝑢𝑚𝑠)𝑀⁡⁡⁡⁡ ∀𝑗,𝑑 ∈ 𝐽⎹⁡𝑗 ≠ 𝑑; ⁡ℎ,𝑢 ∈ 𝐾𝑚⎹⁡𝑢 < ℎ,𝑚 ∈ 𝑁𝑠;𝑠 ∈ 𝐿⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡(9) ∑∑ 𝑑𝑟𝑠𝑗 ℎ𝑟𝑠𝑗𝑅𝑟=1 𝐽 𝑗=1 ≤ 𝐵⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡(10) 𝑑𝑟𝑠𝑗 ≥ 𝑆𝑗𝑠+1 − 𝐶𝑗𝑠⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡(11) (𝑛𝑃𝑚𝑠𝑌𝑛𝑚𝑠 − 𝐶𝑗𝑠)𝑋𝑗ℎ𝑚𝑠 ≥ −𝑀⁡(1 − 𝑊)⁡ ∀𝑗 ∈ 𝐽;∀𝑚 ∈ 𝑁𝑠;∀𝑠 ∈ 𝐿; ℎ ∈ 𝐾𝑚;0 ≤ 𝑛 ≤ 𝑉𝑚𝑠⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡(12) (𝐶𝑗𝑠−𝑝𝑗𝑠(𝜋) − 𝑛𝑃𝑚𝑠𝑌𝑛𝑚𝑠 − 𝐷𝑚𝑠)⁡𝑋𝑗ℎ𝑚𝑠 ≥ −𝑀(𝑊)⁡⁡⁡⁡⁡⁡⁡ ⁡⁡⁡⁡0 ≤ 𝑛 ≤ 𝑉𝑚𝑠⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡(13) 𝑃𝑚𝑠 = ∑ ∑ 𝑋𝑗ℎ𝑚𝑠 ∑ Pr(𝜋) 𝑝𝑗𝑠(𝜋)𝛱𝜋𝐾𝑚ℎ=1𝐽𝑗=1 𝑚𝑚𝑠 ∀𝑚 ∈ 𝑁𝑠;∀𝑠 ∈ 𝐿⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡(14) 𝑇𝑚𝑠 ≤ 𝑛𝑃𝑚𝑠⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡(15) 1 − 𝐴𝑠𝑦𝑠(𝑡) ≥ 𝛽⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡(16) 𝐶𝑗𝑠,𝑑𝑟𝑠𝑗 ≥ 0⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡∀𝑗 ∈ 𝐽; ∀𝑠 ∈ 𝐿⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡(17)⁡ 𝑋𝑗ℎ𝑚𝑠 ∈ {0,1}∀𝑗 ∈ 𝐽;∀𝑚 ∈ 𝑁𝑠;∀𝑠 ∈ 𝐿; ℎ ∈ 𝐾𝑚⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡(18) 𝑌𝑛𝑚𝑠 ∈ {0,1}∀𝑚 ∈ 𝑁𝑠;∀𝑠 ∈ 𝐿;𝑛 ∈ 𝑉𝑚𝑠⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡(19) Eq. (4) indicates the objective function. The constraint (5) ensures that the assignment of each job to one and only one machine at each stage. Constraints set (6) determines that the position of a machine sequence is takes by only one job at each stage. Constraints set (7) states that it is not allowed starting processing jobs at the next stage unless they have completed processing at the previous stage. Constraint sets (8) and (9) indicate that no interference should be taken among jobs on a common machine at any stage if the machine is available. On the other words, the difference between the processing times of any two jobs assigned to the same machine should not have any overlap.The constraint set (10) guarantees that the holding costs should be less than the available budget. Constraints set (11) determines the waiting time of jobs in each buffer.Constraints set (12) and (13) ensure that there are no overlap among operations and maintenance task. Constraint sets (14)-(16) are related to control the system availability.Finally, constraint set (17)-(19) controls the decision variable types. Sadigh Raissi et al. / Three Hybrid Metaheuristic Algorithms… 136 Fig. 2. An example of proposed SFFSSP Gantt chart of this study 3. Solution Methodologies As mentioned earlier, the proposed modelof SFFSSP is an NP-hard problem and solving this problem with exact methods in reasonable CPU time is impossible, so we should use heuristic or metaheuristic algorithms to solve it(Pinedo& Chao, 1999). Pinedo (2002) proved that some exact methods have been developed for solvingSFFSSP’s, but they are not suitable for more than 20 jobs and 10 machines. Also, heuristic algorithms are often may be trapped onsome local solutions, but metaheuristic algorithms can be described as a master strategy that guides and modifies subordinate heuristics to explore the solutions beyond the local optimally (Osam et al., 1996).We describe threehybridmetaheuristic algorithms, including combining parallel SA with two types of PSO algorithms (named PSO-PSAІ and PSO-PSAП), and a GA with parallel SA (named GA-PSA) which have good solution quality in the literature(Singh and Mahapatra, 2012; Kianfar et al., 2012; Rabiee et al., 2014; Tang et al., 2016; Sukkerd and Wuttipornpun, 2016). Theproposed metaheuristic algorithms are explained in the nextsubsections on details. 3.1. Hybrid PSO-PSA algorithms PSO algorithm is an evolutional solution method performed on a population of candidate solutions called particles. These particles move around in the search-space according to simple routine. The particles movements are guided by the best found positions in the search-space, which are continually updated as better positions are found by the particles. At each iteration, the in position of a particle (X vector) is updated by calculating the velocity (Vel vector) using the differences between the current position of the particle and the two following vectors (Kennedy and Eberhart, 1995):  The best position experienced by the particle in all previous iterations. This is called the particle best (Pbest).  The best position experienced by all particles in population in all previous iterations. This is called the global best (gbest). Generally, a PSO is a continuous algorithm inherently, while SA is a discrete one. Experiments show that combining PSO with a discrete algorithm such as SA creates better performances. Poli, Kennedy and Blackwell (2007) presented a review on the variation and the hybridization of the PSO. We proposed two types of hybrid PSO-PSA algorithms named PSO-PSAІ and PSO- PSAП to have both advantages of these methods. The basic idea of the hybrid PSO-PSA algorithm is running the PSO algorithm and improving the best results by applying parallel SA (PSA). In order to have variety in the proposed algorithm, we consider every bad, normal and good solutions could be selected with same chances.Combining PSO and PSA decreases the probability to be trapped in the local optimal solutions. Also, by introducing a suitable neighbourhood formation structure, the search process is enhanced and finds the near global optimum solution. The essential components of our proposed PSO-PSAІ and PSO-PSAП are similar, except in generating initial solutions as described below. The pseudo-code for our PSO-PSA algorithm is shown Figure 6. 3.1.1. Initial solution of PSO-PSAІ The random generation policy is used to generate the initial population of PSO-PSAІ. 3.1.2. Initial solution of PSO-PSAП In PSO-PSAП a new procedure is applied to generate initial solutions. First, we consider SFFSSP with relaxed conditions in which the binary constraints of model are Journal of Optimization in Industrial Engineering Vol.12, Issue 2, Summer & Autumn 2019, 131- 147 137 relaxed. In other word the 𝑋𝑗ℎ𝑚𝑠is considered to be a real number between 0 and 1. Next, the relaxed model was solved by CPLEX. The optimum values of decision variables are numbers between 0 and 1 which are looked as a probability distribution. These values are used as primal inputs of initial solutions for PSO-PSAП. For example, if is X2312 = 0.8 then we set X2312 = 1 with the 0.8 and X2312 = 0 with the 0.2. The initial solutions obtained by this methodare distributed in the high quality of solution space. In order to perform a comprehensive search the low quality of solution space must be investigated. Therefore, some initial solutions are generated, unlike the obtained probability distributions. For example, if X2312 = 0.8, then we set X2312 = 1 with the 0.2 and X2312 = 0 with the 0.8, the two groups if initial solutions are merged and formed initial solutions of PSO- PSAП. The input parameters of both proposed PSO-PSAІ and PSO-PSAП are: the population size (Npop), the number of successive iterations in which the best solution does not change (M), the cognition learning factor (C1), The social factor (C2) the inertial weight (w), the population size (Npop), the number of internal loop (In loop), and the temperature decreasing rate (α). The other components of proposed PSO-PSAІ and PSO-PSAП are the same and described as follows: 3.1.3. The solution structure of PSO-PSA The solution is represented by two sub-matrixes, each of them are associatedto a special area of decision variable. The first sub-matrix presents the sequence of jobs contains S*J matrix where S is the number of stations and J is the number of jobs. An enhanced version of random key representation, proposed by Norman and Bean (1999) is used to show the sequence of jobs which is capable to preserve the solution feasibility. In this way, each job at each station is assigned a real number between(1,𝑁𝑠). The integer part is the number determines the machine number to which the job is assigned and the fractional part is used to sort the jobs assigned to that machine. For example, consider a problem with 4 jobs, 5 stations and 4 machines at each station. An example of a solution for this problem is presented in Figure 3. According to Figure 4, in station 2, the jobs 1 and 4 are both processed on machine 3, and the job 2 is processed on machine 1, and the job 3 is processed on machine 2. Also the order of jobs to be scheduled on machine 3 is job 1 followed by job 4.The second sub-matrix is related to the number of maintenance activity on machines. It is a S*M matrix where S is the number of stations and M is the number of machines at each station. Each cell of this matrix is filled with a random number between 1 and the maximum number of allowable maintenance activity. Figure 4 is an example of a problem with 5 stations and 4 machines at each station. As shown in Figure 4, three, one, three and four maintenance activities are performed on machines 1, 2, 3 and 4 in station 1, respectively. Therefore, the decision variables 𝑌111 , 𝑌211 and 𝑌311will be equal to one. Fig.3. An example of representation of solution (Sequence of jobs) Fig. 4. An example of representation of solution (Maintenance activity) 3.1.4. Particle movements For particle movements, the following formulas are used to update the velocity and position vectors of a particle: 𝑉𝑒𝑙𝑖(𝑘 + 1) = 𝑤 ∗ 𝑉𝑒𝑙𝑖(𝑘) + 𝐶1 ∗ 𝑟1 ∗ (𝑃𝑏𝑒𝑠𝑡𝑖 − 𝑋𝑖(𝑘)) +𝐶2 ∗ 𝑟2 ∗ (𝑔𝑏𝑒𝑠𝑡 − 𝑋𝑖(𝑘)) (19) 𝑋𝑖(𝑘 + 1) =⁡𝑋𝑖(𝑘) +⁡𝑉𝑒𝑙𝑖(𝑘 + 1) (20) In equation (19), Veli(k) is the velocity of particle i in the k th iteration and Xi(k) states the position of particle i in iteration k. Also, Pbesti is the vector for the best known position of particle i and gbest is the best position vector of all particles in the whole population. 𝑤 is called the inertia weight that determines the impact of the current velocity of a particle on its velocity at the next iteration. The parameters C1 and C2 are acceleration coefficients which have constant values to determine the impact of Pbest and gbestvalues in defining the velocity, respectively. r1 and r2 are two random numbers uniformly distributed in [0,1]. 3.1.5. Local search improvements One of the challenges of the algorithms is trapped in local optimal solutions. In the other words, it is possible that the near optimum solution which has found yet, is selected and the algorithm accepts this solution as a final optimum solution and stops. The hybrid algorithms are used to combine the base algorithm with different strategies to improve the algorithm performances. We utilized medium radius local search for the proposed hybrid PSO-PSA algorithms. In this method, the local searches are not used on each swarm, and we used the double change technique. If the improvements on the convergence of fitness function are achieved, we replace it to the previous swarm, but, if no improvements have been seen, we accept this by Bultzen probability. Sadigh Raissi et al. / Three Hybrid Metaheuristic Algorithms… 138 3.1.6. Initial temperature A suitable initial temperature is the one that results in an average increase of acceptance probability near to one. The value of initial temperature will clearly depend on the scaling of fitness and, hence, it should be problem- specific. Therefore, we first generate a large set of random solutions, then a standard division of them is calculated and is used to determine the initial temperature in the way that the acceptance probability of primary generations reach to 0.999. Consequently, the initial Tk is set to 1000based on some preliminary examinations. 3.1.7. Cooling Schedule The performance of this algorithm also depends on the cooling schedule, which is essentially the temperature updating function. In the proportional decrement scheme, temperatures at the k and k+1 steps of the outer loop, 𝑇𝑘and 𝑇𝑘+1, are related by: 𝑇𝑘+1 = 𝛼𝑇𝑘⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡(21) Where 𝛼 is cooling rate and is obtained by some experiment. 3.1.8. Stopping criteria To limit the number of iterations of PSO-PSA algorithms, some convergence experiment was performed and the best criterion was applied as follows: PSO-PSA will be stopped when the best solution does not change after a pre-determined number of successive iterations (M). Also, the PSA is allowed to search in a temperature level, for In-loop iterations. The optimum value of M and In-loop is determined by Taguchi experiments. The pseudo-code of the proposed PSO-PSA is described in Figure 5. Procedure of hybrid PSO-PSA algorithms P: initial particle with random positions and velocities Untiltermination condition is met Do For each particle ido Update the velocity of particle i Update the position of particle i Evaluate particle i Update Pbest and gbest Endfor Start PSA with some of the best and the worst chromosomes Set repetition counter k = 0 Untiltermination condition is met Do Set repetition counter M = 0 UntilM= in_loopDo Generate a neighbour solution: ω2 Calculate Δω1, ω2= f(ω) - f(ω0) If Δω1, ω2 ≤ 0, then ω1 = ω2 If Δω1, ω2 >0, then ω1 = ω2 with probability exp (-Δω1, ω 2 / tk) m = m + 1 End Tk+1= α Tk k = k + 1 End Transfer the improved solution to the particles End Fig. 5. Pseudo-code of hybrid PSO-PSA 3.2. Hybrid GA-PSA algorithm (GA-PSA) Genetic algorithm (GA) has no ability to search effectively to find the best global optimum solution. Also, this algorithm isn’t a capable to complete local searches on solutions. Therefore, we can combine the power of GA in global search with simulated annealing (SA) local searches to address the global optimum solution. This hybrid algorithm which combines GA and SA has both advantages of these algorithms and help to improve solution performances.In this hybrid algorithm, first the GA generates initial solutions with crossover and mutation operators. Next, some of these solutions have been selected for the parallel SA (PSA) as initial solutions. Then, the parallel local search process on the selected solution starts. To have variety in the proposed algorithm, we consider every bad, normal and good solutions could be selected with same chances. The PSA procedure in the same applied in PSO-PSA. The other components of PSO-SA are as follows: 3.2.1. Initialization The input of GA-PSA is the population size (Npop), the number of successive iterations in which the best solution does not change (M), the crossover probability (Pc), and the similarity coefficient (SC), and the mutation probability (Pm) the number of internal loops (In loop), and the temperature decreasing rate (α) are first initialized. Then, to generate an initial population, a random generation policy is utilized in this step. Since the solutions obtained by a metaheuristic algorithm are sensitive to their parameter values, a statistical procedure based on the Taguchi parameter tuning method is used to tune the parameters. 3.2.2. Selection operator One of the most key elements of a GA is the selection operator which is used to select chromosomes (parents) which lead to generate new chromosomes (offspring). The proposed selection operator is roulette wheel selection method in which parent chromosomes are probabilistically selected based on their fitness function value. The better chromosomes are selected with the highest probability. Using the roulette wheel selection each chromosome in the population occupies a slot with a slot size proportional to the chromosome fitness. When the wheel is randomly spun, the chromosome corresponding to the slot where the wheel stopped is selected as the first parent. This process is repeated to find the second parent. Clearly, since better chromosomes have larger slots, they have better chances to be chosen in the selection process. Journal of Optimization in Industrial Engineering Vol.12, Issue 2, Summer & Autumn 2019, 131- 147 139 3.2.3. The chromosome structure The chromosome structure of the GA-PSA is just like that used in “solution structure” in PSO-PSA. 3.2.4. Chromosome evaluation In order to evaluate chromosomes, each chromosome is simulated for 30 times. Then, the obtained results are averaged and considered as the chromosome fitness function. In the other word, by changing the stochastic input parameters, the fitness function of the chromosome will change. Therefore, the fitness function of each chromosome is not under deterministic value and show the stochastic nature of the model. 3.2.5. The crossover operator As defined in the previous section, the chromosomes have two parts, so a crossover operator is applied in these two parts. For each part of the chromosome a single-point crossover is applied. For the first part of the chromosome, a cross point between (1, J) is generated (J is the number of jobs). Then, the crossover operator is applied according to Figure 8. The same procedure is applied to generate the second of the offspring chromosomes. 3.2.6. The mutation operator The mutation operator is used in only some iterations. The similarity coefficient (SC) determines if mutation is applied or not. We can calculate the SC as follows: 𝑆𝐶𝑎𝑏 = ∑ ∑ ∑ ∑ 𝜕(𝑋𝑗ℎ𝑚𝑠(𝑎), 𝑋𝑗ℎ𝑚𝑠(𝑏))𝐿𝑠=1𝑁𝑠𝑚=1𝐾𝑚ℎ=1𝐽𝑗=1 𝑆 × 𝐽 (22) Where 𝑋𝑗ℎ𝑚𝑠(a)⁡and⁡𝑋𝑗ℎ𝑚𝑠(𝑏) are decision variables in chromosomes a and b. For comparing the similarity between two chromosomes, we consider the similarity of each gene that can be expressed as follows: 𝜕(𝑋𝑗ℎ𝑚𝑠(𝑎), 𝑋𝑗ℎ𝑚𝑠(𝑏)) = {1⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡𝑖𝑓⁡𝑋𝑗ℎ𝑚𝑠(𝑎) = 𝑋𝑗ℎ𝑚𝑠(𝑏)⁡⁡0⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒⁡ (23) The average similarity coefficient of the population is calculated as follows: 𝑆𝐶̅̅̅̅ = ∑ ∑ 𝑆𝐶𝑎𝑏𝑁𝑏=𝑎+1𝑁−1𝑎=1 (𝑁2) (24) WhereN is the number of chromosomes in the population. Considering a pre-defined threshold similarity coefficient, the specified mutation operator will be automatically applied. Two swapping types are used in proposed GA- PSA. Swapping type 1 is used to define the neighbourhood N(s) in local search (PSA) and swapping type 2 is used as mutation operator of GA.  Swapping type 1 The mutation operator is applied on both two parts of chromosomes. To this aim, a column of each chromosome sub matrix is randomly selected and is inversely arranged. Figure 6 shows an example of the mutation operator. Fig.6. An example of mutation operator of swapping type 1 (Sequence of jobs)  Swapping type 2 The swapping type 2 works which select two rowsof each chromosome sub matrix is randomly selected and swapped (see Figure7). Fig. 7. An example of mutation operator of swapping type 2 (Sequence of jobs) Sadigh Raissi et al. / Three Hybrid Metaheuristic Algorithms… 140 Fig. 8. An example of crossover operator (Sequence of jobs) 3.2.7. Stopping criteria The stopping criteria in the GA-PSA is similar to the ones described for PSO-PSA. Figure 9 shows the pseudo-code of the proposed GA-PSA algorithm. Procedure of hybrid GA-PSA algorithm P: generate the initial population (Npop) Untiltermination condition is metDo Evaluate the chromosomes Apply the crossover operator Apply the mutation operator Start PSA with some of the best and the worst chromosomes For all the selected chromosomes Do Set repetition counter k = 0 Untiltermination condition is met Do Set repetition counter M= 0 UntilM= in_loopDo Generate a neighbour solution: ω2 Calculate Δω1, ω2= f(ω) - f(ω0) If Δω1, ω2 ≤ 0, then ω1 = ω2 If Δω1, ω2 >0, then ω1 = ω2 with probability exp (-Δω1, ω 2 / tk) m = m + 1 End Tk+1= α Tk k = k + 1 End End Transfer the improved solutions to the genetic population and combine all the generated populations End Fig. 9. Pseudo-code of GA-PSA 4. Computational Evaluation All experimental results have been carried out on an ASUS laptop with a 2.4 GHz, core i5 processor using a 4GB of RAM. All metaheuristic algorithms have been implemented in MATLAB Software (Version 7.10.0.499, R2010a) and the linear programming models have been solved using CPLEX 12. Also, the Minitab 17 software has been used to apply theTaghuchi design of experiment method for parameters tuning.Two sets of test problems are applied to solve the model. These test problems are generated based on the data given in Table 1. Here, we considered 5 scenarios, each of which has equal probability to be occurring.Three machines workat each station, considering the number of jobs and the numbers of stations, the two test sets aregeneratedin small and large problem sizes. First, 18small sized test problems are generated and the solutions, obtained by the algorithms, are evaluated against global optimum amounts.Furthermore, 60 large sizedtest problems are applied to compare the performance of the algorithms. We considered two kinds of total budget called hereinafter by type 1 and 2. Both budget types are associated to the number of jobs. The budget types 1 and 2 are respectively calculating by multiplying of number of jobs in 1200 and 1000. Therefore, the budget type 2 is more strict than type 1. Table 1 Factor and their levels Factors Levels Number of jobs (j) [4, 5, 6, 20, 50, 60, 100] Number of stations (s) [2, 3, 4, 6, 9] Holding cost ($/sec) Uniform [50, 100] Process times (sec) Uniform [1, 99] The minimum of availability 0.95 The repair rate of machine (sec) Uniform [10, 50] The failure rate of machine (sec) Uniform [150, 450] Number of scenarios 5 4.1. Parameter tuning As mentioned before, the initial parameters of hybrid PSO-PSA algorithms (including both PSO-PSAІ and PSO-PSAП) are the population size (Npop), the number of successive iterations in which the best solution does not change (M), the cognition learning factor (C1), The social factor (C2) the inertial weight (w), the population size (Npop), the number of internal loop (In loop), and the temperature decreasing rate (α). Also, the population size (Npop), the number of successive iterations in which the best solution does not change (M), the crossover probability (Pc), and the similarity coefficient (SC), and the mutation probability (Pm) the number of internal loops (In loop), and the temperature decreasing rate (α) were used in GA-PSA.To investigate the influence of those parameters on the performance of the algorithms, we implement the Taguchi’s method in the design of experiments (DOE) (Montgomery, 2005).In the Taguchi’s method, the results are converted into an estimator called single to noise ratio (S/N). The S/N ratio shows both the mean and the variation in the response variables. To minimize the objective function the S/N ratio is calculated as the following formula: 𝑆 𝑁⁄ = −10log ⁡(1 𝑛⁄ ∑ 1𝑦𝑖2𝑛𝑖=1 )⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡(23) Which, n and yiindicate number of replications and process response value at i’th replication. We chose the Journal of Optimization in Industrial Engineering Vol.12, Issue 2, Summer & Autumn 2019, 133- 149 141 orthogonal array L27 for both PSO-PSA and GA-PSA. In Taguchi’s design of experiments, the relative percentage deviations (RPD’s) transform into S/N ratio, we figure out the response value of each parameter in the proposed algorithms.After testing the different parameter levels and analysing the obtained results, the better initial parameter levels were gained. The initial and optimumparameter levels of each proposed hybrid algorithm are shown in Tables 2.According to the parameter values in Table 2, we illustrate the trend of each factor level of PSO-PSA and GA-PSA are in Figures 10 and 11, respectively. Also, the results of each test problem solved by the proposed algorithms are shown in Figure 12. The statistical results of the objective function (makespan)and CPU time for the small and large sized problems are presented in Tables 5 and 6.Additionally, we utilized the relative percentage deviation (RPD) to test the performances of applied metaheuristics as below: 𝑅𝑃𝐷 = 𝐴𝑙𝑔𝑠𝑜𝑙 − 𝐿𝑠𝑜𝑙𝐿𝑠𝑜𝑙 ⁡⁡⁡⁡⁡⁡⁡⁡𝑖 = 1,…,𝑛⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡(32) Where, Algsolis the objective function value for a given algorithm, Lsolis the best value of the objective function between algorithms and n is the number of small size or large size problems. Table 2 Initial and optimum parameter levels of hybrid algorithms Algorithms Parameters Factor levels Optimum amount 1 2 3 PSO-PSA Npop PSO-PSA 100 200 300 200 M 10 20 30 20 C1 1 1.5 2 1.5 C2 0.7 1 1.5 1 Ω 0.3 0.4 1 0.4 In loop 5 10 15 10 α 0.87 0.91 0.95 0.91 GA-PSA Npop GA-PSA 100 200 300 200 M 10 20 30 20 Pc 0.85 0.9 0.95 0.9 SC 0.6 0.7 0.8 0.7 Pm 0.005 0.01 0.015 0.01 In loop 5 10 15 10 α 0.87 0.91 0.95 0.91 Fig. 10. Factor level trend of PSO-PSAalgorithm Fig. 11. Factor level trend of GA-PSAalgorithm Sadigh Raissi et al. / Three Hybrid Metaheuristic Algorithms… 142 Table 5 Solutions found by the algorithms of the small-sized test problems CPLEX PSO-PSAП PSO-PSAІ GA-PSA Problem # Budget type NOJ (j) NOS (s) Makespan CPU time Gap (%) Makespan CPU time Makespan CPU time Makespan CPU time 1 1 2 4 2 205 210 22 22 0 0 205 210 22 23 205 229 26 37 205 233 29 39 2 1 2 209 217 35 36 0 0 209 218 37 35 216 225 52 47 219 228 55 49 3 1 2 3 214 240 95 91 0 0 214 241 97 94 215 248 109 117 227 253 115 122 4 1 2 218 223 155 150 0 0 218 224 156 152 227 237 160 182 249 250 163 186 5 1 2 4 232 236 185 184 0 0 235 238 191 187 243 246 305 287 247 249 321 328 6 1 2 249 250 235 234 0 0 249 252 238 236 257 264 398 382 260 262 415 402 7 1 2 5 2 237 241 645 647 0 0 238 245 651 650 248 252 1039 1077 251 255 1100 1138 8 1 2 253 262 822 830 0 0 256 268 831 833 265 273 1005 1015 270 276 1176 1165 9 1 2 3 270 280 882 890 0 0 274 282 890 893 280 292 996 1016 284 298 1034 1086 10 1 2 289 304 918 922 0 0 292 308 923 931 302 317 987 1083 305 321 1150 1154 11 1 2 4 335 348 1012 1048 0 0 337 350 1024 1055 346 355 1056 1088 350 369 1149 1047 12 1 2 326 351 943 950 0 0 329 353 951 958 336 363 1008 1032 339 368 1062 1094 13 1 2 6 2 349 362 1083 1116 0 0 352 366 1099 1128 361 374 1125 1176 367 379 1139 1182 14 1 2 351 366 1043 1055 0 0 355 371 1074 1094 363 381 1133 1156 369 385 1273 1258 15 1 2 3 350 372 1103 1144 24.54 28.55 354 376 1127 1185 367 392 1176 1192 369 395 1201 1282 16 1 2 347 360 1078 1123 29.17 31.47 351 363 1112 1175 358 371 1185 1211 363 374 1246 1267 17 1 2 4 361 367 1197 1255 68.35 66.92 364 371 1225 1292 373 379 1312 1332 377 383 1424 1487 18 1 2 358 394 1205 1270 74.18 78.53 363 398 1243 1293 370 384 1294 1345 376 399 1431 1481 Average 1 2 268.28 299.06 703.22 720.39 10.90 51.37 288.61 301.89 716.17 734.11 296.22 310.11 798.11 820.33 301.50 315.39 860.17 875.94 Gap (%) = CPLEX optimality gap Fig. 12. The performance of proposed hybrid algorithmsof the large-sized test problems in terms of average makespan 0 2000 4000 6000 8000 10000 12000 14000 16000 1 4 7 10 13 16 19 22 25 28 31 34 37 40 43 46 49 52 55 58 61 64 67 70 73 76 79 82 85 88 91 94 97 10 0 10 3 10 6 10 9 11 2 11 5 11 8 A v e r a g e m a k e sp a n Test problem PSO-PSAП PSO-PSAІ GA-PSA Journal of Optimization in Industrial Engineering Vol.12, Issue 2, Summer & Autumn 2019, 133- 149 143 Table 6 Solutions found by the algorithms of the large-sized test problems Problem # Budget type NOJ (j) NOS (s) PSO-PSAП PSO-PSAІ GA-PSA Makespan CPU time Makespan CPU time Makespan CPU time 1 1 2 20 3 958 1005 2231 2277 937 996 2265 2309 1009 1055 2319 2393 2 1 2 941 994 2226 2202 959 1002 2261 2278 992 1038 2318 2373 3 1 2 964 1032 2230 2282 988 1042 2278 2266 1031 1086 2302 2285 4 1 2 984 1034 2236 2224 981 1034 2279 2273 1035 1088 2333 2383 5 1 2 952 993 2227 2265 985 1071 2291 2281 993 1039 2329 2382 6 1 2 6 1127 1179 2385 2463 1138 1197 2450 2408 1161 1231 2483 2479 7 1 2 1152 1223 2397 2467 1131 1174 2438 2464 1201 1276 2508 2579 8 1 2 1080 1147 2395 241 1097 1142 2428 2426 1119 1165 2469 2486 9 1 2 1188 1252 2373 2434 1149 1223 2443 2406 1245 1297 2471 2532 10 1 2 1164 1218 2380 2379 1170 1228 2431 2426 1104 1159 2468 2464 11 1 2 9 1174 1210 2417 2490 1178 1234 2458 2524 1216 1278 2529 2531 12 1 2 1120 1190 2416 2490 1155 1225 2471 2450 1178 1249 2521 2480 13 1 2 1197 1260 2433 2476 1233 1304 2443 2441 1257 1307 2471 2457 14 1 2 1143 1215 2421 2437 1145 1213 2465 2495 1169 1233 2500 2549 15 1 2 1187 1250 2417 2412 1233 1293 2476 2485 1236 1309 2511 2521 16 1 2 50 3 2617 2706 2903 2990 2684 2793 3021 3030 2767 2912 3197 3137 17 1 2 2652 2758 2911 2945 2730 2881 3028 3058 2813 2936 3134 3216 18 1 2 2562 2693 2910 2956 2639 2753 3074 3124 2697 2799 3251 3339 19 1 2 2635 2740 2908 2952 2697 2813 3010 3013 2768 2916 3159 3130 20 1 2 2589 2676 2901 2936 2680 2789 3050 3042 2757 2858 3180 3270 21 1 2 6 2839 2992 3010 3097 2920 3035 3141 3085 3055 3218 3300 3378 22 1 2 2815 2922 3009 3085 2918 3036 3115 3075 2994 3106 3240 3233 23 1 2 2620 2725 3010 3103 2698 2801 3116 3146 2770 2873 3255 3319 24 1 2 2843 2996 3031 3019 2933 3035 3199 3286 3050 3169 3352 3444 25 1 2 2568 2687 3022 3073 2637 2727 3184 3260 2750 2860 3300 3269 26 1 2 9 2827 2943 3181 3178 2911 3024 3324 3427 3004 3153 3491 3571 27 1 2 2825 2989 3186 3225 2898 3064 3318 3399 3009 3187 3516 3460 28 1 2 3109 3296 3164 3137 3214 3373 3288 3273 3325 3478 3399 3465 29 1 2 2827 2948 3181 3203 2931 3043 3280 3369 3051 3210 3390 3338 30 1 2 3087 3245 3193 3292 3178 3322 3344 3380 3299 3489 3517 3476 31 1 2 60 3 4117 4388 3578 3690 4525 4793 3834 3926 4900 5726 4070 4084 32 1 2 3861 4059 3596 3616 4194 4466 3837 3763 4468 4823 4035 4125 33 1 2 4026 4327 3594 3470 4304 4684 3847 3786 4752 5026 4097 4179 34 1 2 3828 4043 3586 3670 4192 4452 3876 3968 4507 4910 4134 4243 35 1 4048 3589 4424 3839 4787 4112 Sadigh Raissi et al. / Three Hybrid Metaheuristic Algorithms… 144 2 4412 3694 4675 3840 5185 4207 36 1 2 6 4406 4650 3709 3692 4712 5107 3999 3922 5205 5487 4221 4180 37 1 2 4331 4711 3729 3764 4658 4993 4025 3975 5161 5487 4246 4230 38 1 2 4123 4479 3701 3788 4403 4675 4030 4128 4860 5256 4351 4483 39 1 2 4250 4513 3700 3809 4570 4916 3912 3988 5072 5392 4240 4342 40 1 2 4236 4552 3722 3745 4641 4982 3916 3943 5060 5322 4157 4256 41 1 2 9 4668 4927 3954 4002 5008 5350 4221 4189 5415 5760 4494 4638 42 1 2 4443 4749 4024 3961 4793 5051 4248 4177 5155 5592 4499 4553 43 1 2 4249 4607 3947 4001 4626 4958 4156 4232 4943 5218 4431 4511 44 1 2 4705 4857 3958 3937 5156 5283 4227 4232 5302 5673 4449 4477 45 1 2 4569 4985 3949 4037 4761 4912 4226 4168 5018 5325 4554 4537 46 1 2 100 3 9395 10037 4845 4805 10169 11081 5459 5445 11358 12081 6083 6024 47 1 2 9416 10112 4857 5003 10344 11193 5313 5335 11573 12201 6022 5913 48 1 2 9712 10299 4799 4866 10521 11197 5446 5561 11681 12717 5957 5860 49 1 2 9279 9762 4811 4962 10184 10869 5301 5307 10896 11522 5982 5906 50 1 2 9137 9880 4824 4756 9879 10651 5423 5410 11125 11700 6040 6227 51 1 2 6 9938 10733 4968 5022 11013 11748 5629 5732 11839 12457 6274 6258 52 1 2 9725 10589 4972 4932 10744 11570 5587 5544 11916 12867 6240 6367 53 1 2 9747 10378 4920 5059 10919 11581 5415 5340 11947 12926 6167 6248 54 1 2 10115 10917 4937 4948 11005 11661 5566 5489 12176 13142 6305 6327 55 1 2 9844 10468 5017 4972 10732 11580 5479 5498 11919 12673 5978 6091 56 1 2 9 10733 12093 5428 5366 11540 12239 6568 6677 12597 13491 7390 7630 57 1 2 11007 11613 5799 5817 12329 13058 6558 6643 13316 14096 7855 7953 58 1 2 10873 11434 5849 5997 12069 12826 6557 6438 13537 14717 7462 7577 59 1 2 11654 12466 5892 6020 13002 13705 6694 6782 14018 15217 7583 7467 60 1 2 10650 10987 5824 5742 11537 12587 6502 6697 12880 13879 7893 7834 Average 1 2 4547.18 4842.08 3579.70 3617.55 4900.02 5211.50 3850.98 3867.98 5307.42 5656.53 4158.90 4193.93 NOJ (j) =number of jobs NOS (s) =number of stations After solving all test problems, the results showed that in the all test problems, PSO-PSAПhas better performances in terms of makespan and CPU time. The comparisons of makespan show that the PSO-PSAП provides better solution quality with difference average RPD (ARPD) of 5.24, and 12.02 in comparison to PSO-PSAІ, and GA- PSA, respectively. Also, we can see that in the comparison of CPU times, PSO-PSAП give better results in terms of difference ARPD of 6.04, and 13.42 against PSO-PSAІ, and GA-PSA. Figures 13 and 14 show the mean plot of the CPU time and makespan of the proposed metaheuristic algorithms, respectively.According to computational results, it can be inferred that in small- sized test problems, all the algorithms approximately have the same computational effort. Therefore, it can be proved that all the metaheuristic algorithms are able to reach optimal/near optimal solutions.But, as the sizes of the problems are increased, the difference between algorithms is more revealed so that PSO-PSAП overcomes all other algorithms. However, the PSO-PSAІ, and GA-PSA are highly affected by the problem size so that by increasing it, the CPU time is exponentially increased. Figure 15 depicts the 95% confidence intervals of makespan and Figure 16 showsthe 95% confidence intervals for RPD among the proposed algorithms. Journal of Optimization in Industrial Engineering Vol.12, Issue 2, Summer & Autumn 2019, 133- 149 145 Fig. 13. Means plot of CPU time of hybrid algorithms Fig. 14. Means plot of makespan of hybrid algorithms 5. Conclusion In this study, we developed the stochastic flexible flow shop scheduling problem (SFFSSP) considering fixed interval preventive maintenance (PM) and budget constraint. This new type of problem which is the main contribution of the research is presented an integrated mathematical model which is capable to consider all of the influential factors on makespan. In theproposed SFFSSP, there is a buffer between each stage, if all machines are busy or under PM action, the job can wait in the buffer. The budget constraint controls the holding costs of total jobs in all buffers should not be greater than available budget. The proposed model considers not only the preventive maintenance, but also stochastic processing time and budget constraint. By integrating all of these subjects the model will reflect the real performance of a SFFSSP. Since the problem was strongly NP-hard, three hybrid algorithms, including PSA with two types of PSO algorithms (PSO-PSAІ and PSO-SAП), and a GA (GA- PSA) were proposed to solve the model, which have suitable quality solutions in the literature.To compare the computational results, we tested two types of test problems containing 18 small sized test problems and 60 large sized test problems. As the results showed, the PSO- PSAП algorithm provided better quality solutions in both makespan and CPU time among the test problems. Also, the higher performance of proposed PSO-PSAП algorithm with respect to other algorithms is more revealed in the large sized test problems. The presented model is still open to considering other options, such as sequence dependent setup time, machine random breakdown, and the problem of job availability at time zero. Also, it might be exciting in working on bi- objective SFFSSP’s, which the other objective function could be minimizes the maximum tardiness.Another research direction could be incorporating different transportation types to transport jobs between each stage. Additionaleffort can try to solve model by developing a new solution methodology such as a new hybrid algorithm or a new population-based algorithm can be investigated.We assumed that each job has a same holding cost at each intermediate buffer, but it is different at each stage. Another aspect deserving future efforts is to consider that the holding costs of each job are different at any intermediate buffer. Fig 15. The 95% confidence intervals of makespanof the small- sized test problems Fig 16. The 95% confidence intervals of makespanof the large- sized test problems References Akrami, B., Karimi, B., Hosseini, S.M.M., (2006). Two metaheuristic methods for the common cycle 0 5 10 15 20 25 30 20 50 60 100 R P D NOJ PSO-PSAП PSO-PSAІ GA-PSA 0 5 10 15 20 25 20 50 60 100 R P D NOJ PSO-PSAП PSO-PSAІ GA-PSA 0 50 100 150 200 250 300 350 CPLEX PSO-PSAП PSO-PSAІ GA-PSA R P D 0 2 4 6 8 10 12 14 16 PSO-PSAП PSO-PSAІ GA-PSA R P D Sadigh Raissi et al. / Three Hybrid Metaheuristic Algorithms… 146 economic lot sizing and scheduling in flexible flow shops with limited intermediate buffers: the finite horizon case.Applied Mathematics and Computation, 183, 634- 645. Al-Hinai, N., ElMekkawy, T.Y., (2011). Robust and stable flexible job shop scheduling with random machine breakdowns using a hybrid genetic algorithm.International Journal of Production Economics, 132, 279-291. Almeder, C., Hartl, R.F., (2013). A metaheuristic optimization approach for a real-world stochastic flexible flow shop problem with limited buffer.International Journal of Production Economics, 145, 88-95. Arnaout, J.P., (2014). Rescheduling of parallel machines with stochastic processing and setup times.Journal of Manufacturing Systems, 33 (3), 376-384. Brucker, P., Kramer, B., (1995). Shop scheduling problems with multiprocessor tasks on dedicated processors. Annals of Operations Research, 50, 13– 27. Choi, S.H., Wang, K., (2012). Flexible flow shop scheduling with stochastic processing times: A decomposition-based approach.Computers & Industrial Engineering, 63, 362-373. Fahmy, S.A., Sherif, A., Balakrishnan, S., ElMekkawy, T.Y., (2009). A generic deadlock-free reactive scheduling approach.International Journal of Production Research, 47 (20), 5657-5676. González-Neira, E.M., García-Cáceres, R.G., Cabellero- Villalobos, J.P., Molina-Sánchez, L.P., Montoya- Torres, J.R., (2016). Stochastic flexible flow shop scheduling problem under quantitative and qualitative decision criteria.Computer and Industrial Engineering, 101, 128-144. Gupta J.N.D., (1988). Two stage hybrid flow shop scheduling problem. Journal of Operational Research Society, 39(4), 359-64. Hoogeveen, J.A., Lenstra, J.K., Veltman, B., (1996). Minimizing the makespan in a multiprocessor flow shop is strongly NP-hard. European Journal of Operational Research, 89, 172-175. Kennedy, J.,Eberhart,R.C.,(1995). Particle swarm optimization.Proceeding of IEEE International Conference on Neural Network, Piscataway: IEEE 4, 1942-1948. Kianfar, K., FatemiGhomi, S.M.T., OroojlooyJadid, A., (2012). Study of stochastic sequence-dependent flexible flow shop via developing a dispatching rule and a hybrid GA.Engineering Applications of Artificial Intelligence, 25, 494-506. Koulamas, C., Kyparisis, G.J., (2000). Scheduling on uniform parallel machines to minimize maximum lateness.Operations Research Letters, 24(6),175-179. Li, J.Q., Pan, Q.K., (2015). Solving the large-scale hybrid flow shop scheduling problem with limited buffers by a hybrid artificial bee colony algorithm. Information Sciences, 316:487–502. Lin, SW., Ying, K.C., (2013). Minimizing makespan in a blocking flowshop using a revised artificial immune system algorithm.Omega, 2 (4), 383-389. Lin, J.T., Chen, C.M., (2015). Simulation optimization approach for hybrid flow shop scheduling problem in semiconductor back-end manufacturing.Simulation Modeling Practice and Theory, 51, 100-114. Montgomery, D.C. (2005). Design and analysis of experiments. Arizona, John Wiley and Sons. Norman, B.A., Bean, J.C., (1999). A genetic algorithm methodology for complex scheduling problems. Naval Research. Logistics, 46: 199-211. Osman, I.H., Kelly, J.P., (1996). Metaheuristics: an overview. Metaheuristics: theory and applications. Boston: Kluwer Academic Publisher, 1-21. Pinedo, M., Chao, X., (1999). Operations scheduling with applications inmanufacturing and services. New York: McGraw-Hill. Pinedo, M., (2002). Scheduling theory, algorithms, and systems. Englewood Cliffs, New Jersey: Prentice- Hall. Chapter 2. Poli, R., Kennedy, G., Blackwell, T., (2007). Particle swarm optimization: an overview. Swarm Intelligence, 1(3), 33-57. Rabiee, M., Sadeghi Rad, R., Mazinani, M., Shafaei, R., (2014). An intelligent hybrid metaheuristic for solving a case of no-wait two-stage flexible flow shop scheduling problem with unrelated parallel machine.International Journal of Advanced Manufacturing Technology, 71: 1229-1245. Rahmani, D., Heydari, M., (2014). Robust and stable flow shop scheduling with unexpected arrivals of new jobs and uncertain processing times. Original Research Article, 33 (1), 84-92. Rahmani, D., Ramezanian, R., (2016). A stable reactive approach in dynamic flexible flow shop scheduling with unexpected disruptions: A case study. Computers & Industrial Engineering, 98, 360-372. Sangsawang, C., Sethanan, K., Fujimoto, T., Gen, M., (2015). Metaheuristics optimization approaches for two-stage reentrant flexible flow shop with blocking constraint.Expert Systems with Applications, 42, 2395–2410. Singh, M.R., Mahapatra, S.S., (2012). A swarm optimization approach for flexible flow shop scheduling with multiprocessor tasks.International Journal of Advanced Manufacturing Technology, 62: 267-277. Sukkerd, W., Wuttipornpun, T., (2016). Hybrid genetic algorithm and tabu search for finite capacity material requirement planning system in flexible flow shop with assembly operations. Computers & Industrial Engineering, 97, 157- 169. Tang, D., Dai, M., Salido, M.A, Giret, A., (2016). Energy-efficient dynamic scheduling for a flexible flow shop using an improved particle swarm optimization. Computers in Industry, 81, 82–95. Journal of Optimization in Industrial Engineering Vol.12, Issue 2, Summer & Autumn 2019, 133- 149 147 Tran, T.H., Ng, K.M. (2011). A water flow algorithm for flexible flow shop scheduling with intermediate buffers.Journal of Scheduling, 14, 483-500. Villemeur, A., (1991). Reliability, availability, maintainability and safety assessment. USA: Wiley. [32] Wang, X., Tang, L., (2009). A Tabu search heuristic for the hybrid flowshop scheduling with finite intermediate buffers. Computers & Operations Research, 36 (3), 907–918. Wang, K., Choi, S.H., (2014). A holonic approach to flexible flow shop scheduling under stochastic processing times.Computers & Operations Research, 43, 157-168. Wardono, B., Fathi, Y., (2004). A Tabu search algorithm for multi-stage parallel machine problem with limited buffer capacities. European Journal of Operational Research, 155 (2), 380-401. Zabihzadeh, S.R., Rezaeian, J., (2015). Two metaheuristic algorithms for flexible flow shop scheduling problem with robotic transportation and release time. Applied Soft Computing, 40, 319-330. This article can be cited: Raissi, S., Rooeinfar, R. & Ghezavati, V. R. (2019) Three Hybrid Metaheuristic Algorithms for Stochastic Flexible Flow Shop Scheduling Problem with Preventive Maintenance and Budget Constraint Journal of Optimization in Industrial Engineering. 12 (2), 131-147. http://www.qjie.ir/article_543744.html DOI: 10.22094/JOIE.2018.242.1532 http://www.qjie.ir/article_543744.html