Complete MCS-Based Search: Application to Resource Constrained Project Scheduling

1 Introduction

The resource constrained project scheduling problem (RCPSP) is one of the most general scheduling problem that is extensively studied in the literature. It consists in scheduling a project, which is a set of activities linked with precedence constraints, by means of a set of limited resources while minimizing the total duration of the project. The decision variant of the RCPSP, i.e., the problem of determining whether there exists a feasible project of makespan smaller than a given deadline, is NP-hard in the strong sense. The RCPSP is a very popular and frequently studied NP-hard optimization problem and the last 20 years have witnessed a tremendous improvement of both heuristic and exact solution procedures (cf. e.g. the recent surveys given in [Demeulemeester and Herroelen, 2002; Hartmann and Kolisch, 2005]). The currently best lower bounds on the makespan for the general RCPSP are based on solving linear programs using adequate cutting planes [Brucker and Knust, 2000; Baptiste and Demassey, 2004]. State-of-the-art techniques for upper-bounds rely on meta-heuristics such as Genetic Algorithms, Ant Colony Optimization or Large Neighborhood Search. Many scheduling problems such as job shop, cumulative job shop and open-shop can be modeled as special cases of RCPSP.

In this article, we present a pure constraint programming approach based on the exploration of a complete search tree to prove that the project cannot be achieved within a given deadline or to exhibit a feasible project if one exists. The search procedure is based on the detection and resolution of minimal critical sets (MCS) [Laborie and Ghallab, 1995] at each node of the search. MCSs are carefully chosen using a heuristic that tries to minimize the size of the search space. During the search, strong constraint propagation is enforced using classical scheduling constraint propagation techniques such as time-tabling, edge-finding, precedence energy and balance constraints [Laborie, 2003].

Next section recap the definition of the resource constrained project scheduling problem and introduces some notations. Section 3 describes our basic search procedure as well as the heuristic to select MCSs. Section 4 extends the basic search procedure to perform self-adapting shaving. The last part of the paper consists of experimental results on classical benchmarks (general RCPSP, open-shop and cumulative jobshop problems). For general RCPSP, we show that our approach closes more than 15% of previously open instances and improves more than 31% of best known lower bounds of the famous PSPLIB instances [Kolisch and Sprecher, 1996]. The same approach using exactely the same settings was used to close all the hard open-shop instances of [Gu´eret and Prins, 1999] in less than 5s CPU time and to improve the best known lower bounds and close several instances of cumulative jobshop [Nuijten, 1996].

2 Model and notations

The resource constrained project scheduling problem (RCPSP) can be formally stated as follows. A project is made of a set of activities A linked by precedence constraints. Precedence constraints can be represented by a directed acyclic graph G = (A, E) where each node in A represents an activity and each arc (A,B) ∈ E represents a precedence constraint between A and B. Let d(A) denote the fixed duration of activity A ∈ A and s(A) (resp. e(A)) denote the decision variable representing the start (resp. end) time of activity A. A set of discrete capacity resources R is considered, each resource R ∈ R having a maximal available capacity Q(R) over the entire scheduling horizon. Each activity A ∈ A requires a non-negative quantity q(A,R) of resource R. The problem is to find a feasible instantiation s of the activity start times such that precedence and resource constraints are satisfied and the schedule makespan is minimal. More formally:

where S(t) is the set of activities executing at time t:

A resource requirement of activity A on resource R is a triple u = (A,R, q) where q = q(A,R) > 0. If u = (A, R, q) is a resource requirement, we will denote A(u) = A the activity of u, R(u) = R the required resource, q(u) = q the required quantity, s(u) (resp. e(u)) will denote the start (resp. end) time of the activity of u. We will also denote U(R) = {u/R(u) = R} the set of resource requirements on resource R. If ϕ ⊆ U(R) is a subset of resource requirements on a resource R, we will denote q(ϕ) = ∑u∈ϕ q(u) the global resource consumption of R by activities in ϕ.

3 Search

3.1 Branching scheme

Our branching scheme assumes that a temporal network representing the relations between the time-points of all activities (start and end) using the point algebra of [Vilain and Kautz, 1986] is maintained during the search. In our implementation, this is ensured by the precedence graph constraint of ILOG SCHEDULER [ILOG, 2005]. We denote {∅,≺,⪯,=,≻,⪰, ̸=, ?} the set of qualitative relations between time points. If u and v represent two resource requirements, we will denote u ⪯ v if and only if e(u) ⪯ s(v).

The branching scheme relies on the notion of minimal critical sets (MCS) and their resolvers as introduced in [Laborie and Ghallab, 1995]. AMCS is a minimal set of resource requirements on the same resource R that could be executed simultaneously and would, in this case, over-consume resource R. MCSs are a natural generalization to cumulative scheduling of the pairs of activities conflicting for the same unary resource in disjunctive scheduling.

Definition 1 (Minimal critical set) Aminimal critical set on a resource R is a subset φ ⊆ U(R), such that:

1. Q(R) < q(φ)
2. ∀ϕ ⊊ φ, q(ϕ) ≤ Q(R)
3. ∧(u,v)∈φ×φ s(u) ≺ e(v) is consistent with the current temporal network

Informally, the different ways to resolve a minimal critical set consist in posting a precedence constraint between any two of its activities.

Definition 2 (Resolvers of a minimal critical set) If φ ⊆ U(R) is a MCS, we call resolvers of φ the set of temporal constraints Res(φ) = {u ⪯ v : (u, v) ∈ φ2 , u ̸= v}.

As described in [Laborie and Ghallab, 1995], the set of resolvers Res(φ) of a MCS φ can be simplified so as to remove those resolvers ρ ∈ Res(φ) for which there exists another resolver ρ′ ∈ Res(φ) such that ρ ⇒ ρ′ given the current temporal network. Indeed, in such case, the resolver ρ is redundant. Such a simplification procedure can be achieved in O(k3 ) if k is the size of the MCS using the naive Algorithm 1. Line 7 allows removing resolver u ⪯ v if there exists w such that u ⪯ v ⇒ u ⪯ w or u ⪯ v ⇒ w ⪯ v. In what follows, we assume that the set of resolvers of a MCS has been simplified.

Algorithm 1 Resolver simplification algorithm

1: procedure SIMPLIFY RESOLVERS(φ) 2: Res(φ) ← ∅ 3: for all u in φ do 4: for all v in φ \ {u} do 5: keep uv ← TRUE 6: for all w in φ \ {u, v} do 7: if s(v) ⪯ s(w) or e(w) ⪯ e(u) then 8: keep uv ← FALSE 9: break 10: if keep uv then 11: Res(φ) ← Res(φ) ∪ (u ⪯ v) 12: return Res(φ)

At each search node, our branching scheme consists in selecting a MCS φ and branching on its possible resolvers in the children nodes until there is no more MCS. This approach is clearly complete.

3.2 Heuristics

As all the resolvers consist of temporal constraints of the form x ⪯ y where x and y are two time-points (start or end of an activity), we are interested in an estimation of the size of the search space after posting such a precedence constraint. The fraction of the search space that is preserved when adding a precedence constraint is estimated using the complementary of the commitment measure introduced in [Laborie, 2003].

Let x and y be two time-points with respective lower and upper bound for time value: [xmin, xmax] and [ymin, ymax]. The size of the search space is estimated by the Cartesian product of the domain of the two variables, that is, the area of the rectangle [xmin, xmax], [ymin, ymax]. The size of the search space that is preserved when adding the constraint x ⪯ y is the part of that rectangle above the line x = y as illustrated in Figure 1. The fraction of the search space that is preserved can thus be estimated as follows. Let:

The fraction is then equal to:

In the example of Figure 1, this gives A = 30, B = 56, Cmin = 6, Cmax = 2, and preserved(x ⪯ y) = 24/30.

If ω is the size of the search space below the current search node, the size of the search space after posting a temporal constraint x ⪯ y can be estimated by ω · preserved(x ⪯ y). If φ is the MCS that is selected to be resolved at the current search node, the size of the search space to explore below the current node can thus be estimated as the sum of the sizes of the search space below each child node, that is: ω · ∑ρ∈Res(φ) preserved(ρ). Therefore, preserved(φ) = ∑ρ∈Res(φ) preserved(ρ) estimates the fraction of the search space that is preserved when choosing φ as the next MCS to solve 1 . Our heuristic simply chooses to resolve next the MCS φ that minimizes preserved(φ) that is, the one that minimizes our estimation of the size of the explored search space. Next section provides more details about the MCS selection algorithm.

3.3 MCS selection algorithm

At each search node, the MCS selection procedure develops a tree of partial MCSs where the current partial MCS φ is extended adding one resource requirement in each child node. By definition of preserved(φ), it is clear that φ′ ⊂ φ ⇒ preserved(φ′ ) ≤ preserved(φ). Thus, if φ∗ is the best MCS found so far, once a partial MCS φ has been reached such that preserved(φ∗ ) ≤ preserved(φ), the sub-tree of the MCS selection tree rooted at φ can be abandoned.

The algorithm for selecting and branching on a MCS is more precisely described in Algorithm 2 using the following notations: id(u) is a unique index associated with resource requirement u used to break ties. Unranked(u) represents all the resource requirements v that can possibly overlap u given the current temporal constraints, and that are such that q(v) < q(u) or q(v) = q(u) and id(v) < id(u). φ is the (partial) MCS that is currently being extended. A (partial) MCS is represented by a list of resource requirements: φ = (u1, ..., uk). We denote uk = Last(φ) and define the operator ⊕ as follows: φ⊕u = (u1, ..., uk, u); q = q(φ) is the consumption of the current (partial) MCS; p = preserved(φ)

1 Note that this is of course only a rough estimate and in particular, the estimated fraction preserved(φ) can be greater than 1.

is the preserved search space so far in the current (partial) MCS; φ∗ is the best MCS so far and p∗ = preserved(φ∗ ) is the preserved search space of the best MCS so far. The procedure at line 1 calls the MCS rating and selection process on each resource. At line 6, to rate and select MCSs for a given resource, the procedure first sorts the relevant sets of requirements v at lines 7 and 9 by decreasing order of q(v), using id(v) to break ties in order to ensure that each MCS is scanned only once, starting with the smallest MCSs, that is, the ones containing the most greedy requirements. The procedure at line 23 returns TRUE if and only if a given requirement u can possibly overlap all the requirements of a partial MCS given the current temporal network. The procedure at line 28 computes the incremental increase of preserved space due to the insertion of a new requirement u in the current MCS φ. The value of preserved has been described in section 3.2. The main recursive function for selecting MCSs is described at line 12.

Algorithm 2 MCS selection algorithm

1: procedure SELECT MCS 2: p∗ ←+∞ 3: for R in R do 4: SELECT MCS(R) 5: return φ∗ 6: procedure SELECT MCS(R) 7: Sort U(R) by decreasing q 8: for all u in U(R) do 9: Sort Unranked(u) by decreasing q 10: for u in U(R) do 11: RSELECT MCS(R,(u),q(u),0) 12: procedure RSELECT MCS(R,φ,q,p) 13: if q > Q(R) then ▷ φ is a MCS 14: if p < p∗ then ▷ φ is the best MCS so far 15: p∗ ← p 16: φ∗ ← φ 17: else ▷ φ needs to be extended 18: u ← Last(φ) 19: for all v in Unranked(u) do 20: if IS UNRANKED(v,φ) then 21: dp ← DELTA PRESERVED(v,φ) 22: RSELECT MCS(R,φ⊕v,q+q(v),p+dp) 23: procedure IS UNRANKED(u,φ) 24: for all v in φ do 25: if u ⪯ v or v ⪯ u then 26: return FALSE 27: return TRUE 28: procedure DELTA PRESERVED(u,φ) 29: dp←0 30: for all v in φ do 31: dp ← dp + preserved(u,v) + preserved(v,u) 32: return dp

cedure is selected as the one to be solved at the current search node. This MCS is simplified using Algorithm 1 and the search explores all of its resolvers ρ ∈ Res(φ∗ ) in the child nodes by decreasing order of preserved(ρ). This order has no effect when the schedule is not feasible as in this case the complete search tree needs to be explored but it helps finding a solution quicker when a solution exists.

4 Self-adapting shaving

Shaving techniques [Torres and Lopez, 2000] provide a good framework for strengthening constraint propagation and avoiding late failures to be discovered in the search tree. They are all based on the following principle: if adding a constraint C in the current node of the search leads to a failure of the propagation, then, constraint ¬C can be inferred. Due to the cost of propagating a constraint C and the potential number of constraints C to try to shave on, shaving techniques are in general computationally expensive.

To improve the pruning of the search tree, we implemented the following shaving technique based on MCSs. If a MCS φ with resolvers Res(φ) = {ρ1, ..., ρk, ρk+1} is such that ∀i ∈ [1..k], adding ρi in the current schedule leads to a failure of the propagation, then ρk+1 can be inferred. The complexity for shaving a given MCS φ of size n is thus in O(n2P) where P is the cost of full constraint propagation at the current node. Potentially, there is of course an exponential number of MCSs to shave on at each search node and we can expect that many of those MCS do not allow inferring any precedence constraint. An idea to speed-up the shaving process is thus to only try shaving on a subset of MCSs for which the probability to infer a precedence constraint is greater than a given threshold α. Parameter α is an input of the shaving algorithm. We can roughly estimate that the probability that adding a precedence constraint x ⪯ y in the current schedule will lead to a failure of the propagation is proportional 2 to 1 - preserved(x ⪯ y). If ρm = argmaxρ∈Res(φ) preserved(ρ) is the resolver of the MCS φ with maximal preserved search space, we are interested in the MCSs that get a high probability that all their resolvers but ρm will fail, that is, if we assume all the probabilities are independent, the ones such that Πρ∈Res(φ)\{ρm}(1preserved(ρ)) is greater than a given threshold. For those MCSs, if the threshold is close enough to 1, we can assume that preserved(ρ) is small enough so that the first order approximation Πρ∈Res(φ)\{ρm}(1 - preserved(ρ)) ≈ 1 - ∑ρ∈Res(φ)\{ρm} preserved(ρ) is reasonable.

To summarize, we thus only consider for shaving those MCS scanned by the procedure described in Algorithm 2 that are such that preserved(φ) - preserved(ρm) ≤ β, β being a threshold. The computation of this criterion only adds a negligible overhead related with the maintenance of ρm for each MCS in the MCS selection procedure. Due to the numerous approximations, β is not taken to be constant (theoretically equal to 1 - α). The threshold β is computed by self-adaptation in such a way that on average, among the last h shaving attempts, αh lead to the inference of a new precedence. Whenever the number of successful shaving s among the last h ones deviates from αh, the parameter β is adapted accordingly: if s < αh, β is decreased by ϵβ and if s > αh, β is increased by ϵβ. For all our experiments with shaving, we took h = 20, α = 0.75, ϵ = 0.01 and start with β = 1.

2 Note that this estimation is exact at the extreme points when preserved(x ⪯ y) = 0 (propagation will fail for sure) and when preserved(x ⪯ y) = 1 (propagation cannot fail because x ⪯ y has already been discovered given the current domains of x and y).

5 Experimental evaluation

The approach has been implemented on top of ILOG SCHEDULER 6.1 using the timetable, disjunctive, edgefinder, precedence energy and balance constraints [ILOG, 2005]. All the experiments described in this section were run on a Dell Latitude D600 laptop, 1.4 GHz. Detailed results such as individual lower bounds for each problem instance and new optimal solutions are available at http://scheduler.ilog.fr/.

5.1 Results on general RCPSP

We evaluated our approach on the instances of the PSPLIB [Kolisch and Sprecher, 1996] with 60, 90 and 120 activities (resp. sets J60, J90, J120). For each instance, we solve the feasibility problem of finding a schedule with a makespan lower than T, starting with a legal lower bound for T 3 and incrementing T until either the problem is shown to be feasible (in this case, T is the optimal makespan) or a given time limit for solving the problem with makespan T is exceeded (in this case, T is a legal lower-bound but the search stops without providing any legal upper-bound).

In a first series of experiments, we use the basic search described in section 3 without shaving with a time limit of 300s. The previous best lower and upper bounds we compare with are the ones reported in the PSPLIB together with the recent improvements on the J60 instances reported in [Baptiste and Demassey, 2004]. The results are summarized on Table 1 with the following columns:

#O : number of instances previously open

#I : number of improved lower bounds (% of #O)

AGR: average gap (distance from the lower to the upper bound) reduction when a bound is improved

#C : number of closed instances (% of #O)

Out of the 617 previously open instances, we are able to improve 178 lower-bounds with an average gap reduction of 53% and to close 85 instances.

To show the effect of self-adapting shaving, we run a version of our approach using self-adapting shaving with the same time-limit of 300s. The results are summarized on Table 2. Out of the 617 previously open instances, we are able to improve 182 lower-bounds with an average gap reduction of 53.1% and to close 86 instances. The main conclusion is that, within the same time-limit, the addition of self-adapting shaving slightly increases the performances. The two curves below respectively show, on a particular instance (J60 5 2), the number of search nodes with a given depth in the search tree and, for each node depth, the ratio between the number of selected MCSs that effectively lead to the inference of a new precedence and the total number of MCSs selected for shaving. One clearly see that in the region of the search space where most of the nodes are concentrated (between depths 10 and 35), this ratio is effectively close to the target ratio of 0.75. In this instance, 4667 nodes where explored, 4843 MCSs where selected for shaving and among them, 3291 effectively lead to the inference of a new precedence.

3 For instance the lower-bound of the PERT of temporal constraints

Given that self-adapting shaving slightly improves the performances within the same time limit of 300s, we used this configuration with an extended time-limit of 1800s. The results are summarized on Table 3. Out of the 617 previously open instances, we improve 193 lower-bounds (that is more than 31% of the previously open instances) with an average gap reduction of 61.1% and close 96 instances (that is more than 15% of the previously open instances).

5.2 Results on open-shop problems

Open-shop problems can be represented as a special cases of RCPSP where all resources have a unit capacity and additional unary resources are used to model the fact that activities of the same job cannot not overlap. We tested our approach, with the same settings as in previous section on the open-shop problems instances proposed in [Gu´eret and Prins, 1999]. Those instances are considered to be very hard instances of open-shop problems and serve as classical benchmark in open-shop scheduling (see for instance recent work in [Blum, 2005]). The benchmark consists of 80 instances ranging from 3 jobs × 3 machines problems until 10 jobs × 10 machines problems. Out of these 80 problems, 34 instances are still open. Using our approach, we were able to close all those instances in less than 5s CPU time. The optimal makespan for the 34 previously open instances is summarized on Table 4 where the column UB corresponds to the currently best known upper-bound for which no optimality proof did exist.

We also experimented with the open instances of the benchmark of [Brucker et al., 1997]. We closed 3 of these 6 open instances: namely j8-per0-2 (optimal makespan: 1052), j8-per10-0 (optimal makespan 1017) and j8-per10-1 (optimal makespan 1000).

5.3 Results on cumulative jobshop problems

We tested our approach, with the same settings, on the cumulative job-shop problem benchmark described in [Nuijten, 1996]. These instances are derived from classical jobshop scheduling problems by multiplying the number of jobs (and thus the number of activities) and the capacity of the resources by a given factor (×2 or ×3). Our results are summarized on Table 5 where LB is the lower bound using the consistency checking described in [Nuijten, 1996] and NewLB the new lower bound of our approach. We were able to close the ft06×2 and ft06×3 instances as well as to improve 12 lower bounds out of these 38 open instances.

6 Conclusions

We presented a simple complete search procedure implemented on top of classical constraint propagation algorithms and applied it to resource constrained project scheduling problems. In average, this approach outperforms the best algorithms for finding lower bounds on those scheduling problems, even with a time limit of 300s per optimization step 4 . Using this approach in conjunction with a self-adapting shaving procedure, we were able to close more than 15% of the previously open problems of the PSPLIB and improve more than 31% of the best known lower bounds. What is even more remarkable is that this very same approach allows closing all the hard open-shop instances of [Gu´eret and Prins, 1999] in less than 5s CPU time although the approach was not particularly designed to tackle disjunctive scheduling and does not exploit the open-shop nature of the problems.

The understanding of why our method works so well on the instances of the PSPLIB and on many open-shop problems would deserve a deeper study. From one hand, if the problems are highly cumulative, our approach is clearly limited by the explosion of the number of MCSs to consider. From the other hand, when problems are highly disjunctive, we could expect other approaches dedicated to disjunctive scheduling to work better. A first possible explanation could be a good fit between our approach and the "disjunctivity" degree of the hard instances of the PSPLIB as suggested by some recent work [Garaix et al., 2005]. A result of this study could be some kind of hybridizing of MCS-based search with techniques more adapted to highly cumulative problems, MCSbased search being restricted to the resolution of MCSs with small preserved search space (thus small MCSs) at the top of the search tree. A second direction for future work is the generalization of the notion of self-adapting shaving introduced in this paper to other shaving techniques in scheduling.

4 When this time limit is exceeded at an optimization step, usually, the previous steps where fairly quick so that the overall time for computing the lower bound is close to 300s.