European Journal of Operational Research 245 (2015) 70–80

Contents lists available at ScienceDirect

European Journal of Operational Research
journal homepage: www.elsevier.com/locate/ejor

Discrete Optimization

A branch-price-and-cut algorithm for multi-mode resource leveling ✩Eamonn T. Coughlan a, Marco E. Lübbecke b, Jens Schulz a,∗abInstitut für Mathematik, Technische Universität Berlin, MA 5-1, Straße d. 17. Juni 136, 10623 Berlin, Germany
Operations Research, RWTH Aachen University, Kackertstraße 7, 52072 Aachen, Germanya r t i c l e

i n f oArticle history:
Received 4 October 2013
Accepted 22 February 2015
Available online 5 March 2015
Keywords:
Project scheduling
Resource leveling
Dantzig–Wolfe decomposition
Integer programming
Branch-price-and-cut

a b s t r a c t
We propose a new mixed integer programming formulation and solution algorithm for a multi-mode
resource-constrained project scheduling problem with availability constraints (calendars) and the objective
to minimize the resource availability cost. Our model exploits the problem structure and has an exponential
number of variables, which necessitates a column generation approach. The linear programming relaxationis strengthened by adding valid inequalities that need to be carefully separated in order to show the desired
effect. Integer optimal solutions are obtained by an exact state-of-the-art branch-price-and-cut algorithm.
Classical time-indexed mixed integer programming formulations for similar problems quite often fail alreadyon instances with only 30 jobs, depending on the network complexity and the total freedom of arranging
jobs. A reason is the typically very weak linear programming relaxation. Our model can be interpreted asa non-trivial Dantzig–Wolfe reformulation of a time-indexed formulation. In particular, for larger instances
our reformulation gives signiﬁcantly tighter dual bounds, enabling us to optimally solve instances with 50
multi-mode jobs. This outperforms the classical formulation by far.
© 2015 Elsevier B.V. All rights reserved.

1. Introduction
We consider a complex project scheduling problem that is motivated by an industrial application from chemical engineering. It
was recently introduced in Megow, Möhring, and Schulz (2011) as
turnaround scheduling: For the inspection and renewal of parts inchemical plants, these are shut down and disassembled, work is done,
and plants are ﬁnally rebuilt. This induces a partial ordering of jobs to
be done. Jobs are multi-mode, that is, they can be sped-up to a certain
extent by investing in more workers. The time horizon and the number of workers hired for each job determine production downtime
and working cost. We thus face a time-cost tradeoff, and a solution
approach to turnaround scheduling is to binary search over the time
horizon. For a ﬁxed time horizon—which is assumed throughout thispaper—the problem turns into a resource leveling problem, i.e., we✩
A preliminary version appeared in Coughlan, Lübbecke, and Schulz (2010). Our newcomputational results reveal that the tighter dual bounds are a crucial ingredient forsolving the proposed instances. With new cutting planes and a customized branchingstrategy, we increase the critical size of hard instances from before 30 to 50 jobs.
∗
Corresponding author. BASF SE, Scientiﬁc Computing, Carl-Bosch-Straße 31, 67056
Ludwigshafen, Germany. Tel.: +49 621 60 94996.
E-mail addresses: coughlan@math.tu-berlin.de (E.T. Coughlan), marco.luebbecke@rwth-aachen.de (M.E. Lübbecke), jens.schulz@basf.com, jschulz@math.tu-berlin.de(J. Schulz).

http://dx.doi.org/10.1016/j.ejor.2015.02.043
0377-2217/© 2015 Elsevier B.V. All rights reserved.

need to decide how many workers are hired—this number is thenﬁxed for the entire time horizon.
Workers, or more generally renewable resources, are of different
specialization or type. Each resource type is associated with availability periods that can be thought of as working shifts. We assume thateach job needs (possibly several units of) exactly one resource type
for the following reason. The processing time of a job decreases withan increased resource usage. A model with multiple resource types
per job must encode how the processing time of a job is shortened depending on which resource types are increased. Planners would needto encode all possible combinations of resource assignments to a job,and even worse, estimate the resulting processing times, possibly respecting synergy or interference between resources. It is unrealistic
to assume valid data about this in a practical setting on one hand, andhardly doable, if at all, with standard software like MS ProjectTM on
the other hand. Reliable data are diﬃcult to obtain even for one resource type (the actual work increases with increased resource usage
due to, e.g., communication overhead). Nonetheless, the need of multiple resources was studied by other authors, e.g., Santos and Tereso(2011) and can be modeled via generalized precedence constraints.
These blend easily with our approach, if absolutely needed.
We want to level the resources, i.e., we minimize the maximum
capacity requirement of each resource type, at a minimum total resource cost. Balancing the resource usage, i.e., to ﬂatten the variabilityof the resource usage over time, is not an issue at our higher-level
planning stage.

E.T. Coughlan et al. / European Journal of Operational Research 245 (2015) 70–80

From a methodological point of view, when it comes to optimallysolving larger instances of such scheduling problems, mixed integer
programming (MIP) has not fared very well. This is mainly due tothe weak dual bounds obtained from linear programming (LP) relaxations of classical models. We therefore propose a model which
can be interpreted as a non-trivial Dantzig–Wolfe reformulation of atime-indexed formulation and thus provides much stronger bounds.
Variables represent entire sub-schedules, thus are exponential in
number, and necessitate a dynamic generation of variables to solve
the LP relaxation. Yet, the relation to classical models via Dantzig–
Wolfe reformulation enables us to exploit the literature on valid inequalities for scheduling problems, and we demonstrate how to usethem to strengthen the formulation. Constraint programming (CP) is
usually considered a promising approach for solving scheduling problems, due its expressiveness of logical constraints, and we hybridize
some CP elements in our otherwise MIP based approach.
1.1. Our Contribution
We propose a new model for a multi-mode project schedulingproblem with availability constraints (working shifts) and the objective to minimize the resource availability cost. Our model is genericin the way resource calendars are respected. For its solution we design a full-ﬂedged branch-price-and-cut algorithm. This does not only
produce very strong dual bounds, as did a few previous studies inproject scheduling before, but it is also able to provide optimal primal
solutions. With this it is the ﬁrst approach in this area of complex
project scheduling that “works on both ends.” Our model builds on a
non-trivial Dantzig–Wolfe reformulation of a classical time-indexed
model, and we exploit the relationship between these two in designing branching rules and deriving cutting planes. We ﬁnally conduct
an experimental study which shows how the methodological stateof-the-art in branch-and-price should be instantiated in a complex
scheduling context. Our implementation is able to optimally solve
instances with 50 multi-mode jobs, thus outperforming a state-ofthe-art MIP solver by far that is only able to solve 25 percent of theinstances we solve.
2. Formal problem description
We assume familiarity with resource-constrained project scheduling (RCPSP) in general; for a recent survey refer to Hartmann andBriskorn (2010). We are given a set J of non-preemptable jobs anda set R of renewable resources. Precedence constraints are given as
an acyclic digraph G = (J , E) with (i, j) ∈ E if and only if job i has tobe completed before job j starts. Each job j may be run in exactly oneout of a discrete set Mj of modes. Processing job j in mode m ∈ Mj
takes pjm time units and requires rjmk units of resource k ∈ R. Each71

job needs exactly one resource for execution, so we write rjm if theresource is clear from the context.
All jobs have to be completed before the end of the time horizon T.
Each resource k ∈ R has a set Ik := {[a1 , b1 ], . . . , [aik , bik ]} of ik availability periods, also called shifts, where a1 < b1 < · · · < aik < bik . A job
requiring resource k can only be executed during a time interval I ∈ Ik .
Fig. 1 schematically presents a schedule with two resources and threeavailability periods each that are indicated by the bold lines. Jobs are
represented by boxes, the length of which represents processing time
and the height corresponds to resource usage (for simplicity all identical in this example). The arrows indicate precedence constraints. We
use a parameter δkt which is one if resource k is available at time t,i.e., t ∈ I for some I ∈ Ik , and zero otherwise. Each resource k ∈ R is
associated with a per unit cost ck . For each resource k, the available
capacity, which is a variable, is denoted by Rk .
To represent a solution, we denote by S = (S1 , . . . , Sn ) the vectorof start times of jobs, and by M = (m1 , . . . , mn ) the vector of modes
in which jobs are executed. For a given schedule (S, M), denote byA(S, M, t) := {j ∈ J : Sj ≤ t < Sj + pjmj } the set of jobs active at time t.

The amount rk (S, M, t) := j∈A(S,M,t) rjmj k of resource k used at time t
must never exceed the provided capacity Rk . Thus, we obtain resource
constraints with calendars: rk (S, M, t) ≤ Rk · δkt , ∀k ∈ R, ∀t. In addition to this resource feasibility a feasible schedule must obey precedence feasibility, i.e., Si + pimi ≤ Sj for all (i, j) ∈ E.
To conclude, we study a multi-mode project scheduling with mrenewable resources of unbounded capacity with precedence constraints and working shifts and the objective to minimize the total
resource availability cost, i.e., minimizing k∈R ck · Rk .
2.1. Related work
Turnaround scheduling comprises project scheduling with calendars, multi-mode scheduling, and resource leveling; see again Megow
et al. (2011) for an industrial application. The zoo of project scheduling problems is large, and we mention only the most related problems. Makespan minimization is a classical scheduling goal. Lowerbounding schemes for this objective are presented in Brucker andKnust (2000), where column generation is employed to solve a relaxed problem, allowing preemption and precedence constraints formulated as disjunctions. A variable represents a set of jobs selected torun at a certain point in time. For the case of generalized precedence
constraints, lower bounds are derived in Bianco and Caramia (2011)by relaxing resource constraints for jobs which are not precedence
related. This allows a dynamic programming approach on a modiﬁed
activity-on-nodes network. In contrast to minimizing the makespan,other objective functions that measure the variation of resource utilization, e.g., f (rk (S, t)) are of interest in the pre-planning phase; seee.g., Neumann, Schwindt, and Zimmermann (2003).

Fig. 1. Schematic representation of turnaround scheduling with two resources. Arrows indicate precedence constraints between jobs; bold lines represent availability periods.

72

E.T. Coughlan et al. / European Journal of Operational Research 245 (2015) 70–80

The resource leveling problem with single-modes per job, which
is denoted by PS|temp| ck max rk (S, t) with generalized precedence
constraints has been considered earlier under the name resource
investment problem. The special case without generalized prece
dence constraints PS|prec| ck max rk (S, t) has been considered e.g.,
in Demeulemeester (1995) and Möhring (1984). These papers competed on the same instance set which contained about 16 jobs and four
resources, with a time horizon between 47 and 70 time units. Further
computational studies were done containing 15 to 20 jobs and four resources. In the same setting, lower bound computations are proposed
in Drexl and Kimms (2001), one based on Lagrangian relaxation, andone based on column generation, where variables represent schedules
as in our approach. Twenty jobs with small processing times of about
10 time units can be handled; for 30 jobs the Lagrangian relaxation
wins against the column generation approach.
Multi-mode jobs are a key feature of turnaround scheduling.
Such problems of the form MPS|prec|Cmax have been investigatedwith renewable and non-renewable resources, with limited capacity, and makespan minimization, known as multi-mode RCPSP, seee.g., Demeulemeester and Herroelen (2002) and Hartmann (2001).

The notation MPSm, ∞|prec, shifts| ck · max rk (S, M, t) comes close
to our model, but in the literature in an MPS problem, a job requires
different types of resources.
Calendars have been taken into account in previous algorithms as
well. Scheduling problems with ﬁxed processing times and calendars,
but without resource capacities, are considered in Zhan (1992). The
author provides an exact pseudo-polynomial time algorithm (turned
into a polynomial one in Franck, Neumann, and Schwindt (2001)) forcomputing earliest and latest start times for preemptable as well as
non-preemptable jobs.
For a computational benchmarking of project scheduling problems, different problem sets are available in the PSPLib (Kolisch &Sprecher, 1996), where several variants of the RCPSP and of resource
investment problems can be found. For the RCPSP single-mode case,test sets containing 60 jobs could not be solved in total by a vast
number of researchers. In the multi-mode case, instances with 30
jobs are not solved yet. For the resource investment problem, test
sets containing 10, 20, or 30 jobs are available, but they do not contain working shifts, are in single-mode or include generalized precedence constraints. On the other hand a job may need more than one
resource type. Even though none of these problem variants is immediately suited for a direct comparison, they are similar to ours, and thementioned instances inspired us when generating the test set used inthis study (see Coughlan et al., 2010 and Section 5).
3. Mixed integer programming formulations
For solving large-scale scheduling problems, mixed integer programming is not considered as primary choice because of typicallyweak LP relaxations. Moreover, huge numbers of variables (in particular for time-indexed formulations) and constraints may resultin high computation times and memory failures when solving even
only the LP relaxation. The approach we propose demonstrates that
more sophisticated algorithmic techniques can be a partial remedy
for these issues. In the following, we assume the reader to be familiar
with solving MIPs by branch-and-bound, see e.g., Achterberg (2009),
and the basics of column generation and branch-and-price, see e.g.,
Desrosiers and Lübbecke (2005).
3.1. Shortcomings of time-indexed mixed integer programs for RCPSPOne of the most prominent models for the RCPSP is based on a time
discretization and was introduced in Pritsker, Watters, and Wolfe(1969). Their formulation adapted to our resource leveling problem
reads as follows:min

ck · Rk

(1)

k

s.t.

 
t

 
t

t · xjmt = Sj∀j∈J(t + pim )ximt ≤

 

m∈Mi

 

(2)

∀j∈J

m∈Mj

 
t

xjmt = 1

m∈Mj

t
t


(3)

txjmt

∀ (i, j) ∈ E

(4)

m∈Mj

rjmk · xjmτ ≤ Rk

∀ k, t(5)j∈J m∈Mj τ =t−pjm +1
τ ≥0

Sj ≥ 0, ∀ j ∈ Jxjmt ∈ {0, 1}

Rk ≥ 0,∀k∀ t, j, m ∈ Mj

(6)
(7)

Binary (“start-time”) variables xjmt model whether job j starts attime t in mode m or not. Each job j must start exactly once, and inexactly one mode (2). The start time Sj of job j ∈ J is linked to the binary variables xjmt in Equation (3) and is therefore implicitly integral.
This linking is redundant, but we keep it for instructional purposes.
Also precedence constraints (4), and resource capacity constraints (5)can be expressed as linear inequalities. In this model we decide about
the resource capacities Rk ≥ 0 for each resource k, such that the total
resource availability cost is minimized via the objective function (1).
For notational convenience we do not restrict the (j, m, t) indices of
xjmt variables to their feasible combinations. In fact, certain start times
and certain modes would violate the calendar boundaries. However,this can be eﬃciently taken care of in a preprocessing step, yielding ahuge reduction in the number of variables.
Depending on several factors, such as network complexity (thedensity of G) or the time discretization considered, this formulation
may give very weak lower bounds. Very often, we experience that inan optimal solution to the LP relaxation only few binary variables are
fractional, but the points in time used in the convex combination (3) toyield the actual start time Sj of a job are far apart from one another. A
simple example is xjmt1 = xjmt2 = 0.5 which results in Sj = (t1 + t2 )/2.
When t1 and t2 lie in different shifts it may happen that Sj is outside any shift, resulting in an infeasible start time. Fractions of jobs
may violate precedence constraints even though the aggregated constraint (4) is fulﬁlled for the fractional solution. This holds similarly
for the resource constraint (5). Things get more involved in our setting
as different modes can be fractionally used for each job, thus all job
durations between smallest and largest may appear. We informallycall this a “smearing” of start time variables. This smearing gives us
irrelevant information about the schedule and we lose most if not all
structure in the model. This happens often in MIP models, however,it hits us particularly hard for scheduling problems.
Finally, branching a binary variable xjmt to zero is a very weak
decision if the corresponding job can be scheduled one time unit
earlier or later. That is, the decision essentially has no effect. Not
so the upward branch, which ﬁxes a start time and mode and thus
imposes signiﬁcant structure to the overall solution. It is therefore
to be expected that branching on xjmt variables gives an unbalanced
search tree, which is not desired.
3.2. Master problem: A model based on shift conﬁgurations
In order to reduce the effects of “losing the timing information”
caused by the smearing of variables, we propose a model which exploits the problem structure by decomposing the time horizon into
the availability periods of resources. Based on the calendar for eachresource type, every working shift represents a smaller subproblem
for which sub-schedules, or conﬁgurations, are stated for one resource

E.T. Coughlan et al. / European Journal of Operational Research 245 (2015) 70–80

type only. These conﬁgurations are linked by global constraints ensuring that exactly one conﬁguration is chosen for each working shift,
and that precedence constraints are respected.
Deﬁnition 1. A conﬁguration ξ for an interval I ∈ Ik of resource type k
is a schedule for a subset J ⊆ J of jobs that all require resource
type k. Each job j ∈ J is speciﬁed by its start time Sjξ and its completion time Cjξ (this uniquely determines its mode). The maximum
simultaneous usage of resource type k by these jobs is denoted by r̄ξ .
For each conﬁguration ξ for a particular shift I ∈ Ik we introduce abinary variable xξ which indicates whether conﬁguration ξ is chosen
or not. We use the short hand notation j ∈ ξ to express that job j is
executed in the shift corresponding to conﬁguration ξ . Note that themode of each job is determined by the respective start and completion
times. The so-called master problem reads:min

ck · Rk

(8)

k

s.t.
Sj =

∀ (i, j) ∈ E(9)Sj ξ · x ξ

∀j∈J

(10)

C jξ · xξ

∀j∈J

(11)

∀k

(12)

C i ≤ Sj


Cj =




ξ :j∈ξ

r̄ξ · xξ ≤ Rk

∀ I ∈ Ik

ξ :ξ ∈Ixξ = 1∀j∈J(13)∀ k, j ∈ J(14)ξ :j∈ξ

Rk , Sj , Cj ≥ 0xξ ∈ {0, 1}∀ξcertain jobs infeasible. The objective function (16) reﬂects minimizing
the reduced cost.



max

(15)

Each job is executed exactly once by Equation (13). This constraint
also ensures, that for each job exactly one mode (from exactly one
conﬁguration) can be chosen. The start and completion times for eachjob are computed from the chosen conﬁgurations via the linking constraints (10) and (11). Constraints (9) model the precedence relations
between jobs. These could be directly expressed by substituting Sj
and Cj from the linking constraints, but (10) and (11) are helpful inthe pricing problem (Section 3.3) where they help penalizing or encouraging certain start or completion times of jobs. Constraints (12)link resource consumptions to the capacities and ensure that the total
resource capacity available over the planning horizon is at least the
maximum resource usage needed in each shift.
3.3. Column generation: Pricing problem
Since the number of feasible conﬁgurations is exponential in thenumber of jobs, we solve the LP relaxation of (8)–(15) by column
generation. That is, we start with a very small (e.g., heuristically generated) subset of conﬁguration variables, and dynamically add morevariables to the model until it can be proven that no more promisingvariables exist. This optimality proof is provided—as in the standard
simplex method—via non-negativity of reduced costs of all conﬁguration variables. It is checked with the help of an auxiliary optimization problem, the pricing subproblem, which determines a negative
reduced cost conﬁguration variable, or proves that none exist.
We must solve a pricing problem for every shift I ∈ Ik of eachresource type k. For reduced cost computations, we need the dual
variables of constraints (10)–(13) which are denoted by sj , cj , ρkI , andπj , respectively. For a ﬁxed shift I only the subset J ⊆ J of jobs needs
to be considered that can be scheduled in I. This subset is restricted
to jobs using the proper resource type, but also—as we see later—e.g., branching and other decisions may render certain start times of

π j · Xj −



j

s.t.

 

cj · Cj −jsj · Sj − ρkI · r̄k

(16)

j

xjmt = Xj

∀j∈J

(17)

t∈I m∈Mj t · xjmt = Sj

∀j∈J

(18)

t∈I m∈Mj (t + pjm ) · xjmt = Cj ∀ j ∈ J(19)t∈I m∈Mj t

rjmk · xjmτ ≤ r̄k

∀t ∈I

(20)

j∈J m∈Mj τ =t−pjm +1τ ∈Ir̄k ≥ 0Sj , C j ≥ 0

(21)

∀j∈J

xjmt ∈ {0, 1}
Xj ∈ {0, 1}ξ :j∈ξ

73∀ j ∈ J, m ∈ Mj , t ∈ I
∀j∈J

(22)
(23)
(24)

The Xj , Sj , and Cj variables are “reporting” variables, and again
present for instructional purposes only (as they could be substituted).
Binary variable Xj decides whether job j ∈ J is processed, and if so, at
what start time St and completion time Cj . Note once more, that thisuniquely determines the job’s mode. All three variables for each job
j ∈ J characterize a conﬁguration ξ for the speciﬁc shift I ∈ Ik of resource type k. This pricing problem is again a scheduling problem
with a non-regular objective function. By its mechanism, this is again
a time-indexed formulation, with binary start-time variables xjmt , as
is MIP (1)–(7). However, it has a much smaller time horizon (onlyone shift), and it is not required that each possible job j ∈ J is actually
executed. Constraints (17) model the selection of job j and its mode;
Constraints (18) and (19) encode the start and completion times ofjobs according to the chosen mode assignment. Resource capacity
constraints (20) have to be satisﬁed as before and the resource usage
in this shift is again a decision variable, denoted by r̄k . The objectivefunction value is rewarded by πj if job j is selected into the conﬁguration and by multiples of sj and cj according to the start and completion
times. A unit increase of resource usage penalizes the objective valueby ρkI .
The pricing problem (16)–(24) is NP-hard as it contains a resource
leveling problem as special case. This can be seen as follows: Set all πj
to a value large enough to ensure that each job must be scheduled,and let sj and cj be zero for all j. Then, we need to schedule all jobs
during I at a minimum resource cost ρkI · r̄k . As we are not aware ofother exact approaches to this kind of hard scheduling problem, we
solve the pricing problem as a mixed integer program itself.
The master problem (8)–(15) and the pricing problem (16)–(24)are solved alternately, exchanging dual variable values in one, andconﬁguration variables in the other direction, until optimality of themaster problem is proven. Note that an optimal solution to the masterproblem can be fractional, and integrality of conﬁguration variables xξ
still needs to be ensured. This leads us to a branch-and-price algorithm,
in which the LP relaxation in each node of the branch-and-bound tree
is solved by column generation as above.
3.4. Interpretation as non-trivial Dantzig–Wolfe Reformulation
In a Dantzig–Wolfe reformulation of a MIP min{cT x | Ax ≥ b, Dx ≥d, x ∈ Qn1 × Zn2 }, which is called the original formulation in thiscontext, we change the polyhedral representation for a subset of
constraints, say Dx ≥ d. The convex hull of all the integer points(resp. mixed integer sets) that are feasible for Dx ≥ d is expressed74

E.T. Coughlan et al. / European Journal of Operational Research 245 (2015) 70–80
as convex combination x = q λq xq , q λq = 1, λq ≥ 0 of its ﬁnitely
many extreme points xq (assuming that Dx ≥ d is bounded). This convex combination is substituted in the remaining constraints Ax ≥ b,yielding the master problem. The extreme points are given implicitly as solutions to an optimization problem, the pricing problem.
The reformulation, which results in a best possible strengtheningof constraints Dx ≥ d, is applied in order to obtain stronger bounds
from the relaxation. In many practical examples, Dx ≥ d can be separated into K blocks of disjoint sets of variables and constraints,Dκ xκ ≥ dκ , κ = 1, . . . , K, where the dimensions of vectors and matrices are compatible, and each block gives rise to a reformulation.
We refer to Desrosiers and Lübbecke (2005, 2011) for details.
Considering the original time-indexed formulation (1)–(7), a
Dantzig–Wolfe reformulation is not obvious. Note, however, that variables can be grouped “into blocks” according to shifts for all resource
types. That is, variables xjmt form a group for each ﬁxed resource typek and shift I ∈ Ik , containing all jobs j executable during I in all modes
m. For each such group of variables, only those linking constraints (3)and resource constraints (5) for k and all t ∈ I are relevant. These R|
κ = |k=1|Ik | groups of variables and constraints form independent
subproblems Dκ xκ ≥ dκ which are reformulated as sketched above.
All other constraints (2) and (4) are global and kept in the role ofmaster constraints Ax ≥ b. Given k and I, the solution vectors (“extreme points”) of the resulting pricing problem (16)–(24), consistingof (Xj , Sj , Cj )j∈J and r̄k , characterize a conﬁguration ξ . This information,
except the resource usage, entirely encodes a feasible setting of theoriginal xjmt variables; or in other words, the original variable vector
can be expressed as selecting exactly one solution vector (Xj , Sj , Cj )j∈J
and r̄k . Fractionally, this is a convex combination with multipliers xξ ,and the convexity constraint is stated in (13).
It is not necessary to understand the model (8)–(15) as a Dantzig–
Wolfe reformulation of (1)–(7), as our model has a self-supported
interpretation. However, this interpretation motivates the branchingrules we propose as well as the way we make use of cutting planes.
4. Branch-price-and-cut algorithm
4.1. Branching scheme
A solution to the original problem (1)–(7) is given by the resource
capacities Rk , and an assignment of start times Sj and completion
times Cj for each job j. The mode is given by the unique mode mj ,such that pjmj = Cj − Sj . In a modern branch-and-price context one
tries to branch on these original variables which are expressible, via
the Dantzig–Wolfe reformulation, as a sum of master variables, compare (10)–(12), instead of on single master variables (Desrosiers &Lübbecke, 2005). A main reason is that branching decisions on single master variables are hard to respect in the pricing problems, e.g.,a conﬁguration ξ , forbidden via xξ = 0, must not be re-generated.
Branching decisions on (combinations of) master variables represent
additional constraints that give rise to additional dual variables thatneed to be respected in the pricing problem. This usually implies newvariables and new constraints in the pricing problem. However, whenbranching constraints are formulated on the original variables, thesedual variables only affect a subproblem’s objective function, not its
constraints, see e.g., Desrosiers and Lübbecke (2011).
We impose a priority on branching decisions according to the effect we observed in preliminary experiments (Coughlan et al., 2010).
The order is Rk , Sj , and then Cj , from highest to lowest priority. Only
when all higher priority variables assume integral values, we branch
on variables of the next class. Let a superscript star indicate the variable’s value in the current master LP solution. Branching on resource
capacity values has the largest impact as these have non-zero objective coeﬃcients. We ﬁrst evaluate the down-branch Rk ≤ Rk ,and then the up-branch Rk ≥ Rk + 1. Branching on Rk changes thefeasible space of the pricing problems, so we branch on these variables
ﬁrst even if Rk is not fractional. When there is choice between multiple
resource variables, we select a most fractional one and use the objective function coeﬃcient as a tie-breaker (prefer larger values). Start
and completion time variables are considered as branching candidates
only if any corresponding binary conﬁguration variable in (10) is fractional. An unﬁxed start time variable Sj is selected where the points intime corresponding to positive xξ variables differ most, counteracting the smearing of variables. The node is split into two child nodes
with Sj ≤ Sj , and Sj ≥ Sj + 1, respectively. Completion times arehandled accordingly. This scheme is used in conjunction with propagation rules (described next) and aims at creating a more balanced
search tree.
4.2. Propagation
It is a fundamental principle of constraint programming to tighten
variable domains due to logical implications given by the constraints,
and/or already ﬁxed variables. For example, the precedence and resource constraints bear the logical structure of the problem and canbe used to detect infeasible start times of variables which can therefore be eliminated from the domains. This is called propagation. In
branch-and-price, domain propagation can also be used to establish
consistency of the LP relaxation after branching, i.e., already generated master variables which represent conﬁgurations that contradict
any previous branching decision taken on the path that led to thecurrent node need to be locally ﬁxed to zero (variables cannot simply
be eliminated as they are essential in other parts of the search tree).
We mentioned earlier that we only use xjmt variables in a pricing
problem that correspond to feasible combinations of (j, m, t) indices.
These combinations are restricted after branching and propagation,and maybe after cutting. The preprocessing to eliminate infeasiblecombinations is also done in propagation.
In the area of scheduling problems, a large variety of propagation algorithms is known to detect infeasible start times and perform
variable bound adjustments. Edge-ﬁnding is a constraint programming technique concerned with deriving better bounds for earliest
start and latest completion times of jobs using energy arguments. We
adapt the algorithm proposed in Mercier and Van Hentenryck (2008)
to the multi-mode case by using the minimum energy used by anymode for each job, which naturally seems to give weaker bounds. Thisis balanced by the fact that jobs are not preemptive, may not crossshift-bounds and obey precedence constraints which enables further
propagation of start and completion times.
We use this algorithm in every node of the branch-and-price tree,
prior to calling the pricing problem and cutting plane separation.
4.3. Cutting planes
State-of-the-art MIP solvers heavily rely on additional valid inequalities (“cutting planes”) in order to improve the dual bound andby that prune unpromising nodes of the branch-and-bound tree. In
branch-and-price, formulating cutting planes in the master problem
variables is possible, but raises again (as with branching constraints)
the issue of how to respect the additional dual variables in the pricing
problem. It is technically easier to formulate valid inequalities on theoriginal variables and add their Dantzig–Wolfe reformulation, i.e., in
our case their translation to conﬁguration variables, to the masterproblem. Again, this only changes the objective function of the pricing problem, see again Desrosiers and Lübbecke (2011) for details.
Incidentally, this is a good situation for us as the literature knows
several cutting planes for various scheduling problems, all of them
formulated on variables with a meaning as in the standard MIP (1)–(7). Here, we generalize the precedence inequalities as introducedin Christoﬁdes, Alvarez-Valdes, and Tamarit (1987) where jobs have

E.T. Coughlan et al. / European Journal of Operational Research 245 (2015) 70–80

a ﬁxed processing time:∀ (i, j) ∈ E, τ :5. Computational study



xit +

t≥τ

xjt ≤ 1.

(25)

t<τ +pi

These cuts read as follows: For any point in time τ , in an end-tostart precedence relation between jobs i and j, either job i starts thelatest in τ or job j starts before τ + pi . As multiple modes are present
we compute the minimum processing time pi,min a job i may have
locally, i.e., the minimum over all of its modes. Building on these cuts,
we obtain the following cuts in our master problem:∀ (i, j) ∈ E, τ :xξ +ξ :Siξ ≥τxξ ≤ 1.

(i,j)∈E τ ≥tμτ ij +(k,i)∈E τ <t+pim

μτ ki .

In this section we compare the outcome of our branch-price-andcut algorithm against a preprocessed time-indexed formulation (1)–(7) solved by the MIP solver CPLEX. In particular, we are interested in
the dual bounds that are obtained in the root node by either algorithm
since these bounds indicate the strength of the polyhedral description. Furthermore, we give experimental evidence that the branchingstrategy and on some instances a careful use of precedence inequalities play an important role to eﬃciently solve a large number of thebenchmark problems.

(26)
5.1. Benchmark instancesξ :Sjξ ≤τ +pi,min

By construction, this does not entail any structural changes to thepricing problem, but only the coeﬃcients of the objective functionneed to be updated. For each constraint (26) that is added to therelaxation, we introduce the dual variable μτ ij which integrates in
the reduced cost objective function of the pricing problem. If Si ≥ τor Sj < τ + pi,min , the cost coeﬃcient of variable ximt becomes

 

75

(27)

The cutting planes are added to the master problem after new
conﬁguration variables have been generated and each conﬁguration
variable is introduced into the corresponding constraints. Overall, thisgives a branch-price-and-cut algorithm. In Appendix A, we present
these cuts in the context of generalized precedence constraints.
4.4. Primal heuristics
For the master problem, rounding heuristics for fractional solutions are not promising, since values of binary variables may besmeared over the time horizon and precedence constraints are likely
to be violated in a rounded solution. To improve our upper bounds
we extend a leveling heuristic from Megow et al. (2011) that produces initial solutions prior to search. We implement a generic list
scheduling algorithm (a serial schedule generation scheme) and aready scheduling heuristic (a parallel schedule generation scheme,
see Coughlan et al., 2010) which are used as standalone heuristics during the branching process as well as in the initial levelingprocedure.
The two heuristics are called in each node of the search tree. To
this end, we set the capacity of each resource to the fractional variablevalue rounded up, and ﬁx the earliest and latest start and completion
times for each job to the local bounds of the corresponding variables.
We perform list scheduling using the (fractional) LP solution with jobs
sorted by earliest completion times, where each job’s mode is chosen
according to the largest processing time of at most Cj − Sj . If no feasible
solution is obtained we try ready scheduling. Both of these heuristics
produce solutions that are not necessarily feasible w.r.t. the current
primal bound, since resource capacities are rounded up. Regardlessof this, if a feasible schedule is found new columns representing thatschedule are added to the master problem, in order to reduce the total
number of pricing iterations.

There is no publicly available test set of benchmark instances
reﬂecting the setup of our problem. The existing instances in the
PSPLib (Kolisch & Sprecher, 1996) guided the design of our test
set compilation. We use the same generation scheme as before
in Coughlan et al. (2010), now for larger instances. Our set is composed of two sets of job scenarios, with 50 instances each. Each job
can run in three different modes, using one to three units of its resource, with durations ranging from 5 to 12. The ﬁrst set, denotedby N50E70, contains 50 jobs and 70 precedence constraints, whereas
the second set N50E100 contains 50 jobs and 100 precedence constraints. The maximal width W of the precedence graph is six, whichis achieved by constructing W chains of length |J |/W, and randomly
choosing the remaining edges.
There are two different resources which run in ﬁve calendar conﬁgurations, called C1–C5. These calendars are described schematicallyin Fig. 2. In the top row, calendars C1–C3 are shown. In each of these,
the length of the shifts is 60. In C1 and C3 shift breaks are 60 units
long, in C2 only 20 units long. Both resources are available at the sametime in C1, while in C3 availability periods are complementary. In C2
the second resource is offset at 40 units. Calendars C4 and C5 show
shifts with length 20 and breaks having length ﬁve. In C5 one of theresources is offset by 10. All scenarios are tested with each of the
ﬁve different calendars. Time horizons were chosen by computing aminimal and maximal makespan heuristically using an earliest start
list scheduling policy, and averaging these. The two makespans arecomputed by either scheduling all jobs at their highest resource usage(fastest mode) or by assigning the fewest possible resource usage toeach job (slowest mode).
5.2. Experimental setup
All experiments were done on Intel CoreTM i7-870 PCs(2.93 gigahertz, 8 megabytes cache, 8gigabytes memory) runningLinux 2.6.34 (single thread). Each test run had a time limit of 1800 seconds. Our C++ implementation is based on SCIP 2.0.1 (SCIP, 2011) toperform the branch-price-and-cut process, with custom plug-ins forour heuristics, propagation, branching rules, cutting plane separation,and column generation. For the standard MIP (1)–(7) we used CPLEX
12.2 on the same machine, with default parameter settings, again
single thread. Up to two threads were run in parallel on not entirely
idle machines, so run time differences of 5 percent are regarded as“noise.”Fig. 2. Calendar conﬁgurations C1–C3 (top) and C4 and C5 (bottom) used in our test set. Black bars symbolize the temporal location of shifts.

76

E.T. Coughlan et al. / European Journal of Operational Research 245 (2015) 70–80

5.2.1. Usefulness of the heuristics
Previous results showed that it is important for CPLEX and ourbranch-price-and-cut framework to have a good initial solution,
whereas the execution time of the heuristic throughout search is
negligible, see Coughlan et al. (2010).

5.3. Results
We conﬁrm that a Dantzig–Wolfe reformulation can considerably
improve the dual bound as compared to the standard LP relaxation.
As an even stronger statement, we also obtain a dual bound improvement over a state-of-the-art branch-and-cut algorithm. In Fig. 3 we
see that the average improvement per calendar of our conﬁguration
based formulation (8)–(15) over the LP relaxation of (1)–(7) lies between 20 percent and 65 percent throughout all instances. The displayed average deviations and the minimum and maximum values ofthat improvement show the strength of our relaxation.

400
300
200
100
0
C1

C2

C3

C4

C5

Fig. 3. Percentage of improvement of the lower bound obtained from our conﬁguration
based formulation (8)–(15) (no additional cutting planes) over the standard LP relaxation of (1)–(7) as computed by the state-of-the-art MIP solver CPLEX. The minimum(diamond), maximum (triangle), mean, and standard deviation of the improvement inpercent of setting “no cuts” is shown.

100
50
0
−50
C1

C2

C3

C4

C5

Fig. 4. Average mean and standard deviation of root dual bound improvement of ourbranch-price-and-cut approach as compared to CPLEX in percent. The best bound over
all tolerances has been used. Triangles show the best obtained improvement (e.g., 200
percent in C1) and diamonds show the worst lower bound (e.g., 44 percent worse than
CPLEX in C2).

no. solved instances

5.2.3. Settings
We compare the outcome of a CPLEX run on the standard MIP
model (1)–(7) to different settings of our tailored branch-price-andcut approach. In a ﬁrst experiment, we do not separate cuts, denotedby “nocuts.” Then, we evaluate different separation strategies. First,
we separate cuts already in the root node, and secondly, after resource
capacity variables are ﬁxed (by branching or propagation). Checking
for violation only at a promising guess as point in time is denoted by“τgreedy ,” searching for the largest violation cuts is denoted by “τbest .”Furthermore, these cuts are only separated if a certain threshold “tol.”is exceeded, i.e., if in Equations (26) the sum of the variables exceeds 1 + tol. We use thresholds 10−4 , 0.1, 0.3, 0.5 and 0.8, whichrange from a very small violation (gives more cuts) to a very highviolation (gives less cuts).

150

50
40
30
20
10
0
C1

C2

C3

C4

C5

Fig. 5. Number of solved instances for N50E70 with cuts already separated in the root
node; each set of bars from left to right corresponds to rows in Table 1. A CPLEX run on
the standard MIP is black; our branch-and-price algorithm without cuts is white; andthe two groups of different shades of grey show runs of the full branch-price-and-cut
algorithm with precedence cuts enabled, with settings τgreedy and τbest , respectively,
with increasing tolerances (lighter grey is larger tolerance).

no. solved instances

5.2.2. Separation of cuts
There are m · T, i.e., a pseudo-polynomial number of precedence
cuts (26); these are too many to explicitly add them all to the masterproblem, as this would have a negative impact on LP solving times. For
the identiﬁcation of violated cutting planes, which is called separation,
we pursue two approaches. In our ﬁrst approach we only check for aprecedence pair (i, j) at point in time τgreedy = (Ci + Si )/2 whether
Equation (26) holds or is violated by more than a threshold . In thesecond approach, we ﬁnd an optimum point τbest such that the lefthand side of Equation (26) is maximally violated. This can be done bysorting the summands of Equation (26) for each point in time τ suchthat for each value τ the left-hand side can be computed in linear time
(after sorting). This procedure is pseudo-polynomial in the numberof points in time.

200

50
40
30
20
10
C1

C2

C3

C4

C5

Fig. 6. Number of solved instances for N50E100 with cuts already separated in theroot node. See Fig. 5 for the meaning of the bars.

To evaluate the strength of generic cutting planes added in the root
node of the commercial MIP solver, we compare the root dual bounds
in Fig. 4 after all standard cutting planes have been added by CPLEX
versus the root node dual bound in our approach. Our improvement
is no longer that dominating as in Fig. 3 but still for several of the hard
instances in calendar setting C1, there is a 28 percent improvementon the average of the root dual bound compared to CPLEX.
Next, we compare our branch-price-and-cut framework to a standard MIP approach in terms of absolute number of solved instances
and afterward, we discuss the effect of precedence cuts in our modelby evaluating the number of branch-and-bound nodes and the running time. In Fig. 5 and 6 for each calendar, the ﬁrst bar (black) gives
the number of solved instances obtained by CPLEX, the second (white)bar represents this number for our branch-and-price approach without precedence cuts. The ﬁrst ﬁve grey bars symbolize the results forsetting “τgreedy ” in increasing order of tolerances and the last ﬁve bars
stand for settings “τbest ” in increasing order of tolerances.

no. solved instances

E.T. Coughlan et al. / European Journal of Operational Research 245 (2015) 70–80

50
40
30
20
10
0
C1

C2

C3

C4

C5

Fig. 7. Number of solved instances of test set N50E70 if cuts are separated after resource capacity variables are ﬁxed. See Fig. 5 for the meaning of the bars.

no. solved instances

50
40
30
20
10
C1

C2

C3

C4

C5

Fig. 8. Number of solved instances of test set N50E100 if cuts are separated after
resource capacity variables are ﬁxed. See Fig. 5 for the meaning of the bars.

Fig. 5 shows for the instance set N50E70 that the pure branchand-price algorithm without precedence cuts outperforms CPLEX, andadding precedence cuts seems to be a bad idea. Comparing Fig. 5 and6, the set N50E70 with fewer precedence constraints than N50E100 is
the harder one, as expected. In some cases the pure MIP approach iseven better than the branch-price-and-cut approach with additional
precedence cuts. The reason is that the time spent for separating new
cuts, pricing new variables, and the additional LP iterations lead totoo many timeouts and is therefore not competitive. Hence, not using
any precedence cuts at all is, at ﬁrst sight, a good decision.
Nevertheless, it is possible to increase the root dual bound on
several instances by using precedence cuts (26). During root solving
most jobs are able to slide in their time window and the precedence
cuts result in more binary conﬁguration variables that are smeared
over the time horizon. Hence, these cuts should not be used in the root
node but still might be beneﬁcial to prune certain nodes of the branchand-bound tree. Recall that our branching scheme ﬁrst branches on
the resource capacities Rk and afterwardon the start and completiontime variables. Fixed resource capacities in a node already decide on77

a lot of structure for the scheduling problem, since several modes ofa job may no longer be valid. Thus, this seems to be a good point toseparate precedence cuts (26). Figs. 7 and 8 show that more instances
than before can be solved using the precedence cuts. In several casesa tolerance of 10−4 belongs to the best choices for the separation
procedure. Especially, for the harder instances N50E70 some moreinstances of each calendar test set can be solved in the time limit toproven optimality in contrast to CPLEX or a setting without additional
precedence cuts.
Next, we elaborate on the solution time and on the number ofnodes needed to ﬁnd an optimal solution and prove its optimality.
The results on the total number of solved instances showed that itis beneﬁcial to separate precedence cuts after resource variables are
ﬁxed. Therefore, we will only present the numbers for that case.
Tables 1 and 2 reveal that using additional precedence cuts enables
us to decrease the number of tree nodes by 10 percent to 20 percent
on average. For several hard instances, e.g., in calendar C5 in test set
N50E70 a decrease by even 50 percent is possible. Best results interms of nodes are obtained when τbest is used for the separation ofcuts. Nevertheless, this does not carry over to a reduced running time.
For C1 the running times increase (τbest vs. τgreedy ), whereas for C5
it decreases and the running time is about 10 percent faster than if
no precedence cuts are separated. That is, there is the usual trade-off
between quality and time, and cutting planes may be most interesting
in memory critical applications.
A tolerance between 0.1 and 0.3 gives the overall best results as
higher and lower tolerances usually increase the running times.
Tables 1 and 2 only show a slight improvement considering the
average running time and the number of nodes needed (in the shifted
geometric mean) if precedence inequalities are separated after resource capacities are ﬁxed. The performance proﬁles (Dolan & Moré,
2002) in Fig. 9 give a more detailed comparison by comparing theratios of the number of nodes (running time) needed per instanceversus the best running time by any of the settings listed in the tables. We see that, e.g., for calendar C5 on more than 50 percent of theinstances the precedence inequalities remarkably reduce the numberof nodes needed. On several instances, the setting without precedence
inequalities needs between 10 and 100 times more nodes than thebest setting with precedence inequalities. But we also see that the reduction in terms of running times is much smaller as separating thesecuts and computing the new objective coeﬃcients in the pricing problems is costly, too. Similarly, for calendar C4 we observe a decreasein the number of nodes needed, whereas the reduction in terms of
running time is rather modest and can only be seen on less than 20
instances. We do not elaborate on the results for calendars C1–C3
here. On calendars C1 and C2, results are similar to C5, whereas on C3
not much changes, as on these instances the dual bounds have notbeen a bottleneck.

Table 1Comparison of tree nodes and running time in seconds for N50E70—Cuts are only separated after resources are ﬁxed.
Means are computed over the instances solved by all settings, of which there are 27, 14, 48, 26, and 31 instances forcalendars C1–C5.
C1

C2

C3

C4

C5

Tol.

Nodes

Time

Nodes

Time

Nodes

Time

Nodes

Time

Nodes

Time

τgreedy

0.0
0.1
0.3
0.5
0.8

12.0
10.0
10.2
10.6
11.4
11.4

660.7
667.9
662.7
632.6
661.0
611.5

10.7
9.5
8.8
9.0
10.5
10.7

546.3
553.4
518.7
523.4
538.9
560.7

1.0
1.0
1.0
1.0
1.0
1.0

6.4
6.4
6.4
6.4
6.4
6.3

11.4
8.9
8.3
10.0
11.2
11.1

200.5
201.0
197.2
204.2
216.9
200.4

21.3
14.3
16.7
15.9
19.1
19.7

356.2
346.3
361.8
338.9
347.2
351.5

τbest

0.0
0.1
0.3
0.5
0.8

9.2
8.9
9.2
10.1
11.1

704.6
673.1
653.6
660.7
608.4

8.3
8.4
8.6
9.8
10.3

536.0
559.8
542.7
560.6
528.2

1.0
1.0
1.0
1.0
1.0

6.4
6.4
6.4
6.4
6.3

7.4
7.9
7.7
9.1
12.4

196.4
202.8
197.6
201.1
212.4

12.0
11.6
10.9
15.6
16.6

353.6
341.3
323.3
337.6
338.9

No cuts

78

E.T. Coughlan et al. / European Journal of Operational Research 245 (2015) 70–80
Table 2Comparison of tree nodes and running time in seconds for N50E100—Cuts are only separated after resources are ﬁxed.
Means are computed over the instances solved by all settings, of which there are 42, 44, 50, 48, and 47 instances forcalendars C1–C5 respectively.
C1

C2

C3

C4

C5

Tol.

Nodes

Time

Nodes

Time

Nodes

Time

Nodes

Time

Nodes

Time

τgreedy

0.0
0.1
0.3
0.5
0.8

5.3
5.4
4.8
5.4
5.1
5.4

31.0
32.6
31.1
30.5
30.7
29.7

9.8
8.0
8.2
9.1
9.1
10.3

88.8
83.6
86.9
87.3
85.5
88.1

1.1
1.1
1.1
1.1
1.1
1.1

1.0
1.0
1.0
1.0
1.0
0.9

8.4
6.3
6.4
6.5
6.8
7.4

18.1
16.4
16.5
16.3
16.5
16.4

8.2
7.7
7.3
7.7
8.2
9.0

23.8
25.5
24.3
24.1
24.0
23.6

τbest

0.0
0.1
0.3
0.5
0.8

4.8
4.8
4.7
5.1
5.4

32.3
32.5
30.5
31.4
29.8

7.5
7.2
7.6
7.7
10.0

88.8
86.8
87.3
83.2
92.2

1.1
1.1
1.1
1.1
1.1

1.0
1.0
1.0
1.0
0.9

5.9
5.7
5.8
6.4
7.2

16.5
16.4
16.3
16.5
16.2

6.3
6.2
6.5
7.8
8.5

23.5
23.0
22.9
24.0
23.3

No cuts

100
τbest -0.0
τbest -0.1
τbest -0.3
τbest -0.5
τbest -0.8
τgreedy -0.0
τgreedy -0.1
τgreedy -0.3
τgreedy -0.5
τgreedy -0.8
nopc

75
50
25
0
100

100.2

100.4
ratios% of instances% of instances

100

50
25
0

100.6

100

(a) Comparison of the running times for C5.

101
ratios102(b) Comparison of the nodes for C5.

100
τbest -0.0
τbest -0.1
τbest -0.3
τbest -0.5
τbest -0.8
τgreedy -0.0
τgreedy -0.1
τgreedy -0.3
τgreedy -0.5
τgreedy -0.8
nopc

75
50
25
0
100

100.2

100.4
ratios

100.6

100.8

(c) Comparison of the running times for C4.

% of instances

100% of instancesτbest -0.0
τbest -0.1
τbest -0.3
τbest -0.5
τbest -0.8
τgreedy -0.0
τgreedy -0.1
τgreedy -0.3
τgreedy -0.5
τgreedy -0.8
nopc

75

τbest -0.0
τbest -0.1
τbest -0.3
τbest -0.5
τbest -0.8
τgreedy -0.0
τgreedy -0.1
τgreedy -0.3
τgreedy -0.5
τgreedy -0.8
nopc

75
50
25
0
100

101
ratios102(d) Comparison of the nodes for C4.

Fig. 9. Performance proﬁles: Ratios of nodes and running times for all optimally solved instances from N50E70 for calendar C5 on top and for calendar C4 below.

Table 3 indicates the computational efforts spent in the pricing
problems. For each set, N50E70 and N50E100, and per calendar C1–
C5, we list the average running time (avtime) and number of nodes(avnodes) per instance. Column “ptime/time” shows that solving thepricing problem (setting up, solving the MIP, and generating newcolumns) takes from 58 percent up to 99 percent of the total solving
time. Hence, the pricing routine determines the overall execution
time. Recall that in each pricing round, a MIP for each shift needs tobe solved. Column “# calls per node” shows how often the pricer iscalled on average per node, while the number of variables generated
per call is indicated in the last column. The number of shifts per
instance (which is equal to the number of MIPs to be solved with eachcall of the pricing problem) and the number of pricing rounds whensolving a node inﬂuence this number.
Table 4 shows for each set and per calendar how often per instance that was not solved to optimality (“nsubopt”) CPLEX hit thebest known primal (“nbestprimal”) and dual (“nbestdual”) bounds.

The succeeding columns show the average gap of the best primal(dual) solution versus the best known solution (dual bound) from all
our branch-price-and-cut variants, denoted by “avpgap” (“avdgap”).
The last column presents the average gap that CPLEX was able to arrive
at. From that table we can see that for set N50E70, the primal gap isthe hardest for CPLEX, 13 to 23 percent. In contrast, for set N50E100
the average dual gap is much bigger than the average primal gap.
Here, the average dual gap compared to all settings using branchprice-and-cut ranges from 13 percent to 37 percent, while the primal
gap falls between 13 and 17 percent.
6. Summary and conclusions
We studied a practically motivated, complex multi-mode projectscheduling problem with the objective to minimize the resource availability cost. Our aim was to demonstrate how far one can get when insisting on optimal solutions. The literature on related problems and inE.T. Coughlan et al. / European Journal of Operational Research 245 (2015) 70–80

79

Table 3
Share of time spent in the pricing problem (ptime/time), number of calls per node and number ofvariables generated per call on average over all instances per set. A call to the pricer triggers thepricing for each shift.
Set

Calendaravtimeavnodes

ptime/time# calls/node# vars/call

N50E70

C1
C2
C3
C4
C5

1230.69
1424.4951.72
843.63922.76

14.46
40.77
29.76
313.36255.53

0.99
0.97
0.97
0.88
0.92

8.67
10.08
5.46
5.71
4.30

18.44
10.72
9.01
6.63
7.76

N50E100

C1
C2
C3
C4
C5

312.61
273.46
3.06
74.46
152.70276.7076.16
1.13
491.68519.20

0.76
0.93
0.99
0.64
0.56

3.76
6.65
6.68
2.33
2.70

5.73
8.35
22.84
7.13
3.86

Table 4
Evaluation of primal and dual bounds for all instances that could not be solved to optimality by CPLEX.
Set

Calendarnsubopt

nbestprimal

nbestdual

avpgapavdgapavcplexgap

N50E70

C1
C2
C3
C4
C5

37
47
1
42
33

14
23
0
13
8

17
28
0
16
19

0.18
0.13
0.16
0.23
0.22

0.12
0.09
0.50
0.17
0.06

0.34
0.31
0.56
0.34
0.27

N50E100

C1
C2
C3
C4
C5

28
33
0
38
15

14
14
0
25
6

3
6
0
2
6

0.13
0.17
–
0.10
0.16

0.37
0.23
–
0.36
0.13

0.47
0.34
–
0.42
0.24

particular the PSPLib with up-to-date benchmark results suggest thateven today instances with only 30 jobs are hard to solve to optimality
when a “standard” model is used; in our case we compare against an
adapted classical MIP formulation (Pritsker et al., 1969). Our study
demonstrates that it helps to decompose the problem resource-wise
per availability period. Technically, this means applying a Dantzig–
Wolfe reformulation to the standard formulation (in a non-obviousway), even though this is not necessary to understand our approach.
The resulting model is based on variables that represent entire schedules for each resource and shift, of which there are exponentiallymany. This necessitates a solution approach via a column generation
and branch-and-price algorithm. A usually complicating feature inproject scheduling, the resource availability periods (or shift calendars), can even be turned into an advantage in this context, as thisreduces the size of the variable-generating subproblems in column
generation.
In our computational study we showed that the LP relaxationsof the reformulated model are much stronger than that of the standard approach. The dual bound can even be strengthened by applying
additional valid inequalities which are known for the standard formulation, and which we translate to our model. We believe that this isprototypical, and that more inequalities known from the scheduling
literature could be integrated in our approach. The resulting branchprice-and-cut algorithm enables us to solve a large amount of instances with 50 jobs to proven optimality, where state-of-the-art
MIP solvers suffer from the size and weak relaxation of the standard
model.
Our algorithm is generic and components like the pricing problem
may also be solved via constraint programming algorithms or partiallywith a heuristic. We believe that the general approach is well-suited
for similar problems, in particular, when the objective function is“complicated.” We ﬁnally remark that we believe that our ability toprovide solution algorithms to complex models as ours changes ourway we model with integer variables. While classical IP modelingused variables with “little meaning” our study supports the advise toimpose more structure on a single variable’s meaning.

Appendix A. Cuts for generalized precedence inequalities
We generalize the precedence inequalities as introducedin Christoﬁdes et al. (1987) where jobs have a ﬁxed processing time
to our master problem also for the case of generalized precedence
relations. This is interesting for modeling multi-resource jobs, among
others. The cuts of Christoﬁdes et al. are given as follows:∀ (i, j) ∈ E, τ :



xit +

t≥τ



xjt ≤ 1.

(A.1)

t<τ +pi

These cuts read as follows: For any point in time τ , in an endto-start precedence relation between jobs i and j, either job i starts
the latest in τ or job j starts before τ + pi . Building on these cuts, westudy generalized precedence constraints that are of the form startto-start (GP1 ), end-to-start (GP2 ), start-to-end (GP3 ) and end-to-end(GP4 ). Let δij ∈ Z, then we can include these generalized precedence
relations in our master problem by adding the following inequalities
based on start and completion time variables (Sj , Cj ), and also separate
strengthened inequalities, similar to (A.1), based on the conﬁguration
variables (xξ ).

(i, j) ∈ GP1 ⇔ Si + δij ≤ Sj or:xξ +ξ :Siξ ≥τxξ ≤ 1∀τξ :Sjξ <τ +δij

(i, j) ∈ GP2 ⇔ Ci + δij ≤ Sj or:xξ +ξ :Ciξ ≥τxξ ≤ 1∀τξ :Sjξ <τ +δij

(i, j) ∈ GP3 ⇔ Si + δij ≤ Cj or:xξ +ξ :Siξ ≥τxξ ≤ 1∀τξ :Cjξ <τ +δij

(i, j) ∈ GP4 ⇔ Ci + δij ≤ Cj or:ξ :Ciξ ≥τxξ +ξ :Cjξ <τ +δij

xξ ≤ 1∀τ

80

E.T. Coughlan et al. / European Journal of Operational Research 245 (2015) 70–80

Using these cuts in our master problem, the nature of the pricing problem (i.e., its constraints) does not change, only the coeﬃcients of the objective function need to be updated. For each cut
encoded by (τ , (i, j)) ∈ GP we need to add the value of the dual variable μτ ij to the objective function, for all types of generalized precedence relations GP , = 1, . . . , 4 according to the following rules. The
cost coeﬃcient of variable ximt in the pricing problem is increased
by:
GP1 :





(i,j)∈GP1 τ :t≥τ

GP2 :



μ1τ ij +(i,j)∈GP2 τ :t+pim ≥τ

GP3 :





(i,j)∈GP3 τ :t≥τ

GP4 :





(i,j)∈GP4 τ :t+pim ≥τ(k,i)∈GP1 τ :t<τ +δki

μ2τ ij +μ3τ ij +

μ1τ ki(k,i)∈GP2 τ :t<τ +δki(k,i)∈GP3 τ :t+pim <τ +δki

μ4τ ij +

μ2τ kiμ3τ ki(k,i)∈GP4 τ :t+pim <τ +δki

μ4τ ki

References
Achterberg, T. (2009). SCIP: Solving constraint integer programs. Mathematical Programming Computation, 1(1), 1–41.
Bianco, L., & Caramia, M. (2011). A new lower bound for the resource-constrained
project scheduling problem with generalized precedence relations. Computers andOperations Research, 38(1), 14–20.
Brucker, P., & Knust, S. (2000). A linear programming and constraint propagation-based
lower bound for the RCPSP. European Journal of Operational Research, 127(2), 355–
362.
Christoﬁdes, N., Alvarez-Valdes, R., & Tamarit, J. M. (1987). Project scheduling withresource constraints: A branch and bound approach. European Journal of Operational
Research, 29(3), 262–273.

Coughlan, E. T., Lübbecke, M. E., & Schulz, J. (2010). A branch-and-price algorithm for
multi-mode resource leveling. In P. Festa (Ed.), Experimental algorithms. In Lecture
Notes in Computer Science: 6049 (pp. 226–238). Lecture Notes in Computer Science,
vol. 6049. Springer.
Demeulemeester, E. (1995). Minimizing resource availability costs in time-limited
project networks. Management Science, 41, 1590–1598.
Demeulemeester, E. L., & Herroelen, W. S. (2002). Project scheduling: A research handbook. Kluwer.
Desrosiers, J., & Lübbecke, M. E. (2005). A primer in column generation. In G.
Desaulniers, J. Desrosiers, & M. M. Solomon (Eds.), Column generation (pp. 1–32).
Berlin: Springer-Verlag.
Desrosiers, J., & Lübbecke, M. E. (2011). Branch-price-and-cut algorithms. In
J. J. Cochran (Ed.), Encyclopedia of operations research and management science. Chichester: John Wiley & Sons. http://onlinelibrary.wiley.com/doi/10.1002/
9780470400531.eorms0118/abstract.
Dolan, E. D., & Moré, J. J. (2002). Benchmarking optimization software with performance
proﬁles. Mathematical Programming, 91, 201–213.
Drexl, A., & Kimms, A. (2001). Optimization guided lower and upper bounds for the
resource investment problem. The Journal of the Operational Research Society, 52(3),340–351.
Franck, B., Neumann, K., & Schwindt, C. (2001). Project scheduling with calendars. OR
Spektrum, 23, 325–334.
Hartmann, S. (2001). Project scheduling with multiple modes: A genetic algorithm.
Annals of Operations Research, 102(1–4), 111–135.
Hartmann, S., & Briskorn, D. (2010). A survey of variants and extensions of the resourceconstrained project scheduling problem. European Journal of Operational Research,
207(1), 1–14.
Kolisch, R., & Sprecher, A. (1996). PSPLIB—A project scheduling problem library. European Journal of Operational Research, 96, 205–216. http://129.187.106.231/psplib/,
last accessed 2011/05/26
Megow, N., Möhring, R. H., & Schulz, J. (2011). Decision support and optimization in
shutdown and turnaround scheduling. INFORMS Journal of Computing, 23(2), 189–
204.
Mercier, L., & Van Hentenryck, P. (2008). Edge ﬁnding for cumulative scheduling. INFORMS Journal of Computing, 20(1), 143–153.
Möhring, R. H. (1984). Minimizing costs of resource requirements in project networks
subject to a ﬁxed completion time. Operations Research, 32(1), 89–120.
Neumann, K., Schwindt, C., & Zimmermann, J. (2003). Project scheduling with time windows and scarce resources. Springer.
Pritsker, A. A. B., Watters, L. J., & Wolfe, P. M. (1969). Multi project scheduling withlimited resources: A zero-one programming approach. Management Science, 16,
93–108.
Santos, M. A., & Tereso, A. P. (2011). On the multi-mode, multi-skill resource constrained project scheduling problem a software application. In A. Gaspar-Cunha,
R. Takahashi, G. Schaefer, & L. Costa (Eds.), Advances in Intelligent and Soft Computing: 96. Soft computing in industrial applications (pp. 239–248). Berlin/Heidelberg:
Springer. doi:10.1007/978-3-642-20505-7_21.
SCIP (2011). Solving constraint integer programs. http://scip.zib.de/. Version 2.0.1.
Zhan, J. (1992). Calendarization of time planning in MPM networks. ZOR – Methods andModels for Operations Research, 36(5), 423–438.