Constraint-based ACO for a shared resource constrained scheduling problem

Int. J. Production Economics 141 (2013) 230–242

Contents lists available at SciVerse ScienceDirect

Int. J. Production Economics journal homepage: www.elsevier.com/locate/ijpe

Constraint-based ACO for a shared resource constrained scheduling problem Dhananjay Thiruvady a,n, Gaurav Singh b, Andreas T. Ernst b, Bernd Meyer a a b

Clayton School of Information Technology, Monash University, Australia CSIRO Mathematics, Informatics and Statistics, Australia

a r t i c l e i n f o

a b s t r a c t

Article history: Received 2 December 2011 Accepted 9 June 2012 Available online 19 June 2012

We consider a scheduling problem arising in the mining industry. Ore from several mining sites must be transferred to ports to be loaded on ships in a timely manner. In doing so, several constraints must be met which involve transporting the ore and deadlines. These deadlines are two-fold: there is a preferred deadline by which the ships should be loaded and there is a final deadline by which time the ships must be loaded. Corresponding to the two types of deadlines, each task is associated with a soft and hard due time. The objective is to minimize the cumulative tardiness, measured using the soft due times, across all tasks. This problem can be formulated as a resource constrained job scheduling problem where several tasks must be scheduled on multiple machines satisfying precedence and resource constraints and an objective to minimize total weighted tardiness. For this problem we present hybrids of ant colony optimization, Beam search and constraint programming. These algorithms have previously shown to be effective on similar tightly-constrained combinatorial optimization problems. We show that the hybrid involving all three algorithms provides the best solutions, particularly with respect to feasibility. We also investigate alternative estimates for guiding the Beam search component of our algorithms and show that stochastic sampling is the most effective. & 2012 Elsevier B.V. All rights reserved.

Keywords: Scheduling ACO Beam search Constraint programming Stochastic Sampling

1. Introduction Australia is one of the leading producers and exporters of minerals in the world and not surprisingly, the mining industry contributes heavily to the Australian economy. In particular, the industry contributes $ 121 billion a year to the economy. In terms of export income, it generates $ 138 billion per annum, which represents over half (54%) of total goods and services. The mining industry employs 187,400 people directly, and a further 599,680 in support industries1 Australia not only exports raw materials but is also the world’s top supplier of mining technologies and services. At least 60% of the world’s mines operate with Australian-made and designed software.2 This paper considers a scheduling problem motivated from the Australian mining industry and can be described as follows. Each mine has several mineral processing jobs that need to be completed

n

Corresponding author. E-mail addresses: [email protected] (D. Thiruvady), [email protected] (G. Singh), [email protected] (A.T. Ernst), [email protected] (B. Meyer). 1 Minerals Council of Australia—2011-2012 Pre-budget submission (www. mineralscouncil.com.au/file_upload/files/submissions/MCA_Pre%20Budget_FINAL.pdf). 2 Australian Trade Commission (www.utsavaustralia.in/sectoral-capability/ australia-mining-industry.aspx). 0925-5273/$ - see front matter & 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.ijpe.2012.06.012

in order to meet demand for the ships arriving in the future. The production capacity and maintenance requirements at the mine site enforce an earliest date when these jobs can be started and ship arrivals at different periods enforces precedence relations between the jobs. Moreover, during the execution each job also requires a certain quantity of additional resources (e.g., workforce, water, electricity, train wagons or trucks). On the other hand, this resource is of limited supply and therefore must be shared by all the mines. The loading of ships must be completed by a preferred time and thus the problem may be formulated as one where we aim to minimize tardiness. Furthermore, in order to meet a key performance indicator (KPI) the ships must not be delayed more than a certain number of periods beyond the preferred time. This enforces a hard deadline on the completion time of every job. In what follows, this problem can be presented as a resource constrained scheduling problem where every machine has a collection of jobs that need to be processed. The machines can only process one job at a time and each job also requires an additional resource in order to complete. As there is limited availability of this common resource and therefore at any single time point the amount of this resource consumed cannot exceed its maximum availability. This type of resource has been referred to as a renewable resource (Brucker et al., 1999). Section 2 presents a formal description of the problem. Even though the original motivation for the problem considered in this paper is mining industry, it is equally applicable to other supply

D. Thiruvady et al. / Int. J. Production Economics 141 (2013) 230–242

networks where several players need to share a common resource in order to process their jobs. Resource constrained job scheduling (RCJS) has been widely researched in the last 30–40 years with one of the first instances being Johnson’s PhD thesis (Johnson, 1967). The focus of this research has mainly been to minimize makespan as the objective where the aim is to minimize the largest completion time of any task. However, recently minimizing tardiness has also been given some attention (Ballestin and Trautmann, 2008; Singh and Ernst, 2011; Singh and Weiskircher, 2008, 2011). In Ballestin and Trautmann (2008), the authors present a population-based iterated local search method for a resource constrained project scheduling problem where tasks have prescribed completion times and deviation from these times across all tasks is required to be minimized. In Singh and Ernst (2011), the authors tackle a problem where jobs do not have hard deadlines and need to be scheduled on multiple machines to minimize total weighted tardiness subject to release dates, precedence and resource constraints. For this problem they compare a Lagrangian relaxation-based approach with simulated annealing. In Singh and Weiskircher (2008, 2011), the authors also consider problems consisting of multiple machines with resource constraints and minimizing total weighted tardiness as the objective but under the assumption that the data are decentralized. Both these studies approach the problem with agent-based systems and present algorithms that allow agents to find good schedules with minimal information sharing. There are several variants associated with the problem described above (Brucker et al., 1999; Demeulemeester and Herroelen, 2002; Neumann et al., 2003) which occur in project scheduling. Brucker et al. (1999) provide a classification of project scheduling problem where they describe problems consisting of tasks, shared resources, precedences between tasks and deadlines. Similarly, Demeulemeester and Herroelen (2002) describe various project scheduling problems and methods for solving them. A variant of project scheduling with time windows is also examined by Neumann et al. (2003). We investigate methods based on ant colony optimization (ACO) (Dorigo and Stutzle, 2004), constraint programming (CP) (Marriott and Stuckey, 1998) and Beam search (Pinedo, 2005) to determine if they are effective when minimizing tardiness. Recently, hybrid methods have become increasingly popular with the aim of making use of their relative advantages to solve problems more effectively than the independent methods on their own. While ACO has been applied to project scheduling problems (Chen et al., 2010) and so has CP (Artigues and Roubellat, 2000) independently, hybrids of these methods are likely to be more effective given the nature of the problem (tightly constrained with large search spaces). A hybrid of ACO and CP (CP–ACO) was first shown to be successful on a single machine job scheduling problem by Meyer and Ernst (2004). Since this study, Khichane et al. (2008) applied a similar algorithm to the car sequencing problem. ACO and Beam search (Beam-ACO) have also shown to be effective on a similar open shop scheduling problem (Blum, 2005). Combining these methods, Thiruvady et al. (2009, 2011) showed that a hybrid of ACO, CP and Beam search, CP– Beam–ACO can provide significant improvements over CP–ACO or Beam–ACO. We also implement the simulated annealing algorithm of Singh and Ernst (2011) to determine its effectiveness on the problem considered here. In this study, we investigate these hybrids on a modified version of the RCJS problem suggested by Singh and Ernst (2011). They show their problem to be NP-hard and we consider further complicating due time constraints. The results on this problem show that CP–Beam–ACO provides significant advantages, particularly with respect to feasibility, over the other

231

approaches. Additionally, we also conduct experiments with non-CP approaches using plain ACO and Beam–ACO (Blum, 2005) to compare their performance with the CP-based hybrids. We find that ACO on its own is effective but Beam–ACO does not provide an advantage over ACO. This paper is organized as follows. The problem is stated formally in the form of a resource constrained jobs scheduling problem in Section 2. Section 3 describes the methods and algorithms including ACO, CP and their integrations. In Section 4, details of the experiments conducted are provided followed by an analysis of the results in Section 5. Section 6 discusses the results in depth and this paper concludes with Section 7.

2. Problem definition The resource constrained job scheduling (RCJS) problem can be formally defined as follows. A number of machines, M ¼ f1, . . . ,lg must process a number of jobs J ¼ fj1 , . . . ,jn g. Each job i is associated with the following data:

 ri: earliest start time of the job,  pi: processing time of the job, i.e., the number of time units the job will spend on its machine (i.e., mine),

 di: due time of the job which is the desired time (i.e., when ships would prefer to be loaded),

 d i : hard due time of a job (to meet KPIs),  wi: weight of the job (tardiness penalty for completing after di),

 gi: resource used by the job when executing on a machine (e.g., electricity),

 mi: machine (mine) that the job must be scheduled on. Each job may only be scheduled after its earliest time and it must only be executed on its machine. Each machine is capable of executing only one job at a time. There is no setup time between jobs and once scheduled, no job may be stopped before it completes (i.e., no preemption). In addition to the resource constraints precedence constraints are also enforced between jobs (due to the arrival times of ships) which need to be scheduled on same machine. The set of precedences is denoted by PR. For two jobs i,j A J if i precedes j (i-j), then mi ¼mj and j may only start after i completes. Except for wi all data are integral. A sequence of all jobs is specified by a permutation of the jobs p. Let sðÞ be a mapping from sequences to schedules (see Fig. 2 in Section 3.2). A feasible schedule of p, sðpÞ, assigns start times (S ¼ fs1 , . . . ,sn g) and end times (C ¼ fc1 , . . . ,cn g) to all the jobs such that si Z r i , si þ pi r d i and sj Z si þ pi ¼ ci if i-j. Let Pt be the set of jobs either starting at time t or being processed at time t Pt ¼ fj9sj r t osj þpj ,j A J g

ð1Þ

sðpÞ is considered resource feasible if 8t

X

g k rG

ð2Þ

k A Pt

where gk is the amount of resource required (e.g., electricity or train wagons) by the job scheduled at time point t and G is the maximum amount of resource available across all the machines. This equation specifies that resources requirements of all jobs executing at the same time must not exceed the available resource. The objective is to minimize the total weighted-tardiness of a resource feasible schedule sðpÞ TðsðpÞÞ ¼

n X

wpi  Tðsðpi ÞÞ

i¼1

where Tðsðpi ÞÞ is the tardiness of the job pi , maxðcpi dpi ,0Þ.

ð3Þ

232

D. Thiruvady et al. / Int. J. Production Economics 141 (2013) 230–242

3. Methods In this section, we detail the methods investigated in this study including ACO, CP, CP–ACO and CP–Beam–ACO. We also provide details of the original simulated annealing algorithm. 3.1. Ant colony optimization

ConstructPermutation(): The permutation pj of jobs for the jth solution is constructed by selecting a job for each variable. A solution is considered complete when all the variables have unique jobs assigned to them. The selection of a job proceeds as follows. We generate a random number q A ð0; 1 and compare this with a pre-defined parameter q0. For the ith variable, pi , if q 4q0 , a job k for variable i in p is greedily selected according to k¼

Inspired by real ant colonies, ACO is a meta-heuristic for combinatorial optimization originally suggested by Dorigo (1992). ACO is derived from the foraging behavior of ants where they leave their nests looking for food and mark the successful paths to food sources with pheromone. Other ants which go out looking for food will probabilistically follow the paths based on the levels of pheromone on the paths. Thus, paths with stronger pheromone trails are traversed by more ants and these trails in turn receive more pheromone. This crates a positive feedback loop that allows the colony to converge on better food sources over-time (Camazine et al., 2001). In the context of the problem at hand, we specify a model suggested by den Besten et al. (2000). This study examines an ACO algorithm for the single machine problem with the total weight tardiness objective. Hence, the ACO model in this study is a plausible model to use here.3 The pheromones T consist of pheromone values tij for each job j and variable i. Algorithm 1 shows the ACS implementation for RCJS. Algorithm 1. ACS for the RCJS problem. 1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13: 14: 15: 16: 17: 18: 19: 20: 21: 22:

INPUT:

A RCJS instance, na (number of solutions)

pbs : ¼ null (global best) initialize T while termination conditions not satisfied do Siter : ¼ | for j¼ 1 to na pj : ¼ ConstructPermutation() CompletePlacement(pj ) Siter : ¼ Siter [ fpj g end for // Set the iteration best to the best solution pib : ¼ argminff ðpÞ9p A Siter g

pib : ¼ LocalSearch() Update(pib , pbs ) // Update the pheromone trails with the current best solution PheromoneUpdate(T , pbs ) // Determine if the pheromones have converged cf : ¼ ComputeConvergence(pib ) // reset the pheromones if they have converged if cf ¼true then initialize T end if end while bs OUTPUT: p

We examine two variants of ACO that have been successful in practice: ant colony system (ACS) and Max2Min ant system (MMAS). Our initial results show that ACS was more successful for this problem and we hence used ACS for further integrations. We briefly describe the procedures used by this algorithm and highlight some of the differences to the usual ACS. The parameter settings will be discussed in Section 4. 3 Note that we also briefly tested other learning models which shown not to be as effective.

max

j A J \fp1 ,..., pi1 g

tij  Zj

ð4Þ

otherwise, select k (k A J \fp1 , . . . , pi1 g) from the following distribution: pðpi ¼ kÞ ¼ P

tik  Zk t  Zj Þ

j A J \fp1 ,..., pi1 g ð ij

ð5Þ

Zk is defined as wk =dk . This biases the selection of jobs with large weights and earlier due times. While there are other plausible heuristics such as wk =pk , initial tests showed that the above heuristic was the most effective for this problem. Note that the selection of a job does not take precedences into account and these details will be discussed later. Every selection of a job j to a variable i has an associated update to the pheromones tij ¼ tij  r þ tmin

ð6Þ

where r is a learning rate that is chosen so that the pheromones reduce gradually. This form of update essentially allows the algorithm to diversify to avoid getting stuck in a local optimal. tmin ¼ 0:001 is a small value that ensures that no pheromone value is too small such that it will not be considered in future solution constructions. CompletePlacement(): Once a sequence for the current solution has been specified the schedule sðpÞ is determined. This is done using a placement scheme (see Section 3.2) which satisfies precedence and resource constraints to generate a resource feasible schedule. LocalSearch(): Incorporating local search within the algorithm was motivated by the fact that the SA algorithm described previously is very effective on this problem. This algorithm essentially interleaves two neighborhood moves in various ways. Over here, we consider these moves but apply them sequentially a constant number of times. The local search is applied to the best solution found in an iteration. The two following neighborhood moves are applied sequentially for 100 iterations each:

 Random swapping: this is the most basic move and works as



follows. Two values i,j A 1, . . . ,9J 9 are selected uniformly randomly the corresponding jobs pi and pj are swapped creating a new solution–p^ . The old solution is replaced with the new one if the new solution has improved tardiness, f ðsðp^ ÞÞ o f ðsðpÞÞ ) p ¼ p^ . This move is applied a constant number of times. b-sampling: given the best solution from the above procedure, b-sampling (Valls et al., 2003) is applied to it a constant number of times. The basic idea is to select a sub-sequence of the permutation starting at index i (selected uniformly from all the indexes) and move these jobs to the end of the permutation. All subsequent jobs are moved up to i. Fig. 1 demonstrates b-sampling on a permutation of 10 jobs with sample size 4. We use a sample size of 5 similar to (Singh and Ernst, 2011).

Update(pib , pbs ): This procedure sets pbs to pib if f ðsðpib ÞÞ o f ðsðpbs ÞÞ where f ðsðpib ÞÞ is the cost of the iteration best solution (i.e., the weighted tardiness). For the algorithms that do not use CP we first minimize over the number of constraint

D. Thiruvady et al. / Int. J. Production Economics 141 (2013) 230–242

233

Fig. 1. This is demonstration of b-sampling for a permutation of 10 jobs. Using a sample size of 4 and the starting index i¼ 2, jobs 7, 4, 8 and 5 move to the end of the permutation. They maintain their relative sequence as they are shifted. Jobs of the same color are required to be placed on the same machine.

violations (nðsðpib ÞÞ o nðsðpbs ÞÞ) determined from a sequence followed by weighted tardiness. ApplyPheromoneUpdate(T , pbs ): The algorithm updates the pheromones based on the global best solution. All components (i,j) appearing in the global best solution are used to update the corresponding components in the pheromone matrix:

tij ¼ tij  r þ dpbs

ð7Þ bs

where dpbs ¼ Q =f ðp Þ and Q is determined such that 0:01 r dpbs r 0:1. This is the reward factor and ensures that for all instances the reward factor is of the same order. r is the learning rate as specified earlier and is defined to be 0.1 for this study. The pheromones in the CP based algorithms are not rewarded while a feasible solution is not found but is a global best solution is found it is always rewarded even if the current iteration finds no feasible solution. ComputeConvergence(pib ): Here, we compute a convergence measure to determine if the construction mechanism is repeatedly producing the same solution. For this purpose, a list of the past m solutions, lpib , is maintained. The current iteration best pib is compared to these solutions and if they all have the same cost the pheromones are re-initialized: f ðpib Þ ¼ f ðkÞ,kA lpib ) tij ¼ 0:5 8i,j. 3.2. Placement schemes Singh and Ernst (2011) describe the serial and parallel placement schemes for a sequence of jobs. They showed that the serial scheme is usually superior and we therefore focus on this scheme here. The serial method places jobs on their machines starting at time point t0. Each job is tested for placement at a free point on the time-line. If it satisfies the resource and precedence constraints and a sufficiently large sequence of time points are free, the job is immediately placed (see Fig. 2). Here, an example is presented for the three machine fifteen job problem. Fig. 2(a) shows the current state of the placement where a number of jobs are still left in the permutation (p) to schedule. The waiting list (p^ ) has one job at this point. This list consist of those jobs which arrived early in the sequence, but due to precedence constraints must be scheduled after jobs that are still in the sequence. All such jobs are placed in the list. After any job from the sequence is scheduled the list is examined to determine if any job may be placed. These waiting jobs are placed as soon as possible. This is observed in Fig. 2(b) where job 3 is placed on machine one (m1) and job 4 from the waiting list is placed at the earliest possible time after job 3. In this example, job 4 is placed much later than job 3 and after job 14 since job 14 is uses most of the available resource which job 4 is unable to share. The parallel placement scheme works as follows. Starting at t0 the sequence of jobs is examined to determine if any job can be placed. Once a job is placed the sequence is examined again until no possible job can be placed on any machine (due to the resource constraint) at this time point. The procedure then moves along to the next time point and attempts to place jobs in the same fashion. This procedure continues until the job list is empty. A consequence of this placement scheme is that the jobs need not be sequenced in the order of precedence by the ants. While

Fig. 2. This is an example of a three machine 15 job problem. (a) Jobs 1–5 must be placed on machine 1 (m1, green), 6–10 on machine 2 (m2, yellow) and 6–10 on machine 3 (m3, orange). p is the sequence of jobs remaining to be scheduled and p^ is the waiting list. The height of the machines reflects the total resource available at each time point and the sum of the heights of the jobs scheduled across a time point must be less than or equal to the height of any one of the machines. The region in white shows the available resource in any machine. t0 is the point at which the scheduling starts and tn is the last point on the time-line. In this example job 3 must be placed before job 4. (b) Once job 3 is placed, job 4 is placed immediately in the first available location. Note that job 4 has to wait for job 14 to complete since the resource requirement for both these jobs to execute in parallel is greater than the total available resource. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

the precedences could be maintained when constructing sequences, it was found in initial tests that the solution quality found with this procedure was worse than when the precedences were allowed to be violated and handled at the placement phase. This is despite both procedures permitting the optimal sequence (see Appendix A), although the first procedure allows a much larger number of sequences that map to the optimal solution. 3.3. Simulated annealing For the sake of completeness, we provide the SA algorithm suggested by Singh and Ernst (2011) (see Algorithm 2). The algorithm is briefly described here and we refer the author to the original paper for the complete details. We obtain an initial temperature t, the procedure for which is described in Singh and Ernst (2011). pib is set to a random sequence and the algorithm now repeats until the time limit is reached (line 3). Between lines 8 and 20 the algorithm randomly swaps jobs and accepts the new solution if it is an improvement on the original solution. However, if the new solution quality is worse, it is still accepted with some probability (line 15). Outside this loop the b-Sampling procedure is applied at line 22. These steps are repeated five times (lines 5–24) after which the best solution is mutated (line 25). Note that the solution quality is determined by two factors. The first is the number of violations of the hard constraints of a sequence

234

D. Thiruvady et al. / Int. J. Production Economics 141 (2013) 230–242

and the second is the weighted tardiness. A new solution is considered an improvement if it has a smaller violation count or an equal violation count with smaller weighted tardiness. Algorithm 2. Simulated annealing the RCJS problem. 1: 2: 3: 4: 5: 6: 7: 8: 9:

INPUT A RCJS instance initialize t0, pib ¼ RandomSequence() while termination conditions not satisfied do k ¼ 0,t ¼ t 0 while k o5 do tardiness ¼ 1 prb ¼ pib while i o 5000 do rb swap prb l and pm (l,m A f1, . . . ,9J 9g)

10: 11: 12: 13: 14: 15: 16:

if f ðsðprb ÞÞ o tardiness then tardiness ¼ f ðsðprb ÞÞ pib ¼ prb , i¼0 else D ¼ f ðsðprb ÞÞtardiness if eD=t orandomðÞ then rb swap prb l and pm end if i ¼ i þ 1, t ¼ t  0:95 end if end while if f ðsðpib ÞÞ o f ðsðpbs ÞÞ then pbs ¼ pib end if pib ¼ pbs -b Sampling(); k ¼ k þ1 end while pib ¼ Mutate(pbs ) end while bs OUTPUT: p

17: 18: 19: 20: 21: 22: 23: 24: 25: 26: 27:

Mutate(pbs ): This procedure, with equal probability, perturbs bs p in the following way:

 apply b-samplingðÞ once and return best solution is found,  uniformly randomly swap any two jobs nswaps times4 and return the best solution is found,

 return a new random list.

the variables by explicitly enforcing a constraint provided by the program or implicitly by inferring constraints and imposing these on the domains of the variables. The program can obtain information about a variable’s domain (e.g., set of remaining values) by querying the solver. An attempt to assign a concrete value to a variable is called a labeling step. Such a tentative assignment may be rejected by the constraint solver (because it is not consistent with the constraints) or it may be accepted and trigger further propagation of constraints. For further details of the CP paradigm the interested reader is referred to Artigues and Roubellat (2000). CP models for problems are enhanced if high-level constraints can be used. These constraints are typically ones for which specialized algorithms exists and can be handled efficiently. For example, powerful scheduling constraints such as cumulatives() (Gecode, 2010) can be used to effectively tackle scheduling problems. The CP model devised here uses two high-level constraints in addition to a number of low-level constraints. First, both highlevel constraints are an extension of cumulatives(). Specifically, we have for each machine m A f1, . . . ,Mg, cumulatives(sm, pm, cm) which specifies that all jobs i,j A f1, . . . ,J g, in machine m must m m m m m m m have sm j Zsi þpi ¼ ci or si Z sj þ pj ¼ cj . This constraint does not specify any relationship with the resource constraint and therefore cumulatives(s, p, c, r, g) is also used. Here, the start and completion times (s and c) given the processing times (p) are determined such that at any point in time the sum of the resources for the executing tasks (r) is less than g. The execution times are allowed to overlap if the cumulative resource used during a particular time is bounded by the total resource available. Combining the above constraints enforce that all the start times on a single machine do not overlap and also that the cumulative resources used at a particular time point are always less than the available resource. Additionally, we specify the following constraints: 8m A f1, . . . ,Mg48i A f1, . . . ,J m g : 8m A f1, . . . ,Mg48i,j A f1, . . . ,P m g :

si Z r i 4si þpi ¼ ci 4ci r d i i-j ) sj Z si þ pi

The first set of constraints specifies the relationships between the start and end times of the jobs on a machine m. In particular, the job must start after or at its release time, its completion time is the sum of its start time and processing time and its completion time has to be less than the hard due date. The second set of constraints specifies the precedence between jobs on a machine.

3.4. Constraint programming 3.5. Integrating constraint programming Constraint programming (CP) is a paradigm that was originally suggested in the context of satisfaction problems and have been extended to tackle combinatorial optimization problems (Marriott and Stuckey, 1998). Here we deal with finite domain CP where variables and their domains are restricted to be finite and integral. The basic framework underlying CP is that the problem is modeled by means of variables and their associated domains. Restrictions are attached to these domains thereby forming constraints. A constraint solver keeps track of the variables, their domains and all the constraints. The solver automatically enforces restrictions that a program has requested and in doing so updates the domains of the variables. In particular, it provides at least two services to the program: (1) when a new constraint is added, it is checked for consistency with the existing ones. If the resulting domains are consistent the solver signals this as success to the program and failure otherwise. (2) Given any restrictions, the solver automatically updates the domains of 4

nswaps is randomly chosen from f1, . . . ,ng.

Algorithm 3. Construct solutions using CP. 1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13: 14: 15: 16:

i’0, feasible ’ true while i r n & feasible do i’iþ 1 D¼ domain(pi ) repeat j¼ selectJob(D, t) if PRj in sðpÞ then feasible ¼updateJobs(pi , j, p^ ) else feasible ¼true, append j to p^ end if if not(feasible) then post(pi a j) D ¼ D\j until D a| 3 feasible end while Return p

D. Thiruvady et al. / Int. J. Production Economics 141 (2013) 230–242

The integration here is along similar lines of the integrations seen in Meyer and Ernst (2004) and Thiruvady et al. (2009, 2011). Algorithm 3 shows how a solution is constructed using CP. For each variable in the sequence, we select a job by sampling the pheromones. This selection is posted to the CP solver and we continue to select jobs if the solver state is not in failure. If a selection is inconsistent with the solver this job is discarded from the current candidate list and a new selection is made. If we are able to bind all the variables to values we produce a feasible solution. The difference here compared to Meyer and Ernst (2004) and Thiruvady et al. (2009, 2011) is that we keep a waiting list of jobs (p^ in Algorithm 3). These are jobs that may only be scheduled if their preceding jobs have been scheduled. Therefore, a selection of a job does not automatically invoke the solver. This is only the case if the job is scheduled when selected. On line 7, PRj is the set of preceding jobs of j and j may only be scheduled if its preceding list of jobs have been scheduled. If j can be scheduled, this may allow a number of other jobs to be scheduled (updateJobs()). That is, there may be jobs on the waiting list which are free to be scheduled once j has been placed. The placements for all jobs that are scheduled are posted to the solver, immediately, to determine feasibility. The result of this scheme is that a number of jobs may be sequenced without being scheduled. Since the solver is only invoked when they are placed, it may happen that the new job being placed may already be infeasible as a result of the waiting jobs. Hence, this solver potentially leaves it too late to allow jobs to be sequenced (see Section 3.2). However, given that improved sequences are obtained by not maintaining precedences while sequencing jobs, there is no possible way to post to the solver immediately. In terms of obtaining feasibility, the current scheme still provides an effective way to tackle the problem by immediately invoking the solver for most of the jobs. 3.5.1. CP–Beam–ACS Algorithm 4. Probabilistic Beam Search with CP. 1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13: 14: 15: 16: 17: 18: 19: 20: 21: 22: 23: 24:

ðy, m,T Þ // Initialize the beam with y empty solutions B0 ¼ fp1 ¼ ðÞ, . . . , py ¼ ðÞg t ¼0 while t o n and 9Bt 9 40 do for i A Bt do k’0, D ¼domain(pit þ 1 ) while k o m4D a| do ðp , p^ Þ’pi , feasible¼true j ¼selectJob(D, T ) // If the predecessors of j have been scheduled if PRj in sðp Þ then feasible¼ updateJobs(p , p^ , j) else p^ ’ append(p^ ,j) end if if feasible then Bt þ 1 ¼ Bt þ 1 [ ðp , p^ Þ k’kþ 1, D ¼ D\j end while end for Bt þ 1 ¼ ReduceðBt þ 1 , yÞ t’t þ 1 end while OUTPUT: argmaxff ðpÞ j p A Bn1 g INPUT:

CP–Beam–ACS is implemented in a similar fashion to the implementations in Thiruvady et al. (2009, 2011). At a high-level the algorithm works as follows. For the next variable in the partial

235

solution, a number of children are selected based on the pheromone information biased by a heuristic factor. The total set of solutions now includes y  m candidates. From these the best y solutions are selected in a greedy fashion using f. Note that f is also computed for a job that is sent to the waiting list. This set is further reduced to y using the estimate f. See Algorithm 4. A solution is represented by two sequences ðp , p^ Þ where the first is the current sequence of jobs and the second sequence is the waiting list of jobs. A job is selected in line 9 and this selection is posted to the solver only if its preceding jobs (PRj) have been scheduled. This is done in line 11, updateJobs(p , p^ , j). Here, the sequence p , the waiting list p^ and selected job j are passed as arguments. Once j is scheduled, all those waiting jobs that are now free to be placed are immediately scheduled. If a job is not scheduled it is placed in the waiting list (line 13). The new sequence is placed in the beam if it is either feasible or if the new selection is waiting (line 15). The algorithm then proceeds as usual by using estimates at line 19 to eliminate non-promising solutions.

4. Experimental setting The results for experiments with ACS, Beam–ACS, CP–ACS, CP– MMAS and CP–Beam–ACS are reported here. We also compare these results to the original SA algorithm. In order to select an appropriate ACO variant we first conducted experiments with CP– ACS and compared the results to CP–MMAS. As discussed earlier, this was to determine which implementation is more effective on the hard constrained version of the problem. The datasets were obtained from Singh and Ernst (2011). While their study considers a large number of instances, we restrict our experiments to a limited number of instances per machine size. The instances were selected randomly from amongst the instances in each category. The reason for this is that all the algorithms tested here are stochastic and hence several runs per instance are needed to prove their performance statistically. Conducting a large number of runs on every instance was not feasible. The parameter settings for q0 and r were obtained as follows. A problem instance from each machine size was selected and run for 30 min with a cross product of values from q0 ¼ f0:3,0:5,1:0g and r ¼ f0:1,0:01g to get an idea of which parameters were best suited to this problem. Given these tests, q0 ¼1.0 and r ¼ 0:1 were selected for ACS. There was no clear advantage for CP–MMAS and therefore q0 ¼0.5 and r ¼ 0:1 were selected. The CPACO algorithms constructed 10 solutions per iteration and for consistency, y was set to 10 for CP–Beam–ACS. m ¼ 3:0 was also found to be the best multiplier after tuning by hand. The number of samples (N s ) was set to 5. This small number was chosen since a call to the solver is relatively small for this problem and a cost effective estimate can be obtained relatively quickly if only five samples are computed per variable. Thirty runs were conducted for each application of the algorithm for 60 minutes to every instance. All experiments were conducted on the Monash Sun Grid and the Enterprisegrid using Nimrod (Abramson et al., 2000). 4.1. Hard due times The aim is test our algorithms on problems with a medium level of tightness of constraints. This is because very tightly constrained problems are easily solved by CP and loosely constrained problems would do little to test the constraint propagation part of our algorithm. Hence, we set our hard due dates in a way that we expect to result in problems that are neither too tight nor too loose.

236

D. Thiruvady et al. / Int. J. Production Economics 141 (2013) 230–242

The problem instance is setup with a constraint solver with all the basic constraints and the following steps are repeated until the solver state returned is feasible (see Algorithm 5). Initially, a job i has its hard due time set to ð9M1Þ  ðr i þpi Þ (line 3). The idea here is that the multiplier ð9M91Þ provides a substantial cushion for the job to complete without violating the constraints. This however may not be sufficient for the jobs with precedences and therefore these jobs have their hard due time set to the sum of all preceding jobs release times plus processing times (line 4). At this stage we only modify the hard due times for those jobs with precedences. This state is sent to the solver in line 5, where if it fails, the due times for the jobs with precedences are increased until a successful state is achieved (lines 6–11). Specifically, the algorithm iteratively increases the level (l) until the solver state returns true. That is, we initialize a CP solver with all the constraints discussed in Section 3.4 and update the due dates without initializing any jobs. Eventually, the due times will be sufficiently large such that the solver will return a feasible state. Now we are sure that the solver starts with a feasible state but this does not imply that a feasible solution exists. Therefore, as a final step we further increase l’l  k at line 11 and reset the hard due times. In this study, we vary k depending on the problem size, i.e., 9M9 r 4 ) k ¼ 2, 4 r9M9r 9 ) k ¼ 3 or 10 r 9M9 r12 ) k ¼ 4. Even with this arbitrary increase in due dates, we cannot guarantee that a solution to the problem instance will exist. In fact, determining if a feasible solution exists is the decision version of the optimization problem of minimizing makespan. This problem is known to be NP-complete (Demeulemeester and Herroelen, 1992) and therefore we do not attempt to determine this for every instance. Algorithm 5. Procedure to determine hard due times. 1: 2: 3: 4: 5: 6: 7: 8:

INPUT: A RCJS instance m initialize solver, l ¼ 0:0, status ¼fail

d i ¼ ð9M91Þ  ðr i þ pi Þ,8iA J P d i ’d i þ 1:1  j A PRi ðr j þ pj Þ,8i A J status ¼Solve(m) while status¼fail do l’l þ 0:1

d i ¼ di ,fi9PRi a fgg P d i ’d i þ þ1:1  j A PRi ðr j þ pj Þ,8iA J 10: status¼ Solve(m) 11: end while 12: l’l  k 13: d i ’d i þ l  ðr j þ pj Þ,8i A J ,8j A PRi 9:

14:

OUTPUT:

d^ i ,8i A J

More sophisticated methods may be used to obtain these hard due times, however, the procedure followed here suffices for the purpose of generating hard instances to solve. Moreover, by tweaking the parameter l we are also able to create problems of varying hardness.

 CP–ACS and

 

CP–Beam–ACS are the best performing algorithms: J CP–Beam–ACS is most effective with respect to feasibility, stochastic sampling providing a clear advantage. J CP–ACS is by a small margin best performing algorithm where feasibility is found. ACS is surprisingly effective, but not as competitive as CP– ACS given the same number of labeling steps. Beam–ACS preforms poorly: stochastic sampling does not provide an advantage without CP.

Overall we see that ACS, CP–ACS and CP–Beam–ACS are all effective algorithms when solving the RCJS problem. CP–ACS is effective for the hard constrained version of the problem but is slow due to the solver overhead. This is remedied by CP–Beam– ACS which is as effective as CP–ACS but is able to achieve this in a small amount of time. When considering feasibility, CP–Beam– ACS is the most effective method for solving this problem. In terms of solution quality, ACS is the most effective when solutions are easily found, otherwise CP–Beam–ACS is still the most effective algorithm. The first set of results are presented in Table 1. Here we see the comparison between the two variants of ACO and CP–Beam–ACS. The first column specifies the problem instance where the first index is the number of machines and the second index is an identifier. For each algorithm, we report the best (best) solution found across all runs, the average (mean) solution quality across the runs, the associated standard deviations (sd), the number of times the algorithm fails to find feasibility out of 30 runs as a percentage (fail), the number of CP labeling steps (steps) and the number of iterations (iter.) executed by the algorithm in the allowed time. We report CP labeling steps as opposed to total labeling steps since we compare the algorithms in terms of how many calls to the solver are made in the allowed time. The aim will be to determine if a difference in a large number of CP calls is useful in the context of this problem. The best results obtained in a table for an instance (including mean and feasibility) are marked in bold face if they statistically significant at p ¼0.05. The results show that CP–ACS is the preferred implementation over CP–MMAS for this problem. This is true when considering feasibility and solution quality. Therefore, ACS is used in the implementation of the beam search–CP variant. These results are reported in the last six columns of Table 1. The estimate used for CP–Beam–ACS was stochastic sampling ˜ ez and it is used in a similar manner as it was used by Lo´pez-Iba´n et al. (2009) and Thiruvady et al. (2009, 2011). The basic idea is start with a partial solution and complete it by sampling the pheromone trails using Eqs. (4) and (5). This is done Ns times. The samples provide two estimates to the partial solutions. First, they provide an estimate of the number of violations of the hard due times determined from the complete sequence. These samples are constructed by relaxing the hard due time constraint and the lowest number of violations among the samples is returned as the estimate. The second estimate is the minimum value of the weighted tardiness across the samples. More formally   fi ¼ mins ðnðpj ÞÞ,mins ðf ðpj Þ ð8Þ jAN

5. Results The major results of this study considering optimality and feasibility can be summarized as follows5: 5 Note: CPACO refers to a generic CP with ACO hybrid whereas the implementation of CP with a variant of ACO, e.g., ACS is referred to as CP–ACS. Similarly CP with MMAS is referred to as CP–MMAS.

jAN

for the ith variable and j A Ns . We prioritize by minimizing over nðpÞ as feasibility is necessary before we can optimize f ðpj Þ as we do with the higher level ACS search. Alternative estimates based on known bounds (obtained from the literature) were also examined and were found not to be as effective as stochastic sampling (see Appendix B). The results here show that CP–Beam–ACS is the best method for feasibility and optimality on several problem instances.

D. Thiruvady et al. / Int. J. Production Economics 141 (2013) 230–242

237

Table 1 Results for CP–ACS CP–MMAS and CP–Beam–ACS on a subset of problems selected from Singh and Ernst (2011). Statistically significant results (p ¼0.05) between the two best performing algorithms for an instance are marked in boldface. CP–MMAS Best 3–23 3–53 3–5 4–28 4–42 4–73 5–21 5–62 5–7 6–10 6–28 6–58 7–15 7–23 7–47 8–3 8–53 8–77 9–20 9–22 9–47 10–13 10–77 10–7 11–21 11–56 11–63 12–14 12–36 12–80

Mean

158.08 72.53 511.78 23.99 125.37 73.30 200.61 311.66

CP–ACS sd

158.08 72.53 511.90 25.27 151.90 86.03 201.33 348.50

0.00 0.00 0.39 0.92 10.99 10.05 1.56 19.83

1197.06 1334.96 258.16 270.49 333.27 358.54

64.74 6.61 8.60

671.39

698.86

14.92

507.26 1289.38 1039.58 1260.51 1495.98

525.65 1396.30 1117.15 1337.74 1579.70

11.50 40.60 35.54 35.52 42.32

1248.73 1336.53 2621.59 2927.41

43.78 218.78

Fail

Steps

Iter.

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.00 0.67 0.00 0.07 1.00 0.00 1.00 1.00 0.00 0.00 0.00 0.00 0.00 1.00 1.00 1.00 0.00 0.17 1.00 1.00 1.00 1.00

9.71E6 1.29E7 9.04E6 1.07E7 1.49E7 1.26E7 9.07E6 1.01E7

2.86E3 1.12E4 8.05E3 9.40E3 7.47E3 6.30E3 5.66E3 5.59E3

Best

Mean

sd

158.15 72.53 511.96 24.30 144.44 79.28 205.07 355.84

0.17 0.00 1.02 0.95 16.23 5.88 4.93 29.38

1.16E7 5.28E3 1194.41 1296.24 8.56E6 8.18E3 246.61 275.55 1.29E7 5.82E3 323.83 355.95

61.89 13.81 16.23

7.04E6 3.53E3

688.34

25.96

3.07E3 487.00 530.72 2.41E3 1278.31 1348.06 1.67E3 988.60 1057.67 2.45E3 1242.08 1282.66 1.08E3 1283.00 1354.77

19.47 36.08 26.83 28.56 40.21

4.91E6 5.40E6 4.50E6 4.54E6 3.27E6

158.08 72.53 511.78 23.94 123.95 73.30 194.31 316.73

CP–Beam–ACS

636.33

4.74E6 1.20E3 1172.95 3.96E6 9.28E2 2111.05 3943.91 3288.91

1228.23 32.31 2440.69 188.67 3943.91 0.00 3288.91 0.00

3481.23 3969.04

410.83

Fail

Steps

Iter.

Best

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.00 0.37 0.00 0.00 1.00 0.00 1.00 1.00 0.00 0.00 0.00 0.00 0.00 1.00 1.00 1.00 0.00 0.00 0.97 0.97 1.00 0.87

1.08E7 1.76E7 1.11E7 1.26E7 1.68E7 1.25E7 1.01E7 1.32E7

2.47E3 8.84E3 4.45E3 7.77E3 9.05E3 7.56E3 4.11E3 7.60E3

sd

158.61 72.53 515.34 24.45 149.78 78.11 209.73 357.41

0.80 0.00 2.72 0.86 20.02 3.29 6.03 18.37

1.03E7 4.75E3 1211.07 1327.93 9.50E6 5.38E3 251.68 274.67 1.22E7 4.57E3 332.70 364.20

93.13 12.93 16.37

9.47E6 4.51E3

5.67E6 5.75E6 4.53E6 4.08E6 3.19E6

3.75E3 2.50E3 2.21E3 2.36E3 2.07E3

4.20E6 5.13E6 5.68E6 7.47E6

2.02E3 2.17E3 1.52E3 2.07E3

158.08 72.53 511.78 23.94 123.95 73.38 200.61 335.37

Mean

626.95 1009.81 1246.56 495.99 1330.63 1094.99 1316.80 1589.74 4608.86 1835.65

1524.44 2422.26 3320.50 3697.89 5324.54 4.77E6 1.63E3 3856.62

677.38 30.14 1066.33 62.12 1564.81 213.35 537.55 16.33 1449.28 67.74 1213.07 47.43 1417.96 52.13 1673.19 52.68 5036.25 317.75 1941.00 75.59 1609.91 2663.11 3616.50 4108.67 5647.87 4338.52

52.58 116.35 222.56 159.56 222.21 222.52

Fail

Steps

Iter.

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.00 0.00 0.00 0.00 1.00 0.00 0.83 0.00 0.00 0.00 0.00 0.00 0.00 0.13 0.13 1.00 0.00 0.00 0.43 0.03 0.00 0.03

2.64E6 4.21E6 2.44E6 2.55E6 3.22E6 1.82E6 1.60E6 1.46E6

8.54E2 6.42E1 1.35E3 7.20E2 1.02E3 5.52E2 5.30E2 3.81E2

9.21E5 1.40E6 7.39E5 7.58E5 5.95E5 5.24E5 5.44E5 3.74E5 7.87E5 1.18E6

3.29E2 2.45E2 1.66E2 2.83E2 1.97E2 1.53E2 1.67E2 7.56E1 8.64E1 1.80E2

4.19E5 3.32E5 8.19E5 5.81E5 3.96E5 4.40E5

5.87E1 4.45E1 9.03E1 6.20E1 2.28E1 3.22E1

9

11

7.79E5 2.40E2 1.46E6 5.45E2 1.70E6 4.18E2

35 11 10

30

% difference to best

Cumulative % failures

9 8 7 6 5 4

25 20 CP−ACS CP−Beam−ACS

15 10

3

5

2 1

0

0 ACS

CPACS

CPBeamACS

3

4

5

6 7 8 Machines

12

Algorithms Fig. 3. Cumulative failure results for ACS, CP–ACS and CP–Beam–ACS. ACS and CP–ACS fail on several instances where their performance is similar with CP–ACS performing slightly better. For each algorithm, each block corresponds to one problem instance where the shading/colour indicates the problem instance and the height of each block indicates the number of failures. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Differences between the algorithms are also exaggerated for the larger instances with more machines and jobs. Since both algorithms use CP the improved feasibility can be attributed to the estimate provided by stochastic sampling. It is also worth noting that CP–Beam–ACS compared to CP-ACS completes far fewer

Fig. 4. Results for CP–ACS and CP–Beam–ACS for those instances where at least one feasible solution was found. For each algorithm, the % difference to the best performing algorithm is shown averaged across each problem size.

iterations. This amounts to fewer pheromone updates and slower algorithm convergence which is traded-off for a more complex single iteration with beam search. Additionally, CP–Beam–ACS completes fewer CP labeling steps. This is consistent with the number of iterations, i.e., for feasible solutions CP–Beam–ACS completes three times the number of CP labeling steps (m ¼ 3:0) per iteration but CP–ACS on average completes more than three times the number of iterations.

238

D. Thiruvady et al. / Int. J. Production Economics 141 (2013) 230–242

Table 2 Results for SA, ACS and Beam–ACS on the same problem instances. Statistically significant results (p ¼0.05) between the two best performing algorithms for an instance are marked in boldface. SA Best 3–23 3–53 3–5 4–28 4–42 4–73 5–21 5–62 5–7 6–10 6–28 6–58 7–15 7–23 7–47 8–3 8–53 8–77 9–20 9–22 9–47 10–13 10–77 10–7 11–21 11–56 11–63 12–14 12–36 12–80

ACS Mean

sd

Fail

158.08 72.53 511.78 23.94

158.08 74.31 515.05 23.94

0.00 2.02 4.58 0.00

279.24

305.42

401.42

318.94

394.56

366.94

592.60 1274.00 1410.51 1220.75 1242.53

674.02 1461.05 1443.04 1262.48 1310.90

81.04 173.95 843.39 21.35 32.72

0.00 0.00 0.00 0.00 1.00 1.00 0.90 1.00 1.00 1.00 0.83 1.00 1.00 1.00 1.00 1.00 0.07 0.33 0.90 0.00 0.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00

Best

Beam–ACS Mean

sd

158.08 72.53 511.78 23.94 123.95 73.38 194.31 307.80

158.08 72.53 511.78 24.08 146.52 79.78 202.40 365.55

0.00 0.00 0.00 0.48 18.48 7.48 4.05 53.62

1229.13 249.04 325.49

1264.23 271.85 359.39

244.31 13.92 16.63

615.26

672.12

22.23

1310.75 496.66 1259.33 998.94 1215.56 1231.44

1310.75 522.76 1337.33 1045.97 1270.19 1309.94

705.07 15.90 36.38 31.34 28.47 45.20

1106.13 2075.90

1201.74 2255.02

44.63 88.46

3039.25

3113.78

1029.90

The results are also plotted in Figs. 3 and 4. Fig. 3 plots the cumulative failures for each instance for ACS, CP–ACS and CP– Beam–ACS. This figure shows that ACS and CP–ACS fail on over 30% of the total number of runs ( 410=30) where CP–Beam– ACS only fails on about 15% of the runs. If CP–Beam–ACS fails on every run then so do the other algorithms. Fig. 4 shows the % difference to the best algorithm in terms of solution quality for all the instances where feasibility is found averaged across machine size. CP–ACS is the best performing algorithm in this regard. However, there are several instances which can not be compared since CP–ACS does not find solutions for these instances. 5.1. Investigating CP’s contribution Experiments were also conducted to determine how simulated annealing, ACS and Beam–ACS without any CP component perform on these problem instances. The results are reported in Table 2. First, considering SA, the results show that a large number of problem instances (19/30) are too hard for SA to find feasibility. Even for those problem where feasibility is found SA fails several times (e.g., 6–28 and 9–20). Surprisingly, SA even fails on some of the small instances such as 4–42 where every other algorithm finds feasibility. Clearly, the CP component provides significant feasibility advantages as expected. ACS however, performs competitively with CP–ACS. CP–ACS is superior in terms of finding feasibility on some problems (e.g., 6–10 and 7–47) but is surprisingly worse on one instance—8–3. For this instance, CP–ACS finds no feasible solution whereas ACS almost always fails but does find feasibility. On close examination it can be seen that ACS performs a much larger number of iterations and the larger number pheromone updates are able to lead the algorithm to feasibility like the CP

Fail

Iter.

Best

Mean

sd

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.00 0.87 0.00 0.07 1.00 0.00 1.00 0.97 0.00 0.00 0.00 0.00 0.00 1.00 1.00 1.00 0.00 0.00 1.00 0.93 1.00 1.00

6.76E3 9.81E3 2.90E3 1.92E4 2.78E4 1.54E4 1.83E4 2.95E4

158.08 72.53 511.78 23.94 123.95 73.30 200.41 312.47

158.08 72.53 511.78 23.98 132.86 76.83 202.47 340.13

0.00 0.00 0.00 0.14 11.92 2.95 2.68 15.84

1.62E4 1.72E4 2.22E4

1243.48 251.02 321.66

1266.64 269.25 349.56

635.18 9.79 18.50

1.64E4

628.06

659.15

23.45

1.14E4 1.59E4 1.07E4 7.95E3 8.45E3 7.74E3

1383.05 488.93 1283.86 1004.52 1211.87 1307.13

1514.91 504.04 1334.72 1071.36 1244.03 1373.18

836.09 8.59 26.37 31.86 20.29 28.72

1.08E4 7.06E3

1345.47 2422.06

1420.53 2533.90

35.22 79.88

5.39E3

Fail

Iter.

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.00 0.93 0.00 0.00 1.00 0.00 1.00 0.87 0.00 0.00 0.00 0.00 0.00 1.00 1.00 1.00 0.00 0.00 1.00 1.00 1.00 1.00

2.13E2 9.37E2 2.55E2 4.65E2 1.00E3 3.90E2 3.71E2 4.05E2 2.18E2 3.36E2 3.85E2 2.95E2 1.31E2 2.72E2 1.77E2 1.57E2 1.89E2 8.20E1

7.08E1 4.88E1

Table 3 ACS and CP–ACS run for 3000 iterations. Statistically significant results (p ¼0.05) are marked in boldface. ACS best 3–23 3–53 3–5 4–28 4–42 4–73 5–21 5–62 5–7 6–10 6–28 6–58 7–15 7–23 7–47 8–3 8–53 8–77 9–20 9–22 9–47 10–13 10–77 10–7 11–21 11–56 11–63 12–14 12–36 12–80

158.08 72.53 511.78 23.94 128.65 76.10 200.61 324.11

CP-ACS mean

sd

fail

158.14 72.53 513.54 25.79 154.25 85.47 209.16 433.97

0.15 0.00 2.61 1.83 40.51 12.51 6.20 99.49

1198.59 1198.59 255.43 286.44 368.30 385.27

885.03 16.60 78.98

632.56

696.77

31.20

1456.15 494.86 1283.76 1020.93 1230.37 1254.09

1552.45 531.40 1352.90 1053.39 1280.16 1323.17

626.05 18.40 37.35 26.72 31.29 40.33

1146.78 1229.43 2129.98 2285.47

51.16 129.19

4944.91 4944.91 2002.02

0.00 0.00 0.00 0.00 0.37 0.23 0.00 0.40 1.00 0.97 0.00 0.57 1.00 0.03 1.00 0.90 0.00 0.00 0.00 0.00 0.00 1.00 1.00 1.00 0.00 0.10 1.00 1.00 0.97 1.00

best 158.08 72.53 511.78 23.94 125.29 73.30 201.27 315.20

mean

sd

158.08 72.53 512.95 24.82 157.69 84.16 206.13 374.85

0.00 0.00 2.25 1.25 17.02 10.39 4.31 52.56

1149.24 1280.29 252.45 280.70 328.31 352.92

73.56 15.08 14.36

621.40

679.75

26.88

498.74 1293.70 992.55 1223.70 1237.84

524.02 1349.90 1047.93 1276.36 1309.68

16.47 32.07 34.52 30.53 39.15

1123.84 2155.89 3121.37 3621.91

1191.76 46.56 2286.36 110.84 3121.37 0.00 3644.99 23.08

3225.40 3459.17 151.56

fail 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.00 0.73 0.00 0.00 1.00 0.00 1.00 1.00 0.00 0.00 0.00 0.00 0.00 1.00 1.00 1.00 0.00 0.00 0.97 0.93 1.00 0.77

D. Thiruvady et al. / Int. J. Production Economics 141 (2013) 230–242

component of CP–ACS similarly does. Thus, we conducted experiments with ACS and CP–ACS given both algorithms the same number of iterations (3000 iterations). These results are presented in Table 3. Here it can be observed that CP–ACS is in fact superior on most problems where feasibility is an issue. Again, ACS finds a solution for instance 8–3 on one run while CP– ACS does not. Generally, although, both algorithms mostly fail on this instance. In comparison with the previous study for instances with eight or more machines where CP–ACS complete fewer than 3000 iteration always, CP–Beam–ACS is still more effective than CP–ACS even though the algorithm requires many fewer iterations. Beam–ACS however, does not provide any feasibility advantage over ACS. The estimate used here again is based on stochastic sampling and here the repeated pheromone updates (measured by number of iterations) are as effective in terms of feasibility. In fact, ACS finds feasible solutions on occasions for instance 12–14 whereas Beam–ACS does not find any feasible solutions for this instance. In terms of solution quality, Beam–ACS is more effective on small-medium size instances while ACS is more effective on the larger instances. However, the differences between the solution qualities are relatively small. 5.2. No learning and no local search Additional experiments were conducted to demonstrate that pheromone learning contributes to the search. For this implementation, the pheromone updates are turned off. This amounts to repeated solution construction from uniform distributions across all jobs at each variable. Following this, local search is applied as described above to the best solution found at each iteration. The results show that ACS without learning is not effective on this problem in terms of feasibility and solution quality.

239

A final experiment was to determine if the local search does in fact prove useful or not. Here, we turn off the local search altogether when applying the algorithm to all the instances. As the results show the local search component does indeed contribute to the solution quality. Feasibility is also more easily found for the small instances, however, this effect is not seen with the larger instances. Table 4 presents the results for ACS without learning (ACS-NL) and without local search (ACS-NLS). These results show that learning plays a major part in identifying feasibility and in solution quality. While (ACS-NL) finds feasibility almost always for the smaller instances, the solution quality found is poor compared to ACS with learning. The results without local search show that feasibility is not significantly affected here. Local search is able to improve feasibility for the smaller instances, e.g., instances 4–42 and 4– 73. However, this does not seem to be the case for larger instances ( Z 8 machines). In fact feasibility is found on some instances where feasibility was not found earlier with local search (e.g., 12– 36 and 12–80). This can be attributed to the number of pheromone updates due to the number of iterations which increase by one or two orders of magnitude across all problems. Hence, for larger problem instances, more frequent pheromone updates lead to greater feasibility. In terms of solution quality, the local search component almost always produces improved results ( Z 4 machines).

6. Discussion The results show that CP–Beam–ACS with stochastic sampling is the most effective algorithm finding feasibility on more instances than the other algorithms. If feasibility is found CP– ACS is marginally more effective as it is able to find higher quality

Table 4 ACS without learning and with local search. Statistically significant results (p¼ 0.05) are marked in boldface. ACS-NL Best 3–23 3–53 3–5 4–28 4–42 4–73 5–21 5–62 5–7 6–10 6–28 6–58 7–15 7–23 7–47 8–3 8–53 8–77 9–20 9–22 9–47 10–13 10–77 10–7 11–21 11–56 11–63 12–14 12–36 12–80

ACS-NLS Mean

sd

Fail

Iter.

0.00 0.00 0.00 0.00 0.27 1.00 0.00 1.00 1.00 1.00 0.00 1.00 1.00 0.37 1.00 1.00 0.00 0.00 0.00 0.00 0.00 1.00 1.00 1.00 0.03 1.00 1.00 1.00 1.00 1.00

4.27E4 4.99E3 3.19E4 3.79E4 2.97E4 0.00E0 2.78E4 0.00E0 0.00E0 0.00E0 2.65E4 0.00E0 0.00E0 1.45E4 0.00E0 0.00E0 1.33E4 1.01E4 9.79E3 6.63E3 5.63E3 0.00E0 0.00E0 0.00E0 8.49E3 0.00E0 0.00E0 0.00E0 0.00E0 0.00E0

158.59 72.53 533.53 29.35 204.60

159.73 72.53 549.79 30.09 268.35

1.03 0.00 9.29 0.59 104.92

244.86

262.08

8.07

321.60

357.51

15.64

918.66

1028.90

165.56

572.17 1610.71 1393.34 1473.97 1660.98

603.34 1694.76 1464.79 1501.93 1697.08

12.93 53.87 48.04 13.84 19.43

1679.39

1915.76

171.40

Best

Mean

sd

158.08 72.53 511.78 23.94 125.29 76.10 200.61 325.30

159.97 72.53 527.37 27.07 163.36 88.31 215.67 411.61

2.89 0.00 23.55 2.31 27.65 19.97 8.19 71.18

261.26 355.68

298.42 407.59

28.77 75.88

662.58

728.79

33.38

1496.84 510.87 1324.74 1023.57 1270.71 1313.30

1496.84 549.34 1410.66 1145.55 1328.89 1446.77

1574.60 23.78 55.73 67.44 30.11 71.65

1225.69 2105.31

1326.67 2376.21

68.97 165.86

4421.56 3246.14

4670.30 3297.06

1092.16 649.66

Fail

Iter.

0.00 0.00 0.00 0.00 0.13 0.23 0.00 0.50 1.00 1.00 0.00 0.50 1.00 0.03 1.00 0.97 0.00 0.00 0.00 0.00 0.00 1.00 1.00 1.00 0.00 0.33 1.00 1.00 0.90 0.93

1.63E5 3.84E4 1.11E5 2.12E5 3.05E5 2.66E5 1.57E5 1.56E5 0.00E0 0.00E0 1.69E5 1.80E5 0.00E0 1.47E5 0.00E0 5.64E4 1.25E5 1.23E5 1.09E5 1.10E5 9.32E4 0.00E0 0.00E0 0.00E0 1.05E5 7.52E4 0.00E0 0.00E0 7.53E4 8.44E4

240

D. Thiruvady et al. / Int. J. Production Economics 141 (2013) 230–242

solutions in the same time-frame. However, in practical situations, where ships must be loaded on time, feasibility is of critical importance and CP–Beam–ACS proves to be the most effective algorithm in this respect. Also seen in the results (especially for large problems) are that the CP labeling steps and iterations completed by CP–Beam– ACS are a lot fewer than CP–ACS. This can be explained by the cost of stochastic sampling Oðn2  mÞ where we have n variables and m samples. If a call to the CP solver is inexpensive, then as problem instances get larger, the cost of stochastic sampling is going to outweigh the cost of the solver. The question as to whether this is a useful trade-off to arises. Stochastic sampling proves to be effective for RCJS particularly with respect to feasibility. Given that the CP model is relatively simple, we see that the feasibility estimate provides a distinct advantage leading CP–Beam–ACS into feasible regions. This leads to the question of the usefulness of the CP model. CP– ACS outperforms ACS overall in terms of feasibility demonstrating the need for CP. However, stochastic sampling on its own does not provide a feasibility advantage. Thus the best combination in the context of feasibility is CP combined with stochastic sampling. We considered a relatively simple version of local search here. More complex methods such as variable neighborhood descent may be more effective than the constantly applying the two neighborhood moves suggested above. Furthermore, the simulated annealing heuristic may also be useful here. However, the aim here is to examine the effectiveness of parallel solution construction and while the local search may be improved if it is independent of the solution construction and is not expected to have any effect in this component of the algorithm (e.g., biasing the solution construction in the beam).

7. Conclusion In this study we apply ACS, CP–ACS and CP–Beam–ACS to a resource constrained multiple machine job scheduling problem. We see that CP–Beam–ACS using stochastic sampling is the most effective algorithm here finding feasible solutions where the other algorithms struggle. ACS and CP–ACS are both effective when feasible solutions are found often outperforming CP– Beam–ACS. Furthermore, the simulated annealing method suggested by Singh and Ernst (2011) is not effective for the problem with hard due times considered here. As opposed to Thiruvady et al. (2009, 2011) who examined single machine job scheduling and car sequencing, respectively, we see here that CP–Beam–ACS provides a feasibility advantage compared to an optimality advantage. This aspect of CP–Beam– ACO is particularly interesting since we see that advantages over CP–ACO may be seen in different aspects. The feasibility advantage can be accounted for by considering that stochastic sampling provides a feasibility estimate and in the absence of a strong CP model (as we see for the current problem) stochastic sampling provides guidance into feasible regions. Also seen here is that there are many fewer iterations conducted by CP–Beam–ACS compared to the other algorithms due to the cost of stochastic sampling. This is reflected in the results of solution quality where if feasible solutions are found, CP–ACS is the best algorithm. This study further validates pheromone model for sequences for such scheduling problems in addition to CP models for the same problems. It is already known that such models have been effective on similar problems independently (den Besten et al., 2000; Bauer et al., 2000). Hybridizing these models proves to be useful and furthermore incorporating the beam component further enhances the performance of the algorithm.

The form of local search used here is relatively naı¨ve, i.e., repeated application of random swapping and b-sampling. We have also shown that the local search is useful. Therefore, a more sophisticated local search, such as variable neighborhood descent, might provide improved results and this component of the algorithm warrants testing in the future. Finally, there are likely to be significant gains if CP–Beam–ACS is implemented in a parallel framework. ACO’s construction scheme lends itself such methods and has already been investigated. Beam search constructs solutions which are dependent leading to a less straight-forward parallel implementation than ACO but the basic algorithm provides potential to be parallelized. Given the cost of CP solvers, parallel implementations would provide tremendous benefits to significantly reducing the runtime of these algorithms. We are currently investigating in this direction by examining parallel implementations for machines with several cores.

Appendix A. A proof that demonstrates serial scheduling can generate an optimal schedule The following theorem shows that the serial scheduling scheme can produce the optimal solution to the RCJS problem. Theorem 1. If sn is an optimal schedule for the RCJS problem, then it is possible to construct a list of all the jobs, which produces the same schedule as sn using Serial Algorithm. Proof. Among all optimal schedules for the RCJS problem, let sn be the schedule such that the starting time of all the jobs (fs1 , . . . ,sn g) in sn is as early as possible in comparison with any other optimal schedule. In other words, there is no optimal schedule s^ n for the RCJS problem with start times fs^ 1 , . . . , s^ n g such that s^ j osj

8 jAJ n

ð9Þ

Let L ¼ fj1 , . . . ,jn g be a list of all the jobs arranged in nondecreasing order of the sj-values, with ties broken based on machine number for job j. We now apply Serial Algorithm to the list of jobs Ln and let sL be the schedule constructed by it with start times fs 1 , . . . ,s n g. Suppose the theorem is incorrect and for some j A J , s j asj . Among all such jobs let jk A J be the first job in the list Ln such that s jk asjk . Let us first suppose that s jk o sjk . Note that, job jk is scheduled at time s jk in the schedule sL , therefore the machine of job jk is not processing any other job and the resource available during the interval ½s jk ,s jk þ pjk  is at least the resource requirement for job jk. As the jobs are arranged in non-decreasing order of their sj-values and jk is the first job in this list to contradict the theorem, it is possible to start the job jk at time s jk in the schedule sn . As the job is processed only earlier it will not increase the total weighted tardiness, thus it contradicts the selection of sn as an optimal schedule, which satisfies (Dorigo and Stutzle, 2004). On the other hand, if s jk 4 sjk then as earlier, there are only two reasons for the job jk not to be scheduled at time sjk in the schedule sL . Either, the machine for job jk is busy scheduling some other job or the available resource is not sufficient for job jk to be scheduled at that point. In both cases, keeping in mind that the list L was created in non-decreasing order of sj-values and that jk was the first job in this list to contradict the theorem, it contradicts the fact that sn is a feasible schedule. Hence, there is no job j A J , such that s j a sj . &

D. Thiruvady et al. / Int. J. Production Economics 141 (2013) 230–242

241

Table B1 Experiments with bounds for the beam search component in CP–Beam–ACS; tardiness: based on the current partial solution weighted-tardiness; IMIR: infinite machine infinite resource. Statistically significant differences (p ¼ 0.05) are shown in boldface. Tardiness best 3–23 3–53 3–5 4–28 4–42 4–73 5–21 5–62 5–7 6–10 6–28 6–58 7–15 7–23 7–47 8–3 8–53 8–77 9–20 9–22 9–47 10–13 10–77 10–7 11–21 11–56 11–63 12–14 12–36 12–80

IMIR mean

sd

158.08 72.53 519.06 24.76 166.50 98.35 260.02 553.38

162.97 72.90 544.24 29.38 201.73 134.61 340.28 638.20

4.20 1.81 27.32 2.13 23.11 38.82 48.67 56.86

445.36 391.14

631.83 478.32

104.25 48.95

932.67

1076.03

93.45

638.83 1946.27 1519.47 1573.84 2065.75

715.18 2198.73 1767.55 1724.21 2504.39

48.18 186.70 147.59 94.85 265.03

1956.16 3306.15

2578.42 3756.40

331.02 326.50

fail

iter.

0.00 0.00 0.00 0.00 0.00 0.63 0.00 0.00 1.00 1.00 0.10 0.43 1.00 0.43 1.00 1.00 0.00 0.00 0.00 0.00 0.00 1.00 1.00 1.00 0.00 0.70 1.00 1.00 1.00 1.00

4.27E3 5.75E3 6.86E3 6.99E3 5.28E3 5.02E3 2.88E3 3.03E3

158.08 72.53 522.52 24.10 196.09 148.32 263.76 527.12

159.99 72.62 550.17 28.48 287.88 174.56 342.96 661.39

1.08 0.30 14.47 2.53 61.91 22.39 44.30 81.86

2.87E3 2.75E3

1455.33 369.08 430.04

1640.53 414.23 490.12

128.25 28.65 39.21

2.20E3

780.11

959.90

86.43

1.54E3 1.18E3 9.25E2 9.71E2 4.95E2

574.90

644.75

33.58

1348.10 1473.85 1793.36

1541.57 1590.51 2067.01

71.98 64.97 97.34

6.78E2 6.27E2

1780.09 3053.35

1992.25 3477.53

121.65 286.95

Appendix B. Investigating bounds Here we examine two alternative lower bounds for the RCJS problem.6 The first is a quickly obtainable bound based on the current solution weighted-tardiness from which we greedily select the best solutions (see Table B1). As these results show the bound does not drive the solutions towards feasibility or improved solution quality. This is expected since we are only considering past experience and not the potential of the remaining components of the solution. A second more interesting bound – infinite machine infinite resource (IMIR) – is based on a relaxed version of the problem. The first constraint that is relaxed is that jobs that belong to a machine must be placed on the machine. Hence, this bound effectively allows infinite machines. However, precedences between jobs are kept which requires that jobs with precedences must be placed on the same machine. The resource constraint is also relaxed. This is similar to the scheme used by Stinson et al. (1978). Like stochastic sampling, this bound provides two estimates, the number of violations of the hard due times and the weighted tardiness of the relaxed problem. However, as the results show, the estimate of the violations does not lead the partial solutions towards feasible regions. This result is not surprising since the start times of jobs in the relaxed version of the problem are potentially a lot earlier than their actual start times thereby not violating their hard due time. Therefore, this estimate can be misleading. Considering tardiness, the solution qualities are 6 These lower bounds are used to guide the solution construction by providing estimates. The values in Table B1 refer to the total weighted tardiness. No lower bounds for any instances are published here.

best

mean

sd

fail

iter.

0.00 0.00 0.00 0.00 0.70 0.80 0.23 0.87 1.00 0.50 0.00 0.23 1.00 0.13 1.00 1.00 0.00 1.00 0.00 0.00 0.00 1.00 1.00 1.00 0.00 0.57 1.00 1.00 1.00 1.00

5.83E3 7.39E3 5.13E3 6.35E3 9.09E3 3.57E3 3.34E3 3.23E3 1.17E3 1.76E3 2.65E3 1.65E3

1.33E3 7.03E2 8.08E2 5.73E2

6.18E2 5.07E2

improved compared to the previous weighted-tardiness estimate but are not comparable to stochastic sampling. Again, here the relaxation of the constraints potentially result in a misleading estimate of the weighted-tardiness leading the solutions into non-optimal areas of the search space.

References Abramson, D., Giddy, J., Kotler, L., 2000. High performance parametric modeling with Nimrod/G: killer application for the global grid? In: International Parallel and Distributed Processing Symposium (IPDPS), IEEE Computer Society, Washington, DC, USA, pp. 520–528. Artigues, C., Roubellat, F., 2000. A polynomial activity insertion algorithm in a multi-resource schedule with cumulative constraints and multiple modes. European Journal of Operational Research 127 (2), 297–316. Ballestin, F., Trautmann, N., 2008. An iterated-local-search heuristic for the resource-constrained weighted earliness-tardiness project scheduling problem. International Journal of Production Research 46, 6231–6249. Bauer, A., Bullnheimer, B., Hartl, R.F., Strauss, C., 2000. Minimizing total tardiness on a single machine using ant colony optimization. Central European Journal of Operations Research 8, 125–141. Blum, C., 2005. Beam–ACO: hybridizing ant colony optimization with beam search: an application to open shop scheduling. Computers and Operations Research 32, 1565–1591. Brucker, P., Drexl, A., Mohring, R., Neumann, K., Pesch, E., 1999. Resourceconstrained project scheduling: notation, classification, models, and methods. European Journal of Operational Research 112, 3–41. Camazine, S., Deneubourg, J.L., Franks, N.R., Sneyd, J., Theraulaz, G., Bonabeau, E., 2001. Self-Organization in Biological Systems. Princeton University Press, Princeton, NJ, USA. Chen, W.N., Zhang, J., Chung, H.S.H., Huang, z.R., Liu, O., 2010. Optimizing discounted cash flows in project scheduling - an ant colony optimization approach. IEEE Transactions on Systems, Man, and Cybernetics, Part C 40 (1), 64–77. Demeulemeester, E., Herroelen, W., 1992. A branch-and-bound procedure for the multiple resource-constrained project scheduling problem. Management Science 38 (12), 1803–1818. Demeulemeester, E., Herroelen, W., 2002. Project Scheduling: A Research Handbook. Kluwer, Boston, MA, USA.

242

D. Thiruvady et al. / Int. J. Production Economics 141 (2013) 230–242

den Besten, M., Stutzle, T., Dorigo, M., 2000. Ant colony optimization for the total weighted tardiness problem. In: Lecture Notes in Computer Science 1917. pp. 611–620. Dorigo, M., 1992. Optimization, Learning and Natural Algorithms. Ph.D. Thesis, Dip. Elettronica. ¨ Dorigo, M., Stutzle, T., 2004. Ant Colony Optimization. MIT Press, Cambridge, MA, USA. Gecode, 2010. Gecode: Generic Constraint Development Environment. Available from /http://www.gecode.orgS. Johnson, T., 1967. An Algorithm for the Resource-Constrained Project Scheduling Problem. Ph.D. Thesis. Sloan School of Management, MIT, USA. Khichane, M., Albert, P., Solnon, C., 2008. CP with ACO. In: Lecture Notes in Computer Science 5015. pp. 328–332. ˜ ez, M., Blum, C., Thiruvady, D., Ernst, A.T., Meyer, B., 2009. Beam–ACO Lo´pez-Iba´n Based on Stochastic Sampling for Makespan Optimization Concerning the TSP with Time Windows. In: Lecture Notes in Computer Science 5482. pp. 97–108. Marriott, K., Stuckey, P., 1998. Programming With Constraints. MIT Press, Cambridge, MA, USA. Meyer, B., Ernst, A., 2004. Integrating ACO and Constraint Propagation. In: Lecture Notes in Computer Science: Ant Colony, Optimization and Swarm Intelligence 3172/2004. pp. 166–177. Neumann, K., Schwindt, C., Zimmermann, J., 2003. Project Scheduling with Time Windows and Scarce Resources. Springer, Berlin, Germany.

Pinedo, M.L., 2005. Planning and Scheduling in Manufacturing and Services. Springer, New York, USA. Singh, G., Ernst, A.T., 2011. Resource constraint scheduling with a fractional shared resource. Operations Research Letters 39 (5), 363–368. Singh, G., Weiskircher, R., 2008. Collaborative resource constraint scheduling with a fractional shared resource. In: 2008 IEEE/WIC/ACM International Conference on Web Intelligence and Intelligent Agent Technology, vol. 2. IEEE, pp. 359–365. Singh, G., Weiskircher, R., 2011. A multi-agent system for decentralised fractional shared resource constraint scheduling. Web Intelligence and Agent Systems 9 (2), 99–108. Stinson, J.P., Davis, E.W., Khumawala, B.M., 1978. Multiple resource-constrained scheduling using branch and bound. AIIE Transactions 10 (3), 252–259. Thiruvady, D., Blum, C., Meyer, B., Ernst, A.T., 2009. Hybridizing Beam–ACO with constraint programming for single machine job scheduling. In: Lecture Notes in Computer Science 5818. pp. 30–44. Thiruvady, D.R., Meyer, B., Ernst, A.T., 2011. Car sequencing with constraint-based ACO. In: Proceedings of the 13th Annual Conference on Genetic and Evolutionary Computation. GECCO ’11. ACM, New York, USA, pp. 163–170. ISBN 978-1-4503-0557-0. Valls, V., Quintanilla, S., Ballestin, F., 2003. Resource-constrained project scheduling: a critical activity reordering heuristic. European Journal of Operational Research 149, 282–301.