A Branch-Price-and-Cut algorithm for the Multi-Commodity two-echelon Distribution Problem

Petris, Matteo; Archetti, Claudia; Cattaruzza, Diego; Ogier, Maxime; Semet, Frédéric Article A Branch-Price-and-Cut algorithm for the MultiCommodity two-echelon Distribution Problem EURO Journal on Transportation and Logistics (EJTL) Provided in Cooperation with: Association of European Operational Research Societies (EURO), Fribourg Suggested Citation: Petris, Matteo; Archetti, Claudia; Cattaruzza, Diego; Ogier, Maxime; Semet, Frédéric (2024) : A Branch-Price-and-Cut algorithm for the Multi-Commodity two-echelon Distribution Problem, EURO Journal on Transportation and Logistics (EJTL), ISSN 2192-4384, Elsevier, Amsterdam, Vol. 13, Iss. 1, pp. 1-13, https://doi.org/10.1016/j.ejtl.2024.100139 This Version is available at: https://hdl.handle.net/10419/325210 Standard-Nutzungsbedingungen: Die Dokumente auf EconStor dürfen zu eigenen wissenschaftlichen Zwecken und zum Privatgebrauch gespeichert und kopiert werden. Sie dürfen die Dokumente nicht für öffentliche oder kommerzielle Zwecke vervielfältigen, öffentlich ausstellen, öffentlich zugänglich machen, vertreiben oder anderweitig nutzen. Sofern die Verfasser die Dokumente unter Open-Content-Lizenzen (insbesondere CC-Lizenzen) zur Verfügung gestellt haben sollten, gelten abweichend von diesen Nutzungsbedingungen die in der dort genannten Lizenz gewährten Nutzungsrechte. Terms of use: Documents in EconStor may be saved and copied for your personal and scholarly purposes. You are not to copy documents for public or commercial purposes, to exhibit the documents publicly, to make them publicly available on the internet, or to distribute or otherwise use the documents in public. If the documents have been made available under an Open Content Licence (especially Creative Commons Licences), you may exercise further usage rights as specified in the indicated licence. https://creativecommons.org/licenses/by/4.0/ Contents lists available at ScienceDirect EURO Journal on Transportation and Logistics journal homepage: www.elsevier.com/locate/ejtl A Branch-Price-and-Cut algorithm for the Multi-Commodity two-echelon Distribution Problem Matteo Petris a,∗, Claudia Archetti b, Diego Cattaruzza c, Maxime Ogierc, Frédéric Semet c aUniv. Lille, Inria, CNRS, Centrale Lille, UMR 9189 CRIStAL, F-59000 Lille, France bDepartment of Information Systems, Decision Sciences and Statistics, ESSEC Business School, Cergy-Pontoise, France cUniv. Lille, CNRS, Inria, Centrale Lille, UMR 9189 CRIStAL, F-59000 Lille, France ARTICLE INFO Keywords: Two echelon routing problems Multiple commodities Split delivery Branch-Price-and-Cut ABSTRACT In the Multi-Commodity two-echelon Distribution Problem (MC2DP), multiple commodities are distributed in a two-echelon distribution system involving suppliers, distribution centres and customers. Each supplier may provide different commodities and each customer may request several commodities as well. In the first echelon, capacitated vehicles perform direct trips to transport the commodities from the suppliers to the distribution centres for consolidation purposes. In the second echelon, each distribution centre owns a fleet of capacitated vehicles to deliver the commodities to the customers through multi-stop routes. Commodities are compatible, i.e., they can be mixed in the vehicles. Finally, customer requests can be split by commodities, that is, a customer can be visited by several vehicles, but the total amount of each commodity has to be delivered by a single vehicle. The aim of the MC2DP is to minimize the total transportation cost to satisfy customer demands. We propose a set covering formulation for the MC2DP where the exponential number of variables relates to the routes in the delivery echelon. We develop a Branch-Price-and-Cut algorithm (BPC) to solve the problem. The pricing problem results in solving an Elementary Shortest Path Problem with Resource Constraints (ESPPRC) per distribution centre. We tackle the ESPPRC with a label setting dynamic programming algorithm which incorporates ng-path relaxation and a bidirectional labelling search. Pricing heuristics are invoked to speed up the procedure. In addition, the formulation is strengthened by integrating capacity cuts and two families of valid inequalities specific for the multiple commodities aspect of the problem. Our approach solves to optimality 439 over the 736 benchmark instances from the literature. The optimality gap of the unsolved instances is 2.1%, on average. 1. Introduction In a two-echelon distribution system, goods are transferred from origins (depots, suppliers) to destinations (customers) via intermediate facilities (satellites, distribution centres) (see Guastaroba et al., 2016). In the collection echelon, large vehicles bring goods from the origins to the intermediate facilities where consolidation operations are performed. Whereas, in the delivery echelon, smaller vehicles are in charge of distributing the goods to the final customers. Routing decisions are usually required in both echelons. Two-echelon systems take advantage of consolidating goods at intermediate facilities and using different fleets within each echelon to reduce overall transportation costs. An example of this delivery strategy can be encountered in city logistics (Cattaruzza et al.,2017;Crainic et al.,2023) where the aim is also to grant the access in urban areas only to environmental-friendly vehicles that usually have a small capacity. ∗Corresponding author. E-mail addresses: [email protected] (M. Petris), [email protected] (C. Archetti), [email protected] (D. Cattaruzza), [email protected] (M. Ogier), [email protected] (F. Semet). In this article, we consider the two-echelon distribution problem introduced in Gu et al. (2022), namely the Multi-Commodity two-echelon Distribution Problem (MC2DP). In this context, origins, intermediate facilities and destinations are referred to as suppliers,distribution centres and customers, respectively. There are few vehicle routing problems which explicitly deal with multiple commodities within a two-echelon distribution system. To the best of our knowledge, among these problems, the MC2DP is the only one considering a many-to-many setting. In fact, in the MC2DP, the commodity requested by a customer is not pre-assigned to a specific supplier, so it can be collected at any supplier or subset of suppliers where it is available. The amount of the commodities available at the suppliers is limited. In contrast with the usual setting in the literature, the MC2DP requires routing decisions only in the delivery echelon. Indeed, commodities are collected from https://doi.org/10.1016/j.ejtl.2024.100139 Received 31 August 2023; Received in revised form 5 July 2024; Accepted 11 July 2024 EURO Journal on Transportation and Logistics 13 (2024) 100139 Available online 14 July 2024 2192-4376/© 2024 The Author(s). Published by Elsevier B.V. on behalf of Association of European Operational Research Societies (EURO). This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ ). M. Petris et al. the suppliers and brought to the distribution centres via direct round trips. In the delivery echelon, a fleet of vehicles performing routes starting and ending at the same distribution centre is used to deliver the commodities to the customers. All vehicles involved in the distribution system are capacitated and commodities are compatible, i.e., they can be mixed inside all vehicles. Finally, as in the Commodity constrained Split Delivery Vehicle Routing Problem (C-SDVRP) (Archetti et al.,2016), customers can be visited by multiple vehicles as long as the demand of a single commodity is served by a single vehicle. The aim of the MC2DP is to determine a distribution plan to satisfy customer demands while respecting the capacity of the vehicles and not exceeding the commodity availabilities at the suppliers and such that the total transportation cost is minimized. The MC2DP finds an application in the short and local fresh food supply chains (Berti and Mulligan, 2016) where farmers supply different agricultural products to canteens, restaurants or supermarkets through indirect sales. Commonly, a single decision maker, such as an association of farmers, coordinates both the collection and delivery echelons. In this context, the farmers are less numerous than delivery points since the maximal supply of one farmer can cover the demand of several customers. Hence, the collection from the farmers is usually performed via direct round trips. Then, the distribution centres perform the consolidation operations and the deliveries to the customers which are done by vehicles performing routes. We refer to Gu et al. (2022) for more details about the problem application. The authors in Gu et al. (2022) proposed a compact Mixed Integer Linear Programming (MILP) formulation and a sequential heuristic for the MC2DP. The authors decompose the MC2DP in two subproblems: one for the collection from suppliers, and the other one for the delivery to customers. The collection subproblem is modelled as a MILP and solved with a commercial solver while the delivery subproblem is solved by an Adaptive Large Neighbourhood Search (ALNS) algorithm. The contribution of this paper is to present an exended model and to propose the first ever exact approach based on a Branch-Price-andCut (BPC) algorithm to solve the MC2DP. Similar exact approaches have recently been proposed to deal with two-echelon vehicle routing problems (see e.g. Marques et al.,2020;Li et al.,2022;Mhamedi et al., 2022;Marques et al.,2022). However, our BPC algorithm is designed to take into account explicitly the multi-commodity dimension. Specifically, our algorithm relies on a set covering formulation for the MC2DP where the exponentiallymany number of variables correspond to the routes in the delivery echelon starting and ending at each distribution centre. We also strengthen the formulation by the insertion of capacity cuts, valid inequalities arising from the set covering polytope (Balas and Ng,1989) and a new family of valid inequalities based on the number partitioning problem polytope. While capacity cuts are classical inequalities derived for the Capacitated Vehicle Routing Problem (CVRP) (see Laporte et al.,1985), the other two families of inequalities tackle the multi-commodity aspect of the problem. Finally, several state-of-art speed-up techniques are also incorporated in our BPC algorithm. The remainder of the paper is organized as follows. Section 2 provides a literature review. In Section 3, a formal description of the MC2DP is provided. In Section 4, a set covering formulation is presented along with different families of valid inequalities. Our BranchPrice-and-Cut algorithm is described in Section 5. Finally, in Section 6we analyse the results obtained by the proposed algorithm on the benchmark instances introduced in Gu et al. (2022) to assess its effectiveness. 2. Literature review In this section, we review the existing literature on the two-echelon distribution problems, with particular attention to the ones dealing with multiple commodities. The first two-echelon routing problem introduced by Jacobsen and Madsen (1980) was motivated by a specific application. Newspapers have to be distributed from a printing office to sales points possibly passing through some transfer points whose locations are to be decided. Crainic et al. (2004,2009) proposed a formal description of a rich class of two-echelon routing problems along with some economic insights. However, the seminal problem in this class, namely the two-echelon Capacitated Vehicle Routing Problem (2E-CVRP), was introduced in the literature and studied for the first time in Perboli et al. (2011). In the 2E-CVRP, a single commodity has to be transferred from a single origin to several destinations through some intermediate facilities. Two fleets of capacitated vehicles perform routes in the two echelons to transport the commodity from the origin to the intermediate facilities and from the intermediate facilities to the destinations. The objective of the 2E-CVRP is to minimize the total transportation cost of the distribution system. The authors proposed two math-heuristics to solve the problem, a diving and a sub-MIP heuristic. The 2E-CVRP and related problems have received increasing attention in recent years and many variants have been addressed, e.g., 2ECVRP with (i) time windows (Mhamedi et al.,2022); (ii) mobile satellites (Li et al.,2020); (iii) synchronization (Grangier et al.,2016) and bi-synchronization (Li et al.,2021b); (iv) simultaneous pickup and delivery (Li et al.,2022); (v) electric vehicles (Breunig et al.,2019) and battery swapping stations (Jie et al.,2019); (vi) real-time transshipment capacity varying (Li et al.,2018); (vii) covering options (Enthoven et al.,2020); (viii) delivery options (Zhou et al.,2018); (ix) stochastic demands (Sluijk et al.,2022). The interested reader may refer to Cuda et al. (2015), Li et al. (2021a) and Sluijk et al. (2023) for recent surveys on the subject. According to the existing literature, the vast majority of the twoechelon routing problems deal with the single commodity case. Apart from the MC2DP, which is addressed in this paper, only a few works integrate multiple commodities in a two-echelon routing problem (e.g. Dellaert et al.,2021;Jia et al.,2023;Gu et al.,2022). In Dellaert et al. (2021), the authors extended the 2E-CVRP by introducing multiple origins and multiple commodities. In addition, hard time windows are imposed for the delivery at the destinations. In their problem, customers have a commodity demand from a specific origin, i.e., there is a one-to-one setting. Several mathematical formulations are proposed and a BPC algorithm is devised to solve the problem. In Jia et al. (2023), the problem setting is similar to the one of Dellaert et al. (2021). However, the multi-commodity aspect is handled with more restrictions: only two origins are considered and each destination requires one commodity per origin (one-to-one setting). The authors developed an ALNS algorithm to solve large-scale instances of the problem. The MC2DP introduced in Gu et al. (2022) differs from Dellaert et al. (2021) and Jia et al. (2023) for three reasons: (i) there is a manyto-many setting for the commodities, i.e. any commodity requested by a customer can be served from any supplier; (ii) suppliers provide commodities in limited amounts; (iii) routing decisions are not required in the collection echelon. 3. Problem description In the Multi-Commodity two-Echelon Distribution Problem (MC2DP), a set of commodities is distributed in a system involving a set of suppliers (origins) , a set of distribution centres (intermediate facilities) and a set of customers (destinations) . The system is split in two echelons: the collection echelon where the commodities are collected at the suppliers and brought to the distribution centres, and the delivery echelon where the commodities at the distribution centres are delivered to the customers. More precisely, in the collection echelon, each supplier 𝑖∈provides a maximal amount 𝑃𝑖𝑘 ≥0for each commodity 𝑘∈. Note that a supplier 𝑖∈might not supply a commodity 𝑘∈, and in that case, 𝑃𝑖𝑘 takes value 0. An unlimited fleet of homogeneous vehicles of capacity 𝑄𝑆performs direct round trips from the distribution centres to collect the commodities from the EURO Journal on Transportation and Logistics 13 (2024) 100139 2 M. Petris et al. suppliers. The vehicles can transport any subset of commodities. Due to the limited capacity of the vehicles, direct round trips between a distribution centre 𝑜∈and a supplier 𝑖∈may be performed by several vehicles. The problem associated with the collection operations can be modelled as a Multi-commodity Capacitated fixed-charge Network Design Problem (MCNDP, Magnanti and Wong,1984) with a specific cost structure: there is a step-wise cost function defined by a unitary cost associated with each vehicle used between a distribution centre and a supplier. Differently, the problem of distributing the commodities from the distribution centres to the customers is a multi-depot version of the Commodity constrained Split Delivery Vehicle Routing Problem (C-SDVRP). Each customer 𝑗∈has a demand 𝐷𝑗𝑘 ≥0for all commodities 𝑘∈. The request of customer 𝑗is identified by set 𝑗= {𝑘∈ ∶𝐷𝑗𝑘 >0}. Each distribution centre owns an unlimited fleet of homogeneous and capacitated vehicles of capacity 𝑄𝐷which performs routes to deliver the commodities to the customers. Each vehicle has to end its route at its starting distribution centre. As in the collection echelon, a vehicle can be loaded with any commodities. Without loss of generality, we suppose 𝑄𝐷≥max{∑𝑘∈𝑗𝐷𝑗𝑘 ∶𝑗∈}. Furthermore, customer requests can be split, i.e., different vehicles can serve the same customer. However, the demand of a single commodity cannot be split: it has to be delivered by a single vehicle. Note that direct trips from suppliers to customers and inter-connections between distribution centres are not allowed. Finally, the collection and delivery operations taking place in the two echelons are coordinated at the distribution centres by means of the so-called load synchronization strategy (Drexl,2012): the total amount of each commodity collected at the suppliers by each distribution centre must be sufficient to serve the customer demands of that commodity delivered by a vehicle of that distribution centre. We formulate the MC2DP on a directed weighed graph = (,). Set =∪∪contains a vertex for each supplier, distribution centre and customer. Arc set =𝑆∪𝐷is defined as the union of two sets of arcs which model the possible vehicle travels in the two echelons. Specifically, set 𝑆= (×) ∪ (×)includes the arcs modelling the direct trips from suppliers to distribution centres in the collection echelon, whereas 𝐷= (×)∪(×)∪(×)contains all arcs between customers and between distribution centres and customers. Each arc (𝑖, 𝑗) ∈ is assigned with a non-negative cost 𝐶𝑖𝑗 which represent the transportation cost of a vehicle traversing (𝑖, 𝑗). The arc costs are symmetric and satisfy the triangular inequality. In graph , a route in the delivery echelon is a non-empty circuit starting and ending at a distribution centre 𝑜∈. A route is feasible if the total amount of commodities delivered to the customers visited along the route does not exceed vehicle capacity 𝑄𝐷. The cost of any feasible route 𝑟is 𝐶𝑟=∑(𝑖,𝑗)∈(𝑟)𝐶𝑖𝑗 , where (𝑟)is the set of arcs traversed by the route. Finally, the total transportation cost of the distribution system arising from the MC2DP is the sum of the cost of the direct round trips in the collection echelon and the routing costs in the delivery echelon. The aim of the MC2DP is to determine a distribution plan, i.e., the direct round trips in the collection echelon and the routes in the delivery echelon, which satisfies the customer requests, does not exceed the commodity availabilities at the suppliers, satisfies the vehicle capacities in both echelons and respects the load synchronization constraints while minimizing the total transportation cost. 4. Problem formulation We model the MC2DP by means of a set covering formulation, where the exponentially-many variables are associated with the routes in the delivery echelon. For each distribution centre 𝑜∈, we define 𝑜as the set of all feasible routes starting and ending at 𝑜. The set of all feasible routes is denoted by =⋃𝑜∈𝑜. We define a binary coefficient 𝑎𝑟 𝑗𝑘 with value one if commodity 𝑘∈is delivered to customer 𝑗∈by route 𝑟∈ and zero otherwise. For each supplier 𝑖∈and each distribution centre 𝑜∈, we introduce an integer variable 𝑥𝑖𝑜 to represent the number of vehicles traversing arc (𝑖, 𝑜) ∈ . For each 𝑖∈,𝑜∈and 𝑘∈, we define a non-negative continuous variable 𝑞𝑘 𝑖𝑜 that represents the amount of commodity 𝑘collected at supplier 𝑖by distribution centre 𝑜. Finally, for each route 𝑟∈, we introduce a binary variable 𝜆𝑟taking value one if 𝑟is selected in the solution and zero otherwise. The Set Covering formulation [SC] for the MC2DP reads as follows: [SC] min ∑ (𝑖,𝑜)∈ 2𝐶𝑖𝑜𝑥𝑖𝑜 +∑ 𝑟∈ 𝐶𝑟𝜆𝑟(1) s.t. ∑ 𝑜∈ 𝑞𝑘 𝑖𝑜 ≤𝑃𝑖𝑘 ∀𝑖∈,∀𝑘∈(2) ∑ 𝑘∈ 𝑞𝑘 𝑖𝑜 ≤𝑄𝑆𝑥𝑖𝑜 ∀𝑖∈,∀𝑜∈(3) ∑ 𝑟∈ 𝑎𝑟 𝑗𝑘𝜆𝑟≥1 ∀𝑗∈,∀𝑘∈𝑗(4) ∑ 𝑖∈ 𝑞𝑘 𝑖𝑜 ≥∑ 𝑟∈𝑜∑ 𝑗∈ 𝑎𝑟 𝑗𝑘𝐷𝑗𝑘𝜆𝑟∀𝑜∈,∀𝑘∈(5) 𝑥𝑖𝑜 ∈Z≥0∀𝑖∈,∀𝑜∈(6) 𝑞𝑘 𝑖𝑜 ∈R≥0∀𝑖∈,∀𝑜∈,∀𝑘∈(7) 𝜆𝑟∈ {0,1} ∀𝑟∈(8) Objective function (1) minimizes the total transportation cost. Constraints (2) ensure that the commodity availabilities at each supplier are respected. Constraints (3) guarantee that a sufficient number of vehicles perform the collection operations and that the capacity of these vehicles is not exceeded. Covering Constraints (4) impose that each commodity required by a customer is served by at least one route. In addition, the load synchronization constraint linking the collection and delivery echelons is expressed in constraints (5): the quantity of each commodity collected by each distribution centre has to be large enough to satisfy the demand for that commodity delivered by a route of that distribution centre. Finally, Constraints (6),(7) and (8) define variable domains. 4.1. Valid inequalities In this section, we introduce four families of valid inequalities considered to strengthen formulation [SC]. Two of these inequalities are known in the context of vehicle routing problems, while the other two are tailored to deal with the multi-commodity aspect of the MC2DP. Note that such inequalities are valid for the C-SDVRP, hence for the MC2DP. In what follows, given a subset of customers ′⊆, we define 𝐷(′) = ∑𝑗∈′∑𝑘∈𝑗𝐷𝑗𝑘 to be the total demand requested by the customers in ′. In addition, we introduce a binary coefficient 𝑏𝑟 𝑖𝑗 with value one if route 𝑟∈traverses arc (𝑖, 𝑗) ∈ and zero otherwise. Finally, we define 𝑒𝑟 𝑗=∏𝑘∈𝑎𝑟 𝑗𝑘 to be a binary coefficient equal to one if route 𝑟delivers all the commodities of subset ⊆𝑗to customer 𝑗∈and zero otherwise. Bounds on the number of vehicles The following inequalities set bounds on the number of vehicles in the collection and delivery echelons: ∑ (𝑖,𝑜)∈ 𝑥𝑖𝑜 ≥⌈𝐷() 𝑄𝑆⌉(9) and ∑ 𝑟∈ 𝜆𝑟≥⌈𝑣⌉(10a) ∑ 𝑟∈ 𝜆𝑟≤min{||,2𝑣}.(10b) EURO Journal on Transportation and Logistics 13 (2024) 100139 3 M. Petris et al. In inequalities (10a) and (10b), values 𝑣and 𝑣 are obtained by solving an instance of the Bin Packing Problem (BPP), where bins have size equal to the vehicle capacity 𝑄𝐷, and each customer demand has a corresponding item to be packed with size 𝐷𝑗𝑘. Precisely, we solve an integer program for the BPP on such an instance with a commercial solver within a short time limit: 𝑣and 𝑣 are the obtained lower and upper bounds. If the instance is solved to optimality within the time limit, 𝑣=𝑣 holds. The right hand-side of (10b) is the minimum between twice value 𝑣 (see Federgruen and Simchi-Levi,1995) and the number of customers. Capacity cuts Laporte et al. (1985) introduced the capacity cuts to deal with the Capacitated Vehicle Routing Problem: ∑ 𝑟∈(∑ (𝑖,𝑗)∈𝛿−(′) 𝑏𝑟 𝑖𝑗 )𝜆𝑟≥⌈𝐷(′) 𝑄𝐷⌉∀′⊆,(11) where 𝛿−(′) = {(𝑖, 𝑗) ∈ ∶𝑖∉′, 𝑗 ∈′}is the set of arcs of graph  reaching a vertex in ′. Given a subset of customers ′, inequality (11) states that at least ⌈𝐷(′)∕𝑄𝐷⌉vehicles of the delivery echelon are required to cover the requests of the customers in ′. Set covering polytope We present a family of valid inequalities inspired by the facetdefining inequalities proposed in Balas and Ng (1989) for the set covering polytope. Although these inequalities were proposed several years ago, to the best of our knowledge, they have not yet been used in BPC algorithms for vehicle routing problems. However, they are similar to the strong minimum number of vehicles inequalities introduced by Archetti et al. (2011) in the context of a BPC algorithm for the split delivery vehicle routing problem with time windows. Let us first briefly present a formulation for the set covering problem. Let bet a set of elements to be covered, and be a set of subsets of . We denote by 𝑐𝑗the cost associated to subset 𝑗∈, and 𝑑𝑖𝑗 a binary parameter that takes value one if element 𝑖∈is in subset 𝑗∈, and zero otherwise. Let 𝑥𝑗be a binary decision variable taking value one if subset 𝑗∈is selected, zero otherwise. An integer programming formulation for the set covering problem is min ∑ 𝑗∈ 𝑐𝑗𝑥𝑗 s.t. ∑ 𝑗∈ 𝑑𝑖𝑗 𝑥𝑗≥1 ∀𝑖∈ 𝑥𝑗∈ {0,1} ∀𝑗∈ Given a subset ′⊆, the inequalities introduced in Balas and Ng (1989) reads as follows: 2∑ 𝑗∈′ 𝑥𝑗+∑ 𝑗∈ ′ 𝑥𝑗≥2, where ′= {𝑗∈∶𝑑𝑖𝑗 = 1,∀𝑖∈′}is the set of the elements of  which cover ′and  ′= {𝑗∈∶∑𝑖∈′𝑑𝑖𝑗 ≥1 ∧ ∏𝑖∈′𝑑𝑖𝑗 = 0} is the set of the elements of which contain some, but not all, the elements in ′. The inequalities express how subset ′may be covered: either it suffices to select a unique element in that covers ′, i.e., an element in ′, or at least two elements in that partially cover ′have to be selected, i.e., at least two elements in  ′. Under specific conditions, these constraints are facet defining for the set covering polytope. In what follows, we adapt these inequalities to the MC2DP to express how the subsets of commodities required by a given customer may be covered. For the ease of readability, we introduce the following notation. Let 𝑗∈be a customer and 𝑗⊆𝑗be a subset of the commodities requested by 𝑗. We denote by 𝑗 𝑗⊆the subset of routes delivering all commodities in 𝑗to 𝑗, i.e., 𝑗 𝑗= {𝑟∈∶ 𝑒𝑟 𝑗𝑗= 1}. In addition, we write  𝑗 𝑗⊆for the subset of routes which deliver some of the commodities in 𝑗to 𝑗, but not all of them, i.e.,  𝑗 𝑗= {𝑟∈∶∑𝑘∈𝑗𝑎𝑟 𝑗𝑘 ≥1 ∧ 𝑒𝑟 𝑗𝑗= 0}. The set covering polytope inequalities for the MC2DP are defined as follows: 2∑ 𝑟∈𝑗 𝑗 𝜆𝑟+∑ 𝑟∈ 𝑗 𝑗 𝜆𝑟≥2 ∀𝑗∈,∀𝑗⊆𝑗.(12) Inequalities (12) state that subset of commodities 𝑗⊆𝑗of customer 𝑗∈can be covered either by a single route in 𝑗 𝑗or by at least two routes in  𝑗 𝑗. Note that these inequalities are meaningful only if |𝑗|≥3. Indeed, if |𝑗|= 2, they can be retrieved as an aggregation of Covering Constraints (4). Number partitioning polytope We propose a novel family of valid inequalities which exploits the multi-commodity aspect of the MC2DP. More precisely, given a customer 𝑗∈, these inequalities specify the possible combinations of routes to deliver the set of commodities 𝑗required by customer 𝑗. For each customer 𝑗∈, we denote by 𝑙 𝑗the subset of routes which deliver exactly 𝑙= 1,…,|𝑗|commodities to 𝑗, i.e., 𝑙 𝑗= {𝑟∈ ∶∑𝑘∈𝑗𝑎𝑟 𝑗𝑘 =𝑙}. Equalities |𝑗| ∑ 𝑙=1 𝑙∑ 𝑟∈𝑙 𝑗 𝜆𝑟=|𝑗|∀𝑗∈(13) ensure that the selected routes that serve customer 𝑗will exactly bring |𝑗|commodities to customer 𝑗. As an example, let  𝑗∈be a customer having a demand for three commodities, i.e., | 𝑗|= 3. Equality (13) for customer  𝑗states that the commodities of  𝑗can be covered by (i) a single route of 3  𝑗or (ii) one route of 2  𝑗and a route of 1  𝑗or (iii) three routes of 1  𝑗. Proposition 1. Equalities (13) are valid for the MC2DP. More precisely, inequalities ∑|𝑗| 𝑙=1 𝑙∑𝑟∈𝑙 𝑗 𝜆𝑟≥|𝑗|,∀𝑗∈, are implied by Covering Constraints (4) and inequalities |𝑗| ∑ 𝑙=1 𝑙∑ 𝑟∈𝑙 𝑗 𝜆𝑟≤|𝑗|∀𝑗∈(14) are valid for the MC2DP. Proof. It is straightforward that equalities (14) are valid for the MC2DP. Hence, we only need to show that ∑|𝑗| 𝑙=1 𝑙∑𝑟∈𝑙 𝑗 𝜆𝑟≥|𝑗|,∀𝑗∈ , are implied by Covering Constraints (4). Let 𝑗∈be a customer. By summing up the Covering Constraints (4) associated with 𝑗and swapping the summation order, we obtain ∑ 𝑟∈∑ 𝑘∈𝑗 𝑎𝑟 𝑗𝑘𝜆𝑟≥|𝑗|. Let 𝑟 𝑗denote the subset of commodities delivered to customer 𝑗by route 𝑟. We have ∑𝑘∈𝑗𝑎𝑟 𝑗𝑘 =|𝑟 𝑗|. The proof follows from partitioning the set of routes as =⋃|𝑗| 𝑙=0 𝑙 𝑗, where we denoted by 0 𝑗the subset of routes which do not visit 𝑗. Indeed, it holds |𝑗|≤|𝑗| ∑ 𝑙=0 ∑ 𝑟∈𝑙 𝑗|𝑟 𝑗|𝜆𝑟=|𝑗| ∑ 𝑙=1 𝑙∑ 𝑟∈𝑙 𝑗 𝜆𝑟.□ Remark that if we model the MC2DP by means of a set partitioning formulation, i.e., we impose the equality in Constraints (4), Equalities (13) become trivial. Indeed, they can be retrieved as an aggregation of the partitioning constraints. EURO Journal on Transportation and Logistics 13 (2024) 100139 4 M. Petris et al. Given a customer 𝑗∈and 𝑙= 1,…,|𝑗|, we introduce an auxiliary variable 𝑦𝑙 𝑗∈Z≥0defined as 𝑦𝑙 𝑗∶= ∑𝑟∈𝑙 𝑗 𝜆𝑟. Now, let 𝑗∶= {𝑦𝑗∈Z|𝑗| ≥0∶|𝑗| ∑ 𝑙=1 𝑙𝑦𝑙 𝑗≤|𝑗|} be the set of the integer points which satisfy inequality (14), rewritten in terms of 𝑦𝑙 𝑗variables. Proposition 2. The inequalities defining the convex hull of 𝑗,𝑗∈, are valid for the MC2DP. Determining the external description of a convex set is not an easy task, in particular in large dimensions. However, given that customers require at most three commodities in the benchmark instances of Gu et al. (2022) for the MC2DP, we explicitly derive the external description of the convex hull of sets 𝑗⊆Z3,𝑗∈. If the number of commodities is greater than three, software for polyhedral transformations such as PORTA (Christof and Löbel,2009) or PANDA (Lörwald and Reinelt,2015) can be used to determine the external description of the convex hull of sets 𝑗,𝑗∈. Note that inequalities (14) are meaningful only for customers 𝑗∈ who require at least three commodities, i.e., |𝑗|≥3. The external description of the convex hull of sets 𝑗,𝑗∈such that |𝑗|= 3 reads as follows: ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ 𝑦1 𝑗+ 2𝑦2 𝑗+ 3𝑦3 𝑗≤3(a) 𝑦1 𝑗−𝑦2 𝑗≥0(b) 𝑦2 𝑗≥0(c) 𝑦3 𝑗≥0.(d) (15) Inequalities (15c) and (15d) are trivial, indeed, they are implied by the definition of variables 𝑦𝑙 𝑗. Therefore, inequalities (15a) and (15b) are the only meaningful ones in the case of a customer 𝑗∈requiring three commodities (|𝑗|= 3); in terms of 𝜆variables, they are expressed respectively as |𝑗| ∑ 𝑙=1 𝑙∑ 𝑟∈𝑙 𝑗 𝜆𝑟≤|𝑗|∀𝑗∈∶|𝑗|= 3 (16) ∑ 𝑟∈1 𝑗 𝜆𝑟−∑ 𝑟∈2 𝑗 𝜆𝑟≥0 ∀𝑗∈∶|𝑗|= 3.(17) In conclusion, the number partitioning polytope valid inequalities we consider are (16) and (17). 5. Branch-Price-and-Cut algorithm We solve formulation [SC] by means of a Branch-Price-and-Cut (BPC) algorithm (Barnhart et al.,1998), i.e., a variant of the branchand-bound algorithm which deals with integer programming model with exponentially-many variables. Specifically, at each node of the branch-and-bound tree, the Master Problem (MP), that is the linear relaxation of formulation [SC], is solved by a column generation procedure (Desrosiers and Lübbecke,2005). If the solution of the MP is fractional, violated valid inequalities of Section 4.1 may be inserted and the column generation procedure is repeated while some valid inequalities are violated. Finally, branching rules are applied to ensure the integrality of the solution. We impose a time limit as a termination criterion for our BPC algorithm. In this section, we describe the main components of our BPC algorithm. Specifically, in Section 5.1 we present the column generation scheme applied in our BPC algorithm. In Section 5.2, we detail the management of the valid inequalities, and their impact on the pricing problem. Branching strategies and accelerating techniques are presented in Sections 5.3 and 5.4, respectively. 5.1. Column generation At each node of the branch-and-bound tree, a column generation procedure solves the MP defined on the exponentially-many variables 𝜆𝑟,𝑟∈, which correspond to the routes in the delivery echelon. The starting point is the Restricted Master Problem (RMP). The column generation procedure iteratively solves a Restricted Master Problem (RMP), i.e., the MP restricted to a subset of variables 𝜆𝑟. At each iteration of the procedure, after the RMP is solved, a subproblem, named pricing problem is solved. The aim of the pricing problem is to identify a variable (column) with the smallest reduced cost. If such a column has a negative reduced cost, it is added to the RMP in order to decrease (in a minimization problem) the current value of the solution, and the column generation procedure iterates. The procedure ends when the solution of the pricing problem is a non negative reduced cost column, proving the optimality of the MP. More precisely, the pricing problem is [PP] min{  𝐶𝑟∶𝑟∈} where  𝐶𝑟denotes the reduced cost of 𝜆𝑟variable. Note that set of routes can be partitioned per distribution centre, i.e., =⋃𝑜∈𝑜where 𝑜is the set of routes starting and ending at 𝑜. Hence, solving [PP] can be done by solving sequentially ||independent problems with the same structure: [PP(𝑜)] min{  𝐶𝑟∶𝑟∈𝑜}, 𝑜 ∈. Specifically, the aim of problem [PP(𝑜)] is to determine the most negative reduced cost 𝜆𝑟,𝑟∈𝑜, or to detect that none of them exists. The column generation procedure terminates once all problems [PP(𝑜)], 𝑜∈do not yield any negative reduced cost variable. In the following, we detail how a problem [PP(𝑜)] for 𝑜∈is formulated and solved. By denoting 𝜌𝑗𝑘 ≥0,∀𝑗∈, 𝑘 ∈𝑗and 𝜎𝑜𝑘 ≥0,∀𝑜∈, 𝑘 ∈as the optimal dual prices associated with Constraints (4) and (5), respectively, the reduced cost of a 𝜆𝑟,𝑟∈𝑜 variable is defined as follows:  𝐶𝑟=𝐶𝑟−∑ 𝑗∈∑ 𝑘∈𝑗 𝑎𝑟 𝑗𝑘(𝜌𝑗𝑘 −𝐷𝑗𝑘𝜎𝑜𝑘).(18) As mentioned in Section 3, the delivery echelon is a multi-depot version of the C-SDVRP. Hence, the problem [PP(𝑜)] is the pricing problem arising in Branch-Price-and-Cut approaches for the C-SDVRP (see Archetti et al.,2015;Gschwind et al.,2019) and is formulated as an Elementary Shortest Path Problem with Resource Constraints (ESPPRC) on a multi-graph (𝑜) = ((𝑜),(𝑜)). Such graph is analogous to the one presented in Gschwind et al. (2019) to formulate the ESPPRC in the context of the C-SDVRP. Vertex set (𝑜)contains two copies 𝑜′and 𝑜′′ of distribution centre 𝑜and two copies 𝑗′and 𝑗′′ of each customer 𝑗∈. Each arc of set (𝑜)is associated with two resources: demand  𝐷and cost  𝐶. Arc set (𝑜)contains: 1. an arc (𝑖′′, 𝑗′)for each arc (𝑖, 𝑗) ∈ to model the movement of a vehicle from vertex 𝑖to vertex 𝑗; the demand and cost are set to  𝐷𝑖′′𝑗′∶= 0 and  𝐶𝑖′′𝑗′∶= 𝐶𝑖𝑗 , respectively. 2. an arc (𝑗′, 𝑗′′)𝑗for each customer 𝑗∈and each subset 𝑗⊆𝑗to model the delivery of the commodities of 𝑗to 𝑗; the demand and cost are set to  𝐷𝑗 𝑗′𝑗′′ ∶= ∑𝑘∈𝑗𝐷𝑗𝑘 and  𝐶𝑗 𝑗′𝑗′′ ∶= − ∑𝑘∈𝑗(𝜌𝑗𝑘 −𝐷𝑗𝑘𝜎𝑜𝑘), respectively. Solving problem [PP(𝑜)] results in searching for negative reduced cost elementary paths in (𝑜)from 𝑜′′ to 𝑜′such that the resource consumption (demand) does not exceed the vehicle capacity 𝑄𝐷. To do so, we adopt a two phase procedure: Phase 1 computes the Pareto-optimal (demand, cost) pairs ( 𝐷𝑗 𝑗′𝑗′′ ,  𝐶𝑗 𝑗′𝑗′′ )for each customer 𝑗∈. EURO Journal on Transportation and Logistics 13 (2024) 100139 5 M. Petris et al. Phase 2 solves the ESPPRC on multi-graph (𝑜)which includes all arcs of type (𝑖′′, 𝑗′), and only the Pareto-optimal arcs of type (𝑗′, 𝑗′′)𝑗that have been computed in phase 1. Precisely, the ESPPRC is solved by means of a label setting dynamic programming technique (Feillet et al.,2004) which works with an implicit version of the bidirectional labelling search (see Righini and Salani,2006;Bode and Irnich,2012). The elementarity constraints are the bottleneck of such procedure, hence, we partially relax it by solving the ng-path relaxation (Baldacci et al.,2011) of the ESPPRC. For each customer 𝑗∈, we consider a fixed size ng-neighbourhood which includes the 10 closest customers to 𝑗and 𝑗itself. Remark that such relaxation allows a route to serve the same commodity to the same customer multiple times. Hence, the coefficients of the constraints and valid inequalities need to be updated accordingly: e.g., in the Covering Constraints (4),𝑎𝑟 𝑗𝑘 becomes an integer coefficient expressing the number of times customer 𝑗∈is delivered with commodity 𝑘∈𝑗by route 𝑟∈. The reader may refer to Archetti et al. (2015) and Gschwind et al. (2019) for further details. The resolution of the ESPPRCs is the bottleneck of our algorithm, hence, we heuristically solve the ESPPRC with the objective of quickly finding a negative reduced cost column. Precisely, we apply the same heuristic algorithms to solve the ESPPRC as those used in the column generation approach for the C-SDVRP proposed in Petris et al. (2023). One of these heuristics is a two-phase algorithm which exploits the multi-commodity aspect of the problem, while the others are based on reducing the pricing multi-graph by restricting the set of neighbours of the customers and by limiting the total number of splits in a route. When all the heuristics fail to identify a negative reduced cost column, we solve the ESPPRC exactly. 5.2. Management of the valid inequalities In this section, we first describe how the valid inequalities presented in Section 4.1 are considered in the pricing problem. Then, we present the cutting strategy adopted in our BPC algorithm. Impact of the valid inequalities on the pricing problem First, note that inequality (9) imposes a lower bound on the number of vehicles used in the collection echelon. Therefore, it has no impact on the pricing problem. The other inequalities presented in Section 4.1 are all robust, i.e. they do not change the structure of the pricing problem, and their associated dual prices have to be integrated into the objective function of pricing problems [PP(𝑜)], 𝑜∈, i.e. on the cost of arcs in multi-graph (𝑜). The arc costs in multi-graph (𝑜)are modified in the following way: Inequalities (10a) and (10b).Let 𝜏+≥0and 𝜏−≤0be the optimal dual prices associated with valid inequalities (10a) and (10b) respectively. The value 𝜏+∕2 + 𝜏−∕2 is subtracted from the cost of arcs of type (𝑖′′, 𝑗′), if vertices 𝑖or 𝑗represent distribution centre 𝑜. Inequalities (11).Let 𝜉′≥0be the optimal dual prices associated with the capacity cut (11) defined over the subset of customers ′⊆. Let 𝛿−(′)be the subset of arcs in graph entering in vertices of ′. The value 𝜉′is subtracted from the cost of arcs (𝑖′′, 𝑗′), for all (𝑖, 𝑗) ∈ 𝛿−(′). Inequalities (12).Let 𝛾𝑗𝑗≥0be the optimal dual prices associated with the inequality (12) identified by customer 𝑗∈and commodity subset 𝑗⊆𝑗. The value 2𝛾𝑗𝑗is subtracted from the cost of arcs (𝑗′, 𝑗′′)′ 𝑗, for all ′ 𝑗⊆𝑗that contain at least all the commodities of 𝑗, i.e., 𝑗⊆′ 𝑗. The value 𝛾𝑗𝑗is subtracted from the cost of arcs (𝑗′, 𝑗′′)′ 𝑗for all ′ 𝑗that contain some, but not all, commodities of 𝑗, i.e. ′ 𝑗∩𝑗≠∅and ′ 𝑗∩𝑗≠ 𝑗. Inequalities (16) and (17).Let 𝑗∈be a customer requiring exactly three commodities (|𝑗|= 3) and let 𝜓≥0and 𝜒≤0be the optimal dual prices associated with inequalities (16) and (17) defined on 𝑗. For all 𝑗⊆𝑗, the cost of arc (𝑗′, 𝑗′′)𝑗is modified as follows: value |𝑗|𝜓is subtracted, value 𝜒is added if |𝑗|= 2, and value 𝜒is subtracted if |𝑗|= 1. Management of the valid inequalities in the RMP Valid inequalities on vehicle bounds, namely (9),(10a) and (10b), are included in the formulation from the beginning of the solution procedure. Differently, a cut generation procedure manages the insertion of violated inequalities (11),(12),(16) and (17) in the RMP. Such a procedure is called at each node of the branch-and-bound tree of level at most equal to 5, if the associated solution of the RMP is fractional. Specifically, it separates the inequalities hierarchically according to the sequence : (11),(12),(16), and (17). When the separation of a given inequality fails, we separate the next one in the above order. The separation of inequalities (11) is done using the heuristic algorithms presented in Ralphs et al. (2003), namely the extended shrinking heuristic and the greedy shrinking heuristic. Then, although, inequalities (12) are exponentially-many, the size of the problem instances allows the separation by enumeration. The same separation strategy is applied for inequalities (16) and (17), whose number is linear in the number of customers ||. Finally, we limit the number of inequalities (11) to 100 in each cut generation round. For the other inequalities, we include all the violated inequalities. 5.3. Branching strategies Let (𝑥, 𝑞,  𝜆)be a fractional optimal solution of the MP at a certain node of the branch-and-bound tree. We consider seven branching rules that are hierarchically applied. In addition to these rules, the correctness of the algorithm requires the separation of a family of valid inequalities, namely the strong-degree inequalities. Rules 1 and 3 are specific for the MC2DP, while the other ones and the family of valid inequalities are used to solve the C-SDVRP by Branch-and-Price. The interested reader can refer to Gschwind et al. (2019) for more details about the branching strategy for the C-SDVRP. Rule 1 is on the number of vehicles traversing an arc in the collection echelon, i.e., on value 𝑥𝑖𝑜,𝑖∈,𝑜∈. Since 𝜆𝑟variables are not concerned by this rule, there is no impact on the pricing problem. Rule 2 is on the number of vehicles used at each distribution centre 𝑜∈in the delivery echelon, i.e., on value ∑𝑟∈𝑜  𝜆𝑟. Rule 3 forces the assignment of a delivery to a distribution centre. Specifically, given a distribution centre 𝑜∈, a customer 𝑗∈and a commodity 𝑘∈𝑗, we branch on value 𝑝𝑜 𝑗𝑘 ∶= ∑𝑟∈𝑜𝑎𝑟 𝑗𝑘  𝜆𝑟. The branching decisions related to this rule can be expressed as follows: commodity 𝑘required by customer 𝑗is either delivered from distribution centre 𝑜, i.e. ∑𝑜′∈⧵{𝑜}∑𝑟∈𝑜′ 𝑎𝑟 𝑗𝑘𝜆𝑟= 0; or not delivered from 𝑜, i.e. ∑𝑟∈𝑜𝑎𝑟 𝑗𝑘𝜆𝑟= 0. Note that both decisions entail modifications in the pricing problem. As an example, if the first decision is imposed, then we prevent the pricing problem from generating routes starting and ending at distribution centres 𝑜′∈⧵{𝑜}and delivering commodity 𝑘to customer 𝑗. Arcs of type (𝑗′, 𝑗′′)𝑗, with 𝑗⊆𝑗and 𝑘∈𝑗, are removed from all multi-graphs (𝑜′),𝑜′∈⧵{𝑜}. Rule 4 is on the number of visits to each customer 𝑗∈from a distribution centre 𝑜∈. Rule 5 considers the flow on the edges in the delivery echelon coming from a specific distribution centre. EURO Journal on Transportation and Logistics 13 (2024) 100139 6 M. Petris et al. Rules 6 and 7 implement the Ryan and Foster branching rules (Ryan and Foster,1981) which force the two customer requests to be served by different routes or by the same route. Such rules imply the addition of non-robust constraints in the RMP . The management of the associated dual variables in the labelling algorithm used to solve the pricing problem can be found in Gschwind et al. (2019). These seven rules are sufficient to ensure the correctness of the algorithm when only elementary routes are present in formulation [SC]. Indeed, Rule 1 guarantees the integrality of the variables of the collection echelon. Then, regarding the delivery echelon, Rule 3 assigns the deliveries to a specific distribution centre. Once these assignments are done, Rules 2 and 4–7 are enough to guarantee the correctness of the algorithm. Indeed, the delivery echelon is a multi-depot version of the C-SDVRP and such rules ensure the integrality of a solution for the C-SDVRP (see Gschwind et al.,2019). However, as mentioned in Section 5.1, we relax the elementarity requirement of the routes in the second echelon via the ng-path relaxation when solving the pricing problem. Consequently, formulation [SC] may contain routes that serve the same commodity to the same customer more than once. In such a case, applying only Rules 1–7 might lead to a fractional solution as shown by Gschwind et al. (2019) for the C-SDVRP. Hence, after applying Rule 7, strong-degree inequalities (Contardo et al.,2014) have to be separated to ensure providing an integer solution. The strong-degree inequalities read as: ∑ 𝑟∈ 𝜉𝑟 𝑗𝑘𝜆𝑟≥1 ∀𝑗∈,∀𝑗∈𝑗, where 𝜉𝑟 𝑗𝑘 is a binary coefficient with value one if route 𝑟∈ delivers customer 𝑗∈with commodity 𝑘∈𝑗. In our branching strategy, when none of the seven rules are applicable, we separate these inequalities. As for Rules 6 and 7, these inequalities are non-robust. The management of the associated dual variables in the labelling algorithm invoked to solve the pricing problem is described in Contardo et al. (2014). The branch-and-bound tree is explored according to a best-bound first strategy to favour the improvement of the dual bound. The strategies to select the branching decisions are presented in the following. For rule 2, we branch on the fractional value closest to 0.5. For rules 6 and 7, we branch on the first fractional value that is found. For all the other rules, we consider a two-round strong branching procedure (S., 2012) similar to the one presented in Pessoa et al. (2020). In the first round, we evaluate at most 100 branching candidates according to the product rule (Achterberg,2007). More precisely, each candidate gives rise to two branching decisions 𝑑1and 𝑑2and is evaluated by applying such decisions to the RMP and by solving it without generating columns. Then, each candidate is assigned with a score 𝑠𝑐(𝑑1, 𝑑2) = max{𝜖, 𝛥𝐿𝐵1} × max{𝜖, 𝛥𝐿𝐵2}, where 𝜖= 10−6 and 𝛥𝐿𝐵𝑖 is the increase of the lower bound obtained by applying decision 𝑑𝑖to the RMP. The three candidates with the highest scores are sent to the second round, where the same evaluation criterion is used to select the winning candidate. Differently from the first round, here, 𝐿𝐵1and 𝐿𝐵2 are the values of the RMP after a single column generation iteration where the pricing problem is solved heuristically. The strong branching procedure is employed in nodes of the branchand-bound tree of level at most 5. In the other levels, we evaluate the branching candidates based on the fractional value closest to 0.5 for all the rules. 5.4. Accelerating strategies The BPC algorithm incorporates the following accelerating strategies: Initialization of set .We initialize the set of routes to avoid very large dual prices at the first iterations of the column generation procedure which may slow down the pricing solution (Desaulniers,2010). Specifically, for each distribution centre 𝑜∈, we include round-trips (0-𝑗-0) to each customer 𝑗∈, which deliver the commodities of each subset 𝑗⊆𝑗requested by 𝑗. In addition, we modified the randomized Clarke-Wright algorithm (CW) (Clarke and Wright,1964) proposed in Battarra et al. (2008) to take into account the multicommodity aspect of our problem. The algorithm is run 10 times per distribution centre and the obtained routes are inserted into . Heuristic column generators. Before solving the pricing problem to optimality, we consider heuristic column generators to speed up the solution of problems [PP(𝑜)], 𝑜∈. As mentioned in Section 5.1, each problem [PP(𝑜)] is the pricing problem arising in a BPC algorithm for the C-SDVRP. Hence, we apply the same heuristic scheme used in Petris et al. (2023) which proved to be effective in accelerating such pricing problems. This scheme considers two reduced graph heuristics and the two-phase heuristic introduced in Petris et al. (2023) which proved to be effective in dealing with the multi-commodity aspect of the C-SDVRP. The two reduced graph heuristics reduce the size of multi-graphs (𝑜),𝑜∈by limiting both the possibilities of travelling between customers and of deliveries to customers. In the two-phase heuristic, the aim of the first phase is to compute a set of promising customer sequences by solving the ESPPRC on a modified version of multi-graphs (𝑜) where only one delivery per customer is allowed. Specifically, when visiting a customer, the least consuming commodity is delivered and all the profitable dual prices are collected. In the second phase, for each of the customer sequences generated by the first phase, we solve the ESPPRC on the associated acyclic graphs to obtain all negative reduced cost routes arising from the sequence. We refer to Petris et al. (2023) for more details. Restricted master heuristic. We invoke a restricted master heuristic, which consists in solving the formulation [SC] restricted to the subset of variables generated so far, to obtain good upper bounds. Such a technique helps to reduce the integrality gap (see Archetti et al.,2013). Note that variables 𝜆𝑟are then binary. We call the restricted master heuristic every 1000 explored nodes in the branch-and-bound tree as well as when the time limit of the algorithm is reached. In this latter case, we apply a local search procedure based on an adapted version of the mathematical programming operator proposed for the C-SDVRP in Gu et al. (2019). Specifically, we generalized such an operator to deal with the two-echelon case. When the restricted master heuristic is called during the tree exploration a time limit of 3seconds is imposed, while the time limit is 30 seconds when the algorithm terminates. 6. Computational experiments We implemented the BPC algorithm in C++ and compiled it in release mode under a 64-bit version of MS Visual Studio 2019. IBM CPLEX 12.9.0 (64-bit version) is used as a solver. We performed the experiments on a 64-bit Windows machine equipped with a 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. For each run of the algorithm, we impose one hour time limit and allow a single thread. In this section, first, we describe the characteristics of the benchmark instances for the MC2DP introduced in Gu et al. (2022). Then, we discuss the impact of valid inequalities (12),(16) and (17). Finally, we evaluate the effectiveness of the BPC algorithm against solving the compact formulation for the MC2DP presented in Gu et al. (2022) with a commercial solver and we present the results obtained by the BPC algorithm on the benchmark instances. EURO Journal on Transportation and Logistics 13 (2024) 100139 7 M. Petris et al. Table 1 Characteristics of the sets of instances. Set # Characteristics || || ||Description S64 8 30 2, 3 Base set S 164 8 30 2, 3 Unbalanced supplier locations (6-2) S 264 8 30 2, 3 Unbalanced supplier locations (8-0) S 164 8 30 2, 3 Unbalanced customer locations (5-10, with 𝛿= −5,30) S 264 8 30 2, 3 Unbalanced customer locations (5-10, with 𝛿= 10,30) S 364 8 30 2, 3 Unbalanced customer locations (10-5, with 𝛿= −5,30) S 464 8 30 2, 3 Unbalanced customer locations (10-5, with 𝛿= 10,30) S𝑂32 8 30 2 Unbalanced available amounts at the suppliers S𝑎𝑑𝑑 164 10 30 2, 3 Increased number of suppliers to 10 S𝑎𝑑𝑑 264 12 30 2, 3 Increased number of suppliers to 12 S𝑎𝑑𝑑 164 8 50 2, 3 Increased number of customers to 50 S𝑎𝑑𝑑 264 8 70 2, 3 Increased number of customers to 70 Small 36 4, 6 10, 15, 20, 25 2, 3 Small instances 6.1. Benchmark instances Gu et al. (2022) introduced artificial instances as well as instances arising from a real-world case study in the context of a short and local fresh food supply chain. In the following computational experiments, we only consider the artificial instances. Indeed, the sizes of the instances based on the case study are too large to be tackled efficiently with the BPC algorithm. First, Gu et al. (2022) generated a base set of 64 artificial instances Swith two distribution centres (||= 2), eight suppliers (||= 8) and 30 customers (||= 30). The features of the delivery echelon are based on the 64 small instances proposed in Archetti et al. (2016) for the CSDVRP. Each C-SDVRP instance gives rise to a MC2DP instance where the locations of one distribution centre and 15 customers are the ones of the C-SDVRP instance. Such distribution centre and 15 customers are duplicated and their locations are modified by applying a translation of parameter 𝛿= (30,30) to their coordinates. Customer demands are also as in the C-SDVRP instance. Four suppliers are randomly located around each distribution centre. The availability of each commodity at the suppliers is calculated as a fraction of the total demand. The commodity availabilities are the same for all the suppliers. Then, Gu et al. (2022) produced 12 additional sets of instances by applying modifications to one of the characteristics of the base set, such as the suppliers/customers locations, the number of suppliers/customers or the available quantities at the suppliers. In all sets of instances, the number of distribution centres is fixed at two. In Table 1, we summarize the main characteristics of all sets of instances. Each row of the table represents a set of instances. The columns of the table report: set: the name of the set of instances; #: the number of instances in the set; ||: the number of suppliers; ||: the number of customers; ||: the number of commodities; description: a brief description of the main characteristic of the set. In such an entry, we write 𝑛1−𝑛2 to express the distribution of the suppliers/customers around each distribution centre, meaning that 𝑛1suppliers/customers are located around one distribution centre and 𝑛2are located around the other one. Parameter 𝛿is a translation parameter used to determine the locations of the customers/suppliers around the two distribution centres. We refer to Gu et al. (2022) for further details regarding the generation of the set of instances. 6.2. Impact of valid inequalties In this section, we assess the impact of valid inequalities. To do so, we consider the 32 instances of base set Shaving three commodities. Indeed, as mentioned in Section 4.1, if the number of commodities is equal to two, inequalities (12),(16) and (17) can be retrieved as an aggregation of Covering Constraints (4). We examine the following four variants of the BPC algorithm. BPC: valid inequalities on bounds on the number of vehicles are inserted, and no valid inequalities is separated in the course of the algorithm; BPC+CC: only capacity cuts (valid inequalities (11)) are separated; BPC+SC+NP: only the inequalities arising from the set covering polytope (SC), i.e., inequalities (12), and the ones arising from the number partitioning polytope (NP) are separated, i.e., inequalities (16) and (17), are separated; BPC+CC+SC+NP: all valid inequalities are separated. Each row of Table 2 corresponds to a BPC variant. The first two columns report the average lower bound (avg.LB) and time (avg.t[s]) at the root node of the branch-and-bound-tree. The next four columns show the results at the end of the execution of the corresponding BPC variant: the average number of nodes of the branch-and-bound tree (avg.#nodes), the average lower bound at termination (avg.LB) the average time (avg.t[s]) and the number of instances solved to optimality (#opt./#inst.) over the 32 instances. As expected, BPC yields the worst results solving only six instances out of the 32 considered. Variant BPC+SC+NP solves an additional instance w.r.t. BPC, however, the improvement of the lower bound at the root node is mediocre. The best results are obtained when the well-established capacity cuts are separated, namely with variants BPC+CC and BPC+CC+SC+NP. Both variants solve the same 14 instances to optimality and yield the best lower bounds at the root node, being on average equal to 1000.35 and 1001.31 in BPC+CC and BPC+CC+SC+NP, respectively. The same remark applies to the lower bounds at termination which is on average equal to 1039.61 in BPC+CC and BPC+CC+SC+NP. In both cases, lower bounds at the root node and at termination improve significantly with respect to BPC. We also observe that the addition of inequalities (12),(16) and (17) in BPC+CC+SC+NP slightly improves the results with respect to BPC+CC in terms of lower bounds at the root node, number of explored branch-and-bound nodes and solution time. Hence, we choose BPC+CC+SC+NP as the configuration for our BPC algorithm. 6.3. Evaluation of the BPC algorithm The aim of this section is to evaluate the effectiveness of the BPC algorithm. To do so, we compare the results obtained by the BPC algorithm on the instances of set small with the ones obtained by solving a compact formulation for the MC2DP on the same instances with CPLEX 12.8. The latter results are retrieved from Gu et al. (2022) and were obtained on a machine with Intel (R) Core(TM) i7-4600U EURO Journal on Transportation and Logistics 13 (2024) 100139 8