# A PSO-Based Hybrid Metaheuristic for Permutation Flowshop Scheduling Problems

Authors: Le Zhang; Jinnan Wu
Journal: The Scientific World Journal
 (2014)

---

## Abstract

This paper investigates the permutation flowshop scheduling problem (PFSP) with the objectives of minimizing the makespan and the total flowtime and proposes a hybrid metaheuristic based on the particle swarm optimization (PSO). To enhance the exploration ability of the hybrid metaheuristic, a simulated annealing hybrid with a stochastic variable neighborhood search is incorporated. To improve the search diversification of the hybrid metaheuristic, a solution replacement strategy based on the pathrelinking is presented to replace the particles that have been trapped in local optimum. Computational results on benchmark instances show that the proposed PSO-based hybrid metaheuristic is competitive with other powerful metaheuristics in the literature.


---

## Body

## 1. Introduction

Due to the strong industrial background, the permutation flowshop scheduling problem (PFSP) has attracted considerable attention from researchers all over the world. In this problem, a set of jobsJ
=
{
1,2
,
…
,
n
} needs to be processed through a set of machines M
=
{
1,2
,
…
,
m
}. Each job i
∈
J should be processed through these m machines with the same machine order, that is, starting from machine 1 and finishing on the last machine m. The processing time of each job i
∈
J on machine j
∈
M
(
p
i
j
) is nonnegative and known before scheduling. It is assumed that all the jobs are available before processing and once started the processing cannot be interrupted. It is required that each job can only be processed by only one machine at any time, and at the same time each machine cannot process more than one job. The processing of a job cannot start on the next machine j
+
1 until this job has been completed on the current machine j and machine j
+
1 is idle. The objective is to determine the sequence of these n jobs so that a certain performance measure can be optimized. The most commonly studied performance measures are the minimization of makespan (
C
max
) and the minimization of total flowtime (
TFT
). Let π
=
(
π
(
1
)
,
…
,
π
(
n
)
) denote a permutation of the n jobs, in which π
(
k
) represents the job arranged at the kth position, and the completion time of each job π
(
k
) on each machine j can be calculated as follows:
(1)
C
π
(
1
)
,
1
=
p
π
(
1
)
,
1
,
C
π
(
1
)
,
j
=
C
π
(
1
)
,
j
-
1
+
p
π
(
1
)
,
j
,
j
=
2,3
,
…
,
m
,
C
π
(
k
)
,
1
=
C
π
(
k
-
1
)
,
1
+
p
π
(
k
)
,
1
,
k
=
2,3
,
…
,
n
,
C
π
(
k
)
,
j
=
max
{
C
π
(
k
)
,
j
-
1
,
C
π
(
k
-
1
)
,
j
}
+
p
π
(
k
)
,
j
,
k
=
2,3
k
=
k
=
2,3
,
…
,
n
;
j
=
2,3
,
…
,
m
.
Then the makespan can be defined as C
max
(
π
)
=
C
π
(
n
)
,
m, and the total flowtime can be defined as the sum of completion times of all jobs TFT
(
π
)
=
∑
k
=
1
n
C
π
(
k
)
,
m.Since the first introduction of the PFSP [1], considerable attention of researchers has been paid to this problem and many kinds of algorithms have been proposed in the literature. According to the comprehensive review of Ruiz and Maroto [2] and Framinan et al. [3] for PFSP, these solution methods can be classified into three categories: exact methods, heuristics methods, and metaheuristic methods.Since it has been proven that the PFSP with makespan minimization is NP-complete in the strong sense whenm
≥
3 and the PFSP with total flowtime minimization is NP-complete in the strong sense when m
≥
2 [4], few exact methods have been proposed for PFSP in the literature due to their unacceptable computation time. These exact methods include the mixed integer linear programming method [5] and the branch and bound algorithms ([6–8] for the makespan minimization and [9–11] for the total flowtime minimization). However, these exact methods are feasible for only small size problems because they cannot solve large size problems in reasonable computation time.Heuristic methods can be classified into two categories: constructive heuristics and improvement heuristics. Constructive methods start from an empty solution and try to build a feasible solution in a short time. Johnson’s algorithm [1] is the earliest known heuristic for PFSP, which can obtain optimal solutions for m
=
2. Campbell et al. [12] proposed the CDS heuristic, and Koulamas [13] proposed a two-phase heuristic for PFSP, which were both extensions of Johnson’s algorithm. Palmer [14] proposed a slope index heuristic for PFSP, and then Gupta [15] and Hundal and Rajgopal [16] extended Palmer’s heuristic and proposed two simple heuristics. Nawaz et al. [17] proposed a so-called NEH heuristic based on the idea that jobs with high total processing times on all machines should be scheduled as early as possible, and this NEH heuristic is regarded as the best heuristic for PFSP with makespan minimization (Ruiz and Maroto [2], Taillard [18]). Recently, complex heuristics have been proposed for PFSP, for example, Liu and Reeves [19] and Framinan and Leisten [20]. As far as the solution quality is concerned, the FL heuristic proposed by Framinan and Leisten [20] was the best one among simple heuristics (Framinan et al. [3]). Contrary to constructive heuristics, improvement heuristics start from an existing initial solution and try to improve it by a given procedure, for example, a local search. Fan and Winley [21] proposed a heuristic named the intelligent heuristic search algorithm for the PFSP. Suliman [22] proposed a two-phase improvement heuristic, in which an initial solution is generated by the CDS heuristic in the first phase and then improved by a pair exchange neighborhood search in the second phase. Framinan et al. [23] proposed a new efficient heuristic for the PFSP with no-wait constraint.Metaheuristics are high-level strategies that combine and guide other heuristics in the hope of obtaining a more efficient or more robust procedure so that better solutions can be found. The main procedure of metaheuristics generally starts from an initial solution or a set of solutions generated by heuristics and iterates to improve the initial solution or solutions until a given stopping criterion is reached. The metaheuristics proposed for PFSP are mainly the genetic algorithm (Chang et al. [24], Ruiz et al. [25]), the simulated annealing (Hooda and Dhingra [26], Nouri et al. [27]), the tabu search (Gao et al. [28]), the ant colony algorithm (Rajendran and Ziegler [29]), the iterated greedy algorithm (Ruiz and Stützle [30]), and the particle swarm optimization (PSO) (Tasgetiren et al. [31], Wang and Tang [32]). These metaheuristics use the benchmark problems proposed by Taillard [33] to evaluate their performance. The ant colony algorithms proposed by Rajendran and Ziegler [29], named M-MMAS and PACO, obtained much better solutions than constructive heuristics of Framinan and Leisten [20]. The iterated greedy algorithm proposed by Ruiz and Stützle [30] improved the best known results for some instances of PFSP with makespan minimization. The particle swarm optimization (PSO) named PSOVNS, which incorporates variable neighborhood search (VNS) into PSO, proposed by Tasgetiren et al. [31] improved 57 out of 90 best known solutions reported by Framinan and Leisten [20] and Rajendran and Ziegler [29] for the total flowtime criterion.In this paper, we propose an improved PSO for the PFSP. To enhance the exploration ability of PSO, the path relinking and the hybrid simulated annealing with stochastic VNS are incorporated. To improve the search diversification of PSO, a population update method is applied to replace the nonpromising particles. The rest of this paper is organized as follows. Section2 is devoted to describe the proposed PSO algorithm. The computational results on benchmark problems are presented in Section 3. Finally, Section 4 concludes the paper.


## 2. PSO Algorithm for PFSP

### 2.1. Brief Introduction of PSO

PSO algorithm is a population based metaheuristic method introduced by Kennedy and Eberhart [34, 35] based on the social behavior of bird flocking and fish schooling, as well as the means of information exchange between individuals, to solve optimization problems. In the PSO, a swarm consists of m particles and these particles fly around in an n-dimensional search space. The solution of a problem is represented by the position of a particle; that is, the ith particle at the tth generation is denoted as X
i
t
=
[
x
i
1
t
,
x
i
2
t
,
…
,
x
i
n
t
]. At each generation, the flight of each particle is determined by three factors: the inertia of itself, the best position found by itself (
p
best
), and the best position found by the whole swarm (
g
best
). Generally, p
best and g
best are represented as P
i
t
=
[
p
i
1
t
,
p
i
2
t
,
…
,
p
i
n
t
] and G
t
=
[
g
1
t
,
g
2
t
,
…
,
g
n
t
], respectively. Then the velocity of the particle V
i
t
=
[
v
i
1
t
,
v
i
2
t
,
…
,
v
i
n
t
] for the next generation can be obtained from the following equation:
(2)
v
i
j
t
+
1
=
w
·
v
i
j
t
+
c
1
r
1
·
(
p
i
j
t
-
x
i
j
t
)
+
c
2
r
2
·
(
g
j
t
-
x
i
j
t
)
,
x
i
j
t
+
1
=
x
i
j
t
+
v
i
j
t
+
1
,
where w is called the inertiaparameter, c
1 and c
2 are the cognitiveand social parameters, and r
1, r
2 are random numbers between (
0,1
). Based on the above equations, the particle can fly through search space toward p
best and g
best in a navigated way while still exploring new areas by the stochastic mechanism to escape from local optima.


### 2.2. Solution Representation

Since the PSO operates in the continuous space, a job is represented by a dimension of a particle and then then jobs can be denoted as a particle X
i
t
=
[
x
i
1
t
,
x
i
2
t
,
…
,
x
i
n
t
] in the continuous space. Due to the continuous characters of the position values of particles in the PSO, the smallest position value (SPV) rule proposed by Tasgetiren et al. [31] is adopted to transform a particle with continuous position values into a job permutation. A simple example is provided in Table 1 to show the mechanism of the SPV rule. In this instance (n
=
9), the smallest position value is x
i
2
t
=
-
1.75, so job 2 is assigned to the first position of the job permutation according to the SPV rule; then job 9 is assigned to the second position of the job permutation because it has the second smallest position value x
i
9
t
=
-
1.21. With the same way, other jobs are assigned in their corresponding position of the job permutation according to their position values. Thus, based on the SPV rule, the job permutation is obtained; that is, π
i
t
=
(
2,9
,
4,3
,
5,1
,
8,7
,
6
).Table 1
Solution representation and the corresponding job permutation using SPV rule.
Dimensionj
1
2
3
4
5
6
7
8
9
x
i
j
t
0.35
−1.75
−0.02
−0.21
0.02
1.20
1.03
0.67
−1.21
Job,π
i
j
t
6
1
4
3
5
9
8
7
2


### 2.3. Population Initialization

The population withn
pop solutions is initialized with random solutions according to x
i
j
0
=
x
min
+
rand
×
(
x
max
-
x
min
), where rand is a uniform random number in [
0
,
1
], x
min
=
-
4.0, and x
max
=
4.0. Also, we generate the corresponding velocity of each particle by a similar way: v
i
j
0
=
v
min
+
rand
×
(
v
max
-
v
min
), where v
min
=
-
1.0 and v
max
=
1.0. In addition, another solution generated by the NEH heuristic [18] is added to the initial population and replaces a random selected solution so as to ensure the quality of initial population.


### 2.4. Hybrid Method of Simulated Annealing and Stochastic VNS

In the PSOVNS proposed by Tasgetiren et al. [31], a stochastic VNS, which itself is a variant of VNS (Hansen and Mladenović [36]), is developed as the local search. For a given discrete job permutation π
t, let w and z denote two different random integer numbers generated in [
1
,
n
], and then the two stochastic neighborhoods moves used in the stochastic VNS to generate a neighbor solution π
t
′ are (
1
)  π
t
′
=
insert
(
π
t
,
w
,
z
): remove the job at the wth position and insert it in the zth position; and (
2
)  π
t
′
=
interchange
(
π
t
,
w
,
z
): interchange two jobs arranged at the wth position and the zth position. After a job permutation is changed according to a local search operator such as insert or interchange, the position value of each dimension is adjusted correspondingly to guarantee that the permutation that resulted by the SPV rule for new position values is the same as the permutation that resulted by the local search operator. For example, Table 2 shows the interchange move applied to two jobs 3 and 7, and the corresponding position value changes. It is clear that the interchange of jobs 3 and 7 is corresponding to the interchange of position values −1.02 and 0.23. The position value adjustment for the insert move is similar.Table 2
Interchange move on the job permutation and the corresponding position value adjustment.
Dimensionj
1
2
3
4
5
6
7
8
Before interchange
x
i
j
t
0.54
−0.75
−1.02
−0.41
0.92
−1.20
0.23
0.12
Job,π
i
j
t
6
3
2
4
8
7
1
5
After interchange
x
i
j
t
0.54
−0.75
0.23
−0.41
0.92
–1.20
−1.02
0.12
Job,π
i
j
t
6
7
2
4
8
3
1
5
The bold and italic values are used to show the interchange move applied to jobs 3 and 7.To further improve the exploration ability of the local search, we incorporate the solution acceptance scheme of simulated annealing into the stochastic VNS and thus obtain a hybrid method of simulated annealing and stochastic VNS (denoted as SA_VNS). To reduce the computation time and make the search process focus on the intensification phase, we use a decreasing acceptance threshold to act as the cooling procedure of simulated annealing. The procedure of the proposed SA_VNS algorithm is illustrated in Algorithm1.Algorithm 1:The main procedure of theSA_VNS.
Begin:
Initialization:
Letπ
0 be the input initial solution. Set π
=
π
0, the acceptance threshold T
=
0.05, and
o
u
t
e
r
l
o
o
p
_
c
o
u
n
t
e
r
=
0.
while (o
u
t
e
r
l
o
o
p
_
c
o
u
n
t
e
r
<
n
/
5) do
(1) Generate a random numberr in [0, 1], and two random integer numbers w and z.
If r > 0.5, generate π
′ = insert (π
,
w
,
z); otherwise generate π
′ = interchange (π
,
w
,
z).
(2) Seti
n
n
e
r
_
c
o
u
n
t
e
r
=
0.
(3)while  (
i
n
n
e
r
_
c
o
u
n
t
e
r
<
n
×
(
n
-
1
)
)  do
Set k
=
1.
while (k
≤
2) do
(1) Generate two random integer numbers w and z.
(2) If k
=
1, generate π
′
′ = i
n
s
e
r
t
(
π
′
,
w
,
z
).
(3) If k
=
2, generate π
′
′ = i
n
t
e
r
c
h
a
n
g
e
(
π
′
,
w
,
z
).
(4) If f
(
π
′
′
)
<
f
(
π
′
)
,
then set
π
′
=
π
′
′ and k
=
1; otherwise set k
=
k
+
1.
end while
Set i
n
n
e
r
_
c
o
u
n
t
e
r
=
i
n
n
e
r
_
c
o
u
n
t
e
r
+
1.
end while
(4) Iff
(
π
′
)
<
f
(
π
0
)
,
set
π
0
=
π
′
and
π
=
π
′, and then go to step 6; otherwise go to step 5.
(5) Iff
(
π
′
)
≥
f
(
π
0
)
but
(
f
(
π
′
)
-
f
(
π
0
)
)
/
f
(
π
0
)
≤
T
,
set
π
=
π
′; otherwise, generate a random
numberr in [0, 1], and then set π
=
π
0 if r
>
0.5 and π
=
π
′ if r
≤
0.5.
(6) SetT
=
T
×
0.95 and o
u
t
e
r
l
o
o
p
_
c
o
u
n
t
e
r
=
o
u
t
e
r
l
o
o
p
_
c
o
u
n
t
e
r
+
1.
end while
Report the improved solution  π
0.
EndIn the PSOVNS proposed by Tasgetiren et al. [31], the stochastic VNS is applied on the global best particle found at each iteration. However, a drawback of such application is that the starting point of the stochastic VNS may be the same solution if the global best particle cannot be improved for a number of consecutive iterations, and consequently the exploration ability of the PSO may be decreased. Therefore, for a given population at iteration t, we propose to use the following strategy.Step 1.
Apply the SA_VNS on the promising particles satisfying(
f
(
π
t
)
-
f
(
g
best
)
)
/
f
(
g
best
)
≤
0.02, in which f
(
π
t
) is the objective value of particle π
t.Step 2.
Update the global best particle. If a new global best particle is found, then further improve it using the SA_VNS.


### 2.5. Population Update Method

It is well known that the advantage of PSO is that it has a high convergence speed. However, this advantage may become the disadvantage for complex scheduling problems because the scheduling problems generally have many local optimal regions in the search space. That is, for the PSO applied to PFSP, some particles may always fly around a local region and thus are trapped in local optimum. Therefore, we propose a solution replacement strategy based on the pathrelinking [37] to remove these solutions with new solutions with good quality.In our algorithm, a particle is viewed as being trapped in local optimum if its personal best solutionp
best has not been improved for a number of consecutive generations (i.e., 20). For these particles, we give them the last chance to stay in the population by applying the path relinking algorithm on it to check if its personal best p
best can be improved. If so, this particle can remain in the population; otherwise, we replace this particle with a new random particle.The path relinking is originally proposed by Glover et al. [37] to generate new solutions by exploring a path that connects an initial solution and a guiding solution. In this path, moves are selected that introduce attributes contained in the guiding solution. To present the path relinking algorithm, we define the distance between two particles X
i and X
j as d
(
X
i
,
X
j
)
=
d
(
π
i
,
π
j
)
=
|
h
(
π
i
)
-
h
(
π
j
)
|, where π is the corresponding job permutation (obtained by the SPV rule) of particle X and h
(
π
) is the Hash function value of π
=
(
π
(
1
)
,
…
,
π
(
n
)
) that is calculated as h
(
π
)
=
∑
k
=
1
n
k
×
π
(
k
)
×
π
(
k
). Then the path relinking algorithm can be described as follows.Step 1.
Setπ
i
=
(
π
i
(
1
)
,
π
i
(
2
)
,
…
,
π
i
(
n
)
) as the initial solution. Find the particle that has the largest distance to particle X
i in the current population (denoted as X
j), and set its corresponding job permutation π
j
=
(
π
j
(
1
)
,
π
j
(
2
)
,
…
,
π
j
(
n
)
) as the guiding solution. Set k
=
1 and the local optimum π
opt
=
π
i.Step 2.
Ifπ
i
(
k
)
≠
π
j
(
k
), find the job with index π
j
(
k
) in π
i and swap it with π
i
(
k
) to generate a new job permutation π
i
′. If π
i
′ is better than π
opt, then set π
opt
=
π
i
′.Step 3.
Setk
=
k
+
1. If k
>
n
-
10, stop; otherwise, go to Step 2.It should be noted that the above path relinking stops whenk is larger than n
-
10, because if particle X
j is better than X
i, then X
i will be replaced by X
j if the path relinking stops when k
=
n, which may result in duplicated particles and thus decrease the search diversification.


### 2.6. Complete Procedure of the Proposed PSO

Step 1.
InitializationStep
1
.1. Set initial values for the population size n
pop, the inertia weight, and the cognitive and social parameters. Set t
=
0 and N
O
i
=
0
(
i
=
1,2
,
…
,
n
pop
) for each particle in the population. Create the initial population and the initial velocities for each particle using the method described in Section 2.3.Step
1.2. Generate the job permutation for each particle in the population using the SPV rule, and calculate the objective value of each particle.Step
1.3. Set the personal best of each particle to be the particle itself and the global best to be the best one among the population.Step 2.
Update Particle PositionsStep
2.1. Update iteration counter t
=
t
+
1.Step
2.2. Update inertia weight by w
=
w
×
β.Step
2.3. For each particle, update the velocity and position values according to (2).Step
2.4. Generate the job permutation for each particle in the current population using the SPV rule, and calculate the objective value of each particle.Step 3.
Local Search PhaseStep
3.1. Use the SA_VNS algorithm to improve the promising particles in the current population and then the global best particle found so far according to the adoption strategy of SA_VNS described in Section 2.4.Step
3.2. For each particle in the current population, update its personal best p
best. If p
best of particle i is improved, then set N
O
i
=
0; otherwise set N
O
i
=
N
O
i
+
1.Step 4.
Population Update. For each particle in the current population, use the population update method described in Section 2.5 to update the current population.Step 5.
Stopping Criterion. If t
>
T
max (the maximum iteration number) or the runtime has reached the limit, then stop; otherwise, go to Step 2.


## 2.1. Brief Introduction of PSO

PSO algorithm is a population based metaheuristic method introduced by Kennedy and Eberhart [34, 35] based on the social behavior of bird flocking and fish schooling, as well as the means of information exchange between individuals, to solve optimization problems. In the PSO, a swarm consists of m particles and these particles fly around in an n-dimensional search space. The solution of a problem is represented by the position of a particle; that is, the ith particle at the tth generation is denoted as X
i
t
=
[
x
i
1
t
,
x
i
2
t
,
…
,
x
i
n
t
]. At each generation, the flight of each particle is determined by three factors: the inertia of itself, the best position found by itself (
p
best
), and the best position found by the whole swarm (
g
best
). Generally, p
best and g
best are represented as P
i
t
=
[
p
i
1
t
,
p
i
2
t
,
…
,
p
i
n
t
] and G
t
=
[
g
1
t
,
g
2
t
,
…
,
g
n
t
], respectively. Then the velocity of the particle V
i
t
=
[
v
i
1
t
,
v
i
2
t
,
…
,
v
i
n
t
] for the next generation can be obtained from the following equation:
(2)
v
i
j
t
+
1
=
w
·
v
i
j
t
+
c
1
r
1
·
(
p
i
j
t
-
x
i
j
t
)
+
c
2
r
2
·
(
g
j
t
-
x
i
j
t
)
,
x
i
j
t
+
1
=
x
i
j
t
+
v
i
j
t
+
1
,
where w is called the inertiaparameter, c
1 and c
2 are the cognitiveand social parameters, and r
1, r
2 are random numbers between (
0,1
). Based on the above equations, the particle can fly through search space toward p
best and g
best in a navigated way while still exploring new areas by the stochastic mechanism to escape from local optima.


## 2.2. Solution Representation

Since the PSO operates in the continuous space, a job is represented by a dimension of a particle and then then jobs can be denoted as a particle X
i
t
=
[
x
i
1
t
,
x
i
2
t
,
…
,
x
i
n
t
] in the continuous space. Due to the continuous characters of the position values of particles in the PSO, the smallest position value (SPV) rule proposed by Tasgetiren et al. [31] is adopted to transform a particle with continuous position values into a job permutation. A simple example is provided in Table 1 to show the mechanism of the SPV rule. In this instance (n
=
9), the smallest position value is x
i
2
t
=
-
1.75, so job 2 is assigned to the first position of the job permutation according to the SPV rule; then job 9 is assigned to the second position of the job permutation because it has the second smallest position value x
i
9
t
=
-
1.21. With the same way, other jobs are assigned in their corresponding position of the job permutation according to their position values. Thus, based on the SPV rule, the job permutation is obtained; that is, π
i
t
=
(
2,9
,
4,3
,
5,1
,
8,7
,
6
).Table 1
Solution representation and the corresponding job permutation using SPV rule.
Dimensionj
1
2
3
4
5
6
7
8
9
x
i
j
t
0.35
−1.75
−0.02
−0.21
0.02
1.20
1.03
0.67
−1.21
Job,π
i
j
t
6
1
4
3
5
9
8
7
2


## 2.3. Population Initialization

The population withn
pop solutions is initialized with random solutions according to x
i
j
0
=
x
min
+
rand
×
(
x
max
-
x
min
), where rand is a uniform random number in [
0
,
1
], x
min
=
-
4.0, and x
max
=
4.0. Also, we generate the corresponding velocity of each particle by a similar way: v
i
j
0
=
v
min
+
rand
×
(
v
max
-
v
min
), where v
min
=
-
1.0 and v
max
=
1.0. In addition, another solution generated by the NEH heuristic [18] is added to the initial population and replaces a random selected solution so as to ensure the quality of initial population.


## 2.4. Hybrid Method of Simulated Annealing and Stochastic VNS

In the PSOVNS proposed by Tasgetiren et al. [31], a stochastic VNS, which itself is a variant of VNS (Hansen and Mladenović [36]), is developed as the local search. For a given discrete job permutation π
t, let w and z denote two different random integer numbers generated in [
1
,
n
], and then the two stochastic neighborhoods moves used in the stochastic VNS to generate a neighbor solution π
t
′ are (
1
)  π
t
′
=
insert
(
π
t
,
w
,
z
): remove the job at the wth position and insert it in the zth position; and (
2
)  π
t
′
=
interchange
(
π
t
,
w
,
z
): interchange two jobs arranged at the wth position and the zth position. After a job permutation is changed according to a local search operator such as insert or interchange, the position value of each dimension is adjusted correspondingly to guarantee that the permutation that resulted by the SPV rule for new position values is the same as the permutation that resulted by the local search operator. For example, Table 2 shows the interchange move applied to two jobs 3 and 7, and the corresponding position value changes. It is clear that the interchange of jobs 3 and 7 is corresponding to the interchange of position values −1.02 and 0.23. The position value adjustment for the insert move is similar.Table 2
Interchange move on the job permutation and the corresponding position value adjustment.
Dimensionj
1
2
3
4
5
6
7
8
Before interchange
x
i
j
t
0.54
−0.75
−1.02
−0.41
0.92
−1.20
0.23
0.12
Job,π
i
j
t
6
3
2
4
8
7
1
5
After interchange
x
i
j
t
0.54
−0.75
0.23
−0.41
0.92
–1.20
−1.02
0.12
Job,π
i
j
t
6
7
2
4
8
3
1
5
The bold and italic values are used to show the interchange move applied to jobs 3 and 7.To further improve the exploration ability of the local search, we incorporate the solution acceptance scheme of simulated annealing into the stochastic VNS and thus obtain a hybrid method of simulated annealing and stochastic VNS (denoted as SA_VNS). To reduce the computation time and make the search process focus on the intensification phase, we use a decreasing acceptance threshold to act as the cooling procedure of simulated annealing. The procedure of the proposed SA_VNS algorithm is illustrated in Algorithm1.Algorithm 1:The main procedure of theSA_VNS.
Begin:
Initialization:
Letπ
0 be the input initial solution. Set π
=
π
0, the acceptance threshold T
=
0.05, and
o
u
t
e
r
l
o
o
p
_
c
o
u
n
t
e
r
=
0.
while (o
u
t
e
r
l
o
o
p
_
c
o
u
n
t
e
r
<
n
/
5) do
(1) Generate a random numberr in [0, 1], and two random integer numbers w and z.
If r > 0.5, generate π
′ = insert (π
,
w
,
z); otherwise generate π
′ = interchange (π
,
w
,
z).
(2) Seti
n
n
e
r
_
c
o
u
n
t
e
r
=
0.
(3)while  (
i
n
n
e
r
_
c
o
u
n
t
e
r
<
n
×
(
n
-
1
)
)  do
Set k
=
1.
while (k
≤
2) do
(1) Generate two random integer numbers w and z.
(2) If k
=
1, generate π
′
′ = i
n
s
e
r
t
(
π
′
,
w
,
z
).
(3) If k
=
2, generate π
′
′ = i
n
t
e
r
c
h
a
n
g
e
(
π
′
,
w
,
z
).
(4) If f
(
π
′
′
)
<
f
(
π
′
)
,
then set
π
′
=
π
′
′ and k
=
1; otherwise set k
=
k
+
1.
end while
Set i
n
n
e
r
_
c
o
u
n
t
e
r
=
i
n
n
e
r
_
c
o
u
n
t
e
r
+
1.
end while
(4) Iff
(
π
′
)
<
f
(
π
0
)
,
set
π
0
=
π
′
and
π
=
π
′, and then go to step 6; otherwise go to step 5.
(5) Iff
(
π
′
)
≥
f
(
π
0
)
but
(
f
(
π
′
)
-
f
(
π
0
)
)
/
f
(
π
0
)
≤
T
,
set
π
=
π
′; otherwise, generate a random
numberr in [0, 1], and then set π
=
π
0 if r
>
0.5 and π
=
π
′ if r
≤
0.5.
(6) SetT
=
T
×
0.95 and o
u
t
e
r
l
o
o
p
_
c
o
u
n
t
e
r
=
o
u
t
e
r
l
o
o
p
_
c
o
u
n
t
e
r
+
1.
end while
Report the improved solution  π
0.
EndIn the PSOVNS proposed by Tasgetiren et al. [31], the stochastic VNS is applied on the global best particle found at each iteration. However, a drawback of such application is that the starting point of the stochastic VNS may be the same solution if the global best particle cannot be improved for a number of consecutive iterations, and consequently the exploration ability of the PSO may be decreased. Therefore, for a given population at iteration t, we propose to use the following strategy.Step 1.
Apply the SA_VNS on the promising particles satisfying(
f
(
π
t
)
-
f
(
g
best
)
)
/
f
(
g
best
)
≤
0.02, in which f
(
π
t
) is the objective value of particle π
t.Step 2.
Update the global best particle. If a new global best particle is found, then further improve it using the SA_VNS.


## 2.5. Population Update Method

It is well known that the advantage of PSO is that it has a high convergence speed. However, this advantage may become the disadvantage for complex scheduling problems because the scheduling problems generally have many local optimal regions in the search space. That is, for the PSO applied to PFSP, some particles may always fly around a local region and thus are trapped in local optimum. Therefore, we propose a solution replacement strategy based on the pathrelinking [37] to remove these solutions with new solutions with good quality.In our algorithm, a particle is viewed as being trapped in local optimum if its personal best solutionp
best has not been improved for a number of consecutive generations (i.e., 20). For these particles, we give them the last chance to stay in the population by applying the path relinking algorithm on it to check if its personal best p
best can be improved. If so, this particle can remain in the population; otherwise, we replace this particle with a new random particle.The path relinking is originally proposed by Glover et al. [37] to generate new solutions by exploring a path that connects an initial solution and a guiding solution. In this path, moves are selected that introduce attributes contained in the guiding solution. To present the path relinking algorithm, we define the distance between two particles X
i and X
j as d
(
X
i
,
X
j
)
=
d
(
π
i
,
π
j
)
=
|
h
(
π
i
)
-
h
(
π
j
)
|, where π is the corresponding job permutation (obtained by the SPV rule) of particle X and h
(
π
) is the Hash function value of π
=
(
π
(
1
)
,
…
,
π
(
n
)
) that is calculated as h
(
π
)
=
∑
k
=
1
n
k
×
π
(
k
)
×
π
(
k
). Then the path relinking algorithm can be described as follows.Step 1.
Setπ
i
=
(
π
i
(
1
)
,
π
i
(
2
)
,
…
,
π
i
(
n
)
) as the initial solution. Find the particle that has the largest distance to particle X
i in the current population (denoted as X
j), and set its corresponding job permutation π
j
=
(
π
j
(
1
)
,
π
j
(
2
)
,
…
,
π
j
(
n
)
) as the guiding solution. Set k
=
1 and the local optimum π
opt
=
π
i.Step 2.
Ifπ
i
(
k
)
≠
π
j
(
k
), find the job with index π
j
(
k
) in π
i and swap it with π
i
(
k
) to generate a new job permutation π
i
′. If π
i
′ is better than π
opt, then set π
opt
=
π
i
′.Step 3.
Setk
=
k
+
1. If k
>
n
-
10, stop; otherwise, go to Step 2.It should be noted that the above path relinking stops whenk is larger than n
-
10, because if particle X
j is better than X
i, then X
i will be replaced by X
j if the path relinking stops when k
=
n, which may result in duplicated particles and thus decrease the search diversification.


## 2.6. Complete Procedure of the Proposed PSO

Step 1.
InitializationStep
1
.1. Set initial values for the population size n
pop, the inertia weight, and the cognitive and social parameters. Set t
=
0 and N
O
i
=
0
(
i
=
1,2
,
…
,
n
pop
) for each particle in the population. Create the initial population and the initial velocities for each particle using the method described in Section 2.3.Step
1.2. Generate the job permutation for each particle in the population using the SPV rule, and calculate the objective value of each particle.Step
1.3. Set the personal best of each particle to be the particle itself and the global best to be the best one among the population.Step 2.
Update Particle PositionsStep
2.1. Update iteration counter t
=
t
+
1.Step
2.2. Update inertia weight by w
=
w
×
β.Step
2.3. For each particle, update the velocity and position values according to (2).Step
2.4. Generate the job permutation for each particle in the current population using the SPV rule, and calculate the objective value of each particle.Step 3.
Local Search PhaseStep
3.1. Use the SA_VNS algorithm to improve the promising particles in the current population and then the global best particle found so far according to the adoption strategy of SA_VNS described in Section 2.4.Step
3.2. For each particle in the current population, update its personal best p
best. If p
best of particle i is improved, then set N
O
i
=
0; otherwise set N
O
i
=
N
O
i
+
1.Step 4.
Population Update. For each particle in the current population, use the population update method described in Section 2.5 to update the current population.Step 5.
Stopping Criterion. If t
>
T
max (the maximum iteration number) or the runtime has reached the limit, then stop; otherwise, go to Step 2.


## 3. Computational Experiments

To test the performance of our PSO algorithm (denoted as PSO*), computational experiments were carried out on the well-known standard benchmark set of Taillard [33] that is composed of 110 instances ranging from 20 jobs and 5 machines to 200 jobs and 20 machines. This benchmark set contains some instances proven to be very difficult to solve. In this benchmark set there are 10 instances for each problem size. Our PSO algorithm was implemented using C++ and tested on a personal PC with Pentium IV 3.0 GHz CPU and 512 MB memory. To make a fair comparison with the PSOVNS, we use the same parameter setting proposed by Tasgetiren et al. [31]. That is, the population size is taken as n
pop
=
2
n; the initial inertia weight is set to w
=
0.9 and never less than 0.4; the decrement factor β for w is taken as 0.975; the acceleration coefficients are set to c
1
=
c
2
=
2; the maximum iteration number T
max is taken as 500.


### 3.1. Results for PFSP with Makespan Minimization

For the makespan minimization objective, our PSO algorithm was compared with other powerful methods, for example, the ant colony algorithm named PACO of Rajendran and Ziegler [29], the genetic algorithm named HGA_RMA of Ruiz et al. [25], the iterated greedy algorithm named IG
_
RS
LS of Ruiz and Stützle [30], and the PSOVNS algorithm of Tasgetiren et al. [31]. The solution quality was measured by the average relative percent deviation (denoted as ARPD) over R replicated runs for each instance in makespan with respect to the best known upper bounds. More specifically, ARPD is calculated as ARPD
=
∑
i
=
1
R
(
(
(
H
i
-
U
i
)
×
100
)
/
U
i
)
/
R, in which H
i is the makespan obtained by a certain algorithm, whereas U
i is the best known upper bound value for Taillard’s instances as of April 2004 for the makespan criterion. As done by many researchers, R is set to R
=
10 in our experiments.The comparison results for these algorithms are given in Table3, in which the values are the average performance of the 10 instances for each problem size. As seen in Table 3, our PSO
* algorithm achieves the best average performance and it obtains the best results for instances of 20
×
5, 20
×
10, 50
×
5, 50
×
10, 100
×
5, 100
×
20, and 200
×
20. The IG
_
RS
LS method also performs well with the HGA_RMA method being close. More specifically, the PACO method cannot obtain the lowest ARPD for any group of problem size compared to other rival methods. The HGA_RMA method has the lowest ARPD for instances of 20
×
10, 50
×
5, and 100
×
10. The IG
_
RS
LS method demonstrates the best results for instances of 20
×
20, 50
×
5, 50
×
10, 50
×
20, and 200
×
10. For instances of 20
×
10, both the HGA_RMA method and the PSOVNS method give the best performance. For instances of 50
×
5, all the four methods except for PACO can obtain the lowest ARPD. For instances of 100
×
5, only the two PSO algorithms give the best performance. Our PSO
* algorithm performs better than its rivals in 100
×
20 and 200
×
20 instances, which have been proven more difficult to solve. Therefore, it can be concluded that our PSO
* algorithm is competitive with other powerful methods in the literature.Table 3
ARPD comparison of different methods for makespan criterion.
Problem
PACO
HGA_RMA
IG_RSLS
PSOVNS
PSO*
2
0
×
5
0.18
0.04
0.04
0.03
0.00
20
×
10
0.24
0.02
0.06
0.02
0.02
20
×
20
0.18
0.05
0.03
0.05
0.05
50
×
5
0.05
0.00
0.00
0.00
0.00
50
×
10
0.81
0.72
0.56
0.57
0.56
50
×
20
1.41
0.99
0.94
1.36
0.99
100
×
5
0.02
0.01
0.01
0.00
0.00
100
×
10
0.29
0.16
0.20
0.18
0.20
100
×
20
1.93
1.30
1.30
1.45
1.18
200
×
10
0.23
0.14
0.12
0.18
0.18
200
×
20
1.82
1.26
1.26
1.35
1.16
Average
0.65
0.43
0.41
0.47
0.39
The bold font is used to highlight the better solutions.
*used to denote our algorithm.


### 3.2. Results for PFSP with Total Flowtime Minimization

For the total flowtime minimization objective, our PSO algorithm was compared with other powerful methods, for example, the constructive heuristics of Framinan and Leisten [20], the ant colony algorithm of Rajendran and Ziegler [29], and the PSOVNS of Tasgetiren et al. [31], using the benchmark problems of Taillard [33]. The solution quality was measured by the relative percent deviation (denoted as RPD) of the best solution found among R
(
R
=
5
) replicated runs for each instance in the total flowtime criterion with respect to the best known results. That is, RPD is calculated as RPD
=
min
{
(
(
H
i
-
U
i
)
×
100
)
/
U
i
,
i
∈
R
}, in which H
i is the total flowtime value obtained by a certain algorithm, whereas U
i is the best result obtained among the algorithms of Framinan and Leisten [20] and Rajendran and Ziegler [29] (this best result is denoted as LR and RZ).For the minimization of the total flowtime criterion, the PSOVNS algorithm [31] is demonstrated to be a very powerful PSO algorithm because it improved 57 out of 90 best known solutions reported in [20, 29]. The comparison results between our PSO* and the PSOVNS are given in Table 4. From this table, we can see that the PSOVNS algorithm can obtain the best results for instances of 20
×
5, 20
×
10, 50
×
5, 100
×
5, and 100
×
10, while our PSO* algorithm obtains the best results for the other large size instances of 20
×
20, 50
×
10, 50
×
20, and 100
×
20. On average, our PSO* shows a much better performance in the solution quality and robustness than the PSOVNS algorithm.Table 4
Performance comparison of PSOVNS and PSO* for total flowtime criterion.
Problem
PSOVNS
PSO*
RPD
CPU (s)
RPD
CPU (s)
2
0
×
5
−0.175
3.18
−0.168
8.11
20
×
10
−0.037
7.21
−0.035
4.86
20
×
20
2.758
11.93
−0.068
24.97
50
×
5
−0.603
41.71
−0.531
40.28
50
×
10
−0.819
74.49
−0.892
44.32
50
×
20
0.857
143.32
−0.543
50.67
100
×
5
−0.570
222.28
−0.546
409.75
100
×
10
−0.692
407.88
−0.636
414.14
100
×
20
−0.104
824.41
−0.801
442.37
Average
0.068
192.93
−0.469
159.941
The bold font is used to highlight the better solutions.
*used to denote our algorithm.


## 3.1. Results for PFSP with Makespan Minimization

For the makespan minimization objective, our PSO algorithm was compared with other powerful methods, for example, the ant colony algorithm named PACO of Rajendran and Ziegler [29], the genetic algorithm named HGA_RMA of Ruiz et al. [25], the iterated greedy algorithm named IG
_
RS
LS of Ruiz and Stützle [30], and the PSOVNS algorithm of Tasgetiren et al. [31]. The solution quality was measured by the average relative percent deviation (denoted as ARPD) over R replicated runs for each instance in makespan with respect to the best known upper bounds. More specifically, ARPD is calculated as ARPD
=
∑
i
=
1
R
(
(
(
H
i
-
U
i
)
×
100
)
/
U
i
)
/
R, in which H
i is the makespan obtained by a certain algorithm, whereas U
i is the best known upper bound value for Taillard’s instances as of April 2004 for the makespan criterion. As done by many researchers, R is set to R
=
10 in our experiments.The comparison results for these algorithms are given in Table3, in which the values are the average performance of the 10 instances for each problem size. As seen in Table 3, our PSO
* algorithm achieves the best average performance and it obtains the best results for instances of 20
×
5, 20
×
10, 50
×
5, 50
×
10, 100
×
5, 100
×
20, and 200
×
20. The IG
_
RS
LS method also performs well with the HGA_RMA method being close. More specifically, the PACO method cannot obtain the lowest ARPD for any group of problem size compared to other rival methods. The HGA_RMA method has the lowest ARPD for instances of 20
×
10, 50
×
5, and 100
×
10. The IG
_
RS
LS method demonstrates the best results for instances of 20
×
20, 50
×
5, 50
×
10, 50
×
20, and 200
×
10. For instances of 20
×
10, both the HGA_RMA method and the PSOVNS method give the best performance. For instances of 50
×
5, all the four methods except for PACO can obtain the lowest ARPD. For instances of 100
×
5, only the two PSO algorithms give the best performance. Our PSO
* algorithm performs better than its rivals in 100
×
20 and 200
×
20 instances, which have been proven more difficult to solve. Therefore, it can be concluded that our PSO
* algorithm is competitive with other powerful methods in the literature.Table 3
ARPD comparison of different methods for makespan criterion.
Problem
PACO
HGA_RMA
IG_RSLS
PSOVNS
PSO*
2
0
×
5
0.18
0.04
0.04
0.03
0.00
20
×
10
0.24
0.02
0.06
0.02
0.02
20
×
20
0.18
0.05
0.03
0.05
0.05
50
×
5
0.05
0.00
0.00
0.00
0.00
50
×
10
0.81
0.72
0.56
0.57
0.56
50
×
20
1.41
0.99
0.94
1.36
0.99
100
×
5
0.02
0.01
0.01
0.00
0.00
100
×
10
0.29
0.16
0.20
0.18
0.20
100
×
20
1.93
1.30
1.30
1.45
1.18
200
×
10
0.23
0.14
0.12
0.18
0.18
200
×
20
1.82
1.26
1.26
1.35
1.16
Average
0.65
0.43
0.41
0.47
0.39
The bold font is used to highlight the better solutions.
*used to denote our algorithm.


## 3.2. Results for PFSP with Total Flowtime Minimization

For the total flowtime minimization objective, our PSO algorithm was compared with other powerful methods, for example, the constructive heuristics of Framinan and Leisten [20], the ant colony algorithm of Rajendran and Ziegler [29], and the PSOVNS of Tasgetiren et al. [31], using the benchmark problems of Taillard [33]. The solution quality was measured by the relative percent deviation (denoted as RPD) of the best solution found among R
(
R
=
5
) replicated runs for each instance in the total flowtime criterion with respect to the best known results. That is, RPD is calculated as RPD
=
min
{
(
(
H
i
-
U
i
)
×
100
)
/
U
i
,
i
∈
R
}, in which H
i is the total flowtime value obtained by a certain algorithm, whereas U
i is the best result obtained among the algorithms of Framinan and Leisten [20] and Rajendran and Ziegler [29] (this best result is denoted as LR and RZ).For the minimization of the total flowtime criterion, the PSOVNS algorithm [31] is demonstrated to be a very powerful PSO algorithm because it improved 57 out of 90 best known solutions reported in [20, 29]. The comparison results between our PSO* and the PSOVNS are given in Table 4. From this table, we can see that the PSOVNS algorithm can obtain the best results for instances of 20
×
5, 20
×
10, 50
×
5, 100
×
5, and 100
×
10, while our PSO* algorithm obtains the best results for the other large size instances of 20
×
20, 50
×
10, 50
×
20, and 100
×
20. On average, our PSO* shows a much better performance in the solution quality and robustness than the PSOVNS algorithm.Table 4
Performance comparison of PSOVNS and PSO* for total flowtime criterion.
Problem
PSOVNS
PSO*
RPD
CPU (s)
RPD
CPU (s)
2
0
×
5
−0.175
3.18
−0.168
8.11
20
×
10
−0.037
7.21
−0.035
4.86
20
×
20
2.758
11.93
−0.068
24.97
50
×
5
−0.603
41.71
−0.531
40.28
50
×
10
−0.819
74.49
−0.892
44.32
50
×
20
0.857
143.32
−0.543
50.67
100
×
5
−0.570
222.28
−0.546
409.75
100
×
10
−0.692
407.88
−0.636
414.14
100
×
20
−0.104
824.41
−0.801
442.37
Average
0.068
192.93
−0.469
159.941
The bold font is used to highlight the better solutions.
*used to denote our algorithm.


## 4. Conclusions

This paper presents a PSO-based hybrid metaheuristic for the permutation flowshop problems to minimize the makespan and the total flowtime. In this algorithm, a hybrid method of simulated annealing and stochastic variable neighborhood search is incorporated to improve the exploitation ability, and a solution replacement strategy based on the path relinking method is developed to improve the exploration ability. Computational experiments are carried out to test the performance of the proposed PSO-based hybrid metaheuristic, and the results show that the proposed algorithm is competitive or superior to some other powerful algorithms in the literature for this problem. Future research may lie in the application of this algorithm in practical production scheduling problems.



---

*Source: 902950-2014-01-29.xml*