September 6, 2013

Engineering Optimization

PMSDST˙HDCS˙EO˙92˙arxiv

To appear in Engineering Optimization
Vol. 00, No. 00, Month 20XX, 1–22

RESEARCH ARTICLE
Parallel machine scheduling with step deteriorating jobs and setup
times by a hybrid discrete cuckoo search algorithm

arXiv:1309.1453v1 [math.OC] 2 Sep 2013

Peng Guoa Wenming Cheng

a

and Yi Wangb∗

a

b

School of Mechanical Engineering, Southwest Jiaotong University, Chengdu, China;
Department of Mathematics, Auburn University at Montgomery, Montgomery, AL, USA
(Received 00 Month 20XX; final version received 00 Month 20XX)

This article considers the parallel machine scheduling problem with step-deteriorating jobs
and sequence-dependent setup times. The objective is to minimize the total tardiness by
determining the allocation and sequence of jobs on identical parallel machines. In this problem, the processing time of each job is a step function dependent upon its starting time.
An individual extended time is penalized when the starting time of a job is later than a
specific deterioration date. The possibility of deterioration of a job makes the parallel machine scheduling problem more challenging than ordinary ones. A mixed integer programming
model for the optimal solution is derived. Due to its NP-hard nature, a hybrid discrete cuckoo
search algorithm is proposed to solve this problem. In order to generate a good initial swarm,
a modified heuristic named the MBHG is incorporated into the initialization of population.
Several discrete operators are proposed in the random walk of Lévy Flights and the crossover
search. Moreover, a local search procedure based on variable neighborhood descent is integrated into the algorithm as a hybrid strategy in order to improve the quality of elite
solutions. Computational experiments are executed on two sets of randomly generated test
instances. The results show that the proposed hybrid algorithm can yield better solutions
in comparison with the commercial solver CPLEX with one hour time limit, discrete cuckoo
search algorithm and the existing variable neighborhood search algorithm.

Keywords: parallel machine scheduling; step-deterioration; sequence-dependent setup
time; hybrid cuckoo search; total tardiness

1.

Introduction

Production scheduling as a key decision-making process has a great effect on the performance of manufacturing and service systems. There have been a great number of research
works that are involved in the classical scheduling problem. These studies assumed that
the processing time of each job is constant throughout the scheduling period. However,
this assumption may not be true in some real industrial settings. Examples can be found
in steel production, equipment maintenance, and medical emergency, etc., where any
delay or waiting in starting to process a job may incur extended time for its completion.
Such kinds of jobs are called deteriorating jobs. The corresponding scheduling problem
was firstly introduced by Gupta and Gupta (1988) and Browne and Yechiali (1990). Since
then, many scheduling problems with deteriorating jobs have been extensively studied
∗ Corresponding

author. Email: ywang2@aum.edu

1

September 6, 2013

Engineering Optimization

PMSDST˙HDCS˙EO˙92˙arxiv

from various aspects. For details on this stream of research, Alidaee and Womer (1999),
Cheng, Ding, and Lin (2004) and Gawiejnowicz (2008) provided comprehensive surveys
of different models and problems. More recent literature that has explored scheduling
problems with deteriorating jobs include Toksari and Guener (2008); Wang, Jiang, and
Wang (2009); Bahalke, Yolmeh, and Shahanaghi (2010); Cheng, Lee, and Wu (2011);
Huang and Wang (2011); Bank et al. (2012); Jafari and Moslehi (2012); Yin et al. (2013)
and etc.
Among various types of scheduling problems involving deteriorating jobs, those jobs
with step deteriorating effect are considered in this article. A step-deteriorating job means
that if a job fails to be processed prior to a pre-specified threshold, its processing time will
be extended by adding an extra penalty time. The corresponding single machine scheduling problem may date back to Sundararaghavan and Kunnathur (1994), who firstly gave
some solvable cases for minimizing the sum of the weighted completion times. Subsequently, Mosheiov (1995) studied the makespan minimization problem for both singleand multi-machine cases. He also suggested several heuristics for all these complex problems. Subsequently, Jeng and Lin (2004) considered the problem on a single machine with
multiple deteriorating dates, and designed a branch and bound algorithm for deriving
optimal solutions. Cheng and Ding (2001) proved that the total flowtime minimization
problem is NP-complete. Jeng and Lin (2005) developed a branch and bound algorithm
for the problem. Owing to the intractability of the problem, He, Wu, and Lee (2009)
used weighted values of the basic processing times, the deteriorating thresholds and the
penalties to derive a heuristic solution. Cheng et al. (2012) employed a variable neighborhood search algorithm to minimize the total completion time in a parallel machines
environment. Moreover, Layegh, Jolai, and Amalnik (2009) used a memetic algorithm
to solve the single machine total weighted completion time problem. Guo, Cheng, and
Wang (2013) used a general variable neighborhood search to solve a single machine total
tardiness problem. In particular, batch scheduling with step-deterioration has been researched (Leung, Ng, and Cheng 2008; Mor and Mosheiov 2012). In addition, a similar
model which assumes that the processing time of a job is a piecewise linear function
(rather than a step function) of the waiting time is introduced by Kunnathur and Gupta
(1990). Other related works can be referred to Kubiak and Van de Velde (1998); Cheng
et al. (2003); Ji and Cheng (2007); Wu et al. (2009); Moslehi and Jafari (2010); Farahani
and Hosseini (2013).
Sequence dependent setup is commonly encountered in manufacturing environments
(Allahverdi, Gupta, and Aldowaisan 1999; Allahverdi et al. 2008). The research of
scheduling of deteriorating jobs with consideration of setup times is relatively uncommon
compared with those ordinary problems with setup times. The single machine scheduling
problems with past-sequence-dependent setup times and deteriorating jobs were studied
by Zhao and Tang (2010); Cheng, Lee, and Wu (2011); Lai, Lee, and Chen (2011). In
their works, the actual processing time of a job is formulated as an increasing function
of its scheduled position and the actual processing times of jobs already processed. In
addition, Bahalke, Yolmeh, and Shahanaghi (2010) suggested a genetic algorithm and a
tabu search for a single machine scheduling problem with the objective of minimizing
the makespan. Moreover, the branch and bound technique based on some dominance
properties has been adopted by Cheng et al. (2011); Lee, Lin, and Shiau (2011) for
minimizing the maximum tardiness and the number of late jobs, respectively. To the
best of our knowledge, there is no research that concerns the scheduling problem with
step-deteriorating jobs and sequence dependent setup times.
Parallel machine scheduling problem, as a generalization of the single machine problem,
has been widely studied in the literature because it frequently appears in the industrial
production. In this article, we study the parallel machine scheduling problem in which
the two effects of both step deterioration and sequence dependent setup times are con-

2

September 6, 2013

Engineering Optimization

PMSDST˙HDCS˙EO˙92˙arxiv

sidered. Below the abbreviation “PMSDST” will be used to denote the problem under
consideration. Since the minimization of the total tardiness on uniform parallel machines
with sequence-dependent setup times is NP-hard (Armentano and De Franca Filho 2007),
the PMSDST, as a more general problem, is also NP-hard. Though the PMSDST can be
formulated to a mathematical programming model, the model is impossible to be solved
in a reasonable amount of computational time when the problem size increases. Therefore, we propose a hybrid meta-heuristic algorithm based on the cuckoo search to obtain
a near-optimal schedule in a reasonable time. Computational results and comparisons to
other existing algorithms demonstrate the effectiveness and the efficiency of the proposed
meta-heuristic algorithm.
The rest of the article is organized as follows. The next section describes the problem
in more details, and formulates a mixed integer programming model. Section 3 proposes
the search heuristic suggested in our study together with the methods to generate initial
and neighborhood solutions. Then the computational results are reported and discussed
in section 4. Finally, section 5 gives some conclusions along with remarks for future
research.

2.

Problem formulation

In the PMSDST, a set of n independent jobs have to be processed on m identical parallel
machines. All machines are available at time zero, and each machine can handle at most
one job at a time. Preemption, division, or cancelations are not allowed. For job j, its
actual processing time pj is a non-linear step function of its starting time sj and the
deteriorating date hj . If the starting time sj of job j is less than or equal to the corresponding threshold hj , then job j only requires a basic processing time aj . Otherwise,
an extra penalty bj is incurred. Each job j is given a due dates dj . In addition, a sequence
dependent setup time for a job is also considered in this problem. That is to say, there
is a setup time δij when job i precedes immediately job j on the same machine. As a
matter of fact, δjj = 0. In addition, it is assumed that δ0j = 0, which indicates that no
setup time is needed if job j is processed at the first position on a machine. Without loss
of generality, all parameters are assumed to be integers. The tardiness Tj of a job j can
be determined by Tj = max{0, Cj − dj }, where Cj denotes the completion time of job
j. Finally, the goal of the scheduling problem considered here is to find a schedule that
specifies an assignment of all jobs to various machines and the schedule of those jobs on
all machines so that the total tardiness is minimized.
We next propose a mixed integer programming (MIP) model to minimize the total
tardiness for the PMSDST. The basic structure of the MIP model is motivated by a
network model formulated by Radhakrishnan and Ventura (2000). In order to present
constraints conveniently, two dummy jobs 0 and n + 1 are introduced on each machine
and are scheduled at the head and tail of its job sequence, and its processing time is set
to 0. Additionally, the following notation is defined.
• ukij : the binary decision variable that is equal to 1 if job i precedes immediately job j
on machine k, and 0 otherwise.
P
n
•
j=1 Tj : the total tardiness.
• M : a large positive constant such that M → ∞ as n → ∞. In this article, we choose
n
X
M = max dj +
(aj + bj ).
j=1,...,n

j=1

Applying the above notations and variables, the PMSDST can be formulated as follows.
Objective function:

3

September 6, 2013

Engineering Optimization

PMSDST˙HDCS˙EO˙92˙arxiv

min Z =

n
X

Tj

(1)

j=1

Subject to:

(
aj ,
s j 6 hj
, ∀j = 1, . . . , n
pj =
aj + bj , otherwise
n
X

(2)

uk0i = 1, ∀k = 1, . . . , m

(3)

uki(n+1) = 1, ∀k = 1, . . . , m

(4)

i=1
n
X
i=1
n X
m
X

ukij = 1, ∀j = 1, . . . , n, i 6= j

(5)

ukij = 1, ∀i = 1, . . . , n, i 6= j

(6)

i=0 k=1
n+1
m
XX
j=1 k=1

sj > δ0j + M (uk0j − 1), ∀j = 1, . . . , n, ∀k = 1, . . . , m
sj > Ci + δij + M (ukij − 1), ∀i, j = 1, . . . , n, ∀k = 1, . . . , m

(7)
(8)

Cj > sj + pj , ∀j = 1, . . . , n

(9)

Tj > Cj − dj , ∀j = 1, . . . , n

(10)

ukij ∈ {0, 1}, ∀i, j = 1, . . . , n, ∀k = 1, . . . , m
sj , Cj , Tj > 0, ∀j = 1, . . . , n

(11)
(12)

In the above formulation, equation (1) is the objective function that minimizes the sum
of the total tardiness. Constraints (2) define the actual processing time of a job under
the consideration of step-deteriorating effect. Constraints (3) ensure that the dummy job
0 is positioned at the beginning of the sequence before all the real jobs on each machine.
Constraints (4) guarantee that the dummy job n + 1 is processed just after all real jobs
on each machine. Constraints (5) and (6) ensure the precedence relation of jobs assigned
to the same machine. Constraints (7) state that the starting time for the first job of each
machine must be equal to or greater than its initial setup time. For other jobs j, the
starting time sj is determined by constraints (8) and is at least the sum of the completion
time of job i and the setup time from i to j. If i is not the predecessor of j, the subtraction
of M makes constraints (8) non-restrictive. Constraints (9) calculate the completion time
of each job. For each job j (j = 1, . . . , n), the completion time Cj is equal to its starting
time sj plus its actual processing time pj . Constraints (10) ensure that only the positive
value of lateness can be considered as tardiness, that is, Tj = max{0, Cj −dj }. Constraints
(11) declare ukij as binary variables. Finally, constraints (12) state that the starting time,
completion time and tardiness of job j are non-negative.
The above mathematical programming model is used to find optimal solution for our
problem. However, due to the NP-hardness of the proposed problem, it is difficult to be

4

September 6, 2013

Engineering Optimization

PMSDST˙HDCS˙EO˙92˙arxiv

solved optimally for medium- and large-sized instances. This is the main characteristic
of most complex scheduling problems, which makes necessary to develop some metaheuristic algorithms for solving the problem under study. In the next section, we propose
a hybrid cuckoo search algorithm to obtain near-optimal solutions.

3.

Hybrid discrete cuckoo search algorithm

The Cuckoo Search (CS) algorithm is one of the most new swarm intelligence optimization algorithms introduced by Yang and Deb (2009, 2010), inspired by the intelligent
breeding behavior of some species of cuckoo. These species lay their eggs in the nests of
host birds (of other species) with amazing abilities such as selecting the recently spawned
nests and removing existing eggs so that the hatching probability of their own eggs is
increased. On the other hand, some host birds are able to engage direct contest with
infringing cuckoos. For instance, if the alien eggs are discovered by the host bird, it may
either throw out these eggs or abandon the nest completely. This natural phenomenon
has led to the evolution of cuckoo eggs to mimic the egg appearance of local host birds.
For simplicity in the implementation of the algorithm, the following three idealized
rules are used (Yang and Deb 2009): (1) Each cuckoo lays one egg at a time, and dumps
its egg in a randomly chosen nest; (2) The best nests with high quality of eggs will carry
over to the next generation; (3) The number of available host nests is fixed, and the egg
laid by a cuckoo is discovered by the host bird with a probability ρa ∈ [0, 1]. In this case,
the host bird can either discard the egg or abandon the nest so as to build a new nest
in a new location. The last assumption can be approximated by a fraction ρa of the n
nests that shall be replaced by new ones.
The CS algorithm contains a population of P eggs. Egg ı ∈ P is in nest Xı , which
corresponds to a solution of the underlying problem. The word ‘nest’ and the word
‘solution’ shall be used interchangeably. In the initial phase, P solutions are generated
randomly. Then all of these solutions except the one with the smallest objective value
are replaced by new solutions that shall be generated by the Lévy flights. Specifically, if
t is the current iteration index, and let Φt := {Xjt : j ∈ P} be the set of all solutions at
(t+1)

the t-th iteration, then a candidate solution X̃ı
X̃ı(t+1)

=

Xı(t)

+ α0



X(t)

−

is generated by the equation
Xı(t)



· t−λ ,

(13)
(t)

where, the constant parameter α0 is related to the scale of the underlying problem, X
represents a randomly selected solution at the t-th iteration and λ is a constant such
that 1 < λ ≤ 3.
(t)
(t)
The vector α := α0 (X − Xı ) introduces a difference from the current solution
(t)
Xı . This difference mimics the scenario in which a cuckoo’s egg is more difficult to be
recognized by a host if it is more similar to the host’s egg. The term t−λ accounts for
the heavy-tailed power law distribution Lévy(λ) given by
Lévy(λ) ∼ t−λ .

(14)

The distribution Lévy(λ) has infinite variance and mean when 1 < λ ≤ 3, thus the Lévy
flight essentially provides a random walk process in hope to search for a better solution.
If the host bird identifies a cuckoo’s egg in its nest (with the probability ρa ), either the
alien egg shall be discarded by the host or the host shall abandon the nest. In either case,
the cuckoo’s egg will fail to hatch. Thus the cuckoo has to explore a new nest (solution)

5

September 6, 2013

Engineering Optimization

PMSDST˙HDCS˙EO˙92˙arxiv

in order to successfully hatch its egg. This behavior shall be simulated by the crossover
operation. Specifically, let randı , rand be two uniformly distributed random numbers on
the interval [0, 1]. A new candidate solution Vıt+1 is obtained by the equation
Vıt+1

(
X̃ıt+1 + rand · (X̃γt+1
− X̃γt+1
) if randı > ρa
1
2
=
t+1
X̃ı
otherwise,

(15)

where, γ1 and γ2 are two randomly chosen indices such that ı, γ1 and γ2 are pairwise
distinct. If the candidate solution Vıt+1 has a better quality than the old one Xıt , that is,
with a smaller objective function value, then Vıt+1 will be accepted as Xıt+1 . Therefore,
for ı ∈ P, we have at the (t + 1)th iteration,
Xıt+1


=

Vıt+1 , if f (Vıt+1 ) < f (Xıt )
X̃ıt+1 , otherwise,

(16)

where, f is the objective function of the underlying optimization problem. Afterwards,
the best solution among all solutions found so far up to the present iteration is recorded
to memory. The same procedure repeats until a termination criterion is met.
The CS algorithm provides more robust and precise results than other meta-heuristic
algorithms, such as particle swarm optimization, differential evolution, or artificial bee
colony algorithms in solving the continuous optimization problems (Civicioglu and Besdok 2013). Therefore, the CS algorithm has aroused many interests and has been successfully applied to different kinds of optimization problems (Walton et al. 2011; Gandomi,
Yang, and Alavi 2013; Valian et al. 2013; Yildiz 2013). Recently, Ouaarab, Ahiod, and
Yang (Online First 2013) proposed a discrete CS algorithm to solve the traveling salesman problem. Burnwal and Deb (2013) developed a cuckoo search-based approach for
scheduling optimization of a flexible manufacturing system. Li and Yin (2013) designed
a hybrid algorithm that combines the CS and a fast local search to solve a permutation
flow shop scheduling problem.
In the remainder of this section, we propose a hybrid discrete CS algorithm (HDCS) to
solve the considered PMSDST. In order to maintain the remarkable characteristics of the
standard CS, the Lévy flight equation (13) and the crossover operator (15) are redefined
based on an integer encoding scheme. Moreover, a local search procedure based on variable neighborhood descent is employed to refine elite solutions that are not discarded.
In addition, a heuristic named the MBHG is incorporated into the random initialization
to generate a population of initial solutions with certain quality. Furthermore, to avoid
the algorithm trapping into a local optimum, a restarting strategy is embedded into
the search process. The details of the proposed hybrid discrete algorithm are elaborated
below.

3.1

Solution representation

In a meta-heuristic, the solution performance is highly dependent upon the representation of a solution. The most commonly used solution representation for parallel machine
scheduling problem is having an array of jobs for each machine that represents the processing order of jobs assigned to that machine. Nevertheless, in this article, we represent
a solution by a single sequence of the n jobs in order to maintain the simplicity and
the powerful search performance of the CS algorithm. With such a sequence, an actual
schedule can be built simply by assigning the next unscheduled job to the earliest available machine until all jobs are scheduled. When multiple machines are available for an
unscheduled job, the machine with the smallest index will be chosen. The obtained se-

6

Engineering Optimization

PMSDST˙HDCS˙EO˙92˙arxiv

quence of jobs shall be called a solution vector or a nest vector. Once an actual schedule
is constructed, the total tardiness can be easily calculated subsequently.
To illustrate how the nest vector is used to construct an actual schedule, below a simple
example is given with six jobs and two machines. The detailed data for this problem is
provided in Table 1. Suppose a nest vector for the problem is given by X = [2, 6, 4, 1, 5, 3].
An actual schedule can be obtained through the aforementioned decoding scheme, as
shown in Figure 1. Correspondingly, the total tardiness is 65. In Figure 1, jobs 2, 4 and
5 are assigned to machine 1 and jobs 6, 1 and 3 are assigned to machine 2.
Table 1. The input data of the example with six jobs and two machines

Job(j)

1

2

3

4

5

6

Basic processing time aj
Due date dj
Deteriorating date hj
Penalty bj

78
85
70
18
0
5
2
3
3
5

17
48
4
33
9
0
6
5
10
8

97
229
62
1
9
8
0
10
6
4

93
133
19
17
5
2
5
0
9
8

62
220
58
40
4
2
4
3
0
4

53
75
39
31
6
5
7
9
5
0

Setup time δij

Job

1

2

3

4

5

6

Solution

2

6

4

1

5

3

Machine

1

2

Setup Time

Total Tardiness=65

Processing Machine

September 6, 2013

Penalty

6

M2

M1

1

2

4
40

5+

5
80

3+

3

120

160

200

240

Gantt Chart
Figure 1. The decoding scheme of a solution

3.2

Population initialization

In the original CS algorithm, the initial population of solutions is often filled with P ndimensional vectors that are randomly generated. In order to ensure an initial population
with certain quality and diversity, the BHG heuristic (Biskup, Herrmann, and Gupta

7

September 6, 2013

Engineering Optimization

PMSDST˙HDCS˙EO˙92˙arxiv

2008) is modified to generate an initial solution for our problem, whereas the rest P − 1
solutions are initialized with random job permutations.
Generally, a job j is likely to be scheduled in an earlier position if it has a smaller value
of due date dj and a smaller value of deteriorating date hj for our problem. Therefore,
in the modified BHG heuristic, a predefined list of jobs is obtained in ascending order
of their weighted values ωdj + (1 − ω)hj , where 0 < ω < 1. Moreover, in order to
reduce computation time, we simply assign the first m jobs from the predefined list
simultaneously to the m machines, respectively. The detailed modified BHG (MBHG)
heuristic algorithm is described next.
(1) Obtain a list of all jobs in ascending order of the weighted value of due date and
deteriorating date, i.e. ωdj + (1 − ω)hj 6 ωdj+1 + (1 − ω)hj+1 for j = 1, . . . , n − 1.
(2) The first m jobs from the list are simultaneously assigned to the m machines, respectively.
(3) Let U = {m + 1, . . . , n} be the set of unscheduled jobs.
(4) Determine to which machine the job with the smallest index in U should be assigned.
Starting from machine 1, take the job with the smallest index in U and insert it
before and after each job that is already assigned to machine 1. The insertion into
the schedule of machine 1 may start from the position after the last job and continue
sequentially backward to the position before the first job assigned to machine 1. For
each instance of inserting a job into the current schedule of machine 1, calculate a
new value of the total tardiness. The same procedure is applied to machines 2, 3,
. . . , m afterwards. After all the possible insertions into the schedule of each machine
are done, the algorithm chooses the sequence with the lowest total tardiness which
was calculated first. Remove the first job from the set U .
(5) Repeat Step (4) until U is empty.
(6) Output the final schedule on the m machines and its total tardiness.
It is found that the weight ω has a significant influence on the solution quality of the
MBHG. When the weight is ω = 0.1, the solution provided by the MBHG is 116 for the
above mentioned instance listed in table 1. When the weight is ω = 0.5, the solution
obtained by the MBHG is 65 for the same instance. It is hard to find the best value
of weight ω. This motivates the search with different values of ω. In this study, weight
ω is varied from 0.1 to 0.9. The best solution with the least total tardiness among the
solutions obtained from different weights is selected as the final output of the MBHG
heuristics.

3.3

CS-based search

The search process of the original CS algorithm (equations (13), (15), (16)) may fall
into a local optimal solution, particularly in solving a discrete optimization problem.
Therefore, in order to maintain diversity of the population, the following search strategies
are adopted.
For the Lévy flight, a uniformly distributed random variable ψ is firstly generated in
t
the interval (0, 1). Let Xbest
be the best solution found among all solutions in Φt . If
t .
ψ > 0.5, the generation of a new solutions is directed by the current best solution Xbest
Otherwise, the searching process is around a randomly selected solution Xt from the
current population Φt . Parameter λ is critical in delivering a new solution, as observed
from equation (13). A smaller λ increases the global exploration ability while a large λ
tends to improve the local exploitation ability of solutions in the current search area.
Therefore, to balance the two abilities, λ should be linearly increased from a relatively
small value to a relatively large value in the entire process of the HDCS algorithm. By

8

September 6, 2013

Engineering Optimization

PMSDST˙HDCS˙EO˙92˙arxiv

doing so, the HDCS can possess good global search ability at the beginning phase of
the search while having more local search ability near the end of the iterations. In this
article, parameter λ is chose to vary through the formula
λ = λmin +

λmax − λmin
t,
tmax

(17)

where, λmin and λmax are the initial value and the terminal value of λ respectively, and
tmax is the maximum number of iterations. In our algorithm, λmin = 1.1 and λmax = 3.
Based on the above modification, equation (13) of producing a new solution X̃ıt+1 is
replaced by the following equation

(
X̃ıt+1

=

t ) · t−λ ,
Xıt + α0 · (Xıt − Xbest
t
t
t
Xı + α0 · (Xı − X ) · t−λ ,

ψ > 0.5
otherwise.

(18)

In this article, we set α0 = 1.
Since the standard CS algorithm is originally designed for solving a continuous optimization problem, the operations such as subtraction, multiplication and addition involved in the Lévy flight must be redefined for a discrete scheduling problem. These
operations shall be redefined by some neighborhood search strategies to suite a discrete
problem. Next, the discrete operators used in equation (18) are described in detail.
• Subtraction
Our algorithm involves one subtraction operation between two solutions, say, the minuend X1 and the subtrahend X2 . Since the solution is a job permutation, a comparison
strategy is carried out for producing a new solution chain X 0 . If the given elements of
two solutions are identical, the corresponding element of the new chain is set equal to the
number zero. Otherwise, the element of the new chain is copied from the corresponding
one of X1 . Figure 2 shows the redefined subtraction operation.
X1

6

4

1

5

3

2

X2

1

4

2

6

3

5

X'

6

0

1

5

0

2

Figure 2. The subtraction operation X1 − X2 .

• Multiplication
The multiplication operator carries out the product of the scalar σ := t−λ and a
solution chain X 0 obtained by the subtraction operation. Let X 00 denote the output of
the multiplication. In the redefined manner of the multiplication operation, a uniformly
distributed random number between 0 and 1 is generated for each element in the chain
X 0 . The random number is then compared with σ. If the random number is more than or
equal to σ, the element in the offspring chain X 00 is set equal to the element of the same
position in the chain X 0 ; otherwise, its corresponding position is occupied by the number
zero. Figure 3 represents the process of the multiplication operation on the scalar σ and
a solution chain.
• Addition

9

September 6, 2013

Engineering Optimization

PMSDST˙HDCS˙EO˙92˙arxiv

ı=0.15

6

X'

0

1

5

0

2

rand(0, 1) 0.81 0.91 0.13 0.92 0.63 0.09

6

X''

0

0

5

0

0

Figure 3. The multiplication operation

The addition operation is modified by a number of swap moves. Two objects manipulated in the operation are the current solution X and the chain X 00 produced from the
multiplication. Once the addition operation is finished, a new solution Xnew is generated.
Initially, Xnew is set equal to the incumbent X. Then for each element X 00 (i) of X 00 , if
X 00 (i) is more than 0, the job in position X 00 (i) of Xnew and the job in position Xnew (i)
of Xnew are swapped. Figure 4 illustrates the manner in which the addition performs.
(a) Copying the job sequence
X

1

4

2

6

3

5

X''

6

0

0

5

0

0

Xnew

5

4

2

6

1

3

1

4

2

6

3

5

(b) Exchanging the jobs in positions 1 and 6

5

4

2

6

3

1

(c) Exchanging the jobs in positions 5 and 6
Figure 4. The addition operation

As it can be seen, the Lévy flight is redefined by the above three discrete operators.
Likewise, the crossover operation described by equation (15) is replaced by the order
crossover operator. The order crossover (Davis 1985) has shown robust performance in
the genetic algorithm. For each solution, a random number is generated and compared
with the probability ρa to determine whether the crossover operator should be executed.
If the random number is greater than ρa , the solution is selected as a parent for producing
a new offspring. Another solution is randomly chosen from the population of solutions and
regarded as a second parent. With the two parent solutions chosen, the order crossover
operation is performed as below.
(1) Select randomly a subsequence from one parent, called the first parent X1 .
(2) Produce a proto-child by copying the subsequence into its corresponding position.
(3) Check all the jobs of the other parent, called the second parent X2 , from left to right.
If a job is already assigned in the subsequence from X1 , then skip to the next job in
X2 that is not scheduled yet and assign it to the first available position of the new
solution. The jobs taken from X2 are the ones that the proto-child needs.
The above crossover operator will produce two offspring solutions because each parent
may be chosen as the first parent. The incumbent solution shall be updated by the better
one of the two offsprings in the next generation. The procedure is illustrated in figure 5.

3.4

Local search

To improve the performance, the HDCS implements variable neighborhood descent
(VND) as a local search strategy for a set of elite solutions. The elite set is constructed
by finding the best τ number of solutions from the population in each iteration. The

10

September 6, 2013

Engineering Optimization

PMSDST˙HDCS˙EO˙92˙arxiv

Subsequence

Subsequence

X1

1

5

4

2

3

6

X2

4

6

3

1

2

5

X2

4

6

3

1

2

5

X1

1

5

4

2

3

6

X'1

6

1

4

2

3

5

X'2

5

6

3

1

2

4

Figure 5. The order crossover operation

preset value of τ can be chosen by the following equation
τ = max{3, round[P × (1 − ρa )]}.

(19)

Equation (19) ensures that there are at least 3 elite solutions. If there are too few elite
solutions, the local search strategy can not improve the performance of the algorithm
effectively. On the other hand, too many elite solutions may consume too much computational time. Each elite solution indicates a promising region to search for the optimal
solution. Hence, the VND local search increases the possibility of obtaining an optimal
solution from diverse areas of the solution space. In our implementation, we use three
neighborhood structures, which are based on swap, insert and inverse operations. The
swap operation is done by randomly selecting the ith and jth jobs from an elite solution
and then swapping them directly. The insert operation is done by choosing job i at random from a solution and inserting it into the position immediately preceding a randomly
chosen job j of the solution. The inverse operation is done by inverting the subsection
between two randomly chosen jobs of a solution. Each of the three neighborhood operations is repeated n times when it is executed. An example of each operation is shown in
Figure 6.
1

2

3

4

5

6

1

2

3

4

5

6

1

2

3

4

5

6

swap

insert

inverse

1

5

3

4

2

6

1

3

4

5

2

6

1

4

3

2

5

6

Figure 6. The swap, insert and inverse operations for local search

The proposed variable neighborhood descent is applied to enhance the quality of an
elite solution. Each elite solution is regarded as a seed permutation. The search process
repeats each neighborhood for the best local optimum until no improvement appears. If
the local optimum is better than the seed solution, update the seed solution and continue
the search in the neighborhood; otherwise move to the next neighborhood and continue to
search there. The pseudo-code of the suggested local search is summarized in Algorithm
1.
3.5

Restarting strategy and stopping criterion

To maintain the diversity of the population of solutions, it is common to use restarting
strategy in evolutionary algorithms (Ruiz, Maroto, and Alcaraz 2006). In our proposed

11

September 6, 2013

Engineering Optimization

PMSDST˙HDCS˙EO˙92˙arxiv

Algorithm 1 Local search
Input: Elite solutions for the HDCS
Output: Improved elite solutions for the HDCS
1: for each solution X in the elite set do
2:
κ=1;
3:
while κ 6 3 do
4:
if κ == 1 then
0
5:
Apply Swap operations n times to X to obtain a local optimum X ;
6:
else if κ == 2 then
0
7:
Apply Insert operations n times to X to obtain a local optimum X ;
8:
else if κ == 3 then
0
9:
Apply Inverse operations n times to X to obtain a local optimum X ;
10:
end if
0
11:
if X is better than X then
0
12:
X←X ;
13:
Continue the kth neighborhood search operation;
14:
else
15:
κ = κ + 1;
16:
end if
17:
end while
18: end for
algorithm, partial solutions adopt a restarting mechanism at each iteration. Firstly, the
population is sorted in ascending order of objective function values. Then the first 90%
solutions from the sorted list are retained in the population. The remaining 10% solutions
shall be replaced with newly generated nests by random job permutations.
Regarding the stopping criterion, two commonly used termination conditions are employed in our proposed algorithm. That is to say, if the number of iterations t exceeds
the maximum number of iterations tmax or the best solution Xbest is not improved upon
in tnip successive iterations, the procedure is terminated; otherwise, a new iteration will
start and the iteration counter t is increased by one. When the algorithm terminates, the
best found solution Xbest and its total tardiness are output.

3.6

Outline of the HDCS

The search framework of the proposed HDCS algorithm for solving the PMSDST is
shown in Figure 7. The algorithm starts with an initial population of P solutions, among
which one solution is generated by the MBHG, and the rest is generated by random
permutations of n jobs. Then the solution Xbest with the smallest objective value is
identified by calculating the function value of each solution. Subsequently, the algorithm
iterates search for a better solution until a stopping condition is met. At each iteration,
the first τ elite solutions are determined according to equation (19) by their smaller
objective function values. The remaining solutions are regarded as normal solutions. For
each of the normal solutions, the modified Lévy flight by equation (18) is performed to
search for a new solution either around Xbest achieved in the previous iteration or a
randomly selected neighbor. The elite solutions are improved by a local search based on
the variable neighborhood descent. Then, a fraction ρa of the solutions are replaced by the
new candidate solutions produced by performing the order crossover operation. To avoid
premature convergence, the restarting strategy is utilized for the 10% worst solutions.
All the solutions are evaluated for their objective values and the Xbest is updated if a
new better solution is found. The algorithm is then continued to next iteration. When

12

September 6, 2013

Engineering Optimization

PMSDST˙HDCS˙EO˙92˙arxiv

the algorithm is terminated, the best solution found by the algorithm and its objective
value are output for the problem under study.

Figure 7. The flowchart for the proposed hybrid cuckoo search algorithm

4.

Computational experiments

There is no study on parallel machine scheduling problem with step-deterioration jobs
and setup times in the existing literature. Therefore, our article seems to be the first
study of this problem. In order to evaluate the effectiveness of the proposed algorithm, the
performance of the HDCS is compared with the discrete cuckoo search algorithm without
local search strategy (DCS) and the variable neighborhood search proposed by Cheng
et al. (2012) (VNS). The three algorithms are coded in Matlab 7.11 and executed on an
Intel Pentium dual-core 2.6 GHz PC with 2GB RAM under Windows XP environment.
The performance of each algorithm is measured by the relative percentage deviation

13

September 6, 2013

Engineering Optimization

PMSDST˙HDCS˙EO˙92˙arxiv

(RP D) calculated by the formula
RP D(%) =

Algsol − M insol
× 100
M insol

(20)

for each test instance, where Algsol is the objective function value obtained for a given
algorithm and M insol is the best solution of all experiments for each test problem.
4.1

Data generation of test problems

In order to analyze the performance of our algorithm for the scheduling problem under
consideration, different sizes of test instants are needed. In this section, two sets of
instances are defined. For small-sized instances, the following combinations of number n
of jobs and number m of machines are tested: n = {8, 10, 12, 14} and m = {2, 3, 4}. For
large instances, the combinations are n = {30, 40, 50, 60} and m = {4, 6, 8}. In all cases,
the following rules are used to generate job data.
The basic processing times aj of job j is generated from a discrete uniform distribution on the interval [1, 100]. The penalty bj is assumed to be drawn from the uniform
P
distribution U [1, 100×ξ], where ξ = 0.5 is chosen in our study. Set β = j=n
j=1 aj /m. The
deteriorating dates are generated randomly by a discrete uniform distribution on each of
the intervals H1 := [1, 0.5×β], H2 := [0.5×β, β] and H3 := [1, β]. The setup times are produced according to a uniform distribution in the range [1, 10]. Moreover, the due dates of
jobs are randomly generated from a uniform distribution on the interval [1, C̄max ], where
C̄max is the value of the maximal completion time obtained by the following rule: first of
all, a job permutation is achieved by arranging the jobs in the non-decreasing order of
the ratio aj /bj . Then, a corresponding schedule with no idle times can be constructed by
assigning unscheduled jobs to the earliest available machine. Once a schedule is built, the
maximal completion time C̄max is calculated. For each possible combination, 10 random
replicates are generated. Therefore, there are 360 small instances and 360 large instances
to test the algorithms.
4.2

Parameter calibration

Parameter selection can significantly affect the quality of a algorithm. In this study,
suitable parameter values were selected taking into consideration both solution quality
and computational efficiency. To determine appropriate values of parameters, extensive
computational tests were performed on 27 instances, including both small-sized and
large-size ones. These instances are generated from the above defined rule with the following parameter values: P=15,20,25,30,35,40; ρa =0.6,0.7,0.8,0.9; tmax =100,200,300,400;
tnip =30,50,70. The preliminary tests using these parameter values show that the following set of parameter values seems to provide the best performance within a reasonable
computational time: P = 30, ρa = 0.8, tmax = 200 and tnip = 50. Hence, this set of
parameter values will be adopted for all subsequent experiments in this study.
4.3

Results and discussion

In this section, we perform comparative evaluation of the proposed HDCS, the VNS, the
DCS, and a commercial solver CPLEX 12.5 based on the branch-and-cut algorithm on
the same benchmark of instances explained in Section 4.1. Due to the intractability of
the studied problem, the CPLEX 12.5 is only used to solve the MIP model with small
instances. A 3600-second time limit is imposed. When solved by the CPLEX, a particular

14

September 6, 2013

Engineering Optimization

PMSDST˙HDCS˙EO˙92˙arxiv

run is simply terminated and the best current integer solution is returned if the optimal
solution has not been found or verified in that amount of time.
Owing to the stochastic nature of the evaluated meta-heuristics, each algorithm is run
five times for a given instance to reach reliable results. The computational results of the
best (RPDB ), average (RPDA ) and worst (RPDW ) deviation of each combination of n
(number of jobs) and m (number of machines) obtained by the VNS, the DCS, the HDCS
are summarized in tables 2 and 4. Recall that for each combination of n and m in our
experiment, 10 random replicates are generated and for each instance the algorithm is
run 5 times. Thus corresponding to each combination of n and m, each value of RPDB
and RPDW is found in the following way: firstly the best (respectively, the worst) RPD
value is found among the five runs of an algorithm applied to a particular instance, then
the average value of the 10 best (respectively, the worst) RPD values of those 10 random
instances is recorded as RPDB (respectively, RPDW ) in tables 2 and 4. Likewise, the
value of RPDA is obtained in the following way: firstly the average RPD value is found
for the five runs of an algorithm applied to a particular instance, then the average value
of those 10 average RPD values of the 10 random instances is recorded as RPDA in tables
2 and 4.
Table 2 reports the experimental results of the small instances. Column 1 shows the
combination of n and m for an instance. Columns 2–5 report the computational results
of the four methods, respectively. If the CPLEX can yield the optimal solution, the RPD
value is calculated over the optimal total tardiness. The CPLEX solver is able to deliver
optimal solutions for all instances with eight jobs. For the 10-job instances, the CPLEX
is able to optimally solve all the instances with three and four machines, but only 23 of
30 instances with two machines. Regarding the 12- and 14-job cases, the CPLEX only
gives optimal solutions for a small portion of instances. In total, the CPLEX can yield
optimal solutions for 223 of the 360 small instances.
The results revealed by table 2 is exciting. We see that the HDCS can deliver the best
RPDB values except for the instances 10 × 2 with H1 and 14 × 3 with H2 . Even for these
two instances, the RPDB values given by the HDCS are only slightly larger than the
corresponding best values given by one of other three algorithms. The VNS is efficient
in solving a parallel machine flow time problem with step-deterioration (Cheng et al.
2012), however, it does not perform as well as the HDCS for the PMSDST problem.
This is because the consideration of setup times and the total tardiness criterion in the
PMSCST problem makes the VNS more easily get trapped into a local optimal solution.
In order to further analyze the results, the one-way analysis of variance (ANOVA) is
adopted to check whether the observed difference in the RPDA values for different algorithms are statistically significant when one of the deteriorating intervals (H1 , H2 and H3 )
is applied. The means plots along with the Tukey Honestly Significant Difference (HSD)
intervals at the 95% confidence level are given in figure 8 for each algorithm applied
to each of the three different deteriorating intervals. The means obtained by the four
algorithms for the same deteriorating interval (one of H1 , H2 and H3 ) are connected by
a line segment. If the Tukey HSD intervals of two algorithms for the same deteriorating
interval have overlapping, then the performances of the two algorithms are not statistically significantly different. It is clearly seen from figure 8, that the proposed HDCS is
statistically better than other methods by having its Tukey HSD intervals of the RPD
values lower than those of other methods when applied to each of the three intervals H1 ,
H2 and H3 .
The computation times of the four methods for each combination of n (number of jobs)
and m (number of machines) are also presented in table 3. They are obtained in a similar
manner as to obtain the RPDA . When the size of instance increases, the computation
time of the CPLEX grows significantly. Once the number of jobs exceeds 10, the CPLEX
completely exhausts the given time limit for most instances. The DCS consumes the

15

September 6, 2013

Engineering Optimization

PMSDST˙HDCS˙EO˙92˙arxiv

least time comparing with other algorithms, but its performance is the worst among all
algorithms. Alternatively, the computation times of the VNS and the HDCS are very
modest.
Table 2.

Average Relative Percentage Deviation (RPD) of small instances for the algorithms
VNS

Instance

MIP

8×2
8×3
8×4
10×2
10×3
10×4
H1
12×2
12×3
12×4
14×2
14×3
14×4
Average

DCS

HDCS

RPDB

RPDA

RPDW

RPDB

RPDA

RPDW

RPDB

RPDA

RPDW

0.00
0.00
0.00
0.00
0.00
0.00
0.36
0.00
0.00
12.34
1.13
1.51
1.28

0.00
0.00
0.00
0.44
0.00
0.00
0.00
0.14
0.24
0.00
0.12
0.02
0.08

0.22
0.00
0.00
0.96
0.13
0.34
0.72
1.76
2.58
2.74
1.44
1.14
1.00

0.31
0.00
0.00
1.99
0.45
1.50
2.32
5.83
5.18
8.02
3.13
2.71
2.62

0.31
0.00
0.00
6.43
1.14
0.58
15.05
14.26
56.32
53.13
19.91
19.60
15.56

1.63
1.50
0.72
11.13
5.41
3.83
18.06
17.11
58.86
55.09
23.26
20.75
18.11

4.79
5.72
1.21
14.54
7.35
5.79
19.87
18.97
60.87
55.80
24.60
21.83
20.11

0.00
0.00
0.00
0.05
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00

0.00
0.00
0.00
0.05
0.01
0.00
0.00
0.09
0.02
0.41
0.26
0.02
0.07

0.00
0.00
0.00
0.05
0.08
0.00
0.00
0.22
0.03
0.92
0.57
0.11
0.17

8×2
8×3
8×4
10×2
10×3
10×4
H2
12×2
12×3
12×4
14×2
14×3
14×4
Average

0.00
0.00
0.00
0.00
0.00
0.00
0.87
0.23
0.00
4.58
0.28
0.33
0.52

0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.18
0.00
0.03
0.00
0.85
0.09

0.29
0.00
0.00
0.23
0.15
0.04
2.73
1.97
1.02
1.15
1.17
2.31
0.92

0.59
0.00
0.00
0.65
0.33
0.06
10.07
4.18
1.86
3.13
3.13
4.69
2.39

0.11
0.00
0.00
1.39
7.00
1.03
25.28
8.00
4.87
31.56
30.89
16.57
10.56

0.58
0.04
0.00
1.64
8.24
4.51
26.65
9.11
6.94
36.69
36.37
17.47
12.35

0.75
0.18
0.00
1.99
9.96
6.83
27.42
9.66
8.76
38.33
38.96
18.01
13.40

0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.09
0.00
0.01

0.00
0.00
0.00
0.00
0.00
0.00
0.09
0.39
0.10
0.58
0.49
0.96
0.22

0.00
0.00
0.00
0.00
0.00
0.00
0.23
1.05
0.33
2.29
1.63
2.08
0.63

8×2
8×3
8×4
10×2
10×3
10×4
12×2
H3
12×3
12×4
14×2
14×3
14×4
Average

0.00
0.00
0.00
0.00
0.00
0.00
1.00
0.20
0.11
0.61
1.38
0.29
0.30

0.00
0.00
0.00
0.00
0.00
0.00
0.35
0.20
0.04
2.69
0.00
0.67
0.33

0.00
0.00
0.00
0.42
0.82
0.08
3.78
3.61
1.18
5.18
2.70
1.98
1.65

0.00
0.00
0.00
0.79
2.07
0.25
5.76
8.81
3.09
8.25
5.56
4.04
3.22

0.00
0.00
0.00
7.86
3.93
1.54
54.08
21.39
6.84
38.65
26.37
18.38
14.92

0.97
0.37
0.07
10.43
5.92
2.08
61.27
26.33
10.42
46.73
34.83
23.82
18.60

1.70
0.93
0.28
11.33
9.85
2.77
65.25
29.44
12.23
52.13
39.53
27.01
21.04

0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.04
0.00
0.00
0.00
0.00

0.00
0.00
0.00
0.00
0.00
0.00
0.04
0.00
0.35
1.48
0.10
0.60
0.21

0.00
0.00
0.00
0.00
0.00
0.00
0.14
0.00
1.06
3.14
0.52
0.95
0.48

In order to evaluate the algorithms when the size of the problem is increased, experiments are also performed using large-sized instances. Since the CPLEX already consumes
a lot of computational time in solving small-sized instances, it is excluded in this case.
The results of large-sized instances obtained by all the meta-heuristic algorithms are
listed in table 4. In the table, the best RPD value for a specific instance delivered by an
algorithm is boldfaced. We can see that the average results (RPDA ) of the HDCS are
always better than those of the DCS and the VNS. However, there are 15 instances in
which the quality of the best solutions (RPDB ) obtained by the HDCS is poorer than
that of the VNS. This might be due to the fact that more neighborhood structures are
used in the VNS. Nevertheless, the worst solutions (RPDW ) delivered by the HDCS are
significantly better than those of the DCS and the VNS.

16

September 6, 2013

Engineering Optimization

PMSDST˙HDCS˙EO˙92˙arxiv

Table 3. Computation times of small instances for the algorithms ( in seconds)
H1

Instance

H2

H3

MIP

VNS

DCS

HDCS

MIP

VNS

DCS

HDCS

MIP

VNS

DCS

HDCS

8×2
8×3
8×4
10×2
10×3
10×4
12×2
12×3
12×4
14×2
14×3
14×4

5.22
2.01
0.33
1188.281
209.54
14.91
3559.689
2496.766
1873.835
3600.0010
3600.0010
2576.576

1.06
0.99
0.96
1.94
1.80
1.70
3.00
2.87
2.78
4.73
4.59
4.10

0.52
0.58
0.49
0.83
0.81
0.78
1.01
0.95
0.89
1.01
0.85
0.96

1.56
1.50
1.38
2.62
2.60
2.28
4.34
4.02
3.91
6.74
6.88
5.69

30.91
1.17
0.20
2140.394
423.81
6.29
3259.209
3101.888
1160.383
3600.0010
2791.577
2755.647

1.08
1.00
0.93
1.95
1.77
1.62
2.80
2.93
2.76
4.53
4.62
4.24

0.55
0.51
0.49
0.62
0.64
0.59
0.70
0.75
0.72
0.90
0.82
0.95

1.56
1.49
1.37
2.60
2.46
2.25
4.15
4.20
3.61
6.82
6.61
6.24

16.12
2.83
0.31
1590.562
97.56
48.63
3415.899
2328.345
1826.884
3258.919
2862.137
3039.727

1.09
0.99
0.95
1.97
1.75
1.70
3.02
2.86
2.69
4.81
4.41
4.26

0.63
0.56
0.41
0.67
0.74
0.71
0.80
0.83
0.80
1.00
0.91
0.86

1.55
1.46
1.38
2.72
2.44
2.32
4.58
4.00
3.80
6.59
6.28
6.23

Average

1593.93

2.54

0.81

3.63

1605.95

2.52

0.69

3.61

1540.66

2.54

0.74

3.61

∗` indicates that ` out of 10 instances can not be solved optimally by the CPLEX within the 3600 seconds time limit.

Figure 8. Average RPDA Means plot and the Tukey HSD intervals (at the 95% confidence level) for the tested
algorithms (small instances)

In order to validate the statistical significance of the observed difference in solution
quality by different algorithms, we apply again the one-way ANOVA as in the previous
scenario. Figure 9 shows the means plot with the Tukey HSD intervals at 95% confidence
level for each algorithm under different deteriorating intervals. It is clearly seen that, the
HDCS is statistically better than the other two algorithms, as the HDCS consistently
and statistically possesses lower Tukey interval compared to the other two algorithms for
each of the three deteriorating intervals H1 , H2 and H3 .
In Table 4, from the row of “Average”, it can be concluded that the average time taken

17

September 6, 2013

Engineering Optimization

PMSDST˙HDCS˙EO˙92˙arxiv

by the DCS algorithm is the shortest, while the average time needed for the HDCS algorithm is the longest and the time consumed by the VNS lies between the two algorithms.
This is due to the fact that, there are more than one elite solution being improved by the
variable neighborhood descent in the HDCS. But the computational time of the HDCS
is acceptable for obtaining a good quality solution to the underlying problem.
Finally, the effect of deteriorating intervals on the performance of the proposed algorithm is analyzed. As shown in Figure 10, the performance of the HDCS somewhat
depends on a deteriorating interval, but it is not too much significant. When the deteriorating dates are generated from H3 , the mean RPD values delivered by the HDCS is
the smallest.
Table 4. Computational results of large instances for the algorithms
VNS

DCS

HDCS

Instance

RPDB

RPDA

RPDW

Time(s)

RPDB

RPDA

RPDW

Time(s)

RPDB

RPDA

RPDW

Time(s)

30×4
30×6
30×8
40×4
40×6
40×8
H1
50×4
50×6
50×8
60×4
60×6
60×8
Average

1.49
0.24
0.89
0.86
1.43
2.48
1.99
3.62
3.52
6.85
2.63
3.05
2.42

8.50
6.98
8.05
10.31
13.75
10.42
23.64
19.39
16.27
18.67
12.71
11.21
13.33

19.09
15.99
16.09
24.61
25.76
19.29
29.63
42.70
30.23
30.55
21.50
19.78
24.60

41.63
39.69
39.12
95.31
91.41
78.66
164.51
155.02
151.72
349.18
307.53
316.57
152.53

49.61
43.71
45.99
65.89
75.04
50.46
94.65
53.94
76.00
71.67
63.89
57.48
62.36

49.63
43.99
46.34
67.90
80.10
51.35
96.52
58.40
76.25
76.47
63.89
57.48
64.03

49.63
44.06
46.50
68.64
81.40
51.57
97.01
66.08
76.31
79.71
63.89
57.48
65.19

5.16
5.27
6.48
7.40
7.92
7.32
8.76
8.68
8.44
10.48
10.44
10.20
8.05

0.90
1.57
1.60
1.98
3.79
1.06
1.42
3.75
1.79
1.64
2.53
0.60
1.89

4.65
3.58
4.06
6.09
9.16
5.43
7.40
11.96
9.22
10.92
7.67
6.37
7.21

8.50
6.22
6.58
10.23
14.80
9.75
15.21
20.18
16.07
18.80
13.38
12.23
12.66

50.38
51.10
46.58
112.20
110.94
102.04
221.95
217.16
183.92
382.93
315.90
337.98
177.76

30×4
30×6
30×8
40×4
40×6
40×8
H2
50×4
50×6
50×8
60×4
60×6
60×8
Average

1.70
0.94
1.46
0.81
0.87
0.71
2.83
1.96
1.99
3.04
4.55
3.81
2.06

8.81
6.12
7.58
10.33
10.16
8.66
17.59
11.16
15.06
10.77
15.73
8.96
10.91

16.44
13.54
15.95
22.68
23.08
22.33
38.13
21.57
27.47
18.16
28.03
19.10
22.21

43.43
41.03
36.84
97.64
90.96
82.68
156.74
168.50
166.81
310.91
288.22
262.49
145.52

35.19
25.84
28.31
38.18
37.43
49.37
60.14
26.82
38.17
32.29
44.30
20.31
36.36

35.64
27.36
29.78
38.32
38.43
50.00
61.22
26.82
39.24
32.29
44.30
20.31
36.98

35.79
28.15
30.52
38.37
39.01
50.15
61.15
26.82
40.06
32.29
44.30
20.31
37.24

5.45
5.44
6.68
7.00
6.88
6.64
9.00
8.60
8.84
10.44
10.36
10.24
7.96

0.96
0.40
0.60
2.49
2.05
1.66
1.34
1.97
1.96
0.97
1.07
1.61
1.42

5.20
2.51
2.39
7.81
5.67
5.66
7.95
6.39
7.77
6.79
6.62
5.04
5.82

9.70
4.50
3.88
12.19
10.00
9.81
14.98
11.90
13.31
12.50
11.41
8.68
10.24

46.19
45.30
42.12
104.59
91.51
98.96
184.93
179.63
168.26
316.43
298.84
270.85
153.97

30×4
30×6
30×8
40×4
40×6
40×8
H3
50×4
50×6
50×8
60×4
60×6
60×8
Average

2.23
1.27
1.69
1.92
0.92
0.63
3.12
1.01
0.72
4.31
4.05
5.32
2.27

8.59
9.44
7.24
12.24
9.67
8.63
16.66
10.53
13.99
12.80
15.74
17.55
11.92

16.09
17.66
15.42
22.62
19.79
18.36
34.57
22.87
29.39
20.14
28.93
30.32
23.01

42.96
41.98
36.66
97.18
88.09
84.11
165.67
152.81
148.59
312.29
310.69
268.94
145.83

53.71
51.20
50.18
57.73
47.61
66.75
59.92
45.53
81.57
69.43
62.67
73.55
59.99

54.57
52.43
51.67
57.95
48.74
67.53
62.83
45.72
82.41
69.52
62.67
73.55
60.80

54.94
53.04
52.44
58.01
49.02
67.75
67.96
45.77
82.62
69.54
62.67
73.55
61.44

5.82
6.12
6.16
6.96
6.92
6.72
8.68
8.52
8.52
10.48
10.36
10.20
7.95

0.25
2.08
0.35
0.99
1.88
1.20
3.66
1.68
3.48
1.40
1.83
2.06
1.74

3.50
4.31
1.86
5.89
5.67
5.62
10.39
6.67
10.43
8.26
7.78
10.88
6.77

5.76
7.04
3.15
11.00
10.13
10.45
17.57
11.08
18.53
13.26
12.62
18.89
11.62

50.39
46.47
41.30
115.97
103.95
98.52
199.58
189.81
190.57
335.21
332.11
299.49
166.95

5.

Conclusions

In this article, an identical parallel machine scheduling problem with step-deteriorating
jobs and sequence-dependent setup times was considered. The objective of this problem is to find a schedule for minimizing the total tardiness. To solve this problem, a
mixed-integer programming model was presented to deliver the exact solutions for relatively small-sized instances. Due to the intractability of the problem, a hybrid discrete
cuckoo search algorithm was proposed to obtain a near-optimal solution. In the proposed
algorithm, the operators of the standard CS algorithm are redefined by some discrete

18

September 6, 2013

Engineering Optimization

PMSDST˙HDCS˙EO˙92˙arxiv

Figure 9. Average RPDA means plot and the Tukey HSD intervals (at the 95% confidence level) for the tested
algorithms (large instances)

operations. To improve the quality of the solution, the variable neighborhood descent
as local search approach is used to refine a set of chosen elite solutions. Moreover, the
MBHG heuristic was incorporated into the generation of initial solutions. A set of small
to large instances were produced to test the performance of the proposed algorithm.
The computational results show that the hybrid discrete cuckoo search algorithm in general provides good average results for test problems, especially in large-sized instances.
Indeed, the HDCS algorithm outperforms the CPLEX solver, the VNS and the DCS
algorithms in terms of the average RPD values, shown in the row ‘Average’ of tables 2
and 4. Further studies may concern on applying the proposed algorithm to solve other
scheduling problems, such as unrelated parallel machine scheduling problem, the flow
shop problem and their stochastic versions. In addition, some multi-objective scheduling
problems may be considered in future work.

Acknowledgments
This work was partially supported by the National Natural Science Foundation of China
(No.51175442 and 51205328) and the Fundamental Research Funds for the Central Universities (No. 2010ZT03; SWJTU09CX022).

References
Alidaee, B., and N.K. Womer. 1999. “Scheduling with time dependent processing times: Review
and extensions.” Journal of the Operational Research Society 50 (7): 711–720.

19

September 6, 2013

Engineering Optimization

PMSDST˙HDCS˙EO˙92˙arxiv

Figure 10. The effect of different deteriorating intervals on the performance of the HDCS

Allahverdi, A., J.N.D. Gupta, and T. Aldowaisan. 1999. “A review of scheduling research involving
setup considerations.” Omega 27 (2): 219–239.
Allahverdi, A., C.T. Ng, T.C.E. Cheng, and M.Y. Kovalyov. 2008. “A survey of scheduling problems with setup times or costs.” European Journal of Operational Research 187 (3): 985–1032.
Armentano, V.A., and M.F. De Franca Filho. 2007. “Minimizing total tardiness in parallel machine scheduling with setup times: An adaptive memory-based GRASP approach.” European
Journal of Operational Research 183 (1): 100–114.
Bahalke, U., A.M. Yolmeh, and K. Shahanaghi. 2010. “Meta-heuristics to solve single-machine
scheduling problem with sequence-dependent setup time and deteriorating jobs.” International
Journal of Advanced Manufacturing Technology 50 (5-8): 749–759.
Bank, M., S.M.T. Fatemi Ghomi, F. Jolai, and J. Behnamian. 2012. “Application of particle
swarm optimization and simulated annealing algorithms in flow shop scheduling problem under
linear deterioration.” Advances in Engineering Software 47 (1): 1–6.
Biskup, D., J. Herrmann, and J.N.D. Gupta. 2008. “Scheduling identical parallel machines to
minimize total tardiness.” International Journal of Production Economics 115 (1): 134–142.
Browne, S., and U. Yechiali. 1990. “Scheduling Deteriorating Jobs on a Single Processor.” Operations Research 38 (3): 495–498.
Burnwal, S., and S. Deb. 2013. “Scheduling optimization of flexible manufacturing system using
cuckoo search-based approach.” International Journal of Advanced Manufacturing Technology
64 (5-8): 951–959.
Cheng, T.C.E., and Q. Ding. 2001. “Single machine scheduling with step-deteriorating processing
times.” European Journal of Operational Research 134 (3): 623–630.
Cheng, T.C.E., Q. Ding, M.Y. Kovalyov, A. Bachman, and A. Janiak. 2003. “Scheduling jobs
with piecewise linear decreasing processing times.” Naval Research Logistics 50 (6): 531–554.
Cheng, T.C.E, Q. Ding, and B.M.T Lin. 2004. “A concise survey of scheduling with timedependent processing times.” European Journal of Operational Research 152 (1): 1–13.
Cheng, T.C.E., C.J. Hsu, Y.C. Huang, and W.C. Lee. 2011. “Single-machine scheduling with dete-

20

September 6, 2013

Engineering Optimization

PMSDST˙HDCS˙EO˙92˙arxiv

riorating jobs and setup times to minimize the maximum tardiness.” Computers & Operations
Research 38 (12): 1760–1765.
Cheng, T.C.E., W.C. Lee, and C.C. Wu. 2011. “Single-machine scheduling with deteriorating jobs
and past-sequence-dependent setup times.” Applied Mathematical Modelling 35 (4): 1861–1867.
Cheng, W.M., P. Guo, Z.Q. Zhang, M. Zeng, and J. Liang. 2012. “Variable Neighborhood Search
for Parallel Machines Scheduling Problem with Step Deteriorating Jobs.” Mathematical Problems in Engineering 2012: 928312.
Civicioglu, P., and E. Besdok. 2013. “A conceptual comparison of the Cuckoo-search, particle
swarm optimization, differential evolution and artificial bee colony algorithms.” Artificial Intelligence Review 39 (4): 315–346.
Davis, L. 1985. “Applying adaptive algorithms to epistatic domains.” In Proceedings of the 9th
international joint conference on Artificial intelligence, Los Angeles, California. 162–164. San
Francisco, CA, USA.
Farahani, M.H., and L. Hosseini. 2013. “Minimizing cycle time in single machine scheduling
with start time-dependent processing times.” International Journal of Advanced Manufacturing
Technology 64 (9-12): 1479–1486.
Gandomi, A.H., X.S. Yang, and A.H. Alavi. 2013. “Cuckoo search algorithm: A metaheuristic
approach to solve structural optimization problems.” Engineering with Computers 29 (1): 17–
35.
Gawiejnowicz, S. 2008. Time-Dependent Scheduling. Berlin Heidelberg: Springer.
Guo, P., W.M. Cheng, and Y. Wang. 2013. “A general variable neighborhood search for singlemachine total tardiness scheduling problem with step-deteriorating jobs.” Appear in Journal
of Industrial and Management Optimization arXiv:1301.7134.
Gupta, J.N.D., and S.K. Gupta. 1988. “Single facility scheduling with nonlinear processing times.”
Computers & Industrial Engineering 14 (4): 387–393.
He, C.C., C.C. Wu, and W.C. Lee. 2009. “Branch-and-bound and weight-combination search
algorithms for the total completion time problem with step-deteriorating jobs.” Journal of the
operational research society 60 (12): 1759–1766.
Huang, X., and M.Z. Wang. 2011. “Parallel identical machines scheduling with deteriorating jobs
and total absolute differences penalties.” Applied Mathematical Modelling 35 (3): 1349–1353.
Jafari, A., and G. Moslehi. 2012. “Scheduling linear deteriorating jobs to minimize the number
of tardy jobs.” Journal of Global Optimization 54 (2): 389–404.
Jeng, A.A.K., and B.M.T. Lin. 2004. “Makespan minimization in single-machine scheduling with
step-deterioration of processing times.” Journal of the Operational Research Society 55 (3):
247–256.
Jeng, A.A.K., and B.M.T. Lin. 2005. “Minimizing the total completion time in single-machine
scheduling with step-deteriorating jobs.” Computers & Operations Research 32 (3): 521–536.
Ji, M., and T.C.E. Cheng. 2007. “An FPTAS for scheduling jobs with piecewise linear decreasing
processing times to minimize makespan.” Information Processing Letters 102 (2-3): 41–47.
Kubiak, W., and S. Van de Velde. 1998. “Scheduling deteriorating jobs to minimize makespan.”
Naval Research Logistics 45 (5): 511–523.
Kunnathur, A.S., and S.K. Gupta. 1990. “Minimizing the makespan with late start penalties
added to processing times in a single facility scheduling problem.” European Journal of Operational Research 47 (1): 56–64.
Lai, P.J., W.C. Lee, and H.H. Chen. 2011. “Scheduling with deteriorating jobs and past-sequencedependent setup times.” International Journal of Advanced Manufacturing Technology 54 (5-8):
737–741.
Layegh, J., F. Jolai, and M.S. Amalnik. 2009. “A memetic algorithm for minimizing the total
weighted completion time on a single machine under step-deterioration.” Advances in Engineering Software 40 (10): 1074–1077.
Lee, W.C., J.B. Lin, and Y.R. Shiau. 2011. “Deteriorating job scheduling to minimize the number
of late jobs with setup times.” Computers & Industrial Engineering 61 (3): 782–787.
Leung, J.Y.T., C.T. Ng, and T.C.E. Cheng. 2008. “Minimizing sum of completion times for
batch scheduling of jobs with deteriorating processing times.” European Journal of Operational
Research 187 (3): 1090–1099.
Li, X.T., and M.H. Yin. 2013. “A hybrid cuckoo search via Lévy flights for the permutation flow

21

September 6, 2013

Engineering Optimization

PMSDST˙HDCS˙EO˙92˙arxiv

shop scheduling problem.” International Journal of Production Research 51 (16): 4732–4754.
Mor, B., and G. Mosheiov. 2012. “Batch scheduling with step-deteriorating processing times to
minimize flowtime.” Naval Research Logistics 59 (8): 587–600.
Mosheiov, G. 1995. “Scheduling jobs with step-deterioration; Minimizing makespan on a singleand multi-machine.” Computers & Industrial Engineering 28 (4): 869–879.
Moslehi, G., and A. Jafari. 2010. “Minimizing the number of tardy jobs under piecewise-linear
deterioration.” Computers & Industrial Engineering 59 (4): 573–584.
Ouaarab, A., B. Ahiod, and X.S. Yang. Online First 2013. “Discrete cuckoo search algorithm for
the travelling salesman problem.” Neural Computing and Applications 1–11.
Radhakrishnan, S., and J.A. Ventura. 2000. “Simulated annealing for parallel machine scheduling
with earliness-tardiness penalties and sequence-dependent set-up times.” International Journal
of Production Research 38 (10): 2233–2252.
Ruiz, R., C. Maroto, and J. Alcaraz. 2006. “Two new robust genetic algorithms for the flowshop
scheduling problem.” Omega 34 (5): 461–476.
Sundararaghavan, P.S., and A.S. Kunnathur. 1994. “Single machine scheduling with start time
dependent processing times: Some solvable cases.” European Journal of Operational Research
78 (3): 394–403.
Toksari, M.D., and E. Guener. 2008. “Minimizing the earliness/tardiness costs on parallel machine
with learning effects and deteriorating jobs: a mixed nonlinear integer programming approach.”
International Journal of Advanced Manufacturing Technology 38 (7-8): 801–808.
Valian, E., S. Tavakoli, S. Mohanna, and A. Haghi. 2013. “Improved cuckoo search for reliability
optimization problems.” Computers & Industrial Engineering 64 (1): 459–468.
Walton, S., O. Hassan, K. Morgan, and M.R. Brown. 2011. “Modified cuckoo search: A new
gradient free optimisation algorithm.” Chaos, Solitons & Fractals 44 (9): 710–718.
Wang, J.B., Y. Jiang, and G. Wang. 2009. “Single-machine scheduling with past-sequencedependent setup times and effects of deterioration and learning.” International Journal of
Advanced Manufacturing Technology 41 (11-12): 1221–1226.
Wu, C.C., Y.R. Shiau, L.H. Lee, and W.C. Lee. 2009. “Scheduling deteriorating jobs to minimize the makespan on a single machine.” International Journal of Advanced Manufacturing
Technology 44 (11-12): 1230–1236.
Yang, X.S., and S. Deb. 2009. “Cuckoo search via Lévy flights.” In Proceedings of the world
congress on nature and biologically inspired computing, 210–214. IEEE. Coimbatore, India:
IEEE Publications.
Yang, X.S., and S. Deb. 2010. “Engineering optimisation by cuckoo search.” International Journal
of Mathematical Modelling and Numerical Optimisation 1 (4): 330–343.
Yildiz, A.R. 2013. “Cuckoo search algorithm for the selection of optimal machining parameters
in milling operations.” International Journal of Advanced Manufacturing Technology 64 (1-4):
55–61.
Yin, Y.Q., T.C.E. Cheng, J.Y. Xu, S.R. Cheng, and C.C. Wu. 2013. “Single-machine scheduling
with past-sequence-dependent delivery times and a linear deterioration.” Journal of Industrial
and Management Optimization 9 (2): 323–339.
Zhao, C.L., and H.Y. Tang. 2010. “Single machine scheduling with past-sequence-dependent setup
times and deteriorating jobs.” Computers & Industrial Engineering 59 (4): 663–666.

22

