# Solving a Continent-Scale Inventory Routing Problem at Renault

Louis Bouvier, Guillaume Dalle, Axel Parmentier, Thibaut Vidal

► To cite this version:

Louis Bouvier, Guillaume Dalle, Axel Parmentier, Thibaut Vidal. Solving a Continent-Scale Inventory Routing Problem at Renault. Transportation Science, 2023, 58 (1), 10.1287/trsc.2022.0342. hal-03766779v2

HAL Id: hal-03766779

https://hal.science/hal-03766779v2

Submitted on 2 Mar 2026

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Solving a Continent-Scale Inventory Routing Problem at Renault

Louis Bouvier*¹,², Guillaume Dalle¹, Axel Parmentier¹, and Thibaut Vidal³

¹CERMICS, Ecole des Ponts, France

²Groupe Renault

³Polytechnique Montréal, Canada

October 2, 2023

## Abstract

This paper is the fruit of a partnership with Renault. Their reverse logistic requires solving a continent-scale multi-attribute inventory routing problem (IRP). With an average of 30 commodities, 16 depots, and 600 customers spread across a continent, our instances are orders of magnitude larger than those in the literature. Existing algorithms do not scale, so we propose a large neighborhood search (LNS). To make it work, (1) we generalize existing split delivery vehicle routing problem and IRP neighborhoods to this context, (2) we turn a state-of-the-art matheuristic for medium-scale IRP into a large neighborhood, and (3) we introduce two novel perturbations: the reinsertion of a customer and that of a commodity into the IRP solution. We also derive a new lower bound based on a flow relaxation. In order to stimulate the research on large-scale IRP, we introduce a library of industrial instances. We benchmark our algorithms on these instances and make our code open-source. Extensive numerical experiments highlight the relevance of each component of our LNS.

## 1 Introduction

The inventory routing problem (IRP) arises when a supplier manages the delivery of commodities to its customers on a multiple-day horizon in a centralized manner (Archetti and Speranza 2016). It consists in planning routes to deliver commodities from depots to customers with the objective of minimizing inventory and routing costs. This NP-hard problem has received significant attention in the operations research literature over the past 40 years.

The present paper is motivated by a partnership with Renault, a major European car manufacturer who must routinely solve IRP instances of unprecedented continental scale and complexity as part of their reverse logistic problem. Indeed, they receive car parts from suppliers at their plants in packaging, and reuse the latter, which implies the need for reverse packaging logistics. The goal of our partnership is to redesign their IRP algorithm. This is challenging because of (1) the size of the resulting instances, with 600 customers and 16 depots on average, (2) the 30 different commodities involved, and (3) the specific challenges that arise from the geography and timescale. When depots and customers are scattered across a whole continent, travel times can last up to ten days. Beyond requiring a long horizon, 21 days in our case, this makes the problem more difficult because classic decoupling results on the IRP, which were exploited in previous algorithms, are no longer valid. For instance, changing the order of the customers along a route impacts the arrival day at each customer and therefore the customer inventory levels. Hence, routes with suboptimal routing cost may be better because of inventory cost, which is not usually the case with the IRP. It can be compared with the continuous-time IRP discussed in M. Savelsbergh and Song (2008) or Lagos, Boland, and

<footnote>*Corresponding author: louis.bouvier@enpc.fr</footnote>
M. Savelsbergh (2020). Touzout, Ladier, and Hadj-Hamou (2022) go further in this direction by considering a time-dependent setting in which travelling time depends on the departure time. We do not include this specificity in our work because congestion does not significantly impact long trips. Finally, our partner's supply chain process requires that (4) the solution algorithm should not take more than 90 minutes on our computing cluster.

State-of-the-art exact algorithms rely on branch-and-cut (Archetti, Bertazzi, et al. 2007; Coelho and Laporte 2013; Manousakis et al. 2021) and branch-and-price-and-cut methods (Desaulniers, Rakke, and Coelho 2015) with dedicated valid inequalities. They can optimally solve single-commodity single-depot instances with up to 50 customers, but are not appropriate for our large-scale setting.

Typical heuristics include route-based matheuristics (Fischetti and Fischetti 2016), decomposition matheuristics, and metaheuristics. In this field, Bertazzi et al. (2019) and Archetti, Boland, and Speranza (2017) are route-based matheuristics. The main idea is to reduce the size of the mixed-integer linear program (MILP) formulation of the IRP, by selecting promising routes heuristically. Although Bertazzi et al. (2019) is dedicated to the multi-depot case, neither of the two papers handles the multicommodity aspect we must face, and their largest instances have up to six days horizon, six depots and 50 customers. The methods cannot be applied directly in our context, because the MILP remains too large, even when we restrict ourselves to “promising routes”. We instead adapt their principle to our setting, leading to the “reload fixed-path vehicles” subroutine.

Another common approach is to tackle the IRP through a decomposition (Campbell and M. W. P. Savelsbergh 2004; Cordeau et al. 2015). For instance, first set the quantities to be sent, and then create the routes to respect them. The largest instances solved with this two-step method have a single depot, up to five commodities, a six days horizon and 50 customers. We also adapt it to our setting and use it as an initialization heuristic.

Some metaheuristics have been designed for real-world IRP. For instance, Benoist et al. (2011) introduce a randomized local search to address a large-scale single-commodity IRP with pickups, time windows, driver safety and other constraints that are specific to their use case, but less relevant to ours. Su et al. (2020) address a real-world IRP from the ROADEF IRP-Challenge 2016 available at https://www.roadef.org/challenge/2016/fr/ with a large neighborhood search based on mathematical programming. As for Benoist et al. (2011), the single-commodity formulation with additional constraints is not adapted to our multicommodity context. Other large neighborhoods are introduced in Nolz, Absi, and Feillet (2014) for a single-depot single-commodity stochastic version of the IRP. They are based on perturbation ideas such as the removal of every customer visit on a particular day, followed by a best-insertion policy. This process of removal and insertion is at the core of our work, but we leverage MILP formulations for the insertion. A kernel search heuristic based on a preliminary tabu search is considered in Archetti, Guastaroba, et al. (2021). In this framework, smaller MILPs with increasing size are solved iteratively to improve an initial solution. Single-depot, single-commodity instances with up to six days horizon and 200 customers are solved. The study Coelho, De Maio, and Laganà (2020) fixes a part of the decision variables, this time based on the problem's main “axes” – that is to say the different types of sites involved in the deliveries, the routing and the inventory aspects. It solves reduced MILPs in a variable neighborhood search, and defines the multi-attribute IRP as a multicommodity, multi-depot and multi-vehicle IRP. This constraint structure is the closest to ours. Nonetheless, the instances are smaller, with up to six days horizon, 50 customers, six depots and three commodities. Besides, this study is restricted to one-day routes, whereas we explicitly deal with routes that last multiple days, which creates an additional combinatorial challenge. Another difference is that it incorporates a heterogeneous fleet with three distinct vehicle types, whereas we consider a homogeneous infinite fleet of vehicles. Other approaches follow this idea of fixing a part of the solution. For instance, Chitsaz, Cordeau, and Jans (2019) use a three-phase decomposition matheuristic to address the assembly routing problem. The authors also derive ways to adapt it to the single-depot single-commodity IRP. Similarly, Vadseth, Andersson, and Stålhane (2021) introduce an iterative matheuristic for the single-depot single-commodity IRP. The method alternates between generating a small set of routes and solving a path-flow formulation given the fixed routes. Our approach shares this idea of updating sequentially the set of routes and the quantities delivered. Nonetheless, due to our multicommodity aspect, optimizing the
quantities given the routes remains an *NP*-hard problem in our case. We also introduce additional perturbations. Therefore, to the best of our knowledge, no algorithm is known to properly scale to our instances. In this context, our main contributions are the following:

1. We introduce two new large-scale perturbation neighborhoods designed for the multi-attribute IRP. They are based on well-solved MILP formulations, and enable to escape from local minima.

2. We design an efficient large neighborhood search (LNS) built upon these large neighborhoods. We generalize 12 Traveling Salesman Problem (TSP), and Split Delivery Vehicle Routing Problem (SD-VRP) (see, e.g., Dror and Trudeau (1990) and Archetti and Speranza (2008)) neighborhoods from the literature to our IRP context. We propose a new matheuristic inspired by Archetti, Boland, and Speranza (2017) and Bertazzi et al. (2019) and adapted to our large-scale setting.

3. We compute a new lower bound based on a linear program (LP) relaxation (one flow per commodity). To the best of our knowledge, this relaxation is not considered in the literature. We do not expect the bound to be tight, which is a feature of every relaxation in the IRP literature. But it is useful to compare algorithm performance on instances with distinct scales.

4. We provide a publicly available library of realistic multi-attribute IRP instances of a continent scale, as an incentive for further research on the topic.

5. We give access to our open-source Julia (Bezanson et al. 2017) package that implements the ideas of the present paper. It is available on GitHub (Bouvier et al. 2023c), and we also provide the library of instances (Bouvier et al. 2023a) and the solutions we obtained (Bouvier et al. 2023b).

6. We proceed to extensive numerical experiments. Since no algorithm is known to scale to our context, we compare the adapted route-based matheuristic and our large neighborhood search.

We precisely define the problem we consider in Section 2. We then provide an overview of the different solution processes in Section 3. The three next sections 4, 5 and 6 detail the algorithms. We finally report our numerical experiments in Section 7.

## 2 Problem Description

We consider a rich variant of the inventory routing problem. A supplier manages the delivery of several commodities to its customers on a multiple-day horizon in a centralized manner (Archetti and Speranza 2016). The supplier has to plan routes to deliver commodities from its depots to customers with the objective of minimizing inventory and routing costs. Every day, depots release and customers demand commodities. Routes may last several days due to large travel distances. For each route, the starting depot, the ordered list of customers visited, as well as the quantities of each commodity to be delivered at each stop must be decided. We formalize this problem and include additional details in the remaining of Section 2.

Figure 1 (a) illustrates the routing, bin packing, and inventory aspects on a tiny example instance including one depot (in red with chimneys), two customers (in grey), two commodities (red and blue) and a three-days horizon. For simplicity, neither release nor demand occurs. Only initial inventories and deliveries are factored in the inventory dynamics. Inventory levels are proportional to the sizes of the bars on the left of the sites, and measured in the evening of each day. The space occupied by the commodities loaded in the vehicles on each transport section is proportional to the bars above each arrow. On the first day, two routes start from the depot (one represented by the dotted arrow, the other by the full arrow) and reach the first customer on the second day. On the third day, the second route (full arrow) reaches the second customer. No route starts on the third day. Inventory levels are updated in accordance with the quantities delivered. To illustrate the complexity of our instances, Figure 1 (b) shows the solution of a real European instance provided by our algorithm. Routes are drawn with black lines on the continent. The thicker the line, the more the route is used. We understand well the continuous-time aspect at this scale.
<figure><img src="image_2.png" /><figcaption>Figure 1: Instance and solution. Tiny example with two customers and two commodities (a), real European instance and routes solution in black (b).</figcaption></figure>

## 2.1 Notations and Data

Let $Z^+$ be the set of non-negative integers. For $a \in Z^+$, we denote by $[a]$ the set $\{1, ..., a\}$. Besides, for $x \in \mathbb{R}$, we define $(x)^+ := \max(x, 0)$. We denote by $|S|$ the cardinal of a set or list $S$. When we explicitly consider a vector $\mathbf{x} = (x_1, ..., x_p)$ of dimension $p \in Z^+$, we use the notation $\mathbf{x} \in Z$ or $\mathbf{x} \in \{0, 1\}$ instead of $\mathbf{x} \in Z^{|p|}$ or $\mathbf{x} \in \{0, 1\}^{|p|}$ respectively. Let $M$ be the set of commodities, $D$ the set of depots and $C$ the set of customers, that respectively release and demand commodities $m \in M$. The time horizon is $T \in Z^+$ days. At the beginning, each vertex $v$ (depot or customer) has an initial inventory of commodity $m$ denoted by $I_{mv}^0$. On each day $t \in [T]$, a customer $c$ demands a quantity $b_{mct}^{-}$ of commodity $m$. A depot $d$ releases a quantity $b_{mdt}^{+}$ of commodity $m$. We say that a depot $d \in D$ uses a commodity $m \in M$ if it has a positive initial inventory or a positive release for $m$ at least once over the horizon. We denote by $M_d$ the set of commodities used by depot $d$. We similarly define $M_c$ as the set of commodities used by customer $c \in C$, based on initial inventory and demand. A maximum inventory capacity $\kappa_{mvt}$ is set on the night of each day $t$ per vertex $v$ and commodity $m$. Below this capacity, no inventory cost is paid. Above, a cost is set to $c_{mv}^{exc}$ per unit, where "exc" stands for excess. Besides, a price $c_{mv}^{short}$ is paid per unit of unsatisfied demand for commodity $m$ of customer $c$, where "short" stands for shortage. It corresponds to a soft constraint of non-negativity for the customers' inventories. We approximate commodities and vehicles by one-dimensional objects. We associate a length $l_m$ to each commodity $m \in M$. We consider an infinite fleet of homogeneous vehicles of length $L$, to deliver the commodities from depots to customers. They are not assigned to a particular depot. A 1D bin packing problem must be solved for vehicle loading. The depots and customers are the vertices $V = D \cup C$ of a directed graph $D = (V, A)$ that we name the locations graph. The directed aspect is used to model the fact that transport durations and distances depend on the trip direction. There is an arc $\langle a, b \rangle \in A$ for each vertex $a \in D \cup C$ and $b \in C, b \neq a$. Given a vertex $v$, we denote by $\delta^{+}(v)$ the set of arcs outgoing from $v$, and by $\delta^{-}(v)$ the set of arcs incoming to $v$. We associate a distance $\Delta_a$ (in kilometers) and a transport duration $\tau_a$ (in hours) to each arc. We assume that the distances satisfy the triangular inequality. When planning a route, a cost is paid per vehicle $c^{veh}$, per stop (customer visited) $c^{stop}$, and per kilometer travelled $c^{km}$. The number of stops must not exceed $S_{max}$, which is a practical requirement of the car manufacturer. The limit of driving hours per day is $\tau_{max}$. The IRP consists in building a set of routes (see Section 2.2) to deliver commodities from depots to customers, minimizing the sum of the routing, inventory and shortage costs and respecting feasibility constraints detailed in Sections 2.3 and 2.2.

## 2.2 The Route Structure

An admissible path $P = (v_0, v_1, ..., v_k)$ in the locations graph is an elementary path, i.e., a path with pairwise distinct vertices. It starts from a depot $v_0 \in D$, and visits customers $(v_1, ..., v_k) \in C^k, v_i \neq v_j$. We emphasize it means routes are open, they do not end in their starting depots. We have a limit $S_{max}$ to the number of
customers visited:

$$|P| \leq S_{max} + 1. \qquad (1)$$

This constraint is required for practical reasons by the car manufacturer, because of the complexity of the unloading process. One alternative could be to limit the total duration or the total distance travelled by a route. We highlight the fact that a path does not end at its starting depot. Let $\mathcal{P}$ be the set of admissible paths, and $A(P)$ the set of arcs in a path $P$. A route $r$ is a “timed and loaded path”. It is a tuple $r = (t^r, P^r, q^r)$ where:

* $t^r \in [T]$ is the day of the departure.

* $P^r = (v_0^r, v_1^r, ..., v_k^r) \in \mathcal{P}$ is the admissible path followed.

* $q^r = (q_{ms}^r)_{m \in M, s \in [|P^r|-1]} \in (\mathbb{Z}^+)^{|M| \times (|P^r|-1)}$ are the quantities delivered, for each commodity $m \in M$ and to each customer $v_s^r$ for $s \in [|P^r|-1]$.

The total load must not exceed the vehicle capacity $L$, which can be written as:

$$\sum_{m \in M} \ell_m \left( \sum_{s \in [|P^r|-1]} q_{ms}^r \right) \leq L. \qquad (2)$$

Given a route $r$, and the transport durations $\tau_a$ for $a \in \mathcal{A}$, we can compute the arrival day $t_s^r$ at the customer $v_s^r$ for $s \in [|P^r|-1]$ as follows. We first compute the cumulated transport duration in hours up to customer $v_s^r$, with $\tau_0^r = 0$:

$$\tau_s^r = \tau_{s-1}^r + \tau_{(v_{s-1}^r, v_s^r)}, \quad \forall s \in [|P^r|-1]. \qquad (3)$$

Then, the actual day $t_s^r$ of arrival at customer $v_s^r$ takes pauses into account:

$$t_s^r = t^r + \left\lfloor \frac{\tau_s^r}{\tau_{\max}} \right\rfloor, \quad \forall s \in [|P^r|-1]. \qquad (4)$$

Equation (4) means that when a driver exceeds the driving time limit per day $\tau_{\max}$, a pause is made until the next day. The vehicle then goes on from the location of the pause. Since in practice routes start from depots in the morning and the deliveries are only available in the evening at the customers, it is indeed a floor and not a ceiling function we consider in Equation (4). A route must also visit every stop before the horizon $T$, which can be written as:

$$t_s^r \leq T, \quad \forall s \in [|P^r|-1]. \qquad (5)$$

We henceforth denote by $\mathcal{R}$ the set of admissible routes. A direct route follows a path $P = (d, c)$ from a depot $d \in D$ to a customer $c \in C$ in the locations graph. It has only one arc.

## 2.3 Inventory Routing Formulation

The variables we consider are the following. We denote by $z_{mdt}^{-}$ the quantity of commodity $m$ sent from depot $d$ on day $t$, and by $z_{mct}^{+}$ the quantity of commodity $m$ delivered to customer $c$ on day $t$. Let $I_{mvt}$ be the inventory of commodity $m$ at vertex $v$ on the evening of day $t$. We last denote by $x_r$ the number of vehicles
following route $r$. We then consider the MILP formulation:

$$ \min_{\mathbf{x},\mathbf{z},\mathbf{I}} \quad \sum_{r} x_r \left( c^{\text{veh}} + c^{\text{stop}}(|P^r| - 1) + c^{\text{km}} \sum_{a \in A(P^r)} \Delta_a \right) \quad (\text{multi-attribute-IRP}) \\ + \sum_{d,t,m} c_{md}^{\text{exc}} (I_{mdt} - \kappa_{mdt})^{+} \qquad (6a) $$

$$ + \sum_{c,t,m} c_{mc}^{\text{exc}} (I_{mct} - \kappa_{mct})^{+} + c_{mc}^{\text{short}} (b_{mct}^{-} - I_{mc(t-1)})^{+} \qquad (6b) $$

$$ \text{subject to} \quad z_{mdt}^{-} = \sum_{r, \substack{v_0^r=d, \\ t^r=t}} \sum_{s \in [|P^r|-1]} x_r q_{ms}^r, \quad \forall m \in M, \quad \forall d \in D, \quad \forall t \in [T] \qquad (6c) $$

$$ z_{mct}^{+} = \sum_{r} \sum_{s \in [|P^r|-1], t_r^s=t, v_r^s=c} x_r q_{ms}^r, \quad \forall m \in M, \quad \forall c \in C, \quad \forall t \in [T] \qquad (6d) $$

$$ I_{mdt} = I_{md(t-1)} + b_{mdt}^{+} - z_{mdt}^{-}, \quad \forall m \in M, \quad \forall d \in D, \quad \forall t \in [T] \qquad (6e) $$

$$ I_{md0} = I_{md}^{0}, \quad \forall m \in M, \quad \forall d \in D \qquad (6f) $$

$$ I_{mct} = (I_{mc(t-1)} - b_{mct}^{-})^{+} + z_{mct}^{+}, \quad \forall m \in M, \quad \forall c \in C, \quad \forall t \in [T] \qquad (6g) $$

$$ I_{mc0} = I_{mc}^{0}, \quad \forall m \in M, \quad \forall c \in C \qquad (6h) $$

$$ \mathbf{x} \ge 0, \mathbf{z} \ge 0, \mathbf{I} \ge 0 \qquad (6i) $$

$$ \mathbf{x} \in \mathbb{Z}, \mathbf{z} \in \mathbb{Z}, \mathbf{I} \in \mathbb{Z} \qquad (6j) $$

We notice that, given an IRP instance and the route variables $\mathbf{x}$, we can deduce the quantities sent or received $\mathbf{z}$ and the inventory $\mathbf{I}$. The latter are useful to express the IRP as an MILP.

**Objective function.** The first sum models a cost per vehicle, per stop and per kilometer travelled. The second one is related to the excess inventory during the nights at the depots. The third one has both an excess inventory and a shortage part. The quantity $(b_{mct}^{-} - I_{mc(t-1)})^{+}$ is a substitute that is bought separately when a shortage appears.

**Constraints.** The $\mathbf{x}$ variable is used to count the number of vehicles that follow the admissible routes defined in Section 2.2. Equation (6c) is used to bind the total quantities that are sent from each depot to the route deliveries at each customer. Equation (6d) links the total quantities received per customer to the route deliveries. Constraints (6e)-(6f) define the inventory dynamics at the depots, and (6g)-(6h) at the customers. We highlight we cannot deliver a commodity to a customer that does not need it – in the sense of $C_m$ – because the maximum inventory capacity is set to zero and the excess inventory cost to infinity. The MILP (multi-attribute-IRP) is intractable over our instances, we instead suggest several heuristic approaches in the next section.

# 3 Overview of the Algorithms and General Concepts

We emphasize the main principles of the algorithms in Section 3.1 before going into the details of each of their components. We also introduce generic flow graphs and formulations in Section 3.2, concepts useful in the rest of the present paper.

## 3.1 Overview of the Algorithms

We compare three algorithms with increasing degrees of sophistication and performance: an initialization + local search algorithm to quickly derive non-trivial IRP solutions, a route-based matheuristic, and our LNS. The two first algorithms are adapted from frameworks of the literature, the last one is our main contribution.
### 3.1.1 Subroutines

Our algorithms are illustrated on Figure 2 and combine five subroutines. We call inner iteration an iteration within any subroutine, and outer iteration a path through the four types of neighborhoods in the LNS (see the loop in Figure 2). The first subroutine builds an initial solution. 1) The flow relaxation + bin packing (flow relaxation, bin packing) subroutine solves a flow relaxation – thus an LP – per commodity and deduces direct routes by approximately solving bin packing problems to respect vehicle capacity. The other four subroutines improve or perturb an existing solution, and can be applied any number of times in any order. 2) The routing local search (routing local search) subroutine takes a random subset of routes and applies a local search with TSP and SDVRP neighborhoods (see Section 5.1). 3) The reload fixed-path vehicles (reload fixed path vehicles) subroutine solves an MILP per depot to re-optimize the load of the routes starting from it. The MILP is solved with a very low gap threshold in the route-based matheuristic, and up to a larger gap threshold in the LNS, adding a time limit. 4) The customer reinsertion (customer reinsertion) subroutine removes a customer from every delivery of a solution and solves an MILP to reinsert it in the existing routes, also creating new direct routes. 5) The commodity reinsertion (commodity reinsertion) subroutine removes a commodity from every delivery of a solution and solves an MILP to reinsert it. These MILPs are solved up to a gap and time limit. We see the last four subroutines as local search procedures, and call them with a customizable number of inner iterations per outer iteration. In Section 4, we detail the flow relaxation + bin packing subroutine. The routing local search subroutine is described in Section 5. The three large neighborhoods – reload fixed-path vehicles, customer and commodity reinsertion – are detailed in Sections 6.1, 6.2 and 6.3 respectively.

<figure><img src="image_3.png" /><figcaption>Figure 2: Different algorithms ordered by degree of sophistication. Dark gray corresponds to initialization, gray to descent subroutines, and light gray to perturbation subroutines.</figcaption></figure>

### 3.1.2 Algorithms

**Initialization + local search.** This algorithm simply runs the initialization + bin packing subroutine to build an initial solution, and then applies the routing local search subroutine to improve it. It has the advantage of being fast (about four minutes on our large-scale instances on average) since it is based on an LP. It is detailed in Algorithm 1.
**Algorithm 1:** Initialization + local search

input : $\mathcal{I}$ an IRP instance.

output: A solution $\mathbf{r}$ to the IRP.

$\mathbf{y} = \text{flow relaxation}(\mathcal{I})$;

$\mathbf{r} = \text{bin packing}(\mathcal{I}, \mathbf{y})$;

$\mathbf{r} = \text{routing local search}(\mathcal{I}, \mathbf{r})$;

**Route-based matheuristic.** As discussed in Section 1, Archetti, Boland, and Speranza (2017) and Bertazzi et al. (2019) solve the IRP given a subset of “promising routes” that are defined heuristically. This corresponds to reducing the set of feasible solutions to the IRP, which allows solving an MILP. The routes are either selected among those created during a tabu search and leading to cost improvements, or in a constructive manner. The underlying assumption is that they are likely to appear in a good IRP solution. We adapt this idea to our setting. The main difficulty is that our MILP (multi-attribute-IRP) is intractable, even when we restrict the set of “promising routes” to the set of an initial solution. This is due to the multicommodity aspect, the scale of our instances and the routes that last several days. We solve the restricted MILP heuristically, using a large neighborhood approach: 1) Apply the initialization + local search algorithm to get an initial solution. 2) Take the current solution as set of promising routes, and solve sequentially one MILP per depot, with the corresponding promising routes that start from it. It is detailed in Algorithm 2.

**Algorithm 2:** Route-based matheuristic

input : $\mathcal{I}$ an IRP instance, $\mathbf{r} = (r_k)_{1 \le k \le K}$ the current solution with $K \in \mathbb{Z}^+$ routes, time limit a time limit.

output: The solution $\mathbf{r}$ updated.

$\mathbf{r} = \text{initialization} + \text{local search}(\mathcal{I})$;

for $d \in D$ do

$\left| \begin{array}{l} \mathbf{r} = \text{reload fixed path vehicles}(\mathcal{I}, \mathbf{r}, d); \\ \text{if } \text{time elapsed} \ge \text{time limit} \text{ then} \\ \quad \lfloor \text{break;} \end{array} \right.$ // small gap

**Large neighborhood search.** The large neighborhood search Algorithm 3 first uses the initialization + local search approach to find a good initial solution. It then explores four kinds of neighborhoods. Two of them always improve the solution: the routing local search and reload fixed-path vehicles subroutines. Contrary to the route-based matheuristic, the latter is applied with a greater gap limit and an additional time limit, in order to avoid spending too much time within it and cycle over the neighborhoods instead. The two remaining ones are perturbations, which means they can deteriorate the solution. They fix a part of the current solution and optimize the quantities and routes involving a particular customer or commodity over the entire horizon. They both lead to substantial changes, allowing the search to escape from local minima. The LNS uses both iteratively, selecting every customer and ten commodities at random per outer step. Those hyperparameters are tuned experimentally as shown in Appendix E. The outer LNS iterations are illustrated by the loop arc on Figure 2. The LNS returns the best solution found, comparing after each subroutine the current solution with the best one so far. The main idea behind this LNS is to consider the structure of the IRP, “decompose” it along its major axes, and solve smaller natural problems to explore the solution space. One difference with Coelho, De Maio, and Laganà (2020) or Archetti, Guastaroba, et al. (2021), is that our idea is not to fix a part of the MILP variables (multi-attribute-IRP) and optimize with respect to the remaining ones, but to define new smaller MILPs based on the structure of the IRP. Let us now introduce a concept that helps describing our subroutines.
**Algorithm 3:** Large neighborhood search

**input**: $\mathcal{I}$ an IRP instance, $\mathbf{r} = (r_k)_{1\le k\le K}$ the current solution with $K \in \mathbb{Z}^+$ routes, **time limit** a time limit, $n_{\text{cust}}$ the number of customers to reinsert, $n_{\text{comm}}$ the number of commodities to reinsert.

**output**: The solution $\mathbf{r}$ updated.

$\mathbf{r} = \text{initialization} + \text{local search}(\mathcal{I})$;

while $\text{time elapsed} < \text{time limit}$ do

$\mathbf{r} = \text{routing local search}(\mathcal{I}, \mathbf{r});$
for $d \in D$ do
  $\lfloor \mathbf{r} = \text{reload fixed path vehicles}(\mathcal{I}, \mathbf{r}, d); \rfloor$ // large gap
  $C_{\text{sub}} = \text{sample}(C, n_{\text{cust}})$;
  for $c \in C_{\text{sub}}$ do
    $\lfloor \mathbf{r} = \text{customer reinsertion}(\mathcal{I}, \mathbf{r}, c); \rfloor$ // large gap
    $M_{\text{sub}} = \text{sample}(M, n_{\text{comm}})$;
    for $m \in M_{\text{sub}}$ do
      $\lfloor \mathbf{r} = \text{commodity reinsertion}(\mathcal{I}, \mathbf{r}, m); \rfloor$ // large gap

## 3.2 Flow Graphs and Formulations

Let us consider the generic MILP formulation:

$$ \min_{\mathbf{y},\mathbf{x}} \sum_{m \in M} \mathbf{y}_m^\top \mathbf{c}_m + \sum_{r \in \mathbf{r}} x_r c_r \qquad (\text{generic-flow-MILP}) $$

$$ \text{subject to} \quad \sum_{a \in \delta^+(v)} y_{ma} = \sum_{a \in \delta^-(v)} y_{ma}, \quad \forall m \in M, \forall v \in \mathcal{V}^m \qquad (7a) $$

$$ \mathbf{y}_m^{\min} \le \mathbf{y}_m \le \mathbf{y}_m^{\max}, \quad \forall m \in M \qquad (7b) $$

$$ \sum_{m \in M} y_{ma} l_m \le x_r L, \quad \forall a = (d \to (c, r)), \forall r \in \mathbf{r} \qquad (7c) $$

$$ \mathbf{y} \in \mathbb{Z} \qquad (7d) $$

$$ \mathbf{x} \in \{0, 1\}. \qquad (7e) $$

The variable $y_m$ encodes a flow on a given commodity graph thanks to Equations (7a)-(7b). This flow enables modelling the depot and customer inventory dynamics of commodity $m$ defined by constraints (6e)-(6h), and the quantities sent and received. Variable $\mathbf{x}$ encodes the routes that are used to deliver the commodities. Constraint (7c) indeed enforces the flows to respect the vehicle capacity when a route is used. The generic notation $(d \to (c, r))$ refers to the transport arc from the depot $d$ to the customer $c$ in the route $r$, on its departure day. In order to obtain a specific MILP formulation from this generic one, we must specify which commodity graph is used, and which set of routes $\mathbf{r}$ is considered. We are going to use several distinct commodity graphs, see for instance the figures 4, 5, 6 for detailed $(d \to (c, r))$ arcs. However, they all share a common structure which we describe now:

**Depot subgraphs.** A subgraph per depot $d$ (Figure 3 (b)), which models its inventory dynamics. It is shared by the distinct formulations we introduce in this paper and has the following vertices: (t, d, morning) for $t \in [T+1]$, and (t, d, evening) for $t \in [T]$.

**Customer subgraphs.** A subgraph per customer $c$ (Figure 3 (c)), which models its inventory dynamics. It is also shared by the distinct formulations we introduce in this paper, and it has the following vertices: (t, c, morning) for $t \in [T+1]$, and (t, c, evening) for $t \in [T]$.
<table><caption>Table 1: Arcs of the commodity flow graph shared by our formulations. When not stated in the table, Min is 0, Cost is 0 and Max is $\infty$.</caption><thead><tr><th>Subgraph</th><th>Arc type</th><th>Arc description</th><th>Origin</th><th>Destination</th><th>Min</th><th>Max</th><th>Cost</th></tr></thead><tbody><tr><td>Depot</td><td>Incoming</td><td>Initial inventory depot</td><td>initial inventory</td><td>(1, d, morning)</td><td>I<sup>0</sup><sub>md</sub></td><td>I<sup>0</sup><sub>md</sub></td><td></td></tr><tr><td>Depot</td><td>Outgoing</td><td>Final inventory depot</td><td>final inventory</td><td>(T + 1, d, morning)</td><td></td><td></td><td></td></tr><tr><td>Depot</td><td>Incoming</td><td>Release depot</td><td>release</td><td>(t, d, morning)</td><td>b<sup>+</sup><sub>mdt</sub></td><td>b<sup>+</sup><sub>mdt</sub></td><td></td></tr><tr><td>Customer</td><td>Incoming</td><td>Initial inventory customer</td><td>initial inventory</td><td>(1, c, morning)</td><td>I<sup>0</sup><sub>mc</sub></td><td>I<sup>0</sup><sub>mc</sub></td><td></td></tr><tr><td>Customer</td><td>Outgoing</td><td>Final inventory customer</td><td>final inventory</td><td>(T + 1, c, morning)</td><td></td><td></td><td></td></tr><tr><td>Customer</td><td>Outgoing</td><td>Demand customer</td><td>demand</td><td>(t, c, morning)</td><td>b<sup>-</sup><sub>mct</sub></td><td>b<sup>-</sup><sub>mct</sub></td><td></td></tr><tr><td>Customer</td><td>Incoming</td><td>Shortage customer</td><td>shortage</td><td>(t, c, morning)</td><td></td><td></td><td></td></tr><tr><td>Routes</td><td></td><td>Transport d → c</td><td>Formulation specific</td><td>(see Sections 4.1, 6.1, 6.2, 6.3)</td><td></td><td></td><td></td></tr><tr><td>Depot</td><td>Internal</td><td>Daily inventory depot</td><td>(t, d, morning)</td><td>(t, d, evening)</td><td></td><td></td><td></td></tr><tr><td>Depot</td><td>Internal</td><td>Free night inventory depot</td><td>(t, d, evening)</td><td>(t + 1, d, morning)</td><td></td><td></td><td></td></tr><tr><td>Depot</td><td>Internal</td><td>Excess night inventory depot</td><td>(t, d, evening)</td><td>(t + 1, d, morning)</td><td></td><td></td><td></td></tr><tr><td>Customer</td><td>Internal</td><td>Daily inventory customer</td><td>(t, c, morning)</td><td>(t, c, evening)</td><td></td><td></td><td></td></tr><tr><td>Customer</td><td>Internal</td><td>Free night inventory customer</td><td>(t, c, evening)</td><td>(t + 1, c, morning)</td><td></td><td></td><td></td></tr><tr><td>Customer</td><td>Internal</td><td>Excess night inventory customer</td><td>(t, c, evening)</td><td>(t + 1, c, morning)</td><td></td><td></td><td></td></tr><tr><td>Artificial</td><td></td><td>Circulation</td><td>source</td><td>release</td><td></td><td></td><td></td></tr><tr><td>Artificial</td><td></td><td>Circulation</td><td>source</td><td>initial inventory</td><td></td><td></td><td></td></tr><tr><td>Artificial</td><td></td><td>Circulation</td><td>source</td><td>shortage</td><td></td><td></td><td></td></tr><tr><td>Artificial</td><td></td><td>Circulation</td><td>demand</td><td>sink</td><td></td><td></td><td></td></tr><tr><td>Artificial</td><td></td><td>Circulation</td><td>final inventory</td><td>sink</td><td></td><td></td><td></td></tr><tr><td>Artificial</td><td></td><td>Circulation</td><td>sink</td><td>source</td><td></td><td></td><td></td></tr></tbody></table>

**A route subgraph.** Its specific structure depends on the formulation we consider. It contains paths between vertices of the form $(t, d, \text{morning})$ and vertices of the form $(\tilde{t}, c, \text{evening})$ as shown on Figure 3 (a). Flow variables $y$ on those paths model the quantities sent from depots to customers.

**Artificial vertices.** In order to model commodity flows as circulations over commodity graphs, we add artificial vertices connected to the subgraphs above: `source`, `sink`, `initial inventory`, `final inventory`, `release`, `shortage`, and `demand`.

The details of the arcs of the shared subgraphs defined above are in Table 1. For each arc, we give the following information. The subgraph it belongs to is first given. Then, we distinguish “incoming”, “outgoing” and “internal” arcs with respect to the depots and customer subgraphs. A short description is stated to understand the meaning of the arcs. We specify the origin and destination vertices, as well as the minimum and maximum flow capacities associated to the flow variables on the arcs. Last, the cost corresponding to these variables are also given.

Given a graph $\tilde{D} = (\tilde{V}, \tilde{A})$ with capacities associated to its arcs $(y^{\min}, y^{\max})$, we define the set of circulations as $C(\tilde{D}, y^{\min}, y^{\max}) = \{y \in \mathbb{R}^A, y_{\min} \le y \le y_{\max}, \forall v \in \tilde{V}, \sum_{a \in \tilde{\delta}^{+}(v)} y_a = \sum_{a \in \tilde{\delta}^{-}(v)} y_a\}$. We use this notation instead of constraints (7a)-(7b) in the rest of the paper. The minimum $y^{\min}$ and maximum $y^{\max}$ capacities are defined by the `Min` and `Max` columns in tables 1, 4 and figures 4 (b), 5 (b) and 6 (b).

The (generic-flow-MILP) formulation is not fully specified here as we still need to define the route subgraph to obtain a formulation that can be given to a solver. It is an exact formulation if we enumerate all possible routes, but such a column generation is intractable on our large instances. Therefore, we derive heuristic formulations with specific route subgraphs in Sections 4.1, 6.1, 6.2, and 6.3.

**Remark 3.1** In (generic-flow-MILP), routes are modelled with individual paths in the commodity graphs, bound with indicator variables $x$. Sometimes, we define them as paths in another flow graph. In this case, some vertices and arcs are shared between routes. We detail this aspect in Sections 6.2 and 6.3.

**Remark 3.2** We sparsify the commodity graph. Instead of considering the sets of depots $D$ and customers $C$ in the commodity graph $D^m$, we define the subsets $D_m = \{d \in D, m \in M_d\}$ and $C_m = \{c \in C, m \in M_c\}$ the
<figure><img src="image_4.png" /></figure>

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

<figure><img src="image_6.png" /><figcaption>Figure 3: Commodity flow graph. Overview of the global graph structure (a), and then details of the subgraphs “Depot 1” (b) and “Customer 1” (c). The abbreviations “mo” and “ev” stand for morning and evening respectively.</figcaption></figure>

depots and customers that use commodity m. We then restrict the depots and customer subgraphs of $\mathcal{D}^m$ to the ones of $D_m$ and $C_m$.

# 4 Flow Relaxation + Bin Packing Subroutine

The flow relaxation + bin packing subroutine is a fast heuristic to get an initial solution to (multi-attribute-IRP). It takes as input an IRP instance, and returns an initial IRP solution built from intermediate flow solutions that encode who sends what to whom and when.

## 4.1 Multiple Minimum Cost Flows and Relaxation

**Minimum cost flow formulation.** Let $y = (y_{ma})_{m \in M, a \in A^m}$ be a flow variable and $c$ be the corresponding costs defined in Table 1 for the shared subgraphs and Figure 4 for the specific route subgraph. We consider the following LP:

$$ \begin{array}{ll} \min_{y} & \displaystyle\sum_{m \in M} y_m^\top c_m \qquad (\text{flow-relaxation}) \\ \text{subject to} & y_m \in \mathcal{C}(\mathcal{D}^m, y_m^{\min}, y_m^{\max}), \quad \forall m \in M \qquad (8a) \\ & y_m \ge 0 \quad \forall m \in M \end{array} $$

In this variant of the generic MILP (Section 3.2), we do not introduce a route variable $x$ and the corresponding cost. We instead consider a soft version of constraint (7c) in the commodity cost $c_m$. This leads to a separate flow LP for each commodity in (flow-relaxation). Let us now introduce the details of the commodity graph.
**Details of the commodity graph.** We define one graph per commodity *m*, named $\mathcal{D}^m = (\mathcal{V}^m, \mathcal{A}^m)$. The vertices $\mathcal{V}^m$ are exactly the ones defined in Section 3.2. The arcs detailed in Table 1 are included. We specify the route subgraph in the table of Figure 4. It is made of direct transport and delayed transport $d \to c$ arcs. These arcs are added when the date of arrival at the customer is smaller than the horizon $T$. For each tuple $(t, d, c) \in [T] \times D_m \times C_m$, we add one delayed arc $((t, d) \to (\tilde{t}, c))$ per possible delayed arrival day $\tilde{t}$ induced by an indirect path from $d$ to $c$ (thus visiting any set of other customers before $c$) that respects the route constraints defined by $\mathcal{R}$. Those possible delays are pre-computed, using a breadth-first search algorithm over the locations graph, with maximum depth set to $S_{max}$. Indeed, we can browse the locations graph starting from depots, saving the cumulative delay at any vertex and any depth smaller than $S_{max}$. The flow on these arcs models the quantity of commodity sent from depot $d$ on day $t$ to customer $c$ with arrival on day $\tilde{t}$. The precise structure of the graph is illustrated on Figure 4. In arc annotations, capacities are given between brackets (e.g. $[0, \kappa_{mdt}]$), and costs without (e.g. $c_{md}^{exc}$). Dotted arrows are related to the shared artificial vertices, continuous ones to depots, customers and route subgraphs. When not stated in the table, **Min** is 0, **Cost** is 0 and **Max** is $\infty$. On this figure, only two days, one depot and one customer are shown. Besides, some artificial vertices are omitted for simplicity. The cost $c_{mdc}^{tr}$ is detailed below.

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

<table><caption>(b) Additional commodity graph arcs compared with Table 1</caption><thead><tr><th>Subgraph</th><th>Arc description</th><th>Origin</th><th>Destination</th><th>Min</th><th>Max</th><th>Cost</th></tr></thead><tbody><tr><td>Routes</td><td>Transport $d \to c$</td><td>(t, d, morning)</td><td>(t + <sup>$\frac{rdc}{T_{max}}$</sup>&#x20D7;, c, evening)</td><td></td><td></td><td>$c_{mdc}^{tr}$</td></tr><tr><td>Routes</td><td>Delayed transport $d \to c$</td><td>(t, d, morning)</td><td>(\tilde{t}, c, evening)</td><td></td><td></td><td>$c_{mdc}^{tr}$</td></tr></tbody></table>

Figure 4: Details of the commodity graph $\mathcal{D}^m$ for the flow relaxation problem.

Since the routing price is paid at the vehicle level, we cannot derive a minimum cost commodity flow that takes it into account exactly without adding variables for each individual vehicle. Instead, we want to approximate this cost with transportation arcs between depots and customers naturally involving commodity flow variables. A way to do so is to use a “vehicle fraction” unit per commodity, leading to:

$$c_{mdc}^{tr} = \frac{\ell_m}{L} (c^{\text{veh}} + c^{\text{stop}} + c^{\text{km}} \Delta_{dc}), \quad \forall d \in D_m, \forall c \in C_m. \qquad (9)$$

In Equation (9) the factor $\ell_m/L$ is a way to scale the price paid for the delivery of a unit of commodity $m$
based on the percentage of a vehicle it occupies, hence the “vehicle fraction”.

**Proposition 4.1** The optimization problem (flow-relaxation) based on $|M|$ flows is a relaxation of (multi-attribute-IRP). The optimal value of (flow-relaxation) is a lower bound to the cost of an optimal solution to our initial problem.

We now sketch the proof. Given a feasible solution of (multi-attribute-IRP), we can deduce a feasible solution of (flow-relaxation) by fixing the quantities sent by each depot to each customer per day, delay and commodity. The inventory costs are modelled exactly with (flow-relaxation) thanks to the delayed arcs, thus equal to the ones of (multi-attribute-IRP). The transportation costs are lower bounded with Equation (9). The structure of the graph detailed on Figure 4 only allows geographically direct routes between depots and customers, whereas the route constraints allow up to $S_{max}$ stops. Nonetheless, considering additional arcs $c_1 \rightarrow c_2$ with fraction costs $c_{mc_1c_2}^{tr} = \frac{\ell_{m}}{L}(c^{stop} + c^{km}\Delta_{c_1c_2})$ and delay $\tau_{c_1c_2}$ in an extended flow graph would also produce solutions with only direct routes. Indeed, by triangular inequality, it would always be cheaper to send commodities through geographically direct $(d \rightarrow c)$ (possibly delayed) arcs in this framework of “vehicle fraction” costs, rather than sending quantities to intermediate customer $c_1$ before reaching the destination $c_2$.

## 4.2 Bin Packing

We highlight here how the minimum cost flows can be used to derive an IRP solution, a step further from the lower bound computation. The $|M|$ minimum cost flow solutions resulting from (flow-relaxation) enable us to set the quantities sent by each depot to each customer per day and commodity, but do not directly lead to a set of routes. Indeed, for now, we do not know how quantities are loaded in various vehicles. We highlight the fact that the delayed arcs are not used to build an initial solution. They are only introduced to compute a lower bound. To deduce a set of direct routes, we approximately solve one bin packing problem per tuple $(d, c, t) \in D \times C \times [T]$, using the first-fit-decreasing heuristic. The instance of the bin packing is given by the set of commodities to be sent on day $t$ from depot $d$ to customer $c$, their respective lengths, and the length of one vehicle. The solution to the bin packing problem leads to a low number of vehicles each of length $L$, with corresponding loading made of possibly $|M|$ distinct commodities. At this point, we get a set of direct routes as a first feasible solution to (multi-attribute-IRP).

# 5 Routing Local Search Subroutine

Our solution processes emphasized on Figure 2 rely on the routing local search subroutine of Algorithm 4. We detail here both the local search procedure, and the neighborhoods listed in Table 2.

## 5.1 Neighborhoods

Before introducing the local search procedure, we focus on the TSP and SDVRP neighborhoods in an IRP framework. The routing local search subroutine combines those neighborhoods.

**Routes impact inventories.** We highlight the fact that the 12 neighborhoods detailed below in Table 2 alter the routes of a solution and the inventories at the depots or at the customers (contrary to the SDVRP framework). Therefore, whenever a neighborhood is considered, we evaluate the effects on inventories and routes so as to check feasibility and to estimate the cost change. For instance, the optimal order of a route not only depends on the distances $\Delta$, but also on the delays introduced in the inventory dynamics of the customers involved. Therefore, even elementary neighborhoods require calculation. They can be seen as a generalization of the TSP and SDVRP concepts to the continuous-time IRP.
<table><caption>Table 2: Routing local neighborhoods: single-depot and multi-depot variants are considered.</caption><thead><tr><th>Type</th><th>Name</th><th>Description</th></tr></thead><tbody><tr><td rowspan="3">TSP</td><td>relocate</td><td>change the position of one stop in a route</td></tr><tr><td>swap</td><td>exchange the positions of two stops in a route</td></tr><tr><td>2-opt*</td><td>cut a route into three parts and revert the order of the middle one</td></tr><tr><td rowspan="6">SDVRP single-depot</td><td>insert</td><td>give a stop s from route r<sub>1</sub> to route r<sub>2</sub></td></tr><tr><td>swap single depot</td><td>swap adapted to two routes with same depot</td></tr><tr><td>merge</td><td>merge two routes r<sub>1</sub> and r<sub>2</sub> on the same day</td></tr><tr><td>merge multi day</td><td>merge extended to routes with different start dates</td></tr><tr><td>delete route</td><td>delete a route</td></tr><tr><td>change day</td><td>move a route in time, one day before or after</td></tr><tr><td rowspan="3">SDVRP multi-depot</td><td>insert multi depot</td><td>insert extended to routes with distinct depots</td></tr><tr><td>swap multi depot</td><td>swap extended to routes with distinct depots</td></tr><tr><td>2-opt* multi depot</td><td>cut two routes each into two parts and exchange their end parts</td></tr></tbody></table>

## 5.2 Routing Local Search

From the list of neighborhoods emphasized in Table 2, we design the routing local search, Algorithm 4. Usually, local search procedures start from the smallest neighborhoods and use larger ones to escape from local minima. Our case is different because of the size and complexity of the problem, and the short time given to solve it. The delete route, change day and TSP neighborhoods (that we name single-route neighborhoods) can be explored exhaustively, until no improvement is found. The SDVRP neighborhoods involve a lot of computations and cannot be fully browsed. Besides, we want each route in the solution after the routing local search procedure to be optimized from a TSP viewpoint. Therefore, Algorithm 4 first browses randomly the large SDVRP neighborhoods to create new routes, and then refines thoroughly with the single route neighborhoods. Since the neighborhoods tend to reduce the number of routes, the routing local search algorithm deletes routes at several stages. It only applies feasible moves that improve the cost.

As said above, some features of this local search differ from the common SDVRP local search algorithms. Because of the computations involved for insert, swap single depot, and for the multi-depot neighborhoods, and the number of pairs of routes, the routing local search subroutine only samples a fraction of them. To do so, it explicitly samples a subset of the pairs of routes of the solution according to a uniform distribution on each day. We tune this approach with parameter *p* in Algorithm 4 to find a good cost gain per CPU time ratio. We do so instead of restricting the routes candidates with geographic criteria, because the inventory costs cannot be neglected. Two routes that visit customers that are far from each other may still be suitable candidates for a swap for instance, due to the change in inventory cost. The change day function is applied per route one day forward or backward, until no improvement is found.

<table><caption>**Algorithm 4:** Routing local search</caption><tr><td>input</td><td>: $\mathcal{I}$ an IRP instance, $\mathbf{r} = (r_k)_{1 \le k \le K}$ the current solution with $K \in \mathbb{Z}^+$ routes, $n_{it} \in \mathbb{Z}^+$ a number of iterations, $p$ a percentage.</td></tr><tr><td>output</td><td>: The solution $\mathbf{r}$ updated.</td></tr><tr><td>for</td><td>$i = 1 : n_{it}$ do</td></tr><tr><td></td><td>SDVRP multi-depot neighborhoods over $p\%$ of the pairs of routes at random;<br>delete route per day and depot until no improvement;</td></tr><tr><td></td><td>SDVRP single-depot neighborhoods over $p\%$ of the pairs of routes at random;<br>delete route per day and depot until no improvement;</td></tr><tr><td></td><td>SDVRP single-depot neighborhoods over $p\%$ of the pairs of routes at random;<br>single-route neighborhoods until no improvement;</td></tr></table>
# 6 MILP-Based Neighborhoods and Perturbations

We now introduce three subroutines that are based on optimization problems written as MILPs. They all leverage the commodity graph structure emphasized in Section 3.2.

## 6.1 Reload Fixed-Path Vehicles Subroutine

**Reload neighborhood problem.** Let us define the problem behind this large neighborhood. We consider a subset of routes $\mathbf{r}_{reload}$ of the current IRP solution $\mathbf{r}$, in our case the routes that start from a given depot $d$. We solve the following problem: choose the routes to keep in the solution among $\mathbf{r}_{reload}$, and re-estimate the delivered quantities (for the whole set of commodities) of the routes kept, to minimize the total cost. In our large neighborhood setting, we fix the remaining routes of the current solution $\mathbf{r}$.

We denote by $x_r$ for $r \in \mathbf{r}_{reload}$ the indicator variable for keeping route $r$, and by $\mathbf{y} = (\mathbf{y}_m)_{m \in M}$ the set of commodity flow variables. The commodity flow graphs $(\mathcal{D}^m(\mathbf{r}, \mathbf{r}_{reload}))_{m \in M}$ involved are defined below. We model the problem with the following formulation:

$$ \min_{\mathbf{y},\mathbf{x}} \sum_{m \in M} \mathbf{y}_m^\top \mathbf{c}_m + \sum_{r \in \mathbf{r}_{reload}} x_r \left( c^{\text{veh}} + c^{\text{stop}}(|P^r| - 1) + c^{\text{km}} \sum_{a \in A(P^r)} \Delta_a \right) \quad (\text{Reload-MILP}) $$

$$ \text{subject to} \quad \mathbf{y}_m \in \mathcal{C}(\mathcal{D}^m(\mathbf{r}, \mathbf{r}_{reload}), \mathbf{y}_m^{\min}(\mathbf{r}, \mathbf{r}_{reload}), \mathbf{y}_m^{\max}(\mathbf{r}, \mathbf{r}_{reload})), \quad \forall m \in M \quad (10a) $$

$$ \sum_{m \in M} y_{ma} l_m \le x_r L, \quad \forall a = (d \to (c, r)), \quad \forall r \in \mathbf{r}_{reload} \quad (10b) $$

$$ \mathbf{y} \in \mathbb{Z} \quad (10c) $$

$$ \mathbf{x} \in \{0, 1\} \quad (10d) $$

This formulation is very close to the generic MILP introduced in Section 3.2. The objective function is composed of one flow cost per commodity $m \in M$ (inventory and shortage costs), and of the routing cost of each route kept among $\mathbf{r}_{reload}$. Constraint (10b) ensures that the commodity flows from depots to customers only exist along routes that are kept, and that the capacity of the vehicles is respected. The last two constraints define integer and binary variables. We highlight that this MILP exactly formulates the reloading of a given subset of routes. Solving the problem (Reload-MILP) leads to a new feasible solution with lower cost. Indeed, we exactly model the IRP constraints and costs, and optimize with respect to a subset of variables with (Reload-MILP): the commodity flows involved in the routes $\mathbf{r}_{reload}$, and the use of those routes. The rest of the variables are fixed. The current solution is an admissible solution of this MILP, used for warm-start.

The commodity graph $\mathcal{D}^m(\mathbf{r}, \mathbf{r}_{reload})$ for $m \in M$ depends both on the current solution $\mathbf{r}$, and on the routes to potentially keep and reload $\mathbf{r}_{reload}$. As previously, the backbone structure is the same as in Section 3.2: one subgraph per customer, one per depot, one for the routes, and additional vertices to create circulations. The special route subgraph is detailed on Figure 5. In arc annotations, capacities are given between brackets (e.g. $[0, \kappa_{mdt}$]), and costs without (e.g. $c_{md}^{exc}$). Dotted arrows are related to the shared artificial vertices, continuous ones to depots, customers and route subgraphs. When not stated in the table, Min is 0, Cost is 0 and Max is $\infty$. We explicitly create individual route paths, with vertices of the form $(t_s^r, c_s, r)$. It models the fact that when using route $r$, commodities are delivered to customer $c_s$ on day $t_s^r$ at position $s$ of the route. The details of the arcs can be found in the table of Figure 5. Besides, since we fix the quantities sent by the routes in $\mathbf{r} \setminus \mathbf{r}_{reload}$, we need two additional vertices in the route subgraph that we name fixed deliveries sent, and fixed deliveries received. The former is connected to depot vertices in order to take other quantities sent into account. The second is connected to customer vertices to model other quantities received.

## 6.2 Customer Reinsertion Subroutine

As mentioned in the overview Section 3, our neighborhoods are based on a decomposition of the IRP along its main axes. Previous sections focus on the routes. Here, we design a perturbation based on the customers.
<figure><img src="image_8.png" /></figure>

<table><caption>(b) Additional commodity graph arcs compared with Table 1</caption><thead><tr><th>Subgraph</th><th>Arc description</th><th>Origin</th><th>Destination</th><th>Min</th><th>Max</th><th>Cost</th></tr></thead><tbody><tr><td>Routes</td><td>Fixed deliveries sent</td><td>(t, d, morning)</td><td>fixed deliveries sent</td><td>z<sup>-</sup><sub>mdt</sub></td><td></td><td></td></tr><tr><td>Routes</td><td>Fixed deliveries received</td><td>fixed deliveries received</td><td>(t, c, evening)</td><td>z<sup>-</sup><sub>mct</sub></td><td></td><td></td></tr><tr><td>Routes</td><td>Transport d → (c<sub>1</sub>, r)</td><td>(t<sup>r</sup>, d, morning)</td><td>(t<sup>r</sup><sub>1</sub>, c<sub>1</sub>, r)</td><td></td><td></td><td></td></tr><tr><td>Routes</td><td>Transport (c<sub>s</sub>, r) → (c<sub>s+1</sub>, r)</td><td>(t<sup>r</sup><sub>s</sub>, c<sub>s</sub>, r)</td><td>(t<sup>r</sup><sub>s+1</sub>, c<sub>s+1</sub>, r)</td><td></td><td></td><td></td></tr><tr><td>Routes</td><td>Transport (c<sub>s</sub>, r) → c</td><td>(t<sup>r</sup><sub>s</sub>, c, r)</td><td>(t<sup>r</sup><sub>s</sub>, c, evening)</td><td></td><td></td><td></td></tr><tr><td>Artificial</td><td>Circulation</td><td>source</td><td>fixed deliveries received</td><td></td><td></td><td></td></tr><tr><td>Artificial</td><td>Circulation</td><td>fixed deliveries sent</td><td>sink</td><td></td><td></td><td></td></tr></tbody></table>

<caption>Figure 5: Details of the commodity graph $D^m(\mathbf{r}, \mathbf{r}_{\text{reload}})$ for the reload fixed-path vehicles neighborhood.</caption>

We call it perturbation because it may slightly increase the cost of the IRP solution.

**Customer reinsertion problem.** Let us define the customer reinsertion problem. Once the customer $c \in C$ is removed from the solution – that is to say, removed from the routes that deliver to it, leading to zero delivery $z^+$ in the inventory dynamics (6g) – we need to reinsert it in the solution, using only former routes and new direct routes. This means choosing: 1) The insertion position of customer $c$ in each route of the solution in which it is inserted, keeping the relative order of the other stops unchanged. 2) The quantity of each commodity to be delivered by those former routes where $c$ is inserted. 3) New direct routes (path, timing and quantities) to deliver $c$. It can be formulated as a MILP akin to the generic one defined in Section 3.2. The reload fixed-path vehicles and customer reinsertion formulations are different. In the former, we consider only a set of given fixed-path routes in the route subgraph, the associated customers, and one depot. The decisions are selecting or not each route and fixing the quantities to be delivered to each stop. Both inventory and routing costs induced are exactly modelled. In the latter, a customer is removed from every delivery of the current solution. With the reinsertion MILP, we decide in which former route we insert the customer, at which position, and how much of each commodity we deliver to it. We also decide if we create new direct routes from depots to this particular customer and how much of each commodity we
send through them. New direct routes are approximated by one large bin per depot origin and departure date. Inventory dynamics are exactly modelled, at the customers impacted by the insertion due to delays, and at the depots.

**The customer insertion MILP.** Let $\mathbf{y} = (\mathbf{y}_m)_{m \in M_c}$ be the commodity flow variable. Instead of the indicator route variables in Section 3.2, we use another type of graph to model the vehicles with a flow variable $\mathbf{x}$. Indeed, a flow is a convenient tool to model the fact that we have multiple insertion positions in a given route for customer $c$, and we can choose at most one of them. We link the flow variables on the route subgraphs of the commodity graphs with this vehicle flow. It leads to the MILP:

$$ \min_{\mathbf{y},\mathbf{x}} \sum_{m \in M_c} \mathbf{y}_m^\top \mathbf{c}_m + c_x^\top \mathbf{x} \qquad (\text{Cust-MILP}) $$

$$ \text{subject to} \quad \mathbf{y}_m \in \mathcal{C}(\mathcal{D}_c^m, \mathbf{y}_m^{\min}, \mathbf{y}_m^{\max}), \quad \forall m \in M_c \qquad (11a) $$

$$ \mathbf{x} \in \mathcal{C}(\mathcal{D}_c^{\text{veh}}, \mathbf{x}^{\min}, \mathbf{x}^{\max}) \qquad (11b) $$

$$ \sum_{m \in M_c} y_{ma} l_m \le x_a L, \quad \forall a = (d \to c), \quad \forall d \in D \qquad (11c) $$

$$ \sum_{m \in M_c} y_{ma} l_m \le x_a L_{\text{free}}^r, \quad \forall a = (r \to (r, s)), \quad \forall r \in \mathbf{r}, \quad \forall s \in [|P^r|] \qquad (11d) $$

$$ \mathbf{x} \in \mathbb{Z}, \mathbf{y} \in \mathbb{Z} \qquad (11e) $$

The vehicle flow as well as each commodity flow must respect circulation constraints. In Equation (11c), we force the amount of commodities to be sent through new direct routes from the depots to the customer $c$ not to exceed the total content size of the vehicles involved. Indeed, we consider the total vehicle capacity with this constraint, and not individual vehicles each of capacity $L$. Similarly, Equation (11d) does so for the former routes $r \in \mathbf{r}$ at position $s$ having remaining loading space $L_{\text{free}}^r$.

**Details of the commodity flow graphs.** Based on the generic graph structure Section 3.2, we define new commodity graphs $\mathcal{D}_c^m = (\mathcal{V}_c^m, \mathcal{A}_c^m)$ for each commodity $m \in M_c$ for the reinsertion of customer $c \in C$. We highlight the fact that here only one customer is involved in the problem, so only one customer subgraph is present. The depots and customer subgraphs are introduced in Section 3.2. We now define the route subgraph specific to this customer reinsertion MILP. First, we have a vertex **delivery other customers** in the route subgraph. It is connected to the depot subgraphs to model the quantities sent to other customers. Figure 6 shows the specific commodity graph. In arc annotations, capacities are given between brackets (e.g. $[0, \kappa_{mdt}$], and costs without (e.g. $c_{md}^{exc}$). Dotted arrows are related to the shared artificial vertices, continuous ones to depots, customers and route subgraphs. When not stated in the table, **Min** is 0, **Cost** is 0 and **Max** is $\infty$. In this commodity graph, we consider two types of routes. 1) new direct routes are modelled by direct arcs of the type $(t, d, \text{morning}) \to (t', c, \text{evening})$ and do not involve additional vertices. 2) former routes in which we can insert customer $c$. They are modelled with one vertex $(t, r)$ for route $r$ starting on day $t$, connected to the starting depot morning vertex and to each vertex of the form $(t_s^r, s, r)$. The vertex $(t_s^r, s, r)$ is related to the possible insertion position $s$ in route $r$ leading to an arrival day $t_s^r$ at customer $c$. It is connected to the customer subgraph. The details of the route arcs are in the table of Figure 6. We cannot only add one vertex per former route, since the optimal insertion position depends on the commodity flows, and not only on routing costs. Therefore, the optimal insertion position cannot be pre-computed. Besides, we do not introduce an approximate cost for the transportation arcs. Instead, we define another graph $\mathcal{D}_c^{\text{veh}} = (\mathcal{V}_c^{\text{veh}}, \mathcal{A}_c^{\text{veh}})$ as follows.

**Details of the vehicle flow graph.** We consider the $\mathcal{V}_c^{\text{veh}}$ vertices:

* Artificial vertices labelled: **source** and **sink**.
<figure><img src="image_9.png" /></figure>

<table><caption>(b) Additional commodity graph arcs compared with Table 1</caption><thead><tr><th>Subgraph</th><th>Arc description</th><th>Origin</th><th>Destination</th><th>Min</th><th>Max</th><th>Cost</th></tr></thead><tbody><tr><td>Routes</td><td>Delivery other customers</td><td>(t, d, morning)</td><td>delivery other customers</td><td>z<sub>mdt</sub><sup>-</sup></td><td>z<sub>mdt</sub><sup>-</sup></td><td></td></tr><tr><td>Routes</td><td>Transport d → c new route</td><td>(t, d, morning)</td><td>(t + ⌊<sup>τ<sub>dc</sub></sup>⁄<sub>τ<sub>max</sub></sub>⌋, c, evening)</td><td></td><td></td><td></td></tr><tr><td>Routes</td><td>Transport d → r former route</td><td>(t, d, morning)</td><td>(t, r)</td><td></td><td></td><td></td></tr><tr><td>Routes</td><td>Transport r → (r, s) former route</td><td>(t, r)</td><td>(t<sub>s</sub><sup>r</sup>, s, r)</td><td></td><td></td><td></td></tr><tr><td>Routes</td><td>Transport (r, s) → c former route</td><td>(t<sub>s</sub><sup>r</sup>, s, r)</td><td>(t<sub>s</sub><sup>r</sup>, c, evening)</td><td></td><td></td><td></td></tr><tr><td>Artificial</td><td>Circulation</td><td>delivery other customers</td><td>sink</td><td></td><td></td><td></td></tr></tbody></table>

Figure 6: Details of the commodity graph $D_c^m$ for the customer reinsertion neighborhood.

* Depots vertices labelled: (t, d, morning) for $t \in [T]$, $d \in D$.

* Customer vertices labelled: (t, c, evening) for $t \in [T]$.

* Routes vertices labelled: (t, r) for each former route $r \in r$ that can include an additional stop and that does not reach the vehicle capacity $L$ already.

* Nodes per insertion position in former routes: ($t_s^r$, $s$, $r$) for each position $s$ at which we can insert customer $c$ in $r$ without exceeding the time horizon $T$. The date $t_s^r$ on which the route delivers customer $c$ can be pre-computed.

The arcs of the customer reinsertion vehicle flow graph are defined in Table 3. We see that most of its structure is shared with the commodity flow graph described above. Besides, the costs $c_{sc}^r$ are the ones induced by the insertion of customer $c$ at position $s$ in the stops of route $r$. They involve routing and inventory considerations at the other customers delivered by route $r$, because of the delays. They can be exactly computed considering the former list of $r$'s stops, and enumerating the insertion possibilities. We emphasize the 1-maximum capacity on the arcs of the form (t, d, morning) → (t, r) combined with the circulation constraint enforce that at most one insertion position is chosen in former routes.
<table><caption>Table 3: Arcs of the customer reinsertion vehicle flow graph $D_c^{\text{veh}}$. When not stated, Min is 0, Cost is 0 and Max is $\infty$.</caption><thead><tr><th>Arc description</th><th>Origin</th><th>Destination</th><th>Min</th><th>Max</th><th>Cost</th></tr></thead><tbody><tr><td>Start routes</td><td>source</td><td>(t, d, morning)</td><td></td><td></td><td></td></tr><tr><td>End routes</td><td>(t, c, evening)</td><td>sink</td><td></td><td></td><td></td></tr><tr><td>Transport (d → c) new route</td><td>(t, d, morning)</td><td>(t + ⌊$\frac{\tau_{dc}}{r_{max}}$⌋, c, evening)</td><td></td><td></td><td>$c^{\text{veh}} + c^{\text{stop}} + c^{\text{km}}\Delta_{dc}$</td></tr><tr><td>Transport (d → r) former route</td><td>(t, d, morning)</td><td>(t, r)</td><td>1</td><td></td><td></td></tr><tr><td>Transport (r → (r, s)) former route</td><td>(t, r)</td><td>(t<sub>s</sub><sup>r</sup>, s, r)</td><td></td><td></td><td>$c_{sc}^r$</td></tr><tr><td>Transport ((r, s) → c) former route</td><td>(t<sub>s</sub><sup>r</sup>, s, r)</td><td>(t<sub>s</sub><sup>r</sup>, c, evening)</td><td></td><td></td><td></td></tr><tr><td>Circulation</td><td>source</td><td>sink</td><td></td><td></td><td></td></tr><tr><td>Circulation</td><td>sink</td><td>source</td><td></td><td></td><td></td></tr></tbody></table>

**Proposition 6.1** *The problem (Cust-MILP) is a relaxation of the optimal customer reinsertion problem in the IRP solution, where the new direct routes are aggregated per routing arc ($d \to c$).*

We sketch the proof. The constraints for new direct routes do not exactly model the routing structure: instead of a set of vehicles having each a content size $L$, it is as if we had one large vehicle with the total content size on each arc $d \to c$. More precisely, given $d \in D$, filling $p \in \mathbb{Z}^+$ vehicles of content size $L$ and cost $c^{\text{veh}} + c^{\text{stop}} + c^{\text{km}}\Delta_{dc}$ each enables less loading freedom than filling one large vehicle of content size $pL$ with same total cost $p(c^{\text{veh}} + c^{\text{stop}} + c^{\text{km}}\Delta_{dc})$. In the objective function, the costs $c_m$ for $m \in M_c$ and $c_x$ stem from the arc features stated in Table 1, the table of Figure 6 and in Table 3. Because of their values and because we relax the routing structure, this problem is a relaxation of the customer reinsertion in the current solution of the IRP.

**Rebuilding routes.** All the decisions we make are encapsulated in the commodity flow variable $y$. From $y$, we can easily update the solution $r$, filling former routes with the indicated quantities. Then, a bin packing problem is solved approximately per depot and per day to create the new direct routes from the flows on direct routes arcs. Aggregating per ($d \to c$) arc enables us to derive a MILP of reasonable size, but this relaxation can lead to a potential cost increase. Our LNS accepts the new IRP solution if its true cost is not more than 1% greater than the cost of the previous solution.

## 6.3 Commodity Reinsertion Subroutine

The last subroutine we detail here is the commodity reinsertion. As for the customer reinsertion, the idea is to perturb the solution at a broad scale, possibly increasing the cost.

**Commodity reinsertion problem.** We are interested in the following problem. Once a commodity $m \in M$ is removed from every delivery of the current solution, we want to choose the quantity of $m$ to send through each former route, and through new direct routes. Former routes are the routes of the current solution with the commodity $m$ removed that have a remaining load. The restriction to new direct routes is a choice to quickly compute a solution.
**The commodity insertion MILP.** As for the customer reinsertion, we couple a commodity flow variable y with a vehicle flow x. The resulting MILP is given by:

$$ \min_{\mathbf{y},\mathbf{x}} \quad \mathbf{y}_m^\top \mathbf{c}_m + c_x^\top \mathbf{x} \qquad (\text{Comm-MILP}) $$

$$ \text{subject to} \quad \mathbf{y}_m \in \mathcal{C}(\mathcal{D}_m^m, \mathbf{y}_m^{\min}, \mathbf{y}_m^{\max}) \qquad (12a) $$

$$ \mathbf{x} \in \mathcal{C}(\mathcal{D}_m^{\text{veh}}, \mathbf{x}^{\min}, \mathbf{x}^{\max}) \qquad (12b) $$

$$ y_{ma} \le \left\lfloor \frac{L}{\ell_m} \right\rfloor x_a, \quad \forall a = (d \to c) \quad \forall d \in D_m \quad \forall c \in C_m \qquad (12c) $$

$$ \mathbf{x} \in \mathbb{Z}, \mathbf{y} \in \mathbb{Z} \qquad (12d) $$

We highlight only one commodity flow and one vehicle flow are involved here. Based on the subgraphs introduced in Section 3.2, the inventory constraints and costs are exactly modelled. The routing structure is approximated, with a cost paid per unit of a single “large vehicle with total content size” per arc ($d \to c$) for $c \in C_m$ and $d \in D_m$, modelled by constraint (12c).

**Details of the commodity flow graph.** Let $m \in M$ be the commodity we consider. As previously, we build a commodity graph $\mathcal{D}_m^m = (\mathcal{V}_m^m, \mathcal{A}_m^m)$. It shares the backbone structure with customer and depot subgraphs as well as artificial vertices defined in Section 3.2. We now detail the specific route subgraph for the commodity reinsertion problem. It aims at modelling both former routes and new direct routes. For a former route $r \in \mathbf{r}$, for $c$ visited by $r$, we denote by $s \in [|P^r| - 1]$ its position in the list of visited stops, and by $t_s^r$ the date on which it is delivered. We have a vertex $(t_s^r, c, r)$ in the route subgraph. We then have the arcs of Table 4 to explicitly model the flow of the commodity $m$ through former and new routes, in addition to those defined in Table 1.

<table><caption>Table 4: Additional arcs of the commodity reinsertion commodity flow graph compared to Table 1. When not stated, Min is 0, Cost is 0 and Max is $\infty$.</caption><thead><tr><th>Subgraph</th><th>Arc description</th><th>Origin</th><th>Destination</th><th>Min</th><th>Max</th><th>Cost</th></tr></thead><tbody><tr><td>Routes</td><td>Transport <i>d</i> → <i>c</i> new route</td><td>(<i>t</i>, <i>d</i>, morning)</td><td>(<i>t</i> + <sup><i>t</i><sub><i>d</i></sub></sup>/<sub><i>t</i><sub>max</sub></sub>, <i>c</i>, evening)</td><td></td><td></td><td></td></tr><tr><td>Routes</td><td>Transport (<i>d</i> → <i>c</i><sub>1</sub>) former route <i>r</i></td><td>(<i>t</i><sup>r</sup>, <i>d</i>, morning)</td><td>(<i>t</i><sub>1</sub><sup>r</sup>, <i>c</i><sub>1</sub>, <i>r</i>)</td><td></td><td></td><td></td></tr><tr><td>Routes</td><td>Transport (<i>c</i><sub>s</sub> → <i>c</i><sub>s+1</sub>) former route <i>r</i></td><td>(<i>t</i><sub>s</sub><sup>r</sup>, <i>c</i><sub>s</sub>, <i>r</i>)</td><td>(<i>t</i><sub>s+1</sub><sup>r</sup>, <i>c</i><sub>s+1</sub>, <i>r</i>)</td><td></td><td></td><td></td></tr><tr><td>Routes</td><td>Deliver <i>c</i><sub>s</sub> former route <i>r</i></td><td>(<i>t</i><sub>s</sub><sup>r</sup>, <i>c</i><sub>s</sub>, <i>r</i>)</td><td>(<i>t</i><sub>s</sub><sup>r</sup>, <i>c</i><sub>s</sub>, evening)</td><td></td><td></td><td></td></tr></tbody></table>

**Details of the vehicle flow graph.** This flow graph is only used to model new direct routes. Indeed, all former routes are reused in this problem: we only choose the commodity flow through them. The only routing decision is related to the creation of new direct routes. Therefore, the commodity reinsertion vehicles graph vertices $\mathcal{V}_m^{\text{veh}}$ are the following: 1) Artificial vertices labelled: source and sink. 2) Depot vertices labelled: (<i>t</i>, <i>d</i>, morning) for $t \in [T]$, $d \in D_m$. 3) Customer vertices labelled: (<i>t</i>, <i>c</i>, evening) for $t \in [T]$, $c \in C_m$. Let $\mathcal{A}_m^{\text{veh}}$ be the arcs of the vehicle flow graph. They are detailed in Table 5.

**Rebuilding routes.** The commodity reinsertion subroutine proceeds as follows. From the information written in the y flow variable, it fills the former routes by decoding the corresponding flow. Then, it iteratively creates new direct routes to zero the corresponding commodity flow in y. In contrast with the previous section, restricting ourselves to new direct routes is a limitation. It applies a routing local search to address the potential cost rise induced. The new solution is accepted by the LNS, even if its cost is higher.

# 7 Computational Experiments

First, we want to analyze the quality of the solutions of our LNS. We compare them with the route-based matheuristic in terms of costs. We then perform additional experiments to understand the contribution of
<table><caption>Table 5: Arcs of the commodity reinsertion vehicle flow graph $D_m^{\text{veh}}$. When not stated, Min is 0, Cost is 0 and Max is $\infty$.</caption><thead><tr><th>Arc description</th><th>Origin</th><th>Destination</th><th>Min</th><th>Max</th><th>Cost</th></tr></thead><tbody><tr><td>Start routes</td><td>source</td><td>(t, d, morning)</td><td></td><td></td><td></td></tr><tr><td>Transport (d → c)</td><td>(t, d, morning)</td><td>(t + ⌊<sup>τ<sub>dc</sub></sup>⁄<sub>τ<sub>max</sub></sub>⌋, c, evening)</td><td></td><td></td><td>c<sup>veh</sup> + c<sup>stop</sup> + c<sup>km</sup>Δ<sub>dc</sub></td></tr><tr><td>End routes</td><td>(t, c, evening)</td><td>sink</td><td></td><td></td><td></td></tr><tr><td>Circulation</td><td>source</td><td>sink</td><td></td><td></td><td></td></tr><tr><td>Circulation</td><td>sink</td><td>source</td><td></td><td></td><td></td></tr></tbody></table>

the distinct components of our LNS to its performance.

As previously mentioned, all of our code is available on GitHub (Bouvier et al. 2023c). The library of instances (Bouvier et al. 2023a) and the solutions we obtained (Bouvier et al. 2023b) can be found on Zenodo. Note that both the instances and the solutions are downloaded automatically by our test suite, which is run online following every commit. This is meant to improve the reproducibility and reliability of our experiments.

## 7.1 Experimental Setting

**Instances.** Renault gives us access to 71 IRP instances, each of them corresponding to their return logistics at the European scale and over roughly 20 days. They correspond to the same industrial use case but to different periods. We show the instances' number of depots, number of customers and number of commodities on Figure 7. The maximum number of stops is $S_{max} = 3$ to comply with the car manufacturer's requirements. Additional details on the instances are available in Figure 11 in Appendix C. They highlight why our instances are hard to solve, showing how release and demand are spread respectively among depots and customers, as well as how grouping customers in routes with multiple stops is essential.

<figure><img src="image_10.png" /><figcaption>Figure 7: Histograms of the dimensions of the extracted instances.</figcaption></figure>

**Implementation.** Our code is in the Julia language (Bezanson et al. 2017). We use Gurobi optimizer (Gurobi Optimization, LLC 2021) to solve LPs and MILPs, JuMP (Dunning, Huchette, and Lubin 2017) to model mathematical programs, and Graphs.jl (Fairbanks et al. 2021) to define graph structures. We proceed to a warm start using the current solution for each MILP. Our experiments are run on a computing cluster with 189Go of RAM and 32 processors. Multithreading is enabled for Gurobi only.

**Hyperparameters.** Several hyperparameters are introduced for our LNS: the number of routing local search iterations per LNS step, the percentage of random pairs of routes considered in the routing local search, the number of reload fixed-path vehicles, customer, and commodity reinsertion iterations per outer LNS step. To highlight the effect of our new perturbations, we tune the number of customer and commodity reinsertion steps per outer LNS step, as shown in Appendix E in Table 8. The best configuration we find for the short 90 minutes runs is $|C|$ customer reinsertion (the whole set of customers) and 10 commodity
reinsertion steps. We browse the whole set of depots in the reload fixed-path vehicles subroutine, and proceed to one iteration of routing local search per outer LNS step. Besides, we run our code for 90 minutes on each instance. We include longer run results with 300 minutes to solve instances with $S_{max}$ = 10 in the extended experiments of Section 7.3.4.

**Experimental design.** We run our LNS with the hyperparameters described above on the 71 real instances. We compare solutions with the route-based matheuristic, as well as with some versions of the LNS where a subroutine is removed in an ablation study. Every algorithm has the same time budget of 90 minutes. Nonetheless, since both the initialization + local search and route-based matheuristic are constructive, they end before the 90 minutes time budget. In Appendix E, when tuning the LNS hyperparameters, we show in Table 8 a version of the LNS with zero customer and commodity reinsertion step. This version of the LNS can be seen as an iterative extension of the route-based matheuristic. Our aim is to analyze the performance of the LNS, as well as the effect of each of its components. Longer runs with $S_{max}$ = 10 aim at showing rooms for improvement.

## 7.2 Performance analysis

The cost of the IRP solution after our LNS is the first natural metric to evaluate our approach. We can compare it with the cost after the route-based matheuristic defined in Section 3. On Figure 8, we show the box plots of the cost due to the depots' inventory, the customers' inventory, the customers' shortage, the vehicles, the stops and the kilometers. The orange lines indicate the median over the instances, the ends of the boxes the extreme quartiles ($Q1$ and $Q3$), and the whiskers the range [$Q1-1.5(Q3-Q1)$, $Q3+1.5(Q3-Q1)$]. Outlier points correspond to data outside the whiskers. The average total cost over the instances solved by the initialization + local search algorithm is 2.88M€, compared with 2.48M€ after the route-based matheuristic, and 2.10M€ after the LNS.

**Result 7.1** *The LNS enables 37% and 18% cost savings compared with the initialization + local search and route-based matheuristic respectively.*

On Figure 8 we notice that almost all the components of the cost are lower: the depots inventory costs, the shortage costs, as well as the routing costs (vehicles costs, stop costs and kilometer costs) have quartiles corresponding to smaller values after the LNS.

**Result 7.2** *The LNS improves almost every component of the cost.*

## 7.3 Extended analysis

### 7.3.1 Time per Operator

We proceed to a more detailed analysis of the solution steps illustrated on Figure 2. We first compare the total time spent in the initialization + local search algorithm, the routing local search, the reload fixed-path vehicles large neighborhood, and the customer and commodity reinsertion perturbations. We emphasize that if we just proceed to one pass over every neighborhood and perturbation without restricting the depots, customers, or commodities considered for the MILP-based operators, we can barely do two outer iterations of the LNS in 90 minutes. With the hyperparameters we set, we do 2.5 outer iterations of the LNS on average. Adaptive approaches (Gendreau and Potvin 2010) are thus excluded in our context. On the left of Figure 9 we show that the initialization + local search algorithm is indeed fast, taking 2.89 minutes on average. The total time spent in the routing local search amounts to 5.94 minutes on average. The reload fixed-path vehicles neighborhood, customer reinsertion and commodity reinsertion perturbations account for respectively 19.2, 51.0 and 11.2 minutes. About 90% of the time is spent solving MILPs to modify the solution at a large scale. Therefore, there is not much to gain from improving our implementation further, because the solver itself is out of our control.
<figure><img src="image_11.png" /><figcaption>Figure 8: Box plots of the distribution of the solution cost per origin over instances. Green rectangles are related to the route-based matheuristic, blue ones to the LNS.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 9: Box plots of the time spent, and of the cost gain per CPU time per operator.</figcaption></figure>

**Result 7.3** *There is not much performance to gain by further improving our implementation, since most of the time is spent to solve MILPs.*

### 7.3.2 Cost Gain Over Time

On the right of Figure 9, we display the box plots of the cost gain per CPU time (in €/minute) of our four main LNS components: the routing local search, reload fixed-path vehicles neighborhood, customer reinsertion, and commodity reinsertion. We highlight that the routing local search has the highest gain per CPU time, 25.7k€/minute on average. Its large variance over the instances can be related to the balance between routing and inventory costs, which may vary from one instance to another depending on the demand and release profiles. The unit costs themselves do not change between instances. We thus highlight the routing neighborhoods are efficient and crucial. It is partly due to the structure of our perturbations: a customer or commodity reinsertion step creates new routes that are only direct. The reason for this choice is to avoid spending too much time solving greater perturbation MILPs, for they only involve a portion of the solution. It is therefore useful to combine them with routing neighborhoods that merge or mix routes in different manners to increase route length.

The reload fixed-path vehicles neighborhood, customer and commodity reinsertion perturbations alter both routing and inventory variables, with 7.95k€/minute, 8.75k€/minute and 2.77k€/minute ratios on av-
erage. Therefore, the descent neighborhoods (routing local search and reload fixed-path vehicles subroutines) decrease the cost on average, but so do the perturbations, even though they are designed to escape from local minima. We show further details with an ablation study in Section 7.3.3, and tuning in Appendix E. We emphasize that separating the gains per operator is not totally obvious. It is indeed the mix between the axes of the IRP (from routes to commodities and customers) that enables us to substantially reduce the cost in our LNS.

### 7.3.3 Ablation study

During the ablation study, we run the LNS with the same hyperparameters, but without the customer reinsertion, commodity reinsertion or reload fixed-path vehicles neighborhood respectively. We compare the cumulative distributions of the gaps on Figure 10. We must emphasize our lower bound derived in Section 4.1 is not tight, and therefore leads to overestimated gaps. It seems wiser to use these gaps as a relative metric to compare the solution processes, and not as an absolute indicator of solution quality. The first remark we can make is that the route-based matheuristic (blue curve) performs better than our initialization + local search algorithm (orange curve), but worse than the LNS, even when one of the LNS components is removed. This confirms the conclusion drawn from the cost analysis in Section 7.2, emphasizing the performance of our LNS. Besides, for this short 90 minutes run, we show that the customer reinsertion perturbation and the reload fixed-path vehicles large neighborhood are crucial for the LNS performance. Indeed, removing the customer reinsertion subroutine (red curve) entails a 12% gap increase on average. Removing the reload fixed-path vehicles leads to 6% average increase (purple curve). Those two subroutines are complementary in the sense that they modify distinct parts of the solution. The commodity reinsertion perturbation seems to be less crucial for the performance in the short 90 minutes run. Indeed, removing it from the LNS entails only a 2% average gap increase. It may be due to the fact that changing the commodity flow axis is also at the core of the reload fixed-path vehicles subroutine. Nonetheless, the commodity reinsertion perturbation potentially creates new direct routes. We show additional experiments in Appendix E with longer runs to emphasize its utility.

<figure><img src="image_13.png" /><figcaption>Figure 10: Cumulative distributions of the gap among instances solutions.</figcaption></figure>

**Result 7.4** All the neighborhoods are useful. The commodity reinsertion perturbation only becomes fully efficient for longer runs.
### 7.3.4 Results with Longer Routes and Larger CPU Time

As stated above, limiting the routes to $S_{max} = 3$ stops maximum is restrictive, though this is an important constraint for Renault. To emphasize the benefits induced by longer routes, we solve the instances with $S_{max} \in \{3, 10\}$, and a time limit of 300 minutes. Note that both the initialization + local search and route-based matheuristic end before the time limit. The rest of the hyperparameters are unchanged. We summarize the average cost impact with Table 6. First, increasing the LNS running time from 90 minutes to 300 minutes with $S_{max} = 3$ reduces the average gap by an additional seven percent. More experiments on cost evolution with LNS iterations are analyzed in Appendix F. There we also show that we can reach lower costs with $S_{max} = 10$ with each algorithm, observing a large gain after the initialization + local search algorithm. Numerical experiments also demonstrate that our LNS can handle the combinatorics induced by longer routes.

**Result 7.5** *Increasing the running time enables to find better solutions. The LNS has not converged after 90 minutes.*

<table><caption>Table 6: Comparison of the average costs over instances between $S_{max} = 3, 10$ after 300 minutes.</caption><thead><tr><th></th><th>Initialization +<br>local search</th><th>Route-based<br>matheuristic</th><th>LNS</th><th>Lower bound</th></tr></thead><tbody><tr><td>Average cost $S_{max} = 3$</td><td>2.88M€</td><td>2.48M€</td><td>2.01M€</td><td>1.35M€</td></tr><tr><td>Average cost $S_{max} = 10$</td><td>2.41M€</td><td>2.22M€</td><td>1.81M€</td><td>1.35M€</td></tr></tbody></table>

**Result 7.6** *Increasing the maximum length of a route $S_{max}$ enables additional savings when running time can be increased.*

## 8 Conclusions

In this study motivated by an industrial partnership with Renault, we consider European-scale continuous-time multi-attribute IRP instances. This inherently hard problem has to be solved in a limited time of 90 minutes every day. To do so, we design a large neighborhood search based on the generalization of TSP and SDVRP neighborhoods Section 5, a large neighborhood Section 6.1 inspired by recent matheuristics designed for the IRP (Bertazzi et al. 2019; Archetti, Boland, and Speranza 2017), and two new perturbations Section 6.2 and Section 6.3 based on MILPs. We also derive an initialization + local search algorithm that relies on flows Section 4.1. It allows us to quickly initialize our IRP instances with non-trivial solutions, and to derive a lower bound. To the best of our knowledge, this lower bound is unknown in the IRP literature. We extract and process a dataset of 71 European-scale multi-attribute IRP instances that we make available publicly. Numerical experiments in Section 7 show the results of our LNS. We highlight that it outperforms the route-based matheuristic, and that each component of the LNS brings useful contributions.

We also emphasize some limits and perspectives. The choice of the neighborhoods is hard to define a priori, and we do not have time to browse the space of neighborhoods repeatedly. Parallel computing and multi-threading could be some perspectives to consider, although even small changes imply inventory dynamics over the whole time horizon, leading to overlaps between neighborhoods. Cache use would thus be a challenge. Some techniques of machine learning for operations research could also be considered in this direction, for example in a reinforcement learning (Wu et al. 2021) or structured learning (Parmentier 2022) paradigm. Besides, we suffer from the poor quality of our lower bound to compute gaps. Deriving a better relaxation for this multi-attribute IRP is a challenge that has not been addressed in the literature, to the best of our knowledge. We could add valid inequalities leveraging the research on exact solutions of smaller problems (Manousakis et al. 2021; Desaulniers, Rakke, and Coelho 2015) to improve the quality of our neighborhoods. A version of the LNS we designed is currently used in production at Renault, with
additional statistical treatments due to the noise observed in release and demand forecasts. Our simulations enable us to estimate a reduction of thousands of tons of CO2 and millions of euros per year. Dealing with the stochastic IRP (Nolz, Absi, and Feillet 2014; Coelho, Laporte, and Cordeau 2012) is at the core of our future work.

**Acknowledgments** We are grateful to the Renault supply-chain and IT teams for the partnership we have set and the industrial motivation they provide us with, especially to Alain Nguyen, Thaddeus Leonard, Nicusor-Eugen Plescan, Christian Serrano, Ludovic Doudard and Aimé-Frédéric Rosenzweig. We also would like to thank Vincent Leclère for his advice on the form of the article.

## References

Archetti, C., L. Bertazzi, G. Laporte, and M. G. Speranza (Aug. 2007). "A Branch-and-Cut Algorithm for a Vendor-Managed Inventory-Routing Problem". In: *Transportation Science* 41.3, pp. 382-391. ISSN: 0041-1655. DOI: 10/b9kkm7 (cit. on p. 2).

Archetti, C., N. Boland, and M. G. Speranza (Apr. 2017). "A Matheuristic for the Multivehicle Inventory Routing Problem". In: *INFORMS Journal on Computing* 29.3, pp. 377-387. ISSN: 1091-9856. DOI: 10/gbsn75 (cit. on pp. 2, 3, 8, 25).

Archetti, C., G. Guastaroba, D. L. Huerta-Muñoz, and M. G. Speranza (Nov. 2021). "A Kernel Search Heuristic for the Multivehicle Inventory Routing Problem". In: *International Transactions in Operational Research* 28.6, pp. 2984-3013. ISSN: 0969-6016, 1475-3995. DOI: 10/gkzwqh (cit. on pp. 2, 8).

Archetti, C. and M. G. Speranza (2016). "The Inventory Routing Problem: The Value of Integration". In: *International Transactions in Operational Research* 23.3, pp. 393-407. ISSN: 1475-3995. DOI: 10/gctq3d (cit. on pp. 1, 3).

Archetti, C. and M. G. Speranza (2008). "The Split Delivery Vehicle Routing Problem: A Survey". In: *The Vehicle Routing Problem: Latest Advances and New Challenges*. Ed. by B. Golden, S. Raghavan, and E. Wasil. Operations Research/Computer Science Interfaces. Boston, MA: Springer US, pp. 103-122. ISBN: 978-0-387-77778-8. DOI: 10.1007/978-0-387-77778-8_5 (cit. on p. 3).

Benoist, T., F. Gardi, A. Jeanjean, and B. Estellon (Mar. 2011). "Randomized Local Search for Real-Life Inventory Routing". In: *Transportation Science* 45.3, pp. 381-398. ISSN: 0041-1655. DOI: 10.1287/trsc.1100.0360 (cit. on p. 2).

Bertazzi, L., L. C. Coelho, A. De Maio, and D. Laganà (Feb. 2019). "A Matheuristic Algorithm for the Multi-Depot Inventory Routing Problem". In: *Transportation Research Part E: Logistics and Transportation Review* 122, pp. 524-544. ISSN: 1366-5545. DOI: 10/ggdvkf (cit. on pp. 2, 3, 8, 25).

Bezanson, J., A. Edelman, S. Karpinski, and V. B. Shah (2017). "Julia: A fresh approach to numerical computing". In: *SIAM Review* 59.1, pp. 65-98. DOI: 10.1137/141000671. URL: https://epubs.siam.org/doi/10.1137/141000671 (cit. on pp. 3, 21).

Bouvier, L., G. Dalle, A. Parmentier, and T. Vidal (July 2023a). *Continent-scale inventory routing instances*. Zenodo. DOI: 10.5281/zenodo.8177237. URL: https://doi.org/10.5281/zenodo.8177237 (cit. on pp. 3, 21).

— (July 2023b). *Continent-scale inventory routing solutions*. Zenodo. DOI: 10.5281/zenodo.8177271. URL: https://doi.org/10.5281/zenodo.8177271 (cit. on pp. 3, 21).

— (July 2023c). *InventoryRoutingLNS.jl*. Version v0.1.1. DOI: 10.5281/zenodo.8179161. URL: https://doi.org/10.5281/zenodo.8179161 (cit. on pp. 3, 21).

Campbell, A. M. and M. W. P. Savelsbergh (Nov. 2004). "A Decomposition Approach for the Inventory-Routing Problem". In: *Transportation Science* 38.4, pp. 488-502. ISSN: 0041-1655. DOI: 10/fkn54g (cit. on p. 2).

Chitsaz, M., J.-F. Cordeau, and R. Jans (Feb. 2019). "A Unified Decomposition Matheuristic for Assembly, Production, and Inventory Routing". In: *INFORMS Journal on Computing* 31.1, pp. 134-152. ISSN: 1091-9856, 1526-5528. DOI: 10/gfvn3j (cit. on p. 2).
Coelho, L. C. and G. Laporte (Nov. 2013). "A Branch-and-Cut Algorithm for the Multi-Product Multi-Vehicle Inventory-Routing Problem". In: *International Journal of Production Research* 51.23-24, pp. 7156-7169.
ISSN: 0020-7543. DOI: 10/gjvnnf (cit. on p. 2).

Coelho, L. C., A. De Maio, and D. Laganà (Dec. 2020). "A Variable MIP Neighborhood Descent for the Multi-Attribute Inventory Routing Problem". In: *Transportation Research Part E: Logistics and Transportation Review* 144, p. 102137. ISSN: 13665545. DOI: 10/gm954p (cit. on pp. 2, 8).

Coelho, L. C., G. Laporte, and J.-F. Cordeau (2012). *Dynamic and Stochastic Inventory-Routing*. Tech. rep. CIRRELT-2012-37. CIRRELT Montreal (cit. on p. 26).

Cordeau, J.-F., D. Laganà, R. Musmann, and F. Vocaturo (Mar. 2015). "A Decomposition-Based Heuristic for the Multiple-Product Inventory-Routing Problem". In: *Computers & Operations Research* 55, pp. 153-166. ISSN: 03050548. DOI: 10/f6w8pp (cit. on p. 2).

Desaulniers, G., J. G. Rakke, and L. C. Coelho (Oct. 2015). "A Branch-Price-and-Cut Algorithm for the Inventory-Routing Problem". In: *Transportation Science* 50.3, pp. 1060-1076. ISSN: 0041-1655. DOI: 10/f8zpgc (cit. on pp. 2, 25).

Dror, M. and P. Trudeau (1990). "Split Delivery Routing". In: *Naval Research Logistics (NRL)* 37.3, pp. 383-402. ISSN: 1520-6750. DOI: 10.1002/nav.3800370304 (cit. on p. 3).

Dunning, I., J. Huchette, and M. Lubin (2017). "JuMP: A Modeling Language for Mathematical Optimization". In: *SIAM Review* 59.2, pp. 295-320. DOI: 10.1137/15M1020575 (cit. on p. 21).

Fairbanks, J., M. Besançon, S. Simon, J. Hoffman, N. Eubank, and S. Karpinski (2021). *JuliaGraphs/Graphs.jl: an optimized graphs package for the Julia programming language*. URL: https://github.com/JuliaGraphs/Graphs.jl/ (cit. on p. 21).

Fischetti, M. and M. Fischetti (2016). "Matheuristics". In: *Handbook of Heuristics*. Ed. by R. Martí, P. Panos, and M. G. Resende. Cham: Springer International Publishing, pp. 1-33. ISBN: 978-3-319-07153-4. DOI: 10.1007/978-3-319-07153-4_14-1 (cit. on p. 2).

Gendreau, M. and J.-Y. Potvin, eds. (2010). *Handbook of Metaheuristics*. Vol. 146. International Series in Operations Research & Management Science. Boston, MA: Springer US. ISBN: 978-1-4419-1663-1 978-1-4419-1665-5. DOI: 10.1007/978-1-4419-1665-5 (cit. on p. 22).

Gurobi Optimization, LLC (2021). *Gurobi Optimizer Reference Manual*. URL: https://www.gurobi.com (cit. on p. 21).

Lagos, F., N. Boland, and M. Savelsbergh (Jan. 2020). "The Continuous-Time Inventory-Routing Problem". In: *Transportation Science*, trsc.2019.0902. ISSN: 0041-1655, 1526-5447. DOI: 10.1287/trsc.2019.0902 (cit. on p. 1).

Manousakis, E., P. Repoussis, E. Zachariadis, and C. Tarantilis (May 2021). "Improved Branch-and-Cut for the Inventory Routing Problem Based on a Two-Commodity Flow Formulation". In: *European Journal of Operational Research* 290.3, pp. 870-885. ISSN: 03772217. DOI: 10.1016/j.ejor.2020.08.047 (cit. on pp. 2, 25).

Nolz, P. C., N. Absi, and D. Feillet (Jan. 2014). "A Stochastic Inventory Routing Problem for Infectious Medical Waste Collection". In: *Networks* 63.1, pp. 82-95. ISSN: 00283045. DOI: 10.1002/net.21523 (cit. on pp. 2, 26).

Parmentier, A. (Jan. 2022). "Learning to Approximate Industrial Problems by Operations Research Classic Problems". In: *Operations Research* 70.1, pp. 606-623. ISSN: 0030-364X, 1526-5463. DOI: 10.1287/opre.2020.2094 (cit. on p. 25).

Savelsbergh, M. and J.-H. Song (July 2008). "An Optimization Algorithm for the Inventory Routing Problem with Continuous Moves". In: *Computers & Operations Research* 35.7, pp. 2266-2282. ISSN: 03050548. DOI: 10/ft9dtw (cit. on p. 1).

Su, Z., Z. Lü, Z. Wang, Y. Qi, and U. Benlic (Mar. 2020). "A Matheuristic Algorithm for the Inventory Routing Problem". In: *Transportation Science*, trsc.2019.0930. ISSN: 0041-1655, 1526-5447. DOI: 10/gm954w (cit. on p. 2).

Touzout, F. A., A.-L. Ladier, and K. Hadj-Hamou (Aug. 2022). "An assign-and-route matheuristic for the time-dependent inventory routing problem". en. In: *European Journal of Operational Research* 300.3,
pp. 1081-1097. ISSN: 03772217. DOI: 10.1016/j.ejor.2021.09.025. URL: https://linkinghub.elsevier.com/retrieve/pii/S0377221721008006 (visited on 02/08/2023) (cit. on p. 2).

Vadseth, S. T., H. Andersson, and M. Stålhane (July 2021). “An Iterative Matheuristic for the Inventory Routing Problem”. In: *Computers & Operations Research* 131, p. 105262. ISSN: 03050548. DOI: 10/gkzwqj (cit. on p. 2).

Wu, Y., W. Song, Z. Cao, and J. Zhang (2021). “Learning Large Neighborhood Search Policy for Integer Programming”. In: *Advances in Neural Information Processing Systems*. Ed. by M. Ranzato, A. Beygelzimer, Y. Dauphin, P. Liang, and J. W. Vaughan. Vol. 34. Curran Associates, Inc., pp. 30075-30087. URL: https://proceedings.neurips.cc/paper/2021/file/fc9e62695def29ccdb9eb3fed5b4c8c8-Paper.pdf (cit. on p. 25).

## A Notations

In Table 7, we recap the main concepts and notations for the dimensions, the initial inventory, the demand and release, the free inventory capacities, the locations graph, the commodities and vehicles lengths, and the unit costs.

## B Reinsertion cost

In Table 3, we introduce $c_{sc}^r$, the cost induced by the insertion of customer $c$ in route $r$ at position $s$. We define it explicitly here. Let $r$ be a route, with day of departure $t^r$, path $P^r = (v_0^r, v_1^r, \dots, v_k^r)$ and quantities delivered $q^r = (q_{ms}^r)_{m \in M, s \in [|P^r|-1]} \in (\mathbb{Z}^+)^{|M| \times (|P^r|-1)}$. Let $c$ be a customer which does not belong to $P^r$ and $s$ its insertion position. We denote by $\tilde{r}$ the route with corresponding insertion, $P^{\tilde{r}}$ the associated path, and $\tilde{I}$ the inventory variables with updated delivery times due to delays induced by the new stop (see Section 2.3. of the paper). The cost $c_{sc}^r$ is induced by the insertion of customer $c$ at position $s$ in the stops of route $r$. It involves routing and inventory considerations at the other customers delivered by route $r$, because of the delays. It is computed as:

$$c_{sc}^r = c^{\text{stop}} + c^{\text{km}} \left[ \Delta_{(c_{s-1},c)} + \Delta_{(c,c_{s+1})} - \Delta_{(c_{s-1},c_{s+1})} \right] + \\ \sum_{\substack{c' \in (v_1^r, \dots, v_k^r), \\ t \ge t^r, \\ m, \sum_{s \in [|P^r|-1]} q_{sm}^r > 0}} \left[ c_{mc'}^{\text{exc}} (\tilde{I}_{mc'}t - \kappa_{mc'}t)^+ - c_{mc'}^{\text{exc}} (I_{mc'}t - \kappa_{mc'}t)^+ \right] + \\ \sum_{\substack{c' \in (v_1^r, \dots, v_k^r), \\ t \ge t^r, \\ m, \sum_{s \in [|P^r|-1]} q_{sm}^r > 0}} \left[ c_{mc'}^{\text{short}} (b_{mc'}^{-} - \tilde{I}_{mc'}(t-1))^+ - c_{mc'}^{\text{short}} (b_{mc'}^{-} - I_{mc'}(t-1))^+ \right].$$

It has one term of routing cost, one of excess inventory cost, and one of shortage cost. The two last terms are reduced to the customers visited by the route $r$, and browsing only days and commodities impacted by route $r$. This local implementation is important for code performance. Besides, the value of $c_{sc}^r$ does not depend on the commodity flows at the inserted customer, contrary to the inventory cost at the starting depot of the routes, and the inventory cost at the inserted customer. It can therefore be pre-computed before solving the customer reinsertion MILP.

## C Instance details.

The complexity of the industrial problem for Renault comes from the tight coupling between commodities, between nodes (depots and customers), between inventory and routing (increased by the continuous-time
<table><caption>Table 7: Notations</caption><thead><tr><th>Name</th><th>Description</th></tr></thead><tbody><tr><td>t</td><td>Day</td></tr><tr><td>T</td><td>Horizon</td></tr><tr><td>d</td><td>Depot (facility)</td></tr><tr><td>D</td><td>Set of depots</td></tr><tr><td>c</td><td>Customer</td></tr><tr><td>C</td><td>Set of customers</td></tr><tr><td>m</td><td>Commodity</td></tr><tr><td>M</td><td>Set of commodities</td></tr><tr><td>M<sub>c</sub></td><td>Set of commodities used by customer <i>c</i></td></tr><tr><td>M<sub>d</sub></td><td>Set of commodities used by depot <i>d</i></td></tr><tr><td>I<sup>0</sup><sub>md</sub></td><td>Initial inventory of <i>m</i> at <i>d</i></td></tr><tr><td>I<sup>0</sup><sub>mc</sub></td><td>Initial inventory of <i>m</i> at <i>c</i></td></tr><tr><td>b<sup>+</sup><sub>mdt</sub></td><td>Number of commodity <i>m</i> released on day <i>t</i> by <i>d</i></td></tr><tr><td>b<sup>-</sup><sub>mct</sub></td><td>Number of commodity <i>m</i> demanded on day <i>t</i> by <i>c</i></td></tr><tr><td>&kappa;<sub>mdt</sub></td><td>Free inventory capacity of <i>m</i> at <i>d</i> on the evening of day <i>t</i></td></tr><tr><td>&kappa;<sub>mct</sub></td><td>Free inventory capacity of <i>m</i> at <i>c</i> on the evening of day <i>t</i></td></tr><tr><td>V</td><td>Set <i>D</i> &cup; <i>C</i> of the vertices of the locations graph</td></tr><tr><td>v</td><td>Node of <i>V</i></td></tr><tr><td>A</td><td>Arcs of the locations graph</td></tr><tr><td>a</td><td>Arc</td></tr><tr><td>&Delta;<sub>a</sub></td><td>Distance in kilometers corresponding to arc <i>a</i></td></tr><tr><td>&tau;<sub>a</sub></td><td>Duration in hours corresponding to arc <i>a</i></td></tr><tr><td>P</td><td>Path in the locations graph</td></tr><tr><td>S<sub>max</sub></td><td>Maximum number of stops in a route</td></tr><tr><td>&tau;<sub>max</sub></td><td>Number of transport hours per day</td></tr><tr><td>&ell;<sub>m</sub></td><td>Length of a commodity <i>m</i></td></tr><tr><td>L</td><td>Length of a vehicle (homogeneous)</td></tr><tr><td>c<sup>exc</sup><sub>md</sub></td><td>Unit excess inventory storage cost of <i>m</i> at <i>d</i> per night</td></tr><tr><td>c<sup>exc</sup><sub>mc</sub></td><td>Unit excess inventory storage cost of <i>m</i> at <i>c</i> per night</td></tr><tr><td>c<sup>short</sup><sub>mc</sub></td><td>Unit shortage cost of <i>m</i> at <i>c</i> per day</td></tr><tr><td>c<sup>km</sup></td><td>Per kilometer cost of the routes</td></tr><tr><td>c<sup>veh</sup></td><td>Unit cost for using a vehicle</td></tr><tr><td>c<sup>stop</sup></td><td>Unit cost for making a stop</td></tr></tbody></table>

aspect), as well as the size of the instances. We provide some details about the instances' structure here to back our discussion. In Figure 11, we focus on one particular instance with 15 depots and 499 customers. The code to generate the plots, as well as the instances are publicly available. On the top-left plot, we show the cumulative normalized total release distribution over depots per commodity. For each commodity (each line in the plot), we sort depots per total release. Therefore, the depot index axis is different for each commodity. Our aim here is to show how the release of each commodity is spread over depots. We see that for a few commodities, only two or three depots release almost the whole proportion. But for most of the commodities, the release is shared by about eight depots. We recall the distribution of the number of depots per instance is visible in Figure 7. Similarly, on the top-right plot, we show the cumulative normalized total demand distribution over customers per commodity. We see that no commodity is shared by the 499 customers. One is shared by more than 200 customers, and most of them are shared by 50 to 100 customers. On the bottom-left plot, we show the histogram over customers of the fraction of a vehicle represented by the average daily demand. On each day and for each customer, the commodities are scaled by their lengths to compute the total demand length of the day, and the daily average is computed. In this way, we see
that a huge proportion of the customers have a demand that fits in less than 10% of a vehicle. Only one customer has a demand that requires 80% of a vehicle. This analysis has to be combined with the last plot of Figure 11. On the bottom-right plot, we show the proportion of a vehicle length represented by the total daily demand length. Apart from the weekend, we observe that the total demand length represents between 30 and 40 times the length of a vehicle. The bin packing, release profiles and continuous-time aspects lead to many more vehicles per day in practice, as illustrated in Figure 12. We thus show it is crucial to create routes with several stops to find good solutions.

<figure><img src="image_14.png" /><figcaption>Figure 11: Details on the release and demand features of an instance.</figcaption></figure>

**Result C.1** *Instances are complex because of the bindings between depots, customers, days and commodities.*

## D Solution details.

We focus here on the best solution found by the LNS after a 90 minutes run, and compare it to the initialization + local search solution, for the instance highlighted in the previous appendix. As for the instance analysis, the code to generate the plot is publicly available, as well as the solutions found. On Figure 12, we display depot, start date, number of stops and customer distributions over the routes of the two solutions. Blue color is related to the initialization + local search, red color to the LNS. The first remark we can make is that there are much fewer routes in the solution of the LNS than in the solution of the initialization + local search, which explains the decrease of the total routing cost after the LNS. On the top-left of Figure 12, we see that both the solution of the initialization + local search and the best one found by the LNS have routes starting from every depot. The distribution structures are similar, although the numbers of routes are distinct. This can be linked to the release and demand structures. On the top-middle plot, we show that except for the weekends, the departure dates are spread over days. The LNS solution has less starting date variations than the initialization + local search solution. We expect it to be linked
to a better management of the time axis, with the challenging couplings between inventory dynamics and continuous-time routes. We also observe the effect of the end of the time axis, with a decreasing number of departures when approaching the horizon. On the top-right plot of Figure 12, we see that only a few routes in both solutions are direct, meaning they visit an only one customer. We also observe that the LNS leads to a much larger proportion of routes with three stops than the initialization + local search. Recall the remark made on Figure 11 in the previous appendix: the demand distribution is spread over customers, so there is a need to group them in long routes. Last, the bottom plot on Figure 12 shows the distributions of the number of visits per customer in both the solution of the LNS and of the initialization + local search. We see that visits are widely spread over customers. The LNS seems to cancel some visits to customers with small numbers of visits after the initialization + local search, those visits may be related to very expensive routes.

<figure><img src="image_15.png" /><figcaption>Figure 12: Details on LNS and initialization + local search solutions of an instance.</figcaption></figure>

**Result D.1** *The strong couplings between depots, customers, days and commodities is also found in the solutions structure.*

## E LNS hyperparameters tuning.

In the LNS, we jointly tune two hyperparameters: the number of customer reinsertion steps, and the number of commodity reinsertion steps (both for one iteration of the LNS). Their values can range from zero to the number of customers or number of commodities in the instance solved. We could also tune the number of reload fixed-path vehicles steps, as well as the number of routing local search steps per outer LNS iteration. The space of configurations to consider would become too huge to be explored with reasonable computing resources. We thus focus here on the perturbations hyperparameters. We display results in Table 8: the mean, maximum and minimum gaps over the 71 instances after 90 and 300 minutes runs when varying the
hyperparameters, as well as the average number of outer LNS steps done after 90 and 300 minutes. Table 8 also displays the statistics of the initialization + local search algorithm as a reference. Overall, we see that for 90 minutes runs, increasing the number of customer reinsertion steps is crucial to achieve better performance. As shown in the ablation study in Section 7.3.3, the commodity reinsertion perturbation is less important. Yet, the best results are found when doing ten commodity reinsertion steps, and browsing the whole set of customers for the customer reinsertion subroutine per LNS outer steps. It leads to a 58% average gap in 90 minutes. For longer runs of 300 minutes, the commodity reinsertion perturbation is more important. Indeed, the worst hyperparameters configurations found are those with zero commodity reinsertion step. We thus highlight this perturbation is useful in longer runs, when escaping from local minima becomes crucial.

**Result E.1** *Our hyperparameters are tuned for short runs. Longer runs highlight the commodity reinsertion perturbation is relevant to escape from local minima.*

## F LNS cost evolution.

To extend the analysis made on the longer 300 minutes runs in Section 7.3.4, showing that the LNS has not converged in 90 minutes, we track cost with LNS iterations in long 600 minutes runs over 20 instances. We show in Figure 13 the difference between the current and final costs over subroutine iterations in the LNS. We denote by subroutine iteration the pass in one of the routing local search, reload fixed-path vehicles, customer reinsertion and commodity reinsertion subroutines (possibly making several inner iterations in each of them). Since each subroutine has a distinct behavior on each instance, the same time limit does not lead to the same number of subroutine iterations. Overall, we observe a fast cost reduction phase before 90 minutes, and a slower phase after. On average, we reduce by 4% the cost after 600 minutes compared to the 90 minutes runs.

<figure><img src="image_16.png" /><figcaption>Figure 13: Evolution of the LNS cost over 20 instances during 600 minutes runs.</figcaption></figure>

**Result F.1** *Although the LNS has not converged in 90 minutes, cost reduction is much slower on average after this point.*
<table><caption>Table 8: Gap results when tuning the numbers of customer and commodity reinsertion steps per outer LNS step on the 71 instances.</caption><thead><tr><th>Algorithm</th><th>Commodity reinsertion steps</th><th>Customer reinsertion steps</th><th>Mean gap</th><th>Min gap</th><th>Max gap</th><th>Mean number of LNS steps</th></tr></thead><tbody><tr><td>Initialization + local search</td><td>0</td><td>0</td><td>119</td><td>63</td><td>252</td><td>0</td></tr><tr><td></td><td>0</td><td>0</td><td>70</td><td>35</td><td>160</td><td>12.8</td></tr><tr><td></td><td>5</td><td>0</td><td>69</td><td>35</td><td>164</td><td>8.18</td></tr><tr><td></td><td>10</td><td>0</td><td>70</td><td>35</td><td>183</td><td>6.03</td></tr><tr><td></td><td>15</td><td>0</td><td>70</td><td>36</td><td>179</td><td>4.48</td></tr><tr><td></td><td>20</td><td>0</td><td>71</td><td>36</td><td>197</td><td>4.07</td></tr><tr><td></td><td>25</td><td>0</td><td>71</td><td>36</td><td>183</td><td>3.07</td></tr><tr><td></td><td>max</td><td>0</td><td>72</td><td>37</td><td>190</td><td>2.72</td></tr><tr><td></td><td>0</td><td>50</td><td>63</td><td>32</td><td>139</td><td>10.9</td></tr><tr><td></td><td>5</td><td>50</td><td>64</td><td>34</td><td>143</td><td>7</td></tr><tr><td></td><td>10</td><td>50</td><td>63</td><td>34</td><td>135</td><td>5.44</td></tr><tr><td></td><td>15</td><td>50</td><td>65</td><td>34</td><td>154</td><td>4.32</td></tr><tr><td></td><td>20</td><td>50</td><td>64</td><td>35</td><td>145</td><td>3.7</td></tr><tr><td></td><td>25</td><td>50</td><td>66</td><td>35</td><td>154</td><td>3.06</td></tr><tr><td></td><td>max</td><td>50</td><td>65</td><td>36</td><td>151</td><td>2.87</td></tr><tr><td></td><td>0</td><td>100</td><td>62</td><td>31</td><td>140</td><td>9.16</td></tr><tr><td></td><td>5</td><td>100</td><td>61</td><td>32</td><td>132</td><td>6.44</td></tr><tr><td></td><td>10</td><td>100</td><td>62</td><td>33</td><td>133</td><td>4.83</td></tr><tr><td></td><td>15</td><td>100</td><td>62</td><td>32</td><td>129</td><td>4.04</td></tr><tr><td></td><td>20</td><td>100</td><td>64</td><td>34</td><td>136</td><td>3.42</td></tr><tr><td></td><td>25</td><td>100</td><td>65</td><td>34</td><td>161</td><td>2.75</td></tr><tr><td></td><td>max</td><td>100</td><td>65</td><td>35</td><td>132</td><td>2.55</td></tr><tr><td></td><td>0</td><td>200</td><td>60</td><td>31</td><td>131</td><td>7.37</td></tr><tr><td></td><td>5</td><td>200</td><td>61</td><td>32</td><td>127</td><td>5.44</td></tr><tr><td>LNS</td><td>10</td><td>200</td><td>61</td><td>32</td><td>127</td><td>4.17</td></tr><tr><td>90 minutes</td><td>15</td><td>200</td><td>61</td><td>32</td><td>121</td><td>3.64</td></tr><tr><td></td><td>20</td><td>200</td><td>62</td><td>33</td><td>146</td><td>2.9</td></tr><tr><td></td><td>25</td><td>200</td><td>61</td><td>33</td><td>141</td><td>2.63</td></tr><tr><td></td><td>max</td><td>200</td><td>62</td><td>34</td><td>126</td><td>2.32</td></tr><tr><td></td><td>0</td><td>300</td><td>60</td><td>30</td><td>125</td><td>5.95</td></tr><tr><td></td><td>5</td><td>300</td><td>59</td><td>31</td><td>127</td><td>4.8</td></tr><tr><td></td><td>10</td><td>300</td><td>60</td><td>31</td><td>120</td><td>3.73</td></tr><tr><td></td><td>15</td><td>300</td><td>61</td><td>32</td><td>125</td><td>2.92</td></tr><tr><td></td><td>20</td><td>300</td><td>60</td><td>32</td><td>120</td><td>2.9</td></tr><tr><td></td><td>25</td><td>300</td><td>62</td><td>33</td><td>124</td><td>2.24</td></tr><tr><td></td><td>max</td><td>300</td><td>63</td><td>33</td><td>125</td><td>1.85</td></tr><tr><td></td><td>0</td><td>400</td><td>59</td><td>31</td><td>128</td><td>4.96</td></tr><tr><td></td><td>5</td><td>400</td><td>60</td><td>31</td><td>127</td><td>3.76</td></tr><tr><td></td><td>10</td><td>400</td><td>59</td><td>31</td><td>120</td><td>3.11</td></tr><tr><td></td><td>15</td><td>400</td><td>60</td><td>33</td><td>131</td><td>2.55</td></tr><tr><td></td><td>20</td><td>400</td><td>61</td><td>31</td><td>123</td><td>2.39</td></tr><tr><td></td><td>25</td><td>400</td><td>59</td><td>33</td><td>131</td><td>1.91</td></tr><tr><td></td><td>max</td><td>400</td><td>62</td><td>34</td><td>121</td><td>1.86</td></tr><tr><td></td><td>0</td><td>max</td><td>60</td><td>29</td><td>126</td><td>3.55</td></tr><tr><td></td><td>5</td><td>max</td><td>59</td><td>31</td><td>118</td><td>2.97</td></tr><tr><td></td><td>10</td><td>max</td><td>58</td><td>30</td><td>118</td><td>2.54</td></tr><tr><td></td><td>15</td><td>max</td><td>59</td><td>30</td><td>122</td><td>2.27</td></tr><tr><td></td><td>20</td><td>max</td><td>59</td><td>30</td><td>123</td><td>1.9</td></tr><tr><td></td><td>25</td><td>max</td><td>59</td><td>32</td><td>120</td><td>1.77</td></tr><tr><td></td><td>max</td><td>max</td><td>59</td><td>31</td><td>131</td><td>1.63</td></tr><tr><td></td><td>0</td><td>200</td><td>55</td><td>28</td><td>116</td><td>31.1</td></tr><tr><td></td><td>10</td><td>200</td><td>52</td><td>27</td><td>107</td><td>19.4</td></tr><tr><td></td><td>max</td><td>200</td><td>54</td><td>28</td><td>117</td><td>10.1</td></tr><tr><td></td><td>0</td><td>400</td><td>55</td><td>27</td><td>127</td><td>24.2</td></tr><tr><td>LNS</td><td>10</td><td>400</td><td>51</td><td>26</td><td>104</td><td>15.9</td></tr><tr><td>300 minutes</td><td>max</td><td>400</td><td>53</td><td>27</td><td>113</td><td>9.58</td></tr><tr><td></td><td>0</td><td>max</td><td>55</td><td>27</td><td>117</td><td>21.7</td></tr><tr><td></td><td>10</td><td>max</td><td>52</td><td>26</td><td>110</td><td>12.2</td></tr><tr><td></td><td>max</td><td>max</td><td>51</td><td>26</td><td>115</td><td>8.15</td></tr></tbody></table>