A Multi-Commodity Two-Echelon Capacitated Vehicle Routing Problem with Time Windows: Model Formulations and Solution Approach

This paper studies the multi-commodity two-echelon capacitated vehicle routing problem with time windows. This problem takes customer-specific origin to destination, nonsubstitutable demands into account. The main components are: (i) first echelon tours, (ii) second echelon tours and (iii) movement of commodities from first echelon to second echelon, considering route connection and synchronization. Exploiting the structure of the problem, we propose a decomposition scheme which decouples the first and second echelon routing problems, by deleting the third component, and re-couples them by using constraints ensuring movement of commodities between the two echelons. Based on this decomposition, we present (i) a number of model formulations for the MC-2E-VRPTW, and (ii) develop an exact solution approach for this real-life problem. A computational study on a comprehensive set of instances shows the effectiveness in terms of computational effort and solution quality.


Introduction
An important city logistics solution approach for planning freight vehicles going into cities, each delivering small drop quantities, is to consolidate these volumes outside the city, in so-called urban distribution centers. In this case, a second layer of distribution centers is introduced, with the sole purpose to increase efficiency both stream upward (from the shipping warehouses) and stream downward (towards the customers) (Savelsbergh and Van Woensel, 2016). As suggested by (Cuda et al., 2015), these two-echelon problems receive more and more attention (both in the literature and in practice) but mostly these efforts are focused on basic problem variants. Many important real-life issues still need to be incorporated and analyzed.
More specifically, we study an important real-life variant of the two-echelon vehicle routing problem (2E-VRP) where customer-specific (non-substitutable) goods are transferred to customers by using intermediate facilities. There are three types of facilities that are located in a hierarchical structure: (i) Freight depots (distribution centers) are intermodal logistic sites located at the first level and act as the sources of consolidated goods; (ii) Intermediate facilities are used to consolidate goods coming from depots into same vehicle tours and are located at the second level; (iii) Final customers are located at the last level and are destinations of the goods. Whereas the (freight) depots are quite far from the end customers, the intermediate facilities are located closer to the customers, facilitating more efficient last mile distribution with greener vehicles (Grangier et al. (2016)).
The distinguishing real-life feature added to the two-echelon problem considered, is the inclusion of customer-specific demand (i.e. non-substitutable origin-destination flows, e.g., parcel, package). In this case, each demand starts from a specific depot and is supposed to be delivered to a specific customer. Most of the existing studies in the field of 2E-VRP focused on the basic variant of the problem without customer-specific demand. Clearly, this assumption is not realistic in practice (where parcels are typically non-interchangeable among customers). Although some papers introduced the concept of customer-specific de-mand in the 2E-VRP (e.g., Crainic et al. (2009)), the related routing problems are not treated extensively from an optimization point of view.
Our problem on-hand is referred to as the multi-commodity two-echelon capacitated vehicle routing problem with time windows (denoted as MC-2E-VRPTW in the remainder of this paper). We handle customer-specific demands by defining commodities. A commodity consists of the destination customer, its origin depot, a specific volume, and a time window in which the delivery should take place. Note that every customer directly corresponds to a single commodity.
In what follows, we use the terminology proposed by Crainic et al. (2010): the intermediate facilities are called satellites, the first echelon vehicles are referred to as urban vehicles (e.g., large trucks) and the second echelon vehicles are called city freighters (e.g., vans, bikes, etc.). The MC-2E-VRPTW is an extension of the 2E-VRP, which differs significantly by additionally considering multiple commodities and time windows for customers. It is known from the literature that the 2E-VRP is a complex problem for which only a limited num-ber of algorithms (both exact and heuristic) are proposed. For an overview of two-echelon distribution systems, we refer the interested reader to the survey by Cuda et al. (2015).
The objective of this paper is to present (i) a number of model formulations for the MC-2E-VRPTW, and (ii ) develop an effective exact solution approach. The contributions of this paper are as follows: 1. The MC-2E-VRPTW model formulations involve arc-based and path-based formulations, as well as a combined arc-path based model. Solution approaches involving branch-and-price based algorithms are discussed for the developed model formulations.
In this way, we assess the quality of all model formulations and solution approaches, from an optimal solution perspective. Note that the current literature has only focused on the basic variant of the 2E-VRP, ignoring customer-specific demands.
2. We build an exact solution algorithm to deal with our problem. We extend the solution approach developed by Dellaert et al. (2019). This requires better lower bound estimates for the first echelon tours and new types of capacity constraints in the 2-path formulation.
3. The numerical evaluation through a large set of instances demonstrates the power of the developed models and their respective solution methods. We study the strengths and weaknesses of the proposed solution approaches for the different formulations. We are able to solve instances up to 100 customers and 5 satellites to (near-)optimality.
The remainder of this paper is organized as follows. A literature review of the 2E-VRP is given in Section 2. The problem description is discussed in Section 3. We describe the modeling approaches and the different mathematical formulations in Section 4. Details of the MC2E-2P formulation are given in Section 5. The computational experiments are presented and discussed in Section 6. We conclude and give suggestions for future research in Section 7.

Literature review
The problem addressed in this paper fits within the larger two/multi-echelon routing and location-routing literature, but includes attributes that few papers in that literature do. Cuda et al. (2015) discuses both aspects in their 2E-VRP survey, while the location-routing surveys of Nagy and Salhi (2007) and, more recently, Prodhon and Prins (2014) discuss the 2E-VRP literature. Crainic et al. (2009) introduced a 2E-VRP with origin-destination time-dependent demand, capacitated intermediate facilities, and time-dependent, synchronized, multi-depot, heterogeneous fleets (on all echelons), multi-tours, and time windows. The authors introduce a path-based MIP formulation and discuss solution avenues, but do not present any computational experiments. Perboli and Tadei (2010); Perboli et al. (2011) formally introduced the term 2E-VRP for the basic problem setting where most of the attributes mentioned above are not included (e.g., the problem is static with no synchronization, the demand is single substituable product, fleets are homogeneous). The authors proposed a MIP formulation and derived valid inequalities. Huang et al. (2018) study the delivery cost impact of route choice restrictions in the 2E-VRP setting. A two-echelon multi-trip vehicle routing problem with a dynamic satellite is considered in He and J.Li (2019).
A limited number of papers concentrate on developing exact algorithms for the 2E-VRP. Contardo et al. (2012) propose lower and upper bounds for the 2E-CVRP combining a branch-and-cut algorithm based on addressing a new 2-index formulation of the problem and the ALNS of Hemmelmayr et al. (2012). An exact algorithm is also proposed by Baldacci et al. (2013), where a new mathematical formulation for the problem is introduced and used to derive valid lower bounds. The authors develop an exact algorithm based on a decomposition of the problem into a limited set of multi-depot capacitated vehicle routing problems with side constraints. Gonzalez-Feliu (2008) presents a column generation method to solve the linear programming relaxation of the problem. Jepsen et al. (2013) show that the proposed mixed integer liner programming formulation by Perboli et al. (2011) may not provide a correct upper bound when more than two satellites are selected in the solution. They propose a branch and cut method, based on an edge flow model. Santos et al. (2013) and Santos et al. (2014) present two branch-and-price implementations. Marques et al. (2019) propose a branch-cut-and-price algorithm, including several advances features to handle the 2E-CVRP. Dellaert et al. (2019) worked on branch-and-price based algorithms for the 2E-CVRP considering time windows.
As this brief literature survey indicates, several important attributes are still not addressed in the 2E-VRP literature, in particular the multi-commodity, non-substitutable origin-to-destination demand, in combination with customer time windows.
The rest of the paper presents the formal description, formulations, solution methods, and experimental results we propose for the MC-2E-VRPTW.

Problem description
Consider a directed graph G = (V, A) with the vertex set V and the arc set A. The vertex set V is partitioned into V = P ∪ S ∪ Z where P is the set of depot locations, S is the set of satellite locations and Z is the set of customer locations. The edge set A is partitioned to A = A 1 ∪ A 2 , where A 1 is the set of first echelon arcs from node i ∈ P ∪ S to node j ∈ S and from node i ∈ S to node j ∈ P (returning to depot) and A 2 is the set of second echelon arcs from node i ∈ S ∪ Z to node j ∈ Zand from node i ∈ Z to node j ∈ S (returning to satellite). There is a travel cost c ij associated with each arc (i, j) ∈ A which represents the transportation cost of a vehicle traveling on the arc form node i to node j. Moreover, there is a travel time t ij associated with each arc (i, j) which represents the travel time of a vehicle on the arc.
Each satellite is open and ready to give service. A constant service time s i is considered for each satellite i where the de-consolidation and consolidation occurs. Clearly, in many cases, a variable service costs could be incurred as well, representing the loading and offloading of the vehicles. To keep the models tractable, we leave the inclusion of this variable costs for future research. In practice (van Duin et al., 2010;Kin et al., 2018), most of the satellites (or Urban Consolidation Centers) have ample capacity, as they mainly operate as a cross-dock. Mostly, no inventory is carried in the satellites. More restrictive, could be e.g. the docking capacity for the vehicles (city freighters and urban vehicles). This dimension is left out of scope in this paper and we leave this for future research.
Each customer z ∈ Z has a demand of size d z from the depot p z which should be delivered within hard time window [a z , b z ] (the delivery of the commodity is not allowed either before the start of the window or after the end of this window). There is a service time s z , required to unload the commodity from a city freighter and deliver it to the customer z.
A homogeneous fleet of urban vehicles T is used to transfer commodities from depots to satellites. Each urban vehicle τ ∈ T has a capacity of K (1) and a fixed cost h (1) paid per tour. Similarly, a homogeneous fleet of city freighters Υ is used to transfer freight from satellites to customers. Each city freighter υ ∈ Υ has a capacity of K (2) and associated fixed cost h (2) paid per tour.
Define an itinerary of a commodity as an urban vehicle tour (starting from the origin depot of the commodity) continued by a transshipment operation at a satellite and a city freighter tour used to transfer the commodity from the satellite to its destination customer. The associated urban vehicle tour starts (and ends) at the depot where the commodity is available and visits the associated satellite where the transshipment occurs. Similarly, the associated city freighter tours start (and end) at the satellite and visit the customer. An itinerary is feasible if the associated city freighter departs from the satellite after the associated urban vehicle arrived at the satellite and the time windows of the customer is respected.
Every shipment follows the two echelons connected through the satellites where deconsolidation and consolidation activities take place: • First echelon Consolidated commodities are picked up by urban vehicles (first echelon vehicles) from depots and are delivered to the satellites. Each urban vehicle tour starts from a depot, delivers a subset of commodities to a subset of satellites and ends at the same depot.
• De-consolidation and consolidation Commodities received from one or more urban vehicle tours are de-consolidated first and then consolidated to one or multiple city freighters at satellites. Satellites have no capacity constraints.
• Second echelon Consolidated commodities are distributed to the final customers by city freighter tours (second echelon vehicles). Each city freighter tour starts from a satellite and delivers a subset of commodities to their customers and ends at the same satellite.
Transferring goods from one urban vehicles to city freighters at satellites requires a connection and synchronization between the urban vehicle and the city freighter which are used to transfer the same commodity. (i) Connection: the urban vehicle should deliver the commodity to the starting satellite of the city freighter. (ii) Synchronization: the departure time of the city freighter should be later than the arrival time of the urban vehicle to the satellite (clearly, a city freighter which is transferring a subset of commodities should wait for the arrival of all urban vehicles which bring at least one of these commodities).
A solution for the problem consists of the first and second echelon routing decisions such that these two echelon routes are connected and synchronized. Additionally, for each commodity, its path (as a combination of the first and second echelon routes) through the network is described in full (from origin to destination via a specific satellite). Therefore, an optimal solution of the problem consists of the vehicle tours of both echelons such that each commodity is delivered to its associated customer through a satellite, while the time windows of customers are respected, and the total transportation costs are minimized. Figure 1 shows an example of the problem with 2 depots {A, B}, 3 satellites {i, ii, iii} and 8 customers {1, ..., 8}. A possible feasible solution is shown, which consists of two urban vehicle tours and three city freighter tours. An example of the consolidation activities can be observed at satellite (ii) where the commodities of customers {2, 3, 5, 7} delivered by two urban vehicle tours, originating from depots A and B, are consolidated to one city freighter tours to finally deliver them to the customers.

Mathematical model formulations
In this section, we first discuss some of the modeling choices for the MC-2E-VRPTW. Then, the mathematical formulations are explained.

Modeling discussion
A straightforward arc-based formulation is inspired from the arc-based formulation of the (single commodity) 2E-VRPTW proposed by Dellaert et al. (2019). Dealing with multiple commodities brings the need to introduce extra variables and constraints to show how these commodities are transferred from depots to customers. Specifically, we introduce assignment variables to show which urban vehicle, city freighter and satellite is used to handle each commodity. These variables are used to introduce constraints which assure that the first and second echelon vehicle tours are connected and synchronized. This model has the typical challenges of the arc-based models. It is intractable to all but small size instances of the problem. This poor computational performance mainly arises from the fact that the LPrelaxation of the mathematical formulation of this model provides poor lower bounds for the problem.
Inspired on the framework introduced by Crainic et al. (2009), using an integrated modeling approach, we propose a 3-path model formulation, denoted as the multi-commodity twoechelon three-path formulation, or MC2E-3P. The three types of paths are in this model: (1) the first echelon tours, (2) the second echelon tours and (3) the interconnected two-echelon structures (itineraries) which manage the flow and synchronization for the demands. However, solving this proposed 3-path formulation is not straightforward using column generation (see Section 4.3 for a detailed discussion). The challenge is that the reduced cost formulations rely on unknown dual variables (see Section 4.3). There is some literature to estimate the unknown dual variables (e.g., Rezaei (2016) and Gendreau et al. (2002)). Also note that Muter et al. (2013) proposed a simultaneous column-and-row generation mechanism to deal with a similar issue.
As an alternative, we decouple the two routing elements by eliminating the itinerary paths, but re-coupling them with constraints which ensure the synchronization. This decomposition makes the two routing elements independent. Therefore different modeling approaches can be used for each routing element independent of the other one. This approach leads to two model formulation variations: a 2-path formulation and an arc-path formulation. The 2-path formulation builds on the approach of Baldacci et al. (2013) and Dellaert et al. (2019). The combined arc-path model benefits from advantages of both arcbased and path-based modeling approaches. In both cases, because of the high density of the second echelon graph, a path-based modeling approach is used for the second echelon. Additionally, extra variables and constraints are used to connect and synchronize the first and second echelon vehicle tours.
In the next sections, we describe the different model formulations in detail.

A 2-arc formulation (MC2E-2A)
In the 2-arc formulation, we model all flow decisions as arc variables. Specifically, the decision variables are categorized into three main categories: • First echelon routing (arc-based): -Route definition: r τ ij = 1 if arc (i, j), i ∈ P ∪ S, j ∈ S, is traveled by urban vehicle τ ∈ T ; 0 otherwise; -Vehicle selection: u τ = 1 if vehicle τ ∈ T is used at the first echelon; 0 otherwise.; -Commodity-urban vehicle assignment: y τ z = 1 if demand of customer z is assigned to vehicle τ ∈ T ; -Timing: ω τ s shows the arrival time of the urban vehicle τ at satellite s.
• Second echelon routing (path-based): -Route definition: x υ ij = 1 if arc (i, j), i ∈ S ∪ Z, j ∈ Z, is traveled by city freighter υ ∈ Υ; 0 otherwise; -Vehicle selection: u υ = 1 if city freighter υ ∈ Υ is used at the second echelon; 0 otherwise;(note that each vehicle can make only one trip) -Commodity-city freighter assignment: y υ z = 1 if demand of customer z is assigned to city freighter υ ∈ Υ; -Timing: ω υ z shows the arrival time of the city freighter υ at customer z.
• Connecting the first and second echelon routes: -Commodity-satellite assignment: q s z = 1 if the demand of the customer z is satisfied through satellite s.
This formulation is referred to as the MC2E-2A formulation because the arc flow variables are used to represent the components of the problem at both echelons. The MC2E-2A formulation for the MC-2E-VRPTW can be written as follows: x υ zs = u υ , ∀υ ∈ Υ (10) Objective function (1) minimizes the total cost which consists of first echelon transportation cost, second echelon transportation cost and vehicle usage costs.
Constraint (2) is the flow conservation constraint for the first echelon. Constraint (3) relates the urban vehicle selection and first echelon routing variables. Constraint (4) ensures that if a commodity is assigned to an urban vehicle, then that vehicle departs from the depot where the commodity is available. Constraint (5) is the capacity constraint for urban vehicles. Constraint (6) ensures that each commodity is assigned to one urban vehicle. Constraint (7) relates the arrival time of an urban vehicle to a satellite with the arrival time at the preceding location (either a depot or a satellite). Constraint (8) relates the commodity-satellite assignment variables with the first echelon routing variables. Actually, it assures that if a commodity is assigned to a satellite (q s z = 1) and an urban vehicle (y τ z = 1), then the urban vehicle should bring the commodity to the satellite ( k∈P ∪S r τ sk > 0). Constraint (9) is the flow conservation constraint for the second echelon. Constraint (10) relates the city freighter selection and second echelon routing variables. Constraint (11) assures that if a commodity is assigned to a city freighter and is handled through a satellite, then the city freighter should start from the satellite and visit the destination customer of the commodity. Constraint (12) ensures that if a commodity is assigned to a city freighter, then that vehicle visits the destination customer of the commodity. Constraint (13) is the capacity constraint for city freighters. Constraint (14) assures that each commodity is assigned to a city freighter.
Constraint (15) relates the arrival time of an urban vehicle (to a satellite) and arrival time of a city freighter (to a customer) if they meet at a satellite to transfer a commodity and the customer is visited as the first customer in the city freighter tour. Constraint (16) relates the arrival times of a city freighter to a pair of customers (i and j) if j is visited immediately after i. Constraints (17) and (18)

A 3-Path formulation (MC2E-3P)
Let R denote the set of all urban vehicle tours. Each urban vehicle tour r ∈ R starts (and ends) at a depot and visits a subset of satellites S r with an associated cost c r . This captures the first echelon transportation cost and the vehicle usage cost. In a similar way, let L denote the set of all feasible city freighter tours. Each city freighter tour l ∈ L starts (and ends) at a satellite and visits a subset of customers Z l . The associated second echelon transportation cost and the vehicle usage cost are captured in the city freighter tour cost c l .
Moreover, let B z be the set of all feasible itineraries for the commodity of the customer z. Each itinerary β ∈ B z starts from the depot p z , employing urban vehicle tour r(β), arrives at satellite s β ∈ S r(β) , takes a city freighter tour l(β) together with demands of other customers Z l(β) , departs from the satellite and arrives at final customer within the time window. As arrival and departure times of the urban vehicle tours and city freighter tours are known, it is straightforward to establish the feasibility of an itinerary.
Let us define the following decision variables: • First echelon routing: λ r = 1 if urban vehicle tour r is selected; 0 otherwise; • Second echelon routing: δ l = 1 if city freighter tour l is selected; 0 otherwise; • Commodity itineraries: γ β = 1 if itinerary β is selected; 0 otherwise.
The 3-Path formulation of the MC-2E-VRPTW (denoted as MC2E-3P) becomes: The objective function (26) minimizes the total cost which is the sum of urban vehicle and city freighter services cost. Constraint (27) ensures that for any urban vehicle tour the capacity of the associated urban vehicle is respected. Note that, for each urban vehicle tour r and the commodity of customer z, B r z = {β ∈ B z |r(β) = r, r ∈ R} shows all itineraries of the commodity which use that urban vehicle tour. In a similar way, for each city freighter tour l and the commodity of customer z, B l z = {β ∈ B z |l(β) = l, l ∈ L} shows all itineraries of the commodity which use that city freighter tour. Constraint (28) relates city freighter tours with itineraries. Constraint (29) indicates that each commodity should be satisfied employing only one itinerary. Constraints (30 -32) are the domain constraints.

A 2-Path formulation (MC2E-2P)
Observe again that the itineraries govern the flow of the demands, and impose synchronization requirements. Once the itineraries are eliminated, the two echelon routes become decoupled and thus constraints should re-couple them. Note that, for each commodity an urban vehicle and a city freighter should meet at a satellite such that the departure time of the city freighter is after the arrival time of the urban vehicle and a possible service time.
For the 2-Path formulation (denoted as MC2E-2P), we fix a first level routing, by which the variables λ r are no longer decision variables (Baldacci et al., 2013). Following Dellaert et al. (2019), we first determine good lower bounds for all relevant first echelon routes and consider them as candidates according to their lower bounds. Knowing the arrival times of all first level trips and the latest departure times of all second level trips, we can determine for every column/second level trip whether it is feasible for this first echelon routing and which quantities must be delivered by which combination of first level trips. At this phase, it is not necessary to set the commodity itineraries and assign customers of second level trips to first level trips. It is sufficient to make sure that the aggregate capacity to do so is present. The actual assignment problem has no influence on the total cost and only has to be solved in the end when we have an integer feasible solution for the first and second level routing and the only remaining decision is the itinerary choice γ β .
For the aggregate capacity problem, we consider the demand from the depots independently. Let M be the set of first level trips that we consider and L M the set of feasible second level trips for the set M . Let C M dl denote the subset of M that contains the feasible first level trips leaving from depot d for the second level trip l. It is an empty set when no customers from depot d are involved in trip l. For every depot, the number of possible feasible sets is at most two to the power of the number of first level trips from that depot minus 1. Number all combinations 1, ..., S M , where each set C M dl corresponds with such a combination. We denote the demand for trip l to be delivered by combination s as d sl . The total capacity of each combination is based on the total capacity of the first level trips in the combination, denoted by K M s . Finally, let b lz indicate whether customer z is served in trip l. Now, we give our MC2E-2P formulation: In order to illustrate Constraint (34), consider the following example. Suppose that five first level trips visit (i.e. arrive at) satellite 1. The trips 1,2 and 3 are departing from depot 1 and have arrival times of 50, 100 and 150, respectively. The trips 4 and 5 are departing from depot 2, and have arrival times of 60 and 90.
Suppose six trips are leaving from satellite 1. Table 1, describes the quantities required from depot 1 and from depot 2, and the departure time from the satellite (the satellite has a service time of 10). Then, the last 5 columns present the d sl -values for the specific relevant combinations of first level trips.

A combined arc-path formulation (MC2E-A-P)
Let L s z ⊂ L denote the subset of city freighter tours which start from the satellite s and visit the customer z. Moreover, the constant ξ l ij = 1 indicates that the city freighter l visits the customer j immediately after visiting customer i. Similarly, the constant π l ij = 1 indicates that the city freighter l visits both customers i and j, but the customer i is visited as the first customer in the tour. -Commodity-urban vehicle assignment: y τ z = 1 if demand of customer z is assigned to vehicle τ ∈ T ; -Timing: ω τ s shows the arrival time of the urban vehicle τ at satellite s.
• Second echelon routing (path-based): -City freighter tour: δ l = 1 if city freighter tour l is selected; 0 otherwise; -Timing: ω z shows the arrival time at customer z.
• Connecting the first and second echelon routes: -Commodity-satellite assignment: q s z = 1 if the demand of the customer z is satisfied through satellite s.
The objective function (37) minimizes the total transportation and vehicle usage cost at both echelons. Constraints (38-44) are similar as the MC2E-2A formulation.
Constraint (45) makes sure that each commodity is handled through a satellite. Constraint (46) relates the commodity-satellite assignment variables with the second echelon routing variables. It assures that if a commodity is assigned to a satellite, then one city freighter tour should depart from the satellite and deliver the commodity to the associated customer. Constraint (47) relates the arrival time at each customer with the arrival time of an urban vehicle at a satellite, if the customer is visited as the first customer in a city freighter tour and the urban vehicle feeds the city freighter.
Constraint (48) relates the arrival times at pair of customers, if one of the customers is visited immediately after visiting the other customer by a city freighter. Constraints (49) and (50) are hard time windows constraints. Constraints (51-56) are domain constraints.
The MC2E-A-P formulation encompasses a large number of city freighter tours to be generated. A column generation method combined with a branch-and-bound approach (branchand-price) is used to solve this formulation.

Branch-and-Price based solution approach for the MC2E-2P model formulation
The MC2E-2P problem formulation can be solved by a smart enumeration of the first level combinations and using a branch-and-price approach, using ψ s and ϕ z as the dual variables respectively associated to Constraints (34) and (35). After this reformulation, we can use the approach from Dellaert et al. (2019). In the Appendix, we give a basic solution approach for the MC2E-A-P formulation.
For each of the first echelon configurations, we derive a lower bound for the total costs by considering the exact first echelon costs and a dual-cost-based lower bound estimation for the second echelon, based upon aggregate characteristics of the first echelon configuration. We sort the first echelon configurations in increasing costs and step-wise add them to the branch-and-bound tree, whenever the lower bound of the branch-and-bound tree solution exceeds the lower bound of the next first echelon configuration. At each node of the branchand-bound tree, a relaxed master problem (RMP) is obtained by considering the new set of city freighter tours L at the LP-relaxation of the MC2E-2P formulation (33-36).
An insertion algorithm is used to solve the pricing problem to generate negative reduced cost city-freighter tours. Multiple branching strategies are used to create the branch-andbound tree. We use a best-first strategy, which selects the branch with the best lower bound value. Algorithm 1 summarizes the procedure and each step is described in detail in the following sections. Set the LB equal to the lowest lower bound value of the leaf nodes of the tree T

13:
If an integer solution is found: check assignment and update U B

14:
Add a new first echelon configuration when its lower bound< LB.

Evaluating the LP (Lines (3-7))
The column generation scheme searches for the city freighter tour with the most negative reduced cost by solving a pricing problem. Using ψ s and ϕ z as the dual variables associated to Constraints (34) and (35), respectively, and r as the first echelon configuration of the node, the pricing problem is formulated as follows: City freighter tours with negative reduced cost values are added to L and the RMP with new L is re-solved. This is repeated until no city freighter tour with negative reduced cost exists.
The pricing problem (58) uses all information about the schedule of the city freighter (departure time from the satellite, arrival and departure times at customers). Thus, only tours that are feasible for this specific first level configuration r are considered. By considering all relevant combinations of first level trips and satellites one by one, a search space explosion is avoided. In our insertion procedure we first rank the potential customers in decreasing dual costs. This is also the insertion sequence. We also explicitly calculate the savings potential for every combination of remaining capacity and last inserted customer and using this we can easily prune solutions that cannot lead to the most negative reduced cost solution.
Branching strategy (Lines (9-11)) As usual, we use branching to deal with fractional values of the decision variables. We deal with multiple types of decision variables. Therefore, the branching is done on: • Combination of two commodities in a tour, • Commodity-satellite assignment variables, • Fixing one complete tour.
In our strong branching, we first select good candidates, based on fractional solutions close to 0.5, and for the 12 potentially strongest candidates, we determine the improvement to the objective function by solving each of the corresponding LP-relaxations with the existing columns. Ultimately, the one with the largest improvement is chosen.

Issues with larger problems
When we deal with larger problems (i.e., 50 or more commodities, 5 or more satellites), it is important that we have good lower bounds for the first level configurations, as there are many potential candidates. We improve the lower bounds by explicitly considering aggregate characteristics, like the earliest departure times from the satellites and the potential commodities from other depots in every configuration. It is important to keep a good balance between improving the branch-and-bound tree LB, adding new configurations, and searching for a better upper bound. To find a better upper bound, we diversify the search to consider the node with the lowest number of fractional solutions (rather than lowest costs). Of course, finding optimal solutions is not always possible for larger problem sizes and heuristic solutions become inevitable.

Computational study
We present the results of computational tests on newly generated instances of the problem using the proposed models. All tests were performed using C++ on a PC with Intel Core i7-4770 Processor (3.40 GHz) and 8.00 GB of RAM. CPLEX v12.6.0 was used to solve the arc-based formulation and the LP-relaxations. We report the computation time and the bounds to show the performance of our algorithms.
No instances are available in the literature for the MC-2E-VRPTW. We modify instances of the 2E-VRPTW from the literature to generate MC-2E-VRPTW instances. In order to introduce commodities, a random depot is assigned to each customer (commodity). There is one set of instances for the 2E-VRPTW in the literature proposed by Dellaert et al. (2019). These instances are generated by simulating an urban area.
There are four categories of instances differing in the time windows and demands of customers: Ca, Cb, Cc and Cd. For all instances, urban vehicle capacity, city freighter capacity, urban vehicle cost and city freighter cost are set to K (1) = 200, K (2) = 50, h (1) = 50, h (2) = 25. Moreover, the service time for all customers and satellites is set to 10.
We keep the notation of Dellaert et al. (2019) to address instances, but add "MC" to the start of the instance names.

Performance of the MC2E-2A formulation
We test the 2-arc model on instances with 2 depots, 3 satellites and 15 customers, due to the long CPU time required to solve all other larger instances. Table 2 gives the results for the MC2E-2A formulation using CPLEX for these instances. We impose a time limit of one hour. As expected with arc-based formulations, CPLEX fails to solve any of the instances to optimality and is not able to generate any upper bound for 3 instances.

The MC2E-2P model formulation and solution approach
We consider the MC2E-2P formulation and solution approach to work on some larger problems up to 100 customers (|Z| = {50, 100}).
From Table 3 we learn that most instances with 50 customers can be solved to optimality, some of them even within a few minutes. Most of the cases with 100 customers and 2 depots can be solved close to optimality, but when we have 3 depots and 5 satellites we have a maximum gap of almost 4 percent.
In the last column, labeled FLC we give the number of First Level Combinations that are added to the branch-and-bound tree. We see that this number can be really huge for the case with 5 satellites. In the case of 50 customers, it takes less time to add a first level combination than in the case of 100 customers and, therefore, we can add more combinations in our limited computation time. The large number of combinations that are added to the tree is completely different from the single commodity situation, as considered in Dellaert et al. (2019). When there is only a single commodity, which can be supplied by all depots, most first level tours in the optimal situation only visit one (nearby) satellite, whereas in the multi-commodity situation most of the first level tours visit multiple satellites in order to make better second level tours.

The effect of commodities
In this section, we describe the effect of single versus multiple commodities. We evaluate first in terms of the effect on the solution method and, secondly, in terms of costs.
In general, the solution time increases in the multi-commodity situation. Especially in cases with more depots, the number of first level trips increase, as there is no pooling between depots. In the optimal solution, we see far more multi-satellite routes in the first level. We thus need to consider many more first level combinations because of the bigger gap between initial lower bound and optimal solution. We also noticed this effect in the arc-path formulation, where it was even stronger. For the second level routes, we have to consider more satellite options, with different departure times, although often with fewer customers. The percentage of trips with customers from multiple depots is between 30 and 70 percent. In all situations a minimum number of vehicles is used at the first level, and at the second level only occasionally an extra vehicle is used. All in all, we can say that the multi-commodity situation with 3 or more depots is tougher than the single commodity problem.
In terms of costs, we observe a 13 percent total cost increase for the 2 depot situation and 27 percent for the 3 depot situation. The cost increase is the largest for the first level travel costs. Table 4 obtains the percentage cost increases for the four cost categories.

Conclusions and Future Research
In this paper, we introduced the MC-2E-VRPTW which deals with commodities to manage the indirect shipment of non-substitutable demands from specific origin depots to final customer through satellites.
We modeled this problem in different ways. We designed a 3-path formulation which defines urban vehicle tours, city freighter tours and itineraries as the three types of path variables. This formulation is a compact one but has big challenges when column generation is used to solve LP-relaxations. The reduced cost formulations rely on unknown dual variables which seems to make this formulation not suitable for branch-and-price. We de- composed the problem into the first and second echelon routing problems and replaced the itineraries by constraints securing the synchronization. Alternative formulations are discussed. A 2-arc formulation is proposed where arc-flow variables are used to represent the routes at both echelons. We also introduced a combined arc-path formulation which uses an arc-based modeling approach for the first echelon and a path-based modeling approach for the second echelon and lets constraints connect and synchronize the vehicle routes of the first and second echelons. Finally, we elaborated on a 2 path formulation, where again constraints do the flow synchronization. We designed branchand-price algorithms. The 2-path model and solution approach proved to be efficient and effective. We succeeded to solve instances with up to 5 satellites and 100 customers to near-optimality.
The proposed models are essential for further research. In terms of modeling, extensions of the proposed models could be studied to capture more realistic aspects (e.g., uncertainty in problem parameters, time-dependent travel times, etc.). In terms of solution approach, a branch-and-cut-and-price method could be investigated. Moreover, a branch-and-price method could be developed for the 3-path formulation by estimating the unknown dual variables. Alternatively, a simultaneous column-and-row generation mechanism seems promising as well (Muter et al., 2013). In addition, math-heuristic and meta-heuristic methods could be studied to solve bigger problem instances.