Hybrid Algorithm Based on an Estimation of Distribution Algorithm and Cuckoo Search for the No Idle Permutation Flow Shop Scheduling Problem with the Total Tardiness Criterion Minimization

The no idle permutation flow shop scheduling problem (NIPFSP) is a popular NP-hard combinatorial optimization problem, which exists in several real world production processes. This study proposes a novel hybrid estimation of the distribution algorithm and cuckoo search (CS) algorithm (HEDA_CS) to solve the NIPFSP with the total tardiness criterion minimization. The problem model is built on the basis of the starting and ending time point of each job. A discrete solution representation method is applied in HEDA_CS to increase the operation efficiency. A novel probability matrix build method is also designed within the knowledge of the processing time matrix. The partially-mapped crossover operation works effectively during the CS phase. A suitable knowledge-based local search is also designed in the HEDA_CS to balance the exploitation and exploration. Finally, many simulations based on the new hard Ruiz benchmarks are conducted. Computational results demonstrate the effectiveness of the proposed HEDA_CS.


Introduction
The permutation flow shop scheduling problem (PFSP) has been the focus of many studies for decades.The no-idle constraint in scheduling occurs when two consecutive jobs must be processed on the same machine without any interruptions.Given that this constraint appears in real-world production environments, the no idle permutation flow shop scheduling problem (NIPFSP) has considerable academic and practical significance.Similar to the traditional scheduling problem, the NIPFSP is proven to be NP-hard [1,2].Exact optimization methods have limitations in solving large-scale problems because of the calculation time limitation.Therefore, developing effective and efficient algorithms to solve the NIPFSP is significant.
The NIPFSP was first studied by Adiri and Pohoryles [3], and it has received extensive attention ever since.Woollam [4] investigated a solution procedure for a flow shop problem, in which a machine continuously processes all the jobs that must be processed once it is started.That study described a "no idle time allowed" constraint.Saadani et al. [5] determined that the idle characteristic seriously affected the value of the makespan (Cmax) criterion in a three-stage, no idle flow-shop configuration.Given that the traditional exact optimization algorithm has limitations with the age of data exploration, heuristics has received increasing attention by solving the NIPFSP [6][7][8][9][10].Dong et al. [11] proposed an improved NEH-based heuristic with an initial sequence generated by combining the average job processing time.Their proposed strategy was based on the idea of balancing the utilization among all machines and exhibited suitable performance.Baraz and Mosheiov [12] introduced an efficient (O(n 2 )) greedy algorithm, which was shown numerically to perform better than other published heuristics.Kalczynski and Kamburowski [13] addressed the problem of determining a job sequence that minimized the makespan in m-machine flow shops under the no-idle condition.Pan and Wang [14] proposed a novel discrete differential evolution algorithm to solve NIPFSP.That study presented two simple approaches to calculate the makespan and a speed-up method to insert the neighborhood and thus improve the efficiency of the entire algorithm.Pan et al. [15] then proposed a new and novel referenced local search procedure hybridized with both algorithms to further improve the solution quality.The referenced local search exploited the space based on reference positions taken from a reference solution to determine better job positions when performing the insertion operation.Researchers [16] also solved the single machine total weighted tardiness problem with sequence dependent setup times by a discrete differential evolution algorithm.He and Wang [17] proposed a hybrid algorithm that combines evolutionary computation and constraint-handling techniques.Li and Wang [18] proposed a hybrid quantum-inspired genetic algorithm for the multi-objective PFSP.A new random-key representation was used to convert the Q-bit representation to job permutation in evaluating the objective values of the schedule solution.Deng and Gu [19] proposed a hybrid discrete differential evolution (HDDE) algorithm for the no-idle permutation flow shop scheduling problem with makespan criterion.Tasgetiren et al. [20] investigated the utilization of a continuous algorithm for the NIPFSP with tardiness criterion.A differential evolution algorithm with variable parameter search was developed to solve the NIPFSP.Researchers [21] also designed a variable iterated greedy algorithm with differential evolution to solve the NIPFSP in recent years.Pan and Ruiz [22] first proposed an effective iterated greedy algorithm for a mixed no idle flow shop, where several machines had the no-idle constraint, whereas others were regular machines.The researchers also proposed a formula set to accelerate the insertion calculation that was utilized in heuristics and in local search procedures.A novel nature-inspired cuckoo search (CS) algorithm was developed by Yang and Deb [23] in 2009.CS has a series of successful engineering examples [24][25][26].The improved CS algorithm was proposed for hybrid flow shop scheduling problems by Marichelvam et al. [27], and the algorithm was validated with the data from a leading furniture manufacturing company.A discrete version of the inter-species CS algorithm was proposed and applied to solve two significant types of the PFSP [28].
The current study is related to the complex production process with no idle tight constraints in the real world.Given that the starting and ending time points are a suitable method to describe the job process, a new problem formulation is built on the basis of this condition.Considering the last CS algorithm, this study proposes a novel hybrid estimation of distribution algorithm (HEDA) and CS algorithm (HEDA_CS) to solve the NIPFFSP with the total tardiness criterion minimization.A knowledge-based local search is applied to the proposed HEDA_CS to balance the exploitation and exploration of the HEDA_CS.Several latest Ruiz benchmark instances are adopted to test and improve the HEDA_CS performance.A suitable adjustment to the algorithm shows that the HEDA_CS is effective and highly efficient in solving the NIPFFSP with the total tardiness criterion minimization.
The remainder of the paper is organized as follows.Section 2 describes the problem formulation with a time point method to build the mathematical model.Section 3 presents the novel HEDA_CS to solve the NIPFSP with the total tardiness criterion minimization.Several of the latest Ruiz benchmark instances are employed to test the performance of the proposed HEDA_CS in Section 4. The computational results and analysis are also presented in this section.Finally, the conclusions are discussed in Section 5.

Problem Formulation
The permutation flow shop scheduling problem (PFSP) is illustrated in Figure 1.A total of n jobs J = {J 1 , J 2 , . . ., J i , . . ., J n−1 , J n } must be processed on m machines M = {M 1 , M 2 , . . ., M j , . . ., M m−1 , M m } with the same sequence.Consequently, n products {P 1 , P 2 , . . ., P i , . . ., P n−1 , P n } are attained.Therefore, determining the processing sequence of n jobs over m machines in PFSP to satisfy several objectives is widely applied in actual production, especially in one-piece mass production.The NIPFSP is a highly important branch of PFSP and an NP-hard problem.Additional mathematical descriptions are presented in the following subsections.determining the processing sequence of n jobs over m machines in PFSP to satisfy several objectives is widely applied in actual production, especially in one-piece mass production.The NIPFSP is a highly important branch of PFSP and an NP-hard problem.Additional mathematical descriptions are presented in the following subsections.

Notation i, j
normally utilized as loop variables (i.e., i represents the job number, and j represents the machine number) m machine number n job number Job { J1, J2, …, Jn}; represents the job set to be processed π scheduling solution that is the processing sequence of the job set {J1, J2, …, Jn} Ti,j represents the processing time of the i-th job processed on the j-th machine Tsi,j represents the starting time of the i-th processed job on the j-th machine; Given that all of the jobs are prepared to be processed at time zero, then represents the ending time of the i-th processed job on the j-th machine DifTi,j represents the minimum difference time between the π(i)-th processed job completion time of the j-th machine and (j + 1)-th machine di represents the due date of the i-th job ( ) represents the total tardiness of the schedule π

Mathematical Model
The NIPFSP has the same requirements (i.e., processing sequence) as the PFSP [2,5,14].In particular, n jobs must be processed on m machines with the same sequence.Each machine can only process one job at a time.No interruption occurs between the start and end of each job, and all of the jobs are prepared to be processed at time zero.However, a significant difference exists between the two problems.The NIPFSP has tight time constraints for each machine.In particular, each machine must process all of the jobs without any interruption from the start of processing the first job to the end of the last job.The total tardiness (TTd) of the NIPFSP for this feature can be calculated as follows:  represents the processing time of the i-th job processed on the j-th machine Ts i,j represents the starting time of the i-th processed job on the j-th machine; Given that all of the jobs are prepared to be processed at time zero, then Ts π( 1),1 = 0 Te i,j represents the ending time of the i-th processed job on the j-th machine DifT i,j represents the minimum difference time between the π(i)-th processed job completion time of the j-th machine and (j + 1)-th machine d i represents the due date of the i-th job TTd(π) represents the total tardiness of the schedule π

Mathematical Model
The NIPFSP has the same requirements (i.e., processing sequence) as the PFSP [2,5,14].In particular, n jobs must be processed on m machines with the same sequence.Each machine can only process one job at a time.No interruption occurs between the start and end of each job, and all of the jobs are prepared to be processed at time zero.However, a significant difference exists between the two problems.The NIPFSP has tight time constraints for each machine.In particular, each machine must process all of the jobs without any interruption from the start of processing the first job to the end of the last job.The total tardiness (TTd) of the NIPFSP for this feature can be calculated as follows: f irst machine : The aforementioned equations describe the processing sequence through the variables TS and TE.The relationship between TS and TE shows no idle tight time constraints.Therefore, several jobs must be processed with time delay when they are free from several machines.The variable DifT i,j is employed to ensure the exact delay time.The ending time of the last processed job at each machine can be calculated through the sum DifT i,j with T i,j .The ending time of each job at the last machine can then be attained by forward pass calculation.Finally, the total tardiness of the schedule π can be obtained.The objective of solving NIPFSP is to determine a suitable solution (i.e., π) with the total tardiness TTd(π) minimization.An example of the NIPFSP problem with five jobs and three machines is shown in Figure 2 in describing the manufacturing process.The five jobs are processed on the three machines with the same job process vector [4 1 5 3 2].The no idle constraint exists, whether the jobs are processed on the first machine or the following machines.The starting time and end time of each job are also posted in Figure 2.
Sustainability 2017, 9, 953 4 of 16 (1),1 The aforementioned equations describe the processing sequence through the variables TS and TE.The relationship between TS and TE shows no idle tight time constraints.Therefore, several jobs must be processed with time delay when they are free from several machines.The variable DifTi,j is employed to ensure the exact delay time.The ending time of the last processed job at each machine can be calculated through the sum DifTi,j with Ti,j.The ending time of each job at the last machine can then be attained by forward pass calculation.Finally, the total tardiness of the schedule π can be obtained.The objective of solving NIPFSP is to determine a suitable solution (i.e., π) with the total tardiness ( ) TTd  minimization.An example of the NIPFSP problem with five jobs and three machines is shown in Figure 2 in describing the manufacturing process.The five jobs are processed on the three machines with the same job process vector [4 1 5 3 2].The no idle constraint exists, whether the jobs are processed on the first machine or the following machines.The starting time and end time of each job are also posted in Figure 2.

HEDA_CS for NIPFSP
Significant applications of hybrid algorithms to solve scheduling problems with suitable performance have been explored in the past few decades [29][30][31][32].This section presents the novel Sustainability 2017, 9, 953 5 of 16 HEDA_CS for solving the NIPFSP with the total tardiness minimization.First, a discrete vector solution representation and initialization with the building probability model is introduced.Second, a hybrid strategy with CS, updating mechanism, and knowledge-based local search are described in detail.Finally, the flowchart of the HEDA_CS is shown to better understand the circulation mechanism.

Solution Representation
All of the solutions in the traditional estimation of distribution algorithms (EDAs) are designed for continuous optimization problems.As a typical discrete optimization problem, the NIPFSP has distinctive features in its solutions.Thus, a discrete vector decoding method is employed as the solution representation.In particular, each discrete vector demonstrates a solution for the NIPFSP.For example, the discrete vector π = {6, 4, 5, 1, 3, 2} is a scheduling solution of jobs processed in the order of 6, 4, 5, 1, 3, and 2 for the NIPFSP.Each decoded solution in this decoding method shows the unique scheduling result for the NIPFSP.A hybrid discrete algorithm can also be designed to solve the NIPFSP with more efficiency than continuous algorithms.Given that job i is represented by the element π(i), a solution can be decoded as a processing vector π as illustrated in Figure 3.

HEDA_CS for NIPFSP
Significant applications of hybrid algorithms to solve scheduling problems with suitable performance have been explored in the past few decades [29][30][31][32].This section presents the novel HEDA_CS for solving the NIPFSP with the total tardiness minimization.First, a discrete vector solution representation and initialization with the building probability model is introduced.Second, a hybrid strategy with CS, updating mechanism, and knowledge-based local search are described in detail.Finally, the flowchart of the HEDA_CS is shown to better understand the circulation mechanism.

Solution Representation
All of the solutions in the traditional estimation of distribution algorithms (EDAs) are designed for continuous optimization problems.As a typical discrete optimization problem, the NIPFSP has distinctive features in its solutions.Thus, a discrete vector decoding method is employed as the solution representation.In particular, each discrete vector demonstrates a solution for the NIPFSP.For example, the discrete vector π = {6, 4, 5, 1, 3, 2} is a scheduling solution of jobs processed in the order of 6, 4, 5, 1, 3, and 2 for the NIPFSP.Each decoded solution in this decoding method shows the unique scheduling result for the NIPFSP.A hybrid discrete algorithm can also be designed to solve the NIPFSP with more efficiency than continuous algorithms.Given that job i is represented by the element π(i), a solution can be decoded as a processing vector π as illustrated in Figure 3.

Job processing vector
Solution representation in HEDA_CS.

Initialization and Probability Model
Given that each individual is a solution representation as previously described, individuals are generated randomly in the initial population.The individuals in this method are randomly distributed in the entire solution space and have suitable diversity by uniform design.An individual is constructed by the NEH heuristic during the initialization to guarantee the initial population with a certain quality [11].
The probability model is the core brain of the EDA, which is commonly adopted to describe the distribution of the searching solution space.The probability matrix for describing the probability model is built based on several superior solution individuals during the iteration.New solution individuals are then obtained by sampling the probability matrix in the EDA.Therefore, improving the EDA performance is highly advantageous when the probability matrix contains the knowledge from the problems.This study focuses on the NIPFSP.The optimization objective is to determine the best scheduling solution with total tardiness minimization.Therefore, a modified probability matrix building method is proposed by taking advantage of the knowledge from the processing matrix T and designed solution representation.The designed probability matrix P is related to the job processing vector.
Element pij(l) in the probability matrix P represents the probability that job j appears before or in the i-th position at the l-th iteration.The value of pij(l) refers to the importance of a job when decoding a solution into a schedule.As the processing matrix T can be calculated to obtain the total processing time of each job, the original probability matrix P is designed considering the

Initialization and Probability Model
Given that each individual is a solution representation as previously described, individuals are generated randomly in the initial population.The individuals in this method are randomly distributed in the entire solution space and have suitable diversity by uniform design.An individual is constructed by the NEH heuristic during the initialization to guarantee the initial population with a certain quality [11].
The probability model is the core brain of the EDA, which is commonly adopted to describe the distribution of the searching solution space.The probability matrix for describing the probability model is built based on several superior solution individuals during the iteration.New solution individuals are then obtained by sampling the probability matrix in the EDA.Therefore, improving the EDA performance is highly advantageous when the probability matrix contains the knowledge from the problems.This study focuses on the NIPFSP.The optimization objective is to determine the best scheduling solution with total tardiness minimization.Therefore, a modified probability matrix building method is proposed by taking advantage of the knowledge from the processing matrix T and designed solution representation.The designed probability matrix P is related to the job processing vector.
Element p ij (l) in the probability matrix P represents the probability that job j appears before or in the i-th position at the l-th iteration.The value of p ij (l) refers to the importance of a job when decoding a solution into a schedule.As the processing matrix T can be calculated to obtain the total processing time of each job, the original probability matrix P is designed considering the knowledge of the longest total processing time job priority principle and roulette wheel selection.The original probability matrix P's concrete implementation process is illustrated as follows: T n,j The sum of each row and column in the original probability matrix P is evidently one, which fits the requirement of the job processing vector for the sum of each row.The column element still represents the probability that job j appears before or in the i-th position at the first iteration.The amount of the row elements only contains the knowledge from the processing matrix T.An example for better illustrating the process of generating the original probability matrix P is shown in Figure 4. (1) (1) (1 ; 2 ; , ) 1 The sum of each row and column in the original probability matrix P is evidently one, which fits the requirement of the job processing vector for the sum of each row.The column element still represents the probability that job j appears before or in the i-th position at the first iteration.The amount of the row elements only contains the knowledge from the processing matrix T.An example for better illustrating the process of generating the original probability matrix P is shown in Figure 4.

Lévy Flight Strategy in CS
The CS algorithm was developed by Yang and Deb [23] and is a new population-based metaheuristic algorithm inspired by nature.The CS algorithm simulates the cuckoo process to determine a suitable nest location to build an optimization process.The cuckoos have a special ability to lay their eggs in another host's nest, which was recently laid in by the host.The cuckoos lay their eggs among the eggs that were recently laid by the host, or even throw away the host's eggs to increase the successful probability of their own eggs hatching.However, the host can find extraneous eggs and throw them away, or even rebuild a nest in other places.This scenario is a key process to encourage growth and prepare the body for reproduction.Several eggs are hatched successfully, and new individuals arrive at a new suitable nest location through Lévy flights [33].Lévy flights are a global random walk; thus, this strategy provides the algorithm with the capability to search globally and locally.It then converges to the global optimality by exploring the search space efficiently even by the end of the iteration.
A discrete partially-mapped crossover (PMX) operation is introduced in the CS as a novel approach to conceal the cuckoo's eggs and constitute a well-behaved hybrid algorithm.PMX can be viewed as an extension of a two-point crossover.It is an efficient mechanism to mix the partially optimal information of an individual to obtain a better job processing vector.The hybrid algorithm then becomes suitable for the discrete NIPFSP.Nevertheless, another procedure can legalize the new individuals, which are caused by the simple two-point crossover.The node repetition in the PMX crossover can be avoided by utilizing a mapping function.Therefore, the PMX searches for many new better individuals without increasing the computational complexity.The entire procedure is shown in Figure 5.

Lévy Flight Strategy in CS
The CS algorithm was developed by Yang and Deb [23] and is a new population-based metaheuristic algorithm inspired by nature.The CS algorithm simulates the cuckoo process to determine a suitable nest location to build an optimization process.The cuckoos have a special ability to lay their eggs in another host's nest, which was recently laid in by the host.The cuckoos lay their eggs among the eggs that were recently laid by the host, or even throw away the host's eggs to increase the successful probability of their own eggs hatching.However, the host can find extraneous eggs and throw them away, or even rebuild a nest in other places.This scenario is a key process to encourage growth and prepare the body for reproduction.Several eggs are hatched successfully, and new individuals arrive at a new suitable nest location through Lévy flights [33].Lévy flights are a global random walk; thus, this strategy provides the algorithm with the capability to search globally and locally.It then converges to the global optimality by exploring the search space efficiently even by the end of the iteration.
A discrete partially-mapped crossover (PMX) operation is introduced in the CS as a novel approach to conceal the cuckoo's eggs and constitute a well-behaved hybrid algorithm.PMX can be viewed as an extension of a two-point crossover.It is an efficient mechanism to mix the partially optimal information of an individual to obtain a better job processing vector.The hybrid algorithm then becomes suitable for the discrete NIPFSP.Nevertheless, another procedure can legalize the new individuals, which are caused by the simple two-point crossover.The node repetition in the PMX crossover can be avoided by utilizing a mapping function.Therefore, the PMX searches for many new better individuals without increasing the computational complexity.The entire procedure is shown in Figure 5.

Updating Mechanism
Each optimization algorithm implements the population optimization through the iterative approach for solving the problem.Individuals with high fitness can provide search directions to attain the optimum solution in the NIPFSP.Thus, we obtain several better individuals in the population after a series of operations.We can then utilize the information in these individuals to update the probability model P. The details to update the probability model P are described in Equations ( 5) and ( 6).The probability model is directed toward the space that contains better solutions based on the top 10% of the best individuals.The search procedure must also track the potential searching region.( 1) ( 1) ( ) 1, if job appears before or in position 0, else where (0,1)

 
is the learning rate from the new better individuals, and Ii,j describes whether job j is located before position i.Considering the operation of the SP better individuals in Section 3.5, the value of  can be set to be slightly large.A sufficient strategy can maintain the exploration in the HEDA_CS in the aforementioned operation.As an operation to obtain better knowledge information in relatively better solutions, the  value must be set as extremely small.This scenario is also an updating mechanism to obtain new individuals during the iteration.

Knowledge-Based Local Search
A knowledge-based local search is designed as an operator during the iteration to improve the exploitation capability of HEDA_CS.A critical path for the PFSPs refers to a continuous job-path from the beginning to the end of the solution with no idle condition between any two jobs [34].Thus, this continuous job-path contains full processing knowledge of the NIPFSP with the total tardiness criterion minimization.Considering the variable neighborhood search [35], a suitable knowledge-based local search is designed in the HEDA_CS to solve the NIPFSP with the total tardiness criterion minimization.The insertion operator presents superior performance in the local search strategy.The knowledge-based local search mainly relies on such a partial subsequence insertion operator.In particular, the local search algorithm initially fetches a sequence from the job processing vector (i.e., an individual from SP better individuals).The length of the subsequence is at a maximum of n rounded down to the nearest whole unit.Several job numbers can then be obtained, which are the elements of the fetched sequence.Each fetched job is inserted to all of the possible positions of the remaining sequences.Only the best subsequence that has better fitness

Updating Mechanism
Each optimization algorithm implements the population optimization through the iterative approach for solving the problem.Individuals with high fitness can provide search directions to attain the optimum solution in the NIPFSP.Thus, we obtain several better individuals in the population after a series of operations.We can then utilize the information in these individuals to update the probability model P. The details to update the probability model P are described in Equations ( 5) and ( 6).The probability model is directed toward the space that contains better solutions based on the top 10% of the best individuals.The search procedure must also track the potential searching region.
where α ∈ (0, 1) is the learning rate from the new better individuals, and I i,j describes whether job j is located before position i.Considering the operation of the SP better individuals in Section 3.5, the value of α can be set to be slightly large.A sufficient strategy can maintain the exploration in the HEDA_CS in the aforementioned operation.As an operation to obtain better knowledge information in relatively better solutions, the α value must be set as extremely small.This scenario is also an updating mechanism to obtain new individuals during the iteration.

Knowledge-Based Local Search
A knowledge-based local search is designed as an operator during the iteration to improve the exploitation capability of HEDA_CS.A critical path for the PFSPs refers to a continuous job-path from the beginning to the end of the solution with no idle condition between any two jobs [34].Thus, this continuous job-path contains full processing knowledge of the NIPFSP with the total tardiness criterion minimization.Considering the variable neighborhood search [35], a suitable knowledge-based local search is designed in the HEDA_CS to solve the NIPFSP with the total tardiness criterion minimization.The insertion operator presents superior performance in the local search strategy.The knowledge-based local search mainly relies on such a partial subsequence insertion operator.In particular, the local search algorithm initially fetches a sequence from the job processing vector (i.e., an individual from SP better individuals).The length of the subsequence is at a maximum of √ n rounded down to the nearest whole unit.Several job numbers can then be obtained, which are the elements of the fetched sequence.Each fetched job is inserted to all of the possible positions of the remaining sequences.Only the best subsequence that has better fitness than the others is retained.Finally, a complete job sequence is obtained (i.e., an individual that has a minimal total tardiness criterion in such a job processing vector).
If the final individual is different from the original individual, then the aforementioned iteration steps are repeated until the individual is consistent.The better solutions are updated with high exploration quality in this knowledge-based local search.

Overall Implementation
The flowchart of the HEDA_CS for solving the NIPFSP with the total tardiness criterion minimization is illustrated in Figure 6 with the aforementioned designed procedure.At the beginning of the hybrid algorithm, a population is generated with a better solution provided by the NEH heuristic.The probability matrix P is initialized by obtaining the knowledge of the processing time matrix T, and then a CS operator is implemented.Given that several individuals are optimized in the population, all of the individuals are ranked with high fitness (i.e., low total tardiness).The SP better ranked individuals enhance the exploitation by the knowledge-based local search.The probability matrix P is then updated on the basis of the job processing sequence information represented by these better individuals.A new population can be generated by sampling the probability matrix P based on the updating mechanism of the EDA.Given that the stopping condition is not met in the HEDA_CS, the algorithm is iterated to the max generation (Maxgeneration).
The computational complexity at the HEDA_CS iteration can be roughly analyzed as follows.The CS strategy in the updating process requires computational complexity of O(PopSize/2) by the PMX operator on the population.Computational complexity of O(SP × n 2 ) is also observed.Each new generated individual in the sampling process is generated by the roulette strategy with computational complexity of O(n 2 ).The aforementioned analysis shows that the computation complexity of the proposed HEDA_CS is not excessively large.It can solve the NIPFFSP with the total tardiness criterion minimization within an acceptable range of calculations.
than the others is retained.Finally, a complete job sequence is obtained (i.e., an individual that has a minimal total tardiness criterion in such a job processing vector).If the final individual is different from the original individual, then the aforementioned iteration steps are repeated until the individual is consistent.The better solutions are updated with high exploration quality in this knowledge-based local search.

Overall Implementation
The flowchart of the HEDA_CS for solving the NIPFSP with the total tardiness criterion minimization is illustrated in Figure 6 with the aforementioned designed procedure.At the beginning of the hybrid algorithm, a population is generated with a better solution provided by the NEH heuristic.The probability matrix P is initialized by obtaining the knowledge of the processing time matrix T, and then a CS operator is implemented.Given that several individuals are optimized in the population, all of the individuals are ranked with high fitness (i.e., low total tardiness).The SP better ranked individuals enhance the exploitation by the knowledge-based local search.The probability matrix P is then updated on the basis of the job processing sequence information represented by these better individuals.A new population can be generated by sampling the probability matrix P based on the updating mechanism of the EDA.Given that the stopping condition is not met in the HEDA_CS, the algorithm is iterated to the max generation (Maxgeneration).
The computational complexity at the HEDA_CS iteration can be roughly analyzed as follows.The CS strategy in the updating process requires computational complexity of O(PopSize/2) by the PMX operator on the population.Computational complexity of O(SP × n 2 ) is also observed.Each new generated individual in the sampling process is generated by the roulette strategy with computational complexity of O(n 2 ).The aforementioned analysis shows that the computation complexity of the proposed HEDA_CS is not excessively large.It can solve the NIPFFSP with the total tardiness criterion minimization within an acceptable range of calculations.

Results and Analysis
Many tests are performed by utilizing the novel Ruiz benchmark instances in 2015 to investigate the performance of the proposed HEDA_CS [30].All of the data can be obtained at http://soa.iti.es(accessed on 16 April 2017).These novel benchmark instances are near the real production process and have practical significance.Each Ruiz instance has the same structure with Taillards' instances.Several case studies can be easily conducted utilizing the said algorithm.The due date of the i-th job is calculated as [36], where λ represents the tightness factor, T i,j is the processing time of i-th job processed on the j-th machine, and m is the total machine number.The tightness factor λ is set as 1, 2, and 3 to demonstrate the due date loose, medium, and tight, respectively.For the optimization objective of due date, the tightness factor can be seen as a different relaxation condition.Three values of the tightness factor λ represent three different requirements in the processing environment.
All of the experiment results are evaluated by average relative percentage deviation (ARPD) [37] to better evaluate the HEDA_CS performance.Thus, where G best is the total tardiness of the best solution obtained by all of the compared algorithms, and avg corresponds to the average value of the total tardiness of the solution obtained by a selected algorithm.So the lower value of the ARPD means that better solutions are achieved.The proposed HEDA_CS is coded in C++ (Visual Studio 2012) and run on a PC with an Intel(R) Core(TM) i7-2600 CPU 3.40 GHz and 2.85 GB of available main RAM.The computation results demonstrate the quality of the proposed HEDA_CS.The hybrid strategy shows its suitable performance in solving NIPFSP.The following subsections describe the exhaustive concept of the parameter setting, computational results, and discussion.

Parameter Setting
The three main parameters in the proposed HEDA_CS are as follows: Maxgeneration (iteration number), PopSize (population size), and SP (number of the superior ranked individuals in the population).All of the instances containing 60 jobs and 10 machines from the Ruiz benchmark instances are employed to adjust the said parameters.Another benchmark set is utilized to test the performance of the algorithm in the next section.An orthogonal experimental design method [38] is implemented to investigate the influence of these parameters on the HEDA_CS performance.
Given that the tightness factor λ has three different values to demonstrate the due date loose, medium, and tight, the three main parameters can be set differently under different tightness factors.The three levels of each main parameter are listed in Table 1.The HEDA_CS for each experiment environment run each instance 50 times independently (i.e., 50 × 10 × 9 = 4500 times).The entire orthogonal experiment is listed in Table 2.We can then obtain the trends of the tightness factor λ shown in Figure 7, which shows that the trends of the tightness factor λ is equal to its definition in different conditions.When the value of the tightness factor λ is large, the NIPFSP with the total tardiness criterion can be easily solved.Calculating the results in the orthogonal experiments yields the range and rank of the main parameters in the HEDA_CS as listed in Table 3.The trends of the main parameters in the HEDA_CS are shown in Figure 8. Maxgeneration has the most significant impact in the HEDA_CS, which is followed by PopSize.Given the experiment results, the recommended settings for the main parameters in the HEDA_CS are as follows: Maxgeneration = 1000, PopSize = 50, SP = 10.Such parameter settings can improve the efficiency of the HEDA_CS.

Results and Comparison of the Instances
Six sets of benchmark instances are selected from the Ruiz benchmark instances and are used to test the performance of the proposed HEDA_CS.Each set contains 10 different benchmarks with the same job and machine numbers.The factors and levels of these benchmarks are listed in Table 4.The benchmarks can generally describe the NIPFSP characteristics with the total tardiness criterion minimization.The range of the processing time distribution is U(1, 100).The selected benchmark job is in the range of [40,50,60,100].The number of process machines is in the range of [20,40,60].The proposed HEDA_CS is compared with several existing algorithms, such as GA, IEDA, and CS, by utilizing these instances.For each instance, all of the algorithms are run 20 times each.The computational results are summarized in Tables 5-7 with different values of the tightness factor λ. The HEDA_CS obtained nearly all of the total tardiness minimization of the instances.The convergence curves of the four algorithms that solved the instance VFR100_20_3_Gap are shown in Figure 9.The best Gantt charts obtained by the HEDA_CS for better illustrating the scheduling production process and providing the actual production reference for engineers are shown in Figures 10 and 11. Figure 10 illustrates the best scheduling solution of the instance VFR40_20_1_Gap under tight factor λ = 1.After enough iterations, the scheduling solution [36-7-17-6-20-5-29-1-13-26-19-38-28-24-9-25-3-31-30-14-4-21-2-35-11-34-15-8-32-37-10-39-12-27-18-33-40-16-22-23] is obtained by the proposed HEDA_CS.Following this scheduling solution, the 40 jobs in the instance VFR40_20_1_Gap can be processed with the lowest total tardiness among all of the solutions achieved during the iteration in HEDA_CS. Figure 11

Results and Comparison of the Instances
Six sets of benchmark instances are selected from the Ruiz benchmark instances and are used to test the performance of the proposed HEDA_CS.Each set contains 10 different benchmarks with the same job and machine numbers.The factors and levels of these benchmarks are listed in Table 4.The benchmarks can generally describe the NIPFSP characteristics with the total tardiness criterion minimization.The range of the processing time distribution is U(1, 100).The selected benchmark job is in the range of [40,50,60,100].The number of process machines is in the range of [20,40,60].The proposed HEDA_CS is compared with several existing algorithms, such as GA, IEDA, and CS, by utilizing these instances.For each instance, all of the algorithms are run 20 times each.The computational results are summarized in Tables 5-7 with different values of the tightness factor λ. The HEDA_CS obtained nearly all of the total tardiness minimization of the instances.The convergence curves of the four algorithms that solved the instance VFR100_20_3_Gap are shown in Figure 9.The best Gantt charts obtained by the HEDA_CS for better illustrating the scheduling production process and providing the actual production reference for engineers are shown in Figures 10 and 11. Figure 10 illustrates the best scheduling solution of the instance VFR40_20_1_Gap under tight factor λ = 1.After enough iterations, the scheduling solution [36-7-17-6-20-5-29-1-13-26-19-38-28-24-9-25-3-31-30-14-4-21-2-35-11-34-15-8-32-37-10-39-12-27-18-33-40 -16-22-23] is obtained by the proposed HEDA_CS.Following this scheduling solution, the 40 jobs in the instance VFR40_20_1_Gap can be processed with the lowest total tardiness among all of the solutions achieved during the iteration in HEDA_CS. Figure 11

Discussion of Experimental Results
The performance of the proposed HEDA_CS is tested and compared with several existing algorithms by calculating the group of instances in the aforementioned sections.The comparative results show that the HEDA_CS performs better in solving the benchmark instances than the other algorithms.Tables 5-7, show that the proposed HEDA_CS obtains nearly the best scheduling solution with better total tardiness than the compared algorithms at all tightness factor scenarios.The convergence curves show that the proposed HEDA_CS is more efficient than the other

Discussion of Experimental Results
The performance of the proposed HEDA_CS is tested and compared with several existing algorithms by calculating the group of instances in the aforementioned sections.The comparative results show that the HEDA_CS performs better in solving the benchmark instances than the other algorithms.Tables 5-7, show that the proposed HEDA_CS obtains nearly the best scheduling solution with better total tardiness than the compared algorithms at all tightness factor scenarios.The convergence curves show that the proposed HEDA_CS is more efficient than the other algorithms.The proposed HEDA_CS can still continue to optimize the total tardiness of the NIPFSP after much iteration.The discrete job vector coding method can help optimize the HEDA_CS operator.Hence, this hybrid strategy can better balance the HEDA_CS exploration and exploitation.The Gantt chart arrangement is regular.The optimization results can also improve the production process effectively.

Conclusions and Future Work
This study proposes an effective and efficient HEDA_CS to solve the NIPFSP with the total tardiness minimization.A time point-based problem formulation is demonstrated with an example.As a novel hybrid discrete algorithm, a discrete vector decoding method is utilized in the HEDA_CS.Several operators are designed in the EDA with the knowledge of the process time matrix.The Lévy flight strategy enhances the EDA exploration.The PMX operator is applied in the CS, and the knowledge-based local search ensures the HEDA_CS exploitation.The HEDA_CS effectiveness is shown by utilizing the novel Ruiz benchmark instances and comparing with other algorithms.The HEDA_CS is highly efficient in solving the NIPFSP with the total tardiness minimization by setting the suitable main parameters.The proposed HEDA_CS performs best by comparing with other algorithm in solving the NIPFSP.
Future research can be devoted to other hybrid algorithm strategies.The new strategies which better combine the key characteristics of the problem will be more effective.In addition, the proposed approaches can be considered to extend to other scheduling problems with different objectives, such as total completion time, earliness, and makespan.Further studies can also focus on solving this problem with multi-objectives, since existing research mainly considered the single-criterion problems.

Figure 2 .
Figure 2. Example of a no idle PFSP (NIPFSP) problem with five jobs and three machines.Figure 2. Example of a no idle PFSP (NIPFSP) problem with five jobs and three machines.

Figure 2 .
Figure 2. Example of a no idle PFSP (NIPFSP) problem with five jobs and three machines.Figure 2. Example of a no idle PFSP (NIPFSP) problem with five jobs and three machines.

Figure 4 .
Figure 4. Example of generating the original probability matrix P.

Figure 4 .
Figure 4. Example of generating the original probability matrix P.

Figure 5 .
Figure 5. Demonstration process example of a PMX crossover in HEDA_CS.

Figure 5 .
Figure 5. Demonstration process example of a PMX crossover in HEDA_CS.

Figure 6 .
Figure 6.Flowchart of the HEDA_CS for NIPFSP with the total tardiness criterion minimization.Figure 6. Flowchart of the HEDA_CS for NIPFSP with the total tardiness criterion minimization.

Figure 6 .
Figure 6.Flowchart of the HEDA_CS for NIPFSP with the total tardiness criterion minimization.Figure 6. Flowchart of the HEDA_CS for NIPFSP with the total tardiness criterion minimization.

Figure 7 .
Figure 7. Factor trends of the tightness factor λ in NIPFSP with the total tardiness criterion minimization.

Figure 7 .
Figure 7. Factor trends of the tightness factor λ in NIPFSP with the total tardiness criterion minimization.

Figure 8 .
Figure 8. Factor trends of main parameters in the HEDA_CS.

Figure 8 .
Figure 8. Factor trends of main parameters in the HEDA_CS.

Table 1 .
Factor levels of the three main parameters.

Table 3 .
Average value of the ARPD in different main parameter factor levels.

Table 1 .
Factor levels of the three main parameters.

Table 3 .
Average value of the ARPD in different main parameter factor levels.

Table 4 .
Factors and their levels for the selected Ruiz benchmarks.

Table 4 .
Factors and their levels for the selected Ruiz benchmarks.

Table 5 .
Comparison results on each set of benchmark problems (λ = 1).

Table 5 .
Comparison results on each set of benchmark problems (  = 1).

Table 6 .
Comparison results on each set of benchmark problems (  = 2).

Table 6 .
Comparison results on each set of benchmark problems (λ = 2).

Table 7 .
Comparison results on each set of benchmark problems (λ = 3).