# Last-Mile Shuttle Planning for Improving Bus Commuters’ Travel Time Reliability: A Case Study of Beijing

Authors: Weibin Kou; Jiayu Wang; Yanxi Liu; Chenxu Li
Journal: Journal of Advanced Transportation
 (2022)

---

## Abstract

This study proposes a multiobjective mixed integer nonlinear programming model for a last-mile shuttle service to improve bus commuters’ travel time reliability. The approach aims to assign the routes that pick up the transit passengers located at the different stops by shuttle service. A bilevel optimization model is established: the upper model of route design considers the tradeoff between time cost and fare cost when some of the passengers take the shuttles, and the lower model assigns the demand of transit passengers. The proposed model effectively captures the reliability of travel time because related parameters are estimated by a statistical fitting test with a large number of real-world bus geographic information system (GPS) data. Moreover, dynamic demand diverting from conventional transit to shuttle service and travel time reliability, including passenger in-vehicle time (IVT) and waiting time (WT), are fully considered in this model. Since the task is a nonlinear programming model, a two-stage algorithm combined with linearization processing is presented to find an optimal solution. Finally, from the case study of Zhongguancun Software Park zone in Beijing, it is indicated that when last-mile shuttle service is provided, bus passengers’ travel time reliability of last-mile trips can be improved by 14%. The study can be an important reference for improving the low reliability widely existing in the current transit commuters’ last-mile problem.


---

## Body

## 1. Introduction

In China, traffic congestions and delays are becoming more and more frequent and severe with the rapid increase of automobile ownership. For metropolises like Beijing, the development of public transportation systems is given priority to relieve congestions. However, due to the limitation of urban land resources and heavily congested traffic, poor travel time reliability has become the main restriction for improving the service quality of the public transportation system. In downtown areas and residential communities, a large amount of gathered commuting demand in peak hours will cause heavy congestion and additional delay of in-vehicle time (IVT) and waiting time (WT) [1, 2].The area for last-mile trip typically covers narrow roads and suffers frequent traffic congestion. Therefore, by conventional transit systems with fixed routes and stops, passengers cannot avoid passing through those congested areas, which results in great difficulty to get to a destination on time [3]. Providing reliable and efficient service is essential in transit planning and operation to improve passenger satisfaction.This work proposes a robust optimization model for shuttle service considering travel time reliability by designing shuttle routes and estimating the corresponding number of users. It is organized as follows: in Section2, a literature review is conducted, including traditional qualification methods of travel time reliability and how to improve reliability by the flexible transit design. In Section 3, an offline method for estimating travel time reliability is proposed and a bilevel programming formulation is established. In Section 4, a case study of Beijing is discussed. And lastly, the conclusions are discussed in Section 5.


## 2. Literature Review

The review of existing researches mainly focuses on the following two domains: how to estimate transit service reliability and its influences on operation and how to design the shuttle route and operation schedule to improve reliability.


### 2.1. Estimation of Transit Reliability

To improve the unreliable last-mile trip of the transit system, micromobility, including flexible light shuttle, emerges [4]. For transit network design and operation management, reliability is considered as one of the key factors, which includes transit-oriented indicators (e.g., the percentage of on-time buses, the percentage of regular buses, the headway regularity) and passenger-oriented indicators (e.g., the WT at stops, the cleanliness, the space onboard, and the percentage of passengers receiving a punctual service). The studies show that service reliability is twice as important as frequency and almost seven times more important than information for passengers [5].Estimation of transit service reliability includes the regularity and irregularity analysis of passenger WT, IVT, and walking time. Existing researchers have established quantitative volatility indicators to model and analyze the reliability of passenger travel time [6]. Yu et al. [7] used different methods by support vector machine (SVM), artificial neural network (ANN), k nearest neighbours algorithm (k-NN), and linear regression (LR). Salvo et al. [8] believed that the scheduling reliability of the public transport service was crucial to increase attractiveness against private car use and proposed the estimation method of bus speed for various traffic conditions. Li and Quadrifoglio [9] conducted Lempel–Ziv algorithm to quantify passenger travel time reliability. Elefteriadou and Cui [10] proposed a comprehensive framework under different scenarios to analyze reliability and made application to the system performance in Florida. Li et al. [11] examined reliability-based performance indicators for bus routes using global position system (GPS) bus trajectory data.With the development of automotive facilities, multisource data of automatic vehicle location (AVL) data, geographic information system (GIS) data, and GPS data were widely used to investigate the reliability from the passengers’ and companies’ perspectives. Leod [12] contributed to the general theory of estimating headway variance using incomplete data, focusing on the missing buses or discarding spurious bus headways. Chen et al. [13] analyzed bus service reliability at the stop, route, and network level. The service reliability was defined as the possibility that a bus would adhere to the headway between successive buses at each stop within a given time period. Lin et al. [14] developed running time adherence and headway regularity to reflect bus service reliability. Running time adherence (measured in %) was defined as the average difference between the actual and the scheduled running times and headways relative to the scheduled running time. Ömer et al. [15] analyzed the selection problem of electric vehicles charging stations and proposed the approach for a high degree of uncertainty problems. Barabino et al. proposed an offline framework for the diagnosis of time reliability by analyzing AVL data. The bus stops and the time periods in which reliability was insufficient were accurately characterized [16]. Regularity over all bus stops and time periods was measured [17]. He also used the two conceptions of “fraction of passengers waiting less than the scheduled headway plus a possible irregularity threshold” and “the fraction of passengers who will see a bus serving their stop within an acceptably short interval after they arrive” to reflect transit travel time reliability. Both bus and passenger arrival and departure data were collected and processed [18]. In addition, to explore the more accurate qualification results, the effect of AVL anomalies on headways and schedule deviations and the influence of anomalies on the transit service reliability was analyzed [19].


### 2.2. Optimization of Shuttle Design with Reliability

A shuttle bus can be considered as a new type of public transportation service system that is different from the traditional point-fixed and line-fixed operation mode. Commuting shuttle has several similar features with limited buses (e.g., customized bus and feeder bus). It is regarded as the complementary part of the normal conventional transit system, and it serves specific routes during specific hours for specific types of passengers [20]. Compared with conventional bus, the shuttle bus has characteristics of flexible routes and stop selections, diversified service modes, expensive fares, and individual services. It can be used to solve the last-mile trip as an effective measure to alleviate peak travel demand in hotspots [21]. It has wide applicability in many areas, such as hospital emergency vehicle path planning [22], disabled public transportation service planning [23], and regional community school bus planning [24].In most previous transit design optimizations, mathematical programming and simulation-based approaches were used, in which the averages of travel time were considered [25]. However, a big variance in the travel time could lead to poor reliability in actual conditions. Therefore, for transit design, including shuttles, reliability should be considered accurately. The uncertainty of transport networks was applied for transportation network design [26]. Yao et al. [27] presented a robust optimization model considering travel time uncertainty to satisfy the demand of passengers and provide reliable transit service. Markovic et al. [28] analyzed multimodal public transport journeys including proper demand transfer strategy, to improve the whole trip’s travel time reliability. Bao et al. [29] analyzed the optimization model of airport shuttle that bus routes were constructed by taking operational reliability maximization as the main goal.The review of the literature reveals the following critical issues, which deserves further investigation:(1)
Although serval studies have proposed various integrations of shuttle design, few of them consider the reliability of transit passengers. In recent years, with increasing attention to reliability, providing reliable transit service for the last-mile trip has become more important.(2)
Traditional shuttle bus design only considers the fixed shuttle demand. When conventional transit system becomes unreliable during peak hours, some passengers may prefer to choose the new shuttle service. Quantitively, analysis of the dynamic demand should not be neglected, and integration of conventional and shuttle bus services should be analyzed.(3)
Previous studies analyze travel time reliability by different online and offline methods, but most of them are conducted by numerical experiment or case study in small cities. In big cities like Beijing, during peak hours, the traffic is heavily congested and the amount of boarding and alighting passengers at “hot stop” will rapidly increase, which will cause the transit system to become extremely unreliable. Therefore, field data in big cities should be presented to help design more reliable shuttle services.


## 2.1. Estimation of Transit Reliability

To improve the unreliable last-mile trip of the transit system, micromobility, including flexible light shuttle, emerges [4]. For transit network design and operation management, reliability is considered as one of the key factors, which includes transit-oriented indicators (e.g., the percentage of on-time buses, the percentage of regular buses, the headway regularity) and passenger-oriented indicators (e.g., the WT at stops, the cleanliness, the space onboard, and the percentage of passengers receiving a punctual service). The studies show that service reliability is twice as important as frequency and almost seven times more important than information for passengers [5].Estimation of transit service reliability includes the regularity and irregularity analysis of passenger WT, IVT, and walking time. Existing researchers have established quantitative volatility indicators to model and analyze the reliability of passenger travel time [6]. Yu et al. [7] used different methods by support vector machine (SVM), artificial neural network (ANN), k nearest neighbours algorithm (k-NN), and linear regression (LR). Salvo et al. [8] believed that the scheduling reliability of the public transport service was crucial to increase attractiveness against private car use and proposed the estimation method of bus speed for various traffic conditions. Li and Quadrifoglio [9] conducted Lempel–Ziv algorithm to quantify passenger travel time reliability. Elefteriadou and Cui [10] proposed a comprehensive framework under different scenarios to analyze reliability and made application to the system performance in Florida. Li et al. [11] examined reliability-based performance indicators for bus routes using global position system (GPS) bus trajectory data.With the development of automotive facilities, multisource data of automatic vehicle location (AVL) data, geographic information system (GIS) data, and GPS data were widely used to investigate the reliability from the passengers’ and companies’ perspectives. Leod [12] contributed to the general theory of estimating headway variance using incomplete data, focusing on the missing buses or discarding spurious bus headways. Chen et al. [13] analyzed bus service reliability at the stop, route, and network level. The service reliability was defined as the possibility that a bus would adhere to the headway between successive buses at each stop within a given time period. Lin et al. [14] developed running time adherence and headway regularity to reflect bus service reliability. Running time adherence (measured in %) was defined as the average difference between the actual and the scheduled running times and headways relative to the scheduled running time. Ömer et al. [15] analyzed the selection problem of electric vehicles charging stations and proposed the approach for a high degree of uncertainty problems. Barabino et al. proposed an offline framework for the diagnosis of time reliability by analyzing AVL data. The bus stops and the time periods in which reliability was insufficient were accurately characterized [16]. Regularity over all bus stops and time periods was measured [17]. He also used the two conceptions of “fraction of passengers waiting less than the scheduled headway plus a possible irregularity threshold” and “the fraction of passengers who will see a bus serving their stop within an acceptably short interval after they arrive” to reflect transit travel time reliability. Both bus and passenger arrival and departure data were collected and processed [18]. In addition, to explore the more accurate qualification results, the effect of AVL anomalies on headways and schedule deviations and the influence of anomalies on the transit service reliability was analyzed [19].


## 2.2. Optimization of Shuttle Design with Reliability

A shuttle bus can be considered as a new type of public transportation service system that is different from the traditional point-fixed and line-fixed operation mode. Commuting shuttle has several similar features with limited buses (e.g., customized bus and feeder bus). It is regarded as the complementary part of the normal conventional transit system, and it serves specific routes during specific hours for specific types of passengers [20]. Compared with conventional bus, the shuttle bus has characteristics of flexible routes and stop selections, diversified service modes, expensive fares, and individual services. It can be used to solve the last-mile trip as an effective measure to alleviate peak travel demand in hotspots [21]. It has wide applicability in many areas, such as hospital emergency vehicle path planning [22], disabled public transportation service planning [23], and regional community school bus planning [24].In most previous transit design optimizations, mathematical programming and simulation-based approaches were used, in which the averages of travel time were considered [25]. However, a big variance in the travel time could lead to poor reliability in actual conditions. Therefore, for transit design, including shuttles, reliability should be considered accurately. The uncertainty of transport networks was applied for transportation network design [26]. Yao et al. [27] presented a robust optimization model considering travel time uncertainty to satisfy the demand of passengers and provide reliable transit service. Markovic et al. [28] analyzed multimodal public transport journeys including proper demand transfer strategy, to improve the whole trip’s travel time reliability. Bao et al. [29] analyzed the optimization model of airport shuttle that bus routes were constructed by taking operational reliability maximization as the main goal.The review of the literature reveals the following critical issues, which deserves further investigation:(1)
Although serval studies have proposed various integrations of shuttle design, few of them consider the reliability of transit passengers. In recent years, with increasing attention to reliability, providing reliable transit service for the last-mile trip has become more important.(2)
Traditional shuttle bus design only considers the fixed shuttle demand. When conventional transit system becomes unreliable during peak hours, some passengers may prefer to choose the new shuttle service. Quantitively, analysis of the dynamic demand should not be neglected, and integration of conventional and shuttle bus services should be analyzed.(3)
Previous studies analyze travel time reliability by different online and offline methods, but most of them are conducted by numerical experiment or case study in small cities. In big cities like Beijing, during peak hours, the traffic is heavily congested and the amount of boarding and alighting passengers at “hot stop” will rapidly increase, which will cause the transit system to become extremely unreliable. Therefore, field data in big cities should be presented to help design more reliable shuttle services.


## 3. Methodology

### 3.1. Descriptions

This paper proposes a framework to characterize passenger travel time reliability and design shuttle bus considering reliability quantitively. The overall flowchart is presented in Figure1. Firstly, the formulation of passengers’ travel time reliability is discussed and an offline framework to estimate time reliability is proposed. Then, a bilevel optimization formulation is constructed including model establishment and algorithm design. Finally, a case study in Beijing is conducted to validate the model.Figure 1
Flowchart of the method framework.The optimization model is designed not only to give the route plan of the shuttle service but also can estimate how many passengers should be diverted from the original conventional lines to the new shuttles in order to relieve the “hot” OD with poor reliability to an acceptable level. The model is passenger-oriented and passengers’ travel time, reliability, and fare cost are considered. The last-mile shuttle service is formulated as a bilevel programming model. The upper-level model is a mathematical formulation for programming bus routes and corresponding demand. The lower-level model is to conduct the transit assignment for the conventional network. By the demand assignment result from the lower model, passenger travel cost can be estimated.On the basis of the previous studies [30], the following assumptions are made in the paper:(1)
Passengers arrive at the bus stop randomly(2)
Passengers know the mean and variance of IVT and make the selection that has the minimum travel cost(3)
The values of IVT and WT are assumed to be independent and constant(4)
Only cross-platform transfer is considered in this study


### 3.2. Offline Method for Estimating Passenger Travel Time Reliability

With the continuous updating of data acquisition equipment in recent years, transit passenger travel behavior can be described more accurately by mobile GPS data, which can record the bus trajectory by the second level. The method framework includes error data preprocessing, formulation of IVT and WT, and fitting analysis.Table 1
Sample of GPS data.
Line numberDirectionVehicle numberVehicle tripCollection timeLongitudeLatitudeSpeed870J3797120190114060110116.33539.93619.80870J3797120190114060111116.39139.87913.50870J3797120190114060112116.34739.91516.40870J3797120190114060113116.38539.88517.80870J3797120190114060114116.39639.72318.60870J3797120190114060115116.40539.71519.50……………………


#### 3.2.1. Data Processing Based on GPS Data

The GPS data are used to estimate the parameters of bus running conditions, which are collected from on-vehicle GPS facilities. The original data sets mainly have two problems: data losses and data errors. Considering existing data handling methods [12, 19], the data quality controls are conducted by the following three steps:Step 1: data processing of original data. The available fields are collected from the database, including line number, direction, vehicle number, vehicle trip, collection time, longitude, latitude, and speed, which are presented in Table1.Step 2: the fault data of speed, longitude, and latitude are filled and corrected by linear interpolation, which is described as follows:(1)xtj=tj−tixtk+tk−tjxtitk−ti,wherextj is the corrected data of time tj, (ti, tk) is the neighbour time points of tj, xtk and xti are the original values of ti and tk, respectively.Step 3: the GPS data are matched with a GIS map to analyze the spatial locations of passenger travel trajectory, as Figure2 shows. The second-by-second GPS data are aggregated to calculate IVT and WT for different bus trips and lines. In addition, the matching process is verified based on the station number information from passenger IC card data.Figure 2
GPS data matching with GIS.


#### 3.2.2. Formulation of IVT Reliability

When bus passengers plan a trip, they will add the “buffer time” based on the expected travel time, considering the potential fluctuations. Although they may arrive at the destination early, early arrival is better since the negative impact of the late arrival is large [31]. As a result, most passengers may choose to plan the extra buffer time during peak hours. The buffer time can be the quantitative index to measure passenger travel time reliability [32]. The more buffer time passengers need to plan, the worse travel time reliability of the network is.As is depicted in Figure3, IVT considering reliability can be computed as follows:(2)tij=tj−ti,tijin−vehicle=Etij+ρσtij,Etij=∑n=1NtijnN,σtij=1N−1∑n=1Ntijn−Etij2,where tj and ti are samples of bus arrival time at stop i and j, which can be extracted from the GPS data. tij is sample time from node i to node j. tijin−vehicle is IVT considering reliability from node i to node j. E (tij) is the expectation of IVT from node i to node j for all the samples. σtij is the buffer time to measure the reliability, i.e., the standard deviation(SD) of the link ij.Figure 3
Illustration of IVT.ρ is reliability preference parameter. When ρ = 1, it means that the passenger will plan the equivalent buffer time as expected travel time. Passenger is “risk-neutral for on-time arrival”. When ρ = 0, the passenger does not consider the influence of travel time reliability. When ρ > 1, the passenger is regarded as conservative, in which the passenger does not trust bus service reliability and will plan more extra buffer time compared with the risk-neutral passenger. When ρ < 1, the passenger is understood as “risky for on-time arrival” [30].N is the number of samples on link (i, j).


#### 3.2.3. Formulation of WT Reliability

As is presented in Figure4, assuming that passenger arrives at the bus stop on time t, and the following bus arrives on time h, WT can be calculated as follows:(3)w=h−t,twaiting=Ew+ρσw,Ew=Eh−t=Eh−Et,σw=Dw=Dh+Dt,where w is sample WT, t is passenger arrival time at the bus stop, and h is the following bus arrival interval. twaiting is WT considering reliability. E (w) is the expectation of WT for all the samples. σw is the buffer time to measure the reliability, such as SD, and D(w) is the deviation of sample WT.Figure 4
Illustration of WT.Both arrivals of passengers and bus at the stop are random, which can influence WT. Usually, each bus trip does not have a strict arrival timetable on stop level in Beijing. Passenger arrival time and bus arrival time can be considered as an independent event.Studies usually focus on passenger arrival behavior within one fixed bus arrival interval. However, with the increasing influences of traffic congestion and other uncertain factors, the variance of bus arrival time becomes larger, which should not be neglected. Bus intervalh can be regarded following a certain probability distribution h∼h′x. Based on field GPS data, statical characteristics of E(h) are estimated. The passenger arrival rate is regraded following a Poisson distribution on [0, H] [33], which is written as follows:(4)PT=t=zwHte−zwHt!,Et=∫0HtPTdT=zwH,σt=zw,where zw is the ratio of the expected passenger arrival time/bus arrival interval.The model works for high-frequency service. Considering the condition that passengers arrive randomly and the individual WT are different, the statical fitting analysis by a large amount of field data is presented in the case study. Root mean square error (RMSE) and R-square are used to validate the fitting result, which are presented in section case study.


### 3.3. Bilevel Programming Formulation

#### 3.3.1. Formulation of General Cost

The general costs are combined with IVT cost, WT cost, extra waiting delay cost, transfer cost, and a fare cost. For each arc that passenger passes, the formulation of general cost can be calculated as follows:(1)
IVT cost(5)C=utsin−vehicle,∀s∈Sin−vehicle,whereC is general cost, u is unit travel time cost (yuan/h), s is travel arc, and Sin−vehicle is set of in-vehicle travel arc.(2)
WT cost(6)C=utswaiting,∀s∈Swaiting,whereSwaiting is set of waiting arc.(3)
Extra waiting delay time costDuring peak hours, it is a common phenomenon that there is a large number of passengers boarding and alighting at bus stops, which will cause the speed of boarding and alighting to slow down [34]. As a result, there will be an extra delay which is increased with the rising number of passengers. The extra delay is related to the number of existing passengers at the stop. In this paper, the extra delay is determined as follows:(7)Ds=βsYsKn,∀s∈Swaiting,C=uDs,∀s∈Swaiting,Ds is extra delay for waiting arc s. Ys is volume of waiting passengers. K is rated carrying capacity. βs and n are estimated parameters.(4)
Transfer time costUnder the assumption that most passengers will choose to transfer at the same stop, the passenger transfer time is regarded as one additional WT. Therefore, the transfer time cost is expressed as follows:(8)C=utswaiting,∀s∈Stransfer,whereStransfer is set of transfer arc.(5)
Fare cost(9)C=u′ls,∀s∈Sin−vehicle,u′ is unit ticket cost (yuan/km). ls is travel distances (km).


#### 3.3.2. Upper-Level Programming Model of Shuttle Route Design

Last-mile shuttle bus is designed as a measure to relieve congested demand and provide more reliable transit service for passengers. Therefore, the maximum improvement benefit is regarded as the objective function. For the diverted passengers, they can take the more flexible and reliable shuttle bus, but they should burden the additional fare cost. The model is formulated by the following optimization model.(10)maxC=∑o∑dCTod¯yod−∑o∑d∑i∑jutijin−vehicleyodxijod−∑o∑dcsodyod−∑o∑d∑i∑ju′lijyodxijod.(11)s.t.∑i=oxijod=1,∀o,d∈V,∀i,j∈S,(12)∑j=dxijod=1,∀o,d∈V,∀i,j∈S,(13)∑jxijod=∑jxjiod,∀o,d∈V,∀i,j∈S,(14)qod=dod+yod,∀o,d∈V,(15)kod≥k1,∀o,d∈V,(16)yod≤Kfl,∀o,d∈V,(17)xijod∈0,1,∀o,d∈V,∀i,j∈S,(18)dod≥0,∀o,d∈V,(19)yod≥0,∀o,d∈V,where CTod¯ is the average time cost (yuan) and csod is the average ticket cost of conventional bus (yuan/km). The decision variables are xijod and yod, in which xijod means that whether the shuttle from node O to node D passes through road arc (i, j) and yod reflects passengers that take the shuttles.∑o∑dCTodyod¯ and ∑o∑dcsodyod reflect total time cost and ticket cost when passengers take conventional transit. ∑o∑d∑i∑jutijin−vehicleyodxijod and ∑o∑d∑i∑ju′lijyodxijod are total time cost and ticket cost for shuttle passengers. The objective is to maximize passenger benefit when diverted passengers choose the new shuttle service.Equations (11)–(13) are the flow conservation constraints for each node, in which V is set of a bus stop.Equation (14) means that for each OD pair, the total OD demand should equal the sum of shuttle demand and the conventional bus demand, in which qod reflects overall demand and dod is the remaining demand for conventional transit.Equation (15) reflects the reliability threshold, which means that only when kod is larger than the certain level k1, the shuttles will begin to run. Kod quantifies the impact of time reliability, which is the key threshold for judging whether it is necessary to provide shuttle service. The larger the kod is, the worse the reliability condition is. kod is calculated as follows:(20)kod=CRodCTod,∀o,d∈V,CRod=∑suρσs+Dsysod,∀o,d∈V,∀s∈Sin−vehicle∪​Swaiting,CTod=∑s∈Sin−vehicleutsin−vehicleysod+∑s∈Swaitingutswaiting+Dsysod,∀o,d∈V,where CRod reflects the delay time cost which is produced by the unreliability of travel time and CTod represents the total time cost from O to D. ysod is demand volume on arc s.Equation (16) reflects the vehicle carrying capacity constraint, in which fl is shuttle service frequency of line l.Equations (17)–(19) give the value range for the variables.


#### 3.3.3. Lower-Level Programming Model of Demand Assignment

Lower-level model aims to assign the existing bus demand, estimate travel time and cost on the network arc, and quantify the reliability condition. In the lower model, bus passenger IVT cost, WT cost, extra delay cost, transfer time cost, and fare cost are considered.Based on the mentioned analyses, a lower-level flow assignment model based on user equilibrium is formulated as follows:(21)minC=∑o∑d∑p∫0YsCsYsdYs.(22)s.t.∑p∈Pyodp=qod,∀o∈V,∀d∈V,(23)yodp≥0,∀o∈V,∀d∈V,∀p∈P,where CsYs is the general cost function of arc s. Ys is the total volume of passengers on arc s, which is derived that Ys=∑i∑j∑pyodpδodps. yodp is the volume of passengers from i to j on available path p. δodps is the correlation parameter that if travel arc s is on path p from o to d, δodps  = 1, otherwise, δodps  = 0. P is a set of available paths.The objective function is the sum of the results of integrating the general cost functions. The model has two constraints. For (22), it means that the sum of passengers on all paths is equal to the volume of the demand for each OD pair. Constraint (23) means that the decision variable is nonnegative. The optimal solution of the objective function is equivalent to user equilibrium.


### 3.4. Solution Algorithm

Each of the upper-level solutions represents a scheme of shuttle route planning. For each scheme, demand should be assigned by the lower-level model to estimate whether overall reliability is improved to the acceptable threshold (k1). Then, the result of lower-level will be substituted into the upper level to calculate the objective function. In addition, considering the upper model is a nonlinear programming formulation, the computational complexity will be higher, which can be solved by a linearization algorithm.Based on the above analysis, a two-stage algorithm is proposed, which is depicted in Figure5. The first step is for demand assignment, which estimates kod of existing conventional network and selects unreliable OD, which can be diverted to shuttle service. The second step is the shuttle route optimization model including a linearization algorithm to solve the 0-1 nonlinearization programming model. The optimal shuttle service plan will be calculated finally. The algorithm is described in detail as follows:Stage 1: conduct the lower-level demand assignment model based on the existing local bus network and demand. Estimatekod among all OD in the local bus network. Select the congested OD by reliability threshold k1: if kod >k1, goes for stage 2. Otherwise, the selection process ends.Stage 2: conduct the upper-level model. Linearization transforming is introduced for the model. The shuttle route and diverted demand are output finally by the model. The linearization algorithm is described in detail.Figure 5
Overall solution framework.For upper model, the two decision variables are nonnegative integer variables and 0-1 variables. As equation (10) shows, yodxijod is the nonlinear structure. Auxiliary decision variable Mijod (Mijod≥0) and real number M (M≥0) are introduced to construct auxiliary inequations. Auxiliary constraints are introduced as the following equations show:(24)yod−M1−xijod≤Mijod≤yod,∀o,d∈V,∀i,j∈S,(25)Mijod≤Mxijod,∀o,d∈V,∀i,j∈S.Whenxijod=1， yod≤Mijod≤yod. Therefore, Mijod=yodxijod. When xijod=0, according to inequality (25), Mijod≤0. Therefore Mijod=0. There are following inequalities that Mijod≥0 and yod≥0. Within the range of yod and xijod, Mijod always equals yodxijod. As a result, the nonlinear expression yodxijod is transformed into linear expressions with a series of inequalities.The solving difficulty is significantly reduced. In the paper, the commercial linear programming software Gurobi is used to solve the model. After the linearization, the model can be expressed as follows:(26)maxC=∑o∑dCTod¯yod−∑o∑d∑i∑jutijin vehicleMijod−∑o∑dcsodyod−∑o∑d∑i∑ju′lijMijod,s.t.∑i=oxijod=1,∀o,d∈V,∀i,j∈S,∑j=dxijod=1,∀o,d∈V,∀i,j∈S,∑jxijod=∑jxjiod,∀o,d∈V,∀i,j∈S,qod=dod+yod,∀o,d∈V,kod≥k1,∀o,d∈V,yod≤kfl,∀o,d∈V,yod−M1−xijod≤Mijod≤yod,∀o,d∈V,∀i,j∈S,Mijod≤Mxijod,∀o,d∈V,∀i,j∈S,xijod∈0,1,∀o,d∈V,∀i,j∈S,dod≥0,∀o,d∈V,yod≥0,∀o,d∈V.


## 3.1. Descriptions

This paper proposes a framework to characterize passenger travel time reliability and design shuttle bus considering reliability quantitively. The overall flowchart is presented in Figure1. Firstly, the formulation of passengers’ travel time reliability is discussed and an offline framework to estimate time reliability is proposed. Then, a bilevel optimization formulation is constructed including model establishment and algorithm design. Finally, a case study in Beijing is conducted to validate the model.Figure 1
Flowchart of the method framework.The optimization model is designed not only to give the route plan of the shuttle service but also can estimate how many passengers should be diverted from the original conventional lines to the new shuttles in order to relieve the “hot” OD with poor reliability to an acceptable level. The model is passenger-oriented and passengers’ travel time, reliability, and fare cost are considered. The last-mile shuttle service is formulated as a bilevel programming model. The upper-level model is a mathematical formulation for programming bus routes and corresponding demand. The lower-level model is to conduct the transit assignment for the conventional network. By the demand assignment result from the lower model, passenger travel cost can be estimated.On the basis of the previous studies [30], the following assumptions are made in the paper:(1)
Passengers arrive at the bus stop randomly(2)
Passengers know the mean and variance of IVT and make the selection that has the minimum travel cost(3)
The values of IVT and WT are assumed to be independent and constant(4)
Only cross-platform transfer is considered in this study


## 3.2. Offline Method for Estimating Passenger Travel Time Reliability

With the continuous updating of data acquisition equipment in recent years, transit passenger travel behavior can be described more accurately by mobile GPS data, which can record the bus trajectory by the second level. The method framework includes error data preprocessing, formulation of IVT and WT, and fitting analysis.Table 1
Sample of GPS data.
Line numberDirectionVehicle numberVehicle tripCollection timeLongitudeLatitudeSpeed870J3797120190114060110116.33539.93619.80870J3797120190114060111116.39139.87913.50870J3797120190114060112116.34739.91516.40870J3797120190114060113116.38539.88517.80870J3797120190114060114116.39639.72318.60870J3797120190114060115116.40539.71519.50……………………


### 3.2.1. Data Processing Based on GPS Data

The GPS data are used to estimate the parameters of bus running conditions, which are collected from on-vehicle GPS facilities. The original data sets mainly have two problems: data losses and data errors. Considering existing data handling methods [12, 19], the data quality controls are conducted by the following three steps:Step 1: data processing of original data. The available fields are collected from the database, including line number, direction, vehicle number, vehicle trip, collection time, longitude, latitude, and speed, which are presented in Table1.Step 2: the fault data of speed, longitude, and latitude are filled and corrected by linear interpolation, which is described as follows:(1)xtj=tj−tixtk+tk−tjxtitk−ti,wherextj is the corrected data of time tj, (ti, tk) is the neighbour time points of tj, xtk and xti are the original values of ti and tk, respectively.Step 3: the GPS data are matched with a GIS map to analyze the spatial locations of passenger travel trajectory, as Figure2 shows. The second-by-second GPS data are aggregated to calculate IVT and WT for different bus trips and lines. In addition, the matching process is verified based on the station number information from passenger IC card data.Figure 2
GPS data matching with GIS.


### 3.2.2. Formulation of IVT Reliability

When bus passengers plan a trip, they will add the “buffer time” based on the expected travel time, considering the potential fluctuations. Although they may arrive at the destination early, early arrival is better since the negative impact of the late arrival is large [31]. As a result, most passengers may choose to plan the extra buffer time during peak hours. The buffer time can be the quantitative index to measure passenger travel time reliability [32]. The more buffer time passengers need to plan, the worse travel time reliability of the network is.As is depicted in Figure3, IVT considering reliability can be computed as follows:(2)tij=tj−ti,tijin−vehicle=Etij+ρσtij,Etij=∑n=1NtijnN,σtij=1N−1∑n=1Ntijn−Etij2,where tj and ti are samples of bus arrival time at stop i and j, which can be extracted from the GPS data. tij is sample time from node i to node j. tijin−vehicle is IVT considering reliability from node i to node j. E (tij) is the expectation of IVT from node i to node j for all the samples. σtij is the buffer time to measure the reliability, i.e., the standard deviation(SD) of the link ij.Figure 3
Illustration of IVT.ρ is reliability preference parameter. When ρ = 1, it means that the passenger will plan the equivalent buffer time as expected travel time. Passenger is “risk-neutral for on-time arrival”. When ρ = 0, the passenger does not consider the influence of travel time reliability. When ρ > 1, the passenger is regarded as conservative, in which the passenger does not trust bus service reliability and will plan more extra buffer time compared with the risk-neutral passenger. When ρ < 1, the passenger is understood as “risky for on-time arrival” [30].N is the number of samples on link (i, j).


### 3.2.3. Formulation of WT Reliability

As is presented in Figure4, assuming that passenger arrives at the bus stop on time t, and the following bus arrives on time h, WT can be calculated as follows:(3)w=h−t,twaiting=Ew+ρσw,Ew=Eh−t=Eh−Et,σw=Dw=Dh+Dt,where w is sample WT, t is passenger arrival time at the bus stop, and h is the following bus arrival interval. twaiting is WT considering reliability. E (w) is the expectation of WT for all the samples. σw is the buffer time to measure the reliability, such as SD, and D(w) is the deviation of sample WT.Figure 4
Illustration of WT.Both arrivals of passengers and bus at the stop are random, which can influence WT. Usually, each bus trip does not have a strict arrival timetable on stop level in Beijing. Passenger arrival time and bus arrival time can be considered as an independent event.Studies usually focus on passenger arrival behavior within one fixed bus arrival interval. However, with the increasing influences of traffic congestion and other uncertain factors, the variance of bus arrival time becomes larger, which should not be neglected. Bus intervalh can be regarded following a certain probability distribution h∼h′x. Based on field GPS data, statical characteristics of E(h) are estimated. The passenger arrival rate is regraded following a Poisson distribution on [0, H] [33], which is written as follows:(4)PT=t=zwHte−zwHt!,Et=∫0HtPTdT=zwH,σt=zw,where zw is the ratio of the expected passenger arrival time/bus arrival interval.The model works for high-frequency service. Considering the condition that passengers arrive randomly and the individual WT are different, the statical fitting analysis by a large amount of field data is presented in the case study. Root mean square error (RMSE) and R-square are used to validate the fitting result, which are presented in section case study.


## 3.2.1. Data Processing Based on GPS Data

The GPS data are used to estimate the parameters of bus running conditions, which are collected from on-vehicle GPS facilities. The original data sets mainly have two problems: data losses and data errors. Considering existing data handling methods [12, 19], the data quality controls are conducted by the following three steps:Step 1: data processing of original data. The available fields are collected from the database, including line number, direction, vehicle number, vehicle trip, collection time, longitude, latitude, and speed, which are presented in Table1.Step 2: the fault data of speed, longitude, and latitude are filled and corrected by linear interpolation, which is described as follows:(1)xtj=tj−tixtk+tk−tjxtitk−ti,wherextj is the corrected data of time tj, (ti, tk) is the neighbour time points of tj, xtk and xti are the original values of ti and tk, respectively.Step 3: the GPS data are matched with a GIS map to analyze the spatial locations of passenger travel trajectory, as Figure2 shows. The second-by-second GPS data are aggregated to calculate IVT and WT for different bus trips and lines. In addition, the matching process is verified based on the station number information from passenger IC card data.Figure 2
GPS data matching with GIS.


## 3.2.2. Formulation of IVT Reliability

When bus passengers plan a trip, they will add the “buffer time” based on the expected travel time, considering the potential fluctuations. Although they may arrive at the destination early, early arrival is better since the negative impact of the late arrival is large [31]. As a result, most passengers may choose to plan the extra buffer time during peak hours. The buffer time can be the quantitative index to measure passenger travel time reliability [32]. The more buffer time passengers need to plan, the worse travel time reliability of the network is.As is depicted in Figure3, IVT considering reliability can be computed as follows:(2)tij=tj−ti,tijin−vehicle=Etij+ρσtij,Etij=∑n=1NtijnN,σtij=1N−1∑n=1Ntijn−Etij2,where tj and ti are samples of bus arrival time at stop i and j, which can be extracted from the GPS data. tij is sample time from node i to node j. tijin−vehicle is IVT considering reliability from node i to node j. E (tij) is the expectation of IVT from node i to node j for all the samples. σtij is the buffer time to measure the reliability, i.e., the standard deviation(SD) of the link ij.Figure 3
Illustration of IVT.ρ is reliability preference parameter. When ρ = 1, it means that the passenger will plan the equivalent buffer time as expected travel time. Passenger is “risk-neutral for on-time arrival”. When ρ = 0, the passenger does not consider the influence of travel time reliability. When ρ > 1, the passenger is regarded as conservative, in which the passenger does not trust bus service reliability and will plan more extra buffer time compared with the risk-neutral passenger. When ρ < 1, the passenger is understood as “risky for on-time arrival” [30].N is the number of samples on link (i, j).


## 3.2.3. Formulation of WT Reliability

As is presented in Figure4, assuming that passenger arrives at the bus stop on time t, and the following bus arrives on time h, WT can be calculated as follows:(3)w=h−t,twaiting=Ew+ρσw,Ew=Eh−t=Eh−Et,σw=Dw=Dh+Dt,where w is sample WT, t is passenger arrival time at the bus stop, and h is the following bus arrival interval. twaiting is WT considering reliability. E (w) is the expectation of WT for all the samples. σw is the buffer time to measure the reliability, such as SD, and D(w) is the deviation of sample WT.Figure 4
Illustration of WT.Both arrivals of passengers and bus at the stop are random, which can influence WT. Usually, each bus trip does not have a strict arrival timetable on stop level in Beijing. Passenger arrival time and bus arrival time can be considered as an independent event.Studies usually focus on passenger arrival behavior within one fixed bus arrival interval. However, with the increasing influences of traffic congestion and other uncertain factors, the variance of bus arrival time becomes larger, which should not be neglected. Bus intervalh can be regarded following a certain probability distribution h∼h′x. Based on field GPS data, statical characteristics of E(h) are estimated. The passenger arrival rate is regraded following a Poisson distribution on [0, H] [33], which is written as follows:(4)PT=t=zwHte−zwHt!,Et=∫0HtPTdT=zwH,σt=zw,where zw is the ratio of the expected passenger arrival time/bus arrival interval.The model works for high-frequency service. Considering the condition that passengers arrive randomly and the individual WT are different, the statical fitting analysis by a large amount of field data is presented in the case study. Root mean square error (RMSE) and R-square are used to validate the fitting result, which are presented in section case study.


## 3.3. Bilevel Programming Formulation

### 3.3.1. Formulation of General Cost

The general costs are combined with IVT cost, WT cost, extra waiting delay cost, transfer cost, and a fare cost. For each arc that passenger passes, the formulation of general cost can be calculated as follows:(1)
IVT cost(5)C=utsin−vehicle,∀s∈Sin−vehicle,whereC is general cost, u is unit travel time cost (yuan/h), s is travel arc, and Sin−vehicle is set of in-vehicle travel arc.(2)
WT cost(6)C=utswaiting,∀s∈Swaiting,whereSwaiting is set of waiting arc.(3)
Extra waiting delay time costDuring peak hours, it is a common phenomenon that there is a large number of passengers boarding and alighting at bus stops, which will cause the speed of boarding and alighting to slow down [34]. As a result, there will be an extra delay which is increased with the rising number of passengers. The extra delay is related to the number of existing passengers at the stop. In this paper, the extra delay is determined as follows:(7)Ds=βsYsKn,∀s∈Swaiting,C=uDs,∀s∈Swaiting,Ds is extra delay for waiting arc s. Ys is volume of waiting passengers. K is rated carrying capacity. βs and n are estimated parameters.(4)
Transfer time costUnder the assumption that most passengers will choose to transfer at the same stop, the passenger transfer time is regarded as one additional WT. Therefore, the transfer time cost is expressed as follows:(8)C=utswaiting,∀s∈Stransfer,whereStransfer is set of transfer arc.(5)
Fare cost(9)C=u′ls,∀s∈Sin−vehicle,u′ is unit ticket cost (yuan/km). ls is travel distances (km).


### 3.3.2. Upper-Level Programming Model of Shuttle Route Design

Last-mile shuttle bus is designed as a measure to relieve congested demand and provide more reliable transit service for passengers. Therefore, the maximum improvement benefit is regarded as the objective function. For the diverted passengers, they can take the more flexible and reliable shuttle bus, but they should burden the additional fare cost. The model is formulated by the following optimization model.(10)maxC=∑o∑dCTod¯yod−∑o∑d∑i∑jutijin−vehicleyodxijod−∑o∑dcsodyod−∑o∑d∑i∑ju′lijyodxijod.(11)s.t.∑i=oxijod=1,∀o,d∈V,∀i,j∈S,(12)∑j=dxijod=1,∀o,d∈V,∀i,j∈S,(13)∑jxijod=∑jxjiod,∀o,d∈V,∀i,j∈S,(14)qod=dod+yod,∀o,d∈V,(15)kod≥k1,∀o,d∈V,(16)yod≤Kfl,∀o,d∈V,(17)xijod∈0,1,∀o,d∈V,∀i,j∈S,(18)dod≥0,∀o,d∈V,(19)yod≥0,∀o,d∈V,where CTod¯ is the average time cost (yuan) and csod is the average ticket cost of conventional bus (yuan/km). The decision variables are xijod and yod, in which xijod means that whether the shuttle from node O to node D passes through road arc (i, j) and yod reflects passengers that take the shuttles.∑o∑dCTodyod¯ and ∑o∑dcsodyod reflect total time cost and ticket cost when passengers take conventional transit. ∑o∑d∑i∑jutijin−vehicleyodxijod and ∑o∑d∑i∑ju′lijyodxijod are total time cost and ticket cost for shuttle passengers. The objective is to maximize passenger benefit when diverted passengers choose the new shuttle service.Equations (11)–(13) are the flow conservation constraints for each node, in which V is set of a bus stop.Equation (14) means that for each OD pair, the total OD demand should equal the sum of shuttle demand and the conventional bus demand, in which qod reflects overall demand and dod is the remaining demand for conventional transit.Equation (15) reflects the reliability threshold, which means that only when kod is larger than the certain level k1, the shuttles will begin to run. Kod quantifies the impact of time reliability, which is the key threshold for judging whether it is necessary to provide shuttle service. The larger the kod is, the worse the reliability condition is. kod is calculated as follows:(20)kod=CRodCTod,∀o,d∈V,CRod=∑suρσs+Dsysod,∀o,d∈V,∀s∈Sin−vehicle∪​Swaiting,CTod=∑s∈Sin−vehicleutsin−vehicleysod+∑s∈Swaitingutswaiting+Dsysod,∀o,d∈V,where CRod reflects the delay time cost which is produced by the unreliability of travel time and CTod represents the total time cost from O to D. ysod is demand volume on arc s.Equation (16) reflects the vehicle carrying capacity constraint, in which fl is shuttle service frequency of line l.Equations (17)–(19) give the value range for the variables.


### 3.3.3. Lower-Level Programming Model of Demand Assignment

Lower-level model aims to assign the existing bus demand, estimate travel time and cost on the network arc, and quantify the reliability condition. In the lower model, bus passenger IVT cost, WT cost, extra delay cost, transfer time cost, and fare cost are considered.Based on the mentioned analyses, a lower-level flow assignment model based on user equilibrium is formulated as follows:(21)minC=∑o∑d∑p∫0YsCsYsdYs.(22)s.t.∑p∈Pyodp=qod,∀o∈V,∀d∈V,(23)yodp≥0,∀o∈V,∀d∈V,∀p∈P,where CsYs is the general cost function of arc s. Ys is the total volume of passengers on arc s, which is derived that Ys=∑i∑j∑pyodpδodps. yodp is the volume of passengers from i to j on available path p. δodps is the correlation parameter that if travel arc s is on path p from o to d, δodps  = 1, otherwise, δodps  = 0. P is a set of available paths.The objective function is the sum of the results of integrating the general cost functions. The model has two constraints. For (22), it means that the sum of passengers on all paths is equal to the volume of the demand for each OD pair. Constraint (23) means that the decision variable is nonnegative. The optimal solution of the objective function is equivalent to user equilibrium.


## 3.3.1. Formulation of General Cost

The general costs are combined with IVT cost, WT cost, extra waiting delay cost, transfer cost, and a fare cost. For each arc that passenger passes, the formulation of general cost can be calculated as follows:(1)
IVT cost(5)C=utsin−vehicle,∀s∈Sin−vehicle,whereC is general cost, u is unit travel time cost (yuan/h), s is travel arc, and Sin−vehicle is set of in-vehicle travel arc.(2)
WT cost(6)C=utswaiting,∀s∈Swaiting,whereSwaiting is set of waiting arc.(3)
Extra waiting delay time costDuring peak hours, it is a common phenomenon that there is a large number of passengers boarding and alighting at bus stops, which will cause the speed of boarding and alighting to slow down [34]. As a result, there will be an extra delay which is increased with the rising number of passengers. The extra delay is related to the number of existing passengers at the stop. In this paper, the extra delay is determined as follows:(7)Ds=βsYsKn,∀s∈Swaiting,C=uDs,∀s∈Swaiting,Ds is extra delay for waiting arc s. Ys is volume of waiting passengers. K is rated carrying capacity. βs and n are estimated parameters.(4)
Transfer time costUnder the assumption that most passengers will choose to transfer at the same stop, the passenger transfer time is regarded as one additional WT. Therefore, the transfer time cost is expressed as follows:(8)C=utswaiting,∀s∈Stransfer,whereStransfer is set of transfer arc.(5)
Fare cost(9)C=u′ls,∀s∈Sin−vehicle,u′ is unit ticket cost (yuan/km). ls is travel distances (km).


## 3.3.2. Upper-Level Programming Model of Shuttle Route Design

Last-mile shuttle bus is designed as a measure to relieve congested demand and provide more reliable transit service for passengers. Therefore, the maximum improvement benefit is regarded as the objective function. For the diverted passengers, they can take the more flexible and reliable shuttle bus, but they should burden the additional fare cost. The model is formulated by the following optimization model.(10)maxC=∑o∑dCTod¯yod−∑o∑d∑i∑jutijin−vehicleyodxijod−∑o∑dcsodyod−∑o∑d∑i∑ju′lijyodxijod.(11)s.t.∑i=oxijod=1,∀o,d∈V,∀i,j∈S,(12)∑j=dxijod=1,∀o,d∈V,∀i,j∈S,(13)∑jxijod=∑jxjiod,∀o,d∈V,∀i,j∈S,(14)qod=dod+yod,∀o,d∈V,(15)kod≥k1,∀o,d∈V,(16)yod≤Kfl,∀o,d∈V,(17)xijod∈0,1,∀o,d∈V,∀i,j∈S,(18)dod≥0,∀o,d∈V,(19)yod≥0,∀o,d∈V,where CTod¯ is the average time cost (yuan) and csod is the average ticket cost of conventional bus (yuan/km). The decision variables are xijod and yod, in which xijod means that whether the shuttle from node O to node D passes through road arc (i, j) and yod reflects passengers that take the shuttles.∑o∑dCTodyod¯ and ∑o∑dcsodyod reflect total time cost and ticket cost when passengers take conventional transit. ∑o∑d∑i∑jutijin−vehicleyodxijod and ∑o∑d∑i∑ju′lijyodxijod are total time cost and ticket cost for shuttle passengers. The objective is to maximize passenger benefit when diverted passengers choose the new shuttle service.Equations (11)–(13) are the flow conservation constraints for each node, in which V is set of a bus stop.Equation (14) means that for each OD pair, the total OD demand should equal the sum of shuttle demand and the conventional bus demand, in which qod reflects overall demand and dod is the remaining demand for conventional transit.Equation (15) reflects the reliability threshold, which means that only when kod is larger than the certain level k1, the shuttles will begin to run. Kod quantifies the impact of time reliability, which is the key threshold for judging whether it is necessary to provide shuttle service. The larger the kod is, the worse the reliability condition is. kod is calculated as follows:(20)kod=CRodCTod,∀o,d∈V,CRod=∑suρσs+Dsysod,∀o,d∈V,∀s∈Sin−vehicle∪​Swaiting,CTod=∑s∈Sin−vehicleutsin−vehicleysod+∑s∈Swaitingutswaiting+Dsysod,∀o,d∈V,where CRod reflects the delay time cost which is produced by the unreliability of travel time and CTod represents the total time cost from O to D. ysod is demand volume on arc s.Equation (16) reflects the vehicle carrying capacity constraint, in which fl is shuttle service frequency of line l.Equations (17)–(19) give the value range for the variables.


## 3.3.3. Lower-Level Programming Model of Demand Assignment

Lower-level model aims to assign the existing bus demand, estimate travel time and cost on the network arc, and quantify the reliability condition. In the lower model, bus passenger IVT cost, WT cost, extra delay cost, transfer time cost, and fare cost are considered.Based on the mentioned analyses, a lower-level flow assignment model based on user equilibrium is formulated as follows:(21)minC=∑o∑d∑p∫0YsCsYsdYs.(22)s.t.∑p∈Pyodp=qod,∀o∈V,∀d∈V,(23)yodp≥0,∀o∈V,∀d∈V,∀p∈P,where CsYs is the general cost function of arc s. Ys is the total volume of passengers on arc s, which is derived that Ys=∑i∑j∑pyodpδodps. yodp is the volume of passengers from i to j on available path p. δodps is the correlation parameter that if travel arc s is on path p from o to d, δodps  = 1, otherwise, δodps  = 0. P is a set of available paths.The objective function is the sum of the results of integrating the general cost functions. The model has two constraints. For (22), it means that the sum of passengers on all paths is equal to the volume of the demand for each OD pair. Constraint (23) means that the decision variable is nonnegative. The optimal solution of the objective function is equivalent to user equilibrium.


## 3.4. Solution Algorithm

Each of the upper-level solutions represents a scheme of shuttle route planning. For each scheme, demand should be assigned by the lower-level model to estimate whether overall reliability is improved to the acceptable threshold (k1). Then, the result of lower-level will be substituted into the upper level to calculate the objective function. In addition, considering the upper model is a nonlinear programming formulation, the computational complexity will be higher, which can be solved by a linearization algorithm.Based on the above analysis, a two-stage algorithm is proposed, which is depicted in Figure5. The first step is for demand assignment, which estimates kod of existing conventional network and selects unreliable OD, which can be diverted to shuttle service. The second step is the shuttle route optimization model including a linearization algorithm to solve the 0-1 nonlinearization programming model. The optimal shuttle service plan will be calculated finally. The algorithm is described in detail as follows:Stage 1: conduct the lower-level demand assignment model based on the existing local bus network and demand. Estimatekod among all OD in the local bus network. Select the congested OD by reliability threshold k1: if kod >k1, goes for stage 2. Otherwise, the selection process ends.Stage 2: conduct the upper-level model. Linearization transforming is introduced for the model. The shuttle route and diverted demand are output finally by the model. The linearization algorithm is described in detail.Figure 5
Overall solution framework.For upper model, the two decision variables are nonnegative integer variables and 0-1 variables. As equation (10) shows, yodxijod is the nonlinear structure. Auxiliary decision variable Mijod (Mijod≥0) and real number M (M≥0) are introduced to construct auxiliary inequations. Auxiliary constraints are introduced as the following equations show:(24)yod−M1−xijod≤Mijod≤yod,∀o,d∈V,∀i,j∈S,(25)Mijod≤Mxijod,∀o,d∈V,∀i,j∈S.Whenxijod=1， yod≤Mijod≤yod. Therefore, Mijod=yodxijod. When xijod=0, according to inequality (25), Mijod≤0. Therefore Mijod=0. There are following inequalities that Mijod≥0 and yod≥0. Within the range of yod and xijod, Mijod always equals yodxijod. As a result, the nonlinear expression yodxijod is transformed into linear expressions with a series of inequalities.The solving difficulty is significantly reduced. In the paper, the commercial linear programming software Gurobi is used to solve the model. After the linearization, the model can be expressed as follows:(26)maxC=∑o∑dCTod¯yod−∑o∑d∑i∑jutijin vehicleMijod−∑o∑dcsodyod−∑o∑d∑i∑ju′lijMijod,s.t.∑i=oxijod=1,∀o,d∈V,∀i,j∈S,∑j=dxijod=1,∀o,d∈V,∀i,j∈S,∑jxijod=∑jxjiod,∀o,d∈V,∀i,j∈S,qod=dod+yod,∀o,d∈V,kod≥k1,∀o,d∈V,yod≤kfl,∀o,d∈V,yod−M1−xijod≤Mijod≤yod,∀o,d∈V,∀i,j∈S,Mijod≤Mxijod,∀o,d∈V,∀i,j∈S,xijod∈0,1,∀o,d∈V,∀i,j∈S,dod≥0,∀o,d∈V,yod≥0,∀o,d∈V.


## 4. Case Study

### 4.1. Characteristics of Local Bus Network

The last-mile shuttle optimization model is further validated in the local bus network of the Zhongguancun Software Park zone. Zhongguancun Software Park is located in the Haidian district in Beijing, in which many Internet enterprises such as Baidu, Tencent, and Lenovo are located. Most commuters will choose to take the metro to station Xierqi, then transfer to bus system for their last-mile trips. As a high-density workplace zone, the travel time reliability is always insufficient for the last-mile trip during peak hours.Figure6(a) shows the geographic layout of Zhongguancun Software Park and Figure 6(b) shows the structure of the local bus network. There are 12 existing bus lines in the local area. The distance between bus stops and passenger demand is shown in Tables 2 and 3.Figure 6
(a) Road network of the study area. (b) Conventional bus network.
(a)(b)Table 2
Distance between bus stops.
Bus stopDistance (km)Bus stopDistance (km)Bus stopDistance (km)1⟶41.52⟶30.510⟶110.34⟶81.91⟶9111⟶120.41⟶319⟶120.52⟶90.63⟶40.512⟶130.39⟶251.54⟶6113⟶141.81⟶271.36⟶80.914⟶150.727⟶280.74⟶50.715⟶160.61⟶221.35⟶70.816⟶170.722⟶230.61⟶20.31⟶101.123⟶240.8Table 3
OD matrices.
OD (person)235891012131415161718222324261269620937122702291492327139913331095923461119413132315662215232314The average ticket costs of conventional transit lines and shuttle lines are 0.3 yuan/km and 0.6 yuan/km, which is estimated by the existing ticket rate in Beijing. Unit values of travel time and travel time reliability are 11.34 yuan/h and 19.27 yuan/h [31]. Reliability preference is 2 in the case study. Other parameters for the optimization model are described in Table 4.Table 4
Parameter values.
k (shuttle, person/vehicle)k (conventional bus, person/vehicle)fl (vehicle/h)n408063βsk10.50.75


### 4.2. Estimation of Travel Time Reliability

The reliability of the existing transit system is estimated. Field data are collected and processed for the fitting analysis. In addition, GPS data extracted from existing microcirculation lines are chosen to reflect the operation conditions of the proposed shuttle service. In Beijing, microcirculation transit has been developed in recent years, which provides the service in community or business districts. Although it has a fixed route and stops, compared with the other conventional bus, microcirculation bus has the unique operation characteristics that the route length and stop spacing are short, the operating time is during peak hour, and the light-duty vehicles are used, which are similar to the designed shuttle service. Data of 11 conventional lines and 7 microcirculation lines are chosen in the case study. The GPS data cover the period of 6:00–11:00 AM on Monday, with more than five hundred thousand pieces of basic data.


#### 4.2.1. Reliability of IVT

First, overall IVT spatial and temporal distribution are described in Figure7 as an illustration, where the X-axis represents time, the Y-axis represents the section between neighbouring two stops, and the Z axis stands for IVT. It is indicated that for microcirculation lines, IVT variability is significantly large. In addition, it is evidently greater after 7:00 than that before 7:00. However, there is a significant reduction after 9:00. It means that after peak hours, traffic condition becomes better in the running area of microcirculation lines.Figure 7
Spatial and temporal distribution of sample bus lines. (a) Microcirculation lines. (b) Conventional lines.
(a)(b)On the basis of the IVT samples data, the distributions are fitted to explore its statistical characteristics. Four typical fits are applied to test which kind of statistical distribution can best explain the distribution characters of IVT: Gaussian fit, Weibull fit, Laplace fit, and Lognormal fit. Root mean square error (RMSE) and R-square are used to validate the fitting results. The fitting distribution results are shown in Table5. It is concluded that Lognormal distribution can best explain the distribution characteristics of IVT.Table 5
Fitting results of IVT.
LineGaussianWeibullLaplaceLognormalParameters of best fit320 upwardsR-square0.9580.9520.9340.996a = 1.315RMSE0.0190.0190.0220.006b = 0.476320 downwardsR-square0.9650.9790.9140.987a = 1.339RMSE0.0150.0110.0570.009b = 0.406333 upwardsR-square0.9650.9580.9520.992a = 1.328RMSE0.0150.0160.0180.007b = 0.421333 downwardsR-square0.9210.8790.9350.956a = 1.275RMSE0.0210.0260.0190.015b = 0.483362 upwardsR-square0.9520.9520.9350.995a = 1.12RMSE0.0200.0200.0230.006b = 0.445362 downwardsR-square0.9510.9540.9350.990a = 1.113RMSE0.0190.0180.0210.008b = 0.521509 upwardsR-square0.9430.9340.9340.986a = 1.211RMSE0.0210.0230.0230.011b = 0.449509 downwardsR-square0.9270.9180.9480.972a = 1.430RMSE0.0180.0190.0160.011b = 0.495521 upwardsR-square0.9490.9670.9530.980a = 1.143RMSE0.0170.0150.0170.011b = 0.644521 downwardsR-square0.9420.9440.9540.982a = 1.024RMSE0.0220.0210.0200.012b = 0.562570 upwardsR-square0.9730.9730.9380.996a = 1.145RMSE0.0160.0150.0230.006b = 0.402570 downwardsR-square0.9750.9660.9530.992a = 1.192RMSE0.0160.0180.0210.009b = 0.377963 upwardsR-square0.9330.9450.9460.989a = 1.384RMSE0.0150.0130.0140.006b = 0.576963 downwardsR-square0.9180.9230.9310.981a = 1.319RMSE0.0190.0180.0170.009b = 0.478105MR-square0.9770.9680.9630.986a = 1.219RMSE0.0180.0210.0250.014b = 0.346116MR-square0.9590.9670.9520.992a = 1.260RMSE0.0100.0090.01090.004b = 0.501117MR-square0.9320.9510.9290.983a = 1.317RMSE0.0220.0190.0230.011b = 0.485118MR-square0.9600.9710.9490.986a = 1.242RMSE0.0170.0150.0190.010b = 0.496120MR-square0.8940.9180.8900.965a = 1.433RMSE0.0240.0220.0250.014b = 0.463121MR-square0.8830.9380.8930.983a = 1.117RMSE0.0280.0210.0270.011b = 0.531123MR-square0.9530.9520.9380.990a = 1.339RMSE0.0190.0190.0210.009b = 0.406


#### 4.2.2. Reliability of WT

Bus departure headway at the origin stop can directly affect the fluctuation of bus arrival intervals at each midway stop. For ideal scenarios, the arrival interval at the midway stop equals the departure headway. However, due to various uncertain influencing factors, the arrival interval becomes fluctuated. In actual conditions, the departure headway for different bus trips, lines, and running periods is also different. Therefore, based on the stop-level arrival intervals extracted from GPS data, bus arrival interval distribution is fitted under different departure headways. Fitting results are presented in Table6, where H0 represents bus departure headway. 10 scenarios are designed, in which H0 equals from 1 minute to 10 minutes, respectively. Gaussian fit, Weibull fit, Laplace fit, and Lognormal fit are applied to test the best fit of bus arrival intervals. The result indicates that the Laplace fit shows the optimal fitting result.Table 6
Fitting results of bus arrival intervals.
Departure headwayGaussianWeibullLaplaceLognormalParameters for the best fitH0 = 1 minR-square0.9040.9150.9420.828a = 0.972RMSE0.0450.0420.0340.060b = 0.734H0 = 2 minsR-square0.9360.9340.9730.886a = 1.606RMSE0.0220.0220.0140.029b = 2.171H0 = 3 minsR-square0.9020.9040.9620.867a = 2.669RMSE0.0160.0160.0100.019b = 2.740H0 = 4 minsR-square0.8790.8770.9690.757a = 3.018RMSE0.0160.0160.0080.022b = 3.861H0 = 5 minsR-square0.8880.8700.9680.836a = 2.908RMSE0.0160.0170.0080.019b = 4.795H0 = 6 minsR-square0.9120.8790.9890.851a = 2.770RMSE0.0140.0170.0050.018b = 6.021H0 = 7 minsR-square0.9170.8830.9890.856a = 2.626RMSE0.0140.0180.0050.019b = 7.117H0 = 8 minsR-square0.8980.8620.9890.819a = 3.080RMSE0.0140.0160.0060.018b = 8.121H0 = 9 minsR-square0.8590.8190.9760.894a = 2.903RMSE0.0170.0190.0070.013b = 8.849H0 = 10 minsR-square0.9170.8900.9870.919a = 3.185RMSE0.0120.0130.0040.008b = 10.16Following the above fitting analysis, the expectation and standard deviation of passenger IVT and WT in Zhongguancun Software Park is estimated and analyzing result is shown in Table7.Table 7
Expectation and standard deviation of IVT and WT.
LineIVTWTE (min/km)SD (min/km)E (min)SD (min)9023.73.32.45.73203.72.45.15.63333.72.43.65.63622.92.03.65.65093.32.25.16.85213.72.74.16.45703.01.95.16.89633.82.62.45.76363.22.24.46.2853.01.95.16.8823.22.25.16.8


### 4.3. Results and Discussion

CRnetwork is introduced to estimate the overall reliability and quantify the improving effect. CRnetwork represents the sum of extra time cost due to the unreliability of the travel time for all passengers in the network, which can be calculated as follows:(27)CRnetwork=∑o∑d∑suρσs+DsYs,∀o,d∈V,∀s∈Sin−vehicle∪​Swaiting.Based on the optimization result, in local bus network of Zhongguancun software, stop Metro Xierqi to stop Dongbeiwang is the hot OD pair with insufficient reliability during peak hours. The congested travel demand can be diverted by shuttle service. The detailed plan is described in Table8. The number of passengers by shuttle is 240 person/h, which takes account of 12.2% of total transit demand. The running route of the shuttle is: stop Metro Xierqi-Road Ruanjianyuan-stop Dongbeiwang. For the original local bus network, CRnetwork is 15650.8 yuan. After the operating of shuttle, CRnetwork is 13524.9 yuan, which has a significant decrease by 14%.Table 8
Optimization results.
Last-mile shuttle routeDiverting demandOptimal cost of travel time reliabilityStop metro xierqi—stop road ruanjianyuan-dongbeiwang240 person/hOrigin conventional bus networkAfter shuttle operatedOptimization effect15651 yuan13525 yuan14%In the area of Zhongguancun Software Park, the spatial and temporal aggregation of transit demand is significant. Most demand occurrence points (point O) are located at the stop Metro Xierqi station, and the destination points are concentrated in several key stops near the office area. Meanwhile, commuting demand is concentrated in peak hours, and the roads are of low grade, which is narrow and not suitable for oversize vehicles. Although the working area in the study is small, the results show that the shuttle service can significantly improve the travel time reliability of last-mile trips. The applicable scenarios can be expanded to more similar areas such as residential areas or school areas.


## 4.1. Characteristics of Local Bus Network

The last-mile shuttle optimization model is further validated in the local bus network of the Zhongguancun Software Park zone. Zhongguancun Software Park is located in the Haidian district in Beijing, in which many Internet enterprises such as Baidu, Tencent, and Lenovo are located. Most commuters will choose to take the metro to station Xierqi, then transfer to bus system for their last-mile trips. As a high-density workplace zone, the travel time reliability is always insufficient for the last-mile trip during peak hours.Figure6(a) shows the geographic layout of Zhongguancun Software Park and Figure 6(b) shows the structure of the local bus network. There are 12 existing bus lines in the local area. The distance between bus stops and passenger demand is shown in Tables 2 and 3.Figure 6
(a) Road network of the study area. (b) Conventional bus network.
(a)(b)Table 2
Distance between bus stops.
Bus stopDistance (km)Bus stopDistance (km)Bus stopDistance (km)1⟶41.52⟶30.510⟶110.34⟶81.91⟶9111⟶120.41⟶319⟶120.52⟶90.63⟶40.512⟶130.39⟶251.54⟶6113⟶141.81⟶271.36⟶80.914⟶150.727⟶280.74⟶50.715⟶160.61⟶221.35⟶70.816⟶170.722⟶230.61⟶20.31⟶101.123⟶240.8Table 3
OD matrices.
OD (person)235891012131415161718222324261269620937122702291492327139913331095923461119413132315662215232314The average ticket costs of conventional transit lines and shuttle lines are 0.3 yuan/km and 0.6 yuan/km, which is estimated by the existing ticket rate in Beijing. Unit values of travel time and travel time reliability are 11.34 yuan/h and 19.27 yuan/h [31]. Reliability preference is 2 in the case study. Other parameters for the optimization model are described in Table 4.Table 4
Parameter values.
k (shuttle, person/vehicle)k (conventional bus, person/vehicle)fl (vehicle/h)n408063βsk10.50.75


## 4.2. Estimation of Travel Time Reliability

The reliability of the existing transit system is estimated. Field data are collected and processed for the fitting analysis. In addition, GPS data extracted from existing microcirculation lines are chosen to reflect the operation conditions of the proposed shuttle service. In Beijing, microcirculation transit has been developed in recent years, which provides the service in community or business districts. Although it has a fixed route and stops, compared with the other conventional bus, microcirculation bus has the unique operation characteristics that the route length and stop spacing are short, the operating time is during peak hour, and the light-duty vehicles are used, which are similar to the designed shuttle service. Data of 11 conventional lines and 7 microcirculation lines are chosen in the case study. The GPS data cover the period of 6:00–11:00 AM on Monday, with more than five hundred thousand pieces of basic data.


### 4.2.1. Reliability of IVT

First, overall IVT spatial and temporal distribution are described in Figure7 as an illustration, where the X-axis represents time, the Y-axis represents the section between neighbouring two stops, and the Z axis stands for IVT. It is indicated that for microcirculation lines, IVT variability is significantly large. In addition, it is evidently greater after 7:00 than that before 7:00. However, there is a significant reduction after 9:00. It means that after peak hours, traffic condition becomes better in the running area of microcirculation lines.Figure 7
Spatial and temporal distribution of sample bus lines. (a) Microcirculation lines. (b) Conventional lines.
(a)(b)On the basis of the IVT samples data, the distributions are fitted to explore its statistical characteristics. Four typical fits are applied to test which kind of statistical distribution can best explain the distribution characters of IVT: Gaussian fit, Weibull fit, Laplace fit, and Lognormal fit. Root mean square error (RMSE) and R-square are used to validate the fitting results. The fitting distribution results are shown in Table5. It is concluded that Lognormal distribution can best explain the distribution characteristics of IVT.Table 5
Fitting results of IVT.
LineGaussianWeibullLaplaceLognormalParameters of best fit320 upwardsR-square0.9580.9520.9340.996a = 1.315RMSE0.0190.0190.0220.006b = 0.476320 downwardsR-square0.9650.9790.9140.987a = 1.339RMSE0.0150.0110.0570.009b = 0.406333 upwardsR-square0.9650.9580.9520.992a = 1.328RMSE0.0150.0160.0180.007b = 0.421333 downwardsR-square0.9210.8790.9350.956a = 1.275RMSE0.0210.0260.0190.015b = 0.483362 upwardsR-square0.9520.9520.9350.995a = 1.12RMSE0.0200.0200.0230.006b = 0.445362 downwardsR-square0.9510.9540.9350.990a = 1.113RMSE0.0190.0180.0210.008b = 0.521509 upwardsR-square0.9430.9340.9340.986a = 1.211RMSE0.0210.0230.0230.011b = 0.449509 downwardsR-square0.9270.9180.9480.972a = 1.430RMSE0.0180.0190.0160.011b = 0.495521 upwardsR-square0.9490.9670.9530.980a = 1.143RMSE0.0170.0150.0170.011b = 0.644521 downwardsR-square0.9420.9440.9540.982a = 1.024RMSE0.0220.0210.0200.012b = 0.562570 upwardsR-square0.9730.9730.9380.996a = 1.145RMSE0.0160.0150.0230.006b = 0.402570 downwardsR-square0.9750.9660.9530.992a = 1.192RMSE0.0160.0180.0210.009b = 0.377963 upwardsR-square0.9330.9450.9460.989a = 1.384RMSE0.0150.0130.0140.006b = 0.576963 downwardsR-square0.9180.9230.9310.981a = 1.319RMSE0.0190.0180.0170.009b = 0.478105MR-square0.9770.9680.9630.986a = 1.219RMSE0.0180.0210.0250.014b = 0.346116MR-square0.9590.9670.9520.992a = 1.260RMSE0.0100.0090.01090.004b = 0.501117MR-square0.9320.9510.9290.983a = 1.317RMSE0.0220.0190.0230.011b = 0.485118MR-square0.9600.9710.9490.986a = 1.242RMSE0.0170.0150.0190.010b = 0.496120MR-square0.8940.9180.8900.965a = 1.433RMSE0.0240.0220.0250.014b = 0.463121MR-square0.8830.9380.8930.983a = 1.117RMSE0.0280.0210.0270.011b = 0.531123MR-square0.9530.9520.9380.990a = 1.339RMSE0.0190.0190.0210.009b = 0.406


### 4.2.2. Reliability of WT

Bus departure headway at the origin stop can directly affect the fluctuation of bus arrival intervals at each midway stop. For ideal scenarios, the arrival interval at the midway stop equals the departure headway. However, due to various uncertain influencing factors, the arrival interval becomes fluctuated. In actual conditions, the departure headway for different bus trips, lines, and running periods is also different. Therefore, based on the stop-level arrival intervals extracted from GPS data, bus arrival interval distribution is fitted under different departure headways. Fitting results are presented in Table6, where H0 represents bus departure headway. 10 scenarios are designed, in which H0 equals from 1 minute to 10 minutes, respectively. Gaussian fit, Weibull fit, Laplace fit, and Lognormal fit are applied to test the best fit of bus arrival intervals. The result indicates that the Laplace fit shows the optimal fitting result.Table 6
Fitting results of bus arrival intervals.
Departure headwayGaussianWeibullLaplaceLognormalParameters for the best fitH0 = 1 minR-square0.9040.9150.9420.828a = 0.972RMSE0.0450.0420.0340.060b = 0.734H0 = 2 minsR-square0.9360.9340.9730.886a = 1.606RMSE0.0220.0220.0140.029b = 2.171H0 = 3 minsR-square0.9020.9040.9620.867a = 2.669RMSE0.0160.0160.0100.019b = 2.740H0 = 4 minsR-square0.8790.8770.9690.757a = 3.018RMSE0.0160.0160.0080.022b = 3.861H0 = 5 minsR-square0.8880.8700.9680.836a = 2.908RMSE0.0160.0170.0080.019b = 4.795H0 = 6 minsR-square0.9120.8790.9890.851a = 2.770RMSE0.0140.0170.0050.018b = 6.021H0 = 7 minsR-square0.9170.8830.9890.856a = 2.626RMSE0.0140.0180.0050.019b = 7.117H0 = 8 minsR-square0.8980.8620.9890.819a = 3.080RMSE0.0140.0160.0060.018b = 8.121H0 = 9 minsR-square0.8590.8190.9760.894a = 2.903RMSE0.0170.0190.0070.013b = 8.849H0 = 10 minsR-square0.9170.8900.9870.919a = 3.185RMSE0.0120.0130.0040.008b = 10.16Following the above fitting analysis, the expectation and standard deviation of passenger IVT and WT in Zhongguancun Software Park is estimated and analyzing result is shown in Table7.Table 7
Expectation and standard deviation of IVT and WT.
LineIVTWTE (min/km)SD (min/km)E (min)SD (min)9023.73.32.45.73203.72.45.15.63333.72.43.65.63622.92.03.65.65093.32.25.16.85213.72.74.16.45703.01.95.16.89633.82.62.45.76363.22.24.46.2853.01.95.16.8823.22.25.16.8


## 4.2.1. Reliability of IVT

First, overall IVT spatial and temporal distribution are described in Figure7 as an illustration, where the X-axis represents time, the Y-axis represents the section between neighbouring two stops, and the Z axis stands for IVT. It is indicated that for microcirculation lines, IVT variability is significantly large. In addition, it is evidently greater after 7:00 than that before 7:00. However, there is a significant reduction after 9:00. It means that after peak hours, traffic condition becomes better in the running area of microcirculation lines.Figure 7
Spatial and temporal distribution of sample bus lines. (a) Microcirculation lines. (b) Conventional lines.
(a)(b)On the basis of the IVT samples data, the distributions are fitted to explore its statistical characteristics. Four typical fits are applied to test which kind of statistical distribution can best explain the distribution characters of IVT: Gaussian fit, Weibull fit, Laplace fit, and Lognormal fit. Root mean square error (RMSE) and R-square are used to validate the fitting results. The fitting distribution results are shown in Table5. It is concluded that Lognormal distribution can best explain the distribution characteristics of IVT.Table 5
Fitting results of IVT.
LineGaussianWeibullLaplaceLognormalParameters of best fit320 upwardsR-square0.9580.9520.9340.996a = 1.315RMSE0.0190.0190.0220.006b = 0.476320 downwardsR-square0.9650.9790.9140.987a = 1.339RMSE0.0150.0110.0570.009b = 0.406333 upwardsR-square0.9650.9580.9520.992a = 1.328RMSE0.0150.0160.0180.007b = 0.421333 downwardsR-square0.9210.8790.9350.956a = 1.275RMSE0.0210.0260.0190.015b = 0.483362 upwardsR-square0.9520.9520.9350.995a = 1.12RMSE0.0200.0200.0230.006b = 0.445362 downwardsR-square0.9510.9540.9350.990a = 1.113RMSE0.0190.0180.0210.008b = 0.521509 upwardsR-square0.9430.9340.9340.986a = 1.211RMSE0.0210.0230.0230.011b = 0.449509 downwardsR-square0.9270.9180.9480.972a = 1.430RMSE0.0180.0190.0160.011b = 0.495521 upwardsR-square0.9490.9670.9530.980a = 1.143RMSE0.0170.0150.0170.011b = 0.644521 downwardsR-square0.9420.9440.9540.982a = 1.024RMSE0.0220.0210.0200.012b = 0.562570 upwardsR-square0.9730.9730.9380.996a = 1.145RMSE0.0160.0150.0230.006b = 0.402570 downwardsR-square0.9750.9660.9530.992a = 1.192RMSE0.0160.0180.0210.009b = 0.377963 upwardsR-square0.9330.9450.9460.989a = 1.384RMSE0.0150.0130.0140.006b = 0.576963 downwardsR-square0.9180.9230.9310.981a = 1.319RMSE0.0190.0180.0170.009b = 0.478105MR-square0.9770.9680.9630.986a = 1.219RMSE0.0180.0210.0250.014b = 0.346116MR-square0.9590.9670.9520.992a = 1.260RMSE0.0100.0090.01090.004b = 0.501117MR-square0.9320.9510.9290.983a = 1.317RMSE0.0220.0190.0230.011b = 0.485118MR-square0.9600.9710.9490.986a = 1.242RMSE0.0170.0150.0190.010b = 0.496120MR-square0.8940.9180.8900.965a = 1.433RMSE0.0240.0220.0250.014b = 0.463121MR-square0.8830.9380.8930.983a = 1.117RMSE0.0280.0210.0270.011b = 0.531123MR-square0.9530.9520.9380.990a = 1.339RMSE0.0190.0190.0210.009b = 0.406


## 4.2.2. Reliability of WT

Bus departure headway at the origin stop can directly affect the fluctuation of bus arrival intervals at each midway stop. For ideal scenarios, the arrival interval at the midway stop equals the departure headway. However, due to various uncertain influencing factors, the arrival interval becomes fluctuated. In actual conditions, the departure headway for different bus trips, lines, and running periods is also different. Therefore, based on the stop-level arrival intervals extracted from GPS data, bus arrival interval distribution is fitted under different departure headways. Fitting results are presented in Table6, where H0 represents bus departure headway. 10 scenarios are designed, in which H0 equals from 1 minute to 10 minutes, respectively. Gaussian fit, Weibull fit, Laplace fit, and Lognormal fit are applied to test the best fit of bus arrival intervals. The result indicates that the Laplace fit shows the optimal fitting result.Table 6
Fitting results of bus arrival intervals.
Departure headwayGaussianWeibullLaplaceLognormalParameters for the best fitH0 = 1 minR-square0.9040.9150.9420.828a = 0.972RMSE0.0450.0420.0340.060b = 0.734H0 = 2 minsR-square0.9360.9340.9730.886a = 1.606RMSE0.0220.0220.0140.029b = 2.171H0 = 3 minsR-square0.9020.9040.9620.867a = 2.669RMSE0.0160.0160.0100.019b = 2.740H0 = 4 minsR-square0.8790.8770.9690.757a = 3.018RMSE0.0160.0160.0080.022b = 3.861H0 = 5 minsR-square0.8880.8700.9680.836a = 2.908RMSE0.0160.0170.0080.019b = 4.795H0 = 6 minsR-square0.9120.8790.9890.851a = 2.770RMSE0.0140.0170.0050.018b = 6.021H0 = 7 minsR-square0.9170.8830.9890.856a = 2.626RMSE0.0140.0180.0050.019b = 7.117H0 = 8 minsR-square0.8980.8620.9890.819a = 3.080RMSE0.0140.0160.0060.018b = 8.121H0 = 9 minsR-square0.8590.8190.9760.894a = 2.903RMSE0.0170.0190.0070.013b = 8.849H0 = 10 minsR-square0.9170.8900.9870.919a = 3.185RMSE0.0120.0130.0040.008b = 10.16Following the above fitting analysis, the expectation and standard deviation of passenger IVT and WT in Zhongguancun Software Park is estimated and analyzing result is shown in Table7.Table 7
Expectation and standard deviation of IVT and WT.
LineIVTWTE (min/km)SD (min/km)E (min)SD (min)9023.73.32.45.73203.72.45.15.63333.72.43.65.63622.92.03.65.65093.32.25.16.85213.72.74.16.45703.01.95.16.89633.82.62.45.76363.22.24.46.2853.01.95.16.8823.22.25.16.8


## 4.3. Results and Discussion

CRnetwork is introduced to estimate the overall reliability and quantify the improving effect. CRnetwork represents the sum of extra time cost due to the unreliability of the travel time for all passengers in the network, which can be calculated as follows:(27)CRnetwork=∑o∑d∑suρσs+DsYs,∀o,d∈V,∀s∈Sin−vehicle∪​Swaiting.Based on the optimization result, in local bus network of Zhongguancun software, stop Metro Xierqi to stop Dongbeiwang is the hot OD pair with insufficient reliability during peak hours. The congested travel demand can be diverted by shuttle service. The detailed plan is described in Table8. The number of passengers by shuttle is 240 person/h, which takes account of 12.2% of total transit demand. The running route of the shuttle is: stop Metro Xierqi-Road Ruanjianyuan-stop Dongbeiwang. For the original local bus network, CRnetwork is 15650.8 yuan. After the operating of shuttle, CRnetwork is 13524.9 yuan, which has a significant decrease by 14%.Table 8
Optimization results.
Last-mile shuttle routeDiverting demandOptimal cost of travel time reliabilityStop metro xierqi—stop road ruanjianyuan-dongbeiwang240 person/hOrigin conventional bus networkAfter shuttle operatedOptimization effect15651 yuan13525 yuan14%In the area of Zhongguancun Software Park, the spatial and temporal aggregation of transit demand is significant. Most demand occurrence points (point O) are located at the stop Metro Xierqi station, and the destination points are concentrated in several key stops near the office area. Meanwhile, commuting demand is concentrated in peak hours, and the roads are of low grade, which is narrow and not suitable for oversize vehicles. Although the working area in the study is small, the results show that the shuttle service can significantly improve the travel time reliability of last-mile trips. The applicable scenarios can be expanded to more similar areas such as residential areas or school areas.


## 5. Conclusions

This paper proposes a methodology of last-mile shuttle service for a commercial or residential area. Shuttle route is designed to relieve the congested demand between the unreliable OD and improve commuters’ time reliability for the last-mile trip. The main contributions of this study are as follows:(1)
An offline method for estimating reliability is proposed, including data processing methods and statical fitting models. Buffer time is used to quantify travel time reliability considering IVT and WT. With more than 5 hundred thousand pieces of real-world GPS data, fitting analysis is conducted. The results indicate that lognormal distribution and Laplace distribution are the best fit for the distribution of sample IVT and bus arrival time. Key parameters to reflect reliability, including IVT reliability and WT reliability, are estimated based on the fitting results.(2)
Shuttle demand is regarded as a dynamic decision variable in the bilevel optimization model in this study, which can be optimized with the shuttle route jointly. In order to solve the mixed integer nonlinear programming problem, a two-stage solution algorithm is proposed combined with the linearized solution method. The shuttle route and the amount of demand which are diverted from the conventional transit network are estimated. Reliability cost is defined as the key index to evaluate passenger travel time reliability condition.(3)
To illustrate the robustness and quality of the model, the real-world case study is conducted in Zhongguancun Software Park in Beijing. It is indicated that the last-mile shuttle service can significantly improve transit passenger travel time reliability by 14%. In addition, 12.2% of total transit passengers choose the new shuttle service in the case study. The findings can be promoted that in the area of the core business district or residential community, when the commuting demand is concentrated, the shuttle service will be more effective in solving the poor reliability of the last-mile trip and the optimization result will be better.The paper has limitations and more researches can be improved in the future. The last-mile shuttle service is a new transit mode in Beijing. Therefore, actual operation data are lacking, which are replaced by that of the microcirculation bus in the study. The effect can be better verified when the actual service is provided in the future. In addition, the influence of shared-bike and carpool is not considered. Cooperation of multitransportation modes can be analyzed to improve travel time reliability of last-mile trips.



---

*Source: 5117488-2022-01-18.xml*