Matteo Petris
email: matteo.petris@inria.fr

Claudia Archetti
email: archetti@essec.edu

Diego Cattaruzza
email: diego.cattaruzza@centralelille.fr

Maxime Ogier
email: maxime.ogier@centralelille.fr

Frédéric Semet
email: frederic.semet@centralelille.fr

A heuristic with a performance guarantee for the Commodity constrained Split Delivery Vehicle Routing Problem

Keywords: Vehicle routing problems, Multiple commodities, Split delivery, Column generation, Matheuristic, Pricing heuristic

The Commodity constrained Split Delivery Vehicle Routing Problem (C-SDVRP) is a routing problem where customer demands are composed of multiple commodities. A fleet of capacitated vehicles must serve customer demands in a way that minimizes the total routing costs. Vehicles can transport any set of commodities and customers are allowed to be visited multiple times. However, the demand for a single commodity must be delivered by one vehicle only.

In this work, we developed a heuristic with a performance guarantee to solve the C-SDVRP. The proposed heuristic is based on a set covering formulation, where the exponentially-many variables correspond to routes.

First, a subset of the variables is obtained by solving the linear relaxation of the formulation by means of a column generation approach which embeds a new pricing heuristic aimed to reduce the computational time. Solving the linear relaxation gives a valid lower bound used as a performance guarantee for the heuristic. Then, we devise a restricted master heuristic to provide good upper bounds: the formulation is restricted to the subset of variables found so far and solved as an integer program with a commercial solver. A local search based on a mathematical programming operator is applied to improve the solution.

We test the heuristic algorithm on benchmark instances from the literature. Several new (best-known) solutions are found in reasonable computational time. The comparison with the state of the art heuristics for solving C-SDVRP shows that our approach significantly improves the solution time, while keeping a comparable solution quality.

Introduction

Splitting customer demands has proven to be beneficial in reducing the transportation costs and the number of vehicles [see 2 , 6]. A first work in this direction is the article by [START_REF] Dror | Savings by split delivery routing[END_REF]. The authors introduced the Split Delivery Vehicle Routing Problem (SDVRP), where customer demands are composed of a single commodity and can be split among any number of vehicles. This problem and its variants have been widely studied and exact [e.g. 9, [START_REF] Desaulniers | Branch-and-price-and-cut for the split-delivery vehicle routing problem with time windows[END_REF][START_REF] Archetti | Branch-and-cut algorithms for the split delivery vehicle routing problem[END_REF][START_REF] Munari | Compact formulations for split delivery routing problems[END_REF] and heuristic algorithms [e.g. [START_REF] Archetti | A tabu search algorithm for the split delivery vehicle routing problem[END_REF][START_REF] Silva | An iterated local search heuristic for the split delivery vehicle routing problem[END_REF][START_REF] Chen | A novel approach to solve the split delivery vehicle routing problem[END_REF][START_REF] Bortfeldt | The split delivery vehicle routing problem with three-dimensional loading constraints[END_REF] were proposed. Among these, [START_REF] Chen | The self-learning particle swarm optimization approach for routing pickup and delivery of multiple products with material handling in multiple cross-docks[END_REF] studies a particular case of the SDVRP, where customer demands are discretised a priori.

Although this delivery policy brings remarkable cost savings when compared with the policy where no splits are allowed, it is hardly applicable from a practical point of view. Indeed, customers are usually not keen to accept an unconstrained split delivery [START_REF] Archetti | Multicommodity vs. single-commodity routing[END_REF]. One step in the direction of making the split deliveries more adherent to real-world logistics is made in [START_REF] Gulczynski | The split delivery vehicle routing problem with minimum delivery amounts[END_REF]. The authors proposed a variant of the SDVRP where the quantity delivered to each customer has to be greater than a preset minimum amount.

Finally, under a multi-commodity setting, a delivery policy that might reduce customer inconvenience due to split deliveries is to allow demands to be split by commodity. Whenever a vehicle delivers a commodity to a customer, the entire quantity associated with the commodity has to be provided. This policy was firstly studied in the Discrete Split Delivery Vehicle Routing Problem (DSDVRP) proposed in [START_REF] Nakao | A DP-based heuristic algorithm for the discrete split delivery vehicle routing problem[END_REF] to deal with a real-life case study. The problem was formally introduced in the literature under the name of Commodity constrained Split Delivery Vehicle Routing Problem (C-SDVRP) in [START_REF] Archetti | Multicommodity vs. single-commodity routing[END_REF]. In the C-SDVRP, a minimum cost set of routes have to be determined such that the customer demands, composed of multiple commodities, are met, and the capacity of the vehicles is respected. The authors showed that the C-SDVRP is a relaxation of the Capacitated Vehicle Routing Problem (CVRP) where all commodities of each customer are delivered with a single vehicle, and a restriction of the SDVRP. In addition, they proposed an in-depth analysis to assess the benefits of the C-SDVRP in terms of cost savings and applicability in comparison with the CVRP and the SDVRP. To do so, they introduced the first compact mathematical formulation and devised a branch-and-cut and a heuristic algorithm to solve it. Finally, they introduced a first set of benchmark instances characterised by [START_REF] Clarke | Scheduling of vehicles from a central depot to a number of delivery points[END_REF][START_REF] Federgruen | Analysis of vehicle routing and inventory-routing problems[END_REF][START_REF] Solomon | Algorithms for the vehicle routing and scheduling problems with time window constraints[END_REF]60,80 or 100 customers requiring 2 or 3 commodities. Despite its practical relevance, the literature on the C-SDVRP and its variants is quite limited. An exact approach for the C-SDVRP was proposed in [START_REF] Archetti | A branch-price-and-cut algorithm for the commodity constrained split delivery vehicle routing problem[END_REF]. Specifically, the authors modelled the problem by means of a set covering formulation and devised a first branch-price-and-cut (BPC) algorithm. They formulate the pricing problem as an Elementary Shortest Path Problem with Resource Constraints (ESPPRC) and solve the ng-path relaxation by means of a label setting dynamic programming technique. [START_REF] Gschwind | Stabilized branch-price-and-cut for the commodity-constrained split delivery vehicle routing problem[END_REF] enhanced the performances of the BPC algorithm of [START_REF] Archetti | A branch-price-and-cut algorithm for the commodity constrained split delivery vehicle routing problem[END_REF] by embedding new procedures as the implicit bidirectional labelling search to solve the ESPPRC, the separation of non-robust valid inequalities to strengthen the lower bound, and the stabilization of the column generation procedure via dual-optimal inequalities. The authors extended the test-bed introduced in [START_REF] Archetti | Multicommodity vs. single-commodity routing[END_REF] with 336 new instances with 4, 5 and 6 commodities. The enhanced BPC algorithm outperformed the one of [START_REF] Archetti | A branch-price-and-cut algorithm for the commodity constrained split delivery vehicle routing problem[END_REF], being faster and providing several new optima and better lower bounds.

Conversely, [START_REF] Gu | Adaptive large neighborhood search for the commodity constrained split delivery VRP[END_REF] focused on a heuristic algorithm for the C-SDVRP and proposed an adaptive large neighbourhood search (ALNS) that exploits the inherent characteristics of the problem. Specifically, several existing local search moves were adapted to better deal with the multi-commodity aspect, and a mathematical programming operator was developed to reassign commodities to routes. The authors assess the performance of their ALNS on the test-bed introduced in [START_REF] Archetti | Multicommodity vs. single-commodity routing[END_REF]. The ALNS found the optimal value for 81 out of the 84 instances with 15 and 20 customers, and provided 344 new best-known solutions for the 380 instances with more than 40 customers. In [START_REF] Soleilhac | Optimisation de la distribution de marchandises avec sous-traitance du transport: une problématique chargeur[END_REF] the authors propose a small and large neighbourhood search (SLNS) which is capable of solving different variants of routing problems, among those the C-SDVRP. The SLNS is compared with the ALNS proposed in [START_REF] Gu | Adaptive large neighborhood search for the commodity constrained split delivery VRP[END_REF] on 320 instances with 100 customers and 2 and 3 commodities. The SLNS found 155 new best-known solutions, while the computational time is on average three times that of ALNS.

Finally, variants of the C-SDVRP (or of the DSDVRP) have been studied, see e.g. [START_REF] Salani | Branch and price for the vehicle routing problem with discrete split deliveries and time windows[END_REF][START_REF] Wang | A heuristic approach and a tabu search for the heterogeneous multi-type fleet vehicle routing problem with time windows and an incompatible loading constraint[END_REF][START_REF] Alinaghian | Multi-depot multi-compartment vehicle routing problem, solved by a hybrid adaptive large neighborhood search[END_REF][START_REF] Qiu | A tabu search algorithm for the vehicle routing problem with discrete split deliveries and pickups[END_REF][START_REF] Mirzaei | A branch-and-price algorithm for two multi-compartment vehicle routing problems[END_REF][START_REF] Zbib | The commodity-split multi-compartment capacitated arc routing problem[END_REF][START_REF] Guo | Heuristic approaches for a vehicle routing problem with an incompatible loading constraint and splitting deliveries by order[END_REF].

In [START_REF] Salani | Branch and price for the vehicle routing problem with discrete split deliveries and time windows[END_REF], customer demands are composed of multiple items grouped in orders. Each order can be seen as a commodity required by a customer in the C-SDVRP. In addition to the C-SDVRP, this variant includes time windows for the customers and considers service times that depend on the order delivered. The authors proposed a branch-and-price approach.

An extension of the aforementioned problem in a pickup and delivery context is the Vehicle Routing Problem with Discrete Split Deliveries and Pickups proposed by [START_REF] Qiu | A tabu search algorithm for the vehicle routing problem with discrete split deliveries and pickups[END_REF]. [START_REF] Mirzaei | A branch-and-price algorithm for two multi-compartment vehicle routing problems[END_REF] study the effect of splitting customer demands by commodity in a multi-compartment vehicle routing problem. The authors propose a branch-and-price to solve this multi-compartment C-SDVRP. [START_REF] Zbib | The commodity-split multi-compartment capacitated arc routing problem[END_REF] addressed another multi-compartment C-SDVRP in the context of a capacitated arc routing problem arising in the collection of recyclable waste. Vehicles with multiple compartments may make multiple visits to the same household to collect different recyclables, however, the amount of a single recyclable cannot be split.

In this paper, we consider a set covering formulation for the C-SDVRP, where the exponential number of variables are related to routes. Generating all such variables is intractable. Hence, we propose a restricted master heuristic [START_REF] Sadykov | Primal heuristics for branch and price: The assets of diving methods[END_REF] to solve the problem. This heuristic scheme consists in solving the formulation restricted to a subset of variables as a static integer program. Similar approaches have been successfully applied to deal with vehicle routing problems [see, e.g., [START_REF] Taillard | A heuristic column generation method for the heterogeneous fleet VRP[END_REF][START_REF] Muter | Combination of metaheuristic and exact algorithms for solving set coveringtype optimization problems[END_REF][START_REF] Parragh | Hybrid column generation and large neighborhood search for the dial-a-ride problem[END_REF][START_REF] Briant | An efficient and general approach for the joint order batching and picker routing problem[END_REF]. The main difference in the methodologies proposed in these works is the way the subset of variables is generated. [START_REF] Taillard | A heuristic column generation method for the heterogeneous fleet VRP[END_REF] and [START_REF] Muter | Combination of metaheuristic and exact algorithms for solving set coveringtype optimization problems[END_REF] developed a tabu search heuristic to populate a subset of variables for the solution of a routing problem with a heterogeneous fleet of vehicles and the Vehicle Routing Problem with Time Windows, respectively. In [START_REF] Parragh | Hybrid column generation and large neighborhood search for the dial-a-ride problem[END_REF], the authors devised a restricted master heuristic for the dial-a-ride problem: variables are generated by means of a hybrid column generation procedure where a variable neighbourhood search heuristic is employed to identify negative reduced cost columns. Finally, [START_REF] Briant | An efficient and general approach for the joint order batching and picker routing problem[END_REF] deals with the Joint Order Batching and Picker Routing Problem (JOBPRP). The authors proposed formulation with exponentially many variables and solve its linear relaxation by means of column generation. The objective is twofold: determining a subset of variables to use in a restricted master heuristic and calculating a lower bound on the optimal solution value.

The approach proposed in the current paper follows the strategy used in [START_REF] Briant | An efficient and general approach for the joint order batching and picker routing problem[END_REF]. Unlike the existing literature on the C-SDVRP, our approach is a heuristic that provides lower and upper bounds even for large-scale instances of the problem within reasonable computation times. The lower bound serves as a performance guarantee for our heuristic.

In the column generation phase, the pricing problem reduces to solve an ESPPRC. Efficient handling of the pricing problem is essential in a column generation procedure. Therefore, heuristics are commonly used to address the pricing problem before solving it exactly. In this respect, we devise a new pricing heuristic that exploits the multicommodity aspect of the problem. More precisely, the heuristic articulates in two phases: Phase 1 computes a set of promising customer sequences by solving the ESPPRC on a modified version of the pricing graph; Phase 2 is called for each customer sequence produced by the first phase and determines all the negative reduced cost routes arising from the sequence by solving the Shortest Path Problem with Resource Constraints (SPPRC) on acyclic graphs. It is noteworthy that the first phase of our heuristic also provides a valid lower bound on the value of the pricing problem. After the column generation procedure, upper bounds are identified by the restricted master heuristic. Finally, a local search phase is applied to improve the upper bound. This phase uses the mathematical programming operator proposed in [START_REF] Gu | Adaptive large neighborhood search for the commodity constrained split delivery VRP[END_REF] to reassign commodities to the routes. Computational experiments proved that our approach successfully provides upper bounds of good quality in shorter computational times than the state-of-the-art heuristic approaches. More precisely, it is capable of solving large-size instances with four, five, and six commodities and improves a few best-known solution values from the literature. When compared against the state-of-the-art heuristics of [START_REF] Gu | Adaptive large neighborhood search for the commodity constrained split delivery VRP[END_REF] and [START_REF] Soleilhac | Optimisation de la distribution de marchandises avec sous-traitance du transport: une problématique chargeur[END_REF], our approach improves the solution time with an average speedup ratio of 17.0, while keeping the percentage gap with respect to the upper bounds to 0.55% on average. The remainder of the paper is organized as follows. In Section 2, we give a formal description of the C-SDVRP and introduce the notation. In Section 3, we present a set covering formulation for the problem. In Section 4 we describe the main components of the restricted master heuristic we devised to solve it. Finally, the computational results obtained on the benchmark instances are reported and discussed in Section 5.

Problem description

In the Commodity constrained Split Delivery Vehicle Routing Problem (C-SDVRP) the commodities of a set K = {1, . . . , κ} have to be delivered from a depot 0 to a set of customers N = {1, . . . , n}. The request of a customer j ∈ N may be composed of multiple commodities and is identified by set K j = {k ∈ K : D jk > 0}, where D jk ≥ 0 is the demand of commodity k ∈ K to be delivered to customer j. An unlimited fleet performs the distribution of the commodities to the customers. Each vehicle has a capacity Q and is initially based at the depot. The vehicles can transport any subset of commodities provided that their capacity is not exceeded. We suppose without loss of generality that Q ≥ max{D jk : j ∈ N , k ∈ K j }. When a vehicle visits a customer j ∈ N , a non-empty subset M j ⊆ K j of commodities is delivered to j. Hence, a customer request may be split, and a customer may be visited multiple times. However, when a vehicle visits customer j, the amount of each commodity k ∈ M j delivered by the vehicle to j must be equal to D jk . In other words, the demand for a single commodity cannot be split.

The C-SDVRP can be defined on a directed weighted graph G = (V, A). The vertex set V = {0} ∪ N contains a vertex 0 representing the depot, and the set N of vertices representing the customers. The arc set A = {(i, j) : i, j ∈ V, i ̸ = j} contains arcs modelling each possible vehicle travel between two distinct vertices. Each arc (i, j) ∈ A is associated with a non-negative cost C ij which corresponds to the cost of traversing arc (i, j). We suppose that the arc costs satisfy the triangular inequality. A route in graph G is a non-empty circuit starting and ending at the depot. A route is feasible if the total amount of commodities delivered to the customers visited along the route does not exceed the vehicle capacity Q. The set of feasible routes is denoted by R. The cost of a route r is C r = (i,j)∈A(r) C ij , where A(r) is the set of arcs traversed by the route.

The C-SDVRP aims to find a least-cost set of feasible routes such that all the customer requests are served.

Problem formulation

We consider the set covering formulation proposed in [START_REF] Archetti | A branch-price-and-cut algorithm for the commodity constrained split delivery vehicle routing problem[END_REF]. For each feasible route r ∈ R, we introduce a binary coefficient a r jk with value one if commodity k ∈ K is delivered to customer j ∈ N by route r and zero otherwise. Then, for r ∈ R, we introduce a binary variable λ r taking value one if the route is selected in the solution and zero otherwise. Last, we define an auxiliary variable v to count the number of vehicles in the solution.

The Set Covering formulation [SC] reads as follows:

[SC] min r∈R C r λ r (1)
s.t. r∈R a r jk λ r ≥ 1 ∀j ∈ N , ∀k ∈ K j (2) r∈R λ r = v (3) v ≤ v ≤ v (4) λ r ∈ {0, 1} ∀r ∈ R (5)
Objective function (1) minimizes the total routing costs. Constraints (2), which we refer to as covering constraints, ensure that the customer demands are met. Constraints (3) and (4) define an auxiliary variable v and impose a lower bound and an upper bound on it, i.e., on the number of vehicles used in the solution. Finally, constraints (5) define variables λ r as binary.

Bounds in constraints (4) can be obtained by solving an instance of the Bin Packing Problem (BPP), where bins represent vehicles associated with their capacity, and objects, with the respective weights, are the customer demands. The BPP is formulated as an integer program and solved with a commercial solver within a time limit.

A restricted master heuristic

This section describes the main components of the restricted master heuristic we designed to tackle [SC]. The heuristic scheme articulates in three phases.

In the first phase, the Master Problem (MP), i.e., the linear relaxation of the formulation [SC] is solved using a column generation procedure ( see, e.g., [START_REF] Desrosiers | A primer in column generation[END_REF]) to obtain a subset of variables R ′ ⊆ R and a valid lower bound.

Afterwards, if the solution of the MP is fractional, valid inequalities are possibly included to strengthen the lower bound and enrich the set R ′ . The procedure is then repeated. In the second phase, an upper bound is obtained by solving formulation [SC] defined on the variables of R ′ generated in the first phase. Specifically, [SC] restricted to R ′ is solved as a static integer program with a commercial solver run within a time limit. Note that the set R ′ is preprocessed before solving [SC] to repair all the routes whose total amount of delivered commodities is not tight with respect to the vehicle capacity. For each of these routes, we randomly select commodities to be delivered to the customers they visit in order to generate new routes that fill the vehicle capacity. Finally, in the third phase, the mathematical programming operator proposed in [START_REF] Gu | Adaptive large neighborhood search for the commodity constrained split delivery VRP[END_REF] to reassign commodities to the routes is used to improve the upper bound determined in the second phase.

Column generation

As mentioned above, we developrd a column generation algorithm to solve the MP and populate a subset of variables R ′ ⊆ R. The restriction of the MP to R ′ is referred to as Restricted Master Problem (RMP). At each iteration of the procedure, the RMP and the pricing problem are solved sequentially. The pricing problem aims to either identify negative reduced cost variables (columns) to add to R ′ or to produce a certificate of optimality for the solution of the MP. We consider some heuristic approaches to quickly identify negative reduced cost variables when solving the pricing problem. Among others, we devised a novel pricing algorithm which exploits the multicommodity aspect of the C-SDVRP. When the heuristic column generators do not yield negative reduced cost variables, we solve the pricing problem using an exact algorithm to produce a certificate of optimality. In addition, it allows us to compute the Lagrangian bound, a valid lower bound on the value of [SC] which we use to provide an optimality gap for the solution of the restricted master heuristic.

Pricing problem

In this section, we use the terms path and route interchangeably.

As in [START_REF] Archetti | A branch-price-and-cut algorithm for the commodity constrained split delivery vehicle routing problem[END_REF], at each iteration of the column generation procedure, we price out routing variables λ r , r ∈ R. The reduced cost of λ r is given by:

Cr = C r - j∈N k∈Kj a r jk ρ jk -τ,
where ρ jk and τ are the dual prices associated with constraints ( 2) and (3), respectively.

[5] showed that the pricing problem reduces to an Elementary Shortest Path Problem with Resource Constraints (ESPPRC), where the resource is associated with the vehicle capacity. Following to some extent [START_REF] Gschwind | Stabilized branch-price-and-cut for the commodity-constrained split delivery vehicle routing problem[END_REF], we formulate the ESPPRC on a directed multi-graph G ′ = (V ′ , A ′ ) defined over the original graph G as follows. We include in vertex set V ′ two replica i ′ and i ′′ of each vertex i ∈ V. In arc set A ′ , we include an arc (i ′′ , j ′ ) for each arc (i, j) ∈ A to model the trip of a vehicle from vertex i to vertex j. Finally, for all non-empty subsets of commodities M j ⊆ K j , we introduce an arc (j ′ , j ′′ ) Mj to model the delivery of M j to customer j. The resource consumption D is set to zero ( Di ′′ j ′ := 0) on arcs (i ′′ , j ′ ), whereas on arcs (j ′ , j ′′ ) Mj , it is equal to the demand associated with commodity subset M j , i.e., DMj j ′ j ′′ := k∈Mj D jk . Finally, the cost on arcs

(i ′′ , j ′ ) is Ci ′′ j ′ := C ij if i ′′ ̸ = 0 ′′ and j ′ ̸ = 0 ′ , and
Ci ′′ j ′ := C ij -τ /2, otherwise. The cost on arcs (j ′ , j ′′ ) Mj considers the dual prices of constraints (2) associated with customer j ∈ N and the commodities k ∈ M j , that is CMj j ′ j ′′ := -k∈Mj ρ jk . Solving the pricing problem results in searching for a negative reduced cost elementary path in G ′ from 0 ′′ to 0 ′ such that the resource consumption does not exceed the vehicle capacity Q.

Solution of the pricing problem

Negative reduced cost paths are retrieved in the multi-graph G ′ by solving the ESPPRC by means of a label setting dynamic programming algorithm [see 21]. More precisely, labels identify partial paths in G ′ starting at 0 ′′ and are characterised by the following resources: reduced cost, accumulated demand, set of visited customers, and for each customer, the subset of commodities delivered. The starting point of the procedure is a label associated with vertex 0 ′′ with resources set to zero or empty. Then, labels are propagated from a vertex to another while satisfying the elementarity and capacity constraints: each customer is visited at most once along a partial path, and the accumulated demand cannot exceed the vehicle capacity Q. Dominance rules are applied to prune unpromising labels.

In order to accelerate the solution of the pricing problem, we implemented some state-of-the-art procedures.

Specifically, the first one is the ng-path relaxation [START_REF] Baldacci | New route relaxation and pricing strategies for the vehicle routing problem[END_REF] which partially relaxes the elementarity constraint of the paths: a neighbourhood is pre-assigned to each customer, and cycles are allowed only if the customer visited more than once in a path is not in the neighbourhood of its predecessor in that path. In addition, we incorporate an implicit version of the bidirectional labelling search algorithm proposed in [START_REF] Righini | Symmetry helps: Bounded bi-directional dynamic programming for the elementary shortest path problem with resource constraints[END_REF]. The labels are extended from vertex 0 ′′ to the other vertices of G ′ up to a value of the accumulated demand equal to Q/2. Then, the generated labels are merged to obtain complete paths [see 10, for more details].

However, even by embedding the two techniques mentioned above, the exact resolution of the pricing problem might be time consuming. We therefore proceed as follows. First, the pricing problem is solved with a new heuristic coupled with two reduced graph heuristics similar to those presented in [START_REF] Gschwind | Stabilized branch-price-and-cut for the commodity-constrained split delivery vehicle routing problem[END_REF]. Then, the same reduced graph heuristics are also applied on the multi-graph G ′ . Finally, the exact pricing method is invoked. We switch from one pricing algorithm to the next one when the first fails to produce negative reduced cost paths.

In the following, we give a detailed description the heuristics mentioned above, together with the description of the preprocessing phase to reduce the size of the multi-graph G ′ .

Preprocessing phase

We perform the preprocessing procedure proposed in [START_REF] Archetti | A branch-price-and-cut algorithm for the commodity constrained split delivery vehicle routing problem[END_REF] and [START_REF] Gschwind | Stabilized branch-price-and-cut for the commodity-constrained split delivery vehicle routing problem[END_REF] to reduce the size of the multi-graph G ′ . At each iteration of the column generation procedure, we only consider in G ′ arcs of type (j ′ , j ′′ ) Mj whose associated pair demand-cost ( DMj

j ′ j ′′ , CMj j ′ j ′′ ) is Pareto-optimal.
Since the number of commodities is small in the benchmark instances of the C-SDVRP (|K| ≤ 6), the Pareto-optimal commodity subsets can be computed by enumeration.

The reader may refer to [START_REF] Gschwind | Stabilized branch-price-and-cut for the commodity-constrained split delivery vehicle routing problem[END_REF] for a general procedure, based on the solution of the Shortest Path Problem with Resource Constraints (SPPRC) on acyclic graphs, to determine such subsets when the enumeration strategy is not applicable.

A new two-phase pricing heuristic

The heuristic we propose to solve the pricing problem consists of two phases. Phase 1 aims to compute a set of promising customer sequences. To do so, we solve the ESPPRC on a modified graph of reduced size compared with the multi-graph G ′ . Solving the ESPPRC on such a graph is not only faster than solving it on G ′ , but Phase 1 also permits to derive a valid lower bound on the pricing problem value. Phase 2 aims to determine all negative reduced cost paths arising from each of the customer sequences provided by Phase 1. We solve the SPPRC on an acyclic graph for each customer sequence. The topology and size of graph G ′ allow to perform such operation in negligible time (see [START_REF] Jaumard | A two-phase resource constrained shortest path algorithm for acyclic graphs[END_REF]).

The graph used in Phase 1, denoted by G ′′ = (V ′ , A ′′ ), differs from the multi-graph G ′ in the arcs modelling the deliveries to the customers: in G ′′ a unique subset of commodities can be delivered to each customer. Hence, G ′′ is obtained from G ′ by removing, for each customer j ∈ N , all the arcs of type (j ′ , j ′′ ) Mj but one, which we denote by (j ′ , j ′′ ). The demand and cost are set on these arcs so that whenever a customer j ∈ N is visited the least consuming commodity is delivered and all profitable dual prices related to j are collected. Hence, they are defined as Dj ′ j ′′ := min{D jk : k ∈ K j } and Cj ′ j ′′ := -k∈Kj ρ jk , respectively. This definition of demand and cost permits to derive the following properties. In Phase 2, we determine all negative reduced cost routes arising from each customer sequence (path) generated in Phase 1. We do this by solving the SPPRC on an acyclic multi-graph for each path. Specifically, let p = (j ′′ 0 = 0 ′′ , j ′ 1 , j ′′ 1 , . . . , j ′ l(p)-1 , j ′′ l(p)-1 , j ′ l(p) = 0 ′ ) be a path produced in Phase 1, where l(p) denotes the length of p. The acyclic multi-graph G ′ (p) = (V ′ (p), A ′ (p)) associated with p is defined as follows. V ′ (p) is the vertex st that includes only vertices visited along p, i.e., V

′ (p) = {j ′′ 0 = 0 ′′ , j ′ 1 , j ′′ 1 , . . . , j ′ l(p)-1 , j ′′ l(p)-1 , j ′ l(p) = 0 ′ }. A ′ (p)
contains the arcs of G ′ connecting each vertex in V ′ (p) to its successor in p, i.e. (i) arcs (j ′′ h , j ′ h+1 ), h = 0, . . . , l(p) -1 to model the travel from j h to j h+1 , and (ii) arcs (j ′ h , j ′′ h ) Mj h , h = 1, . . . , l(p) -1, M j h ⊆ K j to model the deliveries of subsets of commodities M j h to customer j h .

The negative reduced cost routes arising from path p correspond to the negative cost paths in G ′ (p) from j 0 = 0 ′′ to j l(p) = 0 ′ , which satisfy the capacity constraint. These paths are determined by solving the SPPRC on the multigraph G ′ (p). Although solving the SPPRC on acyclic graphs is NP-hard (see [START_REF] Puglia Pugliese | A survey of resource constrained shortest path problems: Exact solution approaches[END_REF]), the size and particular topology of multi-graphs G ′ (p) allow to do this operation very efficiently in terms of computational time (see [START_REF] Jaumard | A two-phase resource constrained shortest path algorithm for acyclic graphs[END_REF]).

In the following, we provide an example to illustrate how the proposed pricing heuristic works. Figure 1 shows the pricing multi-graph G ′ = (V ′ , A ′ ) arising from such instance at a certain iteration of the column generation procedure. The consumption and cost on arcs of type (i ′′ , j ′ ) ∈ A ′ modelling the movement of the vehicle from one vertex to another are ( Di ′′ j ′ , Ci ′′ j ′ ) = (0, 1). Differently, the consumption and cost on arcs of type (j ′ , j ′′ ) Mj ∈ A ′ modelling the delivery to the customers are reported in the figure with the following notation:

M j :( DMj j ′ j ′′ , CMj j ′ j ′′ ).
We only consider the Pareto-optimal deliveries to the customers. The graph G ′′ = (V ′′ , A ′′ ) built at phase I is shown in Figure 2a. The consumption and cost on arcs of type (i ′′ , j ′ ) ∈ A ′ are as in G ′ and those on arcs (j ′ , j ′′ ) ∈ A ′ are displayed in the figure with the same convention as in Figure 1. In the first phase,we solve the ESPPRC on the graph G ′′ to obtain all non-dominated negative cost paths (customer sequences) that respect the vehicle capacity. The second phase of the heuristic identifies the negative reduced cost routes arising from each of these sequences by solving the ESPPRC on acyclic multi-graphs. As an example, we show how this is done on the most negative path found in phase one (in red in Figure 2a), i.e., on p = (0 ′′ , 2 ′ , 2 ′′ , 3 ′ , 3 ′′ , 0 ′ ) with consumption and cost equal to 3 and -13, respectively. The acyclic multi-graph G ′ (p) associated with path p is shown in Figure 2b. In G ′ (p), all the possible deliveries to customers 2 and 3 are restored.

Finally, the SPPRC is solved on G ′ (p) to obtain all non-dominated feasible routes with negative reduced costs. We obtain six routes that visit customers 2 and 3 in the order imposed by path p and deliver either subset of commodities {1} or {1, 2} to customer 2, both combined with all the possible deliveries to customer 3. The route with the most negative reduced cost (in red in Figure 2b) delivers {1, 2} to customer 2 and {1, 3} to customer 3. Its consumption and reduced cost are 9 and -11, respectively.

0 ′′ 1 ′ 1 ′′ 2 ′ 2 ′′ 3 ′ 3 ′′ 0 ′ {2}:(2,

Reduced graph heuristics

In this section, we present two classical heuristics to decrease the size of the pricing graph. They are applied to the multi-graph G ′ and to the graph G ′′ of the first phase of heuristic, we just described. We discuss them for the case of G ′ , knowing that the case of G ′′ can be treated similarly.

The first heuristic is inspired by [START_REF] Toth | The granular tabu search and its application to the vehicle-routing problem[END_REF] and limits the possibilities of moving between customers. Specifically, a neighbourhood containing the g closest customers is built for each customer j ∈ N . A partial path ending in j can only be extended to customers belonging to its neighbourhood. This is implemented by removing from G ′ all arcs of type (j ′′ , l ′ ) such that l does not belong to the neighborhood of j. The pricing problem is solved considering a sequence of increasing neighbourhood sizes: g = 3, 6, 10, |N |. The value of g is incremented when the associated pricing problem produces no negative reduced cost path.

The second heuristic is specifically designed to handle the multi-commodity aspect of the C-SDVRP. Indeed, it aims at reducing the delivery possibilities to customers. Specifically, we impose an upper bound b on the number of customers whose demand can be split per path. This strategy is motivated by the analysis carried out in [START_REF] Archetti | A branch-price-and-cut algorithm for the commodity constrained split delivery vehicle routing problem[END_REF] 0 ′′ on optimal solutions of the C-SDVRP. They note that in most of them the number of split deliveries is less than three. To count the number of split customers in a path, we introduce an integer resource s in the label definition.

1 ′ 1 ′′ 2 ′ 2 ′′ 3 ′ 3 ′′
The value of s is initially set to zero and is incremented by one unit along arcs of type (j ′ , j ′′ ) Mj , if M j does not correspond to the full delivery to j (M j ⫋ K j ). On the other arcs, the value of s is simply propagated. Once s reaches the bound b, all the following customers visited in the path are delivered only with subset K j , i.e., only arcs (j ′ , j ′′ ) Kj are considered. As for the previous heuristic, we consider an incremental procedure relying on a sequence of increasing upper bounds: b = 0, 1, 2, 3, ∞.

Valid inequalities

We consider a family of robust valid inequalities, the so-called capacity cuts:

r∈R   (i,j)∈δ -(S) b r ij   λ r ≥ j∈S k∈Kj D jk Q ∀S ⊆ N , (6)
where δ -(S) = {(i, j) ∈ A : i / ∈ S, j ∈ S} is the set of arcs of graph G having their final extremity in S and b r ij is a binary coefficient taking value one if route r ∈ R traverses the arc (i, j) ∈ A.

If the solution of the MP is fractional, we separate the capacity cuts [START_REF] Archetti | Multicommodity vs. single-commodity routing[END_REF]. Since the separation of these inequalities is NP-hard, we do so by means of the heuristic algorithms presented in [START_REF] Ralphs | On the capacitated vehicle routing problem[END_REF], namely the extended shrinking heuristic and the greedy shrinking heuristic. The violated cuts are included in the RMP, and the associated dual prices π S are incorporated in the definition of the reduced cost of the variables λ r , and then considered in the pricing problem solution.

Initialization of the set R ′

We initialize the set of routes R ′ to avoid starting the column generation procedure with large dual prices, which usually slows down the pricing problem resolution. Specifically, for each customer j ∈ N , we include a round trip (0-j-0) delivering the commodities of each subset M j ⊆ K j requested by j, feasible with respect to the capacity Q. Moreover, we insert in R ′ the routes obtained by applying a variant of the Clarke-Wright algorithm (CW) [START_REF] Clarke | Scheduling of vehicles from a central depot to a number of delivery points[END_REF]. Precisely, we modified the randomized CW algorithm proposed in [START_REF] Battarra | Tuning a parametric clarke-wright heuristic via a genetic algorithm[END_REF] to take into account the multi-commodity aspect of the C-SDVRP. We set a limit of 20 runs.

Local search

In this section, we present the local search we implement to improve the C-SDVRP solution provided by the restricted master heuristic. Specifically, we consider the Mathematical Programming Operator (MPO) proposed in [START_REF] Gu | Adaptive large neighborhood search for the commodity constrained split delivery VRP[END_REF] to reassign the commodities of a specific customer j ∈ N to the routes of the solution. We iteratively call the MPO for each customer in N .

Let j ∈ N be a customer. We introduce the following notation. We denote by R-j the set of routes in the current C-SDVRP solution where all the visits to customer j are removed. More precisely, R-j contains the routes of the current C-SDVRP solution which do not visit j and, for the ones that visit j, R-j contains a copy of those routes where j has been removed from the sequence of visited customers. For each r ∈ R-j , Q j r denotes the residual capacity in route r. Finally, we indicate by C j r the cost of the cheapest insertion of customer j in route r ∈ R-j . The MPO consists in solving a Capacitated Facility Location Problem (CFLP) [see 27] where all commodities k ∈ K j of customer j have to be assigned to the routes (facilities) of R-j at minimum insertion costs and such that residual capacities of the routes are not exceeded.

The integer program on which the MPO is based makes use of the following decision variables. For each k ∈ K j and each r ∈ R-j , we introduce a binary variable

y j kr =     
1 if commodity k is delivered to customer j by route r 0 otherwise.

In addition, for each r ∈ R-j we include a binary variable z j r defined as follows

z j r =     
1 if route r delivers to customer j at least one commodity 0 otherwise.

Note that j refers to a specific customer and is not used as an index for the variables.

The integer program is for customer j as follows:

min r∈ R-j C j r z j r (7)
s.t.

r∈ R-j

y j kr = 1 ∀k ∈ K j (8) k∈Kj D jk y j kr ≤ Q j r z j r ∀r ∈ R-j (9)
y j kr ∈ {0, 1} ∀k ∈ K j , ∀r ∈ R-j (10)

z j r ∈ {0, 1} ∀r ∈ R-j (11)
The objective function [START_REF] Baldacci | New route relaxation and pricing strategies for the vehicle routing problem[END_REF] minimizes the total insertion cost. Constraints [START_REF] Battarra | Tuning a parametric clarke-wright heuristic via a genetic algorithm[END_REF] guarantee that all commodities of customer j are covered by exactly one route. Constraints [START_REF] Belenguer | A lower bound for the split delivery vehicle routing problem[END_REF] ensure that if some commodities of K j are added to a route, their demand do not exceed the remaining capacity of the route. Finally, Constraints ( 10) and ( 11) are the binary requirements.

Computational experiments

Our In the following, we denote by LB and U B, respectively, the lower and upper bounds returned by our restricted master heuristic. The percentage optimality gap is defined as 100((U B -LB)/LB). The percentage gap with respect to the best-known solution value U B bk from the literature is computed as 100((U B -U B bk )/U B bk ). More precisely, values U B bk are retrieved from [START_REF] Archetti | Multicommodity vs. single-commodity routing[END_REF][START_REF] Gu | Adaptive large neighborhood search for the commodity constrained split delivery VRP[END_REF][START_REF] Gschwind | Stabilized branch-price-and-cut for the commodity-constrained split delivery vehicle routing problem[END_REF] or [START_REF] Soleilhac | Optimisation de la distribution de marchandises avec sous-traitance du transport: une problématique chargeur[END_REF]. Finally, all the solution times are expressed in seconds.

In this section, we first describe the benchmark instances, then, we measure the impact of the new two-phase pricing heuristic on the performance of the overall solution algorithm. Later, we present the results on the whole testbed. Finally, we compare our algorithm with the two existing heuristic approaches from the literature, providing the majority of the best-known solution values.

Benchmark instances

We tested our restricted master heuristic on the benchmark instances for the C-SDVRP proposed by [START_REF] Archetti | Multicommodity vs. single-commodity routing[END_REF] and [START_REF] Gschwind | Stabilized branch-price-and-cut for the commodity-constrained split delivery vehicle routing problem[END_REF]. In each small and mid-size instance, customers' locations are taken from the C101 and R101 Solomon's instance [START_REF] Solomon | Algorithms for the vehicle routing and scheduling problems with time window constraints[END_REF], whereas, in each large instance, locations are taken from the Solomon's RC101 instances. In addition, the following parameters define the instances: (i) number of commodities |K|; (ii) probability p that a customer requires a commodity with a non-zero demand; (iii) interval ∆ to select the non-zero demand of a commodity required by a customer, expressed as a percentage of vehicle capacity; (iv) percentage α of vehicle capacity with respect to the maximum demand (Q = α max{ k∈Kj D jk : j ∈ N }). Table 1 summarises the values of these parameters which characterise each group of instances. Finally, we report a problem with the name of some of the instances in the benchmark with 100 customers, which are called C101, when the locations of the customers correspond to Solomon's RC101 instance. This problem has been corrected in the new database available at https://hal.inria.fr/hal-03836982v1. In addition, we mention that 40 instances with 100 customers with the locations of RC101, and 2 and 3 commodities are available at https://logistik.bwl.uni-mainz.de/forschung/benchmarks/, but have not been tested in any of the former papers dealing with the problem, i.e, in [START_REF] Archetti | Multicommodity vs. single-commodity routing[END_REF][START_REF] Gu | Adaptive large neighborhood search for the commodity constrained split delivery VRP[END_REF][START_REF] Gschwind | Stabilized branch-price-and-cut for the commodity-constrained split delivery vehicle routing problem[END_REF] and [START_REF] Soleilhac | Optimisation de la distribution de marchandises avec sous-traitance du transport: une problématique chargeur[END_REF]. Thus, we exclude this subset of instances when presenting the comparison with benchmark approaches.

Impact of the novel pricing heuristic

In this section, we measure the impact of the novel pricing heuristic presented in Section 4.3.2 on the performance of the overall solution algorithm. To do so, we define RMH-2P (RMH-N2P) to be the variant of the restricted master heuristic presented in Section 4 where we enable (disable) the two-phase pricing heuristic. We test the two variants on a subset of 180 instances characterised by 100 customers and 4 commodities, which is representative of the whole testbed.

The results obtained by comparing RMH-2P against RMH-N2P are shown in Table 2. We group the considered instances by the value of ∆, i.e., by the interval where the commodity demands are selected. Hence, for each group of instances, we have two columns associated with RMH-2P and RMH-N2P, respectively. The row headings of the First, both variants can provide a lower bound for all instances. Variant RMH-2P yields better results on the instances characterised by ∆ = [START_REF] Solomon | Algorithms for the vehicle routing and scheduling problems with time window constraints[END_REF]60]. Indeed, it allows us to reduce the iterations of the exact pricing algorithm by 38% and, consequently, the solution time on average by 23%. In addition, RMH-2P provides upper bounds of slightly better quality (see row avg. UB ). Differently, both variants behave similarly on the instances with ∆ = [START_REF] Alinaghian | Multi-depot multi-compartment vehicle routing problem, solved by a hybrid adaptive large neighborhood search[END_REF]100].

Although RMH-2P shows a good reduction of 24% of the iterations of the exact pricing algorithm, the reduction of the solution time of RMH-2P with respect to the one of RMH-N2P is rather limited (2%). This means that the two-phase heuristic in RMH-2P is not effective enough compared with the reduced graph heuristics and the exact pricing algorithm.

From these results, we infer that the performance of RMH-2P heavily depends on the interval of commodity demands ∆. The reason lies in how Phase 1 of the two-phase pricing heuristic is designed. The main idea of Phase 1 is to reduce the combinatorics in the solution of the ESPPRC due to the multi-commodity aspect of the C-SDVRP. Indeed, each customer is delivered with its least consuming commodity and all the profitable dual prices associated with the customer are collected. The quality of such an approximation is highly affected by the variability of the commodity demands. If ∆ = [START_REF] Solomon | Algorithms for the vehicle routing and scheduling problems with time window constraints[END_REF]60], Phase 1 provides a reasonable approximation of the benefit of serving a customer, whereas this might not be the case if ∆ = [START_REF] Alinaghian | Multi-depot multi-compartment vehicle routing problem, solved by a hybrid adaptive large neighborhood search[END_REF]100]. Indeed, in the latter case, when a customer is visited, we might collect all the profitable dual prices against a minimal consumption of the resource associated with capacity.

Based on the analysis conducted in this section, we apply variant RMH-2P to obtain results on the whole testbed.

In the following, this variant will be referred to simply as RMH.

Results on the whole testbed

In this section, we discuss the results obtained by the restricted master heuristic (RMH) on the 1380 benchmark instances. Due to the large number of instances, the results are presented in an aggregated form. The detailed instance-wise version can be found at https://hal.inria.fr/hal-03836982v1. We compare the results with the best-known solution values from the state-of-the-art exact and heuristic methods for the C-SDVRP available in the literature, namely [START_REF] Archetti | Multicommodity vs. single-commodity routing[END_REF][START_REF] Gu | Adaptive large neighborhood search for the commodity constrained split delivery VRP[END_REF][START_REF] Gschwind | Stabilized branch-price-and-cut for the commodity-constrained split delivery vehicle routing problem[END_REF] and [START_REF] Soleilhac | Optimisation de la distribution de marchandises avec sous-traitance du transport: une problématique chargeur[END_REF]. Note that in Section 5.4, we compare in more details RMH against [START_REF] Gu | Adaptive large neighborhood search for the commodity constrained split delivery VRP[END_REF] and [START_REF] Soleilhac | Optimisation de la distribution de marchandises avec sous-traitance du transport: une problématique chargeur[END_REF].

Table 3 shows the results obtained by RMH on the small and mid-size instances. Each ). However, the number of instances for which the optimality gap exceeds 5% is limited In general, RMH runs in relatively short computational times on small and mid-size instances. Moreover, Table 3 shows that the average solution time grows with the number of commodities and the number of customers.

In Table 4, we report the results obtained on the large instances (|N | = 100) which are grouped by the number of commodities |K|, the interval ∆ for customer demand and the probability p that a customer requires a commodity. The remaining columns reporting the characteristics of the instance are as in [START_REF] Archetti | Multicommodity vs. single-commodity routing[END_REF], [START_REF] Gu | Adaptive large neighborhood search for the commodity constrained split delivery VRP[END_REF] and [START_REF] Soleilhac | Optimisation de la distribution de marchandises avec sous-traitance du transport: une problématique chargeur[END_REF]. Indeed, the gap against the best-known solution values is on average 0.79% and it is zero for five instances with two commodities. Our approach provides a new solution for the 540 instances with |K| = 4, 5, 6. In addition, RMH manages to improve the value of 55 instances with |K| = 2, 3, 40 of these are new solutions. More insights about the comparison against the existing heuristic methods in the literature are drawn in Section 5.4.

We conclude this section with a discussion of the price to be paid to have a performance guarantee, i.e., to compute good lower bounds. In the case of RMH, providing lower bounds involve invoking the exact pricing algorithm in the column generation procedure and the separation of the capacity cuts (6).

To analyse the impact of the lower bound computation, we record the upper bound and the computational time before the first iteration of the exact pricing algorithm. Note that such upper bounds are obtained by solving the formulation [SC] restricted to the subset of columns found so far and applying the mathematical programming operator. Such results are reported in As expected, the price to be paid for having a performance guarantee is paid in terms of computational time: that of RMH is on average 3.9 times that of RMH-NG. The extra time spent in RMH serves not only to compute good lower bounds but also to improve the quality of the upper bounds. Indeed, several additional columns are generated in RMH and the improvement of the upper bounds of RMH compared with respect to those of RMH-NG is on average [START_REF] Gu | Adaptive large neighborhood search for the commodity constrained split delivery VRP[END_REF] and [START_REF] Soleilhac | Optimisation de la distribution de marchandises avec sous-traitance du transport: une problématique chargeur[END_REF] In this section, we analyse how RMH performs against the two existing heuristic algorithms available in the literature providing the majority of the best-known solution values, namely [START_REF] Gu | Adaptive large neighborhood search for the commodity constrained split delivery VRP[END_REF] and [START_REF] Soleilhac | Optimisation de la distribution de marchandises avec sous-traitance du transport: une problématique chargeur[END_REF]. [START_REF] Gu | Adaptive large neighborhood search for the commodity constrained split delivery VRP[END_REF] consider instances with two and three commodities, as those with more commodities were not yet available. Among the instances used in [START_REF] Gu | Adaptive large neighborhood search for the commodity constrained split delivery VRP[END_REF], [START_REF] Soleilhac | Optimisation de la distribution de marchandises avec sous-traitance du transport: une problématique chargeur[END_REF] considers only the ones with 100 customers. Hence, in this section, the testbed is restricted accordingly.

Tables 6 and7 report the comparison against [START_REF] Gu | Adaptive large neighborhood search for the commodity constrained split delivery VRP[END_REF] and [START_REF] Soleilhac | Optimisation de la distribution de marchandises avec sous-traitance du transport: une problématique chargeur[END_REF], respectively. The rows of the tables correspond to subsets of instances with the same number of customers and commodities. The first three columns report the number of customers (|N |), commodities (|K|) and instances (#inst.) in each subset. The next five columns compare the upper bounds reporting the number of instances where RMH matches (#equal ) or improves (#impr.)

the upper bounds of the competing algorithm and the average/minimal/maximal gap (gap[%] avg./min./max.).

The last three columns indicate the average solution times (avg.t[s] ) of RMH, that of the compared heuristic, and the average ratio between the solution times of the compared heuristic with respect to those of RMH (avg.ratio).

The RMH runs much faster than the other heuristics: the speedup ratio is on average 9.4 and 28.0 against [START_REF] Gu | Adaptive large neighborhood search for the commodity constrained split delivery VRP[END_REF] and [START_REF] Soleilhac | Optimisation de la distribution de marchandises avec sous-traitance du transport: une problématique chargeur[END_REF], respectively. In spite of the significant decrease in the solution time, the quality of the solutions provided by RMH remains comparable with that of the competing algorithms. Indeed, when compared with [START_REF] Gu | Adaptive large neighborhood search for the commodity constrained split delivery VRP[END_REF], RMH matches 71 upper bounds, improves 70 of them by an average of 0.26% and the percentage gap of the remaining ones is on average 0.83%. The results show a similar trend in the comparison with [START_REF] Soleilhac | Optimisation de la distribution de marchandises avec sous-traitance du transport: une problématique chargeur[END_REF]: RMH matches five upper bounds, improves 36 of them by an average of0.35% and the percentage gap of the remaining ones is on average 0.69%.

Recall that our heuristic also provides a performance guarantee contrary to [START_REF] Gu | Adaptive large neighborhood search for the commodity constrained split delivery VRP[END_REF] and [START_REF] Soleilhac | Optimisation de la distribution de marchandises avec sous-traitance du transport: une problématique chargeur[END_REF].

Conclusions

In this paper, we considered the Commodity constrained Split Delivery Vehicle Routing Problem (C-SDVRP), a routing problem where customer demands may be composed of multiple commodities, and split deliveries are allowed as long as the demand of a single commodity is delivered all at once. We presented a heuristic with a performance guarantee to solve the problem. Our heuristic is based on a column generation approach which embeds a new pricing heuristic that exploits the multi-commodity aspect of the problem. Such contribution allowed us to reduce the computational time on instances where the variability of the customer demands is not large. We performed a thorough computational analysis on the 1380 benchmark instances available in the literature. We provide an upper bound for all the considered instances, the majority of those are guaranteed to be of good quality (optimality gap less than 5%). Some new best-known solutions are found. Finally, our approach outperforms the state-of-the-art metaheuristic for the C-SDVRP ( [START_REF] Gu | Adaptive large neighborhood search for the commodity constrained split delivery VRP[END_REF] and [START_REF] Soleilhac | Optimisation de la distribution de marchandises avec sous-traitance du transport: une problématique chargeur[END_REF]) in terms of computational time while maintaining the quality of the

Values d and d

  denote the obtained lower and upper bounds, respectively. Hence, the number of vehicles is bounded from below by the ceil function of d (v := ⌈d⌉) and from above by the minimum between twice value d [see 20] and the number of customers (v := min{2 d, |N |}).

Proposition 1 .

 1 All feasible solutions of the ESPPRC on multi-graph G ′ are feasible solutions for the ESPPRC on graph G ′′ . Proposition 2. The optimal solution of the ESPPRC on graph G ′′ provides a lower bound on the optimal value of the ESPPRC on multi-graph G ′ . Corollary 1. If the optimal value of the ESPPRC on the reduced graph G ′′ is positive then the optimal value of the ESPPRC on multi-graph G ′ is positive as well.

Example 1 .

 1 We consider a C-SDVRP instance with three customers N = {1, 2, 3} and three commodities K = {1, 2, 3}: customer 1 requires the commodity of K 1 = {1}, with D 11 = 2; customer 2 requires the commodities of K 2 = {1, 2, 3}, with D 21 = 2, D 22 = 4 and D 23 = 3 ; customer 3 requires the commodities of K 3 = {1, 3}, with D 31 = 2 and D 33 = 1. We assume the travelling cost from the depot to the customers and between customers to be unitary. The vehicle capacity is set to 10.

Figure 1 :

 1 Figure 1: Pricing multi-graph G ′ for the C-SDVRP instance defined in Example 1.

  Graph G ′′ arising from the first phase of the pricing heuristic in Example 1. Acyclic multi-graph G ′ (p) where p = (0 ′′ , 2 ′ , 2 ′′ , 3 ′ , 3 ′′ , 0 ′ ) arising from the second phase of the pricing heuristic in Example 1.

Figure 2 :

 2 Figure 2: Graphs of the first and second phase of the novel pricing heuristic built in Example 1.

  The instances are divided in three groups: small (|N | = 15), mid-size (|N | = 20, 40, 60, 80) and large (|N | = 100).

  Intel(R) Xeon(R) Silver 4214 processor with 24 cores hyper-threaded to 48 virtual cores, with a base clock frequency of 2.2 GHz, and 96 GB of RAM. A time limit of one hour and a single thread are imposed for each run of the algorithm.

algorithm is implemented in C++ and compiled in release mode under a 64-bit version of MS Visual Studio 2019. CPLEX 12.9.0 (64-bit version) is used to solve the RMP in the column generation procedure and the restricted version of the formulation [SC]. All experiments are carried out on a 64-bit Windows machine equipped with a

Table 1 :

 1 Characteristics of the small, mid-size and large instances.

	Values of the parameters

table are

 are

: avg. exact pricing it.: average number of iterations of the exact pricing algorithm; avg. UB : average upper bound value; avg. t[s] : average solution time in seconds.

Table 2 :

 2 Impact of the two-phase pricing heuristic on the instances with |N | = 100 and |K| = 4.

		∆ = [40, 60]	∆ = [1, 100]
		RMH-2P RMH-N2P RMH-2P RMH-N2P
	avg. exact pricing it.	69.37	111.20	128.39	168.42
	avg. UB	3 170.03 3 175.56 2 359.43 2 353.60
	avg. t[s]	530.17	692.17	1900.58	1942.55

  : average solution time in seconds; #opt: number of optima identified by RMH with respect to the ones identified by[START_REF] Gschwind | Stabilized branch-price-and-cut for the commodity-constrained split delivery vehicle routing problem[END_REF]; #equal : number of times RMH returned the best-known solution values from the literature; #impr.: number of times RMH improved the best-known solution values from the literature (considering also the new solutions, i.e., the cases where no solution was available in the literature); avg. gap[%].: average percentage gap with respect to the best-known solution values; the character '-' indicates that no best-known solution value is available.First, note that RMH provides a lower bound for all the small and mid-size instances. For the instances with |N | = 15 and |N | = 20, the branch-price-and-cut algorithm of[START_REF] Gschwind | Stabilized branch-price-and-cut for the commodity-constrained split delivery vehicle routing problem[END_REF] provides 158 and 71 optima over 160 and 80 instances, respectively. RMH is able to identify 119 of them if |N | = 15 and 41 of them if |N | = 20. In respectively 50 and 6 of such cases, RMH proves the optimality of the obtained solutions (opt. gap = 0). The identified solutions are globally good, the optimality gap being on average equal to 0.69% and larger than 2.5% only for nine instances out of 240. Moreover, the average gap with respect to the best-known solution is 0.18%. Six new best solutions are found (see column #impr.) and the optimality gap referred to these six new best solutions is, on average, 2.09%.For the instances with |N | = 40, 60, 80, the number of commodities has a significant impact on the measure of the quality of the solutions found by RMH: the greater the number of commodities, the larger the average optimality gap (see columns opt. gap[%]

table row corresponds to a subset of instances with the same number of customers and commodities. The first four columns report some information regarding the instance subsets: |N |: number of customers; |K|: number of commodities; avg.#CC : average number of customer-commodities ( j∈N |K j |) per instance; #: number of instances in the subsets. The remaining eight columns of the table summarise the results of the RMH: #LB : number of instances for which a LB is found; opt. gap[%] avg./min./max.: average/minimal/maximal optimality gap expressed as a percentage; avg. t[s]

Table 3 :

 3 Results on the small and mid-size instances.

	Instances

Table 3 .

 3 The columns which summarise the results of the RMH are: #LB : number of instances for which a LB is found; opt. gap[%] avg./min./max.: average/minimal/maximal optimality gap expressed in percentage; avg.t[s] : average solution time in seconds; #impr.: number of times RMH improved the best-known solution values from the literature (counting also the new solutions); gap[%] avg./min./max.: average/minimal/maximal percentage gap with respect to the best-known solution values. A character '-' indicates that no solution value was previously known.First, we note that RMH manages to provide a lower bound for 874 of the 900 instances. The 26 instances where it does not succeed are characterised by |K| = 5, 6 and 25 of them also by p = 1 and ∆ = [1, 100], i.e., they have 500 or 600 non-zero customer demands. When a lower bound is found, the optimality gap mirrors the behaviour observed in Table3: it deteriorates as the number of commodities increases. More precisely, it is generally larger

when p = 1, and especially when ∆ =

[START_REF] Alinaghian | Multi-depot multi-compartment vehicle routing problem, solved by a hybrid adaptive large neighborhood search[END_REF] 100]

. In addition, the optimality gap is larger than 5% for 241 instances and, among those, 191 are characterised by five and six commodities. For the instances with fewer commodities, the average optimality gap is 2.71%, meaning that RMH offers a good performance guarantee. For the ones with |K| = 2, 3, RMH shows a behaviour comparable with the approaches from the literature identifying the best-known solution values, namely

  Table 5 in the columns with the heading RMH-NG. The rows of the table represent subsets of instances with the same number of customers. The columns are: |N |: number of customers;

#inst.: number of instances; avg.UB : average upper bounds of RMH-NG and RMH; avg.gap[%] : average percentage gap between the upper bounds of RMH with respect to the ones of RMH-NG; avg.t[s] : average solution time in seconds of RMH-NG and RMH; avg.ratio: ratio between the solution time of RMH with respect to the one of RMH-NG;

Table 4 :

 4 Results on the large instances.

									RMH results
			Instances							best known
							opt. gap[%]				gap[%]
	|K|	∆	p	avg.#CC #inst. #LB avg. min. max. avg.t[s] #impr. avg. min. max.
			0.6	136.40	40	40	1.36 0.03	3.39	30.49	0	0.43	0.00	1.74
		[40,60]
			1	200.00	50	50	1.39 0.09	3.12	68.28	11	0.39 -0.05	1.10
	2
			0.6	136.40	40	40	2.40 1.06	5.24	183.20	0	0.95	0.04	2.87
		[1,100]
			1	200.00	50	50	1.94 0.54	3.95	250.62	16	0.48 -0.23	2.27
			0.6	188.40	40	40	1.77 0.70	3.22	155.42	1	0.45 -0.15	1.79
		[40,60]
			1	300.00	50	50	2.46 0.74	7.22	247.15	12	1.09 -0.60	5.17
	3
			0.6	188.40	40	40	2.46 0.67	5.54	459.32	2	0.89 -0.11	3.49
		[1,100]
			1	300.00	50	50	3.51 1.50	6.71	995.53	13	1.62 -0.45	4.48
			0.6	242.95	40	40	2.42 1.00	5.00	379.13	40	-	-	-
		[40,60]
			1	400.00	50	50	3.91 0.93	9.11	651.01	50	-	-	-
	4
			0.6	243.30	40	40	3.29 0.80	7.33	1100.49	40	-	-	-
		[1,100]
			1	400.00	50	50	5.10 2.05	9.61	2540.64	50	-	-	-
			0.6	299.13	40	40	3.66 1.65	6.19	957.39	40	-	-	-
		[40,60]
			1	500.00	50	50	5.19 2.54 10.08 1609.24	50	-	-	-
	5
			0.6	301.48	40	40	4.82 1.66 13.01 2196.31	40	-	-	-
		[1,100]
			1	500.00	50	45	5.85 2.99 11.50 3479.58	50	-	-	-
			0.6	359.65	40	40	4.78 2.03	8.05	1764.76	40	-	-	-
		[40,60]
			1	600.00	50	50	6.43 2.95 11.98 2641.67	50	-	-	-
	6
			0.6	357.73	40	39	6.04 3.14 15.24 2972.40	40	-	-	-
		[1,100]
			1	600.00	50	30	7.37 4.99 11.89 3622.95	50	-	-	-

Table 5 :

 5 Impact of the performance guarantee.

	Instances			RMH-NG vs. RMH
				avg. UB			avg.t[s]
	|N | #inst. RMH-NG	RMH	avg. gap[%] RMH-NG	RMH	avg. ratio
	15	160	389.35	386.69	-0.65	9.92	32.69	3.93
	20	80	703.78	693.76	-1.41	41.47	160.65	4.07
	40	80	1152.61 1139.68	-1.14	135.62	669.74	5.67
	60	80	1659.90 1644.11	-0.96	254.19 1260.69	4.72
	80	80	2126.85 2113.77	-0.66	361.21 1678.23	4.45
	100	900	2793.71 2780.10	-0.38	382.29 1348.10	3.64
	5.4. Comparison with

Table 6 :

 6 Comparison with[START_REF] Gu | Adaptive large neighborhood search for the commodity constrained split delivery VRP[END_REF] on the instances with |K| = 2, 3 and customers' locations from C101 and R101.

							RMH vs. [23]
		Instances			UB
							gap[%]			avg.t[s]
	|N | |K| #inst. #equal #impr. avg. min. max.	RMH	[23]	avg.ratio
	15	2	32	30	0	0.02	0.00	0.28	0.32	8.91	35.12
	15	3	32	22	0	0.17	0.00	2.47	3.98	15.63	21.05
	20	3	20	12	0	0.18	0.00	0.98	4.52	56.78	14.28
	40	3	20	2	2	0.36 -0.09	1.01	58.76	117.14	2.72
	60	3	20	0	0	0.70	0.05	2.15	170.76 281.83	1.94
	80	3	20	0	2	0.98 -0.26	2.83	354.51 511.00	1.47
		2	160	5	42	0.34 -0.67	2.66	136.07 445.09	9.46
	100
		3	160	0	24	0.91 -0.96	5.17	481.80 822.60	4.02

Table 7 :

 7 Comparison with[START_REF] Soleilhac | Optimisation de la distribution de marchandises avec sous-traitance du transport: une problématique chargeur[END_REF] on the instances with |N | = 100 and |K| = 2, 3 and customers' locations from C101 and R101.

					RMH vs. [39]
	Instances			UB
					gap[%]			avg.t[s]
	|N | |K| #inst. #equal #impr. avg. min. max.	RMH	[39]	avg.ratio
	2	160	5	12	0.49 -0.34	2.84	136.07 1800	44.62
	100
	3	160	0	24	0.63 -1.45	4.03	481.80 1800	11.47

The gap with respect to the best-known solution values is on average equal to 0.70% and larger than 2.5% only in

0.6%. These results show that RMH can be easily adapted according to the needs: in the case where a solution is needed in a short computing time, the RMH-NG can be used without sacrificing solution quality to a large extent.

upper bounds comparable.

Future research may be devoted to the inclusion of additional families of valid inequalities (possibly robust) to improve both the lower and upper bounds. Finally, our approach may be adapted to solve variants of the problem where the additional constraints can be easily handled in the pricing problem (e.g. a multi-compartment C-SDVRP).