E ﬀ ective Heuristic Algorithms Solving the Jobshop Scheduling Problem with Release Dates

: Manufacturing industry reﬂects a country’s productivity level and occupies an important share in the national economy of developed countries in the world. Jobshop scheduling (JSS) model originates from modern manufacturing, in which a number of tasks are executed individually on a series of processors following their preset processing routes. This study addresses a JSS problem with the criterion of minimizing total quadratic completion time (TQCT), where each task is available at its own release date. Constructive heuristic and meta-heuristic algorithms are introduced to handle di ﬀ erent scale instances as the problem is NP-hard. Given that the shortest-processing-time (SPT)-based heuristic and dense scheduling rule are e ﬀ ective for the TQCT criterion and the JSS problem, respectively, an innovative heuristic combining SPT and dense scheduling rule is put forward to provide feasible solutions for large-scale instances. A preemptive single-machine-based lower bound is designed to estimate the optimal schedule and reveal the performance of the heuristic. Di ﬀ erential evolution algorithm is a global search algorithm on the basis of population, which has the advantages of simple structure, strong robustness, fast convergence, and easy implementation. Therefore, a hybrid discrete di ﬀ erential evolution (HDDE) algorithm is presented to obtain near-optimal solutions for medium-scale instances, where multi-point insertion and a local search scheme enhance the quality of ﬁnal solutions. The superiority of the HDDE algorithm is highlighted by contrast experiments with population-based meta-heuristics, i.e., ant colony optimization (ACO), particle swarm optimization (PSO) and genetic algorithm (GA). Average gaps 45.62, 63.38 and 188.46 between HDDE with ACO, PSO and GA, respectively, are demonstrated by the numerical results with benchmark data, which reveals the domination of the proposed HDDE algorithm.


Introduction
Smart manufacturing has become a new competitive advantage in many countries. Therefore, governments have issued smart-manufacturing-related strategies, such as 'Made in China 2025' in China and 'Industrial 4.0' in Germany. One goal of smart manufacturing is to improve customer satisfaction by meeting the requirement of customization, which is a common customer demand pattern in the smart era. Jobshop is a manufacturing system that makes products with bespoke productive processes. For example, a gear is a key part in automatic equipment such as transmission. Precision gears processing time of job j on machine i; C i,j completion time of job j on machine i; C j maximum completion time of job j; a j,i,l = 1, if job j is processed on machine i before machine l; 0, otherwise. z i,j,k = 1, if job k precedes job j immediately on machine i; 0, otherwise. y i,j = 1, if the first operation of job j is processed on machine i; 0, otherwise. M a large positive number.
Therefore, the following MIP model is established.

SPT-DS Heuristic and Lower Bound
Given that the JSS problem is NP-completed, a heuristic algorithm can effectively achieve feasible solutions in a very short time for massive production setting, in which continuous production is more important than optimal solution. Townsend [28] reported the optimality of the shortest processing first (SPT) rule for the single-machine TQCT problem. An SPT-based heuristic is proposed to handle the JSS TQCT problem with release dates. A preemptive single-machine lower bound based on the optimality of the shortest remaining processing time (SRPT) rule for the preemptive single-machine TQCT problem is presented to evaluate the performance of the heuristic [29].

SPT-DS Heuristic
Matrix T stores the processing route of each job. According to the precedence of operations in matrix T, all current operations are stored in set A, where current operation means an immediate successor of the operation that has just been processed for a given job. Matrix S and C store the start time and completion time of all operations, respectively, where s i,j denotes the start time of job j on machine i, i.e., Q i,j , and c i,j denotes the completion time of job j on machine i, ∀j ∈ n, ∀i ∈ m, s i,j ∈ S, c i,j ∈ C. Each job j has a non-negative release date r j and processing time p i,j denotes the execution duration of Q i,j .
Step 1: Set the start time s i,j equals to release date r j and let c i,j = 0, ∀j ∈ n, ∀i ∈ m, Go to Step 2.
Step 2: Assign operation Q i,j with the earliest start time s i,j in set A. If there is more than one job with the earliest start time, then schedule the job with the shortest processing time. Go to Step 3.
Step 3: Update the completion time c i,j = s i,j + p i,j . Delete operation Q i,j in set A. Check the elements in matrix T to determine the machine on which the next execution of job j is processed. If machine k is the target, then add operation Q k,j in set A. Go to Step 4. Step 4: Update the start time of the remaining operations of job j, let s k,j = max (s k,j , c i,j ), 1 ≤ k ≤ m, and the start processing time of machine i, s i,h = max (s i,h , c i,j ), 1 ≤ h ≤ n. Go to Step 5.
Step 5: If set A Ø, turn to Step 2. Otherwise, stop the algorithm and calculate the objective value. A numerical example is provided below to better explain the SPT-DS heuristic.

Example 1.
A JSS problem with three jobs {J 1 , J 2 , J 3 } and three machines {M 1 , M 2 , M 3 } is involved. The input data including processing times, release dates and processing routes are presented in Table 2.
The processing times and the processing routes of three jobs in Table 2 are restored in Matrices P and T. Matrices S and C record the start and completion times of each operation Q i,j , 1 ≤ i ≤ 3, 1 ≤ j ≤ 3. In initial status, the completion times of all operations are set to zero, and the start times of jobs on all machines are equal to their release dates. The first raw in matrix T is [2,3,1] which means the first execution of job j, 1 ≤ j ≤ 3, can be processed on machine M 2 , M 3 or M 1 . Therefore, the data stored in set A is [Q 2,1 , Q 3,2 , Q 1,3 ]. The elements in each matrix are shown as follows.
Next, repeat the above steps to select an operation from set A to schedule until set A is empty. The final operation schedules on machines M1, M2, and M3 are {J3, J2, J1,}, {J1, J2, J3} and {J2, J1, J3}, respectively. The Gantt chart of the SPT-DS schedule is provided in Figure 1. The objective value of the schedule is 12 2 + 10 2 + 15 2 = 469.

Well-Designed Lower Bound
An effective lower bound is proposed for the JSS TQCT problem to estimate the optimal value for large-scale instances. The basic idea of designing the lower bound is to relax the constraints of the relationship between machines and non-preemption for each job. Letting , denote that operation , i k O is scheduled according to the SRPT rule, lower bound is shown as follows.
Given the NP-hardness of the JSS TQCT problem, i.e., it cannot be solved in polynomial time for large-scale instances, an effective lower bound is usually served as a substitute of the optimal schedule to evaluate approximation algorithms. Generally, the lower bound is the optimal solution of a relaxation version of the original problem. An m-machine JSS TQCT problem is reduced to m preemptive versions of the single-machine scheduling (SMS) problem to minimize the TQCT criterion with release dates. Bai [29] proved that the SRPT rule is optimal for the preemptive SMS TQCT problem. Therefore, each of these SRPT schedules is a lower bound of the JSS TQCT problem. As a larger lower bound can provide better performance for minimization criteria, the dominative one among the m lower bounds is selected as the final lower bound for the JSS TQCT problem.
To better explain the lower bound, a numerical example is provided as follows.

Well-Designed Lower Bound
An effective lower bound is proposed for the JSS TQCT problem to estimate the optimal value for large-scale instances. The basic idea of designing the lower bound is to relax the constraints of the relationship between machines and non-preemption for each job. Letting p SRPT i,k denote that operation O i,k is scheduled according to the SRPT rule, lower bound Z LB is shown as follows. where Given the NP-hardness of the JSS TQCT problem, i.e., it cannot be solved in polynomial time for large-scale instances, an effective lower bound is usually served as a substitute of the optimal schedule to evaluate approximation algorithms. Generally, the lower bound is the optimal solution of a relaxation version of the original problem. An m-machine JSS TQCT problem is reduced to m preemptive versions of the single-machine scheduling (SMS) problem to minimize the TQCT criterion with release dates. Bai [29] proved that the SRPT rule is optimal for the preemptive SMS TQCT problem. Therefore, each of these SRPT schedules is a lower bound of the JSS TQCT problem. As a larger lower bound can provide better performance for minimization criteria, the dominative one among the m lower bounds is selected as the final lower bound for the JSS TQCT problem.
To better explain the lower bound, a numerical example is provided as follows.
Example 2. The input data are similar to those in Example 1. The lower bound schedule with the SRPT rule on the three machines is shown in Figure 2. The scheduling on machine M 2 is explained in detail for better understanding. At time t = 0, only job J 2 is available. At time t = 1, job J 1 is released and the remaining time of job J 2 is 5 > p 2,1 = 4. On the basis of the SRPT rule, job J 2 is preempted by job J 1 . At time t = 2, similarly, job J 1 is preempted by job J 3 . Given that no extra job is released after the completion of job J 3 , the remaining parts of jobs J 1 and J 2 are scheduled with the SRPT rule. The full schedule for the lower bound is presented in Figure 2. The associated lower bound value is Z LB = max 2 2 + 4 2 + 6 2 , 11 2 + 3 2 + 6 2 , 9 2 + 4 2 + 14 2 = 293. jobs J1 and J2 are scheduled with the SRPT rule. The full schedule for the lower bound is presented in Figure 2. The associated lower bound value is Z = max 2 + 4 + 6 , 11 + 3 + 6 , 9 + 4 + 14 = 293.

Effective HDDE Algorithm
The SPT-DS heuristic algorithm can quickly output a feasible solution for large-scale instances. However, the quality of the SPT-DS schedule weakens for medium-scale instances. For example, the mean gap is approximately 185% for 20 random instances with the combination m × n = 3 × 30. Relatively, meta-heuristic is a type of higher level heuristic that controls the entire search process. Thus, the global optimal solutions can be achieved systematically and efficiently. The differential evolution (DE) algorithm is a population-based evolutionary algorithm for global optimization in a continuous search space. The DE algorithm improves candidates by executing mutation and crossover operations, and renews the population through greedy one-to-one selection. Compared with the classical evolutionary-based metaheuristics, the superiority of the DE algorithm involves easy implementation, simple structure, robust, and fast convergence. The traditional DE algorithm was originally used to solve the global optimization problem in the continuous search space, in which individuals were encoded by floating point numbers that are invalid for discrete variables. This section transforms individuals into operation-based sequences and proposes the HDDE algorithm to solve the JSS problem.

Encoding and Decoding
Given the diversity of process routes for each job in a JSS model, a job-number-based (JNB) representation is introduced to encode and decode between individuals and feasible schedules in the DDE algorithm. The encoding scheme linearizes the operations in a feasible schedule according to their starting times and represents an individual by the corresponding job number. The decoding scheme restores a JNB individual to a feasible schedule, where each operation is assigned in turn to the available machine following the given machine sequence. Matrices × (operation sequences) and × (processing times of operations) are presented as follows.
A numerical example is provided below to further describe the encoding and decoding processes.

Effective HDDE Algorithm
The SPT-DS heuristic algorithm can quickly output a feasible solution for large-scale instances. However, the quality of the SPT-DS schedule weakens for medium-scale instances. For example, the mean gap is approximately 185% for 20 random instances with the combination m × n = 3 × 30. Relatively, meta-heuristic is a type of higher level heuristic that controls the entire search process. Thus, the global optimal solutions can be achieved systematically and efficiently. The differential evolution (DE) algorithm is a population-based evolutionary algorithm for global optimization in a continuous search space. The DE algorithm improves candidates by executing mutation and crossover operations, and renews the population through greedy one-to-one selection. Compared with the classical evolutionary-based meta-heuristics, the superiority of the DE algorithm involves easy implementation, simple structure, robust, and fast convergence. The traditional DE algorithm was originally used to solve the global optimization problem in the continuous search space, in which individuals were encoded by floating point numbers that are invalid for discrete variables. This section transforms individuals into operation-based sequences and proposes the HDDE algorithm to solve the JSS problem.

Encoding and Decoding
Given the diversity of process routes for each job in a JSS model, a job-number-based (JNB) representation is introduced to encode and decode between individuals and feasible schedules in the DDE algorithm. The encoding scheme linearizes the operations in a feasible schedule according to their starting times and represents an individual by the corresponding job number. The decoding scheme restores a JNB individual to a feasible schedule, where each operation is assigned in turn to the available machine following the given machine sequence. Matrices T m×n (operation sequences) and P m×n (processing times of operations) are presented as follows.
A numerical example is provided below to further describe the encoding and decoding processes.

Example 3.
A JSS model with three jobs {J 1 , J 2 , J 3 } and three machines {M 1 , M 2 , M 3 } is provided to explain the procedure in detail. The input data are similar to those in Example 1. The data in Matrices T 3×3 andP 3×3 are presented as follows.

Initialization
The main task of initialization is to generate the initial population and set the parameters. The initial population is randomly generated by swapping adjacent genes in a chromosome. Orthogonal tests are designed to determine four parameters, including population size Λ, the maximum number of iterations τ max , mutant factor Z, and crossover factor Y, thereby enhancing the performance of the algorithm. The pseudo-code for initialization is provided in Procedure 1. Input: n,Λ,m //n is the number of jobs, m is the number of machines, Λ is the size of the population

Mutation and Crossover
The mutation operator is executed on all individuals in the current population to generate the } are randomly selected from the (k−1)th population to implement the mutation operation with current optimal individual X k−1 best (the individual with the optimal objective value in the (k−1)th population). The mutation operator is expressed as follows.
where rand (•) is a random number assigned to each element of the difference individual and x = 1, 2.
Operator ⊕ in Equality (11) is executed as follows.
where mod (•) is the modular operator, which guarantees that each vector in a mutation individual represents a valid job number. The pseudo-code for the mutation operator is provided in Procedure 2.

Procedure 2: Mutation operator
/Operation of operator ⊗ */ 9 X k−1 αx, j ←An individual randomly selected from the population; 10 X k−1 βx, j ←Another non-repeating individual randomly selected from the population; 11 for (j from 1 to m × n) do //n is the number of jobs and m is the number of machines 12 The crossover operator selectively inserts a mutation individual into a target individual X k−1 h to generate a trial individual U k h . A series of random numbers is assigned to the vectors in a mutant individual V k h . For each element in the mutation individual, if the random number rand (•) < Y, then the element is retained; otherwise, it is removed, where Y ∈ (0, 1) is a crossover factor. Tri-point insertion is executed to obtain a temporary individual, where the remaining elements of a mutation individual are divided into three parts and inserted in three randomly selected positions on target individual X k−1 h . The trial individual U k h is achieved by deleting the extra elements (repeat more than m times) in the temporary individual from left to right. The pseudo-code for crossover operator is presented in Procedure 3.

Procedure 3: Crossover operator
for (j from 0 to m × n − 1) do //n is the number of jobs and m is the number of machines 5 r ← random (0,1); Generate three random insertion points of U k h ; 12 The divided three parts of V k h is inserted into the three random insertion points of U k h in sequence; 13 A numerical example is provided below to further describe the mutation and crossover operations. The procedure of mutation operation is illustrated in Tables 3 and 4, while that of the crossover operation is illustrated in Table 5 and Figure 3.

Hill-Climbing-Based Improvement Strategy
A hill-climbing-based improvement strategy is introduced to enhance the quality of the final solution, aiming to balance diversification and intensification. Hill-climbing is a local search algorithm, which starts with an arbitrary solution of a given problem, then attempts to seek a better solution by making an incremental change to the solution. The change procedure executes continuously until no further improvements can be found.
However, the hill-climbing consumes considerable running time to obtain the local optimal solution in practical settings. For saving computing resource and improving algorithm efficiency, the improvement strategy is presented as follows. If the local optimal solution is found within the specified iterations, the solution is recorded; otherwise, the current best solution is recorded. The dominative parent individual is denoted as . Ten iterations are conducted for the improvement strategy. Each iteration generates 20 neighborhood solutions by exchanging two elements in the parent individual. The optimal one among the 20 neighborhood solutions is compared with the current optimal individual (COI). If the new optimal solution is dominated, then the COI is updated by it and the next iteration is executed. Otherwise, the improvement strategy will continue to be dominated by the original COI until it reaches the maximum number of iterations. The improvement strategy is performed with probability θ to save computing time. The pseudo-codes for the improvement strategy are shown in Procedure 4.

Hill-Climbing-Based Improvement Strategy
A hill-climbing-based improvement strategy is introduced to enhance the quality of the final solution, aiming to balance diversification and intensification. Hill-climbing is a local search algorithm, which starts with an arbitrary solution of a given problem, then attempts to seek a better solution by making an incremental change to the solution. The change procedure executes continuously until no further improvements can be found.
However, the hill-climbing consumes considerable running time to obtain the local optimal solution in practical settings. For saving computing resource and improving algorithm efficiency, the improvement strategy is presented as follows. If the local optimal solution is found within the specified iterations, the solution is recorded; otherwise, the current best solution is recorded. The dominative parent individual is denoted as U k h . Ten iterations are conducted for the improvement strategy. Each iteration generates 20 neighborhood solutions by exchanging two elements in the parent individual. The optimal one among the 20 neighborhood solutions is compared with the current optimal individual (COI). If the new optimal solution is dominated, then the COI is updated by it and the next iteration is executed. Otherwise, the improvement strategy will continue to be dominated by the original COI until it reaches the maximum number of iterations. The improvement strategy is performed with probability θ to save computing time. The pseudo-codes for the improvement strategy are shown in Procedure 4.

Selection
Population updating simulates the law of natural selection: survival of the fittest. Trial individual U k h is decoded as a feasible schedule. If the objective value of U k h dominates that of the target individual X k−1 h , then the individual in the current population will be updated; otherwise, the target individual will be retrained. Executing the mutation, crossover and selection operations for all target individuals generate a new population. The iteration procedure is repeated until the termination condition is satisfied.

Framework of the HDDE Algorithm
With the previous procedures combined, the entire framework of the HDDE algorithm is provided in Algorithm 1.

Numerical Simulation Experiment
Numerous simulation experiments were executed to evaluate the effectiveness of the HDDE algorithm and SPT-DS heuristic for medium-and large-scale instances, respectively. The proposed algorithms were implemented in C++ language, and simulations were run on a Dell computer with an Intel Core i5-8400 (2.8 GHz × 6) CPU and 8.00 GB RAM. The input data were generated because no standard benchmark was available for the JSP problem with release dates. Release date r j was randomly generated from a uniform distribution U[0, 3n]. Without loss of generality, at least one job was released at time zero. Processing time p i,j was randomly generated from a uniform distribution U [1,10]. For operation O i,j , processing time p i,j = 0 indicates that job j does not pass through machine i. The processing routes of each job were randomly generated, and the job might not pass through all machines. For each machine-job combination, 10 random trials were conducted and the average values are reported in the following tables. In Table 6, the GAPs are stable for the constant machine number. For the five-machine example, the GAP values float between 1.0664 and 1.2134 as the problem scale increases from 100 to 500, where the fluctuation range is about 7.5%. Furthermore, the growth trend of GAP values indicates that the SPT-DS heuristic deteriorates with the increase in m for a fixed problem scale. For the 300-job example, the GAP value increases from 0.7349 to 1.4221 as the machine scale increases from 3 to 8. This phenomenon might be attributed to the large machine scale increasing the idle time between adjacent operations in a feasible schedule, thereby weakening the performance of the heuristic.

Improvement of the HDDE Algorithm
To highlight the improvement schemes of the tri-insertion and local search, comparative experiments were executed between the standard DDE (SDDE) and the HDDE algorithms. The machine-job combinations m = {3, 5, 8} with n = {50, 100, 150} were tested under the input settings presented at the beginning of this section. The performance of a meta-heuristic mainly depends on the parameter settings. Thus, a series of orthogonal experiments was conducted to determine the parameters of the algorithm (Appendix A), such as mutation factor Z, crossover factor Y, population size Λ, and local search probability θ (Table A2). The parameters are set as Z = 0.2, Y = 0.1, Λ = 200, τ max = 300, and θ = 0.8.
Relative difference percentage RDP = Z H −Z * Z * × 100% was used for evaluation of the effectiveness of the HDDE algorithm, where Z H is the final objective value obtained by a meta-heuristic for each test, and Z * is the minimum value among all Z H values. Table 7 shows the average RDP (ARDP), minimum RDP (MinRDP), maximum RDP (MaxRDP), and standard deviation (SD) for SDDE and HDDE. ARDP is a performance measure that evaluates an algorithm numerically. The ARDPs of HDDE evidently dominate those of SDDE. For example, the mean values of ARDPs obtained by HDDE and SDDE are 0.6011 and 10.1811, respectively. SD aims to reveal the stability performance of an algorithm. The SD values of HDDE and SDDE are stable in [0, 1.21] and [2.91, 4.34], respectively, except for the maximum values, indicating the robustness of the former relative to the latter. In the minimum objective values of the 90 numerical tests, about 94% (85 out of 90) of the data were obtained by HDDE and about 6% of the data (5 out of 90) were obtained by SDDE. These final experimental outcomes show that the improvement schemes considerably enhance the performance of the HDDE algorithm.

Comparison between HDDE and ACO
Ant Colony Optimization (ACO) algorithm is selected as a contrast objective to indicate the superiority of the HDDE algorithm. The machine-job combinations and input parameters are identical with those in Section 6.2. The RDP measure is used in this section. The parameters used in ACO algorithm are presented as follows. P ij k is the probability of selecting operation j when ant k climbs to position i (i.e., the moment when the current operation is selected).
where Ω k is the tabu list for ant k, which stores operations that have been processed. τ i,j is the pheromone concentration of operation j. ŋ i,j is the current heuristic factor of operation j at present.
where d ij is the processing time of operation j. Parameters α and β indicate the relative importance of pheromones and heuristics.
where Q is a constant, and C k is the objective values of the operations in which ant k has been processed.
∆τ ij is the total increment of the pheromone of operation j.
τ ij is the current pheromone concentration of operation j, ∆τ ij is the pheromone concentration of operation j in the previous iteration stage, and ρ is the evaporation coefficient.
After previous preparation, a basic framework of the ACO algorithm is presented as follows.
Step 1. Set parameters for the algorithm, such as the ant colony size Λ, the maximum number of iterations t max , constant Q, evaporation coefficient ρ, and parameters α and β. The initial pheromone concentration of each operation is set to a certain constant. Set the current number of iterations t = 0.
Step 2. Generate the initial colony. With Equation (14), each ant uses roulette method to continuously select the operations that need to be processed. If all operations are included in the sequence of each ant, the initial ant colony is formed. Use Equations (16)- (18) to update the pheromone of each operation.
Step 4. The current ant uses the roulette method to select an operation with Equation (14). Once an operation is selected, the pheromone increment of the operation is obtained by Equation (16).
Step 5. If h < Λ, set h = h + 1, and then go to Step 4; otherwise go to Step 6.
Step 6. If t < t max , set t = t + 1. The pheromone increment of each operation is obtained by Equation (16). Update the pheromone of each operation in current iteration with Equation (17), and then go to step 3; otherwise go to step 7.
Step 7. Terminate procedure if the termination condition is satisfied. The optimal individual of the current ant colony is output as the final solution.
Similarly, the parameters of ant colony algorithm are set by orthogonal experiment, as shown in Appendix B. The parameters are set as α = 1.0, β = 0.5, ρ = 0.9, Q = 0.5, Λ = 200 and t max = 300. Table 8 shows the average RDP (ARDP), minimum RDP (MinRDP), maximum RDP (MaxRDP), and standard deviation (SD) for ACO and HDDE. The ARDPs of HDDE evidently dominate those of ACO. For example, the mean values of ARDPs obtained by HDDE and ACO are 0.3789 and 22.4422, respectively. The SD values of HDDE and ACO are stable in [0, 1.78] and in [3.12, 9.69], respectively, except for the maximum values, indicating the robustness of the former relative to the latter. In the minimum objective values of the 90 numerical tests, about 92% (83 out of 90) of the data were obtained by HDDE and about 8% of the data (7 out of 90) were obtained by ACO. These numerical results reveal that the HDDE algorithm completely dominates the ACO algorithm.

. Comparison between HDDE and PSO
To highlight the dominance of the HDDE algorithm, a series of comparison experiments are conducted between the HDDE and PSO algorithms. The concrete procedure of the PSO algorithm is referred to [30,31]. The parameters of the HDDE algorithm are identical with those in Section 6.3.1. The parameters of the PSO algorithm are set by orthogonal experiment (shown in Appendix C), i.e., population size Λ = 200, maximum number of iterations τ max = 300, weight w = 5, maximum speed V max = 4, maximum position S max = 4, cognitive coefficient c 1 = 3, and social coefficient c 2 = 3.
The data in Table 9 are the values of ARDP, MinRDP, MaxRDP, and SD for PSO and HDDE. The results of HDDE are all 0 for the 90 trials. This phenomenon indicates that the HDDE algorithm completely dominates the PSO algorithm.

Comparison between HDDE and GA
This section conducted comparison experiments on the HDDE and GA algorithms. The experimental parameters of the GA algorithm are set as follows: mutation probability Z 1 = 0.5, crossover probability Y 1 = 0.5, population size Λ = 200, and maximum number of iterations τ max = 300. The detailed parameter setting process is presented in Appendix D.
The data of ARDP, MinRDP, MaxRDP, and SD for GA and HDDE are shown in Table 10. It is obvious that HDDE's ARDPs is superior to GA's ARDPs. For example, the mean values of ARDPs obtained by HDDE and GA are 0.158 and 23.769, respectively. The SD values of GA are stable in [5.70, 12.91], except for the maximum values, and the SD value of HDDE is only equal to 3.04 at one scale (3 × 150), but equal to zero at all other scales, indicating the robustness of the latter relative to the former. In the minimum objective values of the 90 numerical tests, about 97.8% (88 out of 90) of the data were obtained by HDDE and about 2.2% of the data (2 out of 90) were obtained by GA. These numerical results reveal that the HDDE algorithm completely dominates the GA algorithm.

Comparison under JSS Problem Benchmarks
Currently, no standard benchmark is provided for the JSS problem with release dates. Therefore, the basic benchmarks proposed by Taillard, E. [32] is combined with generated data of release dates from a uniform distribution [0, 3n] to verify the performance of the proposed algorithm.
To control the meta-heuristics running in an appropriate time, the data of processing times with machine-job combinations m × n = 15 × 50, 20 × 50, and 20 × 100 were selected from the benchmarks. The data in Table 11 were the final objective values of the HDDE, ACO, PSO, and GA algorithms, where gap 1 = (Z ACO − Z HDDE )/Z HDDE × 100%, gap 2 = (Z PSO − Z HDDE )/Z HDDE × 100%, gap 3 = (Z GA − Z HDDE )/Z HDDE × 100%, Z HDDE , Z ACO , Z PSO , and Z GA are the final objective values obtained by HDDE, ACO, PSO, and GA, respectively. It is obvious that the HDDE algorithm is superior to other population-based meta-heuristics. However, the tested meta-heuristics consumed different CPU times. For an m × n = 20 × 100 instance, the CPU time expended by the HDDE, ACO, PSO, and GA algorithms are 17,890, 22,055, 837.853, and 238.448 s, respectively. The results reveals that better performance of an algorithm will cost more CPU time. As a trade-off between solution quality and running time, the HDDE algorithm dominates the other meta-heuristics.

Conclusions
This study investigates the JSS problem to minimize the total quadratic completion time where each job is available at its release date. No polynomial algorithm can solve the problem because of its NP-hardness. For large-sized instances, the SPT-DS heuristic algorithm is proposed to achieve approximate solutions in a short computing time. For medium-sized instances, the HDDE algorithm is provided to achieve high-quality solutions, where well-designed tri-insertion crossover and lower search schemes significantly enhance the performance of the meta-heuristic. A series of random experiments demonstrates the effectiveness of the proposed algorithms.
Future researches mainly focus on two aspects. On the one hand, the acceleration scheme will be proposed to save running time for the HDDE algorithm. On the other hand, the model will be generalized to a multi-objective version that is more popular in a production scheduling environment. A non-dominated sorting-based meta-heuristic will be presented to obtain Pareto optimal solutions. A high-efficiency improvement scheme, such as variable neighborhood search, will be combined with the meta-heuristic for JSS problems.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
Orthogonal experiments are designed to determine the appropriate parameters of the DDE algorithm, including population size Λ, mutation factor Z, crossover factor Y, maximum iteration τ max and local search probability θ, because the performance of the DDE algorithm is mainly dependent on parameter setting. The DDE algorithm parameters table is shown in Table A1.
To save computing time, the test scale is fixed at m = 3 and n = 60. The mean improvement percent MIP = Z INI −Z FIN Z FIN × 100% is used to measure the effect of the parameters, where Z INI is the minimum objective value among the initial population, and Z FIN is the final objective value of the output solution at the termination of algorithm. The orthogonal experiment results for the parameters of the HDDE algorithm are shown in Table A2. The main effect graph of the mean value for the orthogonal experiment results is shown in Figure A1. Referring to the peak point of each parameter (2-1-1-5-2) in Figure A1, it is obvious that the best parameter combinations are Z = 0.2, Y = 0.1, Λ = 200, τ max = 300 and θ = 0.8.
X.-y.W. and P.Z.; conceptualization-problem modeling, S.-R.C.; writing-original draft preparation, M.Z.; writingreview and editing, C.-C.W. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
Orthogonal experiments are designed to determine the appropriate parameters of the DDE algorithm, including population size Λ, mutation factor Z, crossover factor Y, maximum iteration τmax and local search probability θ, because the performance of the DDE algorithm is mainly dependent on parameter setting. The DDE algorithm parameters table is shown in Table A1.
To save computing time, the test scale is fixed at m = 3 and n = 60. The mean improvement percent  Table A2. The main effect graph of the mean value for the orthogonal experiment results is shown in Figure A1. Referring to the peak point of each parameter (2-1-1-5-2) in Figure A1, it is obvious that the best parameter combinations are Z = 0.2, Y = 0.1, Λ = 200, τmax = 300 and θ = 0.8. Figure A1. The main effect graph of HDDE. Figure A1. The main effect graph of HDDE.

Appendix B
Similarly, orthogonal experiments are designed to determine the appropriate parameters of the ACO algorithm, including ant colony size Λ, pheromone quantity Q, decay coefficient ρ, maximum iteration t max , accumulated information factor α, and inspiration information factor β. The ACO algorithm parameters level table is shown in Table A3.
To save computing time, the test scale is fixed as m = 3 and n = 60. MIP is used to measure the effect of the parameters. The orthogonal experiment results for the parameters of the ACO algorithm are shown in Table A4. The main effect graph of mean values for the orthogonal experiment results is shown in Figure A2. Referring to the peak point of each parameter (5-4-1-5-1-5) in Figure A2, it is obvious that the best parameter combinations are α = 1.0, β = 0.5, ρ = 0.9, Q = 0.5, Λ = 200, and t max = 300.  Figure A3. The main effect graph of PSO.

Appendix D
The parameters of GA mainly include mutation probability Z1, crossover probability Y1, population size Λ, and maximum iteration τmax. Similarly, an orthogonal experiment with four factors and five levels was performed to determine parameters. The parameters level table of GA is shown in Table A7. The results for the parameters of GA are shown in Table A8. The main effect graph of mean values for the orthogonal experiment results is shown in figure D1. Referring to the peak point of each parameter 5-5-1-5) in Figure A4, it is obvious that the best parameter combinations are Z1 = 0.5, Y1 = 0.5, Λ = 200, τmax = 300.

Appendix D
The parameters of GA mainly include mutation probability Z 1 , crossover probability Y 1 , population size Λ, and maximum iteration τ max . Similarly, an orthogonal experiment with four factors and five levels was performed to determine parameters. The parameters level table of GA is shown in Table A7. The results for the parameters of GA are shown in Table A8. The main effect graph of mean values for the orthogonal experiment results is shown in Table A4. Referring to the peak point of each parameter 5-5-1-5) in Figure A4, it is obvious that the best parameter combinations are Z 1 = 0.5, Y 1 = 0.5, Λ = 200, τ max = 300.
The parameters of GA mainly include mutation probability Z1, crossover probability Y1, population size Λ, and maximum iteration τmax. Similarly, an orthogonal experiment with four factors and five levels was performed to determine parameters. The parameters level table of GA is shown in Table A7. The results for the parameters of GA are shown in Table A8. The main effect graph of mean values for the orthogonal experiment results is shown in figure D1. Referring to the peak point of each parameter 5-5-1-5) in Figure A4, it is obvious that the best parameter combinations are Z1 = 0.5, Y1 = 0.5, Λ = 200, τmax = 300.  0.5 0.5 600 300 Figure A4. The main effect graph of GA.