Fleet Sizing and Allocation for On-demand Last-Mile Transportation Systems

• A1.

  The delivery fleet consists of at most $M$ vehicles, each with integer capacity $c$;

• A2.

  The set of last-mile stops in each service region is finite;

• A3.

  The set of feasible routes for LMTS vehicles in each service region is finite and preselected based on geometry, historical demand patterns, and some practical constraints—e.g., limits on the maximum number of last-mile stops on a single route or the route’s maximum travel distance or travel time;

• A4.

  The inter-arrival time (headway) between arrival trains is deterministic and equal to $h$.

Given a fleet of at most $M$ vehicles and the sets of last-mile stops and preselected routes for LMTS vehicles, we aim to identify: (1) the number of vehicles to allocate for each service region, (2) a routing plan for the allocated vehicles (i.e., route selection) in each service region, and (3) the assignment of passengers to vehicles for different realizations of demand patterns. Decisions (1)–(2) are first-stage decisions that we make before observing demand patterns. If a route is selected, a vehicle should visit all of the last-mile stops specified on this route. The assignment decisions (3) represent the recourse actions in response to the first-stage decisions and the realizations of demand patterns. The quality of fleet-allocation decisions is a function of the vehicle fixed cost (which may include vehicle rental or purchase cost, etc.) plus a weighted sum of passenger waiting time before boarding a vehicle and in-vehicle riding time.

Let integer decision variable $m_{s}$ represent the number of vehicles allocated to service region $s$. The feasible region $\mathcal{M}$ of variables $m$ is defined in (2) such that the number of allocated vehicles is less than or equal to the total number of vehicles in the fleet:

Let variable $w_{i,k,s}$ represent the number of trips on route $k$ starting right after the arrival of train $i$ at service region $s$, for all $i,k,s$. Let variable $v_{i,s}$ represent the number of vehicles waiting at the metro station after arrival of the train (demand batch) $i$ in service region $s$, for all $i,k,s$. Feasible region $\mathcal{W}$ of $(w,v)$ in (6) defines and constrains the number of vehicles waiting at each metro station $s$ after the arrival of each train $i$ (demand batch) and its corresponding service assignment.

Given a feasible $m\in\mathcal{M}$, $w\in\mathcal{W}$, and a joint realization of uncertain parameters $\xi:=(n_{i,j,s})$, we can compute: (1) the number of passengers $z_{i,j,k,s}$ with a destination of last-mile stop $j$ assigned to route $k$ right after the arrival of train $i$ at service region $s$, and (2) the number of unserved passengers $u_{i,j,s}$ with a destinations of each last-mile stop $j$ waiting at the metro station after the arrival of train (demand batch) $i$ at service region $s$ and its corresponding service assignment using the following linear program (see Table 1 for notation):

The objective function (7a) minimizes a linear cost function of the total waiting time and riding time for passengers. Constraints (7b) and (7c) are passenger flow constraints—i.e., they define and constrain the number of unserved passengers with a destinations of each last-mile stop who are waiting at the metro station after the arrival of train $i_{s}$ and its corresponding service assignment. Constraint (7d) ensures the vehicle service capacity restriction. Finally, constraint (7e) specifies the feasible ranges of the decision variables. Accordingly, we formulate the stochastic fleet sizing and allocation of fleet problem as follows:

The SP formulation in (8) searches for vehicle sizing, allocation, routing, and scheduling decisions that minimize the vehicle fixed cost plus a weighted expected sum of passenger waiting time and riding time, where the expectation is taken with respect to a known joint probability distribution $\mathbb{P}$ of $\xi$, where $\xi$ is the a vector of demand (i.e., a vector of all $(n_{i,j,s}),$ for all $s\in S,i\in I_{s}$ and $j\in J_{s}$). The formulation generalizes the deterministic LMTS routing and scheduling formulation of Wang (2019) by (1) considering multiple service regions, (2) considering fleet sizing and allocation decisions, and (3) incorporating the uncertainty of passengers demand for LMTS. In Proposition 1, we show that the number of trips on routes that have at least two last-mile stops is at most one after each train arrival in the optimal solution.

Proposition 1 implies that, at any vehicle dispatch decision after each train $i$ arrives, to reduce passenger riding time, we would never dispatch multiple identical routes (i.e., $w_{i,k,s}\geq 2$) that visit multiple stops (i.e., $\sum_{j\in J_{s}}\phi_{j,k}\geq 2$). In other words, if a route $k$ visits multiple stops (i.e., $\sum_{j\in J_{s}}\phi_{j,k}\geq 2$), we would dispatch at most one vehicle to serve this route $k$ (i.e., $w_{i,k,s}\leq 1$) after each train $i$ arrives. This proposition will render many integer decision variables $w_{i,k,s}$ to binary decision variables, which can help to improve the computational efficiency of the optimization model.

In addition, we let $\mu:=\mathbb{E}[\xi]$ represent the mean (expected) value of $n$. Then we consider the following mean-support ambiguity set $\mathcal{F}(\mathcal{R},\mu)$:

where $\mathcal{P}(\mathcal{R})$ in $\mathcal{F}(\mathcal{R},\mu)$ represents the set of probability distributions supported on $\mathcal{R}$, and each distribution matches the mean values of $n$. Using the ambiguity set $\mathcal{F}(\mathcal{R},\mu)$, we formulate the fleet sizing and allocation as the following min-max problem:

The formulation (12) seeks to identify vehicle sizing, allocation, routing, and scheduling decisions that minimize the worst-case expected cost of passenger waiting time and riding time over a family of distributions of random parameters residing in the ambiguity set $\mathcal{F}(\mathcal{R},\mu)$ for the demand. $Q(m,w,\xi)$ is the recourse problem defined in (7).

Note that the recourse problem $Q(m,w,\xi)$ is a minimization problem. Thus, in (14) we have an inner max-min problem. Next, we use $Q(m,w,\xi)$ properties to derive an equivalent single minimization reformulation of (14). For fixed $m\in\mathcal{M}$, $w\in\mathcal{W}$, and a realized value of $\xi=n$, $Q(m,w,\xi)$ is a linear program. We formulate $Q(m,w,\xi)$ in its dual form as

where $\Gamma$ and $\psi$ are the dual variables associated with constraints (7b)–(7c) and (7d), respectively. Given the dual formulation of $Q(m,w,\xi)$ in (17) and a fixed and feasible ($\boldsymbol{\rho,w}$), we can rewrite the inner maximization problem $\max\limits_{n\in\mathcal{R}}\{Q(m,w,\xi)+\sum\limits_{s\in S}\sum\limits_{i\in I_{s}}\sum\limits_{j\in J_{s}}-n_{i,j,s}\rho_{i,j,s}\}$ in (14) as follows

As we show in the proof of Proposition 3 in A, for fixed $(m,w,\xi)$ and $\rho$, problem (18) is equivalent to the minimization problem in (19) .

Combining the inner problem in the form of (19) with the outer minimization problems in (14) and (12), we derive the following equivalent MILP reformulation of the DR model in (12) (see A)

To generate demand profile for the instance 1-4, we use a similar magnitude of the number of passengers in the LMTS literature (e.g., Wang and Odoni (2016); Wang (2019)), as well as the same random parameter generation procedures in the distributionally robust scheduling and optimization literature (e.g., Jiang et al. (2017); Shehadeh (2021); Mak et al. (2014)). For instance 1–4, we randomly generate the mean values $\mu$ of $n$ from a uniform distribution ${U[1,4]}$ and standard deviation $\sigma=0.5\mu$. We randomly generate $N$ in-sample realizations $(n_{i,j,s}^{\mbox{\tiny n}})$ by following lognormal (LogN) distributions with the generated $\mu_{i,j,s}$ and $\sigma_{i,j}$, for all $s\in S,i\in I_{s},j\in J_{s},n\in N$. LogN is a standard distribution of model customers’ demand and service times in a wide range of applications. Kamath and Pakkala (2002) results suggest that the LogN is a suitable distribution for modeling stochastic demands in an economic context. For the NYC instance, in Table 2 in C, we present the empirical $\mu$ and $\sigma$ of batch demand $n_{i,j,s}$.

According to Gomez-Ibanez et al. (1999), for work trips in San Francisco, the monetary value of a unit of transfer waiting time is 195$\%$ of the user’s after-tax wages, and the monetary value of a unit of in-vehicle riding time is 76$\%$ of the user’s after-tax wages. In general, we should have $\beta^{\mbox{\tiny w}}>\beta^{\mbox{\tiny r}}$ in the objective function. Following similar parameter selections as in Wang (2019), we normalize $\beta^{\mbox{\tiny r}}=1$ and $\beta^{\mbox{\tiny w}}\in[2,3]$ (we perform sensitivity analysis in Section 4.4). As for the vehicle fixed cost $f_{s}$, if the service provider already has a fleet of vehicles, $f_{s}$ could be very small; if the service provider rents vehicles, $f_{s}$ should include the rental fee and operating cost (e.g., fuel and gas); if the service provider needs to purchase a new fleet of vehicles, $f_{s}$ may be very large and a complex depreciation should be considered. Unless stated otherwise, we test two values of $f_{s}$: (1) $f_{s}=0$ (ignoring the fixed cost and assuming an existing vehicle fleet), and (2) $f_{s}=4,000,$ from a range of [2,000, 6,000] (renting vehicles to serve passengers who are less sensitive to riding and waiting times)²²2The range is approximated considering the general after-tax wage for city residents, the rental price of vehicles with capacity 4-10, vehicle price and possible depreciation, fuel cost, etc. (considering vehicle rental fee and/or depreciation and operating cost; see Appendix D). We investigate the impact of $f_{s}$ in Section 4.4.

As in prior applied distributionally robust literature (see, e.g., Jiang et al. (2017), Shehadeh and Sanci (2021), Shehadeh and Padman (2021), and references therein), we respectively use the $20\%$-quantile and $80\%$-quantile values of the $N$ in-sample data to approximate the lower $\underline{n}$ and upper $\overline{n}$ bounds of $n$. We optimize the SP model by using all of the $N$ data points, and the DR model with the corresponding mean, lower bounds, and upper bounds. We implemented the models using AMPL2016 programming language calling CPLEX V12.6.2 as a solver with default settings. We ran all experiments on a computer with an Intel Core i7 processor, 2.5 GHz CPU, and 16 GB (1600MHz DDR3) of memory, and imposed a solver time limit of 1 hour.

Our choice of the sample size $N$ to solve the SAA of SP was motivated by the trade-off between the computational effort required to solve the resulting mixed-integer linear programs (MILPs) and the quality of approximation of the expected value objective of SP by its SAA approximation. On the one hand, the sizes of MILP instances increase with $N$ and their solution times also increase. On the other hand, optimal solutions of SAA instances with larger values of $N$ are likely to be closer to optimality compared with the expected value objective.

The literature on the SAA method provides theoretical insights and guidance for selecting a sample size from this perspective. We implemented the so-called Monte Carlo Optimization (MCO) procedure to compute statistical lower and upper bounds on the optimal value of SP based on an optimal solution to its SAA approximation, which in turn provides a statistical estimate of the relative approximation gap between the optimal value of SP and its SAA approximation (see, e.g., Homem-de Mello and Bayraksan (2014) and Linderoth et al. (2006) for a thorough discussion of MCO). Applying the MCO procedure to our SP model with $N=100$, we estimate the relative approximation gaps for the SP instances described in Table 2 to range between 1% and 5%. In contrast, larger sample sizes resulted in longer solution times without consistent and significant improvements in the relative approximation gaps. Based on these considerations, we selected $N=100$ for our computational experiments.

Table 3 presents solution times (in seconds) using the SP and DR models for different values of $M$ (maximum number of vehicles). We first observe that the DR can quickly solve all instances under the two ranges of the average number of passengers and all values of $f_{s}$ and $M$ and significantly faster than the SP model. In contrast, the SP fails to solve all of the SAA instances corresponding to instance 4 to optimality within the time limit, and terminates with either a large relative MIP (relMIP) gap (relMip:=$\frac{UB-LB}{LB}$, where UB is the best upper bound and LB is the linear programming relaxation-based lower bound obtained at termination) or without any feasible MIP solution, and thus no upper bound. Additionally, the SP’s solution times differ under the two ranges of the average number of passengers and values of $f_{s}$ and $M$.

Under $f_{s}=0$, the two models allocate all vehicles (i.e., $\sum_{s\in S}m_{s}^{*}=M)$. The SP’s solution times decrease as $M$ increases from 40 to 80. Consider instance 1, for example: SP’s solution times decrease from 197 sec and 2,477 sec to 31 sec and 230 sec, respectively, under uncertain demand [1,4] and [3,7]. A possible explanation for this is that when $M$ is large, we can satisfy passenger demand with any vehicle allocations, routing, and scheduling decisions. That is, we can allocate a larger number of vehicles in each region that may be sufficient to transport passengers to their last-mile stops via the direct routes to those steps (e.g., each route serves only one stop). When $M$ is small, vehicle allocations, routing, and scheduling decisions are subtle and difficult to optimize.

Under $f_{s}=4,000$, the two models allocate a subset of the $M$ vehicles to minimize the total cost function (see Table 4). The larger SP solution times indicate that SP’s routing and scheduling decisions are subtle and difficult to optimize in this case.