Joint optimization of production and maintenance scheduling for unrelated parallel machine using hybrid discrete spider monkey optimization algorithm

This paper considers an unrelated parallel machine scheduling problem with variable maintenance based on machine reliability to minimize the maximum completion time. To obtain the optimal solution of small-scale problems, we firstly establish a mixed integer programming model. To solve the medium and large-scale problems efficiently and effectively, we develop a hybrid discrete spider monkey optimization algorithm (HDSMO), which combines discrete spider monkey optimization (DSMO) with genetic algorithm (GA). A few additional features are embedded in the HDSMO: a three-phase constructive heuristic is proposed to generate better initial solution, and an individual updating method considering the inertia weight is used to balance the exploration and exploitation capabilities. Moreover, a problem-oriented neighborhood search method is designed to improve the search efficiency. Experiments are conducted on a set of randomly generated instances. The performance of the proposed HDSMO algorithm is investigated and compared with that of other existing algorithms. The detailed results show that the proposed HDSMO algorithm can obtain significantly better solutions than the DSMO and GA algorithms.


Introduction
The production scheduling problem encompasses various branches, and one of the significant ones is the parallel machine scheduling problems (PMSPs). In the literature, PMSPs are typically classified to three groups: identical parallel machines (Pm), uniform parallel machines (Qm) and unrelated parallel machines (Rm) (Cheng & Sin, 1990). The category of unrelated PMSPs (UPMSPs) encompasses a broader scope than the other two categories, as it involves various machines carrying out identical tasks but with varying processing abilities or capacities. With the development of intelligent manufacturing industries, such as computerized numerical control machine tools, industrial robots, unrelated parallel machine production has become the most common operating mode of enterprises. Solving real-life UPMSPs is a major challenge for industrial experts and researchers, not only because they are mostly NP-hard but also, more importantly, because of the special characteristics or requirements they have in practice.
It is noted that all the above studies assume that machines are continuously available over the scheduling horizon. However, in a real manufacturing environment, machines may not be available because of maintenance operations , tool replacement (Dang et al., 2021), machine breakdowns (Kim & Kim, 2020), etc. UPMSPs with maintenance have been extensively studied in recent years for various maintenance activities and performance measures. For the UPMSP with at most one maintenance activity on each machine, Cheng et al. (2011) proved that the two problems of minimizing the total completion time and the total machine load can be optimally solved in polynomial time, Hsu et al. (2013) proposed models for the same problem as in reference (Cheng et al.,2011) considering three basic types of ageing effects and proved that the models could be solved optimally in polynomial time, and Lu et al. (2018) proposed an ABC-TS algorithm to find an approximately optimal solution in a reasonable time. For the UPMSP considering multiple maintenance activities and ageing effects simultaneously, Yang et al. (2012) applied the group balance principle to determine the optimal positions of the maintenance activities and the number of jobs in each group in the scheduling sequence on each machine. Tavana et al. (2015) presented an integrated three-stage model with the fuzzy analytic hierarchy process, the technique for order of preference by similarity to ideal solution, and goal programming. Gara-Ali et al. (2016) proposed a general model for the UPMSP with different maintenance systems and several performance criteria. For the UPMSP with multiple maintenance activities and sequence-dependent setup times, Avalos-Rosales et al. (2018) proposed a mathematical formulation with valid inequalities to obtain optimal solutions for small to medium instances and an efficient meta-heuristic algorithm based on a multistart strategy for solving larger instances. Lei and Yang (2022) proposed a multisubcolony ABC algorithm to simultaneously minimize the makespan and total tardiness. Lei and Yi (2021) presented a differentiated shuffled frog leaping algorithm and used its strong exploration ability to minimize the makespan. Wang & Pan (2019) presented a novel ICA with an estimation of distribution algorithm to simultaneously minimize the makespan and total tardiness. In addition, the distributed UPMSP with preventive maintenance (Lei & Liu, 2020), UPMSP with additional resources and maintenance (Lei & He, 2022) and UPMSP with release times and maintenance activities (Pang, Tsai, & Chou, 2021) have been researched.
Regarding the above works on the UPMSP with maintenance, three types of maintenance activities are generally assumed: periodically fixed, flexible and variable. For scheduling with fixed maintenance activities, the starting time and duration of each maintenance activity are predefined or given beforehand (Avalos-Rosales et al., 2018;Lei & Yang, 2022;Lei & Yi, 2021;Wang & Pan, 2019;Lei & He, 2022). For scheduling with flexible maintenance activities, the maintenance operation must be performed within a preplanned time window (Lei & Liu, 2020;Beldar et al., 2022) or below a specific threshold (Pang, Tsai, & Chou, 2021), and the duration is fixed (Pang, Tsai, & Chou, 2021) or related to the starting time of the maintenance (Beldar et al., 2022). For scheduling with variable or deteriorating maintenance activities, the maintenance starting times are treated as decision variables, and the maintenance duration is assumed to increase with the starting time (Cheng, Hsu, & Yang, 2011;Hsu et al., 2013;Lu et al., 2018) or is fixed (Yang et al., 2012;Tavana et al. 2015). However, in a real production environment, the starting time of a maintenance activity sometimes depends only on the reliability of the machine. To the best of our knowledge, scheduling with this kind of maintenance is very rare in the literature. Therefore, we considered the UPMSP with reliability-based maintenance in this paper, which assumes that the machine's status follows a discrete degradation process, and if the machine's reliability is less than the minimum acceptable level or threshold before starting job processing, the maintenance operation must be performed.
Additionally, a large number of meta-heuristics based on swarm intelligence, such as the GA (Lin, Pfund, & Fowler, 2011) , SA (Wang et al., 2020;Lin & Ying, 2015), ABC (Lei, Yuan, & Cai, 2021;Lei & Yang, 2022;Lei & Liu, 2020;Lei & He, 2022), TS (Wang & Alidaee, 2019), ICA , and hybrid algorithm ABC-TS (Lu et al., 2018), have been developed to deal with UPMSPs with maintenance. However, the application of meta-heuristics such as spider monkey optimization (SMO) has not been fully investigated. The SMO algorithm proposed by Bansal et al. (2014) is a swarm intelligence optimization algorithm inspired by the intelligent foraging behaviour of fission-fusion social structure-based animals. Since SMO has ability to trade-off between exploration and exploitation, SMO and its variants have been widely used to solve complex real-world optimization problems; these variants include numerical optimization and continuous constrained optimization (Sharma et al., 2016;Gupta et al., 2017). Cheruku et al. (2017) proposed a SMO-based rule miner (SM-RuleMiner) by incorporating a unique fitness function based on diabetes classification. In 2017, Sharma et al. suggested a method for determining the ideal placement and size of capacitors using a combination of a SMO approach based on a limaçon curve and a local search strategy inspired by the same curve. In recent years, SMO has been successfully applied to solve discrete optimization problems. Mumtaz et al. (2020) proposed a hybrid SMO algorithm for multilevel planning and scheduling problems of assembly lines. Yue et al., (2023) proposed a hybrid Pareto spider monkey optimisation algorithm for a two-stage flexible printed circuit board flow shop to minimize TWC and energy consumption simultaneously. Xia et al., (2021) introduced discrete SMO as a solution method for a vehicle routing problem involving uncertain demands. Their findings demonstrate that SMO exhibits strong global search capabilities. This paper presents the HDSMO, a hybrid DSMO algorithm that employs a combination of discrete spider monkey optimization and GA techniques to effectively tackle the UPMSP problem with variable maintenance. To generate feasible initial solution, a three-phase constructive heuristic is proposed. To balance the exploration and exploitation capabilities of DSMO, an individual updating method considering the inertia weight is used. To enhance the search efficiency, a problemoriented neighbourhood search method is applied. Experiments are conducted to compare the performance based on computational time and solution quality of HDSMO, DSMO, and GA on three different scales of the problem.
The research contributions of this paper can be summarized as follows. A DSMO algorithm composed of the spider monkey algorithm and GA is proposed. A hybrid DSMO algorithm is proposed with an initial solution generation method, discrete individual updating method and neighbourhood search method. In addition, the proposed methods have been evaluated and compared with the SMO and GA algorithms in a set of instances.
The subsequent sections of this paper are structured as follows. In Section 2, the UPMSP with maintenance considered in this study is presented, along with its mathematical model. The DSMO, and a proposed HDSMO that includes initial population generation, individual update, and neighbourhood search are presented in Section 3. The experimental results are presented in Section 4, where the proposed algorithm's validity is analyzed. In Section 5, conclusion and future research directions are provided.

Problem description
A set = { , ⋯ , ⋯ , } of jobs is to be processed on unrelated parallel machines , = 1, ⋯ , m. Let denote the number of jobs assigned to , and let ∑ = n. We assume, as in most practical situations, that < n. The jobs are all available for processing at time zero. In this paper, we assume that a machine's reliability follows an exponential distribution, and we let L represent the cumulative processing time of the job being continuously processed or the age of the machine. Then, the reliability of the machine is equal to , and is the machine failure rate. If the reliability falls below the threshold before starting a job's processing, a maintenance operation must be performed to restore the machine to its original condition. Because both the frequency and location of the maintenance activities are decision variables, we define this maintenance activity as the variable maintenance like reference (Beldar et al., 2022). Using the three-field notation α| | introduced by Graham et al. (1979), we denote our problem by Rm/nr, VM/C , where denotes that the jobs are nonresumable, VM denotes variable maintenance, and the objective is to minimize the maximum completion time. The decision is to determine the allocation and sequence of jobs on machines and the maintenance activity arrangements. Since problem Rm/nr, VM/ is NP-hard (Lu et al., 2018), approximate methods are needed to solve real-size instances.
To obtain the near-optimal solution of the studied problem, two principles need to be followed. One is to allocate job to the machine with the smallest processing time as often as possible, where the total processing times of the jobs allocated to the machines should be as close as possible. The second is to sequence the job on machine to minimize the number of maintenance operations. Assume that all jobs between two maintenance operations are grouped into a batch; if the machine's reliability is equal to the threshold before starting the last job in the batch, we define the batch as fully loaded. The decision regarding job sequencing in machine is equivalent to the job-grouping batch decision.

Property:
The batch is almost fully loaded, and the longer the processing time of the last job, the less maintenance is needed, meaning that a better solution can be obtained.
Proof: Suppose batch on machine , consisting of n jobs, can be divided into two sub-batches, job and the remaining n − 1 jobs. Assume that the processing time of job is = max , ∈ and the total processing time of the n − 1 jobs is = ∑ , where > and > . If job is processed after the other n − 1 jobs, the completion time of batch is = ∑ , as shown in Fig. 1(a). If we swap job with any job before it, the completion time of batch is likely to be = ∑ + * , and the number of maintenance operations is 0, as shown in Fig.1(b). Because the continuous processing time of the first n − 1 jobs becomes longer after the swap, the machine reliability before starting job number n in the batch is probably lower than the threshold , and maintenance is needed. Thus, − 0. (1) , ∀ = 1,2, … , = 2,3, … , , ∀ = 1,2, … , In this model, the objective is to minimize , as shown in Eq. (1). Constraints (2) and (3) ensure that each job can only be processed at one position of one machine and that each position of one machine can only be occupied by one job. Constraint (4) specifies the start time of each machine at the first position. Constraint (5) specifies the age of the machine at the first position before processing the job. Constraint (6) specifies the completion time of each machine at the first position. The start and completion times for each machine at each position are defined by constraints (7) and (8). Constraints (9) and (10) define the age of the machine. Constraints (11) and (12) define the reliability function and threshold of the machine. Constraint (13) defines the maximum completion time of each machine.

Encoding and decoding of individuals
According to the definition of the proposed problem, the encoding only needs to specify which machine each job will be processed on; thus, the chromosome can be depicted as a vector that has a length of n + m − 1, where n denotes the total number of jobs and m indicates the total number of machines. For instance: , ⋯ , , 0, ⋯ ,0, ⋯ ,0, , ⋯ , , where the sequence of jobs processed on a single machine is represented by a string of numbers that are divided by 0, while the allocation of jobs to various machines is indicated by 0. An encoding example for a solution of a problem with 15 jobs and 4 machines is illustrated in Fig. 2.

Proposed DSMO algorithm
Empirical studies have demonstrated that the SMO algorithm exhibits strong performance when applied to continuous optimization (Bansal et al., 2014;Gupta et al., 2017). However, the problem / , / is a combination optimization problem. The proposed DSMO algorithm incorporates distinct update techniques for discrete individuals during the local leader phase, global leader phase, and local leader decision phase. The following process outlines how these methods can facilitate the inheritance of genes from both the global optimal individual and the local optimal individual by an individual.

Local leader phase
In this stage, each individual updates its position with the information of , which is the local optima in group , and an individual , which is different from in group , as shown in Eq. (14). The individual is updated in two steps. First, the selected parent individual is attracted toward the local leader by Eq. (15). Second, the generated individual , after Step 1, is updated with the randomly selected individual to avoid premature stagnation, as shown in Eq. (16).
where and are uniformly distributed random numbers in the range [0,1]. and are the crossover rates, and , ∈ 0,1 ; = + 0.4/ . We set the initial value to 0.1 in this paper, and is the maximum number of iterations.
, and , represent the crossover operations between individuals, and two kinds of crossover operators are designed. One retains the genes of the two parent individuals and , and the rest of the genes are randomly sequenced, as shown in Fig. 3 (a). The other operation is for two parent individuals without the same genes; a random interval in 1, + − 1 is selected for crossover by mapping, as shown in Fig. 3 (b). The fitness value of the generated individual is calculated, and if the fitness of the new individual is better than that of the old one, then is replaced with the new individual .
= 0.9 + 0.1 = 0.9 + 0.1 where and are uniformly distributed random numbers in the range [0,1]. is the crossover rate, and ∈ 0,1 . , and , represent the crossover operations between two individuals. The crossover operators are the same as in the local leader phase.

Global leader learning phase and local leader learning phase
The updates its position by using a greedy selection process, and , having the best fitness among all the spider monkeys, is selected as the new position of the . If the position of the remains unchanged, then let the global limit count be = + 1. Like the global leader learning phase, the positions of the s in all the groups are updated, selecting the with the best fitness in each group. If the position of the remains unchanged, then let the local limit count be = + 1.

Local leader decision phase and global leader decision phase
If the position of a group is not updated for a predetermined number of iterations, i.e., > , then the positions of the spider monkeys are updated by using information from both the and based on the probability through Eq. (21).
Similar to the local leader decision phase, if the global limit count is larger than the global leader limit, that is, > , then the population is split into subgroups until the number of groups reaches the maximum allowed number of groups , and then they are combined to form a single group again.

HDSMO algorithm
An HDSMO algorithm integrating the merits of DSMO with three improvements is presented in this section. First, a threephase heuristic is proposed to generate a better initial solution. Second, the position update method considering the inertial weight is proposed to balance the exploration and exploitation capabilities of DSMO. Third, a problem-oriented neighbourhood search method with jump and swap operations is designed to improve the search efficiency of the HDSMO algorithm. The flowchart of HDSMO is shown in Fig. 4.

Initial population generation method
Based on the two-phase scheduling heuristic (Lin, Pfund, & Fowler, 2011) and the properties of the addressed problem, a three-phase heuristic is designed to generate the initial solutions.
Step 1: Apply the following linear programming relaxation model to generate a partial schedule. If = 1, job is assigned to machine ; otherwise, job is not assigned to any machine.
Step 3: A group of jobs are assigned to machine in batches considering the machine's reliability threshold to obtain the solution of the proposed problem. The total processing time of batch is estimated by = + ∑ . Then, the lower bound of the batch number is estimated by Eq. (26).
Assuming is the set of unscheduled jobs on machine sequenced by the longest processing time rule, is the first jobs in the set , and is the remaining − jobs. The process of job-grouping batch is shown in Fig.5.   Fig.5. The process of job-grouping batch

Local leader and global leader update with the inertia weight
The balance between global and local search throughout the course of a run is critical to the success of an evolutionary algorithm (Nickabadi et al., 2011). For the basic DSMO algorithm, the old individual position is totally inherited, which is not good for the solution search. The inertia weight is proposed to balance the exploration and exploitation characteristics of DSMO. The method of updating individuals is improved by considering the inertia weight given in Eqs. (27) and (28), and the individual position update with inertia weight is shown in Eq. (29), which represents the inheritance judgement of the individual to the previous position.
where is a random number in the interval [0,1] and when < , operation is performed by mutation. For the individuals in the first and last 50% of the generation population, the mutation operations of order reversal and two-point exchange were adopted, respectively, as shown in Fig.6.  Fig.6. Two kinds of mutation operations: (a) order reversal and (b) swap

Neighbourhood search
For permutation problems, insert or jump and pairwise swap moves are widely used (Chen et al., 2021;Zhang & Chen, 2022) to improve the search efficiency. Two types of neighbourhood operators are defined in this paper to efficiently explore the solution space: ⑴ move one job from one batch to another batch of a machine or from one machine to another machine if it decrease the completion time, which is called a jump; (2) exchange two jobs from different batches of a machine or different machines if it decrease the completion time, which is called a swap.
For a spider monkey individual , intra-machine or inter-machine neighbourhood search is applied according to the deviation rate of the makespan. It is calculated from the maximum and minimum completion times of the parallel machine = . When 0.1 , the intra-machine neighbourhood search operation is applied; otherwise, the intermachine neighbourhood search operation is performed. If the objective function value of the neighbouring solution is better than that of the current solution, the current solution is updated. Otherwise, the current solution is retained.
Define as the number of batches on machine and n as the number of jobs in batch . * and # represent the machines with the longest completion time and the shortest completion time, respectively. The jump and swap operations are performed sequentially for an individual , and a limit ∆ is set for stopping the neighbourhood search between machines. The process of the jump and swap operations is shown in Fig. 7.

Data generation
In this section, we present the computational experiments on the performance of the proposed algorithm. To evaluate the effectiveness of the proposed HDSMO algorithm, HDSMO is compared with the DSMO algorithm and GA. The algorithms were programmed by using MATLAB R2020a software, and the MIP model was implemented on the LINGO 18.0 x64 platform. All programs were run on an Intel(R) Xeon(R) Gold 6242R @3.10 GHz CPU, 64.0 GB RAM workstation.

Fig.7. The main processes of the jump and swap operations
The machine parameters, job parameters, and maintenance activity parameters are utilized in the execution of the experiments. The parameters of the machine consist of the number of machines , whereas the parameters of the job comprise the quantity of jobs , and the time for processing ; and the maintenance activity parameters include the threshold = 0.4 , maintenance time , and machine failure rate, λ. The instances and experimental parameters are shown in Table 1, as were also shown in reference (Avalos-Rosales et al., 2018). We used 4 machines with 20 jobs (n20m4) to represent the problem instances. For each combination of problem instance sizes, 10 instances were randomly generated.

Parametric tuning
Proper adjustment of the parameters is crucial for optimizing the algorithm's performance, hence employing suitable techniques to tune these parameters is essential. A set of 15 typical examples are selected from the problems in Table 1, and the Taguchi method is employed to establish the algorithm's parameters for problems of varying scales (Yue et al., 2019). Taking the HDSMO algorithm as an example, the parameters of the algorithm include the population size N, the maximum number of iterations , crossover rates and and the inertia weight . The parameters and levels of the HDSMO algorithm determined by the pre-experiment are shown in Table 2.  Taking the makespan as the response variable, each parameter combination experiment is run 10 times. The number of orthogonal experiments under 5 parameters and 4 levels is L 4 , and the total number of experiments required is 16*15*10=2400. The orthogonal experimental design was generated in Minitab, and the signal-to-noise (S/N) ratio analysis result of the HDSMO algorithm under the three instance sets is shown in Fig. 8.   Fig. 8. The mean S/N results for the makespan of the HDSMO algorithm For each parameter, the level with the largest S/N ratio is selected. According to the same method, the Table 3 illustrates the determined parameter values for various algorithms, obtained through experimental analysis of the parameters in the GA and DSMO algorithms.

Computational experiments and discussion
The computational results are shown in Tables 4-6 for all the test instances. The objective value of , and the relative percentage deviation (PD) is applied to compare the three evaluated methods, together with the average computation time (CT) in seconds. The PD values were calculated by = * * . * is the optimal solution obtained from the MIP model of the small-scale problems and the best solution obtained by different algorithms for the medium-and large-scale problems. Due to the low efficiency of the MIP models, the medium-and large-scale problems only include the computation times of the heuristics and intelligent algorithms.  It can be concluded from Table 4 that the MIP model can obtain the optimal solutions of all the small-scale test problems, while CT greatly increases with the increase in the number of machines and jobs. HDSMO and GA can obtain the optimal solutions for 3 of 8 instances, and the average PD of HDSMO is smaller than GA. HDSMO outperforms DSMO and GA in terms of and PD, while DSMO and GA outperform HDSMO in terms of CT. Since HDSMO, SMO and GA can almost obtain near-optimal solutions of the test problems, the difference in the and PD values of the three algorithms is not large. Tables 5 & 6 show that HDSMO clearly outperforms DSMO and GA in terms of and PD, and GA clearly outperforms DSMO in terms of CT for the medium-and large-scale problems. It can be inferred that HDSMO is the most promising solution method for solving medium-and large-scale problems, and it can obtain better and PD within a relatively longer CT than GA, and much shorter CT than SMO. The reason lies in the fact that the HDSMO algorithm balances the exploration and exploitation capabilities of search through three improvements compared with DSMO. To verify whether there is a significant difference in the performance of the four evaluated algorithms, an ANOVA test on the mean was applied to compare the solution approaches, and the results are shown in Table 7. Since the significance is less than 0.05, it indicates that there was a statistically significant difference in the performance of the different solution methods. The boxplot diagram at the 95% confidence level for and PD are shown in Fig. 9. MIP is just for the small problems. According to the ANOVA test and the boxplot diagram, we can find that the performance of HDSMO is significantly superior to that of the other two algorithms in the stability and near optimality of the solution, especially for the medium-and largescale problems.
(a) (b) PD Fig. 9. Boxplot diagram of different algorithms As stated in Section 3.3.3, the neighbourhood search (NS) is designed to balance the exploration and exploitation capabilities of the algorithm. The convergence curves of the five algorithms, GA, DSMO, HDSMO and two hybrid approaches with NS, that is, GA+NS and DSMO+NS, are analyzed to further investigate the effectiveness of NS. The results of classical instances for different scale problems are shown in Fig.10.
It can be seen from the Fig.10 that HDSMO is clearly superior to the other solution approaches in terms of both the and convergence speed for all the test problems, and the performance of DSMO and GA is relatively similar with respect to the quality of the solution for medium-and large-scale problems. In addition, DSMO+NS (or GA+NS) outperforms DSMO in terms of both the solution quality and convergence. Generally, it can be inferred that NS can balance the exploration and exploitation capabilities of the DSMO algorithm for the problem proposed in this paper, and HDSMO is superior to the other solution methods with respect to both the solution quality and convergence speed.

Conclusions
To jointly optimize production and maintenance scheduling for the unrelated parallel machines problem, in this paper, a MIP model of the problem is constructed, and Lingo @ is used to solve the small-scale problem. According to the properties of the addressed problem and the "job-allocating and batch-grouping" decision-making method, a hybrid discrete SMO algorithm is proposed to solve medium-and large-scale problems. The experimental results show that the HDSMO algorithm outperforms the GA and DSMO algorithms in terms of the makespan, with slightly more computation time for the mediumand large-scale problems. In addition, the statistical analysis shows that HDSMO is robust to the size of the problems. Meanwhile, a convergence curves comparison of the different algorithms shows that HDSMO is a promising solution method for solving medium-and large-scale problems, and the neighbourhood search method proposed in this paper is effective in improving the convergence of the algorithm.
UPMSPs with maintenance are very common in real-life production environments, where the practical constraints and performance requirements are diverse and complex. Future research will focus on a hybrid DSMO algorithm study for UPMSPs with dynamic event constraints, such as job release times, urgent jobs, sequence-dependent setup times and multiobjective optimization problems. In addition, smart scheduling algorithms integrating the HDSMO and reinforcement learning methods will be investigated to ensure a quick and reasonable response for dynamic events.