2.1. MOFJSP Description

The MOFJSP is described as follows [11]: there are n jobs Ji,i∈1,2,…,n and m machines Mk,k∈1,2,…,m. Job Ji consists of one or more operations Oij,j∈1,2,…,ni, ni is the total number of operations for Ji. Operation Oij is allowed to be executed by one machine from a given set Mij, which is consisted of all the capable machines of the jth operation for job Ji. If Mij=Mk, it is called total flexible, otherwise it is partially flexible. The task of scheduling is to assign each operation to a suitable machine and to sequence the operations on all the machines, so as to optimize some objectives. In addition, some restrictions must be met:

• One machine can only process one job at a time, and one job can only be processed only on one machine at a time.

• There is a priority restriction between every two different operations of the same job, i.e., only after the previous operation is completed, the next one can be processed.

• The operations of different jobs do not have precedence constraints.

• All jobs have equal priorities, all machines are independent and each machine is ready at zero time.

• Once one operation is started, it will not be interrupted until the process is finished.

• Moving time between operations and setting up time of machines are negligible.

The notation used in this paper is summarized in the follows:

In this paper, three objectives will be optimized simultaneously, as follows [12]:

Subject to:

where Eqs. (1–3) are the three optimization objectives. Eq. (1) is minimizing the makespan. Eq. (2) is minimizing the critical machine workload. Eq. (3) is minimizing the total workload of machines. These objectives are in conflict to some extent and have been widely used in many literatures [13]. Eq. (6) is the operation priority restriction. Eq. (7) is the conflict constraint ensuring that each machine can process only one operation at a time. Eqs. (8) and (9) ensure that each operation can be processed only one time.

2.2. Basic Concepts of Multi-Objective Optimization

Multi-objective optimization problem is usually defined in the following form [14]:

s.t. where x is a decision vector, X is the decision space and q is the number of sub-objectives. The objective vector fx is composed of multiple sub-objectives fix.

Existing research approaches for multi-objective optimization problem can be categorized into three: Pareto dominance-based, aggregation-based and lexicographical order-based, and most of the research are on the first two. Due to the superiority of the Pareto dominance-based approach in the number of solutions, and no need to consider the weight assignment problem among multiple sub-objectives, this paper chooses Pareto dominance-based approach to solve the multi-objective optimization problem.

The goal of Pareto dominance-based approach is to obtain all the nondominated solutions, which are referred to as Pareto-optimal solutions. The definition of Pareto dominance and Pareto optimality are given as follows:

1. Nondominated solutions

   Solution a (it is also called nondominated solution) is said to dominate solution b, denoted by a≺b, if and only if

   fia≤fib,∀i∈1,2,…,qfia<fib,∃i∈1,2,…,q.(12)

2. Pareto optimality

   A solution x∗ is called Pareto optimal if and only if

   ∃x∈X,x≺x∗.(13)

3.1. Encoding and Decoding Scheme

The proposed algorithm optimizes assigning sub-problem and sequencing sub-problem simultaneously, and thus the information of operations and related machines should be encoded in the same chromosome. Two encoding schemes are commonly used to solve MOFJSP [15]: A-B string [8,16,17] and 3-tuple scheme [18–20]. Both schemes represent the same information, and the difference is only in the implementation. In 3-tuple scheme, a gene includes three elements which are the job number, the operation number and the assigned machine. However, in A-B string, two genes are needed to represent the same information. So, 3-tuple scheme is simpler and more intuitive than A-B string in the manner of representation, and the chromosome length of 3-tuple is only half of A-B strings. In this paper, a new encoding scheme, called 2-tuple, is designed for MOFJSP, which is improved on the basis of 3-tuple, and it is more concise than 3-tuple because a gene goes from three elements to two.

Each chromosome is a sequence of genes. A gene is a 2-tuple i,k, in which i and k denote job number and assigned machine separately. The jth occurrence of i, from left to right of chromosome, represents the jth operation of job i. Taking Figure 2 as an example, this chromosome encodes the sequence of six operations, three of job1 and three of job2. The first gene 2,1 of chromosome represents that the 1st operation of job2 will be executed on machine M1, the third gene 2,3 denotes that the 2nd operation of job2 will be processed on machine M3, and the rest can be interpreted in the same manner.

Figure 2

Chromosome encoding.

Decoding is the process of converting each chromosome to scheduling solution. Schedules are grouped into three types by [21]: non-delay schedule, semi-schedule and active schedule. It has been proved that active schedule contains the optimal schedule. In an active schedule, no operation can be processed in advance except putting off another operation's start time or changing the order of operations [22]. In order to obtain active schedule, this paper introduces the left-shift greedy decoding scheme [23]. The advantage of this scheme is that it has a self-repairing mechanism, always searching a sufficient and suitable idle time interval for each operation on the processing machine without delaying any other operations, so as to shift the operation to the left as compact as possible. The detailed procedures are shown as follows:

Step 1: Read gene sequence of a chromosome from left to right successively. Assuming a gene i,k is currently obtained, it denotes operation Oij will be executed on machine k.

Step 2: Search the first idle time interval tkS,tkE on machine k, where tkS and tkE represents the beginning time and the ending time of the interval separately.

Step 3: Identify the beginning time of operation Oij. Due to the precedence constrains among operations, operation Oij can be started only after its immediate job predecessor Oij−1 being completed. The starting time Sij of Oij can be calculated through Eq. (14), where Cij−1 is the completion time of Oij−1. If the operation has no predecessor, Cij−1 is replaced with 0.

Step 4: Evaluate the condition of left-shift. If there is enough time from Sij to tkE, i.e., Sij+pijk≤tkE, where pijk is the processing time of Oij on machine k, then the interval tkS,tkE is available for Oij. That is, Oij can be left shifted, and the completion time of Oij is Sij+pijk. Otherwise, search the next idle time interval. If there is no interval satisfy the condition of left-shift, then Oij will be appended in the rear of the sequence of operations that has already been scheduled before Oij on machine k.

Step 5: Update chromosome. Due to the employ of the left-shift greedy decoding scheme, operation r may be processed earlier than another operation v, which appears before r in the chromosome. Hence an update operation should be performed according to the result of decoding instead of the original chromosome, so as to keep the uniformity of encoding and decoding.

3.2. Population Initialization

The quality of initial population is very important for accelerating convergence speed and avoiding premature convergence. It has been proved that utilizing a hybrid multi-approach to generate initial population is better than using single one in the performance [20]. Thus, in this paper, an approach of integrating strategies is employed for generating initial population. It includes five assignment rules and four sequencing rules.

1. Assignment rules

   • Random rule: It assigns each operation to a capable machine randomly.

   • Minimum processing time rule [20]: For every operation, the machine with the minimum processing time will be selected in the corresponding operation's candidate machine set. If there are multi-machines with the same minimum processing time, select one from them randomly, then assign the operation to the selected machine.

   • Global minimum processing time rule [20]: From the processing time table, find the minimum processing time among all operations, fix the assignment, then add the selected processing time to every other entry in the same column, and conduct the next time selecting, until all operations are assigned. If there are multi-global minimum processing time, select one randomly.

   • Local minimum processing time rule [23]: Find the minimum processing time from the first operation of the first job to the last operation of the last job successively. For each operation, the selected time will be added to other entry in the same column for conducting next time selecting, and each operation will be assigned to the corresponding selected machine. If there are multi-local minimum processing time, select one randomly.

   • Permutation rule: Shuffle the order of operations in processing time table randomly. Take the operation in the first row as the starting operation to select the machine with minimum processing time, fix the assignment, then add the selected processing time to every other entry in the same column for conducting next time's selecting. If there are multi-machines have the same minimum processing time, randomly select one.

2. Sequencing rules

   • Random rule: It places the operations into a chromosome in a random order.

   • Most work remaining rule [6]: One operation is selected according to the remaining work time of all jobs and be placed into a chromosome. The larger the remaining work time is, the earlier the operation of corresponding job can be selected. If there are multi-jobs have the same remaining work time, select one randomly.

   • Most number of operations remaining rule [20]: This rule select operation according to the remaining number of operations, the more the remaining number of operations is, the earlier the operation of corresponding job can be selected. If there are multi-jobs have the same remaining number of operations, select one randomly.

   • Shortest processing time rule [6]: The first operations of the remaining operations from each job are compared, the operation with the minimum processing time will be selected preferentially.

In the researches that utilizing hybrid multi-strategies to generate initial population, commonly a percentage-based approach, allocating different percentages for each approach to generate the required initial population, is adopted. In this paper, a new loop-based method is designed to generate initial population, and the difference is percentage allocation is not needed for each strategy compared with some other excellent algorithm, e.g., EPABC proposed in [13], AIA proposed in [18], etc. The utilization of loop-based method helps to increase the diversity of population and to extend the search space. The detailed procedure is described as follows:

Step 1: Place 20 5×4 combinations of the assignment and sequencing rules in the set CH.

Step 2: A combination is successively selected in the manner of loop from CH to generate an individual. Assuming that the combination selected is A,S. First, employ strategy A to complete operation assignment. Second, utilize strategy S to complete operation sequencing. Finally, combine the result of operation assignment and operation sequencing to generate a chromosome.

Step 3: Repeat perform step 2 until the number of initial individuals reaching Npop.

3.3. Local Search I (LS-I)

Local search is an effective method to exploit better solution around a current solution. It has been widely used in MOFJSP. Inspired by the idea of left-shift greedy decoding scheme (mentioned in Section 3.1), also for keeping initial population's distribution characteristics, a new local search method is proposed. The performing process of LS-I is described as follows:

Step 1: Take out one chromosome from population, decode the chromosome, obtain three objective's value F1,F2 and F3 and the terminal completion time Ck,k∈1,2,…,m of each machine.

Step 2: Decode again. Due to different operations of the same job have precedence constrains and parent population commonly is not the optimal solutions, there are some idle time interval existing in a scheduling solution. Read a gene q from chromosome, find the first interval TuE,TvS which is located between adjacent operations u and v on the same machine, where TuE is the end processing time of u, TvS is the start processing time of v. Let TqeS represents the earliest start processing time of operation q after moving q to the interval TuE,TvS. If TqeS=maxTq−1E,TuE and TqeS<TvS, q has qualification to be moved to the interval TuE,TvS and step 3 is performed. If q cannot be moved, next interval is examined. If all intervals do not satisfy the condition, decode the next gene. In case q has no job predecessor, Tq−1E is replaced with “0”.

Step 3: Judgment. If the moving of operation q results the solution become better, i.e., the new solution dominance the old one, update the chromosome and decode next gene. If the new solution is equal to the old one, but at least one machine's maximum completion time Ck become short, update the chromosome and decode next gene, so as to make the scheduling as compact as possible. Otherwise, q cannot be moved.

Step 4: Repeat step 2 to step 3, until all genes on the chromosome are decoded. Repeat step 1 to step 3, until all the individuals in initial population are optimized.

In order to better illustrate LS-I, taking Figure 3 as a simple example, moving operation O12 to the interval between O41 and O33 will shorten the processing time of machine M3. If it is critical operation, the maximum completion time in all machines will be shorten, if it is not critical operation, the maximum completion time of current machine will be shortened and the scheduling will be more compact.

Figure 3

Local search before evolution operation.

3.4. Crossover Operators

The goal of crossover is to obtain two offsprings which inherit the gene information from their parent's chromosomes and generate new superior gene sequences. In this algorithm, crossover includes two parts: crossover for operation sequence and crossover for machine assignment. These two parts are implemented separately, and the terminal condition is that the amount of offspring equals Npop, where Npop is the amount of the initial population.

1. Crossover for operation sequence

   Many effective approaches are developed for operation sequence's crossover in the past decades, such as order crossover employed by [16], HCO proposed by [12], IPOX designed by [11]. In this paper, considering the chromosome structure, IPOX is adopted. In addition, a new crossover approach, called RPOX, is proposed on the basis of IPOX, so the proposed algorithm can adaptively select one from the two crossover approaches according to the quality of the parent's chromosomes to generate good offspring. The details of IPOX and RPOX are shown as follows:

   Step 1: Fast nondominated sorting [24] is performed for initial population which is divided into different ranks Rαα∈1,2,…,l,…,r. If ∑α=1l−1NRα<0.7Npop and ∑α=1lNRα≥0.7Npop, the individuals in R1 to Rl are called excellent population and the remained individuals are called inferior population, where NRα represents the amount of individuals in rank Rα.

   Step 2: Randomly select two individuals P1 and P2 from initial population as parent's chromosomes, and let C1 and C2 be their children's chromosomes by crossing. All the jobs are randomly divided into two sets J1 and J2.

   Step 3: Copy the elements of P1 that are included in J1 to C1 in the same position and copy the elements of P2 that are included in J2 to C2 in the same position.

   Step 4: If the elements of P2 that are included in J2's rest position in the same order and the elements of P1 that are included in J1's rest position in the same order, this method is called IPOX. If the elements of P2 that are included in J2's rest position in the reverse order and the elements of P1 that are included in J1's rest position in the reverse order, this method is called RPOX. At this moment, keep the machine associated with operation is empty. A simple example about these two crossover methods is shown in Figure 4.

   Step 5: If the parent's chromosomes are in the inferior population, RPOX is adopted. Otherwise, IPOX is adopted. Because the chromosomes selected from the inferior population always possess few good gene information, only with the large change, the offspring can possibly generate useful gene sequence. On the contrary, the chromosomes selected from the excellent population possess good gene sequence, so small change should be made to keep the good gene information to the next generation.

2. Crossover machine assignment

   After operation sequence's crossover is completed, crossover for machine assignment should be performed to determine the machine for each operation in children's chromosomes C1 and C2. Two-point crossover method, called ITP, is used, and the details are described as follows:

   Step 1: Two integers I1 and I2 are selected from the range of 2,L−1, where L is the length of one chromosome.

   Step 2: The genes of C1 between position I1 and I2 inherit the operation's machine information from P1, the remined genes inherit from P2. The genes of C2 between position I1 and I2 inherit the operation's machine information from P2, the remined genes inherit from P1. A simple example is shown in Figure 5.

Figure 4

IPOX and RPOX crossover operation for the operation sequence.

Figure 5

ITP crossover operation for the machine assignment.

3.5. Mutation Operators

Mutation is one of the key technologies for avoiding premature and increasing the diversity of population. A substitution method is adopted in this paper. For operation sequence's mutation, two genes that do not belong to the same job are selected randomly, then their position will be swapped and the corresponding machine information will be updated.

For machine assignment's mutation, the machine Mi in a randomly selected gene is replaced by a new one that is reselected from its candidate machine set. This reselection must meet a rule. If a subset, consisted by the machines with the processing time shorter than Mi, exists, select the new machine from this subset. Otherwise, select randomly from the candidate set but except Mi. The terminal time is controlled by mutation rate β.

3.6. Environment Selection

Environment selection is to select certain amounts of individuals, called new initial population or survived individuals, from the population of parent and children. The survived individuals will be participated in the next time's genetic operation. For environment selection, fitness-based [7] and dominance-based [25] selection methods have been developed and widely used in scheduling problem. Due to the advantage of dominance-based method in term of considering the tradeoff among multi-objectives, it is adopted in this paper.

In dominance-based method, crowding distance [25] is commonly used to select certain individuals from the subpopulation with the same rank. It is effective to the optimization problem with two objectives, but with three objectives it is insufficient to some extent. In order to keep the diversity of the new initial population, a partition selection method is designed to replace crowding distance. The details of environment selection are described as follows:

Step 1: Fast nondominated sorting is performed for all individuals and divided them into different levels Rαα∈1,2,…,l,…r.

Step 2: If ∑α=1l−1NRα<Npop and ∑α=1lNRα>Npop, go to step 3. If ∑α=1lNRα=Npop, all individuals in R1 to Rl will survive to the next generation, and go to step 4. Here NRα is the number of individuals in rank Rα.

Step 3: All individuals in R1 to Rl−1 survive to the next generation directly. The task of environment selection is to determine the survived individuals in Rl. Let NRl represents the number of individuals need to be selected from NRl, where NRl=Npop−∑α=1l−1NRα. First, if two objectives are called a sort operator, six sort operators can be consisted by three objectives F1,F2 and F3, and it is denoted by number 1 to 6 separately. For example, F1,F2 is denoted by number 1. Second, a number is selected randomly from 1 to 6. Assuming the selected number is 1 and it represents the sort operator F1,F2. Sort all individuals of Rl in ascending order based on F1, if F1 of multi-individuals are the same, sort the multi-individuals based on F2. Third, according to the equation nl=∧NRl/NRl, divide NRl individuals into NRl parts, each part except the last one has nl individuals, where ∧ presents taking lower bound value. For example, if NRl/NRl=2.6, nl is 2. The individual number of the last part may be larger than nl. The individual in the first position of each part will be selected to the next generation.

Step 4: An external AS is imbedded in the proposed algorithm. The individuals in R1 will be added into the AS when the environment selection procedure is performed at first time. Otherwise, this step can be negligible.

3.7. Local Search II (LS-II)

The purpose of LS-I mentioned in Section 3.3 is to make the chromosome carry more superior genes in a short time. However, a new local search method, called LS-II, based on the critical path is purposed to find more promising solutions. The emphasis of the two methods are different, LS-I pursues speed, LS-II pursues quality.

In this paper, a concept, called full sequence, is proposed. If no time interval exists between every two adjacent operations on the same machine, this state is called a full sequence, e.g., M2 is under a full sequence in Figure 3. If a machine is under a full sequence and makespan belongs to this machine, it cannot reduce the makespan no matter how to shift the critical operation, unless reassign some operations to other machines.

In order to optimize three objectives simultaneous, the local search method based on critical path is extended. It includes two parts: reassignment and shifting critical operations. Compared with some methods that only shifting critical operations, LS-II not only reduce the makespan (F1), but also balance the critical machine workload and the total workload F2&F3. The details of LS-II are described as follows:

Step 1: According to the result of fast nondominated sorting in Section 3.6, the individuals in R1 will be local searched. The goal of this is to reduce the load of calculation.

Step 2: Take out a chromosome c, exam whether the chromosome has the full sequence and the full sequence belongs to the machine with makespan. If no, go to step 3, else, go to step 4.

Step 3: Shifting critical operation. First, identify the critical paths and the critical blocks on the chromosome. Let Oij represents one critical operation of a critical block, and Oij cannot be the head operation of this critical block. Left-shift Oij to the front of every operation in the critical block successively and calculate the corresponding objective value. If the solution is better, the chromosome after shift and this time's shift will be recorded. Repeat this process, until every critical operation except the head one in every critical block are performed. Multi-critical paths may possess some public critical blocks, so if recorded shift appears, it will be negligible. After shifting critical operations is performed, a new population Ns based on c is formed, the best one in Ns will be selected to replace c. This step can realize reducing the makespan (F1).

Step 4: Reassignment. In all full sequences, if an operation has at least one candidate machine and the candidate machine's processing time is not larger than the current one, this operation is called reassignment operation and these candidate machines compose the replacement machine subset. Select one reassignment operation randomly, select the machine that has the minimum workload from the replacement machine subset, and move the selected operation to the selected machine. If some machines have the same minimum total workload, select one randomly. Repeat performing reassignment until the number of performing equals to the total sum of reassignment operations. If the new chromosome is better than the old one, the chromosome will be recorded. The best one in all new chromosomes will replace the old chromosome. This step can realize reducing the makespan (F1), balancing the critical machine workload and the total workload F2&F3).

Step 5: Repeat step 1 to step 4, until all the individuals in R1 are performed. Compare the individuals in R1 to AS. If some individuals in AS are dominated by one individual in R1, they will be removed and the better individual will be added. If one individual in R1 and individuals in AS are not in the dominance relationship, the individual in R1 will be added into AS.

3.8. Main Algorithm

The framework of the proposed algorithm for MOFJSP can be shown in Figure 6.

Figure 6

The flowchart of the proposed algorithm.

Step 1: Initialization. Loop-based method (Section 3.2) is used to generate initial population Npop.

Step 2: Apply LS-I (Section 3.3) to optimize all individuals in Npop. Update Npop to Newpop.

Step 3: Perform crossover operators (Section 3.4) and mutation operators (Section 3.5) to generate offspring populations Ncpop and Nmpop.

Step 4: Merge Newpop, Ncpop and Nmpop. Remove duplicate individuals. Perform environment selection (Section 3.6) to generate the survive individuals Nspop.

Step 5: Apply LS-II (Section 3.7) to search more promising solutions, update Nspop and AS.

Step 6: If the terminating condition is satisfied, the algorithm ends. else, go to step 2.

4.1. Parameter Setup

The process of determining the feasible parameter is complex in some researches. To simplify the process, in our algorithm, very few parameters are needed, and an adaptive parameter setup method is designed. The amount of initial population Npop is determined by Eq. (15), where n is the amount of jobs. The mutation rate β is calculated by Eq. (16), where iter is the number of iterations. The terminal condition is that the AS members remain unchanged after 30 consecutive times iteration.

4.2. Result and Comparison

To test the performance of the proposed algorithm, the algorithm procedure was implemented in python and run through the PyCharm software on a PC with AMD A-8 6410 2.0 GHz CPU 4.0 G RAM and Windows 10 operating system. Four representative instance sets based on the practical data have been selected. The first one is taken from [26], which includes 10 small size problem instances (SFJS01:10) and 10 medium and large size problem instances (MFJS01:10). The second one is taken from [27] without release time, which includes five problem instances (problem 4×5, 8×8, 10×7, 10×10, 15×10), where problem 4×5 denotes four jobs and five machines, other four problem instances are interpreted in the same manner. The third one is taken from [6], which includes ten problem instances (Mk01-10). The last one is taken from [28], which includes eighteen problem instances (01a–18a). These four instance sets are the most commonly adopted benchmark instances in the researches on MOFJSP.

The algorithm's advantage is that very few parameters setup is needed. Table 1 gives a comparison of several state-of-the-art algorithm's number of parameters, including Li et al. [17,29], Wang et al. [11] and Xing et al. [30]. From this table, the number of parameters needed in the proposed algoithom is the least, which is three.

In order to test the effect of improved coding, crossover, mutation, population generation and environment selection methods on the performance of the proposed algorithm, the two-stage local search was removed from the proposed algorithm, and the remaining part (IEEA) was tested together with the traditional EA through the same case. The proposed algorithm is called IEA. The test case set uses SFJS01-10 and MFJS01-10. The EA is selected from the literature [15], which uses a variety of combination strategies, and has proved to be more superior than other classic EAs. In this section, metric C is used. CA,B measures the fractions of members of B that are dominated by members of A, and Eq. (17) is the calculating equation of CA,B. CA,B=0 means that all solutions in B are not dominated by any solution in A, CA,B=1 means that each solution in B is dominated by some solutions in A. The test results are shown in Table 2.

It can be seen from Table 2 that in the 20 problem instances, IEEA is superior to EA except MFJS03 and MFJS07. Compared with IEEA, IEA does not show obvious advantage in 10 small size problem instances, indicating that for small size problem instances IEA can find the optimal solution without local search, while for 10 medium and large size problem instances, IEA has shown obvious advantage. In addition, the frontier of EA, IEEA and IEA on MFJS02, MFJS03, MFJS05 and MFJS10 are given, as shown in Figure 7. It can also be seen from the figure that IEEA is better than EA, and IEA is better than IEEA.

Figure 7

The frontier of EA, IEEA and IEA on problem instances MFJS02, MFJS03, MFJS05 and MFJS10.

In order to test the overall performance of the proposed algorithm, it was tested on other three instance sets with several state-of-the-art algorithms.

For five problem instances, four state-of-the-art algorithms in Table 1 are also used to compare with the proposed algorithm. The reason of doing this is that the four algorithms are powerful, and it has been proved that it is better than other well-known algorithms (e.g., PSO+SA proposed by [7], AL+CGA proposed by [31], etc.). Hence, if the proposed algorithm is better than the four powerful algorithms, it is also better than PSO+SA and AL+CGA, etc.

Before comparison, some terms should be clarified. Optimal solutions in this paper refer to the nondominated solutions obtained by each algorithm. Pareto front is constructed by the optimal solutions of all compared algorithms and Pareto solutions refer to the solutions in Pareto front. All solutions in the Pareto front are called overall Pareto solutions. If some optimal solutions in a compared algorithm's optimal solutions are the Pareto solutions simultaneously, they are called partial Pareto solutions.

First, compare the amount of optimal solutions. Table 3 is the result of optimal solutions' amount. From this table, the proposed algorithm can obtain the best results for all problem instances. In problem 4×5, the amount of optimal solutions obtained by the proposed algorithm is more than other four. In problem 8×8, the proposed algorithm obtains more optimal solutions than others, and it equals to Xing. However, one sloution (16, 13, 73) in the proposed algorithm's result dominates one solution (17, 13, 73) in Xing's. In problem 10×7, the proposed algorihtm obtains more optimal solutions only than Li, equals to other two. But one solution (11, 11, 61) in the proposed algortithm dominates two solutions (12, 11, 61) and (11, 11, 63) of Li and Pan. In problem 10×10, Li and Pan, Li and Xing all obtain three optimal solutions, the proposed algorithm can achieve four solutions. Compared with Wang, though the optimal solutions' amount is the same, one solution (7, 5, 45) in Wang is dominated by one solution (7, 5, 43) in the proposed algorithm. In problem 15×10, the proposed algorithm's optimal solutions amount is only more than Xing, is fewer than Wang, and equals to Li and Pan and Li. But one solution (11, 11, 91) in the proposed algorithm dominates all solutions in Li and Pan, and one solution (11, 10, 93) dominates two solutions (12, 10, 95) and (11, 10, 98) in Wang. Therefore, it can be concluded that our algorithm has a strong search ability on MOFJSP.

Second, compare the quality of optimal solutions. Table 4 gives an overall perspective of optimal solutions' quality. In this table, H denotes the number of the overall Pareto solutions, N represents the number of optimal solutions, M represents the number of partial Pareto solutions. Q, calculated by Eq. (18), represents the percentage of partial Pareto solutions M in optimal solutions N. R, calculated by Eq. (19), represents the percentage of partial Pareto solutions M in overall Pareto solutions H. For example, in problem 8×8, the proposed algorithm obtains four optimal solutions (i.e., N=4), and the Pareto front includes five solutions (i.e., H=5), all optimal solutions are Pareto solutions (i.e., M=N=4), the percentage Q of partial Pareto solutions M in the optimal solutions N is 100%, and the percentage R of partial Pareto solutions M in overall Pareto solutions H is 80%.

From Table 4, all the optimal solutions obtained by the proposed algorithm in the five problem instances are Pareto solutions, and, in each problem instance except problem 8×8, the number of the optimal solutions are equal to the number of overall Pareto solutions. Compared with other four algorithms, only Xing in problem 10×7 and Li in problem 15×10 achieve the same effectiveness. The last row gives a statistic result in the term of H, N, M, Q and R. In order to best show the statistics result, corresponding value is made into Figure 8. From Figure 8 the proposed algorithm can obtain the best result both in solutions' amount and solutions' quality.

Figure 8

Statistic result on N, M, Q and R.

Table 5 gives a result from a local perspective, i.e., taking every one in the four state-of-the-art algorithms to compare with our algorithm. In this table, metric C is used. The proposed algorithm is denoted by B.

From this table, all optimal solutions obtained by the proposed algorithm are not dominated by any optimal solution obtained by other compared algorithms in the five problem instances. However, at least one optimal solution obtained by the proposed algorithm can dominate some optimal solutions of Li and Pan (A1), Wang (A3) and Xing (A4) in at least one problem instance. Though, the proposed algorithm's optimal solutions are nondominated with Li's, but from Tables 3 and 4, the proposed algorithm can achieve more optimal solutions. So, it can be also concluded that the proposed algorithm is more superior than others compared algorithm.

In addition, the distribution and convergence of optimal solutions is analyzed. One metric called inverted generational distance (IGD) is used. It is calculated by Eq. (20), where P∗ denotes the Pareto front, which consists of optimal solutions, |P∗| denotes the number of Pareto solutions, dx,y is the Euclidean distance between points x and y. The smaller IGDA,P∗ is, the closer A is to the Pareto front. Besides that, the objectives of all chromosomes are normalized by Eq. (21), where fimin and fimax are the minimum and maximal values of the ith objective among all solutions. Table 6 is the calculating result. From this table, the proposed algorithm is closer to Pareto front than other compared algorithms.

For MK01-10, two state-of-the-art algorithms, Bagheri et al. [18] and Xing et al. [30] are selected to compare with the proposed algorithm. The reason of selecting these two algorithms is that they have strong searching capability and give three objective values rather than other popular algorithms' only giving one value. Because the number of nondominated solutions obtained by each algorithm is usually large, it is difficult to judge the algorithms' superiority by the number of solutions. In order to verify the superiority of the proposed algorithm, like the methods used in other papers, this article only compares optimal solution with minimum makespan.

Table 7 lists the optimal solutions with the minimum makespan of each algorithm. In this table, LB is the best-known lower bound of makespan [32]. Dev is the relative deviation, calculated by Eq. (22), where Cpro is the best makespan obtained by the proposed algorithm, Cpar is the best makespan of the algorithm that compares to the proposed algorithm.

According to Dev of this table, the proposed algorithm's result is better than Bagheri et al.'s in five problem instances (MK05, MK06, MK07, MK09 and MK10) and is better than Xing et al.'s in eight problem instances (M01-MK10 except MK03 and MK08). For all optimal solutions of each algorithm in ten problem instances, metric C is used to evaluate the performance of the proposed algorithm, and the calculating result is shown in Table 8, where P represents the proposed algorithm.

From Table 8, eight optimal solutions obtained by Bagheri et al. are dominated by some optimal solutions obtained by the proposed algorithm, and no optimal solution of the proposed algorithm is dominated by Bagheri et al.'s. Three optimal solutions obtained by Xing et al. are dominated by the proposed algorithm's, but only two of the proposed algortihm's are dominated by Xing et al.'s.

In order to further test the performance of the proposed algorithm, the 01a–18a problem instances were used. Two state-of-the-art algorithms, MG designed in [11] and BEG proposed in [33], are selected to compare with the proposed algorithm, because these two algorithms list all solutions in their literature. Metric C is used to evaluate the performance of the proposed algorithm and the calculation results are shown in Table 9.

Compared with MG, IEA has 10 instances better than MG, and MG has 3 better ones in all 18 instances. Compared with BEG, IEA has 16 better ones, while BEG has none. Therefore, it is also concluded that the proposed algorithm has superior solution performance.