Designing zonal-based flexible bus services under stochastic demand
11 DESIGNING ZONAL-BASED FLEXIBLE BUS SERVICES UNDER STOCHASTIC DEMAND
Enoch LEE a , Xuekai CEN* b , Hong K. LO a , and Ka Fai Ng a a Department of Civil and Environmental Engineering, The Hong Kong University of Science and Technology, Hong Kong, China b School of Traffic & Transportation Engineering, Central South University, China * Corresponding author Email addresses: [email protected] (Enoch LEE), [email protected] (Hong K. LO), [email protected] (Xuekai CEN), [email protected] (Ka Fai Ng)
ABSTRACT
In this paper, we develop a zonal-based flexible bus services (ZBFBS) by considering both passenger demands’ spatial (origin-destination or OD) and volume stochastic variations. Service requests are grouped by zonal OD pairs and number of passengers per request, and aggregated into demand categories which follow certain probability distributions. A two-stage stochastic program is formulated to minimize the expected operating cost of ZBFBS, in which the zonal visit sequences of vehicles are determined in Stage-1, whereas in Stage-2, service requests are assigned to either regular routes determined in Stage-1 or ad hoc services that incur additional costs. Demand volume reliability and detour time reliability are introduced to ensure quality of the services and separate the problem into two phases for efficient solutions. In phase-1, given the reliability requirements, we minimize the cost of operating the regular services. In phase-2, we optimize the passenger assignment to vehicles to minimize the expected ad hoc service cost. The reliabilities are then optimized by a gradient-based approach to minimize the sum of the regular service operating cost and expected ad hoc service cost. We conduct numerical studies on vehicle capacity, detour time limit and demand volume to demonstrate the potential of ZBFBS, and apply the model to Chengdu, China, based on real data to illustrate its applicability. Keywords: flexible bus, demand responsive transit, stochastic programming, reliability INTRODUCTION
Public transport (PT), if successful, is instrumental in addressing the transportation needs of many metropolitan areas. Traditionally, without real-time demand information, PT services are operated on fixed routes with fixed schedules (FRFS), regardless of the actual demand, which are neither efficient nor cost-effective. As for demand responsive services, such as taxi, while flexible, they are too costly to serve as a major PT mode due to lack of economies of scale. With the emergence of information technology, travelers can request for ride services in real-time. Transportation network companies (TNCs), such as Lyft, Uber, and Didi Chuxing, can match travelers to vehicles almost instantaneously to provide for ride services, typically through small vehicles that give high flexibility but may worsen traffic congestion in metropolitan cities (Schaller, 2018). TNCs open up the possibility of demand responsive ride sharing services that can be operated at lower costs due to economies of scale, and yet yield much more flexibility than FRFS services, e.g., Didi Bus. The key is to devise a system to allow passengers sharing the same vehicles, which will achieve economies of scale, and at the same time restrict the detour time to pick up multiple passengers to an acceptable level. The concept of demand responsive ride sharing services is not new, which is also referred to as dial-a-ride (DAR) services where customers dial calls for service requests, and the associated problem is named as the dial-a-ride problem (DARP). Such DAR services are available today but limited in scopes, such as airport trips, or trips for the elderly, people with disabilities or special needs. What is new now is that such service requests, enabled by service platforms such as Uber and Didi, are gaining popularity to the extent that demand responsive transit (DRT) may become a substantive PT mode, especially in view of the widely recognized Mohring effect: increasing passenger volumes allows for improved services, initiating a virtuous cycle further attracting new users. In fact, such demand responsive services are gaining recognition in China, referred to as customized bus service, serving people between central business area and residential area especially in peak hours (Liu and Ceder, 2015). DARPs are mainly about designing vehicle routes and schedules for service requests between origins and destinations, and can be considered as a generalization of the pick-up and delivery vehicle routing problem (PDVRP) and the vehicle routing problem with time windows (VRPTW) (Cordeau and Laporte, 2007). Important features of DARPs include vehicle capacity and time window or detour time for picking up and dropping off passengers in a reasonable limit. Ho et al. (2018) divided DARPs into deterministic and stochastic DARPs, which depend on whether the information (such as demand) is known with certainty or not, and stated that few studies focused on static and stochastic DARPs. In this paper, we aim to develop a novel zonal-based flexible bus system, which is a demand responsive transit service considering stochasticity of zonal demand volume and detour time within each zone, vehicle capacity, and detour constraint for each bus. Flexible bus is a class of DRT that provides door-to-door service with a relatively low price, which can bring a higher social welfare to an area with high population density (Zhang et al., 2018). Some researchers studied DRT with analytical models, which typically assumed that the demand is uniformly distributed over space (e.g., Chen and Nie, 2017a, 2017b; Diana et al., 2006; Kim and Schonfeld, 2014, 2015; Quadrifoglio and Li, 2009; Zhang et al., 2017; Daganzo, 1978). Such assumptions of demand distribution are essential for the analytical models, even though they may not reflect reality well, since actual demands tend to aggregate at specific locations, e.g., train stations and shopping centers. To relax this assumption, some researchers utilized historical or actual demand to design flexible bus service (e.g., Eiró et al., 2011; Ma et al., 2017a, 2017b; Martínez et al., 2015). Amirgholy and Gonzales (2016) employed an analytical model to approximate the operating cost for running a DRT system with dynamic uniform distribution of demand. Zhang et al. (2018) developed an analytical model to investigate flexible bus service for a linear corridor connecting city center and suburb residential area. Tong et al. (2017) optimized customized bus service by considering time-space lim its of the requests’ origins and destinations, and solved the associated pickup and delivery vehicle routing problem with time windows (PDVRPTW) by combining various heuristics methods. R. Guo et al. (2018) solved a customized route design problem that characterized the passengers-to-vehicle assignment and compared the solutions from genetic algorithm and branch-and-cut algorithm. As an on-demand service, the flexible bus service system should also consider the dynamics and stochasticity of demand, rather than only deterministic and historical demand. Santos and Xavier (2015) solved the dynamic DARP by matching the taxi and ride sharing service with dynamic passenger demand, while ensuring that the constraints on vehicle capacity, maximum trip cost of each passenger, and maximum delay are satisfied. Ho and Haugland (2011) endogenously considered stochastic demand and formulated the associated DARP as a two-stage stochastic integer linear program with recourse, in which a priori set of vehicle routes was determined in stage 1 to accommodate a certain set of service requests, and in stage 2, the vehicle routes determined in stage 1 were modified by skipping unrealized requests. Hyytiä et al. (2010) devised a dial-a-ride model with stochastic demand arrival and conducted analytical analyses. Q. W. Guo et al. (2018) studied stochastic dynamic switching between fixed and flexible transit services as the demand density evolves continuously over time. The above studies, however, only considered stochasticity of demand volume and neglected the fact that the service requests are spatially stochastic and travel time varies. As flexible bus services detour to serve each service request door-to-door, the spatial stochasticity of demand and variability of travel time should not be neglected. The spatial stochasticity of demand is dependent on the clustering strategies, i.e., how to group neighborhood service requests into different zones and design the zonal flexible bus routes. In fact, a number of mathematical models were proposed to optimize flexible transit service by implementing a zonal strategy (e.g. Aldaihani et al., 2004; Li and Quadrifoglio, 2009; Pan et al., 2015; Quadrifoglio et al., 2006). Chang and Schonfeld (1991) developed optimization models for zonal service operation, which included service zone area and bus capacity as decision variables. Karabuk (2009) introduced a cluster-first route-second zonal heuristic that clusters customers to receive service together and arranges vehicle routes based on the zonal clustering. Recently, there were some studies of DRT that divided the service zone before planning, such as Lu et al. (2017) and Wang et al. (2018). In the former paper, each zone is served by one operator, and vehicles cannot traverse boundaries unless they pick up or drop off customers. In the latter research, a zone is served by only one bus, and passengers are picked up at a common terminal. By this decentralization approach, the problems were divided into several smaller sub-problems, yielding higher efficiency. However, to our best knowledge, none of these studies on DRT have jointly addressed the two aspects of demand stochasticity, namely, volume and spatial stochastic variations. This paper aims to develop a zonal-based flexible bus system by considering both volume and spatial stochasticities of demand. The spatial stochasticity is represented by the assumption that the zonal detour times for picking up and dropping off passengers in each zone follow certain probability distributions that can be inferred from historical service requests. The detour time probability distribution inherits the variation of travel time throughout the day. Similarly, the numbers of service requests also follow certain probability distributions to reflect the volume stochasticity. The frequency-based zonal flexible bus is scheduled to minimize the operating cost, including the total regular service cost and ad hoc service cost, by considering the vehicle capacity, detour constraints within each zone, and spatial and volume stochastic variations. Naturally, it can be formulated as a two-stage stochastic problem with recourse, in which stage 1 determines the zonal visit sequences of the bus fleet and stage 2 assigns buses or ad hoc services to serve the realized service requests. To solve this intrinsically intractable problem, we apply a service reliability (SR)-based gradient approach, which separates the original problem into two phases for solution efficiency and has been adopted for designing ferry service (An and Lo, 2014) and transit network (An and Lo, 2016), managing traffic signals with stochastic demand (Li et al., 2018), and controlling network signal with equilibrium constraint (Huang et al., 2018). Reliabilities on both volume and spatial stochasticity are introduced and optimized. The volume reliability is the probability that all service requests are picked up by the bus fleet, while the detour reliability defines the probability of the time required for a vehicle to pick up passengers. Basically, this approach consists of two phases. In phase-1, given a specific service reliability, the number of requests to be served and the respective detour times are determined and fixed based on the historical probability distributions. The problem essentially becomes dete rministic. We will then determine the fleet and vehicles’ sequences of zonal visits to minimize the operating cost so as to serve the service requests. In phase-2, with the fleet and vehicles’ sequences of zonal visits determined in phase-1, service requests will be assigned to vehicles subject to the detour time and vehicle capacity constraints. If certain requests are not accommodated by the regular fleet due to vehicle capacity or detour time constraints, those requests will be served by ad hoc service, which can also be considered as the cost of lost or unserved demand. The sum of phase-1 operating cost and expected phase-2 ad hoc service cost will be minimized via varying the SR by the gradient based approach. The feasible region of the two-stage stochastic problem with recourse is proven to be equivalent to the feasible region of the proposed SR-based gradient approach. In other words, any optimal solution in the original two-stage stochastic problem with recourse can be found by searching the feasible region of the SR-based formulation. This system optimizes the operating cost for operators, and the abovementioned flexible bus system is suitable for operators that provide flexible bus and taxi service, such as Didi, Uber or Grab. The contributions of this research include (1) the concept of zonal network with intra-zonal detour is introduced and the sequences of zonal visits are optimized; (2) the formulation incorporates heterogeneous service requests, and considers the volume and spatial stochastic variations of the service requests, as well as the variation of travel time endogenously; (3) the service reliability (SR)-based gradient approach is introduced and applied to optimize the sequences of zonal visits; (4) the numerical results show that medium-size vehicles are preferred to minimize the cost while preserving the quality of service in detour time.; and (5) the formulation and solution algorithm are applied to Chengdu, China, based on real data to illustrate their applicablility. The remainder of this paper is organized as follows. In Section 2, the two-stage stochastic formulation of the zonal-based flexible bus service is presented, and the SR-based gradient solution approach is introduced. After that, numerical examples are conducted, and policy implications are presented in Section 3. Section 4 draws the conclusion and provides future research directions. Model Formulation
The following notation is utilized and summarized for ease of reference.
Sets Z Set of zones V Set of vehicles v Z Set of zones traversed by vehicle v RS v Z Set of zones from zone R to zone S traverses by vehicle v P Set of candidate routes E Set of demand categories Set of demand scenarios Set of OD pairs to be served D / D Set of service requests / Set of service requests in scenario Parameters z t Maximum allowable detour time in zone z , z Z t Maximum allowable detour time for the whole trip RS t Maximum allowable detour time from zone R to S, (R,S) z T A matrix of detour time of the vehicle in zone z that considers correlations. cap Vehicle capacity p c Operating cost of route p , p P p m The number of zones traversed by route p ijp Operating cost to provide space for a passenger from zone i to zone j by vehicle of route p , ( , ) , p Pi j d c The ad hoc service cost of serving request d , d D p Probability of scenario , zd Detour time required in zone z to serve service request d , d D , z Z zdb − Half of detour time reduction in zone z if service request d and b are served simultaneously zd + Binary parameter indicating service request d starts in zone z , d D , z Z zd − Binary parameter indicating service request d ends in zone z , d D , z Z d n Number of passengers of service request d , d D ( ,..., ) p pn z z Zonal visit sequence for route p , ,., .., p pn P z Zp z { } prszp b = B A converting matrix to express the number of in-vehicle passengers when leaving zone z for route p , where rs are possible OD pairs, ,( , ,) Z p Pr s z e Random variable of number of service requests of demand category e , e E z Random variable of additional detour time to serve one more service request in zone , z z Z z Detour time required from the boundary to the first point or the last point in zone z II z Detour time per service request served in zone , z z Z e n Number of passengers of demand category e , e E ze + Binary parameter indicating demand category e starts in zone z , e E , z Z ze − Binary parameter indicating demand category e ends in zone z , e E , z Z , M M
Very large positive numbers
Decision variables { } pv x = X Binary variables indicating whether vehicle v is assigned to route p , , p P v V { } p = Χ' Number of vehicles assigned to route p , p P { } dv w = W Binary variables indicating whether service request d is assigned to vehicle v , , d D v V v w The vector of dv w for all service requests for vehicle v RS, { } v v = ζ Number of passengers of OD pair RS picked up by vehicle v in scenario , RS, : (R,S) v v = ζ , ,( ,S) ,R v V RS p Number of passengers from OD pair RS picked up by vehicle of route p , (R,S) , p P { } ev y = Y Number of service request of type e served by vehicle v , , e E v V { } ep y = Y Number of service request of type e served by vehicle with route p , , e E p P vz y Number of pickup or drop off point of vehicle v in zone z , , v V z Z I I { } e = ρ Volume reliability, probability that all requests of service requests of type e are picked up by the fleets, e E II II { } z + = ρ Detour reliability for service request in zone z . The II z -th percentile of detour time is considered in phase-1 problem, z Z e Number of service requests of type e to be picked up in phase-1, e E f C The fixed operating cost determined in stage 1 and phase-1 Q The expected ad hoc service cost t o t a l C The total cost, including of operating cost and ad hoc cost
Zonal-based flexible bus service 2.1.1
Definitions and characteristics of zone, service request and route
The service area is divided into zones geographically, as illustrated in Figure 1, which is a zonal-based network constructed with zone centroids, represented as nodes and roads connecting zone centroids as links. The zone shapes and sizes are considered as given. Indeed, the division of zones, which will affect the planning of the flexible bus service, can be optimized based on the historical demand patterns, and is our current research. The intra-zone service requests are omitted, as ZBFBS is intended for service requests whose origins and destinations fall into different zones. Service requests are characterized by origin-destination (OD) zones, detour times, and the number of passengers per request. The calculation of detour time, considering the correlation between service requests, is explained in Section 2.3. There are two definitions of routes; first, route p in this paper refers to a zonal route , defined as a sequence of zonal visits, and belongs to the candidate route set P , i.e., p P . Each zonal route has a cost p c that is dependent on the length of the route. The ZBFBS problem is generic for any set of zonal routes. Second, an actual route , which specifies the route a flexible bus actually travels including the detours measured in Euclidean distance for simplicity, which can be converted back to street network distance via certain multipliers (Yang et al., 2018). Similar to door-to-door services around the globe, drivers are flexible in determining the actual route based on the traffic condition while visiting the service requests assigned. We minimize the operator’s cost of the ZBFBS by determining the route of every vehicle and assigning service requests to vehicles while ensuring the vehicle capacity constraint and detour constraint for each vehicle are satisfied. Figure 1. An illustration of zonal-based flexible bus service system Assumptions
The following assumptions are made in this paper: A1. The marginal detour time decreases with the number of service requests served in the zone. A2. Door-to-door services are provided without transfers. A3. One vehicle type is used for the regular flexible bus service. A4. Each OD pair has at least one route traversing it. A5. The intra-zone service requests are omitted.
Two-stage stochastic formulation
We aim to design a flexible zonal-based bus system tackling while addressing the uncertainties of demand over space and volume, which follow certain probability distributions. Two types of services are provided: regular service operated with designated frequency, and ad hoc service to accommodate the demand unserved by the regular service. The stochastic formulation is divided into two stages. The first stage decides the regular service provided, where vehicles are either assigned to different routes or idle. Based on the routes determined in the first stage, the second stage matches the realized service requests with the regular service vehicles. Service requests must be served by either regular service or ad hoc service, in which the latter incurs an additional cost to the system. Note the ad hoc service provides one smaller vehicle per service request, with a predefined cost per service request. Regardless of the number of passengers to be carried per request, the entire group of passengers of each request must be served by the same vehicle. The two-stage stochastic problem with recourse (P0) can be stated as,
Stage-1 : min ( ) ppVv pvP c x Q + X,W X (1) subject to pvp P Vx v (2) {0,1}, , pv p P v Vx (3) where, Stage-2 : ( ) min [ (1 )] dDd Vd vv Q p c w = − W X (4) subject to , ,, zv vd zD t z Z v V w T w (5) SRS, R (R,, S) , , v d dD d dvd v Vn w + − = (6) ,, , p pv p v m x cap v V p P − B ζ 1 (7) RS, 1 , (R,S) , , v pvp P
M x v V (8) ,1, V dvv dw D (9) dv w d D v V (10) The goal of this problem is to minimize the expected total system cost through scheduling the regular flexible bus services and introducing ad hoc services if needed in order to accommodate the realized demand. Stage 1 determines the regular zonal-based flexible bus routing. The objective function (1) minimizes the operating cost of the regular flexible bus service and expected Stage 2 cost, where p c is the operating cost of route p P for a vehicle and P is the set of candidate routes. pv x is the binary decision variable which equals 1 if vehicle v V is assigned to route p , otherwise 0, forming the binary matrix { } pv x = X . V is the set of vehicles. Vehicle can only operate on at most one route, as represented in constraint (2). Stage 2 minimizes the expectation of the ad hoc services cost by allocating the realized demand to regular flexible buses with route determined by X in stage 1. Following Lo et al. (2013), we directly formulate the problem in the context of scenario simulation. Each scenario has a probability of p with p = , and a service request set D . { } dv w = W is the binary decision variable which equals 1 if service request d is served by vehicle v ; 0 otherwise. If a service request d is not served, i.e., dvv V w = , an ad hoc service cost d c will be incurred, as denoted in objective function (4) in Stage 2. v w is the vector form of dv w . Constraint (5) enforces that the total detour time of each vehicle in any zone z should not exceed the time limit z t . v w is the transpose of vector v w . z T is a D D square matrix to process the detour time considering the correlation between the service requests, a detailed explanation is in the next subsection. Equation (6) calculates RS, v , the number of passengers from OD pair (R,S) served by vehicle v in scenario , where d n is the number of passengers of request d . R d + and S d − indicate request d starts from zone R and ends in zone S, respectively. is the set of OD pairs to be served. Constraint (7) denotes the vehicle capacity constraint, which ensures that each vehicle v will not exceed its maximum capacity cap over the entire trip. To ensure that vehicles are within the capacity limit in the zone, vehicles first drop off some in-vehicle passengers before picking up new passengers. Vehicle capacity constraints are set for all the zones they traverse except for the end zones where no new passengers are picked up, i.e., vehicle v should have p m − zonal capacity constraints, where p m is the number of zones in route p . v ζ is the vector of RS, v , and RS { } pzp b = B is a ( 1) p m − matrix converting RS, v to the number of in-vehicle passengers when vehicle v leaves each zone, where is the number of OD pairs. p m − is the column vector with p m − items equal to 1, indicating the p m − zonal capacity constraints. To construct the matrix p B , we first rule out all the OD pairs that are not on route p and set all the corresponding RS { } pz b equal to a large number M , ensuring that the service requests from OD pair RS cannot be served by any vehicle on route p not traversing OD pair RS. If zone z is on the path from R to S, then RS { } pz b is equal to 1, meaning that any passengers from R to S are onboard when vehicle v leaves zone z , and 0 otherwise. This is inspired by the capacity constraints in the classical pickup and delivery problems (Savelsbergh and Sol, 1995). Constraint (8) ensures every vehicle serving passengers has a route. Constraint (9) ensures a request to be served by at most one flexible bus. Constraint (10) defines dv w W to be binary, i.e., the requests have to be served entirely regardless of the number of passengers. An illustrative example of the detour constraint, converting matrix, and vehicle capacity constraint is in Appendix A. Note that although the detour time constraint (5) restricts the total detour time per zone, it can be replaced or supplemented by the following constraints that restrict the detour time per trip and/or per OD pair, respectively. ,, v zZ v vz d D t v V w T w (11) RS RS ,, (R,S) , v zv vz dZ D t v V w T w (12) v Z represents all the zones traversed by vehicle v , and RS v Z represents the pair from zone R to zone S traversed by vehicle v . After determining X in stage 1, v Z and RS v Z can be determined from the vehicle route. Detour time quantification
Detour is quantified by the time required to travel between the corresponding zone centroid and the pickup or drop off location, which implicitly captures the traffic conditions. After the demand is realized, the flexible bus does not need to go back to the zone centroid after each pick-up or drop-off, but goes from one service request location to another, hence the detour time between service requests depends on the proximity between their service locations. In general, the total detour time within a zone is shorter than the sum of individual detour times of serving multiple requests (as measured from the individual request location to the zone centroid). This reduction in the total detour time within a zone is determined from the proximity between the pickup or drop-off locations. 0 To facilitate the calculation of detour time, after the realization of service requests,
D D symmetric detour time matrices z T are introduced for every zone z . Specifically, the d th element in the diagonal of z T is the detour time required in zone z to pick-up or drop-off request d (as measured from its location to the zone centroid), denoted by zd , while the other elements are the reduction of detour time when two nearby service requests are served by same vehicle (without having to travel back to the zone centroid for each request), denoted by non-positive zdb − at entry ( , ) d b . Therefore, z zdb bd − − = . To ensure that the total detour time within a zone always increases when an additional service request is served, the sum of potential reductions of detour time for service request d in zone z must be smaller than the detour time of d at z , formally expressed as, T z zdb d z zd db d D z Z − − − (13) where cap is the vehicle capacity, and zd is a set of the cap -smallest zdb − of d for all b D , includes the smallest and the next ( 1) cap − -smallest values of zdb − . If there are less than cap service requests with either origin or destination in zone z , all the negative zdb − are included in the set T zd . After the service requests are realized, zd can be obtained from the historical travel time between zone centroid and location of the service request while zdb − can be calculated based on the proximity of the locations between the service requests. Note that the detour time is often overestimated as the condition (13) limits the reduction of the detour time that may be actually higher, especially when there are a lot of service requests in close proximity within a zone. An example of the detour time matrix is given in Appendix A. Calculations of detour time and detour time reduction from real data are demonstrated in Appendix G. Service reliability-based formulation and solution algorithm P0 is a stochastic binary program that is hard to solve, especially for large numbers of demand requests or a large vehicle fleet. For the sake of computational efficiency, we decompose the problem by introducing service reliabilities, similar to An and Lo (2014) and Lo et al. (2013). The random demand volume is separated into two parts: the first part is based on the volume reliability I ρ , which will be served by regular services; the other part is the demand above that level, which will be served by ad hoc services if needed. Similarly, the random detour time is separated by the detour time reliability II ρ , which is considered when planning the regular services. Upon realization of the actual demand, if the total vehicle detour time exceeds the maximum allowable detour time, some service requests will be served by ad hoc services. The volume and detour time reliabilities separate the original problem P0 into two sub-problems P1 and P2 , corresponding to the phase-1 problem formulated under specified reliability measures, and the phase-2 problem formulated to assign requests to vehicles upon realization of the actual demands. P1 tackles the aggregate service requests, defined by demand categories, while P2 deals with individual service requests. In phase-1, we determine the schedule of regular zonal flexible bus services to cover the stochastic demand up to the 1 specified volume reliability and detour time reliability. The number of passengers onboard should not exceed the vehicle capacity cap in any zone, and the total detour time per vehicle should not exceed the maximum allowable detour time I t designated for every zone I. In phase-2, we determine the allocation of the realized service requests to the existing flexible bus determined in phase-1 and the deployment of ad hoc services. In this problem setting, the formulation is more complex than the previous studies of An and Lo (2014) and Lo et al. (2013), as it not only considers two types of stochastic variations in phase-1, namely volume and spatial stochasticities, but it also matches the service requests with specific vehicles in phase-2, while the matching process was not considered in An and Lo (2014) and Lo et al. (2013). The key to solve the problem is to determine the optimal values of the two demand reliabilities, i.e., volume reliability I ρ representing volume stochasticity and detour time reliability II ρ representing spatial stochasticity, such that the expected total operating cost is minimized. In addition, demand categories e E are introduced to represent the aggregate historical record of service requests d according to OD pairs and passenger groupings d n . A demand category e is specified by its origin and destination, number of passengers per service request de n n = and volume stochasticity given by the random variable of number of requests e . Detour time approximation in phase-1
In the phase-1 problem, the calculation of detour time is separated into two parts as illustrated by arrows with different colors in Figure 2. The first part, indicated by the green arrow, is function z that estimates the detour time between the boundary of zone z and the closest service location, which is assumed to be convex decreasing with respect to vz y , i.e., the number of service requests served by vehicle v in zone z . The function is decreasing because it is more likely to have some service locations closer to the zone boundary with more service requests, hence a shorter detour. By equation (5) of P0 , with more service requests being served, the accumulated reduction of detour time increases. Similarly, function z captures the reduction of detour time with the number of service requests served. The total detour time is shorter than that calculated without correlation in phase-1. The second part is the detour times between locations of service requests, given by the product of II z and the number of trip segments, as indicated by the orange arrows between two points. II z is the detour time between two points specified by the detour time reliability. Similar to the diagonal values of z T in equation (5), this part calculates the detour time per additional service request served. Moreover, the value of II z accounts for the proximity of the location of service requests being served, which is captured by detour time reductions in (5). For example, service request locations that are close will result in a lower II z in phase-1, which is equivalent to higher detour time reduction in z T . 2 Figure 2. An example route of vehicle in a zone By considering the detour time between the boundary and the closest service location and the detour time between the locations of two service requests, vz t , the detour time of vehicle v in zone z is approximated as, ( ) ( ) II ,2 1 , z vz vzvz z y v V z Zt y = + − (14) where vz y , the number of service requests served by vehicle v in zone z , is defined as, ( ) , , z zvz e e eEe v y v V z Zy + − + = (15) ev y is the number of service requests of demand category e assigned to vehicle v , which is one of the decision variables in the phase-1 problem. In the two-stage formulation, the detour times zd , diagonal elements of z T , capture the additional detour times for service requests, which are approximated in phase-1 by II z based on the detour time reliability in each zone. The detour time reductions zdb − that capture the reductions of marginal detour times for the number of service requests served and proximity between the service requests are approximated in phase-1 by II z and z , respectively. Hence, this approximation of detour time is consistent with the two-stage formulation. Formulation of the phase-1 problem (P1)
The phase-1 problem P1 can be formulated as, I II 1, )min ( , f P pv pp vV
C c x = X Y ρ ρ (16) subject to (14), (15) pvp P Vx v (17) I rinf ,| P ( ) e ee e E = (18) , ev ev V Ey e = (19) 3 R SRS ( ,, R,S) v E e e e eve v Vn y + − = (20) , , p pv p v m cap v V p Px − ζ (21) IIII | Pr( )inf , z z z z Z = (22) , , vz z t v V z Zt (23) , v pvp P M x v V (24) ,1 , ,0 pv Vx p P v (25) Y and integral (26) The Phase-1 problem determines the regular flexible bus service to carry the passenger demand up the levels as specified by the volume and detour time reliability measures. Specifically, { } pv x = X , the assignment of vehicle v to route p , and { } ev y = Y , the number of requests of demand category e assigned to vehicle v , are optimized by fixing the volume reliability I ρ in Equation (18) and detour time reliability II ρ in (22). Equation (18) calculates the number of service requests e which are to be served by the regular flexible bus service according to the volume reliability measure I e and the random variable e for each demand category e E . Similar to (18), (22) defines the detour time between two points II z in each zone z without requiring II z to be integers. z is the detour time distribution between two service points in the zone. The objective function (16) minimizes the total regular service operating cost. Constraints (20), (21) and (23) are similar to (5), (6) and (7), which denote the vehicle capacity constraints and zonal detour constraints with demand category e rather than service request d as expressed in (5), (6) and (7). To ensure the service requests can be served, it is reasonably assumed that z z t as the operator should allow enough detour time per zone to pick up or drop off at least one service request. Constraints (17) and (25) are identical with constraints (2) and (3), indicating vehicle v can only be assigned at most one route p P . Equation (19) enforces that all the requests generated from (18) should be served, while they can be served separately by any vehicle v . Constraint (24) ensures that only vehicle assigned to a route, i.e. pvp P x = , can serve customers. M is a large positive number. Constraint (26) is similar to (10) that ensures every request is served entirely. The LHS of constraint (21) is bilinear, so (21) can be reformulated using the linearization technique (e.g., Cen et al., 2018) as follows. ) ,(1 , p p p v m pv m M p P v Vcap M x cap − − + − B ζ 1 1 (27) When pv x = , i.e., vehicle v is not assigned to route p , (27) will always hold. If pv x = , we have p p v m cap − ζ 1 B , which is the designated capacity constraint. M is a positive number 4 much larger than M used in p B . Furthermore, as ( ) z vz y is a convex decreasing function and vz y is a natural number, (14) can be linearized by equation (28) via replacing ( ) z vz y by several linear constraints (Ng et al., 2020). ( 1) ( )) ) ( }( ( , , V, {0, 2) , 4,..., 2 z z z vz z i zyi i i Z v i cap + − − + (28) Through these adjustments, P1 becomes a mixed-integer linear programming (MILP) problem that can be solved by commercial solvers such as IBM CPLEX. Similar to (11) and (12), the detour time constraint (23) can be replaced or supplemented by the following constraints to account for the maximum allowable detour time in a trip, , v vzz Z Vt t v (29) RS RS , (R,S) , v vzz Z t v Vt (30) Phase-2 problem (P2) min (1 ) dvVd D d v Q wc = − W (31) , , v vz zDd t z Z v V w T w (32) R SRS, , (R,S) , v d d dvdd D v Vn w + − = (33) ζ , , p pv p v m cap v V p Px − (34) RS, 1 , (R,S) , v pvp P
M x v V (35) V dvv d Dw (36) dv w d D v V (37) The phase-2 problem is close to the Stage-2 problem in P0 , except that the regular flexible bus deployment { } pv x = X is already fixed after solving the phase-1 problem; hence constraint (34) is now linear in v ζ . The service requests set D in phase-2 is then generated from the random variable of the number of requests e across every demand category e E . Equations (32) to (37) are identical to equations (5) to (10) with only one scenario. (11) and (12) can also be used here if appropriate. Note that the calculation of detour time explained in Section 2.3 also applies here, as the locations of service requests are known in this stage. Based on the vehicle route assignment X solved in P1 , the set of realized service requests 5 D is to be assigned to vehicles, or an ad hoc service cost is incurred. It is a special case of the Stage-2 in P0 with one scenario. P2 is solved for each scenario generated, denote I II , )( = ρ ρ ρ the expected ad hoc cost is given by the following function of ρ , ( ) Q p Q = ρ (38) Finally, the total cost is the sum of fixed cost and ad hoc cost, ( ) ( ) )( total f C C Q = + ρ ρ ρ (39) With fixed volume reliability I ρ and detour time reliability II ρ , f C denotes the regular flexible bus operating cost, which is the objective value in P1 , while Q is the expected ad hoc cost solved in P2 . The objective is to find the optimal reliability I II ( , ) ρ ρ that minimizes the total cost. The basic idea to solve the problem, as will be discussed later, is to develop an SR-based gradient procedure over
I II ( , ) ρ ρ to solve P1 and then repeatedly solve P2 until a certain stopping criterion is reached. Interchangeability between the two-stage and SR-based formulations
To ensure the equivalence between the two-stage and SR-based stochastic formulations, the following additional assumptions are made. It is not necessary to have these assumptions to solve the problem, but the interchangeability between these two formulations is not guaranteed without them. B1. The possible route set P contains one and only one route from each OD pair in , which is taken to be the shortest path. B2. The route operating cost is additive, i.e., c c c = + (cost from zone A to zone B to zone C is equivalent to the sum of cost from zone A to zone B and zone B to zone C). B3. Each OD pair in has at least one one-person service request. The SR-based formulation P1-P2 is introduced by cutting the two-stage stochastic problem with recourse P0 by the reliability measures ρ . The interchangeability of P0 and the SR-based formulations is shown by the fact that any feasible solution of P1 - P2 is feasible in P0 ; and any feasible solution of P0 is feasible in P1-P2 by adjusting the reliabilities. Propositions 1 and 2 below prove the result for the case of one demand category per OD, and Proposition 3 further generalizes the interchangeability to the case with multiple demand categories per OD.
Proposition 1.
Any feasible solution of
P1-P2 is also feasible to the original two-stage stochastic formulation P0 . Proof.
See Appendix B. □ Proposition 2 . For problems with the stochastic volume described by a probability distribution for any integer greater than 0, and the stochastic detour time described by a probability distribution over the range [0, ) + , if there exists only one demand category with one passenger per OD, any optimal solution of the original 2-stage stochastic formulation can be generated by the SR-based stochastic formulation. Proof.
See Appendix C. □ Proposition 3.
For problems similar to those stated in Proposition 2 but with multiple demand categories per OD, including single passenger service requests, any optimal solution of the original 2-stage stochastic formulation can be generated by the SR-based stochastic formulation.
Proof.
See Appendix C. □ SR-based solution procedure
In the above subsections, the two-phase stochastic program P1 - P2 was introduced to solve the original problem P0 . According to Proposition 1 and Proposition 3, the same global optimal solution to the original two-stage problem can be found by searching the entire feasible region of the SR-based formulation. Nevertheless, given that the problem is non-convex, this may not be practical except for very small problems. The proposed solution procedure here aims to find the descent direction in each step, thus decreases the total cost gradually to achieve a local optimal solution of reliability
I II , )( = ρ ρ ρ . The ability to find local optimal solution, when the derivative of the objective function with respect to the reliability is zero, distinguishes our algorithm from heuristic methods in which solution optimality is not clearly defined. This SR-based gradient solution procedure searches for the descent direction through the service reliability vector ρ . It does not involve cutting the feasible region which will lead to increases in the problem or constraint size from iteration to iteration, as in L-shaped or Multi-cut methods. Firstly, for a set of given reliability measures ρ , we obtain the respective optimal route assignment X by solving P1 . Secondly, with X fixed in the phase-1 problem, scenarios are generated to solve P2 and estimate the expected total cost of the phase-2 problem. Then, we approximate the partial derivative of the expected total cost with respect to the reliability measures by sensitivity analysis. The partial derivative is then used to determine the descent direction and update the reliability measures. The process is performed repeatedly until the stopping criteria are satisfied. Optimizing the step size for sensitivity analysis
In the sensitivity analysis, to approximate the partial derivative of the expected total cost with respect to the reliability measures, it is necessary to find reliability measures that produce a different phase-1 optimal solution * X . Otherwise, the expected total cost in phase-2 will not change. One simple way to estimate this partial derivative by perturbation analysis is to increase the reliability measure gradually by a small fixed value, such as 0.05, iteratively until the optimal solution in phase-1 changes (e.g. An and Lo, 2014; Li et al., 2018; Lo et al., 2013). However, doing so will require P1 to be solved many times. To reduce the number of times that P1 needs to be solved, we find the demand volumes and detour times that will make some constraints in P1 active or binding, so that the optimal cost or * X in phase-1 may change. In this section, we establish the necessary conditions that the optimal route assignment * X and objective value in phase-1 may change. 7 An algorithm is proposed to find I e , the maximum increment of the demand volume without violating the constraints, which is derived in Appendix D, and justified in Proposition 4. The equation (40) below calculates the maximum number of additional service requests of type e can be taken by vehicle v that does not affect the phase-1 solution. RI R SII ISR S I min( )min , , ev iv veev i vZ cap t ttn t − − − = (40) i v is the number of passengers onboard the i th zone along route p for vehicle v . ev Z Z is the set of visited zones of vehicle v (ending zone excluded) for any service request of e . *R v y is calculated from * Y in the previous solution by (15). The complete procedure to find I e is formally described in Algorithm 1. Algorithm 1. Finding the maximum demand volume increment
Input: Current optimal route assignment * X , passenger assignment * Y , route set P , vehicle set V , and demand category e E . Step 0: Initialize the maximum demand volume increment I e . Step 1: Find the set of all routes e P P that is possible to serve e . Step 2: For every route e p P Step 2.1:
For every vehicle v V such that pv x = Step 2.1.1: Calculate the bound I ev by equation (40) Step 2.1.2: I I I e e ev +
Output: I e Proposition 4.
The objective value of P1 will not change if the volume of demand category e , I e , is increased by any I I e e , where I e is the maximum demand volume increment given by Algorithm 1. Proof.
See Appendix E. □ Similarly, Proposition 5 shows that the objective function will not change if the travel time between requests is increased by any non-negative value less than II z for zone z given passenger to vehicle assignment * * { } ev y = Y as defined below. The superscript II is the same indicator as detour reliability II , denoting detour time. is a small number to avoid division by zero. z zII * min vv vzz V yt t −= + (41) Proposition 5.
Given passenger assignments to vehicle * Y , for any demand category e , the objective value of P1 does not change if the travel time between requests of zone z is increased by any II II z z
Proof.
See Appendix E. □
8 Algorithm 1 and equation (41) determine the step size in the perturbation analysis of partial derivative more effectively, especially when the variance of the volume or detour time probability distributions are small such that slight changes of reliability measures do not change the optimal solution in phase-1. However, although they are necessary conditions for invoking change in the objective function, they are not sufficient conditions as they only detect when the constraints will become active with fixed * X and * Y . Moreover, even rendering these constraints active will change * X and/or * Y , the objective value may still not change. Nevertheless, it improves the efficiency of the algorithm, especially when it takes a long time to solve P1 . It is used in Step 5.1.1 and 5.1.4 in Algorithm 2. Detailed solution procedure
Recently, various methods for stochastic gradient descent optimization have been developed, such as AdaGrad (Duchi et al., 2011), RMSProp (Hinton et al., 2012), Adam (Kingma and Ba, 2014), etc., as reviewed in Ruder (2016). In this paper, we solve the problem by a mini-batch gradient descent optimization, combining the gradient solution developed in An and Lo (2014), and the Adam method. The Adam method will be used only if the objective value does not improve in three consecutive iterations, as indicated by A = if Adam method is applied, and A = if the analytical method is used. To further enhance the computation efficiency, a decomposition approach to reduce the problem size is introduced in Appendix F. Algorithm 2. Service reliability based gradient descent optimization
Step 1: Set k =1, initialize I II , )( k k k = ρ ρ ρ , where subscript k denotes the iteration indicator, and initialize a binary variable A as 0, indicating that the Adam method is not used. Step 2: Based on k ρ , solve P1 to find the optimal , X Y . Step 2.1: If P1 has no feasible solution, the value of k ρ is reduced until a feasible solution is found in P1 . Step 3: Given the vehicle routes X , generate service request scenarios, and assign the service requests to vehicles by P2 . Calculate the expected ad hoc cost by (38) and the total cost )( k C ρ by (39); set the optimal objective value * C as )( k C ρ if * C does not exist or * ( ) k C C ρ . Step 3.1: If the percentage change of objective value is less than a threshold, say 1%, terminate. Step 4: If the optimal objective value * C does not improve in three consecutive iterations and A is 0, set A as 1, and Adam method will be used in Step 5.2 to determine the step size. Step 5: Optimize ρ . Step 5.1: Calculate the partial derivative of )( C ρ with respect to I ρ and II ρ through the perturbation analysis. As X is binary, it is impossible to analytically find the partial derivatives of )( C ρ with respect to ρ . Therefore, we use perturbation analysis to find the partial derivative as shown in the procedure below. Step 5.1.0: Given k ρ and its corresponding solution ( , ( )) k k C X ρ Step 5.1.1: (a)
Run Algorithm 1 to obtain the maximum demand volume increment I e such that the constraints are not violated. Set I' e k by inverting e such that 9 I e e e + + , other reliabilities being the same as k ρ , denoted as I ' e k ρ . (b) Solve P1 again to obtain the route assignment I ' e k X , if the objective value of P1 changes, go to Step 5.1.2; otherwise, go back to Step 5.1.1(a). Step 5.1.2: Check if I ' e k X is in the solution set or not and get the objective value C directly if it does; otherwise solve P2 repeatedly to obtain C , calculate the sensitivity of I ek : I 'I I I ( ) ( ρ ( ) ) e ke k kk e k ek CC C − − ρ ρ ρ . Save I I '' )( , ( ) k ee k C ρ X into the solution set for future usage. Step 5.1.3: Repeat step 5.1.1 to 5.1.2 for every demand category e in I ρ to obtain sensitivities over I k ρ , denoted by the vector I )( k C ρ . Step 5.1.4: (a) Obtain the detour time increment II z such that the constraints are not violated by equation (41). Set II' z k such that II II II z z z + + by inverting z , other reliability measures being the same as k ρ , denote it as II' z k ρ . is a very small value to violate some constraint. (b) Solve P1 again to obtain the route assignment II' z k X , if the objective value changes, go to Step 5.1.5; otherwise, go back to Step 5.1.4(a). Step 5.1.5: Check if II' z k X is in the solution set or not, get the objective value C directly if it does; otherwise solve P2 repeatedly to obtain II' )( z k C ρ , calculate the sensitivity of II zk : II'II II II ( ) ( ρ )) ( z kzk kk kz zk C C C − − ρ ρ ρ . Save I II'I ' )( , ( ) z kz k C X ρ into the solution set for future usage. Step 5.1.6: Repeat Step 5.1.4 to 5.1.5 for every zone z in II ρ to obtain sensitivities of ( ) k C ρ over II k ρ , denoted by the vector II ( ) k C ρ . Step 5.1.7: Join the vector I )( k C ρ and II ( ) k C ρ to become a vector of all reliabilities ( ) k C ρ Step 5.2: Determine the step size and update the reliability. If A is 0, go to step 5.2.1. Otherwise, go to Step 5.2.2. Step 5.2.1: Update the reliability by the same way as in An and Lo (2014). Step 5.2.1.1: Take the negative sensitivity vector )( k C − ρ as the descent direction. The step size k is chosen as: *2 (( ) ) kk k k CC C − ρ ρ . Step 5.2.1.2: Update the reliability ρ in next iteration. Go to Step 5.2.3. Step 5.2.2: Apply Adam method (Kingma and Ba, 2014) to optimize reliability. Step 5.2.2.1: If it is initialized, go to Step 5.2.2.2. Otherwise, set two parameters k − m 0 and k − v 0 , is set to be a very small number, say − . Step 5.2.2.2: Calculate (1 ) )( k k k C − + − m ρ m Step 5.2.2.3: Calculate
22 1 (1 ) ) )( ( kk kk
C C − − + v v ρ ρ Step 5.2.2.4: Calculate ) ˆ / (1 kk k − m m Step 5.2.2.5: Calculate ) ˆ / (1 kk k − v v Step 5.2.2.6: Update the reliability ρ in next iteration by 0 ˆ ˆ )( k kk k k + = − + ρ ρ m v , an elementwise operation of each reliability. Step 5.2.3: If any element in k + ρ or any element in k + ρ , project that element onto [0,1), set k k = + . Go to Step 2. k , , k , and are parameters determined by trial-and-error. k and k are the k th power of and , respectively. denotes element-wise multiplication (also called the Hadamard product), and denotes the element-wise division of two vectors (also called the Hadamard division). The solution approach outlined above will stop if the stopping criteria is matched, resulting in an optimal X and total cost * C . Numerical examples P1 , P2 and the algorithms are implemented by Python 3.7 PICOS API, and then solved by IBM CPLEX 12.10. The formulation will first be demonstrated by a small example in Section 3.1. Then, a five-zone scenario is introduced in Section 3.2, and solved by Algorithm 2 introduced in Section 2.4.5. Subsequently, sensitivity analyses on the vehicle capacity and detour limit are conducted based on the five-zone scenario in Section 3.3. The computation times under different problem sizes are analyzed in Section 3.4. Finally, a case study based on real data in Chengdu, China, is illustrated in Section 3.5. Small illustrative example on stochastic volume and detour time
We illustrate the ZBFBS formulation in a small area with three zones, with their stochastic demand characteristics shown in Table 1. There is one and only one demand category in OD pair AC. Truncated normal distribution is used to describe the volume and detour time distributions, denoted as ( , , ), TN a b that a normal distribution of mean and standard deviation are bounded by [ , ] a b . As the truncated normal distribution is a continuous distribution, the number of service requests generated is rounded off to the nearest integer. Note that the proposed formulation and algorithm can handle any discrete distribution of integer for the demand distribution, and any continuous distribution for the detour time distribution. Moreover, there are 25 vehicles available, each with vehicle capacity cap = . The ad hoc cost d c to serve a request is 90% of the operating cost of one vehicle (i.e., if service request d is from OD pair IJ, d c c = ). c is the operating cost of the direct route from I to J for a flexible bus and d c is the ad hoc cost of service request d . The ad hoc cost per service request is much higher than the regular service operating cost per service request, as vehicles are shared in regular service, but not in ad hoc service. The maximum allowable detour time is 8 minutes per zone. The detour times from the boundaries of zone A and C are given by A ( ) 0.7 0.03 vz vz y y = − and C ( ) 0.6 0.03 vz vz y y = − respectively, vz y is the number of service requests served by vehicle v in zone z . In phase-2, the detour time reduction for the detour matrix z T is calculated by, ( ) min / 24, z z zdb d b − = − (42) zd and zb are the detour times required in zone z to serve service request d and b respectively, generated from the detour time distribution z . The denominator 24 is to fulfill condition (13) as the maximum vehicle capacity in the numerical examples is 12. 1 Table 1. Demand categories and characteristics Demand category id ( e ) OD pair Demand volume distribution ( e ) Zone A detour time distribution ( A ) Zone C detour time distribution ( C ) Number of passengers per request ( e n ) 1 AC )(16, 6, 0, TN + )(1.5,1, 0, TN + )(1,1, 0, TN +
1 The shortest route is ABC, with an operating cost of HK$10. As the problem size is small, the P1 solutions can be worked out by enumeration. Specifically, reliability measures from 0 to 1, discretized by 0.05 intervals, are enumerated for the stage 1 problem. The expected ad hoc cost in P2 is taken as the average of optimized ad hoc cost of 150 randomly generated scenarios based on the above demand profile. It is observed that there are multiple optimal reliability measures leading to the same best route design, such as I II II1 A C ( , , ) (0.3, 0.25, 0.3) = = ρ or (0.4, 0.2, 0.25) = ρ , in which the expected total cost is HK$24.1 and two vehicles are assigned to route ABC. The former solution gives a higher reliability measure to origin detour, but lower reliability measures to demand volume and destination detour, while the latter solution produces a higher reliability measure in demand volume and destination detour. The second-best deployment has higher reliability measures, e.g., (0.4, 0.45, 0.4) = ρ which assigns three vehicles to route ABC with the expected total cost equal to HK$30.0. This enumeration also shows the trade-off between the regular service operating cost and the ad hoc service cost. Reliability measures that are too high provide excessive regular services that will increase the fixed operating cost f C , but reliability measures that are too low will increase the ad hoc cost Q to serve the excessive demand. The best solution comes from a reliability measure that balances the trade-off between operating cost and variable ad hoc cost, as illustrated in Figure 3. Figure 3. The costs with respect to number of vehicles in this example The number of solutions is limited due to the integral constraints of vehicle to route assignment X . In other words, for some different reliability measures, the same X is 2 generated. Therefore, there are multiple reliability measures that yield the same global optimal solution. Moreover, the total cost C is non-convex with respect to ρ . To show the non-convexity, we enumerate the volume reliability measure I from 0 to 1, discretized by 0.001 intervals, and solve P1-P2 for each volume reliability measure enumerated, given detour time reliability measures
IIA and II C as 0.7. The result, depicted in Figure 4, illustrates the non-convexity of the total cost with respect to volume reliability. It is shown that the total cost increases rapidly at some points where one more vehicle is required to serve the service requests. Therefore, convex optimization approaches that are well-studied in the literature are not applicable to this problem. Figure 4. The score with respect to I , where II IIA C = =
Numerical examples in a five-zone scenario
To demonstrate the capability of the SR-based gradient solution approach, we divide the service area into five zones as shown in Figure 5. The blue dots denote their centroids, and the blue lines are the links connecting the zones, where the blue numbers are the link costs per 10-seater vehicles in HK$. The examples are performed on a desktop computer with Intel i7-3770 3.4Ghz CPU and 32Gb ram. The decomposition outlined in Appendix F is applied to this and subsequent examples to reduce the solving time without loss of generality. Figure 5. The five-zone network topology 3 The demand categories, with their demand volume distributions and numbers of passengers, are given in Table 2. The probability distributions of the travel times between two service locations are given in Table 3. These probability distributions are also used to generate diagonal values of z T , the matrix of detour time for each zone z . The correlations of the detour time, or detour time reductions, are calculated by (42) as in the first example, at the beginning of Section 3.1. Table 2. Volume stochasticity for the five-zone network OD Demand volume distribution ( e ) Number of passengers ( e n ) OD Demand volume distribution ( e ) Number of passengers ( e n ) AB )(3, 4, 0, TN +
1 DE )(5, 4, 0, TN + )(3, 4, 0, TN +
2 ED )(4,9, 0, TN + )(3, 4, 0, TN +
1 CE )(2, 4, 0, TN + )(3, 4, 0, TN +
3 EC )(2,1, 0, TN + )(4, 4, 0, TN +
2 EB )(2,9, 0, TN + )(4, 4, 0, TN +
1 DB )(6,1, 0, TN + )(4,9, 0, TN +
1 BD )(4, 4, 0, TN + )(2,1, 0, TN +
3 CB )(5,9, 0, TN + )(6,16, 0, TN +
2 EA )(2,1, 0, TN + )(3, 4, 0, TN +
1 Table 3. The detour time distributions used for zones Zone Detour time for the segment between two points ( z ) Detour time from the boundary ( ( ) z vz y ) A )(1.5,1, 0, TN + /12 vz y e − B )(1,1, 0, TN + /12 vz y e − C )(1,1, 0, TN + /12 vz y e − D )(1,1, 0, TN + /12 vz y e − E )(1.5,1, 0, TN + /12 vz y e − Under assumption B1, only the shortest route serving each OD is considered, so the candidate route set P includes routes AB, BA, AEC, BC, AD, DA, AE, BE, CD, DE, ED, CE, EC, EB, DEB, BED, CB, and EA. The operating cost is given by summing link costs that a route traverses. A fleet of 80 10-seater vehicles and maximum detour of 8 minutes per zone per vehicle are considered. The problem is solved by the SR-based gradient solution approach outlined in Section 2.4.5. The stopping criterion is set to be 1% of the relative difference between the objective values. k = , = , = , k = and = . The ad hoc cost is 90% of the operating cost of one vehicle, the same as the small illustrative example. Convergence of SR-based gradient solution approach
As this stochastic formulation is inherently non-convex, different initial values of ρ to start the solution procedure may yield different solutions. To compare the effects of different initial points, we initialize the solution approach with different service reliability measures. Every entry in the initialized ρ is enumerated by roughly an interval of 0.1. The optimal 4 solutions resulted from different initial points are plotted in Figure 6. The objective values are represented by the blue lines as the iteration proceeds. Vertical grid lines separate different initial points of service reliabilities, as indicated at the top of each box. The width between each vertical grid line represents the number of iterations, also indicated by the number at the end of each blue line. For instance, for the initial point of = ρ , it takes 5 iterations to converge and the local minimum is HK$952.1. The best local minimum is HK$919.5 with the initial point = ρ . Figure 6. Objective function values versus iteration numbers for different initial points Comparison between deterministic and reliability-based approaches
In this section, we compare with adopting the deterministic approach to find the optimal solution by simply making use of the mean and median of the distributions as point estimates to plan for vehicle zonal routings and assignments, while ignoring the stochastic nature of demand. We then evaluate the deterministic solution by simulating its performance on 200 randomly generated demand scenarios and record its computation time and ad hoc service cost. The objective functions are evaluated and compared with the optimal solutions of the stochastic formulation with different starting points as in Section 3.2.1. Given the current volume and detour time distributions, the computing time is reduced by using deterministic approach and the resultant total cost obtained by the reliability-based approach is lower by about 5%. Nevertheless, if the distribution is not truncated normal distribution, which is sort of symmetric, but follows a lognormal distribution as follows, which is more common, the results are very different:
1( ) 0log ( / 4)20.8 exp 1. 82 , f xx xx = − (43) The computation times and total costs according to the deterministic and reliability-based approaches under lognormal detour time distributions are shown in Table 4. Table 4. Computation time and cost using deterministic and reliability-based approach Mean result Median Optimal reliability-5 result based result Computation time (s) 581 1138 3682 Number of iterations - - 4 Vehicle used 40 40 18 Operating cost (HK$) 1820 1790 820 Ad hoc service cost (HK$) 970 952 1359 Total cost (HK$) 2790 2742 2179 The results demonstrate the capability of the reliability-based solution approach in finding good solutions for all sorts of distributions, while the deterministic approach may not. In the case of the demand following the lognormal distribution, as shown in Table 4, the reliability-based approach reduces the total cost by 20.5%, from HK$2742 to HK$2179, as compared with the deterministic approach.
Sensitivity analysis on vehicle capacity and maximum allowable detour time
In the planning stage, there is a trade-off between economies of scale of capacity and expected occupancy. A larger vehicle yields a lower operating cost per space due to the same fixed costs such as driver cost, vehicle registration fee, but seats are wasted if insufficient passengers are picked up. This section conducts a sensitivity analysis on the vehicle capacity and maximum allowable detour time. Four scenarios with vehicle capacity of 7, 8, 10, 12 are considered, with operating cost factored by 85%, 90%, 100% and 110%, respectively. If a ride request is not served by regular service, an ad hoc service, with a cost of 80% of the link cost shown in Figure 5, is assigned to serve that request. The maximum allowable detour time is 8 minutes. Analyses in this section are conducted by repetitively solving P2
200 times for the best routing X obtained in the SR-based gradient solution approach. The average detour times and vehicle occupancy across all scenarios are taken, as shown in Table 5 and illustrated in Figure 7. Table 5. Results under different vehicle capacities and operating costs Scenario 1 Scenario 2 Scenario 3 Scenario 4 Vehicle capacity cap
7 8 10 12
Factor of operating cost 0.85 0.9 1 1.1 Total cost (HK$) 980.7 965.7 939.6 980.8 Vehicle operating cost f C (HK$) 799 657 810 671 Ad hoc cost (HK$) 181.7 308.7 129.6 309.8 Average volume reliability I ρ II ρ customers’ perspective, due to their lowest detour time per zonal visit (3.97 minutes) as compared with 8-seaters (4.51 minutes), 10-seaters (4.42 minutes), and 12-seaters (4.95 minutes). The vehicle occupancy is generally higher with a smaller-size vehicle. The total cost is higher for vehicle capacity of 7 and 12 and is the lowest for the vehicle capacity of 10, as it provides flexibility while saving the operating cost. In terms of occupancy, vehicles of capacity 10 and 12 yield lower occupancy than vehicles of capacity 7 and 8. The maximum detour times for all vehicle types are far lower than the detour time constraint (set to be 8 minutes) in all three scenarios, indicating that detour time is usually not a binding constraint. The minimum slack in detour time, i.e. the maximum allowable detour time minus the average maximum detour time, is 2.23 minutes across all scenarios at Scenario 4. The maximum detour minute decreases from 5.77 to 4.92 with the usage of smaller vehicle since fewer passengers are picked up per zone on average. Note that the occupancy and total detour time are non-monotonic due to the amount of vehicles used. For example, as there are 18 flexible buses assigned in Scenario 3, i.e., more seats are provided, so that the occupancy is lower than the case when 12-seaters are used. The discussion above poses another question: How will the service patterns change by varying the maximum allowable detour time constraint, from 6, 7, 9, to 10 minutes? The results are summarized in Table 6. Total cost, vehicle occupancy, total detour time and detour time in zone versus detour time limit are illustrated in Figure 8. 7 Table 6. Detour and occupancy under different detour time limits Scenario 5 Scenario 6 Scenario 3 Scenario 7 Scenario 8 Detour time limit (mins)
6 7 8 9 10
Vehicle capacity
10 10 10 10 10
Total cost (HK$) 1119.3 988.2 939.6 887.0 882.9 Vehicle operating cost f C (HK$) 760 730 810 670 670 Ad hoc cost (HK$) 326.3 258.2 129.6 217.0 212.9 Average volume reliability I ρ II ρ Results and computation time for larger instance
As the number of passengers increases, more variables are introduced to
P2, which increases the problem size and potentially increases the solution time. Following the setting of Scenario 3, this example evaluates the performance of the flexible bus system and computation time of the service reliability-based gradient solution approach based on different demand volume level. Two scenarios that multiply the parameter of demand volume distribution by 1.2 and 1.5 respectively are simulated. The detailed result is shown in Table 7. Table 7. Output parameters under different factor of service request Scenario 3 Scenario 9 Scenario 10 Demand volume multiplier
1 1.2 1.5
Total cost (HK$) 939.6 1164.8 1388.7 Vehicle operating cost f C (HK$) 810 730 1130 Ad hoc cost (HK$) 129.6 434.8 258.7 Average volume reliability I ρ II ρ Application for Chengdu, China, with real data
Finally, to illustrate the formulation for a real scenario, we apply it to Chengdu, China, by using the empirical service request data from Didi Chuxing from 1 st November to 30 th November 2016. The flexible bus service is planned based on the data from 9:00am to 9:15am as representative of the morning peak. The area within the second ring road is divided into nine zones as shown in Figure 9. The amount and detour times of service requests are derived from the real data, but the volume is scaled down to 20%, as detailed in Appendix G. The cost to travel from a zone to an adjacent zone with shared border (e.g. from zone A to zone D) is RMB¥5, while to an adjacent zone with shared corner only (e.g. from zone A to 9 zone E) is RMB¥7. To cope with the larger zone size, the maximum allowable detour time within a zone is 15 minutes, the vehicle capacity is 10 passengers. Figure 9. The study area in Chengdu, zone divisions and labels (reproduced from OpenStreetMap) The solution shows that 45 vehicles are used. Specifically, each of AE, BC, BD, BEH, CB, CE, CFI, DBC, EA, ED, EG, EH, FBA, FB, FH, GH, HEB, HEC, HD, HF, HG, IE, IF, IH, CFH and BFI is assigned a vehicle, two vehicles are assigned to route DEF, DHI, EB, FC and FED, and three vehicles are assigned to EC, EF and EI. The total cost is RMB¥496.11, where the regular service cost and ad hoc service cost are RMB¥343.83 and RMB¥152.28, respectively. 88.2% of the passengers are served by the flexible bus, the average number of service locations visited per zone is 4.20, and the average detour time per zone is 10.93 minutes. As the average demand volume is large, the computation time is 165 hours in 6 iterations. Conclusion
In this paper, we proposed a zonal-based flexible bus operation strategy. Demands were clustered geographically into zonal OD pairs and number of passengers. The volume and spatial stochastic variations of demand were estimated by historical record. A two-stage stochastic problem with recourse was formulated to plan the flexible bus routes. Volume 0 reliability measure and detour time reliability measure were introduced to divide the otherwise intractable two-stage stochastic problem into two phases. In phase-1, the zonal visit sequences of vehicles were optimized to serve customer demand up to a certain level of reliability, while satisfying the vehicle capacity and maximum allowable detour time. Upon the realization of demand, the phase-2 optimization matched service requests to vehicles. A service-reliability based gradient solution approach was then implemented to solve the problem by finding optimal reliability measures that minimized the sum of fixed operating cost and ad hoc service cost incurred for customers who cannot be served by regular service. Given the gradient direction, the step sizes were calculated by either a procedure or the Adam method. In our numerical experiments, this hybrid method converged more consistently than either method. It was proved that any optimal solutions obtained in the original two-stage stochastic problem with recourse P0 can be obtained by the SR-based gradient solution approach P1 - P2 , and vice versa. The numerical studies illustrated that the SR-based gradient solution approach consistently produce good solutions under different probability distributions. In a case study, it was demonstrated that vehicles with capacity 10 were preferred, based on the trade-off among the total operating cost, detour time, and occupancy. The results provided insights for operators to decide the routing, vehicle capacity and fleet size based on the aforementioned criteria. The result for the numerical example with a larger demand showed that the computation time increased rapidly with the problem size. Applying the model to Chengdu, China, based on real sample data required over 165 hours of computation time. How to develop a more efficient solution approach for a large demand is an important topic for further studies. The proposed model provides a novel way to formulate and optimize a flexible bus system. The zonal visit sequences can either be treated as traditional long-term strategic planning that is optimized when the demand pattern changes substantially, or short-term planning that is responsive to demand changes due to stochastic demand characteristics. Operationally, the actual route is determined by solving the passenger-to-vehicle assignment problem in real-time, drivers then have the flexiblility on routing as long as all the required pickup or drop off points are visited. Currently, the objective function is only based on the benefit of the mobility service providers, and the benefits of other stakeholders are not considered. Future researches can be conducted to include different stakeholders’ benefits, such as user costs, waiting time and profit. The other possible extensions include improving existing or developing alternative solution methods to solve larger problem sizes more efficiently, incorporating time windows for pick-up and drop-off, and optimizing for dynamic demand realization. Appendix A. Illustration of converting matrix p B , detour correlation matrix z T and capacity constraint To illustrate detour constraint (5) and vehicle capacity constraint (7), we consider demand realization with 4 different service requests, with the detailed information shown in Table 8. The OD of passengers onboard in each portion of route ABC is shown in Figure 10. 1 Figure 10. The OD of passengers onboard from A to B and from B to C. Table 8. Detailed information of 4 service requests Service request no. d Number of passengers d n Origin (detour time d + , minutes) Destination (detour time d − , minutes) Ad hoc cost d c (HK$) 1 2 A (1) C (2) 6 2 3 A (2) C (1) 6 3 2 A (3) B (2) 3 4 1 B (2) C (2) 2 From the above table, we have A A A B C C B C1 2 3 4 1 2 3 4 + + + + − − − − = = = = = = = = indicates the origins and destinations of each service request, and the other ( ) zd + − values are zero for every zone z and service request d . Moreover, the detour time matrices are given by, A B C − − − − − − − − = = = − − − − − −
T T T
The diagonal of matrix z T is the detour time required to serve respective service requests. The other elements are the reduction of detour time by picking up or dropping off more than one service requests at a zone simultaneously. This matrix satisfies condition (13) as the absolute value of the sum of the detour time reduction in every row is smaller than half of the additional detour time. For example, in the first row of A T , the absolute value of sum of the detour time reduction is 0.4, which is smaller than half of the additional detour time 1 at that row. Thus, adding any service requests to a vehicle will only increase the detour time. With the correlation matrices, the detour constraint (5) is stated as, v v v v v v vvv vv vv v v v v v vvv v t v Vt v Vt v Vw w w w w w w w ww w w ww w w w w w w w w + + + − − −+ + −+ − − − (44) where dv w is the binary variable indicating whether request d is served by vehicle . v Equation (6) calculates the number of passengers for each OD pair served by every vehicle v V as follows. 2 AC 1 2AB 3BC 4 v v vv vv v === + (45) We use route
ABC p = , i.e. from zone A to B, then from B to C, to illustrate the converting matrix p B , matrix ABC B is given by,
1( ) 0( ) 1( )0( ) 1( ) 1( ) b b b M M Mb b b M M M = = == = = = B (46) The rows represent the zones where passengers can be picked up; in this case, zone A and B, correspond to the first and the second row respectively. The columns show all the OD pairs in the following order: AB, BC, AC, CB, BA, and CA. Precisely, the passengers of OD pairs AB and AC are onboard in zone A, with corresponding ABA b and ACA b equal to 1 in the first row. ABB b is zero because passengers of OD pair AB will not onboard when the vehicle leaves zone B. Moreover, the entries for OD pairs CB, BA, and CA are represented by a large positive number M , which ensures that every vehicle of route ABC does not serve any service requests from these OD pairs that route ABC does not traverse. With the converting matrix ABC B , if ABC, v x = , the capacity constraint (7) for vehicle v can be expressed as: AB AC CB BA CA1 1 1BC AC CB BA CA1 1 1 for zone Afor zone B v v v v vv v v v v
M M MM M capcapM + + + + + + + + (47) Again the large positive number M ensures that vehicle v of route ABC will not serve any passengers other than the OD pairs in the route ABC, i.e., CB v , BA v and CA v must be zero if ABC, v x is equal to 1. If there are two vehicles in the fleet, i.e., V = , v = for the first vehicle and v = for the another, and cost of route ABC HK$10 c = , combining constraints (44), (45) and (47), the optimal solution is to deploy one vehicle to route ABC if the maximum allowable detour time of each zone BA C , 4.2, t t t and the vehicle capacity cap , and the total operating cost is HK$10. If A t and cap , the optimal solution is to deploy one vehicle to route ABC to serve request 1, 2, 4 and allocate ad hoc service for request 3, i.e., ,1 ,2
1, 0
ABC ABC x x = = , , , , 1,1, 0,1 w w w w = and the total operating cost, including the HK$3 additional ad hoc service cost is HK$13. There is a difference since the sum of detour time in zone A is 4.2 if we serve service requests 1, 2 and 3 simultaneously. Note that as there are multiple optimal solutions by swapping routes of vehicles, for instance, the objective function above is indifferent with ,1 ,2
0, 1
ABC ABC x x = = , , , , 1,1, 0,1 w w w w = . 3 Appendix B. Proof of Proposition 1
Proposition 1.
Any feasible solution of
P1-P2 is also feasible to the original two-stage stochastic formulation P0 . Proof.
We consider a solution ( , , )
X Y W after solving P1 and P2 . Suppose , X Y is a solution of P1 for a given service reliability ρ , and W is a solution of P2 for demand realization . Constraints (2) and (3) in P0 are satisfied as they are the same as (17) and (25) in P1 . Moreover, constraints (5) – (10) are satisfied as W is feasible in P2 and (10) is equal to (37). In short, for any feasible solution ( , , ) X Y W of P1 and P2 , constraints (2), (3) and (5)-(10) are fulfilled, so it is feasible in the two-stage stochastic formulation P0 . □ Appendix C. Derivation and Proof of Proposition 2 and Proposition 3
To show that optimal solutions of P0 can always be obtained by solving P1 for some reliability measures, with the condition that there is one demand category per OD, we show that all solutions of P0 are reproducible from P1 by setting appropriate reliability measures. To facilitate the proof below, some notation and definitions are introduced. Firstly, by assumption B1, there is always one shortest path (defined as direct route ) starting exactly from zone I and ending exactly in zone J for every OD pair (I, J) , denoted by IJ P , which has an operating cost denoted by c . Secondly, we introduce the notion of active detour time constraint. Given any demand category of OD pair IJ and any vehicle v V with route IJ P , the capacity constraint (21) ensures that the number of service requests of type e picked up by vehicle v , i.e., ev y , is bounded above by e capn , e n is the number of passengers in e . Moreover, from the detour time equations (14), (15) and (23), ev y is also bound by ( ) ( ) ( ) ( ) II I J Jmax,time II II J ev y t t = + − + − . The detour constraint is said to be active if max,time eve c yapn , which means that the maximum number of passengers picked up is constrained by the maximum detour bound, but not by the capacity constraint. In the Lemma 1 and Proposition 2 below, by setting appropriate detour reliabilities, the detour constraint is not active, thus the number of requests of service request e , with OD pair IJ, picked up by any vehicle serving route IJ P is bounded by, ev e ncapy (48) Lemma 1.
Given that only the demand category of one passenger exists for each OD pair, if the detour constraints are inactive, and vehicle capacity cap is a factor of demand volume , e e E , serving all demands by their respective direct routes minimizes the operating cost. (i.e. better than or equal to other combinations). Proof.
Define
ABIJ as the cost of providing a space for one passenger from zone I to zone J 4 by route AB P . The spaces may or may not be occupied by a passenger, but the number of spaces provided must be higher than the number of passengers of that OD served. Based on the objective function (16), assigning a vehicle to route IJ P costs c . On the other hand, equation (48) states that it can provide space for at most cap service requests of OD pair IJ as all demand categories have one passenger ( e n = ). Therefore, for any pair of zone IJ, we can express the cost of space provided by its direct route as, JIJIJ 1I ccap = (49) Apart from the direct route IJ P , service requests of OD pair IJ may be served by vehicles with three types of routes traversing zone I and zone J: (1) route starts in zone I, passes through zone J and ends in some other zones; (2) route starts somewhere neither I nor J, passes through zone I and ends in zone J; and (3) route starts somewhere neither I nor J , passes through zone I and zone J respectively, and ends elsewhere. We will show that serving the service requests by these routes does not reduce the operating cost. In case (3), there is a route AB P starting at specific zone A that is neither I nor J, passes through zone I and then zone J, and ends at specific zone B that is neither I nor J. Following the vehicle capacity constraint (21), if a vehicle is serving n service requests from zone I to zone J, it can serve at most cap n − service requests from A to B and n service requests from A to I and from J to B, as illustrated in Figure 11. Serving one customer from I to J will occupy one space for customer from A to B, vice versa. The space from A to I, and J to B obtained by serving customers from I to J can either be filled in by customers from A to I or left vacant. Therefore, for any n cap , we have the following relationship between ABAI , ABIJ , ABJB and ABAB . BAB AB AB ABAI IJ JB AB A )( n n n cap n c + + + − = (50) Figure 11. Seat allocation of route AIJB. From assumption B2, we have AB AI IJ JB c c c c = + + , by substitution to (50),
AB AI IJ JBAB AB AB ABAI IJ JB AB ( ) cap nn n ncap c c c ccap can p cap cnn apn −+ + + − = + + + (51) Since AB P is the (shortest path) direct route to serve demand from A to B, from (49), we have BABAB A ccap = , by substitution to (51) and division by n , 5 AB AB ABAI IJ JB AI IJ JB c c ccap + + = + + (52) Equation (52) establishes the relationship between costs of providing space
ABAI , ABIJ , ABJB by route AB P and route operating costs AI c , IJ c and JB c . It is equal to the cost of providing space with the direct route, as the cost of providing a space for passengers of OD pair AI, IJ and JB by direct route are IAIAI A ccap = , JIJIJ I ccap = and
BJBJB J ccap = respectively. As vehicle capacity is a factor of the demand volume, serving every service request by the respective direct routes left no vacant space. The other allocation of service requests
AI IJAB , ,
P P P and JB P that also waste no seats have the same cost as shown in the above as all spaces are used up. The same conclusion can be drawn for case (1) and case (2) mentioned above by either removing terms with AI or terms with JB. To conclude, as the vehicle capacity is a factor of demand volume, serving the service requests by direct route leaves no empty seats, thus minimizes the cost of providing space. Therefore, the operating cost is minimized by serving the service requests by the direct route. □ Proposition 2 . For problems with the stochastic volume described by a probability distribution for any integer greater than 0, and the stochastic detour time described by a probability distribution over the range [0, ) + , if there exists only one demand category with one passenger per OD, any optimal solution of the original 2-stage stochastic formulation can be generated by the SR-based stochastic formulation proposed. Proof.
The general idea of the proof is that for any optimal solution * * ( , )
X W in P0 , there exists some I ρ and II ρ such that the resultant solution from P1 is * X . Suppose P0 has an optimal solution * * ( , ) X W , we can calculate * p m , the number of vehicles assigned to route p (i.e. * * p vv V p m x = ). For every route p where * p m , we find the corresponding demand category e E with identical origin zone I and destination zone J with route p . Then, for that demand category e , with the condition that the demand volume equal to the amount of vehicle spaces, * pe m cap = , and detour constraints are inactive (e.g. A reliability that makes every z = ), we can always obtain I II , e z by inverting the probability distributions given in Equation (18) and Equation (22). Combining reliabilities from all demand categories, we have I ρ and II ρ . By defining reliability as I ρ and II ρ , we have cap as a factor of e , and detour constraints are not active since they do not further restrict the number of service requests can be served. By Lemma 1, in the above condition, serving all service request of category e by vehicles with direct route p yields the least cost. Therefore, * X is an optimal solution to P1 , given * pe m cap = , II z and by designating I ρ and II ρ . In short, for any optimal solution * X of the original two-stage formulation, it is always possible to identify some I ρ and II ρ in P1 to reproduce * X . □ Proposition 3.
For the same problem in Proposition 2, with multiple demand categories per OD, including single passenger service requests, any optimal solution of the original 2-stage stochastic formulation can be generated by the SR-based stochastic formulation.
Proof.
This proposition extends Proposition 2 to multiple types of demand per OD. Given an optimal solution ( *, *)
X W in P0 , despite there are multiple types of demand per OD, we can always generate * X by varying the volume reliabilities of those single person requests and set the reliability measures of multiple people requests as 0. First, set the reliability I e = for demand category e such that e n , so the service provision for these demand categories is 0 by (18). Therefore, the corresponding ev y is always zero, and the problem is reduced to one type of demand with one passenger per OD by condition (20) and (23). Follow Proposition 2, for each route p with vehicle assigned to, that OD is (I, J) , we find a demand category e with OD pair (I,J) and one passenger (i.e. e n = ), and set each I e , II z such that * pe m cap = and II z = by inverting the probability distributions. □ Appendix D. Algorithms for step size determination
For every demand category e , there is a set of route e P P that can serve that demand category. The spaces left for vehicles assigned to those routes are checked, and the number of additional service requests from demand category e can be served is calculated. First, we consider the capacity constraint (21), which uses p B to map the number of passengers per OD v ζ to a vector of passengers onboard when leaving each zone along the route. We express p v B ζ by the following vector, ( ) T1 2 1, ... p p v v v m v − = ζ B (53) where iv is the number of passengers on vehicle v in the i -th zone along route p . The vector has only p m − entries as we do not need to consider the capacity constraint for the end zone. The capacity constraint (21) limits , , iv cap i v V . If * X does not change, the maximum number of additional service requests from demand category e can be served by vehicle v , denoted as I ev , is bounded by the minimum slack capacity along the route, written as, I min( ) ev ivZev ei capn − (54) e n is the number of passengers of demand category e , ev Z Z is the set of visited zones of vehicle v (ending zone excluded) for any service request of e . Next, we consider the detour time constraints. I ev can be calculated by rearranging the terms in the detour constraint (14), (15) and (23), where zone R and zone S are the origin and destination of demand category e . Therefore, we have equation (40) used in Section 2.4.5.1. 7 RI R SII ISR S I min( )min , , ev iv veev i vZ cap t ttn t − − − = (40) Appendix E. Proof of Proposition 4 and Proposition 5
Proposition 4.
The objective value of P1 will not change if the volume of demand category e , I e , is increased by any I I e e , where I e is the maximum demand volume increment given by Algorithm 1. Proof.
We prove by showing that for each demand category e considered separately, all the constraint will not be violated by increasing I e by any value I e less than I e . Based on I ev specified for each vehicle v in Step 2.1.1 in Algorithm 1, we define the increment of requests picked up by vehicle v as I ev , that satisfies, I I I I ev e e evv V v V = =
Furthermore, from Step 2.1.1 of Algorithm 1, we have,
I I , , ivev ev eve cap i Z v Vn −
Rearranging the term, we get I , , iv ev e ev cap i Z v Vn + As the passengers of demand category e are only onboard when leaving zones in ev Z , we retain the inequality p p v m cap − ζ 1 B for pv x = . Therefore, the constraint (21) is fulfilled. On the other hand, we have RI I IIR R ev ev v t t − at origin zone R of e . Rearranging the term, we have, ( ) ( ) ( )( ) ( )( ) ( )( ) I IIR R R RI* IIR R RI I IIR* * *R R R R* *R R* *RR RI R v v ev vv ev vv ev v evev vv v t y y yy yy ytt y = ++ − + ++ + + += (55) Therefore, detour constraint (23) is fulfilled at origin zone R. Similarly,
SI I II SS e v ve t t − fulfills detour constraint (23) at destination zone S. Therefore, every constraint in P1 is fulfilled and the objective value will not be changed. □ Proposition 5.
Given passenger assignments to vehicle * Y , for any demand category e , the objective value of P1 does not change if the travel time between requests of zone z is increased by any II II z z Proof.
For every vehicle v and every zone z, we have II *II ( ) / , z vz vz zz t Vyt v −
Rearranging the term, we have
I* I t , zvz vz z t v Vy + . As the new detour time between two requests in zone z will become II II z z + , the detour constraint (23) is fulfilled. □ Appendix F: Decomposition of the problem
Without loss of generality, we can decompose the problem to reduce the computational time if there exist a set of service requests i D D or demand category i E E and a set of routes i P P , where each service request i d D or demand category i e E can only be served by routes i p P , and each route i p P can only serve service requests in i D or demand categories in i E . Then, as the set of routes and set of service request neither affecting nor being affected by other routes and demand categories, P1 and P2 can be solved for a smaller set. For example, in the six-zone scenario depicted in Figure 12, if there is demand from OD pairs AB, AE, CD and BF, and the possible routes are AB, ACDE and BCDF, P1 can be separated into two parts. The first part is OD pair AB with possible route AB, the second part is of OD pairs AE, CD and BF, and possible routes ACDE and BCDF. Formally, { } P AB = , A B0 d d
D d D + − = , { , } P ACDE BCDF = and A E C D B F1 | + + = }{ 1 d d d d d d
D Dd + − + − + − = . This decomposition by OD can reduce computational time significantly as both P1 and P2 are NP-hard integer programming problem which the computational times increase exponentially with the size of demand and size of routes. Figure 12. The six-zone scenario Appendix G: Usage of real data: a case study in Chengdu
In this study, the data from Tuesdays, Wednesdays, and Thursdays between 1 st and 30 th of November from 9:00am to 9:15am are used. As the coordinate system of the original data is at GCJ-02 coordinate system, it is first transformed to a Cartesian coordinate system in meter unit. The Euclidean distance is then calculated based on the transformed coordinate system between zone centroid and the location. After that, the detour time is obtained by multiplying the distance by 0.003min/m, which came from a road detour ratio of 1.25 given in the study of Yang et al. (2018) and an assumption that the average speed is 25km/h in Chengdu city area. The demand volume distribution adapted the empirical distribution of number of ride requests between zones in the study period. The function z between the boundary of a zone and the closest service request is inferred by the following algorithm. As the zone is divided by squares, the boundary can be represented by 2 horizontal lines and 2 vertical lines which form a square between them. Algorithm 3. Approximate the detour time function between boundary and the closest service request from the data Input: Sets of origin and destination location at each zone, the boundary of each zone max min max min ( , , , ) z z z z x x y y
Output:
Approximation values of ( ) z i for each zone z Z and meaningful i . Algorithm: For each zone z Z For each number of service request considered i Conduct the below calculation 1000 times Randomly sample i locations in zone z as P , each with x, y coordinate denoted by ( , ) x y Denote max min max min , , ,
X X Y Y
P P P P as the maximal x coordinates, minimal x coordinate, maximal y coordinates, minimal y coordinates of locations in set P Distance max min max minmax min max min z z z zX X Y Y P x x P P y y Pd − + − + − + −
Set the function value ( ) z i by the average of d times 0.003 End for End for End
Based on the data point approximated by the above algorithm, a parameter fitting for the exponential function bx ae c − + was conducted for each zone as ( ) z i used in the phase-1 explained in Section 2.4.1. The result parameters are given in Table 9 below, Table 9. The parameters fitted by zone Parameter Zone A B C D E F G H I a 2.901 3.276 3.135 3.393 3.150 3.399 2.895 3.333 2.811 b 0.308 0.300 0.265 0.350 0.277 0.348 0.233 0.315 0.273 c 0.969 0.564 0.681 0.480 0.669 0.474 0.933 0.510 1.041 To generate a scenario based on historical data, for each zonal OD pair, the demand volume is generated by randomly choosing one of the historical volumes in the study period. The origin and destination locations are then randomly chosen from the historical ride requests with the same zonal OD without replacement. Next, the detour time zd for each zone z to serve request d are calculated from the distance between the coordinate and zone centroid, and the detour time reductions introduced in Section 2.3 are calculated from the proximity of the service request first given by, ( (min( 81 ) ), ) 3 30.28 z bz zd dz z zd zbdb b x yx y − = − − − + − where zd x , zd y are the x- and y-coordinate of service request location in zone z , and 3830.28 is the maximum Euclidean distance within a zone in meter scale, given the 3x3 zonal division. Finally, all detour time reduction zdb − is normalized by a factor to ensure the sum of detour time reduction is less than the detour time, i.e. satisfying equation (13). In short, this appendix explained how the data are chosen, calculation of demand volume, 0 detour time, detour time reduction, and the approximation of detour time between boundary and the closest service request from real data. The same methodology can be expanded to any set of ride request data with origin and destination location to conduct the ZBFBS analysis. Acknowledgements
The authors acknowledge Didi Chuxing providing ride request data for this study (data source: http://gaia.didichuxing.com). This study is supported by the General Research Fund
References
Aldaihani, M.M., Quadrifoglio, L., Dessouky, M.M., Hall, R., 2004. Network design for a grid hybrid transit service.
Transp. Res. Part A Policy Pract.
38, 511 – Transp. Res. Part B Methodol.
92, 234 – Transp. Res. Part B Methodol.
84, 157 – Transp. Res. Part B Methodol.
66, 70 –
89. Cen, X., Lo, H.K., Li, L., Lee, E., 2018. Modeling electric vehicles adoption for urban commute trips.
Transp. Res. Part B Methodol. – Transp. Sci.
25, 281 – Transp. Res. Part B Methodol. –
54. Chen, P. (Will), Nie, Y. (Marco), 2017b. Connecting e-hailing to mass transit platform: Analysis of relative spatial position.
Transp. Res. Part C Emerg. Technol.
77, 444 – Ann. Oper. Res. –
46. Diana, M., Dessouky, M.M., Xia, N., 2006. A model for the fleet sizing of demand responsive transportation services with time windows.
Transp. Res. Part B Methodol.
40, 651 – J. Mach. Learn. Res.
12, 2121 – Transp. Res. Rec. J. Transp. Res. Board – Transp. Res. Part C Emerg. Technol.
94, 288 – J. Transp. Eng. Part A Syst.
COURSERA Neural Networks Mach. Learn.
Ho, S.C., Haugland, D., 2011. Local search heuristics for the probabilistic dial-a-ride problem.
OR Spectr.
33, 961 – Transp. Res. Part B Methodol. – Transp. Res. Part C Emerg. Technol.
95, 394 – Oper. Res. Lett.
38, 432 – Transp. Res. Part B Methodol.
43, 448 – Transp. Res. Part B Methodol.
68, 76 –
97. Kim, M.E., Schonfeld, P., 2015. Maximizing net benefits for conventional and flexible bus services.
Transp. Res. Part A Policy Pract.
80, 116 – AIP Conf. Proc. –
62. Li, L., Huang, W., Lo, H.K., 2018. Adaptive coordinated traffic control for stochastic demand.
Transp. Res. Part C Emerg. Technol.
88, 31 –
51. Li, X., Quadrifoglio, L., 2009. Optimal Zone Design for Feeder Transit Services.
Transp. Res. Rec. J. Transp. Res. Board – Transp. Policy
39, 63 –
76. Lo, H.K., An, K., Lin, W., 2013. Ferry service network design under demand uncertainty.
Transp. Res. Part E Logist. Transp. Rev.
59, 48 –
70. Lu, W., Quadrifoglio, L., Petrelli, M., 2017. Reliability analysis of centralized versus decentralized zoning strategies for paratransit services.
Transp. Res. Procedia
25, 4100 – J. Adv. Transp.
PLoS One
12, 1 –
28. Martínez, L.M., Viegas, J.M., Eiró, T., 2015. Formulating a New Express Minibus Service Design Problem as a Clustering Problem.
Transp. Sci.
49, 85 –
98. Ng, K.F., Lo, H.K., Cen, X., Lee, E., 2020. Global Optimal Solution Algorithm of Mathematical Programs with Equilibrium Constraints. Working paper. Pan, S., Yu, J., Yang, X., Liu, Y., Zou, N., 2015. Designing a Flexible Feeder Transit System Serving Irregularly Shaped and Gated Communities: Determining Service Area and Feeder Route Planning.
J. Urban Plan. Dev.
Transp. Sci.
40, 351 – Transp. Res. Part B
43, 922 – Biometrics
55, 591 – Expert Syst. Appl.
42, 6728 – Transp. Sci.
29, 17 –
29. Schaller, B., 2018. The New Automobility: Lyft, Uber and the Future of American Cities. 2 Brooklyn, New York City, US. Tong, L. (Carol), Zhou, L., Liu, J., Zhou, X., 2017. Customized bus service design for jointly optimizing passenger-to-vehicle assignment and vehicle routing.
Transp. Res. Part C Emerg. Technol.
85, 451 – Int. J. Urban Sci.
0, 1 –
16. Yang, H., Ke, J., Ye, J., 2018. A universal distribution law of network detour ratios.
Transp. Res. Part C Emerg. Technol.
96, 22 –
37. Zhang, J., Wang, D.Z.W., Meng, M., 2018. Which service is better on a linear travel corridor: Park & ride or on-demand public bus?
Transp. Res. Part A Policy Pract. – J. Transp. Eng. Part a-Systems ––