Meta-heuristic algorithms for parallel identical machines scheduling problem with weighted late work criterion and common due date

To our knowledge, this paper investigates the first application of meta-heuristic algorithms to tackle the parallel machines scheduling problem with weighted late work criterion and common due date (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P|{{d}_{j}}=d|{{Y}_{w}}$$\end{document}P|dj=d|Yw). Late work criterion is one of the performance measures of scheduling problems which considers the length of late parts of particular jobs when evaluating the quality of scheduling. Since this problem is known to be NP-hard, three meta-heuristic algorithms, namely ant colony system, genetic algorithm, and simulated annealing are designed and implemented, respectively. We also propose a novel algorithm named LDF (largest density first) which is improved from LPT (longest processing time first). The computational experiments compared these meta-heuristic algorithms with LDF, LPT and LS (list scheduling), and the experimental results show that SA performs the best in most cases. However, LDF is better than SA in some conditions, moreover, the running time of LDF is much shorter than SA.

phrase "late work" was first proposed by Potts and Van Wassenhove (1992), who applied it to the single machine cases. Other researchers such as Hochbaum and Shamir (1990), Hariri et al. (1995), Kovalyov and Kubiak (1999), Kethley and Alidaee (2002) also studied the single machine cases. Then, the late work performance measures were applied to the shop scheduling problems (Błażewicz et al. 2000). Similarly, Blazewicz et al. (2004) focused on the open shop scheduling problem and Blazewicz et al. (2005b) studied the two-machine flow-shop problem. A number of scholars used meta-heuristic approaches to solve the scheduling problem with late work criterion in dedicated machine case. In Blazewicz et al. (2005a) and (2008), three meta-heuristic approaches (tabu search, simulated annealing and variable neighborhood search) were applied to the F 2|d j = d|Y w problem, and genetic algorithm (GA) was used to solve the problem of F |r j |Y (Sterna et al. 2007).
However, meta-heuristic algorithms have not been used in the context of parallel machines case. We propose, for the first time, three meta-heuristic algorithms for the parallel identical machines scheduling problem with weighted late work criterion and common due date which is denoted by P|d j = d|Y w . Three meta-heuristic algorithms include ant colony system algorithm (ACS), genetic algorithm (GA), and simulated annealing (SA).
Two typical scheduling algorithms were also implemented: longest processing time first (LPT) and list scheduling (LS). The LPT should sort all jobs based on their processing time first, and assign each job to a machine with minimum load. While in LS, a job will be directly assigned to the machine with minimum load. Additionally, we proposed a novel algorithm named LDF (largest density first) which is an improved algorithm of LPT. In the simulation experiments, six algorithms mentioned above were compared with each other under different values of parameters.
The rest of the paper is organized as follows. The problem statement of P|d j = d|Y w is introduced in "Problem statement" section. The design and application of meta-heuristic algorithms for P|d j = d|Y w are described in detail in "Meta-heuristic algorithms for P|d j = d|Y w " section. "Largest density first algorithm" section gives a new method which is improved from LPT. "Computational experiments" section analyses the computational experiment results performed by these algorithms. In "Conclusions" section, some conclusions are given.

Problem statement
Generally, the scheduling problem can be defined as assigning a set of jobs to a set of machines under given constrained conditions (Sterna 2011). Hence, the parallel identical machines scheduling problem with weighted late work criterion, could be defined as follows. Given n jobs J = {J 1 , J 2 , . . . , J j , . . . , J n } and m identical machines M = {M 1 , M 2 , . . . , M i , . . . , M m }, each job J j (j = 1, 2, . . . , n) is mainly described by its processing time p j , due date d j and weight w j . Let d j represents the preferred completion time for this job, and weight w j represents the relative importance of this job. Each job can be executed on one of the machines and each machine can execute only one job at a time. In this paper, we focus on a common due date d for all jobs (i.e. d j = d) and we look for a non-preemption schedule that minimizes the total weighted late work Y w .
The completion time of job J j is denoted by C j , the late work Y j for this job is given by the following formula (cf. Fig. 1): According to this definition, we can conclude that the late work criterion evaluates the quality of a feasible solution based on the length of late parts of particular jobs. Late work concentrates the advantages of two kinds of parameters: tardiness and the number of tardy jobs.
The total weighted late work can be defined as the weighted sum of late works of all jobs, is given by the following formula:

Meta-heuristic algorithms for P|d j = d|Y w
In this section, we design three meta-heuristic algorithms to solve the problem P|d j = d|Y w , respectively. The parameters setting of these algorithms will be given in "Computational experiments" section.

ACS for P|d j = d|Y w
ACS algorithm is one of the successful variations of ant colony optimization (ACO) (Dorigo and Gambardella 1997). Recently, ACO methods have been widely used in various scheduling problems, such as single machine case considered by Merkle and Middendorf (2003) and Bauer et al. (1999), job-shop problem (Dorigo et al. 1996;Colorni et al. 1994) and flow-shop case (Stutzle 1998). These research literatures show that ACO methods perform much better than other heuristic algorithms in finding the optimal solution of some benchmark problems (Merkle and Middendorf 2003).

Definition of heuristic information
When assigning a job, an ant should first choose a machine M i as the processing machine randomly, and then add a job J j to the scheduling sequence of machine M i based on heuristic information and pheromone. The heuristic information is denoted by η ij , which indicates the expectation of selecting job J j to assign to machine M i . The value Late work definition schematic diagram. In this figure, d is the common due date, p j denotes the processing time of each job J j (j = 1, 2, . . . , n), and C j denotes the completion time of job J j , then the late work Y j for this job is the minimum of p j and max{0, C j − d} of η ij is calculated according to the heuristic rule named variation of modified due date rule (VMDD) (Merkle and Middendorf 2003), i.e., where T i indicates the total processing time of all jobs which are already assigned to machine M i .

Pheromone initialization
Similar to the definition of heuristic information, the pheromone is denoted by τ ij , indicates the pheromone value of choosing a job J j to assign to machine M i . We initialize the pheromone value as τ 0 = 1/nL H , where L H is the objective function value obtained by the initial solution.

Solution construction and pheromone updating
Initially, the algorithm generates an initial solution randomly and obtains the initial pheromone value τ 0 . Each ant should finish the solution construction process and pheromone updating process. An ant k located on machine M i chooses a next job J j based on the pseudo random proportional rule, i.e., If q > q 0 , the probability of choosing job J j as the next job could be formulated as follows: where q is a random uniform variable in [0,1], and q 0 (0 ≤ q 0 ≤ 1) is a parameter to determine the relative importance between exploitation of a priori knowledge and exploration of a new edge. α and β are factors whose value determines the relative influence of pheromone and heuristic information, respectively. N k is the set of jobs that have not been visited by ant k so far.
To avoid premature convergence, an ant should conduct local pheromone updating after selecting a job J j , according to the following rule: where ξ(0 < ξ < 1) represents the local pheromone evaporation rate. After all of the ants have constructed their valid solutions, the pheromone of the best solution found so far should be updated according to the following rule: where S best represents the best solution found so far and �τ best ij denotes the inverse of the objective function value obtained by S best . ρ is used to control the global pheromone evaporation rate. The algorithm repeats these steps above until meeting the terminal condition, then the S best will be the ultimate solution obtained by ACS. In this paper, we set the terminal condition as the maximum number of iteration.

GA for P|d j = d|Y w
The genetic algorithm (GA) is an optimization technique inspired by natural evolution, and it is widely used in various realistic scheduling problem such as berth allocation (Pratap et al. 2015). In the GA, the individuals (called chromosomes) in a population are encoded to present candidate solutions to an optimization problem. The evolution usually starts with an initial population generated randomly, then the fitness of each individual in the population will be evaluated and a selection mechanism is applied according to the fitness. Then the crossover and mutation operators are applied to generate the offspring. The algorithm repeats these steps until reach the maximum number of iteration.

Chromosomes encoding
In this paper, we adopt a two-dimensional encoding method. Each gene is constituted by a tuple: (x j , y j ), 1 ≤ j ≤ n, where x j denotes the number of the machine which job J j is assigned to, and y j represents that J j is the y j -th job in the job sequence on the current machine. For example, given five jobs and two machines, then we can construct a chromosome with five genes to represent a candidate solution: (1,2), (2,1), (2,2), (1,3), (1,1). The above chromosome means job J 1 is assigned to m 1 and it is the second job on that machine, and J 2 is assigned to m 2 and it is the first job to be executed on that machine, and so on.

Initial population and fitness function
The initial population is generated randomly. The fitness function evaluates the quality of an individual in the population. Since it is a minimization problem, the fitness function can be designed as follows: where Y w (i) is the objective function value (weighted late work) obtained by individual i, and M is a constant which guarantees that the fitness value is positive.

Selection
In this algorithm, the selection operator is based on roulette rule. The probability of selecting the i-th individual can be defined as: where popsize is the population size, and f(i) denotes the fitness value of the i-th individual.

Crossover
In order to generate the offspring, information between two selected parents should be exchanged. In this paper, we consider the single point crossover method. Assuming that the crossover rate is p c , a crossover process can be described as follows: 1. Generate a random number p ∈ (0, 1), if p > p c , then go to (2), else exit this process. 2. Generate a random integer number l ∈ [1, n] , where n is the number of all jobs, also the length of a chromosome. Then, exchange the genes from l to n between two selected parents. Go to (3). 3. Repeat (1) and (2) for all adjacent chromosomes.

Mutation
Mutation operator can help avoid getting trapped in local optimum. In this process, the procedure traverses all individuals in a population, for each individual, the algorithm judges whether it should perform the mutation operator according to the mutation rate p m . When an individual is selected to mutate, we randomly pick a job from it, and reassign this job to another random position which can be on the current machine or on other machines.

SA for P|d j = d|Y w
Simulated annealing (SA) is a random optimization algorithm based on Mente Carlo simulation, which was inspired by the similarity between the annealing process of solid matters and combinational optimization problems. The algorithm was widely applied to the combinational optimization problems (Kirkpatrick et al. 1983;Van Laarhoven and Aarts 1987;Hazir et al. 2008).

Initialization and termination condition
In this paper, we generate an initial solution randomly with an initial temperature. The algorithm is terminated after reaching a minimum temperature or exceeding a given number of iterations.

Neighborhood structure
In order to explore the solution space, two operators are given to generate a neighbor: job move (N 1 ) and jobs interchange (N 2 ).
In the job move operator, a neighbor is generated by selecting a job from a current solution randomly and moving it to another random position. Jobs interchange operator generates a new neighbor through swapping two jobs selected randomly from different machines. The SA selects an operator randomly at each iteration, and performs it ⌊m/2⌋ times, and then obtain a new neighbor.

Cooling scheme
In this paper, we consider the geometric cooling scheme, which was a typical scheme that widely used in many scheduling problems (Hasani et al. 2014). In this scheme, the temperature is updated according to: where 0 < θ < 1 is the cooling factor. T k is reduced after running I iter iterations.

Proposed algorithm procedure
The pseudo-code of the proposed algorithm is described as follows :

Largest density first algorithm
We propose a novel algorithm named LDF (largest density first). This algorithm is improved from LPT (longest processing time first). Considered the weight of each job, we define the density of a job J j as w j /p j . The main idea of this algorithm is assigning the jobs to machines based on their density. The job with the largest density will be scheduled first and assigned to the machine with minimum load. The pseudo-code of LDF is described as follows:
The parameters setting of GA are as follows: set the size of population PopSize = 30, the maximum number of iteration G max = 300, and p c = 0.8, p m = 0.01.
In the SA, set T 0 =5 · n i=1 p i , θ = 0.975, and I iter = 20. Note that a solution with Y w = 0 must be an optimal solution, so the algorithms will "early exit" when generate a solution with Y w =0.
In order to simulate different experimental cases, we also need to set five input parameters, which are the number of machines m, the number of jobs n, the processing time of each job p, the weight of each job w, and the common due date d. The settings of these input parameters are shown in Table 1. The value of p is an integer and it will be rounded down in Poisson distribution. µ is a parameter that belongs to the set 0.8, 0.9, 1.1, 1.2, which controls the value of d. For each group of the five input parameters, 20 different experimental instances are generated to run each algorithm 20 times. For example, set m = 2 and n = 3m (i.e. n = 6), assume p is generated randomly from U[1,10n] (i.e. U [1,60]), w is generated from U[1,100], set µ = 0.8 and we can obtain the corresponding d. Then, 20 different experimental instances are generated according to this group of parameters. The experiments are conducted on a server with win server operating system, E5-2407 CPU and 32G memory.

Results and analysis
For each group of experiments, six different algorithms are all carried out. Because the orders of magnitudes of Y w are quite different under different parameter values, we use Y w / n i=1 p i w i to evaluate and compare the quality of a solution. The average values of Y w / n i=1 p i w i (in percent, denoted by Y ′ Avg) and the average running time TimeAvg (in seconds) of every 20 test instances are computed. Figs. 2, 3, 4, 5 and 6 show the comparison results in Y ′ Avg * , which are merged by Y ′ Avg with the same m, n, p, w and µ, respectively, due to the space limitations.
As we can see from Figs. 2 and 3, SA performs best in most cases, and LDF is slightly better than SA when the problem scale is large, i.e., m = 10 or n = 15m. In Fig. 2, the performance of ACS is significantly influenced by the parameter m which related to the problem scale. ACS is getting worse and worse with the increasing number of machines. Fig. 3 shows that LDF is greatly influenced by the parameter n. LDF is worse than three meta-heuristic algorithms when n = 3m but better than ACS and GA with the increasing ratio from 5 to 15.
In Fig. 4, we can see that the distributions of p have a greater impact on LDF compared with others. LDF is much worse than SA when p obeys U [1,10n] or U [1,20n]. However, it outperforms SA when p obeys Poisson distribution or U[100−5,100+5]. That means LDF can get better performance when the job procession time is concentrated.
The distributions of w gives no evident influence on all these algorithms and SA performs the best in all cases according to Fig. 5. Figure 6 shows the comparison results merged by different common due dates. When µ = 1.1 and µ = 1.2, almost all of these algorithms can obtain the optimal solutions (Y ′ Avg = 0 %), and LDF is a little worse than three meta-heuristic algorithms. When µ = 0.8 and µ = 0.9, the performance of proposed algorithms from high to low are SA, LDF, GA and ACS.
According to Figs. 2, 3, 4, 5 and 6, we can conclude that SA performs best in most cases. LDF also shows excellent performances and it is even better than SA in certain conditions. Though GA and ACS are not so good, they are greatly better than LPT and LS.
The running times of these algorithms are mainly affected by the scale of problem. Note that in meta-heuristic algorithms, we adopt an "early exit" mechanism proposed  Table 4, the values of meta-heuristic algorithms in µ = 0.8 and µ = 0.9 present the average running time merged by µ without "early exit" mechanism and SA is the best.
To make a comparison among the meta-heuristic algorithms in detail, the average number of better solutions is counted. The comparison results of meta-heuristic algorithms in the percent of better solutions are shown in Table 5, which is merged by µ.
The first column of Table 5 presents the values of µ. The rest columns are divided into three groups. In each group, two meta-heuristic algorithms are compared. The first Fig. 5 Comparison results in Y ′ Avg * merged by w. This figure shows the comparison results of six algorithms in Y ′ Avg * which are merged by Y ′ Avg with the same w. Y ′ Avg is denoted as the average value of Y w / n i=1 p i w i Fig. 6 Comparison results in Y ′ Avg * merged by µ. This figure shows the comparison results of six algorithms in Y ′ Avg * which are merged by Y ′ Avg with the same µ. Y ′ Avg is denoted as the average value of Y w / n i=1 p i w i column and the second column in each group show the average amount of better solutions obtained by these two algorithms, respectively. For example, the value "12.64" at the first row and the third column means that when µ = 0.8 and traversing all other parameters, the GA generated an average number of 12.64 solutions which are better than that of ACS. At the same situation, the number of solutions which have equal quality in ACS and GA can be calculated by 20 − 6.32 − 12.64 = 1.04. As we can see from Table 5, SA can usually generate better solutions, apart from the solutions with equal quality to other algorithms.
Combining the above analysis, we can conclude that SA performs best in most cases while its time cost is also acceptable. ACS and GA are also much better than the typical

Conclusions
In this paper, we proposed three meta-heuristic algorithms (i.e., ACS, GA and SA) and a novel LDF algorithm to solve the scheduling problem P|d j = d|Y w for the first time.
The proposed algorithms were compared to two typical scheduling algorithms LPT and LS. The computational experiments demonstrated that SA performed best in terms of finding optimal solutions compared to others in most cases, while the runtime of SA was acceptable in practical problems. The quality of solutions obtained by LDF is better than SA when the problem scale is large or the job processing time is concentrated. LDF also have the advantages in short running time. GA is better than ACS both in the quality of solutions and the running time of algorithm. In terms of the quality of solutions, the overall order in performance from high to low is SA, LDF, GA, ACS, LPT and LS.