An improved memetic algorithm for the flexible job shop scheduling problem with transportation times

In the practical production, the transportation of jobs is existed between different machines. These transportation operations directly affect the production cycle and the production efficiency. In this study, an improved memetic algorithm is proposed to solve the flexible job shop scheduling problem with transportation times, and the optimization objective is minimizing the makespan. In the improved memetic algorithm, an effective simulated annealing algorithm is adopted in the local search process, which combines the elite library and mutation operation. All the feasible solutions are divided into general solutions and local optimal solutions according to the elite library. The general solutions are executed by the simulated annealing algorithm to improve the quality, and the local optimal solutions are executed by the mutation operation to increase the diversity of the solution set. Comparison experiments with the improved genetic algorithm show that the improved memetic algorithm has better search performance and stability.


Introduction
The flexible job shop scheduling problem (FJSP) is an extension of the classic job shop scheduling problem (JSP), which is more in line with the characteristics of discrete manufacture environment. With the development of customization and product quality, the production process will be more flexible. The control of various time factors also needs to be more accurate. In actual production environments, the time factors include not only the processing time, but also the transportation time. These time factors affect the completion time and make the problem more complicated. If this part of the time is not taken into account, it may occur that the processing time is short but the transportation time is long during implementation, resulting in reduced scheduling efficiency. Therefore, a flexible job shop scheduling problem with transportation times (TT-FJSP) is studied in this paper.
The FJSP was first addressed by Brucker and Schile, they solved the problem by polynomial graphics algorithm. 1 Nevertheless, many experiments show that the exact algorithm is not very effective even in solving moderate-size FJSP instances. 2 Currently, the general method for solving the FJSP is to first establish a optimization model and then solve it by the intelligent optimization algorithm. The mixed integer model is one of the common used optimization model. Fattahi and colleagues 3,4 established a mixed integer linear programming model with location attributes for the FJSP with two or more jobs and the FJSP with some constraints between the operations. Nowadays, manufacturing has become refined and relies on intelligent solutions. Scheduling problem exists widely in actual production and manufacturing, and many solving methods are applicable to actual production to provide efficient and feasible solutions for production. With the development of problems, relevant optimization models and theories emerge in endlessly, 5 and the value of practical application research is growing gradually. [6][7][8] The transportation time is an important factor of the actual schedule, which affects the final scheduling results. Naderi et al. studied the displacement flow workshop that takes into account the transportation time of semifinished products, and they established different mixed integer programming models for multi-transporter and single-transporter shops. 9 For the manufacturing system with resource flexibility, sequence-dependent setup, and transportation times, Rossi proposed a swarm intelligence approach to schedule. 10 Karimi et al. 11 proposed an imperialist competitive algorithm to solve the TT-FJSP. Zhang and colleagues 12,13 successively used improved genetic algorithms (GA) solved the TT-FJSP and multi-objective FJSP with multiple time constraints.
According to the search scope, the intelligent optimization algorithms are mainly divided into two categories: global search and local search. The former has a large search range, such as GA, [14][15][16] particle swarm optimization algorithm, 17 ant colony optimization algorithm, 18 and artificial bee colony algorithm. 19 The local search algorithm is the neighborhood search algorithm for solutions, such as the simulated annealing algorithm (SA) 20 and tabu search. 21 Nonetheless, some experimental studies showed that the hybrid algorithms have better results. 22,23 The memetic algorithm (MA) is proposed by Moscato 24 to simulate the evolution of human culture. It is essentially a combination of the global population search and local search. Generally, the MA is similar to the GA, but the mutation part is replaced by a local search. Some scholars use the MA to solve the FJSP. Phu-ang and Thammano 25 presented a honeybee marriage-based MA to avoid falling into local optimal solutions. Gong et al. 26 used an MA for solving the multi-objective FJSP with worker flexibility.
In this paper, an improved memetic algorithm (IMA) is presented for solving the TT-FJSP. In the IMA, an effective SA with the elite library and mutation operation is adopted to form the local search method with parallel structure. To avoid repetitive searches and improve efficiency, the elite library is adopted to maintain outstanding solutions. Feasible solutions, according to the mutation probability, are divided into general solutions and local optimal solutions. General solutions are executed by the SA, and local optimal solutions are executed by the mutation operation. The proposed algorithm is evaluated and compared with the improved genetic algorithm (IGA), 12 and the comparison shows that the IMA is effective for the TT-FJSP.
The rest of this paper is organized as follows. The problem description and model are given in ''Problem formulation'' section. ''The improved memetic algorithm'' section formulates the IMA, and the experimental study results are given in ''Experimental results and comparisons'' section. The conclusions and some ideas for the future research are summarized in the final section.

Problem formulation
The FJSP can be described as follows: There are some jobs and machines. Every job has its own operation set, and the operations are processed under their process constraints. Each operation can be processed by multiple machines, and the processing time on different machines is different. The FJSP is generally divided into the fully flexible FJSP and partially flexible FJSP. 27 In the fully flexible FJSP, each operation can be processed on all machines. In the partially flexible FJSP, each operation can be processed on specific optional machines. A partially flexible FJSP instance is shown in Table 1.
The TT-FJSP has one more transportation times constraint: When one operation is completed, the finish time is not equal to the start time of the next operation, but it needs to determine whether the job needs to be transferred to another machine. If necessary, the transportation time of the job between the two machines should be considered, otherwise the transportation time is 0. The objective is the minimization of the makespan when the processing times and transportation times are considered simultaneously. In the TT-FJSP, several constraints and assumptions are made: Jobs and machines are independent and available at time 0. Job preemption is not allowed, and each machine can handle only 1 job at a time. The operations belonging to the same job cannot be processed simultaneously. Processing time of each operation corresponds to its processing machine, and all processing times are given. Once the operation is completed, the job is transported to the machine of the next operation immediately, and the corresponding transportation time is spent. Transportation time is related only to the selected processing machine and direction, and all transportation times are given. Transportation resources are sufficient, which do not affect the transportation process. Setup times are included in the processing times.
The objective is to minimize the makespan C max . Equation (2) represents the processing sequence constraint of operations. Equation (3) expresses that the makespan is not less than the completion time of each job. Similarly, equation (4) means that the completion time of J i is not less than the finish time of O i,h on M k . Equation (5) indicates that the finish time equals the sum of the start time and processing time of the operation. Equation (6) denotes that the processing of jobs processing starts at time 0. For example, we used the data in Table 1 in our tests. The transportation times are given in Table 2, which shows the time taken to transport from the left corresponding machine to the next corresponding machine. The test result with transportation times is shown in Figure 1(a), and the test result without transportation times is shown in Figure  1(b). In Figure 1(a), the two light yellow rectangles represent the transportation times of J 1 and J 2 , and it takes MT 2,1,4 to transport J 2 from machine M 1 to M 4 , so S 2,3,4 is equal to F 2,2,1 plus MT 2,1,4 . However, S 2,3,4 is equal to F 2,2,1 in Figure 1

The improved memetic algorithm
The JSP has been proved to be an NP-hard problem, 28 while the TT-FJSP is an extension of the JSP and is also an NP-hard problem. This means that it may be difficult to get an optimal solution in finite time through exact solution methods. The MA is a meta-heuristic algorithm with both exploration and exploitation, which is suitable for this problem. The main processes of the IMA are as follows: Step 1: Generate initial population according to encoding rules.
Step 2: Generate a new population for the next generation by the selection and crossover.
Step 3: Apply the local search and save the high-quality individuals.
Step 4: Check the termination criteria. If the criteria are not satisfied, then stop the algorithm and output the optimal individual in the population. Otherwise, go to Step 2.
The framework of the IMA is given in Algorithm 1. The encode and decode, initial population, selection, Table 2. Transportation times of the instance. crossover, and local search are introduced in the following subsections.

Encode and decode
Encode is the first important step in the algorithm. Encode is the process of expressing the solutions of a problem in the form of chromosomes, or as individuals.
Then the genetic operations, such as the selection and crossover, can be implemented. The TT-FJSP contains two subproblems: machine selection and operation sequence. We use the twosection integer encode method to address these issues. The chromosome structure is divided into two sections: machine selection and operation sequencing. The machine selection determines the machines of all operations, and the operation sequencing determines the processing order. The lengths of these two sections are equal and equal to the total number of operations Machine selection. In the machine selection, each gene position corresponds to each operation in sequence, and an integer is used to represent the sequence number of the optional processing machine of the operation, instead of the processing machine number. For example, according to Table 1, if the gene string is 4-1-2-3-3, the first gene position 4 indicates that the operation O 1,1 is processed on the fourth machine in the candidate set of O 1,1 , that is M 5 . The machine of O 1,2 is the first machine in the candidate set of O 1,2 , that is M 2 .
Operation sequencing. The operation sequencing is similar to the chromosomes in the JSP. The encode of this section is based on the processing order of the job. Each gene is the job number and the number of each operation in a job is indicated by the number of occurrences in the operation sequencing section. If the gene string is 2-1-1-2-2, then the corresponding operation

Decode
There are three types of scheduling including active scheduling, semi-active scheduling and non-delay scheduling. 27 It has been proved that the optimal scheduling exists in the active scheduling set. That is, when chromosomes are decoded as actively scheduled, the result is closest to the optimal solution. Therefore, the plug-in decode method is used to decode chromosomes into the active scheduling. Each chromosome contains two sections, so it needs to be decoded separately. . Jm(i,h) and T(i,h) have a 1-to-1 correspondence. According to the instance above, the first operation of J 1 is processed on M 5 , the second operation of J 1 is processed on M 2 (as shown in equations (7) and (8)), and the processing time is 5. Input: parameters for the IMA and an TT-FJSP instance Output: the optimal individual S opt 1 Iter = 0; 2 Set = generating initial population by the initialization method; 4 cal(Iter) = calculate the objective function value of the population in Set; 3 while Failure to satisfy the termination condition 5 Set = selection for individuals in Set according to cal(Iter); 6 Set = crossover for individuals in Set; 7 Set = local search for individuals in Set; 8 Iter = Iter + 1; 9 cal(Iter) = calculate the makespan of the individuals in Set; 10 S(Iter) = select the optimal individual in Set according to cal(Iter); 11 end 12 S opt = select the optimal individual in S; the S i,h,k satisfying the processing sequence constraint (it can be obtained from equation (10). Otherwise, we can get the S i,h,k with equation (11). When there are several idle time periods on a machine, the insertion is attempted sequentially from the beginning to the end.

Population initialization method
The initial population can be obtained through a variety of methods. In general, the initial population method affects the convergence speed and the quality of the population. We use the random method to highlight the performance of other operations of the algorithm. The main steps are as follows: Step 1: Randomly generate a chromosome to join the initial population: In the machine selection part, randomly select one from the optional machine at each gene position, and in the operation sequencing part, all operations are randomly arranged.
Step 2: Repeat Step 1 until the initial population is filled.

Selection operator
Selection is for good individuals to survive with a greater probability and avoid the disruption of good genes by operations such as the crossover and mutation. The tournament selection method was adopted as the selection method for the proposed algorithm. This method randomly selects some individuals from the population for comparison of fitness each time, then puts the highest fitness individual into the crossover pool and waits for the crossover, until the crossover pool is full.

Crossover operator
An excellent crossover can preserve the good genes in the paternal generation by exchanging information among the paternal individuals to produce good new individuals. It can also reduce a blind search and achieve a simple and efficient search. We propose a crossover method for chromosomes in the FJSP, which includes two different crossover modes.
Machine selection. A multi-point crossover is adopted to ensure that the chromosomes of solutions generated by crossover are still feasible. That is to say, the random selection of multiple intersection sites occurs, and then the exchange of genes on the same location of two parents occurs.
Operation sequencing. The multi-job crossover is adopted to avoid generating infeasible solutions. This method selects all the genes of several jobs in the parent chromosomes and crosses them in sequence. The pseudocode of the multi-job crossover is given in Algorithm 2.
In fact, the crossover of the two parts is done simultaneously. When the crossover operator is executed, two individuals are first randomly selected from the population as parents. Then the machine selection and the operation sequencing of the offspring are generated by the multi-point crossover and the multi-job crossover, respectively. This process is repeated until the number of the offspring generated reaches the population size.

Local search
The IMA applies an independent local search (LS) to try to improve the quality of the solutions. The LS finds the local optimal solutions by searching the neighborhood solutions. The LS can be repeated to improve the overall quality of the whole population. In this research, an effective SA was adopted to form a new LS method by parallel structure, which adds the elite library and mutation operation of the GA. According to mutation probability P m and the elite library, the feasible solutions are divided into local optimal solutions and general solutions. When a feasible solution is searched locally, it is judged whether it is in the elite library or not. If the solution is in the elite library or P m . random(0,1), the solution is considered to be a non-potential solution, which is executed by the mutation operation. Otherwise, the solution is considered to be a potential solution, which is executed by the SA to avoid falling into the local optimum and searching for a better solution. The LS process is shown in Figure 2. The symbols in the Figure 2 are defined as follows: mutation probability P m , initial temperature T 0 , current temperature T, termination temperature T f , temperature coefficient T c , cooling factor a, maximum number of perturbation L k . The probability exp(-My/ (T c T)) is the acceptance probability of the worse solution. When My is a small positive number and current temperature is high, the acceptance probability is close to 1, when current temperature is closer to the Algorithm 2. The pseudo-code of the multi-job crossover.
Input: 2 parent individuals S 1 and S 2 Output: 2 offspring S 1 ' and S 2 ' 1 j = the number of jobs in S 1 ; 2 jc = a set of randomly selected jobs; 3 S c1 = the operation sequencing of jc in S 1 ; 4 S c2 = the operation sequencing of jc in S 2 ; 5 S r1 = the remaining gene of S 1 2 S c1 in the original position of S 1 , and the rest is 0; 6 S r2 = the remaining gene of S 2 -S c2 in the original position of S 2 , and the rest is 0; 7 S 1 ' = replace 0 in S r1 with the value in S c2 in order; 8 S 2 ' = replace 0 in S r2 with the value in S c1 in order; termination temperature, the acceptance probability is closer to 0.
Elite library. The elite library is introduced as a distribution strategy in the LS to avoid searching for the neighborhood solutions of the local optimal individual, which may lead to the output individual being the original individual itself or worse than the original individual. The elite library is a limited set and empty at first. Then the individual output by the SA is regarded as a local optimal individual and stored in the elite library as chromosome. When it is found that the individual already exists in the elite library, the mutation operation will be implemented instead of using the SA. The storage space of elite library (the number of local optimal individuals) is adjustable; if the space is too large, it will increase the time to check whether the individual exists in the library, on the contrary, if the space is too small, it will not play the role of avoiding repeated search.
Mutation operator. Common mutation operations include the interchange mutation, reverse sequence mutation, insertion mutation, and neighborhood-based search mutation. We use the random mutation for the machine selection of the chromosomes.
Step 1: Select several locations for the machine part of the mutated chromosome.
Step 2: Choose 1 machine to randomly replace the original machine.
Perturbation operator. The perturbation operator searches the neighborhood of the current individual. Since the maximum number of the perturbation is set to L k , the L k neighborhood solutions of solution S can be obtained and the optimal individuals can be output. The move strategy of operation proposed by Mastrolilli and Gambardella 29 is used to generate the neighborhood solution of the current individual.
Due to the additional transportation time constraint in this paper, some equations and definitions need to be modified. For each operation O i,h , the release time r Oi,h is shown in equation (12), the tail time t Oi,h is shown in equation (13). An operation O i,h is critical if and only if r Oi,h + T i,h,k + t Oi,h = C max .
For a critical operation O i1,h1 , it needs to be inserted to a suit place to minimize the makespan. Let Q k be the sequence containing all the operations processed on machine k (except v if it processed on machine k). Then for the two subsequences of Q k are shown in equations (14) and (15). It has been proved that the set obtained inserting O i1,h1 after all the operations of P k -R k and before all the operations of R k -P k contains at least an optimal insertion of O i1,h1 on k. 29 The perturbation operation jumps out the optimal individual from all neighborhood solutions of the current individual in each loop iteration, and then updates the optimal individual to the current individual.

Termination condition
The termination condition of the IMA is that Iter equals the maximum number of iterations or that the generations of optimal individual unchanged equals n generations.

Experimental results and comparisons
All algorithms were performed in MATLAB R2017a on a PC with an Inter(R) Core(TM) i7-4785 T CPU at 2.2 GHz with 8 GB of RAM. The experiment tested the performance of the IMA for the TT-FJSP. The main parameters of the IMA, these parameter values are determined by extensive preliminary experiments, are shown in Table 3.

Experiment 1
To evaluate the effectiveness of IMA, we selected an instance (8 3 5) from Zhang et al., 12 who proposed the IGA to solve the TT-FJSP. We compared the results of the IGA and IMA after we had run the 2 algorithms 10 times. The Gantt charts of the optimal individuals from the IGA are shown in Figures 3 and 4. The Gantt charts of the optimal individuals from the IMA are shown in Figures 5 and 6. The convergence curve of the solutions of the TT-FJSP from the IGA is shown in Figure 7.
The convergence curve of the solutions of the TT-FJSP from the IMA is shown in Figure 8. As shown in Figures 4 and 6, the makespan of the optimal individual from the IMA is 5 units faster than the makespan of the optimal individual from the IGA. In the Gantt chart of the optimal individual from the IMA in Figure 6, the operations are evenly distributed, and there are no large idle time periods on each    machine. In contrast, in the Gantt chart of the optimal individual from the IGA in Figure 4, there are 2 idle time periods on M 2 . By comparing Figures 5 and 7, we found that the makespan of the optimal individual from the IMA is 4 units faster than the makespan of the optimal individual from the IGA, and there is 3 idle time period on M 4 in Figure 5. By comparing Figures 8 and  9, we find that the value of the optimal individual from the IMA increases with the mean value of the population. The mean value of the population from the IGA does not increase, and the value of the optimal individual easily falls into the local optimum. Therefore, the IMA is more effective than the IGA.

Experiment 2
Because of the lack of experimental data on the TT-FJSP, we chose the IGA as the object of comparison. The experimental data sets were Brandimarte's 21 data sets, and the randomly generated transportation time are shown in Table 4. To obtain meaningful results, we ran our algorithm 10 times on the same instance. The results are shown in Table 5, the opt is the makespan of the optimal individual, the m is the average value, the s is the standard deviation, the Cv is the coefficient of variation. The coefficient of variation is defined as the ratio of the standard deviation to the average value and eliminates the influence of dimension, as shown in equation (16). It can express the stability of algorithm more accurately. MAPE is the mean absolute percentage error; it shows that the percentage of the IMA's optimal individual is higher than that of the IGA's optimal individual, as shown in equation (17). Table 5 shows that the proposed IMA has a better performance than the IGA when solving the 10 problems of Mk01-Mk10 with the transportation time, including all the optimal and average values. The optimal individual from the IMA is 26.67% better than that from the IGA. The stability of the IMA and IGA is compared using the coefficient of variation shown in Figure 9. The average value of Cv IMA is 2.21% and the average value of Cv IGA is 2.66%. Therefore, the robustness of the IMA is better than that of the IGA.
Let us take the Mk01 problem as an example. The Gantt chart of the optimal individual from the IMA is shown in Figure 10. The Gantt chart of the optimal individual from the IMA is shown in Figure 11. As shown in Figures 10 and 11, the makespan of the optimal individual from the IMA is 6.9 units faster than the makespan of the optimal individual from the IGA. However, in both the Gantt charts, there are relatively

Conclusion and future research
In this paper, we studied the TT-FJSP and proposed an IMA. In this algorithm, feasible solutions are treated differently by the local search, general solutions are executed by the SA to improve the quality of the solutions, and the local optimal solutions are executed by the mutation operation to increase the diversity of the solutions. And the high-quality solutions are not lost as the population evolves because of the elite library. The IMA is evaluated and compared with an improved GA. The experimental results show that the optimal individual of each problem obtained from the IMA is better than that obtained from the IGA, and the convergence curve of the solutions is more reasonable with less uncertainty. Therefore, the IMA is effective for solving the TT-FJSP.
Our future work will address two main areas. First, we will further investigate the proposed algorithm. The SA has good search accuracy, but its search speed is not very fast, and it can be further improved. Therefore, the combination of multiple global search algorithms and local search algorithms is also worth trying. In addition, we will continue to explore the FJSP. The transportation constraints are not only the time constraints, we will also consider the limitations of routes and transportation resources.