A hybrid particle swarm optimization for job shop scheduling problem

(1)

A hybrid particle swarm optimization for job shop

scheduling problem

D.Y. Sha

a,b,*

, Cheng-Yu Hsu

b

a

Department of Business Administration, Asia University, 500 Liufeng Road, Wufong, Taichung 413, Taiwan, ROC

b

Department of Industrial Engineering and Management, National Chiao Tung University, 1001 Ta Hsueh Road, Hsinchu 300, Taiwan, ROC

Received 14 January 2006; received in revised form 14 September 2006; accepted 14 September 2006 Available online 27 October 2006

Abstract

A hybrid particle swarm optimization (PSO) for the job shop problem (JSP) is proposed in this paper. In previous research, PSO particles search solutions in a continuous solution space. Since the solution space of the JSP is discrete, we modified the particle position representation, particle movement, and particle velocity to better suit PSO for the JSP. We modified the particle position based on preference list-based representation, particle movement based on swap operator, and particle velocity based on the tabu list concept in our algorithm. Giffler and Thompson’s heuristic is used to decode a particle position into a schedule. Furthermore, we applied tabu search to improve the solution quality. The computational results show that the modified PSO performs better than the original design, and that the hybrid PSO is better than other traditional metaheuristics.

Keywords: Job shop problem; Scheduling; Particle swarm optimization

1. Introduction

The job shop scheduling problem (JSP) is one of the most difficult combinatorial optimization problems. The JSP can be briefly stated as follows (French, 1982; Gen & Cheng, 1997). There are n jobs to be processed through m machines. We shall suppose that each job must pass through each machine once and once only. Each job should be processed through the machines in a particular order, and there are no precedence con-straints among different job operations. Each machine can process only one job at a time, and it cannot be interrupted. Furthermore, the processing time is fixed and known. The problem is to find a schedule to min-imize the makespan (Cmax), that is, the time required to complete all jobs.

Garey, Johnson, and Sethi (1976)demonstrated that JSP is NP-hard, so it cannot be exactly solved in a

reasonable computation time. Many approximate methods have been developed in the last decade to solve

* _{Corresponding author. Tel.: +886 4 23323456x1936; fax: +886 4 2331 6699.}

E-mail addresses:yjsha@mail.nctu.edu.tw(D.Y. Sha),cyhsu.iem92g@nctu.edu.tw(C.-Y. Hsu).

(2)

JSP, such as simulated annealing (SA) (Lourenc¸o, 1995), tabu search (TS) (Nowicki & Smutnicki, 1996; Pezz-ella & Merelli, 2000; Sun, Batta, & Lin, 1995), and genetic algorithm (GA) (Bean, 1994; Gonc¸alves, Mendes, &

Resende, 2005; Kobayashi, Ono, & Yamamura, 1995; Wang & Zheng, 2001). We applied a new evolutionary

search technique – particle swarm optimization (PSO) – to solve the JSP in this paper.

The optimal JSP solution should be an active schedule. In an active schedule the processing sequence is such that no operation can be started any earlier without delaying some other operation (French, 1982). To reduce the search solution space, the tabu search proposed bySun et al. (1995)searches solutions within the set of active schedules. In our algorithm, we appliedGiﬄer and Thompson’s (1960) heuristic to decode a particle position into a schedule. Furthermore, we applied a tabu search to improve the solution quality.

The background of particle swarm optimization (PSO) is introduced in the next section. In Section3, we propose a hybrid PSO for the JSP. In Section4, we test the hybrid PSO onFisher and Thompson (1963) and

Lawrence (1984) and Taillard (1993) test problems. Finally, conclusions and remarks for further works are

given in Section5.

2. The background of particle swarm optimization

Particle swarm optimization (PSO) was developed byKennedy and Eberhart (1995). PSO is a population-based optimization algorithm. Each particle is an individual and the swarm is composed of particles. The problem solution space is formulated as a search space. Each position in the search space is a correlated solu-tion of the problem. Particles cooperate to ﬁnd out the best posisolu-tion (best solusolu-tion) in the search space (solu-tion space).

Particles move toward the pbest position and gbest position with each iteration. The pbest position is the best position found by each particle so far. Each particle has its own pbest position. The gbest position is the best position found by the swarm so far. The particle moves itself according to its velocity. The velocities are randomly generated toward pbest and gbest positions. For each particle k and dimension j, the velocity and position of particles can be updated by the following equations:

vkj w vkjþ c1 rand1 ðpbestkj xkjÞ þ c2 rand2 ðgbestj xkjÞ ð1Þ

xkj xkjþ vkj ð2Þ

In Eqs.(1)and(2), vkjis the velocity of particle k on dimension j, and xkjis the position of particle k on

dimen-sion j. The pbestkj is the pbest position of particle k on dimension j, and gbestj is the gbest position of the

swarm on dimension j. The inertia weight w was ﬁrst proposed byShi and Eberhart (1998a, 1998b), and is used to control exploration and exploitation. The particles maintain high velocities with a larger w, and low velocities with a smaller w. A larger w can prevent particles from becoming trapped in local optima, and a smaller w encourages particles exploiting the same search space area. The constants c1and c2are used

to decide whether particles prefer moving toward a pbest position or gbest position. The rand1and rand2are

random variables between 0 and 1. The process for PSO is as follows:

Step 1: Initialize a population of particles with random positions and velocities on d dimensions in the search space.

Step 2: Update the velocity of each particle, according to Eq.(1). Step 3: Update the position of each particle, according to Eq.(2).

Step 4: Map the position of each particle into solution space and evaluate its fitness value according to the desired optimization fitness function. At the same time, update pbest and gbest position if necessary. Step 5: Loop to step 2 until a criterion is met, usually a sufficiently good fitness or a maximum number of

iterations.

The original PSO design is suited to a continuous solution space. For better suiting to combinatorial opti-mization problems, we have to modify PSO position representation, particle velocity, and particle movement.

Zhang, Li, Li, and Huang (2005)proposed a PSO for resource-constrained project scheduling, and compared

(3)

prior-ity values), and (2) permutation-based representation (particle position represented by a sequential order of activities).Zhang’s results (2005)showed that permutation-based representation is better than priority-based representation.

We modiﬁed the particle position based on preference list-based representation (Davis, 1985) and the par-ticle movement based on a swap operator in this paper. These will be discussed in Section3.

3. A hybrid particle swarm optimization

In this section, we will ﬁrst describe how to associate a particle position into a schedule with two diﬀerent position representations, respectively, the priority-based representation and preference list-based representa-tion. If we implement the priority-based representation, the particle position consists of continuous variables, and it is suited to the original PSO design, as described in Section2. When we implement the preference list-based representation, we have to modify the particle velocity and particle movement, as described in Sections

3.2and3.3. Besides, we propose a diversiﬁcation strategy and a local search procedure for better performance. 3.1. Position representation

3.1.1. Priority-based representation

When we implement the original PSO design, as described in Section2(i.e., the particles search solutions in a continuous solution space), each value of a particle position represents the associated operation priority. For an n-job m-machine problem, we can represent the particle k position by an m· n matrix, i.e.

Xk_¼ xk 11 xk12 xk1n xk 21 x k 22 x k 2n .. . xk m1 xkm2 xkmn 2 6 6 6 6 6 6 4 3 7 7 7 7 7 7 5 where xk

ij denotes the priority of operation oijand oijis the operation of job j that needs to be processed on

machine i. We can map (or decode) a particle position into an active schedule usingGiﬄer and Thompson’s

(1960)heuristic. We brieﬂy describe the G&T algorithm as follows:

Notation:

(i, j): the operation of job j that needs to be processed on machine i. S: the partial schedule that contains scheduled operations.

X: the set of schedulable operations.

s(i,j): the earliest time at which operation (i, j)2 X could be started.

p(i,j): the processing time of operation (i, j).

f(i,j): the earliest time at which operation (i, j)2 X could be ﬁnished, f(i,j)= s(i,j)+ p(i,j).

G&T algorithm:

Step 1: Initialize S = /; X is initialized to contain all operations without predecessors. Step 2: Determine f*_{= min}

(i,j)2X{f(i,j)} and the machine m*on which f*could be realized.

Step 3: (1) Identify the operation set (i0_{, j}0₎_{2 X such that (i}0_{, j}0_{) requires machine m}*_{, and s}

(i0_,j0₎< f*. (2) Choose (i, j) from the operation set identiﬁed in (1) with the largest priority.

(3) Add (i, j) to S.

(4)

Step 4: If a complete schedule has been generated, stop. Else, delete (i, j) from X and include its immediate successor in X, then go to Step 2.

For example, there are two jobs and two machines, as shown onTable 1, and the position of particle k is Xk _¼ 0:6 1:3

0:8 0:5

:

We can use the G&T algorithm to decode Xkinto a schedule in the following steps: Initialization

Step 1: S = /; X = {(1, 1), (2, 2)}. Iteration 1

Step 2: s(1,1)= 0, s(2,2)= 0, f(1,1)= 5, f(2,2)= 4; f*= min{f(1,1), f(2,2)} = 4, m*= 2.

Step 3: Identify the operation set {(2, 2)}; choose operation (2, 2), which has the largest priority, and add it into schedule S, as illustrated inFig. 1(a).

Step 4: Update X = {(1, 1), (1, 2)}. Iteration 2

Step 3: Identify the operation set {(1, 1), (1, 2)}; choose operation (1, 2), which has the largest priority, and add it into schedule S, as illustrated inFig. 1(b).

Step 4: Update X = {(1, 1)}. Iteration 3

Step 2: s(1,1)= 7, f(1,1)= 12; f*= min{f(1,1)} = 12, m*= 1.

Step 3: Identify the operation set {(1, 1)}; choose operation (1, 1), which has the largest priority, and add it into schedule S, as illustrated inFig. 1(c).

Step 4: Update X = {(2, 1)} Iteration 4

Step 2: s(2,1)= 12, f(2,1)= 16; f*= min{f(2,1)} = 16, m*= 2.

Step 3: Identify the operation set {(2, 1)}; choose operation (2, 1), which has the largest priority, and add it into schedule S, as illustrated inFig. 1(d).

Step 4: A complete schedule has been generated, and then stops.

However, there is a shortcoming of priority-based representation. The schedules of two particles may be quite diﬀerent even though their positions are very close to each other. For example, if there are six operations to be sorted on a machine, and there are two positions of two particles as follows:

position 1 : ½0:25; 0:27; 0:21; 0:24; 0:26; 0:23

Table 1 A 2· 2 example

Jobs Machine sequence Processing times

1 1, 2 p(1,1)= 5, p(2,1)= 4

(5)

position 2 : ½0:22; 0:25; 0:23; 0:26; 0:24; 0:21

then we sort the operations according to the decreasing order of their position values as follows: permutation 1 : ½ 2 5 1 4 6 3

permutation 2 : ½ 4 2 5 3 1 6

We can ﬁnd that these two permutations are quite diﬀerent even though the particle positions are very close to each other. This is because the location in the permutation of one operation depends on the position values of other operations.

3.1.2. Preference list-based representation

In the preference list-based representation, there is a preference list for each machine. For an n-job m-machine problem, we can also represent the particle k position by an m· n matrix, and the ith row is the pref-erence list of machine i, i.e.

Xk¼ xk 11 xk12 xk1n xk 21 xk22 xk2n .. . xk m1 xkm2 xkmn 2 6 6 6 6 4 3 7 7 7 7 5 where xk

ij 2 f1; 2; . . . ; ng denotes the job on location j in the preference list of machine i. We can also useGiﬄer

and Thompson’s (1960) heuristic to map a particle position into an active schedule. The same example, as

shown inTable 1, and the position of particle k is

Xk_¼ 2 1 1 2 : Initialization M1 M2 (2, 2) time 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 Partial schedule after the operation (2, 2) scheduled.

M1 (1, 2) M2 (2, 2)

time 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 Partial schedule after the operation (1, 2) scheduled.

M1 (1, 2) (1, 1) M2 (2, 2)

time 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 Partial schedule after the operation (1, 1) scheduled.

M1 (1, 2) (1, 1)

M2 (2, 2) (2, 1)

time 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 A complete schedule after the operation (2, 1) scheduled.

a

b

c

d

(6)

Step 1: S = /; X = {(1, 1), (2, 2)}. Iteration 1

Step 2: s(1,1)= 0, s(2,2)= 0, f(1,1)= 5, f(2,2)= 4; f*= min{f(1,1),f(2,2)} = 4, m*= 2.

Step 3: Identify the operation set {(2, 2)}; choose operation (2, 2), which is ahead of others in the preference list of machine 2, and add it into schedule S, as illustrated inFig. 1(a).

Step 4: Update X = {(1, 1), (1, 2)}. Iteration 2

Step 3: Identify the operation set {(1, 1), (1, 2)}; choose operation (1, 2), which is ahead of others in the pref-erence list of machine 1, and add it into schedule S, as illustrated inFig. 1(b).

Step 4: Update X = {(1, 1)}. Iteration 3

Step 2: s(1,1)= 7, f(1,1)= 12; f*= min{f(1,1)} = 12, m*= 1.

Step 3: Identify the operation set {(1, 1)}; choose operation (1, 1), which is ahead of others in the preference list of machine 1, and add it into schedule S, as illustrated inFig. 1(c).

Step 4: Update X = {(2, 1)} Iteration 4

Step 2: s(2,1)= 12, f(2,1)= 16; f*= min{f(2,1)} = 16, m*= 2.

Step 3: Identify the operation set {(2, 1)}; choose operation (2, 1), which is ahead of others in the preference list of machine 2, and add it into schedule S, as illustrated inFig. 1(d).

Step 4: A complete schedule has been generated, and then stops.

The preference list-based PSO we proposed diﬀers from the original PSO design in that the pbest solutions and gbest solution do not record the best positions found so far, but rather the best schedules generated by the G&T algorithm. For the above example, we do not record the particle position Xk, but record the schedule

Sk ¼ 2 1

2 1

" #

into pbest and gbest solutions if necessary. Because the particle position representation differs from the origi-nal design, we also modified the movement based on the swap operator. This will be discussed in Section3.3. 3.2. Modified particle velocity

When a particle moves in a continuous solution space, due to inertia, the particle velocity not only moves the particle to a better position, but also prevents the particle from moving back to the current position. The velocity can be controlled by inertia weight w in Eq.(1). The larger the inertia weight, the harder the particle backs to the current position.

If we implement preference list-based representation, the velocity of operation oijof particle k is denoted by

vk

ij, vkij2 f0; 1g, where oijis the operation of job j that needs to be processed on machine i. When vkijequals 1, it

means that operation oijin the preference list of particle k (the position matrix, Xk) has just been moved to the

current location, and we should not move it in this iteration. On the contrary, if operation oijis moved to a

new location in this iteration, we set vk

(7)

been moved in the next few iterations. The particle velocity can prevent recently moved operations from mov-ing back to the original location in the next iterations.

Just as the original PSO is applied to a continuous solution space, inertia weight w is used to control particle velocities. We randomly update velocities at the beginning of the iteration. For each particle k and operation oij, if vkij equals 1, vkij will be set to 0 with probability (1 w). This means that if operation oijis ﬁxed on the

current location in the preference list of particle k, oijis allowed to move in this iteration with probability

(1 w). The newly moved operations will then be fixed for more iterations with larger inertia weight, and fixed for less iterations with smaller inertia weight. The pseudo code for updating velocities is given inFig. 2. 3.3. Modified particle movement

The modiﬁed particle movement is based on the swap operator. If vk

ij¼ 0, the job j on x k

i will be moved to

the corresponding location of pbestk_i with probability c1, and will be moved to the corresponding location of

gbestiwith probability c2. Where xki is the preference list of machine i of particle k, pbest k

i is the preference list

of machine i of kth pbest solution, gbestiis the preference list of machine i of gbest solution, c1and c2are

con-stant between 0 and 1, and c1+ c261. The process is described as follows:

Step 1: Randomly choose a location l in xk i.

Step 2: Denote the job on location l in xk i by J1.

Step 3: Find out the location of J1in pbestki with probability c1, or ﬁnd out the location of J1in gbestiwith

probability c2. Denote the location that has been found in pbestki or gbestiby l0, and denote the job in

location l0 _{in x}k i by J2.

Step 4: If J2has been denoted, vkiJ1 ¼ 0, and v

k

iJ2¼ 0, then swap J1and J2in x

k

i, and set v k iJ1 1. Step 5: If all the locations in xk

i have been considered, then stop. Otherwise, if l < n, then set l l + 1, else

l 1, and go to Step 2, where n is the number of jobs. For example, there is a ﬁve-job problem, and xk

i, pbest k

i, gbesti, and vki are shown inFig. 3(a). We set c1= 0.5

and c2= 0.3 in this instance.

In Step 1, we randomly choose a location l = 3. In Step 2, the job in the 3rd location in xk

i is job 4, i.e.

J1= 4. In Step 3, we generate a random variable rand between 0 and 1, and the generated random variable

rand is 0.6. Since c1< rand 6 c1+ c2, we ﬁnd out the location of J1 in gbesti. The location l0= 5, and the

job in the 5th location in xk

i is job 5, i.e., J2= 5. Steps 1–3 are shown in Fig. 3(b). In Step 4, since vki4¼ 0

and vk

i5¼ 0, swap jobs 4 and 5 in xki and set vki4 1 are shown inFig. 3(c). In Step 5, set l 4, and go to Step

2. Repeat the procedure until all the locations in xk

i have been considered.

We also adopt a mutation operator in our algorithm. After a particle moves to a new position, we randomly choose a machine and two jobs on the machine, and then swap these two jobs, disregarding vk

ij. The particle

movement pseudo code is given inFig. 4. 3.4. The diversiﬁcation strategy

If all the particles have the same pbest solutions, they will be trapped into local optima. To prevent such a situation, we proposed a diversiﬁcation strategy to keep the pbest solutions diﬀerent (i.e., keeps the makespans

do foreachparticlekandoperationoij

rand ~ U(0,1) if =1) ( ) ← ≥ (vk and rand w ij then 0 k ij v end if end for

(8)

of pbest solutions diﬀerent). In the diversiﬁcation strategy, the pbest solution of each particle is not the best solution found by the particle itself, but one of the best N solutions found by the swarm so far where N is the size of the swarm. Once any particle generates a new solution, the pbest and gbest solutions will be updated in these three situations:

1. If the particle’s ﬁtness value is better than the ﬁtness value of the gbest solution, set the worst pbest solution equal to the current gbest solution, and set the gbest solution equal to the particle solution.

2. If the particle’s ﬁtness value is worse than the gbest solution, but better then the worst pbest solution and not equal to any gbest or pbest solution, set the worst pbest solution equal to the particle solution. 3. If the particle’s ﬁtness value is equal to any pbest or gbest solution, replace the pbest or gbest solution

(whose ﬁtness value is equal to the particle ﬁtness value) with the particle solution.

The pseudo code for updating the pbest solution and gbest solution with diversiﬁcation strategy is given in

Fig. 5.

3.5. Local search

The tabu search is a metaheuristic approach and a strong local search mechanism. In the tabu search, the algorithm starts from an initial solution and improves it iteratively to ﬁnd a near-optimal solution. This method was proposed and formalized primarily byGlover (1986, 1989, 1990). We applied the tabu search pro-posed byNowicki and Smutnicki (1996) but without back jump tracking. We brieﬂy describe Nowicki and Smutnicki’s method as follows:

a

b

c

(9)

3.5.1. The neighborhood structure

Nowicki and Smutnicki’s method randomly chooses a critical path in the current schedule, and then rep-resents the critical path in terms of blocks. The neighborhood exchanges the ﬁrst two and the last two oper-ations in every block, but excludes the ﬁrst and last operoper-ations in the critical path. The research of Jain,

Rangaswamy, and Meeran (2000)shows that the strategy used to generate the critical path does not materially

aﬀect the ﬁnal solution. Therefore, in this paper, we randomly choose one critical path if there is more than one critical path. For example, there is a schedule for a four-job, three-machine problem, as shown in

Fig. 6(a). We can ﬁnd that there are two critical paths: CP1= {o31, o11, o13, o33} and CP2= {o31, o32, o22,

o21, o24, o14}, where oijis the operation of job j that needs to be processed on machine i. If we randomly choose

CP2, we can represent CP2in terms of blocks: {o31, o32}, {o22, o21, o24}, and {o14}. The possible moves in this

schedule are exchanging {o22, o21} or {o21, o24} (seeFig. 6(b)).

(10)

3.5.2. Tabu list

The tabu list consists of maxt operation pairs that have been moved in the last maxt moves in the tabu search. If a move {oiJ1, oiJ2} has been performed, this move replaces the oldest move in the tabu list, and mov-ing these same two operations is not permitted while the move is recorded in the tabu list.

Fig. 5. Pseudo code of updating pbest solution and gbest solution with diversiﬁcation strategy.

o11 o13 o12 o14 o23 o22 o21 o24 o31 o32 o34 o33 M1 M2 M3

An instance of job shop schedule.

o11 o13 o12 o14 o23 o22 o21 o24 o31 o32 o34 o33 M1 M2 M3

Neighborhood defined by Nowicki & Smutnicki (1996).

a

b

(11)

Table 2

Computational result of FT and LA test problems Problem Size

(n· m)

Best Known Solution (BKS)

Shifting bottleneck Tabu Search

SBI SBII SB-RGLS1 SB-RGLS2 ACM TSAB TSSB

Adams et al. (1988)

Balas and Vazacopoulos (1998) Sun et al. (1995) Nowicki and Smutnicki (1996) Pezzella and Merelli (2000) FT06 6· 6 55 55 55 – – – 55 55 FT10 10· 10 930 1015 930 930 930 930 930 930 FT20 20· 5 1165 1290 1178 – – – 1165 1165 LA01 10· 5 666 666 666 – – – 666 666 LA02 10· 5 655 720 669 655 655 – 655 655 LA03 10· 5 597 623 605 – – – 597 597 LA04 10· 5 590 597 593 – – – 590 590 LA05 10· 5 593 593 593 – – – 593 593 LA06 15· 5 926 926 926 – – – 926 926 LA07 15· 5 890 890 890 – – – 890 890 LA08 15· 5 863 868 863 – – – 863 863 LA09 15· 5 951 951 951 – – – 951 951 LA10 15· 5 958 959 959 – – – 958 958 LA11 20· 5 1222 1222 1222 – – – 1222 1222 LA12 20· 5 1039 1039 1039 – – – 1039 1039 LA13 20· 5 1150 1150 1150 – – – 1150 1150 LA14 20· 5 1292 1292 1292 – – – 1292 1292 LA15 20· 5 1207 1207 1207 – – – 1207 1207 LA16 10· 10 945 1021 978 – – 975 945 945 LA17 10· 10 784 796 787 – – 784 784 784 LA18 10· 10 848 891 859 – – 848 848 848 LA19 10· 10 842 875 860 842 842 842 842 842 LA20 10· 10 902 924 914 – – 902 902 902 LA21 15· 10 1046 1172 1084 1048 1046 1074 1047 1046 LA22 15· 10 927 1040 944 – – 941 927 927 LA23 15· 10 1032 1061 1032 – – 1032 1032 1032 LA24 15· 10 935 1000 976 937 935 954 939 938 LA25 15· 10 977 1048 1017 977 977 1010 977 979 LA26 20· 10 1218 1304 1224 – – 1218 1218 1218 LA27 20· 10 1235 1325 1291 1235 1235 1277 1236 1235 LA28 20· 10 1216 1256 1250 – – 1245 1216 1216 LA29 20· 10 1157 1294 1239 1164 1164 1234 1160 1168 LA30 20· 10 1355 1403 1355 – – 1355 1355 1355 LA31 30· 10 1784 1784 1784 – – 1784 1784 1784 LA32 30· 10 1850 1850 1850 – – 1850 1850 1850 LA33 30· 10 1719 1719 1719 – – 1719 1719 1719 LA34 30· 10 1721 1721 1721 – – 1721 1721 1721 LA35 30· 10 1888 1888 1888 – – 1888 1888 1888 LA36 15· 15 1268 1351 1305 1268 1268 1303 1268 1268 LA37 15· 15 1397 1485 1423 1397 1397 1422 1407 1411 LA38 15· 15 1196 1280 1255 1198 1196 1245 1196 1201 LA39 15· 15 1233 1321 1273 1233 1233 1269 1233 1240 LA40 15· 15 1222 1326 1269 1226 1224 1255 1229 1233 Average gap 3.8796% 1.3838% 0.1157% 0.0591% 1.5184% 0.0501% 0.1015% No. of instance 43 43 13 13 26 43 43 No. of BKS obtained 16 20 8 11 13 37 36

(12)

3.5.3. Back jump tracking

When ﬁnding a new best solution, store the current state (the new best solution, set of moves, and tabu list) in a list L. After the tabu search algorithm performs maxiter_tabu iterations, restart the tabu search algorithm from the latest recorded state, and repeat it until the list L is empty. We did not implement the back jump tracking in our algorithm to reduce computation time.

Table 2 (continued)

Genetic algorithm Particle Swarm Optimization

GASA HGA-Param PSO-priority based PSO-permutation based HPSO

Wang and Zheng (2001) Gonc¸alves et al. (2005) Best solution Average Best solution Average Best solution Average 55 55 55 58.9 55 55.0 55 55.0 930 930 1007 1086.0 937 965.2 930 932.0 1165 1165 1242 1296.7 1165 1178.8 1165 1165.0 666 666 681 705.0 666 666.0 666 666.0 – 655 694 729.7 655 662.1 655 655.0 – 597 633 657.5 597 602.3 597 597.0 – 590 611 648.1 590 592.9 590 590.0 – 593 593 601.1 593 593.0 593 593.0 926 926 926 940.2 926 926.0 926 926.0 – 890 890 941.0 890 890.0 890 890.0 – 863 863 896.6 863 863.0 863 863.0 – 951 953 991.8 951 951.0 951 951.0 – 958 958 976.1 958 958.0 958 958.0 1222 1222 1222 1235.3 1222 1222.0 1222 1222.0 – 1039 1039 1058.4 1039 1039.0 1039 1039.0 – 1150 1150 1179.0 1150 1150.0 1150 1150.0 – 1292 1292 1292.2 1292 1292.0 1292 1292.0 – 1207 1232 1271.7 1207 1207.0 1207 1207.0 945 945 1006 1033.5 945 969.8 945 945.2 – 784 833 883.5 784 787.1 784 784.0 – 848 901 959.9 848 856.8 848 848.0 – 842 895 945.8 842 851.5 842 842.0 – 907 963 1014.0 907 913.3 902 902.3 1058 1046 1201 1247.5 1055 1085.5 1046 1049.8 – 935 1046 1142.5 935 950.5 927 927.0 – 1032 1146 1205.1 1032 1032.0 1032 1032.0 – 953 1082 1140.9 937 967.8 935 937.9 – 986 1107 1176.6 983 1005.9 977 978.2 1218 1218 1409 1468.0 1218 1219.7 1218 1218.0 – 1256 1437 1495.4 1252 1269.1 1235 1251.4 – 1232 1434 1487.4 1216 1241.7 1216 1216.0 – 1196 1359 1429.8 1179 1215.8 1163 1168.8 – 1355 1517 1557.0 1355 1355.0 1355 1355.0 1784 1784 1886 1942.5 1784 1784.0 1784 1784.0 – 1850 2000 2065.6 1850 1850.0 1850 1850.0 – 1719 1832 1896.8 1719 1719.0 1719 1719.0 – 1721 1876 1953.5 1721 1721.0 1721 1721.0 – 1888 2027 2074.5 1888 1888.0 1888 1888.0 1292 1279 1437 1541.0 1291 1317.5 1268 1271.3 – 1408 1539 1628.0 1442 1475.1 1397 1401.6 – 1219 1370 1445.1 1228 1251.1 1196 1200.5 – 1246 1436 1499.4 1233 1285.6 1233 1233.0 – 1241 1380 1457.4 1236 1258.0 1224 1226.2 0.2764% 0.3916% 7.4021% 12.0940% 0.3719% 1.3491% 0.0159% 0.1091% 11 43 43 43 43 9 31 10 31 41

(13)

We implement a tabu search procedure after a particle generates a new solution for further improved solu-tion quality. The tabu search will be stopped after 100 moves that do not improve the solusolu-tion. The research of

Jain et al. (2000)shows that the solution quality of tabu search (Nowicki & Smutnicki, 1996) is mainly aﬀected

by its initial solution. Therefore, in the hybrid PSO, the purpose of the PSO process is to provide good and diverse initial solutions to the tabu search.

4. Computational results

There are three PSOs we tested: (1) priority-based PSO, of which the particle position is represented by the priorities of operations, and implements the original PSO design as described in Section2; (2) preference list-based PSO, of which the particle position is represented by a preference list of machines; (3) hybrid PSO (HPSO), which is the preference list-based PSO with a local search mechanism. The PSOs were tested on

Fish-er and Thompson (1963) (FT06, FT10, and FT20), Lawrence (1984)(LA01 to LA40) and Taillard (1993)

(TA01 to TA80) test problems. These problems are available on the OR-Library web site (Beasley, 1990) (URL:http://people.brunel.ac.uk/~mastjjb/jeb/info.html) and Taillard’s web site (URL:http://ina2.eivd.ch/

Collaborateurs/etd/problemes.dir/ordonnancement.dir/ordonnancement.html).

In the preliminary experiment, four swarm sizes N (10, 20, 30, 50) were tested, where N = 30 was superior and used for all further studies. The other parameters of the priority-based PSO were set to the same common settings as most of the previous research: c1= 2.0, c2= 2.0, the inertia weight w is decreased linearly from 0.9

to 0.4 during a run, and the maximum value ofjxijj and jvijj, Xmaxand Vmaxare equal to the number of jobs n

and n/5, respectively.

The parameters of the preference list-based PSO are determined experimentally. The parameters c1and c2

were tested between 0.1 and 0.5 in increments of 0.1, and the parameter w was tested between 0 and 0.9 in increments of 0.1. The settings c1= 0.5, c2= 0.3 and w = 0.5 were superior. The length of the tabu list maxt

was set to 8 where the value is derived fromNowicki and Smutnicki (1996). The tabu search will be stopped after 100 moves that do not improve the solution. The priority-based PSO and the preference list-based PSO will be terminated after 105iterations, and HPSO will be terminated after 103iterations. The number of iter-ations is determined by the computation time compared withPezzella and Merelli (2000) and Gonc¸alves et al. (2005).

The program was coded in Visual C++, optimized by speed, and run on an AMD Athlon 1700+ PC 20 times for each of the 123 problems. The proposed algorithm is compared with Shifting Bottleneck (Adams,

Table 3

Computation time of FT and LA test problems (in CPU seconds)

Problem Size

(n· m)

HGA-Param

Gonc¸alves et al. (2005)a

Particle swarm optimizationb

PSO-priority based PSO-permutation

based HPSO Best solution time Total time Best solution time Total time Best solution time Total time FT06 6· 6 13 0.0 34 0.0 32 0.0 28 FT10 10· 10 292 1.0 112 21.7 91 4.1 157 FT20 20· 5 204 3.0 180 19.2 138 19.8 219 LA01-05 10· 5 40 0.4 60 5.3 50 0.5 38 LA06-10 15· 5 94 1.0 114 0.1 92 0.1 61 LA11-15 20· 5 192 3.5 177 0.5 143 0.1 100 LA16-20 10· 10 227 0.6 109 15.5 90 19.9 139 LA21-25 15· 10 602 4.8 208 37.2 164 59.6 295 LA26-30 20· 10 1303 12.6 325 103.1 259 90.5 579 LA31-35 30· 10 3691 46.9 652 31.4 520 3.0 1462 LA36-40 15· 15 1920 7.4 331 68.4 254 105.2 471

a _{Run on an AMD Thunderbird 1.333 GHz PC.} b_{Run on an AMD Athlon 1700+ PC.}

(14)

Table 4

Computational result of TA test problems

Problem Size

(n· m)

Optimal solution (or upper bound)

TSAB TSSB HPSO Nowicki and Smutnicki (1996) Pezzella and Merelli (2000)

Best solution Average

TA01 15· 15 1231 1241 1231 1236 TA02 15· 15 1244 1244 1244 1244 1245 TA03 15· 15 1218 1222 1222 1218 1224 TA04 15· 15 1175 1175 1175 1180 TA05 15· 15 1224 1233 1229 1224 1233 TA06 15· 15 1238 1245 1238 1248 TA07 15· 15 1227 1228 1228 1229 TA08 15· 15 1217 1220 1220 1217 1220 TA09 15· 15 1274 1282 1291 1274 1283 TA10 15· 15 1241 1259 1250 1249 1264 TA11 20· 15 (1359) 1371 1366 1386 TA12 20· 15 (1367) 1377 1379 1370 1380 TA13 20· 15 (1342) 1362 1350 1364 TA14 20· 15 1345 1345 1345 1345 1350 TA15 20· 15 (1339) 1360 1350 1364 TA16 20· 15 (1360) 1370 1368 1377 TA17 20· 15 1462 1481 1473 1480 TA18 20· 15 (1396) 1413 1426 1407 1425 TA19 20· 15 (1335) 1352 1351 1335 1353 TA20 20· 15 (1348) 1362 1366 1358 1373 TA21 20· 20 (1644) 1659 1658 1679 TA22 20· 20 (1600) 1623 1614 1625 TA23 20· 20 (1557) 1573 1559 1578 TA24 20· 20 (1646) 1659 1654 1664 TA25 20· 20 (1595) 1606 1616 1632 TA26 20· 20 (1645) 1657 1666 1662 1679 TA27 20· 20 (1680) 1697 1690 1712 TA28 20· 20 (1603) 1622 1617 1627 TA29 20· 20 (1625) 1629 1635 1634 1645 TA30 20· 20 (1584) 1614 1589 1613 TA31 30· 15 1764 1766 1771 1766 1772 TA32 30· 15 (1795) 1841 1840 1823 1848 TA33 30· 15 (1791) 1832 1833 1818 1834 TA34 30· 15 (1829) 1846 1844 1879 TA35 30· 15 2007 2007 2007 2010 TA36 30· 15 1819 1825 1825 1843 TA37 30· 15 1771 1815 1813 1795 1808 TA38 30· 15 1673 1700 1697 1681 1701 TA39 30· 15 1795 1811 1815 1796 1810 TA40 30· 15 (1674) 1720 1725 1698 1714 TA41 30· 20 (2018) 2045 2047 2071 TA42 30· 20 (1949) 1979 1970 1984 TA43 30· 20 (1858) 1898 1899 1928 TA44 30· 20 (1983) 2036 2019 2039 TA45 30· 20 (2000) 2021 2010 2032 TA46 30· 20 (2015) 2047 2041 2070 TA47 30· 20 (1903) 1938 1935 1958 TA48 30· 20 (1949) 2001 1996 1994 2022 TA49 30· 20 (1967) 2013 1992 2015 TA50 30· 20 (1926) 1975 1975 1998 TA51 50· 15 2760 2760 2760 2760

(15)

Balas, & Zawack, 1988; Balas & Vazacopoulos, 1998), Tabu Search (Nowicki & Smutnicki, 1996; Pezzella &

Merelli, 2000; Sun et al., 1995), and Genetic Algorithm (Gonc¸alves et al., 2005; Wang & Zheng, 2001).

The computational results of FT and LA test problems are shown inTable 2. The results show that the preference list-based PSO we proposed is much better than the original design, the priority-based PSO. Since the number of instances tested by each method is diﬀerent, we cannot compare the result by average gap directly. Nevertheless, the result obtained by HPSO is better then other algorithms that tested all of the 43 instances, and the HPSO obtained the best-known solution for 41 of the 43 instances.

Table 3shows the average computation time on FT and LA test problems in CPU seconds. The

‘best-so-lution time’ is the average time that the algorithm takes to first reach the final best so‘best-so-lution, and the ‘total time’ is the average total computation time that the algorithm takes during a run. In HPSO, there is about 99% computation time spent on local search process. As mentioned in Section 3.5, the solution quality of tabu search (Nowicki & Smutnicki, 1996) is mainly affected by its initial solution, and the main purpose of the PSO process is to provide good and diverse initial solutions to tabu search. Therefore, the computational results show that the hybrid method, HPSO, performs better than both TSAB and PSO, and its average gap is 0.356% less than PSO.

Table 4 (continued)

Problem Size

(n· m)

Optimal solution (or upper bound)

TSAB TSSB HPSO Nowicki and Smutnicki (1996) Pezzella and Merelli (2000)

Best solution Average

TA52 50· 15 2756 2756 2756 2758 TA53 50· 15 2717 2717 2717 2717 TA54 50· 15 2839 2839 2839 2840 TA55 50· 15 2679 2679 2684 2679 2694 TA56 50· 15 2781 2781 2781 2785 TA57 50· 15 2943 2943 2943 2943 TA58 50· 15 2885 2885 2885 2885 TA59 50· 15 2655 2655 2655 2666 TA60 50· 15 2723 2723 2723 2732 TA61 50· 20 2868 2868 2868 2868 2896 TA62 50· 20 2869 2902 2942 2930 2958 TA63 50· 20 2755 2755 2755 2755 2774 TA64 50· 20 2702 2702 2702 2702 2718 TA65 50· 20 2725 2725 2725 2735 2759 TA66 50· 20 2845 2845 2845 2848 2869 TA67 50· 20 2825 2841 2865 2840 2861 TA68 50· 20 2784 2784 2784 2784 2802 TA69 50· 20 3071 3071 3071 3071 3096 TA70 50· 20 2995 2995 2995 3005 3041 TA71 100· 20 5464 5464 5519 5595 TA72 100· 20 5181 5181 5211 5305 TA73 100· 20 5568 5568 5581 5655 TA74 100· 20 5339 5339 5355 5412 TA75 100· 20 5392 5392 5466 5563 TA76 100· 20 5342 5342 5396 5504 TA77 100· 20 5436 5436 5444 5493 TA78 100· 20 5394 5394 5394 5476 TA79 100· 20 5358 5358 5363 5434 TA80 100· 20 5183 5183 5183 5209 5364 Average Gap 0.7792% 0.8122% 0.5659% 1.4651% # of instance 33 80 80 # of BKS obtained 12 31 27

(16)

We further tested HPSO on TA test problems (Taillard, 1993). The computational results are shown in

Table 4, and we particularly compared HPSO with TSSB (Pezzella & Merelli, 2000) inTable 5. Since the

max-imum computation time of TSSB is about 3· 104s and our machine is about ten times faster then TSSB (

Pezz-ella & Merelli, 2000), we limited the maximum computation time of HPSO in 3· 103s. As mentioned above,

99% of the computation time is spent on the local search process in HPSO. Therefore, we do not reduce the computation time by decreasing the number of iterations, but decreasing the percentage of particles that per-form a local search procedure. The HPSO will also be terminated after 103iterations, but there are only 34.6% of particles randomly chosen to perform the local search procedure in each iteration on TA51 to TA60 test problems, 26.6% on TA61 to TA70 test problems, and 6.4% on TA71 to TA80 test problems.

Table 5shows the comparison with TSSB (Pezzella & Merelli, 2000). The HPSO performs better than TSSB

on 7 of 8 problem sizes, and only worse than TSSB on the 100· 20 problem size. In the 100 · 20 problem sizes, the ﬁnal best solutions are obtained after 890 iterations of the average (so the best-solution time is very close to the total time). Since the HPSO only performs 103iterations for each run, it shows that the particles of HPSO did not converge in 103iterations, and can further improve the solutions by increasing the maximum iteration. However, since we want to compare HPSO with TSSB, we do not consider increasing the maximum iteration because it takes too much computation time.

5. Conclusions

We have presented a hybrid particle swarm optimization (HPSO) for job shop scheduling problems in this paper. We modiﬁed the representation of particle position, particle movement, and particle velocity to better suit it for JSP. We also applied Tabu Search to improve solution quality. The computational results show that HPSO can obtain better solutions than other methods.

For further research, if the HPSO we proposed is implemented to other sequential ordering problems, there are two aspects for discussion: (1) Modify particle position representation for better suitability to the problem. In the original PSO design, the particles search solutions in a continuous solution space. Although most sequential ordering problems can be represented by the priority-based representation, it may not suit the sequential ordering problems that we illustrated in Section3.1.1. Preference list-based representation or other representations will better suit the algorithm for sequential ordering problems. (2) Design other particle move-ment methods and particle velocity for the modiﬁed particle position representation. Besides, which particle movement method or particle velocity is better could be a further research topic.

Appendix A

A pseudo code of the HPSO for JSP is given below: Table 5

Comparison with TSSB (Pezzella & Merelli, 2000) on TA test problems

Problem Size (n· m) TSSBa _HPSOb

Average gap (%) Total time Average gap (%) Time to get best solution Total time

TA01-10 15· 15 0.4502 2175 0.0726 99 514 TA11-20 20· 15 1.1537 2526 0.5023 345 855 TA21-30 20· 20 1.0840 34910 0.7029 401 1238 TA31-40 30· 15 1.4475 14133 0.7654 1185 2026 TA41-50 30· 20 1.9474 11512 1.6133 1734 2769 TA51-60 50· 15 0.0187 421 0.0000 565 2909c TA61-70 50· 20 0.3960 6342 0.3463 2322 2862c TA71-80 100· 20 0.0000 231 0.5244 2797 3137c

Total average gap 0.8122 0.5659

a _{Run on a Pentium 133 MHz PC.} b

Run on an AMD Athlon 1700 + PC.

c

(17)

initialize a population of particles with random positions. for each particle k do

apply G&T algorithm to decode Xk(the position of particle k) into a schedule Sk. set the kth pbest solution (pbestk) equal to Sk, pbestk Sk.

end for

set gbest solution equal to the best pbestk. repeat

update velocities according toFig. 2. for each particle k do

move particle k according toFig. 4.

apply G&T algorithm to decode xkinto Sk.

update pbest solutions and gbest solution according toFig. 5. apply tabu search on Sk.

update pbest solutions and gbest solution according toFig. 5. end for

until maximum iterations is attained.

References

Adams, J., Balas, E., & Zawack, D. (1988). The shifting bottleneck procedure for job shop scheduling. Management Science, 34(3), 391–401.

Balas, E., & Vazacopoulos, A. (1998). Guided local search with shifting bottleneck for job shop scheduling. Management Science, 44(2), 262–275.

Bean, J. (1994). Genetic algorithms and random keys for sequencing and optimization. Operations Research Society of America (ORSA) Journal on Computing, 6, 154–160.

Beasley, J. E. (1990). OR-Library: Distributing test problems by electronic mail. Journal of the Operational Research Society, 14, 1069–1072.

Davis, L. (1985). Job shop scheduling with genetic algorithm. In J. J. Grefenstette (Ed.), Proceedings of the ﬁrst international conference on genetic algorithms (pp. 140–163). Hillsdale, NJ: Lawrence Erlbaum Associates.

Fisher, H., & Thompson, G. L. (1963). Industrial scheduling. Englewood Cliﬀs, NJ: Prentice-Hall.

French, S. (1982). Sequencing and scheduling: An introduction to the mathematics of the job-shop. UK: Horwood.

Garey, M. R., Johnson, D. S., & Sethi, R. (1976). The complexity of ﬂowshop and jobshop scheduling. Mathematics of Operations Research, 1, 117–129.

Gen, M., & Cheng, R. (1997). Genetic algorithms and engineering design. New York: Wiley.

Giﬄer, J., & Thompson, G. L. (1960). Algorithms for solving production scheduling problems. Operations Research, 8, 487–503. Glover, F. (1986). Future paths for integer programming and links to artiﬁcial intelligence. Computers & Operations Research, 13, 533–549. Glover, F. (1989). Tabu search: Part I. ORSA Journal on Computing, 1, 190–206.

Glover, F. (1990). Tabu search: Part II. ORSA Journal on Computing, 2, 4–32.

Gonc¸alves, J. F., Mendes, J. J. M., & Resende, M. G. C. (2005). A hybrid genetic algorithm for the job shop scheduling problem. European Journal of Operational Research, 167(1), 77–95.

Jain, A. S., Rangaswamy, B., & Meeran, S. (2000). New and ‘‘stronger’’ job-shop neighbourhoods: A focus on the method of Nowicki and Smutnicki (1996). Journal of Heuristics, 6, 457–480.

Kennedy, J., & Eberhart, R. C. (1995). Particle swarm optimization. Proceedings of the 1995 IEEE international conference on neural networks (Vol. 4, pp. 1942–1948). Piscataway, NJ: IEEE Press.

Kobayashi, S., Ono, I., & Yamamura, M. (1995). An eﬃcient genetic algorithm for job shop scheduling problems. In L. J. Eshelman (Ed.), Proceedings of the sixth international conference on genetic algorithms (pp. 506–511). San Francisco, CA: Morgan Kaufman Publishers. Lawrence, S., (1984). Resource constrained project scheduling: An experimental investigation of heuristic scheduling techniques. Graduate

School of Industrial Administration (GSIA), Carnegie Mellon University, Pittsburgh, PA.

Lourenc¸o, H. R. (1995). Local optimization and the job-shop scheduling problem. European Journal of Operational Research, 83, 347–364.

Nowicki, E., & Smutnicki, C. (1996). A fast taboo search algorithm for the job shop problem. Management Science, 42(6), 797–813. Pezzella, F., & Merelli, E. (2000). A tabu search method guided by shifting bottleneck for the job shop scheduling problem. European

Journal of Operational Research, 120(2), 297–310.

Shi, Y., & Eberhart, R. C. (1998a). Parameter selection in particle swarm optimization. In V. W. Porto, N. Saravanan, D. Waagen, & A. E. Eiben (Eds.), Proceedings of the 7th international conference on evolutionary programming (pp. 591–600). New York: Springer. Shi, Y., & Eberhart, R. C. (1998b). A modiﬁed particle swarm optimizer. In D. Fogel (Ed.), Proceedings of the 1998 IEEE international

(18)

Sun, D., Batta, R., & Lin, L. (1995). Eﬀective job shop scheduling through active chain manipulation. Computers & Operations Research, 22(2), 159–172.

Taillard, E. D. (1993). Benchmarks for basic scheduling problems. European Journal of Operational Research, 64, 278–285.

Wang, L., & Zheng, D. (2001). An eﬀective hybrid optimization strategy for job-shop scheduling problems. Computers & Operations Research, 28, 585–596.

Zhang, H., Li, X., Li, H., & Huang, F. (2005). Particle swarm optimization-based schemes for resource-constrained project scheduling. Automation in Construction, 14, 393–404.