A compact reformulation of the two-stage robust resource-constrained project scheduling problem

This paper considers the resource-constrained project scheduling problem with uncertain activity durations. We assume that activity durations lie in a budgeted uncertainty set, and follow a robust two-stage approach, where a decision maker must resolve resource conflicts subject to the problem uncertainty, but can determine activity start times after the uncertain activity durations become known. We introduce a new reformulation of the second-stage problem, which enables us to derive a compact robust counterpart to the full two-stage adjustable robust optimisation problem. Computational experiments show that this compact robust counterpart can be solved using standard optimisation software significantly faster than the current state-of-the-art algorithm for solving this problem, reaching optimality for almost 50% more instances on the same benchmark set.


Introduction
The resource-constrained project scheduling problem (RCPSP) consists of scheduling a set of activities, subject to precedence constraints and limited resource availability, with the objective of minimising the overall project duration, known as the makespan. Given its practical relevance to a number of industries, including construction (Kim, 2013), manufacturing (Gourgand et al., 2008), R&D (Vanhoucke, 2006), and personnel scheduling (Drezet and Billaut, 2008), the RCPSP and many of its variants have been widely studied since a first model was introduced by Pritsker et al. (1969). The vast majority of this research, however, has examined the RCPSP under the assumption that the model parameters are known deterministically (for a survey of the deterministic RCPSP, see Artigues et al. (2008)), but clearly, in practice, large projects are subject to non-trivial uncertainties. For instance, poor weather might delay construction times, uncertain delivery times of parts may delay manufacturing activities, and the duration of research activities are inherently uncertain. As a result, in recent years, increasing attention has been given to the uncertain RCPSP, where scheduling decisions must be made whilst activity durations are unknown.
There exist two main approaches to solving the uncertain RCPSP. The first is to view the problem as a dynamic optimisation problem where scheduling decisions are made each time new information becomes available according to a scheduling policy (Igelmund and Radermacher, 1983a,b;Möhring and Stork, 2000). Most recently, Li and Womer (2015) use approximate dynamic programming to find an adaptive closed-loop scheduling policy for the uncertain RCPSP.
The second approach aims to proactively develop a robust baseline schedule that protects against delays in the activity durations. Zhu et al. (2007) present a two-stage stochastic programming formulation for building baseline schedules for projects with a single resource. Bruni et al. (2015) present a chance-constraint-based heuristic for constructing robust baseline schedules and Lamas and Demeulemeester (2016) introduce a procedure for generating robust baseline schedules that does is independent of later reactive scheduling procedures. For a review of both dynamic and proactive project scheduling, see Herroelen and Leus (2005).
Although frequently referred to as robust, none of the scheduling methods described above make use of robust optimisation as defined in Ben-Tal and Nemirovski (1998, 1999. Over the last 20 years, robust optimisation has emerged as an effective framework for modelling uncertain optimisation problems. Unlike stochastic programming, robust optimisation does not require probabilistic knowledge of the uncertain data. Instead the robust optimisation approach only assumes that the uncertain data lie somewhere in a given uncertainty set, and then aims to find solutions that are robust for all scenarios that can arise from that uncertainty set. The applicability of robust optimisation as a method for solving uncertain optimisation problems has increased following the introduction of adjustable robust optimisation (Ben-Tal et al., 2004;Yanıkoglu et al., 2019). Adjustable robust optimisation extends static robust optimisation into a dynamic setting, where a subset of the decision variables must be determined under uncertainty, whilst other variables can be adjusted following observations of the uncertain data. As well as accurately modelling the decision process undertaken by many real-world decision-makers, adjustable robust optimisation overcomes the overconservativeness that restricts the applicability of static robust optimisation models. For extensive surveys on robust optimisation, see Ben-Tal et al. (2009);Bertsimas et al. (2011);Gorissen et al. (2015); Goerigk and Schöbel (2016).
Despite the successful application of robust optimisation in many different fields (see Bertsimas et al. (2011)), so far robust optimisation has been little used in project scheduling. To the best of our knowledge, to date, only three papers have directly applied robust optimisation in the construction of robust baseline project schedules. Artigues et al. (2013) present an iterative scenario-relaxation algorithm for the uncertain RCPSP with the objective of minimising the worst-case absolute regret (Kouvelis and Yu, 1997). Bruni et al. (2017) introduce a two-stage adjustable robust optimisation model with the objective of minimising the worst-case makespan. For the case of budgeted uncertainty, this model is solved using a Benders'-style decomposition approach (Benders, 1962). Most recently, Bruni et al. (2018) present a computational study of solution methods for solving the two-stage adjustable RCPSP. An additional Benders'-style algorithm is compared against a primal decomposition algorithm, as well as the algorithm presented in Bruni et al. (2017). The primal decomposition algorithm is shown to be the best performing algorithm for solving the two-stage adjustable RCPSP. This paper presents a new compact reformulation of the two-stage adjustable robust RCPSP with budgeted uncertainty. Computational experiments show that this compact reformulation can be solved using standard optimisation software significantly faster, and for a much greater number of instances than the current best algorithm for solving this problem.
The remainder of this paper is organised as follows: Section 2 introduces the two-stage adjustable robust RCPSP in detail, before Section 3 derives a compact reformulation of this problem and computational experiments are presented in Section 4. Concluding remarks are made in Section 5.

The two-stage robust RCPSP
A project consists of a set V = {0, 1, . . . , n, n + 1} of non-preemptive activities, where 0 and n + 1 are dummy source and sink activities with duration 0. Each activity i ∈ V requires an amount r ik ≥ 0 of resource k ∈ K, where K is the set of project resource types. Each resource k ∈ K has a finite availability R k in each time period. Each activity i ∈ V has a nominal duration given byθ i , and a worst-case duration given byθ i +θ i , whereθ i is its maximum deviation. In addition to resource constraints, the project activities must be scheduled in a manner that respects a set E of strict finish-to-start precedence constraints.
A project can be represented on a directed graph G(V, E). An example project involving seven non-dummy activities and a single resource is shown in Figure 1. We assume that the duration of each activity i ∈ V lies somewhere between its nominal valueθ i and its worst-case valueθ i +θ i . Additionally, we follow Bertsimas and Sim (2004) and assume that only a subset of all activities can simultaneously attain their worst-case values. Hence, the set in which we assume durations can lie, known as the uncertainty set, is given by where Γ determines the robustness of the solution by controlling the number of activities that are allowed to reach their worst-case duration simultaneously. For Γ = 0, each activity takes its nominal duration and the problem reduces to the deterministic RCPSP. At the other extreme, when Γ = n, every activity can take its worst-case duration, and this uncertainty set becomes equivalent to interval uncertainty.
The robust RCPSP lends itself naturally to a two-stage decision making process, where resource allocation decisions need to be made at the start of the project, before the uncertain activity durations become known, but the activity start times can be decided following the realisation of the activity durations. Hence, resource allocation decisions constitute the set of first-stage decisions, whilst the activity start times constitute the set of second-stage decisions.
More specifically, the first-stage resource allocation decisions consist of determining a feasible extension of the project precedence relationships E so that all resource conflicts are resolved. A forbidden set (Igelmund and Radermacher, 1983a) is any subset F ⊆ V of non-precedence-related activities such that i∈F r ik > R k for at least one k ∈ K, i.e. the activities of F cannot be processed simultaneously without violating a resource constraint.
A minimal forbidden set is a forbidden set that does not contain any other forbidden set as a subset. We denote the set of minimal forbidden sets by F. For the example project in to the project network. Bartusch et al. (1988) show that solving the RCPSP is equivalent to finding an optimal choice of additional precedence relationships X ⊆ V 2 \ E, such that the extended project network G (V, E ∪ X) is acyclic and contains no forbidden sets. Such an extension X to the project precedence network is referred to as a sufficient selection.
Hence, a solution to the first-stage problem corresponds to the choice of a sufficient selection X. Figure 2 shows the extended project network for a sufficient selection to the example project shown in Figure 1 (arcs in X are dashed).
Given the extended project network resulting from the choice of sufficient selection made in the first stage, the second stage problem consists of determining activity start times in order to minimise the worst-case makespan in this extended network. Since all resource conflicts have been resolved in the first-stage problem, the second stage problem contains no resource constraints. Hence, the two-stage robust RCPSP under budgeted uncertainty is given by: where X is the set of sufficient selections and S(X, θ) denotes the set of feasible activity start times, depending on activity durations θ ∈ U(Γ), as well as on the choice of sufficient selection X. We have To solve this problem we propose a mixed-integer programming formulation, outlined in the following section.

A compact reformulation
In this section, we present a compact reformulation of the two-stage robust RCPSP. We begin by first examining the adversarial sub-problem of maximising the worst-case makespan for a given sufficient selection.

The adversarial sub-problem
Suppose the solution to the first-stage problem provides a sufficient selection X ∈ X , and is given by a vector y ∈ {0, 1} V ×V where The second-stage sub-problem that arises can be considered from the point of view of an adversary who wishes to choose the worst-case scenario of delays for the given first-stage solution. Following the adversary's choice of delays, we can determine the start time of each activity in order to minimise this worst-case makespan.
Let us assume a fixed scenario θ ∈ U(Γ) given by the vector δ ∈ [0, 1] |V | . In this case, the inner minimisation problem becomes where M is some number greater than or equal to the maximum possible minimum makespan.
Taking the dual of this inner minimisation problem, we can find the following non-linear mixed-integer programming formulation for the adversarial sub-problem, first introduced in Bruni et al. (2017): This can be viewed as a longest-path problem, where up to Γ units of delay can be distributed among activities by the adversary.
As shown by Bruni et al. (2017), this model can be linearised as follows: It is claimed in Proposition 4 of Bruni et al. (2017) that this problem is equivalent to its linear relaxation, where α ij ∈ [0, 1] for all (i, j) ∈ V 2 . This, however, is not the case, as the following counter-example demonstrates. Figure 3 shows a project with three non-dummy activities, each with a nominal duration ofθ i = 1, and a maximum deviation ofθ i = 1, i = 1, 2, 3. Suppose a feasible first-stage solution has been found, resulting in the network shown in Figure 3. We consider this problem from the point of view of the adversary, who wishes to distribute up to Γ = 1 units of delay, in order to maximise the minimum makespan. If (13)-(22) is equivalent to its linear relaxation, then the adversary gains no advantage by choosing α ∈ (0, 1) and splitting the unit flow on its route from the source-node 0 to the sink-node 4. However, as can be seen with this example, the adversary does in fact obtain an advantage.
Note that Bruni et al. (2017) attempt to prove that model (13)-(22) is equivalent to its linear relaxation, and therefore polynomially solvable, by showing that the corresponding constraint matrix is totally unimodular. In Appendix A, we identify an error with this proof and show that the constraint matrix is not totally unimodular. This result is consistent with the above counter-example.
Since problem (13)-(22) is not equivalent to its linear relaxation, we cannot apply strongduality to get an equivalent minimisation problem. Therefore, in order to obtain a compact reformulation of the two-stage robust RCPSP, an alternative reformulation of the adversarial sub-problem is required.
A dynamic programming procedure for solving problem (13)-(22) when Γ ∈ Z is presented in Bruni et al. (2017). This procedure works by considering Γ + 1 paths from the source node 0 to node i, for each i ∈ V , where each path π γ i , γ = 0, . . . , Γ, is characterised by the inclusion of exactly γ delayed activities. Given a path π γ i , its extension to each successor node j ∈ Succ i is evaluated by considering two possibilities: either the succes-    (22) is not equivalent to its linear relaxation.
sor activity j is delayed, resulting in the path π γ+1 j , or it is not delayed, resulting in the path π γ j . Hence, the dynamic programming algorithm has a state ST (j, γ) for each node j at level γ, and the value of each state V (ST (j, γ)) is computed through the following recursion: This dynamic programming algorithm can be viewed as finding the critical path through the augmented project network built from Γ + 1 copies of the original project network (an example of such a network is shown in Figure 4). The inclusion of an inter-level arc, e.g. a dashed arc in Figure 4, in the critical path corresponds to the delay of the activity at the origin of that arc.
Since the second stage problem is simply a longest-path problem on this augmented network, it can be recast into the following mixed-integer linear program: where α ijγ is the flow from node i to node j in level γ and β ijγ is the flow from node i in level γ − 1 to node j in level γ. The constraints model a unit flow through the augmented network from node 0 in level 0 (Constraint (30)) to node n + 1 in level Γ (Constraint (31)).
Constraints (27) are flow-conservation constraints that ensure that for node each in level γ = 1, . . . , Γ − 1, the incoming flow from levels γ and γ − 1 must be equal to the outgoing flow to levels γ and γ + 1. Constraints (28) and (29) conserve flow over the nodes in the special cases of the first and last level, respectively. Note that this model includes more α ijγ and β ijγ variables than indicated in Figure 4, with the edges shown in Figure 4 corresponding to the edges for which y ij = 1. The edges that are not shown are penalised by constant M in the objective (26) when y ij = 0. To ensure that it is always possible to find a path from node 0 in level 0 to node n + 1 in level Γ in the augmented network (if Γ is larger than the number of activities included in the longest path from node 0 to node n + 1 in the original project network, such a path may not be possible), the final sink nodes of each layer are connected by enforcing y n+1,n+1 = 1 (see dotted arcs in Figure 4). Sinceθ n+1 +θ n+1 = 0 these additional edges can be traversed at no extra cost to reach node n + 1 in level Γ.
The first-stage problem aims determine a sufficient selection X ∈ X that minimises the objective value of the second-stage objective value. This first-stage problem can be modelled with a flow-based formulation, as proposed by Artigues et al. (2003). This formulation makes use of continuous resource flow variables f ijk , which determine the amount of resource type k ∈ K that is transferred upon the completion of activity i to activity j.
Additionally, binary variables y ij capture the choice of sufficient selection by representing precedence relationships of the extended project network.
Thus, having dualised the second-stage problem (26)-(33) into a minimisation problem, the first and second-stages can be combined to obtain the following compact reformulation of the full two-stage robust RCPSP with budgeted uncertainty: s.t. S 00 = 0 (35) where M , as before, is chosen to be greater than or equal to the maximum possible minimum makespan, and N k is some number greater than or equal to R k . Constraints (35) Constraints (39) ensure that resource flow respects precedence relationships, and constraints (40) and (41)

Computational experiments
This section compares results obtained by solving the compact robust counterpart (34)-(44), and three slight extensions to this model, with the current state-of-the-art approach to solving the two-stage robust RCPSP proposed in Bruni et al. (2018). Before outlining the proposed extensions to the basic model detailed in the previous section, we introduce the test instances used in this computational study.

Instances
The test instances used in this computational study have been converted from deterministic RCPSP instances involving 30 activities, taken from the PSPLIB (Kolisch and Sprecher (1997), http://www.om-db.wi.tum.de/psplib/). The difficulty of these instances is measured and controlled by the following three parameters:

Implementations
The following section compares the performance of model (34)-(44) with that of three slight extensions to this model. Here, we outline these extensions and clarify details regarding the practical implementation of these models.
The first variant of the basic model (34)-(44) includes the following transitivity constraints on the y-variables: As explained in Section 3.2, these transitivity of the y-variables is not strictly necessary to ensure the feasibility of the activity start-times. We include them as an extension to model (34)-(44) in order to assess their impact on the computational performance of the model.
The second extension involves the provision of a heuristic warm-start solution to the solver software. This heuristic solution is obtained with the following procedure: 1. Given an uncertain RCPSP instance, a heuristic solution is found to the corresponding deterministic instance using the latest-finish-time (LFT) priority-rule heuristic (Kolisch, 1996).
2. From this solution, a feasible set of y-variables is obtained by setting where s j is the start time of activity j, and f i is the finish time of activity i.
3. These y variables are passed to the basic model (34)-(44), which is solved to provide a feasible warm-start solution.
This warm-start solution can be used to tighten the big-M constraints (36) and (37), and thereby further improve the basic model. This is achieved by setting M ij = LF * i − ES j for each (i, j) ∈ V 2 , where ES j is the earliest start time of activity j, and LF * i is the latest finish time of activity i, calculated relative to the makespan of the warm-start solution.
These values are computed recursively via a forward-pass and backward-pass of the project network, respectively.
Note that, although the S-variables of model (34)-(44) are in general continuous, for the purposes of this computational study, the S-variables have been set to be integer. Sinceθ = θ /2 ∈ Z for the instances solved in this study, the correctness of the model is unaffected by this specification.
In summary, the following section presents results from the following five solution ap-  Bruni et al. (2018). This is the strongest existing approach for solving the two-stage robust RCPSP.
All the models proposed in this paper have been solved using Gurobi 9.0.1, running on 4 cores of a 2.30GHz Intel Xeon CPU, limited to 16GB RAM. Note that the specifications of this machine have been chosen to be as similar as possible to that of the CPU used in the experiments performed in Bruni et al. (2017) and Bruni et al. (2018). A limit of 20 minutes was imposed on the solution time of each model, the same as used for the experiments performed in Bruni et al. (2017) and Bruni et al. (2018). Results for the primal method have been reproduced from Bruni et al. (2018).

Results
In this section, we first present and analyse results from solving model (34)-(44) and the three variants proposed in the previous section, before we compare these results with those from the current best iterative algorithm presented in Bruni et al. (2018).
We start by considering the performance profile (Dolan and Moré, 2002) plot shown in Figure 5. The performance profile uses the performance ratio as a measure by which the different models can be compared. The performance ratio of model m ∈ M for problem instance i ∈ I is defined to be where t im is the time required to solve instance i using model m. If model m is unable to solve instance i to optimality within the 20 minute time-limit, then p im = P , where P ≥ max i,m r im . The performance profile of model m ∈ M is defined to be the function i.e. the probability that the performance ratio of model m is within a factor τ of the best performance ratio. The performance profile in Figure 5 has been plotted on the log scale for clarity.
It is clear from Figure 5 that the provision of a heuristic warm-start solution improves solution time, with the models that make use of a warm-start solution being faster to solve for a greater proportion of instances that their respective models without a warm-start.
It can also be seen that the models that make use of transitivity constraints are slower to solve to optimality for a greater proportion of instances than their respective models that do not use transitivity constraints. However, the inclusion of transitivity constraints does increase the proportion of instances that can be solved to optimality, by 5.3% for the basic model, and by 5.2% for the model with warm-start. Figure 6 plots the cumulative percentage of instances solved to within a given optimality gap within the 20 minute time-limit. Note that the left-hand y-intercept of this figure gives the same information as the right-hand y-intercept in Figure 5, that is, the proportion of instances solved to optimality using each model. Looking at Figure 6, it can be seen that as well as increasing the proportion of instances that can be solved to optimality, the inclusion of transitivity constraints increases the proportion of instances that can be solved to within a given optimality gap. Of the 255 instances for which an optimal solution was unable to be found with any model, but for which a feasible solution was found using all models, the average optimality gap was 24.53% for the basic model, 22.80% with the inclusion of transitivity constraints, 24.71% with the inclusion of a warm-start solution, and 22.36%  In Table 1, we now compare the performance of the basic model (34)-(44) and its strongest extension, with the results of the strongest existing algorithm for the two-stage robust RCPSP, the primal method (Bruni et al., 2018). For each set of test instances, J301, . . . , J3048, Table 1 reports instance parameters (NC, RF, RS), the average CPU time required to solve the instances that were solved to optimality (time), the average optimality gap for those instances which were not solved to optimality but for which a feasible solution was obtained (gap), and the number of instances solved to optimality (#solv).
Note that in four of the most challenging instance sets, J3013, J3021, J3029, J3041, the   Table 1 show that the primal method solves one or two instances to optimality, sometimes outperforming the model proposed in this paper over these instance sets. However, these optimal solutions are obtained whilst simultaneously reaching the maximum time-limit of 1200 seconds, and it is therefore unclear whether or not this is a numerical inaccuracy in the results presented in Bruni et al. (2018).
The results in this table show that the models proposed in this paper solve almost 50% more instances than the primal method, and do so in a considerably shorter computation time. These results confirm the strength of the new model proposed in this paper.

Conclusion
This paper has introduced a new mixed-integer linear programming formulation for the robust counterpart to the two-stage adjustable robust RCPSP. This new compact formulation has been derived by considering a reformulation of the second-stage adversarial sub-problem of maximising the worst-case delayed makespan for a project without resource conflicts.
The reformulation of this sub-problem is equivalent to a longest-path problem over an augmented project network made from multiple copies of the original project network. Hence, the dual of this longest-path problem can be inserted into the first-stage resource allocation problem to obtain a compact minimisation problem for the full two-stage robust RCPSP.
The performance of this new formulation has been examined over 1440 instances of varying characteristics and difficulty, and results show that the proposed formulation can be solved by standard optimisation software significantly faster than the current best algorithm for solving this problem, and can be solved to optimality for almost 50% more instances.
Regarding future research on the two-stage robust RCPSP, the development of heuristic approaches for solving larger and more-challenging instances of this problem would seem to be a natural and worthwhile objective.
where A is a |V | × |E| arc-node incidence matrix, B is a |E| × |V | matrix where B (i,j),i = 1 for each (i, j) ∈ E, and 0 otherwise, I V and I E are identity matrices of dimension |V | and |E| respectively, and e T V is a |V | × 1 vectors of 1's. The rows of this matrix have been grouped according to the constraints that they represent, and similarly, the columns have been grouped by the variables that they represent. Ghouila-Houri (1962) showed that a matrix A is totally unimodular if and only if for every subset of rows R, there exists a partition of R into two disjoint subsets R 1 and R 2 such that i∈R1 a ij − i∈R2 a ij ∈ {−1, 0, 1}, ∀j = 1, . . . , n.
Therefore, finding a subset of rows of matrix C for which this condition cannot hold will prove that C is not totally unimodular.
Consider the constraint matrix of the example shown in Figure  Take R to be the subset of rows consisting of the first row of Group 1, and first two rows of Groups 2 and 3. We will refer to these rows as R1, . . . , R5. To ensure that the sum of column C9 is in {−1, 0, 1}, R2 and R3 must be assigned opposite signs. R4 and R5 must have opposite signs to R2 and R3, respectively to ensure that the sum of columns C5 and C6 are in {−1, 0, 1}. Then, whatever the choice of sign for R1, the sum of column C1 and the sum of column C2 cannot both be in {−1, 0, 1}. Hence, there exists a subset of rows for which the Ghouila-Houri characterisation of total unimodularity does not hold, thus proving that matrix C is not totally unimodular, and that model (13)-(22) is not equivalent to its linear relaxation.