Proactive Algorithms for Scheduling with Probabilistic Durations

1 Introduction

Classical models of scheduling assume all information is known and static. When faced with uncertainty, proactive scheduling techniques seek to incorporate models of the uncertainty in the off-line problem solving and to build solutions that can achieve a satisfactory level of performance at execution time. [Beck and Wilson, 2004] addressed the problem of job shop scheduling with probabilistic durations with two styles of algorithm both using Monte Carlo simulation to evaluate the probabilistic makespan of solutions. This paper builds on that work, making the following contributions: (i) A novel constraint-based search technique is introduced that performs a series of complete branch-and-bound searches on highly constrained problem models. The models are gradually relaxed with the cost of the best solution model found in the previous searches being used as an upper bound. Experimental results show that the algorithm performs as well as the existing branch-and-bound style of search on small problems and signi�cantly out-performs it on larger problems. (ii) A number of novel deterministic �ltering algorithms are presented. On larger problems such algorithms out-perform the other algorithms tested. Their performance is affected by the following two factors: the quality of the deterministic solutions found and the correlation between the deterministic and

∗This material is based upon works supported by the Science Foundation Ireland under Grant No. 00/PI.1/C075 and by ILOG, SA.

Nic Wilson

Cork Constraint Computation Centre Department of Computer Science University College Cork, Ireland n.wilson@4c.ucc.ie

probabilistic makespan of solutions. (iii) It is demonstrated that the method introduced by Beck & Wilson for the incorporation of uncertainty in deterministic problems can increase the correlation between deterministic and probabilistic makespan, providing an explanation for its strong performance.

In the next section we brie�y review previous work. The six search algorithms investigated in this paper are de�ned in Section 3 and Section 4 presents our experiments and results. In Section 5, we present our conclusions.

2 Background

The job shop scheduling problem (JSP) consists of a set A of activities each with a positive integer-valued duration, di. A is partitioned into jobs, and with each job is associated a total ordering on that set of activities. Each activity speci�es a resource on which it must execute without interruption. No activities that require the same resource can overlap in their execution. We represent this formally by a partition of A into resource sets. A solution consists of a total ordering on each resource set such that the union of the resource and job orderings is an acyclic relation on A. The makespan of a solution is the difference between the minimum start time and the maximum end time of all activities. Finding a solution with minimum makespan is NP-hard [Garey and Johnson, 1979]. An independent probabilistic job shop scheduling problem is a JSP where the duration di associated with an activity Ai ∈ A is a positive integer-valued random variable. These random variables are fully independent. di has expected value µi = E[di] and variance σ2 i = Var[di]. The makespan make(s) of solution s is therefore also a random variable. We �x a degree of con�dence p, e.g., we set p = 0.95 in the experiments. The minimum probabilistic makespan, D(s), of a solution s is the smallest value D such that the probability that all jobs will �nish before D is at least p. That is, D(s) = min{D : Pr(make(s) ≤ D) ≥ p}. An optimal probabilistic makespan is the minimumD(s) over all solutions s.

Beck & Wilson introduce an analytical model and develop two main results. First, it is demonstrated that Monte Carlo simulation can be used to �nd the probabilistic makespan for a solution and to �nd a lower bound on the probabilistic makespan of a partial solution. Second, it is shown that a deterministic JSP instance can be constructed from a probabilistic JSP instance by setting the duration di of each activity to µi+qσi and that, for the appropriate choice of non-negative q, the optimal makespan of the deterministic instance is a lower bound on the optimal probabilistic makespan of the corresponding probabilistic JSP. Three algorithms are de�ned: a branch-and-bound approach where Monte Carlo simulation is used at each node to �nd a lower bound (at internal nodes) or an upper bound (at a leaf node) on the probabilistic makespan, and two deterministic �ltering algorithm consisting of using either constraint-based search or tabu search to �nd a series of increasingly good deterministic solutions, each of which is simulated. Empirical results indicated that the �rst algorithm was able to �nd optimal solutions only for very small problems, that for medium-sized problems the constraint-based �ltering algorithm found the best solutions, and for the largest problems the tabu-based �ltering algorithm performed best. Choosing q > 0 led to stronger performance of the �ltering algorithms. No empirically founded explanations were provided to explain the differences in problem solving performance.

3 Algorithms for Probabilistic JSP

In this section, we describe six algorithms for solving the probabilistic job shop scheduling problem.

3.1 Approximately Complete Branch-and-Bound

We implemented a branch-and-bound algorithm identical to that of Beck & Wilson except that we use stronger heuristics that make decisions on the sequence of activities in each resource set (see Section 3.5). We refer to this algorithm as the B&B-N as it performs Branch-and-Bound with simulation at each Node.

3.2 Iterative Descending Bound Search

By setting activity durations based on a q value that ensures a lower bound on the probabilistic makespan, we can use standard constraint propagation on the deterministic durations to implicitly calculate a lower bound at each internal node. At each leaf node, simulation is used as in B&B-N to �nd the probabilistic makespan. Rather than parameterize the algorithm with a �xed q value, we perform repeated tree searches with a descending q value. Starting with a high q value, we begin a tree search. When the �rst leaf, e, is reached, simulation is used to �nd D(se). It is likely that Dq(se), the minimum makespan of se in the corresponding deterministic problem, will be larger than D(se). Since we enforce Dq(se) ≤ D(se), estimating D(se) through simulation causes the search to backtrack to a interior node, i, very high in the tree. At node i, Dq(Si) ≤ D(se): Si is the set of solutions in the subtree below node i and Dq(Si) is a lower bound on the deterministic makespan of those solutions based on standard constraint propagation. With high q values, we commonly observe that there are very few nodes that meet this criterion and the search space is quickly exhausted. We then decrement the q value (e.g., by 0.05) and start a new tree search. Eventually, and often very quickly, we reach a q value such that there exists a full solution where Dq(se) ≤ D(se). The probabilistic makespan of the best solution that has been

B&B-DQ-L():

Returns the solution with lowest probabilistic makespan

1 (s∗ ,D∗ ) ←�ndFirstB&BSimLeaves(∞, 0) 2 q ← q0 3 while q ≥ 0 AND not timed-out do 4 (s, D) ← �ndOptB&BSimLeaves(D∗, q) 5 if s ̸= NIL then 6 s∗ ← s; D∗ ← D end 7 q ← q - qdec end

8 return s∗

found by any point in the overall search is used, as in B&B-N, as an upper bound on subsequent search.

Algorithm 1 presents pseudocode for the basic algorithm. We make use of two functions not de�ned using pseudocode: �ndFirstB&BSimLeaves(c, q) creates a JSP with activity durations de�ned based on the q value and conducts a branch-and-bound search where Monte Carlo simulation is used for each leaf node and standard constraint propagation is used at interior nodes. The �rst solution that is found whose probabilistic makespan is less than c is returned. �ndOptB&BSimLeaves(c, q) does the same as �ndFirstB&BSimLeaves(c, q) except the solution with lowest probabilistic makespan is returned rather than the �rst one. If no solution is found, a NIL value is returned. Unless the q value is low enough that the deterministic makespan is a lower bound on the probabilistic makespan, this function does not necessarily return the globally optimal solution. We �nd a starting solution with q = 0 to serve as an initial upper bound on the probabilistic makespan. In practice, B&B-DQ-L is run with a maximum CPU time. If q = 0 is reached within the time limit, this algorithm is approximately complete.

We refer to this algorithm as B&B-DQ-L as it does a series of Branch-and-Bound searches with Descending Q values and with simulation is used at the Leaves of the tree.

3.3 Branch-and-Bound Filtering Algorithms

Preliminary experiments with the branch-and-bound �ltering algorithm presented by Beck & Wilson showed that a signi�- cant amount of simulation was being done early in the search when the upper bound on the deterministic makespan was relatively high. Such solutions would seem to have little chance of being the optimal probabilistic solution as the deterministic solution is very poor. We present two constraint-based �ltering algorithms that seek to avoid these wasted simulations.

Branch-and-Bound Timed Better Solution A simple method is to spend a �xed amount of CPU time, tinitial to �nd a solution, s∗ , with a low deterministic makespan Dinitial using a �xed q value and constraint-based search. Then search can be restarted using the same q value and using Dinitial as an upper bound on the deterministic makespan. Whenever a solution, s, is found such that Dq(s) ≤ Dinitial , a simulation is run to evaluate the probabilistic makespan. The solution

B&B-TBS(q):

Returns the solution with lowest probabilistic makespan 1 (s∗ ,Dinitial) ←�ndOptB&B(∞, q, tinitial) 2 D∗ ←∞ 3 while termination criteria unmet AND not timed-out do 4 (s, D) ← �ndNextB&B(Dinitial, q, tremaining ) 5 D′ ← simulate(s) 6 if D′ < D∗ then 7 s∗ ← s; D∗ ← D′ end end 8 return s∗

with the lowest probabilistic makespan is returned when either the entire tree has been searched or when the maximum allowed CPU time has expired. Algorithm 2 contains pseudocode for this approach.

As above, we make use of a number of functions not de�ned with pseudocode: �ndOptB&B(c, q, t) creates a JSP with activity durations de�ned based on q and conducts a deterministic branch-and-bound search for up to t CPU seconds using c as an upper bound on the makespan. When the search tree is exhausted or the time-limit is reached, the best solution found, together with its makespan, are returned. No simulation is done. �ndNextB&B(c, q, t) produces a sequence of solutions (one solution each time it is called) whose deterministic makespan is less than or equal to c. The problem is de�ned using the q value and t is the CPU time limit. The solutions produced are the leaves of the B&B search tree in the order encountered by the algorithm. The c value does not change. Given enough CPU time, the algorithm will evaluate the probabilistic makespan of all solutions whose deterministic makespan is less than or equal to Dinitial. Finally, simulate(s) runs Monte Carlo simulation on solution s and the probabilistic makespan is returned.

The algorithm is called B&B-TBS for Branch-and-Bound Timed Better Solution. The algorithm is not complete, as there is no guarantee that the optimal probabilistic solution will have a deterministic makespan less than Dinitial .

Branch-and-Bound Iterative Best Solution A more extreme �ltering algorithm �rst �nds an optimal deterministic solution and uses it as a �lter for choosing the solutions to simulate. Using a �xed q value, a solution with an optimal deterministic makespan, D∗ q , is found and simulated. Then we execute a series of iterations beginning with i = 0. For each iteration, all solutions, se, with deterministic makespan Dq(se) ≤ D∗ q × (1 + i/100) are simulated and the one with the lowest probabilistic makespan is returned. If an optimal deterministic solution cannot be found within the CPU time limit, the best deterministic solution found is simulated and returned (i.e., only one simulation is done). Algorithm 3 presents the pseudocode.

The algorithm is theoretically complete. When i is large enough that the cost bound is greater than the deterministic makespan of all activity permutations, they will all be sim-

B&B-I-BS(q):

Returns the solution with smallest probabilistic makespan 1 (s∗ ,Dinitial) ←�ndOptB&B(∞, q, t0 -1) 2 D∗ ← simulate(s∗) 3 i ← 0 4 while not timed-out do 5 while termination criteria unmet do 6 (s, Dq ) ← �ndNextB&B( Dinitial × (1 + i/100), q, tremaining) 7 D ← simulate(s) 8 if D < D∗ then 9 s∗ ← s; D∗ ← D end end 10 i ← i + 1 end 11 return s∗

ulated. However, i would have to grow unreasonably large and therefore we treat this algorithm as, practically, incomplete. We refer to this algorithm as B&B-I-BS for Branchand-Bound-Iterative-Best Solution.

3.4 Local Search Filtering Algorithms

We also experiment with two local search �ltering algorithms analogous to the constraint-based algorithms above except that tabu search is used in place of constraint-based search.

Tabu Timed Better Solution For a �xed q and for a �xed CPU time, tinitial, a solution with the lowest deterministic makespan, Dinitial, possible is sought. Search is then restarted and whenever a solution, s, is found that has a deterministic makespan Dq(s) ≤ Dinitial, simulation is used to evaluate the probabilistic makespan. The solution with the lowest probabilistic makespan is returned. The pseudocode for Tabu-TBS is the same as that for B&B-TBS (Algorithm 2) with the following changes. The function �ndNextTabu(c, q, t) replaces �ndNextB&B(c, q, t); it produces a sequence of solutions whose deterministic makespan is less than c. The �ndBestTabu(c, q, t) function replaces �ndOptB&B(c, q, t); tabu search is run for up to t CPU seconds and the solution with the lowest deterministic makespan that is less than c is returned. We call this algorithmTabu-TBS for Tabu-Timed Better Solution.

Tabu Iterative Best Solution The core tabu search implementation does not necessarily use the entire CPU time. We can therefore create an iterative tabu-based solver for the probabilistic JSP similar to B&B-I-BS. In the �rst phase, using the overall time limit less one second, tabu search is used to �nd a good deterministic solution, based on a �xed q value. That solution is then simulated. Any remaining time is spent generating solutions with a deterministic makespan within a �xed percentage of the initial solution's deterministic makespan. As with B&B-I-BS, iterations are run with increasing i value starting with i = 0. In each iteration, solutions found by the tabu search whose deterministic makespan are within 1 + i/100 of the initial deterministic makespan are simulated. The solution with the lowest probabilistic makespan is returned. The algorithm is termed Tabu-I-BS for Tabu-Iterative-Best Search. The pseudocode for this algorithm is the same as that for B&B-I-BS (Algorithm 3) with the same substitutions that were made for Tabu-TBS.

3.5 Core AlgorithmDetails

The branch-and-bound algorithms described use texturebased heuristics [Beck and Fox, 2000] to decide on the the pair of activities to sequence and which sequence to try �rst. When constraint propagation is used (i.e., all B&B algorithms except B&B-N), we use temporal propagation, timetables [Le Pape et al., 1994], edge-�nder [Nuijten, 1994], and the balance constraint [Laborie, 2003].

The tabu search is the same one used by Beck & Wilson: the TSAB algorithm of [Nowicki and Smutnicki, 1996].

4 Empirical Investigations

Our empirical investigations aim to evaluate and investigate the performance of the algorithms. In particular, we are interested in algorithm scaling with problem size and uncertainty level and with the usefulness of representing uncertainty information using q values.

4.1 Experimental Details

We use four sets of probabilistic JSPs of size {4×4, 6×6, 10× 10, 20 × 20} with three uncertainty levels uj ∈ {0.1, 0.5, 1}. A deterministic problem is generated using an existing generator [Watson et al., 2002] with integer durations drawn uniformly from the interval [1, 99]. Three probabilistic instances are then produced by setting the mean durations µi to be the deterministic durations and by uniformly drawing the standard deviation σi from the interval [0, ujµi]. The distribution of each duration is approximately normal. For each problem size, we generate 10 deterministic problems which are transformed into 30 probabilistic instances.

Given the stochastic nature of the simulation and the tabu search algorithm, each algorithm is run 10 times on each problem instance with different random seeds. Each run has a time limit of 600 CPU seconds. Each Monte Carlo simulation uses 1000 trials. For B&B-TBS and Tabu-TBS, tinitial is 60 CPU seconds.

For the heuristic techniques, a deterministic duration is assigned to each of the activities based on the q value. We experiment with the same four q values used by Beck & Wilson as displayed in Table 1. The B&B-DQ-L algorithm requires an initial value of q, q0, and a decrement, qdec. For all except the largest problems, we used q0 = 1.25 and qdec = 0.05. Preliminary experiments with the 20×20 problems indicated that q0 = 1.25 resulted in long runs with no solutions due to a large search space with few if any solutions. To avoid such runs, we used q0 = 0.9, for which solutions could be found.

The probabilistic makespans are based on con�dence level p = 0.95. Our primary evaluation criterion is the mean normalized probabilistic makespan (MNPM) that each algorithm achieves: MNPM(a, L) = 1 |L| ∑l∈L D(a,l) Dlb(l) . L is a set of problem instances, D(a, l) is the mean probabilistic makespan found by algorithm a on l over 10 runs, Dlb(l) is the lower bound on the probabilistic makespan for l. For all problems except 20 × 20, the Dlb is found by solving the deterministic problems to optimality using the q1 durations (see Section 3 of [Beck and Wilson, 2004]). For the 20×20 problems, the optimal solutions could not be found and so the Dlb values represent the best deterministic solutions found. These values, therefore, are not a true lower bound.

The hardware used for the experiments is a 1.8GHz Pentium 4 with 512 M of main memory running Linux RedHat 9. All algorithms were implemented using ILOG Scheduler 5.3.

4.2 Results

Table 2 presents the results of our experiment. Each cell is the mean over 10 independent runs of 10 problems. The observed mean standard deviations for each cell are small: the maximum over all cells is less than 10.3% of the corresponding mean makespan, while the mean over all cells is less than 0.8%. We did not observe a large difference among the q values provided q > 0 and therefore present the results for q2.

For the 4 × 4 problems, a number of algorithms �nd lower mean probabilistic makespans than B&B-N despite the fact that B&B-N terminates before the limit on CPU time and therefore �nds approximately optimal solutions. This is an artifact resulting from algorithms that simulate the same solution multiple times. In B&B-N, a particular solution is only simulated once. In other algorithms, the same solution may be simulated multiple times leading to a bias to lower probabilistic makespan values. Based on the theoretic model of Beck & Wilson, these values still correspond to approximately optimal solutions.

4.3 Analysis

Overall the empirical results agree with those presented by Beck & Wilson. Here, we will focus on investigations of B&B-DQ-L and of the behavior of the heuristic algorithms.

Branch-and-Bound Descending-Q Leaves

B&B-DQ-L out-performs B&B-N across all problem sets even when it is not able to reach q = 0. As an indication of the number of iterations, the mean minimum values q for each problem size are: 4×4 : 0.007, 6×6 : 0.56, 10×10 : 1.00, 20×20 : 0.9.

For B&B-DQ-L, the deterministic durations de�ned by the q value serve to guide and prune the search for each iteration and, therefore, as with the heuristic algorithms (see below), the search is heuristically guided to the extent that solutions with low deterministic makespans also have low probabilistic makespans. However, the characteristics of the solutions found by the search are unclear. Ideally, we would like the search with high q to �nd solutions with very good probabilistic makespans both because we wish to �nd good solutions quickly and because the probabilistic makespan values are used to prune subsequent search. Analysis indicates that the �rst solution found by B&B-DQ-L is consistently better than the �rst solution found by B&B-N and for the larger problems (10×10 and 20×20) the relative increase in solution quality over subsequent iterations is greater for B&B-DQ-L.

Heuristic Algorithms

To investigate the performance of the heuristic algorithms, we look at the performance differences coming from using the q values, and the extent to which the performance of the heuristic algorithms depends on two factors: their ability to �nd good deterministic solutions and the underlying correlation between deterministic and probabilistic makespans in the experimental problems.

The Effect of the q Values Using a non-zero q value almost always results in better performance. This appears to be particularly true when either the uncertainty or the problem size increases. Table 3 shows the difference between mean normalized probabilistic makespan when using q0 and q2. On a few problem sets, negative values indicate that the runs with q0 deliver a better mean solution. However, by far the majority of the values are greater than 0 and their magnitude increases with both uncertainty and problem size, indicating increasingly good relative performance on the q2 problems. We show below that one explanation for this performance is that problem instances with q > 0 have a greater correlation between deterministic and probabilistic makespans.

Finding Good Deterministic Makespans We hypothesize that the performance of the heuristic techniques (and B&BDQ-L) is partially dependent upon the ability to �nd solutions with low deterministic makespans. We looked at two measures of the performance of B&B-I-BS and Tabu-I-BS: the quality of the best deterministic solutions found and the number of iterations. Because B&B-I-BS performs systematic search of all deterministic solutions with makespan below a given threshold, we hypothesize that on problem sets where it outperforms Tabu-I-BS, it will also have found better deterministic solutions and/or be able to systematically search solutions up to a higher threshold. In contrast, for problem sets where Tabu-I-BS outperforms B&B-I-BS, it will be due to �nding better deterministic solutions.

Table 4 presents data using q = q2. The mean normalized deterministic makespan (MNDM) is the calculated as follows: MNDM(a, L) = 1 |L| ∑l∈L Dq (a,l) Dq,min(l,B&B-I-BS) . The notation is as above with the addition of Dq,min(l, B&B-I-BS), the lowest deterministic makespan found by the B&B-I-BS algorithm over all runs on problem l. MNDM, therefore, provides a relative measure of the mean deterministic makespans from the two algorithms: the higher the value, the worse is the average makespan found relative to B&B-I-BS. The second measure (Iters) is the mean number of iterations performed by each algorithm.

Table 4 is consistent with our hypotheses. For those problems sets with large performance differences (i.e., 10 × 10 and 20×20) in probabilistic makespan, the better performing algorithm does �nd better deterministic makespans.

The Correlation Between Deterministic and Probabilistic Makespan The explanatory power of an algorithm's ability to �nd good deterministic makespans is, by itself, insuf�cient unless there is a correlation between the good deterministic and good probabilistic makespans. It is reasonable to expect that the level of uncertainty has an impact on any correlation: at low uncertainty, the variations in duration are small, so we can expect the probabilistic makespan to be relatively close to the deterministic makespan. When the uncertainty level is high, the distribution of probabilistic makespans for a single solution will be wider, resulting in a lower correlation. We hypothesize that this impact of uncertainty level contributes to the observed performance degradation (see Table 2) of the heuristic techniques with higher uncertainty levels.

We generated 100 new 10 × 10 deterministic JSP problem instances with the same generator and parameters used above. The standard deviation for the duration of each activity in the 100 instances were generated independently for each of �ve uncertainty levels uj ∈ {0.1, 0.5, 1, 2, 3} resulting in a total of 500 problem instances. For each instance and each of the four q values above, we randomly generated and simulated 100 deterministic solutions. Using R [R Development Core Team, 2004], we measured the correlation coef�cient, r, for each problem set.

Table 5 supports our hypothesis. As the uncertainty level increases, the correlation between the deterministic makespan and corresponding probabilistic makespan lessens. The strength of the correlation is somewhat surprising: even for the highest uncertainty level, r is more than 0.94 for q2 and q3. This suggests that the heuristic algorithms can be successfully applied to problem with higher uncertainty levels provided the appropriate q value can be found.

The correlation results in Table 5 also provide an explanation for the effect of the q values: they increase the correlation between deterministic and probabilistic makespans.

5 Conclusion

In this paper, we presented and conducted empirical analysis of a number of algorithms to proactively solve the job shop scheduling problem with probabilistic durations. Our empirical studies have demonstrated that a novel algorithm, B&B-DQ-L, based on iteratively reducing a parameter that determines the validity of the lower bound results in equal performance on small problems and much better performance on larger problems when compared to the previous branchand-bound technique. However, highest performance was observed in heuristic algorithms based on using the deterministic makespan to �lter the solutions which are to be simulated. It was argued, using detailed analysis of the experimental data and results from an ancillary experiment, that the performance of these algorithms is affected by their ability to �nd good deterministic makespans and the correlation in the problems between the quality of deterministic and probabilistic makespans. The correlation data on problems with higher uncertainty suggest that such algorithms may scale well to such problems. Central to the success of the heuristic algorithms was the use of a q value that governs the extent to which uncertainty (i.e., the variance) was represented in the durations of activities in deterministic problems. It was shown that such an incorporation of uncertainty data leads to a stronger correlation between deterministic and probabilistic makespans, and suggests a corresponding ability to �nd better probabilistic makespans.