Three Hybrid Scatter Search Algorithms for Multi-Objective Job Shop Scheduling Problem

: The Job Shop Scheduling Problem (JSSP) consists of finding the best scheduling for a set of jobs that should be processed in a specific order using a set of machines. This problem belongs to the NP-hard class problems and has enormous industrial applicability. In the manufacturing area, decision-makers consider several criteria to elaborate their production schedules. These cases are studied in multi-objective optimization. However, few works are addressed from this multiobjective perspective. The literature shows that multi-objective evolutionary algorithms can solve these problems efficiently; nevertheless, multi-objective algorithms have slow convergence to the Pareto optimal front. This paper proposes three multi-objective Scatter Search hybrid algorithms that improve the convergence speed evolving on a reduced set of solutions. These algorithms are: Scatter Search/Local Search (SS/LS), Scatter Search/Chaotic Multi-Objective Threshold Accepting (SS/CMOTA), and Scatter Search/Chaotic Multi-Objective Simulated Annealing (SS/CMOSA). The proposed algorithms are compared with the state-of-the-art algorithms IMOEA/D, CMOSA, and CMOTA, using the MID, Spacing, HV, Spread, and IGD metrics; according to the experimental results, the proposed algorithms achieved the best performance. Notably, they obtained a 47% reduction in the convergence time to reach the optimal Pareto front.


Introduction
The Job Shop Scheduling Problem (JSSP) consists of a set of jobs, formed by operations, which must be processed in a set of machines subject to restrictions of precedence and resource capacity. For a job to be completed, all of its operations must be processed in a given sequence. This problem belongs to the NP-hard class [1], is challenging of solving it, and has significant industrial applicability [2]. In JSSP, we must determine the order or sequence for processing a set of jobs through several machines minimizing one or more objective functions. An essential function of JSSP is the coordination and control of complex activities, both optimum resource allocation and the sequence in performing those activities [3].
In a real operation context, it is common to consider more than one criterion simultaneously, which defines a multi-objective optimization problem whose solution involves generating a set of non-dominated solutions. This set provides the decisionmaker with several alternatives to choose the one according to the needs of the manufacturing process [4][5][6].
A detailed analysis of the state-of-the-art for multi-objective JSSP shows that few works approach the problem from the multi-objective perspective, in which at least three

Related Literature
The Pareto Archived Simulated Annealing (PASA) algorithm was applied to a JSSP with two objectives: the makespan and the mean flow time [12]. This algorithm was evaluated with 82 instances from the literature. The results obtained are evaluated in terms of the number of non-dominated schedules generated by the algorithm and the proximity of the obtained non-dominated front to the Pareto front.
A successful algorithm based on Simulated Annealing (SA) for several objectives named AMOSA was proposed [13]; also, it was reported that this algorithm performed better than some MOEA algorithms, such as the NSGA-II [14].
A two-stage genetic algorithm (2S-GA) was proposed for JSSP. This algorithm is applied to minimize three objectives makespan, total weighted earliness, and total weighted tardiness [15]. This algorithm is composed of two Stages: Stage 1 applies parallel GA to find the best solution for each individual objective function, and in Stage 2 the populations are combined. The performance of the algorithm is tested with benchmark instances and compared with results from other published papers. The experimental results show that 2S-GA is efficient in solving instances of the job shop scheduling problem in terms of the quality solution.
The Contemporary Design and Integrated Manufacturing Technology (CDMIT) laboratory proposed the Improved Multiobjective Evolutionary Algorithm based on Decomposition (IMOEA/D) to minimize the makespan, tardiness, and total flow time [16]. This algorithm was evaluated with 58 instances using the performance metrics Coverage [17] and Mean Ideal Media (MID) [18] to the evaluation. This algorithm stands out from the rest of the literature because it considered three objectives, applied two performance metrics, and reported good results.
In 2017, was proposed a hybrid algorithm with NSGA-II and a linear programming approach [19] to solve the FT10 instance [20]. In this paper, the objectives are to minimize weighted tardiness and energy costs. Furthermore, the authors applied the Hypervolume metric to evaluate the performance. Furthermore, in 2019, was proposed MOMARLA; this is a new algorithm based on Q-Learning to solve MOJSSP [21]. In this work, each agent represents a specific objective and uses two action selection strategies to find a diverse and accurate Pareto front.
The Contemporary Design and Integrated Manufacturing Technology (CDMIT) laboratory proposed the Improved Multiobjective Evolutionary Algorithm based on Decomposition (IMOEA/D) to minimize the makespan, tardiness, and total flow time [16]. This algorithm was evaluated with 58 instances using the performance metrics Coverage [17] and Mean Ideal Media (MID) [18] to the evaluation. This algorithm stands out from the rest of the literature because it considered three objectives, applied two performance metrics, and reported good results. Furthermore, in 2021, two multi-objective algorithms for minimizing makespan, total tardiness, and total flow time, were published. These algorithms are Chaotic Multi-Objective Simulated Annealing (CMOSA) and Chaotic Multi-Objective Threshold Accepting (CMOTA) [11]. They incorporate dominance criteria and a chaotic perturbation strategy to improve their performance. Experimental evaluation results indicated that they overpassed the state-of-the-art algorithms MOMARLA, MOPSO, CMOEA, and SPEA [21]. The algorithms proposed in the present paper seek to enhance the best algorithms of this group.
Finally, the scheduling will probably be directed to increasingly automated and use intelligent systems in the future. Under the Industry 4.0 environment, the computational workload could be greatly reduced and the systems probably will become more flexible and agile [22].

Background
This section describes basic concepts and algorithms in the multi-objective area which are related to this work. Furthermore, we present the multi-objective Job Shop Scheduling formulation and the main performance metrics used in this work.

Multiobjective Optimization Concepts
The Multi-Objective optimization algorithms use the concept of domination where two solutions are compared to determine if one solution dominates the other or not. Key concepts for Multi-Objective optimization are described below.
Pareto Dominance: For any optimization problem, solution A dominates another solution B if the following conditions are met: A is strictly better than B on at least one objective, and A is not worse than B in all objectives [23].
Non-dominated set: Among a set of P solutions, the subset of non-dominated solutions P1 is integrated by solutions that accomplish the following conditions: • Any pair of P1 solutions must be non-dominated (one regarding the other) • Any solution that does not belong to P1 is dominated by at least one member of P1 [23].
Pareto front: It is the graphical representation of the non-dominated solutions in the space of the objectives of the multi-objective optimization problem [24].

Performance Metrics
In the case of Multi-Objective Optimization, defining quality is complicated because two or more conflicting objective functions could exist. Then in an experimental comparison of different optimization algorithms, it is necessary to have the notion of performance. Some of the performance metrics are shown in Table 1.
A large number of performance metrics or quality indicators can be found. These metrics consider mainly the following three aspects of a non-dominated solution set [25]: • Convergence: that is a feature related to the closeness to the theoretical Pareto optimal front.
• Diversity: this feature for any distribution of non-dominated solutions is measured by Spread and Spacing.

•
The number of solutions.
It is difficult to find a single performance metric that encompasses all of the above criteria. However, according to the characteristics they measure, the metrics can be grouped as [25]:

Metric Type Formula
Mean Ideal Distance Accuracy Hypervolume Inverted Generational Distance Accuracy/Diversity The MID metric is calculated with Equation (1). This metric calculates the closeness of the calculated Pareto front (PFcalc) solutions with an ideal point [18]. In this equations, Q is the number of solutions in the PFcalc, = √( 1 ) 2 + ( 2 ) 2 + ( 3 ) 2 , and 1 , 2 , 3 are the values of the i-th non-dominated solution for their first, second, and third objective function, respectively. In Equation (2) is showed the formula to calculate S; this metric evaluates the distribution of the non-dominated solutions in the PFcalc. The algorithm with the smallest S value is the best [26]; di measures the distance in the space of the objectives functions between the i-th solution and its nearest neighbor; that is the j-th solution in the PFcalc of the algorithm, ̅ is the average of di, that is ̅ = ∑ ⁄ =1 and = (| 1 ( ) − 1 ( )| + | 2 ( ) − 2 ( )| + ⋯ + | ( ) − ( )|) , where 1 , 2 are the values of the j-th non-dominated solution for their first and second objective function, respectively. Furthermore, M is the number of objective functions and i,j=1,…,Q. Equation (3) shows the formula to calculate HV. In this formula represents a hypercube which is constructed with a reference point W (this can be found constructing a vector with the worst values of the objective function) and the solution i as the diagonal corners of the hypercube [23]. This metric calculates the volume in the objective space that is covered by all solutions in the non-dominated set [27]. Therefore, an algorithm that obtains the highest HV value is the best. The data should be normalized to calculate the HV by transforming the value in the range [0, 1] for each objective separately.
The Spread metric is calculated using Equation (4), which considers the distance to the extreme points of the True Pareto front (PFtrue) and was proposed to have a more precise coverage value [14]. Where measures the distance between the extreme point of the PFtrue for the k-th objective function and the nearest point of the PFcalc.
Finally, Equation (5) shows the formula to calculate IGD; this metric finds the average distance between the points of the PFtrue to the PFcalc [28]. Where = { 1 , 2 , … | | } that is, the solutions in the PFtrue and |T| is the cardinality of T, p is an integer parameter, in this case, p = 2 and ̂ is the Euclidian distance from tj to its nearest objective vector q in Q. In where fm(t) is the m-th objective function value of the t-th member of T. Another important metric is the number of non-dominated solutions generated by the algorithm. The greater the number of solutions, the greater the number of alternatives the decision maker will have to choose the desired solution [4][5][6].

MOJSSP Formulation
In JSSP, there are a set of n jobs, consisting of operations, which must be processed in m different machines. There are a set of precedence constraints for these operations, and there is a resource capacity constraint for ensuring that each machine should process only one operation at the same time. The processing time of each operation is known in advance.
The objective of JSSP is to determine the sequence of the operations in each machine (the start and finish time of each operation) to minimize certain objective functions. The most common objective is the makespan, which is the total time in which all the problem operations are processed. Nevertheless, real scheduling problems are multi-objective, and several objectives should be considered simultaneously.
This work tries to optimize three objectives simultaneously, makespan, total tardiness, and total flow time. The formal MO-JSSP model can be formulated as follows [29,30]: where q is the number of objectives, x is the vector of decision variables, and S represents the feasible region, defined by the next precedence and capacity constraints, respectively: The objective functions of makespan, total tardiness, and total flow time, are defined by Equations (7)-(9), respectively.
where Cj is the makespan of job j.
where Tj = max (0, Cj − Dj) is the tardiness of job j, and Dj is the due date of job j and is calculated with = ∑ , =1 [31], where pj,i is the time required to process the job j in the machine i. In this case, the due date of the j job is the sum of the processing time of all its operations on all machines, multiplied by a narrowing factor ( ), which is in the range 1.5 ≤ t ≤ 2.0 [31,32].

Proposed Hybrid Algorithms MO
Three hybrid algorithms based on Scatter Search (SS) are proposed for the solution of the MO-JSSP problem. Figure 1 shows the Scatter Search framework showing the three different process which distinguish the three proposed hybrid algorithms. Algorithm 1 contains the pseudocode of our Scatter Search Algorithm, which is described in detail in the next subsection. Notice that in line five, one of the algorithms CMOSA, CMOTA , or LS can be executed. The goal of Algorithm 1 is to improve the solutions of the reference set.

Scatter Search General Framework
Scatter Search (SS) is an algorithm proposed by Glover [10], and it is composed of the following methods: • Generator of diverse solutions, in which a set P of diverse solutions of size 30 is generated.

•
Updater and creator of refSet, from the P solutions, the first three non-dominated and three most diverse are selected, using the Euclidean distance, to form the reference set (RefSet) of size 6.

A Chaotic Threshold Accepting Algorithm (CMOTA). Threshold Accepting
(TA) is an algorithm proposed in [33]. In this enhanced method CMOTA, a version adapted to JSSP is used [11]. 3. Chaotic Simulated Annealing Algorithm (CMOSA), SA was originally proposed in [34], and in 2021 a new version is implemented under the multiobjective approach [11].
In improvement processes 2) and 3), an analytical tuning process is performed for the algorithm parameters [35]. A regular perturbation is also applied to generate a new solution that is compared to the current one. From the previous comparison, the dominant one is selected, and the non-dominated is discarded. If both are not dominated, one is saved in the set of non-dominated solutions, and the other continues the search. When new non-dominated solutions are not found in both algorithms, a chaotic perturbation is applied. This perturbation consists of using the equation of the logistic maps [36] as a mechanism to escape from stagnation and search diversity in the solutions. A reheating process is also applied, which consists of raising the value of the current temperature parameter to be able to carry out a new scan. The improvement algorithms implemented are described in more detail in the following sections.

Hybrid Scatter Search with Local Search (SS/LS)
In this algorithm, a Local Search (LS) is applied in the solution improvement phase. Algorithm 2 shows the local search algorithm used. In this algorithm, a set of iterations is performed. In each of them, a regular perturbation is applied to the solution (exchange of two operations) to generate a new one. Dominance is verified between the current solution and the new one created with the perturbation in each iteration. In this verification, three possible cases are generated: • Case A. If the new solution dominates the current one, then the new solution is saved and replaces the current one to continue the search process. • Case B. If the current solution dominates the new one, then the new solution is discarded.
• Case C. If the current and the new solution are non-dominated, the current solution is saved, and the new one replaces the current one to continue the search.

Hybrid Scatter Search with Chaotic Simulated Annealing (SS/CMOSA)
The SS/CMOSA algorithm is based on Chaotic Simulated Annealing Algorithm (CMOSA) [11]. This hybrid algorithm, SS/CMOSA, receives the solutions obtained by a combination process as shown in Figure 1, while CMOSA is shown in Algorithm 3.: one that controls the stop condition and the other internal (Metropolis cycle) that controls the number of iterations carried out for each temperature parameter value. In Algorithm 2, a perturbation is made to the current solution in the internal cycle to generate a new solution. Dominance is verified between the current solution and the new one in each iteration. In this verification, the same three possible previous cases are generated: • Cases A and C are evaluated similarly to the SS/LS version. • Case B. According to the Boltzmann probability distribution, when the current solution dominates the new one, the latter could be replaced by the former.
Additionally, a chaotic perturbation and a reheating are applied when stagnation occurs, consisting of a predetermined number of iterations without finding nondominated solutions. The chaotic perturbation uses the equation of the logistic maps [37], whose main characteristic is that it generates different outputs even in small changes in its input data. Then chaos or chaotic perturbation is a process carried out to restart the search from another point in the space of solutions. Reheating is the process by which the current temperature parameter of the SA algorithm is high; this helps to perform reprocessing that allows reinitializing the search process.

Hybrid Scatter Search with Chaotic Threshold Accepting (SS/CMOTA)
The Chaotic Threshold Accepting (CMOTA) algorithm is used as an improvement method. In CMOTA (Algorithm 4), there are also two cycles such as CMOSA, one that controls the stop condition (temperature) and the other internal (Metropolis cycle) that controls the number of iterations that are carried out for each value of the temperature parameter.
The same three possible previous cases are generated. Cases A and C are evaluated in the same way as the SS/CMOSA version. In Case B, if the current solution dominates the new solution, the latter can replace the current one in the searching process by using a threshold established in the algorithm.
Similar to CMOSA, this algorithm also applies chaotic perturbation and a reheating process.

Computational Experimentation
This section describes the dataset used, the conditions of the experimentation as well as the results obtained.

Datasets
To perform the experimental evaluation of the algorithm, a benchmark of 70 instances of the problem are used. All the solutions of this dataset have different sizes and degrees of complexity. The instances are divided in six sets: •  TA01, TA11, TA21, TA31, TA41, TA51, TA61, and TA71 taken fromTaillard [41].
The instance sizes of this dataset range from six jobs on six machines (the FT06 instance) to 100 jobs on 20 machines (the instance TA71).

Experiment Description
Two experiments were carried out with a different number of instances to evaluate the performance of the proposed algorithms.
The first experimentation was carried out with 58 common instances with only three algorithms (IMOEA/D [16], CMOSA, and CMOTA [11]) of the state-of-the-art that used the MID performance metric and the same three objectives. The second experimentation was carried out with 70 instances to compare the results with the CMOSA and CMOTA algorithms proposed in [11]. Each instance was executed 30 times using 30 initial solutions. The set of non-dominated solutions was obtained from the total solutions generated by the 30 executions The performance metrics used are MID, Spacing, Spread, HV, IGD, Runtime, and the number of non-dominated solutions. Finally, we applied these metrics to the set of nondominated solutions obtained by the algorithms at the end of their processes.
The execution of the proposed algorithms was carried out in a terminal of the Ehécatl cluster of the Technological Institute of Ciudad Madero, with the following characteristics: Intel ® Xeon ® processor at 2.30 GHz, Memory: 64 Gb (4 × 16 Gb) ddr4-2133, Linux operating system CentOS. C language was used for the implementation, and GCC compiler. Table 2 shows a comparison with the average results obtained by the MID metric for the 58 instances used in the algorithms IMOEA/D [16], CMOSA [11], CMOTA [11], and the proposed SS algorithms. We observed that our hybrid algorithms SS/LS, SS/CMOTA, and SS/CMOSA obtained the best results. Furthermore, the hybrid SS/CMOSA obtained the best result surpassing IMOEA/D by 17%.  Table 3 shows the results obtained for each of the 70 instances by the three proposed algorithms, comparing them with the best of the state-of-the-art (CMOSA, CMOTA), taking as a reference the value obtained by the MID metric. The last row shows the average value of the MID metric for each of the algorithms. In this table, we observed that the SS/LS algorithm obtained a better performance since the value of the MID metric is smaller than the other algorithms analyzed.  Table 4 summarizes the experiment results with 70 instances and five metrics; it contains the average values obtained by executing the algorithms 30 times. The first column has the name of the evaluated metric. The next two columns show the results of the two best state-of-the-art algorithms (CMOSA and CMOTA). Finally, the last three columns show the results for the three hybrid proposed algorithms (SS/LS, SS/CMOSA, and SS/CMOTA). In Table 4, the best values are highlighted and marked with an asterisk (*). An approximate Pareto front is used in the case of metrics in which it is necessary to use a True Pareto front [24]. The approximate Pareto front is generated by previous executions of developed algorithms throughout the study of this problem.

Comparative Results
We can observe that SS/LS obtains the best result for MID and IGD metrics. This algorithm uses the shortest processing time, which means that it has the best convergence and generates the highest number of non-dominated solutions of the three proposed hybrid algorithms. The results indicate that the solutions found by SS/LS are closer to the origin point (0,0,0); they are closer on average to the approximate front and were achieved with the lowest amount of execution time. SS/CMOSA algorithm obtains the best Spread, which means that the generated solutions are very well distributed on the non-dominated solutions front. On the other hand, SS/CMOTA achieved the best Spacing and HV values, which indicate that this algorithm-has a more uniform Spacing and the best solution space coverage. In other words, the hybrid algorithms proposed in this paper obtained the best results for these datasets.
Finally, the non-dominated solution fronts obtained by the proposed algorithms for the 70 instances are included in Appendix A.

Conclusions
This paper presents three Hybrid Multi-Objective algorithms for JSSP, named SS/LS, SS/CMOSA, and SS/CMOTA. Three objectives are considered: makespan, total tardiness, and total flow time. Furthermore, we present an experimental evaluation applying six performance metrics.
Regarding the results from the comparison, we observe that SS/LS generates solutions closer to the origin point and PFapprox. It provides more solutions than the others' algorithms and uses the minimum runtime. SS/CMOSA generates solutions with a better distribution concerning PFapprox. Furthermore, SS/CMOTA generates better-distributed solutions in the PFcalc and better coverage, as shown by the HV metric. The results obtained by the proposed algorithms SS/LS, SS/CMOSA, and SS/CMOTA compared with some of the best algorithms in the literature show that they are among the best in the area. We highlight that our proposed SS/LS algorithm reduces processing time by 43% compared to the fastest in the state-of-the-art.

Informed Consent Statement: Not applicable.
Data Availability Statement: In this section, Data supporting are located in this paper.

Acknowledgments:
We thank LanTI Lab for letting us use its computers and Conacyt for the scholarship to some of us while we studied our Ph.D.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Non-Dominated Solutions Obtained
The non-dominated solutions obtained by the proposed SS/LS, SS/CMOSA, and SS/CMOTA algorithms for the 70 instances used in this paper are shown in Tables A1-A6,  Tables A7-A12 and Tables A13-A18, respectively.
In these tables, MKS is the makespan, TDS is the total tardiness, and FLT is the total flow time. For each instance, the best value for each objective function is highlighted with an asterisk (*) and in bold type.  Table A2. Non-dominated solutions obtained by SS/LS for the JSSP instances proposed by [37].  Table A3. Non-dominated solutions obtained by SS/LS for the JSSP instances proposed by [38]. LA02  LA03  LA04  LA05  MKS  TDS  FLT  MKS  TDS  FLT  MKS  TDS  FLT  MKS  TDS  FLT  MKS  TDS Table A4. Non-dominated solutions obtained by SS/LS for the JSSP instances proposed by [39]. ABZ6  ABZ7  ABZ8  ABZ9  MKS TDS  FLT MKS TDS FLT MKS TDS  FLT  MKS TDS  FLT MKS TDS Table A6. Non-dominated solutions obtained by SS/LS for the JSSP instances proposed by [41]. TA01  TA11  TA21  TA31  TA41  MKS TDS  FLT MKS TDS  FLT MKS TDS  FLT MKS TDS  FLT MKS TDS Table A8. Non-dominated solutions obtained by SS/CMOSA for the JSSP instances proposed by [37]. ORB2  ORB3  ORB4  ORB5  MKS TDS  FLT MKS TDS FLT MKS TDS  FLT MKS TDS  FLT MKS TDS FLT Table A9. Non-dominated solutions obtained by SS/CMOSA for the JSSP instances proposed by [38].