Dual-driven path elimination for vehicle routing with idle times and arrival-time consistency

Highlights Dual-driven path elimination for vehicle routing with idle times and arrival-time consistency Jorge Riera-Ledesma, Inmaculada Rodríguez-Martín, Hipólito Hernández-Pérez •A dual-driven technique for infeasible path elimination is introduced. •The method targets vehicle routing with idle times and arrival-time consistency. •It integrates infeasibility detection into a branch-and-cut framework. •Three formulations are proposed: a compact model, an SEC-enhanced formulation, and an IPEC-based model. •Instances with up to 100 customers and 5 periods are solved to optimality. Dual-driven path elimination for vehicle routing with idle times and arrival-time consistency Jorge Riera-Ledesmaa,c,∗, Inmaculada Rodríguez-Martínb,c, Hipólito Hernández-Pérezb,c aDepartamento de Ingeniería Informática y de Sistemas, Universidad de La Laguna, La Laguna, Spain bDepartamento de Matemáticas, Estadística e Investigación Operativa, Universidad de La Laguna, La Laguna, Spain cInstituto de Matemáticas y Aplicaciones, Universidad de La Laguna, La Laguna, Spain Abstract We present a simple dual-driven methodology for generating infeasible path elimination constraints within branch-and-cut algorithms for vehicle routing problems that incorporate idle times and arrival-time consistency requirements. By leveraging dual information from a feasibility-checking subproblem, the approach systematically identifies the combinatorial sources of infeasibility and uses them to generate and strengthen valid inequalities. We apply the method to the Consistent Traveling Salesperson Problem with idling, which enforces temporal consistency across multiple service days while allowing idle time between tasks. This problem, defined by basic routing and synchronization constraints, serves as an ideal case study to demonstrate the method’s effectiveness. Computational experiments on a benchmark set of 756 instances, based on multi-period extensions of classical TSPLIB datasets, show that the approach solves 536 instances to proven optimality, including cases with up to 100 customers and a five-day planning horizon, all within a two-hour time limit. Keywords: exact algorithms, vehicle routing, synchronization, branch-and-cut, path elimination, time consistency 1. Introduction An infeasible Path Elimination Constraint (iPEC) is a valid inequality that is dynamically introduced to exclude combinations of subpaths that are infeasible in any solution due to synchronization constraints, arrival-time consistency, or other problem-specific requirements arising in vehicle routing problems. These constraints play a central role in the bounding phase of branch-and-cut algorithms, where a master problem determines routing decisions while a subproblem verifies feasibility conditions to eliminate inconsistent partial solutions. IPECs are especially valuable when certain conditions are difficult to express linearly or when their explicit modeling would incur significant computational cost. Instead of relying on large penalty constants (the so-called Big-Mapproach), these inequalities provide a more efficient mechanism to rule out subroutes that implicitly violate global constraints, without resorting to more complex or less tractable formulations. Identifying the elements required to construct an iPEC is not always straightforward. In some cases, this task relies on a feasibility checking module that only reports whether a candidate solution is feasible, ∗Corresponding author Email addresses: [email protected] (Jorge Riera-Ledesma), [email protected] (Inmaculada Rodríguez-Martín), [email protected] (Hipólito Hernández-Pérez) Preprint submitted to Computers and Operations Research November 8, 2025 without revealing the specific cause of infeasibility. As a result, one only knows that the violation stems from the full combination of arcs defining the solution, with no further insight into the structure of the conflict. Moreover, since iPECs eliminate infeasible solutions by removing a single arc in the sequence responsible for the infeasibility, their effect on the quality of the Linear Programming (LP) relaxation is often weak. For this reason, it is common practice to strengthen such constraints to enhance algorithmic performance. This strengthening may involve identifying minimally infeasible subsets of arcs within the original solution, adding additional arcs to the left-hand side of the constraint, or combining both strategies. To do this effectively, it is essential to understand the underlying combinatorial structure of the arcs involved. Although narrower in scope than more general frameworks (such as the Combinatorial Benders Cuts approach introduced by Codato and Fischetti [10]), the approach presented in this work provides a tailored and efficient alternative for constructing exact algorithms for a class of vehicle routing problems involving synchronization constraints. Specifically, our methodology targets a practically motivated subclass centered on operation synchronization, which introduces complex temporal interdependencies between tasks. It not only enables the effective generation of iPECs, but also provides the necessary tools for their systematic strengthening. We now review the modeling foundations of this class of problems and how they integrate with broader temporal consistency requirements. 1.1. Overview of vehicle routing problems with synchronization Many practical routing applications feature synchronization constraints that introduce complex interdependencies between tasks. Our methodology targets a specific yet widely applicable class of such problems, centered on operation synchronization. We begin by reviewing its modeling foundations and how it integrates with broader temporal consistency requirements. Operation synchronization In many vehicle routing problems, operation synchronization constraints arise when two or more tasks, executed by different vehicles, must be temporally coordinated to carry out a joint operation. It represents a practical and increasingly relevant feature in modern logistics, especially in applications involving cooperative handling, transfers, or team-based services. The aim is to ensure that such tasks, although performed on separate routes, occur within a specified maximum allowed time window width (as maxTime Width throughout this paper) to maintain time consistency. Enforcing operation synchronization in vehicle routing problems, as discussed in Soares et al. [38], usually requires computing customer arrival times based on routing decisions. This is commonly achieved through precedence constraints that encode conditional relationships, often of the if-then type, between arrival time sequences. While these constraints can be linearized, as in the classic Miller-Tucker-Zemlin (MTZ) constraints for the Traveling Salesperson Problem, such formulations often rely on auxiliary variables and large constants. If not carefully calibrated, these constructs may weaken the LP relaxation, especially when routing variables take fractional values, thus shifting the computational burden to the branching phase. The study of synchronization in vehicle routing problems was first systematically formalized by Drexl [16], who proposed a classification including five types: task, operation, movement, load, and resource synchronization. More recently, Soares et al. [38] revisited this taxonomy and argued for a simplified distinction 2 between operation synchronization and movement synchronization. In this work, we focus exclusively on the former. In operation synchronization, the core requirement is that the arrival times of two related tasks must lie within a given temporal offset. These tasks are typically associated with different vehicles and routes, yet are logically part of the same operation. Movement synchronization, by contrast, involves the coordinated displacement of vehicles, which is beyond the scope of this study. A critical modeling issue in this context is the distinction between locations,tasks, and operations. Synchronization constraints often involve multiple tasks at the same location, which requires treating locations and tasks as separate modeling elements. Likewise, operations are understood as abstractions of task pairs that are interdependent. In our framework, any task carried out collaboratively by more than one vehicle is considered part of an operation composed of several tasks. We formalize an operation, here involving two tasks, as the set {i= (c, k), j = (c, k′)}, where iand j are tasks performed at the same location c, but assigned to different routes kand k′. In our setting, for all practical purposes, operation and location can be treated as equivalent concepts. Synchronization constraints The definition of synchronization constraints requires computing the arrival times of tasks, which in our context are represented by continuous variables si, where idenotes a task. It is useful to distinguish three types of constraints typically involved in this computation: lower-bound constraints,upper-bound constraints, and those that combine both. Consider two tasks i= (c, k)and j= (c′, k)assigned to the same route kbut associated with different locations cand c′. If task iprecedes task jin the route, a lower-bound constraint imposes a minimum on the arrival time of task j. This lower bound is defined relative to the arrival time of task iand is determined by a temporal offset tij , yielding the constraint: si+tij ≤sj. By contrast, an upper-bound constraint restricts the arrival time of task jrelative to task i, leading to: si+tij ≥sj. Both types of constraints are typically introduced by logical (nonlinear) formulations, e.g., “if task iprecedes task j, then . . . ”, which must be linearized in order to be used in a Mixed Integer Linear Programming (MILP) framework. The literature on vehicle routing has examined the role of idle time (intentional waiting periods between tasks) from multiple modeling perspectives, ranging from its explicit inclusion in the objective function to its treatment as an intrinsic element of the problem definition (see Soares et al. [38] for a comprehensive review). While incorporating idle time typically increases the total service duration (and may require additional constraints to limit route length) it also significantly enhances solution feasibility in the presence of strict arrival-time consistency requirements. 3 Vehicle routing problems with synchronization can generally be divided into two categories: those that allow idle time and those that do not. In problems that allow it, synchronization is typically enforced using only lower-bound constraints. In contrast, problems that prohibit idle time use a combination of lowerand upper-bound constraints to prevent any temporal slack between tasks. In this work, we focus exclusively on lower-bound precedence constraints, which are widely used in models that permit idle time to ensure temporal feasibility. Arrival-time consistency has emerged as a key factor in customer satisfaction, particularly in serviceoriented delivery systems where minimizing temporal variability across service days is essential. This requirement is especially critical for companies operating regular delivery schedules, such as food and beverage suppliers, grocery chains, and other retailers, where predictable service times can significantly enhance the customer experience (see Kovacs et al. [26]). To operationalize this concept, arrival-time consistency is typically enforced by limiting the variability in arrival times associated with synchronized tasks, that is, tasks belonging to the same operation but scheduled on different service days. Consider two such tasks, i= (c, k)and j= (c, k′), executed on days kand k′for the same customer c, with arrival times siand sj. The temporal offset between them is given by |sj−si|. While some studies consider both lower and upper bounds on this offset, in our setting we focus exclusively on imposing a maxTimeWidth, i.e., a maximum allowable width for arrival-time variation. This requirement, combined with the flexibility to insert idle time between consecutive tasks, gives rise to a class of vehicle routing problems with synchronization constraints that we refer to as vehicle routing problems with idling and arrival-time consistency (VRP-IATC). These problems capture the key temporal dimensions relevant to our study and provide the structural foundation for the proposed methodology. In the next section, we introduce representative problems within this class that satisfy the necessary conditions for its application. Idle times and arrival-time consistency in vehicle routing problems Several vehicle routing problems in the literature incorporate synchronization requirements stemming from idle times and arrival-time consistency. In this subsection, we review representative VRP-IATCs that exemplify these temporal features and provide a suitable ground for applying our approach. A prominent example, and the focus of this work, is the Consistent Traveling Salesperson Problem (ConTSP), chosen for its conceptual clarity and suitability to illustrate the core aspects of our technique. The ConTSP generalizes the Periodic Traveling Salesperson Problem (PTSP) [9], which assigns customers to service days within a fixed horizon. The ConTSP adds temporal consistency in the form of maxTimeWidth, requiring customers visited on multiple days to be served at similar times, thus promoting predictability. The goal is to minimize routing cost while maintaining consistent arrivals. Subramanyam and Gounaris [40,41] study two variants: ConTSP1, which forbids idling and minimizes travel cost subject to arrival-time deviations, but relies on precedence constraints incompatible with our approach; and ConTSP2, which permits idling and enforces a maxTimeWidth, enlarging the solution space. The latter is solved through a decomposition scheme and will serve as the case study in this paper. DíazRíos and Salazar-González [17] further compare four formulations for both CTSP variants, introducing a 4 multi-commodity flow model that enables efficient Benders decomposition. A related problem is the Selective Routing Problem with Synchronization [36], motivated by telescope scheduling, where multiple sensors must be synchronized. It can be modeled as a periodic Orienteering Problem with strict arrival-time consistency, and is solved via a branch-and-price-and-cut scheme. The Consistent Vehicle Routing Problem (ConVRP) [21] extends the ConTSP to capacitated fleets and driver consistency. Most formulations prohibit idling, though variants such as Goeke et al. [19] allow flexible departures and waiting at customers, solved through column generation. A further extension is the VRP with Synchronization at Delivery Locations [37], where idle times are permitted and customers define their own maxTimeWidths; an adaptive large neighborhood search is proposed. In food distribution, the Consistent VRP with Time Windows [30] enforces customer–driver consistency with hard time windows and irregular demand, solved with a heuristic minimizing the number of drivers. Collaboration settings are addressed by the Collaborative ConVRP with Workload Balance [32], where multiple providers cooperate under maxTimeWidths constraints, idle times, and workload balance, maximizing coalition profit. Finally, the Consistent Electric VRP with Backhauls and Charging [33] involves EV routing with deliveries and pickups under consistency and energy constraints. Idle times are strategically exploited for battery recharging, and the problem is solved using a hybrid ALNS with constraint and quadratic programming components. These variants illustrate the breadth of VRPs with synchronization, from conceptual models such as the ConTSP to application-driven problems in logistics, astronomy, and electromobility. 1.2. Exact methods for vehicle routing problems with synchronization This section reviews exact methods developed for solving vehicle routing problems with synchronization constraints, with particular emphasis on those that model idle times and arrival-time consistency, though not exclusively. All such methods, in one way or another, rely on the well-established branch-and-bound paradigm, which supports many exact algorithms for combinatorial optimization problems. As is well known, this technique systematically explores a search tree in which each node corresponds to a subproblem defined by additional constraints that reduce the feasible region. The objective is to decompose the original problem into more tractable instances while pruning portions of the search space that cannot contain optimal solutions. The key ingredients of any branch-and-bound method are a branching strategy, an effective bounding mechanism, and a node selection policy. We review both general-purpose approaches based on MILP formulations, typically solved by commercial solvers whose core is a branch-and-bound engine, and tailored branch-and-bound algorithms specifically designed for certain classes of synchronization problems. Direct MILP-based solution methods A first approach in the development of exact methods for vehicle routing problems with synchronization is to formulate the problem using a compact model that can be directly solved with a MILP solver such 5 as IBM ILOG CPLEX1, Gurobi2, or SCIP3. These formulations are typically based on extensions of the classical vehicle routing model, enhanced with additional variables and constraints that capture the temporal dependencies between synchronized tasks. A common modeling strategy involves using MTZ-type constraints to represent the sequence of visits and compute customer arrival times. Although these constraints are linear, they require the introduction of auxiliary variables to model service times and often rely on large constants to enforce conditional validity. If not properly tuned, these constants can significantly degrade the quality of the linear relaxation, thus impairing solver performance. While this modeling approach offers advantages in terms of simplicity and ease of implementation, it is limited in scalability. The linear relaxations of such formulations are often weak, leading to long solution times even for medium-sized instances. Nevertheless, compact MILP models remain a valuable tool for earlystage model testing and for solving small instances to global optimality. In fact, they are frequently used to provide a complete formal description of the problem, even when heuristic or alternative exact methods are later proposed for large-scale solution. This is the case in works such as Sarasola and Doerner [37], Kovacs et al. [28], Stavropoulou [39], Mancini et al. [32], and Nolz et al. [33], among others, who introduce compact formulations for routing problems with synchronization, but do not develop specific exact algorithms for their resolution. Similarly, the survey by Soares et al. [38] uses compact models as a framework to formally describe a broad family of synchronization based routing problems, without addressing their solution methodology. Some authors, however, have sought to strengthen these formulations to improve their computational effectiveness. For example, Díaz-Ríos and Salazar-González [17] propose a number of enhancements that significantly improve the quality of the LP relaxation and, consequently, the performance of MILP-based solution algorithms. From a computational perspective, commercial MILP solvers internally implement branch-and-bound algorithms, branching primarily on the values of the integer (typically arc binary) decision variables. During the bounding phase, these solvers automatically apply various techniques to tighten the linear relaxation, such as general-purpose cutting planes (e.g., Gomory cuts), and employ primal heuristics to quickly find feasible integer solutions. These internal processes are largely transparent to the user but contribute significantly to solver efficiency. In summary, compact MILP formulations provide a useful formal framework for representing vehicle routing problems with synchronization, and are well suited for solving small instances. However, due to their limited scalability stemming from weak LP relaxations, more advanced methods are often required to tackle larger and more complex problems. Branch-and-cut strategies Branch-and-cut algorithms extend the classical branch-and-bound approach by dynamically adding valid inequalities, commonly referred to as cuts, to strengthen the linear relaxation of the problem. These cuts 1https://www.ibm.com/products/ilog-cplex-optimization-studio 2https://www.gurobi.com/ 3https://www.scipopt.org/ 6 exclude fractional regions of the solution space without removing any feasible integer solutions, significantly improving the bounding power of the relaxed model. Two main categories of cuts are commonly employed. Feasibility cuts are introduced to eliminate fractional solutions that violate structural constraints of the problem, such as subtour elimination, capacity restrictions, or synchronization requirements. In contrast, reinforcement cuts aim to strengthen the formulation by tightening the linear relaxation, regardless of whether they are associated with explicit violations. The terminology used in the literature to refer to these cutting-plane strategies is not entirely consistent. Some authors reserve the term branch-and-cut for algorithms that incorporate only integrality reinforcement cuts, while approaches involving feasibility cuts, particularly those added dynamically, are referred to as branch-and-benders methods or lazy constraint generation schemes. In this work, we adopt a unified perspective and refer to all such methods under the umbrella of branchand-cut, regardless of the nature or purpose of the added cuts. A traditional approach for deriving problem-specific cuts is to analyze the underlying combinatorial structure and identify valid inequalities accordingly. For relatively simple problems, polyhedral analysis can lead to strong inequalities such as facet-defining cuts (see, e.g., Laporte et al. [29]). However, this becomes impractical as the complexity of the problem increases, particularly in the presence of synchronization constraints. Alternatively, valid cuts can be derived through decomposition techniques such as Benders decomposition (see, Benders [7]), where master and subproblem variables are separated and cuts are generated based on the dual solution of the subproblem. This approach is especially useful for problems with block structures or naturally decomposable formulations. Moreover, it is effective even when little is known about the subproblem’s structure. For instance, Riera-Ledesma and Salazar-González [35] propose a branch-and-cut algorithm where Benders cuts, derived from an opaque subproblem, are strengthened via left-hand-side rounding. Nevertheless, Benders cuts are often weaker than those obtained through direct combinatorial analysis, as they focus solely on subproblem feasibility without exploiting the structure of the original problem. Their form is typically not known in advance, which makes them difficult to strengthen, and their implementation can be computationally demanding, especially when Big-Mconstants are required to model feasibility conditions. As noted by Codato and Fischetti [10], “the classical Benders’ approach can be viewed as a tool to accelerate the solution of the LP relaxation, but not to improve its quality”. A complementary strategy for generating problem-specific cuts is the use of iPECs, which aim to eliminate invalid routes by identifying incompatible subpaths. This technique has been widely adopted and shown to be effective across various routing problems. For instance, Ascheuer et al. [4] apply it to the Asymmetric Traveling Salesperson Problem with Time Windows, Cordeau [11] to the Dial-a-Ride Problem, Luo et al. [31] to a manpower routing problem with synchronization, Alba Martínez et al. [1] to a vehicle routing problem with loading constraints, Subramanyam and Gounaris [40] to the ConTSP1, and Riera-Ledesma and Salazar-González [36] to the Selective Routing Problem with Synchronization. In many cases, the key challenge is to identify incompatible subpaths and generate cuts that remove the corresponding invalid solutions. These cuts are initially weak, as they typically eliminate only a single arc from the conflicting path. However, in some cases, they can be strengthened by analyzing the subpath’s 7 structure and adding extra arcs to the inequality’s left-hand side. This strengthening process is essential for improving relaxation quality and, consequently, algorithmic performance. A noteworthy example of this is the work of Ascheuer et al. [4]. Thus, the method consists of two main steps: identifying conflicting subpaths and generating reinforced cuts to eliminate them. While in some cases the identification of conflicting subpaths is straightforward, such as in the Asymmetric Traveling Salesman Problem with Time Windows [4] or the Dial-a-Ride Problem [11], in others, the process can be considerably more complex. For example, Alba Martínez et al. [1] and Iori and Riera-Ledesma [25] explore route-load compatibility in TSP variants with multiple stacks to identify infeasible loading and unloading sequences that induce subpath conflicts. Similarly, Subramanyam and Gounaris [40] detect pairs of routes that violate temporal consistency constraints under the current scheduling plan. In the context of the Selective Routing Problem with Synchronization, Riera-Ledesma and SalazarGonzález [36] address infeasibility arising from idle times by solving a dedicated feasibility subproblem. Since the VRP-IATC belongs to the class of problems in which identifying iPECs is particularly challenging, the present work builds upon and formalizes this latter strategy. Specifically, it introduces a systematic framework for detecting conflicting subpaths and generating strengthened iPECs tailored to the structural features of the ConTSP2. The novelty of this approach lies in the identification of conflicting subpaths by leveraging the underlying combinatorial structure revealed in the dual formulation of the temporal consistency verifier subproblem. This dual perspective allows us to isolate violations in a principled manner, leading to the derivation of effective infeasibility cuts. It is important to emphasize that, despite exploiting the dual of the temporal consistency verifier subproblem, the proposed method is not based on traditional Benders cut generation. Benders cuts are derived directly from the dual but typically do not capture any specific combinatorial structure. In contrast, the key contribution of our approach lies in using the dual formulation to uncover the combinatorial configuration responsible for infeasibility, an aspect that, to the best of our knowledge, has not been previously explored in this context. From a practical perspective, modern MILP solvers such as CPLEX and Gurobi allow for the implementation of customized branch-and-cut algorithms via user-defined callbacks. These interfaces enable the detection of infeasible fractional solutions during the search process and the dynamic addition of corresponding cuts. This allows users to focus on the design and validation of problem-specific cuts, while delegating branching, node selection, and other internal mechanisms to the solver. The methodology proposed in this work integrates seamlessly into this framework, enabling the rapid development of effective exact algorithms. In summary, branch-and-cut algorithms combine the structural power of valid inequalities with the flexibility of branching strategies, offering a robust framework to tackle vehicle routing problems with complex constraints such as synchronization and temporal consistency. Branch-and-price-and-cut strategies A more advanced class of exact algorithms for vehicle routing problems with synchronization relies on branch-and-price-and-cut techniques (see, e.g., Desrosiers and Lübbecke [13]), which integrate column 8 three service days, while customers at locations c8and c13 require service on only one day. This demand table indicates that the instance burma14_p3_f50_lH contains a total of 19 delivery tasks, in addition to 3 departure and 3 arrival tasks. Table 1: Demands for each location in set Cfor each service day, for the instance burma14_p3_f50_lH. C K0c2c3c4c5c6c7c8c9c10 c11 c12 c13 c14 1 2 3 Figure 1provides a graphical representation of a feasible solution for the instance burma14_p3_f50_lH. The figure includes two subfigures illustrating different aspects of the same schedule. The first subfigure presents the scheduling of delivery tasks across the three service days in K, while the second subfigure displays the same schedule, organized by customer locations in C′. The scheduling time of each task is represented by a vertical segment, with a label at the bottom indicating its identifier. The number on the left of each row in subfigure 1(a) denotes the corresponding service day. For example, the first vertical segment in row 1 of subfigure 1(a) indicates that delivery task (c9,1), associated with customer 9 on service day 1, is scheduled at time 160. Transitions between tasks are depicted as gray-shaded rectangles, with the transport time between two tasks specified inside. Additionally, subfigure 1(a) highlights a violet-shaded region preceding task (c12,1), following task (c7,1). This violet area represents a waiting time of 146 units before performing task (c12,1), after a travel duration of 163 units from task (c7,1). Although Figure 1does not explicitly display them, each service day kbegins with departure task (0, k) and concludes with arrival task (n+ 1, k). Subfigure 1(b) presents the same scheduling information as subfigure 1(a), but organized by customer. This visualization facilitates verification of constraint (1) compliance. A blue-shaded rectangle, centered at the average scheduling time of a customer’s tasks, represents the maxTimeWidth. This graphical representation confirms that all delivery tasks for a customer meet constraints (1), as they remain within the rectangle’s boundaries. Each row’s leftmost identifier corresponds to the associated customer. In subfigure 1(a), on the first day, the vehicle completes its deliveries in 2660 time units, though the route cost is 2514, considering a waiting time of 146 units. On the second service day, the vehicle consumes 2561 time units, and on the third service day, it consumes 2233 time units. None of the routes exceed the defined maximum service time. The total cost of this solution is 7308. Including waiting times helps achieve feasible solutions by accommodating the maxTimeWidth constraint. Consider the tasks (c12,1) and (c12,2) in subfigure 1(b). Both tasks are scheduled within the maxTimeWidth, allowing the vehicle to wait 146 time units on the first service day before performing task (c12,2). Without 15 342 (c13,3) 124 (c7,3) 168 (c6,3) 400 (c5,3) 491 (c4,3) 310 (c14,3) 398 3 160 (c9,2) 990 (c5,2) 418 (c12,2) 221 (c14,2) 211 (c3,2) 491 (c8,2) 70 2 160 (c9,1) 482 (c7,1) 163 146 (c12,1) (c6,1) 400 (c5,1) 491 (c4,1) 289 (c3,1) 510 1 0 500 1000 1500 2000 2500 (a) Timeline representation of delivery tasks across the 3 service days. Each row corresponds to a service day, as indicated on the left. 0 500 1000 1500 2000 2500 c14 (c14,2) (c14,3) c13 (c13,3) c12 (c12,1) (c12,2) c9 (c9,1) (c9,2) c8 (c8,2) c7 (c7,3) (c7,1) c6 (c6,3) (c6,1) c5 (c5,3) (c5,2) (c5,1) c4 (c4,3) (c4,1) c3 (c3,2) (c3,1) (b) Timeline representation of delivery tasks organized by customer location. Each row corresponds to a customer, as indicated on the left. Figure 1: Graphical depiction of a feasible solution for the ConTSP2 instance burma14_p3_f50_lH. The first subfigure illustrates tasks scheduled per service day, while the second organizes them by customer. this adjustment, the solution would be infeasible. However, incorporating waiting times comes at the cost of increasing the total service time. Additionally, imposing the maximum service time constraint ω0may result in a solution that is feasible in terms of max TimeWidth but infeasible due to exceeding the allowable service duration. 4. Formulations This section introduces and analyzes three MILP formulations for the ConTSP2, each designed to isolate and assess the impact of specific modeling choices on algorithmic performance. All three formulations build upon the combinatorial description of feasible solutions provided in Section 2.3, which serves as the structural foundation for the modeling approach. The first model, referred to as the Compact Formulation, encodes the full set of feasibility constraints using a polynomial number of inequalities. In particular, arrival-time consistency is enforced explicitly via time variables, which are used to model synchronization constraints. This model is fully expressible within 16 off-the-shelf MILP solvers and serves as a natural baseline for evaluating solution quality and computational performance. The second model, called the SEC-Enhanced Formulation, strengthens the compact model by incorporating classical connectivity constraints, which significantly improve the linear relaxation. As in the compact model, synchronization constraints are modeled through time variables, allowing us to evaluate the effect of standard structural reinforcement techniques without altering the overall modeling approach. Finally, the Path Elimination Formulation omits the arrival-time variables altogether and replaces the timing and consistency constraints with a family of infeasible path elimination constraints based solely on routing variables. This model leverages structural insights into the problem to generate strong feasibility cuts, and serves as the basis for the dual-driven methodology proposed in this work. All three formulations are tested and compared in Section 7.2 to evaluate their effectiveness and highlight the performance gains achieved by integrating iPEC-based reasoning. The best-performing formulation is then used in the computational comparison against the state-of-the-art algorithm of Subramanyam and Gounaris [41]. For clarity, we define the following standard notation. Given vertex subsets S, T ⊆V, we let: A(S) := {(i, j)∈A|i, j ∈S}, A(S:T) := {(i, j)∈A|i∈S, j ∈T}, δ−(S) := {(i, j)∈A|i∈V\S, j ∈S}, δ+(S) := {(i, j)∈A|i∈S, j ∈V\S}. If S={i}, we use the shorthand δ−(i)and δ+(i). We define two families of variables. The first is used in all three models; the second appears only in the first two models: •For each arc (i, j)∈A, the binary variable xij takes the value 1 if the vehicle performs task jimmediately after task i. •For each task j∈V, the continuous variable sjrepresents the arrival time at task j. For any subset of arcs F⊆A, we write x(F)to denote P(i,j)∈Fxij. The first two models, both based on arrival-time variables, include three main components: (i) constraints defining a Hamiltonian path for each service day, (ii) precedence constraints encoding arrival times, and (iii) arrival-time consistency constraints. The path elimination model also incorporates Hamiltonian path constraints but replaces the latter two components with iPECs. Each of the three families of constraints is described in what follows. 17 Hamiltonian paths Given the structure of the graph G= (V, A), defined above, the following constraints define |K|Hamiltonian paths, one for each service day, based on the xvariables: x(δ+(j)) = 1 j∈V\0−(2) x(δ−(j)) = 1 j∈V\0+(3) x(A(S)) ≤ |S| − 1S⊂V, |S| ≥ 2(4) xij ∈ {0,1}(i, j)∈A, (5) where 0−:= {(n+ 1, k)|k∈K}and 0+:= {(0, k)|k∈K}are the sets of arrival and departure tasks at the depot, respectively. Constraints (2)–(3) ensure that each task is visited exactly once. Constraint (4) prevents cycles in the service path. Constraint (5) enforces the binary nature of the variables x. Precedence constraints These constraints determine the arrival times sjbased on a selected path defined by variables x. They can be expressed using logical terms as follows: xij = 1 ⇒si+tij ≤sj(i, j)∈A. (6) These logical constraints ensure that if task jfollows i, then jstarts no earlier than the completion of iplus the travel time between them, for each arc (i, j)∈A. They naturally allow for waiting time between tasks. To linearize these constraints, we apply the standard big-Mreformulation: si+tij −M(1 −xij)≤sj(i, j)∈A(7) Here, Mis a large constant that effectively deactivates the constraint when xij = 0. These are known as MTZ constraints (see, e.g., Gouveia and Pires [20]). It should be noted that (7) also prevents subtours, making (4) redundant in integer solutions. A valid choice for Mis ω0, since no sjcan exceed this value due to the maximum service time constraints. Arrival time consistency constraints This final group of constraints enforces the arrival time consistency requirements. Specifically, for any pair of deliveries (i, j)to the same customer c∈C: sj−si≤ωc(i, j)∈B[c], c ∈C. (8) These inequalities ensure that delivery tasks for the same customer occur within the maxTimeWidth. They also implicitly enforce the upper bound ω0on the total duration of each service day. 4.1. Compact formulation The compact formulation serves as a baseline model. It captures all relevant constraints using a polynomial number of inequalities and can be directly handled by any commercial MILP solver. This model 18 combines standard routing constraints with a linearized version of the precedence logic, using big-Mcoefficients to enforce the correct timing of tasks. In addition, arrival time consistency is enforced by bounding the temporal gap between visits to the same customer. Despite its structural simplicity and ease of implementation, this model suffers from weak linear relaxations due to the presence of big-Mterms. These constraints become particularly ineffective when xvariables take fractional values, resulting in poor lower bounds and increased reliance on the branching mechanism. Nevertheless, the compact formulation is useful for analyzing the overall structure of the problem and provides a reference point for evaluating the effectiveness of more advanced models. min X (i,j)∈A dijxij (9) subject to: x(δ+(j)) = 1 j∈V\0−(2) x(δ−(j)) = 1 j∈V\0+(3) xij ∈ {0,1}(i, j)∈A(5) si+tij −M(1 −xij)≤sj(i, j)∈A(7) sj−si≤ωc(i, j)∈B[c], c ∈C. (8) The objective function (9) minimizes the total distance traveled. Time variables sjdo not require explicit upper bounds thanks to (8), and initializing s(0,k)= 0 for each k∈Kis also unnecessary. 4.2. SEC-Enhanced formulation The SEC-Enhanced formulation builds directly upon the compact model by incorporating a classical family of connectivity constraints, which are added dynamically during the branch-and-bound process. These constraints explicitly eliminate infeasible cycles in the routing structure and provide a significant strengthening of the LP relaxation, especially in early nodes of the search tree. This enhancement allows for a more accurate assessment of the solution space and reduces the number of nodes required to prove optimality. Although the use of connectivity constraints introduces an exponential number of constraints, their separation is well understood and efficient in practice. This formulation provides an intermediate step between the purely compact representation and the more sophisticated path elimination approach described next. min X (i,j)∈A dijxij (9) 19 subject to: x(δ+(j)) = 1 j∈V\0−(2) x(δ−(j)) = 1 j∈V\0+(3) x(A(S)) ≤ |S| − 1S⊂V, |S| ≥ 2(4) xij ∈ {0,1}(i, j)∈A(5) si+tij −M(1 −xij)≤sj(i, j)∈A(7) sj−si≤ωc(i, j)∈B[c], c ∈C. (8) Constraints (4) are dynamically added on demand during the solution process to eliminate infeasible cycles and strengthen the LP relaxation. 4.3. Path elimination formulation The third formulation introduces a qualitative shift in modeling: instead of relying on precedence and timing variables, it enforces feasibility through a family of infeasibility path elimination constraints (iPECs) expressed solely in terms of the routing variables x. These constraints are derived from structural incompatibilities between subpaths that cannot coexist in a feasible solution due to synchronization requirements. This approach eliminates the need for big-Mcoefficients and leads to tighter LP relaxations by directly targeting infeasible routing patterns. It forms the core of the dual-driven separation strategy proposed in this work, where infeasibility is identified through subproblem analysis and used to generate strong, problemspecific cuts. The resulting model preserves the routing structure of the SEC-Enhanced formulation, but replaces timing and consistency constraints with iPECs dynamically added during the solution process. As shown in the computational study, this formulation achieves substantial improvements in both solution times and bounds, especially for larger or more tightly constrained instances. min X (i,j)∈A dijxij (9) subject to: x(δ+(j)) = 1 j∈V\0−(2) x(δ−(j)) = 1 j∈V\0+(3) x(A(S)) ≤ |S| − 1S⊂V, |S| ≥ 2(4) X P∈T x(P)≤X P∈T |P| − 1T ∈ T(10) xij ∈ {0,1}(i, j)∈A. (5) Constraints (10) remove those combinations of subpaths that lead to infeasibility because of maxTime Width violation. The set Trepresents the collection of all sets of mutually incompatible paths, where each T ∈ Tdenotes a specific subset of paths that are incompatible with each other. 20 At this point, several questions arise: how can these new inequalities (10) be generated, and how can we ensure they are strong enough to serve as an alternative to the compact model? We address these questions in the next section. 5. Identifying iPECs To identify iPECs that render any infeasible solution, we outline a methodology in this section. First, we establish a method to determine whether a solution described in terms of variables xis feasible for the ConTSP2. Then, we propose a method to generate inequalities that, based on this feasibility check, help avoid infeasible solutions. Finally, we introduce a third step that reinforces the obtained constraints by tightening them, further improving the formulation’s strength. 5.1. Checking feasibility Consider a solution ˜xobtained from the model (2)–(5) and (9), where the integrality constraints (5) have been relaxed. Additionally, branching constraints imposed on the xvariables may be present, as well as some constraints from the family (10) that may have been introduced earlier. The solution ˜xis feasible for the ConTSP2 if there exist values of ssatisfying the following constraints: −sj+ ˜xijsi≤ −˜xijtij (i, j)∈A(αij)(11) sj−si≤ωc(i, j)∈B[c], c ∈C. (βij)(8) Where constraints (11) particularize constraints (6) for the solution ˜x. To determine whether there exists a vector ssatisfying constraints (11) and (8) for a given solution ˜x, we apply duality theory by associating dual variables αand βwith the synchronization constraints (11) and the maxTimeWidth constraints (8), respectively. The feasibility of this model is then directly related to the unboundedness of the following dual linear program: min −X (i,j)∈A ˜xijtij αij +X (i,j)∈B[c],c∈C ωcβij (12) s.t. X (j,i)∈A[k] ˜xjiαji −X (i,j)∈A[k] αij = X (j,i)∈B[c] βji −X (i,j)∈B[c] βij j= (c, k)∈V(13) αij ≥0 (i, j)∈A(14) βij ≥0 (i, j)∈B. (15) The set of constraints (13) enforces flow conservation at each vertex j= (c, k)in the network G= (V, A∪B). Within each temporal layer A[k], for every k∈K, flow is driven by the αij variables over arcs (i, j)∈A[k] that are active in ˜x, while flow across arcs in Bis determined by the βij variables. 21 (0, k) (c1, k) (c2, k) (c3, k) (c4, k) (n+ 1, k) Figure 2: Example of a fractional solution containing two paths with value 0.5 from (0, k)to (n+ 1, k). The first path follows (0, k)→(c1, k)→(c2, k)→(c4, k)→(n+ 1, k), while the second follows (0, k)→(c1, k)→(c3, k)→(c4, k)→(n+ 1, k). Solid arcs indicate edges with ˜x= 1, while dashed arcs correspond to edges with ˜x= 0.5. It is important to note that the solution ˜xsatisfies constraints (2)–(4), with 0≤xij ≤1, for all (i, j)∈A. Therefore, for each period k∈K, the subgraph induced by arcs in A[k]contains only paths, either integer or fractional, connecting (0, k)to (n+ 1, k). As a consequence, there are no directed cycles entirely contained within any individual layer A[k]under ˜x. If a layer A[k]were to contain arcs forming a (fractional) cycle O⊆A[k], for instance, due to the presence of two or more overlapping fractional paths in ˜x(e.g., O= ((c2, k),(c3, k)),((c3, k),(c2, k)), as illustrated in Figure 2), with ˜xij <1, for (i, j)∈O, then the left-hand side of constraint (13) would violate the flow conservation condition for each vertex involved in the cycle. That is, for any task jwithin the cycle, the corresponding constraint would fail to satisfy X (j,i)∈O ˜xjiαji =X (i,j)∈O αij, thereby disrupting the feasibility of a cyclical flow in the dual solution. Constraints (14) and (15) impose non-negativity on the dual variables αand β, respectively. Note that this dual model is always feasible, as it admits the trivial solution (˜α, ˜ β) = (0,0), although our interest lies in its potential unboundedness as a certificate of primal infeasibility. In conjunction with the objective function (12), this model determines whether the solution ˜xis feasible for the ConTSP2. The objective function consists of two components: the first one is proportional to the values of ˜xon the selected arcs and to the transition service times between tasks in A, while the second one is proportional to the maxTimeWidth associated with the arcs in Bthat participate in the cycle. Practically, the objective function compares the total transition service time across the cycle’s arcs in A with the aggregate width of the time windows associated with the arcs in Bthat complete the cycle. If the transition time exceeds the cumulative time-window allowance, the objective becomes unbounded, signaling infeasibility. 22 5.2. Generating iPECs To derive valid infeasible path inequalities from an infeasible solution ˜x, we establish upper bounds on the variables αand β. To this end, we replace constraints (14) and (15) with the following bounds: 0≤αij ≤¯α(i, j)∈A(16) 0≤βij ≤¯ β(i, j)∈B, (17) where ¯αand ¯ βcan take the same value to prevent unbounded solutions. This reformulation ensures that the objective function (12) remains bounded. We now show that any infeasible solution ˜xto the ConTSP2 gives rise to a sequence of mutually incompatible subpaths. This structural insight serves as the foundation for deriving valid inequalities for the ConTSP2, as well as for developing techniques to strengthen them within an exact solution framework. Proposition 1. Let ˜xbe any infeasible solution to the ConTSP2. Then, there exists a finite cyclic structure in the network G= (V, A ∪B), composed of a sequence of subpaths P1, . . . , Pp, with p≥1, each associated with a service day and connected by consistency arcs. This structure, referred to as a consistency cycle, is responsible for the violation of the temporal constraints of the problem. Proof. Assume that ˜xis an infeasible solution to the ConTSP2. Consider the dual model D(˜x), defined by constraints (13), (16), (17), and objective function (12). If ˜xis infeasible, then the dual problem attains a strictly negative optimal value, indicating that the total service time across transitions exceeds the permissible cumulative time defined by the consistency windows. Therefore, there exists an optimal dual solution (˜α, ˜ β)with non-zero objective value. Let ˜ A:= {(i, j)∈A|˜αij >0}and ˜ B:= {(i, j)∈B|˜ βij >0}denote the support of this solution, and define the corresponding support graph ˜ G= (V, ˜ A∪˜ B). Since the flow conservation constraints (13) enforce zero net flow at every vertex, each connected component of ˜ Gmust form a closed structure within G. Furthermore, as shown in Section 5.1, intra-layer cycles within each A[k]are excluded by construction. Consequently, any closed component in ˜ Gmust necessarily include at least one arc from B, linking subpaths across different service periods. Two distinct types of consistency cycles may arise: •A single intra-day path within some layer A[k], closed by an arc ((n+1, k),(0, k)) ∈B, violating the global time window ω0. •A cross-period sequence of subpaths connected by arcs in B[c]for some c∈C, resulting in a violation of the inter-period consistency constraint ωc. Moreover, if the dual solution is obtained using a vertex-based algorithm such as the Simplex method, it corresponds to an extreme point of the feasible polyhedron. Hence, the connected components of the support graph cannot be expressed as convex combinations of simpler structures; they form minimal, elementary cycles that expose the infeasibility of ˜x. In summary, each component of the support graph ˜ Gdefines a consistency cycle, composed of subpaths joined across days by consistency arcs, whose aggregate service time violates the time consistency conditions of the ConTSP2. This completes the proof. 23 342 (c13,3) 1307 (c7,3) 182 (c6,3) 1413 (c5,3) 491 (c4,3) 310 (c14,3) 398 3 70 (c8,2) 2941 (c3,2) 211 (c14,2) 221 (c12,2) 418 (c5,2) 990 (c9,2) 160 2 1784 (c9,1) 482 (c7,1) 163 (c12,1) 400 (c5,1) 491 (c4,1) 289 (c3,1) 510 1 Figure 3: Example of an infeasible partial solution for the instance burma14_p3_f50_lH, illustrating the values of ˜xand (˜α, ˜ β) across service days. Transitions included in ˜xare shown as shaded rectangles; arcs corresponding to ˜αare highlighted in blue, and those corresponding to ˜ βin red. (c9,1) (c7,1) (c7,3) (c6,3) (c5,3) (c4,3) (c14,3) (c14,2) (c12,2) (c5,2) (c9,2) 482 182 1413 491 310 221 418 990 617 617 617 Figure 4: Graph representation of the values (˜α, ˜ β)shown in Figure 3, embedded in graph G. Figure 3provides a graphical representation of the structure of the cycles described by D(˜x)through a concrete example. The figure depicts a solution in terms of the αand βvariables derived from a single infeasible solution ˜x. The input data corresponds to the instance burma14_p3_f50_lH, whose details have been presented in Section 3. The figure displays the values of ˜x, that is, the transitions between task pairs (i, j)included in solution ˜x, as gray-shaded rectangles. In this particular case, they all correspond to ˜xij = 1, since the infeasible solution ˜xobtained in this iteration is integer. Blue-shaded rectangles similarly represent transitions included in ˜x, highlighting the involved tasks. The transportation time between tasks is shown inside each transition rectangle. The figure also illustrates the solution values ˜αand ˜ βfrom the variables αand βfor D(˜x). Blue-shaded rectangles denote transitions with non-zero values of ˜αfor task pairs i= (c, k)and j= (c′, k)where ˜xij >0, for some c, c′∈Cand k∈K. Red arcs represent non-zero values of ˜ βand connect segments from different service days, specifically, tasks i= (c, k)and j= (c, k′), for some c∈Cand k=k′∈K. Figure 4shows the consistency cycle associated with the same ˜αand ˜ βvalues from Figure 3, now embedded in the graph G. In this example, the total transportation time across the arcs in Abelonging to the three subpaths equals 4507 time units. Meanwhile, the total time-window allowance across the three arcs in B, which link those subpaths, adds up to 1851 time units. This results in a violation of 2656 time units. Building on this result, we observe that, for any infeasible solution ˜xto the ConTSP2, the set Tof subpaths extracted from the corresponding dual solution (˜α, ˜ β)defines a valid inequality of the form: X P∈T x(P)≤X P∈T |P| − 1.(18) This inequality eliminates ˜xby explicitly excluding the combination of subpaths that jointly violate the problem’s consistency constraints. 24 Algorithm 2 Separation process for infeasible path elimination constraints. 1: function IdentifyCycle(˜x)▷Returns (T) 2: (˜α, ˜ β)←Solve(D(˜x)) 3: ˜ A← {(i, j)∈A|˜αij >0} 4: ˜ B← {(i, j)∈B|˜ βij >0} 5: T ← GetSubpaths(˜ A, ˜ B) 6: return T 7: end function 8: 9: function IdentifyIpec(T)▷Returns (C) 10: find a partition T1,T2of Tsuch that PP∈T1t(P) + PP∈T2t∗(P)> ω(T) 11: C ← PP∈T1x(∆+(P)) + PP∈T2x(∆ϕ(P)) ≤PP∈T |P| − 1 12: return C 13: end function 14: 15: function SeparateIpec(R)▷Returns (Rip) 16: ˜x←Solve(minx∈R dtx) 17: T ← IdentifyCycle(˜x) 18: C ← IdentifyIpec(T) 19: Rip ← R ∪ C 20: return Rip 21: end function 31 violations are detected. At each iteration, the set of constraints R′is updated with any newly identified inequalities. Once the iterative process concludes, the solution ˜xobtained in the final execution of the SeparateIpec routine is used to compute the dual bound z∗for the current search node, using the expression z∗=dt˜x. Algorithm 3 Global separation process. 1: function Separation(R)▷Returns (R′) 2: R′← R 3: C ← ∅ 4: repeat 5: Rco ←SeparateConn(R′) 6: Rip ←SeparateIpec(Rco) 7: C ← R′\ Rip 8: R′← Rip 9: until C=∅ 10: return R′ 11: end function 6.2. Initial heuristic A good initial upper bound can accelerate convergence by enabling earlier pruning in the branch-andbound process. We implemented two complementary strategies to construct such bounds. The first follows a template-based approach inspired by Subramanyam and Gounaris [41]. A TSP solution over all customers in Cis computed and used to define a global precedence structure. Each daily route adheres to this structure, ensuring that the arrival-time consistency constraints for customer visits are satisfied (though not necessarily the overall route duration constraint). The second strategy solves a series of TSPs, one per service day, each incorporating a precedence pattern updated from the previous iteration. This sequential refinement ensures consistency across days. Both approaches guarantee satisfaction of the consistency constraints for repeated customer visits but may violate time constraints at the depot. For the first strategy, we use the nearest neighbor and cheapest insertion heuristics, followed by 2-opt local search. When no further improvement is observed, diversification is triggered and the best solution is retained. This procedure is repeated for a fixed number of iterations. The second strategy requires solving multiple instances of the Precedence Constrained Traveling Salesperson Problem (PC-TSP), for which we use a custom MILP model. As in the main algorithm, path elimination constraints enforce precedence compliance. A simple constructive heuristic provides initial solutions to help guide the solver. Each instance is limited to 30 seconds, resulting in a worst-case time of 150 seconds for five service days. Results from both strategies are reported in Tables 3and 4(see Section 7). For each instance, we report the total solution time, root node gap, and time spent solving heuristic subproblems. These metrics allow us to assess the relative impact of the heuristic phase on overall computational performance. 32 7. Computational results The computational study presented in this section pursues a dual objective. First, it aims to validate the design of the proposed dual-driven path elimination methodology by comparing the three formulations introduced in Section 4, evaluating their effectiveness and quantifying the performance improvements resulting from the integration of iPEC-based reasoning. Second, it seeks to provide experimental evidence that the dual-driven approach, when embedded within a branch-and-cut framework for solving the ConTSP2, is not only effective but also competitive with, and in some respects superior to, the best-performing exact algorithm currently available in the literature. To this end, we begin by comparing the performance of the three aforementioned formulations: a baseline compact model that uses arrival-time variables, a strengthened variant that augments this model with connectivity constraints, and a third formulation based on infeasible path elimination constraints, which entirely avoids the use of arrival-time variables. This comparison reveals the impact of structural reinforcement and illustrates the advantages of enforcing temporal consistency through routing-based cuts. We then contrast the performance of the most effective of these models, the Path Elimination formulation, with that of the decomposition-based method proposed by Subramanyam and Gounaris [41], which remains the current state of the art for solving the ConTSP2. All implementations were developed in C++ and compiled with gcc version 13.3.0 on a Linux Ubuntu 24.04 LTS system using the -O2 optimization flag. Experiments were run in single-threaded mode on an Intel Core i5-7500 processor with 20 GB of RAM. CPLEX 22.1, accessed through its callable library, was employed to solve the branch-and-cut formulations, with a time limit of 7,200 seconds per instance, in line with the computational settings adopted in Subramanyam and Gounaris [41]. The implementation relies primarily on embedding direct replicas of Algorithm 3within CPLEX’s callback mechanisms. In particular, the procedure is invoked through the LazyConstraintCallbackI and UserCutCallbackI classes, used respectively to separate feasibility cuts and to generate valid inequalities that strengthen the LP relaxation. In addition, the same routine is called within the IncumbentCallbackI to validate candidate integer solutions; whenever a violated cut is identified, the incumbent solution is discarded. To provide an immediate perspective on the scope of the experiments and the effectiveness of the proposed approach, we open this section with a concise summary of the main computational outcomes. Table 2reports the aggregated performance of the main algorithms considered in this study. For each method, it presents the number of benchmark instances solved to proven optimality and the corresponding average runtime, together with the number of unsolved instances and their mean residual optimality gap. The results reveal a clear progression from the SEC-enhanced formulation to the dual-driven path elimination model, which solves the largest number of instances within the shortest average time. For reference, the aggregate results reported by Subramanyam and Gounaris [41] are also included. The remainder of this section is organized as follows. Section 7.1 describes the benchmark instance sets and summarizes their main characteristics. Section 7.2 presents a preliminary comparison among the three proposed formulations, highlighting the advantages of the path elimination model. Section 7.3 contrasts the results of this model with those reported by Subramanyam and Gounaris [41], the best-performing exact 33 Table 2: Aggregated computational performance of the main proposed algorithms. For each method, the table reports the number of benchmark instances (out of 756) solved to proven optimality (#) and their average runtime (seconds), as well as the number of unsolved instances (out of 756) and their mean residual gap (%). Results reported by Subramanyam and Gounaris [41] are included for reference. Proven optimal Residual gap Method # t(sec) # % gap Proposed algorithmsa SEC-Enhanced 401 469.3 355 1.4 Path elimination 536 381.8 220 1.6 Subramanyam and Gounaris [41]b524 492.1 232 2.8 aIntel Core i5-7500 CPU @ 3.40GHz. bIntel Xeon @ 2.80GHz family. method for the ConTSP2. This part includes both a methodological comparison between the two approaches and a concise performance assessment, which reinforces the contribution of our dual-driven strategy in terms of simplicity and ease of implementation. Finally, Section 7.4 provides details on the availability of the complete datasets and results generated by the proposed algorithms. 7.1. Test instances To evaluate our algorithm, we used the two ConTSP benchmark sets proposed in Subramanyam and Gounaris [40,41]. The first set is based on the 23 TSPLIB instances (Reinelt [34]) with up to 51 nodes, which were extended into threeand five-period versions of the ConTSP. Each instance is defined by two parameters: (a) the probability f∈0.5,0.7,0.9that a customer requires service in a given period, and (b) the value ωfor customers in C′, set at 10%, 15%, or 20% of the maximum (across periods) travel time of the corresponding optimal single-period traveling salesman tours. Each original TSPLIB instance was thus transformed into nine three-period and nine five-period ConTSP instances, yielding a total of 414 instances. In these instances, travel costs and travel times coincide, but they are not necessarily symmetric or hold with the triangle inequality. The second set considers 342 instances following the same procedure described in Kovacs et al. [27] and Subramanyam and Gounaris [40], applied to all TSPLIB instances with up to 101 nodes. These instances were also extended into threeand five-period versions of the ConTSP. Subramanyam and Gounaris have defined a value for the duration limit ω0. They have set ω0to 1.1 times the maximum (across periods) travel time of the corresponding optimal single-period TSP tours. With the experimental setup defined and the benchmark instances introduced, we now proceed to evaluate the performance of the proposed solution strategy. We begin by comparing the three algorithmic implementations described above, focusing on their behavior across different instance families. 7.2. Evaluation of the proposed implementations The three formulations presented in Section 4reflect different modeling philosophies: the Compact model emphasizes simplicity and ease of implementation; the SEC-Enhanced model incorporates classical 34 structural reinforcements; and the Path Elimination model leverages problem-specific infeasibility cuts based on incompatible subpaths. Unlike the first two models, which explicitly encode synchronization constraints via arrival-time variables, the Path Elimination formulation avoids these variables altogether, embedding feasibility through structurally derived cuts. To assess their relative performance, we conducted a preliminary computational study on a representative subset of instances. All three models, Compact, SEC-Enhanced, and Path Elimination, were implemented and solved using an identical MILP solver configuration. The results exhibit a clear performance gradient. The Compact formulation performed poorly even on smalland medium-sized instances, primarily due to the weak LP relaxation induced by big-Mconstraints. Adding connectivity constraints in the SEC-Enhanced model significantly improved both bounds and solution times. Finally, the Path Elimination formulation consistently outperformed both alternatives, producing tighter bounds and faster convergence, thanks to the integration of strong, problem-specific infeasibility cuts. Table 3summarizes this initial evaluation. It reports the performance of the three implementations across 42 instance families. Due to the poor performance of the Compact formulation, we report its results only for the first 23 families (the first 414 instances). The main goal of this study is to assess the impact of reinforcing the Compact model with connectivity constraints and replacing synchronization constraints (expressed via arrival-time variables) with subpath elimination cuts. While these enhancements compromise the structural compactness of the model, they substantially improve the number of instances solved to optimality, the root node relaxation quality, and the overall solution time. Table 3is organized into three main blocks, each corresponding to one of the evaluated implementations. For each block, we report the following four performance indicators: •r(%): average root node gap, calculated as ub−lb ub ×100, where lb and ub are the lower and upper bounds after root node processing; •#: number of instances solved to optimality (out of 18) within each family; •t(sec): average solution time (in seconds) for solved instances; •g(%): average optimality gap for unsolved instances at termination. The first two columns of the table provide additional context: •Instance: name of each instance family; •h(sec): average running time (in seconds) of the initial heuristic for that family. The impact of connectivity constraints on the Compact model is particularly notable, even though they are not required for feasibility. These constraints act as strong reinforcements of the LP relaxation. While the Compact formulation alone solves only 184 of the first 414 instances, its connectivity-enhanced counterpart solves 325, and does so with an average time of 566 seconds, compared to the 988 seconds required by the 35 original formulation to solve half as many. The average root gap also improves significantly, from 10.7% in the Compact model to just 1.4% in the SEC-Enhanced version, confirming the value of connectivity as a relaxation strengthening mechanism. The Path Elimination implementation clearly outperforms both alternatives. It solves all 18 instances in 16 of the 42 families and achieves global optimality for 536 of the 756 benchmark instances. Notably, 15 of these 16 families involve relatively small instances with 13-41 customers. However, the approach fails to solve any instances in two large families with 69 and 75 customers, respectively. Despite this, the Path Elimination model exhibits superior runtime performance, even when compared to the SEC-Enhanced formulation, which appears faster only in cases where it solves fewer instances. The performance gap between the SEC-Enhanced and Path Elimination models is evident as early as the first 15 instance families: the latter solves all of them, while the former succeeds with only nine. In many cases, solution times for the SEC-Enhanced model exceed those of the Path Elimination model by an order of magnitude. In the few rows where the SEC-Enhanced version appears faster, it corresponds to instances it failed to solve. Overall, across the 42 instance families, the SEC-Enhanced model solves 131 fewer instances and, on average, requires over 20% more time than the Path Elimination model on those it does solve. A more detailed analysis focused exclusively on the SEC-Enhanced and Path Elimination models, disaggregated by instance parameters, is presented in Table 4. This table also aims to experimentally assess the impact of replacing arrival-time variables and their associated consistency constraints with infeasible path elimination constraints. The results are grouped according to key characteristics such as the number of customers (n), the number of delivery tasks (m), the number of service periods (|K|), and two instancegeneration parameters: the fraction parameter fand the width of ω(see Subsection 7.1). The columns of the table are defined as follows: •Parameters: description of each parameter group, including the number of customers n, delivery tasks m, service periods |K|, the fraction parameter f, and the width of ω; •h (sec): average runtime (in seconds) of the initial heuristic for instances in the group; •#ins: number of instances in the parameter group. For each of the two algorithms implementing the SEC-Enhanced and Path Elimination models, the following performance indicators are reported: •#sol: number of instances solved to optimality within the group; •r (%): average root node gap, computed as ub−lb ub ×100; •t (sec): average runtime (in seconds) for solved instances; •#nod: average number of branch-and-bound nodes explored (for solved instances); •#conn: average number of connectivity cuts added during the solution process (for solved instances); •#ipec:(Path Elimination model only) average number of infeasible path elimination constraints added (for solved instances); 36 Table 3: Performance comparison of three exact implementations for solving the ConTSP2, grouped by instance family. CompactaSEC-Enhanced Path elimination Instance h(sec) r(%) # t(sec) g(%) r(%) # t(sec) g(%) r(%) # t(sec) g(%) burma14 0.5 11.1 16 201.5 2.0 0.6 18 0.2 - 0.6 18 0.1 - ulysses16 0.8 11.3 12 744.0 3.9 0.9 18 0.6 - 0.8 18 0.2 - br17 0.2 45.2 15 90.8 2.3 0.0 18 0.0 - 0.0 18 0.0 - gr17 1.4 8.1 16 755.0 0.6 0.6 18 3.2 - 0.6 18 0.4 - gr21 0.6 0.4 18 0.5 - 0.1 18 0.1 - 0.1 18 0.0 - ulysses22 0.7 13.3 3 2,289.6 5.5 0.1 18 0.4 - 0.1 18 0.1 - gr24 1.2 2.7 18 191.5 - 0.4 18 194.7 - 0.5 18 12.0 - fri26 1.3 5.9 12 711.3 1.6 1.1 17 680.5 0.2 1.0 18 14.9 - bayg29 2.2 3.5 12 289.4 0.8 0.4 18 4.9 - 0.5 18 0.5 - bays29 1.8 3.4 15 491.4 0.8 0.3 18 16.6 - 0.3 18 1.3 - ftv33 2.2 6.1 13 1,060.4 2.1 3.3 14 182.2 1.1 2.5 18 330.8 - ftv35 2.6 6.1 8 1,826.4 1.9 3.6 14 985.5 1.6 3.2 18 429.8 - ftv38 4.3 4.6 9 974.3 2.1 2.8 11 663.0 1.5 2.2 18 1,427.1 - dantzig42 12.9 11.5 0 - 6.9 1.2 16 570.2 0.2 1.2 18 34.0 - swiss42 7.6 5.9 6 1,493.7 0.9 1.5 15 847.0 0.7 1.4 18 158.9 - p43 39.6 64.7 0 - 9.0 0.1 6 60.1 0.1 0.1 9 653.2 0.0 ftv44 4.7 5.9 0 - 2.3 3.7 7 2,641.9 1.7 3.2 9 108.5 1.3 att48 12.1 9.3 0 - 5.1 1.3 10 957.5 0.9 1.4 16 151.8 0.3 ftv47 7.5 4.3 3 397.8 2.3 3.0 9 513.0 2.6 2.6 12 664.5 1.4 gr48 18.0 5.5 5 2,771.7 2.3 1.1 12 644.4 0.6 1.1 14 77.9 0.5 hk48 15.8 5.5 0 - 1.8 1.4 12 1,083.2 0.9 1.6 14 87.7 0.5 ry48p 17.9 8.0 0 - 3.7 1.7 14 1,082.0 0.8 1.3 17 39.9 0.7 eil51 15.2 4.6 3 2,498.1 1.8 2.5 6 1,884.1 1.2 2.2 11 1,193.5 1.2 Aggr. Inst. 1 - 414 7.4 10.7 184 987.5 2.8 1.4 325 565.9 1.0 1.2 372 234.2 0.8 berlin52 40.5 0.5 14 347.3 0.3 0.5 18 23.4 - ft53 49.6 2.7 3 814.5 1.1 2.2 13 232.1 1.9 ftv55 16.4 3.2 8 116.0 1.9 2.6 12 217.4 3.2 brazil58 33.3 0.1 17 28.3 0.2 0.2 17 1.8 0.0 ftv64 54.7 3.8 3 64.4 2.4 3.1 5 2,035.9 1.4 ft70 211.0 4.1 0 - 3.0 4.1 0 - 3.8 st70 36.9 1.2 7 105.5 0.6 1.2 12 272.2 0.7 ftv70 51.4 4.1 0 - 2.4 3.4 5 2,440.0 2.3 eil76 75.5 2.9 0 - 2.2 2.7 3 1,541.7 2.2 pr76 198.4 3.0 0 - 1.9 2.9 0 - 1.9 gr96 162.5 2.0 0 - 0.9 1.7 9 1,333.4 1.0 rat99 178.2 2.2 1 6,207.0 1.4 2.3 7 2,217.6 1.6 kroA100 105.3 1.6 3 61.0 0.8 1.0 10 1,194.5 0.7 kroB100 148.1 1.8 3 1,312.2 1.3 1.6 5 1,523.2 1.1 kroC100 142.1 1.2 5 1,195.9 0.5 0.9 14 1,063.9 0.1 kroD100 96.0 1.1 3 4,894.0 0.4 1.0 11 910.7 0.2 kroE100 154.3 1.7 1 606.7 0.9 1.6 7 717.2 1.0 rd100 137.9 1.6 5 1,002.3 1.2 1.6 10 449.4 1.1 eil101 154.5 2.3 3 227.6 1.6 2.2 6 1,668.5 2.1 Aggr. Inst. 1 - 756 52.8 1.7 401 469.3 1.4 1.6 536 381.8 1.6 aOnly the first 23 families. 37 •g (%): average optimality gap at termination for unsolved instances, computed as ub∗−lb∗ ub∗×100. A final block of the table, titled Improvement, summarizes the impact of replacing arrival-time variables with infeasible path elimination constraints, reporting the increase in the number of instances solved and the average reduction in root node gap. A close examination of Table 4reveals several noteworthy patterns regarding how the performance of both algorithms is influenced by different instance characteristics, and highlights the superiority of the Path Elimination model over the SEC-Enhanced approach. As expected, problem size, measured by the number of customers (n), has a major impact on tractability in both cases. All instances with n≤20 are solved to optimality by the SEC-Enhanced algorithm, and the Path Elimination model maintains full success up to this same threshold. However, the success rate steadily declines as nincreases, dropping sharply for n > 60. This decrease in solvability is accompanied by a significant increase in computational effort: larger instances require longer runtimes, explore more branch-and-bound nodes, and generate more cutting planes. A similar trend is observed with respect to the number of delivery tasks (m): nearly all instances with m≤100 are efficiently solved by the Path Elimination algorithm, whereas performance deteriorates for larger values in both solution rate and computational efficiency. Nevertheless, the termination gaps for unsolved instances remain relatively small, indicating that the algorithm continues to provide high-quality dual bounds even when optimality is not reached within the time limit. The number of service periods (|K|) also affects performance. Both algorithms exhibit the same qualitative behavior in this regard, differing primarily in the number of instances solved. Problems with three periods are consistently easier to solve than those with five, which demand similar computational effort but result in fewer optimal solutions. Regarding the instance-generation parameters, higher values of the fraction parameter f, which lead to denser delivery schedules, appear to facilitate the solution process. In particular, instances with f= 0.9 exhibit the lowest average root node and termination gaps for both algorithms. Finally, the results show that algorithmic performance remains stable across different values of ω, suggesting that the Path Elimination model is robust to variations in arrival-time flexibility. The last block of Table 4summarizes the impact of replacing arrival-time variables with infeasible path elimination constraints. The Path Elimination model consistently outperforms the SEC-Enhanced formulation, yielding lower average root node gaps that translate into marked improvements in the number of instances solved across all parameter settings. In addition, we provide two performance profiles in the sense of Dolan and Moré [15], which offer an intuitive visualization of relative performance across the benchmark set. These plots depict, for each algorithm, the proportion of instances it solves within a given factor of the best observed performance. As such, performance profiles provide a robust and widely accepted framework for comparing the efficiency of optimization algorithms over diverse problem instances. Figure 6displays the profiles of the three proposed algorithms on the reduced set of 414 instances (to allow the inclusion of the compact model), while Figure 7compares the SEC-Enhanced and Dual-driven 38 Table 4: Average computational performance of the implementation of the SEC-Enhanced and Path elimination model as a function of benchmark characteristics SEC-Enhanced Path elimination Improvement Parameters h(sec) #ins #sol r(%) t(sec) #nod #conn g(%) #sol r(%) t(sec) #nod #conn #ipec g(%) #sol r(%) n∈(0,20] 0.7 90 90 0.4 0.8 1,379 49 - 90 0.4 0.1 94 46 45 - 0 0.01 n∈(20,40] 2.0 144 128 1.5 305.5 115,839 299 1.4 144 1.3 277.1 35,521 286 591 - 16 0.21 n∈(40,60] 20.8 252 149 1.7 743.4 209,621 593 1.1 198 1.6 216.8 31,612 470 1,045 1.2 49 0.16 n∈(60,80] 104.6 108 10 3.2 93.2 16,742 415 2.2 25 2.9 1,210.8 90,325 954 1,789 2.2 15 0.26 n∈(80,100] 142.1 162 24 1.7 1,553.7 93,806 877 1.0 79 1.6 1,158.5 98,861 1,117 674 1.1 55 0.19 m∈(0,100] 2.8 261 245 1.1 202.3 87,234 260 0.8 258 1.0 90.4 14,974 188 384 0.1 13 0.12 m∈(100,200] 35.6 288 132 2.1 703.9 181,414 552 1.5 205 1.9 473.3 52,879 591 1,122 1.8 73 0.16 m∈(200,300] 85.9 138 21 2.1 1,985.5 153,848 786 1.3 64 1.9 998.8 89,861 996 906 1.7 43 0.22 m∈(300,400] 195.5 42 0 2.1 - - - 1.3 2 1.8 1,505.5 129,736 1,239 469 1.3 2 0.25 m∈(400,500] 328.3 27 3 1.3 1,331.1 18,040 1,049 0.9 7 0.9 2,479.8 105,228 2,099 858 0.7 4 0.39 Three-period 25.5 378 237 1.6 412.1 103,262 403 1.1 321 1.5 379.8 39,332 475 673 1.4 84 0.18 Five-period 80.1 378 164 1.8 551.9 147,139 371 1.6 215 1.7 384.7 41,047 456 828 1.6 51 0.17 f= 0.544.8 252 154 2.0 433.6 134,449 376 1.7 188 1.8 255.1 35,535 346 677 1.9 34 0.12 f= 0.741.1 252 114 2.0 369.5 136,421 399 1.6 167 1.9 564.3 54,314 539 905 1.8 53 0.17 f= 0.972.4 252 133 1.2 596.0 92,832 398 0.9 181 1.0 345.0 31,489 527 639 0.9 48 0.22 Low ω54.0 252 127 1.8 444.0 122,503 384 1.4 173 1.6 489.5 50,499 499 817 1.5 46 0.20 Med. ω52.6 252 135 1.7 578.4 143,779 408 1.4 178 1.6 287.3 35,503 447 688 1.5 43 0.17 High ω51.8 252 139 1.7 386.4 98,099 377 1.3 185 1.5 372.0 34,566 458 703 1.6 46 0.15 39 path elimination formulations on the complete set of 756 instances. These profiles are consistent with the aggregated picture in Table 2, providing additional evidence of the competitive performance of the pathelimination model. Figure 6illustrates the superior performance of the Path Elimination model compared to the other two proposed algorithms. As reflected in both the plot and Table 3, the Dual-Driven Path Elimination algorithm solves 89% of the first group of 414 instances, outperforming its competitors also in terms of speed (i.e., for smaller values of τ). The SEC-Enhanced model, while significantly improving upon the baseline compact formulation, still lags behind the path-based approach, solving 78% of the instances. In contrast, the Compact model exhibits substantially lower performance, solving only 44% of the instances and generally requiring longer computation times. Figure 7confirms the same performance pattern when comparing the SEC-Enhanced and Path Elimination models over the full benchmark set. However, the inclusion of 342 additional, larger instances leads to a proportional decrease in the percentage of solved cases for both models. 100101102103104105 0 0.2 0.4 0.6 0.8 1 Performance ratio τ Fraction of instances solved within a factor τof the best performance Compact SEC-Enhanced Path elimination Figure 6: Dolan–Moré performance profiles of the three proposed algorithms over the reduced set of 414 instances. The horizontal axis represents the performance ratio τ, while the vertical axis indicates the fraction of instances solved within τ times the best performance. 7.3. Benchmarking against the state of the art The most effective exact algorithm currently available for solving the ConTSP2 is the one proposed by Subramanyam and Gounaris [41], which follows a decomposition based framework, as previously discussed in Subsection 1.2. Their algorithm is capable of solving not only the problem addressed in this work (the ConTSP2), but also its variant that disallows idle times (the ConTSP1). 40 We are especially grateful to Anirudh Subramanyam for his willingness to clarify technical aspects of his work and for his consistently kind and constructive responses. His openness and availability greatly facilitated our comparative analysis. Declaration of generative AI and AI-assisted technologies in the writing process During the preparation of this work the authors used ChatGPT in order to improve the clarity, style, and linguistic accuracy of the manuscript. After using this tool/service, the authors reviewed and edited the content as needed and take full responsibility for the content of the publication. References [1] M. Alba Martínez, J.-F. Cordeau, M. Dell’Amico, M. Iori, A Branch-and-Cut Algorithm for the Double Traveling Salesman Problem with Multiple Stacks, INFORMS J. Comput. 25 (2013) 41–55, URL https: //doi.org/10.1287/ijoc.1110.0489. [2] D. L. Applegate, R. E. Bixby, V. Chvatál, W. J. Cook, Cuts from Blossoms and Blocks, chap. 7, Princeton University Press, 185–198, URL http://www.jstor.org/stable/j.ctt7s8xg.10, 2006. [3] N. Ascheuer, M. Fischetti, M. Grötschel, A polyhedral study of the asymmetric traveling salesman problem with time windows, Networks 36 (2000) 69–79, URL https://doi.org/10.1002/ 1097-0037(200009)36:2<69::AID-NET1>3.0.CO;2-Q. [4] N. Ascheuer, M. Fischetti, M. Grötschel, Solving the Asymmetric Travelling Salesman Problem with time windows by branch-and-cut, Math. Programming 90 (2001) 475–506, URL https://doi.org/10. 1007/PL00011432. [5] E. Balas, M. Fischetti, W. R. Pulleyblank, The precedence-constrained asymmetric traveling salesman polytope, Mathematical Programming 68 (1995) 241–265, URL https://doi.org/10.1007/ BF01585767. [6] R. Baldacci, A. Mingozzi, R. Roberti, New Route Relaxation and Pricing Strategies for the Vehicle Routing Problem, Oper. Res. 59 (2011) 1269–1283, URL https://doi.org/10.1287/opre.1110.0975. [7] J. F. Benders, Partitioning Procedures for Solving Mixed-variables Programming Problems, Numer. Math. 4 (1) (1962) 238–252, URL https://doi.org/10.1007/BF01386316. [8] N. Bélanger, G. Desaulniers, F. Soumis, J. Desrosiers, Periodic airline fleet assignment with time windows, spacing constraints, and time dependent revenues, Eur. J. Oper. Res. 175 (2006) 1754–1766, URL https://doi.org/10.1016/j.ejor.2004.04.051. [9] N. Christofides, J. E. Beasley, The period routing problem, Networks 14 (1984) 237–256, URL https: //doi.org/10.1002/net.3230140205. [10] G. Codato, M. Fischetti, Combinatorial Benders’ Cuts for Mixed-Integer Linear Programming, Oper. Res. 54 (2006) 756–766, URL https://doi.org/10.1287/opre.1060.0286. 47 [11] J. F. Cordeau, A Branch-and-Cut Algorithm for the Dial-a-Ride Problem, Oper. Res. 54 (2006) 573–586, URL https://doi.org/10.1287/opre.1060.0283. [12] G. Dantzig, R. Fulkerson, S. Johnson, Solution of a Large-Scale Traveling-Salesman Problem, Journal of the Operations Research Society of America 2 (1954) 393–410, URL http://www.jstor.org/stable/ 166695. [13] J. Desrosiers, M. E. Lübbecke, Branch-price-and-cut algorithms, Encyclopedia of Operations Research and Management Science. John Wiley & Sons, Chichester (2011) 109–131. [14] A. Dohn, E. Kolind, J. Clausen, The manpower allocation problem with time windows and job-teaming constraints: A branch-and-price approach, Comput. Oper. Res. 36 (2009) 1145–1157, URL https: //doi.org/10.1016/j.cor.2007.12.011. [15] E. D. Dolan, J. J. Moré, Benchmarking optimization software with performance profiles, Math. Programming 91 (2002) 201–213, URL https://doi.org/10.1007/s101070100263. [16] M. Drexl, Synchronization in Vehicle Routing –Survey of VRPs with Multiple Synchronization Constraints, Transportation Sci. 46 (2012) 297–316, URL https://doi.org/10.1287/trsc.1110.0400. [17] D. Díaz-Ríos, J. J. Salazar-González, Mathematical formulations for consistent travelling salesman problems, Eur. J. Oper. Res. 313 (2024) 465–477, URL https://doi.org/10.1016/j.ejor.2023.08. 021. [18] M. Fischetti, J. J. Salazar-González, P. Toth, Solving the Orienteering Problem Through Branch-AndCut, INFORMS J. Comput. 10 (1998) 133–148, URL https://doi.org/10.1287/ijoc.10.2.133. [19] D. Goeke, R. Roberti, M. Schneider, Exact and Heuristic Solution of the Consistent Vehicle-Routing Problem, Transportation Sci. 53 (2019) 1023–1042, URL https://doi.org/10.1287/trsc.2018.0864. [20] L. Gouveia, J. M. Pires, The asymmetric travelling salesman problem and a reformulation of the Miller–Tucker–Zemlin constraints, Eur. J. Oper. Res. 112 (1999) 134–146, URL https://doi.org/ 10.1016/S0377-2217(97)00358-5. [21] C. Groër, B. Golden, E. Wasil, The Consistent Vehicle Routing Problem, Manufacturing Service Oper. Management 11 (2009) 630–643, URL https://doi.org/10.1287/msom.1080.0243. [22] S. Gélinas, M. Desrochers, J. Desrosiers, M. M. Solomon, A new branching strategy for time constrained routing problems with application to backhauling, Ann. Oper. Res. 61 (1995) 91–109, URL https: //doi.org/10.1007/BF02098283. [23] I. Ioachim, J. Desrosiers, F. Soumis, N. Bélanger, Fleet assignment and routing with schedule synchronization constraints, Eur. J. Oper. Res. 119 (1999) 75–90, URL https://doi.org/10.1016/ S0377-2217(98)00343-9. 48 [24] I. Ioachim, S. Gélinas, F. Soumis, J. Desrosiers, A dynamic programming algorithm for the shortest path problem with time windows and linear node costs, Networks 31 (1998) 193–204, URL http: //doi.org/10.1002/(SICI)1097-0037(199805)31:3<193::AID-NET6>3.0.CO;2-A. [25] M. Iori, J. Riera-Ledesma, Exact algorithms for the double vehicle routing problem with multiple stacks, Comput. Oper. Res. 63 (2015) 83–101, URL https://doi.org/10.1016/j.cor.2015.04.016. [26] A. A. Kovacs, B. L. Golden, R. F. Hartl, S. N. Parragh, Vehicle routing problems in which consistency considerations are important: A survey, Networks 64 (2014) 192–213, URL https://doi.org/10.1002/ net.21565. [27] A. A. Kovacs, B. L. Golden, R. F. Hartl, S. N. Parragh, Vehicle routing problems in which consistency considerations are important: A survey, Networks 64 (2014) 192–213, URL https://doi.org/10.1002/ net.21565. [28] A. A. Kovacs, S. N. Parragh, R. F. Hartl, A template-based adaptive large neighborhood search for the consistent vehicle routing problem, Networks 63 (2014) 60–81, URL https://doi.org/10.1002/net. 21522. [29] G. Laporte, J. Riera-Ledesma, J. J. Salazar-González, A Branch-and-Cut Algorithm for the Undirected Traveling Purchaser Problem, Oper. Res. 51 (2003) 940–951, URL https://doi.org/10.1287/opre. 51.6.940.24921. [30] H. Lespay, K. Suchan, A case study of consistent vehicle routing problem with time windows, Int. Trans. Oper. Res. 28 (2021) 1135–1163, URL https://doi.org/10.1111/itor.12885. [31] Z. Luo, H. Qin, W. Zhu, A. Lim, Branch-and-price-and-cut for the manpower routing problem with synchronization constraints, Naval Res. Logist. 63 (2016) 138–171, URL http://doi.org/10.1002/ nav.21683. [32] S. Mancini, M. Gansterer, R. F. Hartl, The collaborative consistent vehicle routing problem with workload balance, Eur. J. Oper. Res. 293 (2021) 955–965, URL https://doi.org/10.1016/j.ejor.2020. 12.064. [33] P. C. Nolz, N. Absi, D. Feillet, C. Seragiotto, The consistent electric-Vehicle routing problem with backhauls and charging management, Eur. J. Oper. Res. 302 (2022) 700–716, URL https://doi.org/ 10.1016/j.ejor.2022.01.024. [34] G. Reinelt, TSPLIB—A Traveling Salesman Problem Library, ORSA J. Comput. 3 (1991) 376–384, URL https://doi.org/10.1287/ijoc.3.4.376. [35] J. Riera-Ledesma, J. J. Salazar-González, A branch-and-cut algorithm for the continuous error localization problem in data cleaning, Comput. Oper. Res. 34 (2007) 2790–2804, URL https://doi.org/ 10.1016/j.cor.2005.10.016. 49 [36] J. Riera-Ledesma, J. J. Salazar-González, Selective routing problem with synchronization, Comput. Oper. Res. 135 (2021) 105465, URL https://doi.org/10.1016/j.cor.2021.105465. [37] B. Sarasola, K. F. Doerner, Adaptive large neighborhood search for the vehicle routing problem with synchronization constraints at the delivery location, Networks 75 (2020) 64–85, URL https://doi. org/10.1002/net.21905. [38] R. Soares, A. Marques, P. Amorim, S. N. Parragh, Synchronisation in vehicle routing: Classification schema, modelling framework and literature review, Eur. J. Oper. Res. 313 (2024) 817–840, URL https: //doi.org/10.1016/j.ejor.2023.04.007. [39] F. Stavropoulou, The Consistent Vehicle Routing Problem with heterogeneous fleet, Comput. Oper. Res. 140 (2022) 105644, URL https://doi.org/10.1016/j.cor.2021.105644. [40] A. Subramanyam, C. E. Gounaris, A branch-and-cut framework for the consistent traveling salesman problem, Eur. J. Oper. Res. 248 (2016) 384–395, URL https://doi.org/10.1016/j.ejor.2015.07. 030. [41] A. Subramanyam, C. E. Gounaris, A Decomposition Algorithm for the Consistent Traveling Salesman Problem with Vehicle Idling, Transportation Sci. 52 (2018) 386–401, URL https://doi.org/10.1287/ trsc.2017.0741. 50