# New exact and heuristic algorithms to solve the prize-collecting job sequencing problem with one common and multiple secondary resources

Aurélien Froger, Ruslan Sadykov

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Aurélien Froger, Ruslan Sadykov. New exact and heuristic algorithms to solve the prize-collecting job sequencing problem with one common and multiple secondary resources. European Journal of Operational Research, In press, 10.1016/j.ejor.2022.07.012. hal-03287769v3

HAL Id: hal-03287769

https://inria.hal.science/hal-03287769v3

Submitted on 9 Sep 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# New exact and heuristic algorithms to solve the prize-collecting job sequencing problem with one common and multiple secondary resources

Aurélien Froger¹,²

Ruslan Sadykov²

¹ Université de Bordeaux, CNRS, Bordeaux INP, IMB, UMR 5251, 33405 Talence, France

² Centre Inria de l'université de Bordeaux, 33405 Talence, France

## Abstract

We study the prize-collecting job sequencing problem with one common and multiple secondary resources. In this problem, a set of jobs is given, each with a profit, multiple time windows for its execution, and a duration during which it requires the main resource. Each job also requires a preassigned secondary resource before, during, and after its use of the main resource. The goal is to select and schedule the subset of jobs that maximize the total profit. We present a new mixed integer linear programming formulation of the problem and a branch-cut-and-price algorithm as an exact solution method. We also introduce a heuristic algorithm to tackle larger instances. Extensive numerical experiments show that our exact algorithm can solve to optimality literature instances with up to 500 jobs for a particular dataset and up to 250 jobs for another dataset with different characteristics. Our heuristic builds high-quality solutions in a small computational time. It computes new best-known solutions for most of the larger instances.

**Keywords** — scheduling; sequencing; labeling algorithm; branch-cut-and-price; iterated local search

## 1 Introduction

This work focuses on solving the prize-collecting job sequencing problem with one common and multiple secondary resources (PC-JSOCMSR). This problem was introduced by Horn et al. (2018). It has application in particle therapy scheduling (Maschler and Raidl, 2020) and in pre RUNTIME scheduling of avionic systems. We refer to the work Horn et al. (2021a) for details about these real-world applications (see §3.3 of the previously cited work).

The PC-JSOCMSR is defined as follows. We are given a set $\mathcal{J} = \llbracket 1, n \rrbracket$ of non-preemptive jobs, meaning that the execution of a job cannot be interrupted and resumed later (hereafter, $\llbracket a_1, a_2 \rrbracket = [a_1, a_2] \cap \mathbb{Z}$ for every $a_1, a_2 \in \mathbb{Z}$ such that $a_1 \le a_2$). Each job $j \in \mathcal{J}$ should be entirely scheduled (from its start to its end) within one of the time windows of set $W_j = \{T_{jk} | k \in \llbracket 1, w_j \rrbracket\}$ with $T_{jk} = [t_{jk}^1, t_{jk}^2]$. Time windows are indexed in chronological order of the $t_{jk}^1$ values. Denoting $t^1 = \min_{j \in \mathcal{J}} \{t_{j1}^1\}$ and $t^2 = \max_{j \in \mathcal{J}} \{t_{jw_j}^2\}$, the planning horizon is defined as the interval $[t^1, t^2]$.

We are also given a set of (renewable) disjunctive resources $\mathcal{R}_0 = \{0\} \cup \mathcal{R}$. We refer to resource 0 as the main resource and to $\mathcal{R} = \llbracket 1, m \rrbracket$ as the set of secondary resources. Each job $j \in \mathcal{J}$ requires a single secondary resource $r_j \in \mathcal{R}$ (preassigned to j) for $p_j$ consecutive time periods. If the execution of job j starts at time

Email addresses: aurelien.froger@math.u-bordeaux.fr (Aurélien Froger), ruslan.sadykov@inria.fr (Ruslan Sadykov)
<table><caption>(a) Instance</caption><thead><tr><th rowspan="2">Job j</th><th rowspan="2">r<sub>j</sub></th><th rowspan="2">b<sub>j</sub></th><th colspan="3">p<sub>j</sub></th><th colspan="4">W<sub>j</sub></th></tr><tr><th>p&#773;<sub>j</sub></th><th>p<sup>0</sup><sub>j</sub></th><th>p<sup>+</sup><sub>j</sub></th><th>t<sup>1</sup><sub>j1</sub></th><th>t<sup>2</sup><sub>j1</sub></th><th>t<sup>1</sup><sub>j2</sub></th><th>t<sup>2</sup><sub>j2</sub></th></tr></thead><tbody><tr><td>J<sub>1</sub></td><td>2</td><td>4</td><td>2</td><td>2</td><td>1</td><td>0</td><td>9</td><td>10</td><td>15</td></tr><tr><td>J<sub>2</sub></td><td>1</td><td>1</td><td>3</td><td>3</td><td>2</td><td>3</td><td>12</td><td>-</td><td>-</td></tr><tr><td>J<sub>3</sub></td><td>2</td><td>2</td><td>1</td><td>1</td><td>1</td><td>1</td><td>12</td><td>-</td><td>-</td></tr><tr><td>J<sub>4</sub></td><td>1</td><td>2</td><td>2</td><td>3</td><td>2</td><td>4</td><td>14</td><td>-</td><td>-</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>(b) Solution</figcaption></figure>

<caption>Figure 1: An example of an instance and a feasible solution</caption>

period $t$, then this job requires the main resource (in addition to the resource $r_j$) from time period $t + p_j^-$ for $p_j^0$ consecutive time periods (with $p_j^- + p_j^0 \le p_j$). We denote $\hat{p}_r^- = \max_{j \in \mathcal{J}: r_j=r} \{p_j^-\}$ and $\hat{p}_r^+ = \max_{j \in \mathcal{J}: r_j=r} \{p_j^+\}$ where $p_j^+ = p_j - (p_j^- + p_j^0)$. We have $p_j^- \ge 0$ and $p_j^+ \ge 0$ for every job $j$. The order of processing starting times on every secondary resource is a suborder of the order of processing starting times on the main resource.

Scheduling a job $j \in \mathcal{J}$ generates a profit $b_j \ge 0$, whereas not scheduling it generates no profit. The problem is to schedule the jobs while maximizing the total profit. A feasible solution to the problem is a schedule that satisfies the following constraints: 1) each job is scheduled at most once, 2) if a job is scheduled, it is scheduled uninterruptedly and entirely (from its start to its end) within one of its time windows, and 3) at any time, every resource (main or secondary) is used by at most one job. Hereafter, we refer to constraints 1), 2), and 3) respectively as *elementarity*, *time window*, and *resource usage* constraints. Figure 1 displays a feasible solution to an instance with 4 jobs $J_1, J_2, J_3$, and $J_4$. On the right, the vertical axis lists the resources and the horizontal axis represents time. Each rectangle represents the use of a particular resource by a job during a time interval. In the displayed solution, job $J_2$ is not processed. The starting time of jobs $J_1, J_3$, and $J_4$ (i.e., the time at which the preassigned secondary resource starts being used) is equal to 0, 8, and 4 respectively. As $p_{J_4}^- = 2$, job $J_4$ starts using the main resource at time $6(= 4 + 2)$ for $p_{J_4}^0 = 3$ time periods. The same is true for the other jobs. The solution value is equal to $b_{J_1} + b_{J_2} + b_{J_3} = 8$.

We now outline the main original contributions of this paper:

* We introduce dominance rules that define properties satisfied by at least one optimal schedule.

* We present a pseudo-polynomial time-indexed mixed integer linear programming (MILP) model for the PC-JSOCMSR. Although indexing variables on time may limit a priori its interest, its numerical performance on existing benchmark instances is largely better than an existing compact order-based MILP model introduced by Horn et al. (2021b), and a significant number of instances are solved to proven optimality in less than two hours.

* We formulate the PC-JSOCMSR as a resource constrained shortest path problem (RCSPP) and describe labeling algorithms to solve it. We discuss different modeling approaches. We present how we use the dominance rules.

* We present a branch-cut-and-price (BCP) algorithm to solve an extended MILP model we introduce for the PC-JSOCMSR. The elementarity constraints are handled in the master problem, which allows to efficiently solve the pricing problem as a RCSPP. Our BCP includes many ingredients found in state-of-the-art BCP algorithms: limited arc memory rank-1 cuts, automatic dual price smoothing stabilization, a procedure for enumerating elementary paths, and a multiphase strong branching procedure.

* We present a heuristic algorithm to initialize the BCP algorithm and to compute high-quality solutions for large-sized instances that remain beyond the reach of an exact algorithm.

Extensive computational experiments on benchmark instances from the literature show that the new exact approaches we present to tackle the PC-JSOCMSR solve to proven optimality instances with up to 250 jobs and
up to 500 jobs for some specific instances. They also show that our heuristic algorithm computes high-quality
solutions within short running time and outperforms existing heuristics.

The paper is organized as follows. Section 2 presents a brief overview of related literature. Section 3 presents two dominance rules. Section 4 presents two MILP models for the PC-JSOCMSR, one of which is new. Section 5 presents how we model the problem as a RCSPP. Section 6 describes the BCP algorithm. Section 7 presents the heuristic algorithm. Section 8 shows the results of our computational experiments. Section 9 contains final remarks.

## 2 Related literature

The PC-JSOCMSR is NP-hard (Horn et al., 2018). It is a recently introduced problem that has been little studied. Particular attention has been paid to the use of multivalued decision diagrams (MDDs). In this context, an MDD is a directed weighted acyclic multigraph (usually layered) storing solutions of the PC-JSOCMSR and their values as paths between a starting node and an ending node. An MDD is said relaxed (resp. restricted) when it encodes a superset (resp. subset) of all feasible solutions. In the top-down compilation procedure, an MDD is built using a dynamic programming formulation of the PC-JSOCMSR with nodes to merge identified using an attached dynamic programming state. If additional nodes are merged to limit the size of the graph, the procedure builds a relaxed MDD. Removing some nodes during the procedure leads to a restricted MDD. In the incremental refinement procedure, a relaxed MDD is first built with one node per layer and is then incrementally refined by removing arcs and splitting nodes to eliminate infeasible solutions. We refer to (Castro et al., 2022) for a survey on the use of MDDs in discrete optimization. Maschler and Raidl (2018) presented top-down compilation (TDC) and incremental refinement procedures to build relaxed multivalued decision diagrams (MDDs) used to obtain dual bounds. To compute heuristic solutions, they built a restricted MDD with a TDC procedure and introduced a general variable neighborhood search (GVNS). The former method dominated the latter in providing quality solutions for instances with less than 100 jobs, while the opposite was more often true for instances with 100 to 300 jobs. To our knowledge, the state-of-the art algorithm for the PC-JSOCMSR is presented in Horn et al. (2021a). It relies on relaxed and restricted MDDs build with a TDC procedure that uses a A* search principle and cost-based filtering procedures to remove arcs in the MDDs. Some instances with 50 jobs are solved to optimality within the 15 minutes time limit. For larger instances, the gap between the value of the best computed solution and the best dual bound is on average between 8 and 15%.

As mentioned by previous works on the PC-JSOCMSR, a related problem is the single-machine scheduling problem introduced by van der Veen et al. (1998), in which a change-over time between the processing of two consecutive tasks must be respected. The value of this change-over time depends on whether the two jobs belongs to the same group. Contrary to this problem where the change-over time is only related to the previously scheduled job, the change-over time incurred before starting a job *j* in the PC-JSOCMSR depends on two jobs (when different): the last job scheduled on the main resource and the last job scheduled on the secondary resource *r*<sub>j</sub>. There exist additional differences with the PC-JSOCMSR: no time windows and an objective function focused on minimizing the total sequencing time.

More generally, the PC-JSOCMSR shares similarities with single-machine scheduling problems with sequence-
dependent setup times. Indeed, a solution to the PC-JSOCMSR can be represented as a sequence of jobs using
one after another the main resource. As pointed out before, the setup time prior to running a job is dependent
not only on the last scheduled job on the main resource, but also on another job scheduled before. To a cer-
tain extent, the problem can be classified as a scheduling problem with past-sequence-dependent setup times,
concept introduced by Koulamas and Kyparisis (2008). To our knowledge, the literature on this problem is
nonetheless essentially focused on potential learning or deterioration effects related to the past jobs. We refer
readers to the survey of Allahverdi (2015). Moreover, another major difference with such scheduling problems
is the presence of several time windows in the PC-JSOCMSR.

As mentioned by Horn et al. (2021a), the PC-JSOCMSR is a resource-constrained project scheduling problem (RCPSP) with fixed time lags between some of the tasks and the objective of minimizing a weighted sum of tardy
jobs. To understand this, each job $j \in \mathcal{J}$ is split into three tasks with the adequate precedence constraints. The three tasks require the secondary resource $r_j$ and have duration $p_j^-, p_j^0$, and $p_j^+$. The second task also requires the main resource. The weight associated with these tasks are 0 for the first and second task, and $b_j$ for the third one. While making this connection explicit is interesting, it does not seem likely that adapting methods introduced for RCPSPs would be suitable for the PC-JSOCMSR due to the existence of multiple time windows for the tasks in the latter.

Finally, as mentioned by Horn et al. (2021a), the PC-JSOCMSR shares similarities with orienteering problems with time windows. The two problems share the same objective function and time windows should be respected for the start time of the service of a customer in orienteering problems and for the starting time of a job in the PC-JSOCMSR. The two main differences are that there is a travel time between two successive visits to clients in orienteering problems, and that the waiting time before scheduling a task on the main resource in the PC-JSOCMSR depends not only on its time window, but also on the completion time of a job that is not the last one scheduled on the main resource. Despite these differences, taking inspiration from solution methods to orienteering problems with time windows (Vansteenwegen et al., 2011; Gunawan et al., 2016) is of interest for solving the PC-JSOCMSR.

## 3 Dominance rules

With the aim of reducing the search space, we introduce two dominance rules, directly deduced from the problem definition, that describe properties satisfied by at least one optimal schedule to the PC-JSOCMSR. Specifically, in Proposition 1 and 2 we provide sufficient conditions for restricting the order in which two jobs use one after the other the main resource. These conditions guarantee that a better feasible schedule always exists by reversing the order of the two jobs (Proposition 1), by inserting another job in between or by replacing one of the two jobs (Proposition 2). For each proposition, its precise mathematical description and its proof are given in §S.1 of the supplementary material. Solution methods of the PC-JSOCMSR can take advantage of these results.

**Proposition 1.** Let $j, j' \in \mathcal{J}$ such that $j \neq j' \land r_j = r_{j'}, k \in [1, w_j]$, and $k' \in [1, w_{j'}]$. The following conditions are sufficient to ensure that there exists an optimal schedule where the main resource is not used by job $j'$ during time window $T_{j'k'}$ and then used consecutively by job $j$ during time window $T_{jk}$:

$$t_{jk}^1 \le t_{j'k'}^1 \land t_{jk}^2 \le t_{j'k'}^2 \land p_j^- \ge p_{j'}^- \land p_j^+ \le p_{j'}^+$$

**Proposition 2.** Assume that a job $j' \in \mathcal{J}$ is scheduled during time window $T_{j'k'}$ with $k' \in [1, w_{j'}]$. There exists an optimal schedule in which:

1. Job $j'$ does not use the main resource just after its use by another job $j'' \in \mathcal{J}\{j'\}$ during time window $T_{j''k''}$ with $k'' \in [1, w_{j''}]$ (resp. $j'$ is not the first job using the main resource) if there exists a job $j \in \mathcal{J}\{j', j''\}$ (resp. $j \in \mathcal{J}\{j'\}$) that cannot use the main resource before $j''$ and after $j'$ (resp. after $j'$) and that satisfies at least one of the two following conditions:

a) Job $j$ uses the same secondary resource as $j'$, is at least as profitable as $j'$, and can always be scheduled in the place of $j'$ regardless of when this latter is scheduled during $T_{j'k'}$.

b) Job $j$ can use the main resource i) just after $j''$ without preventing $j''$ to end its execution at time $t_{j''k''}^2$ and imposing restriction on the processing of precedent jobs, and ii) just before $j'$ without preventing $j'$ to start its execution at time $t_{j'k'}^1$ and imposing restriction on the processing of subsequent jobs (resp. ii))

2. Job $j'$ does not use the main resource just before its use by another job $j'' \in \mathcal{J}\{j'\}$ during time window $T_{j''k''}$ with $k'' \in [1, w_{j''}]$ (resp. $j'$ is not the last job using the main resource) if there exists a job $j \in \mathcal{J}\{j', j''\}$ (resp. $j \in \mathcal{J}\{j'\}$) that cannot use the main resource before $j'$ and after $j''$ (resp. before $j'$) and that satisfies at least one of the two following conditions:
a) Job j uses the same secondary resource as j', is at least as profitable as j', and can always be scheduled in the place of j' regardless of when this latter is scheduled during $T_{j'k'}$.

b) Job j can use the main resource i) just after j' without preventing j' to end its execution at time $t_{j'k'}^2$ and imposing restriction on the processing of precedent jobs, and ii) just before j'' without preventing j'' to start its execution at time $t_{j''k''}^1$ and imposing restriction on the processing of subsequent jobs (resp. i))

# 4 Mixed integer linear programming models

The PC-JSOCMSR can be formulated using mixed linear integer programming. We present two models. The first formulation is a compact MILP model introduced in Horn et al. (2021b). The second formulation is a pseudo-polynomial MILP model that relies on time discretization.

## 4.1 Compact order-based MILP model of Horn et al. (2021b)

The compact model introduced by Horn et al. (2021b) is an ordered-based model. It relies on the following decision variables. Binary variable $y_j$ is equal to 1 if job $j \in \mathcal{J}$ is scheduled, 0 otherwise. Binary variable $z_{jk}$ is equal to 1 if job $j \in \mathcal{J}$ is scheduled within its time window $T_{jk}$. Continuous variable $s_j$ is the starting time period of job $j$. Binary variable $u_{jj'}$ is equal to 1 if job $j$ finishes using the main resource before its use (not necessarily consecutive) by job $j'$. Note that when $u_{jj'} = 1$, job $j'$ is not necessarily the next immediate job using the main resource after $j$. We introduce $\delta_{jj'}$ as the minimum time required between the start of job $j$ and the beginning of another job $j'$ to avoid that these jobs overlap on the main and secondary resources when they are scheduled in this order, but not necessarily consecutively. The value of $\delta_{jj'}$ is equal to $p_j$ if $j$ and $j'$ use both the same secondary resource (i.e., $r_j = r_{j'}$) and to $p_j - p_{j'}$ if they do not (i.e., $r_j \neq r_{j'}$). The model, denoted [F1], is as follows:

$$ [F1] \quad \max \sum_{j \in \mathcal{J}} b_j y_j \qquad (1) $$

$$ \text{s.t.} \quad \sum_{k=1}^{w_j} z_{jk} = y_j \qquad j \in \mathcal{J} \qquad (2) $$

$$ u_{jj'} + u_{j'j} \ge y_j + y_{j'} - 1 \qquad j, j' \in \mathcal{J}, j \ne j' \qquad (3) $$

$$ s_{j'} \ge s_j + \delta_{jj'} - (t_{jw_j}^2 - p_j - t_{j'1}^1 + \delta_{jj'}) (1 - u_{jj'}) \qquad j, j' \in \mathcal{J}, j \ne j' \qquad (4) $$

$$ s_j \ge t_{j1}^1 + \sum_{k=1}^{w_j} (t_{jk}^1 - t_{j1}^1) z_{jk} \qquad j \in \mathcal{J} \qquad (5) $$

$$ s_j \le t_{jw_j}^2 - p_j + \sum_{k=1}^{w_j} (t_{jk}^2 - t_{jw_j}^2) z_{jk} \qquad j \in \mathcal{J} \qquad (6) $$

$$ y_j \in \{0, 1\} \qquad j \in \mathcal{J} \qquad (7) $$

$$ z_{jk} \in \{0, 1\} \qquad j \in \mathcal{J}, k \in \{1, \dots, w_j\} \qquad (8) $$

$$ s_j \in [t_{j1}^1, t_{jw_j}^2 - p_j] \qquad j \in \mathcal{J} \qquad (9) $$

$$ u_{jj'} \in \{0, 1\} \qquad j, j' \in \mathcal{J}, j \ne j' \qquad (10) $$

The objective function (1) maximizes the total profit for the scheduled jobs. Constraints (2) state that a job is deemed scheduled if it is scheduled within one of its time window. Constraints (3) impose a linear order for the use of the main resource for every couple of jobs that are scheduled. Constraints (4) ensure that jobs are not using the same resources simultaneously. Specifically, according to the order in which they use the main resource, the constraints guarantee that the starting times of two scheduled jobs is sufficient to avoid overlapping. Constraints (5) and (6) impose that jobs are scheduled within their time windows. Finally, constraints (7)–(10) set the domain of the decision variables.
## 4.2 Pseudo-polynomial time-indexed MILP model

Model [F1] is a continuous-time model as the start time of every jobs is a continuous variable that can take an infinite number of values provided that the constraints of the problem are respected. We introduce a pseudo-polynomial formulation that relies on time discretization. If the data is not integer, provided that the granularity of the time discretization is fine enough, a discrete-time formulation has the same optimal value as a continuous-time formulation of the same problem. Available benchmark instances introduced in previous studies use integer data (see Section 8).

The model, denoted [F2], uses the following decision variables. As in [F1], binary variable $y_j$ is equal to 1 if job $j \in \mathcal{J}$ is scheduled, 0 otherwise. We denote $\mathcal{T}_j$ as the set of possible starting time periods for job $j \in \mathcal{J}$. Specifically we have $\mathcal{T}_j = \bigcup_{k \in \{1, \dots, w_j\}} \{t_{jk}^1, t_{jk}^1 + 1, \dots, t_{jk}^2 - p_j\}$. We denote as $\mathcal{T}$ the set of time periods (i.e., $\mathcal{T} = \{t^1, t^1 + 1, \dots, t^2\}$). Binary variable $x_{jt}$ is equal to 1 if job $j \in \mathcal{J}$ starts being executed at time period $t \in \mathcal{T}_j$, 0 otherwise. The model is as follows:

$$ [F2] \quad \max \sum_{j \in \mathcal{J}} b_j y_j \qquad (11) $$

$$ \text{s.t.} \quad \sum_{t \in \mathcal{T}_j} x_{jt} = y_j \qquad j \in \mathcal{J} \qquad (12) $$

$$ \sum_{j \in \mathcal{J}} \sum_{t' \in \mathcal{T}_j \cap \{t-p_j^- - p_j^0 + 1, \dots, t-p_j^-\}} x_{jt'} \le 1 \qquad t \in \mathcal{T} \qquad (13) $$

$$ \sum_{\substack{j \in \mathcal{J} \\ \wedge r_j = r}} \sum_{t' \in \mathcal{T}_j \cap \{t-p_j^0 + 1, \dots, t\}} x_{jt'} \le 1 \qquad t \in \mathcal{T}, r \in \mathcal{R} \qquad (14) $$

$$ y_j \in \{0, 1\} \qquad j \in \mathcal{J} \qquad (15) $$

$$ x_{jt} \in \{0, 1\} \qquad j \in \mathcal{J}, t \in \mathcal{T}_j \qquad (16) $$

The objective function (11) is identical to the one of model [F1] and maximizes the total profit for the scheduled jobs. Constraints (12) state that a job is deemed scheduled if its start is scheduled at any time period that is adequate with respect to its time windows. Constraints (13) and (14) ensure that jobs are not using the same resources simultaneously. Specifically, the number of jobs using the main resource and each secondary resource at any time $t \in \mathcal{T}$ is less than or equal to one. Constraints (15) and (16) set the domains of the decision variables.

# 5 Graph-based model

Previous studies on the PC-JSOCMSR using MDDs introduce very large acyclic graphs where every arc represents the processing of a job with a weight equal to its profit and every path between a starting node and an ending node represents a solution to the PC-JSOCMSR (Maschler and Raidl, 2018; Horn et al., 2021a). An optimal solution is then computed by solving a shortest path problem. If the MDD is relaxed, this solution may be infeasible by potentially violating elementarity, time window, and resource usage constraints.

In this section, we reformulate the PC-JSOCMSR as a resource-constrained shortest path problem (RCSPP) on a small graph that contains cycles. Every arc represents the processing of a job or a transition between the use of the main resource by two jobs. It has a cost and a vector of resource consumption. A path between the source and the sink of this graph represents a feasible solution to the PC-JSOCMSR (i.e., a schedule) if and only if at each node of this path the bounds on the resource consumption accumulated along the path are not violated. An optimal schedule for the PC-JSOCMSR can then be obtained by solving the RCSPP on this small graph by the labeling algorithm presented below. We start with the graph model which describes all feasible schedules, before adapting it to the case when only semi-active schedules are admitted.

We formulate the PC-JSOCMSR as an RCSPP on a directed graph $G = (V, A)$ defined from a set $V$ of nodes and a set $A$ of arcs between these nodes. There are three special nodes $v_{source}$, $v_*$ and $v_{sink}$ in $V$ that are used to model, respectively, the beginning of a schedule, any intermediate moment between the use of the
main resource by two tasks, and the end of a schedule. The set $V$ also contains a node $v_{j,k}$ for every job $j$ and time window $k \in \{1, \dots, w_j\}$. We denote by $V_j = \{v_{j,k} : k \in \{1, \dots, w_j\}\}$ the subset of nodes associated with job $j$. The set of arcs $A$ contains the arcs $(v_{source}, v_\star)$ and $(v_\star, v_{sink})$, as well as arcs $(v_\star, v_{j,k})$ and $(v_{j,k}, v_\star)$ for every pair $(j,k)$. Traversing the arc $(v_{j,k}, v_\star)$ in a path corresponds to processing job $j$ during its time window $T_{jk}$ in the corresponding schedule. Each path from $v_{source}$ to $v_{sink}$ in $G$ defines the jobs to be scheduled and the order in which these jobs use the main resource. There also exists an idle (loop) arc $(v_\star, v_\star)$ to model idles times that are required due to the presence of time windows and sometimes necessary when the earliest starting time of a job inserted at the end of a partial schedule is different when we calculate it independently for the main resource and its secondary resource. We denote as $\Gamma^-(v)$ and $\Gamma^+(v)$ the set of arcs entering in and leaving every node $v$, respectively. A path starting at $v_{source}$ is called a forward path, while a path starting at $v_{sink}$ and composed of arcs traversed in their opposite direction is called a backward path. Any forward (resp. backward) path defines a partial schedule that is built from time $t^1$ (resp. built backward from time $t^2$).

To ensure elementarity, time window and resource usage constraints, we associate a resource consumption to each arc, these consumed resources being in one-to-one correspondence with the resources of set $\mathcal{R}_0$. The accumulated consumption of resource 0 represents time. The accumulated consumption of resource $r \in \mathcal{R}$ represents the time to be added to the current accumulated consumption of resource 0 to obtain the time when $r$ becomes available for jobs. Let $Q_0$ be the accumulated consumption of resource 0 (i.e., the current time) and $Q_r$ be the accumulated consumption of a secondary resource $r$. In a forward (resp. backward) path, this means that the resource 0 is available during the time interval $[Q_0, t^2)$ (resp. $[t^1, Q_0)$), and that resource $r$ is available during the time interval $[Q_0 + Q_r, t^2)$ (resp. $[t^1, Q_0 + Q_r)$). In a forward (resp. backward) path, a negative (resp. positive) value of the accumulated consumption of a secondary resource means that the resource is already available for some time, whereas a positive (resp. negative) value means that the resource is currently used. The consumption $q_{a,r}$ of resource $r$ when traversing arc $a \in A$ is defined as follows:

$$ q_{a,r} = \begin{cases} 0 & \text{if } a = (v_{source}, v_\star) \lor a \in \Gamma^+(v_\star) \setminus (v_\star, v_\star) \\ p_j^0 & \text{if } a = (v_{j,k}, v_\star) \text{ and } r = 0 \\ -p_j^0 & \text{if } a = (v_{j,k}, v_\star) \text{ and } r \in \mathcal{R} \setminus \{r_j\} \\ p_j^- + p_j^+ & \text{if } a = (v_{j,k}, v_\star) \text{ and } r = r_j \\ 1 & \text{if } a = (v_\star, v_\star) \text{ and } r = 0 \\ -1 & \text{if } a = (v_\star, v_\star) \text{ and } r \in \mathcal{R} \end{cases} $$

As a path defines the order in which jobs use the main resource, we increase the accumulated consumption of resource 0 (i.e., we increase the current time) by $p_j^0$ when scheduling job $j$, which requires subtracting $p_j^0$ from the accumulated consumption of the secondary resources different from $r_j$ because the time intervals during which they are continuously available are not modified. If job $j$ starts using the main resource at the current time, it means that it started using $r_j$ for $p_j^-$ time units. After $j$ stops using the main resource, it still uses resource $r_j$ for $p_j^+$ time units, which explains the arc consumption for this resource. As traversing the idle arc $(v_\star, v_\star)$ means going to next time period, we increase the current time by one, which requires subtracting one for the accumulated consumption of every secondary resource (because the time intervals during which they are continuously available are not modified).

Each node $v \in V$ is given an interval $[l_{v,r}, u_{v,r}]$ for the accumulated consumption of each resource $r \in \mathcal{R}_0$. The intervals are defined as follows:

$$ [l_{v,r}, u_{v,r}] = \begin{cases} [t^1, t^1] & \text{if } v = v_{source} \text{ and } r = 0 \\ [0, 0] & \text{if } v = v_{source} \text{ and } r \in \mathcal{R} \\ [t^1, t^2] & \text{if } v = v_\star \text{ and } r = 0 \\ [-\hat{p}_r^-, \hat{p}_r^+] & \text{if } v = v_\star \text{ and } r \in \mathcal{R} \\ [t^2, t^2] & \text{if } v = v_{sink} \text{ and } r = 0 \\ [0, 0] & \text{if } v = v_{sink} \text{ and } r \in \mathcal{R} \\ [t_{jk}^1 + p_j^-, t_{jk}^2 - (p_j^0 + p_j^+)] & \text{if } v = v_{j,k} \text{ and } r = 0 \\ [-p_j^-, -p_j^-] & \text{if } v = v_{j,k} \text{ and } r = r_j \\ [-\hat{p}_r^-, \hat{p}_r^+] & \text{if } v = v_{j,k} \text{ and } r \in \mathcal{R} \setminus \{r_j\} \end{cases} $$
The node $v$ can be traversed by a forward path (resp. by a backward path) only if the resource consumption accumulated up to $v$ belongs to $[l_{v,0}, u_{v,0}]$ for resource 0 and is less than or equal to $u_{v,r}$ (resp. larger than or equal to $l_{v,r}$) for the secondary resources. For a node $v_{j,k}$, the interval $[t_{jk}^1 + p_j^-, t_{jk}^2 - (p_j^0 + p_j^+)]$ for resource 0 guarantees that job $j$ is executed during time window $T_{jk}$. The interval $[-p_j^-, -p_j^-]$ for secondary resource $r_j$ makes sure that this resource is available for long enough when job $j$ starts using the main resource. Indeed, job $j$ can start using resource 0 in a forward (resp. backward) path at the current time $Q_0$ only if the resource $r_j$ is available for a time interval that contains $[Q_0 - p_j^-, t^2)$ (resp. $[t^1, Q_0 - p_j^-)$). The intervals $[-\hat{p}_r^-, \hat{p}_r^+]$ for the secondary resources different from $r_j$ are not restrictive in the sense they do not impact the respect of past and future accumulated consumption intervals. Specifically, in a forward (resp. backward) path, every secondary resource $r$ is available in time interval $[Q_0 + \hat{p}_r^+, t^2)$ (resp. $[t^1, Q_0 - \hat{p}_r^+)$) if the current time is $Q_0$.

The intervals for the accumulated consumption of secondary resources also have the following characteristic which is essential for building a feasible schedule. An accumulated consumption in a forward (resp. backward) path that is smaller than $l_{v,r}$ (resp. larger than $u_{v,r}$) should be updated to $l_{v,r}$ (resp. $u_{v,r}$). This is necessary in the following cases. Let $Q_0$ be the accumulated consumption of resource 0 (i.e., the current time) and $Q_r$ be the accumulated consumption of every secondary resource $r \in \mathcal{R}$. If $Q_0 + Q_{r_j} < Q_0 - p_j^-$ in a forward path, then job $j$ can start using resource 0 at time $Q_0$, which means it starts using resource $r_j$ at time $Q_0 - p_j^-$. An increase of $Q_{r_j}$ to $-p_j^-$ is therefore necessary to account for the starting time of job $j$ and to ensure that the accumulated consumption of resource $r_j$ becomes equal to $p_j^+$ after traversing the arc $(v_{j,k}, v_\star)$. Similarly, if $Q_0 + Q_{r_j} > Q_0 - p_j^-$ in a backward path, then job $j$ can start using resource 0 at time $Q_0$ (time is considered in backward direction), which means that $j$ starts using resource $r_j$ at time $Q_0 - p_j^-$. A decrease of $Q_{r_j}$ to $-p_j^-$ is therefore necessary to account for the fact that $r_j$ is then available only during the time interval $[t^1, Q_0 - p_j^-)$. If, in a forward (resp. backward) path, a resource $r$ different from $r_j$ is available at time $Q_0$ during a time interval that contains $[Q_0 - \hat{p}_r^-, t^2)$ (resp. $[t^1, Q_0 + \hat{p}_r^+)$), meaning that $Q_r < -\hat{p}_r^-$ (resp. $Q_r > \hat{p}_r^+$), then increasing (resp. decreasing) $Q_r$ to $-\hat{p}_r^-$ (resp. $\hat{p}_r^+$) is valid because no job that could start using resource 0 at time $Q_0$ would require the use of $r$ before time $Q_0 - \hat{p}_r^-$ (resp. after time $Q_0 + \hat{p}_r^+$).

We finally define the cost $c_a$ of each arc $a \in A$ as:

$$c_a = \begin{cases} b_j & \text{if } a \in \Gamma^+(v_{j,k}) \\ 0 & \text{otherwise} \end{cases}$$

<figure><img src="image_4.png" /><figcaption>Figure 2: Graph G associated with the instance of Figure 1</figcaption></figure>
To illustrate the graph formulation of the problem, we show in Figure 2 the graph G associated with the instance of Figure 1. The solution of Figure 1 is represented in this graph by the following path: $(v_{source}, v_\star)$, $(v_\star, v_\star)$, $(v_\star, v_\star)$, $(v_\star, v_{11})$, $(v_{11}, v_\star)$, $(v_\star, v_\star)$, $(v_\star, v_\star)$, $(v_\star, v_{41})$, $(v_{41}, v_\star)$, $(v_\star, v_{31})$, $(v_{31}, v_\star)$, $(v_\star, v_\star)$, $(v_\star, v_\star)$, $(v_\star, v_\star)$, $(v_\star, v_\star)$, $(v_\star, v_\star)$, $(v_\star, v_{sink})$. This path can be computed by the labeling algorithms introduced in §5.1.

The PC-JSOCMSR is then equivalent to an RCSPP whose goal is to find the path from $v_{source}$ and $v_{sink}$ in G with maximum cost that satisfies all the resource consumption interval and elementarity constraints. The RCSPP with elementarity constraints is NP-hard for the case of a general graph with cycles (Dror, 1994). We describe in §5.1 the labeling algorithms used to solve this RCSPP. We present in §5.2 how we can modify the label extension to consider only semi-active schedules. We then present in §5.3 an alternative definition of graph G that allows us to take into account the dominance rules introduced in Section 3.

## 5.1 Forward and backward labeling algorithms

Enumeration of all feasible (resource-constrained) paths in graph G can be performed by a labeling algorithm applying bi-directional search (i.e., a forward search from $v_{source}$ and a backward search from $v_{sink}$).

**Forward labeling algorithm** The forward labeling algorithm uses a forward label $\vec{L} = (c^{\vec{L}}, v^{\vec{L}}, Q^{\vec{L}}, J^{\vec{L}})$ to represent a path starting at $v_{source}$ and finishing at node $v^{\vec{L}} \in V$ with accumulated resource consumption $Q^{\vec{L}}$ and a cost equal to $c^{\vec{L}}$. Forward label $\vec{L}$ is also characterized by a set $J^{\vec{L}}$ of scheduled jobs. The algorithm starts with one unextended label $(0, v_{source}, (l_{v_{source},r})_r, \emptyset)$. At each iteration, it selects an unextended label $\vec{L}$ such that $v^{\vec{L}} \neq v_{sink}$, proceeds to its extension along every arc $a = (v^{\vec{L}}, v) \in \Gamma^+(v^{\vec{L}})$ such that $v \notin \cup_{j \in J^{\vec{L}}} V_j$, and marks $\vec{L}$ as extended. When extending the label $\vec{L}$ along an arc $a$, it creates an unextended label $\vec{L}'$ with the following characteristics:

$$c^{\vec{L}'} = c^{\vec{L}} + c_a$$

$$v^{\vec{L}'} = v$$

$$Q_r^{\vec{L}'} = \begin{cases} Q_r^{\vec{L}} + q_{a,r} & r = 0 \\ \max\{l_{v,r}, Q_r^{\vec{L}} + q_{a,r}\} & r \in \mathcal{R} \end{cases}$$

$$J^{\vec{L}'} = \begin{cases} J^{\vec{L}} \cup \{j\} & v \in V_j \\ J^{\vec{L}} & \text{otherwise} \end{cases}$$

The computation of $Q_r^{\vec{L}'}$ for $r \in \mathcal{R}$ can be interpreted as follows: if the resource $r$ is available at the current time $Q_0^{\vec{L}'}$ for more than $l_{v,r}$ time periods (meaning that the $Q_r^{\vec{L}} + q_{a,r} < l_{v,r}$), then $Q_r^{\vec{L}'}$ can always be at least set equal to $l_{v,r}$ as no job that would start using resource 0 at $Q_0^{\vec{L}'}$ could use $r$ before time $Q_0^{\vec{L}'} - l_{v,r}$ (according to the values set for $l_{v,r}$). The label $\vec{L}'$ is deleted if $Q_0^{\vec{L}'} \notin [l_{v,0}, u_{v,0}]$ or there exists $r \in \mathcal{R}$ such that $Q_r^{\vec{L}'} > u_{v,r}$. The algorithm continues as long as there are still unextended labels at a node $v \in V \setminus \{v_{sink}\}$. After the algorithm termination, the set of labels at $v_{sink}$ corresponds to the set of feasible resource-constrained paths in G from $v_{source}$ to $v_{sink}$.

**Backward labeling algorithm** The backward labeling algorithm uses a backward label $\vec{L} = (c^{\vec{L}}, v^{\vec{L}}, Q^{\vec{L}}, J^{\vec{L}})$ to represent a (backward) path starting at $v_{sink}$ and finishing at node $v^{\vec{L}} \in V$. It starts with one unextended label $(0, v_{sink}, (u_{v_{sink},r})_r, \emptyset)$. At each iteration, it selects an unextended label $\vec{L}$ such that $v^{\vec{L}} \neq v_{source}$ and extend it along every arc $a = (v, v^{\vec{L}}) \in \Gamma^-(v^{\vec{L}})$ (considering $a$ in its opposite direction) such that $v \notin \cup_{j \in J^{\vec{L}}} V_j$, and marks $\vec{L}$ as extended. When extending the label $\vec{L}$ along an arc $a$, we create an unextended label $\vec{L}'$ with the following characteristics:

$$c^{\vec{L}'} = c^{\vec{L}} + c_a$$

$$v^{\vec{L}'} = v$$
$$Q_r^{\tilde{L}'} = \begin{cases} Q_r^{\tilde{L}} - q_{a,r} & r = 0 \\ \min\{u_{v,r}, Q_r^{\tilde{L}} - q_{a,r}\} & r \in \mathbb{R} \end{cases}$$

$$J^{\tilde{L}'} = \begin{cases} J^{\tilde{L}} \cup \{j\} & v \in V_j \\ J^{\tilde{L}} & \text{otherwise} \end{cases}$$

The label $\tilde{L}'$ is deleted if $Q^{\tilde{L}'} \notin [l_{v,0}, u_{v,0}]$ or if there exists $r \in \mathbb{R}$ such that $Q_r^{\tilde{L}'} < l_{v,r}$. The algorithm continues as long as there are still unextended labels at a node $v \in V \setminus \{v_{source}\}$.

In the forward and backward labeling algorithms, we extend a forward label $\tilde{L}$ (resp. backward label $\tilde{L}$) in non-decreasing (resp. non-increasing) order of the vector $Q_r^{\tilde{L}}$ (resp. $Q_r^{\tilde{L}}$), with two labels being compared using the lexicographical order on the Cartesian product $\mathbb{Z}^{m+1}$. Since $G$ only contains cycles where the consumption of resource 0 increases (resp. decreases) by at least 1, these algorithms are label-setting algorithms, meaning that they only extend labels that can never be dominated by another label created after that extension. The maximum number of labels in these algorithms is bounded above by $1 + 2^n((t^2 - t^1) \prod_{r \in \mathbb{R}} (\hat{p}_r^+ - \hat{p}_r^-) + \sum_{j \in \mathcal{J}} \sum_{k=1}^{w_j} ((t_{jk}^2 - t_{jk}^1 - p_j) \prod_{r \in \mathbb{R} \setminus \{r_j\}} (\hat{p}_r^+ - \hat{p}_r^-)))$. In the worst case, each label is extended for every possible addition of a job to the partial schedule The time and space complexity of the labeling algorithms are therefore exponential in the worst case.

**Bi-directional labeling algorithm** We run the forward and backward algorithms where only labels $\tilde{L}$ such that $Q_0^{\tilde{L}} \le t^*$ and labels $\tilde{L}$ such that $Q_r^{\tilde{L}} > t^*$ are kept (an extension that leads to a label that does not have this property is deleted). The value of $t^*$ is set equal to $\frac{t_1+t_2}{2}$. We then compute all the paths from $v_{source}$ to $v_{sink}$ by concatenating a forward path, an arc, and a backward path. We can concatenate a (partial) path represented by a forward label $\tilde{L}$ and a (partial) backward path represented by a backward label $\tilde{L}$ along an arc $a = (v^{\tilde{L}}, v^{\tilde{L}}) \in A$ if the following conditions are satisfied:

$$Q_0^{\tilde{L}} + q_{a,0} = Q_0^{\tilde{L}} \\ Q_r^{\tilde{L}} + q_{a,r} \le Q_r^{\tilde{L}} \quad r \in \mathbb{R} \\ J^{\tilde{L}} \cap J^{\tilde{L}} = \emptyset$$

The cost of the resulting path is equal to $c^{\tilde{L}} + c_a + c^{\tilde{L}}$.

The forward (resp. backward) labels computed when traversing forward (resp. backward) the path associated with the solution of Figure 1 are as follows:

<figure><img src="image_5.png" /></figure>

## 5.2 Semi-active schedules, modification of label extension and dominance checks

The set of semi-active schedules (i.e., schedules that do not admit local left shifts for every job) is dominant for the PC-JSOCMSR, as any schedule can be transformed into a semi-active one. We can therefore restrict the search of an optimal solution to the PC-JSOCMSR to these type of schedules. Due to the idle arc $(v_*, v_*)$, schedules associated to feasible solutions of the RCSPP are not necessary semi-active. To limit the enumeration to semi-active schedules, we remove this idle arc from $A$ and adapt the RCSPP definition. The resource consumption accumulated in a forward (resp. backward) path up to $v$ should be less than or equal to $u_{v,0}$ (resp. larger than or equal to $l_{v,0}$). The update of the accumulated consumption of the secondary resources according
to the intervals defined at each node is done in a similar way as described before, but an additional update is introduced to handle the insertion of idle times. Let $Q_0$ be the accumulated consumption of resource 0 (i.e., the current time) and $Q_r$ be the accumulated consumption of a secondary resource $r$. In a forward path up to a node $v$, if $Q_0$ is less than $l_{v,0}$ (the job cannot start using resource 0), then the time should increase at least by $l_{v,0} - Q_0$. Similarly, if $Q_r \ge u_{v,r}$ (resource $r$ is not available for long enough at the current time), then the time should increase at least by $Q_r - u_{v,r}$ as this results in subtracting the opposite of this value to $Q_r$. The path corresponds to a feasible schedule if the maximum of these minimum time increases (corresponding to a number of idle times) is such that the resulting time is less than or equal to $u_{v,0}$. Similarly, in a backward path up to a node $v$, if $Q_0$ is larger than $u_{v,0}$, then the time should decrease at least by $Q_0 - u_{v,0}$. If $Q_r \le l_{v,r}$, then the time should decrease at least by $l_{v,r} - Q_r$. The backward path corresponds to a feasible schedule if the maximum of these minimum time decreases is such that the resulting time is larger than or equal to $l_{v,0}$. Unlike forward paths that are associated with semi-active schedules, a backward path is associated with a schedule that does not admit local right shifts for every job. Further in the text, when concatenating a forward path with a backward path, we derive a semi-active schedule by local left shifting of the jobs of the latter path.

We describe below the modifications that are required in the forward, backward, and bi-directional algorithms. The time and space complexity of these algorithms in the worst case remains exponential. Pseudo-code of these algorithms are available in §S.2 of the supplementary material.

### Forward labeling algorithm

The value of the accumulated resource consumption of a label $\vec{L}'$, result of the extension of $\vec{L}$ along an arc $a = (v^{\vec{L}}, v)$, becomes:

$$Q_r^{\vec{L}'} = \begin{cases} Q_r^{\vec{L}} + q_{a,r} + \bar{\Delta}(\vec{L}, a) & r = 0 \\ \max\{l_{a,r}, Q_r^{\vec{L}} + q_{a,r} - \bar{\Delta}(\vec{L}, a)\} & r \in \mathcal{R} \end{cases}$$

where

$$\bar{\Delta}(\vec{L}, a) = \max \left\{ 0, l_{v,0} - Q_0^{\vec{L}} - q_{a,0}, \max_{r \in \mathcal{R}} \left\{ Q_r^{\vec{L}} + q_{a,r} - u_{v,r} \right\} \right\}$$

The value of $\bar{\Delta}(\vec{L}, a)$ corresponds to the minimum idle time before the job associated with node $v$ starts using the main resource. The value $l_{v,0} - Q_0^{\vec{L}} - q_{a,0}$ corresponds to the idle time before the beginning of the time window of this job, while the value $\max_{r \in \mathcal{R}} \{Q_r^{\vec{L}} + q_{a,r} - u_{v,r}\}$ corresponds to the idle time necessary to ensure that the secondary resources are available for long enough. The value of $Q_0^{\vec{L}'}$ ultimately represents the time at which the job (possibly dummy) associated with node $v$ starts using the main resource. The label $\vec{L}'$ is deleted if $Q_0^{\vec{L}'} > u_{v,0}$.

### Backward labeling algorithm

The value of accumulated resource consumption of a label $\vec{L}'$, result of the extension of $\vec{L}$ along an arc $a = (v, v^{\vec{L}})$ (traversed backward), becomes:

$$Q_r^{\vec{L}'} = \begin{cases} Q_r^{\vec{L}} - q_{a,r} - \bar{\Delta}(\vec{L}, a) & r = 0 \\ \min\{u_{v,r}, Q_r^{\vec{L}} - q_{a,r} + \bar{\Delta}(\vec{L}, a)\} & r \in \mathcal{R} \end{cases}$$

where

$$\bar{\Delta}(\vec{L}, a) = \max \left\{ 0, Q_0^{\vec{L}} - q_{a,0} - u_{v,0}, \max_{r \in \mathcal{R}} \left\{ l_{v,r} - (Q_r^{\vec{L}} - q_{a,r}) \right\} \right\}$$

The value of $\bar{\Delta}(\vec{L}, a)$ represents a number of idle time periods for the main resource. The label $\vec{L}'$ is deleted if $Q_0^{\vec{L}'} < l_{v,0}$.

### Bi-directional labeling algorithm

The two first conditions for concatenating (a path represented by) a forward label $\vec{L}$ and (a path represented by) a backward label $\vec{L}$ along an arc $a = (v^{\vec{L}}, v^{\vec{L}})$ become:

$$\begin{aligned} Q_0^{\vec{L}} + q_{a,0} &\le Q_0^{\vec{L}} \\ Q_r^{\vec{L}} + q_{a,r} - (Q_0^{\vec{L}} - (Q_0^{\vec{L}} + q_{a,0})) &\le Q_r^{\vec{L}} \quad r \in \mathcal{R} \end{aligned}$$
The solution of Figure 1 is a semi-active schedule and can therefore be represented (due to the modification of the labeling algorithms) as the following path in the graph G without idle arc: $(v_{source}, v_{\star})$, $(v_{\star}, v_{11})$, $(v_{11}, v_{\star})$, $(v_{\star}, v_{41})$, $(v_{41}, v_{\star})$, $(v_{\star}, v_{31})$, $(v_{31}, v_{\star})$, $(v_{\star}, v_{sink})$. The forward (resp. backward) labels computed when traversing the path forward (resp. backward) are as follows:

$$ (\text{Forward}) \quad (0, v_{\text{source}}, 0, 0, 0, \{\}) \to (0, v_{\star}, 0, 0, 0, \{\}) \to (0, v_{11}, 2, -2, -2, \{\}) \to (4, v_{\star}, 4, -3, 1, \{J_1\}) \to (4, v_{41}, 6, -2, -1, \{J_1\}) \to \\ (6, v_{\star}, 9, 2, -2, \{J_1, J_4\}) \to (6, v_{31}, 9, 2, -1, \{J_1, J_4\}) \to (8, v_{\star}, 10, 1, 1, \{J_1, J_3, J_4\}) \to (8, v_{\text{sink}}, 15, 0, 0, \{J_1, J_3, J_4\}) $$

$$ (\text{Backward}) \quad (0, v_{\text{sink}}, 15, 0, 0, \{\}) \to (0, v_{\star}, 15, 0, 0, \{\}) \to (2, v_{31}, 10, 2, -1, \{J_3\}) \to (2, v_{\star}, 10, 2, -1, \{J_3\}) \to (4, v_{41}, 7, -2, 1, \{J_3, J_4\}) \to \\ (4, v_{\star}, 7, -2, 1, \{J_3, J_4\}) \to (8, v_{11}, 5, 0, -2, \{J_1, J_3, J_4\}) \to (8, v_{\star}, 5, 0, -2, \{J_1, J_3, J_4\}) \to (8, v_{\text{source}}, 0, 0, 0, \{J_1, J_3, J_4\}) $$

To eliminate labels and accelerate the labeling algorithm during the forward and backward search, we perform dominance checks. The following conditions (with at least one inequality or inclusion strictly satisfied) are sufficient to state that a forward label $\vec{L}_1$ dominates a forward label $\vec{L}_2$:

$$ \begin{array}{l} c^{\vec{L}_1} \ge c^{\vec{L}_2} \\ v^{\vec{L}_1} = v^{\vec{L}_2} \\ Q_0^{\vec{L}_1} \le Q_0^{\vec{L}_2} \\ Q_r^{\vec{L}_1} - (Q_0^{\vec{L}_2} - Q_0^{\vec{L}_1}) \le Q_r^{\vec{L}_2} \quad r \in \mathcal{R} \\ J^{\vec{L}_1} \subseteq J^{\vec{L}_2} \end{array} $$

With these conditions, every feasible completion of the partial path associated with $\vec{L}_2$ into a path $\vec{P}_2$ from $v_{source}$ to $v_{sink}$ is also a feasible completion of the partial path associated with $\vec{L}_1$ into a full path $\vec{P}_1$, and the cost of $\vec{P}_2$ is always lower than the cost of $\vec{P}_1$. The third condition considers the possibility of waiting the difference in time between the time associated with $\vec{L}_1$ and the time associated with $\vec{L}_2$ and checks if the secondary resources can be available for a longer time.

Similarly to the forward search, the following conditions are sufficient to state that a backward label $\vec{L}_1$ dominates a backward label $\vec{L}_2$:

$$ \begin{array}{l} c^{\vec{L}_1} \ge c^{\vec{L}_2} \\ v^{\vec{L}_1} = v^{\vec{L}_2} \\ Q_0^{\vec{L}_1} \ge Q_0^{\vec{L}_2} \\ Q_r^{\vec{L}_1} + (Q_0^{\vec{L}_1} - Q_0^{\vec{L}_2}) \ge Q_r^{\vec{L}_2} \quad r \in \mathcal{R} \\ J^{\vec{L}_1} \subseteq J^{\vec{L}_2} \end{array} $$

## 5.3 Alternative definition of graph G

Dominance checks are performed only between labels ending at the same node. An alternative definition of graph G has therefore a potential impact on the efficiency of solving the RCSPP. We introduce a new graph by removing node $v_{\star}$ and concatenating every couple of arcs $((v_{j,k}, v_{\star})(v_{\star}, v_{j',k'}))$ into a single arc. It has the interest of forbidding two consecutive occurrences of the same job directly by the structure of the graph. The introduction of this new graph allows to take into account the dominance rules introduced in Section 3.

We denote $G' = (V', A')$ the new graph with $V' = V \setminus \{v_{\star}\}$. The set $A'$ is different from A and does not contain any idle arc. For every job $j \in \mathcal{J}$, it contains the arc $(v_{source}, v_{j,k})$ for $k = 1$ and the arc $(v_{j,k}, v_{sink})$ for every $k$. It also contains an arc $(v_{j,k}, v_{j',k'})$ if job j can be executed in time window $T_{jk}$ and use the main resource immediately before its use by job $j' \neq j$ executed in time window $T_{j'k'}$ (i.e, $t_{jk}^1 + p_j \le t_{j'k'}^2 - p_{j'}$ if $r_j = r_{j'}$ or $t_{jk}^1 + p_j^- + p_j^0 \le t_{j'k'}^2 - p_{j'} + p_{j'}^-$ if $r_j \neq r_{j'}$). The consumption $q_{a,r}$ of resource r in arc a is defined as follows:

$$ q_{a,r} = \begin{cases} 0 & \text{if } a \in \Gamma^+(v_{source}) \\ p_j^0 & \text{if } a \in \Gamma^+(v_{j,k}) \text{ and } r = 0 \\ -p_j^0 & \text{if } a \in \Gamma^+(v_{j,k}) \text{ and } r \in \mathcal{R} \setminus \{r_j\} \\ p_j^- + p_j^+ & \text{if } a \in \Gamma^+(v_{j,k}) \text{ and } r = r_j \end{cases} $$
The accumulated resource consumption interval $[l_{v,r}, u_{v,r}]$ for each node $v \in V'$ and each resource $r \in R_0$ is identical to the one described for graph $G$.

The labeling algorithms used to enumerate all feasible paths from $v_{source}$ to $v_{sink}$ in $G'$ are strictly identical to those already presented.

Using the dominance rules introduced in Section 3, we now describe how we reduce the number of arcs of $G'$ while ensuring that at least one of the optimal solutions to the RCSPP is retained. We recall that using an arc $(v_{j,k}, v_{j',k'})$ in a path of graph $G'$ translates into a schedule where the main resource is used by job $j$ during time window $T_{jk}$ and then used by job $j'$ during time window $T_{j'k'}$. The arc $(v_{j',k'}, v_{j,k})$ can be safely removed from $G'$ if the sufficient conditions of Proposition 1 are satisfied. The arc $(v_{j'',k''}, v_{j',k'})$ (resp. $(v_{source}, v_{j',k'})$) can be safely removed if the sufficient conditions stated in item 1 of Proposition 2 are satisfied. Similarly, the arc $(v_{j',k'}, v_{j'',k''})$ (resp. $(v_{j',k'}, v_{sink})$) can be safely removed if the sufficient conditions stated in item 2 of Proposition 2 are satisfied.

We show in Figure 3 the graph $G'$ associated with the instance of Figure 1. The solution of Figure 1 is represented in this graph by the following path: $(v_{source}, v_{11}), (v_{11}, v_{41}), (v_{41}, v_{31}), (v_{31}, v_{sink})$. Using Proposition 1 with $j = J_1$, $j' = J_3$, $k = 1$ and $k' = 1$, we remove the dotted orange arc $(v_{31}, v_{11})$. We remove the dashed red arcs using Proposition 2 with $j = J_4$ and $j', j'', k'$ and $k''$ defined according to the arc to remove. The arcs in gray are removed due to time window incompatibility.

The forward (resp. backward) labels computed when traversing the path forward (resp. backward) are as follows:

$$ \text{(Forward)} \quad (0, v_{source}, 0, 0, 0, \{\}) \rightarrow (0, v_{11}, 2, -2, -2, \{J_1\}) \rightarrow (4, v_{41}, 6, -2, -1, \{J_1, J_4\}) \rightarrow (6, v_{31}, 9, 2, -1, \{J_1, J_3, J_4\}) \rightarrow \\ \qquad (8, v_{sink}, 15, 0, 0, \{J_1, J_3, J_4\}) $$

$$ \text{(Backward)} \quad (0, v_{sink}, 15, 0, 0, \{\}) \rightarrow (2, v_{31}, 10, 2, -1, \{J_3\}) \rightarrow (4, v_{41}, 7, -2, 1, \{J_3, J_4\}) \rightarrow (8, v_{11}, 5, 0, -2, \{J_1, J_3, J_4\}) \rightarrow \\ \qquad (8, v_{source}, 0, 0, 0, \{J_1, J_3, J_4\}) $$

<figure><img src="image_6.png" /><figcaption>Figure 3: Graph $G'$ associated with the instance of Figure 1</figcaption></figure>
# 6 A branch-cut-and-price algorithm

Due to the well-known curse of dimensionality, solving the PC-JSOCMSR as a RCSPP using the labeling algorithm described in Section 5 is not practical. The main reason explaining the explosion of the number of labels despite the dominance checks is related to dealing with the elementarity constraints. To overcome this issue, we propose to take these constraints into account in a different way.

We first introduce an extended MILP model for the PC-JSOCMSR. Let $S$ be the set of schedules satisfying time window restrictions and resource usage constraints. We denote by $b_S$ the profit associated with a schedule $S \in S$ (sum of the individual profit of the jobs scheduled in $S$) and by $\tilde{y}_j^S \in \mathbb{N}$ the number of times job $j$ is scheduled in $S$ ($\tilde{y}_j^S$ can be larger than 1 since $S$ does not necessarily satisfy elementarity constraints). Our new model, where the binary variable $\lambda_S$ equal to 1 if and only if schedule $S$ is selected, is as follows:

$$ [F3] \quad \max \sum_{S \in S} b_S \lambda_S \qquad (17) $$

$$ \text{s.t.} \quad \sum_{S \in S} \tilde{y}_j^S \lambda_S \le 1 \qquad j \in \mathcal{J} \qquad (18) $$

$$ \sum_{S \in S} \lambda_S = 1 \qquad (19) $$

$$ \lambda_S \in \{0, 1\} \qquad S \in S \qquad (20) $$

The objective function (17) maximizes the total profit associated with the selected schedules. Constraints (18) impose the elementarity constraints. Constraint (19) imposes the selection of a single schedule. Constraints (20) set the domains of the decision variables.

We introduce a branch-cut-and-price (BCP) algorithm to solve formulation [F3]. Since the number of variables is exponential in the number of jobs, we apply column generation to solve its linear relaxation. The variables $\lambda$ are dynamically generated by solving the pricing problem described in §6.2 and defined as a RCSPP on graph $G$ or $G'$ (with some adjustments). To tighten the dual bound, we generate cuts derived from Chvátal-Gomory rank-1 cuts (21). This is explained in §6.1. The whole column and cut generation procedure is detailed in §6.3. We finally describe in §6.4 the branching strategy we use to discard non-integer solutions.

## 6.1 Chvátal-Gomory rank-1 cuts

Let $J \subseteq \mathcal{J}$ and $(\mu_j)_{j \in J}$ be a vector of strictly positive multipliers, a Chvátal-Gomory rounding of set packing constraints (18) results in the following Chvátal-Gomory rank-1 cut (R1C)

$$ \sum_{S \in S} \left| \sum_{j \in J} \mu_j \tilde{y}_j^S \right| \lambda_S \le \left| \sum_{j \in J} \mu_j \right| \qquad (21) $$

that is a valid inequality for [F3]. R1Cs are non-robust in the sense they impact the structure of the pricing problem (due to the floor rounding operation). However, their separation can improve the quality of the Lagrangian dual bound. As a solution to [F3] is a single schedule, we only consider R1Cs such that $\left| \sum_{j \in J} \mu_j \right| < 1$ (i.e., $|J| = 1$) as other R1Cs are practically never violated due to constraint (19). Specifically, we separate the following R1Cs:

$$ \sum_{S \in S} |0.5 \tilde{y}_j^S| \lambda_S \le 0 \quad j \in \mathcal{J} \qquad (22) $$

These cuts can be violated by the linear relaxation solution since the schedules in $S$ do not necessarily satisfy the elementarity constraints.

To mitigate the impact of these cuts on the pricing problem difficulty, we consider the limited memory version of R1Cs (22) introduced by Pecin et al. (2017b). Each R1C is associated with an arc memory set $AM \subseteq A$ as described in Pecin et al. (2017a). With the arc memory set, the coefficient $|0.5 \tilde{y}_j^S|$ of variable $\lambda_S$ is replaced by a smaller coefficient $\varphi(j, AM, S)$. The limited-arc-memory R1C cut associated with job $j \in \mathcal{J}$ reads
as follows:

$$ \sum_{S \in \mathcal{S}} \varphi(j, AM, S) \lambda_S \le 0 \qquad (23) $$

Algorithm 1 describes the computation of $\varphi(j, AM, S)$. The variable $e$ is referred to as the state variable. Depending on its value, earlier occurrences of job $j$ may be forgotten when an arc does not belong to the memory.

<caption>**Algorithm 1:** Computing $\varphi(j, AM, S)$</caption>

```
\varphi \leftarrow 0, e \leftarrow 0
    \text{foreach } a \text{ arc of schedule } S \text{ do}
        \text{if } a \notin AM \text{ then } e \leftarrow 0
        \text{if } a \in \bigcup_{v \in V_j} \Gamma^-(v) \text{ then}
            \text{e} \leftarrow e + 0.5
            \text{if } e = 1 \text{ then } \varphi \leftarrow \varphi + 1, e \leftarrow 0
        \text{end}
    \text{end}
    \text{return } \varphi
```

The arc memory set $AM$ is computed as the minimal subset of arcs that ensure that the coefficient in (23) of each variable $\lambda_S$ with a positive value in the current linear relaxation solution to [F3] is equal to its coefficient in (22) (see Pecin et al. (2017b) for details). If an arc $(v, v')$ is added to $AM$, then the arc $(v', v)$ is also added to $AM$.

Using a limited memory version of R1Cs implies that a cut of type (22) for a job $j$ may be violated by the current linear relaxation solution although a limited memory version of this cut exists. This means that a different arc memory set $AM'$ must be used to cut off the solution. In that case, we create a new cut combining $AM'$ and the arc memory sets of existing cuts associated with $j$.

## 6.2 Pricing problem

Let $(\pi_j)_{j \in \mathcal{J}}$ be the dual variables associated with constraints (18), $\zeta$ be the dual variable associated with constraints (19), and $(\sigma_g)_{g \in \mathcal{G}}$ be the dual variables associated with active R1Cs. We denote by $j(g)$ the job associated with R1C $g \in \mathcal{G}$ and $AM(g)$ its arc memory set. The reduced cost of a schedule $S$ is equal to:

$$ \bar{c}_S = \sum_{j \in \mathcal{S}} (b_j - \pi_j) - \sum_{g \in \mathcal{G}} \varphi(j(g), AM(g), S) \sigma_g - \zeta \qquad (24) $$

The pricing problem can be stated as: $\max_{S \in \mathcal{S}} \{\bar{c}_S\}$.

We define the pricing subproblem as a RCSPP on the directed graph $G = (V, A)$ (or indifferently on $G' = (V', A')$) introduced in Section 5. The cost $c_a$ of an arc $a$ is replaced by $\bar{c}_a$ its reduced cost defined as:

$$ \bar{c}_a = \begin{cases} b_j - \pi_j & \text{if } a \in \Gamma^+(v_{j,k}) \\ 0 & \text{otherwise} \end{cases} $$

Introducing $\tilde{x}_a^S$ as a parameter equal to the number of times arc $a$ belongs to the path in one-to-one correspondence with schedule $S$, we can rewrite (24) as follows:

$$ \bar{c}_S = \sum_{a \in A} \tilde{x}_a^S \bar{c}_a - \sum_{g \in \mathcal{G}} \varphi(j(g), AM(g), S) \sigma_g - \zeta \qquad (25) $$

We solve the pricing problem applying the bi-directional algorithm described in §5.1 with the adaptions of §5.2 for generating paths only associated with semi-active schedules. We now describe additional modifications that are required, in particular in the definition of the labels, their extension, and the dominance checks. Pseudo-code of these algorithms are available in §S.2 of the supplementary material.
**Forward labeling algorithm** The forward labeling algorithm uses a forward label $\vec{L} = (\vec{c}^L, v^L, Q^L, E^L)$ to represent a path starting at $v_{source}$ and finishing at node $v^L \in V$ with accumulated resource consumption $Q^L$ and a reduced cost of $\vec{c}^L$ (in the place of $c^L$). Forward label $\vec{L}$ is also characterized by a vector of states $E^L$ (see the state variable introduced in §6.1) associated with R1Cs active in the current linear relaxation solution (the set $J^L$ is dropped). When extending the label $\vec{L}$ along an arc $a = (v^L, v)$, the resulting label $\vec{L}'$ has the following characteristics (only differences are described):

$$ \vec{c}^{\vec{L}'} = \vec{c}^{\vec{L}} + \bar{c}_a - \sum_{g \in \mathcal{G}} [a \in AM(g) \wedge v \in V_{j(g)} \wedge E_g^{\vec{L}} + 0.5 = 1] \sigma_g \\ E_g^{\vec{L}'} = \begin{cases} 0 & (a \notin AM(g) \wedge v \notin V_{j(g)}) \vee (AM(g) \wedge v \in V_{j(g)} \wedge E_g^{\vec{L}} + 0.5 = 1) \\ 0.5 & a \notin AM(g) \wedge v \in V_{j(g)} \\ E_g^{\vec{L}} & a \in AM(g) \wedge v \notin V_{j(g)} \\ E_g^{\vec{L}} + 0.5 & a \in AM(g) \wedge v \in V_{j(g)} \wedge E_g^{\vec{L}} + 0.5 < 1 \end{cases} \quad g \in \mathcal{G} $$

where the square brackets are Iverson brackets ($[P] = 1$ if statement P is true, and $[P] = 0$ otherwise). When checking if a forward label $\vec{L}_1$ dominates another forward label $\vec{L}_2$, we replace the first condition (related to the label costs) by:

$$ \vec{c}^{\vec{L}_1} \ge \vec{c}^{\vec{L}_2} + \sum_{g \in \mathcal{G}: E_g^{\vec{L}_1} > E_g^{\vec{L}_2}} \sigma_g $$

**Backward labeling algorithm** The backward labeling algorithm uses a backward label $\vec{L} = (\vec{c}^L, v^L, Q^L, E^L)$ to represent a (backward) path starting at $v_{sink}$ and finishing at node $v^L \in V$. When extending the label $\vec{L}$ along an arc $a = (v, v^L)$, the resulting label $\vec{L}'$ has the following new characteristics:

$$ \vec{c}^{\vec{L}'} = \vec{c}^{\vec{L}} + \bar{c}_a - \sum_{g \in \mathcal{G}} [a \in AM(g) \wedge v \in V_{j(g)} \wedge E_g^{\vec{L}} + 0.5 = 1] \sigma_g \\ E_g^{\vec{L}'} = \begin{cases} 0 & (a \notin AM(g) \wedge v \notin V_{j(g)}) \vee (AM(g) \wedge v \in V_{j(g)} \wedge E_g^{\vec{L}} + 0.5 = 1) \\ 0.5 & a \notin AM(g) \wedge v \in V_{j(g)} \\ E_g^{\vec{L}} & a \in AM(g) \wedge v \notin V_{j(g)} \\ E_g^{\vec{L}} + 0.5 & a \in AM(g) \wedge v \in V_{j(g)} \wedge E_g^{\vec{L}} + 0.5 < 1 \end{cases} \quad g \in \mathcal{G} $$

When checking if a backward label $\vec{L}_1$ dominates another backward label $\vec{L}_2$, we replace the first condition (related to the label costs) by:

$$ \vec{c}^{\vec{L}_1} \ge \vec{c}^{\vec{L}_2} + \sum_{g \in \mathcal{G}: E_g^{\vec{L}_1} > E_g^{\vec{L}_2}} \sigma_g $$

Using the same reasoning as in §5.1, the time and space complexity of the labeling algorithms is pseudo-polynomial in the worst case if $\mathcal{G} = \emptyset$, while otherwise it is exponential.

**Bi-directional labeling algorithm** The last condition (related to elementarity constraints) for concatenating a forward path with a backward path is dropped. The reduced cost of a path obtained as the concatenation of a (partial) path represented by a forward label $\vec{L}$ and a (partial) backward path represented by a backward label $\vec{L}$ along an arc $a = (v^L, v)$ is equal to $\vec{c}^L + \bar{c}_a + \vec{c}^L - \sum_{g \in \mathcal{G}: E_g^{\vec{L}} + E_g^{\vec{L}} \ge 1} \sigma_g$.

**Bucket graph** We use a bucket graph implementation of the labeling algorithms as proposed in (Sadykov et al., 2021). Labels are grouped and stored in buckets to improve the efficiency of the dominance checks. At each node $v$, we create $u_{v,0} - l_{v,0} + 1$ buckets, each bucket being associated with a specific time. Every label corresponding to a partial path ending up at the same node and with the same accumulated resource consumption of the main resource (i.e., with the same time) are stored in the same bucket. Using buckets helps reducing the number of dominance checks and the computational time needed for solving the pricing problem.
Dominance checks are performed only with labels within the same bucket when testing the insertion of a new label. Before extending a label, dominance checks are performed by considering the labels stored in all buckets created at a node and associated with a lower time. Let us denote each bucket as a couple $(v, t)$ of a node $v \in V$ and a time $t \in [l_{v,0}, u_{v,0}]$. In the forward (resp. backward) bucket graph implementation of $G$, for every bucket $(v, t)$ and every arc $a = (v, v')$ (resp. $a = (v', v)$), we create an arc between bucket $(v, t)$ and bucket $(v', \max\{l_{v',0}, t+q_{a,0}\})$ (resp. an arc between bucket $(v, t)$ and bucket $(v', \min\{u_{v',0}, t-q_{a,0}\})$. The bucket graph implementation of the graphs $G$ and $G'$ displayed in Figures 2 and 3 are given in §S.3 of the supplementary material. The bucket graph is used for the reduced cost-based arc elimination procedure described in §6.3.

## 6.3 Column and cut generation algorithm

We use an adaptation of the BCP algorithm introduced by Sadykov et al. (2021). To make the paper self-contained, we describe the main components of this algorithm.

For each node of the branch-and-bound tree, we compute the linear relaxation of the corresponding restricted master program using a three-stage column generation procedure. The first two stages use heuristic pricing and generates at most 30 columns per call, whereas the last stage uses exact pricing and generates at most 150 columns per call. During the first stage, the pricing problem is solved by keeping only in every bucket the eight labels with the largest reduced cost. During the second phase, the dominance checks do not include the conditions on the R1Cs. The last stage performs full dominance checks. After each pricing in this stage, the threshold value $t^*$ for the bidirectional labeling algorithm (initially set equal to $(t^1 + t^2)/2$), is automatically adjusted as follows: if the number of non-dominated labels in the forward (resp. backward) algorithm is larger than the number of non-dominated labels in the backward (resp. forward) algorithm, then $t^*$ is increased (resp. decreased) by 5%.

To improve the column generation convergence, we use the automatic dual price smoothing technique described in Pessoa et al. (2018) in each stage.

When the number of columns in the restricted master program exceeds 10000, we remove 34% of the columns by selecting the non-basic columns whose reduced cost is the lowest.

When column generation converges, we eliminate arcs of graph $G$ (or $G'$) based on reduced costs. If the maximum reduced cost of a path using an arc $a \in A$ is smaller than the difference between the primal bound and the dual bound, then $a$ can be removed from $G$ as it does not belong to any improving solution to [F3] (see Irnich et al. (2010)). The value of a feasible heuristic solution to the PC-JSOCMSR is given as input and initializes the value of the best primal bound. To test if an arc $a = (v, v')$ should be removed, we test all the combinations of a forward label $\vec{L}$ such that $v^{\vec{L}} = v$ with arc $a$ and a backward label $\bar{L}$ such that $v^{\bar{L}} = v'$ that lead to a feasible path. As the labeling algorithms described in §6.2 use the bucket graph implementation of Sadykov et al. (2021), we apply the *bucket arc elimination procedure* described in this latter work. Specifically, using the condition previously described, we remove arcs in the bucket graph implementation of $G$ (or $G'$) rather than in $G$ (or $G'$). For the validity of the algorithm, we need to introduce “jump” arcs between the buckets of the same node associated with successive times. A jump arc $a$ between bucket $(v, t)$ and bucket $(v, t')$ is such that $c_a = 0$, $q_{a,0} = t' - t$, and $q_{a,r} = t - t'$ for every $r \in \mathbb{R}$. Labels are also extended along these arcs. The bucket arc elimination procedure is stronger than a classical arc elimination procedure in $G$, since an arc of $G$ can only be eliminated based on reduced cost only if all the arcs associated with it in the bucket graph implementation of $G$ can be eliminated. The arc elimination procedure is performed after the initial convergence and also each time the current dual-primal gap decreased by at least 10% since it was last applied.

We separate the R1Cs (22) by enumeration. In each cut generation round, we consider at most 50 R1Cs (the most violated in priority) and add the limited memory version (23) of these cuts.

Column and cut generation is stopped either when after 3 rounds of cut generation the dual-primal gap decreases by less than 2% per round (tailing-off condition), or when the time spent to solve the pricing subproblem exceeds 5 seconds (time-efficiency condition).

If the current dual-primal gap is less than 1%, we try to solve the pricing problem using the bidirectional labeling algorithm described in §5.2 and the additional constraint that the reduced cost should be larger than
the difference between the primal bound and the dual bound (reduced costs are considered like in §6.2). For efficiency consideration, we stop the algorithm if more than 1,000,000 forward (backward) labels are generated. If the algorithm succeeds (i.e., it has not been stopped), the pricing subproblem is solved by inspection in future column generation iterations. If the total number of enumerated paths is less than 10,000, all the paths are added to S and formulation [F3] is solved with a MILP solver. An optimal solution at the current node is obtained, and this node can be pruned (after updating the current best primal bound if necessary). Otherwise, we proceed to branching. The branch-and-bound search tree is explored using a depth-first search strategy. Among the nodes with same depth, the algorithm selects the one with the largest dual bound.

## 6.4 Branching

To reach integrality, we perform branching on the value of the following aggregated variables:

* the number of times a job j is performed:

$$ \sum_{S \in S} \tilde{y}_j^S \lambda_S = \sum_{S \in S} \sum_{a \in \bigcup_{v \in V_j} \Gamma^+(v)} \tilde{x}_a^S \lambda_S $$

* the number of times a job j is performed during a particular time window $k \in [1, w_j]$:

$$ \sum_{S \in S} \sum_{a \in \Gamma^+(v_{j,k})} \tilde{x}_a^S \lambda_S $$

* the number of times arc $a \in A$ is used:

$$ \sum_{S \in S} \tilde{x}_a^S \lambda_S $$

Adding a branching constraint to the (restricted) master linear program does not make the pricing problem harder to solve. The dual variable $\psi$ associated with a generic branching constraint written as $\sum_{S \in S} \sum_{a \in A} h_a \tilde{x}_a^S \lambda_S \le H$ is considered in the pricing problem by subtracting $h_a \psi$ from $\bar{c}_a$ for every $a \in A$.

After performing the previous branching strategy, we may end up with a fractional solution. In that case, we branch over accumulated consumption of resource 0 (Pessoa et al., 2021). Let $Q_a^S$ be the accumulated consumption of resource 0 at the origin of arc $a$ in schedule S (i.e., the time at which the current occurrence of the job associated with $a$ starts using the main resource). We compute for each job j the value $t_j^*$ such that $\sum_{S \in S} \sum_{a \in \bigcup_{v \in V_j} \Gamma^+(v)} [Q_a^S \le t_j^*] \lambda_S$ is the most fractional. If there exists a threshold $t_j^*$ for a job j leading to a fractional value, then we create two branches: in the left branch we set $u_{v,0} = t_j^*$ for every $v \in V_j$ such that $u_{v,0}$ is larger than $t_j^*$, and in the right branch we set $l_{v,0} = t_j^*$ for every $v \in V_j$ such that $l_{v,0}$ is less than $t_j^*$.

We use the three-phase evaluation procedure described by Sadykov et al. (2021) to choose the most promising branching candidate among the aggregated variable (presented above) with a fractional value. This branching procedure is inspired by strong branching. First, it chooses at most 30 branching candidates: half are chosen based on the history (if not empty) of previous calls to the branching procedure, and the remaining candidates are chosen according to the distance between their fractional value and the nearest integer. Second, it solves the LP relaxation of the restricted master problem modified with the corresponding branching constraint. Third, the procedure selects the three candidates with the largest product of dual bound improvements in the two branches, and performs column generation using heuristic pricing. The candidate with the largest product of dual bound improvements in the two branches is selected for branching.

## 7 A heuristic algorithm

We introduce a new heuristic algorithm for the PC-JSOCMSR with two goals. First, we aim at providing a good initial solution quickly to the BCP algorithm as it is essential to speed-up its convergence (for arc elimination
and nodes pruning). Second, we aim at computing high quality solutions for large-sized instances that remain beyond the reach of an exact algorithm.

Our heuristic is an iterated local search (ILS) algorithm. ILS is a metaheuristic that iteratively applies a local search phase to produce local optima, and a perturbation mechanism to escape from them. We use a large neighborhood search (LNS) for the local search phase. Algorithm 2 outlines the general structure of our heuristic. First, we build a solution *s* using the LNS repair (insertion) operator (see §7.1). The algorithm then enters an iterative process. At every iteration, the LNS phase (see §7.1) starts from solution *s* and accepts only improving solutions. It stops when either the number of iterations reaches $\delta^{max}$ or the number of iterations without improvements reaches $\theta^{max}$. If appropriate, the algorithm updates the current best solution $s^*$. Then, the algorithm performs the perturbation phase (see §7.2) to avoid being trapped in a local optimum that cannot be escaped (or escaped in a reasonable number of iterations) by successive application of the destroy and reparator operators. It modifies $s^*$ to produce a new solution *s* for the next LNS phase with the aim of exploring a new region of the search space. We only accept the solution produced by the perturbation operator if its relative gap with respect to the best solution computed so far is less than or equal to $\epsilon$. This is to avoid too much re-optimization and to increase the possibility to find improving solutions with the LNS. The algorithm stops when a total of $\delta^{max}$ iterations of the LNS have been performed out of all iterations. The remainder of this section describes the LNS and perturbation phases.

<caption>**Algorithm 2: ILS**</caption>

input: two positive integer values $\delta^{max}$ and $\theta^{max}$, a float value $\epsilon$
Build $s^0$ using the LNS repair operator
$\delta \leftarrow 0, s^* \leftarrow s^0, s \leftarrow s^0$
while $\delta < \delta^{max}$ do
  $\theta \leftarrow 0$
  while $\delta < \delta^{max} \land \theta < \theta^{max}$ do
    $s' \leftarrow repair(destroy(s))$
    if $f(s') > f(s)$ then
      $s \leftarrow s'$
    else
      $\theta \leftarrow \theta + 1$
    end
    $\delta \leftarrow \delta + 1$
  end
  if $f(s') > f(s^*)$ then $s^* \leftarrow s$;
  repeat
    $s \leftarrow perturb(s^*)$
  until $1 - f(s)/f(s^*) \le \epsilon$;
end
return $s^*$

$\triangleright$ Initialization (see §7.1)

$\triangleright$ LNS phase (see §7.1)

$\triangleright$ $f(s)$ denotes the value of the objective function for a solution $s$

$\triangleright$ Perturbation phase (see §7.2)

## 7.1 The LNS phase

LNS is a local search algorithm where the current solution is successively partially destroyed and then repaired. Our LNS uses a single destroy and repair operator. It explores the solution space without allowing any violation of the time constraints (nor other constraints). It has two stopping conditions: a maximum total number of iterations and a maximum number of iterations without improvement. Its implementation borrows ideas from the iterated local search introduced by Vansteenwegen et al. (2009) to solve the TOPTW.

Our LNS relies on the modeling of a solution to the PC-JSOCMSR as a sequence of jobs using the main resource. To reduce the complexity and increase the time efficiency of testing the feasibility of operations modifying the solution, we maintain predecessor and successor information on each resource for each scheduled job. Let $\mathcal{J}^{sch}$ the set of jobs scheduled during the planning horizon in the current solution. We associate with each job $j \in \mathcal{J}^{sch}$ its predecessor $pred_j^r$ and successor $succ_j^r$ on resource $r \in \mathcal{R}$. We introduce two dummy jobs denoted by *head* and *tail* and define them as the predecessor of the first scheduled job and successor of the last scheduled job on every resource. The processing time attributes associated with these dummy jobs
are set to 0. We associate with each scheduled job its starting time on the main resource (denoted by $start_j^0$) and an attribute (denoted by $shift_j$) representing the maximum value by which $start_j^0$ can be shifted while still maintaining time window feasibility for the subsequent jobs. We also set $start_{head}^0 = t^1$, $start_{tail}^0 = t^2$, $shift_{head} = shift_{tail} = 0$. For better readability and to ease the presentation of the algorithm, we hereafter denote by $start_j$, $end_j$, and $end_j^0$ the starting time of job j, its ending time on the secondary resource and its ending time on the main resource, respectively. Specifically, we have $start_j = start_j^0 - p_j^-$, $end_j = start_j + p_j$, and $end_j^0 = start_j^0 + p_j^0$. Note that we only work with semi-active schedules in the algorithm, meaning that a job is scheduled after the end of the previous job on resource 0 as early as possible with respect to the time windows. The attributes are computed as follows and updated throughout the algorithm:

$$ 
\begin{aligned}
start_j^0 &= \begin{cases} t_{j1}^1 + p_j^- & \text{if } pred_j^0 = \text{head} \\ \phi_j^{es}(end_{pred_j^0}^0, end_{pred_j^r_j}) + p_j^- & \text{if } pred_j^0 \neq \text{head} \end{cases} \\
shift_j &= \begin{cases} t_{jw_j}^2 - (start_j^0 + p_j^0 + p_j^+) & \text{if } succ_j^0 = \text{tail} \\ \phi_j^{ls}(start_{succ_j^0}^0 + shift_{succ_j^0}, start_{succ_j^r_j} + shift_{succ_j^r_j}) & \text{if } succ_j^0 \neq \text{tail} \end{cases}
\end{aligned}
$$

where $\phi_j^{es}(t, t')$ denotes the earliest starting time possible for j if the resources 0 and $r_j$ become available at time t and $t'$, and $\phi_j^{ls}(t, t')$ denotes the latest starting time possible for j if the resources 0 and $r_j$ should become available at time t and $t'$ (see their mathematical definition in §S.1.2 of the supplementary material).

The heuristic relies on a destroy operator that removes jobs from the current solution (i.e., from $\mathcal{J}^{sch}$). This destroy operator works with a position $l \in [1, |\mathcal{J}^{sch}|] = \{c \in \mathbb{N} : 1 \le c \le |\mathcal{J}^{sch}|\}$ in the sequence of jobs scheduled in the current solution and a value $\beta \in [\beta^{min}, \beta^{max}]$ representing the number of jobs to remove from that position. If there are less than $\beta$ jobs scheduled after position $l$, then it removes jobs from the beginning of the time horizon. The rationale behind it is that contiguous jobs and contiguous free space increase the possibility for insertion of unscheduled jobs. The value of $l$ is incremented by one every time an improving solution has not been found at the last iteration. If $l$ reaches the end of the sequence (i.e., $l = |\mathcal{J}^{sch}|$), then $l$ is reset to 1 and $\beta$ is incremented by one if it is below its maximum possible value $\beta^{max}$.

At the end of the destroy phase, we also apply with a probability of 0.5 a random perturbation to the solution. Specifically, we replace $\eta$ tasks by non-scheduled tasks from $\mathcal{J} \setminus \mathcal{J}^{sch}$. The selection among feasible replacements is performed using a uniform probability distribution.

After jobs removal and a potential perturbation, we use a repair operator to build a new solution. This operator considers two operations: replacement and insertion. A replacement consists of replacing a scheduled job by an unscheduled job, while an insertion consists of inserting an unscheduled job into the solution. Replacing job $j_1$ by job $j_2$ is feasible if the following conditions are satisfied:

$$ 
\begin{cases}
\phi_j^{es}(end_{pred_{j1}^0}^0, end_{pred_{j1}^{r_{j2}}}^0) + p_{j2}^- + p_{j2}^0 \le start_{succ_{j1}^0}^0 + shift_{succ_{j1}^0} \\
\phi_j^{es}(end_{pred_{j1}^0}^0, end_{pred_{j1}^{r_{j2}}}^0) + p_{j2} \le start_{succ_{j1}^{r_{j2}}}^{r_{j2}} + shift_{succ_{j1}^{r_{j2}}}
\end{cases}
$$

We check if the shift induced by the operation to the subsequent jobs on the main and secondary resources is feasible without any time violation. If this is the case, we set $start_{j2}^0 = \phi_j^{es}(end_{pred_{j1}^0}^0, end_{pred_{j1}^{r_{j2}}}^0)$ and update $\mathcal{J}^{sch}$ to $(\mathcal{J}^{sch} \setminus \{j_1\}) \cup \{j_2\}$. Inserting job $j_2$ after job $j_1$ is feasible if the following conditions (similar to the conditions for a replacement) are satisfied:

$$ 
\begin{cases}
\phi_j^{es}(end_{j1}^0, end_{pred_{succ_{j1}^0}^{r_{j2}}}^0) + p_{j2}^- + p_{j2}^0 \le start_{succ_{j1}^0}^0 + shift_{succ_{j1}^0} \\
\phi_j^{es}(end_{j1}^0, end_{pred_{succ_{j1}^0}^{r_{j2}}}^0) + p_{j2} \le start_{succ_{j1}^{r_{j2}}}^{r_{j2}} + shift_{succ_{j1}^{r_{j2}}}
\end{cases}
$$

If this is the case, we set $start_{j2}^0 = \phi_j^{es}(end_{j1}^0, end_{pred_{succ_{j1}^0}^{r_{j2}}}^0)$ and update $\mathcal{J}^{sch}$ to $\mathcal{J}^{sch} \cup \{j_2\}$.

We compute a score for each feasible operation. The score $w_o$ of an operation $o$ corresponding to the insertion of job $j$ after a currently scheduled job $j'$ is equal to

$$ w_o = \frac{m}{m+1} \times b_j \times \left( \frac{1}{p_j^0 + \alpha.wait^0} + \frac{1}{m(p_j + \alpha.wait)} \right) $$
with $wait_0 = \min \{ \max\{0, start_{succ_j}^0 - end_j^0\}, \max\{0, start_j^0 - end_{pred_j}^0\} \}$ and $wait = \min \{ \max\{0, start_{succ_j} - end_j\}, \max\{0, start_j - end_{pred_j}\} \}$. The score $w_o$ of an operation $o$ corresponding to the replacement of job $j'$ at its current position by job $j$ is equal to

$$w_o = \frac{m}{m+1} \times (b_j - b_{j'}) \times \left( \frac{1}{p_j^0 + \alpha.wait^0} + \frac{1}{m(p_j + \alpha.wait)} \right)$$

The score of an operation increases as the inserted job has a large profit by time unit of use of the main resource and the secondary resource. The use of the secondary resource is considered with a factor equal to $1/m$. Inserting a job that leads to an important waiting time with respect to both its predecessor and successor is penalized with a factor equal to $\alpha \in [0, 1]$. The algorithm selects the operation to apply among the $\gamma$ operations with the highest score using a uniform probability distribution. It continues as long as there exist an operation with a strictly positive score. This means that a replacement can only be applied if it improves the objective function of the solution. This avoids a cycling effect, but may restrict the search space. The random perturbation applied before the reparation phase is meant to alleviate this limitation, since non-improving replacements can be applied. Note that the removed jobs are candidates for an insertion at a different position in the solution. These jobs may end up at their initial position at the end of the repair procedure, but only if it is profitable for the objective function.

After any modification of the current solution, we recompute the solution attributes requiring an update. Specifically, if $j'$ is the first job scheduled after a sequence of consecutive jobs that has been removed or after the job that has been inserted (either by an insertion or a replacement operation), we update $start_j^0$ to $start_j^0$ and $shift_j$ to $shift_j$ as follows:

$$\overline{start}_j^0 = \begin{cases} start_j^0 & \text{if } j \in \{\text{head}, \text{tail}\} \lor j \text{ scheduled before } j' \text{ on the main resource} \\ \phi_j^{es}(\overline{end}_{pred_j}^0, \overline{end}_{pred_j}^0) & \text{otherwise} \end{cases}$$

$$\overline{shift}_j = \begin{cases} shift_j - (start_j^0 - start_j^0) & \text{if } j = j' \lor (j \neq tail \land j \text{ scheduled after } j' \text{ on the main resource}) \\ \phi_j^{ls}(start_{succ_j}^0 + shift_{succ_j}, start_{succ_j} + shift_{succ_j}) - start_j^0 & \text{otherwise} \end{cases}$$

We start by computing the attributes forward from $j'$ to tail and then backward from $pred_j^0$ to head. We stop the update in each direction as soon as the starting time of the last job considered by the update procedure has not changed on every resource.

## 7.2 The perturbation phase

To escape from a local optimum, we perturb the current best solution $s^*$ by removing $\beta^{max}$ jobs using a uniform probability distribution. Contrary to the repair operator, the jobs removed are not necessarily consecutive. Next, we replace randomly $\eta$ scheduled jobs by non-scheduled ones. We finally use the LNS repair operator described in the previous section to improve the solution. We want the resulting solution to have similarities to $s^*$, but enough differences to explore a different part of the search space in the next LNS phase.

# 8 Computational experiments

We coded all the algorithms using C++ and compiled them with GCC 9.2.0. We used the BaPCod C++ library v.053 (Sadykov and Vanderbeck, 2021) for the implementation of the BCP algorithm. For the bucket graph implementation of the labeling algorithm, path enumeration, and the separation of limited-memory R1Cs, we used the C++ code developed by Sadykov et al. (2021). Modifications to these generic codes have been required to handle specific characteristics of the PC-JSOCMSR, especially for the labeling algorithms. We used IBM CPLEX Optimizer version 12.10.0 with the default parameters as the LP solver in column generation, and as the MILP solver for the formulations presented in Section 4. The experiments were run on a 2 Deca-core Ivy-Bridge Haswell Intel Xeon E5-2680 v3 server running at 2.50 GHz. The 128 GB of available RAM was
shared between six copies of the algorithm running in parallel on the server. Each instance is solved by one
copy of the algorithm using a single thread.

We considered the benchmark instances<fnref n="1" /> used in Horn et al. (2021a). The instances are split into two sets.
The total number of jobs *n* in these instances belongs to {50, 100, 150, 200, 250, 300, 400, 500} (we excluded the
instances with 350 and 450 to reduce the number of tests). The set denoted *A* contains instances related to
the pre-runtime scheduling of avionic systems. It has the following characteristics: the length of the planning
horizon is 1000, the number of secondary resources *m* is equal to three or four, the number of jobs that can
be scheduled within their time windows is independent of the total number of jobs, “regular” jobs require the
main resource either at the beginning or at the end of their processing and for between one fifth and one eighth
of their total processing time (that is between 200 and 320 for “regular” jobs), “communication” jobs have a
processing time on the main and secondary resource equal to 40, between 20 to 30 jobs can be scheduled during
the planning horizon, and the prize of each job refers to a priority. The set denoted *P* are instances related to
particle therapy scheduling. Their characteristics are the following: the length of the planning horizon (from
100 to 800) depends on the number of jobs, the number of secondary resources *m* is equal to two or three, the
number of jobs that can be scheduled within their time windows grows approximately linearly with the total
number of jobs, the jobs require the main resource on average during the second third of their processing, the
average processing time of a job is 40, around 30% of the jobs can be scheduled during the planning horizon,
and the prize of each job is correlated to its processing time on the main resource. In all the instances of set
*A* and *P*, the secondary resources are equally distributed among the jobs. Jobs have between one and three
time windows (except the communication jobs in set *A* which can have up to six time windows). There are 30
instances in each set for each combination of *n* and *m*, which makes a total of 960 instances. We refer to the
Appendix D in Horn et al. (2021a) for details about the generation procedure.

The first aim of our computational experiments is to assess the performance of a commercial MILP solver
for solving instances of the PC-JSOCMSR using the compact and pseudo-polynomial formulations introduced
in section 4. These results are presented in §8.1. The second aim of our experiments is to assess the effectiveness
of our heuristic algorithm introduced in section 7 to build high-quality solutions to the PC-JSOCMSR. The
results and comparison with the literature are presented in §8.2. The third and most important aim of our
computational experiments is to assess the quality of our BCP algorithm as an exact solution method for the
PC-JSOCMSR. We also compare with results from the literature. We report in §S.4 of the supplementary
material the results for each combination of instance parameters, meaning that each line is the average of the
results of 30 instances (60 instances when reporting the results either according to the number of jobs or the
number of secondary resources). The detailed results for every single tested instance are available online<fnref n="2" />. In
all the tables, the CPU time is given in seconds.

8.1 Results for [F1] and [F2] with a MILP solver

We empirically tested the performance of solving formulations [F1] and [F2] with a MILP solver. We restrict
our analysis to the instances with 50 and 100 jobs. We ran the two MILP formulations with a two-hour CPU
time limit (single thread). For each formulation, Table 1 reports the number of instances with a solution proven
to be optimal (#Opt), the average CPU time in seconds (Time), the average gap (Gap) with respect to the
upper bound retrieved by the solver. For a better comparison, we also report the average gap (Gap^[F2]^) with
respect to the upper bound retrieved by the solver on the model [F2] (this bound is always lower than the
bound obtained with the model [F1]). We compute the gap with respect to a value z̄ as (z̄ - z)/z̄, where z
is the objective function value of the best integer solution returned by the solver. For unsolved instances, the
computational time is set to the time limit.

We observe that the MILP solver provides proven optimal solutions for all the instances when using model
[F2], whereas it usually fails to return optimal solutions within the time limit when using model [F1]. In the

<footnote n="1">¹Available from https://www.ac.tuwien.ac.at/research/problem-instances/ (last accessed: 2021/07/13)</footnote>

<footnote n="2">²https://www.math.u-bordeaux.fr/~afroger001/documents/Results_PCJSOCMSR_Article.tar.xz</footnote>
<table><caption>Table 1: Results for the MILP formulations [F1] and [F2] according to the instance characteristics</caption><thead><tr><th rowspan="2">Set</th><th rowspan="2">n</th><th rowspan="2">m</th><th colspan="4">[F1]</th><th colspan="4">[F2]</th></tr><tr><th>#Opt</th><th>Time</th><th>Gap</th><th>Gap<sup>[F2]</sup></th><th>#Opt</th><th>Time</th><th>Gap</th></tr></thead><tbody><tr><td rowspan="4">A</td><td rowspan="2">50</td><td>3</td><td>0/30</td><td>7200</td><td>15.2%</td><td>0.5%</td><td>30/30</td><td>6</td><td>0.0%</td></tr><tr><td>4</td><td>0/30</td><td>7200</td><td>19.2%</td><td>0.8%</td><td>30/30</td><td>5</td><td>0.0%</td></tr><tr><td rowspan="2">100</td><td>3</td><td>0/30</td><td>7200</td><td>50.7%</td><td>11.1%</td><td>30/30</td><td>12</td><td>0.0%</td></tr><tr><td>4</td><td>0/30</td><td>7200</td><td>51.5%</td><td>12.2%</td><td>30/30</td><td>11</td><td>0.0%</td></tr><tr><td rowspan="4">P</td><td rowspan="2">50</td><td>2</td><td>22/30</td><td>3782</td><td>2.8%</td><td>0.1%</td><td>30/30</td><td>0</td><td>0.0%</td></tr><tr><td>3</td><td>13/30</td><td>5615</td><td>8.3%</td><td>0.3%</td><td>30/30</td><td>0</td><td>0.0%</td></tr><tr><td rowspan="2">100</td><td>2</td><td>0/30</td><td>7200</td><td>57.6%</td><td>15.3%</td><td>30/30</td><td>57</td><td>0.0%</td></tr><tr><td>3</td><td>0/30</td><td>7200</td><td>52.3%</td><td>11.2%</td><td>30/30</td><td>45</td><td>0.0%</td></tr></tbody></table>

latter case, the gap is very high. It seems to come from a very weak upper bound resulting from solving the linear relaxation of model [F1] and, to a lesser extent, from a difficulty for finding high-quality primal solutions (cf. column Gap$^{[F2]}$). We conclude that time discretization allows us to build a model that is more interesting from a computational point of view.

We empirically tested the performance of solving [F2] for larger instances. We report the results in Table 2. The solver returns an optimal solution for all instances of set A, whereas it fails to provide an optimal solution for instances of set P within the time limit for most of the instances with 150 jobs or more. The gap for the unsolved instances of set P can still be important after two hours of computational time.

Our results suggest that instances for which the profit associated with each job is correlated with the processing time on the main resource (i.e., instances of set P) are the most difficult. These instances have jobs with a shorter duration than instances of set A and more jobs can be scheduled. Using Proposition 2, we are able for the instances of set A to detect a significant number of jobs that are dominated by other jobs (a job j dominates a job j' if it is only profitable to schedule j' when j is scheduled). We can also detect some jobs (around 5% on average) that are never profitable to schedule. We also observe that instances with the smallest number of secondary resources or a large number of tasks are more challenging.

<table><thead><tr><th>Set</th><th>n</th><th>m</th><th>#Opt</th><th>Time</th><th>Gap</th></tr></thead><tbody><tr><td rowspan="10">A</td><td rowspan="3">150</td><td>3</td><td>30/30</td><td>18</td><td>0.0%</td></tr><tr><td>4</td><td>30/30</td><td>17</td><td>0.0%</td></tr><tr><td>3</td><td>30/30</td><td>23</td><td>0.0%</td></tr><tr><td rowspan="3">200</td><td>4</td><td>30/30</td><td>25</td><td>0.0%</td></tr><tr><td>3</td><td>30/30</td><td>28</td><td>0.0%</td></tr><tr><td>4</td><td>30/30</td><td>29</td><td>0.0%</td></tr><tr><td rowspan="6">300</td><td>3</td><td>30/30</td><td>35</td><td>0.0%</td></tr><tr><td>4</td><td>30/30</td><td>36</td><td>0.0%</td></tr><tr><td>3</td><td>30/30</td><td>48</td><td>0.0%</td></tr><tr><td>4</td><td>30/30</td><td>49</td><td>0.0%</td></tr><tr><td>3</td><td>30/30</td><td>59</td><td>0.0%</td></tr><tr><td>4</td><td>30/30</td><td>63</td><td>0.0%</td></tr></tbody></table>

<table><caption>Table 2: Additional results for the MILP formulation [F2] according to the instance characteristics</caption><thead><tr><th>Set</th><th>n</th><th>m</th><th>#Opt</th><th>Time</th><th>Gap</th></tr></thead><tbody><tr><td rowspan="4">P</td><td rowspan="2">150</td><td>2</td><td>0/30</td><td>7200</td><td>2.3%</td></tr><tr><td>3</td><td>24/30</td><td>3367</td><td>0.1%</td></tr><tr><td rowspan="2">200</td><td>2</td><td>0/30</td><td>7200</td><td>3.8%</td></tr><tr><td>3</td><td>1/30</td><td>6995</td><td>0.9%</td></tr><tr><td rowspan="2">250</td><td>2</td><td>0/30</td><td>7200</td><td>4.1%</td></tr><tr><td>3</td><td>2/30</td><td>6986</td><td>0.9%</td></tr></tbody></table>

## 8.2 Results for the heuristic algorithm

For these experiments, we considered all the instances and performed 10 runs for each instance. We ran the heuristic algorithm with the following parameters set after performing preliminary experiments: $\delta^{max} = 200000$, $\theta^{max} = 2000$, $\epsilon = 0.02$, $\beta^{min} = 4$, $\beta^{max} = \max\{4, 0.2|\mathcal{J}^{sch}|\}$, $\alpha = 0.25$, $\eta = 2$, $\gamma = 5$.

For comparison purposes, we have computed the best known solution (BKS) for each instance from the results of Maschler and Raidl (2018), Horn et al. (2021a), our heuristic algorithm, the MILP formulation [F2]
and the BCP algorithm applied to solve [F3] (the two latter methods are not considered for instances of set P with more than 250 jobs). As a result, all BKS for set A and almost all BKS for instances of set P with less than 250 jobs are optimal.

Table 3 reports the average over 10 runs of the CPU time in seconds (Time) and the average gap (Gap) with respect to the BKSs. We compute the gap with respect to a value $\bar{z}$ as $(\bar{z} - z)/\bar{z}$, where z is the mean (Mean), minimum (Worst) or maximum (Best) objective function value of the best solution computed by the heuristic over the 10 runs. We report results after 20000 and 200000 iterations. To judge the efficiency of the heuristic to find good-quality solutions in a short computational time, we also tested solving formulation [F2] in a heuristic way. Specifically, for each instance we set a time limit equal to the average computational time of the heuristic for 200000 iterations. We report the gap with respect to the BKSs ([F2] - Heur.). We also report in the two last columns the gap with respect to the BKS for the GVNS introduced by Maschler and Raidl (2018) (GVNS) and the TDC with A* search and cost-based filtering described in Horn et al. (2021a) (TDC-A*-f). These two last methods were run for 900 seconds.

First, if we look at the results for 20000 iterations, we can see that the heuristic usually finds very quickly good solutions, especially for the instances of set A. Nonetheless, we see that the remaining iterations are necessary to find better solutions. The computational time remains reasonably low. The heuristic exhibits a stable behavior for instances of set A, while the difference between the best and the worst solution found over the 10 runs is on average 0.5% for instances of set P. On the tested instances, results are slightly better when $m$ is larger. We clearly see that our heuristic outperforms existing methods in terms of solution quality. It computes a solution with an objective strictly better than or identical to the best reported in the literature for all but one of the 960 instances of sets A and P if we consider the best result over the 10 runs for each instance (and all but 61 out of the 960 instances if we consider the worst result). The improvement increases with the number of jobs and are primarily for instances of set P. We also observe that the results found by the MILP solver with model [F2] are better for the instances of set A. It finds the optimal solutions in less than one minute on average and clearly outperforms our heuristic. For the instances of set P, the solver is not competitive for finding good primal solutions for instances with 150 jobs and more.

To have a better overview, Figure 4 shows the number of instances whose computed solution has a gap with respect to the BKS smaller than a certain value. As a whole, we see that our heuristic outperforms the other methods.

<figure><img src="image_7.png" /><figcaption>Figure 4: Performance profiles of several heuristic methods on the 960 instances of set A and set P</figcaption></figure>

We tested the LNS without the ILS framework (it corresponds to setting $\theta^{max} = \delta^{max}$). Results are worse on average. We observe a 0.5% difference between the average gap with respect to the BKSs with and without ILS for the instances of set P with 200 and 400 jobs.

We finally tested removing the random perturbation performed before calling the repair operator. We observe worst (or almost equal) results on average for the instances of set A and for the instances of set P with a number of jobs less than or equal to 200. For the instances of set P with a number of jobs larger than 200
<table><caption>Table 3: Comparison of the results obtained by different heuristic methods</caption><thead><tr><th rowspan="2">Set</th><th rowspan="2">n</th><th colspan="3">Heuristic / δ = 20000</th><th colspan="3">Heuristic / δ = 200000</th><th rowspan="2">[F2] -<br>Heur.<br>Gap</th><th colspan="2">Literature</th></tr><tr><th>Worst</th><th>Mean</th><th>Best</th><th>Worst</th><th>Mean</th><th>Best</th><th>GVNS<br>Gap</th><th>TDC-A*-f<br>Gap</th></tr></thead><tbody><tr><td rowspan="10">A</td><td>50</td><td>0.16%</td><td>0.03%</td><td>0.00%</td><td>2</td><td>0.01%</td><td>0.00%</td><td>0.00%</td><td>0.00%</td><td>0.00%</td></tr><tr><td>100</td><td>0.13%</td><td>0.03%</td><td>0.00%</td><td>7</td><td>0.00%</td><td>0.00%</td><td>0.00%</td><td>0.00%</td><td>0.02%</td></tr><tr><td>150</td><td>0.22%</td><td>0.05%</td><td>0.00%</td><td>12</td><td>0.01%</td><td>0.00%</td><td>0.00%</td><td>0.00%</td><td>0.04%</td></tr><tr><td>200</td><td>0.15%</td><td>0.05%</td><td>0.00%</td><td>16</td><td>0.03%</td><td>0.01%</td><td>0.00%</td><td>0.00%</td><td>0.06%</td></tr><tr><td>250</td><td>0.10%</td><td>0.04%</td><td>0.00%</td><td>21</td><td>0.01%</td><td>0.00%</td><td>0.00%</td><td>0.00%</td><td>0.10%</td></tr><tr><td>300</td><td>0.15%</td><td>0.06%</td><td>0.01%</td><td>25</td><td>0.02%</td><td>0.01%</td><td>0.00%</td><td>0.00%</td><td>0.23%</td></tr><tr><td>400</td><td>0.22%</td><td>0.06%</td><td>0.00%</td><td>34</td><td>0.03%</td><td>0.01%</td><td>0.00%</td><td>0.00%</td><td>0.37%</td></tr><tr><td>500</td><td>0.24%</td><td>0.07%</td><td>0.01%</td><td>42</td><td>0.06%</td><td>0.02%</td><td>0.01%</td><td>0.00%</td><td>0.26%</td></tr><tr><td>All</td><td></td><td>0.17%</td><td>0.05%</td><td>0.00%</td><td>20</td><td>0.02%</td><td>0.01%</td><td>0.00%</td><td>0.00%</td><td>0.13%</td></tr><tr><td rowspan="10">P</td><td>50</td><td>0.36%</td><td>0.06%</td><td>0.00%</td><td>2</td><td>0.00%</td><td>0.00%</td><td>0.00%</td><td>16</td><td>0.00%</td></tr><tr><td>100</td><td>1.11%</td><td>0.42%</td><td>0.02%</td><td>6</td><td>0.28%</td><td>0.07%</td><td>0.00%</td><td>60</td><td>0.02%</td></tr><tr><td>150</td><td>1.69%</td><td>1.05%</td><td>0.47%</td><td>13</td><td>0.76%</td><td>0.38%</td><td>0.13%</td><td>129</td><td>0.61%</td></tr><tr><td>200</td><td>2.09%</td><td>1.54%</td><td>1.04%</td><td>23</td><td>1.27%</td><td>0.86%</td><td>0.48%</td><td>235</td><td>1.19%</td></tr><tr><td>250</td><td>1.70%</td><td>1.24%</td><td>0.75%</td><td>36</td><td>0.91%</td><td>0.59%</td><td>0.25%</td><td>367</td><td>2.81%</td></tr><tr><td>300</td><td>1.36%</td><td>0.91%</td><td>0.48%</td><td>51</td><td>0.61%</td><td>0.31%</td><td>0.01%</td><td>517</td><td>8.16%</td></tr><tr><td>400</td><td>1.23%</td><td>0.81%</td><td>0.40%</td><td>84</td><td>0.49%</td><td>0.26%</td><td>0.00%</td><td>860</td><td>14.64%</td></tr><tr><td>500</td><td>1.18%</td><td>0.77%</td><td>0.40%</td><td>125</td><td>0.48%</td><td>0.24%</td><td>0.00%</td><td>1282</td><td>17.86%</td></tr><tr><td>All</td><td></td><td>1.33%</td><td>0.84%</td><td>0.44%</td><td>42</td><td>0.60%</td><td>0.34%</td><td>0.11%</td><td>433</td><td>5.66%</td></tr><tr><td rowspan="4">Set</td><td>m</td><td colspan="3"></td><td colspan="3"></td><td colspan="2"></td></tr><tr><td>3</td><td>0.18%</td><td>0.05%</td><td>0.01%</td><td>20</td><td>0.02%</td><td>0.01%</td><td>0.00%</td><td>196</td><td>0.00%</td></tr><tr><td>4</td><td>0.16%</td><td>0.05%</td><td>0.00%</td><td>20</td><td>0.02%</td><td>0.01%</td><td>0.00%</td><td>201</td><td>0.00%</td></tr><tr><td>2</td><td>1.54%</td><td>1.02%</td><td>0.56%</td><td>47</td><td>0.74%</td><td>0.43%</td><td>0.14%</td><td>483</td><td>6.17%</td></tr><tr><td></td><td>3</td><td>1.11%</td><td>0.67%</td><td>0.32%</td><td>38</td><td>0.46%</td><td>0.25%</td><td>0.08%</td><td>383</td><td>5.16%</td></tr></tbody></table>

instances, we observe slightly better results with a difference between the average gap with respect to the BKS approximately equal to 0.05% for the solutions computed by the heuristic with and without this perturbation.

## 8.3 Results for solving [F3] with a BCP algorithm

We tested the performance of solving [F3] with our BCP algorithm on the instances with up to 250 jobs (set P) and up to 500 jobs (set A). We set a two-hour CPU time limit (single thread). We used as input the value of the best solution found by our heuristic algorithm with δ = 200000 described in Section 7 minus one (the total computational time we report includes the time to compute this solution). We tested the following BCP configurations:

* **G**: pricing problem defined on graph *G*

* **G + R1Cs**: configuration *G* + separation of R1Cs and use of the limited memory version of these cuts

* **G'**: pricing problem defined on graph *G'*

* **G' + R1Cs**: configuration *G'* + separation of R1Cs and use of the limited memory version of these cuts

* **G' + R1Cs + AR**: configuration *G'* + R1Cs + reduction of the number of arcs in *G'* according to the dominance rules introduced in §3.

Table 4 gives the comparison of four BCP configurations. It shows the average values for the root gap (Root Gap), the root CPU time (Root Time), the gap for the computed solution with respect to the best dual bound (Gap), the number of branch-and-bound nodes (TN), the pricing problem solution time (SP Time), and the total CPU time (Time). It also reports the number of instances optimally solved within the time limit (#Opt). We compute the root gap considering the best (primal) solution found before starting branching and the linear relaxation value. For unsolved instances, the computational time is set to the time limit.

Solutions to instances of set *P* contain a lot more jobs than those of set *A*, which explained why they are easier to solve with a BCP algorithm. We see that results obtained for the BCP configurations with a
<table><caption>Table 4: Results for the BCP configurations according to the pricing problem definition and the use of R1Cs</caption><thead><tr><th rowspan="2">Set</th><th colspan="6">G</th><th colspan="6">G + R1Cs</th></tr><tr><th>Root Gap</th><th>Root Time</th><th>Gap</th><th>#Opt</th><th>TN</th><th>SP Time</th><th>Time</th><th>Root Gap</th><th>Root Time</th><th>Gap</th><th>#Opt</th><th>TN</th><th>SP Time</th><th>Time</th></tr></thead><tbody><tr><td>A</td><td>0.10%</td><td>19</td><td>0.01%</td><td>475/480</td><td>26</td><td>55</td><td>320</td><td>0.09%</td><td>20</td><td>0.00%</td><td>480/480</td><td>1</td><td>17</td><td>219</td></tr><tr><td>P</td><td>1.05%</td><td>79</td><td>0.53%</td><td>180/300</td><td>71</td><td>1745</td><td>3253</td><td>0.77%</td><td>656</td><td>0.39%</td><td>188/300</td><td>7</td><td>2511</td><td>2968</td></tr></tbody></table><table><thead><tr><th rowspan="2">Set</th><th colspan="6">G'</th><th colspan="6">G' + R1Cs</th></tr><tr><th>Root Gap</th><th>Root Time</th><th>Gap</th><th>#Opt</th><th>TN</th><th>SP Time</th><th>Time</th><th>Root Gap</th><th>Root Time</th><th>Gap</th><th>#Opt</th><th>TN</th><th>SP Time</th><th>Time</th></tr></thead><tbody><tr><td>A</td><td>0.09%</td><td>289</td><td>0.00%</td><td>478/480</td><td>9</td><td>267</td><td>514</td><td>0.09%</td><td>292</td><td>0.00%</td><td>480/480</td><td>1</td><td>254</td><td>493</td></tr><tr><td>P</td><td>1.01%</td><td>244</td><td>0.43%</td><td>197/300</td><td>49</td><td>2058</td><td>2931</td><td>0.74%</td><td>993</td><td>0.23%</td><td>230/300</td><td>4</td><td>1913</td><td>2289</td></tr></tbody></table>

pricing problem defined on $G$ are better than for those with a pricing problem defined on $G'$ for instances of set $A$, whereas it is the opposite for instances of set $P$ (with a greater difference). When using graph $G$, dominance performed at node $v_*$ reduces the total number of generated labels, which reduces the pricing problem solution time. However, for instances of set $P$, we remove a larger proportion of (bucket) arcs with the reduced cost-based filtering procedures when using $G'$. Moreover, the set of schedules considered with graph $G'$ is smaller than the one considered with graph $G$. Indeed, it is forbidden by construction of $G'$ to have the same job using the main resource more than once in a row in a particular time window. This does not seem to have a significant impact on the root gap for the tested instances, but this makes finding integer feasible solution easier when starting branching. When using R1Cs with graph $G$, we observe that the number of non-dominated labels when solving the pricing problem explodes. This can be explained by larger arc memory sets than with graph $G'$. This also confirms that solutions to the pricing problem contain more jobs with multiple occurrences for BCP configurations $G$ than for BCP configurations $G'$. Our results also show that using R1Cs improves the convergence of the BCP algorithm as it is able to compute more proven optimal solutions within the same computational time. The root gap decreases by approximately 0.3% for instances of set $P$ and the branch-and-bound tree is smaller.

Table 5 shows a comparison of the results obtained when reducing the number of arcs in $G'$. Although we only see a slight improvement in the computational time (due to a faster resolution of the pricing problem), the arc reduction procedure nevertheless allows to solve to optimality four additional instances.

<table><caption>Table 5: Results for the BCP algorithm depending on using the arc reduction results</caption><thead><tr><th rowspan="2">Set</th><th colspan="6">G' + R1Cs</th><th colspan="6">G' + R1Cs + AR</th></tr><tr><th>Root Gap</th><th>Root Time</th><th>Gap</th><th>#Opt</th><th>TN</th><th>SP Time</th><th>Time</th><th>Root Gap</th><th>Root Time</th><th>Gap</th><th>#Opt</th><th>TN</th><th>SP Time</th><th>Time</th></tr></thead><tbody><tr><td>A</td><td>0.09%</td><td>292</td><td>0.00%</td><td>480/480</td><td>1</td><td>254</td><td>493</td><td>0.09%</td><td>170</td><td>0.00%</td><td>480/480</td><td>1</td><td>151</td><td>370</td></tr><tr><td>P</td><td>0.74%</td><td>993</td><td>0.23%</td><td>230/300</td><td>4</td><td>1913</td><td>2289</td><td>0.74%</td><td>907</td><td>0.22%</td><td>234/300</td><td>4</td><td>1831</td><td>2204</td></tr></tbody></table>

Table 6 compares solving [F3] with our BCP algorithm to solving model [F2] with a MILP solver. Our BCP configurations are not competitive to solve the easiest instances of set $A$. For the more difficult and more challenging instances of set $P$, they outperform the MILP solver, which demonstrates their interest. Instances with more than 200 jobs are challenging to solve as there are still a lot of them currently open. Instances with more secondary resources are also more difficult to solve by the BCP algorithm, which is explained by weaker dominance rules.
<table><caption>Table 6: Results for the best BCP configurations and for formulation [F2] depending on the instance characteristics</caption><thead><tr><th rowspan="2">Set</th><th rowspan="2">n</th><th colspan="3">G + R1Cs</th><th colspan="3">G' + R1Cs + AR</th><th colspan="3">[F2]</th></tr><tr><th>Opt</th><th>Gap</th><th>Time</th><th>Opt</th><th>Gap</th><th>Time</th><th>Opt</th><th>Gap</th><th>Time</th></tr></thead><tbody><tr><td rowspan="8">A</td><td>50</td><td>60/60</td><td>0.00%</td><td>24</td><td>60/60</td><td>0.00%</td><td>24</td><td>60/60</td><td>0.00%</td><td>6</td></tr><tr><td>100</td><td>60/60</td><td>0.00%</td><td>71</td><td>60/60</td><td>0.00%</td><td>78</td><td>60/60</td><td>0.00%</td><td>12</td></tr><tr><td>150</td><td>60/60</td><td>0.00%</td><td>122</td><td>60/60</td><td>0.00%</td><td>139</td><td>60/60</td><td>0.00%</td><td>18</td></tr><tr><td>200</td><td>60/60</td><td>0.00%</td><td>173</td><td>60/60</td><td>0.00%</td><td>206</td><td>60/60</td><td>0.00%</td><td>24</td></tr><tr><td>250</td><td>60/60</td><td>0.00%</td><td>222</td><td>60/60</td><td>0.00%</td><td>293</td><td>60/60</td><td>0.00%</td><td>29</td></tr><tr><td>300</td><td>60/60</td><td>0.00%</td><td>278</td><td>60/60</td><td>0.00%</td><td>406</td><td>60/60</td><td>0.00%</td><td>36</td></tr><tr><td>400</td><td>60/60</td><td>0.00%</td><td>382</td><td>60/60</td><td>0.00%</td><td>701</td><td>60/60</td><td>0.00%</td><td>49</td></tr><tr><td>500</td><td>60/60</td><td>0.00%</td><td>483</td><td>60/60</td><td>0.00%</td><td>1109</td><td>60/60</td><td>0.00%</td><td>61</td></tr><tr><td>All</td><td>480/480</td><td>0.00%</td><td>219</td><td>480/480</td><td>0.00%</td><td>370</td><td>480/480</td><td>0.00%</td><td>29</td></tr><tr><td rowspan="5">P</td><td>50</td><td>60/60</td><td>0.00%</td><td>16</td><td>60/60</td><td>0.00%</td><td>16</td><td>60/60</td><td>0.00%</td><td>0</td></tr><tr><td>100</td><td>60/60</td><td>0.00%</td><td>63</td><td>60/60</td><td>0.00%</td><td>64</td><td>60/60</td><td>0.00%</td><td>51</td></tr><tr><td>150</td><td>58/60</td><td>0.02%</td><td>577</td><td>60/60</td><td>0.00%</td><td>280</td><td>24/60</td><td>1.20%</td><td>5284</td></tr><tr><td>200</td><td>10/60</td><td>0.69%</td><td>6662</td><td>44/60</td><td>0.20%</td><td>3405</td><td>1/60</td><td>2.35%</td><td>7097</td></tr><tr><td>250</td><td>0/60</td><td>1.24%</td><td>7522</td><td>10/60</td><td>0.90%</td><td>7253</td><td>2/60</td><td>2.50%</td><td>7093</td></tr><tr><td>All</td><td>188/300</td><td>0.39%</td><td>2968</td><td>234/300</td><td>0.22%</td><td>2204</td><td>147/300</td><td>1.21%</td><td>3267</td></tr><tr><td rowspan="3">Set A</td><td>m</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>3</td><td>240/240</td><td>0.00%</td><td>209</td><td>240/240</td><td>0.00%</td><td>273</td><td>240/240</td><td>0.00%</td><td>29</td></tr><tr><td>4</td><td>240/240</td><td>0.00%</td><td>230</td><td>240/240</td><td>0.00%</td><td>467</td><td>240/240</td><td>0.00%</td><td>29</td></tr><tr><td rowspan="2">Set P</td><td>m</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>2</td><td>97/150</td><td>0.46%</td><td>2819</td><td>130/150</td><td>0.18%</td><td>1726</td><td>60/150</td><td>2.04%</td><td>4331</td></tr><tr><td></td><td>3</td><td>91/150</td><td>0.32%</td><td>3117</td><td>104/150</td><td>0.26%</td><td>2681</td><td>87/150</td><td>0.38%</td><td>3479</td></tr></tbody></table>

Last, to provide a global overview and finer comparison of solution time, Figure 5 shows the number of instances optimally solved according to the solution time for the different exact methods presented in this paper. We come to the same conclusion as before. The MILP solver performs better for instances of set A and that the BCP configuration G+R1Cs+AR performs better for instances of set P.
<figure><img src="image_8.png" /><figcaption>Figure 5: Performance charts of the best exact algorithms on the instances of set A (top) and set P (bottom)</figcaption></figure>

9 Conclusion

We have proposed a new pseudo-polynomial time-indexed formulation ([F2]) for the PC-JSOCMSR. Computational experiments have shown that a MILP solver finds significantly better solutions with this formulation in comparison with an existing compact order-based formulation.

As an exact solution method, we have introduced a BCP algorithm to solve an extended MILP formulation of the problem based on the generation of schedules that may contain several occurrences of the same job. We have presented two alternative ways of defining the pricing problem as a resource-constrained shortest path problem. We have also presented dominance rules to reduce the search space. While it only has a slight impact on the efficiency of the method, these rules are general enough to be interesting for other algorithms. Our BCP algorithm is able to compute proven optimal solutions for more instances in comparison with solving [F2] with a MILP solver. Using R1Cs with limited memory to strengthen the linear relaxation provides better results. We have optimally solved instances with up to 250 jobs (500 jobs on a particular set of instances), while in the literature only a few instances with 50 jobs were optimally solved.

We have developed an ILS to heuristically solve the PC-JSOCMSR. Our experiments have shown this heuristic algorithm gives equal or better results than existing heuristic methods of the literature for more than 93% of the tested instances in the worst case of our 10 runs (99% in the best case) and the gap between the solution value returned by our heuristic and that of the literature is on average for the most difficult set of instances between approximately 1.5% and 2%.

As a perspective, refining relaxations of the time-indexed formulation [F2] may be a path worth investigating for solving exactly larger instances. A possible extension of this work concerns the case of multiple main resources with every secondary resource being associated with a dedicated main resource (e.g., the case of patients that
can choose to have their cancer treatment in one location out of several).

## Acknowledgments

The experiments presented in this paper were carried out using the Federative Platform for Research in Computer Science and Mathematics (PlaFRIM), created under the Inria PlaFRIM development action with support from Institut Polytechnique de Bordeaux, Laboratoire Bordelais de Recherche en Informatique, Institut de Mathématiques de Bordeaux, Conseil Régional d'Aquitaine, Université de Bordeaux, Centre National de la Recherche Scientifique, and Agence Nationale de la Recherche in accordance with the "Programme d'Investissements d'Avenir" (see www.plafrim.fr/en/home).

## References

* Allahverdi, A. (2015). The third comprehensive survey on scheduling problems with setup times/costs. *European Journal of Operational Research*, 246(2):345-378.

* Castro, M. P., Cire, A. A., and Beck, J. C. (2022). Decision diagrams for discrete optimization: A survey of recent advances. arXiv preprint arXiv:2201.11536, https://arxiv.org/abs/2201.11536.

* Dror, M. (1994). Note on the complexity of the shortest path models for column generation in vrptw. *Operations Research*, 42(5):977-978.

* Gunawan, A., Lau, H. C., and Vansteenwegen, P. (2016). Orienteering Problem: A survey of recent variants, solution approaches and applications. *European Journal of Operational Research*, 255(2):315-332.

* Horn, M., Maschler, J., Raidl, G. R., and Rönnberg, E. (2021a). A*-based construction of decision diagrams for a prize-collecting scheduling problem. *Computers and Operations Research*, 126:105125.

* Horn, M., Raidl, G. R., and Rönnberg, E. (2018). An A* algorithm for solving a prize-collecting sequencing problem with one common and multiple secondary resources and time windows. In *Proceedings of PATAT*, pages 235-256.

* Horn, M., Raidl, G. R., and Rönnberg, E. (2021b). A* search for prize-collecting job sequencing with one common and multiple secondary resources. *Annals of Operations Research*, 302:477-505.

* Irnich, S., Desaulniers, G., Desrosiers, J., and Hadjar, A. (2010). Path-Reduced Costs for Eliminating Arcs in Routing and Scheduling. *INFORMS Journal on Computing*, 22(2):297-313.

* Koulamas, C. and Kyparisis, G. J. (2008). Single-machine scheduling problems with past-sequence-dependent setup times. *European Journal of Operational Research*, 187(3):1045-1049.

* Maschler, J. and Raidl, G. R. (2018). Multivalued decision diagrams for a prize-collecting sequencing problem. In *Proceedings of the 12th International Conference of the Practice and Theory of Automated Timetabling, PATAT 2018*. pp. 375-397.

* Maschler, J. and Raidl, G. R. (2020). Particle therapy patient scheduling with limited starting time variations of daily treatments. *International Transactions in Operational Research*, 27(1):458-479.

* Pecin, D., Contardo, C., Desaulniers, G., and Uchoa, E. (2017a). New Enhancements for the Exact Solution of the Vehicle Routing Problem with Time Windows. *INFORMS Journal on Computing*, 29(3):489-502.

* Pecin, D., Pessoa, A., Poggi, M., and Uchoa, E. (2017b). Improved branch-cut-and-price for capacitated vehicle routing. *Mathematical Programming Computation*, 9(1):61-100.
Pessoa, A., Sadykov, R., and Uchoa, E. (2021). Solving bin packing problems using vrpsolver models. *Operations Research Forum*, 2(2):20.

Pessoa, A., Sadykov, R., Uchoa, E., and Vanderbeck, F. (2018). Automation and combination of linear-programming based stabilization techniques in column generation. *INFORMS Journal on Computing*, 30(2):339-360.

Sadykov, R., Uchoa, E., and Pessoa, A. (2021). A bucket graph-based labeling algorithm with application to vehicle routing. *Transportation Science*, 55(1):4-28.

Sadykov, R. and Vanderbeck, F. (2021). BaPCod - a generic branch-and-price code. Technical report, Inria Bordeaux Sud-Ouest.

van der Veen, J. A. A., Woeginger, G. J., and Zhang, S. (1998). Sequencing jobs that require common resources on a single machine: A solvable case of the TSP. *Mathematical Programming*, 82(1-2):235-254.

Vansteenwegen, P., Souffriau, W., and Oudheusden, D. V. (2011). The orienteering problem: A survey. *European Journal of Operational Research*, 209(1):1-10.

Vansteenwegen, P., Souffriau, W., Vanden Berghe, G., and Van Oudheusden, D. (2009). Iterated local search for the team orienteering problem with time windows. *Computers & Operations Research*, 36(12):3281-3290.
# S.1 Proofs

## S.1.1 Proof of Proposition 1

*Proof.* Let $j, j' \in \mathcal{J}$ such that $j \neq j'$, $k \in [1, w_j]$, and $k' \in [1, w_{j'}]$, and assume that $t_{jk}^1 \le t_{j'k'}^1$, $t_{jk}^2 \le t_{j'k'}^2$, $p_j^- \ge p_{j'}^-$, and $p_j^+ \le p_{j'}^+$.

Let $S_{j' \to j}$ be a schedule with job $j'$ using the main resource just before job $j$ uses it. Let $start_j$ and $start_{j'}$ be the starting time of job $j$ and $j'$ in this schedule. As $t_{jk}^1 \le t_{j'k'}^1$, we can always assume that $start_j = start_{j'} + p_{j'}$. In $S_{j' \to j}$, the main resource is either available or used by the two jobs from time $start_{j'} + p_{j'}$ to time $start_j + p_j^- + p_j^0$. The secondary resource $r_j$ is used by the two jobs from time $start_{j'}$ to time $start_j + p_j$.

Let $S_{j' \to j}$ be the schedule obtained from $S_{j' \to j}$ by switching the processing order of job $j$ and $j'$ and with $j$ starting at time $start_{j'}$ and $j'$ starting at time $start_{j'} + p_j$. As $t_{jk}^1 \le t_{j'k'}^1$, job $j$ can start at time $start_{j'}$. As $t_{jk}^2 \le t_{j'k'}^2$, job $j'$ can start at time $start_{j'} + p_j$. This means that in $S_{j' \to j}$ the secondary resource $r_j$ is still used by the two jobs from time $start_{j'}$ to time $start_{j'} + p_j + p_{j'} = start_j + p_j$. In $S_{j' \to j}$, the main resource is used from time $start_{j'} + p_j$ to time $start_j + p_{j'} + p_{j'}^0$. As $p_j^- \ge p_{j'}^-$, we have $start_{j'} + p_j^- \ge start_{j'} + p_{j'}$. As $p_j^+ \le p_{j'}^+$, we also have $start_j + p_{j'}^- + p_{j'}^0 \le start_{j'} + p_j + p_{j'}^- + p_{j'}^0$. This means that the main resource is required in $S_{j' \to j}$ during a sub-interval of the interval during which it is required in $S_{j' \to j}$.

We can therefore always convert a schedule in which $j'$ uses the main resource just before $j$ to a schedule in which $j$ uses the main resource just before $j'$. □

## S.1.2 Proof of Proposition 2

We first introduce the functions $\phi_{jk}^{es}$ and $\phi_{jk}^{ls}$ defined for every $j \in \mathcal{J}$ and $k \in [1, w_j]$:

* $\phi_{jk}^{es}(t, t')$ denotes the earliest starting time possible for a job $j$ in time window $T_{jk}$ if the resources 0 and $r_j$ become available at time $t$ and $t' (\infty)$ if not existing):

$$ \phi_{jk}^{es}(t, t') = \begin{cases} t_{jk}^1 & \text{if } \max\{t - p_j^-, t'\} \le t_{jk}^1 \\ \max\{t - p_j^-, t'\} & \text{if } t_{jk}^1 < \max\{t - p_j^-, t'\} \le t_{jk}^2 - p_j \\ \infty & \text{if } t_{jk}^2 - p_j < \max\{t - p_j^-, t'\} \end{cases} $$

* $\phi_{jk}^{ls}(t, t')$ denotes the latest starting time possible for a job $j$ in time window $T_{jk}$ if the resources 0 and $r_j$, after their use by $j$, should become available at time $t$ and $t' (-\infty)$ if not existing):

$$ \phi_{jk}^{ls}(t, t') = \begin{cases} t_{jk}^2 - p_j & \text{if } \min\{t - p_j^-, p_{j'}^0, t' - p_j\} \ge t_{jk}^2 - p_j \\ \max\{t - p_j^-, t'\} & \text{if } t_{jk}^1 \le \min\{t - p_j^-, p_{j'}^0, t' - p_j\} < t_{jk}^2 - p_j \\ -\infty & \text{if } \min\{t - p_j^-, p_{j'}^0, t' - p_j\} < t_{jk}^1 \end{cases} $$

From the functions $\phi_{jk}^{es}$ and $\phi_{jk}^{ls}$ defined for every $j \in \mathcal{J}$ and $k \in [1, w_j]$, we mathematically define the functions $\phi_j^{es}$ and $\phi_j^{ls}$ as follows:

$$ \phi_j^{es}(t, t') = \min_{k \in [1, w_j]} \{\phi_{jk}^{es}(t, t')\} $$

$$ \phi_j^{ls}(t, t') = \max_{k \in [1, w_j]} \{\phi_{jk}^{ls}(t, t')\} $$

We then mathematically translate the content of Proposition 2.

* Let $j' \in \mathcal{J}$, $k' \in [1, w_{j'}]$, and $j \in \mathcal{J} \setminus \{j'\}$. Job $j$ cannot use the main resource before job $j'$ scheduled during time window $T_{j'k'}$ if:

$$ \left\{ \begin{array}{ll} t_{j1}^1 + p_j > t_{j'k'}^2 - p_{j'} & \text{if } r_j = r_{j'} \\ t_{j1}^1 + p_j^- + p_j^0 > t_{j'k'}^2 - p_{j'} + p_{j'}^- & \text{if } r_j \neq r_{j'} \end{array} \right. \qquad (26) $$
Conditions (26) state that scheduling *j* at its earliest starting time implies that *j'* cannot be scheduled during time window *T*<sub>*j*′*k*′</sub>.

* Let $j' \in \mathcal{J}$, $k' \in [1, w_{j'}]$, $j'' \in \mathcal{J} \setminus \{j', j''\}$, $k'' \in [1, w_{j''}]$, and $j \in \mathcal{J} \setminus \{j', j''\}$. Job $j$ cannot use the main resource before this resource is used by job $j'$ scheduled during time window $T_{j'k'}$ and then by job $j''$ scheduled during time window $T_{j''k''}$ if:

$$
\left\{
\begin{array}{ll}
t_{j1}^1 + p_j > \phi_{j'k'}^{ls}(t_{j''k''}^2 - p_{j''} + p_{j''} - t_{j''k''}^2 - p_{j''}) & \text{if } r_j = r_{j'} = r_{j''} \\
t_{j1}^1 + p_j > \phi_{j'k'}^{ls}(t_{j''k''}^2 - p_{j''} + p_{j''} - t_{j''k''}^2) & \text{if } r_j = r_{j'} \wedge r_{j''} \neq r_j \\
t_{j1}^1 + p_j^- + p_j^0 > \phi_{j'k'}^{ls}(t_{j''k''}^2 - p_{j''} + p_{j''} - t_{j''k''}^2 + p_{j''} - t_{j1}^1 + p_j > t_{j''k''}^2 - p_{j''}) & \text{if } r_j = r_{j''} \wedge r_{j'} \neq r_j \\
t_{j1}^1 + p_j^- + p_j^0 > \phi_{j'k'}^{ls}(t_{j''k''}^2 - p_{j''} + p_{j''} - t_{j''k''}^2 - p_{j''}) + p_{j''} - t_{j1}^1 + p_j^0 > \phi_{j'k'}^{ls}(t_{j''k''}^2 - p_{j''} + p_{j''} - t_{j''k''}^2 + p_{j''} - t_{j1}^1 + p_j) & \text{if } r_j \neq r_{j'} \wedge r_{j'} = r_{j''} \\
t_{j1}^1 + p_j^- + p_j^0 > \phi_{j'k'}^{ls}(t_{j''k''}^2 - p_{j''} + p_{j''} - t_{j''k''}^2 - p_{j''}) + p_{j''} - t_{j1}^1 + p_j^0 > \phi_{j'k'}^{ls}(t_{j''k''}^2 - p_{j''} + p_{j''} - t_{j''k''}^2 - p_{j''}) + p_{j''} - t_{j1}^1 & \text{if } r_j \neq r_{j'} \wedge r_j \neq r_{j''} \wedge r_{j'} \neq r_{j''}
\end{array}
\right.
\tag{27}
$$

Conditions (27) state that scheduling *j* at its earliest starting time implies that *j'* and *j''* cannot use the main resource in this order after *j*.

* Let $j' \in \mathcal{J}$, $k' \in [1, w_{j'}]$, and $j \in \mathcal{J} \setminus \{j'\}$. Job $j$ cannot use the main resource after job $j'$ scheduled during time window $T_{j'k'}$ if:

$$
\left\{
\begin{array}{ll}
t_{j'k'}^1 + p_{j'} > t_{j'w_j}^2 - p_j & \text{if } r_j = r_{j'} \\
t_{j'k'}^1 + p_{j'}^- + p_{j'}^0 > t_{j'w_j}^2 - p_j + p_{j'}^- & \text{if } r_j \neq r_{j'}
\end{array}
\right.
\qquad (28)
$$

Conditions (28) state that scheduling *j* at its latest starting time implies that *j'* cannot be scheduled
during time window $T_{j'k'}$.

* Let $j' \in \mathcal{J}$, $k' \in [1, w_{j'}]$, $j'' \in \mathcal{J} \setminus \{j'\}$, $k'' \in [1, w_{j''}]$, and $j \in \mathcal{J} \setminus \{j', j''\}$. Job $j$ cannot use the main resource after this resource is used by job $j'$ scheduled during time window $T_{j'k'}$ and then by job $j''$ scheduled during time window $T_{j''k''}$ if:

$$
\left\{
\begin{array}{ll}
\phi_{j''k''}^{es}(t_{j'k'}^1 + p_{j'}^- + p_{j'}^+, t_{j'k'}^1 + p_{j'}^-) + p_{j''} > t_{j'w_j}^2 - p_j & \text{if } r_j = r_{j'} = r_{j''} \\
\phi_{j''k''}^{es}(t_{j'k'}^1 + p_{j'}^- + p_{j'}^+, 0) + p_{j''} + p_{j''}^0 > t_{j'w_j}^2 - p_j + p_{j'}^- \lor t_{j'k'}^1 + p_{j'}^- > t_{j'w_j}^2 - p_j & \text{if } r_j = r_{j'} \land r_{j''} \neq r_j \\
\phi_{j''k''}^{es}(t_{j'k'}^1 + p_{j'}^- + p_{j'}^+, 0) + p_{j''} > t_{j'w_j}^2 - p_j & \text{if } r_j = r_{j''} \land r_{j'} \neq r_j \\
\phi_{j''k''}^{es}(t_{j'k'}^1 + p_{j'}^- + p_{j'}^+, t_{j'k'}^1 + p_{j'}^-) + p_{j''} + p_{j''}^0 > t_{j'w_j}^2 - p_j + p_{j'}^- & \text{if } r_j \neq r_{j'} \land r_{j'} = r_{j''} \\
\phi_{j''k''}^{es}(t_{j'k'}^1 + p_{j'}^- + p_{j'}^+, 0) + p_{j''} + p_{j''}^0 > t_{j'w_j}^2 - p_j + p_{j'}^- & \text{if } r_j \neq r_{j'} \land r_j \neq r_{j''} \land r_{j'} \neq r_{j''}
\end{array}
\right.
\tag{29}
$$

Conditions (29) state that scheduling *j* at its latest starting time implies that *j'* and *j''* cannot use the main resource in this order before *j*.

* Let $j' \in \mathcal{J}$, $k' \in [1, w_{j'}]$, and $j \in \mathcal{J} \setminus \{j'\}$. Job $j$ can use the main resource just before job $j'$ scheduled during $T_{j'k'}$ without preventing $j'$ to start its execution at time $t_{j'k'}^1$ and imposing restriction on the processing of subsequent jobs:

$$
\left\{
\begin{array}{ll}
t_{j1}^1 + p_j \le t_{j'k'}^1 & \text{if } r_{j'} = r_j \\
t_{j1}^1 + p_j^- + p_j^0 \le t_{j'k'}^1 + p_{j'}^- \wedge t_{j1}^1 + p_j \le t_{j'k'}^1 + p_{j'}^- + p_{j'}^0 - \hat{p}_{r_j} & \text{if } r_{j'} \ne r_j
\end{array}
\right.
\qquad (30)
$$

Conditions (30) state that there exists a starting time for *j* such that *j'* can start at any time during time
window $T_{j'k'}$ and that the next job using the secondary resource $r_j$ (possibly $j'$) can start at any time. If
$r_j = r_{j'}$, only the secondary resource needs consideration. If $r_j \neq r_{j'}$, observe that the main resource (after
its use by $j'$) can never be available for subsequent jobs before time $t_{j'k'}^1 + p_{j'}^-$ and $p_{j'}^0$. The earliest possible
starting time of a task $j'''$ such that $r_{j''} = r_j$ is therefore equal to $t_{j'k'}^1 + p_{j'}^- + p_{j'}^0 - p_{j'''} \ge t_{j'k'}^1 + p_{j'}^- + p_{j'}^0 - \hat{p}_{r_j}$.

* Let $j' \in \mathcal{J}$, $k' \in [1, w_{j'}]$, and $j \in \mathcal{J} \setminus \{j'\}$. Job $j$ can use the main resource just after job $j'$ scheduled during $T_{j'k'}$ without preventing $j'$ to start its execution at time $t_{j'k'}^2 - p_{j'}$ and imposing restriction on the
processing of precedent jobs if:

$$
\left\{
\begin{array}{ll}
t_{j'k'}^2 + p_j \le t_{jw_j}^2 & \text{if } r_{j'} = r_j \\
t_{j'k'}^2 - p_{j'}^+ + p_j^0 + p_j^+ \le t_{jw_j}^2 \wedge t_{j'k'}^2 - p_{j'} + p_{j'}^- + \hat{p}_{r_j}^+ + p_j \le t_{jw_j}^2 & \text{if } r_{j'} \ne r_j
\end{array}
\right.
\qquad (31)
$$

Conditions (31) state that there exists a starting time for *j* such that *j'* can start at any time during time window *T*<sub>j'*k'*</sub> and that the previous job using the secondary resource *r*<sub>j</sub> (possibly *j'*) can start at any time. As before, only the secondary resource needs consideration if *r*<sub>j</sub> = *r*<sub>j'*</sub>. If *r*<sub>j</sub> ≠ *r*<sub>j'*</sub>, observe that the main resource (before its use by *j'* ) should be available not later than time *t*<sub>j'*k'</sub><sup>2</sup> − *p*<sub>j'</sub> + *p*<sub>j'</sub><sup>−</sup>. The latest possible completion time of a task *j'''* such that *r*<sub>j'''</sub> = *r*<sub>j</sub> is therefore equal to *t*<sub>j'*k'</sub><sup>2</sup> − *p*<sub>j'</sub> + *p*<sub>j'</sub><sup>−</sup> + *p*<sub>r<sub>j</sub></sub><sup>+</sup>.

* Let *j'* ∈ 𝒩, *k'* ∈ [1, *w*<sub>j'</sub>], *j''* ∈ 𝒩\{*j'*}, *k''* ∈ [1, *w*<sub>j''</sub>], and *j* ∈ 𝒩\{*j'* , *j''*}. Job *j* can use the main resource just after *j'* and just before *j''* without preventing *j''* to start its execution at time *t*<sub>j''k''</sub><sup>1</sup> and *j'* to end at time *t*<sub>j'*k'</sub><sup>2</sup> and imposing restriction on the processing of precedent and subsequent jobs if:

$$
\left\{
\begin{array}{ll}
\phi_j^{es}(t_{j'k'}^2 - p_{j'}^+, t_{j'k'}^2) + p_j \le t_{j''k''}^1 & \text{if } r_j = r_{j''} = r_{j'} \\
\phi_j^{es}(t_{j'k'}^2 - p_{j'}^+, t_{j'k'}^2 - p_{j'}^+ - p_{j'}^0 + \hat{p}_{r_j}^+) + p_j \le t_{j''k''}^1 & \text{if } r_j = r_{j''} \wedge r_j \ne r_{j'} \\
(\phi_j^{es}(t_{j'k'}^2 - p_{j'}^+, t_{j'k'}^2) + p_{j'}^- + p_j^0 \le t_{j''k''}^1 + p_{j''}^- + p_{j''}^0 - \hat{p}_{r_j}^-) & \text{if } r_j = r_{j'} \wedge r_j \ne r_{j''} \\
(\phi_j^{es}(t_{j'k'}^2 - p_{j'}^+, t_{j'k'}^2 - p_{j'}^+ - p_{j'}^0 + \hat{p}_{r_j}^+) + p_{j'}^- + p_j^0 \le t_{j''k''}^1 + p_{j''}^- + p_{j''}^0 - \hat{p}_{r_j}^-) & \text{if } r_j \notin \{r_{j'}, r_{j''}\}
\end{array}
\right.
\tag{32}
$$

Conditions (32) state that there exists a starting time for *j* such that *j'* and *j''* can start at any time during time window *T*<sub>j'*k'*</sub> and *T*<sub>j''k''</sub> and such that the previous and next job using the secondary resource *r*<sub>j</sub> (possibly *j'* and *j''*) can start at any time. We use the same arguments as above to derive the sufficient condition. Since *j* is scheduled between *j'* and *j''*, the use of function $\phi_j^{es}(\cdot, \cdot)$ is required.

In Lemma 1, we identify when a job *j* can always be scheduled during time window *T*<sub>*jk*</sub> in the place of a
job *j'* scheduled during time window *T*<sub>*j'*k'*</sub>.

**Lemma 1.** Let $j, j' \in \mathcal{J}$ such that $j \neq j' \wedge r_j = r_{j'}, k \in [1, w_j]$, and $k' \in [1, w_{j'}]$. The following conditions are sufficient to ensure that job *j* can always be scheduled during time window $T_{jk}$ in the place of job *j'* scheduled during time window $T_{j'*k'*}$.

$$
\exists d \in [\max\{0, p_{j'}^{-} - p_{j}^{-}\}, \min\{p_{j'}^{-} - p_{j}^{-} + p_{j'}^{0} - p_{j}^{0}, p_{j'}^{-} - p_{j}^{-}\}] : t_{jk}^{1} \le t_{j'*k'*}^{1} + d \wedge t_{jk}^{2} + p_{j'}^{-} - p_{j}^{-} - d \ge t_{j'*k'*}^{2} \quad (33)
$$

Proof. Let $S_{j'}$ be a schedule with job $j'$ and $start_{j'}$ the starting time of this job. In $S_{j'}$, the main resource is used by $j'$ from time $start_{j'} + p_{j'}^{-}$ to time $start_{j'} + p_{j'}^{-} + p_{j'}^{0}$. The secondary resource $r_j$ is used by $j'$ from time $start_{j'}$ to time $start_{j'} + p_{j'}^{-}$.

Let $S_j$ be the schedule obtained from $S_{j'}$ by replacing $j'$ by $j$ and setting $start_j = start_{j'} + d$. As $start_{j'} \ge t_{j'*k'*}^1$ and $t_{jk}^1 \le t_{j'*k'*}^1 + d$, we have $start_j \ge t_{jk}^1$. As $start_{j'} + p_{j'}^- \le t_{j'*k'*}^2$ and $t_{jk}^2 + p_{j'}^- - p_j - d \ge t_{j'*k'*}^2$, we have $start_j + p_j^- \le t_{jk}^2$.

In $S_j$, the main resource is used from time $start_j + p_j^-$ to time $start_j + p_j^+ + p_j^0$. Since $d \ge p_j^- - p_j^-$, we have $start_j + p_j^- \ge start_{j'} + p_{j'}^-$. Since $d \le p_j^- - p_j^- + p_{j'}^0 - p_j^0$, we have $start_j + p_j^- + p_j^0 \le start_{j'} + p_{j'}^- + p_{j'}^0$. This means that the main resource is required in $S_j$ during a sub-interval of the interval during which it is required in $S_{j'}$.

In $S_j$, the secondary resource is used from time $start_j$ to time $start_j + p_j$. Since $d \ge 0$, we have $start_j \ge start_{j'}$. Since $d \le p_{j'} - p_j$, we have $start_j + p_j \le start_{j'} + p_{j'}$. This means that the secondary resource is required in $S_j$ during a sub-interval of the interval during which it is required in $S_{j'}$.

We can therefore always replace job *j'* scheduled in time window *T*<sub>*j'*k'*</sub> by job *j*.

<figure><img src="image_9.png" /></figure>
We finally reformulate Proposition 2 using the conditions previously introduced.

**Proposition 2.** Assume that a job $j' \in \mathcal{J}$ is scheduled during time window $T_{j'k'}$ with $k' \in [1, w_{j'}]$. There exists an optimal schedule in which:

1. Job $j'$ does not use the main resource just after its use by another job $j'' \in \mathcal{J}\{j'\}$ during time window $T_{j''k''}$ with $k'' \in [1, w_{j''}]$ (resp. $j'$ is not the first job using the main resource) if there exists a job $j \in \mathcal{J}\{j', j''\}$ (resp. $j \in \mathcal{J}\{j'\}$) such that $(j'', k'', j', k', j)$ satisfies (27) and (29) (resp. $(j', k', j)$ satisfies (28)) and that satisfies at least one of the two following conditions:

a) $r_j = r_{j'}, b_j \ge b_{j'}$ and $\exists k \in [1, w_j]$ such that $(j, j', k, k')$ satisfies (33)

b) $(j'', k'', j)$ satisfies (31) and $(j', k', j)$ satisfies (30) (resp. $(j', k', j)$ satisfies (30))

2. Job $j'$ does not use the main resource just before its use by another job $j'' \in \mathcal{J}\{j'\}$ during time window $T_{j''k''}$ with $k'' \in [1, w_{j''}]$ (resp. $j'$ is not the last job using the main resource) if there exists a job $j \in \mathcal{J}\{j', j''\}$ (resp. $j \in \mathcal{J}\{j'\}$) such that $(j', k', j'', k'', j)$ satisfies (27) and (29) (resp. $(j', k', j)$ satisfies (26)) and that satisfies at least one of the two following conditions:

a) $r_j = r_{j'}, b_j \ge b_{j'}$ and $\exists k \in [1, w_j]$ such that $(j, j', k, k')$ satisfies (33)

b) $(j', k', j)$ satisfies (31) and $(j'', k'', j)$ satisfies (30) (resp. $(j', k', j)$ satisfies (31))

*Proof.* Case 1: Let $S$ be a schedule where the main resource is consecutively used by $j''$ and $j'$ (in that order). This means that $j$ is not scheduled in $S$.

a) It is always possible to build a more profitable schedule $S'$ that is identical to $S$ with the exception of either job $j$ scheduled such that it uses the main resource between $j''$ and $j'$ (resp. before $j'$), or job $j$ scheduled in the place of $j'$.

b) It is always possible to build a more profitable schedule $S'$ that is identical to $S$ with the exception of job $j$ scheduled such that it uses the main resource between $j''$ and $j'$ (resp. before $j'$)

Case 2: Let now $S$ be a schedule where the main resource is consecutively used by $j'$ and $j''$ (in that order). This means that $j$ is not scheduled in $S$.

a) It is always possible to build a more profitable schedule $S'$ that is identical to $S$ with the exception of either job $j$ scheduled such that it uses the main resource between $j'$ and $j''$ (resp. after $j'$), or job $j$ scheduled in the place of $j'$.

b) It is always possible to build a more profitable schedule $S'$ that is identical to $S$ with the exception of job $j$ scheduled such that it uses the main resource between $j'$ and $j''$ (resp. after $j'$)

□
## S.2 Labeling algorithms

<caption>**Algorithm 3: runMonodirectionalAlgorithm (dir)**</caption>

input : dir ∈ {forward, backward}
output: a set of extended non-dominated labels
if dir = forward then else
L ← createInitialLabel (dir)
N ← {L}, E ← ∅
while N ≠ ∅ do
L ← getNextLabelToExtend (N, dir) ▷ See Function 6
N ← N \ {L}
for L' ∈ E such that vL' = vL and QL' 0 < QL' 0 do
| if Dominates (L', L, dir) then go to line 3 ▷ See Function 8
end
E ← E ∪ {L}
if dir = forward then
| Γ ← Γ+(vL, QL' 0 ) ▷ Set of outgoing arcs from vL at time QL' 0 in the bucket implementation of the graph
else
| Γ ← Γ−(vL, QL' 0 ) ▷ Set of incoming arcs at vL at time QL' 0 in the bucket implementation of the graph
end
for a ∈ Γ do
| if Extend (L, a, L', dir) and not IsBorderCrossed (L', dir) then ▷ See Functions 7 and 9
| for L'' ∈ N ∪ E such that vL'' = vL' and QL'' 0 = QL' 0 do
| | if Dominates (L'', L', dir) then continue the loop of line 3 with the next arc
| end
| for L'' ∈ N such that vL'' = vL' and QL'' 0 = QL' 0 do
| | if Dominates (L', L'', dir) then N ← N \ {L''}
| end
N ← N ∪ {L'}
end
end
end
return E
Algorithm 4: runBidirectionalAlgorithm()

<figure><img src="image_10.png" /></figure>

Function 5: CreateInitialLabel

<figure><img src="image_11.png" /></figure>

Function 6: getNextLabelToExtend

<figure><img src="image_12.png" /></figure>
## Function 7: Extend

**input**: L is a label, a is an arc, dir ∈ {forward, backward}, L' is the result of the extension of L along a in direction dir

**output**: true if the extension succeeds, false otherwise

**Function Extend(L, a, L', dir):**

if dir = forward then

$$v^{L'} \leftarrow \text{head of } a \\ \Delta \leftarrow \max \{0, l_{vL',0} - Q_0^L - q_{a,0}, \max_{r \in \mathcal{R}} \{Q_r^L + q_{a,r} - u_{vL',r}\}\}$$

$$Q_0^{L'} \leftarrow Q_0^L + q_{a,r} + \Delta$$

if $Q_0^{L'} > u_{vL',0}$ then return false

for $r \in \mathcal{R}$ do

$$Q_r^{L'} \leftarrow \max \{l_{a,r}, Q_r^L + q_{a,r} - \Delta\}$$

end

else

$$v^{L'} \leftarrow \text{queue of } a \\ \Delta \leftarrow \max \{0, Q_0^L - q_{a,0} - u_{vL',0}, \max_{r \in \mathcal{R}} \{l_{vL',r} - (Q_r^L - q_{a,r})\}\}$$

$$Q_0^{L'} \leftarrow Q_0^L - q_{a,r} - \Delta$$

if $Q_0^{L'} < l_{vL',0}$ then return false

for $r \in \mathcal{R}$ do

$$Q_r^{L'} \leftarrow \min \{u_{vL',r}, Q_r^L - q_{a,r} + \Delta\}$$

end

end

if $v^{L'} \in \bigcup_{j \in J} \bar{V}_j$ then return false

$$c^{L'} \leftarrow c^L + c_a \qquad \triangleright \text{Only for the algorithm of §5.2}$$

$$\text{if } \exists j \in \mathcal{J} \text{ such that } v^{L'} \in \bar{V}_j \text{ then } J^{L'} \leftarrow J^L \cup \{j\} \text{ else } J^{L'} \leftarrow J^L \qquad \triangleright \text{Only for the algorithm of §5.2}$$

$$\bar{c}^{L'} \leftarrow \bar{c}^L + \bar{c}_a \qquad \triangleright \text{Only for the algorithm of §6.2}$$

$$\text{for } g \in \mathcal{G} \text{ do} \qquad \triangleright \text{Only for the algorithm of §6.2}$$

$$\text{if } a \notin AM(g) \text{ then } E_g^{L'} \leftarrow E_g^L \text{ else } E_g^{L'} \leftarrow 0$$

if $v \in V_{j(g)}$ then

$$E_g^{L'} \leftarrow E_g^{L'} + 0.5$$

if $E_g^{L'} = 1$ then $E_g^{L'} \leftarrow 0, \bar{c}^{L'} \leftarrow \bar{c}^{L'} - \sigma_g$

end

end

return true
## Function 8: Dominates

input : L and L' are two labels, dir ∈ {forward, backward}

output: true if L dominates L', false otherwise

Function Dominates(L, L', dir):

if $v^L \neq v^{L'}$ then return false

if $J^L \not\subseteq J^{L'}$ then return false

if $c^L < c^{L'}$ then return false

if $\bar{c}^L < \bar{c}^{L'} + \sum_{g \in G: E_g^L > E_g^{L'}} \sigma_g$ then return false

if dir = forward then

if $Q_0^L > Q_0^{L'}$ then return false

for $r \in \mathcal{R}$ do

if $Q_r^L - (Q_0^{L'} - Q_0^L) > Q_r^{L'}$ then return false

end

else

if $Q_0^L < Q_0^{L'}$ then return false

for $r \in \mathcal{R}$ do

if $Q_r^L + (Q_0^L - Q_0^{L'}) < Q_r^{L'}$ then return false

end

end

return true

▷ Only for the algorithm of §5.2

▷ Only for the algorithm of §5.2

▷ Only for the algorithm of §6.2

## Function 9: IsBorderCrossed

input : L is a label, dir ∈ {forward, backward}

output: true if the border in direction dir is crossed by L, false otherwise

Function IsBorderCrossed(L, dir):

if dir = forward then

return $Q_0^L > t^*$

else

return $Q_0^L \le t^*$

end

## Function 10: ConcatenationCost

input : $\vec{L}$ is a forward label, a an arc, and $\vec{L}$ a backward label

output: $-\infty$ if the two labels cannot be concatenated,
the cost of the solution resulting from the concatenation otherwise

Function ConcatenationCost($\vec{L}$, a, $\vec{L}$):

if $Q_0^{\vec{L}} + q_{a,0} > Q_0^{\vec{L}}$ then return $-\infty$

for $r \in \mathcal{R}$ do

if $Q_r^{\vec{L}} + q_{a,r} - (Q_0^{\vec{L}} - (Q_0^{\vec{L}} + q_{a,0})) > Q_r^{\vec{L}}$ then return $-\infty$

end

if $J^{\vec{L}} \cap J^{\vec{L}} \neq \emptyset$ then return $-\infty$

return $c^{\vec{L}} + c_a + c^{\vec{L}}$

return $\bar{c}^{\vec{L}} + \bar{c}_a + \bar{c}^{\vec{L}} - \sum_{g \in G: E_g^{\vec{L}} + E_g^{\vec{L}} \ge 1} \sigma_g$

▷ Only for the algorithm of §5.2

▷ Only for the algorithm of §6.2
## S.3 Example of the bucket graph implementation

In the forward (resp. backward) bucket graph implementation, the arcs incoming to $v_{sink}$ (resp. going out of $v_{source}$) are represented as dotted arcs to make the graphs easier to read.

A label $L'$, obtained by extension from label $L$ along an arc $a = ((v, t), (v', t'))$ of the bucket graph implementation of $G$ (or $G'$) is not necessarily contained in the bucket $(v', t')$. It is contained in the bucket $(v', Q_0^{L'})$.

<figure><img src="image_13.png" /><figcaption>Figure 6: The forward bucket graph implementation of the graph $G$ displayed in Figure 1</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 7: The backward bucket graph implementation of the graph $G$ displayed in Figure 2</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 8: The forward bucket graph implementation of $G'$ displayed in Figure 3</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 9: The backward bucket graph implementation of $G'$ displayed in Figure 3</figcaption></figure>
## S.4 Detailed computational results

<table><caption>Table 7: Detailed computational results for the heuristic methods</caption><thead><tr><th rowspan="2">Set</th><th rowspan="2">n</th><th rowspan="2">m</th><th colspan="4">Heuristic / δ = 20000</th><th colspan="4">Heuristic / δ = 200000</th><th rowspan="2">[F2] -<br>Heur.<br>Gap</th><th colspan="2">Literature</th></tr><tr><th>Worst</th><th>Mean</th><th>Best</th><th>Time</th><th>Worst</th><th>Mean</th><th>Best</th><th>Time</th><th>GVNS<br>Gap</th><th>TDC-A*-f<br>Gap</th></tr></thead><tbody><tr><td rowspan="12">A</td><td rowspan="3">50</td><td>3</td><td>0.24%</td><td>0.04%</td><td>0.00%</td><td>2</td><td>0.00%</td><td>0.00%</td><td>0.00%</td><td>21</td><td>0.00%</td><td>0.00%</td><td>0.00%</td></tr><tr><td>4</td><td>0.08%</td><td>0.03%</td><td>0.00%</td><td>2</td><td>0.01%</td><td>0.00%</td><td>0.00%</td><td>22</td><td>0.00%</td><td>0.00%</td><td>0.06%</td></tr><tr><td>100</td><td>3</td><td>0.14%</td><td>0.03%</td><td>0.00%</td><td>7</td><td>0.01%</td><td>0.00%</td><td>0.00%</td><td>66</td><td>0.00%</td><td>0.00%</td><td>0.03%</td></tr><tr><td rowspan="3">150</td><td>4</td><td>0.11%</td><td>0.03%</td><td>0.00%</td><td>7</td><td>0.00%</td><td>0.00%</td><td>0.00%</td><td>70</td><td>0.00%</td><td>0.00%</td><td>0.01%</td></tr><tr><td>3</td><td>0.16%</td><td>0.05%</td><td>0.00%</td><td>11</td><td>0.01%</td><td>0.00%</td><td>0.00%</td><td>113</td><td>0.00%</td><td>0.05%</td><td>0.06%</td></tr><tr><td>4</td><td>0.27%</td><td>0.05%</td><td>0.01%</td><td>12</td><td>0.02%</td><td>0.01%</td><td>0.00%</td><td>119</td><td>0.00%</td><td>0.03%</td><td>0.01%</td></tr><tr><td rowspan="3">200</td><td>3</td><td>0.20%</td><td>0.07%</td><td>0.00%</td><td>16</td><td>0.04%</td><td>0.01%</td><td>0.00%</td><td>161</td><td>0.00%</td><td>0.07%</td><td>0.10%</td></tr><tr><td>4</td><td>0.10%</td><td>0.04%</td><td>0.00%</td><td>16</td><td>0.03%</td><td>0.01%</td><td>0.00%</td><td>163</td><td>0.00%</td><td>0.05%</td><td>0.05%</td></tr><tr><td>250</td><td>3</td><td>0.12%</td><td>0.04%</td><td>0.00%</td><td>21</td><td>0.00%</td><td>0.00%</td><td>0.00%</td><td>206</td><td>0.00%</td><td>0.03%</td><td>0.08%</td></tr><tr><td rowspan="6">P</td><td rowspan="3">100</td><td>3</td><td>0.08%</td><td>0.03%</td><td>0.00%</td><td>21</td><td>0.02%</td><td>0.00%</td><td>0.00%</td><td>211</td><td>0.00%</td><td>0.16%</td><td>0.17%</td></tr><tr><td>2</td><td>0.97%</td><td>0.35%</td><td>0.00%</td><td>6</td><td>0.22%</td><td>0.06%</td><td>0.00%</td><td>59</td><td>0.03%</td><td>0.59%</td><td>0.24%</td></tr><tr><td>3</td><td>1.20%</td><td>0.48%</td><td>0.03%</td><td>6</td><td>0.33%</td><td>0.08%</td><td>0.00%</td><td>61</td><td>0.01%</td><td>0.83%</td><td>0.33%</td></tr><tr><td rowspan="3">150</td><td>2</td><td>1.83%</td><td>1.10%</td><td>0.46%</td><td>13</td><td>0.86%</td><td>0.40%</td><td>0.12%</td><td>132</td><td>0.86%</td><td>2.13%</td><td>0.99%</td></tr><tr><td>3</td><td>1.51%</td><td>0.96%</td><td>0.46%</td><td>13</td><td>0.64%</td><td>0.36%</td><td>0.13%</td><td>126</td><td>0.36%</td><td>1.72%</td><td>1.63%</td></tr><tr><td>200</td><td>2</td><td>2.52%</td><td>1.90%</td><td>1.32%</td><td>25</td><td>1.55%</td><td>1.07%</td><td>0.59%</td><td>253</td><td>1.57%</td><td>3.11%</td><td>2.59%</td></tr><tr><td rowspan="6">P</td><td rowspan="3">250</td><td>3</td><td>1.61%</td><td>1.16%</td><td>0.75%</td><td>21</td><td>0.96%</td><td>0.64%</td><td>0.35%</td><td>218</td><td>0.81%</td><td>2.41%</td><td>2.60%</td></tr><tr><td>2</td><td>2.18%</td><td>1.65%</td><td>1.07%</td><td>40</td><td>1.24%</td><td>0.83%</td><td>0.37%</td><td>405</td><td>3.33%</td><td>3.26%</td><td>2.31%</td></tr><tr><td>3</td><td>1.20%</td><td>0.82%</td><td>0.42%</td><td>32</td><td>0.58%</td><td>0.35%</td><td>0.12%</td><td>328</td><td>2.29%</td><td>2.20%</td><td>2.90%</td></tr><tr><td rowspan="3">300</td><td>2</td><td>1.67%</td><td>1.18%</td><td>0.63%</td><td>57</td><td>0.79%</td><td>0.42%</td><td>0.01%</td><td>587</td><td>11.30%</td><td>3.07%</td><td>2.17%</td></tr><tr><td>3</td><td>1.04%</td><td>0.64%</td><td>0.32%</td><td>44</td><td>0.42%</td><td>0.21%</td><td>0.00%</td><td>448</td><td>5.01%</td><td>2.07%</td><td>3.62%</td></tr><tr><td>400</td><td>2</td><td>1.51%</td><td>1.03%</td><td>0.52%</td><td>95</td><td>0.63%</td><td>0.34%</td><td>0.00%</td><td>975</td><td>14.69%</td><td>3.88%</td><td>3.29%</td></tr><tr><td rowspan="3">500</td><td>3</td><td>0.95%</td><td>0.60%</td><td>0.28%</td><td>73</td><td>0.35%</td><td>0.18%</td><td>0.00%</td><td>744</td><td>14.59%</td><td>3.08%</td><td>4.59%</td></tr><tr><td>2</td><td>1.35%</td><td>0.91%</td><td>0.50%</td><td>139</td><td>0.59%</td><td>0.30%</td><td>0.00%</td><td>1439</td><td>17.54%</td><td>4.62%</td><td>3.55%</td></tr><tr><td>3</td><td>0.98%</td><td>0.62%</td><td>0.30%</td><td>110</td><td>0.37%</td><td>0.18%</td><td>0.00%</td><td>1125</td><td>18.18%</td><td>3.80%</td><td>5.17%</td></tr></tbody></table>
<table><caption>Table 8: Detailed computational results for BCP configuration G</caption><thead><tr><th>Set</th><th>n</th><th>m</th><th>Root Gap</th><th>Root Time</th><th>Gap</th><th>#Opt</th><th>TN</th><th>SP Time</th><th>Time</th></tr></thead><tbody><tr><td rowspan="12">A</td><td rowspan="3">50</td><td>3</td><td>0.15%</td><td>2</td><td>0.03%</td><td>29</td><td>91</td><td>105</td><td>263</td></tr><tr><td>4</td><td>0.09%</td><td>1</td><td>0.00%</td><td>30</td><td>1</td><td>1</td><td>24</td></tr><tr><td>100</td><td>3</td><td>0.11%</td><td>3</td><td>0.03%</td><td>29</td><td>83</td><td>100</td><td>308</td></tr><tr><td rowspan="3">150</td><td>3</td><td>0.08%</td><td>4</td><td>0.00%</td><td>30</td><td>2</td><td>3</td><td>74</td></tr><tr><td>4</td><td>0.12%</td><td>5</td><td>0.00%</td><td>30</td><td>6</td><td>10</td><td>135</td></tr><tr><td>200</td><td>3</td><td>0.09%</td><td>7</td><td>0.01%</td><td>29</td><td>48</td><td>116</td><td>366</td></tr><tr><td rowspan="3">250</td><td>3</td><td>0.09%</td><td>7</td><td>0.00%</td><td>30</td><td>54</td><td>74</td><td>400</td></tr><tr><td>4</td><td>0.12%</td><td>14</td><td>0.00%</td><td>30</td><td>41</td><td>82</td><td>351</td></tr><tr><td>4</td><td>0.08%</td><td>7</td><td>0.00%</td><td>30</td><td>1</td><td>5</td><td>214</td></tr><tr><td rowspan="3">300</td><td>3</td><td>0.08%</td><td>19</td><td>0.00%</td><td>30</td><td>1</td><td>17</td><td>230</td></tr><tr><td>4</td><td>0.08%</td><td>12</td><td>0.00%</td><td>30</td><td>2</td><td>17</td><td>276</td></tr><tr><td>4</td><td>0.08%</td><td>30</td><td>0.00%</td><td>30</td><td>1</td><td>27</td><td>286</td></tr><tr><td rowspan="3">400</td><td>3</td><td>0.09%</td><td>20</td><td>0.01%</td><td>29</td><td>38</td><td>90</td><td>611</td></tr><tr><td>4</td><td>0.08%</td><td>51</td><td>0.00%</td><td>30</td><td>1</td><td>47</td><td>393</td></tr><tr><td>500</td><td>3</td><td>0.09%</td><td>21</td><td>0.00%</td><td>30</td><td>1</td><td>16</td><td>437</td></tr><tr><td rowspan="12">P</td><td rowspan="3">50</td><td>2</td><td>0.83%</td><td>0</td><td>0.00%</td><td>30</td><td>1</td><td>0</td><td>16</td></tr><tr><td>3</td><td>0.72%</td><td>0</td><td>0.00%</td><td>30</td><td>1</td><td>0</td><td>17</td></tr><tr><td>100</td><td>2</td><td>0.44%</td><td>2</td><td>0.00%</td><td>30</td><td>1</td><td>1</td><td>61</td></tr><tr><td rowspan="3">150</td><td>2</td><td>0.47%</td><td>5</td><td>0.00%</td><td>30</td><td>1</td><td>3</td><td>65</td></tr><tr><td>3</td><td>1.27%</td><td>35</td><td>0.03%</td><td>29</td><td>88</td><td>368</td><td>1185</td></tr><tr><td>200</td><td>2</td><td>0.75%</td><td>63</td><td>0.08%</td><td>28</td><td>60</td><td>571</td><td>1292</td></tr><tr><td rowspan="3">250</td><td>2</td><td>1.84%</td><td>64</td><td>1.52%</td><td>1</td><td>268</td><td>2991</td><td>7356</td></tr><tr><td>3</td><td>0.99%</td><td>187</td><td>0.70%</td><td>2</td><td>116</td><td>4443</td><td>7213</td></tr><tr><td>2</td><td>2.19%</td><td>112</td><td>2.04%</td><td>0</td><td>125</td><td>3459</td><td>7610</td></tr><tr><td rowspan="3">3</td><td>2</td><td>1.02%</td><td>319</td><td>0.93%</td><td>0</td><td>46</td><td>5612</td><td>7717</td></tr></tbody></table>
<table><caption>Table 9: Detailed computational results for BCP configuration G + R1Cs</caption><thead><tr><th>Set</th><th>n</th><th>m</th><th>Root Gap</th><th>Root Time</th><th>Gap</th><th>#Opt</th><th>TN</th><th>SP Time</th><th>Time</th></tr></thead><tbody><tr><td rowspan="12">A</td><td rowspan="3">50</td><td>3</td><td>0.13%</td><td>3</td><td>0.00%</td><td>30</td><td>1</td><td>1</td><td>24</td></tr><tr><td>4</td><td>0.09%</td><td>1</td><td>0.00%</td><td>30</td><td>1</td><td>1</td><td>24</td></tr><tr><td>100</td><td>3</td><td>0.11%</td><td>4</td><td>0.00%</td><td>30</td><td>1</td><td>2</td><td>69</td></tr><tr><td rowspan="3">150</td><td>3</td><td>0.08%</td><td>4</td><td>0.00%</td><td>30</td><td>1</td><td>3</td><td>73</td></tr><tr><td>4</td><td>0.08%</td><td>5</td><td>0.00%</td><td>30</td><td>1</td><td>3</td><td>118</td></tr><tr><td>200</td><td>3</td><td>0.09%</td><td>7</td><td>0.00%</td><td>30</td><td>1</td><td>6</td><td>127</td></tr><tr><td rowspan="3">250</td><td>3</td><td>0.10%</td><td>14</td><td>0.00%</td><td>30</td><td>1</td><td>9</td><td>177</td></tr><tr><td>4</td><td>0.08%</td><td>7</td><td>0.00%</td><td>30</td><td>1</td><td>5</td><td>214</td></tr><tr><td>300</td><td>3</td><td>0.08%</td><td>19</td><td>0.00%</td><td>30</td><td>1</td><td>17</td><td>230</td></tr><tr><td rowspan="3">400</td><td>3</td><td>0.08%</td><td>16</td><td>0.00%</td><td>30</td><td>1</td><td>12</td><td>270</td></tr><tr><td>4</td><td>0.08%</td><td>30</td><td>0.00%</td><td>30</td><td>1</td><td>27</td><td>286</td></tr><tr><td>500</td><td>3</td><td>0.08%</td><td>29</td><td>0.00%</td><td>30</td><td>1</td><td>24</td><td>370</td></tr><tr><td rowspan="3">P</td><td>3</td><td>0.08%</td><td>51</td><td>0.00%</td><td>30</td><td>1</td><td>47</td><td>393</td></tr><tr><td>4</td><td>0.09%</td><td>21</td><td>0.00%</td><td>30</td><td>1</td><td>16</td><td>437</td></tr><tr><td>50</td><td>2</td><td>0.09%</td><td>103</td><td>0.00%</td><td>30</td><td>1</td><td>91</td><td>528</td></tr><tr><td rowspan="12">P</td><td rowspan="3">100</td><td>2</td><td>0.83%</td><td>0</td><td>0.00%</td><td>30</td><td>1</td><td>0</td><td>16</td></tr><tr><td>3</td><td>0.72%</td><td>0</td><td>0.00%</td><td>30</td><td>1</td><td>0</td><td>17</td></tr><tr><td>150</td><td>2</td><td>0.38%</td><td>2</td><td>0.00%</td><td>30</td><td>1</td><td>1</td><td>61</td></tr><tr><td rowspan="3">200</td><td>2</td><td>0.35%</td><td>4</td><td>0.00%</td><td>30</td><td>1</td><td>2</td><td>65</td></tr><tr><td>3</td><td>0.37%</td><td>69</td><td>0.00%</td><td>30</td><td>1</td><td>29</td><td>202</td></tr><tr><td>250</td><td>2</td><td>0.40%</td><td>345</td><td>0.04%</td><td>28</td><td>3</td><td>725</td><td>952</td></tr><tr><td rowspan="3">250</td><td>2</td><td>1.06%</td><td>1297</td><td>0.71%</td><td>7</td><td>16</td><td>5301</td><td>6177</td></tr><tr><td>3</td><td>0.88%</td><td>1424</td><td>0.67%</td><td>3</td><td>16</td><td>6293</td><td>7147</td></tr><tr><td>3</td><td>1.70%</td><td>1328</td><td>1.57%</td><td>0</td><td>21</td><td>6474</td><td>7639</td></tr><tr><td rowspan="3">250</td><td>2</td><td>0.97%</td><td>2087</td><td>0.91%</td><td>0</td><td>14</td><td>6289</td><td>7405</td></tr><tr><td>3</td><td>0.97%</td><td>2087</td><td>0.91%</td><td>0</td><td>14</td><td>6289</td><td>7405</td></tr></tbody></table>
<table><caption>Table 10: Detailed computational results for BCP configuration G'</caption><thead><tr><th>Set</th><th>n</th><th>m</th><th>Root Gap</th><th>Root Time</th><th>Gap</th><th>#Opt</th><th>TN</th><th>SP Time</th><th>Time</th></tr></thead><tbody><tr><td rowspan="12">A</td><td rowspan="3">50</td><td>3</td><td>0.15%</td><td>4</td><td>0.03%</td><td>29</td><td>119</td><td>117</td><td>264</td></tr><tr><td>4</td><td>0.09%</td><td>3</td><td>0.00%</td><td>30</td><td>1</td><td>2</td><td>26</td></tr><tr><td>100</td><td>3</td><td>0.08%</td><td>12</td><td>0.00%</td><td>30</td><td>1</td><td>9</td><td>78</td></tr><tr><td rowspan="3">150</td><td>3</td><td>0.08%</td><td>16</td><td>0.00%</td><td>30</td><td>1</td><td>13</td><td>86</td></tr><tr><td>4</td><td>0.08%</td><td>31</td><td>0.00%</td><td>30</td><td>1</td><td>24</td><td>144</td></tr><tr><td>200</td><td>3</td><td>0.08%</td><td>47</td><td>0.00%</td><td>30</td><td>1</td><td>39</td><td>166</td></tr><tr><td rowspan="3">250</td><td>4</td><td>0.11%</td><td>94</td><td>0.00%</td><td>30</td><td>2</td><td>90</td><td>280</td></tr><tr><td>3</td><td>0.08%</td><td>89</td><td>0.00%</td><td>30</td><td>1</td><td>68</td><td>295</td></tr><tr><td>4</td><td>0.08%</td><td>211</td><td>0.00%</td><td>30</td><td>1</td><td>180</td><td>422</td></tr><tr><td rowspan="3">300</td><td>3</td><td>0.08%</td><td>157</td><td>0.00%</td><td>30</td><td>1</td><td>127</td><td>410</td></tr><tr><td>4</td><td>0.08%</td><td>410</td><td>0.00%</td><td>30</td><td>1</td><td>360</td><td>665</td></tr><tr><td>400</td><td>3</td><td>0.08%</td><td>275</td><td>0.00%</td><td>30</td><td>1</td><td>220</td><td>607</td></tr><tr><td rowspan="3">500</td><td>3</td><td>0.08%</td><td>843</td><td>0.00%</td><td>30</td><td>1</td><td>746</td><td>1185</td></tr><tr><td>4</td><td>0.09%</td><td>485</td><td>0.00%</td><td>30</td><td>1</td><td>397</td><td>900</td></tr><tr><td>4</td><td>0.09%</td><td>1896</td><td>0.01%</td><td>29</td><td>6</td><td>1840</td><td>2471</td></tr><tr><td rowspan="12">P</td><td rowspan="3">50</td><td>2</td><td>0.86%</td><td>0</td><td>0.00%</td><td>30</td><td>1</td><td>0</td><td>16</td></tr><tr><td>3</td><td>0.72%</td><td>0</td><td>0.00%</td><td>30</td><td>1</td><td>0</td><td>17</td></tr><tr><td>100</td><td>2</td><td>0.39%</td><td>4</td><td>0.00%</td><td>30</td><td>1</td><td>3</td><td>63</td></tr><tr><td rowspan="3">150</td><td>2</td><td>0.42%</td><td>8</td><td>0.00%</td><td>30</td><td>1</td><td>7</td><td>69</td></tr><tr><td>3</td><td>1.09%</td><td>42</td><td>0.00%</td><td>30</td><td>13</td><td>75</td><td>263</td></tr><tr><td>200</td><td>2</td><td>0.72%</td><td>96</td><td>0.00%</td><td>30</td><td>35</td><td>374</td><td>760</td></tr><tr><td rowspan="3">250</td><td>2</td><td>1.79%</td><td>167</td><td>0.70%</td><td>14</td><td>244</td><td>3252</td><td>5876</td></tr><tr><td>3</td><td>0.98%</td><td>494</td><td>0.71%</td><td>3</td><td>80</td><td>5231</td><td>7079</td></tr><tr><td>2</td><td>2.15%</td><td>456</td><td>1.97%</td><td>0</td><td>86</td><td>5338</td><td>7610</td></tr><tr><td rowspan="3">3</td><td>2</td><td>1.02%</td><td>1169</td><td>0.94%</td><td>0</td><td>27</td><td>6298</td><td>7562</td></tr></tbody></table>
<table><caption>Table 11: Detailed computational results for BCP configuration $G' + R1Cs$</caption><thead><tr><th>Set</th><th>n</th><th>m</th><th>Root Gap</th><th>Root Time</th><th>Gap</th><th>#Opt</th><th>TN</th><th>SP Time</th><th>Time</th></tr></thead><tbody><tr><td rowspan="12">A</td><td rowspan="3">50</td><td>3</td><td>0.15%</td><td>5</td><td>0.00%</td><td>30</td><td>1</td><td>3</td><td>27</td></tr><tr><td>4</td><td>0.09%</td><td>3</td><td>0.00%</td><td>30</td><td>1</td><td>2</td><td>26</td></tr><tr><td>100</td><td>3</td><td>0.08%</td><td>12</td><td>0.00%</td><td>30</td><td>1</td><td>9</td><td>78</td></tr><tr><td rowspan="3">150</td><td>3</td><td>0.08%</td><td>16</td><td>0.00%</td><td>30</td><td>1</td><td>13</td><td>86</td></tr><tr><td>4</td><td>0.08%</td><td>31</td><td>0.00%</td><td>30</td><td>1</td><td>24</td><td>144</td></tr><tr><td>200</td><td>3</td><td>0.08%</td><td>47</td><td>0.00%</td><td>30</td><td>1</td><td>39</td><td>166</td></tr><tr><td rowspan="3">250</td><td>3</td><td>0.11%</td><td>107</td><td>0.00%</td><td>30</td><td>1</td><td>97</td><td>281</td></tr><tr><td>4</td><td>0.08%</td><td>89</td><td>0.00%</td><td>30</td><td>1</td><td>68</td><td>295</td></tr><tr><td>300</td><td>3</td><td>0.08%</td><td>210</td><td>0.00%</td><td>30</td><td>1</td><td>179</td><td>421</td></tr><tr><td rowspan="3">400</td><td>3</td><td>0.08%</td><td>157</td><td>0.00%</td><td>30</td><td>1</td><td>126</td><td>411</td></tr><tr><td>4</td><td>0.08%</td><td>409</td><td>0.00%</td><td>30</td><td>1</td><td>359</td><td>665</td></tr><tr><td>500</td><td>3</td><td>0.08%</td><td>270</td><td>0.00%</td><td>30</td><td>1</td><td>214</td><td>602</td></tr><tr><td rowspan="6">P</td><td rowspan="3">50</td><td>2</td><td>0.86%</td><td>0</td><td>0.00%</td><td>30</td><td>1</td><td>0</td><td>16</td></tr><tr><td>3</td><td>0.72%</td><td>0</td><td>0.00%</td><td>30</td><td>1</td><td>0</td><td>17</td></tr><tr><td>100</td><td>2</td><td>0.39%</td><td>4</td><td>0.00%</td><td>30</td><td>1</td><td>3</td><td>63</td></tr><tr><td rowspan="3">150</td><td>2</td><td>0.35%</td><td>8</td><td>0.00%</td><td>30</td><td>1</td><td>6</td><td>69</td></tr><tr><td>3</td><td>0.37%</td><td>52</td><td>0.00%</td><td>30</td><td>1</td><td>33</td><td>184</td></tr><tr><td>200</td><td>2</td><td>0.42%</td><td>185</td><td>0.00%</td><td>30</td><td>2</td><td>178</td><td>364</td></tr><tr><td rowspan="3">250</td><td>2</td><td>0.98%</td><td>828</td><td>0.01%</td><td>29</td><td>7</td><td>1304</td><td>1838</td></tr><tr><td>3</td><td>0.84%</td><td>1724</td><td>0.32%</td><td>15</td><td>9</td><td>4928</td><td>5621</td></tr><tr><td>3</td><td>1.53%</td><td>2522</td><td>1.05%</td><td>5</td><td>9</td><td>6082</td><td>7145</td></tr><tr><td></td><td></td><td></td><td>0.97%</td><td>4611</td><td>0.88%</td><td>1</td><td>5</td><td>6591</td><td>7570</td></tr></tbody></table>
<table><caption>Table 12: Detailed computational results for BCP configuration G' + R1Cs + AR</caption><thead><tr><th>Set</th><th>n</th><th>m</th><th>Root Gap</th><th>Root Time</th><th>Gap</th><th>#Opt</th><th>TN</th><th>SP Time</th><th>Time</th></tr></thead><tbody><tr><td rowspan="12">A</td><td rowspan="3">50</td><td>3</td><td>0.13%</td><td>3</td><td>0.00%</td><td>30</td><td>1</td><td>2</td><td>24</td></tr><tr><td>4</td><td>0.09%</td><td>2</td><td>0.00%</td><td>30</td><td>1</td><td>1</td><td>25</td></tr><tr><td>100</td><td>3</td><td>0.08%</td><td>7</td><td>0.00%</td><td>30</td><td>1</td><td>5</td><td>73</td></tr><tr><td rowspan="3">150</td><td>3</td><td>0.08%</td><td>13</td><td>0.00%</td><td>30</td><td>1</td><td>10</td><td>82</td></tr><tr><td>4</td><td>0.08%</td><td>16</td><td>0.00%</td><td>30</td><td>1</td><td>12</td><td>129</td></tr><tr><td>200</td><td>3</td><td>0.08%</td><td>29</td><td>0.00%</td><td>30</td><td>1</td><td>24</td><td>148</td></tr><tr><td rowspan="3">250</td><td>3</td><td>0.08%</td><td>31</td><td>0.00%</td><td>30</td><td>1</td><td>23</td><td>192</td></tr><tr><td>4</td><td>0.08%</td><td>57</td><td>0.00%</td><td>30</td><td>1</td><td>48</td><td>221</td></tr><tr><td>4</td><td>0.08%</td><td>42</td><td>0.00%</td><td>30</td><td>1</td><td>32</td><td>249</td></tr><tr><td rowspan="3">300</td><td>3</td><td>0.08%</td><td>126</td><td>0.00%</td><td>30</td><td>1</td><td>110</td><td>338</td></tr><tr><td>4</td><td>0.08%</td><td>73</td><td>0.00%</td><td>30</td><td>1</td><td>57</td><td>326</td></tr><tr><td>4</td><td>0.08%</td><td>231</td><td>0.00%</td><td>30</td><td>1</td><td>206</td><td>486</td></tr><tr><td rowspan="3">400</td><td>3</td><td>0.08%</td><td>206</td><td>0.00%</td><td>30</td><td>1</td><td>192</td><td>558</td></tr><tr><td>4</td><td>0.08%</td><td>502</td><td>0.00%</td><td>30</td><td>1</td><td>456</td><td>844</td></tr><tr><td>500</td><td>3</td><td>0.09%</td><td>214</td><td>0.00%</td><td>30</td><td>1</td><td>171</td><td>630</td></tr><tr><td rowspan="12">P</td><td rowspan="3">50</td><td>2</td><td>0.86%</td><td>0</td><td>0.00%</td><td>30</td><td>1</td><td>0</td><td>16</td></tr><tr><td>3</td><td>0.72%</td><td>0</td><td>0.00%</td><td>30</td><td>1</td><td>0</td><td>17</td></tr><tr><td>100</td><td>2</td><td>0.39%</td><td>2</td><td>0.00%</td><td>30</td><td>1</td><td>1</td><td>61</td></tr><tr><td rowspan="3">150</td><td>2</td><td>0.35%</td><td>5</td><td>0.00%</td><td>30</td><td>1</td><td>4</td><td>66</td></tr><tr><td>3</td><td>0.38%</td><td>39</td><td>0.00%</td><td>30</td><td>1</td><td>19</td><td>171</td></tr><tr><td>200</td><td>2</td><td>0.41%</td><td>151</td><td>0.00%</td><td>30</td><td>3</td><td>197</td><td>388</td></tr><tr><td rowspan="3">250</td><td>2</td><td>0.98%</td><td>687</td><td>0.00%</td><td>30</td><td>6</td><td>1067</td><td>1588</td></tr><tr><td>3</td><td>0.84%</td><td>1554</td><td>0.39%</td><td>14</td><td>8</td><td>4546</td><td>5223</td></tr><tr><td>2</td><td>1.53%</td><td>2181</td><td>0.92%</td><td>10</td><td>9</td><td>5752</td><td>6794</td></tr><tr><td></td><td>3</td><td>0.97%</td><td>4450</td><td>0.89%</td><td>0</td><td>5</td><td>6723</td><td>7712</td></tr></tbody></table>