The multi-mode resource investment problem: a benchmark library and a computational study of lower and upper bounds

The multi-mode resource investment problem (MRIP) is the multi-mode extension of the resource investment problem, which is also known under the name resource availability cost problem. It is a project scheduling problem with a given due date as well as precedence and resource constraints. The goal is to find a precedence feasible schedule that minimises the resource costs of the allocated resources. To compute these costs, the maximum resource peak is considered regarding renewable resource types, whereas the sum of allocated nonrenewable resource units is used in the case of nonrenewable resources. Many practical and complex project scheduling settings can be modelled with this type of problem. Especially with the usage of different processing modes, time and cost compromises can be utilised by the project manager. In the literature, some procedures for the MRIP have been investigated; however, the computational experiments in these studies have not been carried out on common benchmark instances. This makes a fair comparison of methods difficult. Therefore, we generated novel instances specifically designed for this problem and published them on the website https://riplib.hsu-hh.de. On this website, the instances as well as best-known solution values are available and researchers can also contribute their findings. We investigate these novel instances by proposing and evaluating lower bounds for the MRIP. Additionally, we analyse the proposed instances by applying heuristic as well as exact methods. These experiments suggest that most of the instances are challenging and further research is needed. Finally, we show some computational complexity properties of the NP-hard MRIP.


Introduction
Project scheduling is an important part of project management. Often, projects have a deadline or a due date and the tasks of the project compete against each other in the use of scarce resources. In the resource investment problem (RIP), the decision maker can determine how many units of each resource are allocated to the project in order to find a schedule that meets the due date. This schedule needs to respect the allocated resources as well as the precedence relations among the activities of the project. The goal is to minimise the project costs that are defined by the allocated units of resources. In this work, we consider two types of resources: renewable and nonrenewable resources. The amount of available renewable resource units is replenished after each time period. So, for the cost calculation of the renewable resource costs, the maximum resource usage peak is multiplied with a cost factor. This cost factor defines the costs of adding an extra resource unit to the project that is available in each period of the project horizon. Nonrenewable resource units are consumed for the entire project horizon and their resource costs are calculated by considering the total sum of consumed resource units.
The resource investment problem was firstly introduced by Möhring (1984) and is also known as the resource availability cost problem (RACP) in the literature. Möhring (1984) describes it as a "problem of scarce time" in contrast to the wellknown and related resource-constrained project scheduling problem (RCPSP) which is a "problem of scarce resources". With the RIP, it is possible to model a wide range of applications such as the construction or dismantling of buildings or software development projects, just to name a few (e.g. Bartels 2009). Several extensions of the problem were introduced in the literature. In this work, we consider the multi-mode resource investment problem (MRIP), which was firstly introduced by Hsu and Kim (2005). Here, each activity can be processed in multiple modes which vary in the resource consumption and the processing duration. For the RCPSP and its multi-mode extension, several benchmark instance libraries, such as the Boctor datasets (Boctor 1993), the PSPLIB (Kolisch and Sprecher 1997), or the MMLIB (Van Peteghem and Vanhoucke 2014), exist. Yet, for the MRIP, no benchmark instances are available to the public which makes the comparison of solution procedures hard. For this reason, we designed a set of benchmark instances and made them available to the public on the website https ://ripli b.hsu-hh.de.
The contributions of this work are the following: • We show that it is sufficient to consider only one nonrenewable resource in the MRIP and provide a transformation for instances with multiple nonrenewable resources. • We describe and compare different approaches to compute lower bounds for the MRIP. • We have set-up and maintain a website https ://ripli b.hsu-hh.de. Here, researchers can access a new set of benchmark instances for the MRIP and share their results with others.

3
The multi-mode resource investment problem: a benchmark library… • We use different exact and heuristic procedures to solve the novel benchmark instances. This enables us to gain insights into which characteristics lead to "easy" and "hard" instances.
This paper is organised as follows: In Sect. 2, we give a formal definition of the problem and show that we can aggregate multiple nonrenewable resources into a single one. Section 3 provides an overview of the existing literature concerning the RIP or extensions of it. Next, we present several lower bound procedures as well as different optimisation techniques to generate upper bounds for the MRIP in Sect. 4. We applied mixed-integer programming (MIP), constraint programming (CP), simulated annealing (SA) and a multi-start local search (MLS). In Sect. 5, we explain which parameters we used for generating benchmark instances and what difficulties we encountered in the process of doing so. Sect. 6 contains a computational study with our findings considering the performance of the procedures presented in Sect. 4. Finally, in Sect. 7, we finish this work with a critical appraisal and an outlook for further research.

Problem statement
An instance of the MRIP is defined by the following properties: a set of activities A = {0, … , n + 1} , a set of precedence relations E ⊆ A × A , a set of renewable resources R and a set of nonrenewable resources R n . For each activity, i ∈ A , a set of modes M i exists and for each mode m ∈ M i , the duration d im ∈ ℤ + 0 is given. Activity 0 and n + 1 are dummy activities that mark the beginning and end of the project and their duration and resource requirements are equal to 0. In this work, we consider only finish-to-start precedence constraints among activities which are represented by the set E. So, if (i, j) ∈ E for activities i, j ∈ A this means that there is a minimum time lag of 0 between the finish of activity i and the start of activity j, i.e. j can only start after i is finished.
The due date D ∈ ℤ + is fixed and defines the maximum project duration. The resource requirement r imk ∈ ℤ + 0 ( r n imk ∈ ℤ + 0 ) of each non-dummy activity i depends on the mode m and the renewable resource k ∈ R (nonrenewable resource k ∈ R n ). For each renewable resource k ∈ R (nonrenewable resource k ∈ R n ) a resource cost factor c k ∈ ℤ + 0 ( c n k ∈ ℤ + 0 ) defines the price of allocating one unit of the resource to the project. For renewable resources, the allocated amount of the resource is replenished after each time period. A schedule is resource feasible with respect to the renewable resources if the amount of allocated resource units is larger than or equal to the maximum resource peak usage. That means that for each period of the project horizon, the sum of resource requirements of activities that are in process in this period has to be smaller than or equal to the allocated amount. To calculate the renewable resource costs, we multiply the allocated amount that corresponds with the maximum resource peak with the resource cost factor. The renewable resource type is useful to model for example workers or machines.
The amount of allocated nonrenewable resources, however, is consumed by the activities over the whole project and they are not replenished. To calculate the nonrenewable resource costs, we sum up the nonrenewable resource requirements of all activities and multiply them with the respective cost factor (i.e. no peak usage is considered as the resource units do not replenish here). Resources of this type can represent budgets or rare materials. The nonrenewable resources are also useful to model outsourcing certain activities to external contractors.
A small example of an MRIP instance is depicted in Fig. 1. Here, we have five activities with activities 0 and 4 being the dummy start and dummy end activity of the project, respectively. The arcs in the network represent the precedence relations among activities. So, e.g. activity 2 can only start after activity 1 is finished. For each activity, there is only one mode available except for activity 3. Fig. 1b depicts the duration of each mode as well as the renewable resource consumption r i,m,1 and the nonrenewable resource consumption r n i,m,1 . Hence, we consider only one renewable and one nonrenewable resource in this example. The unit cost factors are c 1 = 2 and c n 1 = 1 and there is no upper bound on the capacity of each resource. The due date of the project is D = 4 . In Fig. 2a, we show a Gantt chart, where each activity is processed in mode 1. Since there is a renewable resource usage of 3 units in period 1 and 2, we need to allocate 3 units of the renewable resource to the project. Because of the mode choice, 0 units of the nonrenewable resource are utilised, and thus, the schedule has a cost value of 6. The precedence relations as well as the due date are respected. By switching the mode of activity 3 (mode 2 instead of mode 1), we obtain a better solution (Fig. 2b). Here, only 2 units of the renewable resource are allocated to the project since the peak resource  0  1  0  0  0  1  1  2  1  0  2  1  1  2  0  3  1  2  2  0  2  3  1  1  4 1 0 0 0 (b) Duration and resource consumption The multi-mode resource investment problem: a benchmark library… usage is now 2 instead of 3, but also 1 nonrenewable resource unit is needed. Therefore, the cost of the schedule depicted in Fig. 2b is 5. In equations (1)- (7), we display a mathematical model formulation for the MRIP, which is also used in Sect. 4 in a MIP procedure. It is an adaptation of the model of Talbot (1982) for the multi-mode extension of the RCPSP. We use two types of decision variables in the mathematical model: real-valued variables a k and a n k for the resource allocation of resource k ∈ R ∪ R n and binary variables x imt for the mode and scheduling choice. The variable x imt takes a value of 1 if and only if activity i is processed in mode m ∈ M i and starts at period t. For the start period, we compute lower (and upper) bounds called ES i ( LS i ) using forward and backward calculation (FBC, similar to Kelley 1963). Here, we use the minimum duration of each activity and the due date D acts as an upper bound of the latest start time of the project dummy end activity n + 1 (i.e. LS n+1 = D ). We also compute an upper bound LF i on the latest finish period with respect to the due date D using backward calculation, and hence, LF i − d im is the latest possible start of activity i if it is processed in mode m.
In the objective function (1), the resource costs are minimised. Constraints (2) enforce that for each activity i exactly one mode and one start period is chosen. The inequalities (3) represent the precedence relations: if (i, j) ∈ E , then the finish period of activity i (left side of the inequality) has to be lower than or equal to the start period of activity j (right side). Inequalities (4) and (5) make sure that the amount of allocated resources a n k and a k is as least as high as the consumption of the nonrenewable and renewable resources, respectively. Finally, in terms (6) and (7), the two different types of decision variables are depicted. According to , the binary variables are so-called pulse variables over discrete time periods and inequalities (3) are the socalled aggregated precedence constraints. Hence, we call the formulation displayed in (1)-(7) the "aggregated discrete time formulation based on pulse variables" (PDT). In Sect. 4, we will use this formulation and others (using the disaggregated version of the precedence constraints and/or other decision variable types) in a MIP.
Since the resource investment problem with a single mode per activity is NP-hard (cf. Möhring 1984), the multi-mode extension is also NP-hard. This is still true if we extend the problem setting with nonrenewable resources. However, it is sufficient to consider only a single nonrenewable resource since we can aggregate multiple nonrenewable resources into a single one. We compute the novel resource requirement of an activity i and a mode m to be the sum of the old resource requirements times the respective resource cost factor: The cost factor of this single nonrenewable resource is equal to 1. Note that this is possible since there is no upper bound on the resource use and the resource consumption is not time dependent for nonrenewable resources. For this aggregated nonrenewable resource, we can set the allocation to be a n 1 = ∑ k∈R n c n k ⋅ a n k based on allocations a n k of the former resources. Obviously, the objective value does not change by this aggregation and also all constraints such as precedence relations and renewable resource constraints are not changed. This transformation is the reason why we consider only instances with a single nonrenewable resource in this study and it can be used to convert instances with multiple nonrenewable resources into the single resource case.

Literature review
The RIP is closely related to other project scheduling problems such as the RCPSP or the resource levelling problem (RLP). However, the goal of the RCPSP and its multi-mode extension (MRCPSP) is the minimisation of the makespan with fixed resource availabilities. Several heuristic and exact procedures have been proposed for the MRCPSP (e.g. Geiger 2017 andSchnell andHartl 2016, respectively). In the RLP, also a due date for the latest project completion is given, yet the objective function often differs. Several resource levelling objective functions are known in the literature such as the total squared utilisation cost or the total overload cost (Rieck and Zimmermann 2015). Bianco et al. (2016) extended the RLP setting with generalised precedence relations and variable intensities in the execution of the activities. Another closely related problem to the RIP is the time-constrained project scheduling problem (TCPSP). It can be seen as a combination of the RIP and the RCPSP since there is a given due date as well as resource capacities. However, additional

3
The multi-mode resource investment problem: a benchmark library… capacities can be temporarily allocated to the project at a certain cost. So, it has to be decided in which periods extra resource units are added. The goal is to minimise the cost for additional resources under the given due date. This problem was firstly proposed by Deckro and Hebert (1989), and lately, Verbeeck et al. (2017) proposed an artificial immune system (AIS) implementation for the TCPSP. Next, we give a brief literature overview on the existing work for the RIP and some of its extensions. Several exact solution methods for the RIP exist. Möhring (1984) was the first to use the RIP to model a bridge construction project. He proposed an exact method using graph theoretical algorithms to solve the problem. Another exact algorithm, called minimum bounding algorithm (MBA), was introduced by Demeulemeester (1995), where a branch-and-bound procedure for the RCPSP is applied iteratively. Improving this procedure, Yamashita (2010, 2015) proposed the modified minimum bounding algorithm by utilising a feasible initial solution that is found heuristically. Lately, Kreter et al. (2018) studied the RIP as well as an extension with general temporal constraints and one with calendar constraints. They provided mixed-integer linear programming formulations as well as constraint programming (CP) implementations for the problems. With their CP procedure, they were able to close all available instances from the literature for these single-mode problems. The authors also give an overview on previous work on the RIP with generalised precedence constraints (RIP/max). On their website, the authors provide instances for the single-mode RIP as well as the RIP/max with up to 500 activities per project. Drexl and Kimms (2001) studied the computation of lower bounds (LB) for the RIP. They proposed two procedures: one using Lagrangian relaxation and the other based on column generation techniques. Both procedures also yield feasible solutions for the problem as a by-product. The authors conducted computational experiments on project instances with up to 30 activities and up to 8 resources.
The initial application of a metaheuristic for the RIP was proposed by Yamashita et al. (2006). They implemented a scatter search (SS) which is a population-based metaheuristic and it outperformed two simple multi-start heuristics as well as the upper bounds obtained by Drexl and Kimms (2001) on instances with up to 120 activities. Ranjbar et al. (2008) also implemented population-based metaheuristics: path relinking (PR) and genetic algorithm (GA). Here, the PR achieved slightly better results than the GA on the instances generated in Yamashita et al. (2006), but due to different hardware they did not compare their results directly to the SS mentioned above. An implementation of the AIS metaheuristic as well as a benchmark set with 30 activities and 4 resources (called RACP30) is presented by Van Peteghem and Vanhoucke (2013). They show that their approach outperforms the GA for the tardiness permitted extension presented by Shadrokh and Kianfar (2007) and the proposed instances are also used in Van Peteghem and Vanhoucke (2015). Meng et al. (2016) proposed a hybrid metaheuristic by combining the tabu search (TS) metaheuristic with the SS metaheuristic. They tested their procedure, called tabued scatter search (TSS), on 48 adapted RCPSP instances from the PSPLIB. Experiments showed that on average TSS is able to obtain better solutions than the SS of Yamashita et al. (2006), but by spending a higher computational time. A novel heuristic approach called multi-start iterative search heuristic (MSIS) is proposed by Zhu et al. (2017). Here, the authors combine local search techniques with path relinking and show that their method outperforms the PR, the GA and the SS procedures presented earlier. The computational experiments were performed on adapted PSPLIB instances and newly generated ones. Shadrokh and Kianfar (2007) began studying the extension of the RIP, where exceeding the due date (tardiness) is permitted but penalised in the objective function. The penalty costs rise for each time period that the project finishes later than its specified due date by a fixed tardiness cost factor. It is called resource investment problem with tardiness (RIPT) and the authors applied a GA to tackle this novel problem. Van Peteghem and Vanhoucke (2015) apply a metaheuristic procedure called invasive weed optimisation algorithm (IWO) to the RIP as well as the RIPT. Computational experiments performed on the RACP30 instance set as well as on 30 activity projects adapted from the PSPLIB indicate that the IWO outperforms both the AIS from Van Peteghem and Vanhoucke (2013) (2007).
The multi-mode extension of the RIP was firstly studied Hsu and Kim (2005). They proposed a heuristic procedure that combines two priority rules: one regarding the increase in costs when adding an activity to a partial schedule and the other considering how the finish time of the current activity affects its successors' remaining start times. They tested different weight combinations of this combined priority rule heuristic against heuristics that applied different priority rules sequentially (one to select the next activity and another one to select the mode and start time). For the experiments, they used MRCPSP instances from the PSPLIB with 12 to 30 activities, but treated the nonrenewable resources as renewable resources. Another heuristic procedure for the MRIP was presented by Qi et al. (2015). The authors proposed a novel schedule generation scheme as well as the modified particle swarm optimisation (MPSO) metaheuristic, which is a combination of particle swarm optimisation (PSO) and SS. They also adapted PSPLIB instances for the MRCPSP with 12-30 activities to test their MPSO heuristic against a PSO implementation and an adaptation of the GA of Ranjbar et al. (2008). Colak and Azizoglu (2014) proposed a heuristic procedure for the special case of the MRIP with a single renewable resource. Their approach uses different construction procedures and tries to improve a solution using several neighbourhood search strategies. The authors performed experiments on instances with up to 100 activities and up to 10 modes per activity.
For the MRIP, two exact procedures are known. Yamashita and Morabito (2009) combined the MBA of Demeulemeester (1995) with an exact branch-and-bound algorithm for the MRCPSP of Sprecher and Drexl (1998). In order to compute time/ cost trade-off curves, they investigated several different due dates per instance. Since the procedure relies on solving multiple MRCPSP exactly, they solved only small instances with 15 activities per project. The other exact procedure was presented by Coughlan et al. (2015). The authors also added calendar constraints to the MRIP setting, making some resources not available at certain times. They proposed a Dantzig-Wolfe reformulation in combination with a column generation and a branch-and-price algorithm, where they heavily utilise the calendar structure of the resource availability. In computational experiments, Coughlan et al. (2015) compared their approach to a standard MIP implementation in CPLEX 12.2. Their approach was able to close several 50 activity instances, while a normal MIP implementation in CPLEX often failed to solve the instance.
Another extension of the MRIP with generalised temporal constraints and socalled cumulative resources was proposed by Bartels (2009). Cumulative resources are similar to renewable resources but are used to model storage. The start and finish of an activity can change the resource level in a positive or negative way. Bartels proposed two application cases where one considers the dismantling of nuclear power plants (RBP). There, at most two modes are available for each activity but only the resource usage can vary and not the activity processing time. The author used projects with 20 − 60 activities and 6 renewable resources in the computational study. The objective is to minimise the net present value associated with modedependent costs (which is similar to the nonrenewable resource in the MRIP setting but time dependent). The second area of application models is the scheduling of testing in the automotive industry (VTP). Here, tests (modelled as activities) are assigned to experimental vehicles (represented by cumulative resources) through different modes and the goal is to minimise the total number of utilised experimental vehicles. Again, the processing time does not vary with the mode choice as only resource utilisation is affected by the mode. In Zimmermann (2009, 2015), this problem type is further analysed and a priority rule-based schedule generation scheme as well as a GA are utilised. The authors applied these procedures to instances with 20 and 600 activities. Each instance incorporates one renewable resource to model the construction of the experimental vehicles and a cumulative resource for each experimental vehicle.
The problem extension with both a tardiness penalty and multiple modes (MRIPT) was studied by Gerhards and Stürck (2018). They proposed a hybrid large neighbourhood search (LNS) procedure that uses MIP techniques to solve subproblems exactly. They carried out experiments on adapted 30 activity MRCPSP instances of the PSPLIB.
Several other extensions of the resource investment problem exist, where the authors adapted the objective function. In the work of Najafi and Niaki (2006), a RIP extension with discounted cash flows and net present value maximisation is proposed and a GA was applied to the problem.  extend this setting further by adding generalised precedence constraints to the problem and also applied a GA. A RIP extension with net present value maximisation and tardiness penalty was studied by Najafi and Azimi (2009) and a priority rule-based heuristic was proposed. Yamashita et al. (2007) added uncertainty to the activity durations in the RIP. The authors applied a SS with PR to different scenarios to obtain a robust solution.
In Table 1, we give an overview of the instance properties used in multi-mode RIP studies. Here, |M| is the maximum number of modes per activity. It is clear that most of the studies only address small projects with rather few activities. The majority of the instances are based on MRCPSP instances and only the work of Gerhards and Stürck (2018) used nonrenewable resources (when nonrenewable resources were available in the original MRCPSP instance, the other authors treated them as renewable resources). Hence, it is desirable to have a more challenging and publicly available benchmark set of instances for future studies.

Lower and upper bounds
Next, we present several procedures to obtain lower (Sect. 4.1) and upper bounds (Sect. 4.2) for the MRIP. We apply these methods in Sect. 6 to evaluate the novel set of benchmarks instances.

Lower bounds
In order to obtain lower bounds, we use four formulae as well as the linear programming relaxations of different mathematical formulations and the destructive improvement approach. Two rather simple lower bounds for the minimum resource level for the renewable resources in the RIP presented by Drexl and Kimms (2001) can be adapted for the multi-mode setting ( a 1 k and a 2 k ). We calculate how many units for a resource k ∈ R are needed so that every activity can be performed in its least resource consuming mode (with r min ik = min m∈M i r imk we denote the minimum resource demand of activity i for resource k). Hence, the minimum resource level for the first method is as follows.
As a second way to compute a lower bound on the minimal consumption, we can distribute the required resources equally over the planning horizon. Hence, we divide the sum of the minimal products of the resource consumption and the duration of all activities by the due date D. 1 3 The multi-mode resource investment problem: a benchmark library… Next, we utilise so-called core times to compute bounds on the minimal consumption (cf. Klein and Scholl 1999). The core time of an activity i is the periods where the activity has to be processed, i.e. after its latest start LS i and before it earliest finish EF i = ES i + min{d im } . So, this concept uses the precedence relations as well as the due date as an upper bound on the project makespan. However, if the due date is large or the minimum activity duration is small, this core time interval can be empty. The bound is computed by checking each period and summing up the minimal resource consumptions of activities that have a core time in this period.
As a fourth way to compute minimal resource consumptions, we want to identify triplets of activities that cannot be scheduled sequentially, and therefore, at least two of them have to be processed with some overlap. We look at triplets of activities i, j, l ∈ A with no precedence relations among them. There are six potential order- we check whether at least one of them is feasible with respect to the due date. If not, at least two of the three activities have to be carried out in intersecting time intervals and we can add the minimal resource usages. For example, let us assume we investigate the ordering i ≺ j ≺ l , i.e. activity i has to be finished before j can start and j has to be finished before l can start. This means that we can compute a new earliest start time ES j = max{ES j , ES i + d min i } and latest start time LS j = min{LS j , LS l − d min j } . The ordering is not feasible if LS j < ES j . If all six potential orderings are infeasible with respect to the due date, then we add the triplet to the set of infeasible triplets IT. For each triplet in IT, at least two activities have to be scheduled with overlap and we get the following lower bound for the resource usage: For a nonrenewable resource, the only simple lower bound R n k is the sum of the minimal consumptions of all activities.
By multiplying these minimal resource levels with the corresponding resource cost factors and summing them, we get the following simple lower bound for the MRIP.
We compare LB 0 with bounds obtained by solving the linear programming (LP) relaxation of the MIP. Therefore, we use the mathematical formulation displayed in (1)-(7) (PDT) and we refer to the objective value of this LP relaxation as LB 1 . Furthermore, we also use a mathematical formulation where the precedence relations are modelled in a disaggregated way. The constraints are displayed in (15). Here, for each pair (i, j) ∈ E and for each time period t of the planning horizon a constraint is added that enforces the precedence relations. When the right side of (15) is equal to 1, i.e. the successor j starts before or in period t, then it forces the left side to take a value of 1 as well. Hence, the predecessor i has to be started before t − d im depending on the mode m.
Artigues (2017) reported that this disaggregated formulation is stronger (w.r.t to the LP relaxation) than the aggregated formulation displayed in constraints (3) which can be seen since constraints (3) are implied by (2) and (15). The "disaggregated discrete time formulation based on pulse variables" (PDDT) is defined by equations (1), (2), (4)- (7) and (15). With LB 2 , we refer to the lower bound obtained by solving the LP relaxation of the PDDT formulation. Note that depending on the value of D, the PDDT formulation can contain considerable more constraints than the PDT formulation, and thus, setting up the mathematical model in the solver most likely takes a longer time.
Another method for computing lower bounds is the destructive improvement strategy introduced by Klein and Scholl (1999) for the RCPSP. It is an iterative approach that is started with LB 0 as a starting lower bound B . In each iteration, we try to proof that if we take B as an upper bound on the objective value, then no feasible solution can exist. If we succeed with the proof, we know that B + 1 is a valid lower bound for the instance and we can use this value in the next iteration. This process is repeated until the proof of infeasibility fails. For the proof, we either use the PDT formulation and a MIP solver ( LB 3 ) or the CP formulation (displayed in "Appendix 1.3") and the IBM ILOG CP Optimizer solver ( LB 4 ). In each iteration, we add an extra constraint that bounds the objective value by the current upper bound B to the respective problem formulation. To limit the overall run-time of the procedure, we allow the solver a run-time of 60 s for each iteration. If no infeasibility can be detected after that time, the procedure stops and returns the best-known lower bound. However, if the respective solver finds a feasible solution in an iteration, this solution has to be optimal and the procedure terminates as well.

Upper bounds
To obtain upper bounds, i.e. the costs of feasible solutions, we use different approaches: On the one hand, we implement heuristic procedures such as a multi-start local search (MLS), a simulated annealing (SA) procedure and an adaptation of the priority rule heuristic (PRH) of Hsu and Kim (2005). On the other hand, we also apply The multi-mode resource investment problem: a benchmark library… exact methods such as MIP and CP solvers that are able to detect optimality for some instances (although with a run-time restriction this is not always the case). Before we explain how the MLS and the SA work in particular, we explain the schedule generation scheme (SGS) that is utilised by both methods. In our case, we use a serial SGS that takes a scheduling sequence and a mode as an input and tries to build a feasible solution (i.e. a schedule with start and finish times and resource allocations). The SGS that we apply works similar to Algorithm 3.7.1 introduced by Neumann et al. (2003) but also has some slight differences. We schedule the activities one at a time and in the order specified by the sequence (in Neumann et al. 2003 the SGS schedules critical activities first and then selects the next activity according to a priority rule; one could interpret our scheduling sequence as a special kind of priority rule). So, the SGS ends after n + 2 iterations and we expect the scheduling sequence to respect the precedence constraints. In iteration l of the SGS, it tries to schedule activity i = l in mode m = i at the least cost increasing feasible time period that respects the precedence relations. Thus, the SGS checks all feasible start times t ∈ {ES i , … , LF i − d im } and computes the cost increase in the partial schedule with the activity added at each particular start time. If there are multiple start times with the same cost increment, the earliest one of them is assigned. It is possible that the SGS may not find a feasible start time with respect to the due date since the durations of the chosen modes in are too high. If the SGS detects such an infeasibility, it fails and does not return a feasible schedule.
First, we explain the concept of the multi-start local search method that we used. Algorithm 1 depicts the outline of the MLS. We initialise the current best mode vector b by choosing the mode with minimal duration for each activity (ties are resolved arbitrarily). This initial mode vector is feasible with respect to the due date D. Next, we determine an initial scheduling sequence b that respects the precedence relations. This means, that a successor j cannot occur at an earlier position in the sequence b than its predecessor i if (i, j) ∈ E . We generate such a sequence by adding the activities one at a time to the sequence. An activity can only be added to the sequence if it is not part of the sequence yet and if all of its predecessors are already in the sequence. If more than one activity is eligible, we either choose one arbitrarily among the eligible activities (random generation) or choose the activity with the highest value of some priority rule. The priority rules that we used in this process are the following: minimum LST, minimum LFT, minimum slack (i.e. LFT − EFT ), minimum number of successors, maximum number of successors and maximum rank positional weight (cf. Kolisch 1996). We generate a sequence for each priority rule and one using the random selection technique and use the SGS and the mode vector b to translate it into a schedule. We keep the scheduling sequence with the lowest cost value and store it as b .
In the main loop of Algorithm 1, we start by assigning a perturbed mode vector, called p in line 5. Here, we select a random number of activities and change their mode from the one stored in b to an arbitrary one. For the remaining activities (that were not selected), we assign in p the same mode as in b . The current best sequence b is also modified (by extracting a random number of activities from the sequence and inserting them at a random, precedence feasible position in the modified sequence) and then copied to p . We use these two steps to obtain a new start point for the local search. In order to limit the degree of diversification from the current best solution p and p , we limit the number of changes by n ⋅ max M and n ⋅ max S , respectively. Here, max M , max S ∈ [0, 1] are parameters of the MLS. Then, we perform a local search starting at p and p in line 7. This local search procedure aims to improve the solution costs by altering p and p locally. First, the LS procedure checks if altering the chosen mode of an activity leads to lower costs by applying the SGS to the altered mode vector and the unchanged sequence. After checking all possible mode changes of an activity, the LS explores all feasible (with respect to precedence relations) swaps in the scheduling sequence. So, it exchanges the activities at two positions of the scheduling sequence and if the sequence still respects the precedence relations, it checks if the application of the SGS leads to a lower cost value. If an alteration of p or p leads to an improvement, we apply the change and repeat the local search procedure with the updated mode vector and scheduling sequence. Hence, we apply a first improving move policy. In lines 8 − 12 , we update the current best mode vector and scheduling sequence if the LS found a strictly better cost value. The MLS repeats until the specified time limit is reached.
We briefly explain the priority rule heuristic of Hsu and Kim (2005): First, we initialise the investment upper bound IUB = LB 1 . Then, the PRH iteratively tries to find a feasible schedule with costs lower than or equal to IUB . If we cannot find such a schedule, we increase IUB by one at the end of the iteration. To obtain a schedule, we schedule one activity at a time. Therefore, we compute for each activity i that is eligible for scheduling (i.e. all of its predecessors are already scheduled) the mode m * and the start time t * that results in the lowest priority value v(i, m * , t * ) . Then, the activity i * with the worst minimum value is 1 3 The multi-mode resource investment problem: a benchmark library… selected for scheduling. The priority value v is a linear combination of two priority functions v 1 and v 2 . The parameter controls the portion of the two priority functions.
Here, v 1 is called the transformed slack priority function and measures how the selected start time and mode influence the set of remaining start times for the successors of this activity (it favours early finish times such that many possible start times remain for the successors). Priority function v 2 , called transformed investment priority function, calculates how much it costs to schedule an activity in a specific mode and a start time with respect to the increase in costs of the partial solution (in contrast to the original work, we also consider nonrenewable resources and their costs). If the costs excel the given bound IUB , the v 2 value is ∞ . The scheduling procedure stops, if either a start time for each activity was determined or the costs of the (partial) solution excelled the given upper bound IUB . The whole priority rule heuristic stops once a feasible solution is found.
As a second metaheuristic approach, we adapted the simulated annealing procedure with reheating presented by Józefowska et al. (2001) for the MRCPSP to fit our MRIP setting. Here, in each iteration, we generate a solution candidate by perturbing the currently best-known solution. It is accepted if it has a better cost value or with some probability that depends on the delta of the cost values as well as the so-called temperature. The temperature is used to control the acceptance rate of worse solutions during the search and decreases constantly by a cooling factor . However, we also use reheating to escape local optima when the search is stuck. The total number of reheats is one of the algorithm parameters of the SA approach as well as the cooling factor .
Next, we explain the set-up of the MIP experiments. We implemented six formulations: the two based on pulse variables PDT (displayed in (1)- (7)) and PDDT (displayed in (1), (2), (4)- (7) and (15)). Furthermore, we adapted two formulations based on step variables (SDT and SDDT) and two based on on/off variables (OODT and OODDT) which are displayed in more detail in "Appendix 1.1" and 1.2, respectively. In order to model and solve the MIP, we used Gurobi Optimizer 9.0 via the C# API. For large instances, setting up the model can require some time (in our experiments, the maximal measured set-up time was 8.16 s), while for the smaller instances, no big difference in set-up times between the two formulations was measured.
Another exact approach that becomes more and more relevant in the field of project scheduling is constraint programming. Constraint programming solvers became more efficient in recent years and they are applied to a growing number of scheduling problems (e.g. Kreter et al. 2018 andSchnell andHartl 2016). The CP formulation of the MRIP used in this study is displayed in "Appendix 1.3". We implemented the CP model displayed in (47)-(54).

Instance generation and benchmark library
In this section, we explain how we generated the benchmark instances for the MRIP. Our main concerns with the instances are diversity of the instances with respect to certain characteristics and that the feasible mode space cannot be reduced easily. We used the following characteristics to generate different instances: number of activities n, maximum number of modes per activity |M|, number of renewable resources |R|, due date factor , order strength OS and resource factor RF. Here, the parameter order strength is a measure of the number of precedence relations (Mastor 1970). It is defined as the fraction of precedence relations in E compared to the total number of possible relations and, hence, is an indicator whether the precedence structure of the project is more parallel or more serial (Schwindt 1998 showed that OS is a good approximation for the restrictiveness of the precedence relations). The resource factor value is the average portion of different resources required for the processing of an activity in a mode (Kolisch and Sprecher 1997). Note that it only varies for the renewable resources. For the single nonrenewable resource, RF = 1 for all instances since each mode has a positive resource requirement for the nonrenewable resource except the modes of the dummy start and finish activities. The due date factor is used to compute the due date of the project based on the earliest start time ES n+1 of the dummy end activity (calculated by forward calculation). We calculate the due date as follows: The round function in equation (17) applies the rounding to the nearest even number strategy in the case of midpoint values. Table 2 displays the values that we used. For each parameter combination, we generated 5 instances, giving us in total 4950 instances. As shown in Sect. 2, it is sufficient to consider only one nonrenewable resource.
Concerning the mode space, the instances should neither contain inefficient nor infeasible modes. We call a mode m ∈ M i of activity i inefficient if there is another mode m � ∈ M i for that activity such that the other mode has a shorter or equal duration and it has lower or equal resource requirements for all resources (cf. Kolisch et al. 1995). It would never be beneficial to include an inefficient mode m into a solution since with m ′ the resource usage would be lower, and hence, we can omit mode m from the mode space. It is possible to check if a mode is inefficient in polynomial 1 3 The multi-mode resource investment problem: a benchmark library… time. We call a mode m ∈ M i infeasible if its duration d im is too long to finish the project before the due date. In mathematical terms, this happens if the difference between the latest finish time (LF) of the activity and the earliest start time (ES) is smaller than the duration, i.e.: Especially when using project data that was generated for the MRCPSP and when the due date is very close to ES n+1 , many modes are infeasible and can be omitted from the mode space. Again, we can check in polynomial time if a mode is infeasible since we can calculate the ES and the LF in polynomial time with the FBC. Gerhards and Stürck (2018) adapted MRCPSP instances from the PSPLIB with 30 activities per project. They also computed the due date as in (17) and used values of ∈ {1.0, 1.1, … , 1.5} . However, they did not investigate if the duration of each mode is short enough for the resulting due dates. In Table 3, we present the average number of infeasible modes for the different values of as well as the maximum and minimum proportion of infeasible modes for single instances (labelled min and max in the table). Especially for = 1.0 and = 1.1 , many modes are infeasible. Furthermore, we also computed how many modes become inefficient when we try to make the infeasible modes feasible. We altered the duration of all infeasible modes by setting d im = LF i − ES i . Again, for low values of many modes can be omitted from the mode space since they are inefficient (the altered modes have a shorter duration and dominate unchanged ones). For this reason, we consider in this study slightly larger values for the parameter and check during the generation of the instances for infeasible and inefficient modes.
In the following, we describe how we computed each instance with the desired properties. First, we used the network generator RanGen Demeulemeester et al. (2003) to generate an activity-on-the-node network with the specified number of activities and OS value (this gives us the set E of precedence relations). For each activity i ∈ A , the cardinality of the mode set is set to the desired number, i.e. |M i | = |M| (except for the dummy activities which have only one mode). We draw for each activity i and all its modes m ∈ M i the duration d im as a discrete uniformly distributed random number U{1, 10} . Similarly, the resource requirements r imk ( r n imk ) for each resource k ∈ R ( k ∈ R n ) are also taken from U{1, 10} . If the value of RF < 1 , then we arbitrarily set sufficiently many of the renewable resource requirements of each mode to 0 such that the resource factor is achieved. After we determined all the resource requirements and the durations of an activity, we check if there are inefficient modes. If an inefficient mode occurs, we take new random values for the modes that cause the inefficiency. This is repeated until each activity has no inefficient modes. Then, we use the forward calculation to compute ES and determine the due date D as in (17). Based on the due date, we can compute latest finish times LF for all activities. Next, we check each activity for infeasible modes. If we encounter an infeasible mode with some activity, we start over and redraw all values (durations and resource requirements) of all activity modes of the instance. By taking new random values for all of the modes and not just the infeasible ones, we want to avoid that infeasible modes are set to a similar duration (for example, setting the duration to d im = LF i − ES i would repair the infeasibility, but could also result in inefficient modes and biased duration values). We repeat this until all modes of all activities are feasible and none of them is dominated. Finally, we determine the random cost factors c k for all renewable resources k ∈ R from U{1, 10} . Although, we only investigate the MRIP without tardiness permitted in this study, we also added a value for the tardiness cost factor (also from U{1, 10} ) that can be useful in future work. In the next section, we will investigate the novel instances. For most instances (3725 of 4950), we drew feasible mode durations on the first try. However, on average we had to redraw 385 times per instance until all modes were feasible. The maximum number of redraws was 67,653 for instance MRIP_30_79. We analysed if the redrawing has an effect on the distribution of the mode duration values. Therefore, we performed a Pearson's chi-squared test of goodness of fit (Cochran 1952) with significance level = 0.01 to check whether the duration data fits to a discrete uniform distribution. Only for 47 of the 4950 instances, we were able to reject this hypothesis, and hence, we suspect that there is no strong bias in the duration values for most of the generated instances.

Computational study
To test the "hardness" of the instances, we compute both lower and upper bounds with the procedures presented in Sect. 4 and compare them. To compute upper bounds, we apply metaheuristic search procedures as well as different MIP and CP implementations. The goal is to examine how hard the proposed instances are for these general purpose methods. For the lower bounds, we propose some rather simple approaches and compare them to linear programming relaxations of different mathematical formulations as well as so-called destructive improvement methods. We want to investigate if drawing the cost factors from a uniform distribution has an effect on the "hardness" of the instances. Therefore, we perform experiments with the random cost factor instances and additionally with the same instances, but with all cost factors set to 1 (called equal resource cost factors). We conducted all following computational experiments on an Intel Xeon Silver 4214 CPU running at 2.20 GHz and all procedures were implemented in C#. We utilised Gurobi Optimizer 9.0 to solve LP relaxations and the MIPs. For the CP-based destructive improvement 1 3 The multi-mode resource investment problem: a benchmark library… procedure and the standalone CP implementation, we used IBM ILOG CP Optimizer 12.9.0 (cf. Laborie et al. 2018).
The five lower bound procedures proposed in Sect. 4.1 are applied on the new benchmark instances. In Table 4, we depict the average of the relative improvement over the worst lower bound LB min for all five lower bounds, i.e. LB i −LB min LB min , i = 0, … , 4 . If the value of ∅ improvement is close to 0% , then the respective lower bound value was close to the minimum value among the calculated bounds. Hence, the higher the ∅ improvement the better. We used random resource cost factors for the results depicted in Table 4. For instances with unit resource cost factors, we get a similar trend.
Clearly, LB 0 (the bound based on minimum resource level formulae) is always the worst lower bound. Only for instances with a small number of modes ( |M| = 3 ) or a high resource factor ( RF = 1 ), the gap between the more advanced procedures gets smaller. The LP relaxation-based lower bound with disaggregated precedence constraints LB 2 performs best on average. Since the disaggregated formulation is stronger than the aggregated one (cf. Artigues 2017), LB 1 (based on the aggregated precedence constraints) performs slightly worse than LB 2 regarding the relative improvement. The destructive improvement methods ( LB 3 and LB 4 ) work especially well on small-to medium-sized instances but are outperformed by the LP-based approaches on the n = 100 instances. However, LB 3 finds optimal solution for 31 instances, while LB 4 can solve 4 instances to optimality. In absolute terms, LB 3 gives the best lower bound for 3736 instances (and 3403 of those bounds were solely found by this procedure). Yet, the performance of the MIP solver seems to be significantly worse for the instances with 100 activities, as we can see a drop in the average improvement there. Furthermore, in Table 5, we display the average computation time (in s). As expected, LB 0 has a significantly lower running time than the other procedures. The PDT LP relaxation ( LB 1 ) is solved about 10 times faster than the disaggregated version. While the MIP-based destructive improvement procedure takes the longest computation time, it also closes most instances and performs best on the instances with 30 and 50 activities. The CP-based variant is not able to find as many optimal solutions as its MIP-based counterpart, yet, it terminates faster. For instances with 30 and 50 activities, LB 2 seems to have the best quality per computation time ratio, while on the larger 100 activities instances, LB 1 seems most efficient. However, we solved the LP relaxations ( LB 1 and LB 2 ) much faster when the mode count was higher. To investigate this unexpected behaviour, we modified the instances

3
The multi-mode resource investment problem: a benchmark library… with three modes per activity by duplicating each mode of each activity and solving the modified instance again as an LP relaxation. Note that the objective value does not change when solving this modified instance. However, the average computation time decreased from 18.44 to 15.96 s for the aggregated variant and from 234.38 to 163.49 s for the disaggregated version. So, basically the same instance was solved up to 30% faster just by adding the identical modes. Why the LP solver is so much faster with these additional decision variables is still an open problem and further research is needed here. Next, we analyse the upper bound procedures presented in Sect. 4.2. Therefore, let us first examine how many instances the MIP and CP procedures solved to optimality within a time limit of one hour. When we used equal resource cost factors, the PDT formulation MIP solved 17.4% of the instances to optimality, whereas the CP solver found optimal solutions for 30.4% of the instances. Using random cost factors, PDT solved only 9.3% and CP 28.4% of the instances to optimality. In total, the exact procedures solved 1430 instances to optimality, and hence, we found optimal solutions for less than 29% of the instances with random cost factors. The random cost factors make the instances more challenging, especially for the MIP-based approaches. Therefore, we take a closer look at these instances.
In Table 6, we can observe how the exact approaches performed given different run-time limits (60, 600, and 3600 s). Surprisingly, the CP implementation finds more optimal solutions in 60 s than any of the MIP approaches in one hour. The aggregated versions of the step-based and the on/off-based formulations cannot find feasible solutions for all of the instances, even with a time limit of one hour. The disaggregated variants find at least one feasible solution for each instance except for PDDT and OODDT with the 60 second time limit. PDT can solve the most instances to optimality among the MIP-based approaches. However, the CP approach seems to be much better at proving optimality of a solution and, hence, achieves also faster average run-times.   In Table 7, we can observe how many instances were solved optimally by CP or one of the MIP formulations after one hour of computation depending on the respective instance parameter. We see that for "small" instances, i.e. the instances with only 30 activities, CP outperforms all other approaches and is the only procedure that finds optimal solutions for instances with 100 activities. In general, instances with 30 activities, a small RF = 0.25 or a small due date factor are more likely to be solved to optimality by one of the approaches. A similar phenomenon as with the LP relaxations occurs with the PDDT MIP approach where more optimal solutions are found for instances with higher mode count as the MIP solver has to solve several LPs during the search.
Next, we compare the results of the exact procedures to the MLS, the SA approach and the PRH. To calibrate the algorithm parameters of the MLS, PRH and SA, we used the software package López-Ibáñez et al. (2016) and a set of 990 training instances (with the same instance parameters). After performing 1000 experiments with each algorithm, the best parameter values on the training instance set are max S = 0.92 , max M = 0.2 and = 0.28 for the MLS and the PRH. For the SA approach, the number of reheats was 23, 880 and 5280 and the temperature cooling factor was 0.52, 0.9 and 0.9 for a run-time limit of 60, 600 and 3600 s, respectively.

3
The multi-mode resource investment problem: a benchmark library… We see that using mostly the transformed investment priority function performs better for PRH (in Hsu and Kim 2005, an value of 0.4 performed best). On average, the PRH took 12.5 s to compute a feasible solution. The fastest PRH running time took 0.08 s, while the maximum computation time in our experiments was 258.5 s for one instance.
Metaheuristic procedures are in general not able to detect if a solution is optimal or not. So, we compare the procedures based on the average relative deviation ( ∅ RD, cf. (19)) from the respective best-known solution value ( UB min ) found by the four procedures. So, if this average is 0% , this means that the procedure always found a solution with lowest objective function value for each instance.
We also compare them by the average relative deviation from the best-known lower bound ( ∅RDLB , cf. (20)).
Tables 8 and 9 depict the average relative deviations ∅RD and ∅RDLB of the tested procedures, respectively. Note, that SDT and OODT did not find a feasible solution for each instance and we omitted instances with no solution in the calculation of the average for these two formulations. This means that the true average for those procedures may be higher than displayed and a comparison with the other methods is difficult. We observe that the PDT formulation reaches the lowest average deviation over all instances although CP was able to find more optimal solutions. One could suspect that the stronger PDDT formulation should also have a better performance than the weaker aggregated PDT variant. But we need to keep in mind, that solving the disaggregated LPs takes longer and the extra constraints also need more memory, and hence, it is not a priori clear which MIP formulation performs better (cf. Artigues 2017). However, for the project instances with 30 activities, the CP implementation achieves the best results. For instances with more activities, again PDT achieves the best results and also the gap to the metaheuristic procedures grows. Among the heuristic procedures, MLS works best for small-to medium-sized project instances but SA performs better on the 100 activity instances. However, the heuristic approaches are not competitive with the exact methods and there is a lot of potential for improvement.
In Table 9, we see that PDDT, SDDT and OODDT have a lower deviation for instances with a higher number of modes. We think that this phenomenon is related to the one experienced with the lower bounds. When solving the MIP, also several LPs are solved, and hence, we think this is related to the better performance of the disaggregated MIP procedures on instances with 6 modes per activity. A low resource factor of 0.25 seems to make it easier for the CP and PDT formulation which is in line with the findings in Table 7 where the most optimally solved instances also had a RF = 0.25 . The choice of the due date factor has no strong impact on the solution quality of the CP solver. The MIP approaches, especially SDDT and OODDT, perform worse for 1 3 The multi-mode resource investment problem: a benchmark library… higher due date factors. OODDT performs best over all procedures on the = 1.2 and worst on the = 2 instances. In total, most best-known solutions (BKS) were found by CP (2228) followed by OODDT (1362), PDDT (1268) and PDT (1239). The MLS found only 15 BKS and the PRH only a single one, while the SA approach did not find any BKS at all.

Summary and conclusions
In this paper, we investigated the multi-mode resource investment problem. It is a prominent project scheduling problem where each activity is processed in one of multiple modes that vary in the activities' duration and resource consumption. We have shown that it is sufficient to consider a single nonrenewable resource with a transformation procedure. Furthermore, we propose a novel set of benchmark instances for this problem as no common set of instances was used and known in the literature so far. In total, we computed 4950 instances with a diverse set of instance characteristics. We especially assured that none of the modes can be reduced because of infeasibility or inefficiency reasons. By maintaining the website https ://ripli b.hsu-hh.de, we make these instances available to the public. In addition, researchers can access best-known solution values for the instances as well as share and validate their results on this website. We encourage researchers to test their solution procedures on the benchmark dataset and compare their results.
In extensive computational experiments, we examined lower and upper bounds for the MRIP. For the lower bounds, our experiments revealed that using the LP relaxation of the so-called disaggregated discrete time indexed formulation yields better lower bounds for most instances at hand. However, we also proposed destructive improvement methods that yielded good results for the small-and medium-sized instances and even provided optimal solutions in some cases.
We also tested several procedures to obtain good upper bounds for the MRIP. The metaheuristic procedures, a multi-start local search, simulated annealing, and a priority rule heuristic from the literature, were not able to compete with the MIP and CP implementations. In total, the exact procedures were able to prove the optimality of 1340 of the 4950 instances. That means that over 60% of the instances are still open and we encourage researchers to investigate them further. Our experiments also indicated that the instances are more challenging when we use random cost factors.
For future research, we advise the application of more advanced metaheuristic procedures. In addition, further extensions such as general temporal constraints or the tardiness penalty in the objective function could be an interesting addition to take some important aspects of project scheduling into account.
The multi-mode resource investment problem: a benchmark library… are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/. This formulation is called the aggregated version of the discrete time formulation based on step variables (SDT). The objective function (22) is the same as in the PDT formulation and minimises the resource costs. The equation in (23) ensures that exactly one mode and one step (i.e. a start) is chosen for each activity. The inequalities in (24) let the step variables grow, while the equality in (26) makes sure that the start time of each activity is scheduled before the mode specific latest start time LF i − d im . In (27), we display the aggregated version of the precedence constraints. The next inequalities (28) model the nonrenewable resources. (29) displays the renewable resource constraints. Finally, in (30) and (31), the domains of the decision variables are defined. There is also a version with disaggregated precedence constraints which are displayed in (32).

3
The multi-mode resource investment problem: a benchmark library… By combining (21) and (33), we get the following formula for the start time of activity i The model using on/off variables and aggregated precedence constraints (OODT) is displayed in (35)-(45). Since we are in a multi-mode setting, we need to introduce an additional binary variable e im to represent the mode choice. This decision variable is equal to 1 if and only if activity i is processed in mode m.