Multimodal Optimization of Permutation Flow-Shop Scheduling Problems Using a Clustering-Genetic-Algorithm-Based Approach

: Though many techniques were proposed for the optimization of Permutation Flow-Shop Scheduling Problem (PFSSP), current techniques only provide a single optimal schedule. Therefore, a new algorithm is proposed, by combining the k-means clustering algorithm and Genetic Algorithm (GA), for the multimodal optimization of PFSSP. In the proposed algorithm, the k-means clustering algorithm is ﬁrst utilized to cluster the individuals of every generation into different clusters, based on some machine-sequence-related features. Next, the operators of GA are applied to the individuals belonging to the same cluster to ﬁnd multiple global optima. Unlike standard GA, where all individuals belong to the same cluster, in the proposed approach, these are split into multiple clusters and the crossover operator is restricted to the individuals belonging to the same cluster. Doing so, enabled the proposed algorithm to potentially ﬁnd multiple global optima in each cluster. The performance of the proposed algorithm was evaluated by its application to the multimodal optimization of benchmark PFSSP. The results obtained were also compared to the results obtained when other niching techniques such as clearing method, sharing ﬁtness, and a hybrid of the proposed approach and sharing ﬁtness were used. The results of the case studies showed that the proposed algorithm was able to consistently converge to better optimal solutions than the other three algorithms.


Introduction
Scheduling problems such distribution, transportation, management, construction, engineering and manufacturing, are omnipresent in real-world applications [1]. These scheduling problems are complex and very difficult to solve using conventional techniques. The complexity of these scheduling problems arises due to the presence of a large number of constraints. Due to these reasons, there was ample research conducted in the past few decades that focused on developing different techniques for the optimization of these scheduling problems. Permutation flow-shop scheduling problem (PFSSP), a variation of the classical flow-shop scheduling problem (FSSP), is a problem in operation research in which n ideal jobs, J1, J2, . . . , Jn, of varying processing times are assigned to m resources (usually machines). The objective is to allocate the available resources to the jobs to maximize (or minimize) performance indicators, such as makespan, tardiness, etc. In classical FSSP, for n jobs on m machines, there are (n!)m different alternatives for sequencing jobs on machines, while in permutation problems, the search space is reduced to n!. Like any scheduling problem, there are constraints and assumptions associated with PFSSP. In this paper, these constraints and assumptions are as follows.
(1) On the shop floor, all machines are independent from each other. (2) The operation of all jobs have deterministic processing times. (3) No pre-emption of job operation is permitted.
(4) The set-up times are not considered in these problems. (5) No machine down-time or maintenance time is considered. (6) The machine used for each operation of a job is known and fixed. (7) A job cannot revisit a machine group for more than one operation. (8) Operation constraints must always be followed. (9) All jobs must enter the machines in the same order.
Since PFSSP are commonly encountered in various industries, researchers developed many methods to optimize the key performance indicators (KPIs) of PFSSP. Liu and Liu [2] used a hybrid discrete artificial bee colony algorithm for minimizing the makespan in PFSSP. They utilized Greedy Randomized Adaptive Search Procedure (GRASP) to create the initial population and operators such as insert, swap, path relinking to generate new solutions. The authors also improved the local search by further optimizing the best solution. When applied to PFFP, their algorithm had superior performance compared to other algorithms like particle sawm optimization (PSO) algorithm, hybrid GA (HGA), etc. Govindan et al. [3] combined decision tree (DT) and scatter search (SS) algorithms to solve PFSSP. DT was used initially to convert the problem into a tree structure using the entropy function, followed by SS, to do an extensive investigation of the solution space. Simulation results and statistical test comparisons showed the advantage of the authors' proposed algorithm. Ancău [4] proposed two algorithms for solving PFSSP. (1) A constructive greedy heuristic (CG) and (2) stochastic greedy heuristic (SG). The CG was based on a greedy selection, while the SG was a modified version of CG with iterative stochastic start. The authors validated the proposed algorithms by using them to solve benchmark problems, and the results showed that the algorithm was able to come within 6% of the best-known solution. Zobolas et al. [5] created a hybrid algorithm by combining greedy heuristic, GA, and variable neighborhood search (VNS) algorithm. The hybrid algorithm was able to take advantage of both GA and VNS, and obtain the best-known makespan for the benchmark problems in a short computational time.
Zhang et al. [6] proposed an enhanced GA to solve the distributed assembly permutation flow shop scheduling problems (DAPFSP). The authors used a greedy mating pool to select better candidates. They also updated the crossover strategy and incorporated local search strategies to improve local convergence. Andrade et al. [7] minimized the total flowtime of PFFSP by using a Biased Random-Key Genetic Algorithm (BRKGA). Their proposed technique used a shaking strategy to perturb individuals from the elite population, while resetting the remainder of the population. The authors validated their proposed approach by its application to 120 benchmark problems and comparing the results to those obtained by Iterative Greedy Search and Iterated Local Search. The results showed that the proposed approach was more efficient and was able to identify new bounds as well as optimal solutions for the problems. Ruiz et al. [8] used the iterated Greedy (IG) search algorithm to solve the distributed permutation flow shop scheduling problems. The authors improved the initialization, construction, and deconstruction procedures, and improved the local search (ILS) to obtain better results. A big advantage of the proposed technique was that it required very little problem-specific knowledge to be implemented. Pan et al. [9] proposed three constructive heuristics, DLR (inspired by LR heuristic), DNEH (inspired by NEH heuristic), a hybrid DLR-DNEH heuristic, and four metaheuristics-Discrete Artificial Bee Colony (DABC), Scatter Search (SS), ILS, & IG-for minimizing the total flowtime criterion of the distributed permutation flow shop scheduling problem (DPFSP). Both the heuristics and metaheuristic were able to significantly improve upon the existing results in literature. Their results also showed that the trajectory-based metaheuristics were considerably better than the population-based metaheuristics.
Though ample work is done in solving PFSSP, the current techniques only provide a single optimal schedule upon execution. Furthermore, real-life operators and schedulers have preferences that cannot be modeled and added to the fitness function of optimization algorithms. The single optimal schedule found using these techniques might theoretically satisfy the objective function. However, it might not be able to satisfy the aforementioned complexities. Although multimodal optimization provides multiple similar solutions for the same problem under consideration, it cannot handle complex situations such as machine breakdowns, delayed supply chain, etc., without changing the fitness function and constraints. However, having multiple similar solutions for the same problem could make it possible to provide opportunity for the operators to have the flexibility to handle emergency using their experience, while waiting to solve the underlying problem. Therefore, it is important to obtain multiple optimal schedules that satisfy the objective function, i.e., perform multimodal optimization (MMO). Over the years, many different techniques were developed to accomplish MMO, such as niching techniques [10][11][12][13][14][15][16], differential evolutionary strategies [17][18][19][20][21][22][23][24], use of GA [25][26][27][28] etc. However, these techniques were only tested on benchmark mathematical problems. Furthermore, the aim of these techniques was to find both global and local optima, which increases their computational complexity. Additionally, it is possible for local optima to have much worse objective value than global optima, which is undesirable in scheduling problems.
Therefore, in this study, a new algorithm is proposed for the MMO of PFSSP. This is accomplished by using the k-means clustering algorithm with GA. In each iteration of the proposed algorithm, first, solutions are clustered together, using k-means clustering algorithm. Next, the operators of GA are utilized within each cluster to generate solutions for the next generation. Various techniques used by Perez et al. were also modified and utilized for the MMO of PFSSP and the performance of the algorithms was tested on benchmark PFSSP. The rest of the paper is organized as follows. Section 2 presents information about the modified GA operators used in the proposed algorithm to represent a solution and create new solutions during the iterations of the algorithm. Section 3 provides detailed information about the proposed algorithm, the clearing method (CM), and the sharing fitness methods used for the MMO of PFSSP, as well as the features used to cluster the solutions. Section 4 presents information about the benchmark problems used to test the proposed algorithm and presents the results obtained by utilizing various algorithms. Section 5 draws upon conclusions based on the results obtained from the case studies and mentions the future direction of the proposed technique.

Encoding Method and Initial Solution Creation
In this study, the optimization objective was to find the sequence of operations on each machine that would lead to a minimization of the makespan. An important task in any optimization problem is to determine the encoding method. Over the years, many different methods were used to encode the solutions of scheduling problems, such as the direct method [29], the binary method [30], the circular method [31], and the permutation with repetition [32] method. As the jobs must be processed in the same order on each machine in PFSSP, a solution is created using permutation encoding for a single machine. The length of the chromosome (solution) can be determined by the size of the problem. For example, for a 9 × 3 problem (9 jobs with 3 operations each) the chromosome length would be 9. An example of permutation for the 9 × 3 problem is shown in Figure 1. According to Figure 1, the order of operations on any machine will be job 3, 1, 2, 9, 5, 4, 6, 8 and 7. By utilizing permutation encoding, only feasible initial solutions would be created, which would reduce the computational cost of the algorithm.

Adaption of Genetic Operators
In the optimization process, GA operators are utilized to produce the initial population, as well as the population belonging to the next generations. To create only feasible solutions for the iterations of the algorithm, appropriate crossover and mutation operator must be developed. In the crossover operator, half the number of jobs of the problem under consideration were randomly selected and the corresponding genes were switched on in both the parents. An example of the crossover process is shown in Figure 2 for a PFSSP with nine jobs. In the example, jobs 2, 3, 4, 5, 6, and 9 were randomly chosen and the corresponding genes (3,2,4,9,5,6) in Parent A and (2,3,6,9,5,4) in Parent B were switched on. The crossover operator was limited to solutions belonging to the same cluster, in order to preserve some features of the cluster. Similarly, in the mutation operator, half the jobs were randomly chosen and their corresponding location within the solution was switched on. An example of the mutation operator is shown in Figure 3 for a PFSSP with nine jobs. As can be seen from Figures 2 and 3, using the crossover and mutation operator only creates feasible solution.

Multimodal Optimization Methodologies
As stated earlier, an algorithm is proposed for the MMO of PFSSP. To accomplish this, the GA is combined with the k-means clustering algorithm, in order to converge to multiple optimal solutions. Additionally, GA is also combined with niching techniques, such as sharing fitness and CM, and the performance of the three algorithms is compared. Sharing fitness and CM were used as these were shown to have good performance when utilized for the MMO of JSSP [33,34]. The steps of the three different algorithms are outlined below.

Clustering-Optimization Approach
In traditional GA, the elite population survives each generation as they represent the best solution. Majority of the offspring are generated using crossover which has a high probability of selecting individuals with better fitness values. Though these operations ensure that the population would converge towards an optimal solution, they also reduce the diversity of the population significantly. This strategy is applicable when the objective is to obtain a single optimal solution. However, it is severely lacking when multiple optimal or near optimal solutions are required.
To achieve MMO, the diversity of the elite population and the remainder of the population must be maintained throughout the iterations of GA. In the proposed approach, this is achieved by splitting the population of GA in multiple niches, based on similarity and dissimilarities. Since cluster algorithms maximize the distance between different clusters and by assuming that different optimal solutions have high dissimilarity, we can utilize the cluster algorithms to split the search space into different clusters and find multiple optimal solution in each cluster. By limiting crossover and optimization within each niche, we can maintain the high diversity of the population, as we can avoid exposure from other clusters. Finally, to search the solution space more effectively, mutation is performed after every iteration. The basic procedure of the proposed algorithm is given below and is also shown in Figure 4.

1.
Define the fitness function of the specific case study. In this study, the maximum completion time (makespan) was used as the fitness value.

2.
Define the parameters for the k-means clustering algorithm (k value) as well as population size, crossover rate, mutation rate, generation limit, and initialize the best solution set. 3.
Generate the initial population. 4.
Evaluate the current population, calculate the feature matrix of every solution, and update the best solution set. 5.
Using the k-means clustering algorithm, divide the current population into k clusters, based on the features calculated in step 4. 6.
Within each cluster, select elite individuals that would be a part of the next generations population. 7.
Within each cluster, select individuals for crossover using roulette wheel selection to produce offspring for the next generation. 8.
Randomly select individuals (based on mutation probability) and mutate them. 9.
Repeat Steps 4-8 until the generation limit is reached. Output the final population and the best solution set.
To cluster solutions together, a feature matrix must be defined. Typically, in MMO, the distance between two individuals is used as a feature to cluster the solutions together. However, in PFSSP, a distance cannot be defined using the proposed encoding method, therefore, the feature matrix used to cluster the solutions together in the proposed algorithm is defined as the sequence of jobs on a particular machine. For example, the feature matrix for Parent A shown in Figure 2 would be (1, 3,2,4,9,5,8,6,7). An example of applying the proposed algorithm is shown in Figure 5. It should be noted that clustering would not improve the algorithm efficiency. However, that is not the top priority in this study. The focus is to obtain multiple and different solutions with the same objective value.

Sharing Fitness
In the sharing fitness method, the fitness of an individual decreases in accordance with the number of similar individuals in the population. Figure 6 shows the overall flow chart of the sharing fitness method. After evaluating the individuals of the current population, their shared fitness is calculated according to Equation (1): where, f i is the fitness of individual i, and SH(d(i,j)) is the sharing function. Here, the sharing function depends on the distance (d(i,j)) between individuals i and j. The sharing function can be calculated using Equation (2).
where σ share is the maximum distance of a niche and α is the slope of the sharing fitness function. A value of 1 is most commonly used for the slope, while σ share depends on the problem under consideration. Next, based on the new fitness f * i , parents are selected for crossover and mutation, using the roulette-wheel selection method. Once offspring are generated using modified crossover and mutation operators, their fitness is evaluated. Next, the fitness of the parents and offspring are updated to f * * i using Equation (3). where The distance in Equations (1)-(4) is problem dependent and for the sharing fitness and CM, it was defined as the number of operations situated in different places for each machine. For example, the distance between Parents A and B shown in Figure 2 would be 8 (only job 1 is located in the same place for the two solutions). The above definition of distance was used for CM and the sharing fitness method, as it was shown to provide a large number of optimal solutions when used by Perez et al. for the MMO of JSSP.

Clearing Method
CM is similar to sharing fitness but is based on the concept of limited resources of the environment. CM is applied after evaluating all individuals but before the selection process. In CM, the individuals are ordered according to their fitness, from best to worst. The first individual is called the dominant (and belongs to the first niche) and it is compared with the remaining n−1 individuals. The comparison is done to determine whether the remaining individuals belong to the first niche. Only k individuals from each move forward to the next generation as elite individuals, while the rest of the individuals belonging to that niche have their fitness reset to zero. The above procedure is repeated for individuals whose fitness is greater than zero. A pseudocode of the CM is given below, and a flowchart is given in Figure 7.

Experimental Study and Performance Benchmarking
To measure the performance of the algorithms outlined above, these were used for the MMO of multiple PFSSP. Eighteen benchmark problems from various datasets were taken for experiments. The first nine problems were proposed by Carlier [35], the next eight were proposed by Reeves [36] and the last one was proposed by Heller [37]. To better assess the algorithms, 10 independent simulations were performed for each test problems. The performance of the algorithms used was then measured on the basis of four indicators. (1) The best solution found in the 10 simulations. (2) The mean value and the standard deviation of the optimal solution found for the 10 simulations. (3) The number of times the algorithms converged to the best optimal solution. (4) The average number of best optimal solutions found for the 10 simulations. For indicator number 4, only simulations in which the algorithm converged to the best solution from (1) were taken into consideration i.e., if the algorithm only converged to the best solution in 4 out of 10 simulations, then the average number of optimal solutions found was calculated using those 4 simulations. The parameters used for the three algorithms were determined after numerous simulations and are given in Table 1. Additionally, the performance of the three algorithms was also compared to the performance of Improved Efficient Genetic Algorithm (IEGA) and the Hybrid IEGA (HIEGA) proposed by Basset et al. [38], when applicable. These algorithms were selected as they showed promising results when used for the optimization of benchmark PFSSPs.

The Best Solution Found
As mentioned above, 10 independent simulations were ran for each test problem using the three algorithms. The first method used to compare the performance of these algorithms included observing the best solution obtained during 10 simulations. The results obtained by these algorithms were also compared to the best-known solution, as well as the best-known solutions obtained using different algorithms in the available literature. These results are shown in Table 2. According to Table 2, the proposed algorithm, sharing fitness algorithm, and HIEGA were able to find the best-known solutions for all Carlier benchmark problems. IEGA was able to converge to the best-known solution for all Carlier problems besides Car-03, while CM was able to find the best-known solution for Car-01, 05, 06, 07, and 08. For the Reeves problems, the HIEGA had the best performance, as it was able to converge to the best-known solution for 6 out of the 9 problems. The proposed algorithm was able to obtain the best-known solution for Rec-07, Rec-09, and Rec-11, and was able to obtain a better solution than HEIGA for . It was also able to converge within 1% of the best-known solution for the remainder of the problems. Sharing fitness, IEGA, and CM were unable to converge to any of the best-known solutions for the Reeves problems and were within 4% and 15% of the best-known solutions. Similarly, only HIEGA converged to the best-known solution of the Heller problem. However, the proposed algorithm was within 0.7% of the best solution, while sharing fitness, IEGA, and CM were unable to do so.

Mean Value and Standard Deviation of the Optimal Solution
Next, to evaluate the stability of the algorithms, the mean value and the standard deviation of the optimal solution found in the 10 simulations was calculated. Since these were 10 independent simulations, the starting points were randomized in each simulation. These results are shown in Table 3. Table 3. Mean value and standard deviation of the optimal solution found after 10 independent simulations using the algorithms. Even though HIEGA had a better performance than the proposed algorithm when converging to the best-known solution, the proposed algorithm was able to consistently converge to a better solution than HEIGA, as indicated by the mean value and standard deviation. For the Reeves problems, the proposed algorithm had a better performance than HEIGA. For 7 out of 9 problems, its average optimal value was much closer to the bestknown solution than that of HEIGA, IEGA, sharing fitness, and CM. The standard deviation varied between 0.00 and 16.54 for the proposed algorithm, while it varied between 4.64 and 35.10 for HEIGA, 11.74 and 29.53 for IEGA, 15.61 and 39.63 for sharing fitness and 20.34 and 30.15 for CM. The proposed algorithm had the best results for all Carlier problems, as it always converged to the best-known solution. HEIGA had a 100% convergence rate for 4 out of the 8 problems, sharing fitness had a 100% convergence rate to the best-known solutions of problems 1, 2, and 8, while CM and IEGA had a 100% convergence rate to the best-known solutions for 2 out of 8 problems. Lastly, for Heller, the proposed algorithm consistently converged to within 0.7% of the best-known solution. However, HIEGA converged to a better average. These results show the robustness of the proposed algorithm i.e., its ability to consistently converge to a better optimal solution.

Number of Times the Algorithm Converged to the Optimal Solution
To ensure that the algorithms could find multiple solutions consistently, the number of simulations in which the algorithm converged to the best-known solution (N C ) was also recorded. Since none of the algorithms converged to majority of the best-known solution for the Reeves problems, the best optimal solution found was used instead. For example, the best-known solution for Rec-01 was 1247, while the best optimal solution found by the algorithms was 1249, therefore, 1249 was used to calculate N C . The results for N C are given in Table 4. A complete analysis for HEIGA and IEGA could not be done for this metric as these results were not reported by the authors (this is indicated by 'NA' in Table 4). However, some inferences could be made based on the results reported. For example, we could assume that both HEIGA and IEGA always converged to the best-known solution for Car-01 and Car-02, as the mean values were equal to the best-known solution and the standard deviations were zero. Table 4. Number of simulations in which the algorithms were able to converge to the best optimal solution (N C ). Since CM had the worst optimal solution of the three algorithms for the Reeves problems, its N C 's were 0. Similarly, sharing fitness was unable to find the best optimal solution for all but problem 1, and therefore, its N C for those problems was 0. The proposed algorithm had the best performance of all algorithms (where results were available). It had a N C of 10 for problems 7 and 9, 1 for problems 5, 13, and 17, 2 for problems 1 and 15, 4 for problem 11, and 8 for problem 3.

Problem
The proposed algorithm had an N C of 10 for all Carlier problems, indicating that it had a 100% converge rate to the best-known solution for those problems. Sharing fitness also had an N C of 10 for problems, 1, 2, and 8, and had an N C of 4, 9, 3, 3, and 8 for the remaining problems. CM had the worst performance as it was unable to converge to the best-known solution for problems 2, 3, and 4. It had an N C of 2, 6, and 9 for problems 1, 5, and 6, respectively, and had an N C of 10 for problems 7 and 8. Lastly, for the Heller problem, the proposed algorithm had an N C of 10 while the other two algorithms had an N C of 0. These results further solidified the argument that the proposed algorithm had the best performance of the three algorithms, as it was able to consistently converge to the best optimal solution.

Average Number of Best Optimal Solutions Found
Lastly, the average number of best optimal solutions found for each test problem using the different algorithms was compared. For this index, only simulations that had an N C greater than 0 were used. For example, since the proposed algorithm had an N C of two for Rec-01, N O was calculated using the results of those two simulations. N O for the different test problems using the different algorithms are given in Table 5. Since the purpose of HIEGA and IEGA was to only find a single optimal solution, these were excluded for the purpose of this metric. Like the previous three indicators, the proposed algorithm had a significantly better performance than the sharing fitness and CM. Though the sharing fitness found 18 solutions for Rec-01, it was unable to find multiple best optimal solutions for rest of the Reeves problems. The proposed algorithm was able to find 2-6.8 best optimal solutions for all Reeves. CM was unable to find multiple best optimal solutions for any of the Reeves problems, as it never converged to the best optimal solution.
Thought both CM and sharing fitness found more optimal solutions compared to the proposed algorithm for Car-01, their performance was worse for problems 2, 3, 4, and 5. The proposed algorithm found 30.40, 6.20, 13.50, and 1.90 different best-known solutions on average for those problems, while sharing fitness only obtained 9.89, 1.25, 8.67, and 1.00 and the CM obtained 0.00 for all. All algorithms were only able to find 1.00 different optimal solutions for problems 6-8, indicating that the global optimal solution is a unique one. Lastly, for the Heller problem, the proposed algorithm again had the best performance, as it was able to obtain 6.30 different optimal solutions on average, while both sharing fitness and CM obtained 0.00. These results indicate that the proposed algorithm was able to consistently converge to the best optimal solution for all test problems, while finding multiple best optimal solutions.
An attempt was also made to combine the proposed algorithm with sharing fitness, to take advantage of both algorithms, i.e., the ability of the proposed algorithm to consistently converge to the best optimal solution and the ability of sharing fitness to find multiple best optimal solutions. The above task was accomplished by calculating the distance used in the sharing fitness method using the feature matrix of the proposed method. This distance could be calculated using Equation (5).
where F i is the feature matrix of individual i and F j is the feature matrix of individual j.
For the hybrid algorithm, the value of σ share was changed to 10. The performance of the hybrid algorithm was then compared to the performances of sharing fitness. The proposed algorithm and the results for the four indicators are given in Tables 6-9. Table 6. Comparison of the best optimal solution found using the new hybrid algorithm and the previous three algorithms.

Problem Best Optimal Solution Found Using the New Hybrid Algorithm
Best Optimal Solution Found Using Previous Three Algorithms Table 7. Comparison of the mean value and standard deviation of the optimal solutions found using the new hybrid algorithm and from the previous three algorithms.

Best Mean Value and Standard
Deviation of the Optimal Solution Found Using the New Hybrid Algorithm Best Mean Value and Standard Deviation of the Optimal Solution Found Using the Previous Three Algorithms Hel-02 0 10 As can be seen from the results above, the new hybrid algorithm was able to achieve the same best optimal makespan for the Carlier problems but not for the Reeves and Heller problems. In terms of the mean value and standard deviation of the optimal solution found, the new hybrid algorithm had a similar performance to the proposed algorithm for the Carlier problems 1, 4, and 8, while it had a higher mean value and standard deviation for rest of the benchmark problems. When comparing N c , the new hybrid algorithm was able to converge to the best optimal solutions in all 10 simulations for Carlier problems 1, 4, and 8, but was unable to achieve similar results for the remaining problems. In terms of N O , the new hybrid algorithm was unable to obtain more optimal solutions for any benchmark problem. These results indicate that the new hybrid algorithm did not have a superior performance compared to sharing fitness or the proposed algorithm. The reason that the proposed algorithm was able to achieve better performance than sharing fitness, CM, and the hybrid algorithm, was its utilization of a feature vector to separate the individuals into different clusters rather than a single value.
In this section, the performance of the three algorithms was measured by their application to 18 benchmark PFSSP. Four indicators were used to measure the performance-(1) best optimal found, (2) mean value and standard deviation of the optimal solution found in 10 simulations, (3) number of simulations in which the algorithm converged to the best optimal solution, and (4) average number of best optimal solution found. The results indicated that the proposed algorithm had a better performance than the sharing fitness and CM, and also a hybrid algorithm created by combining the proposed algorithm and sharing fitness.

Conclusions
In this study, a new algorithm is proposed for the MMO of PFSSP, i.e., provide multiple solutions with the same objective value. The proposed algorithm utilizes GA for the optimization procedure with the k-means clustering algorithm being used to cluster the individuals of each generation. Unlike standard GA, where all individuals belong to the same cluster, the individuals in the proposed approach were split into multiple clusters and crossover operator is restricted to individuals belonging to the same cluster. Doing so enabled the possibility of the proposed algorithm to find multiple global optima in each cluster. Sharing fitness and CM utilized by Perez et al. for the MMO of JSSP were also modified and used for the MMO of PFSSP. The performance of the three algorithms was tested on various benchmark PFSSP and the results indicated that the proposed algorithm had a better performance than CM and sharing fitness, in terms of its ability to find the best-known optimal solution, consistently provide a better optimal solution, and converge to multiple optimal solutions.
As the sharing fitness method had a better performance on a few problems in terms of finding more optimal solutions, an attempt was made to combine the sharing fitness method with the proposed algorithm. This was done by using the feature matrix of the proposed algorithm in the sharing fitness method. When applied for MMO of benchmark PFSSP, the new hybrid algorithm was unable to outperform or equal the proposed algorithm.
Though the proposed algorithm had the best performance of the algorithms tested, further research needs to be performed in the following areas: 1.
Introduction of more rigorous similarity measure metrics for a more efficient clustering.

2.
Comparison of the performance of the proposed approach to algorithms specifically developed for optimization of PFSSP.

3.
Embedding the proposed approach in algorithms specifically developed for optimization of PFSSP and use it to achieve the MMO of PFFSP.

4.
Applying the proposed approach to real-life operational environment, testing the robustness of the proposed algorithms and its parameters, and measuring its performance under different criteria.
Author Contributions: P.Z. and M.R., analysis of the data, development of code, and writing of the manuscript. S.Y.L., guideline for the proposed approach and comments for the paper. All authors have read and agreed to the published version of the manuscript.

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