A New Class of Compact Formulations for Vehicle Routing Problems

A VRPTW problem instance is described using the following terms: (a) a depot located in space; (b) a set of customers located in space, each of which has an integer demand, service time, and time window for when service starts; and (c) an unlimited number of homogeneous vehicles with integer capacity.

A solution to VRPTW assigns vehicles to routes, where each route satisfies the following properties: (1) The route starts and ends at the depot. (2) The total demand of customers serviced on the route does not exceed the vehicle’s capacity. (3) The route cost is the total distance traveled. (4) The vehicle leaves a customer immediately after servicing that customer. (5) Service starts at a given customer within the time window of that customer. (6) A customer is serviced at most once in a route.

The VRPTW solver selects a set of routes with the goal of minimizing the total distance traveled while ensuring that each customer is serviced exactly once across all selected routes.

We use $N$ to denote the set of customers, which we index by $u$. We use $N^{+}$ to denote $N$ augmented with the starting/ending depot (which are co-located), and denoted $\alpha,\bar{\alpha}$ respectively. Each customer $u\in N$ has a demand denoted $d_{u}$. We define the demand at the depots as $0$, which we write formally as $d_{\alpha}=d_{\bar{\alpha}}=0$. We are given an unlimited number of homogeneous vehicles with capacity $d_{0}$ that start at the depot at time $t_{0}$ where $t_{0}$ is the length of the time horizon of optimization. Note that this is without loss of generality as we can easily apply a cost (distance) term to each departure/arrival from the depot to restrict the number of vehicles deployed. We use $t^{+}_{u}$ and $t^{-}_{u}$ to denote the earliest/latest times for the service at customer $u$ to begin. The service windows for the start/end depot are denoted by $t_{\alpha}^{+}=t_{\alpha}^{-}=t_{\bar{\alpha}}^{+}=t_{0}$ and $t_{\bar{\alpha}}^{-}=0$. Servicing a given customer $u\in N$ takes $t_{u}^{*}$ units of time, where $t_{\alpha}^{*}=t_{{\bar{\alpha}}}^{*}=0$.

For each $u\in N^{+},v\in N^{+}$ we use $c_{uv}=t_{uv}$ to denote the cost and time required to travel from $u$ to $v$. For simplicity, $t_{uv}$ is altered to include the service time at customer $u$, which is $t_{u}^{*}$. Thus, we set $t_{uv}\leftarrow t_{uv}+t_{u}^{*}$, then set $c_{uv}\leftarrow t_{uv}$, which offsets the cost of a solution by a constant.

We now describe the two-index formulation for the VRPTW of interest. We use the following decision variables. We set binary decision variable $x_{uv}=1$ to indicate that a vehicle departs $u$ and travels immediately to $v$. $x_{uv}$ is defined for every pair of customers/depots, where $u$ is not the ending depot and $v$ is not $u$ or the starting depot. The term $x_{uv}$ is also not defined for cases where the route that proceeds from beginning to end $[\alpha,u,v,\bar{\alpha}]$ is infeasible due to time or capacity. The set of existing $x_{uv}$ terms are denoted $E^{*}$.

We use continuous decision variable $\tau_{u}\in\mathbb{R}_{0+}$ to denote the amount of time remaining when service starts at customer $u$. We enforce the bounds on the start of the service with the following inequalities: $t^{+}_{u}\geq\tau_{u}\geq t^{-}_{u}\quad\forall u\in N$.

We use continuous decision variable $\delta_{u}\in\mathbb{R}_{0+}$ to denote the capacity remaining in the vehicle immediately before servicing customer $u$. We enforce that sufficient capacity remains to service customer $u$ with the following inequalities $d_{0}\geq\delta_{u}\geq d_{u}\quad\forall u\in N$.

We now write the two-index form for the VRPTW as a compact MILP :

In (1a) we minimize the total distance traveled by all vehicles. In (1b), we enforce that each customer is serviced exactly once. In (1c), we enforce that a vehicle must leave each customer once. In (1d), we enforce that if $u$ is immediately preceded by $v$ in a route, then the amount of demand remaining immediately before servicing $u$ (denoted $\delta_{u}$) is no greater than that of $v$ (denoted $\delta_{v}$) minus $d_{v}$; and is otherwise unrestricted. Observe that in (1d) that $(d_{0}+d_{v})(1-x_{vu})$ operates as a big-M term, which makes the inequality inactive when $x_{vu}=0$. In (1e), we enforce that if $u$ is preceded by $v$ in a route, then the amount of time remaining immediately before servicing $u$ (denoted $\tau_{u}$) is no greater than that of $v$ (denoted $\tau_{v}$) minus $t_{vu}$; and is otherwise unrestricted. Observe that in (1e) that $(t^{+}_{u}+t_{vu})(1-x_{vu})$ operates as a big-M term, which makes the inequality inactive when $x_{vu}=0$.

In (1f) we enforce that a minimum number of vehicles exit the depot; equal to a lower bound on the number of vehicles required to service all demand. This lower bound is the sum of the demands divided by the vehicle capacity, then rounded up to the nearest integer. This is an optional constraint that is not common in the CG literature.

It is well know that this formulation can lead to weak LP bounds because of the following issues. (1) It does not explicitly eliminate sub-tours in the LP solution. (2) Capacity and time feasibility are not well enforced.

In this section, we provide additional variables and associated constraints using Local Area (LA) arcs, which were introduced in the arXiv paper (Mandal et al., 2022). LA-arcs are designed to eliminate cycles of customers localized in space as described by the $x$ variables used in (1). As a precursor, we provide some notation for describing LA-arcs, which we illustrate in Fig 1. For each $u\in N$ we use $N^{\odot}_{u}$ to denote the set of customers nearby $u$, not including $u$, called the LA-neighbors of $u$. Specifically $N^{\odot}_{u}$ is the set of $K$ closest customers to $u$ in space that are reachable from $u$ considering time and capacity; for user-defined $K$.¹¹1Other definitions can be used such as excluding from $N^{\odot}_{u}$ customers for which the route ordered as follows $\alpha,u,v,u,\bar{\alpha}$ is infeasible by time/capacity. For ease of notation we define $N^{\odot}_{+u}=\{N^{\odot}_{u}\cup u\}$. We use $N^{\odot\rightarrow}_{u}$ to denote the subset of $N^{+}$ not in $N^{\odot}_{+u}$ or the starting depot meaning $N^{\odot\rightarrow}_{u}=N^{+}-(N^{\odot}_{+u}\cup\alpha)$. For each $u\in N,v\in N^{\odot\rightarrow}_{u},\hat{N}\subseteq N^{\odot}_{u}$, we use $R^{+}_{u\hat{N}v}$ to denote the subset of all elementary sequences of customers (called orderings) starting at $u$; ending at $v$; including as intermediate customers $\hat{N}$; that are feasible with regard to time/capacity. Not all sets of intermediate customers $\hat{N}\subseteq N_{u}^{\odot}$ are feasible given $u,v$ due to time/capacity constraints. For any given $u,v$ and set of intermediate customers in $\hat{N}\subseteq N_{u}$, we define $R_{u\hat{N}v}\subseteq R^{+}_{u\hat{N}v}$ to contain only those orderings on the efficient frontier with respect to (a) travel distance, (b) latest feasible departure time from $u$ and (c) earliest feasible departure time from $v$. Feasibility of (b),(c) is a consequence of the time windows for the customers. Thus, orderings are preferred with (a) lower cost, (b) that can be started later, or (c) that can be finished earlier. In practice, $R_{u\hat{N}v}$ is a small subset of all possible orderings $R^{+}_{u\hat{N}v}$ as most do not lie on the efficient frontier. We use $R_{u}$ to denote the union of $R_{u\hat{N}v}$ over all $v,\hat{N}$ where we clip off the final term in the ordering (meaning $v$). Observe that any optimal integer solution to VRPTW must select a route containing an ordering in $R_{u}$ followed immediately by some $v\in N^{\odot\rightarrow}_{u}$ for each $u\in N$. In practice, $R_{u}$ is smaller than the sum of the sizes of the $R_{u\hat{N}v}$ terms as many terms are replicated between $R_{u\hat{N}v}$ sets. We provide a dynamic programming approach for generating $R_{u}$ in Section 6.

We describe an ordering $r\in R_{u}$ as follows. For any $w\in N^{\odot}_{+u}$, $v\in N^{\odot}_{u}$ we define $a_{wvr}=1$ if $w$ immediately precedes $v$ in ordering $r$, and otherwise define $a_{wvr}=0$. For any $w\in N^{\odot}_{+u}$ we define $a_{w*r}=1$ if $w$ is the final customer in ordering $r$, and otherwise define $a_{w*r}=0$.

We use decision variables $y_{r}\in\mathbb{R}_{0+}\;\forall r\in R_{u},\;u\in N$ to describe a selection of orderings to use in our solution. We set $y_{r}=1$ for $r\in R_{u}$ to indicate that after servicing customer $u$ the sequence of customers described by ordering $r$ is used after, which a member of $N_{u}^{\odot\rightarrow}$ follows immediately, and otherwise set $y_{r}=0$. We use $E^{*\odot}_{u}\subseteq E^{*}$ to denote the set of edges that can be contained in an ordering starting at $u$; thus $E^{*\odot}_{u}=\{wv\in E^{*},\mbox{ s.t. }w\in N^{\odot}_{+u},v\in N^{\odot}_{u}-w\}$. We use $E^{*\odot}_{uw}$ to denote the set of edges in $E^{*}$ that can succeed the final customer $w$ in an ordering starting at $u$; thus $E^{*\odot}_{uw}=\{wv\in E^{*};w\in N^{\odot}_{+u},v\in N^{\odot\rightarrow}_{u}\}$.

We now add the following equations to (1), which adds minimization over $y$ to enforce that the solution to the $x$ variables is consistent with the solution $y$, for which we provide an exposition below the equations.

In (2a) we enforce that one ordering in $R_{u}$ is selected to describe the activities succeeding $u$ for each $u\in N$. In (2b) we enforce that if an ordering $r\in R_{u}$ is selected in which $w$ is immediately followed by $v$; then it must be the case that $x_{wv}$ is selected meaning that $w$ is succeeded by $v$ in the solution over $x$. In (2c) we enforce that if an ordering $r\in R_{u}$ is selected in which $w$ is the final customer in the ordering then, $w$ must be followed by a customer (or depot) that is not an LA-neighbor of $u$ (and not $u$) as described by $x$. Observe that if $x$ is binary then $y$ must be binary as well. The terms $y_{r}$ are not present in the objective for optimization. The use of LA-arcs to remove cycles, is a key contribution as it allows the user to grow the size of the formulation as the computational capabilities of MILP solvers improve over time.

It is possible to show (see Section 4.5.2) that the number of LA-neighbors required can be less than that of $N^{\odot}_{u}$ for some $u$. To this end, we introduce $N_{u}\subseteq N_{u}^{\odot}$ to denote this reduced set. We use $N_{+u}$ to denote $N_{u}\cup u$. We use $N^{\rightarrow}_{u}$ to denoted the subset of $N^{+}$ that does not lie in $N_{+u}\cup\alpha$. We use $E^{*}_{u}\subseteq E^{*}$ to denote the set of edges that can be contained in an ordering starting at $u$, thus $E^{*}_{u}=\{wv\in E^{*},\mbox{ s.t. }w\in N_{+u},v\in N_{u}-w,\}$. We use $E^{*}_{uw}$ to denote the set of edges in $E^{*}$ that can succeed the final customer $w$ in an ordering starting at $u$, thus $E^{*}_{uw}=\{wv\in E^{*};v\in N^{\rightarrow}_{u}\}$. We write (2) using LA-neighbors $N_{u}$ below.

Note that $R_{u}$ is determined by $N^{\odot}_{u}$ not $N_{u}$; while the set of constraints enforced in (3b)-(3c) is defined by $N_{u}$. Observe that this can create redundant $y$ terms. These are simply removed before optimization.

This section considers using capacity/time discretization to tighten the LP relaxation in (1). As mentioned earlier, this discretization is inspired by the work of (Boland et al., 2017, 2019). Earlier discretization methods can be found in (Appelgren, 1969, 1971; Wang and Regan, 2002).

Consider that we have a partition of capacity associated with each customer $u$ denoted $D_{u}$ into continuous “buckets,” e.g. $D_{u}=[[d_{u},d_{u}+3],[d_{u}+4,d_{u}+5],[d_{u}+6,d_{u}+9]...[d_{0}-4,d_{0}]]$. We index the buckets in $D_{u}$ from smallest to greatest remaining demand with $k$; where $k=1$ is the bucket that contains $d_{u}$ and $k=|D_{u}|$ contains $d_{0}$. The lower/upper bound values of the $k^{\prime}th$ bucket for customer $u$ are referred to as $d^{-}_{(u,k)},d^{+}_{(u,k)}$ respectively.

Consider a directed unweighted graph $G^{D}$ (which we also use for node-set) with edge set $E^{D}$. For each $u\in N$, $k\in D_{u}$ we create one node $i=(u,k)$ for which we alternatively write as $(u_{i},d^{-}_{i},d^{+}_{i})$ where $u_{i}\leftarrow u,d^{-}_{i}\leftarrow d^{-}_{k},d^{+}_{i}\leftarrow d^{+}_{k}$. There is a single node $(\alpha,d_{0},d_{0})$ also denoted $(\alpha,1)$ associated with the starting depot and a single node $(\bar{\alpha},0,d_{0})$ also denoted $(\bar{\alpha},1)$ associated with the ending depot. We now define the edges $E^{D}$.

• •

  For each $u$ we create a directed edge from $(\alpha,1)$ to $(u,|D_{u}|)$. Traversing this edge indicates that a route starts with $u$ as its first customer.

• •

  For each $u$ we create a directed edge from $(u,1)$ to $(\bar{\alpha},1)$. Traversing this edge indicates that a route ends after having $u$ as its final customer.

• •

  For each $i=(u,k)$ for $k>1$ we connect $i$ to $j=(u,k-1)$. Traversing this edge indicates that the vehicle servicing $u$ will not use all possible capacity. These edges can be thought of as dumping extra capacity.

• •

  For each pair of nodes $i=(u,k),j=(v,m)$ we connect $i$ to $j$ if the following are satisfied: $u\neq v$, $d^{+}_{i}-d_{u}\geq d^{-}_{j}$ and if $k>1$ we also require that $d^{-}_{j}>d^{+}_{(u,k-1)}-d_{u}$. Traversing this arc indicates that a vehicle has between $d^{-}_{i}$ and $d^{+}_{i}$ units of capacity remaining before servicing $u$ and, upon servicing, $u$ travels directly to $v$ where it has between $d^{-}_{j}$ and $d^{+}_{j}$ units of capacity remaining before servicing $v$. The second condition $d^{-}_{j}>d^{+}_{(u,k-1)}-d_{u}$ ensures that redundant edges are not generated.

We now describe flow formulation on the discretized capacity. This flow governs the amount of the resource capacity available at a given node. The following constraints are written below using decision variable $z^{D}_{ij}\in\mathbb{R}_{0+}$ to describe the edges traversed .

In (4a), we enforce that the solution $z^{D}$ is feasible with regard to capacity by enforcing that the flow incoming to a node is equal to the flow outgoing. In (4b), we enforce that the solution $x$ is consistent with the solution capacity graph $z^{D}$.

We now describe flow formulation on discretized time as we previously did for capacity. Time buckets and associated graph and edge sets are denoted $T_{u},G^{T},E^{T}$ respectively. The constraints are defined analogously to capacity (with minimum/maximum values of $t^{-}_{u}$ and $t^{+}_{u}$ respectively in place of $d_{u}$ and $d_{0}$). We write the corresponding constraints below. We use non-negative decision variables $z^{T}_{ij}$ to describe the movement of vehicles. We set $z^{T}_{ij}=1$ to indicate that a vehicle leaves $u_{i}$ with time remaining in between $t^{-}_{i}$ and $t^{+}_{i}$ and leaves at $u_{j}$ with time remaining between $t^{-}_{j}$ and $t^{+}_{j}$. We write the constraints for time analogous to (4a) and (4b) below.

We refer to the finest possible granularity of capacity and time buckets for each $u\in N$ as $D^{\odot}_{u},T^{\odot}_{u}$.

Given a set of LA-neighbors, time buckets, and capacity buckets for each $u\in N$ denoted $N_{u},T_{u},D_{u}$ respectively, we write the MILP enforcing no-localized cycles in space (via (3)), capacity bucket feasibility (via (4)), and time bucket feasibility (via (5)) as $\bar{\Psi}(\{N_{u},T_{u},D_{u}\},\forall u\in N)$ below, with its LP relaxation denoted $\Psi(\{N_{u},T_{u},D_{u}\},\forall u\in N)$.

To be able to use an off-the-shelf MILP solver to solve for (6) efficiently, we must have parameterization $(\{N_{u},T_{u},D_{u}\}\forall u\in N)$ that is as tight as possible and small as possible in terms of using fewer variables and constraints. Since the tightest LP expressible is $\Psi(\{N^{\odot}_{u},T^{\odot}_{u},D^{\odot}_{u}\}\forall u\in N)$ we want $\Psi(\{N_{u},T_{u},D_{u}\}\forall u\in N)$ to be equal to $\Psi(\{N^{\odot}_{u},T^{\odot}_{u},D^{\odot}_{u}\}\forall u\in N)$ and refer to such parameterizations as sufficient.

We want the LP $\Psi(\{N_{u},T_{u},D_{u}\}\forall u\in N)$ and hence each node in the MILP $\bar{\Psi}(\{N_{u},T_{u},D_{u}\}\forall u\in N)$ to be easily solved. One key condition for such parameterizations is that of being parsimonious; any further contraction of the sets in the parameterization decreases the LP value of the optimal solution. Another way to say this is that the parameterization has the widest buckets and smallest LA-neighborhoods possible, given fixed LP tightness, leading to maximally efficient solutions to the LPs at each node of the branch & bound tree of the corresponding MILP. We refer to LPs that contain these two properties that as sufficient and parsimonious parameterization (SPP) and in Section 5 consider the construction of such LPs.

We now consider some sufficient (but not always necessary) conditions of sufficient parameterizations. To this end we consider a candidate optimal LP solution $(x^{*},y^{*},z^{*})$ given parameterization $(\{N^{*}_{u},T^{*}_{u},D^{*}_{u}\}\forall u\in N)$.

We now consider some necessary but not always sufficient properties of parsimonious parameterizations.

Recall from Section 4.5.1 that for any $i,j\in E^{D}$ for which $u_{i}=u_{j}$ if $\pi^{D}_{i}=\pi^{D}_{j}$ then we can merge nodes $i,j$ without loosening the bound and remove $d^{+}_{j}$ from $W^{D}_{u_{j}}$. The same logic applies to time. For any $i,j\in E^{T}$ for which $u_{i}=u_{j}$ if $\pi^{T}_{i}=\pi^{T}_{j}$, we can merge nodes $i,j$ without loosening the bound and remove $t^{+}_{j}$ from $W^{T}_{u_{j}}$.

To enforce bucket feasibility with respect to capacity we do the following. For each $i,j$ s.t. $z^{D}_{ij}>0$ where $i=(u_{i},d^{-}_{i},d^{+}_{i})$ and $j=(u_{j},d^{-}_{j},d^{+}_{j})$ and $j\neq\bar{\alpha}$ and $u_{i}\neq u_{j}$: We add to $W^{D}_{u_{j}}$ the term $d^{+}_{i}-d_{u_{i}}$ (if this term is not already present in $W^{D}_{u_{j}}$ ). Observe that this operation can never loosen the relaxation as undoing it simply corresponds to requiring that $\pi^{D}_{i}=\pi^{D}_{j}$ in the dual.

The same approach applies to enforcing bucket feasibility with respect to time. For each $i,j$ s.t. $z^{D}_{ij}>0$ where $i=(u_{i},t^{-}_{i},t^{+}_{i})$ and $j=(u_{j},t^{-}_{j},t^{+}_{j})$ and $j\neq\bar{\alpha}$: We add to $W^{T}_{u_{j}}$ the term $\min(t^{+}_{i}-t_{u_{i},u_{j}},t^{+}_{j})$ (if this term is not already present in $W^{T}_{u_{j}}$ ). Observe that this operation can never loosen the relaxation as undoing it simply corresponds to requiring that $\pi^{T}_{i}=\pi^{T}_{j}$ in the dual.

Recall from (9) that we can set the number of LA-neighbors to the smallest number for which there is an associated non-zero dual value without loosening the bound. Thus our contraction operation on LA-neighborhoods is written as follows.

To expand the neighborhoods, we can iterate over $u\in N$, then, for each $u$, identify the smallest number of LA-neighbors required to induce a violation of the current LP. To determine if a given number of LA-neighbors $k$ induces a violation of the current LP we simply check to see if a fractional solution $y_{r}$ exists to satisfy (6i), (6g) and (6h) where we fix the $x$ terms to their current values and $N_{u}$ to $N^{k}_{u}$.

A simpler approach is as follows. We can periodically simply reset the number of LA-neighbors to the maximum number $|N^{\odot}_{u}|$ for each $u\in N$. Next, we solve the LP and apply a contraction operation described in (10). Our experiments employ this approach when expansion operations on capacity/time buckets have not recently tightened the LP bound.

We now consider the construction of an SPP by iteratively solving the LP relaxation and tightening the LP relaxation to aim toward sufficiency; each time we tighten the relaxation, we decrease the size of the parameterization with parsimony. We provide our algorithm in Alg 1 with exposition below.