- Email: [email protected]
Comparison of new metaheuristics, for the solution of an integrated jobs-maintenance scheduling problem
Accepted Manuscript
Comparison of new metaheuristics, for the solution of an integrated jobs-maintenance scheduling problem Massimo Bertolini , Davide Mezzogori , Francesco Zammori PII: DOI: Reference:
S0957-4174(18)30799-1 https://doi.org/10.1016/j.eswa.2018.12.034 ESWA 12375
To appear in:
Expert Systems With Applications
Received date: Revised date: Accepted date:
17 July 2018 8 November 2018 18 December 2018
Please cite this article as: Massimo Bertolini , Davide Mezzogori , Francesco Zammori , Comparison of new metaheuristics, for the solution of an integrated jobs-maintenance scheduling problem, Expert Systems With Applications (2018), doi: https://doi.org/10.1016/j.eswa.2018.12.034
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
HIGHLIGHTS The paper focuses on a single machine Jobs-Maintenance scheduling problem.
We enhanced the MHS by Zammori et. al (2014).
We also design a Cuckoo Search (CS) metaheuristics based on a new Lévy Flight.
Performance of the MHS improved, but this metaheuristic remains sub-optimal.
The newly developed CS is, instead, highly effective and robust.
AC
CE
PT
ED
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
Comparison of new metaheuristics, for the solution of an integrated jobs-maintenance scheduling problem Massimo Bertolini, Davide Mezzogori* and Francesco Zammori
CR IP T
Department of Engineering and Architecture, University of Parma, Italy Parco Aree delle Scienze, 181/A, 43124, Parma
Co-Author Prof. Massimo Bertolini
Phone: +39 0521 905705
*Corresponding Author Dr. Davide Mezzogori Email: [email protected]
Prof. Francesco Zammori
ED
Co-Author
M
Phone: +39 0521 905887
PT
Email: [email protected]
AC
CE
Phone: +39 0521 905887
AN US
Email: [email protected]
ACCEPTED MANUSCRIPT
Comparison of new metaheuristics, for the solution of an integrated jobs-maintenance scheduling problem ABSTRACT This paper presents and compares new metaheuristics to solve an integrated jobs-maintenance scheduling problem, on a single machine subjected to aging and failures. The problem, introduced by Zammori et al. (2014), was originally solved using the Modified Harmony Search (MHS) metaheuristic. However, an
CR IP T
extensive numerical analysis brought to light some structural limits of the MHS, as the analysis revealed that the MHS is outperformed by the simpler Simulated Annealing by Ishibuchi et al. (1995). Aiming to solve the problem in a more effective way, we integrated the MHS with local minima escaping procedures and we also developed a new Cuckoo Search metaheuristic, based on an innovative Lévy Flight. A thorough comparison confirmed the superiority of the newly developed Cuckoo Search, which is capable to find better solutions in
AN US
a smaller amount of time. This an important result, both for academics and practitioners, since the integrated job-maintenance scheduling problem has a high operational relevance, but it is known to be extremely hard to be solved, especially in a reasonable amount of time. Also, the developed Cuckoo Search has been designed in an extremely flexible way and it can be easily readapted and applied to a wide range of combinatorial problems.
M
Keywords - Cuckoo Search, Harmony Search, Metaheuristics, Scheduling, Single Machine, Planned Maintenance.
ED
1. INTRODUCTION
Scheduling is an important tool for manufacturers and production engineers, as it deals with the way in which human and machinery resources should be allocated to meet customers‟ demand. Specifically, given a
PT
set of production orders (or jobs), scheduling allows one to decide when and in which order jobs should be processed, with which resources and on which equipment. The objective is to control and to optimize
CE
workloads, to minimize production time and costs, respecting all production constraints (such as jobs‟ routing and jobs‟ releasing time) and, possibly, delivery due dates. Therefore, an optimized schedule may have a major impact on the performance of a manufacturing system and can be source of competitive
AC
advantage, as it allows reducing work-in-process, queuing times and industrial costs. It is worth noting that scheduling problems have always been of great interest also for the academia because, being very hard to be solved, they offer a challenging test bench to assess performance of sophisticated algorithms and/or metaheuristics (Allahverdi, 2015). As a matter of fact, scheduling models proposed in technical literature are becoming more and more complex, since current research‟s trends are toward integrating new and practical constraints in traditional scheduling models (El Khoukhi et al., 2017). In this regard, an issue of recent interest concerns the possibility to increase machines‟ availability and to reduce the overall production costs, by including Planned Maintenance tasks (PMs) into the jobs‟ schedule (Salmasnia et al., 2017). Traditionally, in fact, machines and equipment have been considered as fault free and the focus was limited to the way in which jobs should be sequenced and allotted to the available resources. Yet, this assumption is
ACCEPTED MANUSCRIPT rather heroic because machines can fail and, due to failures and reactive maintenance tasks (needed to restore the operating state of the fault machines) processing and waiting times gets longer, routings could need to be changed accordingly and, consequently, the original schedule may become non-optimal and should be reevaluated. Also, and perhaps more important, since breakdowns reduce productivity and increase costs, PMs should also be scheduled, as additional tasks, as an effective way to improve machines‟ availability (Safari and Sadjadi, 2011). The opportunity to define the scheduling sequence, taking into account possible failures and the need to perform PMs is becoming even more engaging, now that we are entering in the so-called
CR IP T
forth industrial revolution (Macchi et al., 2017). Indeed, due to the recent advances in manufacturing industry, which is leading toward the systematical deployment of Cyber-Physical Systems, it is now possible not only to remotely monitoring the system‟s state and to take actions based on real time data, but also to synchronize the physical factory floor and the cyber computational space (Lee et al., 2014), where optimization algorithms can be effectively used to dynamically solve the scheduling problem, in an acceptable amount of time.
AN US
Among the papers dealing with the Integrated Job Maintenance Scheduling problem (IJMS), this work focuses on the one by Zammori et al. (2014), who faced the problem of scheduling n jobs on a single machine subjected to aging and failures and proposed a Modified-Harmony-Search (MHS) metaheuristic to generate a quasi-optimal solution, in terms of earliness-tardiness penalties. The rational aspect of this choice is twofold:
The underlying hypotheses of the model are very pertinent with the industrial reality, as the single
M
machine is characterized by a time dependent (and increasing) failure rate and both reactive and
ED
planned maintenance tasks are executed to restore the functionality of the system. The model has a high practical relevance: although the single machine is the simplest form of shop
2013).
PT
floor configuration, it is also the one with the largest number of practical applications (Quan and Xu,
Concerning the last point, we note that the single machine can be used to adequately represent traditional
CE
continuous process plants (such as refining), but it can be easily readapted to model batch processes such as those of bakeries, biotechnology, chemical and pharmaceutical firms. Not only these production systems are
AC
very common, but it is exactly in these production environments that scheduling decisions are extremely important (to boost productivity) and where failures may have a dramatic impact, as they could interrupt the whole manufacturing process. Additionally, the single machine can also be used to schedule production in a pure flow shop or even in a job shop system, provided that the single machine is used to represent the real bottleneck of the manufacturing system. We also note that our interest (in this model) is mostly operational and computational and the aim is to introduce alternative solution approaches, more performing than the original MHS. The need to find better solutions was urged by the results of the analysis made by Zammori et al. (2014), who comparted the Modified Harmony Search (MHS) to the Hybrid Kangaroo Simulated Annealing by Soltani et al. (2010), a metaheuristic known to be particularly effective with NP-hard scheduling problems. Although MHS resulted
ACCEPTED MANUSCRIPT as the best approach, the gap was not as wide as one could expect, when comparing a sophisticated algorithm as the MHS with a simpler one. To shed light on this unexpected fact, we repeated the experimental campaign on a wider problem landscape and we observed that the solutions‟ space is rather flat (i.e., the average gradient is low), but it is disseminated of many local minima, with values close to the global one. As noted by Yang (2010), in these peculiar conditions performances of the Harmony Search (HS) get worse and the probability to get stuck, almost immediately, in a local minimum becomes high. Conversely, when the solution space is flat, only a search strategy assuring a high degree of randomness and little to none influence
CR IP T
to the initial starting point, may find the global optimum. To verify this assumption, we solved all the original problems‟ set using the very famous and general-purpose Simulated Annealing (SA) by Ishibuchi et al. (1995) that, due to a relatively simple iterated greedy strategy, fully satisfies the above-mentioned requirements. Obtained results were quite astonishing as the SA outperformed the MHS, in any operating condition.
Due to the above-mentioned issues, the speculative and operational goal of this paper is to ascertain if and
AN US
how this specific instance of the IJMS problem could be optimally solved also by a population-based metaheuristic. To this aim: (i) we modified the MHS, by introducing new local minima escaping techniques and (ii) we developed a new population-based metaheuristic that implements a searching strategy like the one of the Cuckoo Search (CS) by Yang and Deb (2009). Specifically, the CS has been chosen because, with respect to other population-based heuristics, it is less prone to early convergence, as cross-influence among
M
saved solutions plays a minor role in the evolutionary process.
The rest of the paper is organized as follows. Literature review on the IJMS problem is given in Section
ED
2, and the specific single machine IJMS problem considered in this paper is detailed in Section 3. New searching strategies and metaheuristics are introduced in Section 4 and 5 and their parameters are fine-tuned in Section 6. Furthermore, in Section 6 a comprehensive experimental campaign is made to assess and
PT
compare performance of the proposed metaheuristics. Lastly, conclusions and directions for future
CE
researches are drawn in Section 7.
2. LITERATURE REVIEW ON THE IJMS PROBLEM Scheduling problems have been generally treated under the assumption of full machines‟ availability.
AC
Only in the last decade this simplifying hypothesis has been questioned by several researchers, who postulated that, due to failures, the overall productivity of the manufacturing system may benefit from the inclusion of PMs between jobs. However, and quite surprisingly, almost the entirety of papers that start from this hypothesis do not fully consider the effect that PMs may have on machines‟ availability. Indeed, most of the researchers implicitly assume that, if properly performed, PMs will ensure sufficient machine availability to complete productions‟ orders in the allotted time. Basically, the problem is simplified as it follows: jobs‟ characteristics (in terms of processing times, set up times, release dates, routings, etc.) and the optimal number (or alternatively the optimal frequency) of the required PMs are considered as known inputs, and the scheduling problem reduces to the insertion of jobs between two consecutive PMs, in a way that optimizes a certain objective (such as the minimization of makespan or of the maximal tardiness), without violating (at
ACCEPTED MANUSCRIPT least within a certain tolerance limit) the expected PMs frequency. Yet, the possibility that faults may still occur is neglected, as well as the fact that failures‟ probability is strictly related to the time in which planned maintenance was performed. Among the most recent papers belonging to this category, the following ones can be mentioned. Mosheiov and Sarig (2009) extended the well-known scheduling problem of the “single machine with duewindows assignment”, by including in the schedule one additional maintenance activity. This maintenance activity is optional and, to be performed, it requires a fixed amount of time, during which production must be
CR IP T
halted. Although this maintenance task implies a waste of processing time, it may be beneficial because, once it has been performed, the single machine is supposed to become more efficient, and the successive workings will require a lower processing time. The authors also presented a polynomial time algorithm that generates the optimal schedule that minimizes costs related to jobs‟ earliness and tardiness and to the size and starting time of the due-windows.
Low et al. (2010) considered the problem of scheduling jobs on a single machine subject to periodic
AN US
maintenance activities. In this case, both the planned start time and the duration of each maintenance activity are known in advance, and the problems is to insert jobs between consecutive PMs, trying to minimize the makespan. To this aim the authors proposed an innovative version of the Particle Swarm algorithm that makes use of: (i) a “job-to-position” representation for the particles and (ii) the largest positional value rule to generate a particle sequence.
M
Zhu and Zhou (2015) also considered the single machine scheduling problem, but the peculiarity is that they introduced a deteriorating and resource-dependent “rate-modifying” planned maintenance. As usual, the
ED
hypothesis made is that planned maintenance alters (i.e., increases) the productivity of the machine. However, on the other one hand, planned maintenance is considered “deteriorating”, as its duration depends both on the resources allocated to it and on its starting time, i.e., the lateer it is performed, the worse the
PT
conditions of the machine become, and the longer the maintenance duration becomes. The authors also demonstrated that this problem can be optimally solved in a polynomial time with respect to many objective
date).
CE
functions (i.e., makespan, flowtime, maximum tardiness and combination of earliness, tardiness and due-
A similar problem was faced by Luo et al. (2015), who assumed that PM must start before a given deadline
AC
and that its duration increases with its starting time. The authors also provided polynomial-time algorithms to solve the problem with respect to the same objective functions used by Zhu and Zhou (2015). Similarly, Mor and Mosheiov (2015) studied the single machine due-windows assignment problem, in case of a deteriorating maintenance task. Specifically, two different types of deterioration were considered: timedependent (where the maintenance time increases as a function of its starting time), and position-dependent (where it is a function of its position in the sequence). In both scenarios a polynomial-time O(n4) optimisation algorithm is proposed to minimize costs related to jobs‟ earliness and tardiness and to the size and starting time of the due-windows.
ACCEPTED MANUSCRIPT Yu and Seif (2016) considered a flow shop scheduling problem, with the peculiarity that the operating conditions of each machine is expressed by a set of different Maintenance Levels (MLs). When a job is processed on a machine, each MLs is decreased and, as soon as a ML drops to zero, PM is performed. With the objective of minimizing tardiness and maintenance costs, the problem becomes NP-hard and it is solved with a lower-bound based genetic algorithm. Lee and Wang (2017) considered a single machine that processes two types of jobs, and the objective is to minimize the makespan of Type #1 jobs, provided that the maximum tardiness of Type #2 jobs does not
CR IP T
violate a predefined threshold limit. Additionally, PMs must be completed within specific time windows. Also, in this case, the number of PMs and their (deterministic) duration are known in advance and the problem is solved using branch and bound (based on a set of dominance rules) and via genetic algorithms. Similarly, also Touat et al. (2017) considered the problem of scheduling jobs on a single machine that requires PMs. In this case the specificity is that PMs must be performed inside predefined time windows, under human resources competences and availability constraints. To minimize the total tardiness, the authors
AN US
proposed two fuzzy genetic algorithms, differing in the way in which jobs, PMs and human resources are scheduled.
Cui et al. (2017) considered the single machine scheduling problem with jobs‟ release dates. The single machine requires flexible PMs that are activated following a “run-based” policy, i.e., a maintenance task must be performed any time the machine working age exceeds a predefined threshold value. To solve the
M
problem minimizing the makespan, a performing heuristic (based on some properties of the optimal solution) and a branch and bound (that utilizes several dominance rules) are proposed for large and medium sized
ED
problems, respectively.
Yazdani et al. (2017) considered a single machine scheduling problem characterized by multiple unavailability periods, occurring at known and constant intervals, during which overhauls are performed.
PT
Aiming to minimize the sum of earliness-tardiness penalties, the authors proposed an algorithm characterized by a variable neighbourhood search that, in addition to local search, also implements a “shake” procedure,
CE
which diversifies solutions, by performing a leap toward another local neighbourhood. The algorithm is also integrated with a “knowledge module” that keeps track of the structure of good solutions and that takes advantage of it during the search.
AC
El-Khoukhi et al. (2017) investigated the flexible job shop scheduling problem, when a predefined number of PMs, of known and deterministic length, are needed. To minimize makespan, they solved the problem using an innovative Dual-Ants Colony Optimisation. The metaheuristic is also integrated with a dynamic memory, which takes advantage of local search strategies and several dispatching rules. Woo and Kim (2018) considered the parallel machines scheduling problem characterized by “ratemodifying” PMs activities. Specifically, each machine is characterized by a starting processing time that increases through time, accordingly to a certain deterioration rate. However, if PMs is performed on a machine, the processing time is restored to its original value. Given a certain number of time slots where PMs can be performed and considering as known and deterministic the PMs‟ time, the main decisions to be
ACCEPTED MANUSCRIPT taken concern: (i) the number and position of PMs and (ii) the ordered assignment of jobs on parallel machines. Specifically, to minimize the makespan, the authors proposed a specific heuristic that takes advantage of Simulated Annealing and of Genetic Algorithm to decompose the problems into sub-problems than can be analytically and optimally solved. As we have already said, all these papers consider PMs as a prerequisite needed to assure a minimum level of availability, as required by the manufacturing system. Therefore, the number and the duration of the PMs are assumed to be known and the positive effect that PMs may have on the hazard rate of the machines
CR IP T
is disregarded. Very few papers reduced these drawbacks, by considering a time-dependent hazard rate and by including in the analysis also the additional losses that may be generated by faults and by reactive maintenance tasks. As far as the knowledge of the authors, in addition to the work by Zammori et al. (2014), which will be discussed in depth in Section 3, this approach was followed only by eight recent papers that share the following common hypotheses: (i) Machines‟ deterioration increases through time and it is described by a time dependent hazard rate, (ii) when performed, planned maintenance restores the original
AN US
operating conditions of the machine (i.e., perfect repair) or, at least, it reduces the current degradation state, (iii) although PMs are performed, failures may still occur. In this case reactive maintenance tasks are needed to restore the machines, nonetheless their degradation remains unchanged (i.e., minimal repair model). Safari et al. (2011) are, probably, the first authors that fully considered the effects of PMs on the hazard rate. Specifically, the authors deal with a flow-shop configuration under the assumption of condition-based
M
maintenance to minimize expected makespan. More precisely, every T time units, each machine is inspected and, depending on its current degradation (modelled with a Poisson distribution) a reparation may be
ED
performed. Also, failures can occur in between two consecutive inspections, according to a nonhomogeneous Poisson process of parameter
. A hybrid heuristic based on genetic algorithm and
simulated annealing is also proposed to minimize makespan.
PT
Cui et al. (2014) tackles the problem of scheduling jobs and PMs on a single machine subject to unexpected failures. The objective function is evaluated both in terms of solution robustness and of solution quality.
CE
Specifically, the first metric is computed as the total deviation between the planned start time and the actual start time of the jobs, whereas the second metric is measured as the maximum tardiness of the jobs. Also, to generate optimal solutions, a three-phased heuristic is proposed. At first the jobs‟ sequence and the relative
AC
position of PMs is determined; next some extra idle time (or “time buffer”) is included in the schedule to protect due dates from unpredictable failures. Lastly, in the third and final phase, a pair-wise-swap-based local search is performed on the jobs‟ sequence, aiming to improve the initial solution. Lu et al. (2015) tackled the problem of finding a robust and stable schedule for a single machine that suffers of unexpected breakdowns. According to the authors, the best way to maximise robustness is to: (i) add time buffers to protect the initial schedule and (ii) execute PMs anytime machine‟s reliability drops below a threshold limit. Owing to these issues, the authors solved the problem with a genetic algorithm based on a double chromosome: the first chromosome codifies (with natural number) the jobs‟ sequence, the second one codifies (in binary code) the desired reliability limit.
ACCEPTED MANUSCRIPT Mokhtari et al. (2015) considered a flexible flow-shop scheduling problem, in which machines‟ failure rates are time-varying and proposed a heuristic, based on Simulated Annealing and on Monte Carlo simulation, to minimize the number of tardy jobs. The job shop is considered flexible because any operation can be performed by more than one machine and failure rate is considered time-varying as it depends on environmental situations like temperature, light, humidity, etc. Quite interestingly, this model is based on the hypothesis that the time needed to complete a batch of orders is much smaller than the average time between two planned maintenance tasks. Consequently, the use of PMs is not contemplated and, indeed, only reactive
CR IP T
maintenance tasks are executed anytime a failure occurs. Nouiri et al. (2017) focused on the flexible job shop scheduling problems under machines breakdowns. Makespan minimization is considered as the main objective, but, as secondary objective, the obtained solution must guarantee a stable schedule even in case of breakdowns. To this aim a two-stages particle swarm optimization was proposed: at first breakdowns are disregarded and an initial sequence is found, next breakdowns are randomly generated (supposing that machines with the heavier workload are more likely to
AN US
fail) and the original solution is modified to assure robustness and stability. As in the work by Mokhtari et al. (2015), also in this case PMs are not considered, and only reactive maintenance is used to counteract machines‟ failures.
Salmasnia and Mirabadi-Dastjerd (2017) considered the scheduling problem of a degraded single machine subjected to failures. Specifically, they considered different levels of PMs: although planned maintenance
M
always leads to machine‟s rejuvenation, only if the top and most expensive kind of planned maintenance is performed the machine is restored to its original state (i.e., perfect repair). The objective is to jointly
ED
minimize the total maintenance cost and the total jobs completion time. Since there is a clear trade-off between these objectives, the Pareto Set of the non-dominated solutions was generated using a genetic algorithm.
PT
Rahmati et al. (2017) extended the flexible job shop scheduling problem presented by Yu and Seif (2016), by considering machines characterized by an increasing hazard rate. Briefly, to schedule jobs and PMs, a
CE
simulation-based optimization similar to the one by Nouiri et al. (2017) was proposed. At first, Harmony Search is used to find both the jobs‟ sequence and the machines‟ assignments that minimize the makespan; next the optimal time between PMs is defined via simulation.
AC
Cui et al. (2018) dealt with the integration of production scheduling and maintenance planning to optimize robustness and solution quality (evaluated as in Cui et al. (2014)), for a flow shop with failure uncertainty. As in Lu et al. (2015), buffer times are proactively added into the initial schedule and an initial solution is obtained following a constructive mathematical model. Next, the solution is improved using a two loops algorithm that takes advantage of numerical simulation and Monte Carlo sampling methods. More precisely, the outer loop performs a local search to optimize jobs‟ sequence and PMs‟ positions; the inner loop implements a genetic algorithm that tries to optimize the idle times when jobs‟ sequence and PMs‟ positions are fixed. Results of the inner loop are returned as an evaluation criterion to guide the search direction of the outer loop.
ACCEPTED MANUSCRIPT For the sake of clarity, results of the literature review on IJMS problems published starting from 2014 are summarized in Table I. Although the considered models are case specific and cannot be directly compared,
AC
CE
PT
ED
M
AN US
CR IP T
we tried to summarize, in the table, pros and cons of each one of them.
CR IP T
ACCEPTED MANUSCRIPT
Table I Main papers dealing with the IJMS problem from 2014 to now Paper
System
Objective Function
Features of PMs
Solution Approach
Cui et al. (2014)
Single Machine
PMs executed at variable time intervals to increase machine‟s availability.
Schedule Robustness and Quality
Zammori et al. (2014)
Single Machine.
PMs executed at variable time intervals to increase machine‟s availability.
Earliness and Tardiness Penalties.
Luo et al. (2015)
Single Machine
PM is used to improve machine‟s productivity; PM is “deteriorating”, as is duration increases over time.
Makespan; Other classic objectives.
Lu et al. (2015)
Single Machine
PM is executed any time drops below a threshold limit, to restore its initial state.
Robust Schedule Stable Schedule
Mokhtari et al. (2015)
Flexible Flow Shop
Machines are characterized by a time-varying failure rate. When faults occur, maintenance is needed to reactivate the system.
Number of tardy jobs
Mor and Mosheiov (2015)
Single Machine
PM is used to improve machine‟s productivity; PM is “deteriorating”, i.e. duration is time and position dependent.
Costs minimization.
Optimal Polynomial Time Algorithm.
Zhu and Zhou (2015)
Single Machine
Makespan; Other classic objectives.
Optimal Polynomial Time Algorithm.
Yu and Seif (2016)
Flow Shop
Tardiness and Maintenance costs
Genetic Algorithm
Cui et al. (2017)
Single Machine.
Makespan.
Specific Heuristic; Branch and Bound procedure.
Hazard rate is considered. Robustness is obtained using time buffers and the number of PMs is analytically optimized. Hazard rate is considered. Failures and reactive maintenance are considered. The number of PMs is analytically optimized. A deteriorating rate, affecting PM‟s time is considered; Optimal solution can be obtained in polynomial time Hazard rate is considered; Failures are considered; The optimal threshold limit for is optimized.
AN US Modified Harmony Search. Optimal Polynomial Time Algorithm.
Genetic Algorithm with double chromosome Simulated Annealing and Monte Carlo simulation
M
ED
AC
CE
PT
PM is used to improve machine‟s productivity; PM is “deteriorating”, as is duration increases over time. Machines degrade over time and their productivity decrease. PMs are needed as soon as the Maintenance Level drop to zero. The machine degrades and its productivity decreases. PMs are executed, anytime the degradation exceeds a threshold.
Three-phased heuristic
Pros
Time varying hazard rate depending on environmental conditions. Two kinds of deteriorating rates are considered, and several costs are considered. Optimal solution can be obtained in polynomial time. A deteriorating rate, affecting PM‟s time, is considered; Optimal solution can be obtained in polynomial time. Different Maintenance Levels can be used. The number and the kind of PMs to be executed is an output of the genetic algorithm. Introduces a “degradation” concept. PMs executed accordingly to a runbased policy to reduce degradation.
Cons
A single level of PM (i.e., perfect repair) is considered.
A single level of PM (i.e., perfect repair) is considered. Solution approach is not very performing. Hazard rate is not considered; Faults are not considered; A single maintenance task must be included in the schedule. A single level of PM (i.e., perfect repair) is considered.
Only faults and reactive maintenance tasks are considered. PMs are not considered as a way to improve availability. Hazard rate is not considered; Faults are not considered; At most one maintenance task can be included in the schedule. Hazard rate is not considered; Faults are not considered; At most one maintenance task can be included in the schedule. Hazard rate is not considered (i.e., it is substituted by a degradation rate); Faults are not considered; Hazard rate is not considered (i.e., it is substituted by a degradation rate); Faults are not considered;
ElKhoukhi et al. (2017)
Flexible Flow Shop.
Lee and Wang (2017)
Single Machine.
Nouiri et al. (2017)
PMs are needed to assure the required availability and do not alter machine productivity.
Dual Ant Colony with dynamic memory
PMs are needed to assure the required availability and must be performed within predefined time windows.
Makespan; Tardiness.
Genetic Algorithm.
Flexible Job Shop.
Fault probability depends on the work load; anytime a fault occurs, reactive maintenance is needed.
Makespan
Rahmati et al. (2017)
Flexible Job Shop.
PMs executed at variable time intervals to increase machine‟s availability.
Makespan.
Touat et al. (2017)
Single Machine.
PMs are needed to assure the required availability and must be performed within predefined time windows.
Total Tardiness.
Salmasnia and Mirabadi (2017)
Single Machine
PMs executed at variable time intervals to increase machine‟s availability.
Total Cost Total Compl. Time
Yazdani et al. (2017)
Single Machine
Overhauls are performed (stopping production) at predefine time intervals.
Earliness and Tardiness Penalties
Variable Neighborhood Search
Cui et al. (2018)
Flow Shop
PMs executed at variable time intervals to increase machine‟s availability.
Schedule Robustness and Quality
Dispatching Rules, Simulations and Monte Carlo Sampling
Woo and Kim (2018)
Parallel Machines
The machine degrades and its productivity decreases. PMs must be executed, within time windows to reduce degradation.
Makespan
Simulated Annealing and Genetic Algorithm
Optimization via simulation
Harmony Search and Simulation
Fuzzy Genetic Algorithms.
Genetic Algorithm with double chromosome
M
ED
PT
CE
AC
The model is standard, the innovative part lies in the solution approach that is very performing. Number and duration of PMs are known, but their starting time can vary. Jobs are classified into families with different importance.
AN US
Makespan
CR IP T
ACCEPTED MANUSCRIPT
Failure rate is time varying and it also depends on work load. Hazard rate is considered. Failures and reactive maintenance are considered. The number of PMs is analytically optimized. PMs must be placed within predefined time windows, that also depends on human resources competences and availability constraints. Hazard rate , failures and reactive maintenance are considered. Different levels of PM are considered, and their number is analytically optimized. The model is standard, the innovative part lies in the solution approach that is very performing. Hazard rate is considered. Failures and reactive maintenance are considered. The number of PMs is analytically optimized. Time buffers are used to increase robustness A machines‟ deteriorating rate is considered. The optimal number of time windows (and their duration) is optimized.
Hazard rate is not considered; Faults are not considered; Number and duration of PMs are known in advance. Hazard rate is not considered; Faults are not considered; Number and duration of PMs are known in advance. Only faults and reactive maintenance tasks are considered. PMs are not considered as a way to improve availability. The simulation model is case specific A single level of PM (i.e., perfect repair) is considered. Optimization is obtained via simulation and it is case specific. Hazard rate is not considered; Faults are not considered; Number and duration of PMs are known in advance. Due to tradeoff among the objectives, a global optimum cannot be found; only a Pareto Optimal frontier can be generated. Hazard rate is not considered; Faults are not considered; Number and duration of PMs are known in advance.
A single level of PM (i.e., perfect repair) is considered. Optimization is obtained via simulation and it is case specific. Hazard rate λ(t) is not considered (i.e., it is substituted by a degradation rate); Faults are not considered;
ACCEPTED MANUSCRIPT
3. THE 3.1. In this paper, we focus on the single machine IJMS problem introduced by Zammori et al. (2014), whose features are listed below:
The single machine must process a certain number of jobs, all available at the starting time
For each job, its processing time, due date and sequence dependent setup time are known;
The machine is subjected to failures and production may be halted both to perform reactive and
;
planned maintenance tasks; A time dependent hazard rate
is considered and minimal and perfect repair models are used to
CR IP T
describe machine‟s restoration after the execution of a reactive and of a planned maintenance task, respectively;
The goal is to find the integrated sequence of jobs and planned maintenance tasks that minimizes the sum of Earliness-Tardiness Penalties (ETP).
AN US
In short, using the Kendall‟s notation, the problem can be classified as: (1| Due Dates, Sequence Dependent Set-Up Times, Planned Maintenance| Earliness-Tardiness Penalties). 3.2. Computing the expected completion time of the jobs Let
and
be, respectively, the number of jobs and of PMs. Then, any possible solution can be
expressed as a combination of two vector and. The vector possible permutations of
{
jobs and each one of its element
M
the
(
be processed in position . Similarly,
represents one of } codifies the job that will
) represents one of the (
) possible ways
ED
in which PMs can be inserted between the first and the last job; in this case, each element of index of the job (within ) that precedes the j-th PM. For instance, with } would be denoted as:
,
PT
{
[
]
CE
Using this notation, the expected completion time for each job *
AC
∑(
[
[ {
)
(∑
and
codifies the , the schedule
.
+ can be expressed as in Eq. (1):
) (1)
]
∑
(
)
(
]
∑
(
)
(
∑
∑
)
)
ACCEPTED MANUSCRIPT
where: {
} is the set of jobs;
is the vector of the jobs‟ processing times; (
) is the
matrix of the sequence dependent setup times; is the vector of the due dates;
is the average time of a reactive maintenance task;
is the average time of a planned maintenance task;
is the expected number of failures in the interval [
].
CR IP T
According to Eq. (1) the expected completion time of a job can be computed using a recursive formula that takes into consideration the position of the job within the scheduling sequence, relatively to the nearest PM. More precisely, the whole schedule is seen as an aggregate of sub sequences of jobs divided by
AN US
consecutive PMs. Also, due to the perfect repair hypothesis, all these sequences can be considered independent because any time PM is performed, the hazard rate consequently, the interval [
is reset to its original value and,
], on which the expected number of failures
is computed, shrinks
accordingly. Owing to these issues, the completion time of a job belonging to a generic subsequence equals the sum of: (i) the completion time of the last subsequence
and (iii) the expected down time, due to both failures and
M
the preceding jobs in the same subsequence
, (ii) the processing and setup times of all
maintenance, from the start of subsequence up to the considered job.
becomes: (2)
can be easily obtained as:
CE
Hence,
shape parameter , the hazard rate
PT
with characteristic life
ED
Also, under the reasonable and common assumption of the time to failure being distributed as a Weibull,
( )
(3)
AC
Based on the expected completion times, ETP can be finally evaluated as in (4), where earliness and tardiness penalties (per unit of time) relative to job ∑(
and
are the
. )
(4)
where: ,
-
(5)
ACCEPTED MANUSCRIPT
|
,
-|
(6)
3.3. Solution approach Zammori et al. (2014) also proposed a Modified Harmony Search (MHS) to solve the (1| Due Dates, Sequence Dependent Set-Up Times, Planned Maintenance| ETP) scheduling problem. Basically, the MHS is an advancement of a standard Harmony Search (HS), a metaheuristic inspired by the creativity process of
CR IP T
jazz musicians that, during jam sessions, expand their repertoire using pure inspiration and/or reworking on old harmonies (Geem et al., 2001).
Formally, the HS has its essence in the Harmony Memory (HM), a matrix that represents a musician‟s repertoire. This matrix contains
row vectors
[
, where the superscript
] indicates the position in
the HM, each one coding a feasible solution and its Objective Function (OF). For the present problem, and
are joined to form an
whose generic element
(with
(
-dimensional vector
AN US
vectors
),
) contains the index of the i-th task of the integrated job-
maintenance schedule and the element
contains the ETP of the schedule. Since tasks may be either
jobs or maintenance, natural numbers in the interval [ jobs and PMs, respectively.
] and in the interval [
] denote
M
Randomly initialized, solutions of the HM follow an evolutionary process based, mainly, on the Harmony Memory Consideration (HMC). Specifically, new solutions
taken from some “promising” solutions stored in the HM. To avoid rapid
ED
a time, the elements
are generated by recombining, one element at
convergence, selection is made using a roulette wheel selection rule so that, although rarely, even the elements of a solution
with a poor objective function may be selected. For the same aim, two additional
PT
strategies, namely Random Selection (RS) and Pitch Adjustment (PA), are used to convey a certain degree of randomness to the generative process. According to the first strategy, to create a new solution, some may be randomly generated from the values of the domain [
] to which they belong
CE
elements
to. Similarly, PA is used to slightly perturb one or few elements of the new solution. For instance, in the
AC
MHS, PA is implemented as a pairwise interchange generation scheme, that is two positions and new solution
are randomly selected and the corresponding elements
and
of the
are swapped. Lastly, a new
solution is accepted if its objective function is better than that of the worst solution stored in the HM; in this case, the worst solution is replaced by the new one. The peculiarity of the MHS is the implementation of an additional local minima escaping strategy, called Large Portion Recovery (LPR). Anytime LPR is triggered, the HM is completely regenerated with new solutions created by recombining, through 2-points cross-over, entire sub-arrays of the precedent solutions contained in the HM. LPR is activated with a rate (LPRR) that depends on the Saturation‟s Level of the HM, a measure of similarity of the solutions of the HM, computed as in Eq. (7):
ACCEPTED MANUSCRIPT
(7) Where
and
are, respectively, the number of copies and the total number of solutions of the HM.
3.4. Problem peculiarities As anticipated in Section 1, when the ETP is used as objective function, regardless of others problem settings, the solution space becomes rather flat, but it is scattered by many points of local minimum. For the
CR IP T
sake of brevity, we will not report the detailed results of this preliminary analysis, but we will show a simple example, which illustrates this concept well. Specifically, Figure 1.a shows the ETP of all the 3600 solutions of three randomly generated instances of a 6-Jobs, 1-PM scheduling problem. In each problem, the average job‟s processing time equals 100-time units, planned and reactive maintenance equals 40 and 20-time units, respectively. What differentiates each case is the value of the shape parameter
of the Weibull probability
AN US
distribution, which equals 90%, 60% and 30% of the average jobs‟ processing time. ------------------------------Insert Figure 1.a here -------------------------------All contiguous points on the graph (i.e., points labelled as
,
) are neighbour solutions generated as
M
follows. At first, starting from the basic sequence {1, 2, 3, 4, 5, 6}, its 720 permutations were generated using the shifting procedure i.e., two locations if
and
are randomly selected and job
is moved before
and vice-versa. Next, all these permutations were replicated for five times and PM was added
ED
after job in position 1 for the first group of 720 permutations, in position 2 for the next 720 permutations, and
CE
PT
so on, as shown in the generation scheme of Table II:
,
… …
… ,
,
…
AC
,
Table II The neighbors‟ generation scheme
As it can be seen from Figure 1.a, the general pattern does not significantly change: neighbor solutions
tend to scatter, almost randomly, around an average value that does not display a significant trend. The fact that the solution space is rather flat is even better highlighted by Figure 1.b, which shows the ETP curves obtained after sorting (in ascending order) all the 3600 solutions of the problem. ------------------------------Insert Figure 1.b here --------------------------------
ACCEPTED MANUSCRIPT Lastly, we note that also the number of local minima is high, as shown by Table III, which reports the main features of the above-mentioned example.
Table III Basic features of the problems
= 60
= 30
Minimum Cost
10,094
13,455
19,838
Maximum Cost
13,067
17,861
26,608
Solutions within 10% from the optimum
404
703
698
Solutions within 5% from the optimum
82
165
170
AN US
3.5. Computational drawbacks
CR IP T
= 90
As noted by Yang (2010), when the solution space is so flat, neither Random Selection nor Pitch Adjustment can convey enough randomness to effectively escape an early convergence process: all the solutions of the Harmony Memory become similar, if not even identical, and the probability to get immediately stuck in a local minimum becomes high. Although the MHS implements a third local minimum
M
escaping strategy (i.e., the Large Portion Recovery), it is licit to ask if this approach is sufficient to obviate the early convergence problem. To answer this question, we compared MHS with the very famous and general-purpose Simulated Annealing (SA) by Ishibuchi et al. (1995), whose main features are detailed in
ED
Appendix A. The rational of this choice can be motivated as follows: The metaheuristic implements a relatively simple iterated greedy strategy that assures a high degree
PT
of randomness during the exploration of the solution space; The metaheuristic is flexible, as its performance does not depend, significantly, neither on the
CE
specific nature of the combinatorial problem under analysis, nor on the initial condition used as the starting point;
The number of operating parameters to be fine-tuned is very low.
AC
We also note that we opted for a general-purpose algorithm rather than for one specifically designed to tackle IJMS problems, as it was not possible to find an algorithm suitable for the specific IJMS under consideration. Indeed, as clearly discussed in Section 2 and showed in Table I, most of the approaches proposed in the literature have the limit not to consider the effect of PMs on the hazard rate of the machines, which is, instead, one of the peculiar feature of the model herein considered. Even the few works that considered the hazard rate, cannot be used as benchmark, since they are very specific and cannot be effectively readapted to the framework of the problem under consideration. Obtained results were astonishing, as SA outperformed MHS, even in case of 10 randomly generated problems of small size (with parameters defined as indicated in Section 6), characterized by 14-Jobs and 2-
ACCEPTED MANUSCRIPT PMs and for which the optimum was found through an exhaustive search. More precisely, each problem instance was solved 8 times and obtained results where quantified in terms of Average Percentage Error (APE) and Average Percentage Deviation (APD), computed as in Eq. (8) and Eq. (9), respectively. Note that quantifies the average difference between the best solution found by the j-th algorithm and the real optimal one. Similarly,
measures how much, on average, the j-th algorithm diverges from the real
optimum. (
Where:
(8)
(
)
(9)
AN US
∑∑
)
CR IP T
∑
is the k-th solution generated by the j-th algorithm for the i-th problem; in this case [
],
] and
[
];
is the best solution for problem found by algorithm i.e.,
{
};
is a priori known global optimum for problem .
M
[
In case of SA, both APE and APD were equal to zero; conversely, for MHS an APE of 0.5% and an APD of 1.8% were found. These values demonstrated that, whereas SA consistently found the global minimum,
ED
MHS always converged to some local minimum rather distant from the global one.
PT
4. ENHANCEMENTS OF THE MHS
Due to the preliminary results discussed in Section 3, we tried to improve the MHS through the addition
CE
of four additional local minima escaping strategies that will be detailed in the following Sub-Sections. 4.1.
Block Shifting: a new Pitch Adjustment procedure
As explained in Section 3, in the MHS PA has been implemented as a simple pairwise interchange
AC
operator. Consequently, due to the small number of PMs compared to the number of jobs, it is quite rare for the PA to swap a job with a PM. In general, this is not an issue and the evolving process works fine, as PMs repositioning is assured by the Harmony Memory Consideration. Yet, when the Harmony Memory begins to saturate, Harmony Memory Consideration loses effectiveness and the exploration of the neighborhood proceeds mainly thanks to PA. To exploit this search, it is thus fundamental to add to the original pairwise interchange, which modifies jobs‟ sequencing leaving almost unaltered the position of PMs, a dual PA operator that, conversely, acts on PMs leaving almost unaltered jobs‟ sequencing. Owing to these considerations, we defined the Block Shifting (BS) operator, which works as follows:
A whole subsequence of jobs (i.e., jobs placed between two PMs) is randomly selected;
ACCEPTED MANUSCRIPT
Let
be the selected subsequence, then a new neighbor solution is generated by moving the block
,
- either after the first job of the subsequence
subsequence
.
{
For instance, in case of sequence subsequence
}, if
} was selected and block {
{
or before the last jobs of the
{
} was moved forward, the sequence } would be finally obtained. Please note that
could have been generate also by means of pairwise interchange; yet this would be extremely improbable
CR IP T
because, considering that sequences with 2 PM in adjacent positions are not allowed, it would require a minimum of five consecutive moves. 4.2.
A second crossover operator
Aiming to improve the regeneration of the HM through Large Portion Recovery, in addition to the 2points crossover, we also considered and compared two alternative types of crossover: (i) the order crossover
AN US
and (ii) the double 2-1-points crossover, readapted from the recent one proposed by Salmasnia and MirabadiDastjerd (2017).
The ordered crossover was chosen to foster the exploration of large portions of the solution space, as it is known to have the ability to generate offspring solutions sensibly different from the parent ones (Gen et Cheng, 2000). Briefly, we recall that the order crossover is based on the following steps: Select a random swath of consecutive genes from Parent 1;
Drop the swath down to Child 1 and mark out these genes in Parent 2;
Starting on the right side of the swath, take genes from Parent 2 and insert them in Child 1, at the
ED
M
right edge of the swath;
Generate Child 2 by flipping Parent 1 and Parent 2 and repeating all the other steps.
highlighted in Italics:
{
AC
CE
{
PT
To clarify this process, a simple example is given next, where the swath of genes selected from Parent 1 is }, },
{
},
{
}.
The integrated 2-1-points crossover is the only alternative crossover operator that was found in the
literature, among the subset of papers that deal with the IJMS problem and that fully consider the effect of PM on machines‟ hazard rate (see Table II in Section 2, for further details). In its original version, 2-1-points crossover operates on two distinct chromosomes that codify jobs and different levels of PM, respectively. However, in the present problem there is a single kind of PM, and so we used the second gene to codify the position of PMs within the jobs‟ sequence. In virtue of this slight modification, the 2-1-points crossover operates as it follows:
Jobs and PMs are separated and coded with the
-
notation of Section 3.1;
ACCEPTED MANUSCRIPT
2-points crossover and 1-point crossover are made on the
vector (i.e., jobs part) and on the
vector (i.e., maintenance part), respectively;
The two constructed parts are recombined together to generate a new solution.
To clarify this process, the integrated 2-1-points crossover is applied to the parents used in the example above; note that both vectors have been recoded using the
notation and that the symbol ˈ|ˈ has been
-
used to highlight the cutting points. {
}
|
|
|
,
{
}
|
|
|
,
CR IP T
,
{
},
{
}.
Lastly, aiming to select the best crossover operator we: (i) developed three standard genetic algorithms,
AN US
differentiated only by the type of crossover and (ii) we tested them on the same problems‟ set described in Section 3.4. At the end of the analysis the choice fell on the ordered crossover, as it was the only one able to find the best solution for every instance of the problems set. 4.3.
Maximizing the Harmony Memory degree of differentiation
In the original MHS, anytime the Large Portion Recovery is invoked, crossover is activated to regenerate
M
the HM. Ideally, at the end of this process, the HM should be completely new, yet, as clearly shown by Figure 2, the true picture is rather different.
ED
------------------------------Insert Figure 2 here
PT
--------------------------------
As it can be seen, anytime the Saturation Level (SL) reaches a threshold limit (in this case set to 0.3), the Large Portion Recovery is triggered and, immediately after, the SL drops down; yet, unlike one could expect,
CE
it is always greater than zero. This fact implies that, due to an excessive degree of similarity among parents, duplicates are created in the offspring. Conversely, to maximize the degree of differentiations of the HM,
AC
without losing good solutions, we decided to limit the regenerative process only to the duplicated solutions, which are substituted with newly generated ones. More precisely, letting
,
and
be,
respectively, the number of solutions, of copies and of original (non-duplicated) solutions contained in the HM, the regenerative process works as follows:
Identify all different solutions contained within HM;
Use a roulette wheel selection (without replacement if select
, with replacement otherwise) to
parents among the original solutions and apply order-crossover, to generate
offspring
solutions
Selectively replace all the copies contained in the HM using the
generated offspring solutions.
ACCEPTED MANUSCRIPT It is worth noting that, the regeneration of the HM may introduce solutions that are much worse than the current best, as indicated by the spikes in the “worst value curve” (in red) shown in Figure 3. Anyhow, as indicated by the “current optimum curve” (in green), the regenerative process strongly improves the exploration of the solution space and, indeed, improvements often occur immediately after a regenerative process has taken place. ------------------------------Insert Figure 3 here
4.4.
A dynamic Large Portion Recovery Rate
As a last modification, instead of using a fixed
, we decided to dynamically update this threshold
limit, using the cubic function of Eq. (10). In this way, rapidly rises as
CR IP T
--------------------------------
remains low for small values of
increases. This behavior is rationale because when
, but it
is high almost all the solutions of
AN US
the HM are identical and only crossover can refresh the HM. For this reason, also
must approach one,
to activate almost certainly the Large Portion Recovery.
{
ED
Where:
is the value of
PT
is the minimum value of
CE
(10)
⁄(
;
below which
is the value of
)
M
[
above which
takes its minimum value; becomes 1;
).
AC
We conclude this section noting that all the above-mentioned improvements take advantage of the specific structure of the IJMS problem herein considered. Consequently, from here on, we will refer to this metaheuristic as the Harmony Search for Integrated Job Maintenance Scheduling, indicated as HS/4/IJMS. For further details the interested readers are referred to Appendix B, where the full pseudo-code can be found. 5.
CUCKOO SEARCH VIA LÉVY FLIGHT AND POSSIBLE ENHANCMENTS We also decided to investigate the performances of the Cuckoo Search via Lévy Flight (CS-LF), a
popular population-based metaheuristic introduced by Yang and Deb in 2009. CS-LF was originally designed to solve continuous problems, but its preformances are excellent also in case of discrete
ACCEPTED MANUSCRIPT combinatorial problems (see for example Ouaarab et al., 2014; Komaki et al., 2017). Also, and perhaps more important, since its search strategy assures little or none cross influence among saved solutions, CS-LF is a very good candidate to tackle the IJMS problem. 5.1.
Main Features of the standard CS-LF
CS-LF is a nature inspired metaheuristic, which mimics the brood parasitism of some species of cuckoo, a bird which does not brood its eggs but lays them in other nests. Each cuckoo explores its surroundings, by means of a Lévy Flight (LF), in search of a good nest to lay its egg. Evaluating the options, the cuckoo will
CR IP T
select the best one and it will lay its egg. Once an egg has been laid, it can be detected and purged by the owner of the nest, or it can successfully hatch. If so, the new cuckoo will continue the exploration of the surrounding area, looking for new nests to colonize.
Briefly, the key element of the CS-LF is that, at each step, a pairwise interchange scheme is iteratively applied to generate a neighbor
from an initial solution
, then a neighbor
and so on (i.e.,
, a parameter that represents the
AN US
) for, at most, a number of times equal to
from
distance travelled by the cuckoo. Also, this travel is modelled as a Lévy Flight that is, essentially, a random walk in which the length of the random step is drawn from a probability distribution that is heavy-tailed i.e., for which
, with
, so that the variance is infinite. For
, the Pareto
distribution (see Eq. (11)) is one of the distributions satisfying these requirements. (
M
{
)
(11)
ED
So, using a Pareto, in virtue of the inverse transformation techniques, the step length
is a random number uniformly distributed in the interval [
(12) ].
CE
where
PT
Eq. (12):
can be generated as in
Also, to model the possibility that some of the eggs are detected and smashed, a
of the worst solutions are
substituted with randomly generated ones. For the full algorithm, the interested readers are referred to Yang
AC
and Deb (2009); an extended and modified version is detailed in Sections 5.2 and 5.3. 5.2.
Cuckoo Search with Migration and Nests Location Memory (CS-MNM)
We also altered the searching strategy implemented by the CS-LF to explore the solution space in an even
better way. The key of such re-design is based on a modification of the standard Lévy Flight and on the introduction of three different features, namely: (i) Nests Location Memory, (ii) Risk Aversion and (iii) Clan Migration. 5.2.1
Double Shaped Lévy Flight
As shown in Figure 4, two flight modes are possible:
Mode #1 – A flight of progressive departure from the starting point
;
ACCEPTED MANUSCRIPT
Mode #2 – A peripheral flight centered on an intermediate point .
Note that, starting from
(highligthed in black in Figure 4), the overall flight (from 0 to 7) can be
decomposed into a series of intermediate points
(each one representing a feasible solution) where the
cuckoo pauses to assess the goodness of the surranding area. As long as the quality improve, with respect to the last visited point, the cuckoo continues to move away (as for points { , (as for points { ,
,
} and { ,
}). Otherwise
}), it explores the neighborhood of the last point of improvement (as for point
). The
latter case can be seen as the cuckoo coming back and then restarting its search or, simply, as the cuckoo
number of visited points
CR IP T
performing a sort of circular flight centered on the last point of improvement. Anyhow, the maximum is defined by the Lévy parameter of Eq. (12) and, as in the standard CS-LF, the
search is immediately interrupted if a point
improves the overall best solution found so far. In this case the
cuckoo stops its search and lays an egg. -------------------------------
AN US
Insert Figure 4 here -------------------------------Operatively, starting from
, all the visited points
are generated applying the pairwise interchange
operator to the last solution that produced an improvement. More precisely the first point neighbor of
; next, if
is better than
, provided that
M
is generated as a neighbor of
, the second point
is generated as a
is generated as a neighbor of
, otherwise
. This process is iteratively repeated until at most
point have been explored. For instance, in the case of Figure 4, the generating scheme is the following one:
5.2.2
ED
{
}.
Nest Location Memory and Risk Aversion (RA)
PT
Most of the times, the cuckoo starts a new flight from the last point of improvement found during its last flight (i.e.,
in the example of Figure 4). This always happens in case of a successful
CE
mission when a new host nest has been identified and the cuckoo has laid its egg. This also happens in case of an unsuccessful mission, provided that Mode #1 has been used at least once, during the Double Shaped Lévy Flight. Conversely, if only Mode #2 has been used, the cuckoo may try to exploit all the information
AC
obtained during its peripheral flight. If so, it might decide not to come back to the initial starting point rather to restart from the best intermediate point
, but
found during the flight, although this point was judged
inadequate to lay an egg and it is probably a bad move in the short-term. Clearly this choice carries a risk that is accepted by the cuckoo with a certain probability (namely the Non-Return Probability NRP). Also, the NRP is not fixed, but, as in Simulated Annealing, it is exponentially decreased over time, to better represent the risk aversion profile of the cuckoo, which is supposed to become more and more circumspect as it gets older. This is shown in Eq. (13).
ACCEPTED MANUSCRIPT
{ Where cuckoo
.
is the change in the objective function and
/}
, is the parameter that links the risk aversion of
with the number of performed iterations . To this aim
from an initial value
, until it reaches an ending value
is progressively and linearly reduced
at the last iteration , as shown in Eq. (14): )
(14)
CR IP T
( Concerning
(13)
, as suggested by Kirkpatrick et al. (1983), this value can be chosen to assure a high initial
Non-Return Probability
as in Eq. (15): ̅̅̅̅
(15)
Eq. (15) can be used to compute 5.2.3
AN US
where ̅̅̅̅ is the average change of the objective function over a sample of neighbor solutions. Also note that too, by substituting
Clan Migration (CM)
with a low final Non-Return Probability
.
We also introduced a migration strategy, to mimic the ability of birds to share part of their past experience. Specifically, since members of a clan tend to follow the behavior of the most successful
M
individuals, we hypothesized that, after a certain number of iterations (i.e., an epoch), some of the cuckoos that ended up in areas that offer low opportunities for improvement, migrate in the surrounding of the best
ED
and most promising areas so far discovered by them and or by other cuckoos. Operatively, this behavior has been implemented through the following steps:
Solutions are sorted in descending orders, in terms of the consecutive number of flights that did not
PT
produce any further improvement (i.e., solutions corresponding to cuckoos which are more likely to be stuck in a local minimum are placed at the top of the list); The first
cuckoos are preselected and, among these one, all those who exceeded a predefined
CE
threshold limit (expressed as number of consecutive failed flights) are finally selected for
AC
migration;
The starting points of migrating cuckoos are substituted with new points, generated as neighbors of the
5.3.
best solutions found so far.
CS-MNM full implementation
The implementation is based on two sets and its nest
and
storing, for each cuckoo , its current starting point
i.e., the best solution found so far by cuckoo . At each iteration, both sets are updated
accordingly to the outcome of the flight. Let us consider a generic cuckoo {
}, with
placed at point
. Also, let
, be the points visited by the cuckoo during its flight and let
be the best
among them. In accordance with the generation scheme of the set
(i.e., Double Shaped Lévy Flight),
ACCEPTED MANUSCRIPT anytime
proves to be a point of improvement, it replaces
following neighbor
is obtained. Also note that, as clarified in Sub Section 5.2.1, the flight is
immediately interrupted as soon as a coincides with both vectors
as the new starting point, from which the
improves the current best
and it is also better than
and
been unsuccessful,
. That is, whenever
. In this case the mission has been successful, and thus
are updated so that
. Conversely, if
is not updated and, typically, the cuckoo will go back to
changed during the flight (Mode #2 was employed exclusively), even if
, the mission has
. However, if
is worse than
has not
, the cuckoo might
and resume its next exploration from that point, with a certain probability (as explained in
Sub-Section 5.2.2). In this case
is updated so that
CR IP T
choose to go in
, the last point
.
For the sake of clarity, the full pseudocode of the CS-MNM is shown in Table IV. Briefly, the Double Shaped Lévy Flight is implemented from Line #14 to #24 and the part of code between Line #20 to #21 is responsible to switch between the two flight modes, as described in section 5.2.1. Cuckoo‟s Risk Aversion is
AN US
implemented by Line #25 to #28 and, as it can be seen, it is enabled anytime the Double Shaped Lévy Flight proved unsuccessful. Lastly, Clan Migration is detailed by Line #29 to #32. Table IV CS-MNM pseudocode
PT
ED
M
Set , Set P% Set , Set N Set Randomly generate the set N of host nests EN u 0N k0 Do While Or (Stop Criterion) Randomly select Cuckoo i: set , j0
14
Do While
CE
13
15
Risk Aversio n
16 17 18 19 20 21 22 23 24 25
From generate new neighbor If Then Break If Then
jj+1 Loop If
# Exploration set, impose initial condition E ≡ N # Vector of consecutive unsuccessful flights
accordingly
Accept
and store it in set # Improvement over # Replace both starting point and last host nest # Reset count of unsuccessful flights # Halt Double Shaped Lévy Flight # Improvement over # Replace only the starting point # Increase count of unsuccessful flights
# No improvements over ( )
26 27
# Levy Flight lower limit and Power Law Exponent # Epoch migration population ratio # Initial and Final Non-Return Probability # Number of iterations # Threshold of consecutive unsuccessful flights
Generate travel distance d
AC
Double Shaped Lévy Flight
1 2 3 4 6 7 8 5 9 10 11 12
accordingly to Eq. (13)
# Best of generated neighbors
Clan Migration
ACCEPTED MANUSCRIPT 28 29 30 31 32 33
AC
CE
PT
ED
M
AN US
CR IP T
34
If Accepted ( Then: # Accept a worsening exploration starting point If An Epoch has Elapsed Then B P% of elements of N, sorted in ascending order accordingly to W P% of elements of E, sorted in descending order accordingly to , s.t. Replace each vector solution of W with a neighbor of a vector solution of B kk+1 Loop
ACCEPTED MANUSCRIPT 6. PERFORMED TESTS All the tests discussed in this Section were performed on a test set of randomly generated problems characterized by the following features (Zammori et al., 2014):
Processing Times are Normally distributed with mean
and standard deviation
of 100 and 25
time-units, respectively;
Setup Times are Uniformly distributed in the interval [
Due Dates are Uniformly distributed in the interval [
], where
and
is the Tardiness Factor, representing the quota of jobs likely to be tardy,
is the relative Range of Due Dates (Tan et al. 2000);
The earliness and tardiness penalties have been set, respectively, at
Time-to-Failure follows a Weibull distribution with scale parameter ;
and
;
and shape parameter
The average time to execute Planned and Reactive Maintenance are set, respectively, to and
AN US
is the
CR IP T
number of jobs,
];
.
6.1. Parameters optimization
To assure a fair comparison of the metaheuristics introduced in Sections 3 to 5, their operating parameters were fine-tuned using a test set of 50 small size problems, characterized by 14-Jobs and 2-PMs. Problems of
M
this size were chosen because: (i) their complexity is high enough to prevent the metaheuristics from finding, almost always, the optimal solution and, (ii) they can be optimally solved in a reasonable amount of time
ED
(i.e., around one and a half hour on a computer with a 3.2 GHz processor and 16 GB of RAM memory), through an exhaustive search. Also, parameters that can be reasonably considered independent were finetuned one by one, using a sequential optimization strategy. Conversely, parameters for which interaction
More precisely:
MHS was not considered in the fine-tuned process, as we used the parameters (reported in Table V)
CE
PT
effects cannot be excluded a priori, were pre-analyzed and optimized using Design-Of-Experiment (DOE).
suggested by Zammori et al. (2014), which had been already optimized on the same problem set; SA and CS-LF (i.e., the original metaheuristics) were optimized first. Parameters were assumed
AC
independent and fine-tuned following the sequential approach;
HS/4/IJMS and CS-MNM were optimized next. In this case, parameters shared with the original metaheuristics were directly inherited from the previous optimization step; conversely, specific parameters were optimized using DOE, to better investigate possible interactions among them.
6.1.1. Sequential optimization Let us consider a generic algorithm , and let
{
} be the set of its parameters, ordered
from the most to the least important one. Sequential optimization proceeds as follows:
Initialize each parameter pi using a reasonable starting value. To this aim we used values suggested in the literature;
ACCEPTED MANUSCRIPT
Given these initial settings, all problems of the test set are solved and the accuracy of
is
computed. To this aim, we measured accuracy in terms of APE, as in Eq. (8), using 8 repetitions of each problem instance;
The parameter
, having the higher influence on the algorithm performance, is selected and its
value is slightly modified. All the other parameters {
} are kept unaltered and all the
problems of the test set are solved again. This process is iteratively repeated, proceeding at small step, until the accuracy of
continues to grow; is frozen and the second parameter
is
CR IP T
When accuracy does not increase anymore, the value of
selected for possible improvements and it is fine-tuned using the same approach described above;
The procedure is repeated for all the other parameters of the set .
Obtained results are shown in Table V, where optimized parameters are listed in order of importance. The optimal value is shown in the last column, whereas the range and the step size used in the optimization are
AN US
reported in the third and fourth column, respectively. Table V Sequential optimization results Algorithm
Parameter Random Selection Rate Pitch Adjustment Rate Harmony Memory Size Start Probability End Probability Minimum Flight dmin Pareto Shape parameter Purge percentage P% Population Size
Step Size
Optimal Value
[0.05; 0.50] [0.001; 0.03] [1; 10] [1; 2] [0.01; 0.25] [0.5NJ; 2NJ]
0.05 0.005 1 0.3 0.05 0.25
0.05 0.1 0.5NJ 0.1 0.02 7 1.9 0.1 NJ
M
MHS
Range
ED
SA
CS-LF
PT
Where NJ is the total number of Jobs
6.1.2. Design of Experiment
CE
Given the higher complexity of HS/4/IJMS and CS-MNM, a 2-levels full factorial design was made (using the levels shows in Table VI), limited to the following parameters: Concerning HS/4/IJMS, we considered
AC
,
and
, because these parameters
conjointly govern the rate at which the Large Portion Recovery is triggered;
Concerning CS-MNM we considered
,
, and the Epoch‟s Length, because these parameters
conjointly influence the evolution of the solutions contained within sets
and .
More precisely, the 10 hardest problem instances (among the 50 previously generated ones) were selected and each one was solved 8 times, both with HS/4/IJMS and with CS-MNM. Testing
combinations of
parameters levels, a total of 640 solutions (i.e., 10·8·23) were generated with each metaheuristic. It is worth noting that, to assess interaction effects through 3-ways ANOVA, there is the need to collect more than a
ACCEPTED MANUSCRIPT single value of the response variable for each one of the 10·23 combinations of problems and parameters setting. Thus, we replaced APE with the Relative Percentage Error defined in Eq. (16): (16) Where, attempt
[
is the (a priori known) optimal solution of problem [
] and
] by the used metaheuristic, with parameters configuration
[
is the solution found at ].
Operating in this way, two remarkable results were found, as DOE and the related ANOVA tests revealed
CR IP T
that:
Concerning possible interaction effects, all p-values were higher than 0.25 (definitely greater than the usual threshold of 0.05), and so parameters can be considered independent;
Concerning the effect of the operating parameters on the algorithm performance, all p-values were higher than 0.1, and so performances are stable and little influenced by specific values of the
AN US
operating parameters.
Since both results confirmed the robustness of the proposed metaheuristics, as well as the absence of any interactions among parameters, we completed the optimization process following the same sequential approach detailed above. Obtained results are shown in Table VI.
CS-MNM
Levels
Optimal Value
Low saturation threshold SLow
[0.1; 0.3]
0.3
High saturation threshold SHigh Large Portion Recovery Rate LPRRLow Starting Risk Aversion SRA Ending Risk Aversion ERA Epoch E
[0.5; 0.7]
0.7
[0.05; 0.1] [0.1; 0.5] [0.01; 0.05] [50; 100]
0.05 0.1 0.01 75
PT
HS/4/IJMS
Parameter
ED
Algorithm
M
Table VI DOE results
CE
6.2. Numerical Assessment
Lastly, SA, CS-LF, HS/4/IJMS and CS-MNM were compared in relative terms. Please note that the
AC
original MHS was not included in this analysis, because, as detailed in Section 3.6, it is outperformed by the SA, even with problems of relatively small dimensions. We also note that, to get a preliminary idea of the quality of the modifications made on the MHS, we tested its improved version (i.e., the HS/4/IJMS) against the same test set of 10 problems with 14-Jobs and 2-PMs. In this case, results were encouraging as HS/4/IJMS behaved almost like the SA and it consistently found the global optimum (i.e., APE = 0 and APD = 0.01). Comparisons were made relatively to two test sets with a different level of complexity: (i) 10 hard problems with 25-Jobs and 3-PMs, and (ii) 10 very hard problems with 50-Jobs and 8-PMs. In both cases, the number of PMs coincides with the optimal value, calculated accordingly to the procedure proposed by Zammori et al. (2014), to which the interested reader is referred for further insights.
ACCEPTED MANUSCRIPT Each problem instance was solved 8 times, by each metaheuristic, starting from a different initial population. Also, we let each metaheuristic run until convergence, so that computation was halted when the best solution did not improve for a sufficiently high number of consecutive iterations. Lastly, performances were assessed using APE and APD, as defined in Eq. (8) and Eq. (9), respectively. Please note that, in this case, the real optimum were not obtained through an exhaustive search; thus, in both equation the parameter substituted by the best solution found for problem , i.e.
{
is
}.
Table VII Final comparison 25-Jobs 3-PMs
CR IP T
Obtained results are reported in Table VII.
50-Jobs 8-PMs
APE
APD
Max Error
# Wins
APE
APD
HS/4/IJMS
0.007
0.023
0.0491
0%
0.0580
0.2141
0.3562
0%
SA
0.001
0.011
0.0400
20%
0.0030
0.0213
0.0647
40%
CS-LF
0.087
0.104
0.1416
0%
[NC]
[NC]
[NC]
[NC]
CS-MNM
0
0.004
0.0109
40%
0.0008
0.0074
0.0187
60%
AN US
Algorithm
Max Error
# Wins
As expected, the SA works very well and, even with 25-Jobs and 3-PMs it outperforms both CS-LF and HS/4/IJMS. Specifically, whereas performance of HS/4/IJMS are still acceptable, in its standard form the
M
Cuckoo Search is definitely the worst metaheuristics, both in terms of solutions‟ quality and stability. For these reasons CS-LF was not even considered in the second test set (as denoted by the [NC] values in
ED
columns 6 to 9 of Table VII). Frankly speaking, the very poor performance of the CS-LF was not expected. Probably, although solutions evolve independently and cross influence among them is minimized, the solution space is explored neither efficiently nor effectively. Alone, the Lévy Flight does not seem to be
PT
enough to assure the right random behavior required for a thorough exploration of the solution space. Also, purging the worst solutions (at every iteration) may limit the neighborhood search, thus further decreasing
CE
the overall effectiveness of the approach. Conversely, in its improved version, the CS-MNM turned out to be the best approach. CS-MNM, indeed, shows the minimum values of APE, APD and Max Error; particularly, a value of APE equal to zero ascertains that, for each problem instance, at least one of the eight solutions
AC
generated by the CS-MNM was the overall best one (relatively to the 32 generated by all the tested metaheuristics). CS-MNM also proved to be very precise in founding the best solutions as certified by a very low value of both APD and Max Error and by the highest value of percentage wins. Relatively to this last metric, it is important to note that summing the percentage wins of column 5 a total of 60% is obtained. This means that 40% of the times SA and CS-MNM converged to the same, presumably optimal, solution. The superiority of the CS-MNM is even better outlined by the outcomes of the very hard problems‟ set. Not only CS-MNM is the metaheuristic that located the best solution with the highest frequency but, as certified by the lowest APD value, it proved to be extremely robust because, even when it missed the best solution, it still found a solution very close to it. All these outcomes indicate that the introduced modifications (i.e., nests location memory, cuckoo‟s risk aversion, clan migration and double shaped Lévy
ACCEPTED MANUSCRIPT Flight), made it possible to overcome the drawbacks of the standard approach, assuring an efficient and effective exploration of the solution space. As a final note, we want to highlight the following fact. Since convergence was used as stop criterion, one could argue that comparison between CS-MNM and SA was not fair. Indeed, being a population-based metaheuristic, CS-MNM might have found better solutions only because it ran for a total number of iterations higher than that of the SA. More precisely, given a fixed number of iterations , SA evolves the same solution exactly
of the set , at iteration each solution will have been evolved
times, on average. It is thus
CR IP T
among the
times; conversely, considering that at each iteration CS-MNM selects a solution
plausible to expect that, to find the best solution, CS-MNM needs to run for a number of iterations
times
bigger than that of the SA. To ascertain this fact, we compared the evolution path of SA with that of CSMNM, for each generated solution. Quite surprisingly, the test showed that 75% of times, when SA reached convergence, CS-MNM had already found a better solution and its evolution process was still on progress; a fact that further certifies the superiority of the CS-MNM, which is also faster than the SA in finding the
AN US
optimal solution. This result is graphically displayed by Figure 5.a and Figure 5.b that show an example of the evolution path followed by the two metaheuristics.
------------------------------Insert Figure 5.a here
--------------------------------
M
------------------------------Insert Figure 5.b here
6.3.
ED
--------------------------------
Operational Implications
PT
Once the superiority of the CSM-MNM algorithm was verified, a last check was made to quantify the operational advantage that could be achieved using the CSM-MNM algorithm to schedule jobs and PMs in
CE
an integrated way. To this aim we compared the cost of the best schedule found by the CSM-MNM algorithm with the cost of the presumably best schedule that could be generated by a production manager
AC
equipped with a standard scheduling system. The latter one was generated using the following naïve algorithm:
At first only jobs are considered (i.e., the single machine is considered as fault free) and the best schedule S is generated. This schedule, that using Kendall‟s notation can be classified as (1| | Earliness-Tardiness Penalties), corresponds to the one that could be generated by the standard scheduling system and, in the present case is obtained evaluating all the (n!) permutations of the jobs‟ sequence. Next, the optimal number M of PMs is inserted in S, without altering the ordering of the jobs. This is made in the following reasonable way: - S is split into (M +1) parts of almost equal duration;
ACCEPTED MANUSCRIPT - PMs are inserted at these subdivision points, to assure the same level of availability in each subsection of the schedule. The comparison has been made testing the exhaustive solver and the CS-MNM against 5 different instances of the IJMS problem characterized by 50-Jobs and 8-PMs. The CS-MNM always outperformed the naïve exhaustive, and solutions generated by the CS-MNM were on average (23 ± 1.1) % more efficient, in terms of expected operating cost, than that obtained by the naïve exhaustive. These results stress even more the
CR IP T
relevancy and applicability of the proposed approach.
7. CONCLUSIONS AND FUTURE WORKS
In this paper, we focused on an integrated jobs-maintenance scheduling problem, relevant to many industrial realities, where failures and/or reduced availability may jeopardize productivity. Specifically, we considered a single machine subjected to failures and aging with the objective to minimize earliness-
AN US
tardiness penalties. These features make the problem particularly challenging, because the solution space becomes very flat, but disseminated of many local minima.
Considering the practical relevance of the problem and its inherent complexity, we developed and tested different metaheuristics, aiming to find robust and quasi-optimal schedules in a reasonable amount of time. Specifically, we benchmarked the famous SA by Ishibuchi et al. (1995), with two population-based
M
metaheuristics, namely the Modified Harmony Search by Zammori et al. (2014) and the Cuckoo Search by Lévy Flight by Yang and Deb (2009). MHS was chosen because it was specifically designed to tackle this problem, CS-LF was chosen as it implements a very effective local minimum escaping strategy. In addition,
ED
we also developed two alternative approaches: (i) HS/4/IJMS, a hybrid Harmony Search and Genetic Algorithm, specifically tied to this problem and (ii) CS-MNM, an innovating modification of the Cuckoo
PT
Search via Lévy Flight.
After a fine-tuning procedure, a comprehensive comparison was made. Outcomes showed that, in these
CE
peculiar conditions, simpler heuristics characterized by a high degree of randomness may perform better than sophisticated population-based metaheuristics. Indeed SA (i.e., the simpler algorithm) outperformed MHS, CS-LF and HS/4/IJMS, especially in case of very hard problems, characterized by 50-Jobs and 8 PMs. Only
AC
the newly developed CS-MNM subverted this situation, as it was the best approach both in term of precision and accuracy of the obtained solutions. We also showed that CS-MNS offers the following operating advantages:
It has relatively few parameters to optimize; Parameters are independent (there are no interaction effects among them) and can be optimized independently; It performs an excellent exploration of the solution space but, at the same time, it assures a fast convergence and a low computation time, also for large problems.
All these features make the CS-MNS algorithm a good candidate for a possible inclusion in the optimizer module of a Cyber Physical System, fully compliant with the principles of the Industry 4.0 paradigm.
ACCEPTED MANUSCRIPT Frankly speaking we must admit that, at present, our model requires the number of PMs that must be inserted in the schedule as input; although this value can be analytically estimated, a better approach would be to use a two-level optimization to generate, not only the optimal job-maintenance sequence, but also the optimal number of PMs. Another limit of the model is that a single level of planned maintenance is considered because, anytime PMs are executed, the hazard rate of the machine is brought back to its original level (i.e., perfect repair model). However, a more realistic scenario could include different levels of maintenance tasks, each one characterized by different costs and times and capable to reduce (without
CR IP T
necessarily reset) the aging and/or the degradation of the machine. Lastly, we note that, although the single machine scheduling problem is a rather flexible model, as it can be applied to a variety of common production systems, it could also be extended to other common layouts, such as parallel machines, flow shops or even job shops. All these issues could be topics for future researches.
AN US
Authors contributions section
F.Zammori and D.Mezzogori conceived of presented idea and development of the CS-MNM algorithm. M.Bertolini implemented the rest of the algorithms used as benchmark. D.Mezzogori executed the Design of Experiment section and F.Zammori performed final numerical comparison. M.Bertolini finally supervised
M
and checked the findings of this work. All authors discussed the results and contributed to the final
ED
manuscript.
PT
8. REFERENCES
AC
CE
Allahverdi, A. (2015). The third comprehensive survey on scheduling problems with setup times/costs, European Journal of Operational Research, 377 (2), 345-378. El-Khoukhi, F., Boukachour, J., & Alaoui, A.L.H. (2017). The „„Dual-Ants Colony”: A novel hybrid approach for the flexible job shop scheduling problem with preventive maintenance, Computers & Industrial Engineering, 106, 236– 255. Cui, W., Lu, Z., & Pan, E. (2014). Integrated production scheduling and maintenance policy for robustness in a single machine, Computers & Operations Research, 47, 81-91. Cui, W., & Lu, Z. (2017). Minimizing the makespan on a single machine with flexible maintenances and jobs‟ release dates, Computers & Operations Research, 80, 11-22. Cui, W., Lu, Z., Li, C., & Han, X. (2018). A proactive approach to solve integrated production scheduling and maintenance planning problem in flow shops, Computers & Industrial Engineering, 115, 342–353. Geem, Z., Kim, J., & Loganathan, G. (2001). A new heuristic optimization algorithm: harmony search, Simulation, 76 (2), 60-68. Gen, M., & Cheng, R. (2000). Genetic algorithms and engineering optimization. John Wiley & Sons, New York, USA. Ishibuchi, H., Misaki, S., & Tanaka, H. (1995). Modified simulated annealing algorithms for the flow shop sequencing problem, European Journal of Operational Research, 81 (2), 388–398. Kirkpatrick, S., Gelatt, C., & Vecchi, M. (1983). Optimization by simulated annealing, Science, 220 (4598), 671-680. Komaki, G., Teymourian, E., Kayvanfar, V., & Booyavi, Z. (2017). Improved discrete cuckoo optimization algorithm for the three-stage assembly flowshop scheduling problem, Computers & Industrial Engineering, 105, 158–173. Lee, J., Bagheri, B., & Kao, H. A. (2015). A cyber-physical systems architecture for industry 4.0-based manufacturing systems. Manufacturing Letters, 3, 18-23. Lee W.C., & Wang J.Y. (2017). A three-agent scheduling problem for minimizing the makespan on a single machine. Computers & Industrial Engineering, 106, 147-160.
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
AN US
CR IP T
Low C., Hsu C., & Su C., (2010) A modified particle swarm optimization algorithm for a single-machine scheduling problem with periodic maintenance. Expert System with Applications, 37, 6429–6434. Lu Z., Cui, W., & Han X., (2015). Integrated production and preventive maintenance scheduling for a single machine with failure uncertainty. Computers & Industrial Engineering, 80, 236–244. Luo, W., Cheng, T.C.E., Ji, M., (2015). Single-machine scheduling with a variable maintenance activity. Computers & Industrial Engineering, 79, 168 – 174. Macchi M., Roda I., Fumagalli L., 2017, 'On the Advancement of Maintenance Management Towards Smart Maintenance in Manufacturing', IFIP Advances in Information and Communication Technology, 383-390 Metropolis, N., Rosenbluth, A., Rosenbluth, M., Teller, A., & Teller, E. (1953). Equation of State Calculations by Fast Computing Machines. The Journal of Chemical Physics, 21 (6), 1087-1092. Mokhtari, H., Dadgar, M., (2015). Scheduling optimization of a stochastic flexible job-shop system with time-varying machine failure rate, Computers & Operations Research, 61, 31-45. Mor, B., & Mosheiov, G., (2015). Scheduling a deteriorating maintenance activity and due-window assignment. Computers & Operations Research, 57, 33-40. Mosheiov, G., & Sarig, A., (2009). Scheduling a maintenance activity and due-window assignment on a single machine. Computers & Operations Research, 36, 2541-2545 Nouiri M., Bekrar, A., Jemai, A., Trentesaux, D., Ammari, C.A., & Niar, S., (2017). Two stage particle swarm optimization to solve the flexible job shop predictive scheduling problem considering possible machine breakdowns, Computers & Industrial Engineering , 112, 595-606. Ouaarab, A., Ahiod, B., & Yang, X.Y., (2014). Discrete cuckoo search algorithm for the travelling salesman problem. Neural Computing & Application, 24, 1659–1669 Quan, O. Y., & Xu, H. Y. (2013). The review of the single machine scheduling problem and its solving methods. In Applied Mechanics and Materials (Vol. 411, pp. 2081-2084). Trans Tech Publications. Rahmati, S.H.A., Ahmadi, A., & Govindan, K., (2017). A novel integrated condition-based maintenance and stochastic flexible job shop scheduling problem: simulation-based optimization approach, Annals of Operations Research, 139. Safari, E., & Sadjadi S., (2011). A hybrid method for flowshops scheduling with condition-based maintenance constraint and machines breakdown. Expert System with Applications, 38, 2020-2029. Salmasnia, A., & Mirabadi-Dastjerd, D., (2017). Joint production and preventive maintenance scheduling for a single degraded machine by considering machine failures. TOP, 25, 544–578. Soltani, R., Jolai, F., & Zandieh, M. (2010). Two robust meta-heuristics for scheduling multiple job classes on a single machine with multiple criteria. Expert Systems with Applications, 37 (8), 5951–5959. Tan, K. C. T., Narasimhan, R., Rubin, R. A., Ragatz, & G. L., (2000). A comparison of four methods for minimizing total tardiness on a single processor with sequence dependent setup times, Omega, 28, 313-326. Touat, M., Bouzidi-Hassini, S., Benbouzid-Sitayeb, F., & Benhamou, B. (2017). A hybridization of genetic algorithms and fuzzy logic for the single-machine scheduling with flexible maintenance problem under human resource constraints, Applied Soft Computing, 59, 556-573. Woo, Y., Kim, B., (2018). Methaeuristic approaches for parallel machine scheduling problem with time-dependent deterioration and multiple rate-modifying activities, Computers & Operations Research, 95, 97-112. Yang, X.-S. (2010). Nature-inspired metaheuristic algorithms 2nd Ed., Luniver Press, Frome, UK. Yang, X.-S., & Deb, S. (2009). Cuckoo search via Lévy flights, Proceeding of World Congress on Nature & Biologically Inspired Computing, December 2009, India, IEEE Publications, USA, 210-214. Yazdani, M., Khalili, S.M., Babagolzadeh,M., & Jolai, F., (2017). A single machine scheduling problem with multiple unavailability constraints: A mathematical model and an enhanced variable neighborhood search approach, Journal of Computational Design and Engineering, 4, 46–59. Yu, A.J., & Seif J.,(2016). Minimizing tardiness and maintenance costs in flow shop scheduling by a lower-boundbased GA. Computers & Industrial Engineering, 97, 26 – 40. Zammori, F., Braglia, M., & Castellano, D. (2014). Harmony search algorithm for single-machine scheduling problem with planned maintenance, Computers & Industrial Engineering, 76, 333–346. Zhu, H., Li, M., Zhou, J. (2015). Machine scheduling with deteriorating and resource-dependent maintenance activity, Computers & Industrial Engineering, 88, 479-486.
ACCEPTED MANUSCRIPT Appendix A. Simulated Annealing (SA) Introduced by Metropolis et al. (1953) and made popular by Kirkpatrick et al. (1983), SA is a probabilistic heuristic, whose inspiration come from annealing in metallurgy, a slow cooling process used to bring melted metals to a uniform, low-energy crystallization state. Briefly, at each one of neighborhood of the current solution
iterations, SA generates a new solution
from the
and accepts it as the new starting point, with a probability that, as iterations go
on, decreases following the negative exponential law show in Eq. (a.1), valid for a minimization problem.
Where
(
)}
(a.1)
CR IP T
{
is the change in the objective function and
, the so-called temperature at iteration , is a
parameter that accounts for the number of completed iterations. To this aim, the temperature is progressively and linearly reduced from an initial value
, as in Eq. (a.2): (
)
AN US
For the sake of clarity, the full SA Pseudo-Code is reported in Table A.I
(a.2)
Table A.I SA by Ishibuchi et al. (1995)
# Number of Iterations, Initial and Ending Temp. # Number of Solutions explored for each level of T
# Current and Best Solution # Best and Current Objective Function
ED
M
Set N, T0, TN Set K Randomly generate an initial solution y x, x* y x, x* f(y) i0 Do While T0 < TN Generate K new neighbors of x {y1,…, yK}
# Best OF, among the newly generated solutions # Best of the K generated solutions
PT
y Accept accordingly to Eq. (a.1) If Accepted ( Then: x y, x y If y < x* Then: x* y, x* y ii+K Update Ti as in Eq. (a.2) Loop Return x*, x*
# Temperature is linearly reduced
CE
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
AC
Also, neighbor solutions are generated using the pairwise interchange operator i.e., two locations and are randomly selected and the corresponding elements
and
are swapped:
and
Concerning the initial and ending temperatures, their values are set accordingly to the classical approach by Kirkpatrick et al. (1983). Specifically,
is chosen to assure a desired (high) Initial Acceptance Probability (IAP) computed as in
Eq. (a.3), where ̅̅̅̅ is the average variation in the objective function when a neighbor is generated. ̅̅̅̅
Please note that Eq. (a.3) can be used to compute Acceptance Probability (EAP).
(a.3)
too, by simply substituting IAP with the desired (low) Ending
ACCEPTED MANUSCRIPT
# Novelty
AN US
CR IP T
Set # Harmony Memory Consideration Rate Set # Pitch Adjustment Rate Set # Block Shifting Rate Set # Number of iterations Set # Number of job and maintenance tasks Randomly generate the Do While Or (Stop Criterion) Compute accordingly to Eq. (10) If Then Generate empty solution For To Do If Then: Get an element via Random Selection (RS) Else Get an element via Harmony Memory Consideration (HMC) # Add element to candidate solution Next If Then: Perform Pitch Adjustment (PA) on Else If Then: Perform Block Shifting (BS) on
# Novelty
# Novelty
ED
M
If Then: # Overwrite worst solution Else Let be the number of copies within unique solutions within Initialize set of offsprings Do While : Select two solutions from using Wheel Selection Rule Perform Order Crossover (OC) on the selected solutions Store generated solutions into For Each Offspring In Do If Then Perform Pitch Adjustment (PA) on the Offspring Replace a copy within the with the Offspring Loop
PT
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34
# Novelty
CE
Large Portion Recovery
HM Consideration
APPENDIX B. HS/4/IJMS Pseudo-Code
All modifications (made to the original MHS algorithm) are clearly defined in the HS/4/IJMS pseudocode.
AC
Specifically:
Line #9 is used to dynamically compute the value of LPRR, based on the saturation of the Harmony Memory, which will trigger either the original Harmony Memory Consideration or the Large Portion Recovery;
Line #18 introduces the Block Shifting operator;
Lines from #21 to #32 specify how the Large Portion Recovery works. Specifically, lines #29 to #32aim to minimize the saturation of the Harmony Memory, by replacing all repeated solutions.
Note that standard HS corresponds to lines #10 and line #20, except for line #18.
ACCEPTED MANUSCRIPT Fig. 1.a – Neighbour solutions of three instances of a 6-Jobs 1-Planned Maintenance problem Fig. 1.b – Sorted solutions of three instances of a 6-Jobs, 1-Maintenance problem Fig. 2 – Typical trend of the Saturation Index during a computational run Fig. 3. – Typical trend of the objective function during a computational run Fig. 4. – Cuckoo‟s flight modes Fig. 5.a – Typical evolution in the initial phases Fig. 5.b – Typical evolution before convergence is reached
AC
CE
PT
ED
M
AN US
CR IP T
1. 2. 3. 4. 5. 6. 7.
AN US
CR IP T
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
Fig. 1.a – Neighbour solutions of three instances of a 6-Jobs 1-Planned Maintenance problem
CR IP T
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
AN US
Fig. 1.b – Sorted solutions of three instances of a 6-Jobs, 1-Maintenance problem
CR IP T
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
AN US
Fig. 2 – Typical trend of the Saturation Index during a computational run
CR IP T
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
AN US
Fig. 3. – Typical trend of the objective function during a computational run
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
AN US
CR IP T
Fig. 4. – Cuckoo‟s flight modes
CR IP T
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
AN US
Fig. 5.a – Typical evolution in the initial phases
CR IP T
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
AN US
Fig. 5.b – Typical evolution before convergence is reached
Comparison of new metaheuristics, for the solution of an integrated jobs-maintenance scheduling problem Massimo Bertolini , Davide Mezzogori , Francesco Zammori PII: DOI: Reference:
S0957-4174(18)30799-1 https://doi.org/10.1016/j.eswa.2018.12.034 ESWA 12375
To appear in:
Expert Systems With Applications
Received date: Revised date: Accepted date:
17 July 2018 8 November 2018 18 December 2018
Please cite this article as: Massimo Bertolini , Davide Mezzogori , Francesco Zammori , Comparison of new metaheuristics, for the solution of an integrated jobs-maintenance scheduling problem, Expert Systems With Applications (2018), doi: https://doi.org/10.1016/j.eswa.2018.12.034
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
HIGHLIGHTS The paper focuses on a single machine Jobs-Maintenance scheduling problem.
We enhanced the MHS by Zammori et. al (2014).
We also design a Cuckoo Search (CS) metaheuristics based on a new Lévy Flight.
Performance of the MHS improved, but this metaheuristic remains sub-optimal.
The newly developed CS is, instead, highly effective and robust.
AC
CE
PT
ED
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
Comparison of new metaheuristics, for the solution of an integrated jobs-maintenance scheduling problem Massimo Bertolini, Davide Mezzogori* and Francesco Zammori
CR IP T
Department of Engineering and Architecture, University of Parma, Italy Parco Aree delle Scienze, 181/A, 43124, Parma
Co-Author Prof. Massimo Bertolini
Phone: +39 0521 905705
*Corresponding Author Dr. Davide Mezzogori Email: [email protected]
Prof. Francesco Zammori
ED
Co-Author
M
Phone: +39 0521 905887
PT
Email: [email protected]
AC
CE
Phone: +39 0521 905887
AN US
Email: [email protected]
ACCEPTED MANUSCRIPT
Comparison of new metaheuristics, for the solution of an integrated jobs-maintenance scheduling problem ABSTRACT This paper presents and compares new metaheuristics to solve an integrated jobs-maintenance scheduling problem, on a single machine subjected to aging and failures. The problem, introduced by Zammori et al. (2014), was originally solved using the Modified Harmony Search (MHS) metaheuristic. However, an
CR IP T
extensive numerical analysis brought to light some structural limits of the MHS, as the analysis revealed that the MHS is outperformed by the simpler Simulated Annealing by Ishibuchi et al. (1995). Aiming to solve the problem in a more effective way, we integrated the MHS with local minima escaping procedures and we also developed a new Cuckoo Search metaheuristic, based on an innovative Lévy Flight. A thorough comparison confirmed the superiority of the newly developed Cuckoo Search, which is capable to find better solutions in
AN US
a smaller amount of time. This an important result, both for academics and practitioners, since the integrated job-maintenance scheduling problem has a high operational relevance, but it is known to be extremely hard to be solved, especially in a reasonable amount of time. Also, the developed Cuckoo Search has been designed in an extremely flexible way and it can be easily readapted and applied to a wide range of combinatorial problems.
M
Keywords - Cuckoo Search, Harmony Search, Metaheuristics, Scheduling, Single Machine, Planned Maintenance.
ED
1. INTRODUCTION
Scheduling is an important tool for manufacturers and production engineers, as it deals with the way in which human and machinery resources should be allocated to meet customers‟ demand. Specifically, given a
PT
set of production orders (or jobs), scheduling allows one to decide when and in which order jobs should be processed, with which resources and on which equipment. The objective is to control and to optimize
CE
workloads, to minimize production time and costs, respecting all production constraints (such as jobs‟ routing and jobs‟ releasing time) and, possibly, delivery due dates. Therefore, an optimized schedule may have a major impact on the performance of a manufacturing system and can be source of competitive
AC
advantage, as it allows reducing work-in-process, queuing times and industrial costs. It is worth noting that scheduling problems have always been of great interest also for the academia because, being very hard to be solved, they offer a challenging test bench to assess performance of sophisticated algorithms and/or metaheuristics (Allahverdi, 2015). As a matter of fact, scheduling models proposed in technical literature are becoming more and more complex, since current research‟s trends are toward integrating new and practical constraints in traditional scheduling models (El Khoukhi et al., 2017). In this regard, an issue of recent interest concerns the possibility to increase machines‟ availability and to reduce the overall production costs, by including Planned Maintenance tasks (PMs) into the jobs‟ schedule (Salmasnia et al., 2017). Traditionally, in fact, machines and equipment have been considered as fault free and the focus was limited to the way in which jobs should be sequenced and allotted to the available resources. Yet, this assumption is
ACCEPTED MANUSCRIPT rather heroic because machines can fail and, due to failures and reactive maintenance tasks (needed to restore the operating state of the fault machines) processing and waiting times gets longer, routings could need to be changed accordingly and, consequently, the original schedule may become non-optimal and should be reevaluated. Also, and perhaps more important, since breakdowns reduce productivity and increase costs, PMs should also be scheduled, as additional tasks, as an effective way to improve machines‟ availability (Safari and Sadjadi, 2011). The opportunity to define the scheduling sequence, taking into account possible failures and the need to perform PMs is becoming even more engaging, now that we are entering in the so-called
CR IP T
forth industrial revolution (Macchi et al., 2017). Indeed, due to the recent advances in manufacturing industry, which is leading toward the systematical deployment of Cyber-Physical Systems, it is now possible not only to remotely monitoring the system‟s state and to take actions based on real time data, but also to synchronize the physical factory floor and the cyber computational space (Lee et al., 2014), where optimization algorithms can be effectively used to dynamically solve the scheduling problem, in an acceptable amount of time.
AN US
Among the papers dealing with the Integrated Job Maintenance Scheduling problem (IJMS), this work focuses on the one by Zammori et al. (2014), who faced the problem of scheduling n jobs on a single machine subjected to aging and failures and proposed a Modified-Harmony-Search (MHS) metaheuristic to generate a quasi-optimal solution, in terms of earliness-tardiness penalties. The rational aspect of this choice is twofold:
The underlying hypotheses of the model are very pertinent with the industrial reality, as the single
M
machine is characterized by a time dependent (and increasing) failure rate and both reactive and
ED
planned maintenance tasks are executed to restore the functionality of the system. The model has a high practical relevance: although the single machine is the simplest form of shop
2013).
PT
floor configuration, it is also the one with the largest number of practical applications (Quan and Xu,
Concerning the last point, we note that the single machine can be used to adequately represent traditional
CE
continuous process plants (such as refining), but it can be easily readapted to model batch processes such as those of bakeries, biotechnology, chemical and pharmaceutical firms. Not only these production systems are
AC
very common, but it is exactly in these production environments that scheduling decisions are extremely important (to boost productivity) and where failures may have a dramatic impact, as they could interrupt the whole manufacturing process. Additionally, the single machine can also be used to schedule production in a pure flow shop or even in a job shop system, provided that the single machine is used to represent the real bottleneck of the manufacturing system. We also note that our interest (in this model) is mostly operational and computational and the aim is to introduce alternative solution approaches, more performing than the original MHS. The need to find better solutions was urged by the results of the analysis made by Zammori et al. (2014), who comparted the Modified Harmony Search (MHS) to the Hybrid Kangaroo Simulated Annealing by Soltani et al. (2010), a metaheuristic known to be particularly effective with NP-hard scheduling problems. Although MHS resulted
ACCEPTED MANUSCRIPT as the best approach, the gap was not as wide as one could expect, when comparing a sophisticated algorithm as the MHS with a simpler one. To shed light on this unexpected fact, we repeated the experimental campaign on a wider problem landscape and we observed that the solutions‟ space is rather flat (i.e., the average gradient is low), but it is disseminated of many local minima, with values close to the global one. As noted by Yang (2010), in these peculiar conditions performances of the Harmony Search (HS) get worse and the probability to get stuck, almost immediately, in a local minimum becomes high. Conversely, when the solution space is flat, only a search strategy assuring a high degree of randomness and little to none influence
CR IP T
to the initial starting point, may find the global optimum. To verify this assumption, we solved all the original problems‟ set using the very famous and general-purpose Simulated Annealing (SA) by Ishibuchi et al. (1995) that, due to a relatively simple iterated greedy strategy, fully satisfies the above-mentioned requirements. Obtained results were quite astonishing as the SA outperformed the MHS, in any operating condition.
Due to the above-mentioned issues, the speculative and operational goal of this paper is to ascertain if and
AN US
how this specific instance of the IJMS problem could be optimally solved also by a population-based metaheuristic. To this aim: (i) we modified the MHS, by introducing new local minima escaping techniques and (ii) we developed a new population-based metaheuristic that implements a searching strategy like the one of the Cuckoo Search (CS) by Yang and Deb (2009). Specifically, the CS has been chosen because, with respect to other population-based heuristics, it is less prone to early convergence, as cross-influence among
M
saved solutions plays a minor role in the evolutionary process.
The rest of the paper is organized as follows. Literature review on the IJMS problem is given in Section
ED
2, and the specific single machine IJMS problem considered in this paper is detailed in Section 3. New searching strategies and metaheuristics are introduced in Section 4 and 5 and their parameters are fine-tuned in Section 6. Furthermore, in Section 6 a comprehensive experimental campaign is made to assess and
PT
compare performance of the proposed metaheuristics. Lastly, conclusions and directions for future
CE
researches are drawn in Section 7.
2. LITERATURE REVIEW ON THE IJMS PROBLEM Scheduling problems have been generally treated under the assumption of full machines‟ availability.
AC
Only in the last decade this simplifying hypothesis has been questioned by several researchers, who postulated that, due to failures, the overall productivity of the manufacturing system may benefit from the inclusion of PMs between jobs. However, and quite surprisingly, almost the entirety of papers that start from this hypothesis do not fully consider the effect that PMs may have on machines‟ availability. Indeed, most of the researchers implicitly assume that, if properly performed, PMs will ensure sufficient machine availability to complete productions‟ orders in the allotted time. Basically, the problem is simplified as it follows: jobs‟ characteristics (in terms of processing times, set up times, release dates, routings, etc.) and the optimal number (or alternatively the optimal frequency) of the required PMs are considered as known inputs, and the scheduling problem reduces to the insertion of jobs between two consecutive PMs, in a way that optimizes a certain objective (such as the minimization of makespan or of the maximal tardiness), without violating (at
ACCEPTED MANUSCRIPT least within a certain tolerance limit) the expected PMs frequency. Yet, the possibility that faults may still occur is neglected, as well as the fact that failures‟ probability is strictly related to the time in which planned maintenance was performed. Among the most recent papers belonging to this category, the following ones can be mentioned. Mosheiov and Sarig (2009) extended the well-known scheduling problem of the “single machine with duewindows assignment”, by including in the schedule one additional maintenance activity. This maintenance activity is optional and, to be performed, it requires a fixed amount of time, during which production must be
CR IP T
halted. Although this maintenance task implies a waste of processing time, it may be beneficial because, once it has been performed, the single machine is supposed to become more efficient, and the successive workings will require a lower processing time. The authors also presented a polynomial time algorithm that generates the optimal schedule that minimizes costs related to jobs‟ earliness and tardiness and to the size and starting time of the due-windows.
Low et al. (2010) considered the problem of scheduling jobs on a single machine subject to periodic
AN US
maintenance activities. In this case, both the planned start time and the duration of each maintenance activity are known in advance, and the problems is to insert jobs between consecutive PMs, trying to minimize the makespan. To this aim the authors proposed an innovative version of the Particle Swarm algorithm that makes use of: (i) a “job-to-position” representation for the particles and (ii) the largest positional value rule to generate a particle sequence.
M
Zhu and Zhou (2015) also considered the single machine scheduling problem, but the peculiarity is that they introduced a deteriorating and resource-dependent “rate-modifying” planned maintenance. As usual, the
ED
hypothesis made is that planned maintenance alters (i.e., increases) the productivity of the machine. However, on the other one hand, planned maintenance is considered “deteriorating”, as its duration depends both on the resources allocated to it and on its starting time, i.e., the lateer it is performed, the worse the
PT
conditions of the machine become, and the longer the maintenance duration becomes. The authors also demonstrated that this problem can be optimally solved in a polynomial time with respect to many objective
date).
CE
functions (i.e., makespan, flowtime, maximum tardiness and combination of earliness, tardiness and due-
A similar problem was faced by Luo et al. (2015), who assumed that PM must start before a given deadline
AC
and that its duration increases with its starting time. The authors also provided polynomial-time algorithms to solve the problem with respect to the same objective functions used by Zhu and Zhou (2015). Similarly, Mor and Mosheiov (2015) studied the single machine due-windows assignment problem, in case of a deteriorating maintenance task. Specifically, two different types of deterioration were considered: timedependent (where the maintenance time increases as a function of its starting time), and position-dependent (where it is a function of its position in the sequence). In both scenarios a polynomial-time O(n4) optimisation algorithm is proposed to minimize costs related to jobs‟ earliness and tardiness and to the size and starting time of the due-windows.
ACCEPTED MANUSCRIPT Yu and Seif (2016) considered a flow shop scheduling problem, with the peculiarity that the operating conditions of each machine is expressed by a set of different Maintenance Levels (MLs). When a job is processed on a machine, each MLs is decreased and, as soon as a ML drops to zero, PM is performed. With the objective of minimizing tardiness and maintenance costs, the problem becomes NP-hard and it is solved with a lower-bound based genetic algorithm. Lee and Wang (2017) considered a single machine that processes two types of jobs, and the objective is to minimize the makespan of Type #1 jobs, provided that the maximum tardiness of Type #2 jobs does not
CR IP T
violate a predefined threshold limit. Additionally, PMs must be completed within specific time windows. Also, in this case, the number of PMs and their (deterministic) duration are known in advance and the problem is solved using branch and bound (based on a set of dominance rules) and via genetic algorithms. Similarly, also Touat et al. (2017) considered the problem of scheduling jobs on a single machine that requires PMs. In this case the specificity is that PMs must be performed inside predefined time windows, under human resources competences and availability constraints. To minimize the total tardiness, the authors
AN US
proposed two fuzzy genetic algorithms, differing in the way in which jobs, PMs and human resources are scheduled.
Cui et al. (2017) considered the single machine scheduling problem with jobs‟ release dates. The single machine requires flexible PMs that are activated following a “run-based” policy, i.e., a maintenance task must be performed any time the machine working age exceeds a predefined threshold value. To solve the
M
problem minimizing the makespan, a performing heuristic (based on some properties of the optimal solution) and a branch and bound (that utilizes several dominance rules) are proposed for large and medium sized
ED
problems, respectively.
Yazdani et al. (2017) considered a single machine scheduling problem characterized by multiple unavailability periods, occurring at known and constant intervals, during which overhauls are performed.
PT
Aiming to minimize the sum of earliness-tardiness penalties, the authors proposed an algorithm characterized by a variable neighbourhood search that, in addition to local search, also implements a “shake” procedure,
CE
which diversifies solutions, by performing a leap toward another local neighbourhood. The algorithm is also integrated with a “knowledge module” that keeps track of the structure of good solutions and that takes advantage of it during the search.
AC
El-Khoukhi et al. (2017) investigated the flexible job shop scheduling problem, when a predefined number of PMs, of known and deterministic length, are needed. To minimize makespan, they solved the problem using an innovative Dual-Ants Colony Optimisation. The metaheuristic is also integrated with a dynamic memory, which takes advantage of local search strategies and several dispatching rules. Woo and Kim (2018) considered the parallel machines scheduling problem characterized by “ratemodifying” PMs activities. Specifically, each machine is characterized by a starting processing time that increases through time, accordingly to a certain deterioration rate. However, if PMs is performed on a machine, the processing time is restored to its original value. Given a certain number of time slots where PMs can be performed and considering as known and deterministic the PMs‟ time, the main decisions to be
ACCEPTED MANUSCRIPT taken concern: (i) the number and position of PMs and (ii) the ordered assignment of jobs on parallel machines. Specifically, to minimize the makespan, the authors proposed a specific heuristic that takes advantage of Simulated Annealing and of Genetic Algorithm to decompose the problems into sub-problems than can be analytically and optimally solved. As we have already said, all these papers consider PMs as a prerequisite needed to assure a minimum level of availability, as required by the manufacturing system. Therefore, the number and the duration of the PMs are assumed to be known and the positive effect that PMs may have on the hazard rate of the machines
CR IP T
is disregarded. Very few papers reduced these drawbacks, by considering a time-dependent hazard rate and by including in the analysis also the additional losses that may be generated by faults and by reactive maintenance tasks. As far as the knowledge of the authors, in addition to the work by Zammori et al. (2014), which will be discussed in depth in Section 3, this approach was followed only by eight recent papers that share the following common hypotheses: (i) Machines‟ deterioration increases through time and it is described by a time dependent hazard rate, (ii) when performed, planned maintenance restores the original
AN US
operating conditions of the machine (i.e., perfect repair) or, at least, it reduces the current degradation state, (iii) although PMs are performed, failures may still occur. In this case reactive maintenance tasks are needed to restore the machines, nonetheless their degradation remains unchanged (i.e., minimal repair model). Safari et al. (2011) are, probably, the first authors that fully considered the effects of PMs on the hazard rate. Specifically, the authors deal with a flow-shop configuration under the assumption of condition-based
M
maintenance to minimize expected makespan. More precisely, every T time units, each machine is inspected and, depending on its current degradation (modelled with a Poisson distribution) a reparation may be
ED
performed. Also, failures can occur in between two consecutive inspections, according to a nonhomogeneous Poisson process of parameter
. A hybrid heuristic based on genetic algorithm and
simulated annealing is also proposed to minimize makespan.
PT
Cui et al. (2014) tackles the problem of scheduling jobs and PMs on a single machine subject to unexpected failures. The objective function is evaluated both in terms of solution robustness and of solution quality.
CE
Specifically, the first metric is computed as the total deviation between the planned start time and the actual start time of the jobs, whereas the second metric is measured as the maximum tardiness of the jobs. Also, to generate optimal solutions, a three-phased heuristic is proposed. At first the jobs‟ sequence and the relative
AC
position of PMs is determined; next some extra idle time (or “time buffer”) is included in the schedule to protect due dates from unpredictable failures. Lastly, in the third and final phase, a pair-wise-swap-based local search is performed on the jobs‟ sequence, aiming to improve the initial solution. Lu et al. (2015) tackled the problem of finding a robust and stable schedule for a single machine that suffers of unexpected breakdowns. According to the authors, the best way to maximise robustness is to: (i) add time buffers to protect the initial schedule and (ii) execute PMs anytime machine‟s reliability drops below a threshold limit. Owing to these issues, the authors solved the problem with a genetic algorithm based on a double chromosome: the first chromosome codifies (with natural number) the jobs‟ sequence, the second one codifies (in binary code) the desired reliability limit.
ACCEPTED MANUSCRIPT Mokhtari et al. (2015) considered a flexible flow-shop scheduling problem, in which machines‟ failure rates are time-varying and proposed a heuristic, based on Simulated Annealing and on Monte Carlo simulation, to minimize the number of tardy jobs. The job shop is considered flexible because any operation can be performed by more than one machine and failure rate is considered time-varying as it depends on environmental situations like temperature, light, humidity, etc. Quite interestingly, this model is based on the hypothesis that the time needed to complete a batch of orders is much smaller than the average time between two planned maintenance tasks. Consequently, the use of PMs is not contemplated and, indeed, only reactive
CR IP T
maintenance tasks are executed anytime a failure occurs. Nouiri et al. (2017) focused on the flexible job shop scheduling problems under machines breakdowns. Makespan minimization is considered as the main objective, but, as secondary objective, the obtained solution must guarantee a stable schedule even in case of breakdowns. To this aim a two-stages particle swarm optimization was proposed: at first breakdowns are disregarded and an initial sequence is found, next breakdowns are randomly generated (supposing that machines with the heavier workload are more likely to
AN US
fail) and the original solution is modified to assure robustness and stability. As in the work by Mokhtari et al. (2015), also in this case PMs are not considered, and only reactive maintenance is used to counteract machines‟ failures.
Salmasnia and Mirabadi-Dastjerd (2017) considered the scheduling problem of a degraded single machine subjected to failures. Specifically, they considered different levels of PMs: although planned maintenance
M
always leads to machine‟s rejuvenation, only if the top and most expensive kind of planned maintenance is performed the machine is restored to its original state (i.e., perfect repair). The objective is to jointly
ED
minimize the total maintenance cost and the total jobs completion time. Since there is a clear trade-off between these objectives, the Pareto Set of the non-dominated solutions was generated using a genetic algorithm.
PT
Rahmati et al. (2017) extended the flexible job shop scheduling problem presented by Yu and Seif (2016), by considering machines characterized by an increasing hazard rate. Briefly, to schedule jobs and PMs, a
CE
simulation-based optimization similar to the one by Nouiri et al. (2017) was proposed. At first, Harmony Search is used to find both the jobs‟ sequence and the machines‟ assignments that minimize the makespan; next the optimal time between PMs is defined via simulation.
AC
Cui et al. (2018) dealt with the integration of production scheduling and maintenance planning to optimize robustness and solution quality (evaluated as in Cui et al. (2014)), for a flow shop with failure uncertainty. As in Lu et al. (2015), buffer times are proactively added into the initial schedule and an initial solution is obtained following a constructive mathematical model. Next, the solution is improved using a two loops algorithm that takes advantage of numerical simulation and Monte Carlo sampling methods. More precisely, the outer loop performs a local search to optimize jobs‟ sequence and PMs‟ positions; the inner loop implements a genetic algorithm that tries to optimize the idle times when jobs‟ sequence and PMs‟ positions are fixed. Results of the inner loop are returned as an evaluation criterion to guide the search direction of the outer loop.
ACCEPTED MANUSCRIPT For the sake of clarity, results of the literature review on IJMS problems published starting from 2014 are summarized in Table I. Although the considered models are case specific and cannot be directly compared,
AC
CE
PT
ED
M
AN US
CR IP T
we tried to summarize, in the table, pros and cons of each one of them.
CR IP T
ACCEPTED MANUSCRIPT
Table I Main papers dealing with the IJMS problem from 2014 to now Paper
System
Objective Function
Features of PMs
Solution Approach
Cui et al. (2014)
Single Machine
PMs executed at variable time intervals to increase machine‟s availability.
Schedule Robustness and Quality
Zammori et al. (2014)
Single Machine.
PMs executed at variable time intervals to increase machine‟s availability.
Earliness and Tardiness Penalties.
Luo et al. (2015)
Single Machine
PM is used to improve machine‟s productivity; PM is “deteriorating”, as is duration increases over time.
Makespan; Other classic objectives.
Lu et al. (2015)
Single Machine
PM is executed any time drops below a threshold limit, to restore its initial state.
Robust Schedule Stable Schedule
Mokhtari et al. (2015)
Flexible Flow Shop
Machines are characterized by a time-varying failure rate. When faults occur, maintenance is needed to reactivate the system.
Number of tardy jobs
Mor and Mosheiov (2015)
Single Machine
PM is used to improve machine‟s productivity; PM is “deteriorating”, i.e. duration is time and position dependent.
Costs minimization.
Optimal Polynomial Time Algorithm.
Zhu and Zhou (2015)
Single Machine
Makespan; Other classic objectives.
Optimal Polynomial Time Algorithm.
Yu and Seif (2016)
Flow Shop
Tardiness and Maintenance costs
Genetic Algorithm
Cui et al. (2017)
Single Machine.
Makespan.
Specific Heuristic; Branch and Bound procedure.
Hazard rate is considered. Robustness is obtained using time buffers and the number of PMs is analytically optimized. Hazard rate is considered. Failures and reactive maintenance are considered. The number of PMs is analytically optimized. A deteriorating rate, affecting PM‟s time is considered; Optimal solution can be obtained in polynomial time Hazard rate is considered; Failures are considered; The optimal threshold limit for is optimized.
AN US Modified Harmony Search. Optimal Polynomial Time Algorithm.
Genetic Algorithm with double chromosome Simulated Annealing and Monte Carlo simulation
M
ED
AC
CE
PT
PM is used to improve machine‟s productivity; PM is “deteriorating”, as is duration increases over time. Machines degrade over time and their productivity decrease. PMs are needed as soon as the Maintenance Level drop to zero. The machine degrades and its productivity decreases. PMs are executed, anytime the degradation exceeds a threshold.
Three-phased heuristic
Pros
Time varying hazard rate depending on environmental conditions. Two kinds of deteriorating rates are considered, and several costs are considered. Optimal solution can be obtained in polynomial time. A deteriorating rate, affecting PM‟s time, is considered; Optimal solution can be obtained in polynomial time. Different Maintenance Levels can be used. The number and the kind of PMs to be executed is an output of the genetic algorithm. Introduces a “degradation” concept. PMs executed accordingly to a runbased policy to reduce degradation.
Cons
A single level of PM (i.e., perfect repair) is considered.
A single level of PM (i.e., perfect repair) is considered. Solution approach is not very performing. Hazard rate is not considered; Faults are not considered; A single maintenance task must be included in the schedule. A single level of PM (i.e., perfect repair) is considered.
Only faults and reactive maintenance tasks are considered. PMs are not considered as a way to improve availability. Hazard rate is not considered; Faults are not considered; At most one maintenance task can be included in the schedule. Hazard rate is not considered; Faults are not considered; At most one maintenance task can be included in the schedule. Hazard rate is not considered (i.e., it is substituted by a degradation rate); Faults are not considered; Hazard rate is not considered (i.e., it is substituted by a degradation rate); Faults are not considered;
ElKhoukhi et al. (2017)
Flexible Flow Shop.
Lee and Wang (2017)
Single Machine.
Nouiri et al. (2017)
PMs are needed to assure the required availability and do not alter machine productivity.
Dual Ant Colony with dynamic memory
PMs are needed to assure the required availability and must be performed within predefined time windows.
Makespan; Tardiness.
Genetic Algorithm.
Flexible Job Shop.
Fault probability depends on the work load; anytime a fault occurs, reactive maintenance is needed.
Makespan
Rahmati et al. (2017)
Flexible Job Shop.
PMs executed at variable time intervals to increase machine‟s availability.
Makespan.
Touat et al. (2017)
Single Machine.
PMs are needed to assure the required availability and must be performed within predefined time windows.
Total Tardiness.
Salmasnia and Mirabadi (2017)
Single Machine
PMs executed at variable time intervals to increase machine‟s availability.
Total Cost Total Compl. Time
Yazdani et al. (2017)
Single Machine
Overhauls are performed (stopping production) at predefine time intervals.
Earliness and Tardiness Penalties
Variable Neighborhood Search
Cui et al. (2018)
Flow Shop
PMs executed at variable time intervals to increase machine‟s availability.
Schedule Robustness and Quality
Dispatching Rules, Simulations and Monte Carlo Sampling
Woo and Kim (2018)
Parallel Machines
The machine degrades and its productivity decreases. PMs must be executed, within time windows to reduce degradation.
Makespan
Simulated Annealing and Genetic Algorithm
Optimization via simulation
Harmony Search and Simulation
Fuzzy Genetic Algorithms.
Genetic Algorithm with double chromosome
M
ED
PT
CE
AC
The model is standard, the innovative part lies in the solution approach that is very performing. Number and duration of PMs are known, but their starting time can vary. Jobs are classified into families with different importance.
AN US
Makespan
CR IP T
ACCEPTED MANUSCRIPT
Failure rate is time varying and it also depends on work load. Hazard rate is considered. Failures and reactive maintenance are considered. The number of PMs is analytically optimized. PMs must be placed within predefined time windows, that also depends on human resources competences and availability constraints. Hazard rate , failures and reactive maintenance are considered. Different levels of PM are considered, and their number is analytically optimized. The model is standard, the innovative part lies in the solution approach that is very performing. Hazard rate is considered. Failures and reactive maintenance are considered. The number of PMs is analytically optimized. Time buffers are used to increase robustness A machines‟ deteriorating rate is considered. The optimal number of time windows (and their duration) is optimized.
Hazard rate is not considered; Faults are not considered; Number and duration of PMs are known in advance. Hazard rate is not considered; Faults are not considered; Number and duration of PMs are known in advance. Only faults and reactive maintenance tasks are considered. PMs are not considered as a way to improve availability. The simulation model is case specific A single level of PM (i.e., perfect repair) is considered. Optimization is obtained via simulation and it is case specific. Hazard rate is not considered; Faults are not considered; Number and duration of PMs are known in advance. Due to tradeoff among the objectives, a global optimum cannot be found; only a Pareto Optimal frontier can be generated. Hazard rate is not considered; Faults are not considered; Number and duration of PMs are known in advance.
A single level of PM (i.e., perfect repair) is considered. Optimization is obtained via simulation and it is case specific. Hazard rate λ(t) is not considered (i.e., it is substituted by a degradation rate); Faults are not considered;
ACCEPTED MANUSCRIPT
3. THE 3.1. In this paper, we focus on the single machine IJMS problem introduced by Zammori et al. (2014), whose features are listed below:
The single machine must process a certain number of jobs, all available at the starting time
For each job, its processing time, due date and sequence dependent setup time are known;
The machine is subjected to failures and production may be halted both to perform reactive and
;
planned maintenance tasks; A time dependent hazard rate
is considered and minimal and perfect repair models are used to
CR IP T
describe machine‟s restoration after the execution of a reactive and of a planned maintenance task, respectively;
The goal is to find the integrated sequence of jobs and planned maintenance tasks that minimizes the sum of Earliness-Tardiness Penalties (ETP).
AN US
In short, using the Kendall‟s notation, the problem can be classified as: (1| Due Dates, Sequence Dependent Set-Up Times, Planned Maintenance| Earliness-Tardiness Penalties). 3.2. Computing the expected completion time of the jobs Let
and
be, respectively, the number of jobs and of PMs. Then, any possible solution can be
expressed as a combination of two vector and. The vector possible permutations of
{
jobs and each one of its element
M
the
(
be processed in position . Similarly,
represents one of } codifies the job that will
) represents one of the (
) possible ways
ED
in which PMs can be inserted between the first and the last job; in this case, each element of index of the job (within ) that precedes the j-th PM. For instance, with } would be denoted as:
,
PT
{
[
]
CE
Using this notation, the expected completion time for each job *
AC
∑(
[
[ {
)
(∑
and
codifies the , the schedule
.
+ can be expressed as in Eq. (1):
) (1)
]
∑
(
)
(
]
∑
(
)
(
∑
∑
)
)
ACCEPTED MANUSCRIPT
where: {
} is the set of jobs;
is the vector of the jobs‟ processing times; (
) is the
matrix of the sequence dependent setup times; is the vector of the due dates;
is the average time of a reactive maintenance task;
is the average time of a planned maintenance task;
is the expected number of failures in the interval [
].
CR IP T
According to Eq. (1) the expected completion time of a job can be computed using a recursive formula that takes into consideration the position of the job within the scheduling sequence, relatively to the nearest PM. More precisely, the whole schedule is seen as an aggregate of sub sequences of jobs divided by
AN US
consecutive PMs. Also, due to the perfect repair hypothesis, all these sequences can be considered independent because any time PM is performed, the hazard rate consequently, the interval [
is reset to its original value and,
], on which the expected number of failures
is computed, shrinks
accordingly. Owing to these issues, the completion time of a job belonging to a generic subsequence equals the sum of: (i) the completion time of the last subsequence
and (iii) the expected down time, due to both failures and
M
the preceding jobs in the same subsequence
, (ii) the processing and setup times of all
maintenance, from the start of subsequence up to the considered job.
becomes: (2)
can be easily obtained as:
CE
Hence,
shape parameter , the hazard rate
PT
with characteristic life
ED
Also, under the reasonable and common assumption of the time to failure being distributed as a Weibull,
( )
(3)
AC
Based on the expected completion times, ETP can be finally evaluated as in (4), where earliness and tardiness penalties (per unit of time) relative to job ∑(
and
are the
. )
(4)
where: ,
-
(5)
ACCEPTED MANUSCRIPT
|
,
-|
(6)
3.3. Solution approach Zammori et al. (2014) also proposed a Modified Harmony Search (MHS) to solve the (1| Due Dates, Sequence Dependent Set-Up Times, Planned Maintenance| ETP) scheduling problem. Basically, the MHS is an advancement of a standard Harmony Search (HS), a metaheuristic inspired by the creativity process of
CR IP T
jazz musicians that, during jam sessions, expand their repertoire using pure inspiration and/or reworking on old harmonies (Geem et al., 2001).
Formally, the HS has its essence in the Harmony Memory (HM), a matrix that represents a musician‟s repertoire. This matrix contains
row vectors
[
, where the superscript
] indicates the position in
the HM, each one coding a feasible solution and its Objective Function (OF). For the present problem, and
are joined to form an
whose generic element
(with
(
-dimensional vector
AN US
vectors
),
) contains the index of the i-th task of the integrated job-
maintenance schedule and the element
contains the ETP of the schedule. Since tasks may be either
jobs or maintenance, natural numbers in the interval [ jobs and PMs, respectively.
] and in the interval [
] denote
M
Randomly initialized, solutions of the HM follow an evolutionary process based, mainly, on the Harmony Memory Consideration (HMC). Specifically, new solutions
taken from some “promising” solutions stored in the HM. To avoid rapid
ED
a time, the elements
are generated by recombining, one element at
convergence, selection is made using a roulette wheel selection rule so that, although rarely, even the elements of a solution
with a poor objective function may be selected. For the same aim, two additional
PT
strategies, namely Random Selection (RS) and Pitch Adjustment (PA), are used to convey a certain degree of randomness to the generative process. According to the first strategy, to create a new solution, some may be randomly generated from the values of the domain [
] to which they belong
CE
elements
to. Similarly, PA is used to slightly perturb one or few elements of the new solution. For instance, in the
AC
MHS, PA is implemented as a pairwise interchange generation scheme, that is two positions and new solution
are randomly selected and the corresponding elements
and
of the
are swapped. Lastly, a new
solution is accepted if its objective function is better than that of the worst solution stored in the HM; in this case, the worst solution is replaced by the new one. The peculiarity of the MHS is the implementation of an additional local minima escaping strategy, called Large Portion Recovery (LPR). Anytime LPR is triggered, the HM is completely regenerated with new solutions created by recombining, through 2-points cross-over, entire sub-arrays of the precedent solutions contained in the HM. LPR is activated with a rate (LPRR) that depends on the Saturation‟s Level of the HM, a measure of similarity of the solutions of the HM, computed as in Eq. (7):
ACCEPTED MANUSCRIPT
(7) Where
and
are, respectively, the number of copies and the total number of solutions of the HM.
3.4. Problem peculiarities As anticipated in Section 1, when the ETP is used as objective function, regardless of others problem settings, the solution space becomes rather flat, but it is scattered by many points of local minimum. For the
CR IP T
sake of brevity, we will not report the detailed results of this preliminary analysis, but we will show a simple example, which illustrates this concept well. Specifically, Figure 1.a shows the ETP of all the 3600 solutions of three randomly generated instances of a 6-Jobs, 1-PM scheduling problem. In each problem, the average job‟s processing time equals 100-time units, planned and reactive maintenance equals 40 and 20-time units, respectively. What differentiates each case is the value of the shape parameter
of the Weibull probability
AN US
distribution, which equals 90%, 60% and 30% of the average jobs‟ processing time. ------------------------------Insert Figure 1.a here -------------------------------All contiguous points on the graph (i.e., points labelled as
,
) are neighbour solutions generated as
M
follows. At first, starting from the basic sequence {1, 2, 3, 4, 5, 6}, its 720 permutations were generated using the shifting procedure i.e., two locations if
and
are randomly selected and job
is moved before
and vice-versa. Next, all these permutations were replicated for five times and PM was added
ED
after job in position 1 for the first group of 720 permutations, in position 2 for the next 720 permutations, and
CE
PT
so on, as shown in the generation scheme of Table II:
,
… …
… ,
,
…
AC
,
Table II The neighbors‟ generation scheme
As it can be seen from Figure 1.a, the general pattern does not significantly change: neighbor solutions
tend to scatter, almost randomly, around an average value that does not display a significant trend. The fact that the solution space is rather flat is even better highlighted by Figure 1.b, which shows the ETP curves obtained after sorting (in ascending order) all the 3600 solutions of the problem. ------------------------------Insert Figure 1.b here --------------------------------
ACCEPTED MANUSCRIPT Lastly, we note that also the number of local minima is high, as shown by Table III, which reports the main features of the above-mentioned example.
Table III Basic features of the problems
= 60
= 30
Minimum Cost
10,094
13,455
19,838
Maximum Cost
13,067
17,861
26,608
Solutions within 10% from the optimum
404
703
698
Solutions within 5% from the optimum
82
165
170
AN US
3.5. Computational drawbacks
CR IP T
= 90
As noted by Yang (2010), when the solution space is so flat, neither Random Selection nor Pitch Adjustment can convey enough randomness to effectively escape an early convergence process: all the solutions of the Harmony Memory become similar, if not even identical, and the probability to get immediately stuck in a local minimum becomes high. Although the MHS implements a third local minimum
M
escaping strategy (i.e., the Large Portion Recovery), it is licit to ask if this approach is sufficient to obviate the early convergence problem. To answer this question, we compared MHS with the very famous and general-purpose Simulated Annealing (SA) by Ishibuchi et al. (1995), whose main features are detailed in
ED
Appendix A. The rational of this choice can be motivated as follows: The metaheuristic implements a relatively simple iterated greedy strategy that assures a high degree
PT
of randomness during the exploration of the solution space; The metaheuristic is flexible, as its performance does not depend, significantly, neither on the
CE
specific nature of the combinatorial problem under analysis, nor on the initial condition used as the starting point;
The number of operating parameters to be fine-tuned is very low.
AC
We also note that we opted for a general-purpose algorithm rather than for one specifically designed to tackle IJMS problems, as it was not possible to find an algorithm suitable for the specific IJMS under consideration. Indeed, as clearly discussed in Section 2 and showed in Table I, most of the approaches proposed in the literature have the limit not to consider the effect of PMs on the hazard rate of the machines, which is, instead, one of the peculiar feature of the model herein considered. Even the few works that considered the hazard rate, cannot be used as benchmark, since they are very specific and cannot be effectively readapted to the framework of the problem under consideration. Obtained results were astonishing, as SA outperformed MHS, even in case of 10 randomly generated problems of small size (with parameters defined as indicated in Section 6), characterized by 14-Jobs and 2-
ACCEPTED MANUSCRIPT PMs and for which the optimum was found through an exhaustive search. More precisely, each problem instance was solved 8 times and obtained results where quantified in terms of Average Percentage Error (APE) and Average Percentage Deviation (APD), computed as in Eq. (8) and Eq. (9), respectively. Note that quantifies the average difference between the best solution found by the j-th algorithm and the real optimal one. Similarly,
measures how much, on average, the j-th algorithm diverges from the real
optimum. (
Where:
(8)
(
)
(9)
AN US
∑∑
)
CR IP T
∑
is the k-th solution generated by the j-th algorithm for the i-th problem; in this case [
],
] and
[
];
is the best solution for problem found by algorithm i.e.,
{
};
is a priori known global optimum for problem .
M
[
In case of SA, both APE and APD were equal to zero; conversely, for MHS an APE of 0.5% and an APD of 1.8% were found. These values demonstrated that, whereas SA consistently found the global minimum,
ED
MHS always converged to some local minimum rather distant from the global one.
PT
4. ENHANCEMENTS OF THE MHS
Due to the preliminary results discussed in Section 3, we tried to improve the MHS through the addition
CE
of four additional local minima escaping strategies that will be detailed in the following Sub-Sections. 4.1.
Block Shifting: a new Pitch Adjustment procedure
As explained in Section 3, in the MHS PA has been implemented as a simple pairwise interchange
AC
operator. Consequently, due to the small number of PMs compared to the number of jobs, it is quite rare for the PA to swap a job with a PM. In general, this is not an issue and the evolving process works fine, as PMs repositioning is assured by the Harmony Memory Consideration. Yet, when the Harmony Memory begins to saturate, Harmony Memory Consideration loses effectiveness and the exploration of the neighborhood proceeds mainly thanks to PA. To exploit this search, it is thus fundamental to add to the original pairwise interchange, which modifies jobs‟ sequencing leaving almost unaltered the position of PMs, a dual PA operator that, conversely, acts on PMs leaving almost unaltered jobs‟ sequencing. Owing to these considerations, we defined the Block Shifting (BS) operator, which works as follows:
A whole subsequence of jobs (i.e., jobs placed between two PMs) is randomly selected;
ACCEPTED MANUSCRIPT
Let
be the selected subsequence, then a new neighbor solution is generated by moving the block
,
- either after the first job of the subsequence
subsequence
.
{
For instance, in case of sequence subsequence
}, if
} was selected and block {
{
or before the last jobs of the
{
} was moved forward, the sequence } would be finally obtained. Please note that
could have been generate also by means of pairwise interchange; yet this would be extremely improbable
CR IP T
because, considering that sequences with 2 PM in adjacent positions are not allowed, it would require a minimum of five consecutive moves. 4.2.
A second crossover operator
Aiming to improve the regeneration of the HM through Large Portion Recovery, in addition to the 2points crossover, we also considered and compared two alternative types of crossover: (i) the order crossover
AN US
and (ii) the double 2-1-points crossover, readapted from the recent one proposed by Salmasnia and MirabadiDastjerd (2017).
The ordered crossover was chosen to foster the exploration of large portions of the solution space, as it is known to have the ability to generate offspring solutions sensibly different from the parent ones (Gen et Cheng, 2000). Briefly, we recall that the order crossover is based on the following steps: Select a random swath of consecutive genes from Parent 1;
Drop the swath down to Child 1 and mark out these genes in Parent 2;
Starting on the right side of the swath, take genes from Parent 2 and insert them in Child 1, at the
ED
M
right edge of the swath;
Generate Child 2 by flipping Parent 1 and Parent 2 and repeating all the other steps.
highlighted in Italics:
{
AC
CE
{
PT
To clarify this process, a simple example is given next, where the swath of genes selected from Parent 1 is }, },
{
},
{
}.
The integrated 2-1-points crossover is the only alternative crossover operator that was found in the
literature, among the subset of papers that deal with the IJMS problem and that fully consider the effect of PM on machines‟ hazard rate (see Table II in Section 2, for further details). In its original version, 2-1-points crossover operates on two distinct chromosomes that codify jobs and different levels of PM, respectively. However, in the present problem there is a single kind of PM, and so we used the second gene to codify the position of PMs within the jobs‟ sequence. In virtue of this slight modification, the 2-1-points crossover operates as it follows:
Jobs and PMs are separated and coded with the
-
notation of Section 3.1;
ACCEPTED MANUSCRIPT
2-points crossover and 1-point crossover are made on the
vector (i.e., jobs part) and on the
vector (i.e., maintenance part), respectively;
The two constructed parts are recombined together to generate a new solution.
To clarify this process, the integrated 2-1-points crossover is applied to the parents used in the example above; note that both vectors have been recoded using the
notation and that the symbol ˈ|ˈ has been
-
used to highlight the cutting points. {
}
|
|
|
,
{
}
|
|
|
,
CR IP T
,
{
},
{
}.
Lastly, aiming to select the best crossover operator we: (i) developed three standard genetic algorithms,
AN US
differentiated only by the type of crossover and (ii) we tested them on the same problems‟ set described in Section 3.4. At the end of the analysis the choice fell on the ordered crossover, as it was the only one able to find the best solution for every instance of the problems set. 4.3.
Maximizing the Harmony Memory degree of differentiation
In the original MHS, anytime the Large Portion Recovery is invoked, crossover is activated to regenerate
M
the HM. Ideally, at the end of this process, the HM should be completely new, yet, as clearly shown by Figure 2, the true picture is rather different.
ED
------------------------------Insert Figure 2 here
PT
--------------------------------
As it can be seen, anytime the Saturation Level (SL) reaches a threshold limit (in this case set to 0.3), the Large Portion Recovery is triggered and, immediately after, the SL drops down; yet, unlike one could expect,
CE
it is always greater than zero. This fact implies that, due to an excessive degree of similarity among parents, duplicates are created in the offspring. Conversely, to maximize the degree of differentiations of the HM,
AC
without losing good solutions, we decided to limit the regenerative process only to the duplicated solutions, which are substituted with newly generated ones. More precisely, letting
,
and
be,
respectively, the number of solutions, of copies and of original (non-duplicated) solutions contained in the HM, the regenerative process works as follows:
Identify all different solutions contained within HM;
Use a roulette wheel selection (without replacement if select
, with replacement otherwise) to
parents among the original solutions and apply order-crossover, to generate
offspring
solutions
Selectively replace all the copies contained in the HM using the
generated offspring solutions.
ACCEPTED MANUSCRIPT It is worth noting that, the regeneration of the HM may introduce solutions that are much worse than the current best, as indicated by the spikes in the “worst value curve” (in red) shown in Figure 3. Anyhow, as indicated by the “current optimum curve” (in green), the regenerative process strongly improves the exploration of the solution space and, indeed, improvements often occur immediately after a regenerative process has taken place. ------------------------------Insert Figure 3 here
4.4.
A dynamic Large Portion Recovery Rate
As a last modification, instead of using a fixed
, we decided to dynamically update this threshold
limit, using the cubic function of Eq. (10). In this way, rapidly rises as
CR IP T
--------------------------------
remains low for small values of
increases. This behavior is rationale because when
, but it
is high almost all the solutions of
AN US
the HM are identical and only crossover can refresh the HM. For this reason, also
must approach one,
to activate almost certainly the Large Portion Recovery.
{
ED
Where:
is the value of
PT
is the minimum value of
CE
(10)
⁄(
;
below which
is the value of
)
M
[
above which
takes its minimum value; becomes 1;
).
AC
We conclude this section noting that all the above-mentioned improvements take advantage of the specific structure of the IJMS problem herein considered. Consequently, from here on, we will refer to this metaheuristic as the Harmony Search for Integrated Job Maintenance Scheduling, indicated as HS/4/IJMS. For further details the interested readers are referred to Appendix B, where the full pseudo-code can be found. 5.
CUCKOO SEARCH VIA LÉVY FLIGHT AND POSSIBLE ENHANCMENTS We also decided to investigate the performances of the Cuckoo Search via Lévy Flight (CS-LF), a
popular population-based metaheuristic introduced by Yang and Deb in 2009. CS-LF was originally designed to solve continuous problems, but its preformances are excellent also in case of discrete
ACCEPTED MANUSCRIPT combinatorial problems (see for example Ouaarab et al., 2014; Komaki et al., 2017). Also, and perhaps more important, since its search strategy assures little or none cross influence among saved solutions, CS-LF is a very good candidate to tackle the IJMS problem. 5.1.
Main Features of the standard CS-LF
CS-LF is a nature inspired metaheuristic, which mimics the brood parasitism of some species of cuckoo, a bird which does not brood its eggs but lays them in other nests. Each cuckoo explores its surroundings, by means of a Lévy Flight (LF), in search of a good nest to lay its egg. Evaluating the options, the cuckoo will
CR IP T
select the best one and it will lay its egg. Once an egg has been laid, it can be detected and purged by the owner of the nest, or it can successfully hatch. If so, the new cuckoo will continue the exploration of the surrounding area, looking for new nests to colonize.
Briefly, the key element of the CS-LF is that, at each step, a pairwise interchange scheme is iteratively applied to generate a neighbor
from an initial solution
, then a neighbor
and so on (i.e.,
, a parameter that represents the
AN US
) for, at most, a number of times equal to
from
distance travelled by the cuckoo. Also, this travel is modelled as a Lévy Flight that is, essentially, a random walk in which the length of the random step is drawn from a probability distribution that is heavy-tailed i.e., for which
, with
, so that the variance is infinite. For
, the Pareto
distribution (see Eq. (11)) is one of the distributions satisfying these requirements. (
M
{
)
(11)
ED
So, using a Pareto, in virtue of the inverse transformation techniques, the step length
is a random number uniformly distributed in the interval [
(12) ].
CE
where
PT
Eq. (12):
can be generated as in
Also, to model the possibility that some of the eggs are detected and smashed, a
of the worst solutions are
substituted with randomly generated ones. For the full algorithm, the interested readers are referred to Yang
AC
and Deb (2009); an extended and modified version is detailed in Sections 5.2 and 5.3. 5.2.
Cuckoo Search with Migration and Nests Location Memory (CS-MNM)
We also altered the searching strategy implemented by the CS-LF to explore the solution space in an even
better way. The key of such re-design is based on a modification of the standard Lévy Flight and on the introduction of three different features, namely: (i) Nests Location Memory, (ii) Risk Aversion and (iii) Clan Migration. 5.2.1
Double Shaped Lévy Flight
As shown in Figure 4, two flight modes are possible:
Mode #1 – A flight of progressive departure from the starting point
;
ACCEPTED MANUSCRIPT
Mode #2 – A peripheral flight centered on an intermediate point .
Note that, starting from
(highligthed in black in Figure 4), the overall flight (from 0 to 7) can be
decomposed into a series of intermediate points
(each one representing a feasible solution) where the
cuckoo pauses to assess the goodness of the surranding area. As long as the quality improve, with respect to the last visited point, the cuckoo continues to move away (as for points { , (as for points { ,
,
} and { ,
}). Otherwise
}), it explores the neighborhood of the last point of improvement (as for point
). The
latter case can be seen as the cuckoo coming back and then restarting its search or, simply, as the cuckoo
number of visited points
CR IP T
performing a sort of circular flight centered on the last point of improvement. Anyhow, the maximum is defined by the Lévy parameter of Eq. (12) and, as in the standard CS-LF, the
search is immediately interrupted if a point
improves the overall best solution found so far. In this case the
cuckoo stops its search and lays an egg. -------------------------------
AN US
Insert Figure 4 here -------------------------------Operatively, starting from
, all the visited points
are generated applying the pairwise interchange
operator to the last solution that produced an improvement. More precisely the first point neighbor of
; next, if
is better than
, provided that
M
is generated as a neighbor of
, the second point
is generated as a
is generated as a neighbor of
, otherwise
. This process is iteratively repeated until at most
point have been explored. For instance, in the case of Figure 4, the generating scheme is the following one:
5.2.2
ED
{
}.
Nest Location Memory and Risk Aversion (RA)
PT
Most of the times, the cuckoo starts a new flight from the last point of improvement found during its last flight (i.e.,
in the example of Figure 4). This always happens in case of a successful
CE
mission when a new host nest has been identified and the cuckoo has laid its egg. This also happens in case of an unsuccessful mission, provided that Mode #1 has been used at least once, during the Double Shaped Lévy Flight. Conversely, if only Mode #2 has been used, the cuckoo may try to exploit all the information
AC
obtained during its peripheral flight. If so, it might decide not to come back to the initial starting point rather to restart from the best intermediate point
, but
found during the flight, although this point was judged
inadequate to lay an egg and it is probably a bad move in the short-term. Clearly this choice carries a risk that is accepted by the cuckoo with a certain probability (namely the Non-Return Probability NRP). Also, the NRP is not fixed, but, as in Simulated Annealing, it is exponentially decreased over time, to better represent the risk aversion profile of the cuckoo, which is supposed to become more and more circumspect as it gets older. This is shown in Eq. (13).
ACCEPTED MANUSCRIPT
{ Where cuckoo
.
is the change in the objective function and
/}
, is the parameter that links the risk aversion of
with the number of performed iterations . To this aim
from an initial value
, until it reaches an ending value
is progressively and linearly reduced
at the last iteration , as shown in Eq. (14): )
(14)
CR IP T
( Concerning
(13)
, as suggested by Kirkpatrick et al. (1983), this value can be chosen to assure a high initial
Non-Return Probability
as in Eq. (15): ̅̅̅̅
(15)
Eq. (15) can be used to compute 5.2.3
AN US
where ̅̅̅̅ is the average change of the objective function over a sample of neighbor solutions. Also note that too, by substituting
Clan Migration (CM)
with a low final Non-Return Probability
.
We also introduced a migration strategy, to mimic the ability of birds to share part of their past experience. Specifically, since members of a clan tend to follow the behavior of the most successful
M
individuals, we hypothesized that, after a certain number of iterations (i.e., an epoch), some of the cuckoos that ended up in areas that offer low opportunities for improvement, migrate in the surrounding of the best
ED
and most promising areas so far discovered by them and or by other cuckoos. Operatively, this behavior has been implemented through the following steps:
Solutions are sorted in descending orders, in terms of the consecutive number of flights that did not
PT
produce any further improvement (i.e., solutions corresponding to cuckoos which are more likely to be stuck in a local minimum are placed at the top of the list); The first
cuckoos are preselected and, among these one, all those who exceeded a predefined
CE
threshold limit (expressed as number of consecutive failed flights) are finally selected for
AC
migration;
The starting points of migrating cuckoos are substituted with new points, generated as neighbors of the
5.3.
best solutions found so far.
CS-MNM full implementation
The implementation is based on two sets and its nest
and
storing, for each cuckoo , its current starting point
i.e., the best solution found so far by cuckoo . At each iteration, both sets are updated
accordingly to the outcome of the flight. Let us consider a generic cuckoo {
}, with
placed at point
. Also, let
, be the points visited by the cuckoo during its flight and let
be the best
among them. In accordance with the generation scheme of the set
(i.e., Double Shaped Lévy Flight),
ACCEPTED MANUSCRIPT anytime
proves to be a point of improvement, it replaces
following neighbor
is obtained. Also note that, as clarified in Sub Section 5.2.1, the flight is
immediately interrupted as soon as a coincides with both vectors
as the new starting point, from which the
improves the current best
and it is also better than
and
been unsuccessful,
. That is, whenever
. In this case the mission has been successful, and thus
are updated so that
. Conversely, if
is not updated and, typically, the cuckoo will go back to
changed during the flight (Mode #2 was employed exclusively), even if
, the mission has
. However, if
is worse than
has not
, the cuckoo might
and resume its next exploration from that point, with a certain probability (as explained in
Sub-Section 5.2.2). In this case
is updated so that
CR IP T
choose to go in
, the last point
.
For the sake of clarity, the full pseudocode of the CS-MNM is shown in Table IV. Briefly, the Double Shaped Lévy Flight is implemented from Line #14 to #24 and the part of code between Line #20 to #21 is responsible to switch between the two flight modes, as described in section 5.2.1. Cuckoo‟s Risk Aversion is
AN US
implemented by Line #25 to #28 and, as it can be seen, it is enabled anytime the Double Shaped Lévy Flight proved unsuccessful. Lastly, Clan Migration is detailed by Line #29 to #32. Table IV CS-MNM pseudocode
PT
ED
M
Set , Set P% Set , Set N Set Randomly generate the set N of host nests EN u 0N k0 Do While Or (Stop Criterion) Randomly select Cuckoo i: set , j0
14
Do While
CE
13
15
Risk Aversio n
16 17 18 19 20 21 22 23 24 25
From generate new neighbor If Then Break If Then
jj+1 Loop If
# Exploration set, impose initial condition E ≡ N # Vector of consecutive unsuccessful flights
accordingly
Accept
and store it in set # Improvement over # Replace both starting point and last host nest # Reset count of unsuccessful flights # Halt Double Shaped Lévy Flight # Improvement over # Replace only the starting point # Increase count of unsuccessful flights
# No improvements over ( )
26 27
# Levy Flight lower limit and Power Law Exponent # Epoch migration population ratio # Initial and Final Non-Return Probability # Number of iterations # Threshold of consecutive unsuccessful flights
Generate travel distance d
AC
Double Shaped Lévy Flight
1 2 3 4 6 7 8 5 9 10 11 12
accordingly to Eq. (13)
# Best of generated neighbors
Clan Migration
ACCEPTED MANUSCRIPT 28 29 30 31 32 33
AC
CE
PT
ED
M
AN US
CR IP T
34
If Accepted ( Then: # Accept a worsening exploration starting point If An Epoch has Elapsed Then B P% of elements of N, sorted in ascending order accordingly to W P% of elements of E, sorted in descending order accordingly to , s.t. Replace each vector solution of W with a neighbor of a vector solution of B kk+1 Loop
ACCEPTED MANUSCRIPT 6. PERFORMED TESTS All the tests discussed in this Section were performed on a test set of randomly generated problems characterized by the following features (Zammori et al., 2014):
Processing Times are Normally distributed with mean
and standard deviation
of 100 and 25
time-units, respectively;
Setup Times are Uniformly distributed in the interval [
Due Dates are Uniformly distributed in the interval [
], where
and
is the Tardiness Factor, representing the quota of jobs likely to be tardy,
is the relative Range of Due Dates (Tan et al. 2000);
The earliness and tardiness penalties have been set, respectively, at
Time-to-Failure follows a Weibull distribution with scale parameter ;
and
;
and shape parameter
The average time to execute Planned and Reactive Maintenance are set, respectively, to and
AN US
is the
CR IP T
number of jobs,
];
.
6.1. Parameters optimization
To assure a fair comparison of the metaheuristics introduced in Sections 3 to 5, their operating parameters were fine-tuned using a test set of 50 small size problems, characterized by 14-Jobs and 2-PMs. Problems of
M
this size were chosen because: (i) their complexity is high enough to prevent the metaheuristics from finding, almost always, the optimal solution and, (ii) they can be optimally solved in a reasonable amount of time
ED
(i.e., around one and a half hour on a computer with a 3.2 GHz processor and 16 GB of RAM memory), through an exhaustive search. Also, parameters that can be reasonably considered independent were finetuned one by one, using a sequential optimization strategy. Conversely, parameters for which interaction
More precisely:
MHS was not considered in the fine-tuned process, as we used the parameters (reported in Table V)
CE
PT
effects cannot be excluded a priori, were pre-analyzed and optimized using Design-Of-Experiment (DOE).
suggested by Zammori et al. (2014), which had been already optimized on the same problem set; SA and CS-LF (i.e., the original metaheuristics) were optimized first. Parameters were assumed
AC
independent and fine-tuned following the sequential approach;
HS/4/IJMS and CS-MNM were optimized next. In this case, parameters shared with the original metaheuristics were directly inherited from the previous optimization step; conversely, specific parameters were optimized using DOE, to better investigate possible interactions among them.
6.1.1. Sequential optimization Let us consider a generic algorithm , and let
{
} be the set of its parameters, ordered
from the most to the least important one. Sequential optimization proceeds as follows:
Initialize each parameter pi using a reasonable starting value. To this aim we used values suggested in the literature;
ACCEPTED MANUSCRIPT
Given these initial settings, all problems of the test set are solved and the accuracy of
is
computed. To this aim, we measured accuracy in terms of APE, as in Eq. (8), using 8 repetitions of each problem instance;
The parameter
, having the higher influence on the algorithm performance, is selected and its
value is slightly modified. All the other parameters {
} are kept unaltered and all the
problems of the test set are solved again. This process is iteratively repeated, proceeding at small step, until the accuracy of
continues to grow; is frozen and the second parameter
is
CR IP T
When accuracy does not increase anymore, the value of
selected for possible improvements and it is fine-tuned using the same approach described above;
The procedure is repeated for all the other parameters of the set .
Obtained results are shown in Table V, where optimized parameters are listed in order of importance. The optimal value is shown in the last column, whereas the range and the step size used in the optimization are
AN US
reported in the third and fourth column, respectively. Table V Sequential optimization results Algorithm
Parameter Random Selection Rate Pitch Adjustment Rate Harmony Memory Size Start Probability End Probability Minimum Flight dmin Pareto Shape parameter Purge percentage P% Population Size
Step Size
Optimal Value
[0.05; 0.50] [0.001; 0.03] [1; 10] [1; 2] [0.01; 0.25] [0.5NJ; 2NJ]
0.05 0.005 1 0.3 0.05 0.25
0.05 0.1 0.5NJ 0.1 0.02 7 1.9 0.1 NJ
M
MHS
Range
ED
SA
CS-LF
PT
Where NJ is the total number of Jobs
6.1.2. Design of Experiment
CE
Given the higher complexity of HS/4/IJMS and CS-MNM, a 2-levels full factorial design was made (using the levels shows in Table VI), limited to the following parameters: Concerning HS/4/IJMS, we considered
AC
,
and
, because these parameters
conjointly govern the rate at which the Large Portion Recovery is triggered;
Concerning CS-MNM we considered
,
, and the Epoch‟s Length, because these parameters
conjointly influence the evolution of the solutions contained within sets
and .
More precisely, the 10 hardest problem instances (among the 50 previously generated ones) were selected and each one was solved 8 times, both with HS/4/IJMS and with CS-MNM. Testing
combinations of
parameters levels, a total of 640 solutions (i.e., 10·8·23) were generated with each metaheuristic. It is worth noting that, to assess interaction effects through 3-ways ANOVA, there is the need to collect more than a
ACCEPTED MANUSCRIPT single value of the response variable for each one of the 10·23 combinations of problems and parameters setting. Thus, we replaced APE with the Relative Percentage Error defined in Eq. (16): (16) Where, attempt
[
is the (a priori known) optimal solution of problem [
] and
] by the used metaheuristic, with parameters configuration
[
is the solution found at ].
Operating in this way, two remarkable results were found, as DOE and the related ANOVA tests revealed
CR IP T
that:
Concerning possible interaction effects, all p-values were higher than 0.25 (definitely greater than the usual threshold of 0.05), and so parameters can be considered independent;
Concerning the effect of the operating parameters on the algorithm performance, all p-values were higher than 0.1, and so performances are stable and little influenced by specific values of the
AN US
operating parameters.
Since both results confirmed the robustness of the proposed metaheuristics, as well as the absence of any interactions among parameters, we completed the optimization process following the same sequential approach detailed above. Obtained results are shown in Table VI.
CS-MNM
Levels
Optimal Value
Low saturation threshold SLow
[0.1; 0.3]
0.3
High saturation threshold SHigh Large Portion Recovery Rate LPRRLow Starting Risk Aversion SRA Ending Risk Aversion ERA Epoch E
[0.5; 0.7]
0.7
[0.05; 0.1] [0.1; 0.5] [0.01; 0.05] [50; 100]
0.05 0.1 0.01 75
PT
HS/4/IJMS
Parameter
ED
Algorithm
M
Table VI DOE results
CE
6.2. Numerical Assessment
Lastly, SA, CS-LF, HS/4/IJMS and CS-MNM were compared in relative terms. Please note that the
AC
original MHS was not included in this analysis, because, as detailed in Section 3.6, it is outperformed by the SA, even with problems of relatively small dimensions. We also note that, to get a preliminary idea of the quality of the modifications made on the MHS, we tested its improved version (i.e., the HS/4/IJMS) against the same test set of 10 problems with 14-Jobs and 2-PMs. In this case, results were encouraging as HS/4/IJMS behaved almost like the SA and it consistently found the global optimum (i.e., APE = 0 and APD = 0.01). Comparisons were made relatively to two test sets with a different level of complexity: (i) 10 hard problems with 25-Jobs and 3-PMs, and (ii) 10 very hard problems with 50-Jobs and 8-PMs. In both cases, the number of PMs coincides with the optimal value, calculated accordingly to the procedure proposed by Zammori et al. (2014), to which the interested reader is referred for further insights.
ACCEPTED MANUSCRIPT Each problem instance was solved 8 times, by each metaheuristic, starting from a different initial population. Also, we let each metaheuristic run until convergence, so that computation was halted when the best solution did not improve for a sufficiently high number of consecutive iterations. Lastly, performances were assessed using APE and APD, as defined in Eq. (8) and Eq. (9), respectively. Please note that, in this case, the real optimum were not obtained through an exhaustive search; thus, in both equation the parameter substituted by the best solution found for problem , i.e.
{
is
}.
Table VII Final comparison 25-Jobs 3-PMs
CR IP T
Obtained results are reported in Table VII.
50-Jobs 8-PMs
APE
APD
Max Error
# Wins
APE
APD
HS/4/IJMS
0.007
0.023
0.0491
0%
0.0580
0.2141
0.3562
0%
SA
0.001
0.011
0.0400
20%
0.0030
0.0213
0.0647
40%
CS-LF
0.087
0.104
0.1416
0%
[NC]
[NC]
[NC]
[NC]
CS-MNM
0
0.004
0.0109
40%
0.0008
0.0074
0.0187
60%
AN US
Algorithm
Max Error
# Wins
As expected, the SA works very well and, even with 25-Jobs and 3-PMs it outperforms both CS-LF and HS/4/IJMS. Specifically, whereas performance of HS/4/IJMS are still acceptable, in its standard form the
M
Cuckoo Search is definitely the worst metaheuristics, both in terms of solutions‟ quality and stability. For these reasons CS-LF was not even considered in the second test set (as denoted by the [NC] values in
ED
columns 6 to 9 of Table VII). Frankly speaking, the very poor performance of the CS-LF was not expected. Probably, although solutions evolve independently and cross influence among them is minimized, the solution space is explored neither efficiently nor effectively. Alone, the Lévy Flight does not seem to be
PT
enough to assure the right random behavior required for a thorough exploration of the solution space. Also, purging the worst solutions (at every iteration) may limit the neighborhood search, thus further decreasing
CE
the overall effectiveness of the approach. Conversely, in its improved version, the CS-MNM turned out to be the best approach. CS-MNM, indeed, shows the minimum values of APE, APD and Max Error; particularly, a value of APE equal to zero ascertains that, for each problem instance, at least one of the eight solutions
AC
generated by the CS-MNM was the overall best one (relatively to the 32 generated by all the tested metaheuristics). CS-MNM also proved to be very precise in founding the best solutions as certified by a very low value of both APD and Max Error and by the highest value of percentage wins. Relatively to this last metric, it is important to note that summing the percentage wins of column 5 a total of 60% is obtained. This means that 40% of the times SA and CS-MNM converged to the same, presumably optimal, solution. The superiority of the CS-MNM is even better outlined by the outcomes of the very hard problems‟ set. Not only CS-MNM is the metaheuristic that located the best solution with the highest frequency but, as certified by the lowest APD value, it proved to be extremely robust because, even when it missed the best solution, it still found a solution very close to it. All these outcomes indicate that the introduced modifications (i.e., nests location memory, cuckoo‟s risk aversion, clan migration and double shaped Lévy
ACCEPTED MANUSCRIPT Flight), made it possible to overcome the drawbacks of the standard approach, assuring an efficient and effective exploration of the solution space. As a final note, we want to highlight the following fact. Since convergence was used as stop criterion, one could argue that comparison between CS-MNM and SA was not fair. Indeed, being a population-based metaheuristic, CS-MNM might have found better solutions only because it ran for a total number of iterations higher than that of the SA. More precisely, given a fixed number of iterations , SA evolves the same solution exactly
of the set , at iteration each solution will have been evolved
times, on average. It is thus
CR IP T
among the
times; conversely, considering that at each iteration CS-MNM selects a solution
plausible to expect that, to find the best solution, CS-MNM needs to run for a number of iterations
times
bigger than that of the SA. To ascertain this fact, we compared the evolution path of SA with that of CSMNM, for each generated solution. Quite surprisingly, the test showed that 75% of times, when SA reached convergence, CS-MNM had already found a better solution and its evolution process was still on progress; a fact that further certifies the superiority of the CS-MNM, which is also faster than the SA in finding the
AN US
optimal solution. This result is graphically displayed by Figure 5.a and Figure 5.b that show an example of the evolution path followed by the two metaheuristics.
------------------------------Insert Figure 5.a here
--------------------------------
M
------------------------------Insert Figure 5.b here
6.3.
ED
--------------------------------
Operational Implications
PT
Once the superiority of the CSM-MNM algorithm was verified, a last check was made to quantify the operational advantage that could be achieved using the CSM-MNM algorithm to schedule jobs and PMs in
CE
an integrated way. To this aim we compared the cost of the best schedule found by the CSM-MNM algorithm with the cost of the presumably best schedule that could be generated by a production manager
AC
equipped with a standard scheduling system. The latter one was generated using the following naïve algorithm:
At first only jobs are considered (i.e., the single machine is considered as fault free) and the best schedule S is generated. This schedule, that using Kendall‟s notation can be classified as (1| | Earliness-Tardiness Penalties), corresponds to the one that could be generated by the standard scheduling system and, in the present case is obtained evaluating all the (n!) permutations of the jobs‟ sequence. Next, the optimal number M of PMs is inserted in S, without altering the ordering of the jobs. This is made in the following reasonable way: - S is split into (M +1) parts of almost equal duration;
ACCEPTED MANUSCRIPT - PMs are inserted at these subdivision points, to assure the same level of availability in each subsection of the schedule. The comparison has been made testing the exhaustive solver and the CS-MNM against 5 different instances of the IJMS problem characterized by 50-Jobs and 8-PMs. The CS-MNM always outperformed the naïve exhaustive, and solutions generated by the CS-MNM were on average (23 ± 1.1) % more efficient, in terms of expected operating cost, than that obtained by the naïve exhaustive. These results stress even more the
CR IP T
relevancy and applicability of the proposed approach.
7. CONCLUSIONS AND FUTURE WORKS
In this paper, we focused on an integrated jobs-maintenance scheduling problem, relevant to many industrial realities, where failures and/or reduced availability may jeopardize productivity. Specifically, we considered a single machine subjected to failures and aging with the objective to minimize earliness-
AN US
tardiness penalties. These features make the problem particularly challenging, because the solution space becomes very flat, but disseminated of many local minima.
Considering the practical relevance of the problem and its inherent complexity, we developed and tested different metaheuristics, aiming to find robust and quasi-optimal schedules in a reasonable amount of time. Specifically, we benchmarked the famous SA by Ishibuchi et al. (1995), with two population-based
M
metaheuristics, namely the Modified Harmony Search by Zammori et al. (2014) and the Cuckoo Search by Lévy Flight by Yang and Deb (2009). MHS was chosen because it was specifically designed to tackle this problem, CS-LF was chosen as it implements a very effective local minimum escaping strategy. In addition,
ED
we also developed two alternative approaches: (i) HS/4/IJMS, a hybrid Harmony Search and Genetic Algorithm, specifically tied to this problem and (ii) CS-MNM, an innovating modification of the Cuckoo
PT
Search via Lévy Flight.
After a fine-tuning procedure, a comprehensive comparison was made. Outcomes showed that, in these
CE
peculiar conditions, simpler heuristics characterized by a high degree of randomness may perform better than sophisticated population-based metaheuristics. Indeed SA (i.e., the simpler algorithm) outperformed MHS, CS-LF and HS/4/IJMS, especially in case of very hard problems, characterized by 50-Jobs and 8 PMs. Only
AC
the newly developed CS-MNM subverted this situation, as it was the best approach both in term of precision and accuracy of the obtained solutions. We also showed that CS-MNS offers the following operating advantages:
It has relatively few parameters to optimize; Parameters are independent (there are no interaction effects among them) and can be optimized independently; It performs an excellent exploration of the solution space but, at the same time, it assures a fast convergence and a low computation time, also for large problems.
All these features make the CS-MNS algorithm a good candidate for a possible inclusion in the optimizer module of a Cyber Physical System, fully compliant with the principles of the Industry 4.0 paradigm.
ACCEPTED MANUSCRIPT Frankly speaking we must admit that, at present, our model requires the number of PMs that must be inserted in the schedule as input; although this value can be analytically estimated, a better approach would be to use a two-level optimization to generate, not only the optimal job-maintenance sequence, but also the optimal number of PMs. Another limit of the model is that a single level of planned maintenance is considered because, anytime PMs are executed, the hazard rate of the machine is brought back to its original level (i.e., perfect repair model). However, a more realistic scenario could include different levels of maintenance tasks, each one characterized by different costs and times and capable to reduce (without
CR IP T
necessarily reset) the aging and/or the degradation of the machine. Lastly, we note that, although the single machine scheduling problem is a rather flexible model, as it can be applied to a variety of common production systems, it could also be extended to other common layouts, such as parallel machines, flow shops or even job shops. All these issues could be topics for future researches.
AN US
Authors contributions section
F.Zammori and D.Mezzogori conceived of presented idea and development of the CS-MNM algorithm. M.Bertolini implemented the rest of the algorithms used as benchmark. D.Mezzogori executed the Design of Experiment section and F.Zammori performed final numerical comparison. M.Bertolini finally supervised
M
and checked the findings of this work. All authors discussed the results and contributed to the final
ED
manuscript.
PT
8. REFERENCES
AC
CE
Allahverdi, A. (2015). The third comprehensive survey on scheduling problems with setup times/costs, European Journal of Operational Research, 377 (2), 345-378. El-Khoukhi, F., Boukachour, J., & Alaoui, A.L.H. (2017). The „„Dual-Ants Colony”: A novel hybrid approach for the flexible job shop scheduling problem with preventive maintenance, Computers & Industrial Engineering, 106, 236– 255. Cui, W., Lu, Z., & Pan, E. (2014). Integrated production scheduling and maintenance policy for robustness in a single machine, Computers & Operations Research, 47, 81-91. Cui, W., & Lu, Z. (2017). Minimizing the makespan on a single machine with flexible maintenances and jobs‟ release dates, Computers & Operations Research, 80, 11-22. Cui, W., Lu, Z., Li, C., & Han, X. (2018). A proactive approach to solve integrated production scheduling and maintenance planning problem in flow shops, Computers & Industrial Engineering, 115, 342–353. Geem, Z., Kim, J., & Loganathan, G. (2001). A new heuristic optimization algorithm: harmony search, Simulation, 76 (2), 60-68. Gen, M., & Cheng, R. (2000). Genetic algorithms and engineering optimization. John Wiley & Sons, New York, USA. Ishibuchi, H., Misaki, S., & Tanaka, H. (1995). Modified simulated annealing algorithms for the flow shop sequencing problem, European Journal of Operational Research, 81 (2), 388–398. Kirkpatrick, S., Gelatt, C., & Vecchi, M. (1983). Optimization by simulated annealing, Science, 220 (4598), 671-680. Komaki, G., Teymourian, E., Kayvanfar, V., & Booyavi, Z. (2017). Improved discrete cuckoo optimization algorithm for the three-stage assembly flowshop scheduling problem, Computers & Industrial Engineering, 105, 158–173. Lee, J., Bagheri, B., & Kao, H. A. (2015). A cyber-physical systems architecture for industry 4.0-based manufacturing systems. Manufacturing Letters, 3, 18-23. Lee W.C., & Wang J.Y. (2017). A three-agent scheduling problem for minimizing the makespan on a single machine. Computers & Industrial Engineering, 106, 147-160.
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
AN US
CR IP T
Low C., Hsu C., & Su C., (2010) A modified particle swarm optimization algorithm for a single-machine scheduling problem with periodic maintenance. Expert System with Applications, 37, 6429–6434. Lu Z., Cui, W., & Han X., (2015). Integrated production and preventive maintenance scheduling for a single machine with failure uncertainty. Computers & Industrial Engineering, 80, 236–244. Luo, W., Cheng, T.C.E., Ji, M., (2015). Single-machine scheduling with a variable maintenance activity. Computers & Industrial Engineering, 79, 168 – 174. Macchi M., Roda I., Fumagalli L., 2017, 'On the Advancement of Maintenance Management Towards Smart Maintenance in Manufacturing', IFIP Advances in Information and Communication Technology, 383-390 Metropolis, N., Rosenbluth, A., Rosenbluth, M., Teller, A., & Teller, E. (1953). Equation of State Calculations by Fast Computing Machines. The Journal of Chemical Physics, 21 (6), 1087-1092. Mokhtari, H., Dadgar, M., (2015). Scheduling optimization of a stochastic flexible job-shop system with time-varying machine failure rate, Computers & Operations Research, 61, 31-45. Mor, B., & Mosheiov, G., (2015). Scheduling a deteriorating maintenance activity and due-window assignment. Computers & Operations Research, 57, 33-40. Mosheiov, G., & Sarig, A., (2009). Scheduling a maintenance activity and due-window assignment on a single machine. Computers & Operations Research, 36, 2541-2545 Nouiri M., Bekrar, A., Jemai, A., Trentesaux, D., Ammari, C.A., & Niar, S., (2017). Two stage particle swarm optimization to solve the flexible job shop predictive scheduling problem considering possible machine breakdowns, Computers & Industrial Engineering , 112, 595-606. Ouaarab, A., Ahiod, B., & Yang, X.Y., (2014). Discrete cuckoo search algorithm for the travelling salesman problem. Neural Computing & Application, 24, 1659–1669 Quan, O. Y., & Xu, H. Y. (2013). The review of the single machine scheduling problem and its solving methods. In Applied Mechanics and Materials (Vol. 411, pp. 2081-2084). Trans Tech Publications. Rahmati, S.H.A., Ahmadi, A., & Govindan, K., (2017). A novel integrated condition-based maintenance and stochastic flexible job shop scheduling problem: simulation-based optimization approach, Annals of Operations Research, 139. Safari, E., & Sadjadi S., (2011). A hybrid method for flowshops scheduling with condition-based maintenance constraint and machines breakdown. Expert System with Applications, 38, 2020-2029. Salmasnia, A., & Mirabadi-Dastjerd, D., (2017). Joint production and preventive maintenance scheduling for a single degraded machine by considering machine failures. TOP, 25, 544–578. Soltani, R., Jolai, F., & Zandieh, M. (2010). Two robust meta-heuristics for scheduling multiple job classes on a single machine with multiple criteria. Expert Systems with Applications, 37 (8), 5951–5959. Tan, K. C. T., Narasimhan, R., Rubin, R. A., Ragatz, & G. L., (2000). A comparison of four methods for minimizing total tardiness on a single processor with sequence dependent setup times, Omega, 28, 313-326. Touat, M., Bouzidi-Hassini, S., Benbouzid-Sitayeb, F., & Benhamou, B. (2017). A hybridization of genetic algorithms and fuzzy logic for the single-machine scheduling with flexible maintenance problem under human resource constraints, Applied Soft Computing, 59, 556-573. Woo, Y., Kim, B., (2018). Methaeuristic approaches for parallel machine scheduling problem with time-dependent deterioration and multiple rate-modifying activities, Computers & Operations Research, 95, 97-112. Yang, X.-S. (2010). Nature-inspired metaheuristic algorithms 2nd Ed., Luniver Press, Frome, UK. Yang, X.-S., & Deb, S. (2009). Cuckoo search via Lévy flights, Proceeding of World Congress on Nature & Biologically Inspired Computing, December 2009, India, IEEE Publications, USA, 210-214. Yazdani, M., Khalili, S.M., Babagolzadeh,M., & Jolai, F., (2017). A single machine scheduling problem with multiple unavailability constraints: A mathematical model and an enhanced variable neighborhood search approach, Journal of Computational Design and Engineering, 4, 46–59. Yu, A.J., & Seif J.,(2016). Minimizing tardiness and maintenance costs in flow shop scheduling by a lower-boundbased GA. Computers & Industrial Engineering, 97, 26 – 40. Zammori, F., Braglia, M., & Castellano, D. (2014). Harmony search algorithm for single-machine scheduling problem with planned maintenance, Computers & Industrial Engineering, 76, 333–346. Zhu, H., Li, M., Zhou, J. (2015). Machine scheduling with deteriorating and resource-dependent maintenance activity, Computers & Industrial Engineering, 88, 479-486.
ACCEPTED MANUSCRIPT Appendix A. Simulated Annealing (SA) Introduced by Metropolis et al. (1953) and made popular by Kirkpatrick et al. (1983), SA is a probabilistic heuristic, whose inspiration come from annealing in metallurgy, a slow cooling process used to bring melted metals to a uniform, low-energy crystallization state. Briefly, at each one of neighborhood of the current solution
iterations, SA generates a new solution
from the
and accepts it as the new starting point, with a probability that, as iterations go
on, decreases following the negative exponential law show in Eq. (a.1), valid for a minimization problem.
Where
(
)}
(a.1)
CR IP T
{
is the change in the objective function and
, the so-called temperature at iteration , is a
parameter that accounts for the number of completed iterations. To this aim, the temperature is progressively and linearly reduced from an initial value
, as in Eq. (a.2): (
)
AN US
For the sake of clarity, the full SA Pseudo-Code is reported in Table A.I
(a.2)
Table A.I SA by Ishibuchi et al. (1995)
# Number of Iterations, Initial and Ending Temp. # Number of Solutions explored for each level of T
# Current and Best Solution # Best and Current Objective Function
ED
M
Set N, T0, TN Set K Randomly generate an initial solution y x, x* y x, x* f(y) i0 Do While T0 < TN Generate K new neighbors of x {y1,…, yK}
# Best OF, among the newly generated solutions # Best of the K generated solutions
PT
y Accept accordingly to Eq. (a.1) If Accepted ( Then: x y, x y If y < x* Then: x* y, x* y ii+K Update Ti as in Eq. (a.2) Loop Return x*, x*
# Temperature is linearly reduced
CE
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
AC
Also, neighbor solutions are generated using the pairwise interchange operator i.e., two locations and are randomly selected and the corresponding elements
and
are swapped:
and
Concerning the initial and ending temperatures, their values are set accordingly to the classical approach by Kirkpatrick et al. (1983). Specifically,
is chosen to assure a desired (high) Initial Acceptance Probability (IAP) computed as in
Eq. (a.3), where ̅̅̅̅ is the average variation in the objective function when a neighbor is generated. ̅̅̅̅
Please note that Eq. (a.3) can be used to compute Acceptance Probability (EAP).
(a.3)
too, by simply substituting IAP with the desired (low) Ending
ACCEPTED MANUSCRIPT
# Novelty
AN US
CR IP T
Set # Harmony Memory Consideration Rate Set # Pitch Adjustment Rate Set # Block Shifting Rate Set # Number of iterations Set # Number of job and maintenance tasks Randomly generate the Do While Or (Stop Criterion) Compute accordingly to Eq. (10) If Then Generate empty solution For To Do If Then: Get an element via Random Selection (RS) Else Get an element via Harmony Memory Consideration (HMC) # Add element to candidate solution Next If Then: Perform Pitch Adjustment (PA) on Else If Then: Perform Block Shifting (BS) on
# Novelty
# Novelty
ED
M
If Then: # Overwrite worst solution Else Let be the number of copies within unique solutions within Initialize set of offsprings Do While : Select two solutions from using Wheel Selection Rule Perform Order Crossover (OC) on the selected solutions Store generated solutions into For Each Offspring In Do If Then Perform Pitch Adjustment (PA) on the Offspring Replace a copy within the with the Offspring Loop
PT
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34
# Novelty
CE
Large Portion Recovery
HM Consideration
APPENDIX B. HS/4/IJMS Pseudo-Code
All modifications (made to the original MHS algorithm) are clearly defined in the HS/4/IJMS pseudocode.
AC
Specifically:
Line #9 is used to dynamically compute the value of LPRR, based on the saturation of the Harmony Memory, which will trigger either the original Harmony Memory Consideration or the Large Portion Recovery;
Line #18 introduces the Block Shifting operator;
Lines from #21 to #32 specify how the Large Portion Recovery works. Specifically, lines #29 to #32aim to minimize the saturation of the Harmony Memory, by replacing all repeated solutions.
Note that standard HS corresponds to lines #10 and line #20, except for line #18.
ACCEPTED MANUSCRIPT Fig. 1.a – Neighbour solutions of three instances of a 6-Jobs 1-Planned Maintenance problem Fig. 1.b – Sorted solutions of three instances of a 6-Jobs, 1-Maintenance problem Fig. 2 – Typical trend of the Saturation Index during a computational run Fig. 3. – Typical trend of the objective function during a computational run Fig. 4. – Cuckoo‟s flight modes Fig. 5.a – Typical evolution in the initial phases Fig. 5.b – Typical evolution before convergence is reached
AC
CE
PT
ED
M
AN US
CR IP T
1. 2. 3. 4. 5. 6. 7.
AN US
CR IP T
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
Fig. 1.a – Neighbour solutions of three instances of a 6-Jobs 1-Planned Maintenance problem
CR IP T
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
AN US
Fig. 1.b – Sorted solutions of three instances of a 6-Jobs, 1-Maintenance problem
CR IP T
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
AN US
Fig. 2 – Typical trend of the Saturation Index during a computational run
CR IP T
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
AN US
Fig. 3. – Typical trend of the objective function during a computational run
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
AN US
CR IP T
Fig. 4. – Cuckoo‟s flight modes
CR IP T
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
AN US
Fig. 5.a – Typical evolution in the initial phases
CR IP T
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
AN US
Fig. 5.b – Typical evolution before convergence is reached