A Hybrid Local Search Algorithm for the Sequence Dependent Setup Times Flowshop Scheduling Problem with Makespan Criterion

This paper focuses on the flowshop scheduling problem with sequence dependent setup times (FSSP-SDST), which has been an investigated object for decades as one of the most popular scheduling problems in manufacturing systems. A novel hybrid local search algorithm called HLS is presented to solve the flowshop scheduling problem with sequence dependent setup times with the criterion of minimizing the makespan. Firstly, the population is initialized by the Nawaz-Enscore-Hoam based problem-specific method (NEHBPS) to generate high quality individuals of the current population. Then, a global search embedded with a light perturbation is designed to produce a new population. After that, to improve the quality of individuals in the current population, a single insertion-based local search is applied. Meanwhile, a further local search strategy based on the insertion-based local search is used to find better solutions for the individuals which are non-improved. Finally, the heavy perturbation is used to explore potential solutions in the neighbor region. To validate the performance of HLS, we compare our proposed algorithm with other competitive algorithms on Taillard benchmark problems. From the experimental results, it can be concluded that the proposed algorithm outperforms the benchmark algorithms.


Introduction
Permutation flowshop scheduling is one of the classical combinatorial optimization problems since the seminal work of Johnson Seminal [1], which is widely applied in real life, such as computer work, industrial engineering, mathematics. In particular, the permutation flowshop scheduling problem with sequence dependent setup times has attracted the attention of many researchers. For this problem, the sequence dependent means that the setup times of the processing jobs depend both on the preceding job and the job to be processed. It is different from the sequence independent which depends only on the job to be processed [2,3]. Meanwhile, the setup time is regarded as negligible in a lot of cases thus it is often ignored or dealt with a part of processing time [4][5][6][7].
Recently, many heuristic methods have been introduced to solve this problem since exact methods [8] may take too much computation cost. However, the exact methods are applied widely in the hybridization algorithms for many combinational optimization problems. Such as, Nishi and Hiranaka employed the Lagrangian relaxation and cut generation technique for solving sequence-dependent setup time flowshop scheduling problems to minimise the total weighted tardiness effectively [9], and in the book about the use of metaheuristics in optimization [10], the chapters [11,12] surveyed the integer linear programming techniques and metaheuristics derived from branch and bound for the combinatorial optimization respectively. Furthermore, Andreagiovanni and Nardin [13] proposed an improved ant-colony algorithm for designing sensor networks in 2015, where linear programming relaxations were employed to take better variable fixing decisions. Andreagiovanni and et al. also raised a very effective and fast algorithm [14] integrated an ant colony-like algorithm with an approximation algorithm and linear relaxations for solving multiperiod network design. Besides, Gambardella and et al. [15] presented how ant colony algorithms combined with solution improvement phases following some simple but effective rules can lead to good improvements in the quality of obtained solutions. In terms of the heuristic algorithms, the advantage of them is that some excellent solutions can be gained quickly, contrarily the disadvantage is that the quality of the solutions can not be guaranteed in general. There are three types of heuristics including constructive heuristic, meta-heuristic and hybrid heuristic respectively to deal with FSSP-SDST in the literature. A constructive heuristic comprises a stochastic construction and a greedy construction, and constructs feasible solutions by adding jobs one by one to the partial permutation. NEH algorithm was proposed by Nawaz et al. [16] and proved effective in finding the good solution for FSSP-SDST. Based on it, Rios-Mercado and Bard [17] proposed an extension constructive deterministic heuristic (NEHT_RMB) with high effectiveness for FSSP-SDST with the makespan objective. Meanwhile, a greedy randomized adaptive search procedure (GRASP) was also presented to improve the fitness of solutions. Ruiz et al. [18] used two algorithms namely a novel genetic algorithm (GA) and a memetic algorithm (MA) based on the classical genetic algorithm and they performed better than that of Rios-Mercado and Bard [19]. Rajendran and Ziegler [20] proposed an ant colony optimization (PACO) to enhance the solutions for FSSP-SDST with the objective of makespan and total flowtime of jobs respectively using NEHT_RMB method. Another ant colony optimization algorithm (ACO) was developed by Gajpal et al. [21] for FSSP-SDST with makespan criterion. Tseng et al. [22] presented an inventive heuristic algorithm and compared it with an existing index heuristic algorithm. The migrating birds optimization meta-heuristic was proposed by Benkalai et al. [23] and its good performance was proved by related experiments. The last but not least type of the heuristic is the hybrid algorithm; it also contains many algorithms blended with local search methods. Rios-Mercado and Bard [19] developed a heuristic algorithm called the hybrid algorithm which was the modification of the four heuristics compared with three benchmarks presented by Simons [24]. Besides, a new heuristic method called iterative greedy algorithm (IG) was proposed by Jacobs and Brusco [25], which was applied widely in other heuristics and had shown better performance for FSSP-SDST in the algorithm developed by Ruiz et al. [26]. Based on the basic IG (IG_RS), an effective local search was added to it in the destruction and construction phases for FSSP-SDST. In the same way, that effective local search was also added to the original MA which was called MA LS to solve the same problem. By the experiments, it has to be noted that the above two hybrid algorithms, IG_RS LS and MA LS , outperformed others including earlier heuristics and that of Rajendran and Ziegler [27] in getting better solutions. In addition, another effective iterated local search algorithm (ILS) developed by Wang et al. was shown its high efficiency with the objective of makespan in FSSP-SDST [28]. However, most of the above algorithms still cannot find solutions with high quality in a reasonable computational time. In this paper, we propose a hybrid local search algorithm for the sequence dependent setup times flowshop scheduling problem with makespan criterion. The difference of our proposed algorithm is that in most of the above algorithms, the initial individuals are generated randomly to keep the diversity of the population and most of them are not of high quality and this leads to solutions of low quality in next generations. In HLS, a NEH based problem-specific method is applied to guarantee the diversity and quality of the initial solutions at the same time. Besides, some of the above methods only focus on the global search capability and on contrary close attention to the local search ability is paid in some algorithms such as ILS and ACO. Hence, we combine the global search with the local search to enhance the solution quality in HLS. Furthermore, without effective perturbation methods, some algorithms such as IG_RS and IG_RS LS , are easy to be trapped into a local optimum after a number of iterations of local search without improving the quality of solutions further. Thus, in HLS, an efficient perturbation method is carried out to guide the search towards another area of the solutions and achieve better exploration.
Based on the above analyses, the hybrid local search algorithm can be summarized as follows: first of all, an overall framework of HLS based on local search and perturbation operators is presented to solve FSSP-SDST. During the search, the NEH based problem-specific method is used to initialize the population. Then to improve the performance of current population, a global search with a perturbation is applied. Next, we use an insertion-based local search to enhance the solutions for FSSP-SDST. What's more, a further local search method is used to exploit neighbors with good quality if the individuals cannot be improved by single local search method. At last, an update method is applied to update the current population. In order to show the performance of our proposed algorithm, six state-of-the-art algorithms namely GA, MA, MA LS , IG_RS, IG_RS LS and PACO are compared to test its high effectiveness. For the purpose of experimentation, we use the Taillard benchmark problems for FSSP-SDST at four levels including SDST10, SDST50, SDST100 and SDST125 with 12 sizes of data to analyze the performance of all the above seven algorithms. The statistical results show the high efficiency and competitive performance of HLS.
The rest of this paper is organized as follows: Section 2 describes the flowshop scheduling problem with sequence dependent setup times. The detail of the proposed algorithm is given in Section 3. In Section 4, the computational experiments are conducted. Finally, we draw the conclusion and the future work in Section 5.

The Flowshop Scheduling Problem with Sequence Dependent Setup Times
For the flowshop scheduling problem with sequence dependent setup times, a set of n jobs or tasks in J = {1, 2, 3, ..., n} has to be processed on a set M = {1, 2, 3, ..., m} of m machines sequentially, such as first on machine 1, then on machine 2, till the end on machine m. The objective of this problem is to find a processing sequence with the minimized makespan (C max ). For all the machines, the operating order of each processing sequence is identical. For a job permutation sequence π = {π 1 , π 2 , ..., π n }, the value of the setup times depends on the current job processing on the certain machine and the next job in the same permutation. Each job can be processed on one machine every time and meanwhile each machine can process only one job. All n × m operations are non-preemptive and none of them on each machine has the priority. The definition of the flowshop scheduling problem with sequence dependent setup times can be described as follows: where S j,π i ,π k is denoted as the setup time when producing job π k on the machine j after having processed job π i . p π i ,j denotes the processing time of job π i on machine j. The C(π i , j) is the complete time of job π i on machine j. In this paper, we assume that S j,0,π i = 0, C π i ,0 = 0 and C π 0 ,j = 0 (j ∈ M, i, k ∈ J).
Next, in order to show the problem simply, a problem with three jobs and three machines are introduced. The setup times and processing times are presented in Tables 1 and 2 respectively. A scheduling chart is outlined in Figure 1. From this figure, we can calculate the value of makespan is 20.

Overall Framework
In this section, the hybrid local search framework for solving FSSP-SDST with the objective of makespan is proposed by integrating some problem-specific methods. Each operator plays an important role in improving the efficiency of HLS.
At the beginning of the proposed algorithm, a NEH based problem-specific heuristic is utilized to initialize the population with N individuals. In other words, N permutation sequences for FSSP-SDST are formed. Then an improved population with good neighbor individuals is developed by the global search method. In addition, an effective single insertion-based local search method is applied to each individual for exploiting good individuals around it. If a better individual is found, the previous one will be updated. Otherwise, the multiple insertion-based local search is executed deeply as a further search regarding the individuals that are trapped into the local optima. Then, the solution will be stored in PotentialInd to maintain the good individual completely. In order to explore individuals transversely, the heavy perturbation approach is used to partly guide the trend of search. At last, we apply an effective update method to update the current population. The detail schemes of HLS are given in Algorithm 1.

Input:
The population size N. The number of jobs n. The maximum times of further search FSNum max .

Output:
The best individual s best .

Initialization
Regarding the initialization stage, it is divided into three parts including the population initialization, PotentialInd initialization and the initialization of s best . First, for the population initialization, it is generated randomly. Then, in order to produce some high quality solutions, the NEH based problem-specific method is applied to initialize the population which is composed of N random individuals. After that, HLS can provide some high quality solutions at the beginning of the search process. NEH method can be presented as follows in detail: Select a solution s from P c randomly, the first two jobs are extracted and two partial possible schedules with these two jobs are evaluated. The better one based on makespan is chosen as the current sequence. (2) For each unscheduled job j in s, put it into all the possible positions of the current scheduled sequence to generate all the possible partial sequences. The best one is selected as the current sequence for scheduling the next job.
A new individual s new is formed after scheduling all jobs, if s new is better than s, then it will replace s.
With respect to NEHBPS, for N times of iterations, a random individual is chosen to carry out NEH. Not all the individuals are processed by NEH, hence, the diversity of the initial population can be ensured. In summary, a population P c is generated with the feature of quality and diversity.
Secondly, HLS will initialize PotentialInd by assigning it a determinate sequence to store some problem-specific heuristic information. It is promised that a good individual is stored completely to get ready for finding the best individual during the search process.
Lastly, we initialize s best using the first individual in the population. This individual is updated via the iterations of search process in HLS. At the end of HLS, s best will be outputted with the minimum makespan for FSSP-SDST.

Global Search Method
The goal of this method is to gain an improved population by perturbing the good individuals and then updating them. In this subsection, the top 25% individuals are chosen as a new elite population to carry out the light perturbation. This setting aims to explore large search space using the appropriate solutions. Furthermore, we update the current population by using the individual which has undergone the light perturbation for enhancing the fitness of the population overall. The framework of GlobalSearch is shown in Algorithm 2.

Input:
The number of jobs n. The current population P c . The maximal number of GlobalSearch GSNum max . The mutation probability MP.

Output:
The improved population P b . 1: i = 0, j = 0; 2: Sort individuals in P c by an ascending makespan order. 3: while i ≤ GSNum max do 4: for each individual s in the first 25% range of P c do 5: while j < MP × n do 6: s = LightPerturbation(s); 7: P c = U pdateInds(s ); 8: if C max (s ) < C max (s) then 9: s = s ; 10: end if 11: j = j + 1; 12: end while 13: end for 14: Sort individuals in P c by an ascending makespan order; 15: i = i + 1; 16: end while 17: P b = P c ; 18: return P b This method allows the perturbed solution to maintain some characteristics of the previous solution and find other characteristics in new solutions for improving the diversity of the population. It also helps HLS to find other optimal solutions in the neighborhoods. The light perturbation strength is so sufficient to lead the search trajectory to another neighbor region which can result in a different solution. In conclusion, GlobalSearch can both enhance the diversity of the population and improve the convergence rate of HLS.

Update Method
In this section, two algorithms, namely Algorithms 3 and 4, are presented. In Algorithm 3, we use the specified individual to update the individual randomly in the current population. At the same time, the update population method is presented in Algorithm 4. After achieving it, an elite population including some individuals with top fitness has acquired.

Input:
The current population P c . An individual to be used for updating s.

Output:
The current population with updated individuals P c . 1: TempP c = P c ; 2: while TempP c = ∅ do 3: Randomly select an individual s in P c ; 4: if C max (s) < C max (s ) then Randomly select a position k in t; 5: t = t \ k; 6: if C max (s i ) < C max (s k ) then 7: s k = s i ; break; 8: end if 9: iter = iter + 1; 10: end while 11: return P c Except for the above two update methods, another update method is used for PotentialInd. As Algorithm 1 dedicates, if the current solution has improved successfully in the single insertion-based local search or the further search section, PotentialInd will be updated.

Perturbation and Local Search Methods
In the single insertion-based local search and the further search, if the current individual has been replaced with a better individual, then the better one is stored. Otherwise, it means that this solution has been trapped into a local optimum, the heavy perturbation will be employed to the current solution. Compared with the light perturbation, the heavy perturbation has larger interference strength to provide enough perturbation for the high quality solutions. The heavy perturbation is used to obtain some characteristics of the previous solution, but it isn't apt for all the local search methods. Such as in GlobalSearch, the light perturbation method is more suitable as it has a better result for disturbing the good solution with the relatively fitting strength. The light perturbation and the heavy perturbation operations can be defined as follows: The heavy perturbation: Input a solution s. Three different positions p 1 , p 2 , p 3 of s are randomly selected, where Let S 1 represent the partial sequence between p 1 and p 2 and S 2 denote the other partial sequence between p 2 and p 3 (not including p 2 ), then exchange S 1 and S 2 to generate a new solution s new to be outputted.
The light perturbation: (1) Input a solution s. Two different positions p 1 , p 2 of s are randomly selected, where p 1 < p 2 .
Let S represent the partial sequence between p 1 and p 2 (not including p 1 ), then move the job in position p 1 behind S to generate a new solution s new to be outputted.
On the first phase, a single insertion-based local search is employed to HLS for exploring better solutions. In [29][30][31][32][33][34][35][36][37], this local search has shown its high effectiveness in obtaining the better quality solutions than the local search based on swap operators. The curial idea of this method is to compare the previous makespan with the new makespan of the sequence, which has inserted a job into all positions of the sequence. Then the better one will be accepted. But if a solution cannot be improved, it will experience the further search to make it rope out of the local optimum on the second phase. Moreover, the further search with the purpose of exploiting good individuals in the neighborhood can increase the convergence speed of HLS. The effective insertion-based local search is described as follows: Input a solution s and a position r.
Let s best be the best solution based on the minimal makespan among the n − 1 neighborhood solutions. (4) Output s best .

How to Balance between Exploitation and Exploration
For designing the hybrid heuristic algorithm, it is known to all that diversity and convergence are the two basic issues. In terms of diversity, exploitation is executed in HLS for improving the convergence rate. Otherwise, explorations to other search directions on breadth is applied to HLS for maintaining or enhancing the diversity of the population. HLS can trade off the convergence and diversity efficiently where the diversity is pursued by the exploration and the convergence is increased in the manner of the exploitation. (1) Regarding the exploitation: It means that a further search in the large and deep search space can be executed. For the purpose of searching better solutions, firstly we have used the single insertion-based local search to begin the evolution of population. Furthermore, if the current individual cannot be replaced with the new one after carrying out the single local search, thus secondly we choose the further search based on the insertion-based local search to search the solutions deeply. After above methods, the rate of convergence can be improved quickly. (2) Regarding the exploration: In other words, the central idea of it is to maintain the diversity of the population. In the search process, we have applied the heavy perturbation method to the individuals which are trapped in the local optima. This way effectively avoids the current population being trapped into the local optima and urges that population to an anticipant direction for next generation. It is particularly worth mentioning here that the light perturbation in GlobalSearch keeps the diversity of the improved population, which has a high convergence for searching good solutions at the same time.
As we expect, a good population with high quality and diversity will be formed. Moreover, at the foundation of the above analyses, HLS can balance the exploitation with exploration in the search process successfully.

Environmental Setup
HLS is performed on a PC with 4 GB RAM and a CPU of 3.40 GHz on Windows 7 OS. It is programmed in C++ by Microsoft Visual Studio 2013. It is obviously known that the algorithm with the more running time can bring out the better result. To compare with other algorithms fairly, the stop criterion of HLS is CPU time limit given by (n × m/2) × f milliseconds. Besides, f is set as 30, 60, 90 respectively, which is the same as Ruiz used in [38].

Benchmark Problem Instances and Benchmark Algorithms
To verify the performance of the proposed algorithm, Taillard-based sets from Vallada et al. (2003) including 4 sets and 480 benchmark instances are used. The instances in each set range from 20 jobs with 5 machines to 500 jobs with 20 machines, consisting of 20 × 5, 20 × 10, 20 × 20, 50 × 5, 50 × 10, 50 × 20, 100 × 5, 100 × 10, 100 × 20, 200 × 10, 200 × 20 and 500 × 20 respectively. Besides, each size has 10 specific examples. It's remarkable that each set is different in processing times and setup times. The first two sets are SDST10 and SDST50, in which the setup times are 10% and 50% of the average processing times p π i ,j . In other words, the setup times are generated equably from the distribution ranges [1,9] and [1,49] because p π i ,j in this benchmark are uniformly generated in the distribution range [0, 99]. Similarly, the last two sets are SDST100 and SDST125, in which the setup times are generated uniformly from the distribution ranges [1,99] and [1,124], respectively 100% and 125% of p π i ,j .
Then, to show the effectiveness of HLS, six efficient algorithms including GA, MA, MA LS , PACO, IG_RS, and IG_RS LS on FSSP-SDST [18,38] are used to compare with our proposed algorithm. Furthermore, a response variable called the relative percentage deviation (RPD) is applied to show the increase between a specific solution produced from a certain algorithm and the best bound found in http://soa.iti.es/problem-instances.

Experimental Parameter Settings
HLS has five control parameters, comprised of N, MP, GSNum max , UpdateNum max and FSNum max . Instances Ta071-Ta080 (100 × 10) in SDST50 are used to calibrate the parameters as the base case due to space limitations in the following experiments. Each experiment runs ten times independently with CPU time equaling to (n × m/2) × 30. The good results relatively are marked in bold.

The Influence of the Population Size N
It is obvious that the population size is an essential parameter to heuristic algorithms. Selecting an appropriate size has been known as a challenging and puzzling task. If the size is too small, it is easy for the algorithm to have a fast convergence which leads to local optima prematurely. However, using too large size will bring about additional computation costs. To get the influence of the population size on the search performance distinctly, the population size has been varied from 10 to 40 with the RPD results summarized in Table 3. From the average results of RPD in Table 3, the population size 20 can provide better average RPD (ARPD) in some instances including Ta071, Ta075, Ta078 and Ta080. In Ta078, a population size 10 provides the same least ARPD as a population size 20 with the value of 0.41. Besides a population size 10 and 30 can bring out the same least ARPD (0.56) in Ta079. Moreover, HLS with a population size 30 can gain a better RPD in instances Ta074, Ta077 and Ta079. Figure 2 gives an empirical insight of influence on each population size using the confidence intervals at the 95% confidence level of ARPD for Ta071-Ta080 and it is observed that the population size 20 beats other settings with the best ARPD. Therefore, it is concluded that the population size 20 outperforms other settings in producing better ARPD results.

The Influence of the Maximal Number of GlobalSearch GSNum max
In this experiment, we will discuss the effect of GSNum max . In GlobalSearch, GSNum max indicates the maximum number of iterations that conduct the search on the selected individuals who can be further interfered. To demonstrate the effect of GSNum max on the proposed algorithm, Ta071-Ta080 are also applied to present the performance measured as RPD with GSNum max chosen from set {3, 4, 5, 6}. Other parameters temporarily set as the previous experiment in this experiment. Table 4 summaries RPD results of each GSNum max . From Table 4, it is easily observed that GSNum max = 3 gains only one best result in Ta071, while GSNum max = 4 can perform better in Ta075, Ta076, Ta078, Ta079 and Ta080. What's more, GSNum max = 5 provides better results on three instances including Ta073, Ta074 and Ta077, however only in Ta072, GSNum max = 6 beats other settings. In general, GSNum max = 4 gives the best average RPD on instance set Ta071-Ta080. Figure 3 provides the 95% confidence interval error graph in terms of ARPD for Ta071-Ta080 with different GSNum max . We can find that GSNum max = 4 has a strong adaptability with higher convergence ability and it is good in solving FSSP-SDST with the best ARPD among all the settings. Based on the above analyses, it is summarized that GSNum max = 4 results in the better performance for most instances in terms of RPD. Therefore, GSNum max = 4 is the best choice for HLS.

The Effect of the Mutation Probability MP
As HLS introduced before, in GlobalSearch, a mutation probability was used to control iteration times to implement the perturbation and update process. It is so important that too large MP will bring a move at angles to the space with better results and small values of MP will guide the search to the parallel axes in the exploring space. Hence, a set of MP ranging from 0.1 to 0.4 is used to present the relative performance in getting better results with Ta071-Ta080. Table 5 details the results with different values of MP.
It is shown in Table 5 that MP = 0.1 can gain best RPD only in Ta075, while MP = 0.2 provides better RPD in instances set Ta071, Ta074. Besides, MP = 0.3 can give better results in most of the instances including Ta072, Ta073, Ta076, Ta079 and Ta080 with the smallest average RPD. For Ta077 and Ta078, MP = 0.4 performs better results among all the settings. ARPD concluded in average 95% CI is presented in Figure 4. Generally speaking, the experiments implies that MP = 0.3 might be a good setting to solve FSSP-SDST with the better ARPD.

The Influence of the Maximal Number of Update UpdateNum max
In the update method UpdateP c of HLS, we use UpdateNum max to present the maximal iteration number of update taking advantage of an individual separately. To maintain the diversity of the population, the value of it can't be too large and it should be less than the population size N according to the framework of UpdateP c , however if it was too small, the update method has little effect on gaining better result in the search space. A certain number is not apt for all instances. Here UpdateNum max is set as 4, 6, 8, 10 to test its suitability. Other settings are as the preliminary experiments. Table 6 shows the detail results.
From Table 6, UpdateNum max = 4 is suitable for Ta072, Ta073, Ta079 and Ta080. For Ta071, Ta074 and Ta077, UpdateNum max = 6 beats others in getting better RPD. By the same way, Ta072 and Ta075 can give better results when UpdateNum max = 8 while UpdateNum max = 10 provides better RPD in instances Ta076 and Ta078. Herein UpdateNum max = 4 and UpdateNum max = 8 present the same results in Ta072. As can be seen in the average RPD, UpdateNum max = 6 outperforms other numbers of UpdateNum max . It can be found in Figure 5, an average RPD plot with the confidence intervals at the 95% confidence level is shown and it is easy to find out that UpdateNum max = 6 performs best among all the settings. There can be concluded that different numbers of UpdateNum max will affect the performance of the proposed algorithm variously and we set UpdateNum max to 6 for its good achivement.   If some individuals cannot be improved in the single insertion-based local search, then a further search method will be applied to them. In this section, the effect of FSNum max used to define the maximal number of the further exploitation will be discussed. We have taken advantage of Ta071-Ta080 measured as RPD to present the effect. The value of it varies through the set {10, 15, 20, 25} and other parameter settings are the same as the previous experiments. The results of the settings are summarized in Table 7.
From Table 7, FSNum max = 10 provides good RPD on Ta072, however for Ta078 and Ta073, FSNum max = 20 and FSNum max = 25 perform better than other settings respectively. For the rest instances including Ta071, Ta074, Ta075, Ta076, Ta077, Ta079 and Ta080, FSNum max = 15 is a good setting beating other values in gaining better RPD. It is obvious that FSNum max = 15 outperforms others on most of the instances. Besides, Figure 6 demonstrates ARPD with 95% confidence interval of all the ten instances. It can be easily seen that FSNum max = 15 has the best result for solving the instances in HLS.

Effect of NEH Based Problem-Specific Heuristic
In the initial stage of the proposed algorithm, a NEH based problem-specific heuristic method is used to initialize the population which is produced by a random method. It plays an important part in enhancing the results of the algorithm. Without it, the search of getting good individuals will spread out in a quite broad direction. We respectively use the instances Ta071-Ta080 (100 × 10) in SDST50 and SDST125 ran ten times to demonstrate the vital effect of it. Besides, the rest frameworks are the same in order to compare fairly. Table 8 shows the results of instances in SDST50 and SDST125 clearly. The better RPD results are marked in bold.
From Table 8, with the exception of Ta071, RPD of other instances in SDST50 using the method with NEHBPS are all better than those without NEHBPS. Besides NEHBPS on all the instances of SDST125 behaves the same as that on SDST50 except for Ta073. It means that RPD of the instances in SDST125 taking advantage of the method with NEHBPS outperforms the method without NEHBPS except Ta071 and Ta073. ARPD plot with 95% confidence interval of SDST50 and SDST125 is directly shown in Figure 7. It has demonstrated that if NEHBPS is not adopted, the performance of HLS seriously degrades. In summary, the NEH based problem-specific heuristic method has affected the proposed algorithm to a great extent. HLS that making use of it to initialize the population is a great strategy.

Effect of Different Perturbation Operators
In GlobalSearch, a perturbation has applied to disturb some individuals which trapped in local optima of the current population in terms of balancing the exploration and exploitation. It has an important role in finding individuals with high quality and reducing the further local search. However, there are many kinds of perturbation operators with different effects on the proposed algorithm. Here, two operators consisting of light and heavy perturbation mentioned before are compared by the experiments using instances Ta071-Ta080 in SDST50 and SDST125. In order to compare fairly, it is necessary to make sure the rest frameworks are the same. The effect of the two perturbations details in Table 9 which the better results are marked in bold.
From Table 9, for the instances in SDST50, the heavy perturbation method performs better in Ta072, Ta074, Ta078 and Ta080. In contrast, the light perturbation behaves better in the other seven instances. In terms of the instances in SDST125, the light perturbation also performs better in most of the instances including Ta072, Ta073, Ta074, Ta075, Ta076, Ta077, Ta078 and Ta080. However, the heavy perturbation beats the light one only in instances Ta071 and Ta079. Figure 8 can provide the direct comparison between the two methods. It is proved that the light perturbation exhibits the superior performance in gaining better overall ARPD compared with the heavy perturbation. To conclude, the heavy perturbation is a quite efficient method to help the individuals quit the local optimal condition.

Effectiveness Evaluation of Different Local Search Operators
In the beginning search process and the further search of the proposed algorithm, an insertion-based local search local search was executed to find individuals with better fitness in the neighbour search space. Different types of local search can guide the individuals to different directions no matter they are good or bad. In this section, two other local search operators are used to compare with the effective insertion-based local search by instances from Ta071 to Ta080 (100 × 10) in SDST50 and SDST125 running ten times. The compared two operators can be described as follows.
(1) exchange_based local search: For the positions ranging from 1 to n − 1 of the scheduling permutation, exchange the job in it with the job in the adjacent next position. (2) swap_based local search: For every position of the scheduling permutation (from 1 to n), swap the job in it with the job in the given position.
It is noted that the other frameworks are the same in this comparison to make the comparison fair. Table 10 demonstrates the details of the three operators. The better results are marked in bold. There are definite differences among these three local search operators. For the standard deviation (SD) of these instances in SDST50 and SDST125, the insertion-based local search has the smallest standard deviation 0.26 and 0.47 respectively among these local search methods. It means that the insertion-based local search is much stabler in getting high-qualified solutions. Furthermore, ARPD that obtained by the insertion-based local search of HLS is much smaller. Hence, it is obvious that the insertion-based local search which can generate high quality solutions performs better than other two methods on all the instances in SDST50 and SDST125. Figure 9 demonstrates that the insertion-based local search in the proposed algorithm can yield better overall mean namely ARPD. In conclusion, the insertion-based local search exhibits its superior performance with high efficiency in the proposed algorithm.

Comparison Results with Some State-Of-The-Art Approaches
In this subsection, the proposed algorithm with the best algorithm settings yielded from preceding sections is tested with six compared benchmark algorithms containing GA, MA, MA LS , IG_RS, IG_RS LS and PACO which can be found in the literature [38] to present a comprehensive performance comparison. Among them IG_RS LS and MA LS are integrated IG_RS and MA respectively with an effective local search. To examine the relative ranking of all the algorithms which are dependent on the computation time with the comparative performances, three series of experiments setting the f aforesaid at 30, 60 and 90 are carried out. For each benchmark instance in the following experiments, HLS runs independently ten times and thirty times. HLS(10) denotes the running time of HLS is ten, while HLS (30) represents that HLS is running thirty times. RPD is calculated to compare the performance of different algorithms for each instance. Besides, the standard deviations for each type of instances are calculated to show the robustness of the compared algorithms. The comparative results are shown from Tables 11-16 grouped by instances type and size. Because of the limitation of space, only the average RPD of problems in 12 different scales is presented in the sets including SDST10, SDST50, SDST100 and SDST125. It means that ARPD values for each scale with 10 instances, ARPD and SD of all the scales in each type of instance sets are shown in the corresponding tables. The best results of ARPD for the specific algorithm are marked at bold. Besides the negative values represent that the results found by HLS are better than the current best solution provided in http://soa.iti.es/problem-instances. It has to be noted that in the following analyses the ARPDs are achieving by HLS running ten times and the improvement of ARPD means the decrease of that value compared with a certain value because the lower value of ARPD represents the better performance of the algorithm. From ARPD through Table 11, which the computation time is (n × m/2) × 30, the algorithms except for HLS, the best cross ARPD are all gained by IG_RS LS . It is said that the IG_RS LS owns the optimal performance in the compared algorithms. In SDST10, HLS can provide slightly better results for instances in the scales from 20 × 5 to 200 × 20 than the algorithm IG_RS LS . Among them, the maximum improvement compared with IG_RS LS is 0.41 in the scale of 100 × 20 and the minimum improvement is 0.01 for the scale of 20 × 10. As for the instances in large scale 500 × 20, IG_RS LS gives a better result of ARPD which is 0.46 than the ARPD of our proposed algorithm 0.72. For the rest algorithms, HLS performs better than all of them in obtaining better ARPD which has the lower value. Overall HLS is better than IG_RS LS that is the suboptimal algorithm with an ARPD improvement of 0.17. In SDST50, the promotion of ARPD which is developed by HLS has a bit improvement than the type of SDST10. Except for the scale of 500 × 20, HLS presents the best performance on each scale among all the algorithms. But in terms of ARPD in the scale 500 × 20, HLS is just next to IG_RS LS and is superior to the rest algorithms owing the second-best performance. The least improvement of ARPD is 0.04 and the most improvement is 0.99 compared with the suboptimal RPD obtained by IG_RS LS . To the overall ARPD, HLS beats all the other algorithms on average with the value of 0.45, better than ARPD 1.00 of IG_RS LS . The ranking order of the algorithms as the best performing to the worst performing are HLS, IG_RS LS , MA LS , MA, IG_RS, PACO and GA. From Table 12, for the large sets SDST100 and SDST125, HLS can generate better results compared to the small sets including SDST10 and SDST50. For SDST100, the least difference which is better than IG_RS LS is 0.19 which is the difference between 0.27 and 0.08 in the scale of 20 × 20. The above data records the most difference is 1. 30. The values of ARPD in each scale of instances gaining by HLS are all the best with the exception of the scale 500 × 20. Considering the cross average in this table, the ranking of HLS is best. For SDST125, the least and the most difference which are improved compared with IG_RS LS are 0.19 and 1.58 respectively. Seen from the overall ARPD, there is an ARPD 0.62 developed by HLS better than the other algorithms and it has the value of 0.82 in terms of the improvement compared with the suboptimal ARPD in the algorithm IG_RS LS . It is noted that it has a much higher average percentage promotion than other five algorithms. It also implies that HLS can provide better performance with the increase of ratio between the setup times and the processing times.
From Tables 13 and 14, which are ARPD results of all the types in all the instances when we set the computation time as (n × m/2) × 60, the same pattern are obtained in the four data sets. With the exception of the instances in 500 × 20, all other ARPD are improved rather better gradually with the increase of the setup times in the range of SDST10, SDST50, SDST100 and SDST125. It is clear that HLS exceeds all other algorithms because of its good performance on the least overall ARPD.  Tables 15 and 16 when the computation time is (n × m/2) × 90. For SDST50, except for the scale of 500 × 20, ARPD of each algorithm in each scale decreases as the order of GA, PACO, IG_RS, MA, MA LS , IG_RS LS , HLS, which shows the best performance of HLS among all the compared algorithms. For SDST100 and SDST125, it is shown the same pattern. However, the tendency of increasing softens in SDST10 with the exceptions of ARPD in the scale of 20 × 10 and 500 × 20. For brevity, the remaining analyses of the minimum and maximum of ARPD between the relatively good two algorithms which are HLS and IG_RS LS in each scale are omitted. But the trend of performance for each algorithm is the same as the previous tables has shown. As we can see from the overall ARPD, it is evident that there is the least ARPD in HLS which surpasses by other algorithms attributed to the high effectiveness of each component in HLS. Furthermore, via the standard deviations from Tables 11-16 of each compared algorithm, it is demonstrated that the high stability of HLS in gaining good solutions.  In addition, it is obvious that each algorithm is benefited from additional computation time in all the types of instances including SDST10, SDST50, SDST100 and SDST125. Figures 10-13 are used to prove the importance of the computation time in different scales of instances by ARPD evaluated with HLS running ten times. The trajectory of each curse in the figure overall declines and becomes closer to the X axis along with the increase of f which represents the better cross ARPD. Besides, it is also found that the average RPD shows an increasing trend because the increase of the instance size in the number of machines classified by the number of jobs including 20, 50, 100, 200 and 500. Via the above figures, what can be summarized is that the value of the cross ARPD is decreasing and becomes better clearly with the increase of f regardless of the instance type and instance size.

It is the same as
Besides, the interactions between the ratio of the setup times in the processing times and the scales have to be studied to verify the effectiveness of HLS. It is shown in Figures 14-16 to present the different types of instances segregated by f = 30, 60 and 90 for HLS running ten times respectively. As a matter of fact, the average RPD has reduced rather more acutely with a larger difference in each algorithm for the large type of SDST100 and SDST125 compared with the type sets of SDST10 and SDST50. Moreover, the trend of each curve in Figures 15 and 16 demonstrates the performance ranking of these compared algorithms are as the ascendent order of GA, PACO, IG_RS, MA, MA LS , IG_RS LS , HLS visually. However, Figure 14 has presented the descendent performance order of algorithms is HLS, IG_RS LS , MA, IG_RS, MA LS , PACO and GA. It is summarized that HLS gives higher priority to the increase of the computation time and the increase ratio between the setup time and the processing time.  The next tables from Tables 17-23 present the newbounds of makespan discovered by HLS running ten times and the aforementioned best-known bounds grouped by different scales in different types of instances. The differences between the two values are presented in the tables which shows the high effectiveness in getting high-qualified solutions of HLS. Besides, they become larger as the increase of the ratio of the setup times in the processing times and the scales of instances. By the same pattern, the number of newbounds has raised progressively. For an example, in the scale 50 × 5, the difference in Ta039 of SDST10 is −2, but for Ta038 of SDST100, the difference is −27. Table 18 provides the newbounds, the best-known bounds and the differences between them for the instances in the scale 50 × 10. Hence, it is evident that the differences are bigger in SDST50 and SDST100 overall which mean that HLS is much more effective. The tables are shown that HLS is statistically superior to other well-performing algorithms. For the illustration purpose, what follows it is Table 24 that shows some permutations of jobs chosen from the scale 50 × 5 in the instances of SDST10 and SDST100 because of space limitations. The permutation sequence begins from 0 and ends in 49 in each instance which is the best solution with the newbound of makespan developed by HLS.   To illustrate the robustness of HLS, the instances in the set SDST50 are chosen to be the base examples. Besides, Figures 17-19 use ARPD in confidence intervals at the 95% confidence level taking the CPU time of (n × m/2) × 30, (n × m/2) × 60 and (n × m/2) × 90 respectively to analyse the robustness. The ARPDs here are developed by HLS running ten times. The length of the error bar represents the stability degree of each algorithm. It means if the length of the bar for one algorithm is short, then the algorithm gives the better robustness. It is clearly shown that HLS has the strongest robustness with the shortest interval among all the compared algorithms no matter how long the CPU time is. The robustness order of other algorithms are PACO, GA, IG_RS, MA, MA LS , IG_RS LS from worst to best as the above figures have provided. Moreover, the above figures have presented the best effectiveness of HLS with a smallest value of the average RPD. As a result, the average RPD also shows that the effectiveness ranking of other algorithms is GA, PACO, IG_RS, MA, MA LS , IG_RS LS from the last rank to the first rank. Although PACO has the better effectiveness which has a lower ARPD compared with GA, in terms of robustness GA is stabler than PACO. Through synthetical consideration of the effectiveness and the robustness for each algorithm, HLS is resulted to the best algorithm in all the compared algorithms.        36 6 18 20 42 40 5 22 32 44 2 43 11 28 47 31 12 10 41 3 As we can see from the above tables and figures, it can be concluded that it is GA, PACO, IG_RS, MA, MA LS , IG_RS LS , HLS from the worst to the best performance for the comprehensive consideration particularly because ARPD of all the 480 instances is decreasing in line. Although in the instances of the scale 500 × 20, IG_RS LS can provide a better ARPD, HLS can present a better overall ARPD regardless of instance type. Besides, HLS performs much better as the increase ratio of the setup times in the processing times. Based on the above results, it can be said that HLS performs better than other six algorithms. It is verified that HLS is an effective and robust algorithm for solving the flowshop scheduling problem with sequence dependent setup times.

Comparison Results with Some Recent Algorithms
In order to further illustrate the effectiveness of the proposed algorithm, we reconstruct an experiment to compare with two recent algorithms, namely the adaptive hybrid algorithm (AHA) [39], and the enhanced migrating birds optimization (EMBO) [40] for FSSP-SDST with the makespan criterion.
In AHA, each job is assigned an inheriting factor. For dynamically updating the factor, a novel operator is constructed. Therefore, both good and bad genes can be explored. A new crossover operator is proposed by inheriting good genes to the offspring and destroying the bad genes with a high probability. Hence, the offspring is integrated with more and more good genes generation after some generations. It helps to improve the effectiveness of AHA in obtaining high-qualified solutions. The algorithm AHA is chosen to compare with HLS. Besides, for the migrating birds optimization (MBO) [41], it is a metaheuristic inspired from the flight of migrating birds. To save the energy, the known V-flight shape3 is developed. Thereafter, a basic migrating birds optimization (BMBO) [23] is proposed to solve FSSP-SDST with the makespan criterion. In terms of EMBO, since the performance of BMBO keeps decreasing with the increase in the size of instances, EMBO is based on the BMBO with an or-opt neighbourhood which was designed for the travelling salesman problem (TSP) first and the well-known heuristics to generate its leader bird. The or-opt neighbourhood operator is moving a block of one, two, three jobs or four jobs and inserting it elsewhere in the sequence.
HLS runs 30 times to compare with the above algorithms. The RPD results for each scale data in SDST10, SDST50, SDST100, SDST125 of each compared algorithm is presented in Table 25. The relatively good results are marked in bold. From Table 25, it is seen that HLS is more suitable than the other compared algorithms for achieving good results. In detail, HLS has the best values of ARPD 0.31 and 0.40 in SDST10 and SDST50 respectively among the ARPDs gaining by the compared algorithms. For SDST100 and SDST125, the values of the overall ARPD are 0.47 and 0.52 which are clearly better than the other ARPD values of EMBO and AHA. In addition, Figure 20 is presented to show the good effectiveness of HLS in obtaining high-qualified results intuitively.
In conclusion, HLS is very competitive to the good performing algorithms for solving FSSP-SDST with makespan criterion.

Conclusions and Future Work
In this paper, the flowshop scheduling problem with sequence dependent setup times is addressed with minimizing the makespan. A hybrid local search algorithm based on novel local search methods is presented to deal with this problem. First, to initialize the population, we apply an effective NEH based problem-specific method to the initialize the population. Second, the global search embedded with a light perturbation is used to generate the better population. Then, to find good individuals in the current population, the insertion-based local search is adopted. Next, a further local search is applied to individuals which are trapped into local optima. Last, a heavy perturbation is used to explore better neighbours in the new research regions.
In order to demonstrate the performance of the proposed algorithm, extensive experiments are conducted on all the 120 instances of four different scales. It is obviously shown in the compared experiments that HLS gives better results than other six state-of-the-art algorithms, including GA, MA, MA LS , IG_RS, IG_RS LS and PACO. The experimental results confirm that it is appropriate in using the combination of local search and perturbation methods in HLS.
In the future, the proposed HLS algorithm is expected to be used to solve other combinational problems such as no-wait flowshop scheduling problems (NWFSSP), hybrid flowshop scheduling problems (HFSSP), blocking flowshop scheduling problems (BFSSP) under makespan criterion. Different values of parameters for HLS can be experimented comprehensively to enhance the solution quality. In addition, different neighbourhood structures can be developed for improving the effectiveness of HLS. Meanwhile, the interactions between the exploitation and the exploration strategies should be investigated in depth. In other words, other techniques on different ways of local search and perturbation which is more suitable for FSSP-SDST can be studied to guide the search to other extensive spaces and enhance the performance of HLS. Moreover, it helps to construct other metaheuristics with some effective ingredients adding to the basic HLS. For example, the restart strategy that a new population is generated randomly with the current best individual not improving after a given number of iterations and the elite strategy that keeps some elite individuals in current population for a number of iterations.
Furthermore, we can integrate HLS with other heuristics, such as the variable neighbourhood search and the tabu search. This would be more effective to keep balance between the intensification and the diversification and improve the quality of solutions.