# The Robust Steiner Team Orienteering Problem with Decreasing Priorities under Budgeted Uncertainty

Lucas Assunção, Andréa Cynthia Santos

► To cite this version:

Lucas Assunção, Andréa Cynthia Santos. The Robust Steiner Team Orienteering Problem with Decreasing Priorities under Budgeted Uncertainty. Operations Research Perspectives, In press, 15, pp.100344. 10.1016/j.orp.2025.100344. hal-05129981

HAL Id: hal-05129981

https://hal.science/hal-05129981v1

Submitted on 15 Jul 2025

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The Robust Steiner Team Orienteering Problem with Decreasing Priorities under budgeted uncertainty

Lucas Assunção ID*, Andréa Cynthia Santos

LITIS, ISEL, Université Le Havre Normandie, 25 rue Philippe Lebon, Le Havre, 76600, France

## ARTICLE INFO

**Keywords:**

Humanitarian logistics
OR in disaster relief
Robust optimization
Cut-and-branch
Metaheuristics

## ABSTRACT

Post-disaster relief operations have gained attention over the past decade, focusing on enhancing resilience in labor and social environments. This work introduces the Robust Steiner Team Orienteering Problem with Decreasing Priorities (R-STOP-DP) to model emergency rescue operations where several locations might need relief shuttles, but exact demands cannot be foreseen. R-STOP-DP is a variation of the vehicle routing problem with location priorities that applies robust optimization to model the variability on service times incurred by visiting locations. Locations are sub-divided into mandatory and optional, being the latter linked to priorities that linearly decrease over time. The goal is to find robust feasible routes maximizing the priorities collected, while considering the worst-case conditions of service times within an *uncertainty budget* and a routes' duration limit. We propose two compact formulations – reinforced by valid inequalities adapted from the literature – and solve them in a cut-and-branch fashion. In addition, we propose a *kernel search mat-heuristic* and a *simulated annealing* heuristic. Computational experiments suggest the strict dominance of one formulation, improving dual bounds by 12.29% on average over the 360 instances tested. The cut-and-branch algorithm based on the stronger model also stands out, solving 20 more instances than the other. The simulated annealing heuristic obtains a remarkable performance by improving over and/or reaching the best-known bounds for the complete benchmark, within an average execution time of 2.52 s. In turn, the kernel search mat-heuristic reaches or improves the bounds for 81% of the instances within 4.5 min of average running time.

## 1. Introduction

Post disaster relief operations and risk management studies have drawn significant attention over the past decade, aiming at increasing resilience in labor and social environments through prevention and systematic response to emergency events [1–3]. Within this context, Operational Research (OR) plays a key role in devising mathematical methods to assess data, guide decisions and automate processes [4–6]. In disaster relief and risk management (as in other “real world” operations), demands, availability of resources and costs may vary according to the circumstances, and the quality of the decisions made is directly impacted by changes in the current scenario. Thus, considering data uncertainty is crucial to devising more accurate and robust solutions.

Different approaches can be applied to model problems under uncertain data, such as stochastic programming and Robust Optimization (RO). In the former, uncertainty is modeled through a probabilistic description and is mostly applied whenever the probability law associated with the uncertain data is known in advance. On the other hand, RO has gained particular attention over the last decades due to its power in modeling uncertainty while still applying the standard tooling of

classical/deterministic optimization. This characteristic is particularly advantageous for modeling disaster relief situations, where devising solutions with a short-time response is crucial, being even more important than the optimality of the plans in some cases. Not surprisingly, a recent trend is to apply RO to model these crisis situations, such as deciding where to position field hospitals after large-scale disasters (e.g., the earthquakes that hit Kathmandu, Nepal in April 2015) [7], as well as pre-positioning relief resources to mitigate supply shortage and surplus of relief items [8].

In this work, we introduce a new RO version of the Team Orienteering Problem with Decreasing Priorities (TOP-DP) [9], namely the Robust Steiner Team Orienteering Problem with Decreasing Priorities (R-STOP-DP). In the deterministic TOP-DP, a set of candidate locations to be visited is given, and each of them is linked to a priority value that linearly decreases over time. Thus, any feasible route must satisfy the deadline constraints at locations (i.e., all the locations selected must be served before their priorities reach zero), as well as a total route time limit. The goal is to maximize the total sum of priorities incurred by visiting locations, taking into account their decrease over time.

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: lucas.assuncao-de-almeida@univ-lehavre.fr (L. Assunção), andrea-cynthia.duhamel@univ-lehavre.fr (A.C. Santos).</footnote>
In R-STOP-DP, the variation we address in this work, we additionally consider locations with top priorities, which must necessarily be visited in the solution routes, and uncertain servicing times for visiting locations. Precisely, we treat the service times on locations as independent random variables which range according to nominal values and maximum (worst-case) deviations. Moreover, an *uncertainty budget* [10] parameter is applied to manage the level of uncertainty allowed by limiting the number of locations whose service times deviate to the worst-case scenario. This RO criterion has been chosen to model the data uncertainty in this work due to its potential of controlling the conservatism of the solutions (through tuning the budget value), yet avoiding the non-linearity/intractability that other criteria (e.g., the *ellipsoidal uncertainty*) can bring [10]. This same approach was adopted to model uncertainty in the R-STOP-DP version with no mandatory locations, namely Robust TOP-DP (R-TOP-DP) [11], and other works in disaster relief management [12,13].

R-STOP-DP is particularly suitable to model rescue and recovery operations in emergencies caused by natural disasters and incidents in industrial sites. In this sense, consider that, within a region affected by a disaster, several locations are potentially in need of a rescue shuttle or relief resources (such as food, water, medicine or other supply). Accordingly, R-STOP-DP can be used to devise rescue/supply routes – using land vehicles – for visiting/supplying first-aid tents and other selected locations where there might be victims, taking into account limitations in available resources (e.g., restricted autonomy of the vehicles, amount of relief resources, etc.). In these situations, priority must be given to locations that are more likely to have survivors close-by, also considering that the chances of finding victims still alive directly relies on a quick response, thus decreasing over time. Values of priority can be estimated based on factors as population density or historical data related to similar occurrences (in case of floods or mudslides, for instance) in the region or even by the existence of pre-established points of gathering in case of incidents in industrial sites. In addition, locations previously selected for first-aid tents, which are known to have victims, are set as mandatory visiting points in the solution routes. Moreover, it is expected that the servicing time at each possible location might be uncertain, due to unpredictable number of victims to be rescued or even uncertain waiting times caused by queuing, for example.

R-STOP-DP also generalizes the variation of the Team Orienteering Problem (TOP) [14] with mandatory vertices, namely the Steiner Team Orienteering Problem (STOP) [15,16]. STOP (and, by consequence, R-STOP-DP) is especially challenging, as solely finding a feasible solution (or proving its infeasibility) is already NP-hard [15]. This difficulty in generating a viable initial solution adds a new layer of complexity to the problem as opposed to its version without mandatory vertices – R-TOP-DP [11] –, which can be particularly problematic in disaster relief applications, where response speed plays a crucial role, being even more important than the optimality of the plans in some cases.

In this work, we propose two compact formulations for R-STOP-DP and further reinforce them with two families of valid inequalities adapted from the literature, namely arc-vertex inference cuts and clique conflict cuts. The first formulation incorporates the restriction of mandatory vertices to the flow/commodity model of [11] for R-TOP-DP, which has the drawback of defining the uncertainty-related space by means of “big-M” constraints. In the second formulation, we avoid these complicating constraints by describing the commodity variation in an alternative way: instead of increasing along the route, the commodity is consumed in each of the selected arcs. We experimentally show that the latter formulation tightens the dual limits provided by the former, and, in particular, gives strictly better bounds for all the instances tested. The reinforcement implied by the valid inequalities is also quantified for both formulations. We also solve the integral models by means of a cut-and-branch algorithm – incorporating the valid inequalities as root user cuts – to assess how challenging the original problem is in practice.

In addition, we propose a mat-heuristic based on the *kernel search* framework [17], as well as a *simulated annealing* [18] heuristic. The primal bounds obtained by the latter are used to warm-start the cut-and-branch algorithm, favoring its convergence.

The experiments indicate that the commodity-based formulation, coupled with the valid inequalities, can significantly improve the dual bounds (on average, 17.12%) when compared to the weaker model. Moreover, the stronger model, when solved to (integer) optimality, was able to provide 19 more optimality certificates. The cutting-plane algorithm achieved a limited solution improvement when compared to the plain models: it solves one more instance to optimality, while decreasing the average optimality gap by almost 4%. The results suggest that solving R-STOP-DP exactly can be already challenging for small to medium-sized instances of 25 and 50 locations.

With respect to the heuristics proposed, the simulated annealing clearly stands out both in terms of solution quality and execution times. Precisely, its multi-threading implementation finds, for all the instances tested, equal or better bounds than the best ones provided by the exact models within a tiny average execution time of 2.52 s. The kernel search heuristic, in turn, achieves moderate success, improving or matching the best-known bounds for 81% of the instances. With an average execution time of 4.5 min, the mat-heuristic remains a good solution option due to its simplicity and low implementation effort.

The remainder of this study is organized as follows. Related works are briefly discussed in Section 2, and the formal definition of R-STOP-DP is presented in Section 3. The mathematical formulations proposed are described in Section 4, and the cut-and-branch algorithm used to solve the formulations exactly is detailed in Section 5, along with the families of valid inequalities used to reinforce the original models. The kernel search mat-heuristic is described in Section 6, followed by the simulated annealing heuristic proposed (Section 7). Some implementation choices are pinpointed in Section 8, and Section 9 is devoted to the computational results. Concluding remarks are drawn in the last section.

## 2. Related works

R-STOP-DP is closely related to the well-studied deterministic orienteering problems [19], which are defined on a graph, usually complete and undirected, where a value of priority is associated with each vertex, and a traverse time is associated with each edge (or arc). They aim at finding one or more routes (in case of single or multi-vehicle variations) from an origin vertex to a destination one, while satisfying a total traverse time constraint and maximizing the cumulative priorities of the vertices visited. As in the original orienteering sport, a vertex/location can only be serviced once. Contrary to the classical traveling salesman problem [20], solutions for orienteering problems do not necessarily visit all the vertices of the graph, but prioritize the ones associated with higher priority values.

Several variations of orienteering problems have been addressed in the literature. The original problem, namely the Orienteering Problem (OP) [19], considers a single vehicle. When the origin and destination vertices coincide, the problem is known as the selective traveling salesman problem [21]. Moreover, when a team of vehicles working together is considered, the problem becomes the Team Orienteering Problem (TOP) [14]. Time windows constraints [22,23], priorities decreasing over time [9], precedence order in the visitation of vertices [24] and heterogeneous vehicles [25] are also addressed in the literature [22,23]. When some of the vertices to be visited are set as mandatory, the Steiner TOP (STOP) arises [15,16].

Currently, STOP finds its main applications in devising the itinerary of the delivery of goods by shipping companies and in the planning of technical visits, where a priority value – which may rely on factors such as the urgency of the request and the customer’s rank/category – is associated with visiting each customer/location, and deliveries/visits with top priority (e.g., those whose deadlines are expiring) are set as
mandatory. Quite recently, the potential of orienteering problems for modeling post disaster relief operations was explored for the first time in an application of rescuing victims from industrial sites in Italy during emergency situations [1].

An RO version of OP was also proposed to model the planning of information gathering by Unmanned Aerial Vehicles (UAV) in military missions [26]. In the work, uncertainty is considered both in the resource consumption of the routes and in the priorities of the locations to be visited. The authors propose a mathematical formulation and solve instances with up to 20 vertices directly through a commercial solver. Several other non-deterministic variations of OP have already been studied, but modeling the uncertainty through stochastic programming. In these problems, the uncertainty is considered either in the priority values [27–29] or in time-related parameters (travel time, service time, etc.) [30–32], or both [33]. We refer to [11,34] for an extensive survey of classical and non-deterministic variations of OP.

To the best of our knowledge, R-STOP-DP has not been addressed in the literature before, except for its specific case with no mandatory vertices, namely R-TOP-DP [11]. In the study, the authors propose a compact formulation and a route-based one (with an exponential number of variables), which is solved through a branch-and-price algorithm. They also present a tabu search heuristic to solve large-sized instances with up to 100 vertices. The addition of the constraint on mandatory vertex visitation makes R-STOP-DP more challenging than R-TOP-DP both in theory and in practice. While a trivial solution of empty routes is always feasible for R-TOP-DP, finding an initial solution for R-STOP is known to be NP-hard [15], and consists of solving a multiple traveling salesman problem [35] in, at least, the worst-case scenario for service times. Our computational results on solving the compact formulation extended from [11] show, in practice, the additional difficulty of R-STOP-DP compared to R-TOP-DP for instances with the same dimensions.

In addition to dealing with a more challenging version of the problem, we propose a stronger bound reformulation obtained by replacing the “big-M” constraints of the previously introduced model and further strengthen it through the separation of two families of valid inequalities. Another contribution of our work is the proposal of two heuristic approaches. The first one is based on the kernel search framework [17]. To our knowledge, this is only the second time the framework is applied to address an OP variation [24] and the first time for an RO problem. The second heuristic we propose follows a more classical path, using well-known local searches from the literature of vehicle routing within a simulated annealing metaheuristic. A key distinction from prior work lies in its implementation choices, being tailored to fully benefit from multi-threading and, as a result, it finds solutions attaining and improving over the best-known bounds at remarkably tiny execution times for the whole benchmark of instances tested.

### 3. Problem definition and notation

R-STOP-DP is defined on a digraph $G = (N, A)$, where $N$ is the vertex set, and $A$ is the arc set. Let $0 \in N$ be the origin, $S \subseteq N \setminus \{0\}$ be the subset of mandatory vertices, and $P \subseteq N \setminus \{0\}$ be the set of optional vertices, such that $S \cap P = \emptyset$ and $N = S \cup P \cup \{0\}$. A priority $p_i \in \mathbb{Z}^+$ is associated with each vertex $i \in P$, and a traverse time $t_{ij} \in \mathbb{R}^+$ is associated with each arc $(i, j) \in A$. A maximum route duration limit $T$ is given, and the arc traverse times are assumed to satisfy the triangle inequality. Moreover, we assume that each priority $p_i, i \in P$, decreases linearly over the time with a decay ratio $d_i$. Accordingly, a vertex can only be visited while its priority exists, i.e., before its deadline, given by $D_i = p_i/d_i$ for all $i \in P$. We consider a fixed $D_i = T$ for all $i \in S \cup \{0\}$.

Moreover, at each vertex $i \in N$, the nominal service time is given by $\bar{s}_i$, and $\hat{s}_i$ is its largest deviation. Naturally, $\bar{s}_0 = \hat{s}_0 = 0$. Each vehicle of the given homogeneous fleet $F$ can run for no more than the time limit $T$. Then, R-STOP-DP aims at finding up to $|F|$ routes that start and end at the same vertex 0, such that every mandatory vertex in $S$ belongs to

exactly one route and the total sum of priorities of the optional vertices visited is maximized. Here, each optional vertex in $P$ can be visited by at most one vehicle, thus avoiding the multiple accounting of a same priority. In R-STOP-DP, the feasible region is subject to uncertainty, and a solution for the problem must comply with an uncertainty budget applied to the uncertain service times at vertices, as detailed in the sequel.

#### 3.1. Uncertainty characterization

Here, the uncertain service times at vertices are treated as independent random variables that can assume values in the range

$$ \bar{s}_i = \bar{s}_i + \xi_i \hat{s}_i, \xi_i \in [-1, 1], \forall i \in N \setminus \{0\}, \quad (1) $$

where $\xi_i$ are parameters used to simulate the uncertainty. Since the worst-case scenario for service times only occur for positive deviations from the nominal values, we assume $\xi_i \in [0, 1]$ instead for simplicity. Then, as in [11], we apply the budgeted uncertainty [10] as RO criterion to control the randomness associated with $\bar{s}$. Accordingly, the budget uncertainty set is defined as

$$ \chi = \left\{ \bar{s} \in \mathbb{R}^{|N \setminus \{0\}|} : \bar{s}_i = \bar{s}_i + \xi_i \hat{s}_i, \xi_i \in [0, 1], \forall i \in N \setminus \{0\}, \sum_{i \in N \setminus \{0\}} \xi_i \le \Gamma \right\}, \quad (2) $$

where $\Gamma \in \mathbb{Z}$, $0 \le \Gamma \le |N \setminus \{0\}|$, is a parameter called uncertainty budget, which determines the level of conservatism of the solution by limiting the number of vertices whose service times reach the worst-case scenario $\bar{s} + \hat{s}$ in a route. Notice that, when $\Gamma = 0$, the uncertainty is completely eliminated, and the deterministic counterpart of the problem arises, with $\bar{s} = \hat{s}$. The uncertainty characterization described above implicitly defines the set of possible duration times of (not necessarily feasible) routes for R-STOP-DP.

**Theorem 1.** The (worst-case) scenarios of interest in (2) can be fully represented by $\xi$ in its integer bounds, i.e., $\xi \in \{0, 1\}^{|N \setminus \{0\}|}$.

**Proof.** Given a solution route for R-STOP-DP, its total time duration is given by the sum of (i) the (deterministic) time incurred by traversing its arcs, (ii) the nominal service times of the vertices visited, and (iii) the time values that deviate from the nominal service times. Since (i) and (ii) are fixed, we only have to compute the maximum value of (iii) to determine the worst-case scenario for the route. Notice that this task corresponds to the continuous knapsack problem: each vertex $i \in N \setminus \{0\}$ in the route is treated as an item of weight $\bar{s}_i$, the knapsack capacity is given by the uncertainty budget $\Gamma$, and the solution aimed corresponds to an assignment of $\xi \in [0, 1]^{|N \setminus \{0\}|}$ — the (fractional) selection of items. Since $\Gamma$ is an integer and all $\xi$ values linked to the vertices are at most one, there is always an optimal solution with $\xi$ integer [36]. $\square$

Considering Theorem 1, we model the uncertainty characterization in a compact manner by means of the continuous decision variables $a_{ij}$, which indicate, for all $i \in N$, $\gamma \in \{0, ..., \Gamma\}$, the maximum arrival time at a vertex $i$ when up to $\gamma$ predecessors vertices (in that same route) have reached their worst-case service times. These arrival times can be efficiently computed though a set of recursive equations based on dynamic programming [11,37–39] as follows. Given a route $r = (v_0 = 0, v_1, ..., v_m, v_{m+1} = 0)$ visiting $m$ vertices of $N \setminus \{0\}$,

$$ a_{v_i, \gamma} = \begin{cases} 0, & i = 0, 0 \le \gamma \le \Gamma, \\ t_{0v_i}, & i = 1, 0 \le \gamma \le \Gamma, \end{cases} \quad (3a) $$

$$ a_{v_{i-1}, \gamma} + \bar{s}_{v_{i-1}} + t_{v_{i-1}, v_i}, \quad i > 1, \gamma = 0, \quad (3b) $$

$$ \max \left\{ a_{v_{i-1}, \gamma} + \bar{s}_{v_{i-1}}, a_{v_{i-1}, \gamma-1} + \bar{s}_{v_{i-1}} + \hat{s}_{v_{i-1}} \right\} + t_{v_{i-1}, v_i}, \quad i > 1, 0 < \gamma \le \Gamma. \quad (3c) $$

$$ \max \left\{ a_{v_{i-1}, \gamma} + \bar{s}_{v_{i-1}} + t_{v_{i-1}, v_i}, a_{v_{i-1}, \gamma-1} + \bar{s}_{v_{i-1}} + \hat{s}_{v_{i-1}} \right\} + t_{v_{i-1}, v_i}, \quad i > 1, 0 < \gamma \le \Gamma. \quad (3d) $$
Eqs. (3a) and (3b) set the starting time in the origin to zero and make sure that the arrival time at the first vertex of the route corresponds to the time needed to traverse the arc from the origin to it. Eqs. (3c) define the arrival times at vertices $v_i$, $i > 1$, when no vertex in the route is allowed to reach the worst-case scenario (i.e., $\gamma = 0$). In this case, $a_{v_i\gamma} = a_{v_i0}$ is computed by simply adding to $a_{v_{i-1}0}$ – the arrival time at the predecessor vertex $v_{i-1}$ – the nominal service time of $v_{i-1}$, plus the traverse time of the arc $(v_{i-1}, v_i)$. When at least one vertex is allowed to deviate from the nominal values (i.e., $0 < \gamma \le \Gamma$), the arrival time at a vertex $v_i$, $i > 1$, is the maximum value among two situations: (i) $\gamma$ predecessors have already been set to their maximum service times before reaching $v_{i-1}$, or (ii) only up to $\gamma - 1$ vertices in the route have attained their worst-case scenarios before reaching $v_{i-1}$, and, thus, there is still budget to include $v_{i-1}$ at its maximum deviation time $\bar{s}_{v_{i-1}} + \hat{s}_{v_{i-1}}$. These two situations are described in (3d).

Finally, a route $r = (v_0 = 0, v_1, \dots, v_m, v_{m+1} = 0)$ is said to be (robust) feasible with respect to the budgeted uncertainty if the arrival times (when up to $\Gamma$ vertices are at their worst-case) at all of its optional vertices is at most their deadlines (i.e., $a_{v_i\Gamma} \le D_{v_i}$, for all $1 \le i \le m$) and, obviously, the total time limit of the route is satisfied (i.e., $a_{v_{m+1}\Gamma} \le T$, $v_{m+1} = v_0 = 0$).

# 4. Mathematical formulations

In this section, we describe the two compact formulations for R-STOP-DP addressed in this work: the first one was obtained by incorporating the mandatory vertices restriction to the flow-based formulation of [11] for R-TOP-DP, where flow is used to represent the increasing arrival times at vertices. The second formulation is also based on flow/commodity, but, now, each vehicle is expected to leave the origin with $T$ units of commodity to be consumed while traversing each arc and vertex.

## 4.1. Vertex-flow compact formulation ($\mathcal{F}_1$)

Consider the decision variables $y$ on the choice of vertices belonging or not to the solution routes, such that $y_i = 1$ if the vertex $i \in N \setminus \{0\}$ is visited by a vehicle of the fleet, and $y_i = 0$, otherwise. Likewise, let the binary variables $x$ identify the solution routes themselves: $x_{ij} = 1$ if the arc $(i, j) \in A$ is traversed in the solution; $x_{ij} = 0$, otherwise. Additionally, the continuous variables $a_{ij\gamma}$, for all $i \in N$, $\gamma \in \{0, \dots, \Gamma\}$, represent the time a vehicle arrives at a vertex $i$ when at most $\gamma$ predecessors' service times have reached their worst-case values, and the slack variable $\varphi$ keeps the number of vehicles that are not used in the solution. Given a vertex $i \in N$, we define the sets of arcs leaving and entering $i$ as $\delta^+(i) = \{j \in N : (i, j) \in A\}$ and $\delta^-(i) = \{j \in N : (j, i) \in A\}$. The first compact formulation is defined as follows.

$$ (\mathcal{F}_1) \max \sum_{j \in P} (p_j y_j - d_j a_{j\Gamma}), \quad (4) $$

$$ \text{s.t.} \quad y_i = 1 \quad \forall i \in S \cup \{0\}, \quad (5) $$

$$ \sum_{j \in \delta^+(i)} x_{ij} = y_i \quad \forall i \in S \cup P, \quad (6) $$

$$ \sum_{j \in \delta^+(0)} x_{0j} + \varphi = |F|, \quad (7) $$

$$ \sum_{j \in \delta^+(i)} x_{ij} - \sum_{j \in \delta^-(i)} x_{ji} = 0 \quad \forall i \in N, \quad (8) $$

$$ a_{i\gamma} \ge t_{0i} y_i \quad \forall i \in \delta^+(0), \forall \gamma \in \{0, \dots, \Gamma\}, \quad (9) $$

$$ a_{i\gamma} + \bar{s}_i + t_{ij} \le a_{j\gamma} + M'_{ij}(1 - x_{ij}) \quad \forall i \in S \cup P, j \in \delta^+(i), \forall \gamma \in \{0, \dots, \Gamma\}, \quad (10) $$

$$ a_{i\gamma-1} + \bar{s}_i + \hat{s}_i + t_{ij} \le a_{j\gamma} + M''_{ij}(1 - x_{ij}) \quad \forall i \in S \cup P, j \in \delta^+(i), \forall \gamma \in \{1, \dots, \Gamma\}, \quad (11) $$

$$ a_{i\gamma} \le D_i y_i \quad \forall i \in N, \forall \gamma \in \{0, \dots, \Gamma\}, \quad (12) $$

$$ x_{ij} \in \{0, 1\} \quad \forall (i, j) \in A, \quad (13) $$

$$ y_i \in \{0, 1\} \quad \forall i \in N, \quad (14) $$

$$ a_{i\gamma} \ge 0 \quad \forall i \in N, \forall \gamma \in \{0, \dots, \Gamma\}, \quad (15) $$

$$ 0 \le \varphi \le |F|. \quad (16) $$

The objective function in (4) gives the sum of priorities incurred by visiting optional vertices. Notice that the loss in the priority values (over time) linked to the decay ratio is taken into account. In particular, for each vertex $j \in P$, the corresponding decrease rate is multiplied by the vertex's maximum arrival time in its worst-case scenario ($a_{j\Gamma}$), i.e., when up to $\Gamma$ vertices reach their maximum service times in the route from the origin to $j$. Here, the other $a$ variables indexed in $\{0, \dots, \Gamma - 1\}$ are used to calculate $a_{j\Gamma}$ in a bottom-up way – through constraints (9)–(11) – as discussed in the sequel. Constraints (5) impose that all mandatory vertices are selected, while constraints (6) ensure that each vertex in $S \cup P$ is visited at most once. Restrictions (7) ensure that at most $|F|$ vehicles leave the origin 0, whereas constraints (8), along with constraints (5) and (6), guarantee that, if a vehicle visits a vertex $i \in N$, then it must enter and leave this vertex exactly once.

Constraints (9)–(11) represent the recursive Eqs. (3b)–(3d), which compute the possible arrival times at vertices when considering the uncertainty budget $\Gamma$. In particular, (9) are related to (3b). Here, the equality is relaxed to simplify the model without losing its correctness. We also point out that (3a) were not represented, as they are implicitly inferred by (9). Constraints (10) and (11) depict (3c) and (3d) by means of Miller–Tucker–Zemlin-like constraints [40] that are activated/deactivated through sufficiently large constants $M'_{ij} = D_i + \bar{s}_i + t_{ij}$ and $M''_{ij} = D_i + \hat{s}_i + \bar{s}_i + t_{ij}$ for all $(i, j) \in A$, $i \neq 0$. Here, the management of the flow associated with the variables $a$ also avoids the existence of sub-tours in the solutions.

Constraints (12) force the solution routes to meet the vertex deadlines and the total duration of the route $T$. Recall that a deadline $T$ is assumed for the mandatory vertices and the origin. These constraints also link the flow variables $a$ to $y$ by making sure that, if the arrival time at a given vertex is other than zero (regardless of the amount of uncertainty budget used), then the vertex is selected in the solution. The remaining restrictions define the domain of the decision variables.

## 4.2. Decreasing arc-flow compact formulation ($\mathcal{F}_2$)

In this formulation, inspired by the one of [15] for the deterministic STOP, time is treated as a commodity to be spent by the vehicles when traversing each arc in their routes, such that every vehicle departs from 0 with an initial amount of $T$ units of commodity, the time limit. Consider the same $x$, $y$ and $\varphi$ decision variables of $\mathcal{F}_1$ and let the flow variables $f_{ij\gamma}$, for all $(i, j) \in A$, $\gamma \in \{0, \dots, \Gamma\}$, represent the amount of time still available for a vehicle after traversing the arc $(i, j)$ when up to $\gamma$ predecessors' service times (including vertex $i$) have reached their worst-case values. Note that the maximum possible values for $f_{ij}$, $(i, j) \in A$, correspond to the minimum possible arrival time $a_j$ of $\mathcal{F}_1$. The second compact formulation is defined as follows.

$$ (\mathcal{F}_2) \max \sum_{j \in P} (p_j y_j - d_j (Ty_j - \sum_{i \in \delta^-(j)} f_{ij} \Gamma)), \quad (17) $$

s.t. Constraints (5)–(8),

$$ f_{ij\gamma} \le (T - \Delta_{0i} - \bar{s}_i - t_{ij}) x_{ij} \quad \forall (i, j) \in A, \forall \gamma \in \{0, \dots, \Gamma\}, \quad (18) $$

$$ f_{ij\gamma} \le \sum_{k \in \delta^-(i)} f_{k\gamma} - (\bar{s}_i + t_{ij}) x_{ij} \quad \forall (i, j) \in A, i \neq 0, \forall \gamma \in \{0, \dots, \Gamma\}, \quad (19) $$

$$ f_{ij\gamma} \le \sum_{k \in \delta^-(i)} f_{k\gamma-1} - (\bar{s}_i + \hat{s}_i + t_{ij}) x_{ij} \quad \forall (i, j) \in A, i \neq 0, \forall \gamma \in \{1, \dots, \Gamma\}, \quad (20) $$
$$f_{ij\gamma} \geq \max\{T - D_j, \bar{s}_j + \Delta_{j0}\}x_{ij} \quad \forall j \in N \setminus \{0\}, i \in \delta^{-}(j), \forall \gamma \in \{0, \dots, \Gamma\}, \quad (21)$$

$$x_{ij} \in \{0, 1\} \quad \forall (i, j) \in A, \quad (22)$$

$$y_i \in \{0, 1\} \quad \forall i \in N, \quad (23)$$

$$f_{ij\gamma} \geq 0 \quad \forall (i, j) \in A, \forall \gamma \in \{0, \dots, \Gamma\}, \quad (24)$$

$$0 \le \varphi \le |F|, \quad (25)$$

where $\Delta_{ij}$, $i, j \in N$, is the minimum time of a route from $i$ to $j$, including the vertices nominal service times, except those of $i$ and $j$. Naturally, $\Delta_{ii} = 0$ for all $i \in N$.

The objective function represents the amount of priorities linked to the vertices visited, still considering their decrease over time. Constraints (18)–(20) depict Eqs. (3b)–(3d) in an alternative way, without the need of the “big-M” complicating constants of $F_1$. While (19) and (20) simulate the time consumption on arcs – based on (3c) and (3d) –, (18) ensure (3b). Precisely, when $i = 0$, (18) state that the time available after traversing any arc $(0, j) \in A$ is at most the initial amount $T$, minus $t_{0j}$ (recall that $\Delta_{00} = \bar{s}_0 = 0$). The remaining (18) constraints give upper bounds on the $f$ values, while (21) provide lower bounds. Restrictions (22)–(25) define the variables’ domain.

Although both $F_1$ and $F_2$ are compact, their size can grow significantly due to restrictions (9)–(12), (18)–(21) and variables $a$ and $f$. In particular, they define the feasible space subject to uncertainty, which is directly proportional to the level of uncertainty controlled by the $\Gamma$ parameter.

## 5. Cut-and-branch algorithm

In this section, we describe the cut-and-branch algorithm used to solve formulations $F_1$ and $F_2$ exactly. To this end, we first discuss two families of valid inequalities originally devised for STOP [41], but extensible to R-STOP-DP. As detailed by the end of this section, these inequalities are separated within the cut-and-branch algorithm.

### 5.1. Arc-vertex inference cuts (AVICs)

Consider the set $E = \{\{i, j\} : (i, j) \in A \text{ and } (j, i) \in A, i, j \neq 0\}$. AVICs are defined as

$$|y_i - y_j| \le 1 - (x_{ij} + x_{ji}) \quad \forall \{i, j\} \in E, \quad (26)$$

which can be linearized as

$$y_i - y_j \le 1 - (x_{ij} + x_{ji}) \quad \forall \{i, j\} \in E, \text{ and} \quad (27)$$

$$y_j - y_i \le 1 - (x_{ij} + x_{ji}) \quad \forall \{i, j\} \in E. \quad (28)$$

The validity of AVICs relies on two trivial properties that are also inherent to feasible solutions for $F_1$ and $F_2$: (i) for all $\{i, j\} \in E$, arcs $(i, j)$ and $(j, i)$ cannot be simultaneously selected (i.e., $x_{ij} + x_{ji} \le 1$), and (ii) once an arc $(i, j) \in A$ is selected in a solution (i.e., $x_{ij} = 1$), we must have $y_i = y_j = 1$. Notice that we excluded the origin 0 from the definition of $E$, as condition (i) does not hold for trivial routes $r = (v_0, v_1, v_2)$, with $v_0 = v_2 = 0$ and $v_1 \in S \cup P$. In this work, AVICs (27) and (28) are separated by complete enumeration with $O(|A|)$ time complexity, as the total number of inequalities is equal to $2 \times |E|$, which is at most $|A|$.

### 5.2. Clique conflict cuts (CCCs)

To define CCCs, consider the undirected conflict graph $G_c = (N \setminus \{0\}, K)$, where $K$ is the set of vertex pairs which cannot be simultaneously in a same valid route. Precisely, for every pair $\{i, j\} \in K$, with $i, j \in N \setminus \{0\}$, any route that visits both $i$ and $j$ (in any order) exceeds the limit $T$ or the deadline arriving time constraints related to

D. Also, let $Y$ be the set of all the cliques of $G_c$, which are referred to as conflict cliques. CCCs are defined as

$$\sum_{e \in \delta^{-}(V)} x_e \ge \sum_{i \in \sigma} y_i \quad \forall \sigma \in Y, \forall V \subseteq N \setminus \{0\}, \sigma \subseteq V, \quad (29)$$

$$\sum_{e \in \delta^{+}(V)} x_e \ge \sum_{i \in \sigma} y_i \quad \forall \sigma \in Y, \forall V \subseteq N \setminus \{0\}, \sigma \subseteq V. \quad (30)$$

These inequalities work by giving a lower bound on the number of routes leaving/entering a conflict clique. This bound corresponds to the number of vertices of that clique that are visited in a feasible solution, thus preventing conflicting vertices from belonging to a same route. Although the number of conflict cliques in a graph might be huge, it has been shown [41] that CCCs related to non-maximal conflict cliques can be discarded, as they are necessarily dominated by maximal ones. We highlight that isolated vertices in the conflict graph also count as maximal conflict cliques and, in this case, the corresponding CCCs turn into classical general connectivity constraints [42]. For separating CCCs of interest, we adapt the heuristic procedure of [41], in which a sequence of maximum flow problems are solved on the residual graph induced by the solution of the linear relaxation of formulations $F_1$ and $F_2$. The details are given in the following.

#### 5.2.1. Separation of CCCs

We start by computing the set $K$ of conflicting vertices. Precisely, two vertices $i, j \in N \setminus \{0\}$, $i \neq j$, are set as conflicting if two conditions hold: (i) there exists no path from 0 to 0 on $G$ that traverses both $i$ and $j$ (in this order) and satisfies the total time limit $T$ and the deadline of vertex $j$, and (ii) there is no path from 0 to 0 that traverses $j$ and $i$ (in this order) and satisfies both $T$ and the deadline of $i$. For simplicity, we only consider a subset $\tilde{K} \subseteq K$ of conflicting vertex pairs, such that, for all $i, j \in N \setminus \{0\}, i \neq j$,

$$\{i, j\} \in \tilde{K} \text{ iff } \begin{cases} (i) \Delta_{0i} + \bar{s}_i + \Delta_{ij} + \bar{s}_j + \Delta_{j0} > T \text{ or } \Delta_{0i} + \bar{s}_i + \Delta_{ij} > D_j, \text{ and} \\ (ii) \Delta_{0j} + \bar{s}_j + \Delta_{ji} + \bar{s}_i + \Delta_{i0} > T \text{ or } \Delta_{0j} + \bar{s}_j + \Delta_{ji} > D_j. \end{cases}$$

Notice that $\tilde{K}$ is not necessarily equal to the actual set of conflicting vertices $K$, since the criteria used to compute the latter simply aggregate sub-paths in $G$, thus yielding possible non-elementary routes. Once the corresponding conflict graph $\tilde{G}_c = (N \setminus \{0\}, \tilde{K})$ is built, we compute a subset $\tilde{Y} \subseteq Y$ of conflict cliques by finding all the maximal cliques of $\tilde{G}_c$, including the ones of size one. To this end, we apply the depth-first search algorithm of [43].

Let $(\bar{x}, \bar{y})$ be a possibly fractional solution referring to the linear relaxation of either $F_1$ and $F_2$, and define $\tilde{G} = (N, \tilde{A})$ as the residual graph induced by this solution, such that each arc $(i, j) \in A$ belongs to $\tilde{A}$ if, and only if, $\bar{x}_{ij} > 0$. In this graph, each arc $(i, j) \in \tilde{A}$ is given a capacity $c[i, j] = \bar{x}_{ij}$. Hereafter, given the residual graph $\tilde{G}$, we denote by $max\text{-}flow_{i \to j}(\tilde{G})$ the problem of finding the maximum flow $F_{i \to j}$ from $i$ to $j$, where $i, j \in N$. The corresponding minimum cut is denoted by $\theta_{i \to j}$, with $i \in \theta_{i \to j}$. Then, we iteratively look for violated CCCs by solving maximum flow problems on auxiliary graphs, which are built based on the residual graph of the current fractional solution. Precisely, for each conflict clique $\sigma' \in \tilde{Y}$, two auxiliary graphs are built, one for each type of CCC. The first one, denoted by $\tilde{G}'_1$, consists of the residual graph $\tilde{G}$ with an extra artificial vertex $\beta_1$ and $|\sigma'|$ arcs $(i, \beta_1)$, one for each $i \in \sigma'$. The second one, denoted by $\tilde{G}'_2$, is built by reversing all the arcs of $\tilde{G}$ and, then, adding an artificial vertex $\beta_2$, as well as an arc $(i, \beta_2)$ for each $i \in \sigma'$.

For short, $\tilde{G}'_1 = (\tilde{N}'_1, \tilde{A}'_1)$, where $\tilde{N}'_1 = N \cup \{\beta_1\}$ and $\tilde{A}'_1 = \tilde{A} \cup \{(i, \beta_1) : i \in \sigma'\}$, and $\tilde{G}'_2 = (\tilde{N}'_2, \tilde{A}'_2)$, where $\tilde{N}'_2 = N \cup \{\beta_2\}$ and $\tilde{A}'_2 = \{(v, u) : (u, v) \in \tilde{A}\} \cup \{(i, \beta_2) : i \in \sigma'\}$. The original capacities of the arcs of $\tilde{G}$ are preserved in $\tilde{G}'_1$ and $\tilde{G}'_2$. Moreover, all the additional arcs have a same capacity value, which is equal to a sufficiently large number (say $|F|$, the number of vehicles). Fig. 1 illustrates the construction of the residual and the auxiliary graphs described above.

Recall that the capacity of any arc in the residual graph is at most one, due to the bounds of the binary variables $x$. Moreover, due to
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Example of a residual graph and its corresponding auxiliary graphs built when considering a maximal conflict clique {1, 3} and a fractional solution $\tilde{x}$, with $\tilde{x}_{01} = \tilde{x}_{12} = 0.4$, $\tilde{x}_{03} = \tilde{x}_{32} = 0.6$ and $\tilde{x}_{20} = 1$. Here, the values associated with the arcs are their corresponding capacities, and the infinity symbol stands for a sufficiently large value.</figcaption></figure>

**Input:** A fractional solution $(\tilde{x}, \tilde{y})$, its corresponding residual graph $\tilde{G} = (N, \tilde{A})$ and the subset $\tilde{Y} \subseteq Y$ of maximal conflict cliques.

**Output:** A set of CCCs violated by $(\tilde{x}, \tilde{y})$.

1. for all $(\sigma' \in \tilde{Y})$ do

2. Step I. *Building the auxiliary graphs*

3. Build $\tilde{G}'_1$ and $\tilde{G}'_2$;

4. Step II. *Looking for a violated CCC (29)*

5. $(F_{0 \to \beta_1}, \theta_{0 \to \beta_1}) \leftarrow \text{max-flow}_{0 \to \beta_1}(\tilde{G}'_1)$;

6. if $(F_{0 \to \beta_1} < \sum_{i \in \sigma'} \bar{y}_i)$ then

7. Found violated CCC (29) with $V = N\theta_{0 \to \beta_1}$ and $\sigma = \sigma'$;

8. else

9. Step III. *Looking for a violated CCC (30)*

10. $(F_{0 \to \beta_2}, \theta_{0 \to \beta_2}) \leftarrow \text{max-flow}_{0 \to \beta_2}(\tilde{G}'_2)$;

11. if $(F_{0 \to \beta_2} < \sum_{i \in \sigma'} \bar{y}_i)$ then

12. Found violated CCC (30) with $V = N\theta_{0 \to \beta_2}$ and $\sigma = \sigma'$;

13. end-if

14. end-if-else

15. end-for;

<caption>**Fig. 2.** Algorithm used to separate violated CCCs.</caption>

constraints (6), the sum of flow entering/leaving every vertex other than the origin is at most one. Therefore, since the capacities of the arcs connecting the cliques to the artificial vertices is sufficiently large, in any minimum cut from $0$ to $\beta_1$ (or $\beta_2$), the clique is disconnected from $0$ and connected to the artificial vertex.

Given the current conflict clique $\sigma' \in \tilde{Y}$ considered and its corresponding auxiliary graphs, the algorithm first looks for a violated CCC (29) by computing the maximum flow from $0$ to $\beta_1$ on $\tilde{G}'_1$. Accordingly, let $(F_{0 \to \beta_1}, \theta_{0 \to \beta_1})$ be the solution of $\max - \text{flow}_{0 \to \beta_1}(\tilde{G}'_1)$. If $F_{0 \to \beta_1}$ is smaller than $\sum_{i \in \sigma'} \bar{y}_i$, a violated CCC (29) is identified, with $V = N\theta_{0 \to \beta_1}$. If the separation algorithm does not find a violated CCC (29) for the current clique, it then looks for a violated CCC of type (30) by computing $(F_{0 \to \beta_2}, \theta_{0 \to \beta_2}) = \max - \text{flow}_{0 \to \beta_2}(\tilde{G}'_2)$. Similarly, if $F_{0 \to \beta_2}$ is smaller than $\sum_{i \in \sigma'} \bar{y}_i$, then a violated CCC (30) is found, with $V = N\theta_{0 \to \beta_2}$. The whole separation procedure described above is summarized in Fig. 2. Since the number of maximal conflict cliques and CCCs can be large, we adopt some filtering criteria during the search of violated CCCs. First, every time a violated CCC is found for a given clique $\sigma$, we discard all the yet not explored cliques that contain at least one vertex of $\sigma$. Likewise, we only keep the most violated CCC found (if any) and the ones that are sufficiently orthogonal to it, as a way of balancing the strength and diversity of the cuts separated, while limiting the model size. More details are given in Section 8.

### 5.3. Main cut-and-branch algorithm

The algorithm here adopted to solve either $F_1$ or $F_2$ exactly follows a cutting-plane scheme. First, the linear relaxation of the original problem is reinforced by AVICs, which are separated by complete enumeration. Then, at each iteration of a cutting-plane phase, the algorithm looks for CCCs that are violated by the solution of the current LP model, as described in Section 5.2.1. This iterative procedure stops when either no more violated cuts are found or the bound improvement obtained by the current model (when compared to the one from the previous iteration) is no greater than a tolerance $\epsilon$. Once the cutting-plane phase is over, all the AVICs and CCCs separated are turned into user cuts, the integrality of the variables $x$ and $y$ is restored, and the resulting reinforced model is solved to optimality.

### 6. Kernel search (KS) mat-heuristic

Kernel search consists of a simple, yet effective, model-based heuristic framework (mat-heuristic) originally devised for the multi-dimensional knapsack problem [17] and later generalized to solve any mixed integer linear programming problem [44]. It relies on the idea that the decision variables can be ordered according to their more or less likelihood of appearing in an optimal solution. Precisely, given a priority criterion applied to the decision variables (like their actual values in the linear relaxation solution or their corresponding reduced costs), the “best” variables are kept in a kernel, and the
remaining ones are distributed into buckets. Then, iteratively, the original integral problem is solved considering only the current kernel's variables and the ones of a given candidate bucket, one at a time. The kernel incorporates the bucket's variables which compose a new incumbent solution every time a better quality solution is found.

In this work, we solve formulation $\mathcal{F}_2$ in a *Kernel Search* fashion by considering the $y$ variables (related to the selection of vertices) as the candidate variables to compose the kernel, as detailed in Fig. 3. The same strategy could be applied to solve formulation $\mathcal{F}_1$, but, here, we chose to test it only for the formulation with strongest bounds, as later shown in Section 9.2.

Consider the following input parameters of the heuristic: maximum size $\kappa$ of the buckets, the range $[\delta_{min}, \delta_{max}]$ to which the execution time of each sub-problem is subject, and the rate $\mu$ at which the execution time limit of each sub-problem decreases linearly. The algorithm starts by building the initial kernel $K$ and the buckets $\lambda \in \Lambda$ of candidate variables (lines 1–4). To this end, the linear relaxation of the formulation is solved, and the $y$ variables are ordered according to the following criterion. For all $i \in N$, let $\bar{y}_i$ and $rc(y_i)$ be the value of $y_i$ in the linear relaxation solution and its corresponding reduced cost, respectively. The origin and the mandatory vertices are given higher priority in the ordering, followed by the vertices whose $\bar{y}$ values are other than zero (in non-increasing order) and, at last, by the remaining vertices, in non-increasing order of reduced costs. The reduced costs (in non-increasing order) are also used to decide the cases of tied relaxation values. Then, considering this ordering, the initial kernel $K$ is filled with the first $\max(|S| + 1, \kappa)$ vertices, while the buckets are composed with the remaining ones, with up to $\kappa$ vertices per bucket. Note that the size of the initial kernel is defined as to guarantee the availability of the origin and the mandatory vertices in any sub-problem, thus avoiding infeasibility.

Then, the iterative local search procedure takes place (lines 5–13). Let $\mathcal{F}(N', v, \delta)$ denote the sub-problem of solving either $\mathcal{F}_1$ or $\mathcal{F}_2$ when only the vertices in $N' \subseteq N$ are eligible, an initial feasible solution $v$ is given, and only up to $\delta$ seconds are available for searching for a solution. The current best (possibly feasible) solution $v^*$ is initialized by solving $\mathcal{F}(K, \emptyset, \delta_{max})$, i.e., the sub-problem restricted to the initial variables in the kernel (lines 5–7). Here, no initial solution is considered, and the time limit is setup to $\delta \leftarrow \delta_{max}$ seconds. The search proceeds by solving, for each bucket $\lambda \in \Lambda$ in the sequence, the sub-problem $\mathcal{F}(K \cup \lambda, v^*, \delta)$ obtained by temporarily adding the bucket's variables to the current kernel and setting the current best solution as a starting point. At each iteration, the execution time limit is linearly decreased according to the rate $\kappa$, down to $\delta_{min}$ (line 10), as to devote less computational effort in the last buckets, which are less likely to have vertices of interest. When a new best bound solution is found by the end of an iteration, the new vertices (if any) that compose this solution are incorporated into the kernel (line 12) and considered in the subsequent sub-problems. At the end, the best solution $v^*$ found is returned.

## 7. Simulated annealing (SA) metaheuristic

In this section, we propose an SA heuristic for R-STOP-DP. As in the original SA framework [18], our heuristic works by sequentially searching for better solutions in the neighborhood of a given initial solution. In contrast to population-based metaheuristics, a single solution is kept in each iteration and used as starting point of a next iteration of the algorithm. Although the final goal is to find a best possible quality solution, degenerative moves (i.e., the ones that compromise solution quality) are allowed in the first iterations as a way to favor diversity. Accordingly, a temperature parameter which gradually cools down over time is used as the abstraction to mark the transition from initial iterations exploring larger neighborhoods (diversification phase) to more restricted searches for local optimality (intensification phase). The heuristic we propose also incorporates destroying and repairing

procedures more commonly used in large neighborhood search heuristics [45], coupled with classical local searches from the literature of routing problems in general. In particular, the latter procedures, which are described in Sections 7.2 and 7.4, are also present in most of the successful heuristics proposed to solve the deterministic TOP and STOP (e.g., [41,46–51]).

We describe, in Fig. 4, the general algorithm of the SA heuristic proposed. The heuristic receives as input a temperature decrease rate $\tau \in (0, 1)$ and an initial feasible solution $Y$, which is found by directly solving formulation $\mathcal{F}_2$ for the sub-problem considering only the mandatory vertices. This solution is improved (line 1, Fig. 4) through the local search procedures later detailed in Section 7.2 and, then, used to compute the initial temperature $t \leftarrow t_0$ (line 2), as described in Section 7.1.

At this point, an iterative procedure is performed until the temperature reaches zero (lines 3–14). First, a copy of the current solution $Y$ is stored as $Y'$. Then, $Y'$ is partially destroyed by the removal of some vertices, as later described in Section 7.3. In the sequel, the algorithm successively tries to improve $Y'$ (lines 6–8). To this end, the local searches of Section 7.2 are performed, followed by attempts of inter-route vertex swap perturbations, which are detailed in Section 7.4. If the latter attempts succeed (i.e., there exists a vertex swap that improves the current solution), the algorithm resumes to another round of local searches. Otherwise, based on the current temperature and the current quality of $Y'$, $Y$ is either updated to $Y'$ or left untouched. Precisely, let $cost(Y)$ denote the priority sum of the vertices visited in a solution $Y$. We update $Y \leftarrow Y'$ if either $cost(Y') > cost(Y)$ or $p \le \exp(-\frac{(cost(Y)-cost(Y'))}{\tau})$, where $p$ is randomly chosen in the interval $[0, 1]$ (lines 9–12$^t$).

In practical terms, according to the behavior of the exponential function, a deteriorated $-$ in terms of priority sum $-$ solution $Y'$ has higher chances of being accepted for subsequent iterations whenever we have greater temperatures (thus, in initial iterations) and/or smaller solution degradations, given by $\max\{0, cost(Y)-cost(Y')\}$, since R-STOP-DP is a maximization problem. In turn, improved solutions will always be accepted as new starting points for later iterations, regardless of the current temperature.

By the end of each iteration of the main loop, the temperature is updated according to the decrease factor $\tau$ (line 13). When the temperature reaches zero, the algorithm terminates, and the solution with the greatest priority sum found is returned.

### 7.1. Computing the initial temperature

We compute the initial temperature for the SA heuristic through the method proposed by [52], based on sampling. The procedure, which is detailed in Fig. 5, receives as input the sampling size $\phi$, the target acceptance probability $\hat{\chi}_0$ and a feasible solution $Y$. Intuitively speaking, the algorithm estimates the initial temperature for which the solution obtained in a given iteration of the SA would have an acceptance probability of $\hat{\chi}_0$, i.e., the solution would be accepted as starting point of a next iteration in $100 \cdot \hat{\chi}_0\%$ of the cases. Recall that the acceptance criteria are summarized in lines 9–12, Fig. 4.

The procedure starts by computing the sampling set $\Phi$ (line 1, Fig. 5). Precisely, starting from the feasible solution $Y$, it iteratively runs lines 4–8 of the algorithm in Fig. 4, and checks if the resulting $Y'$ solution degrades $Y$ in terms of priority sum, i.e., $cost(Y') < cost(Y)$. If so, the transition $(Y \rightarrow Y')$ is added to $\Phi$. At every iteration, the current $Y'$ is used as starting point of the next iteration, until $\phi$ degrading transitions are stored in $\Phi$. After the sampling, an initial candidate temperature $T$ is computed (line 2), and the sampling set is explored to further adjust this temperature until its estimated acceptance probability $\hat{\chi}(T)$ is sufficiently close to the target probability $\hat{\chi}_0$ (lines 3–10). We refer to [52] for a detailed explanation of the procedure, as well as the theory and the implementation details that guarantee its convergence.
**Input:** maximum size $κ$ of the buckets, time range [$\delta_{min}, \delta_{max}$] and decreasing rate $μ$ of the sub-problem's execution time.

**Output:** a heuristic solution $v^*$.

1. // Step I. Build the initial kernel and buckets.

2. Solve the linear relaxation of $F_1$;

3. Order the y variables prioritizing origin and mandatory vertices, then non-increasing order of relaxation values and reduced costs;

4. Select $\max\{|S| + 1, κ\}$ vertices to compose the kernel $K$, and distribute the remaining vertices into buckets $λ \in Λ$ of size up to $κ$, following the ordering;

5. // Step II. Search in the initial kernel.

6. $δ \leftarrow δ_{max}$;

7. $v^* \leftarrow F(K, \emptyset, δ)$;

8. // Step III. Search in candidate buckets.

9. **for all** ($λ \in Λ$) **do**

10. $δ \leftarrow \max\{\delta_{min}, μδ\}$;

11. $v^* \leftarrow F(K \cup λ, v^*, δ)$;

12. $K \leftarrow K \cup v^*$;

13. **end-for**;

<caption>Fig. 3. Kernel Search mat-heuristic used to solve the proposed formulations.</caption>

**Input:** A temperature decrease rate $τ \in (0, 1)$ and an initial feasible solution $Y$

1. Improve $Y$ through local searches (Section 7.2);

2. Compute and set initial temperature $t \leftarrow t_0$ (Section 7.1) based on $Y$;

3. **while** ($t \approx 0$) **do**

4. $Y' \leftarrow$ copy of $Y$;

5. Partially destroy $Y'$ by removing vertices (Section 7.3);

6. **do**

7. Improve $Y'$ through local searches (Section 7.2);

8. **while** (perform inter-route vertex swap perturbations on $Y'$) (Section 7.4);

9. $p \leftarrow$ random number in $[0, 1]$;

10. **if** ($cost(Y') > cost(Y)$ or $p \le \exp(-\frac{cost(Y)-cost(Y')}{t})$) **then**

11. $Y \leftarrow Y'$;

12. **end-if**;

13. $t \leftarrow t\tau$;

14. **end-while**;

<caption>Fig. 4. Description of the general SA algorithm.</caption>

**Input:** Sampling size $φ$, a target acceptance probability $\hat{\chi}_0$ and a feasible solution $Y$.

**Output:** Temperature $T$ with estimated acceptance probability close to $\hat{\chi}_0$.

1. Generate $t = 1, \dots, φ$ random degrading transitions $(Y_t \to Y'_t)$, $cost(Y'_t) < cost(Y_t)$, and store them in $\Phi$;

2. Set the temperature $T \leftarrow \frac{\sum_{t \in \Phi} (cost(Y_t) - cost(Y'_t))}{\Phi \ln(\hat{\chi}_0)}$;

3. **do**

4. Compute $\hat{\chi}(T) \leftarrow \frac{\sum_{t \in \Phi} \exp(-\frac{cost(Y_t)}{T})}{\sum_{t \in \Phi} \exp(-\frac{cost(Y'_t)}{T})}$;

5. **if** ($\hat{\chi}(T) \approx \hat{\chi}_0$) **then**

6. **break**;

7. **else**

8. $T \leftarrow T \cdot \left( \frac{\ln(\hat{\chi}(T))}{\ln(\hat{\chi}_0)} \right)$;

9. **end-if-else**;

10. **while** ($\hat{\chi}(T) \approx \hat{\chi}_0$);

11. **return** $T$;

<caption>Fig. 5. Procedure for computing the initial temperature of the SA algorithm of Fig. 4.</caption>
**Input:** An initial feasible solution.

1. **do**

2. Do inter and intra-route improvements (Section 7.2.1);

3. Try to insert unvisited vertices (Section 7.2.2);

4. Do vertex replacements (Section 7.2.3);

5. Try to insert unvisited vertices (Section 7.2.2);

6. **while** (did any improvement);

<figure><img src="image_3.png" /><figcaption>**Fig. 7.** Outline of the sequence of the inter and intra-route local searches.</figcaption></figure>

## 7.2. Local searches

The local search methods applied aim to enhance the quality of a given feasible solution, either by (i) increasing the total sum of priorities or (ii) reducing the overall route times while maintaining the sum of priorities. In fact, in all the procedures described next, a move will only be allowed if it results in one of these two improvements. The main local search algorithm systematically applies inter and intra-route improvements (Section 7.2.1), vertex insertion (Section 7.2.2) attempts and vertex replacements (Section 7.2.3), as summarized in Fig. 6. The algorithm terminates when a locally optimal solution is found, meaning that no further improvement can be made through the neighborhoods defined.

### 7.2.1. Inter and intra-route improvements

In each iteration of the algorithm outlined in Fig. 7, the current feasible solution is first subjected to inter-route improvements, which consist of exchanging vertices between different routes (lines 2–9, Fig. 7). Specifically, three types of exchanges are performed for all route combinations: (i) 1–1, where one vertex from a route is swapped with one from another route, (ii) 1–0, where a vertex is moved from one route to another, and (iii) 2–1, where two adjacent vertices from one route are exchanged with a single vertex from another. In all cases, an exchange is only executed if it maintains the solution's feasibility and satisfies one of the improvement goals (Section 7.2). Each vertex exchange procedure performs only the first possible valid exchange, examining the routes sequentially from beginning to end.

The inter-route improvement process continues until no further exchanges are possible. Afterwards, the algorithm moves on to intra-route improvements (lines 10–13, Fig. 7), where the classical 2-opt operator [53] is applied to optimize each route of the solution. If any improvement is made, the algorithm returns to the main loop, resuming the inter-route improvement process. If no further improvements are found, the algorithm terminates.

**Input:** An initial feasible solution.

1. **do**

2. Do 1-1 unvisited vertex exchange;

3. Do 2-1 unvisited vertex exchange;

4. **while** (did any improvement);

<caption>**Fig. 8.** Description of the sequence of vertex replacements.</caption>

### 7.2.2. Vertex insertion

Each time the insertion procedure is invoked, one of two priority orders for the unvisited vertices is randomly selected: either non-increasing or non-decreasing based on their original priority values. According to the chosen order, the unvisited vertices are tested individually for possible insertion into the current solution. If a vertex can be added without rendering the routes infeasible, it is inserted into the route and position that yield the smallest increase in the total time duration of the routes. If the vertex cannot be inserted, it remains unvisited, and the next vertex in the sequence is considered for insertion.

### 7.2.3. Vertex replacements

The vertex replacement procedure enhances a feasible solution by replacing visited vertices with unvisited ones (Fig. 8). The selection of unvisited vertices for potential insertion follows the same criteria described in the insertion procedure of Section 7.2.2. Here, two types of replacements are considered: (i) 1–1, where one visited vertex is replaced by an unvisited one, and (ii) 2–1, where two visited vertices are replaced by a single unvisited vertex. A replacement is only permitted if it maintains the feasibility and results in one of the solution improvements aimed (Section 7.2). For both types of replacements, the algorithm executes only the first one found, analyzing the routes sequentially from beginning to end.

## 7.3. Route (partial) destruction

The route destruction process involves eliminating a subset of optional vertices (while preserving all the mandatory ones) from a given feasible solution. To do this, a fixed parameter, denoted by removal_percentage ∈ [0,1], is introduced. Initially, we calculate an upper limit on the number of vertices that can be removed, referred to as max_number_of_removals. This value is randomly chosen from an open interval between zero and the product of the removal_percentage and the total number of optional vertices in the current solution. Then, the procedure makes up to max_number_of_removals attempts to remove optional vertices, according to one out of three possible orderings, chosen randomly: non-increasing order of priority values, non-decreasing order of priority values, or completely random order. Here, a vertex is only removed if it does not compromise the feasibility of the resulting route.

## 7.4. Inter-route vertex swap perturbation

The inter-route swap perturbation works by individually testing if each of the visited vertices can be moved to any other route. This is done from the first to the last route, from the beginning to the end of each route, in the order originally displayed in the solution. Additionally, each vertex can only be moved once. For every route, the algorithm attempts to relocate each of its vertices (whether mandatory or optional) to a different route within the solution. The goal is to select a destination route/position that results in the maximum possible increase in the solution's priority sum. A vertex may be placed in any position within the destination route, provided that both origin and destination routes maintain their feasibility.
<figure><img src="image_4.png" /><figcaption>Fig. 9. Average solution values obtained by the Simulated Annealing Heuristic (SAH) segmented by instance size, route duration limit (T) and uncertainty budget (Γ). Results segmented by instance size $|N\setminus\{0\}| \in \{25, 25, 100\}$, route duration limit (T) quartile in {0.8, 1, 2, 3, 4} and uncertainty budget $\Gamma \in \{0, 1, 5\}$.</figcaption></figure>

If at least one vertex is successfully relocated during the vertex movement attempt from a specific route, the algorithm continues by trying to insert unvisited vertices into that route. This insertion follows the same procedure outlined in Section 7.2.2, but restricted to the route currently under consideration.

## 8. Implementation details and parameter configurations

All the code was developed in C++, along with the optimization solver ILOG CPLEX 22.11¹ under default configurations. Regarding the separation of CCCs, we solved the maximum flow sub-problems with the implementation of the preflow push-relabel algorithm of [54] provided by the open-source Library for Efficient Modeling and Optimization in Networks — LEMON [55]. Additionally, we adopted the absolute violation criterion to determine which cuts are violated by a given solution during the separation of CCCs (see Section 5.2.1). In turn, to rank the CCCs found in terms of violation, we adopted the distance criterion. Given an n-dimensional column vector $w$ of binary variables, a point $\bar{w} \in \mathbb{R}^n$ and an inequality of the general form $a^T w \le b$, with $a \in \mathbb{R}^n$ and $b \in \mathbb{R}$, the absolute violation of this inequality with respect to $\bar{w}$ is simply given by $a^T \bar{w} - b$. Moreover, the distance

from $a^T w \le b$ to $\bar{w}$ corresponds to the Euclidean distance between the hyper-plane $a^T w = b$ and $\bar{w}$, which is equal to $\frac{(a^T \bar{w} - b)}{\|a\|}$, where $\|a\|$ is the Euclidean norm of $a$.

In our implementation of the cutting-plane algorithm, we set two extra parameters: a precision one (namely absolute violation precision) used to determine if a CCC is violated with respect to the absolute violation criterion, and another one to discard cuts that are not sufficiently orthogonal to the most violated ones. The latter parameter determines the minimum angle that an arbitrary cut must form with the most violated cut, as not to be discarded. In practice, this parameter establishes the maximum acceptable inner product between the arbitrary cut and the most violated one. Accordingly, we call it the maximum inner product. In the case of two inequalities $a_1^T w \le b_1$ and $a_2^T w \le b_2$, with $a_1, a_2 \in \mathbb{R}^n$ and $b_1, b_2 \in \mathbb{R}$, the inner product between them is given by $\frac{(a_1^T a_2)}{\|a_1\|\|a_2\|}$ and corresponds to the cosine of the angle defined by them. Here, the absolute violation precision and the maximum inner product were set to 0.01 and 0.03 (which corresponds to an approximate angle of 88.2°), respectively, based on pilot experiments. Moreover, the tolerance input value $\epsilon$ was set to $10^{-3}$.

Regarding the KS heuristic, we set the input parameters to: $[\delta_{min}, \delta_{max}] = [25, 72]$, $\kappa = 5$, $\mu = 0.86$. In turn, the SA heuristic was run under the following configurations: $\tau = 0.98$, $\phi = 20$, $\chi_0 = 0.9$ and removal_percentage = 0.75. All the parameter choices were made with

<footnote n="1">¹ https://www.ibm.com/analytics/cplex-optimizer</footnote>
the assistance of the automatic parameter tuning tool called Irace [56]. We also highlight that the implementation of the SA heuristic was done as to allow multi-threading execution. Precisely, the iterative code block from lines 3–14, Fig. 4 (in particular, the operations of lines 4 and 13), was implemented in a thread-safe way, as to be executed in concurrent threads.

## 9. Computational results

The computational experiments were performed on a 64-bit AMD EPYC-Milan virtual machine with 100.0 GB of RAM and running at 2.4 GHz under the Linux operating system. The machine has 16 cores, with hyper-threading turned off.

### 9.1. Benchmark of instances

In our experiments, we adapted the R-TOP-DP instances of $|N\setminus\{0\}| = \{25, 50, 100\}$ vertices of [11] by randomly selecting a subset of $[5\%|N\setminus\{0\}|]$ vertices as mandatory. The small percentage of mandatory vertices was chosen so as to avoid infeasible instances with respect to route duration. The instances consider three types of data sets, namely C, R and RC, where vertices are, respectively, clustered, randomly distributed and in a mix of the previous two cases. The first node is considered the origin, and each vertex is given a 2-dimensional position coordinate and a priority. Accordingly, the travel time of an arc is set to the Euclidean distance between the corresponding vertices. The routes' time limit $T$ assumes the values of the first, second, third and fourth quartiles of all vertices' deadlines of a given instance. We also considered a value of $T$ set to 80% of the first quartile to assess how the bound improvement incurred by the valid inequalities would behave for a lower routes' duration limit. The number of vehicles $|F|$ takes values from $\{2, 4\}$. At last, the worst-case service time deviation of a vertex $i \in P$ is set to $\hat{s}_i = \omega\bar{s}_i$, with $\omega \in \{10\%, 25\%\}$, and the uncertainty budget $\Gamma$ assumes values in $\{1, 5\}$. In total, considering all the possible parameter realizations mentioned above, we used 360 instances in our experiments.

### 9.2. Strength of the formulations and bound improvement by valid inequalities

In the first experiment, we compared the strength of formulations $\mathcal{F}_1$ and $\mathcal{F}_2$ and analyzed the impact of the inequalities discussed in Section 5. To this end, we computed the dual (upper) bounds obtained from solving the linear relaxations of these formulations, namely $\mathcal{L}_1$ and $\mathcal{L}_2$, with and without adding the AVICs and CCCs. Accordingly, for the models reinforced by the cuts (referred to as $\mathcal{L}_1$+cuts and $\mathcal{L}_2$+cuts), we solved the cutting-plane phase described in Section 5. $\mathcal{L}_1$ is used as baseline, as it directly extends the previously known model in the literature of R-TOP-DP [11].

In these experiments, we considered the whole benchmark of 360 instances. The instance sets are named $\langle t \rangle(|N\setminus\{0\}|)Q(q)$, with instance type $t \in \{C, R, RC\}$ and the quartile of vertex deadlines $q \in \{0.8, 1, 2, 3, 4\}$. Recall that the 0.8 case stands for 80% of the value of the first quartile. For each set, we considered all possible combinations of number of vehicles, service time deviation and uncertainty budget. The results are summarized in Table 1.

The first, second and third columns display, for each instance set, its name, the size of its instances and its cardinality. Then, for each (re-)inforced model, we give the average and the standard deviation (over all the instances in each set) of the percentage bound improvements obtained when compared to the baseline $\mathcal{L}_1$. Without loss of generality, given an instance, its percentage improvement in a configuration $i \in \{\mathcal{L}_1$+cuts, $\mathcal{L}_2$, $\mathcal{L}_2$+cuts\} is given by $100 \cdot \frac{UB_{\mathcal{L}_1} - UB_i}{UB_{\mathcal{L}_1}}$, where $UB_{\mathcal{L}_1}$ denotes the bound provided by $\mathcal{L}_1$, and $UB_i$ stands for the bound obtained for model $i$. The results are summarized in the “sub-total” row for each

instance size, and the last row displays the numerical results while considering the whole benchmark of instances. Bold entries highlight, for each instance set, the best model(s) in terms of average bound improvement.

The results clearly indicate the bound dominance of formulation $\mathcal{F}_2$ for the instances tested. In addition, the valid inequalities proposed were also capable of further strengthening both $\mathcal{F}_1$ and $\mathcal{F}_2$, especially for the instance sets with lower route time limits $T$ (see, e.g., $q = 0.8$). This behavior is expected, as these instances are more likely to have conflicting vertices (i.e., vertices that, when on a same route, exceed the maximum route duration limit $T$), thus yielding more CCCs. On average, the bound improvement achieved by $\mathcal{L}_2$ when compared to $\mathcal{L}_1$ was of 12.29%. Moreover, the average bound improvement of $\mathcal{L}_1$+cuts over $\mathcal{L}_1$ was of 4.30%, while that of $\mathcal{L}_2$+cuts over $\mathcal{L}_2$ was of 5.67%. We also highlight that, in a instance-per-instance analysis, the results show that both variations of $\mathcal{F}_2$ always provide strictly better bounds than those of $\mathcal{F}_1$ alone.

### 9.3. Solving the integer models

In order to assess how challenging the original problem is in practice, we set CPLEX to solve formulations $\mathcal{F}_1$ and $\mathcal{F}_2$ under default configuration, with a time limit of 3600 s, wall-clock time, for each instance. The results are shown in Table 2.

The first and second columns display the name of each instance set and the size of its instances, respectively. Then, for each algorithm, we give four result values described as follows. The first one corresponds to the number of instances solved at optimality (or to proven infeasibility) out of the complete instance set. The second one is the average wall-clock processing time (in seconds) spent. The last couple of result values provides the average and the standard deviation of the relative optimality gaps obtained by the algorithms. These gaps are given by $\frac{UB-LB}{UB}$, where $LB$ and $UB$ are, respectively, the best lower and upper bounds obtained by the corresponding algorithm for a given instance within one hour of execution. If, for a given instance, no feasible solution is found within the time limit and its infeasibility is also not proven, the corresponding optimality gap is assumed to be 100%. Likewise, this gap is set to 0% whenever the instance is proven to be infeasible. For the cases of $UB = LB = 0$, the optimality gaps are also set to 0%. The results are summarized in the “sub-total” row for each instance size, and the last row gives the overall results considering the 360 instances tested. Bold entries highlight, for each instance set, the best algorithm(s) in terms of number of instances solved to optimality and/or optimality gaps.

From the results, $\mathcal{F}_2$ outperforms $\mathcal{F}_1$ for all instance sets both in terms of number of instances solved to optimality (or proven infeasibility) as in optimality gaps. This behavior is in accordance with the fact that $\mathcal{F}_2$ gives much stronger bounds than $\mathcal{F}_1$ (as shown in Table 1). In summary, $\mathcal{F}_1$ and $\mathcal{F}_2$ solved to optimality, respectively, 16 and 35 out of the 360 instances, and obtained average gaps of 67.19% and 58.94%, respectively. In terms of execution times, the average values for the two formulations are comparable, with a same order of magnitude.

### 9.4. Heuristics

We also performed experiments to assess the performance of the proposed heuristics, namely Kernel Pump Heuristic (KPH) and Simulated Annealing Heuristic (SAH). KPH was run once for each instance, while SAH considers a unique set of 10 seeds used in independent executions, as it is built under randomized choices. Both algorithms were run in multi-threading mode, and the results are summarized in Table 3.

For each algorithm and instance set, we report five values: (i) the size of the instances in the set, (ii) the number of instances for which the heuristic was able to reach or improve over the best-known primal bounds, (iii) the average wall-clock execution times – in the case of
<caption>**Table 1**</caption>

Average percentage bound improvements obtained by the linear relaxation of (i) the vertex-flow formulation reinforced by the cuts of Section 5 ($\mathcal{L}_1$+cuts), (ii) the decreasing arc-flow formulation alone ($\mathcal{L}_2$) and with the additional cuts ($\mathcal{L}_2$+cuts). For each instance set, it is displayed the size of its instances and the total number of instances considered (#). $\mathcal{L}_1$ is used as bound baseline and the bold entries highlight, for each instance set, the best model(s) in terms of average bound improvement.

<table><thead><tr><th rowspan="2">Set</th><th rowspan="2">|N\{0}|</th><th rowspan="2"># inst.</th><th colspan="2">L<sub>1</sub>+cuts</th><th colspan="2">L<sub>2</sub></th><th colspan="2">L<sub>2</sub>+cuts</th></tr><tr><th>Avg (%)</th><th>StDev (%)</th><th>Avg (%)</th><th>StDev (%)</th><th>Avg (%)</th><th>StDev (%)</th></tr></thead><tbody><tr><td>C25Q0.8</td><td>25</td><td>8</td><td>4.53</td><td>8.39</td><td>10.27</td><td>0.88</td><td>18.74</td><td>9.92</td></tr><tr><td>C25Q1</td><td>25</td><td>8</td><td>0.00</td><td>0.00</td><td>10.18</td><td>0.82</td><td>10.20</td><td>0.82</td></tr><tr><td>C25Q2</td><td>25</td><td>8</td><td>0.00</td><td>0.00</td><td>10.16</td><td>0.81</td><td>10.20</td><td>0.82</td></tr><tr><td>C25Q3</td><td>25</td><td>8</td><td>0.00</td><td>0.00</td><td>10.15</td><td>0.81</td><td>10.19</td><td>0.82</td></tr><tr><td>C25Q4</td><td>25</td><td>8</td><td>0.00</td><td>0.00</td><td>10.15</td><td>0.81</td><td>10.19</td><td>0.82</td></tr><tr><td>R25Q0.8</td><td>25</td><td>8</td><td>14.09</td><td>12.91</td><td>13.71</td><td>1.06</td><td>37.78</td><td>11.60</td></tr><tr><td>R25Q1</td><td>25</td><td>8</td><td>0.00</td><td>0.00</td><td>13.47</td><td>1.00</td><td>13.48</td><td>1.01</td></tr><tr><td>R25Q2</td><td>25</td><td>8</td><td>0.00</td><td>0.00</td><td>13.44</td><td>1.00</td><td>13.46</td><td>1.00</td></tr><tr><td>R25Q3</td><td>25</td><td>8</td><td>0.00</td><td>0.00</td><td>13.43</td><td>1.00</td><td>13.46</td><td>1.00</td></tr><tr><td>R25Q4</td><td>25</td><td>8</td><td>0.00</td><td>0.00</td><td>13.42</td><td>0.99</td><td>13.45</td><td>1.00</td></tr><tr><td>RC25Q0.8</td><td>25</td><td>8</td><td>0.00</td><td>0.00</td><td>10.49</td><td>0.79</td><td>10.52</td><td>0.78</td></tr><tr><td>RC25Q1</td><td>25</td><td>8</td><td>0.00</td><td>0.00</td><td>10.48</td><td>0.79</td><td>10.51</td><td>0.78</td></tr><tr><td>RC25Q2</td><td>25</td><td>8</td><td>0.00</td><td>0.00</td><td>10.47</td><td>0.79</td><td>10.51</td><td>0.78</td></tr><tr><td>RC25Q3</td><td>25</td><td>8</td><td>0.00</td><td>0.00</td><td>10.47</td><td>0.79</td><td>10.51</td><td>0.78</td></tr><tr><td>RC25Q4</td><td>25</td><td>8</td><td>0.00</td><td>0.00</td><td>10.46</td><td>0.79</td><td>10.50</td><td>0.79</td></tr><tr><td>Subtotal</td><td></td><td>120</td><td>1.24</td><td>5.21</td><td>11.38</td><td>1.72</td><td>13.58</td><td>7.85</td></tr><tr><td>C50Q0.8</td><td>50</td><td>8</td><td>6.25</td><td>6.68</td><td>11.23</td><td>0.45</td><td>23.55</td><td>5.98</td></tr><tr><td>C50Q1</td><td>50</td><td>8</td><td>0.00</td><td>0.00</td><td>11.21</td><td>0.44</td><td>11.22</td><td>0.44</td></tr><tr><td>C50Q2</td><td>50</td><td>8</td><td>0.00</td><td>0.00</td><td>11.20</td><td>0.44</td><td>11.21</td><td>0.44</td></tr><tr><td>C50Q3</td><td>50</td><td>8</td><td>0.00</td><td>0.00</td><td>11.19</td><td>0.44</td><td>11.20</td><td>0.44</td></tr><tr><td>C50Q4</td><td>50</td><td>8</td><td>0.00</td><td>0.00</td><td>11.19</td><td>0.44</td><td>11.20</td><td>0.44</td></tr><tr><td>R50Q0.8</td><td>50</td><td>8</td><td>35.78</td><td>6.44</td><td>14.48</td><td>1.20</td><td>50.49</td><td>5.69</td></tr><tr><td>R50Q1</td><td>50</td><td>8</td><td>2.04</td><td>3.86</td><td>13.40</td><td>0.52</td><td>16.39</td><td>3.69</td></tr><tr><td>R50Q2</td><td>50</td><td>8</td><td>0.00</td><td>0.00</td><td>13.37</td><td>0.52</td><td>13.40</td><td>0.52</td></tr><tr><td>R50Q3</td><td>50</td><td>8</td><td>0.00</td><td>0.00</td><td>13.36</td><td>0.51</td><td>13.39</td><td>0.52</td></tr><tr><td>R50Q4</td><td>50</td><td>8</td><td>0.00</td><td>0.00</td><td>13.35</td><td>0.51</td><td>13.39</td><td>0.52</td></tr><tr><td>RC50Q0.8</td><td>50</td><td>8</td><td>9.99</td><td>11.03</td><td>11.85</td><td>0.46</td><td>30.22</td><td>7.45</td></tr><tr><td>RC50Q1</td><td>50</td><td>8</td><td>0.00</td><td>0.00</td><td>11.81</td><td>0.44</td><td>11.83</td><td>0.45</td></tr><tr><td>RC50Q2</td><td>50</td><td>8</td><td>0.00</td><td>0.00</td><td>11.78</td><td>0.43</td><td>11.80</td><td>0.44</td></tr><tr><td>RC50Q3</td><td>50</td><td>8</td><td>0.00</td><td>0.00</td><td>11.77</td><td>0.43</td><td>11.80</td><td>0.44</td></tr><tr><td>RC50Q4</td><td>50</td><td>8</td><td>0.00</td><td>0.00</td><td>11.76</td><td>0.43</td><td>11.79</td><td>0.44</td></tr><tr><td>Subtotal</td><td></td><td>120</td><td>3.60</td><td>9.78</td><td>12.20</td><td>1.17</td><td>16.86</td><td>10.82</td></tr><tr><td>C100Q0.8</td><td>100</td><td>8</td><td>31.36</td><td>8.06</td><td>11.83</td><td>0.25</td><td>41.67</td><td>4.35</td></tr><tr><td>C100Q1</td><td>100</td><td>8</td><td>0.64</td><td>0.89</td><td>11.75</td><td>0.24</td><td>13.63</td><td>2.23</td></tr><tr><td>C100Q2</td><td>100</td><td>8</td><td>0.00</td><td>0.00</td><td>11.74</td><td>0.24</td><td>11.77</td><td>0.25</td></tr><tr><td>C100Q3</td><td>100</td><td>8</td><td>0.00</td><td>0.00</td><td>11.73</td><td>0.24</td><td>11.75</td><td>0.24</td></tr><tr><td>C100Q4</td><td>100</td><td>8</td><td>0.00</td><td>0.00</td><td>11.72</td><td>0.24</td><td>11.73</td><td>0.24</td></tr><tr><td>R100Q0.8</td><td>100</td><td>8</td><td>45.62</td><td>6.08</td><td>24.00</td><td>2.36</td><td>58.89</td><td>2.81</td></tr><tr><td>R100Q1</td><td>100</td><td>8</td><td>7.54</td><td>8.21</td><td>13.70</td><td>0.29</td><td>23.24</td><td>5.11</td></tr><tr><td>R100Q2</td><td>100</td><td>8</td><td>0.00</td><td>0.00</td><td>13.67</td><td>0.29</td><td>13.69</td><td>0.29</td></tr><tr><td>R100Q3</td><td>100</td><td>8</td><td>0.01</td><td>0.03</td><td>13.66</td><td>0.29</td><td>13.68</td><td>0.29</td></tr><tr><td>R100Q4</td><td>100</td><td>8</td><td>0.01</td><td>0.03</td><td>13.65</td><td>0.29</td><td>13.67</td><td>0.29</td></tr><tr><td>RC100Q0.8</td><td>100</td><td>8</td><td>31.67</td><td>5.93</td><td>12.52</td><td>0.27</td><td>45.34</td><td>5.98</td></tr><tr><td>RC100Q1</td><td>100</td><td>8</td><td>4.08</td><td>4.53</td><td>12.41</td><td>0.26</td><td>17.49</td><td>3.60</td></tr><tr><td>RC100Q2</td><td>100</td><td>8</td><td>0.00</td><td>0.00</td><td>12.39</td><td>0.25</td><td>12.42</td><td>0.26</td></tr><tr><td>RC100Q3</td><td>100</td><td>8</td><td>0.00</td><td>0.00</td><td>12.38</td><td>0.25</td><td>12.42</td><td>0.26</td></tr><tr><td>RC100Q4</td><td>100</td><td>8</td><td>0.00</td><td>0.00</td><td>12.38</td><td>0.25</td><td>12.42</td><td>0.26</td></tr><tr><td>Subtotal</td><td></td><td>120</td><td>8.06</td><td>15.04</td><td>13.30</td><td>3.03</td><td>20.92</td><td>14.80</td></tr><tr><td>Total</td><td></td><td>360</td><td>4.30</td><td>11.12</td><td>12.29</td><td>2.26</td><td>17.12</td><td>11.87</td></tr></tbody></table>

SAH, for each instance, it is considered the average time over the 10 executions –, (iv) the average and (v) the standard deviation of the relative bound improvements given by 100 · $\frac{LB-LB^*}{LB^*}$, where $LB$ is the primal bound (cost of the heuristic solution) obtained by the corresponding heuristic. In the case of SAH, $LB$ corresponds to the average (over the 10 executions) primal bound. Moreover, $LB^*$ is the best-known bound for the instance (obtained by either $F_1$ or $F_2$ within 3600 s of execution). For the cases where, for a given instance and heuristic, no feasible solution is provided by the exact formulations, but at least one feasible solution is obtained by the heuristic in any execution (or proven to be infeasible), the due improvements of (iv) are set to 100%. By its turn, if the heuristic fails to find a feasible solution (or prove the instance's infeasibility) while the exact formulations succeed, the corresponding improvements are set to -100%. In fact, in

any other case, negative improvements indicate that the heuristic gives worse primal bounds than the best-known ones. Moreover, when the instance is proven to be infeasible both by the exact formulations and the heuristic, the improvement is set to 0%, as well as when $LB^* = LB = 0$. Entries in bold discriminate the cases where each heuristic obtained the best performance in terms of number of best bounds found or quality of the solutions.

From the results, SAH clearly outperforms KSH in terms of solution quality and time execution for almost all the instance sets tested. In fact, SAH was able to reach or improve over the best-known bounds for all instances, obtaining an average bound improvement of 15.68% with a tiny average execution time of 2.52 s. We highlight that, although the performance of KSH was poorer than that of SAH, KSH still obtained decent results, as it improves or reaches the best-known solutions in
<caption>**Table 2**</caption>

Results for formulations $\mathcal{F}_1$ and $\mathcal{F}_2$ while directly solved by CPLEX with a time limit of 3600 s. Bold entries highlight, for each instance set, the best algorithm(s) in terms of number of instances solved to optimality and/or optimality gaps.

<table><thead><tr><th rowspan="2">Set</th><th rowspan="2">|N\{0}|</th><th colspan="3">F<sub>1</sub></th><th colspan="3">F<sub>2</sub></th></tr><tr><th>#opt/total</th><th>Time (s)</th><th>Gap (%)</th><th>#opt/total</th><th>Time (s)</th><th>Gap (%)</th></tr><tr><td></td><td></td><td></td><td></td><td>Avg</td><td></td><td></td><td>Avg</td></tr><tr><td></td><td></td><td></td><td></td><td>StDev</td><td></td><td></td><td>StDev</td></tr></thead><tbody><tr><td>C25Q0.8</td><td>25</td><td>0/8</td><td>3600.00</td><td>77.54</td><td>15.10</td><td>2/8</td><td>2858.11</td><td>48.02</td><td>31.46</td></tr><tr><td>C25Q1</td><td>25</td><td>0/8</td><td>3600.00</td><td>63.07</td><td>17.51</td><td>0/8</td><td>3600.00</td><td>59.87</td><td>19.15</td></tr><tr><td>C25Q2</td><td>25</td><td>0/8</td><td>3600.00</td><td>40.17</td><td>14.52</td><td>0/8</td><td>3600.00</td><td>35.30</td><td>15.99</td></tr><tr><td>C25Q3</td><td>25</td><td>0/8</td><td>3600.00</td><td>35.29</td><td>9.67</td><td>0/8</td><td>3600.00</td><td>29.54</td><td>9.88</td></tr><tr><td>C25Q4</td><td>25</td><td>0/8</td><td>3600.00</td><td>35.20</td><td>9.72</td><td>0/8</td><td>3600.00</td><td>29.46</td><td>9.90</td></tr><tr><td>R25Q0.8</td><td>25</td><td>1/8</td><td>3269.10</td><td>63.69</td><td>27.79</td><td>8/8</td><td>26.92</td><td>0.00</td><td>0.00</td></tr><tr><td>R25Q1</td><td>25</td><td>0/8</td><td>3600.00</td><td>66.50</td><td>14.05</td><td>4/8</td><td>2547.78</td><td>21.08</td><td>22.96</td></tr><tr><td>R25Q2</td><td>25</td><td>0/8</td><td>3600.00</td><td>43.29</td><td>16.29</td><td>0/8</td><td>3600.00</td><td>37.09</td><td>17.87</td></tr><tr><td>R25Q3</td><td>25</td><td>0/8</td><td>3600.00</td><td>39.37</td><td>12.65</td><td>0/8</td><td>3600.00</td><td>31.91</td><td>13.60</td></tr><tr><td>R25Q4</td><td>25</td><td>0/8</td><td>3600.00</td><td>38.19</td><td>11.99</td><td>0/8</td><td>3600.00</td><td>31.02</td><td>12.73</td></tr><tr><td>RC25Q0.8</td><td>25</td><td>0/8</td><td>3600.00</td><td>67.09</td><td>15.00</td><td>0/8</td><td>3600.00</td><td>63.13</td><td>16.28</td></tr><tr><td>RC25Q1</td><td>25</td><td>0/8</td><td>3600.00</td><td>54.60</td><td>17.41</td><td>0/8</td><td>3600.00</td><td>50.51</td><td>17.44</td></tr><tr><td>RC25Q2</td><td>25</td><td>0/8</td><td>3600.00</td><td>39.92</td><td>15.37</td><td>0/8</td><td>3600.00</td><td>33.50</td><td>16.44</td></tr><tr><td>RC25Q3</td><td>25</td><td>0/8</td><td>3600.00</td><td>36.75</td><td>12.50</td><td>0/8</td><td>3600.00</td><td>30.24</td><td>12.76</td></tr><tr><td>RC25Q4</td><td>25</td><td>0/8</td><td>3600.00</td><td>35.20</td><td>11.15</td><td>0/8</td><td>3600.00</td><td>29.20</td><td>11.61</td></tr><tr><td>Sub-total</td><td></td><td>1/120</td><td>3577.94</td><td>49.06</td><td>20.24</td><td>14/120</td><td>3242.19</td><td>35.33</td><td>21.70</td></tr><tr><td>C50Q0.8</td><td>50</td><td>0/8</td><td>3600.00</td><td>90.49</td><td>6.23</td><td>0/8</td><td>3600.00</td><td>86.52</td><td>7.00</td></tr><tr><td>C50Q1</td><td>50</td><td>0/8</td><td>3600.00</td><td>86.29</td><td>7.32</td><td>0/8</td><td>3600.00</td><td>81.68</td><td>8.66</td></tr><tr><td>C50Q2</td><td>50</td><td>0/8</td><td>3600.00</td><td>78.59</td><td>9.48</td><td>0/8</td><td>3600.00</td><td>75.26</td><td>11.08</td></tr><tr><td>C50Q3</td><td>50</td><td>0/8</td><td>3600.00</td><td>60.52</td><td>13.12</td><td>0/8</td><td>3600.00</td><td>55.59</td><td>14.70</td></tr><tr><td>C50Q4</td><td>50</td><td>0/8</td><td>3600.00</td><td>56.72</td><td>9.03</td><td>0/8</td><td>3600.00</td><td>51.51</td><td>10.63</td></tr><tr><td>R50Q0.8</td><td>50</td><td>1/8</td><td>3150.10</td><td>81.09</td><td>33.19</td><td>6/8</td><td>1100.08</td><td>15.52</td><td>28.89</td></tr><tr><td>R50Q1</td><td>50</td><td>0/8</td><td>3600.00</td><td>88.55</td><td>7.51</td><td>0/8</td><td>3600.00</td><td>83.88</td><td>8.78</td></tr><tr><td>R50Q2</td><td>50</td><td>0/8</td><td>3600.00</td><td>72.11</td><td>11.30</td><td>0/8</td><td>3600.00</td><td>68.05</td><td>12.14</td></tr><tr><td>R50Q3</td><td>50</td><td>0/8</td><td>3600.00</td><td>64.11</td><td>12.35</td><td>0/8</td><td>3600.00</td><td>60.45</td><td>13.16</td></tr><tr><td>R50Q4</td><td>50</td><td>0/8</td><td>3600.00</td><td>58.53</td><td>9.57</td><td>0/8</td><td>3600.00</td><td>52.12</td><td>11.16</td></tr><tr><td>RC50Q0.8</td><td>50</td><td>0/8</td><td>3600.00</td><td>92.55</td><td>5.08</td><td>1/8</td><td>3463.24</td><td>69.61</td><td>30.35</td></tr><tr><td>RC50Q1</td><td>50</td><td>0/8</td><td>3600.00</td><td>85.76</td><td>7.30</td><td>0/8</td><td>3600.00</td><td>83.15</td><td>8.73</td></tr><tr><td>RC50Q2</td><td>50</td><td>0/8</td><td>3600.00</td><td>69.31</td><td>12.74</td><td>0/8</td><td>3600.00</td><td>63.74</td><td>14.00</td></tr><tr><td>RC50Q3</td><td>50</td><td>0/8</td><td>3600.00</td><td>61.40</td><td>12.18</td><td>0/8</td><td>3600.00</td><td>57.70</td><td>15.15</td></tr><tr><td>RC50Q4</td><td>50</td><td>0/8</td><td>3600.00</td><td>57.34</td><td>10.00</td><td>0/8</td><td>3600.00</td><td>51.46</td><td>11.38</td></tr><tr><td>Sub-total</td><td></td><td>1/120</td><td>3570.01</td><td>73.56</td><td>17.61</td><td>7/120</td><td>3424.22</td><td>63.75</td><td>22.80</td></tr><tr><td>C100Q0.8</td><td>100</td><td>4/8</td><td>1801.40</td><td>48.05</td><td>51.38</td><td>4/8</td><td>1801.28</td><td>47.68</td><td>51.02</td></tr><tr><td>C100Q1</td><td>100</td><td>0/8</td><td>3600.00</td><td>95.96</td><td>3.71</td><td>0/8</td><td>3600.00</td><td>95.19</td><td>5.11</td></tr><tr><td>C100Q2</td><td>100</td><td>0/8</td><td>3600.00</td><td>93.22</td><td>4.43</td><td>0/8</td><td>3600.00</td><td>90.84</td><td>6.18</td></tr><tr><td>C100Q3</td><td>100</td><td>0/8</td><td>3600.00</td><td>84.78</td><td>6.29</td><td>0/8</td><td>3600.00</td><td>83.26</td><td>9.46</td></tr><tr><td>C100Q4</td><td>100</td><td>0/8</td><td>3600.00</td><td>75.41</td><td>7.86</td><td>0/8</td><td>3600.00</td><td>82.34</td><td>12.91</td></tr><tr><td>R100Q0.8</td><td>100</td><td>4/8</td><td>1801.40</td><td>48.62</td><td>51.99</td><td>4/8</td><td>1801.44</td><td>44.15</td><td>47.51</td></tr><tr><td>R100Q1</td><td>100</td><td>1/8</td><td>3150.50</td><td>84.89</td><td>34.41</td><td>1/8</td><td>3150.54</td><td>83.21</td><td>33.91</td></tr><tr><td>R100Q2</td><td>100</td><td>0/8</td><td>3600.00</td><td>92.28</td><td>4.60</td><td>0/8</td><td>3600.00</td><td>90.64</td><td>6.87</td></tr><tr><td>R100Q3</td><td>100</td><td>0/8</td><td>3600.00</td><td>88.05</td><td>6.27</td><td>0/8</td><td>3600.00</td><td>82.84</td><td>8.67</td></tr><tr><td>R100Q4</td><td>100</td><td>0/8</td><td>3600.00</td><td>82.68</td><td>7.73</td><td>0/8</td><td>3600.00</td><td>79.68</td><td>11.00</td></tr><tr><td>RC100Q0.8</td><td>100</td><td>4/8</td><td>1801.49</td><td>48.63</td><td>51.99</td><td>4/8</td><td>1801.37</td><td>48.26</td><td>51.62</td></tr><tr><td>RC100Q1</td><td>100</td><td>1/8</td><td>3150.52</td><td>84.76</td><td>34.36</td><td>1/8</td><td>3150.54</td><td>83.88</td><td>34.16</td></tr><tr><td>RC100Q2</td><td>100</td><td>0/8</td><td>3600.00</td><td>90.49</td><td>6.62</td><td>0/8</td><td>3600.00</td><td>87.12</td><td>7.11</td></tr><tr><td>RC100Q3</td><td>100</td><td>0/8</td><td>3600.00</td><td>85.89</td><td>7.96</td><td>0/8</td><td>3600.00</td><td>83.07</td><td>7.92</td></tr><tr><td>RC100Q4</td><td>100</td><td>0/8</td><td>3600.00</td><td>80.43</td><td>10.39</td><td>0/8</td><td>3600.00</td><td>83.94</td><td>13.21</td></tr><tr><td>Sub-total</td><td></td><td>14/120</td><td>3180.35</td><td>78.94</td><td>29.99</td><td>14/120</td><td>3180.34</td><td>77.74</td><td>29.79</td></tr><tr><td>Total</td><td></td><td>16/360</td><td>3442.77</td><td>67.19</td><td>26.58</td><td>35/360</td><td>3282.25</td><td>58.94</td><td>30.58</td></tr></tbody></table>

294 out of the 360 instances and attains a significant average bound improvement of 13.31% within less than five minutes of average time. Then, taking into account its simplicity in terms of implementation effort, KSH stands as a good option when a slight solution deterioration is not a big issue.

## 9.5. Cut-and-branch algorithms

We also assessed the performance of the cut-and-branch described in Section 5, running for up to 3600 s. The resulting algorithms, referred to as CBA₁ and CBA₂, solve formulations $\mathcal{F}_1$ and $\mathcal{F}_2$, respectively, and

receive as warm-start the primal solutions obtained by the SA heuristic introduced in this work (Section 7). Table 4 summarizes the results for the cut-and-branch algorithms according to the same structure provided in Table 2 for the compact models.

In general, the performances of CBA₁ and CBA₂ follow the same pattern observed by $\mathcal{F}_1$ and $\mathcal{F}_2$, with a limited improvement on the quality of the bounds obtained. In particular, while compared to $\mathcal{F}_2$, CBA₂ was able to decrease the average optimality gap from 58.94% to 55.23% while solving to optimality one more instance. Nevertheless, we point out that the number of instances unsolved, as well as the
<caption>**Table 3**</caption>

Summary of the results obtained by the Kernel Search Heuristic (KSH) and the Simulated Annealing Heuristic (SAH) while running in multi-threading mode. Bold entries highlight, for each instance set, the best algorithm(s) in terms of average bound improvements and/or number of instances that reached or improved over the best-known bounds.

<table><thead><tr><th rowspan="3">Set</th><th rowspan="3">|N\{0}|</th><th colspan="4">Kernel Search Heuristic (KSH)</th><th colspan="4">Simulated Annealing Heuristic (SAH)</th></tr><tr><th rowspan="2">#best/total</th><th rowspan="2">Time (s)</th><th colspan="2">Improv. (%)</th><th rowspan="2">#best/total</th><th rowspan="2">Time (s)</th><th colspan="2">Improv. (%)</th></tr><tr><th>Avg</th><th>StDev</th><th>Avg</th><th>StDev</th></tr></thead><tbody><tr><td>C25Q0.8</td><td>25</td><td>8/8</td><td>44.70</td><td>0.00</td><td>0.00</td><td>8/8</td><td>0.11</td><td>0.00</td><td>0.00</td></tr><tr><td>C25Q1</td><td>25</td><td>8/8</td><td>97.16</td><td>0.27</td><td>0.77</td><td>8/8</td><td>0.20</td><td>0.27</td><td>0.77</td></tr><tr><td>C25Q2</td><td>25</td><td>3/8</td><td>187.23</td><td>-0.21</td><td>0.46</td><td>8/8</td><td>0.98</td><td>0.17</td><td>0.26</td></tr><tr><td>C25Q3</td><td>25</td><td>5/8</td><td>182.91</td><td>-0.03</td><td>0.04</td><td>8/8</td><td>1.58</td><td>0.06</td><td>0.16</td></tr><tr><td>C25Q4</td><td>25</td><td>3/8</td><td>185.22</td><td>-0.13</td><td>0.20</td><td>8/8</td><td>1.74</td><td>0.01</td><td>0.03</td></tr><tr><td>R25Q0.8</td><td>25</td><td>4/8</td><td>5.46</td><td>-2.35</td><td>5.42</td><td>8/8</td><td>0.10</td><td>0.00</td><td>0.00</td></tr><tr><td>R25Q1</td><td>25</td><td>6/8</td><td>55.69</td><td>-1.49</td><td>2.76</td><td>8/8</td><td>0.15</td><td>0.00</td><td>0.01</td></tr><tr><td>R25Q2</td><td>25</td><td>2/8</td><td>185.44</td><td>-0.61</td><td>0.97</td><td>8/8</td><td>0.58</td><td>0.34</td><td>0.78</td></tr><tr><td>R25Q3</td><td>25</td><td>3/8</td><td>193.05</td><td>-0.52</td><td>0.73</td><td>8/8</td><td>0.92</td><td>0.21</td><td>0.20</td></tr><tr><td>R25Q4</td><td>25</td><td>3/8</td><td>195.30</td><td>-0.45</td><td>0.89</td><td>8/8</td><td>1.30</td><td>0.22</td><td>0.36</td></tr><tr><td>RC25Q0.8</td><td>25</td><td>4/8</td><td>78.18</td><td>-0.53</td><td>0.88</td><td>8/8</td><td>0.15</td><td>0.12</td><td>0.28</td></tr><tr><td>RC25Q1</td><td>25</td><td>6/8</td><td>160.31</td><td>0.87</td><td>3.95</td><td>8/8</td><td>0.31</td><td>1.71</td><td>3.50</td></tr><tr><td>RC25Q2</td><td>25</td><td>5/8</td><td>191.48</td><td>0.10</td><td>0.78</td><td>8/8</td><td>0.97</td><td>0.69</td><td>1.05</td></tr><tr><td>RC25Q3</td><td>25</td><td>3/8</td><td>190.69</td><td>-0.96</td><td>1.40</td><td>8/8</td><td>1.28</td><td>0.17</td><td>0.28</td></tr><tr><td>RC25Q4</td><td>25</td><td>4/8</td><td>192.43</td><td>-0.09</td><td>0.14</td><td>8/8</td><td>1.65</td><td>0.01</td><td>0.02</td></tr><tr><td>Sub-total</td><td></td><td>67/120</td><td>143.02</td><td>-0.41</td><td>2.01</td><td>120/120</td><td>0.80</td><td>0.26</td><td>1.03</td></tr><tr><td>C50Q0.8</td><td>50</td><td>6/8</td><td>84.52</td><td>1.64</td><td>10.80</td><td>8/8</td><td>0.56</td><td>5.54</td><td>5.57</td></tr><tr><td>C50Q1</td><td>50</td><td>7/8</td><td>256.24</td><td>11.22</td><td>11.14</td><td>8/8</td><td>0.64</td><td>11.66</td><td>10.29</td></tr><tr><td>C50Q2</td><td>50</td><td>8/8</td><td>307.23</td><td>21.47</td><td>8.53</td><td>8/8</td><td>0.96</td><td>24.21</td><td>6.42</td></tr><tr><td>C50Q3</td><td>50</td><td>7/8</td><td>315.89</td><td>8.27</td><td>8.88</td><td>8/8</td><td>2.68</td><td>9.95</td><td>9.09</td></tr><tr><td>C50Q4</td><td>50</td><td>7/8</td><td>318.13</td><td>1.75</td><td>1.74</td><td>8/8</td><td>4.91</td><td>3.26</td><td>1.99</td></tr><tr><td>R50Q0.8</td><td>50</td><td>5/8</td><td>12.52</td><td>-2.53</td><td>3.79</td><td>8/8</td><td>0.49</td><td>0.02</td><td>0.03</td></tr><tr><td>R50Q1</td><td>50</td><td>8/8</td><td>90.37</td><td>4.47</td><td>3.32</td><td>8/8</td><td>0.54</td><td>4.55</td><td>3.25</td></tr><tr><td>R50Q2</td><td>50</td><td>8/8</td><td>310.32</td><td>11.51</td><td>6.41</td><td>8/8</td><td>1.07</td><td>14.63</td><td>6.61</td></tr><tr><td>R50Q3</td><td>50</td><td>8/8</td><td>314.63</td><td>8.68</td><td>6.54</td><td>8/8</td><td>1.77</td><td>11.20</td><td>6.87</td></tr><tr><td>R50Q4</td><td>50</td><td>7/8</td><td>314.41</td><td>1.52</td><td>1.42</td><td>8/8</td><td>3.89</td><td>2.80</td><td>1.41</td></tr><tr><td>RC50Q0.8</td><td>50</td><td>6/8</td><td>44.82</td><td>0.80</td><td>3.94</td><td>8/8</td><td>0.50</td><td>2.09</td><td>3.89</td></tr><tr><td>RC50Q1</td><td>50</td><td>7/8</td><td>164.97</td><td>4.90</td><td>13.41</td><td>8/8</td><td>0.59</td><td>9.73</td><td>8.51</td></tr><tr><td>RC50Q2</td><td>50</td><td>8/8</td><td>317.47</td><td>14.43</td><td>11.22</td><td>8/8</td><td>1.53</td><td>19.16</td><td>14.62</td></tr><tr><td>RC50Q3</td><td>50</td><td>8/8</td><td>316.11</td><td>8.05</td><td>6.17</td><td>8/8</td><td>2.44</td><td>10.17</td><td>5.45</td></tr><tr><td>RC50Q4</td><td>50</td><td>7/8</td><td>317.74</td><td>1.74</td><td>1.54</td><td>8/8</td><td>4.75</td><td>2.86</td><td>1.33</td></tr><tr><td>Sub-total</td><td></td><td>107/120</td><td>232.36</td><td>6.53</td><td>9.42</td><td>120/120</td><td>1.82</td><td>8.79</td><td>9.21</td></tr><tr><td>C100Q0.8</td><td>100</td><td>8/8</td><td>162.70</td><td>25.54</td><td>46.40</td><td>8/8</td><td>2.80</td><td>32.29</td><td>64.08</td></tr><tr><td>C100Q1</td><td>100</td><td>8/8</td><td>463.96</td><td>17.58</td><td>16.82</td><td>8/8</td><td>2.79</td><td>13.52</td><td>16.54</td></tr><tr><td>C100Q2</td><td>100</td><td>8/8</td><td>581.04</td><td>29.79</td><td>17.48</td><td>8/8</td><td>3.03</td><td>35.17</td><td>17.51</td></tr><tr><td>C100Q3</td><td>100</td><td>8/8</td><td>581.26</td><td>50.26</td><td>16.13</td><td>8/8</td><td>4.88</td><td>59.11</td><td>17.53</td></tr><tr><td>C100Q4</td><td>100</td><td>8/8</td><td>682.97</td><td>20.76</td><td>14.51</td><td>8/8</td><td>12.64</td><td>25.86</td><td>16.01</td></tr><tr><td>R100Q0.8</td><td>100</td><td>8/8</td><td>104.69</td><td>6.67</td><td>15.27</td><td>8/8</td><td>2.83</td><td>8.27</td><td>15.05</td></tr><tr><td>R100Q1</td><td>100</td><td>8/8</td><td>257.91</td><td>14.36</td><td>21.35</td><td>8/8</td><td>2.75</td><td>2.70</td><td>43.44</td></tr><tr><td>R100Q2</td><td>100</td><td>8/8</td><td>601.06</td><td>53.37</td><td>26.53</td><td>8/8</td><td>3.27</td><td>66.69</td><td>32.35</td></tr><tr><td>R100Q3</td><td>100</td><td>8/8</td><td>592.80</td><td>52.68</td><td>26.09</td><td>8/8</td><td>4.27</td><td>62.88</td><td>30.87</td></tr><tr><td>R100Q4</td><td>100</td><td>8/8</td><td>608.11</td><td>56.72</td><td>52.33</td><td>8/8</td><td>10.74</td><td>64.00</td><td>56.28</td></tr><tr><td>RC100Q0.8</td><td>100</td><td>8/8</td><td>121.22</td><td>17.85</td><td>22.78</td><td>8/8</td><td>2.82</td><td>18.33</td><td>23.32</td></tr><tr><td>RC100Q1</td><td>100</td><td>8/8</td><td>333.38</td><td>30.15</td><td>37.19</td><td>8/8</td><td>2.72</td><td>18.33</td><td>25.40</td></tr><tr><td>RC100Q2</td><td>100</td><td>8/8</td><td>594.83</td><td>48.53</td><td>20.61</td><td>8/8</td><td>3.57</td><td>62.61</td><td>23.86</td></tr><tr><td>RC100Q3</td><td>100</td><td>8/8</td><td>585.90</td><td>51.12</td><td>20.67</td><td>8/8</td><td>4.62</td><td>60.61</td><td>24.21</td></tr><tr><td>RC100Q4</td><td>100</td><td>8/8</td><td>580.62</td><td>31.93</td><td>24.32</td><td>8/8</td><td>10.34</td><td>39.30</td><td>28.64</td></tr><tr><td>Sub-total</td><td></td><td>120/120</td><td>456.83</td><td>33.82</td><td>30.57</td><td>120/120</td><td>4.94</td><td>37.98</td><td>37.67</td></tr><tr><td>Total</td><td></td><td>294/360</td><td>277.40</td><td>13.31</td><td>23.65</td><td>360/360</td><td>2.52</td><td>15.68</td><td>27.57</td></tr></tbody></table>

optimality gaps, are rather high, indicating that the problem is indeed challenging in practice, even for small to medium-sized instances of 25 and 50 vertices.

## 9.6. The price of robustness

As a last experiment, we compared the average solution values (in terms of the sum of the collected priorities) segmented by instance size, route duration limit ($T$) and uncertainty budget ($\Gamma$). The goal was to measure the impact on the priority sum incurred by restricting

the duration limit for the routes and/or increasing the uncertainty allowed. To this end, we considered the average solution values provided by our simulated annealing heuristic SAH running for the complete benchmark (10 seeds/executions for each instance, as in Section 9.4), and, additionally, for the corresponding deterministic cases where the uncertainty budget is limited to zero. The summary of the results is shown in Fig. 9.

The results confirm the intuition that looser route duration limits yield solutions with a greater sum of collected priorities, since it is usually possible to visit more vertices. Likewise, increased robustness
<caption>**Table 4**</caption>

Results for the cut-and-branch algorithms CBA₁ and CBA₂ within 3600 s of execution. Bold entries highlight, for each instance set, the best algorithm(s) in terms of number of instances solved to optimality and/or optimality gaps.

<table><thead><tr><th rowspan="3">Set</th><th rowspan="3">|N\{0}|</th><th colspan="3">CBA<sub>1</sub></th><th colspan="3">CBA<sub>2</sub></th></tr><tr><th rowspan="2">#opt/total</th><th rowspan="2">Time (s)</th><th colspan="2">Gap (%)</th><th rowspan="2">#opt/total</th><th rowspan="2">Time (s)</th><th colspan="2">Gap (%)</th></tr><tr><th>Avg</th><th>StDev</th><th>Avg</th><th>StDev</th></tr></thead><tbody><tr><td>C25Q0.8</td><td>25</td><td>0/8</td><td>3600.00</td><td>76.83</td><td>14.73</td><td>2/8</td><td>2841.64</td><td>48.76</td><td>31.22</td></tr><tr><td>C25Q1</td><td>25</td><td>0/8</td><td>3600.00</td><td>62.75</td><td>17.73</td><td>0/8</td><td>3600.00</td><td>59.41</td><td>18.98</td></tr><tr><td>C25Q2</td><td>25</td><td>0/8</td><td>3600.00</td><td>39.88</td><td>14.75</td><td>0/8</td><td>3600.00</td><td>34.57</td><td>15.37</td></tr><tr><td>C25Q3</td><td>25</td><td>0/8</td><td>3600.00</td><td>35.27</td><td>9.80</td><td>0/8</td><td>3600.00</td><td>29.44</td><td>9.94</td></tr><tr><td>C25Q4</td><td>25</td><td>0/8</td><td>3600.00</td><td>35.20</td><td>9.77</td><td>0/8</td><td>3600.00</td><td>29.46</td><td>9.93</td></tr><tr><td>R25Q0.8</td><td>25</td><td>1/8</td><td>3321.30</td><td>60.91</td><td>25.81</td><td>8/8</td><td>28.46</td><td>0.00</td><td>0.00</td></tr><tr><td>R25Q1</td><td>25</td><td>0/8</td><td>3600.00</td><td>66.46</td><td>13.94</td><td>4/8</td><td>2524.31</td><td>23.21</td><td>24.84</td></tr><tr><td>R25Q2</td><td>25</td><td>0/8</td><td>3600.00</td><td>43.09</td><td>16.21</td><td>0/8</td><td>3600.00</td><td>36.09</td><td>17.66</td></tr><tr><td>R25Q3</td><td>25</td><td>0/8</td><td>3600.00</td><td>39.18</td><td>12.61</td><td>0/8</td><td>3600.00</td><td>31.63</td><td>13.43</td></tr><tr><td>R25Q4</td><td>25</td><td>0/8</td><td>3600.00</td><td>38.03</td><td>11.83</td><td>0/8</td><td>3600.00</td><td>30.83</td><td>12.52</td></tr><tr><td>RC25Q0.8</td><td>25</td><td>0/8</td><td>3600.00</td><td>66.39</td><td>15.29</td><td>0/8</td><td>3600.00</td><td>63.08</td><td>16.39</td></tr><tr><td>RC25Q1</td><td>25</td><td>0/8</td><td>3600.00</td><td>53.48</td><td>16.61</td><td>0/8</td><td>3600.00</td><td>49.09</td><td>17.67</td></tr><tr><td>RC25Q2</td><td>25</td><td>0/8</td><td>3600.00</td><td>38.94</td><td>15.22</td><td>0/8</td><td>3600.00</td><td>33.07</td><td>16.02</td></tr><tr><td>RC25Q3</td><td>25</td><td>0/8</td><td>3600.00</td><td>36.26</td><td>12.19</td><td>0/8</td><td>3600.00</td><td>30.09</td><td>12.73</td></tr><tr><td>RC25Q4</td><td>25</td><td>0/8</td><td>3600.00</td><td>35.10</td><td>11.25</td><td>0/8</td><td>3600.00</td><td>29.13</td><td>11.61</td></tr><tr><td>Sub-total</td><td></td><td>1/120</td><td>3581.42</td><td>48.52</td><td>19.84</td><td>14/120</td><td>3239.63</td><td>35.19</td><td>21.60</td></tr><tr><td>C50Q0.8</td><td>50</td><td>0/8</td><td>3600.00</td><td>87.34</td><td>7.03</td><td>0/8</td><td>3600.00</td><td>84.33</td><td>8.63</td></tr><tr><td>C50Q1</td><td>50</td><td>0/8</td><td>3600.00</td><td>81.17</td><td>9.46</td><td>0/8</td><td>3600.00</td><td>78.78</td><td>10.58</td></tr><tr><td>C50Q2</td><td>50</td><td>0/8</td><td>3600.00</td><td>72.21</td><td>11.50</td><td>0/8</td><td>3600.00</td><td>68.72</td><td>12.81</td></tr><tr><td>C50Q3</td><td>50</td><td>0/8</td><td>3600.00</td><td>57.01</td><td>11.20</td><td>0/8</td><td>3600.00</td><td>51.56</td><td>12.43</td></tr><tr><td>C50Q4</td><td>50</td><td>0/8</td><td>3600.00</td><td>55.12</td><td>9.12</td><td>0/8</td><td>3600.00</td><td>49.38</td><td>10.03</td></tr><tr><td>R50Q0.8</td><td>50</td><td>1/8</td><td>3150.11</td><td>75.52</td><td>31.73</td><td>6/8</td><td>1176.71</td><td>15.40</td><td>28.60</td></tr><tr><td>R50Q1</td><td>50</td><td>0/8</td><td>3600.00</td><td>85.69</td><td>8.29</td><td>0/8</td><td>3600.00</td><td>83.19</td><td>9.49</td></tr><tr><td>R50Q2</td><td>50</td><td>0/8</td><td>3600.00</td><td>67.71</td><td>12.87</td><td>0/8</td><td>3600.00</td><td>62.75</td><td>14.76</td></tr><tr><td>R50Q3</td><td>50</td><td>0/8</td><td>3600.00</td><td>60.09</td><td>12.38</td><td>0/8</td><td>3600.00</td><td>53.87</td><td>14.24</td></tr><tr><td>R50Q4</td><td>50</td><td>0/8</td><td>3600.00</td><td>57.12</td><td>10.00</td><td>0/8</td><td>3600.00</td><td>50.49</td><td>11.46</td></tr><tr><td>RC50Q0.8</td><td>50</td><td>0/8</td><td>3600.00</td><td>89.19</td><td>6.47</td><td>2/8</td><td>3415.93</td><td>57.10</td><td>37.24</td></tr><tr><td>RC50Q1</td><td>50</td><td>0/8</td><td>3600.00</td><td>82.88</td><td>9.07</td><td>0/8</td><td>3600.00</td><td>80.61</td><td>10.22</td></tr><tr><td>RC50Q2</td><td>50</td><td>0/8</td><td>3600.00</td><td>62.53</td><td>12.39</td><td>0/8</td><td>3600.00</td><td>57.55</td><td>13.97</td></tr><tr><td>RC50Q3</td><td>50</td><td>0/8</td><td>3600.00</td><td>57.59</td><td>11.60</td><td>0/8</td><td>3600.00</td><td>51.91</td><td>13.15</td></tr><tr><td>RC50Q4</td><td>50</td><td>0/8</td><td>3600.00</td><td>55.68</td><td>10.04</td><td>0/8</td><td>3600.00</td><td>49.83</td><td>11.28</td></tr><tr><td>Sub-total</td><td></td><td>1/120</td><td>3570.01</td><td>69.79</td><td>17.38</td><td>8/120</td><td>3426.18</td><td>59.70</td><td>23.22</td></tr><tr><td>C100Q0.8</td><td>100</td><td>4/8</td><td>1801.55</td><td>45.66</td><td>48.84</td><td>4/8</td><td>1801.41</td><td>44.74</td><td>47.85</td></tr><tr><td>C100Q1</td><td>100</td><td>0/8</td><td>3600.00</td><td>93.92</td><td>5.62</td><td>0/8</td><td>3600.00</td><td>93.08</td><td>6.35</td></tr><tr><td>C100Q2</td><td>100</td><td>0/8</td><td>3600.00</td><td>88.45</td><td>7.17</td><td>0/8</td><td>3600.00</td><td>86.90</td><td>8.11</td></tr><tr><td>C100Q3</td><td>100</td><td>0/8</td><td>3600.00</td><td>75.07</td><td>9.26</td><td>0/8</td><td>3600.00</td><td>71.70</td><td>10.46</td></tr><tr><td>C100Q4</td><td>100</td><td>0/8</td><td>3600.00</td><td>69.33</td><td>7.01</td><td>0/8</td><td>3600.00</td><td>65.20</td><td>7.87</td></tr><tr><td>R100Q0.8</td><td>100</td><td>4/8</td><td>1801.53</td><td>46.09</td><td>49.30</td><td>4/8</td><td>1801.41</td><td>44.64</td><td>47.77</td></tr><tr><td>R100Q1</td><td>100</td><td>1/8</td><td>3150.62</td><td>82.91</td><td>33.81</td><td>1/8</td><td>3150.55</td><td>81.86</td><td>33.56</td></tr><tr><td>R100Q2</td><td>100</td><td>0/8</td><td>3600.00</td><td>84.56</td><td>7.33</td><td>0/8</td><td>3600.00</td><td>82.17</td><td>8.45</td></tr><tr><td>R100Q3</td><td>100</td><td>0/8</td><td>3600.00</td><td>77.26</td><td>8.93</td><td>0/8</td><td>3600.00</td><td>73.69</td><td>10.25</td></tr><tr><td>R100Q4</td><td>100</td><td>0/8</td><td>3600.00</td><td>72.00</td><td>7.26</td><td>0/8</td><td>3600.00</td><td>67.70</td><td>8.32</td></tr><tr><td>RC100Q0.8</td><td>100</td><td>4/8</td><td>1801.38</td><td>46.69</td><td>49.93</td><td>4/8</td><td>1801.46</td><td>45.85</td><td>49.03</td></tr><tr><td>RC100Q1</td><td>100</td><td>1/8</td><td>3150.56</td><td>82.69</td><td>33.75</td><td>1/8</td><td>3150.55</td><td>81.90</td><td>33.56</td></tr><tr><td>RC100Q2</td><td>100</td><td>0/8</td><td>3600.00</td><td>82.52</td><td>8.18</td><td>0/8</td><td>3600.00</td><td>80.05</td><td>9.37</td></tr><tr><td>RC100Q3</td><td>100</td><td>0/8</td><td>3600.00</td><td>76.93</td><td>8.96</td><td>0/8</td><td>3600.00</td><td>73.68</td><td>10.17</td></tr><tr><td>RC100Q4</td><td>100</td><td>0/8</td><td>3600.00</td><td>72.80</td><td>7.24</td><td>0/8</td><td>3600.00</td><td>68.89</td><td>8.23</td></tr><tr><td>Sub-total</td><td></td><td>14/120</td><td>3180.38</td><td>73.13</td><td>28.67</td><td>14/120</td><td>3180.36</td><td>70.80</td><td>28.43</td></tr><tr><td>Total</td><td></td><td>16/360</td><td>3443.93</td><td>63.81</td><td>24.94</td><td>36/360</td><td>3282.05</td><td>55.23</td><td>28.69</td></tr></tbody></table>

in terms of uncertain visit times at vertices comes at the price of poorer priority sums. Therefore, when calibrating the level of robustness targeted, one must take into account not only the computational effort but also the degradation of the objective function to find an adequate trade-off of one's specific needs.

## 10. Concluding remarks

In this work, we proposed a new robust optimization variation of the team orienteering problem that considers mandatory vertices and decreasing priority values for visiting locations/vertices. Service times on visiting locations are subject to uncertainty, and a budget parameter is used to control the level of uncertainty. Two compact

formulations are presented, namely vertex-flow ($F_1$) and decreasing arc-flow ($F_2$), along with two classes of valid inequalities adapted from the literature. Computational results suggest the dominance of the decreasing arc-flow formulation over the other one in terms of strength of the linear relaxation. Precisely, the latter formulation was able to improve in 12.29%, on average, the bounds provided by the vertex-flow for a benchmark of 360 instances. When reinforced by the valid inequalities, the average improvement was 17.12%. Additionally, we tested the proposed models on solving the integer problems and the cut-and-branch algorithms obtained from coupling the formulations with the valid inequalities and a simulated annealing heuristic also developed in this work. The cut-and-branch algorithm based on $F_2$ provided better results than the one using $F_1$, being able to solve, in total, 20 more instances than the latter.
In terms of non-exact methods, we proposed a mat-heuristic based on kernel search (KSH) and a simulated annealing heuristic (SAH). KSH obtained decent results, by improving or reaching the best-known primal bounds (from either $\mathcal{F}_1$ or $\mathcal{F}_2$ running for up to one hour) for 294 out of the 360 instances tested, with an average execution time of less than five minutes. The multi-threading version of SAH achieved even better results, by improving or reaching the best-known bounds for all instances tested, while running for a tiny average time of 2.52 s. KSH stands as an easily reproducible option – as it only relies on the existence of an MILP model – when a slight solution deterioration is not problematic, while the more elaborate SAH should be applied when aiming to achieve the best possible solution quality with a fast response time.

Although, from the experiments, it is clear that the decreasing arc-flow formulation provides stronger bounds, a formal proof of the strength of the formulation is yet to be devised, which is an interesting future research direction. Moreover, the theoretical difficulty of the problem stands out in practice, as, for all the exact algorithms tested, the number of unsolved instances is still high, even for medium-sized instances of 50 vertices. Therefore, innovative decomposition approaches and reformulations should be investigated for the problem. In addition, a more practical and experimental study can be conducted to assess the proposed methods, especially SAH, at solving realistic problem instances built from real data.

So far, we have already put effort in reformulating the proposed models in a Benders' decomposition fashion, by dualizing the continuous feasible space of the complicating constraints (9)–(12) and (18)–(21) from formulations $\mathcal{F}_1$ and $\mathcal{F}_2$, respectively. Unfortunately, the convergence of the resulting Benders' decomposition algorithms was rather poor, even after attempts of selecting cuts in a wiser way [57]. Therefore, we suggest that future research take a complementary direction by exploring column-generation approaches.

## CRediT authorship contribution statement

**Lucas Assunção:** Writing – review & editing, Writing – original draft, Software, Project administration, Methodology, Investigation, Funding acquisition, Conceptualization. **Andréa Cynthia Santos:** Supervision, Resources, Conceptualization.

## Declaration of competing interest

The authors have declared no conflict of interest.

## Acknowledgments

This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 101034329. Recipient of the WINNINGNormandy Program supported by the Normandy Region.

## Data availability

Data will be made available on request.

## References

[1] Baffo I, Carotenuto P, Rondine S. An orienteering-based approach to manage emergency situation. Transp Res Procedia 2017;22:297-304.

[2] Duhamel C, Santos AC, Brasil D, Châtelet E, Birregah B. Connecting a population dynamic model with a multi-period location-allocation problem for post-disaster relief operations. Ann Oper Res 2016;247:693-713.

[3] Zheng Y-J, Chen S-Y, Ling H-F. Evolutionary optimization for disaster relief operations: A survey. Appl Soft Comput 2015;27:553-66.

[4] Bayram V. Optimization models for large scale network evacuation planning and management: A literature review. Surv Oper Res Manag Sci 2016;21(2):63-84.

[5] Caunhye AM, Nie X, Pokharel S. Optimization models in emergency logistics: A literature review. Socio- Econ Plan Sci 2012;46(1):4-13, Special Issue: Disaster Planning and Logistics: Part 1.

[6] Kamyabniya A, Sauré A, Salman FS, Bénichou N, Patrick J. Optimization models for disaster response operations: a literature review. OR Spectrum 2024;46(3):737-83.

[7] Coco AA, Santos AC, Noronha TF. Robust min-max regret covering problems. Comput Optim Appl 2022;83:111-41.

[8] Wang D, Yang K, Yang L, Dong J. Two-stage distributionally robust optimization for disaster relief logistics under option contract and demand ambiguity. Transportation Research Part E: Logistics Transportation Review 2023;170:103025.

[9] Murat Afsar H, Labadie N. Team orienteering problem with decreasing profits. Electron Notes Discrete Math 2013;41:285-93.

[11] Bertsimas D, Sim M. The price of robustness. Oper Res 2004;52(1):35-53.

[12] Yu Q, Cheng C, Zhu N. Robust team orienteering problem with decreasing profits. INFORMS J Comput 2022;34(6):3215-33.

[13] Faiz TI, Vogiatzis C, Liu J, Noor-E-Alam M. A robust optimization framework for two-echelon vehicle and UAV routing for post-disaster humanitarian logistics operations. Networks 2024;84(2):200-19.

[14] Li J, Che A, Chu F. Prepositioning of emergency supplies for predictable disasters using distributionally robust optimization. IFAC- Pap 2022;55(10):3100-5, 10th IFAC Conference on Manufacturing Modelling, Management and Control MIM 2022.

[15] Chao IM, Golden BL, Wasil EA. The team orienteering problem. European J Oper Res 1996;88(3):464-74.

[16] Assunção L, Mateus GR. A cutting-plane algorithm for the Steiner team orienteering problem. Comput Ind Eng 2019;135:922-39.

[17] Letchford AN, Nasiri SD, Theis DO. Compact formulations of the Steiner traveling salesman problem and related problems. European J Oper Res 2013;228(1):83-92.

[18] Angelelli E, Mansini R, Grazia Speranza M. Kernel search: A general heuristic for the multi-dimensional knapsack problem. Comput Oper Res 2010;37(11):2017-26.

[19] Kirkpatrick S, Gelatt CD, Vecchi MP. Optimization by simulated annealing. Science 1983;220(4598):671-80.

[20] Tsiligirides T. Heuristic methods applied to orienteering. J Oper Res Soc 1984;35(9):797-809.

[21] Dantzig GB, Fulkerson DR, Johnson SM. Solution of a large-scale traveling-salesman problem. Oper Res 1954;2:393-410.

[22] Laporte G, Martello S. The selective travelling salesman problem. Discrete Appl Math 1990;26(2-3):193-207.

[23] Kantor MG, Rosenwein MB. The orienteering problem with time windows. J Oper Res Soc 1992;43(6):629-35.

[24] Montemanni R, Gambardella L. An ant colony system for team orienteering problems with time windows. Found Comput Decision Sci 2009;34(4):287-306.

[25] Hanafi S, Mansini R, Zanotti R. The multi-visit team orienteering problem with precedence constraints. European J Oper Res 2020;282(2):515-29.

[26] Chakraborty S, Agarwal YK. Solving the team orienteering problem with nonidentical agents: A Lagrangian approach. Networks 2021;78(4):443-60.

[27] Evers L, Dollevoet T, Barros AI, Monsuur H. Robust UAV mission planning. Ann Oper Res 2014;222:293-315.

[28] Angelelli E, Archetti C, Filippi C, Vindigni M. The probabilistic orienteering problem. Comput Oper Res 2017;81:269-81.

[29] Ilhan T, Iravani S, Daskin M. The orienteering problem with stochastic profits. IIE Trans 2008;40:406-21.

[30] Song Y, Ulmer MW, Thomas BW, Wallace SW. Building trust in home services—Stochastic team-orienteering with consistency constraints. Transp Sci 2020;54(3):823-38.

[31] Campbell AM, Gendreau M, Thomas BW. The orienteering problem with stochastic travel and service times. Ann Oper Res 2011;186:61-81.

[32] Evers L, Glorie K, van der Ster S, Barros AI, Monsuur H. A two-stage approach to the orienteering problem with stochastic weights. Comput Oper Res 2014;43:248-60.

[33] Tang H, Miller-Hooks E. Algorithms for a stochastic selective travelling salesperson problem. J Oper Res Soc 2005;56(4):439-52.

[34] Gunawan A, Lau HC, Vansteenwegen P. Orienteering problem: A survey of recent variants, solution approaches and applications. European J Oper Res 2016;255(2):315-32.

[35] Bektas T. The multiple traveling salesman problem: an overview of formulations and solution procedures. Omega 2006;34(3):209-19.

[36] Korte B, Vygen J. Combinatorial optimization: Theory and algorithms. fifth ed.. Springer Publishing Company, Incorporated; 2012.

[37] Agra A, Christiansen M, Figueiredo R, Hvattum LM, Poss M, Requejo C. The robust vehicle routing problem with time windows. Comput Oper Res 2013;40(3):856-66, Transport Scheduling.
[38] Lee C, Lee K. Robust vehicle routing problem with deadlines and travel time/demand uncertainty. J Oper Res Soc 2012;63:1294-306.

[39] Munari P, Moreno A, De La Vega J, Alem D, Gondzio J, Morabito R. The robust vehicle routing problem with time windows: Compact formulation and branch-price-and-cut method. Transp Sci 2019;53(4):1043-66.

[40] Miller CE, Tucker AW, Zemlin RA. Integer programming formulation of traveling salesman problems. J ACM 1960;7(4):326-9.

[41] Assunção L, Mateus GR. Coupling feasibility pump and large neighborhood search to solve the Steiner team orienteering problem. Comput Oper Res 2021;128:105175.

[42] Bianchessi N, Mansini R, Speranza MG. A branch-and-cut algorithm for the team orienteering problem. Int Trans Oper Res 2018;25(2):627-35.

[43] Tomita E, Tanaka A, Takahashi H. The worst-case time complexity for generating all maximal cliques and computational experiments. Theoret Comput Sci 2006;363(1):28-42.

[44] Guastaroba G, Savelsbergh M, Speranza M. Adaptive Kernel search: A heuristic for solving mixed integer linear programs. European J Oper Res 2017;263(3):789-804.

[45] Shaw P. Using constraint programming and local search methods to solve vehicle routing problems. In: Maher M, Puget J-F, editors. Principles and practice of constraint programming — CP98. Berlin, Heidelberg: Springer Berlin Heidelberg; 1998, p. 417-31.

[46] Dang DC, Guibadj RN, Moukrim A. An effective PSO-inspired algorithm for the team orienteering problem. European J Oper Res 2013;229(2):332-44.

[47] Ke L, Archetti C, Feng Z. Ants can solve the team orienteering problem. Comput Ind Eng 2008;54(3):648-65.

[48] Ke L, Zhai L, Li J, Chan FT. Pareto mimic algorithm: An approach to the team orienteering problem. Omega 2016;61:155-66.

[49] Kim BI, Li H, Johnson AL. An augmented large neighborhood search method for solving the team orienteering problem. Expert Syst Appl 2013;40(8):3065-72.

[50] Souffriau W, Vansteenwegen P, Berghe GV, Oudheusden DV. A path re-linking approach for the team orienteering problem. Comput Oper Res 2010;37(11):1853-9.

[51] Vansteenwegen P, Souffriau W, Berghe GV, Oudheusden DV. A guided local search metaheuristic for the team orienteering problem. European J Oper Res 2009;196(1):118-27.

[52] Ben-Ameur W. Computing the initial temperature of simulated annealing. Comput Optim Appl 2004;29:369-85.

[53] Lin S. Computer solutions of the traveling salesman problem. Bell Syst Tech J 1965;44(10):2245-69.

[54] Goldberg AV, Tarjan RE. A new approach to the maximum-flow problem. J ACM 1988;35(4):921-40.

[55] Dezső B, Jüttner A, Kovács P. LEMON — an open source C++ graph template library. Electron Notes Theor Comput Sci 2011;264(5):23-45, Proceedings of the Second Workshop on Generative Technologies (WGT) 2010.

[56] López-Ibáñez M, Dubois-Lacoste J, Pérez Cáceres L, Birattari M, Stützle T. The irace package: Iterated racing for automatic algorithm configuration. Oper Res Perspect 2016;3:43-58.

[57] Fischetti M, Salvagnin D, Zanette A. A note on the selection of Benders' cuts. Math Program 2010;124:175-82.