A Variable Block Insertion Heuristic for Solving Permutation Flow Shop Scheduling Problem with Makespan Criterion

In this paper, we propose a variable block insertion heuristic (VBIH) algorithm to solve the permutation flow shop scheduling problem (PFSP). The VBIH algorithm removes a block of jobs from the current solution. It applies an insertion local search to the partial solution. Then, it inserts the block into all possible positions in the partial solution sequentially. It chooses the best one amongst those solutions from block insertion moves. Finally, again an insertion local search is applied to the complete solution. If the new solution obtained is better than the current solution, it replaces the current solution with the new one. As long as it improves, it retains the same block size. Otherwise, the block size is incremented by one and a simulated annealing-based acceptance criterion is employed to accept the new solution in order to escape from local minima. This process is repeated until the block size reaches its maximum size. To verify the computational results, mixed integer programming (MIP) and constraint programming (CP) models are developed and solved using very recent small VRF benchmark suite. Optimal solutions are found for 108 out of 240 instances. Extensive computational results on the VRF large benchmark suite show that the proposed algorithm outperforms two variants of the iterated greedy algorithm. 236 out of 240 instances of large VRF benchmark suite are further improved for the first time in this paper. Ultimately, we run Taillard’s benchmark suite and compare the algorithms. In addition to the above, three instances of Taillard’s benchmark suite are also further improved for the first time in this paper since 1993.


Introduction
Sustainability in manufacturing industries is mainly measured by their competitiveness in the market place. Competitiveness is referred to timely product delivery with the best quality, minimum manufacturing time and price to customers. Minimum manufacturing time can be obtained by optimal production sequences that can minimize makespan or total flowtime. Note that a manufacturing company can fail to satisfy production plans although the other production entities such as operators, maintenance, inventory, quality control, etc. are in control due to the lack of optimal or near optimal production sequences in the shop floor. For this reason, seeking optimal or near-optimal production sequences and schedules is vital to manufacturing companies in order to minimize the makespan, which also minimizes idle times on the machines and maximize machine utilization.
The permutation flow shop scheduling problem (PFSP) has been widely studied in the literature and has extensively been applied in the inustry. There are many different fields in real-life where PFSP can be used [1]. It is yet an exceptionally active topic of investigation, especially because flow shop environments are at the center of real-life scheduling problems in various fields of high social or economic impact. In addition, the flow shop layout is a regular configuration in many manufacturing companies. The basic PFSP consists of a set of n jobs which are processed by m machines. These jobs follow the same route and their operations on the machines cannot be interrupted. All the jobs must be processed in the same order on the machines and the aim is to find the best permutation π = {π 1 , π 2 , . . . , π n } of these jobs with respect to the given objective.
In this study, our aim is to maximize the throughput of the system by maximizing the utilization rate of the machines which means minimizing makespan. To compute the makespan, π denotes the given arbitrary solution, where job π i is the job at the ith position of solution π. C i,k is denoted as the completion time of job π i on machine k at position i. Following this notation, completion times of jobs at each machine are computed as in the following Equations (1)- (5), where p π i ,k is the processing time of job π i at the kth machine. The makespan of solution π, denoted as C max (π), is the completion time of the last job (i.e., n) on the last machine (i.e., m). It is simply denoted as C nm and calculated as follows: (1) C i,1 = C i−1,1 + p π i ,1 ∀i = 2, . . . , n C 1,k = C 1,k−1 + p π 1 ,k ∀k = 2, . . . , m C i,k = max C i−1,k , C i,k−1 + p π i ,k ∀i = 2, . . . , n; ∀k = 2, . . . , m C max (π) = C nm .
The PFSP with makespan criterion is denoted as Fm|Permu|Cmax according to the notation of [2] and has been proven to be NP-hard for the makespan criterion [3], so it is challenging to solve it with exact methods. Therefore, metaheuristic algorithms were employed to solve the problem and obtain near-optimal solutions. In recent years, various metaheuristic algorithms have been presented to solve various variants of PFSP with different objectives. One of the state-of-the-art algorithms for PFSPs is the iterated greedy algorithm (IG) presented by [4]. We focused on the recent literature that considers the IG algorithm in their solution approaches.
The IG algorithm was employed to PFSP with makespan criterion in [4][5][6][7][8][9]. In [5], to improve the solution quality, a local search was applied to the partial solution after the destruction step of the IG algorithm, while in [6] sequence depended setup times were employed for the PFSP with makespan criterion. In addition, in [7] the authors studied the PFSP with makespan and proposed a tie-breaking mechanism for the IG algorithm, while in [8] an IG and a discrete differential evolution algorithm were proposed and compared. In this study, we employ new hard VRF instances which are first introduced in [9], and they also applied an IG algorithm. In addition, the same problem was studied in [10] to minimize the makespan over Taillard's benchmark suite.
The IG algorithm was applied to various variants of PFSP such as no-wait flow shops in [11][12][13]; blocking flow shops in [14][15][16][17]; no-idle flow shops in [18][19][20]; energy-efficient PFSP in [21,22]; multi-objective PFSP in [23,24] where both studies presented a restarted iterated Pareto greedy algorithm. In a no-wait variant of PFSP, distributed no-wait flow shop problem [11], tabu-based reconstruction strategy [12], and sequence depended setup times [13] were employed with IG algorithm. In blocking variant of PFSP, IG algorithms were combined with local search algorithms [14], constructive heuristics [15,16], and also embedded in differential evolution framework [17]. In [25] profile fitting and NEH heuristic algorithms were proposed for the same problem. In a no-idle variant of PFSP, iterated reference greedy algorithm [18], and variable IG algorithm [19] were presented. In addition, IG algorithm was employed for the mixed no-idle PFSP [20]. IG algorithm was also applied to PFSP with different objective functions such as total tardiness criterion [26,27]; total flowtime criterion [28]. In [1], they carried out an exhaustive review and computational evaluation of heuristics and metaheuristics published until 2017 for the PFSP to minimize the makespan. Therefore, for the further analysis of the literature of PFSP, the indicated manuscript [1] should be examined.
In traditional search algorithms, swap and insertion neighborhood structures are generally employed. The swap neighborhood exchanges two jobs in a solution, whereas the insertion one removes a job from a solution and inserts it into another position in the solution. Recently, block move-based local search algorithms are presented for the single machine-scheduling problems in the literature [29][30][31][32]. Xu et al. [31] developed a Block Move neighborhood structure in which l consecutive jobs (called a block) are inserted into another position in the solution. They represent a block move by a triplet (i, k, l), where i denotes the position of the first job of the block, k the target position of the block to be inserted and l the size of the block. It is obvious that one edge insertion, two edge-insertion and 3-block insertion are the block move neighborhoods with l = 1, l = 2, and l = 3. Similarly, Gonzales and Vela [32] developed a variable neighborhood descent algorithm with three distinct block move neighborhoods and employed in a memetic algorithm. Then, a memetic algorithm with block insertion heuristic is presented in [29]. Moreover, in [33], a variable block insertion heuristic (VBIH) algorithm was employed to solve the blocking PFSP with makespan criterion.
In IG algorithms, some solutions components are removed from the current solution and reinserted into the partial solutions. In other words, a number, dS, of jobs are removed randomly, which is known as the destruction phase. Then, in the construction phase, these dS jobs are reinserted into the partial solution in the same order they are removed. For each of dS jobs, it makes a number n − dS + 1 of insertions. However, the VBIH algorithm removes a block of jobs π b with size b from the current solution and it makes a number n − b + 1 of block insertions only. That is the difference between IG and VBIH algorithms.
The main contributions of the paper can be outlined as follows. VBIH is employed to solve PFSP with makespan criterion using the new hard VRF benchmark sets [9]. Detailed computational results show that VBIH algorithm outperforms two variants of the iterated greedy algorithm. 236 out of 240 instances of large VRF benchmark suite are further improved for the first time in this paper, while the results of the remaining four instances are found as the same with the current results. In addition, the formulation of two mathematical models is given to solve the small benchmark set in order to verify the results of our proposed VBIH algorithm. One hundred and eight out of 240 small instances are proven to be optimal. Therefore, this paper proposes new lower bounds with the use of an efficient algorithm, which differentiates the study from the current literature. We also show that the speed up method of Taillard is substantially effective when solving the PFSP with makespan criterion.
The remaining part of the paper is organized as follows: Section 2 introduces the formulation of PFSP including mixed integer programming (MIP) model and constraint programming (CP) model whereas Section 3 presents all the heuristic algorithms. Section 4 explains the computational results of the MIP and CP models on small VRF instances to show the solution quality of the heuristic algorithms and the limitations of the models. Section 5 reports the experimental results of the heuristic algorithms and the improvements to the large VRF instances. Finally, Section 6 summarizes the concluding remarks.

Mathematical Model Formulation
This paper proposes MIP and CP models to solve small VRF instances for PFSP with the makespan criterion in order to verify the solution quality of proposed heuristic algorithms. The input parameters used in the models are presented in follows:

Parameters:
n: Total number of jobs, i = 1, . . . , n m: Total number of machines, k = 1, . . . , m p i,k : Processing time of job i on machine k M : A sufficient large constant integer.

The MIP Model
The MIP decision variables, objective function and the constraints are given in the following equations. The MIP formulation of PFSP, which were proposed by Manne [34], is used.

Decision Variables:
C max : Makespan C i,k : Completion time of job i on machine k D i,j : Binary variable: 1 if job i is scheduled before job j; 0, otherwise; i < j
The objective function (6) minimizes the makespan while Constraint (7) calculates the maximum completion time of all jobs on the last machine. In PFSP, all jobs follow the same route through the machines so that their final processes will be done on the last machine. Constraint (8) computes the completion time of each job on machine 1 ensuring that they cannot occur earlier than the duration of their processing time on machine 1 which is the starting machine for all jobs. Constraint (9) ensures that the completion time of each job on each machine cannot be processed before their completion time on the previous machine. Constraints (10) and (11) specify the relationship between the processing of two consecutive jobs on the same machine. Constraint (11) starts that if job i precedes job j in the permutation, then job i should be completed before job j on each machine. Otherwise, job j should precede job i on each machine which is shown by Constraint (10).

The CP Model
CP decision variables, objective function and the constraints are presented in the following equations using the OPL API of CP Optimizer. To express the processing times of the jobs on the machines, the model uses interval variables denoted as JobInt. In addition, sequence variables for the machines are defined in the model as Machine which collects all these interval variables.

Decision Variables:
JobInt i,k : Interval variable for job i on machine k with duration p i,k Machine k : Sequence variable for machine k over JobInt i,k 1 ≤ i ≤ n .
The CP model minimizes the makespan by computing the maximum end date of each job on the last machine (13). Constraint (14) impose the precedence constraints between the consecutive operations of each job on the sequence of machines. Machines are the disjunctive resources and can process only one job at a time, which is expressed by the noOverlap Constraint (15) over the sequence variables associated with machines. The relationship between sequence variables and the interval variables are provided while defining the decision variables. The last constraint sameSequence (16) guarantees that all the jobs are processed in the same order on each machine. Therefore, the permutation of the jobs will be the same for each machine.

Taillard's Speed Up Method for PFSP with Makespan Criterion
Insertion neighborhood structure is very effective for makespan minimization. The size of the insertion neighborhood is (n − 1) 2 . Since each objective function evaluation takes O(nm) time, its computational complexity is O(n 3 m). In [35], a speed-up method is proposed where it reduces the computational complexity from O(n 3 m) to O(n 2 m) for the PFSP with makespan criterion. In order to execute the insertion procedure in time O(nm), this speed-up method can be explained as follows: Suppose that job π i will be inserted in a position l. Then the speed up method can be described below:

3.
Compute the earliest relative completion time f i,k on the lth machine of job π j inserted at the lth position. Completion time of an inserted job on the first machine is zero.

4.
The value of the makespan C max,l when inserting job j at the lth position is: In order to illustrate the speed up the procedure, we give the 7-job 2-machine example. Note that Johnson's algorithm [36] solves this problem to optimality. Hence, in Table 1, we provide the problem instance with the processing times as well as the optimal solution. According to the Johnson's algorithm [36], the optimal solution is {1, 2, 7, 3, 5, 4, 6} with C max = 36. Now, suppose that we remove job 7 and obtain the partial solution, {1, 2, 3, 5, 4, 6}. Suppose that we insert job 7 into position l = 3 of the partial solution to obtain the optimal solution. We follow the speed up method now:
Speed-up calculation of the partial solution is given in Figure 1. Speed-up calculation of the partial solution is given in Figure 1. Speed-up calculation of the complete solution is given in Figure 2.
The value of the makespan C max,l when inserting job π i at the lth position is:  It is clear that the above speed-up method reduces the complexity of the whole insertion neighborhood from ( ) to ( ) . This speed-up method is the key to success for any algorithm for PFSP with makespan criterion. For this reason, we have chosen the Car8 instance from the literature in order to illustrate the speed-up method above in detail. From the literature, we know It is clear that the above speed-up method reduces the complexity of the whole insertion neighborhood from O(n 3 m) to O(n 2 m). This speed-up method is the key to success for any algorithm for PFSP with makespan criterion. For this reason, we have chosen the Car8 instance from the literature in order to illustrate the speed-up method above in detail. From the literature, we know that best or optimal solution is {7, 3, 8, 5, 2, 1, 6, 4} with C max = 8366. In Appendix A, we remove job 2 from the optimal solution and re-insert it into the 5th position again. A detailed implementation of Taillard's speed up method is given in Appendix A in order to ease the understanding of it.

IG Algorithms
IG algorithms mainly have four components; namely, initial solution, destruction-construction (DC) procedure, local search, and acceptance criterion. The traditional IG RS is proposed by [4]. In this algorithm, the initial solution is constructed by the NEH heuristic in [37]. In the destruction step, dS jobs are randomly removed from the solution π without repetition and stored in π D . The remaining jobs are also stored in π P that represents the partial solution. In the construction step, each job in π D is inserted into the partial solution π P , in the order in which they were removed, until a complete solution of n jobs is constructed. Having carried out the destruction and construction procedure, a local search is employed to further enhance solution quality. After a local search, if the solution is better than or equal to the incumbent solution, it is accepted. Otherwise, it is accepted with a simple simulated annealing-type acceptance criterion, which is suggested by [38]: where τP is a parameter to be adjusted. The pseudo-code of the traditional IG RS is given in Algorithm 1, where r is a uniform random number between 0 and 1.
The IG RS algorithm for the PFSP under makespan minimization employs an initial solution generated by the NEH heuristic. In addition, the NEH heuristic was extended to the FRB5 heuristic with a local search on the partial solutions [39]. Both heuristics are simple and very effective for minimizing the makespan, and its pseudo-code is given in Algorithm 2. In the first phase, the sum of the processing times on all machines are calculated for each job. Then, jobs are sorted in decreasing order to obtain δ. In the second phase, the first job in δ is selected to establish a partial solution π 1 . The remaining jobs in δ are inserted in the partial solution one by one. After each iteration, optionally, a local search is applied to the partial solution. Local search is implemented as long as the partial solution is improved. After having inserted all jobs, a complete solution is obtained. Note that the NEH heuristic is denoted as FRB5 heuristic with an optional local search to partial solutions.
//Algorithm 3 f or FRB5 heuristic end f or return π with n jobs and f (π) The IG RS algorithm employs insertion neighborhood structure as a local search after destruction and construction procedure. Insertion neighborhood is very effective with the speed-up method explained in Section 3.1 for makespan minimization. Insertion neighborhood can be deterministic or stochastic depending on the decision of choosing a job from solution to be removed. The deterministic variant is given in Algorithm 3. This procedure removes π i from the solution π and inserts it into all possible positions of the incumbent solution π. When the best-improving insertion position is found, job π i is inserted into that position. These steps are repeated for all jobs. If an improvement is observed, the local search is re-run until no better solution is obtained.
In the stochastic variant given in Algorithm 4, jobs are randomly chosen from solutions to make insertions. In Algorithm 4, job π k at position k is randomly chosen from the solution π without repetition, and partial solution π P is obtained. Then, job π k is inserted into all possible positions of the partial solution π P . When the best-improving insertion position is found, job π k is inserted into that position, and a complete solution π * is obtained. These steps are repeated for all jobs. If an improvement is found, the local search is rerun again until no better solution is obtained.

Algorithm 4: First improvement insertion neighborhood(π)
f or i = 1 to n do π P = Remove job π k f rom solution π randomly and without repetition Recently, a new IG ALL algorithm has been presented in the literature [5] with excellent results for the PFSP with makespan minimization. The difference between IG ALL and IG RS is that IG ALL applies an additional local search to partial solutions after destruction, which substantially enhances solution quality. In the IG RS algorithm, local search is applied to the complete solution after the construction phase to improve the current candidate solution whereas, in IG ALL algorithm, local search is applied to a partial solution after destruction phase. This idea is applied in heuristic algorithms by Reference [39].
They study on vehicle routing problem and apply local search on the routes in the construction phase. Applying local search to the partial solution is more advantageous in terms of computational time and providing different search directions. Due to having a partial solution, a local search is applied to the smaller size of the complete solution so that the search procedure can be conducted quickly. Another difference between IG RS and IG ALL is due to the fact that the initial solution is constructed by FRB5 heuristic. The pseudo code of IG ALL algorithm is presented in Algorithm 5.
return π best and f (π best ) endprocedure Note that Algorithm 3 is used in the FRB5 heuristic in order to construct the initial solution with a single run due to its deterministic property. In both algorithms, Algorithm 4 is employed in applying local search to both partial and complete solutions.

Variable Block Insertion Algorithm
In this paper, we propose a VBIH algorithm as follows. The VBIH algorithm employs the FRB5 heuristic as an initial solution. It has a minimum block size (b min ), and a maximum block size (b max ). It removes a block of jobs (π b ) with size b from the current solution and obtains a partial solution (π P ). Similar to the IG ALL algorithm, it applies the local search in Algorithm 4 to the partial solution. Then, it makes a number, n − b + 1, of block insertion moves sequentially in the partial solution. It chooses the best one amongst those solutions from block insertion moves. Well-known RIS local search in the literature is applied to the complete solution found after block insertion moves. If the new solution obtained after the local search is better than or equal to the current solution, it replaces the current solution. As long as it improves, it retains the same block size (i.e., b = b). Otherwise, the block size is incremented by one (i.e., b = b + 1) and a simulated annealing-based acceptance criterion, similar to IG RS and IG ALL algorithms, is employed to accept the new solution to escape from local minima. This process is repeated until the block size reaches its maximum limit (i.e., b ≤ b max ). The outline of the VBIH algorithm is given in Algorithm 6. Note that π R is the reference sequence; tP is temperature parameter for the acceptance criterion and r is a uniform random number between 0 and 1.
Regarding the local search algorithm that will be applied only to complete solutions, we use a well-known referenced insertion scheme local search, RIS [8,40]. RIS is guided by a reference solution π R , which is the best solution obtained so far during the search process. For instance, if the reference solution is given by π R = {3, 5, 1, 4, 2} and the current solution by π = {1, 2, 3, 4, 5}. The reference solution π R implies that job 3 in the current solution π might not be in a proper position. For this reason, the RIS local search first removes job 3 from the current solution π and inserts it into all possible slots of the partial solution π P . A new solution with the best insertion slot is replaced by the current solution, and the iteration counter is reset to one if any improvement occurs. Otherwise, the iteration counter is incremented by one. Then, it removes job 5 from the current solution π and inserts it into all possible positions of the partial solution π P . This procedure is repeated until the iteration counter is greater than the number of jobs n, and a new solution is obtained. The pseudo-code of the RIS local search is given in Algorithm 7.
After the local search phase, it should be decided if the new solution is accepted as the incumbent solution for the next iteration. A simple simulated annealing-type of acceptance criterion is used with a constant temperature similar to the IG RS and IG ALL algorithms. Note that Taillard's speed-ups are employed wherever possible in our code.

Design of Experiment for Parameter Tuning
In this section, we present a Design of Experiments (DOE) approach [41] for parameter settings of the VBIH algorithm. In order to carry out experiments, we generate random instances with the method proposed in [9]. In other words, random instances are generated for each combination of n ∈ {100, 200, 300, 400, 500, 600, 700, 800} and m ∈ {20, 40, 60}. Five instances are generated for each job and machine combination. Ultimately, we obtained 1200 instances in total. We consider three parameters in the DOE approach. These are maximum block size (bMax), temperature adjustment parameter (τP), and the decision of whether or not to implement the local search to the partial solution after removal of a block of jobs. We have taken the maximum block size with seven levels as bMax ∈ (2, 3, 4, 5, 6, 7, 8); the temperature adjustment parameter with ten levels as τP ∈ (0.1, 0.2, 0.3, 0.4, 0.5); and the decision on the local search to partial solutions as pL ∈ (1, 2). Note that pL = 1 means that the local search is applied to partial solutions whereas pL = 2 does not apply the local search to partial solutions. In the design of VBIH algorithm, there are 7 × 5 × 2 = 70 algorithm configurations, i.e., treatments. The VBIH algorithm is coded in C++ programming language on Microsoft Visual Studio 2013, and a full factorial design of experiments is carried out for each algorithm on a Core i5, 3.40 GHz, 8 GB RAM computer. Each instance is run for 70 treatments with a maximum CPU time T max = 10 × n × m milliseconds. Note that it took 18 days to run the full factorial design. We calculate the relative percent deviation (RPD) for each instance-treatment pair as follows: where CMAX i is the makespan value generated by the VBIH algorithm in each treatment and CMAX min is the minimum makespan value found amongst 70 treatments. For each job size-treatment pair, the average RPD value is calculated by taking the average of five instances in each job size. Then, the response variable (ARPD) of each treatment is obtained by averaging these RPD values of all job sizes. After determining the ARPD values for each treatment as mentioned above, the main effects plots of the parameters are analyzed and given in Figure 3. As it can be seen from Figure 3, the following parameters have better ARPD values than the others: = 2, = 0.5, and = 1. Furthermore, in order to see whether or not there is an interaction effect between parameters, an ANOVA analysis is also given in Table 2.  Table 2 indicates that , , and were statistically significant since higher magnitude of values and -values of parameter interaction effects are less than the significance level = 0.05. High magnitude of value for also suggest that applying local search to partial solutions has a significant impact on the solution quality as mentioned in [5]. In terms of interaction effects, it can be observed that × interaction is not significant because the p-value is much higher than the significance level = 0.05 . However, × and × interactions were significant since their p values are less than the significance level = 0.05. The interaction effects plot for × is given in Figure 4. As it can be seen from Figure 3, the following parameters have better ARPD values than the others: bMax = 2, τP = 0.5, and pL = 1. Furthermore, in order to see whether or not there is an interaction effect between parameters, an ANOVA analysis is also given in Table 2.  Table 2 indicates that bMax, tP, and pL were statistically significant since higher magnitude of F values and p-values of parameter interaction effects are less than the significance level α = 0.05. High magnitude of F value for pL also suggest that applying local search to partial solutions has a significant impact on the solution quality as mentioned in [5]. In terms of interaction effects, it can be observed that bMax × tP interaction is not significant because the p-value is much higher than the significance level α = 0.05. However, bMax × pL and tP × pL interactions were significant since their p values are less than the significance level α = 0.05. The interaction effects plot for bMax × pL is given in Figure 4.   Figure 4 indicates that maximum block size should be taken as = 2 and local search to the partial solution should be applied. Since × interaction is also significant, we provide the interaction plot in Figure 5.   Ultimately, we set the parameters of VBIH algorithm as follows: = 2, = 0.5, and = 1.

Computational Results
In this section, the computational results for the small and large set of VRF benchmark sets are provided. MIP and CP models were written in OPL and run on the IBM ILOG CPLEX 12.8 software  Figure 4 indicates that maximum block size should be taken as bMax = 2 and local search to the partial solution should be applied. Since tP × pL interaction is also significant, we provide the interaction plot in Figure 5.   Figure 4 indicates that maximum block size should be taken as = 2 and local search to the partial solution should be applied. Since × interaction is also significant, we provide the interaction plot in Figure 5.   Ultimately, we set the parameters of VBIH algorithm as follows: = 2, = 0.5, and = 1.

Computational Results
In this section, the computational results for the small and large set of VRF benchmark sets are provided. MIP and CP models were written in OPL and run on the IBM ILOG CPLEX 12.8 software  Ultimately, we set the parameters of VBIH algorithm as follows: bMax = 2, τP = 0.5, and pL = 1.

Computational Results
In this section, the computational results for the small and large set of VRF benchmark sets are provided. MIP and CP models were written in OPL and run on the IBM ILOG CPLEX 12.8 software suite, while all the heuristic algorithms were being written in Visual C++ 13 and carried out on an Intel Core i5, 3.40 GHz, 8 GB RAM computer. The proposed VBIH algorithm is compared to IG RS and IG ALL algorithms. In addition, the results of these algorithms are obtained without the Taillard's speed up method, and they are denoted as IG RS *, IG ALL * and VBIH*. Regarding parameters of them with, destruction size ds, and temperature adjustment factor, tP are taken as ds = 4 and tP = 0.4 for IG RS and IG RS * as suggested in [4]. They are taken as ds = 2 and tP = 0.7 for IG ALL and IG ALL * as indicated in [5]. As explained in the previous section DOE is conducted for the VBIH algorithm and its parameters are determined as follows: bMax = 2, τP = 0.5, and pL = 1, which are also used for the VBIH* algorithm.

MIP Versus CP
Computational results are given in Table 3 for each combination, giving a total of 240 small VRF instances. For each combination, the table summarizes the number of optimal solutions (nOpt) found for ten instances of each job-machine combination (n × m), the average relative percent deviation (ARPD%) from the upper bounds given in [9], the average CPU time for its ten instances, and the optimality gap percentage (GAP%) on termination, which means the gap between best lower and best upper bound. The maximum CPU time is restricted to an hour (3600 s). The result of CP and MIP models are compared for job sizes 10 and 20. While MIP model can find solutions for very small sized instances (10 jobs) in a shorter time than CP model, it becomes hard for MIP to solve large sized problems (20 jobs and more). Both models cannot always find optimal solutions when the machine size becomes greater than 5, but the MIP model has larger gaps than the CP model. The results show that CP is more efficient than MIP on PFSP, except for very small-sized instances. The results of the remaining instances are obtained only from the CP model because of very large gaps by MIP model. CP model always captures optimal solutions when the machine number is five regardless of the number of jobs. Besides, CP can find optimal solutions in some of the instances when the machine size is 10. Overall, within the time limit, the CP model verifies optimality for 108 out of 240 instances.

Comparison of Heuristic Algorithms with Exact Solutions
In order to compare performances of heuristic algorithms with CP exact method, we run all algorithms for five independent replications with different seed numbers. Relative percent deviation values from upper bounds for ten different instances of each job-machine combinations are calculated as follows: where M is the makespan value generated by any heuristic; and M UB is the upper bound provided in [9]. Note that, for each instance, we record the average RFD of five replications for statistical analysis purposes, especially, for interval graphs. The solutions of the CP model are limited to 3600 s and its average CPU times are given in Table 4. IG ALL , IG RS , and VBIH algorithms are run for five replications with three different time limits 15, 30, and 45 × n × m. As expected, the performance of these algorithms is much better than those by CP exact model, and they improve the upper bounds provided in [9], which means that the proposed algorithm and other IG algorithms can find good (optimal in some cases) solutions in a very short time. As the solution time increases, the solution quality of VBIH algorithm increases and according to the RPD, it gives the best solutions amongst all other algorithms. It should be noted that the VBIH algorithm further improves 64 out 240 upper bounds for small VRF instances within a very short time.

Large VRF Instances
Note that both IG ALL and VBIH algorithms employ the FRB5 heuristic for constructing initial solution whereas IG RS uses the traditional NEH heuristic. For the large VRF instances, Table 5 summarizes the ARPD values of heuristic methods such as NEH, NEH without speed-up, denoted as NEH*, and extended NEH heuristic with a local search on partial solutions denoted as FRB5. As shown in Table 5, NEH is very fast with 0.04 s on overall average CPU time. However, its overall average of ARPD is 3.33%. Although FRB5 heuristic is computationally very expensive, which is 43.76 s on overall average CPU time, its average ARPD is only 0.89% from the upper bounds. It is obvious from Table 5 that FRB5 heuristic is substantially better than NEH with a very large margin at the expense of increased CPU time. It is interesting to observe the CPU time performance of the NEH heuristic without the speed-up method of Taillard. Table 5 clearly indicates that the Taillard's speed-up method is substantially effective since the overall average CPU time is jumped from 0.04 s to 4.95 s without the speed-up method of Taillard. In addition to the above, we present the interval graph of both constructive heuristics in Figure 6 in order for justification. Figure 6 indicates that differences in ARPDs are significantly meaningful on the behalf of FRB5 heuristic since their confidence intervals do not coincide.

Computational Results of Metaheuristics
In this section, the performance of VBIH algorithm is compared to the best-performing algorithms, IG RS and IG ALL , from the literature. All algorithms are run five replications to solve the large VRF instances. In Table 6, we present average, minimum and maximum ARPD values for the CPU time limit T max = 15 × n × m milliseconds.

Computational Results of Metaheuristics
In this section, the performance of VBIH algorithm is compared to the best-performing algorithms, IGRS and IGALL, from the literature. All algorithms are run five replications to solve the large VRF instances. In Table 6, we present average, minimum and maximum ARPD values for the CPU time limit = 15 × × milliseconds.    As seen in Table 6, VBIH generated better Avg, Min and Max RPD values on the overall average. On overall average, it was able to further improve the upper bounds up to −0.05%; its best overall performance was −0.18% indicating that 0.18% of 240 instances are further improved and its worst-case performance was 0.08%. In order to see if differences in ARPDs are statistically significant, we provide the 95% confidence interval plot of algorithms in Figure 7, where we can observe that differences in ARPD values are statistically significant on the behalf of VBIH against IG RS and IG ALL algorithms because their confidence intervals do not coincide. As seen in Table 6, VBIH generated better Avg, Min and Max RPD values on the overall average. On overall average, it was able to further improve the upper bounds up to −0.05%; its best overall performance was −0.18% indicating that 0.18% of 240 instances are further improved and its worstcase performance was 0.08%. In order to see if differences in ARPDs are statistically significant, we provide the 95% confidence interval plot of algorithms in Figure 7, where we can observe that differences in ARPD values are statistically significant on the behalf of VBIH against IGRS and IGALL algorithms because their confidence intervals do not coincide. Computational results for Avg, Min and Max ARPD values with the CPU time limit = 30 × × milliseconds are given in Table 7. As seen in Table 7, VBIH was able to generate better Avg, Min and Max ARPD values on the overall average. On overall average, it was able to further improve the upper bounds by −0.11% in Avg value, −0.24% of upper bounds are further improved on Min value and its worst-case performance was 0.02%. However, as CPU times increased, the performance of IGALL algorithm was also remarkable. Briefly, both VBIH and IGALL outperformed IGRS in almost each problem set.   Computational results for Avg, Min and Max ARPD values with the CPU time limit T max = 30 × n × m milliseconds are given in Table 7. As seen in Table 7, VBIH was able to generate better Avg, Min and Max ARPD values on the overall average. On overall average, it was able to further improve the upper bounds by −0.11% in Avg value, −0.24% of upper bounds are further improved on Min value and its worst-case performance was 0.02%. However, as CPU times increased, the performance of IG ALL algorithm was also remarkable. Briefly, both VBIH and IG ALL outperformed IG RS in almost each problem set.

IGRS I G ALL V B I H
In order to see if these results are statistically significant, we provide the 95% confidence interval plot of algorithms in Figure 8, where we can observe that differences in ARPD values are statistically significant on the behalf of both VBIH and IG ALL algorithms against IG RS algorithm because their confidence intervals do not coincide with IG RS . In other words, VBIH and IG ALL algorithms were statistically equivalent but significantly superior to IG RS .  In order to see if these results are statistically significant, we provide the 95% confidence interval plot of algorithms in Figure 8, where we can observe that differences in ARPD values are statistically significant on the behalf of both VBIH and IGALL algorithms against IGRS algorithm because their confidence intervals do not coincide with IGRS. In other words, VBIH and IGALL algorithms were statistically equivalent but significantly superior to IGRS. Computational results for average, minimum and maximum RPD values with the CPU time limit = 45 × × milliseconds are given in Table 8, where VBIH outperformed IGRS and IGALL algorithms with respect to average, minimum and maximum RPD values on the overall average. On overall average, it was able to further improve the upper bounds by −0.25% on the average value, -0.36% on the minimum value, and its worst-case performance was −0.13%. These statistics indicate that VBIH generated much better results than both the IGRS and IGALL algorithms. Computational results for average, minimum and maximum RPD values with the CPU time limit T max = 45 × n × m milliseconds are given in Table 8, where VBIH outperformed IG RS and IG ALL algorithms with respect to average, minimum and maximum RPD values on the overall average. On overall average, it was able to further improve the upper bounds by −0.25% on the average value, −0.36% on the minimum value, and its worst-case performance was −0.13%. These statistics indicate that VBIH generated much better results than both the IG RS and IG ALL algorithms. In order to see if these results are statistically significant, we provide the 95% confidence interval plot of algorithms in Figure 9, where we can observe that differences in ARPD values are statistically significant on the behalf of VBIH algorithm against both IG RS and IG ALL algorithms because their confidence intervals do not coincide. In other words, VBIH algorithm was statistically superior to both IG RS and IG ALL algorithm. In order to see if these results are statistically significant, we provide the 95% confidence interval plot of algorithms in Figure 9, where we can observe that differences in ARPD values are statistically significant on the behalf of VBIH algorithm against both IGRS and IGALL algorithms because their confidence intervals do not coincide. In other words, VBIH algorithm was statistically superior to both IGRS and IGALL algorithm.  In the Supplementary Materials, we summarize all the best-known solutions found for the first time by IG ALL and VBIH algorithms. The VBIH algorithm further improves 230 out of 240 instances. In addition, 173 out of 240 instances are improved by the IG RS algorithm, while the IG ALL algorithm further improves 222 out of 240 instances. The IG ALL algorithm improves six instances that are not improved by VBIH algorithm. Ultimately, 236 out of 240 instances are further improved by all algorithms within 45 × n × m time limits with the remaining four solutions being equal.
As mentioned before, IG ALL algorithm is presented in [5], where they analyzed the performances of IG RS and IG ALL on both Taillard's [42] and large VRF instances. They observed that the results obtained by using Taillard's benchmark set, both algorithms do not present very significant differences with respect to the RPDs obtained. In fact, they have shown that both algorithms did not show any statistically significant differences. However, statistically significant differences between IG RS and IG ALL have been shown when large VRF instances are employed. In order to validate this observation, we have run three algorithms on Taillard's benchmark set with a stopping criterion T max = 45 × n × m milliseconds. Furthermore, we run three algorithms without the Taillard's speed up method and they are denoted as IG RS *, IG ALL * and VBIH*. The computational results are given in Table 9. As seen in Table 9, VBIH produced much better RPDs than IG RS and IG ALL algorithms when the Taillard's speed up method is employed since its overall RPD was 0.17 from the best-known solutions. However, IG RS and IG ALL algorithms do not show so many differences in terms of RPDs. Interval plots of the algorithms in Figure 10 show that differences in RFDs are not statistically significant because their confidence intervals do coincide. This suggests a fact that researches on PFSP and its variants should employ VRF benchmark suite to see differences in algorithms newly presented. Figure 10 also shows that the Taillard's speed up method is significantly effective for all three algorithms. During these runs, we were also able to find 3 new best-known solutions for the Taillard's benchmark suite (ta054 = 3719, ta55 = 3610, ta56 = 3680) and their permutations are also provided in the Supplementary Materials. significant because their confidence intervals do coincide. This suggests a fact that researches on PFSP and its variants should employ VRF benchmark suite to see differences in algorithms newly presented. Figure 10 also shows that the Taillard's speed up method is significantly effective for all three algorithms. During these runs, we were also able to find 3 new best-known solutions for the Taillard's benchmark suite (ta054 = 3719, ta55 = 3610, ta56 = 3680) and their permutations are also provided in the supplementary materials.  Figure 10. Interval plot at the 95% confidence level for Taillard's instances.

Conclusions
This paper presents a variable block insertion heuristic (VBIH) algorithm for solving the permutation flow shop scheduling problem (PFSP) with makespan criterion. In addition, we introduce mixed integer programming (MIP) and constraint programming (CP) models to solve the small benchmark set and to verify the results of our proposed heuristic algorithm. By employing the time limited CP model, we can find optimal solutions for some of small VRF instances for the first time in the literature. Furthermore, all algorithms can generate better solution values than upper those currently exist in the literature. We adapted a well-known speed-up method of Taillard and applied all the necessary parts while coding the heuristic algorithms. The parameters of the proposed VBIH algorithm is tuned through a design of experiments on randomly generated benchmark instances. Extensive computational results on two new VRF benchmark suites show that the VBIH algorithm is superior to the best performing algorithms from the literature.