Design of a two-echelon last-mile delivery model

Pina-Pardo, Juan C. et al. Article Design of a two-echelon last-mile delivery model EURO Journal on Transportation and Logistics (EJTL) Provided in Cooperation with: Association of European Operational Research Societies (EURO), Fribourg Suggested Citation: Pina-Pardo, Juan C. et al. (2022) : Design of a two-echelon last-mile delivery model, EURO Journal on Transportation and Logistics (EJTL), ISSN 2192-4384, Elsevier, Amsterdam, Vol. 11, Iss. 1, pp. 1-13, https://doi.org/10.1016/j.ejtl.2022.100079 This Version is available at: https://hdl.handle.net/10419/325156 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-nc-nd/4.0/ EURO Journal on Transportation and Logistics 11 (2022) 100079 Available online 20 April 2022 2192-4376/© 2022 The Authors. 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-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). Design of a two-echelon last-mile delivery model Juan C. Pina-Pardo a , * , Matheo Moreno b , Miguel Barros c , Alexandre Faria d , Matthias Winkenbach e , Milena Janjevic e a School of Industrial Engineering, Pontificia Universidad Cat´ olica de Valparaíso, Chile b Business School, Universidad San Francisco de Quito, Ecuador c Institute for Educational Development, Fundaç˜ ao Getulio Vargas, Brazil d Executive MBA Program, Insper, Brazil e Center for Transportation & Logistics, Massachusetts Institute of Technology, USA ARTICLE INFO Keywords: Network design Facility location Last-mile Two-stage stochastic program Two-echelon Vehicle routing ABSTRACT Due to high congestion in cities and growing demand for last-mile delivery services, several companies have been implementing two-echelon distribution strategies over the past few years. Notably, the installation of urban transshipment points has gained increasing attention, used by logistics operators to transfer goods from large freight trucks to smaller and more agile vehicles for last-mile delivery. Nevertheless, the main challenge is how to decide the number and location of these facilities under the presence of demand uncertainty. In this paper, we develop a two-stage stochastic program to design two-echelon last-mile delivery networks under demand uncertainty. This approach decomposes the problem into strategic decisions (facility location) and operational decisions (daily distribution of goods). To address large-scale instances, we solve the model through the sample average approximation (SAA) technique and estimate the optimal routing costs (of the SAA counterpart) using a continuous approximation method. Using a real-world case study with more than 1300 customers from New York City, our results provide several managerial insights regarding the mix of transportation modes, facility location, and the impact of allowing the outsourcing of customer demand. We provide extensive validation of the twostage stochastic program results through a simulation-based approach and the calculation of the value of the stochastic solutions. 1. Introduction Retail e-commerce sales worldwide are expected to reach 4.9 trillion U.S. dollars by 2021, representing an increase of 268% from 2014 (Statista, 2020). This is accompanied by the increasing expectation of consumers for fast delivery. According to the World Economic Forum (2020), instant and same-day delivery services are the fastest-growing segments in the last-mile context, with growth rates at nearly 30% per year. This trend puts high pressure on retailers and e-tailers, since the last-mile delivery is the most complex, expensive, and time-consuming stage of the e-commerce logistics chain (Business Insider, 2018). The growth of cities in recent years has also made last-mile logistics more complex (World Economic Forum, 2020). The world population is projected to be 9.8 billion by 2050, and 68% of people are expected to live in cities (Ritchie, 2018). The swelling population of urban areas puts stress on existing infrastructure and augments traffic congestion and air pollution, reducing people’s quality of life. Without any intervention, traffic congestion is estimated to increase by 21% by 2030, which is equivalent to adding 11 minutes of commute time per passenger per day (World Economic Forum, 2020). On a similar note, a recent study estimates that looking for parking accounts for about 30% of the total travel time of delivery vehicles in downtown Seattle due to congestion (Dalla Chiara and Goodchild, 2020). To efficiently respond to the demand for fast deliveries in increasingly complex urban settings, companies have been designing and implementing multi-echelon delivery networks over the past several years (Janjevic and Ndiaye, 2014). As Cuda et al. (2015) mentioned in their work, these networks represent an effective response to congestion, pollution emissions, and noise generated by freight transportation in urban areas. Examples of these include the use of micro-hubs or urban transshipment points (so-called satellite facilities), which are used by logistics companies to transfer goods from freight trucks to small and eco-friendly vehicles for last-mile deliveries (Urban Freight Lab, 2002). In this context, this paper studies the design of two-echelon last-mile * Corresponding author. Avenida Brasil 2241, Valparaíso, 2362807, Chile. E-mail address: [email protected] (J.C. Pina-Pardo). Contents lists available at ScienceDirect EURO Journal on Transportation and Logistics journal homepage: www.sciencedirect.com/journal/euro-journal-on-transportation-and-logistics https://doi.org/10.1016/j.ejtl.2022.100079 Received 6 August 2021; Received in revised form 31 March 2022; Accepted 7 April 2022 EURO Journal on Transportation and Logistics 11 (2022) 100079 2 delivery networks, composed of distribution centers (DCs) and urban transshipment points (UTPs). In these distribution systems, customer products are transported in large vehicles from DCs (typically located in suburban areas) to different UTPs located throughout the city, from where they are collected by smaller and more agile vehicles (e.g., vans or cargo bikes) for delivery. The problem studied consists of determining the locations where DCs and UTPs should be installed to minimize the total installation and transportation costs under demand uncertainty. An illustrative example of the two-echelon distribution scheme studied is shown in Fig. 1. The shaded triangles and boxes represent the installed DCs and UTPs, respectively. Continuous lines show the flow of customer products in the first distribution echelon. In the second distribution echelon, continuous lines represent the shipment of products using vans, while dashed lines show the shipment of products using cargo bikes. These vehicles are just examples of the transportation means that each UTP could use for deliveries. According to their geographical locations, end customers are clustered into several segments. The size of each customer segment is related to the corresponding customer density. The problem studied is known in the literature as the Two-Echelon Location-Routing Problem (2E-LRP) (Cuda et al., 2015; Janjevic et al., 2019). The term is because it integrates two well-known problems in operations research: the Facility Location Problem (FLP) and the Vehicle Routing Problem (VRP). The FLP aims to find the optimal locations where facilities must be opened to serve a set of customers (Farahani and Hekmatfar, 2009), whereas the VRP seeks to determine the optimal set of routes that a fleet of vehicles must perform to deliver customer orders (Toth and Vigo, 2014). Note that, if we assume the existence of a single (previously installed) DC and that the deliveries are made directly from that location to the customers (i.e., without using UTPs), the solution of the problem will be given by performing an optimal VRP set of routes. Therefore, the 2E-LRP is 𝒩 ℘-hard since it has the VRP as a special case. Due to the inherent complexity of the 2E-LRP, we use the continuous approximation function developed by Winkenbach et al. (2016a) to estimate the optimal routing costs to serve each customer segment (see Fig. 1), thus avoiding defining explicit vehicle routes to visit each individual consumer. Furthermore, to consider demand uncertainty, we develop a two-stage stochastic program, where we decompose the problem into facility location and distribution decisions. This model is solved using the sample average approximation technique, which has been successfully applied to network design problems under uncertainty (Li and Zhang, 2018; Santoso et al., 2005). Using a real-world case study with more than 1300 consumers from Manhattan, NY, our results provide several managerial recommendations regarding the mix of transportation modes, facility location, and the impact of allowing the outsourcing of customer demand. Additionally, we show the importance of considering demand uncertainty when designing two-echelon distribution networks by computing the value of the stochastic solutions and through a simulation-based optimization approach. The rest of the paper is structured as follows. Related works to our stochastic 2E-LRP are presented in Section 2. In Section 3, we formally describe the problem studied and the two-stage stochastic program developed. Section 4 presents the SAA algorithm used to solve the twostage stochastic program. The case study is described in Section 5, whereas computational results are presented and analyzed in Section 6. Finally, conclusions and recommendations for further developments are discussed in Section 7. 2. Literature review Location-Routing Problems (LRPs) arise to integrate both facility location and routing decisions, which are commonly made in a sequential or hierarchical fashion (i.e., routing decisions are typically subordinate to location decisions) (Nagy and Salhi, 2007). As Nagy and Salhi (2007) pointed out, the main argument for separating these decisions is that facility location is a long-term strategic decision, while routing is a more tactical or operational decision that is commonly made on a daily basis. Nevertheless, as Salhi and Rand (1989) showed, making these decisions separately may lead to highly suboptimal solutions. Prodhon and Prins (2014) conducted a survey where they presented different exact and heuristic solution approaches for the LRP. Additionally, the effectiveness of the identified metaheuristics is compared on problems that consider capacitated vehicles, as well as two distribution echelons. Finally, these authors described new extensions of the LRP that consider several decision periods, inventory management decisions, and uncertain data. Other recent surveys of the LRP can be found in Drexl and Schneider (2015), and in Schneider and Drexl (2017). A well-known variant of the LRP emerges when considering twoechelon delivery networks, which are employed by logistics companies to distribute their products to end consumers through intermediate facilities (Guastaroba et al., 2016). Cuda et al. (2015) provided a survey on two-echelon routing problems. These authors classified these problems into three categories: the 2E-LRPs, where both the location of facilities and the routing of vehicles must be decided; the two-echelon Fig. 1. Example of the two-echelon delivery network studied. J.C. Pina-Pardo et al. EURO Journal on Transportation and Logistics 11 (2022) 100079 3 VRPs, where the set of facilities is given and the vehicle routing must be specified at each tier of the network; and the truck and trailer routing problems, where there are customers that can be visited by either a truck pulling a trailer or by a truck alone, and others that can exclusively be served by a truck alone (this means that trailers must be parked at some location so that vehicles can perform the second distribution stage alone). For each class of problem, Cuda et al. (2015) described route-based mathematical formulations and organized the existing solution methods in the literature. To solve large-scale 2E-LRP instances, defining explicit routing to serve each individual end consumer becomes computationally intractable (Janjevic et al., 2021). Therefore, to approximate the optimal routing costs, the use of route length estimation methods has gained increasing attention in the literature over the past few years (Nagy and Salhi, 2007; Smilowitz and Daganzo, 2007). Winkenbach et al. (2016a) developed an Augmented Routing Cost Estimation (ARCE) function that extends the approximation formula developed by Smilowitz and Daganzo (2007) to consider multiple vehicle options, vehicle capacity constraints, and a maximum service time restriction. These authors proposed a solution scheme that can be briefly summarized in locating first and routing second, where customers are clustered and the ARCE function is employed to approximate the cost of serving each customer cluster. Using real-world data provided by La Poste, Winkenbach et al. (2016a) showed that this two-stage heuristic approach can produce high-quality solutions while substantially reducing computational times. Their proposed method subsequently informed the strategic redesign of La Poste’s mail and parcel networks (Winkenbach et al., 2016b). Recently, Janjevic et al. (2019) proposed a two-tier delivery network where customer products are first shipped from a single warehouse to different satellite facilities. The products are then collected by last-mile delivery vehicles, which deliver them directly to end consumers or to pick-up locations (such as parcel lockers). They formulated the problem as a non-linear optimization problem, where the routing cost component is approximated using an adaptation of the ARCE function proposed by Winkenbach et al. (2016a). Their solution method is based on enumerating all possible configuration/subset of satellites and then, for each configuration, solving the resulting problem using a constructive heuristic approach where pick-up locations are iteratively added. Similarly, Rautela et al. (2021) developed a 2E-LRP to investigate the impact of collection-and-delivery points, such as parcel lockers, on the overall cost performance of the last-mile delivery system. Finally, regarding problems that consider stochastic demand, Snoeck et al. (2018) proposed a two-stage stochastic program for solving a 2E-LRP with facility and vehicle capacity constraints. In this formulation, the first-stage optimization model aims to determine the strategic decisions – related to facility renting and vehicle acquisition – for all possible demand scenarios (i.e., for all realizations of the random variables that describe the demand), while the second-stage model seeks to find the optimal operational decisions – associated with the daily distribution of customer products – for a given demand scenario with fixed strategic decisions. To address the problem, the authors divided the distribution area into several rectangular pixels, developing an adaptation of the ARCE formula of Winkenbach et al. (2016a) to approximate the cost of serving each pixel. Finally, their solution approach consisted of solving the first-stage model on a representative subset of demand scenarios and then obtaining the network design candidates that minimize the expected total cost. For the interested reader, a similar study can be found in Snoeck and Winkenbach (2020). Additionally, a recent methodology for solving multi-period 2E-LRPs under uncertainty is described in Mohamed et al. (2020), who also introduced a two-stage stochastic program that is solved through Benders decomposition and the Sample Approximation Approach (SAA) technique (see Ahmed and Shapiro (2002) for a detailed description of the SAA method). 2.1. Research gaps The vast majority of the literature has addressed the design of twoechelon distribution networks from a deterministic perspective (Cuda et al., 2015; Guastaroba et al., 2016). The different exact methods that have been developed, however, fail to address realistic-sized instances, primarily because defining explicit routes to serve each end consumer becomes impossible in reasonable times using state-of-the-art commercial solvers. As identified by Janjevic et al. (2021), instances of up to 500 customers have been heuristically addressed. Nevertheless, real-world problems often exceed 100,000 customers (Snoeck and Winkenbach, 2020). Furthermore, to the best of our knowledge, only Snoeck et al. (2018) and Snoeck and Winkenbach (2020) have studied the design of large-scale last-mile networks under uncertainty. Therefore, in light of the aforementioned research gaps, the contributions of this paper are: •We develop a two-stage stochastic program to design two-echelon last-mile delivery networks under demand uncertainty. Unlike Snoeck et al. (2018) and Snoeck and Winkenbach (2020), our stochastic model explicitly considers the installation of distribution centers and the flow of customer goods in the first distribution echelon, considering both the transportation costs (between the installed facilities) and the flow balance constraints involved. •The stochastic model is successfully solved through the SAA technique. We provide an analysis of the solution quality and the computational effort required. •Based on deterministic counterparts of the two-stage stochastic program, an extensive analysis is carried out to validate its results and show the importance of considering demand uncertainty when solving the problem. •We provide several managerial insights by using a real-world case study with more than 1300 consumers from Manhattan, NY. 3. Mathematical formulations In this section, we begin by presenting a Mixed-Integer Linear Programming (MILP) formulation for the deterministic case of the problem, where each customer’s demand is a known parameter. We then describe the two-stage stochastic program, where we decompose the problem into strategic decisions (facility location) and operational decisions (daily distribution of goods). Regarding our modeling considerations and assumptions, we consider that each UTP has a set of available vehicle types (e.g., vans and cargo bikes) for making deliveries (see Fig. 1). The number of each type of vehicle at each UTP is assumed to be unlimited. Furthermore, we assume that all customers order the same product but in different quantities. To ensure the computational tractability of the mathematical optimization models developed, we cluster the customers into several segments and use the ARCE function developed by Winkenbach et al. (2016a) to approximate the optimal routing costs to serve each of them (see Appendix A for a detailed description of the ARCE function). The quality of this function (to obtain optimal or near-optimal values) was exhaustively validated in Winkenbach et al. (2016a) and Janjevic et al. (2019), and therefore outside our scope. 3.1. Deterministic MILP formulation This section presents a deterministic MILP formulation to introduce the mathematical notation and explain the main decisions involved in the problem studied. In this model, we assume that the demand of each customer segment is a known parameter. Using a simulation-based approach, we will employ this MILP model to validate the results of our two-stage stochastic program (see Section 6.2). Let K be the set of customer segments, I the set of locations where a DC can be opened, and J the set of locations where a UTP can be opened. J.C. Pina-Pardo et al. EURO Journal on Transportation and Logistics 11 (2022) 100079 4 Furthermore, let V be the set of available vehicle types to deliver products from UTPs to customer segments. The parameters and decision variables of the model are listed below. In particular, the parameter cv jk is computed using the ARCE function created by Winkenbach et al. (2016a). 3.1.1. Parameters •c ij : Cost for shipping products from DC i ∈I to UTP j ∈J. •cv jk: Cost of serving the demand of customer segment k ∈K from UTP j ∈J using vehicles of type v ∈V. •f i : Fixed cost by opening a DC at location i ∈I. •g j : Fixed cost by opening a UTP at location j ∈J. •d k : Units of products demanded by customer segment k ∈K. •D i : Capacity of DC i ∈I. •S j : Capacity of UTP j ∈J. •Q: Load capacity of first-echelon vehicles. 3.1.2. Decision variables •y i : 1, if a DC is opened at location i ∈I. 0, otherwise. •γ j : 1, if a UTP is opened at location j ∈J. 0, otherwise. •x ij : Number of first-echelon vehicles used to ship products from DC i ∈I to UTP j ∈J. •r ij : Number of products sent from DC i ∈I to UTP j ∈J. •sv jk: 1, if UTP j ∈J serves customer segment k ∈K using vehicles of type v ∈V. 0, otherwise. 3.1.3. MILP formulation The MILP formulation is described as follows. minimize ∑ i∈I fi⋅yi+∑ j∈J gj⋅γj+∑ i∈I∑ j∈J cij⋅xij +∑ j∈J∑ k∈K∑ v∈Vcv jk⋅sv jk,(1) subject to ∑ j∈J rij ≤Di⋅yi,∀i∈I,(2) ∑ i∈I rij ≤Sj⋅γj,∀j∈J,(3) ∑ i∈I rij =∑ k∈K∑ v∈V dk⋅sv jk,∀j∈J,(4) ∑ j∈J∑ v∈V sv jk =1,∀k∈K,(5) rij ≤Q⋅xij,∀i∈I,j∈J,(6) yi∈ {0,1},∀i∈I,(7) γj∈ {0,1},∀j∈J,(8) rij ≥0,∀i∈I,j∈J,(9) xij ∈Z+ ,∀i∈I,j∈J,(10) sv jk ∈ {0,1},∀j∈J,k∈K,v∈V.(11) The objective function (1) minimizes the total installation and transportation costs. Constraints (2) state that a DC can send products as long as it has been opened. Similarly, constraints (3) force each UTP to be open in case of receiving products from DCs. The flow balance constraints of each UTP are included in (4). Constraints (5) force each customer segment to be served by exactly one UTP using exactly one type of vehicle. The variables r ij and x ij are linked in constraints (6). Finally, non-negativity and integrality constraints are included in (7)– (11). Note that, if a customer segment can be served by more than one UTP using more than vehicle type, then constraints (11) should be replaced by sv jk ∈ [0,1],∀j∈J,k∈K,v∈V. In this case, each variable sv jk represents the percentage of demand of the customer segment k that is delivered from UTP j using vehicles of type v. The rest of the MILP model remains the same. 3.2. Stochastic model Let us now consider the case of stochastic demand. We define a realization of the daily demand of all customer segments as a scenario ω ∈Ω, where Ω represents a set that contains all possible scenarios. For each customer segment k ∈K and scenario ω ∈Ω, we denote its daily demand as d k ( ω ). Additionally, we consider that each scenario ω has a probability of occurrence π ( ω ). The approach used to model the stochastic 2E-LRP is inspired by Snoeck and Winkenbach (2020). In this approach, we decompose the problem into the strategic decisions – related to the location of facilities – and the operational decisions, associated with the daily distribution of products in each stage of the delivery network. Strategic decisions are made considering all possible demand scenarios, while operational decisions are determined for each scenario ω ∈Ω individually. This results in a two-stage stochastic program, which is described below. 3.2.1. First-stage MILP model Let y:= (yi)i∈I∈ {0,1}I and γ:= (γj)j∈J∈ {0,1}J be vectors of binary variables that indicate the locations where DCs and UTPs will be opened, respectively. Given y and γ, we define the function C( ω , y, γ) as the minimum operational cost to satisfy the demand described by ω ∈Ω. The first-stage MILP model aims to determine the values of y ∈{0,1} I and γ ∈{0,1} J that minimizes the sum of installation costs and the expected distribution cost for all demand scenarios, as presented in (12). min{f⊤y+g⊤γ+∑ ω ∈Ω π ( ω )⋅C( ω ,y,γ):(y,γ) ∈ {0,1}I× {0,1}J}(12) 3.2.2. Second-stage MILP model Given that the total daily demand may exceed the capacity of all opened UTPs, this model considers that each customer segment can also be served directly from DCs. In case this is not enough, we assume that the company pays a penalty p k for each unit of product not delivered to customer segment k ∈K. (Alternatively, a similar assumption is that the company purchases the remaining products from a third party, which also handles their distribution.) Furthermore, since the cost of serving each customer segment’s demand depends on the particular demand scenario, we update the parameters c ik and cv jk by c ik ( ω ) and cv jk( ω ), respectively. The new decision variables, as well as the redefinition of variables sv jk, are described as follows. 3.2.2.1. Decision variables. • τ ik : Fraction of the demand of customer segment k ∈K sent from DC i ∈I. J.C. Pina-Pardo et al. EURO Journal on Transportation and Logistics 11 (2022) 100079 5 •β k : Fraction of the demand of customer segment k ∈K not delivered or outsourced. •sv jk: Fraction of the demand of the customer segment k ∈K that is delivered from UTP j ∈J using vehicles of type v ∈V. 3.2.2.2. MILP formulation. Given y ∈{0,1} I and γ ∈{0,1} J , as well as a scenario ω ∈Ω, the second-stage MILP model seeks to minimize the distribution and penalty (per lost or outsourced demand) costs, as presented below. C( ω ,y,γ) = min ∑ i∈I∑ j∈J cij⋅xij +∑ i∈I∑ k∈K cik( ω )⋅ τ ik +∑ j∈J∑ k∈K∑ v∈Vcv jk( ω )⋅sv jk +∑ k∈K pk⋅dk( ω )⋅βk, (13) subject to (3), (6), (9), (10), ∑ j∈J rij +∑ k∈K dk( ω )⋅ τ ik ≤Di⋅yi,∀i∈I,(14) ∑ i∈I rij =∑ k∈K∑ v∈V dk( ω )⋅sv jk,∀j∈J,(15) ∑ j∈J∑ v∈V sv jk +∑ i∈I τ ik +βk=1,∀k∈K,(16) sv jk ∈ [0,1],∀j∈J,k∈K,v∈V,(17) τ ik ∈ [0,1],∀i∈I,k∈K,(18) βk∈ [0,1],∀k∈K.(19) Constraints (14) extend constraints (2) by considering that DCs can deliver products directly to customer segments. Similarly, constraints (5) are extended in (16), taking into consideration that the demand of each customer segment is now divided into the fraction delivered by UTPs, the fraction provided by DCs, and the fraction of lost (or outsourced) demand. The flow balance constraints associated with each UTP are presented in (15). Finally, the domain of the new decision variables is specified in (17)–(19). 4. Solution approach To make the two-stage stochastic program tractable and solve it in reasonable computational times, we employ the SAA technique, which has been successfully applied to supply chain network design problems under uncertainty (Li and Zhang, 2018; Santoso et al., 2005). In the following, we provide a theoretical overview of the SAA technique and then present the algorithm used to solve the 2E-LRP studied. 4.1. Sample average approximation Following the SAA technique, the expected distribution cost of Equation (12) is approximated (i.e., estimated) by the average distribution cost of a sample { ω n}N n=1 of N scenarios, as presented in the corresponding SAA problem, min{f⊤y+g⊤γ+1 N∑N n=1C( ω n,y,γ):(y,γ) ∈ {0,1}I× {0,1}J}.(20) Let us denote by o* and oN the corresponding optimal values of (12) and (20), respectively. As discussed in several works, such as Ahmed and Shapiro (2002) and Kleywegt et al. (2002), oN is a consistent estimator of o*. That is, the optimal solution value to the SAA problem (20) converges with probability one to the optimal solution value of (12) as the sample size N increases. Further, it is also well known that E[oN] ≤ o∗ (Mak et al., 1999). Therefore, it is possible to compute a lower bound to the true optimal solution value o* by estimating the expected value of oN, as explained below. 4.1.1. Computing lower bounds Let us consider that M independent samples are generated, each containing N scenarios. Let om N be the optimal solution value and (ym N,γm N) an optimal solution of the corresponding SAA problem, with m ∈{1, …, M}. Since the quantity oM N:=1 M∑M m=1om N(21) is an unbiased estimator of E[oN](Ahmed and Shapiro, 2002), it follows that oM N represents a (statistical) lower bound to o*. 4.1.2. Computing upper bounds Consider now a feasible solution of (12), for example, one of the previously generated solutions (ym N,γm N), with m ∈{1, …, M}. The expected distribution cost of implementing this solution can be estimated by generating a new (independent) sample { ω n}N′ n=1 of N′scenarios. Thereby, the quantity υ m N′:=f⊤ym N+g⊤γm N+1 N′∑N′ n=1C( ω n,ym N,γm N)(22) is an unbiased estimator of the true objective value of implementing (ym N, γm N). Since (ym N,γm N)is a feasible solution of (12), it follows that (22) sets a (statistical) upper bound on o*. Note that computing (22) requires solving N′independent secondstage MILP problems. However, given that the facility location decisions are fixed, in practice, it is possible to use a sample size N′much bigger than N (Santoso et al., 2005). 4.2. SAA scheme for the stochastic 2E-LRP As discussed by Schütz et al. (2009), Santoso et al. (2005), and others, repeatedly solving the SAA problem (20) for M independent samples can be more (computationally) efficient than increasing the sample size N, which we also noticed in our preliminary experiments. Moreover, this procedure also allows computing an estimator of the optimality gap, which can be useful for decision-makers to have information about the quality of a given candidate solution. Thus, the proposed SAA scheme for solving the stochastic 2E-LRP is presented in Algorithm 1. In Steps 1 to 3 of Algorithm 1, M SAA problems are solved, each containing an independent sample of N scenarios. Then, the statistical lower bound oM N is computed in Step 4 (see Section 4.1.1). In Step 5, a new set of N′scenarios is created to obtain the statistical upper bound. To do so, in Steps 6 and 7, the solution of each SAA problem is fixed and then N′independent second-stage MILP models are solved (see Section 4.1.2). Following Santoso et al. (2005), a network design with the smallest upper bound is selected in Step 8. In Step 9, the optimality gap estimate is computed. The algorithm terminates if the optimality gap is small enough; otherwise, we increase M, N, or N′and return to Step 1 (Bidhandi and Patrick, 2017). J.C. Pina-Pardo et al. EURO Journal on Transportation and Logistics 11 (2022) 100079 6 Algorithm 1. The SAA algorithm for the stochastic 2E-LRP 5. Case study This section presents the case study used for evaluating the performance of Algorithm 1. We begin by describing the data and our parameter choices. Then, we present the scenario generation procedure used. 5.1. Data and parameter definition In our experiments, we used a sample dataset containing 5-day demand information from 1385 end consumers from Manhattan, NY (the information about the company and the demanded products is not revealed in this paper to protect confidentiality). For each end consumer, this dataset shows the geographical location (latitude and longitude), zip code, the number of units demanded, and the day their order was placed. We obtained the zip codes information of Manhattan from NYC Open Data (2018) to cluster the end consumers into different segments. For the sake of exposition, we discarded those zip codes whose area is less than 0.1 [km 2 ] (e.g., individual buildings). Fig. 2 shows the resulting 40 zip codes (i.e., customer segments) employed, where the intensity of the color is related to the size of each zip code’s area. Based on the information provided by the company, we considered 3 and 8 candidate locations for opening DCs and UTPs, respectively. Fig. 3 depicts these locations and the customer segments considered (the size of each segment is related to its demand). We used road network information from Open Street Maps to obtain the distances between each node of the underlying graph. We assumed that DCs and UTPs are homogeneous. The fixed cost for opening a DC is f i =1000 cost units and its capacity is D i =3000 parcels, ∀i ∈I. For each UTP j ∈J, we considered a fixed cost of g j =400 cost units and a capacity of S j =1000 parcels. Finally, for the first distribution echelon, we considered a homogeneous fleet of trucks with a load capacity of 250 parcels. We set the fixed cost for using a truck as 200 cost units. In the second echelon, we employed vans and cargo bikes for last-mile deliveries. Based on Winkenbach et al. (2016a), the parameters associated with the first and second echelon vehicles are presented in Table 1, used as input to the ARCE function to approximate the optimal routing cost to serve each customer segment. 5.2. Scenario generation To generate demand scenarios, we assumed that each customer segment’s demand follows a normal distribution. We computed the mean and standard deviation for each segment using the five days of demand presented in our database. With this, we generated a total of 500 scenarios (i.e., 500 realizations of the random daily demand of all customer segments), from which the samples of Steps 2 and 5 of Algorithm 1 are obtained. Regarding Algorithm 1, we solved M =20 SAA problems to obtain the statistical lower bound for our two-stage stochastic program (see Steps 1 to 3). For the SAA problems, we considered sample sizes of N ∈{10, 20, 30} scenarios. Furthermore, to compute the statistical upper bound and thus calculate the estimator for the optimality gap, we used a sample of N′=250 scenarios (see Steps 5 to 9). Similar values were used by Santoso et al. (2005) and Snoeck et al. (2018). 6. Computational experiments and results In this section, we evaluate the performance of the SAA algorithm developed for the two-stage stochastic program. The mathematical models and Algorithm 1 were programmed in Python 3.7 on a workstation Intel(R) Xeon(R) Gold 5118 CPU @2.30 GHz (12 cores) with 64 GB RAM, using Gurobi 9.1 when necessary as a MILP solver. 6.1. Stochastic model results To analyze the computational performance of the SAA algorithm and provide managerial insights, we calculated the following metrics: •Optimality gap reached by Algorithm 1. •CPU time (in seconds) incurred by Algorithm 1. •Total cost (𝒞T), fixed cost (𝒞F), expected first-echelon distribution cost (𝒞1E), and expected second-echelon distribution cost (𝒞2E). Fig. 2. Customer segments (i.e., zip codes of Manhattan, NY) employed. J.C. Pina-Pardo et al. EURO Journal on Transportation and Logistics 11 (2022) 100079 7 •Expected total cost of shipping products from DCs to customer segments (𝒞DC→CS). It represents the expected value of the second term of equation (13). •Expected total cost per lost (or outsourced) demand (𝒞L). It represents the expected value of the last term of equation (13). We begin by describing the computational performance of Algorithm 1. Given the computational complexity of the SAA problem solved in Step 3, we report the results reached by Gurobi in a maximum run-time of one hour. Considering a penalty per unit of lost (or outsourced) demand of p k =25 cost units, ∀k ∈K, Table 2 shows the optimality gaps and CPU times obtained for the different sample sizes (N) considered. As expected, optimality gaps decrease as the sample size increases, although at the expense of the computational effort required. Algorithm 1 reached an optimality gap of 0.39% in 100 seconds for N =10 scenarios, and an optimality gap of 0.22% in 210 seconds for N =20. However, while an optimality gap of 0.11% was reached for a sample size of N =30 scenarios, CPU time increased considerably to 723 seconds. Furthermore, after observing the final network design obtained for each sample size considered, all solutions suggest installing the same facilities: DC 1 and UTP 4. This is not surprising since (1) Algorithm 1 selects the network design with the minimum upper bound estimator (see Step 8) and (2) these locations are close to where demand is concentrated (see Fig. 3). Given the satisfactory results achieved for N =20 scenarios, we retained this value for the rest of the computational experiments. Regarding the cost-based metrics, Table 3 shows that the second distribution echelon is the highest cost driver, responsible for about 71% of total costs. This can be explained because each customer segment contains many individual end consumers (more than 1300 in total) that must be visited by the second-echelon vehicles. Additionally, the low value of the 𝒞DC→CS metric shows the importance of using intermediate facilities: the SAA algorithm prefers to open UTP 4 instead of shipping all products directly to customer segments from DCs. This is mainly because, as in real life, making last-mile deliveries using first-echelon vehicles (i.e., large freight trucks) is expensive (see Table 1). Finally, the expected total cost per lost demand (𝒞L)shows that Algorithm 1 decides to satisfy almost all the customer segments’ demand, which is likely due to the high penalty per unit established (p k =25, ∀k ∈K). 6.1.1. Effects of the facility capacity As mentioned before, Algorithm 1 recommends installing DC 1 and UTP 4 regardless of the sample size considered. Notably, given the low amount of lost demand, the installation of a single DC suggests that the total demand of all customers does not exceed its capacity. This led us to carry out a sensitivity analysis to understand the effects of the facilities’ capacity considered. For this, we analyzed the results of the SAA technique over the base-line values (i.e., D i =3000 parcels, ∀i ∈I, and S j =1000 parcels, ∀j ∈J), half of these capacities, and a quarter of these capacities. Considering a penalty cost per unit of lost (or outsourced) demand of p k =25 cost units, ∀k ∈K, Table 4 shows the different network designs obtained and the capacity utilization of each of the installed facilities (i. e., the total quantity of products shipped divided by their capacity). As Fig. 3. Locations to open logistics facilities. Table 1 First and second echelon vehicle parameters. Truck Van Cargo bike Capacity 250.0 50.0 15.0 Fixed cost 200.0 100.0 25.0 Cost per distance 1.0 0.5 0.1 Cost per time 80.0 40.0 30.0 Speed inter stop 10.0 15.0 7.0 Speed line-haul 15.0 20.0 10.0 Service time 0.1 0.1 0.1 Route setup time 0.7 0.5 0.3 Allowed service time – 10.0 10.0 Table 2 Computational performance of Algorithm 1. N Upper bound Lower bound Gap CPU time 10 11022.8 10980.1 0.39% 100.0 20 11040.9 11017.0 0.22% 210.6 30 11030.7 11018.4 0.11% 723.0 J.C. Pina-Pardo et al. EURO Journal on Transportation and Logistics 11 (2022) 100079 8 expected, as the capacity considered decreases, the number of open facilities increases. Furthermore, facilities are installed based on their proximity to high-demand areas (see Fig. 3). When considering a quarter of the base-line capacities, 100% of the operational capacity of UTP 1, UTP 3, UTP 4, and UTP 5 is used. These results state that these locations are attractive and play an essential role in minimizing total costs. 6.1.2. Effects of the penalty cost We now study the effects of the penalty cost (per unit of lost demand) parameter. For this, we consider that p k takes values of 12.5, 25.0, or 50.0 cost units, ∀k ∈K. Fig. 4 shows the percentage of the total cost attributed to the installation (or fixed) cost, to the expected secondechelon distribution cost, and to the expected cost per lost demand, considering different penalties and facility capacities. (For the sake of exposure, we omit the percentages contributed by the 𝒞1E and 𝒞DC→CS metrics, as they remain relatively stable – at 13.1% and 0.4% on average, respectively – when varying penalties and capacities.) Considering the base-line capacities and half of these, the percentage associated with installation costs varied around 0.5% when the penalty increased from 12.5 to 25.0 per unit of lost demand. However, this percentage drastically increased by 8.8% (from 19.7% to 28.5%) when considering a quarter of the facilities’ capacities. These results provide an interesting insight for managers: when the cost per unit of lost (or outsourced) demand is relatively high, it is more convenient to open more facilities, even if they have a small capacity. The latter can also be observed in the percentage of total costs contributed to the expected cost per lost demand: while it reached about 25% when considering the lowest penalty and a quarter of the base-line capacities, it drastically decreased near to 0% when the penalty parameter increases to p k =25 cost units (∀k ∈K). Finally, the contribution of the expected secondechelon distribution cost quickly increased as the penalty per unit of lost demand increased from 12.5 to 25.0 cost units. Nevertheless, it remained stable when increasing the penalty parameter from 25 to 50 cost units (e.g., the percentage of total costs attributed to the second distribution echelon reached about 69% for both p k =25 and p k =50 when considering the half of the base-line capacities). Regarding the location of the facilities, Table 5 presents the different network designs obtained for all combinations of the penalty and capacity parameters used. This table shows that the installation of DC 2, UTP 2, UTP 6, UTP 7, and UTP 8 is always discarded, primarily because they are farther away from the demand concentration (see Fig. 3). When considering p k =12.5 and the base-line capacities, the solution of Algorithm 1 shows that UTP 3 should be installed. However, when the facilities’ size decreases and the penalty increases to 25 or 50 cost units, the results suggest opening UTP 4 and UTP 5. The latter coincides with the results of the simulation-based approach, as we will show in Section 6.2.2. 6.2. Validation of the two-stage stochastic model results This section provides an extensive validation of the two-stage stochastic program results. We begin by describing the value of the stochastic solutions (Maggioni and Wallace, 2012). We then present a simulation-based approach that employs the deterministic MILP model (1)–(11). 6.2.1. Value of the stochastic solutions A first approach to validate the results of the two-stage stochastic program is to solve the Expected Value Problem (EVP), which is the deterministic counterpart of (12) where all random variables (i.e., customer demand) are replaced by their expected values (Maggioni and Table 3 Cost metrics of the solution of Algorithm 1 for N =20. 𝒞T 𝒞F 𝒞1E 𝒞DC→CS 𝒞2E 𝒞L 11040.9 1400.0 1598.8 63.0 7815.6 163.5 Table 4 Facility location (represented by 1) and their capacity utilization, considering the base-line capacities, half and a quarter of them. Facility Facility location Capacity utilization Base-line Half Quarter Base-line Half Quarter DC 1 1 1 1 45.9% 85.9% 100.0% DC 2 0 0 0 0.0% 0.0% 0.0% DC 3 0 0 1 0.0% 0.0% 36.5% UTP 1 0 0 1 0.0% 0.0% 100.0% UTP 2 0 0 0 0.0% 0.0% 0.0% UTP 3 0 0 1 0.0% 0.0% 100.0% UTP 4 1 1 1 100.0% 100.0% 100.0% UTP 5 0 1 1 0.0% 100.0% 100.0% UTP 6 0 0 0 0.0% 0.0% 0.0% UTP 7 0 0 0 0.0% 0.0% 0.0% UTP 8 0 0 0 0.0% 0.0% 0.0% Fig. 4. Percentage of total cost attributed to the total installation cost, the expected total cost per lost demand, and the expected second-echelon distribution cost, considering p k ∈{12.5, 25.0, 50.0} and the base-line capacities, half and a quarter of them. J.C. Pina-Pardo et al.