Uncertainty analysis of autonomous delivery robot operations for last-mile logistics in European cities

Abstract

Although autonomous delivery robots (ADRs) are widely anticipated to significantly enhance the efficiency of last-mile logistics operations in dense urban environments in the coming years, their impact on logistics service providers’ supply chains has yet to be accurately assessed on a large scale. The primary objective of this article is to quantify the uncertainty in the cost of ADR operations as a function of the stochastic behavior of given input variables. First, ADR operations are modeled using the continuous approximation methodology. The mathematical formulations proposed in this article relate certain ADR and service area input parameters to provide an estimate of the carrier’s last-mile operating costs. An uncertainty analysis based on the Monte Carlo approach is then performed. The numerical results indicate that the implementation of a two-echelon delivery scheme with heavy-duty vehicles cooperating with ADRs through an urban logistics micro-hub would, on average, reduce last-mile operating costs by more than 10% in small and medium-sized European cities, given adequate operational conditions.

Keywords:

• autonomous vehicle
• continuous approximation
• Monte Carlo
• two-echelon
• urban consolidation center

1. Introduction

Freight vehicles account for around 20% of urban traffic (Russo & Comi, 2012). Deliveries are expected to become more frequent and fragmented in the coming years, especially due to the skyrocketing use of e-commerce. The on-demand economy generates much more externalities, as more vehicles have to be used to deliver parcels in a very limited window of time. Last-mile logistics operations in city centers lead to very negative effects: increased traffic congestion, safety problems for pedestrians, cyclists, and delivery staff, as well as air and noise pollution.

Autonomous delivery robots (ADRs) are a specific application of autonomous technologies for last-mile logistics operations. They are expected to significantly improve the efficiency of last-mile distribution operations in the coming years. In addition, the COVID-19 pandemic has certainly had a positive impact on the acceptance of ADRs by end users, as they increase the social distance between individuals. However, the real impact of the pandemic on the use and acceptance of ADRs remains to be accurately quantified (Pani et al., 2020).

From the perspective of logistics service providers (LSPs), the potential of ADRs in terms of reducing externalities and operational costs is widely believed to be high (Figliozzi, 2020; Jennings & Figliozzi, 2019; Milakis et al., 2017). Given this attractiveness, many prototypes of ADRs are either under development or have already been piloted in very specific use cases. Baum et al. (2019) identified 39 automated micro-vehicles for the last mile, which are either laboratory prototypes or already operating in cities. The market is currently in full expansion.

Nevertheless, the magnitude of cost and externality reductions is highly dependent on some specific operational variables of the robots and the service area. In some cases, business-as-usual (BAU) operations with conventional vans or other delivery strategies (Rechavi & Toch, 2020) may be more competitive. In other words, of the 39 ADR-based initiatives identified by Baum et al. (2019), many will result in higher operating costs, depending on the delivery schemes adopted, hindering their adoption by LSPs. The market is not yet mature, and many directions of development are still being explored. A more in-depth assessment of the uncertainty of ADR operations and the optimal use cases for ADR is clearly needed.

Currently, a large part of the literature on last-mile logistics focuses on the development of formulations and algorithms to numerically optimize operations under some constraints (see Table 1; Balcik et al., 2008; Boysen et al., 2018; Simoni et al., 2020; Ulmer & Streng, 2019). These methods are fundamental from an operational perspective. However, our goal in this article is to provide strategic insights to identify the most promising operational schemes that will emerge thanks to ADRs, and discrete numerical methods are not as efficient for this type of task. To evaluate ADR operations in several cities or metropolitan regions with different urban characteristics, many scenarios have to be studied, i.e., many simulations have to be performed, which may be too time-consuming due to the high computational cost.

Table 1. Comparative summary of related literature on ADR operation optimization problem.

Reference	Operation type	Problem setup	Modeling technique
Daganzo (1984)	Delivery vehicles from DC	TSP	CA
del Castillo (1999)	Delivery vehicles from DC	TSP	CA
Estrada and Roca-Riu (2017)	Urban consolidation center	Two-echelon VRP	CA
Figliozzi (2020)	Truck-launched ADRs	TSP with Robot	CA
Jennings and Figliozzi (2019)	Truck-launched ADRs	TSP with Robot	CA
Balcik et al. (2008)	Delivery trucks from DC	Inventory and VRP	Mixed-integer programming
Boysen et al. (2018)	Truck-launched ADRs with micro depots	TSP with Robot	Mixed-integer programming
Simoni et al. (2020)	Truck-launched ADRs	TSP with Robot	Integer programming
Ulmer and Streng (2019)	Pickup stations and ADRs	Dynamic dispatching problem for same-day delivery	Markov decision process

TSP: traveling salesman problem; VRP: vehicle routing problem.

Consequently, in order to position the paper at a strategic level, the proposed mathematical model will be based on the continuous approximation (CA) method (Daganzo, 1984). The main advantage of the CA method is that it provides compact and transparent equations in which the main decision variables and tradeoffs can be identified (Daganzo et al., 2012). In addition, a large number of scenarios and hypotheses can be evaluated quickly and easily because these models do not rely on complex numerical methods.

As we have already explained, a large number of ADR designs and distribution strategies are emerging. Boysen et al. (2018), Figliozzi (2020), and Simoni et al. (2020) analyze the performance of truck-launched ADRs. In this particular delivery scheme, ADRs work in conjunction with “mothership” trucks that can pick them up and drop them off at a few specific locations in the city, with the aim of reducing the total distance traveled by the entire fleet. Unfortunately, this scheme is not necessarily the most optimal, as some regulatory challenges and traffic issues may complicate the operation of these “mothership” trucks in dense and crowded areas (Jennings & Figliozzi, 2019).

To address these challenges that truck-based ADRs may face, we propose to investigate another type of ADR last-mile operations. In this article, it is assumed that a logistics micro-hub has been implemented in the considered service region. This logistics micro-hub is conceived by the authors as a dedicated infrastructure where the items to be delivered can be transferred from heavy-duty vehicles (HDVs) to ADRs. This logistics micro-hub is a “safe space” where the carrier can carry out the operations without disturbing other road users, avoiding, for example, the problem of insufficient parking (Dablanc, 2019; EIT Urban Mobility, 2023; Estrada & Roca-Riu, 2017; Navarro et al., 2016). This approach is very similar to the work of Ulmer and Streng (2019), but they use a discrete numerical method to solve the optimization problem. As explained above, this article is situated at a more strategic level.

The analytical formulation proposed in this article relates some input parameters of ADRs, service regions, and demand density and provides an estimate of the last-mile operational costs. A Monte Carlo approach is used to account for the uncertainty of the input parameters in our model. To the best of our knowledge, this is completely novel in the field of ADR operations assessment. Instead of working with deterministic values, some probabilistic distributions are defined for the different input parameters. Due to the novelty of ADRs and the small number of real implementations in urban environments with different traffic flows, there is still uncertainty in some cost and performance metrics. The total cost of last-mile operations therefore follows its own probabilistic distribution, obtained through the Monte Carlo process. Much more insight can be gained with this methodology as the model output is described by a mean and a standard deviation.

To summarize and conclude this section, the main contributions of this article to the literature are as follows:

• The evaluation of a last-mile delivery operation scheme based on the use of a heterogeneous fleet of vans and ADRs, cooperating in different parts of the city according to their intrinsic characteristics. A comprehensive cost structure is also proposed.

• The optimization of delivery operations is based on two key decision variables, namely the micro-hub service area radius and the distance between the micro-hub and the city center.

• The application of our logistics optimization model to a set of more than 150 small and medium-sized European cities, aggregating accurate data from different sources. This enables us to provide representative logistics use cases for ADRs across Europe.

• Identifying the key operational parameters that would enhance the economically efficient use of ADR-based delivery operations compared to conventional van deliveries.

2. Methodological framework

In this article, we will compare ADR last-mile deliveries with BAU distribution patterns. The service region is modeled as a disk of radius $r_c$ [km]. We assume that it is representative of a core city, as defined by Dijkstra et al. (2019). This definition of a city is independent from local administrative delimitations, i.e., we will be able to compare cities from different European countries. We assume that the carrier daily demand density ${\delta }_c\left(r\right)$ [receivers/carrier/km²/day] probability distribution function (PDF) is continuous and decreases with distance to the city center $r$ [km] because the population density is higher in the city center. We assume that all variables of the problem are independent from the angular coordinate, and the ring-radial metric is used. Nevertheless, we used the L2 metric in the drawings presented in the article for the sake of clarity. The authors are aware that the shape of the routes presented in the different drawings differs from the routes obtained with the ring-radial metric.

In the BAU situation, parcels are directly taken from a distribution center (DC) to final receivers with light commercial vehicles (LCVs) (see Figure 1). The expected distance between the DC and the center of the service region is $l_{\mathit{\text{DC}}}$ [km]. In this configuration, the carrier’s main objective is to minimize the distance traveled by its fleet. This is a particular instance of the so-called vehicle routing problem (VRP).

Figure 1. Business-as-usual delivery scheme.

In the case of using ADRs for last-mile operations, we assume that the carrier’s distribution network is made up of two distinct echelons (see Figure 2a,b).

Figure 2. Two-echelon delivery scheme.

First of all, a logistics micro-hub is located in the city, at a distance $r_h$ [km] from the city center. For all receivers that are located within a radius $r_0$ [km] from the city center (Zone A, see Figure 2a), parcels are first carried with HDVs from the carrier’s DC to the micro-hub. Then, parcels are transhipped and delivered to final customers using ADRs. Implementing this intermediate step with the logistics micro-hub is almost mandatory because ADRs’ range of action is limited and they cannot directly come from the carrier’s DC. ADRs are much more agile and flexible than HDVs, which is a great advantage at a local scale, especially in dense and complex urban environments, but their battery capacity is limited.

If the final receiver is located beyond $r_0$ (Zone B, see Figure 2b), the parcel is delivered with LCVs, directly from the DC without going through the logistics micro-hub. Nevertheless, if necessary, LCVs can travel through the city center to access the receivers more rapidly.

The outputs of this modeling work are the expressions of some performance metrics such as run distance, travel times, and finally total operation costs in both the BAU and two-echelon distribution strategies (see next section).

The analytical formulations of the carrier’s total operation costs depend on input parameters (either from the LCVs, ADRs, HDVs, or considered cities) whose values are highly uncertain. As a consequence, a completely deterministic approach is not enough to quantify the reduction in operation costs induced by the use of ADRs. We will use the Monte Carlo methodology to quantify the model uncertainty.

The first stage of the Monte Carlo approach is to assume some given PDFs for the model input parameters. Then the last-mile operation cost PDF is obtained by computing the outputs of the analytical model for a very large number of input parameter random samplings. One input parameter random sampling returns one output value. By repeating this process many times, the model output PDF can be obtained.

3. Operational modeling

The main objective of this section is to analytically model the total operational costs in both BAU and two-echelon distribution patterns. The estimation of the different tour lengths is largely based on the heuristics developed by del Castillo (1999).

Table 2 presents the notations used in the modeling of logistics operations.

Table 2. Logistics key performance indicator notations.

Notation	Description
$L\left(N,R\right)$	Optimized length of a traveling salesman tour visiting $N$ customers in a disk of radius $R,$ considering a ring-radial metrics and demand decreasing with distance to city center (del Castillo, 1999)
$S(N)$	Number of subsectors of a traveling salesman tour visiting $N$ customers in a disk, considering a ring-radial metrics and demand decreasing with distance to city center (del Castillo, 1999)
$\rho (N)$	Optimum partition ratio of a traveling salesman tour visiting $N$ customers in a disk, considering a ring-radial metrics and demand decreasing with distance to city center (del Castillo, 1999)
$N_i^{\mathit{\text{veh}}}$	Total daily number of recipients visited by vehicle fleet $\mathit{\text{veh}}$ in distribution scheme $i$
$N_{t,c}$	Total daily number of visited recipients visited by carrier $c$ in a given city
${\overline{r}}_{\mathit{\text{BAU}}}$	Expected distance between the city center and final receivers in the BAU distribution scheme
$r_0$	Radius of Zone A in the two-echelon delivery scheme
${\rho }_h$	ADR access distance between the logistics micro-hub and final recipients in Zone A in the two-echelon distribution scheme
$r_h$	Distance between the city center and the urban micro-hub in the two-echelon delivery scheme
${\overline{r}}_z$	Expected distance between the city center and final receivers in Zone $z$ of the two-echelon delivery scheme
$d_{\mathit{\text{BAU}}}^{\mathit{\text{LH}}}$	Expected line-haul distance between the carrier’s DC and final receivers in the BAU distribution scheme
$d_B^{\mathit{\text{LH}}}$	Expected line-haul distance between the carrier’s DC and final receivers located in Zone B of the two-echelon distribution scheme
$l_d^{\mathit{\text{ADR}}}$	Expected distance between two consecutive recipients visited by an ADR in the two-echelon distribution scheme
$t_d^{\mathit{\text{ADR}}}$	Expected ADR total delivery time in the two-echelon distribution scheme, including travel time between two consecutive recipients and parcel hand-over process to recipient
${\mathrm{\Psi }}_{t,i}^{\mathit{\text{veh}}}$	Expected capacity, i.e., total number of recipients visited along a delivery route of vehicle $\mathit{\text{veh}},$ considering time restriction, in distribution scheme $i$
${\mathrm{\Psi }}_{b,i}^{\mathit{\text{veh}}}$	Expected capacity, i.e., total number of recipients visited along a delivery route, considering battery restriction (if applicable) of vehicle $\mathit{\text{veh}}$ in distribution scheme $i$
${\mathrm{\Psi }}_i^{\mathit{\text{veh}}}$	Expected capacity, i.e., total number of recipients visited along a delivery route of vehicle $\mathit{\text{veh}}$ in distribution scheme $i$
$n^{\mathit{\text{ADR}}}$	Expected size of ADR fleet in the two-echelon distribution scheme
$D_i^{\mathit{\text{veh}}}$	Expected total distance traveled daily by vehicle fleet $\mathit{\text{veh}}$ in distribution scheme $i$
$D_B^{\mathit{\text{LCV}}}$	Expected distance traveled daily by the LCV fleet in Zone B of the two-echelon distribution scheme
$T_i^{\mathit{\text{veh}}}$	Expected total daily working time of vehicle fleet $\mathit{\text{veh}}$ in distribution scheme $i$
$E_i^{\mathit{\text{veh}}}$	Expected total daily energy consumption of vehicle fleet $\mathit{\text{veh}}$ in distribution scheme $i$
$Z_i^{\mathit{\text{veh}}}$	Expected total daily operation costs of vehicle fleet $\mathit{\text{veh}}$ in distribution scheme $i$
${\mathrm{\Omega }}_h$	Expected total daily costs of micro-hub in the two-echelon distribution scheme, including infrastructure and personnel costs
$Z_i$	Expected total daily operation costs of delivery scheme $i$

3.1. Tour length estimation

As previously mentioned, we consider a demand density decreasing with distance to the city center and a ring-radial metric. The objective of this subsection is to provide an estimation of the optimized length $L\left(N,R\right)$ [km] of a tour visiting $N$ customers in a disk of radius $R.$

The first step of the tour length estimation is to compute the number of subsectors $S(N).$ (1) $$S\left(N\right)=\mathit{\text{min}}\left({\left[\frac{N}{25}\right]}^{+};{\left[\frac{N}{60}+2\right]}^{+};24\right)$$(1) where $N$ is the total number of receivers that are to be visited along the tour.

The expression of $S\left(N\right)$ was obtained after “exhaustive computer tests” (del Castillo, 1999), and the numerical coefficients presented in Equation (1)(1) $$S\left(N\right)=\mathit{\text{min}}\left({\left[\frac{N}{25}\right]}^{+};{\left[\frac{N}{60}+2\right]}^{+};24\right)$$(1) are in any case related to the capacity of the vehicle operating the route. The formulation presented here generates an approximate solution to the traveling salesman problem, not the VRP (del Castillo, 1999).

The optimum partition ratio $\rho$ is then given by (2) $$\rho \left(N\right)=\left(\frac{2}{3}C-\frac{\theta }{4}-2\right){\left(\frac{1}{6}C\theta +\frac{2}{3}C-\frac{1}{2}\theta \right)}^{-1}$$(2) where $C=\mathit{\text{max}}\left\{\frac{N}{S};7\right\}$ and $\theta =\frac{2\pi }{S}$

Finally, the total tour length $L\left(N,R\right)$ visiting $N$ points in a disk of radius $R$ is estimated as (3) $$L\left(N,R\right)=\hspace{0.17em}\mathit{\text{RS}}[2\rho +\frac{1}{12}{\rho }^{2}C\theta +\frac{1}{2}\left(1+\rho \right)\theta \left(1-\frac{\rho }{2}\right)+\frac{1}{3}C{\left(1-\rho \right)}^2]$$(3)

3.2. Business-as-usual delivery strategy

In this subsection, the carrier’s total operation costs in the BAU delivery strategy are estimated.

First, we compute the total daily number of points $N_{t,c}$ [receivers/day] that are to be visited by the LCV fleet of carrier $c.$ (4) $$N_{t,c}=2\pi \underset{0}{\overset{r_c}{\int }}r{\delta }_c\left(r\right)\mathit{\text{dr}}$$(4) where ${\delta }_c\left(r\right)$ [receivers/km²/day] is the demand density of carrier $c$ and $r_c$ [km] the city radius.

The conditional PDF $f_{\mathit{\text{BAU}}}(r)$ is defined as (5) $$f_{\mathit{\text{BAU}}}\left(r\right)=\frac{{\delta }_c\left(r\right)}{{\int }_0^{r_c}{\delta }_c\left(r\right)\mathit{\text{dr}}}$$(5)

The expected distance ${\overline{r}}_{\mathit{\text{BAU}}}$ [km] between the city center and the final receivers is then computed. (6) $${\overline{r}}_{\mathit{\text{BAU}}}=\underset{0}{\overset{r_c}{\int }}r{f}_{\mathit{\text{BAU}}}\left(r\right)\mathit{\text{dr}}$$(6)

The expected line-haul distance $d_{\mathit{\text{BAU}}}^{\mathit{\text{LH}}}$ [km] between the carrier’s DC and the final receivers is given by Equation (7)(7) $$d_{\mathit{\text{BAU}}}^{\mathit{\text{LH}}}=\mathit{\text{max}}\left({\overline{r}}_{\mathit{\text{BAU}}},l_{\mathit{\text{DC}}}\right)+\frac{\pi -2}{\pi }\mathit{\text{min}}\left({\overline{r}}_{\mathit{\text{BAU}}},l_{\mathit{\text{DC}}}\right)$$(7) . (7) $$d_{\mathit{\text{BAU}}}^{\mathit{\text{LH}}}=\mathit{\text{max}}\left({\overline{r}}_{\mathit{\text{BAU}}},l_{\mathit{\text{DC}}}\right)+\frac{\pi -2}{\pi }\mathit{\text{min}}\left({\overline{r}}_{\mathit{\text{BAU}}},l_{\mathit{\text{DC}}}\right)$$(7)

The reader can refer to the Supplemental Information for more details about the modeling of $d_{\mathit{\text{BAU}}}^{\mathit{\text{LH}}}.$

The next step is to compute the LCV capacity ${\mathrm{\Psi }}_{\mathit{\text{BAU}}}^{\mathit{\text{LCV}}},$ which corresponds to the total number of receivers visited per LCV route. ${\mathrm{\Psi }}_{\mathit{\text{BAU}}}^{\mathit{\text{LCV}}}$ is subject to two restrictions. First, LCVs have a limited volume capacity, i.e., a given maximum number of parcels can be loaded in the vehicle at the beginning of the delivery route. Secondly, we assume that all delivery operations must be done within a given time window $H$ [h], i.e., the LCV route duration cannot exceed $H.$

As a consequence, ${\mathrm{\Psi }}_{\mathit{\text{BAU}}}^{\mathit{\text{LCV}}}$ is estimated as (8) $${\mathrm{\Psi }}_{\mathit{\text{BAU}}}^{\mathit{\text{LCV}}}=\mathit{\text{min}}\left\{N_{t,c};C^{\mathit{\text{LCV}}};\frac{H-\frac{d_{\mathit{\text{BAU}}}^{\mathit{\text{LH}}}}{v_{\mathit{\text{LH}}}^{\mathit{\text{LCV}}}}}{\frac{1}{v_L^{\mathit{\text{LCV}}}}\frac{L\left(N_{t,c};r_c\right)}{N_{t,c}}+{\tau }_d^{\mathit{\text{LCV}}}}\right\}$$(8) where $C^{\mathit{\text{LCV}}}$ [parcels] is the LCV maximum volume capacity, $v_{\mathit{\text{LH}}}^{\mathit{\text{LCV}}}$ [km/h] the LCV expected line-haul commercial speed, $v_L^{\mathit{\text{LCV}}}$ [km/h] the LCV expected commercial speed in the local urban grid, and ${\tau }_d^{\mathit{\text{LCV}}}$ [h] the LCV expected stop time per delivery, including parking and parcel hand-over to the final recipient.

The daily total distance $D_{\mathit{\text{BAU}}}^{\mathit{\text{LCV}}}$ [veh-km/day] traveled by the LCV fleet in this BAU scenario is computed. (9) $$D_{\mathit{\text{BAU}}}^{\mathit{\text{LCV}}}=2d_{\mathit{\text{BAU}}}^{\mathit{\text{LH}}}\frac{N_{t,c}}{{\mathrm{\Psi }}_{\mathit{\text{BAU}}}^{\mathit{\text{LCV}}}}+L\left(N_{t,c},r_c\right)$$(9) where $L\left(N_{t,c},r_c\right)$ is the traveling salesman total tour length to visit $N_{t,c}$ receivers in a disk of radius $r_c,$ as modeled in Subsection 3.1.

The daily total LCV working time $T_{\mathit{\text{BAU}}}^{\mathit{\text{LCV}}}$ [veh-h/day] is estimated as (10) $$T_{\mathit{\text{BAU}}}^{\mathit{\text{LCV}}}=\frac{1}{v_{\mathit{\text{LH}}}^{\mathit{\text{LCV}}}}2d_{\mathit{\text{BAU}}}^{\mathit{\text{LH}}}\frac{N_{t,c}}{{\mathrm{\Psi }}_{\mathit{\text{BAU}}}^{\mathit{\text{LCV}}}}+\frac{1}{v_L^{\mathit{\text{LCV}}}}L\left(N_{t,c},r_c\right)+{\tau }_d^{\mathit{\text{LCV}}}{N}_{t,c}$$(10)

Finally, the daily total operation costs $Z_{\mathit{\text{BAU}}}$ [€/day] in the BAU delivery scenario are (11) $$Z_{\mathit{\text{BAU}}}=c_d^{\mathit{\text{LCV}}}{D}_{\mathit{\text{BAU}}}^{\mathit{\text{LCV}}}+c_t^{\mathit{\text{LCV}}}{T}_{\mathit{\text{BAU}}}^{\mathit{\text{LCV}}}$$(11) where $c_d^{\mathit{\text{LCV}}}$ [€/veh-km] is the LCV expected unit distance cost and $c_t^{\mathit{\text{LCV}}}$ [€/veh-h] the LCV expected unit time cost.

3.3. Two-echelon delivery strategy

In this subsection, we will model the ADR and HDV operation costs for the receivers that are located at less than $r_0$ [km] from the city center (Zone A, see Figure 2a) as well as the operation costs induced by the LCV fleet in the city’s external ring (Zone B, see Figure 2b).

First, we will consider the recipients located in Zone A (see Figure 2a). The total daily number of receivers $N_{\mathit{\text{ech}}}^{\mathit{\text{ADR}}}$ [receivers/day] visited by the ADR fleet of carrier $c$ is given by (12) $$N_{\mathit{\text{ech}}}^{\mathit{\text{ADR}}}=2\pi \underset{0}{\overset{r_0}{\int }}r{\delta }_c\left(r\right)\mathit{\text{dr}}$$(12)

The conditional PDF $f_A\left(r\right)$ in Zone A is then defined. (13) $$f_A\left(r\right)=\frac{{\delta }_c\left(r\right)}{{\int }_0^{r_0}{\delta }_c\left(r\right)\mathit{\text{dr}}}$$(13)

The expected distance ${\overline{r}}_A$ [km] between the city center and the final receivers in Zone A is estimated as (14) $${\overline{r}}_A\left(r_0\right)=\underset{0}{\overset{r_0}{\int }}r{f}_A\left(r\right)\mathit{\text{dr}}$$(14)

The ADR expected access distance ${\rho }_h$ [km] between the logistics micro-hub and the final recipients in Zone A is given by Equation (15)(15) $${\rho }_h=\mathit{\text{max}}\left\{{\overline{r}}_A;r_h\right\}+\frac{\pi -2}{\pi }\mathit{\text{min}}\left\{{\overline{r}}_A;r_h\right\}$$(15) (see Supplemental Information). (15) $${\rho }_h=\mathit{\text{max}}\left\{{\overline{r}}_A;r_h\right\}+\frac{\pi -2}{\pi }\mathit{\text{min}}\left\{{\overline{r}}_A;r_h\right\}$$(15) where $r_h$ [km] is the distance between the city center and the logistics micro-hub (see Figure 2a).

We also need to compute the line-haul distance between the DC and the micro-hub $l_h$ [km] (see Figure 2a). (16) $$l_h=\mathit{\text{min}}\left\{l_{\mathit{\text{DC}}}+r_h;l_{\mathit{\text{DC}}}-r_h+{\theta }_h{r}_h\right\}$$(16) where ${\theta }_h$ [rad] is the angular coordinate of the micro-hub (see Figure 2a).

The expected distance between two consecutive recipients visited by an ADR in the two-echelon distribution scheme $l_d^{\mathit{\text{ADR}}}$ [km] and the expected ADR total delivery time in the two-echelon distribution scheme, including travel time between two consecutive recipients and parcel hand-over process to the recipient, $t_d^{\mathit{\text{ADR}}}$ [h] are estimated as (17a) $$l_d^{\mathit{\text{ADR}}}=\frac{L\left(N_{\mathit{\text{ech}}}^{\mathit{\text{ADR}}};r_0\right)}{N_{\mathit{\text{ech}}}^{\mathit{\text{ADR}}}}$$(17a) (17b) $$t_d^{\mathit{\text{ADR}}}=\frac{1}{v^{\mathit{\text{ADR}}}}\frac{L\left(N_{\mathit{\text{ech}}}^{\mathit{\text{ADR}}};r_0\right)}{N_{\mathit{\text{ech}}}^{\mathit{\text{ADR}}}}+{\tau }_d^{\mathit{\text{ADR}}}$$(17b) where ${\tau }_d^{\mathit{\text{ADR}}}$ [h] is the ADR expected stop time per parcel delivery, including waiting time and parcel hand-over to the final recipient.

As we previously did in the modeling of the BAU delivery scheme, we need to compute the expected number of receivers visited per ADR route ${\mathrm{\Psi }}_{\mathit{\text{ech}}}^{\mathit{\text{ADR}}}.$ As in the case of LCVs, ${\mathrm{\Psi }}_{\mathit{\text{ech}}}^{\mathit{\text{ADR}}}$ is subject to the vehicle’s limited volume capacity and time window restrictions. In addition, ADR operations must comply with another restriction linked to the ADR’s limited battery capacity. The total energy needed to perform an entire ADR delivery route cannot exceed the total ADR battery capacity. This is an important difference between vehicles with internal combustion engines and battery-electric ones. In Equation (18a)(18a) $${\mathrm{\Psi }}_{b,\mathit{\text{ech}}}^{\mathit{\text{ADR}}}=\frac{{\mathit{\text{BC}}}^{\mathit{\text{ADR}}}-2{\rho }_h\left[{\beta }_{\mathit{\text{mech}}}^{\mathit{\text{ADR}}}+\frac{P_e^{\mathit{\text{ADR}}}}{v^{\mathit{\text{ADR}}}}\right]}{l_d^{\mathit{\text{ADR}}}{\beta }_{\mathit{\text{mech}}}^{\mathit{\text{ADR}}}+t_d^{\mathit{\text{ADR}}}{P}_e^{\mathit{\text{ADR}}}}$$(18a) , ${\mathrm{\Psi }}_{b,\mathit{\text{ech}}}^{\mathit{\text{ADR}}}$ [parcels] estimates the expected number of recipients visited per ADR route if only the ADR battery capacity restriction was considered. Similarly, in Equation (18b)(18b) $${\mathrm{\Psi }}_{t,\mathit{\text{ech}}}^{\mathit{\text{ADR}}}=\frac{H-2\frac{l_h}{v_{\mathit{\text{LH}}}^{\mathit{\text{HDV}}}}-\frac{1}{2}{\mathrm{\Psi }}_{\mathit{\text{ech}}}^{\mathit{\text{HDV}}}{\tau }_{\mathit{\text{LU}}}^{\mathit{\text{HDV}}}-2\frac{{\rho }_h}{v^{\mathit{\text{ADR}}}}}{t_d^{\mathit{\text{ADR}}}}$$(18b) , ${\mathrm{\Psi }}_{t,\mathit{\text{ech}}}^{\mathit{\text{ADR}}}$ [parcels] represents the expected number of recipients visited per ADR route if only the delivery time window restriction was considered. (18a) $${\mathrm{\Psi }}_{b,\mathit{\text{ech}}}^{\mathit{\text{ADR}}}=\frac{{\mathit{\text{BC}}}^{\mathit{\text{ADR}}}-2{\rho }_h\left[{\beta }_{\mathit{\text{mech}}}^{\mathit{\text{ADR}}}+\frac{P_e^{\mathit{\text{ADR}}}}{v^{\mathit{\text{ADR}}}}\right]}{l_d^{\mathit{\text{ADR}}}{\beta }_{\mathit{\text{mech}}}^{\mathit{\text{ADR}}}+t_d^{\mathit{\text{ADR}}}{P}_e^{\mathit{\text{ADR}}}}$$(18a) (18b) $${\mathrm{\Psi }}_{t,\mathit{\text{ech}}}^{\mathit{\text{ADR}}}=\frac{H-2\frac{l_h}{v_{\mathit{\text{LH}}}^{\mathit{\text{HDV}}}}-\frac{1}{2}{\mathrm{\Psi }}_{\mathit{\text{ech}}}^{\mathit{\text{HDV}}}{\tau }_{\mathit{\text{LU}}}^{\mathit{\text{HDV}}}-2\frac{{\rho }_h}{v^{\mathit{\text{ADR}}}}}{t_d^{\mathit{\text{ADR}}}}$$(18b) (18c) $${\mathrm{\Psi }}_{\mathit{\text{ech}}}^{\mathit{\text{ADR}}}=\mathit{\text{min}}\left\{N_{\mathit{\text{ech}}}^{\mathit{\text{ADR}}};C^{\mathit{\text{AD}}\mathrm{R}};{\mathrm{\Psi }}_b^{\mathit{\text{ADR}}};{\mathrm{\Psi }}_t^{\mathit{\text{ADR}}}\right\}$$(18c) where $C^{\mathit{\text{ADR}}}$ [parcels] is the ADR maximum volume capacity, $v_{\mathit{\text{LH}}}^{\mathit{\text{HDV}}}$ [km/h] the HDV expected line-haul commercial speed, ${\mathrm{\Psi }}_{\mathit{\text{ech}}}^{\mathit{\text{HDV}}}$ [parcels] the HDV expected capacity, ${\tau }_{\mathit{\text{LU}}}^{\mathit{\text{HDV}}}$ [h] the expected time to load and unload a parcel into the HDV, $v^{\mathit{\text{ADR}}}$ [km/h] the ADR expected commercial speed, $L\left(N_{\mathit{\text{ech}}}^{\mathit{\text{ADR}}},r_0\right)$ the total tour length to visit $N_{\mathit{\text{ech}}}^{\mathit{\text{ADR}}}$ receivers in a disk of radius $r_0,$ ${\mathit{\text{BC}}}^{\mathit{\text{ADR}}}$ [kWh] the ADR battery capacity, ${\beta }_{\mathit{\text{mech}}}^{\mathit{\text{ADR}}}$ [kWh/km] the ADR unit distance energy consumption rate for mechanical propulsion, and $P_e^{\mathit{\text{ADR}}}$ [kW] the ADR electronics and sensor power.

As a consequence, the total daily distance traveled by the ADR fleet $D_{\mathit{\text{ech}}}^{\mathit{\text{ADR}}}$ [veh-km/day] is (19) $$D_{\mathit{\text{ech}}}^{\mathit{\text{ADR}}}\left(r_0\right)=2{\rho }_h\frac{N_{\mathit{\text{ech}}}^{\mathit{\text{ADR}}}}{{\mathrm{\Psi }}_{\mathrm{e}\mathit{\text{ch}}}^{\mathit{\text{ADR}}}}+L\left(N_{\mathit{\text{ech}}}^{\mathit{\text{ADR}}},r_0\right)$$(19)

The total daily ADR working time $T_{\mathit{\text{ech}}}^{\mathit{\text{ADR}}}$ [veh-h/day] and total daily ADR fleet energy consumption $E_{\mathit{\text{ech}}}^{\mathit{\text{ADR}}}$ [kWh/day] in this two-echelon delivery scheme are (20) $$T_{\mathit{\text{ech}}}^{\mathit{\text{ADR}}}=\frac{D_{\mathit{\text{ech}}}^{\mathit{\text{ADR}}}}{v^{\mathit{\text{ADR}}}}+N_{\mathit{\text{ech}}}^{\mathit{\text{ADR}}}{\tau }_d^{\mathit{\text{ADR}}}$$(20) (21) $$E_{\mathit{\text{ech}}}^{\mathit{\text{ADR}}}=D_{\mathit{\text{ech}}}^{\mathit{\text{ADR}}}{\beta }_{\mathit{\text{mech}}}^{\mathit{\text{ADR}}}+T_{\mathit{\text{ech}}}^{\mathit{\text{ADR}}}{P}_e^{\mathit{\text{ADR}}}$$(21)

We then define the total number of ADRs $n^{\mathit{\text{ADR}}}$ [ADRs]. (22) $$n^{\mathit{\text{ADR}}}=\mathit{\text{max}}\left\{\frac{T_{\mathit{\text{ech}}}^{\mathit{\text{ADR}}}}{H-2\frac{l_h}{v_{\mathit{\text{LH}}}^{\mathit{\text{HDV}}}}-\frac{1}{2}{\mathrm{\Psi }}_{\mathit{\text{ech}}}^{\mathit{\text{HDV}}}{\tau }_{\mathit{\text{LU}}}^{\mathit{\text{HDV}}}};\frac{E_{\mathit{\text{ech}}}^{\mathit{\text{ADR}}}}{{\mathit{\text{BC}}}^{\mathit{\text{ADR}}}}\right\}$$(22) where $H-2\frac{l_h}{v_{\mathit{\text{LH}}}^{\mathit{\text{HDV}}}}-\frac{1}{2}{\mathrm{\Psi }}_{\mathit{\text{ech}}}^{\mathit{\text{HDV}}}{\tau }_{\mathit{\text{LU}}}^{\mathit{\text{HDV}}}$ corresponds to the ADR operation time window to deliver all parcels in Zone A.

We assume that ADRs’ batteries are only recharged at night at the logistics micro-hub and that they cannot be swapped.

To ensure that our modeling process is representative, Inequality (23) must be verified (Robusté et al., 1990). (23) $${\mathrm{\Psi }}_{\mathit{\text{ech}}}^{\mathit{\text{ADR}}}\ge 1.5\frac{N_{\mathit{\text{ech}}}^{\mathit{\text{ADR}}}}{{\mathrm{\Psi }}_{\mathit{\text{ech}}}^{\mathit{\text{ADR}}}}$$(23)

As for HDVs, their expected capacity ${\mathrm{\Psi }}_{\mathit{\text{ech}}}^{\mathit{\text{HDV}}}$ [parcels] is estimated as (24) $${\mathrm{\Psi }}_{\mathit{\text{ech}}}^{\mathit{\text{HDV}}}=\mathit{\text{min}}\left\{N_{\mathit{\text{ech}}}^{\mathit{\text{ADR}}};C^{\mathit{\text{HDV}}};\frac{H-2\frac{l_h}{v_{\mathit{\text{LH}}}^{\mathit{\text{HDV}}}}}{{\tau }_{\mathit{\text{LU}}}^{\mathit{\text{HDV}}}}\right\}$$(24) where $C^{\mathit{\text{HDV}}}$ [parcels] is the HDV’s expected maximum volume capacity.

Finally, the total daily distance traveled by the HDV fleet $D_{\mathit{\text{ech}}}^{\mathit{\text{HDV}}}$ [veh-km/day] and total daily working time $T_{\mathit{\text{ech}}}^{\mathit{\text{HDV}}}$ [veh-h/day] are estimated as (25) $$D_{\mathit{\text{ech}}}^{\mathit{\text{HDV}}}\left(r_0\right)=2l_h\frac{N_{\mathit{\text{ech}}}^{\mathit{\text{ADR}}}}{{\mathrm{\Psi }}_{\mathit{\text{ech}}}^{\mathit{\text{HDV}}}}$$(25) (26) $$T_{\mathit{\text{ech}}}^{\mathit{\text{HDV}}}=\frac{1}{v_{\mathit{\text{LH}}}^{\mathit{\text{HDV}}}}{D}_{\mathit{\text{ech}}}^{\mathit{\text{HDV}}}+N_{\mathit{\text{ech}}}^{\mathit{\text{ADR}}}{\tau }_{\mathit{\text{LU}}}^{\mathit{\text{HDV}}}$$(26)

Let us now go through the LCV operations in the city external ring, i.e., the receivers located beyond $r_0$ [km] from the city center (Zone B, see Figure 2b).

The total daily number of receivers $N_{\mathit{\text{ech}}}^{\mathit{\text{LCV}}}$ [receivers/day] that are visited daily by the LCV fleet of carrier $c$ in this second echelon is given by Equation (27)(27) $$N_{\mathit{\text{ech}}}^{\mathit{\text{LCV}}}\left(r_0\right)=2\pi \underset{r_0}{\overset{r_c}{\int }}r{\delta }_c\left(r\right)\mathit{\text{dr}}$$(27) . (27) $$N_{\mathit{\text{ech}}}^{\mathit{\text{LCV}}}\left(r_0\right)=2\pi \underset{r_0}{\overset{r_c}{\int }}r{\delta }_c\left(r\right)\mathit{\text{dr}}$$(27)

The conditional PDF in Zone B $f_B\left(r\right)$ is defined as (28) $$f_B\left(r\right)=\frac{{\delta }_c\left(r\right)}{{\int }_{r_0}^{r_c}{\delta }_c\left(r\right)\mathit{\text{dr}}}$$(28)

The expected distance ${\overline{r}}_B$ [km] between the city center and the receivers located in Zone B is computed as (29) $${\overline{r}}_B\left(r_0\right)=\underset{r_0}{\overset{r_c}{\int }}r{f}_B\left(r\right)\mathit{\text{dr}}$$(29)

This expression of ${\overline{r}}_B$ enables us to compute the expected line-haul distance $d_B^{\mathit{\text{LH}}}$ [km] between the carrier’s DC and the receivers located in Zone B (see Supplemental Information). (30) $$d_B^{\mathit{\text{LH}}}=\mathit{\text{max}}\left(l_{\mathit{\text{DC}}},{\overline{r}}_B\right)+\frac{\pi -2}{\pi }\mathit{\text{min}}\left(l_{\mathit{\text{DC}}},{\overline{r}}_B\right)$$(30)

The total distance $D_B^{\mathit{\text{LCV}}}$ [veh-km/day] traveled daily by the LCV fleet in Zone B is given by Equation (31)(31) $$D_{\mathrm{B}}^{\mathit{\text{LCV}}}\left(r_0\right)=L\left(N_{\mathit{\text{ech}}}^{\mathit{\text{LCV}}}+N_{\mathit{\text{ech}}}^{\mathit{\text{ADR}}};r_c\right)-N_{\mathit{\text{ech}}}^{\mathit{\text{ADR}}}{\overline{r}}_A\frac{\pi }{3S}-2r_{0}S+\pi r_0$$(31) (see Supplemental Information). (31) $$D_{\mathrm{B}}^{\mathit{\text{LCV}}}\left(r_0\right)=L\left(N_{\mathit{\text{ech}}}^{\mathit{\text{LCV}}}+N_{\mathit{\text{ech}}}^{\mathit{\text{ADR}}};r_c\right)-N_{\mathit{\text{ech}}}^{\mathit{\text{ADR}}}{\overline{r}}_A\frac{\pi }{3S}-2r_{0}S+\pi r_0$$(31) where $S$ is the number of subsectors of the traveling salesman tour considering the ring-radial metric and demand decreasing with distance to the city center, as defined by del Castillo (1999).

The LCV expected capacity ${\mathrm{\Psi }}_{\mathit{\text{ech}}}^{\mathit{\text{LCV}}}$ [parcels] is estimated as (32) $${\mathrm{\Psi }}_{\mathit{\text{ech}}}^{\mathit{\text{LCV}}}=\mathit{\text{min}}\left\{N_{\mathit{\text{ech}}}^{\mathit{\text{LCV}}};C^{\mathit{\text{LCV}}};\frac{H-2\frac{d_B^{\mathit{\text{LH}}}}{v_{\mathit{\text{LH}}}^{\mathit{\text{LCV}}}}}{\frac{1}{v_L^{\mathit{\text{LCV}}}}\frac{D_B^{\mathit{\text{LCV}}}}{N_{\mathit{\text{ech}}}^{\mathit{\text{LCV}}}}+{\tau }_d^{\mathit{\text{LCV}}}}\right\}$$(32)

The total distance traveled daily by the LCV fleet in this two-echelon distribution pattern $D_{\mathit{\text{ech}}}^{\mathit{\text{LCV}}}$ [veh-km/day] is then (33) $$D_{\mathit{\text{ech}}}^{\mathit{\text{LCV}}}\left(r_0\right)=D_B^{\mathit{\text{LCV}}}+2d_B^{\mathit{\text{LH}}}\frac{N_{\mathit{\text{ech}}}^{\mathit{\text{LCV}}}}{{\mathrm{\Psi }}_{\mathit{\text{ech}}}^{\mathit{\text{LCV}}}}$$(33)

The total LCV daily working time $T_{\mathit{\text{ech}}}^{\mathit{\text{LCV}}}$ [veh-h/day] in this two-echelon distribution scheme is estimated as (34) $$T_{\mathit{\text{ech}}}^{\mathit{\text{LCV}}}\left(r_0\right)=\frac{1}{v_L^{\mathit{\text{LCV}}}}{D}_B^{\mathit{\text{LCV}}}+\frac{1}{v_{\mathit{\text{LH}}}^{\mathit{\text{LCV}}}}2d_B^{\mathit{\text{LH}}}\frac{N_{\mathit{\text{ech}}}^{\mathit{\text{LCV}}}}{{\mathrm{\Psi }}_{\mathit{\text{ech}}}^{\mathit{\text{LCV}}}}+N_{\mathit{\text{ech}}}^{\mathit{\text{LCV}}}{\tau }_d^{\mathit{\text{LCV}}}$$(34)

Finally, the total daily operation costs $Z_{\mathit{\text{ech}}}$ [€/day] of the two-echelon distribution scheme are computed as (35a) $$Z_{\mathit{\text{ech}}}=Z_{\mathit{\text{ech}}}^{\mathit{\text{ADR}}}+Z_{\mathit{\text{ech}}}^{\mathit{\text{HDV}}}+Z_{\mathit{\text{ech}}}^{\mathit{\text{LCV}}}+{\mathrm{\Omega }}_h$$(35a) (35b) $$Z_{\mathit{\text{ech}}}^{\mathit{\text{ADR}}}=c_d^{\mathit{\text{ADR}}}{D}_{\mathit{\text{ech}}}^{\mathit{\text{ADR}}}+c_t^{\mathit{\text{ADR}}}{T}_{\mathit{\text{ech}}}^{\mathit{\text{ADR}}}$$(35b) (35c) $$c_d^{\mathit{\text{ADR}}}={\gamma }_{\mathit{\text{maintenance}}}^{\mathit{\text{ADR}}}+{\beta }_{{\mathit{\text{mech}}}^e}^{\mathit{\text{ADR}}}$$(35c) (35d) $$c_t^{\mathit{\text{ADR}}}={\epsilon }^{\mathit{\text{ADR}}}+{\gamma }_{\mathit{\text{insurance}}}^{\mathit{\text{ADR}}}+{\gamma }_{\mathit{\text{structure}}}^{\mathit{\text{ADR}}}+P_{e^e}^{\mathit{\text{ADR}}}$$(35d) (35e) $$Z_{\mathit{\text{ech}}}^{\mathit{\text{HDV}}}=c_d^{\mathit{\text{HDV}}}{D}_{\mathit{\text{ech}}}^{\mathit{\text{HDV}}}+c_t^{\mathit{\text{HDV}}}{T}_{\mathit{\text{ech}}}^{\mathit{\text{HDV}}}$$(35e) (35f) $$Z_{\mathit{\text{ech}}}^{\mathit{\text{LCV}}}=c_d^{\mathit{\text{LCV}}}{D}_{\mathit{\text{ech}}}^{\mathit{\text{LCV}}}+c_t^{\mathit{\text{LCV}}}{T}_{\mathit{\text{ech}}}^{\mathit{\text{LCV}}}$$(35f) (35g) $${\mathrm{\Omega }}_h=N_{\mathit{\text{ech}}}^{\mathit{\text{ADR}}}{a}^p{\omega }_{\mathit{\text{inf}}}\left(r_h\right)+\frac{T_{\mathit{\text{ech}}}^{\mathit{\text{ADR}}}}{{\kappa }^{\mathit{\text{ADR}}}}{\omega }_p$$(35g) where $c_d^{\mathit{\text{ADR}}}$[€/veh-km] and $c_t^{\mathit{\text{ADR}}}$ [€/veh-h] are the ADR expected unit distance and unit time costs, ${\gamma }_{\mathit{\text{maintenance}}}^{\mathit{\text{ADR}}}$ [€/veh-km] the ADR unit distance maintenance costs, ${\gamma }_{\mathit{\text{insurance}}}^{\mathit{\text{ADR}}}$ [€/veh-h] the ADR unit time insurance costs, ${\gamma }_{\mathit{\text{structure}}}^{\mathit{\text{ADR}}}$ [€/veh-h] the ADR unit time structure costs (Observatory of road freight transport in Catalonia, 2019), $c_d^{\mathit{\text{HDV}}}$ [€/veh-km] and $c_t^{\mathit{\text{HDV}}}$ [€/veh-h] the HDV expected unit distance and unit time costs, ${\mathrm{\Omega }}_h$ [€/day] the expected daily operation costs of the logistics micro-hub, $a^p$ [m²] the unit area occupied by one parcel within the logistics micro-hub, ${\omega }_{\mathit{\text{inf}}}$ [€/m²/day] the logistics micro-hub daily unit area infrastructure cost, ${\kappa }^{\mathit{\text{ADR}}}$ [ADRs] the number of ADRs an operator can monitor and supervise, ${\omega }_p$ [€/h] the unit time personnel costs, and ${\$}_e$ [€/kWh] the electricity price.

4. Case study

As presented in the modeling framework, the operation costs in the BAU and two-echelon scenarios depend on several input parameters, both from the city and vehicle characteristics. Now that the mathematical formulation has been explained, the objective of the rest of the article is to present some numerical insight.

To have a broader vision of ADR operations and not restrain our study to a unique use case, a Monte Carlo analysis considering several European cities will be conducted. In this section, we will explain the preparatory work we performed to create a representative data sample on which to perform the Monte Carlo analysis.

4.1. City selection

First of all, all European functional urban areas (FUAs) identified by the OECD (2019) were gathered. An important distinction has to be made between the FUA’s core center (also called the city) and its commuting zone. The modeling approach presented in the previous section deals with the optimization of last-mile operations within the FUA’s city, not the commuting zone.

The total area $A_c$ [km²] of each city is computed using QGIS 3.22.4 (QGIS Development Team, 2022). Some cities were removed from the data set because their shape could hardly be modeled with an idealized circular disk (see previous section) without an important loss of representativeness. The cities that were removed were essentially coastal ones because the hypothesis of problem variables independent from the angular coordinate is too far from reality. The filtering process was done manually, following subjective criteria.

Finally, only FUAs with a population inferior to 250,000 inhabitants (i.e., small and medium-sized FUAs; Dijkstra et al., 2019) were considered in our sample. For metropolitan and large metropolitan FUAs (with a total population greater than 250,000 inhabitants), the assumption of having a unique logistics micro-hub in the city center is not representative enough. To serve a large population, several urban distribution hubs would be needed. Figure 3 shows the location of the 164 cities that were finally used in this article. The most represented countries are Poland, Italy, and Germany, with 34, 26, and 16 cities, respectively (see Supplemental Information).

Figure 3. Small and medium European cities considered in the sample.

Figure 4 presents the city population PDF. The city mean population within our sample is 79,850 inhabitants, and the mean city area is 96 km².

Figure 4. (a) City population PDF and (b) city area PDF, per country.

After the filtering process, the equivalent city radius $r_c$ is estimated. We model the city as an idealized disk (see previous section). (36) $$r_c=\sqrt{\frac{A_c}{\pi }}$$(36)

4.2. Carrier demand density estimation

In each city, the total number of deliveries per day generated by B2C e-commerce is estimated considering the population aged between 16 and 74 years (Eurostat, 2022b), combined with the e-commerce demand rate per inhabitant per day (Eurostat, 2021). The parcel generation rate we obtain from Eurostat (2021) is equal to 0.026 parcels per inhabitant per day, on average at the European scale, which corresponds to approximately 6 parcels per inhabitant per year (assuming 250 working days per year), which is in line with the literature (Figliozzi & Unnikrishnan, 2021). In this article, we decided to focus on e-commerce parcel delivery operations. Other use cases, such as food or grocery delivery, are out of scope.

Once the total number of deliveries per day has been computed, the carrier’s competitive market in the city has to be modeled. We assume that the market share $p_k^{\mathit{\text{fua}}}$ of carrier $k$ in the FUA $\mathit{\text{fua}}$ follows a Poisson distribution. Carriers with a market share lower than 0.05 (5%) are not considered in this article. (37) $$p_k^{\mathit{\text{fua}}}=\frac{e^{-\lambda }{\lambda }^k}{k!};\lambda =‐\mathrm{\text{ln}}\left(p_0\right)$$(37) where $p_0$ is the market share of the national post operator in the country where the FUA is located (European Commission et al., 2018).

Once the number of deliveries per carrier, obtained by multiplying the carrier’s market share by the total number of deliveries in the city, has been estimated, the last step is to formulate the spatial distribution of the demand density ${\delta }_c(r).$ As a reminder, we assumed during the modeling process that this demand density decreases with distance to the city center $r.$ We assume that the demand density is proportional to the population density and that it can be modeled with a negative exponential spatial PDF (von Bergmann, 2019). (38) $${\delta }_c\left(r\right)={\delta }_0{e}^{-r/d};\frac{P_c}{P_t}=\frac{d-\left(d+r_c\right)e^{-r_c/d}}{d-\left(d+R\right)e^{-R/d}};{\delta }_0=\frac{N_{t,c}}{2\pi d\left(d-R{e}^{-R/d}-d{e}^{-R/d}\right)}$$(38) where $P_c$ [inhabitants] is the city population, $P_t$ [inhabitants] the FUA total population including the commuting zone, $R$ [km] the FUA radius, $N_{t,c}$ [parcels/day] the daily total number of parcels delivered by carrier $c$ in the whole FUA, and $d$ [km] a decrease rate that characterizes the concentration of recipients in the core city compared to the commuting zone.

Figure 5 presents an illustration of the carrier demand spatial PDF in 3 cities taken from our sample.

Figure 5. Carrier demand spatial PDF. Examples in different European cities.

4.3. Logistics hub infrastructure cost modeling

The objective of this subsection is to present the modeling of the logistics micro-hub unit area daily infrastructure cost ${\omega }_{\mathit{\text{inf}}}(r),$ which depends on the distance to the city center $r.$ We assume that the real estate price could be described as an exponential function of the distance to the city center (Ackland & Wargentin, 2014; d'Acci, 2019; Lukavec & Kadeřábková, 2017; Rosenthal et al., 2022). As a consequence, ${\omega }_{\mathit{\text{inf}}}(r)$ can be described by an exponential function of $r.$ In addition, considering the lack of available data, we assume that the exponential decrease rate of ${\omega }_{\mathit{\text{inf}}}(r)$ is equal to $d$ (exponential decrease rate of the demand density ${\delta }_c(r),$ see previous subsection). (39) $${\omega }_{\mathit{\text{inf}}}\left(r\right)=\frac{{\gamma }_c-{\gamma }_o}{1-e^{-r_c/d}}\left(e^{-r/d}-1\right)+{\gamma }_c$$(39) where ${\gamma }_c$ [€/m²/day] is the micro-hub daily infrastructure cost in the city center and ${\gamma }_o$ [€/m²/day] the micro-hub daily infrastructure cost in the city outskirts (see Figure 6; Numbeo, 2022).

Figure 6. Micro-hub unit area daily infrastructure cost as a function of distance to city center. Examples in different European cities.

To compute ${\gamma }_c$ and ${\gamma }_o,$ we use the “Price per Square Meter to Buy Apartment in City Centre” and the “Price per Square Meter to Buy Apartment Outside of Centre” (Numbeo, 2022) as proxies to evaluate the cost of implementing a logistics micro-hub. To estimate the expected daily price of such an infrastructure, we assume an investment of 10 years (with 250 working days per year).

Finally, concerning the urban micro-hub, we assume that its angular coordinate is null (${\theta }_h$ = 0, see the previous subsection) for all results presented later on. This means that the urban micro-hub is located on a straight line that joins the carrier’s DC and the city center. We make this assumption for two reasons. Firstly, when implementing an urban micro-hub in a given city, a logistics company will look for the optimum location to decrease total operation costs. As we previously assumed that the micro-hub infrastructure costs only depend on the distance to the city center, and not the angular coordinate, placing the micro-hub between the DC and the city center will minimize the HDV total operation costs in the two-echelon delivery scheme. Secondly, as it will be demonstrated in the next section, the value of ${\theta }_h$ has a very limited impact on the two-echelon total delivery costs. Our model is quite insensitive to the value of this parameter. As a consequence, assuming that ${\theta }_h$ = 0 will reduce the complexity of the model while ensuring the representativeness of the obtained results.

4.4. Vehicle operation characteristics

All the vehicle operational input parameters considered in our use case are gathered in Table 3. We assume that the LCV volume capacity $C^{\mathit{\text{LCV}}}$ is equal to 150 parcels (Renault Renault, 2023), and the HDV volume capacity $C^{\mathit{\text{HDV}}}$ is equal to 1,000 parcels (StreetScooter, 2023). As a matter of illustration, 150 parcels with an average volume of 0.03 m³ (30 cm × 30 cm × 30 cm approximately) and an average weight of 3 kg (Perboli & Rosano, 2019) occupy a volume of 4.5 m³ and have a total weight of 450 kg, which is compatible with the characteristics of a delivery van (Renault Renault, 2023). We assume that the expected commercial line-haul speed for both HDVs and LCVs is equal to 50 km/h for all cities. Similarly, the expected commercial speed in the local urban grid for both HDVs and LCVs is assumed to be 20 km/h for all cities. To estimate the expected line-haul and local commercial speeds of LCVs and HDVs, we used data from 5 line-haul and 8 local driving cycles and computed the average speed across both driving cycle types (see Table 4; DieselNet, 2023; NREL, 2023).

Table 3. Parameter and cost assumptions.

Parameter	Value assumption	Reference
Carrier daily demand density ${\delta }_c(r)$ [receivers/km^2/day]	Follows a negative exponential PDF which depends on the considered city and carrier	(European Commission et al., 2018; Eurostat, 2021, 2022b; von Bergmann, 2019)
City radius $r_c$ [km]	Depends on the considered city	(OECD, 2019)
Line-haul distance between the carrier’s DC and city center $l_{\mathit{\text{DC}}}$ [km]	Follows a uniform PDF between 10 and 30	(Heitz & Dablanc, 2015)
Operation time window $H$ [h]	8	(Allen et al., 2018)
Micro-hub daily unit area infrastructure cost ${\omega }_{\mathit{\text{inf}}}$ [€/m^2/day]	Follows a negative exponential PDF which depends on the considered city	(Ackland & Wargentin, 2014; d'Acci, 2019; Lukavec & Kadeřábková, 2017; Numbeo, 2022; Rosenthal et al., 2022)
Unit time personnel cost ${\omega }_p$ [€/h]	Depends on the considered city	(Numbeo, 2022; Observatory of road freight transport in Catalonia, 2019)
Unit area $a^p$ occupied by one parcel at the logistics micro-hub [m^2/parcel]	1	(Navarro et al., 2016)
LCV volume capacity $C^{\mathit{\text{LCV}}}$ [parcels]	150	(Perboli & Rosano, 2019; Renault Renault, 2023)
LCV line-haul speed $v_{\mathit{\text{LH}}}^{\mathit{\text{LCV}}}$ [km/h]	50	(DieselNet, 2023; NREL, 2023)
LCV local speed $v_L^{\mathit{\text{LCV}}}$ [km/h]	20	(DieselNet, 2023; NREL, 2023)
LCV stop time per delivery ${\tau }_d^{\mathit{\text{LCV}}}$ [min]	4	(Allen et al., 2018; Perboli & Rosano, 2019)
LCV unit distance operation cost $c_d^{\mathit{\text{LCV}}}$ [€/veh-km]	Depends on the considered city	(Numbeo, 2022; Observatory of road freight transport in Catalonia, 2019)
LCV unit time operation cost $c_t^{\mathit{\text{LCV}}}$ [€/veh-h]	Depends on the considered city	(Numbeo, 2022; Observatory of road freight transport in Catalonia, 2019)
ADR volume capacity $C^{\mathit{\text{ADR}}}$ [parcels]	15	Own prototype of road ADR
ADR commercial speed $v^{\mathit{\text{ADR}}}$ [km/h]	10	(Lemardelé et al., 2023)
ADR stop time per delivery ${\tau }_d^{\mathit{\text{ADR}}}$ [min]	2	(Perboli & Rosano, 2019)
ADR unit time depreciation cost ${\epsilon }^{\mathit{\text{ADR}}}$ [€/veh-h]	1	Variable value used for sensitivity analysis
ADR insurance costs ${\gamma }_{\mathit{\text{insurance}}}^{\mathit{\text{ADR}}}$ [€/veh-h]	Assumed to be equal to LCV insurance costs in the considered city	(Numbeo, 2022; Observatory of road freight transport in Catalonia, 2019)
ADR maintenance costs ${\gamma }_{\mathit{\text{maintenance}}}^{\mathit{\text{ADR}}}$ [€/veh-km]	Assumed to be equal to LCV maintenance costs in the considered city	(Numbeo, 2022; Observatory of road freight transport in Catalonia, 2019)
ADR structure costs ${\gamma }_{\mathit{\text{structure}}}^{\mathit{\text{ADR}}}$ [€/veh-h]	Assumed to be equal to LCV structure costs in the considered city	(Numbeo, 2022; Observatory of road freight transport in Catalonia, 2019)
ADR unit distance energy consumption rate for mechanical propulsion ${\beta }_{\mathit{\text{mech}}}^{\mathit{\text{ADR}}}$ [Wh/km]	20	(Lemardelé et al., 2023)
ADR electronics and sensor power $P_e^{\mathit{\text{ADR}}}$ [W]	500	(Lemardelé et al., 2023)
ADR battery capacity ${\mathit{\text{BC}}}^{\mathit{\text{ADR}}}$ [kWh]	3.75 (assuming a rate of discharge of 80%)	(Lemardelé et al., 2023)
Number of ADRs an operator is able to monitor and supervise ${\kappa }^{\mathit{\text{ADR}}}$ [ADRs]	5	(National Academies of Sciences, Engineering, and Medicine, 2023)
HDV volume capacity $C^{\mathit{\text{HDV}}}$ [parcels]	1,000	(Perboli & Rosano, 2019; StreetScooter, 2023)
HDV line-haul speed $v_{\mathit{\text{LH}}}^{\mathit{\text{HDV}}}$ [km/h]	50	(DieselNet, 2023; NREL, 2023)
HDV loading/unloading time per parcel ${\tau }_{\mathit{\text{LU}}}^{\mathit{\text{HDV}}}$ [s]	45	(Dell’Amico & Hadjidimitriou, 2012; Hofmann et al., 2017)
HDV unit distance operation cost $c_d^{\mathit{\text{HDV}}}$ [€/veh-km]	Depends on the considered city	(Numbeo, 2022; Observatory of road freight transport in Catalonia, 2019)
HDV unit time operation cost $c_t^{\mathit{\text{HDV}}}$ [€/veh-h]	Depends on the considered city	(Numbeo, 2022; Observatory of road freight transport in Catalonia, 2019)
Price of electricity ${\$}_e$ [€/kWh]	Depends on the country where the city is located	(Eurostat, 2022a)

Table 4. Line-haul and local driving cycles.

Driving cycle name	Driving cycle type	Average speed (km/h)
CADC line-haul	Line-haul	57.5
EPA UDDS line-haul	Line-haul	41.3
NEDC line-haul	Line-haul	66.4
NREL parcel delivery truck	Line-haul	30.7
WLTC line-haul	Line-haul	56.7
Braunschweig city driving cycle	Local	30.1
CADC urban	Local	17.6
CSHVC	Local	22.2
EPA UDDS urban	Local	25.8
NEDC urban	Local	18.1
NREL Class 6 EV	Local	25.1
NYCC	Local	11.4
WLTC urban	Local	18.9

The expected LCV stop time per delivery ${\tau }_d^{\mathit{\text{LCV}}}$ is equal to 4 min (Allen et al., 2018; Perboli & Rosano, 2019). We assume that the expected HDV loading/unloading time per parcel ${\tau }_{\mathit{\text{LU}}}^{\mathit{\text{HDV}}}$ is equal to 45 s (Dell’Amico & Hadjidimitriou, 2012; Hofmann et al., 2017).

As for ADRs, we consider the characteristics of a prototype developed by our team, for which data is available. Its volume capacity is approximately 0.5 m³ so, considering an expected parcel volume of 0.03 m³, the ADR volume capacity $C^{\mathit{\text{ADR}}}$ is equal to 15 parcels. Assuming an average parcel weight of 3 kg, the total weight of 15 parcels (which corresponds to a fully loaded robot) would be equal to 45 kg, which is consistent with the characteristics of our ADR. Nevertheless, as the value of $C^{\mathit{\text{ADR}}}$ is dependent on the actual design configuration of the prototype and also to increase the representativeness of our results, we will perform a sensitivity analysis with respect to the value of this parameter in the next section. In addition, we consider an expected ADR unit distance energy consumption rate for mechanical propulsion ${\beta }_{\mathit{\text{mech}}}^{\mathit{\text{ADR}}}$ of 20 Wh/km (Lemardelé et al., 2023), and an ADR electronics and sensor power $P_e^{\mathit{\text{ADR}}}$ of 500 W (Lemardelé et al., 2023). Finally, the ADR battery capacity ${\mathit{\text{BC}}}^{\mathit{\text{ADR}}}$ is assumed to be equal to 3.75 kWh, including a rate of discharge of 80% (Lemardelé et al., 2023). We assume the ADR commercial speed $v^{\mathit{\text{ADR}}}$ to be equal to 10 km/h (Lemardelé et al., 2023), and its stop time per delivery ${\tau }_d^{\mathit{\text{ADR}}}$ 2 min (Perboli & Rosano, 2019). This prototype is a road ADR. Sidewalk robots are out of scope for this study. Concerning the delivery process, ADRs are not able to directly access the final receivers’ doorsteps. The consignee receives a message through a mobile application stating that the ADR is at its destination, ready for delivery. The consignee must then come out and pick up his or her parcel.

The expected unit time operation costs of non-autonomous vehicles are assumed to be the sum of the personnel costs, vehicle depreciation costs, vehicle insurance costs, diet costs, and carrier’s structure costs. For LCVs and HDVs, the value of these different costs is obtained using data from the Observatory of road freight transport in Catalonia (2019) and Numbeo (2022) (see Supplemental Information). The expected unit distance costs of each vehicle are assumed to be the sum of the fuel and maintenance costs, updated in each city using Numbeo (2022) (see Supplemental Information). In the case of ADRs, because of a lack of available operational data, we assume that their insurance, structure, and maintenance costs are equal to LCVs’. The value of ADR insurance costs in future years is highly uncertain (Li & Kunze, 2023). It is true that ADRs present a lower level of danger (smaller vehicles traveling at a reduced speed when compared to conventional delivery vans), but the novelty of autonomous technologies used in ADRs may make their insurance costs increase. The same phenomenon might occur with maintenance operations. Even if electric vehicles are assumed to present lower maintenance costs than internal combustion engine ones (Propfe et al., 2012), the usage of autonomous technologies may require more maintenance from a software perspective, especially at the beginning when the technology needs more maturity. For all these reasons, we assume that ADR maintenance, insurance, and structure costs (including administrative expenses, IT network costs, etc.) are equal to LCVs. These hypotheses should be validated in future years when more large-scale operational data are available within the research community. At the moment, this information is sensitive, and robot manufacturers and operators are very secretive in this regard.

The impact of temperature on battery performance has not been considered in this article. In addition, we assume that the road gradient is null in all cities. This aspect, which is highly relevant in hilly cities, is left for further research.

We assume that the expected distance $l_{\mathit{\text{DC}}}$ between the carrier’s DC and the city center follows a uniform PDF between 10 and 30 km (Heitz & Dablanc, 2015). We assume that the operation time window $H$ is equal to 8 h (Allen et al., 2018). Finally, the unit area occupied by one parcel within the logistics micro-hub $a^p$ is assumed to be 1 m² per parcel (Navarro et al., 2016).

5. Results

We will now go deeper into numerical results concerning the optimization of ADR last-mile operations. Figures 7–9 present the operation delivery cost ratio $Z_{\mathit{\text{ech}}}/Z_{\mathit{\text{BAU}}}$ as a function of the micro-hub service radius $r_0$ and distance $r_h$ between the micro-hub and the city center. The situation in three cities is presented.

Figure 7. Operation delivery cost ratio $Z_{\mathit{\text{ech}}}/Z_{\mathit{\text{BAU}}}$ as a function of micro-hub delivery zone radius $r_0$ and distance $r_h$ between the logistics micro-hub and city center in Roanne.

Figure 8. Operation delivery cost ratio $Z_{\mathit{\text{ech}}}/Z_{\mathit{\text{BAU}}}$ as a function of micro-hub delivery zone radius $r_0$ and distance $r_h$ between the logistics micro-hub and city center in Legnica.

Figure 9. Operation delivery cost ratio $Z_{\mathit{\text{ech}}}/Z_{\mathit{\text{BAU}}}$ as a function of micro-hub delivery zone radius $r_0$ and distance $r_h$ between the logistics micro-hub and city center in Potenza.

In Figures 7–9, we indicate with $\oplus$ the values of $\left(r_0^{*};r_h^{*}\right)$ for which the ratio $Z_{\mathit{\text{ech}}}/Z_{\mathit{\text{BAU}}}$ is minimum.

It seems reasonable to assume that the total operation costs $Z_{\mathit{\text{ech}}}$ can be optimized depending on the values of $(r_0,$ $r_h).$ When plotting the delivery cost ratio $Z_{\mathit{\text{ech}}}/Z_{\mathit{\text{BAU}}}$ as a function of $(r_0,$ $r_h)$ (see Figures 7–9), a minimum appears, above all because the solution search space $(r_0,$ $r_h)$ is bounded. As a matter of illustration, in Roanne, considering a line-haul distance $l_{\mathit{\text{DC}}}$ = 11 km and 646 deliveries, an optimal value of $Z_{\mathit{\text{ech}}}$ is reached for $r_h$ = 0 km (the logistics micro-hub is located at the city center) and $r_0$ = 3.2 km approximately (see Figure 7). In Legnica, for $l_{\mathit{\text{DC}}}$ = 26 km and 498 deliveries, the minimum two-echelon operation costs are reached for $r_h$ = 1 km and $r_0$ = 3.5 km, approximately. In Potenza, the minimum two-echelon operation costs are reached for $r_h$ = 2 km and $r_0$ = 5 km, i.e., the maximum value of $r_0$ we considered in this study (see Figure 9). Here it is necessary to emphasize that the optimization function used to find the minimum value of $Z_{\mathit{\text{ech}}}$ has some limitations. In the case of Legnica (see Figure 8), the optimal value of $Z_{\mathit{\text{ech}}}$ found by the minimization algorithm and represented by $\oplus$ does not correspond to the actual global minimum, which is located at ($r_0$ = 1 km, $r_h$ = 1 km) approximately, i.e., the numerical algorithm is not able to find the global minimum in this particular case. This may alter the overall results of the article if this underestimation occurs in many circumstances. The relative difference between the global minimum and the actual minimum estimated by the numerical minimization algorithm is less than 8% (see Figure 9). The improvement of the minimization numerical methods is beyond the scope of this article and is left for further research.

As depicted in Figures 7–9, the two-echelon delivery costs $Z_{\mathit{\text{ech}}}$ can be minimized in Roanne, Legnica and Potenza. The same result was observed in numerous evaluation cases, and we assume this result can be generalized to all cities and delivery configurations.

The delivery cost reduction depends on the particular configuration in each city. In Figures 7–9, the ratio $Z_{\mathit{\text{ech}}}/Z_{\mathit{\text{BAU}}}$ is displayed as a function of $r_0$ considering that $r_h$ = 0 km, i.e., assuming the logistics micro-hub is located at the city center. For all cities highlighted in these figures, the same pattern can be observed. While the micro-hub delivery radius $r_0$ increases, the two-echelon LCV costs decrease while the ADR, HDV, and logistics micro-hub costs increase. This pattern was expected a priori. The key element to determining if the two-echelon delivery scheme is able to reduce last-mile operation costs is to analyze if the LCV cost decrease is balanced by the micro-hub, ADR, and HDV cost increases. In the city of Roanne, the LCV cost decrease is not compensated by the hub and ADR cost increases, i.e., the total two-echelon delivery costs (which are the sum of these three elements) decrease when $r_0$ increases. On the contrary, in the case of Legnica (see Figure 8), the LCV cost decrease is fully compensated by the logistics micro-hub and ADR cost increases (the micro-hub cost increase being the main contributor). To explain the tradeoff at stake, many input parameters of diverse nature have to be considered. Especially the comparison between personnel and micro-hub infrastructure costs is expected to play a key role in the two-echelon delivery scheme’s performance. Some explanatory variables will be identified in the following.

In Figures 7–9, the delivery cost ratio $Z_{\mathit{\text{ech}}}/Z_{\mathit{\text{BAU}}}$ is also displayed as a function of $r_h.$ For low values of $r_0$ (i.e., $r_0$ = 1 km), two-echelon delivery costs are not sensitive to the value of $r_h$ (see Supplemental Information) because they are essentially due to LCV operations, which are independent from $r_h.$ For higher values of $r_0,$ two-echelon delivery costs are sensitive to the value of $r_h,$ but the situation depends on the considered city. In the case of Potenza and Legnica, the two-echelon total operation cost minimum is reached for an intermediate value of $r_h$ (around 1 km for Legnica and 2 km for Potenza). For low values of $r_h,$ the micro-hub cost decreases when $r_h$ increases because the micro-hub unit area infrastructure cost decreases when going away from the city center and this is the biggest driver for the micro-hub cost decrease. Nevertheless, if the value of $r_h$ keeps increasing, the micro-hub costs start increasing because the ADR total working time increases (ADRs have to travel a greater distance to access final receivers), i.e., the ADR fleet size increases, and more personnel are needed to supervise operations. The micro-hub personnel cost increase compensates for the infrastructure cost decrease when going away from the city center. In addition, when $r_h$ increases, the total two-echelon delivery costs increase because of the ADR cost increase. Delivery robots have to travel a greater distance to access final recipients for higher values of $r_h.$

For a more in-depth analysis of three representative cities (including the estimation of several delivery key performance indicators), the reader can refer to the Supplemental Information.

In the study, we decided to limit the value of $r_0$ to 5 km. If the city radius $r_c$ is inferior to 5 km, the value of $r_0$ is also limited by $r_c.$ Assuming that the micro-hub is located exactly at the city center, considering that $r_0$ is inferior to 5 km means that each ADR can be reached at any moment within 15 min (assuming a speed of 20 km/h, which is the local urban grid speed we considered previously). This is important for monitoring delivery operations and ensuring safety. If a safety breach is foreseen, for instance, a supervision operator located at the micro-hub can reach any ADR within 15 min, which seems a reasonable intervention time to the authors.

In the rest of the article, the values of $r_0$ and $r_h$ that minimize two-echelon delivery operations will be noted $r_0^{*}$ and $r_h^{*}.$ The values of $r_0^{*}$ and $r_h^{*}$ are obtained using the “minimize” function of the Scipy Python library, considering that $r_0$ is bounded between 1 km and $\mathit{\text{min}}\left\{5;0.99r_c\right\}$ km. The value of $r_h$ is also bounded between 0 and 5 km. In the optimization process, we use a tolerance of 10⁻⁴ and 15 as the maximum number of iterations.

Figure 10 shows the decomposition of the ratio $Z_{\mathit{\text{ech}}}^{*}/Z_{\mathit{\text{BAU}}}$ per vehicle type (either HDVs, LCVs, or ADRs) and micro-hub contribution for a given sample of cities. The BAU delivery cost per parcel delivery is also indicated.

Figure 10. Decomposition of the optimized delivery cost ratio $z_{\mathit{\text{ech}}}^{*}/z_{\mathit{\text{BAU}}}$ per vehicle type and micro-hub contribution for several cities.

In the sample, the BAU operation costs range from 1.4 €/delivery to 5.9 €/delivery. In some cities (see Reus), the delivery cost reduction reaches 23% because city configuration is favorable to ADRs. On the other hand, in some cases (see Hradec Králové, for instance), delivery configuration is not favorable to ADRs, i.e., no significant delivery cost reduction is observed when implementing an urban micro-hub and robotic devices. In the rest of the article, we will give some insight into the configurations that are favorable to ADRs or not.

In the city sample considered in Figure 10, ADR costs represent between approximately 0 (see Hradec Králové) and 20% (see Satu Maru) of the two-echelon total operation costs. In Hradec Králové, implementing a micro-hub and ADRs would not significantly reduce operation costs, i.e., almost all delivery costs are incurred by LCVs. In favorable conditions (see Salamanca, for instance), most parcels are distributed with ADRs, and LCV operations are not needed any more. On the contrary, if ADRs are not able to reduce delivery operations (see Hradec Králové), LCV operation costs represent the vast majority of total operation costs, even in the two-echelon delivery scheme (ADRs are not used at all).

Another main component of two-echelon total operation costs is the micro-hub contribution, which can represent up to 73% of the two-echelon total costs (see the city of Lucerne with 627 parcels per day, for instance). In use cases with numerous ADRs, the number of parcels that are transhipped from HDVs to ADRs is higher, i.e., the micro-hub needs to be bigger, which generates higher micro-hub infrastructure costs. In the case of numerous ADRs operating in the city, the micro-hub personnel costs are also higher because more operators are needed to supervise the robot fleet. As a matter of illustration, in the city of Villingen, considering a total number of 1,200 deliveries in the whole city, approximately 700 parcels would be delivered through a 700-m² micro-hub in the optimal two-echelon configuration (see Supplemental Information). The remaining 500 parcels (in the city’s external ring) would be delivered with LCVs. This 700-m² urban micro-hub is related to the concept of logistics hotels, as implemented in the city of Paris, whose area can range from around 1,000 m² to more than 20,000 m² (Dablanc, 2019). In Villingen, around 20 ADRs are needed to deliver the 700 parcels that are transhipped at the micro-hub, i.e., the micro-hub would be run by 4 operators since we assumed that an operator is able to supervise 5 ADRs.

Another interesting aspect is to analyze the sensitivity of the two-echelon delivery cost $z_{\mathit{\text{ech}}}^{*}$ to the line-haul distance $l_{\mathit{\text{DC}}}.$ In Figure 11, the cost ratio $z_{\mathit{\text{ech}}}^{*}/z_{\mathit{\text{BAU}}}$ is displayed as a function of $l_{\mathit{\text{DC}}}$ for several cities. The first interesting result is that, in most cases (see Piotrków Trybunalski for the exception), two-echelon delivery costs decrease when $l_{\mathit{\text{DC}}}$ increases. For all cities displayed in Figure 11, except Piotrków Trybunalski, increasing $l_{\mathit{\text{DC}}}$ from 10 to 30 km generates a 10% decrease in the cost ratio $z_{\mathit{\text{ech}}}^{*}/z_{\mathit{\text{BAU}}},$ approximately.

Figure 11. Optimized delivery cost ratio $z_{\mathit{\text{ech}}}^{*}/z_{\mathit{\text{BAU}}}$ as a function of the line-haul distance $l_{\mathit{\text{DC}}}$ for several cities.

These results were expected. If the carrier’s DC is located far away from the city center, LCV costs in the BAU delivery scheme are higher because vans have to travel a higher distance to access final recipients, i.e., BAU delivery costs increase. In the meantime, in the two-echelon delivery scheme, full advantage is taken of the HDV economies of scale when the carrier’s DC is located far away from the city center. The combination of these two phenomena explains why cost reduction is higher for large values of $l_{\mathit{\text{DC}}}.$ In Ostrów Wielkopolski, there is a clear threshold of around 15 km. If the value of $l_{\mathit{\text{DC}}}$ is lower, BAU delivery costs are higher. On the contrary, if the value of $l_{\mathit{\text{DC}}}$ is higher, two-echelon delivery costs are lower. This highlights the importance of considering the logistics sprawl when planning last-mile operations in a particular city.

In the case of Piotrków Trybunalski, the cost reduction is independent of $l_{\mathit{\text{DC}}}$ because minimum two-echelon delivery costs are obtained for very low values of $r_0$ (see Supplemental Information). As a consequence, LCV operations are the main contributor to the total two-echelon delivery costs for all values of $l_{\mathit{\text{DC}}}$ and no cost reduction is observed. The BAU and two-echelon delivery schemes are almost equal because city configuration is not favorable to ADR operations.

Finally, we can also confirm that the assumption we made (${\theta }_h$ = 0 rad for all two-echelon delivery configurations) is valid because the error in the estimation of two-echelon delivery costs is lower than 1% (see Supplemental Information). This is also confirmed by Figure 10 because the share of HDV operation costs in two-echelon total delivery costs is lower than 20% and the main source of HDV costs is the truck unloading operations at the micro-hub. Our model is quite insensitive to the value of the micro-hub angular coordinate ${\theta }_h.$

The next step of the study is to perform a Monte Carlo analysis.

Figure 12 presents the results of the Monte Carlo analysis considering 10,000 iterations (computational cost of 1.7 h), an ADR depreciation ${\epsilon }^{\mathit{\text{ADR}}}$ = 1 €/h, and an ADR capacity $C^{\mathit{\text{ADR}}}$ = 15 parcels. We assume that an operator located at the micro-hub supervises 5 robots, i.e., ${\kappa }^{\mathit{\text{ADR}}}$ = 5 (National Academies of Sciences, Engineering, and Medicine, 2023).

Figure 12. Histogram of cost ratio $z_{\mathit{\text{ech}}}^{*}/z_{\mathit{\text{BAU}}}$ as a function of (a) micro-hub infrastructure/personnel cost ratio ${\omega }_{\mathit{\text{inf}}}\left(r_h=0\right)/{\omega }_p,$ (b) demand density $\delta =N_{t,c}/A_c,$ (c) distribution center distance/city radius ratio $l_{\mathit{\text{DC}}}/r_c,$ (d) city area $A_c,$ (e) personnel cost/ADR depreciation ratio ${\omega }_p/{\epsilon }^{\mathit{\text{ADR}}},$ and (f) total number of deliveries $N_{t,c}.$

The BAU mean cost is equal to 2.7 €per delivery, whereas the two-echelon mean delivery cost is 2.3 €per delivery. The mean cost reduction when passing from the BAU model to the two-echelon one is then 16%.

Figure 12 presents the histogram of the cost ratio $z_{\mathit{\text{ech}}}^{*}/z_{\mathit{\text{BAU}}}.$ In the different subplots, a color scale indicates the values of some relevant input parameters. The objective is to identify the main decision variables that would foster or prevent the use of a micro-hub and ADRs.

It seems that the most important input parameters to explain the variability of the results are the demand density $\delta$ (see Figure 12b), the infrastructure/personnel cost ratio ${\omega }_{\mathit{\text{inf}}}\left(r_h=0\right)/{\omega }_p$ (see Figure 12a), and the personnel cost/ADR depreciation ratio ${\omega }_p/{\epsilon }^{\mathit{\text{ADR}}}$ (see Figure 12e). We use ${\omega }_{\mathit{\text{inf}}}(r_h=0)$ (micro-hub unit area daily infrastructure cost at the city center) as a proxy to evaluate the intensity of the logistic micro-hub infrastructure costs. Operation cost reductions of around 40% are achieved in the context of a high demand density ($\delta$ > 11 parcels/km²/day) and values of the ratio ${\omega }_{\mathit{\text{inf}}}\left(r_h=0\right)/{\omega }_p$ below 0.04. This result seems to be consistent with an a priori analysis. If the relative weight of the micro-hub infrastructure cost is high compared to the personnel costs, the two-echelon delivery strategy will not be able to reduce operation costs, and the BAU delivery scheme with a delivery person remains more competitive.

Figure 12e indicates that higher delivery cost reductions are achieved in the context of high personnel costs. For values of the ratio ${\omega }_p/{\epsilon }^{\mathit{\text{ADR}}}$ above 23 (personnel costs 23 times higher than ADR depreciation), delivery costs can be reduced by 40%. City area and distance ratio $l_{\mathit{\text{DC}}}/r_c$ do not seem to have a significant impact on delivery cost reduction (see Figure 12c,d).

We will now perform a sensitivity analysis of the model to the ADR depreciation cost ${\epsilon }^{\mathit{\text{ADR}}},$ volume capacity $C^{\mathit{\text{ADR}}},$ and the number of ADRs monitored by an operator ${\kappa }^{\mathit{\text{ADR}}}.$ A variation in the ADR volume capacity $C^{\mathit{\text{ADR}}}$ would require a modification of the robot design to be able to load more parcels. Performing a sensitivity analysis of the model with respect to this parameter is interesting to validate the current robot design or give insight into future developments of our prototype.

Figure 13a presents the optimized two-echelon operation cost per parcel delivery $z_{\mathit{\text{ech}}}^{*}$ for different values of ${\epsilon }^{\mathit{\text{ADR}}},$ $C^{\mathit{\text{ADR}}},$ and ${\kappa }^{\mathit{\text{ADR}}}.$ Figure 13b shows the variation in last-mile operation costs when implementing a two-echelon delivery scheme. Finally, Figure 13c depicts the probability of having lower two-echelon operation costs than BAU ones $P\left(z_{\mathit{\text{ech}}}^{*}\le z_{\mathit{\text{BAU}}}\right).$ Each matrix cell in Figure 13 was computed using the Monte Carlo approach previously described with 1,000 iterations.

Figure 13. (a) Two-echelon mean operation costs per parcel delivery, (b) operation costs variation when implementing two-echelon delivery scheme, and (c) probability of obtaining lower two-echelon operations costs $P\left(z_{\mathit{\text{ech}}}^{*}\le z_{\mathit{\text{BAU}}}\right).$

Comparing $z_{\mathit{\text{ech}}}^{*}$ and $P\left(z_{\mathit{\text{ech}}}^{*}\le z_{\mathit{\text{BAU}}}\right)$ will help us discriminate between situations in which cost reduction is high but quite unlikely and situations in which cost reduction is lower but more likely to occur.

As a matter of illustration, let us have a look at Figure 14. The two-echelon operation cost PDF 1 has a much lower mean, but in some cases, the two-echelon operation costs are higher than the BAU ones. On the contrary, PDF 2 has a higher mean than PDF 1, but it is less uncertain than PDF 1 because its standard deviation is lower. Studying the mean operation costs as well as the probability of having lower two-echelon operation costs $P\left(z_{\mathit{\text{ech}}}^{*}\le z_{\mathit{\text{BAU}}}\right)$ will help us characterize the two-echelon total operation costs PDF.

Figure 14. Hypothetical two-echelon operation PDFs.

In the case of ${\kappa }^{\mathit{\text{ADR}}}$ = 5 ADRs per operator, depending on the values of ${\epsilon }^{\mathit{\text{ADR}}}$ and $C^{\mathit{\text{ADR}}},$ two-echelon delivery costs range from 2.1 €/delivery to 2.7 €/delivery. Considering that the BAU mean operation cost is equal to 2.7 €/delivery, the cost decrease ranges from −17% to null (see Figure 13b). As it was expected, the highest two-echelon mean operation costs occur for high values of ${\epsilon }^{\mathit{\text{ADR}}}$ and low values of $C^{\mathit{\text{ADR}}},$ i.e., high depreciation costs and a low number of visited receivers per ADR route. For instance, it can be seen in Figure 13a,b that having the combination (${\epsilon }^{\mathit{\text{ADR}}}$ = 0.5 €/h, $C^{\mathit{\text{ADR}}}$ = 5) approximately leads to the same two-echelon operation costs as the combination (${\epsilon }^{\mathit{\text{ADR}}}$= 4 €/h, $C^{\mathit{\text{ADR}}}$ = 20). Globally, in the case of ${\kappa }^{\mathit{\text{ADR}}}$ = 5 ADRs per operator, mean operation costs are reduced by more than 10% for $C^{\mathit{\text{ADR}}}\ge$ 15 receivers per route and ${\epsilon }^{\mathit{\text{ADR}}}\le$ 1 €/h. In the case of ${\kappa }^{\mathit{\text{ADR}}}$ = 10 ADRs per operator, mean delivery operation costs are reduced by more than 10% for $C^{\mathit{\text{ADR}}}\ge$ 10 receivers per route and ${\epsilon }^{\mathit{\text{ADR}}}\le$ 2 €/h. The biggest cost reduction that can be achieved is −24% for $C^{\mathit{\text{ADR}}}=$ 20 receivers per route and ${\epsilon }^{\mathit{\text{ADR}}}=$ 0.5 €/h with an operator monitoring 10 robots.

The conclusions obtained in Figure 13c are quite equivalent. In the case of ${\kappa }^{\mathit{\text{ADR}}}$ = 5 ADRs per operator, for ${\epsilon }^{\mathit{\text{ADR}}}$ = 0.5 €/h, the two-echelon delivery operations are very robust independently from the operation context because the probability of having lower two-echelon operation costs is above 67% for all values of $C^{\mathit{\text{ADR}}}.$ This means that in more than 67% of the configurations considered in the Monte Carlo analysis, two-echelon delivery costs are lower than BAU ones. On the contrary, for instance, if ${\epsilon }^{\mathit{\text{ADR}}}$ = 4 €/h and $C^{\mathit{\text{ADR}}}$ = 5 receivers per route, the probability of having lower operation costs in the two-echelon scenario is 38%, i.e., in almost 2/3 of cities, implementing this type of ADR will generate higher delivery costs. Adequate use cases will have to be chosen carefully.

If ${\kappa }^{\mathit{\text{ADR}}}$ = 5, ${\epsilon }^{\mathit{\text{ADR}}}$ = 0.5 €/h, and $C^{\mathit{\text{ADR}}}$ = 5 receivers per route, the average cost reduction is 6%, and the probability of having $z_{\mathit{\text{ech}}}^{*}$ lower than $z_{\mathit{\text{BAU}}}$ is 67%. If ${\epsilon }^{\mathit{\text{ADR}}}$ = 4 €/h and $C^{\mathit{\text{ADR}}}$ = 15 receivers per route, the average cost reduction is 4%, but the probability of having $z_{\mathit{\text{ech}}}^{*}$ lower than $z_{\mathit{\text{BAU}}}$ is 58%, i.e., lower than the 67% previously mentioned. It seems better to implement small robots with reduced depreciation costs than larger robots (with a large volume capacity) with higher depreciation costs. If ${\epsilon }^{\mathit{\text{ADR}}}$= 2 €/h and $C^{\mathit{\text{ADR}}}$ = 10 receivers per route, the cost reduction is around 10%, but the probability of having $z_{\mathit{\text{ech}}}^{*}$ lower than $z_{\mathit{\text{BAU}}}$ is 73%. In a majority of cities, the two-echelon scheme generates fewer operation costs, but in a few cities, the two-echelon costs are higher, worsening the average value at a European level. In this configuration, logistics operators willing to implement ADR services should make a reasonable choice of the cities where to implement ADR operations so as not to worsen their operational performance globally.

6. Conclusion and further research

In this article, we have developed a heuristic algorithm to model the total operating costs of BAU and two-echelon last-mile delivery operations in given European cities. In the case of the two-echelon delivery pattern, receivers located in the city center are visited by an ADR fleet operating from an urban logistics micro-hub. Recipients located in the outskirts of the city are visited by conventional delivery vans. Data from different European cities were collected and filtered to fit our modeling methodology. This city sample was used to perform a Monte Carlo analysis of BAU and two-stage total operating costs.

The first main result of this article is the possibility of optimizing two-echelon delivery operations depending on the size of the micro-hub delivery zone and the distance between the micro-hub and the city center. LSPs wishing to implement ADR operations should take these parameters into account when designing urban logistics micro-hubs. Micro-hub delivery zones that are too small or too large, too far or too close to the city center, would lead to increased operating costs.

Then, taking into account the assumptions made in the modeling phase and the definition of the input parameters, the Monte Carlo analysis showed that the use of ADRs with a depreciation cost of 1 €/h and a volume capacity of 15 parcels would reduce the cost of last-mile operations by 13% on a European scale, considering that one operator is able to supervise a fleet of 5 robots. This expected cost reduction on a European scale hides a great deal of diversity for each individual city. In some configurations, two-echelon operations can help reduce last-mile delivery costs by more than 40%. On the contrary, in inappropriate frameworks, the implementation of ADRs would not significantly reduce last-mile costs. In extreme cases, the implementation of ADRs could even increase them. As a consequence, it will be necessary for practitioners to carefully assess ADR deployment scenarios in a given urban context in order to avoid cost increases.

The two main input parameters at stake in explaining the variability of two-echelon delivery costs are the relative weight of infrastructure costs compared to labor costs and demand density. In order to benefit from last-mile cost reductions of more than 40% using ADRs, a service area with the following characteristics should be considered under the previously defined hypotheses:

• An expected demand density in the city greater than 10 parcels/km²/day. This demand density is calculated as the total number of final recipients divided by the total area of the city. A higher cost reduction would be achieved in the case of a high concentration of demand in the city center.

• Personnel costs, including salary and food, of more than 20 euros per hour.

• A daily infrastructure cost per unit area for the urban micro-hub of less than 1 €/m²/day.

The figures given above are orders of magnitude but should give a representative overview of the main tradeoffs at stake.

Finally, it was observed that in order to achieve operating cost reductions of more than 10% (still in the case of one operator supervising 5 robots), an ADR depreciation cost of less than 1 €/h combined with a capacity of more than 15 recipients served per ADR route should be considered.

The main limitations of our model are the following. Firstly, cities whose shape cannot be accurately represented by a disk (coastal cities, for illustration) are not analyzed. The representativeness of our study would be improved by generalizing our model to arbitrary city shapes. Secondly, in the case of large cities, it seems reasonable to assume that several urban micro-hubs would be strategically located in the service region, whereas we consider only one micro-hub in our study. Improving our model by considering the implementation of several micro-hubs is a future line of research.

Finally, we assumed that the values of some operational parameters, such as line-haul and local speeds, expected vehicle volume capacity, stop time per delivery, or expected ADR commercial speed, are constant for all cities. However, these values vary from city to city and from country to country, affecting the cost of last-mile operations. A more complete set of data on last-mile operating conditions in each city would be needed to improve the representativeness of the model.

Acknowledgments

The first author would like to personally acknowledge CARNET Barcelona for the funding of this research article, developed in the framework of his PhD thesis. The participation of the second author has taken place under the METROPOLIS project, PLEC2021-007609, funded by MCIN/AEI/10.13039/501100011033 and by the European Union NextGenerationEU/PRTR. The authors also acknowledge the comments of anonymous reviewers, which greatly helped in improving and clarifying the article.

Disclosure statement

The authors report there are no competing interests to declare.

Additional information

Funding

The first author would like to personally acknowledge CARNET Barcelona for the funding of this research article, developed in the framework of his PhD thesis. The participation of the second author has taken place under the METROPOLIS project, PLEC2021-007609, funded by MCIN/AEI/10.13039/501100011033 and by the European Union NextGenerationEU/PRTR. The authors also acknowledge the comments of anonymous reviewers, which greatly helped in improving and clarifying the article.