{'fraction_non_alphanumeric': 0.02379087945100161, 'fraction_numerical': 0.00874183122855475, 'mean_word_length': 4.6607027724049, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Demand-Adaptive Route Planning and Scheduling for Urban Hub-based High-Capacity Mobility-on-Demand Services

In this study, we propose a three-stage framework for the planning and scheduling of high-capacity mobility-on-demand services (e.g., micro transit and flexible transit) at urban activity hubs. The proposed framework consists of (1) the route generation step to and from the activity hub with connectivity to existing transit systems, and (2) the robust route scheduling step which determines the vehicle assignment and route headway under demand uncertainty. Efficient exact and heuristic algorithms are developed for identifying the minimum number of routes that maximize passenger coverage, and a matching scheme is proposed to combine routes to and from the hub into roundtrips optimally. With the generated routes, the robust route scheduling problem is formulated as a two-stage robust optimization problem. Model reformulations are introduced to solve the robust optimization problem into the global optimum. In this regard, the proposed framework presents both algorithmic and analytic solutions for developing the hub-based transit services in response to the varying passenger demand over a short-time period. To validate the effectiveness of the proposed framework, comprehensive numerical experiments are conducted for planning the HHMoD services at the JFK airport in New York City (NYC). The results show the superior performance of the proposed route generation algorithm to maximize the citywide coverage more efficiently. The results also demonstrate the cost-effectiveness of the robust route schedules under normal demand conditions and against worst-case-oriented realizations of passenger demand.


Introduction
High-capacity public transit such as bus and metro serves as an affordable mobility solution to urban commuters. When properly planned and operated, public transit may significantly reduce commuters' dependence on private vehicles and play a vital role in a sustainable mobility system in dense populated urban areas [1]. Unfortunately, the configuration and infrastructure of existing fixed-route public transit systems may no longer be attractive to urban travelers due to their fast-changing mobility needs. One notable evidence is the recent rise of the ride-hailing industry and subsequent loss of public transportation ridership. For instance, in New York City (NYC), the bus ridership declined by 1.3%, 5.1%, and 5.8% from years 2016 to 2018 [2] while the number of for-hire vehicles (FHV) trips had increased by 300% during the time [3]. This leads to more congestion, emissions, and energy consumption [4] in cities, and poses a significant need for a better balance and connectivity between emerging mobility options and public transportation.
One challenge associated with the planning of the transit system is the trade-off between accessibility and service coverage. It is often impractical for high-capacity urban transit such as fixed-route buses to provide door-to-door service to meet the varying mobility needs of the general public. Nevertheless, with flexible transit of adaptive service routes and schedules, it is possible to offer a more effective transit service for serving passengers sharing similar origins and destinations [5]. One appealing and applicable scenario is the public transit rooted in major urban activity hubs such as shopping malls, financial districts, railway stations, and airports. These locations have high passenger volume, and it is likely to find travelers to and from these hubs of similar travel routes. We show in Figure 1 the empirical evidence that supports the applicability of hub-based transit service, where the cumulative pickups and drop-offs of both taxis and FHVs are plotted at the taxi-zone level. In particular, we can verify that around 20% of the taxi zones account for over 60% of total pickups and drop-offs for both taxis and FHV. Moreover, we observe from the data that over 25% of yellow taxi trips and 15% of FHV trips have either the pickup or drop-off at one of the airports in NYC during peak hours. These observations reveal huge opportunities for travel mileage reduction if individual trips can be optimally consolidated, and motivate us to investigate the design of the hub-based high-capacity mobilityon-demand (HHMoD) service as a demand-adaptive transit system at urban activity hubs. The HHMoD is a special case of the ridesharing scheme, and the main idea is to satisfy passenger demand at activity hubs with high capacity vehicles that use optimized routes with minimum deviation from the shortest travel routes. But different from the ad-hoc ridesharing with fewer passengers, the high capacity feature and the higher and more stable passenger demand at activity hubs allow the routes and schedules of the HHMoD to be planned proactively. This avoids the unnecessary waiting for the realization of requests from many passengers and is also more attractive to potential passengers as they will be informed of the scheduled services ahead of time. And the routes and schedules of the HHMoD are adaptive since they can be reoptimized from time to time (e.g., every 30 minutes or one hour) based on the predicted future demand to ensure the planned services are flexible and aligned with the time-varying mobility needs.  To promote the idea of the HHMoD, this study concerns the modeling framework for the route planning and scheduling of the HHMoD over a short time period, which can be regarded as a novel variant of the classical transit route network design problem (TRNDP). There exists broad literature on TRNDP and readers may refer to [6,7,8] for comprehensive reviews of related studies. In general, the TRNDP is intrinsically difficult and the analytical framework cannot be applied to solve real-world problems. In this regard, heuristic approaches constitute the mainstream for solving TRNDP, which compromises sub-problems of the design of operation routes and the scheduling of transit service [7]. There are two main directions to tackle the sub-problems. The first direction divides the process of transit planning into route generation and headway optimization. Most of the TRNDP studies followed this direction and various methods for route generation were proposed. For instance, Lampkin and Saalmans [9] proposed the route generation method that starts with random control points and consecutively adds stops between these points. Silman et al. [10] divided the city into zones and used a heuristic algorithm for generating routes that minimize walking time. Recently, Cipriano et al. [11] used Genetic Algorithms (GA) to generate candidate route sets, which were evolved according to the fitness measure of travel time. Nikolic and Teodorovic [12] developed the route generation process based on Bee Colony Optimization, where random routes got improved by modeling each route as a bee and modifying its route with a pheromone objective function. The second direction follows an iterative route generation and improvement approach where route and headway are jointly determined in an iterative process. Studies that follow the second direction include Drezner et al. [13], Fan and Wei [14], Chakroborty and Partha [15], where GA was widely adopted in these studies to create routes and evolve the routes over iterations, and travel links or control points were chosen as genes which may change during the search of an optimal solution. It can be seen from the literature that, due to the difficulties in solving TRNDP exactly, heuristic and metaheuristic approaches are heavily used to obtain the candidate route set and determine the frequency setting. Nevertheless, there is no guarantee on the solution quality with these approaches and the complicated parameter tuning is often required in order to obtain a satisfactory solution. Moreover, even for a finetuned heuristic method, the performance may vary significantly with respect to different passenger demand realizations, which creates a critical barrier for using such a framework to solve the time-varying flexible route transit design problem. We also note that there are several recent studies that investigated the design of flexible route transit and the demand-adaptive transit, such as the simulation model for demand responsive transit design [16], the demand-responsive system with compulsory and optional stops [17], the flexible transit for low demand areas in grid networks [18], the semi-flexible transit systems [19], the feeder transit with fixed and flexible routes [20], the point-to-point high coverage transit system [21], and the paired-line hybrid transit with radial route structure [22,23] and the flexible line length bus system [24]. There exist three major gaps based on the results from the previous studies. First, as the flexible route service was primarily proposed for sparsely populated areas, the effectiveness of the existing framework remains unknown if implemented in populated urban areas. Second, existing studies rely on the simulation framework [16,21], heuristics [20] or metaheuristics [22] for obtaining the optimal design due to the high computational cost. Finally, the demand-responsive mobility services were mainly investigated under deterministic or revealed demand. But an effective demand-responsive system needs to account for future passenger realizations under uncertainty. To the best of our knowledge, an analytical framework for the route generation and robust route scheduling of urban mobility-on-demand services has yet been established.
In light of the fundamental challenges and the practical usefulness of the HHMoD for serving the mobility needs in populated urban areas, this study establishes the framework for operating HHMoD that is tailored to the spatiotemporally varying mobility needs of urban travelers based on the real-time trip and traffic data. The primary goal of the planned HHMoD is to find the set of routes that maximize passenger coverage and identify the optimal schedules of the routes that minimize the worst-case oriented operation cost under demand uncertainties. Specifically, we take advantage of the hub-based structure and propose both exact and heuristic algorithms that generate the round-trip HHMoD routes with maximum passenger coverage. The solution algorithms are further extended to incorporate connections to existing public transportation systems. Given the candidate routes, a two-stage robust optimization model is then developed to schedule the vehicle assignment and route headway optimally. The computational performances and the effectiveness of the proposed HHMoD solution are demonstrated with large-scale numerical experiments at the JFK airport in NYC with 200 candidate stops and connectivity to the NYC subway system.
The rest of the study is organized as follows. The second section gives an overview of the research framework for developing the HHMoD system. The third section discusses the route generation problems and proposes exact and heuristic algorithms for obtaining the candidate route set. A heuristic shortestpath based route generation mechanism in the literature is also introduced as the benchmark algorithm. Afterward, we present the route combination model for joining planned routes to and from the hub. The fourth section develops the two-stage robust optimization model in light of the worst passenger realization for the optimal scheduling given the generated route set. The fifth presents comprehensive numerical experiments on planning the HHMoD service at JFK airport in NYC. Finally, we conclude our study with major findings and future directions in the last section.

Problem overview
In this study, we follow the transit planning framework as summarized in [7] by first generating the set of candidate routes and then developing the optimal service frequency under demand uncertainty. We propose a three-stage framework for developing the HHMoD system, as shown in Figure 2. Our framework is different from the route generation step in existing studies -we propose an innovative route generation algorithm that is objective-driven (maximizing passenger coverage). One may choose to find a predefined number of routes with maximum passenger coverage or find the minimum number of routes that cover all passenger demand. Besides, the route generation algorithm offers additional flexibility for connecting the HHMoD system to existing transportation facilities. The overall framework takes the set of candidate stop locations, the trip demand level, and road traffic conditions as the inputs. It outputs the set of operating routes, the assigned number of vehicles to each route, and each route's corresponding operating frequency (headway). Here we  Figure 2: HHMoD Planning framework briefly summarize the three major steps for the development of the HHMoD system: 1. In the first stage, we develop the route generation algorithm which builds the set of K feasible candidate routes that maximize the coverage of potential passengers, where K is the user-specified number. A route is considered as feasible if the trip time of all passengers is no larger than λ times the minimum trip time by taking alternative travel modes. Given the candidate stops, trip demand and road traffic condition, the route generation algorithm will first create two sets of feasible routes for incoming and outgoing trips at activity hubs respectively.
2. The second step then identifies the optimal combination of the two sets of routes to and from the hub. The route combination step focuses on balancing the passenger demand between two directions and is subject to the vacant trip distance constraint. The constraint helps to avoid excessive empty miles for connecting two distant routes and will force the vehicle to return to the hub if there are no matching routes in the opposite direction within the distance constraint.
3. The integrated round trip routes will then serve as the input for the third stage problem, which is to optimize the vehicle assignment and route service frequency under demand uncertainty. For the route scheduling optimization problem, our objective is to minimize the weighted combination of the total operation cost, passenger waiting time together with the penalty for unsatisfied passenger demand. A two-stage robust optimization problem is developed to identify the worst-case oriented optimal fleet and headway setting.
As the problem requires the information on candidate stop locations and the associated demand level, such knowledge can be directly mined from large-scale GPS data [25]. And the real-time travel time information can also be obtained from various sources such as Google Maps. With the availability of this information, the proposed method can be used to deploy the HHMoD service over a short-time period (e.g., 30 minutes or an hour interval) based on predicted demand or the distribution of future passenger demand from historical observations. This offers great flexibility to align the available fleet resources with the passengers' timevarying mobility needs, and the framework may, therefore, advance the planning and operation of urban MoD with high capacity vehicles. We next present the details of the three-stage framework.

Route generation 3.1 Problem description
As the first step of the HHMoD development, the route generation problem focuses on identifying the set of K candidate routes that yield the highest possible passenger coverage. We term this target set of K routes as the K maximum coverage routes (K-MCR). We consider the study area consists of |N | candidate stops and we denote G(N, E, Q, W, λ) as the network representation of the study area, with N being the set of stops, e i,j ∈ E denoting the shortest travel path between stop i and stop j, q i ∈ Q being the potential passenger demand (to hub or from hub) at stop i, and w i,j ∈ W being the weight (shortest travel time) on edge e i,j . We make the following reasonable assumptions to conceptualize the routing generation problem: 1. We assume that the HHMoD is a station-based system where the candidate stops and the associated demand level at each stop is given.
2. The travel time between two stops are known (e.g., obtained from historical or real-time traffic data).
3. We consider that passengers arriving at station i are covered by the HHMoD if there exists a route that stops at i.

4.
The HHMoD service will use the shortest travel time path when travelling between two stops.
5. The shortest travel time on different segments between stations should satisfy the triangular inequality constraint.
6. Passengers using the HHMoD are time sensitive.
The first four assumptions set the scope of the route generation problem with candidate stop locations, known passenger demand and travel time information. The fifth assumption imposes a practical constraint on travel time, t i,j + t j,k > t i,k , so that the route travel time should strictly increase with additional intermediate stops.
And the sixth assumption suggests that passengers may use alternate travel mode with trip time t m (such as subway and taxis) if the travel time of the proposed HHMoD service exceeds a certain threshold. This motivates us to take the route travel time into consideration while deriving the MCR, where we can translate the time constraints for passengers who arrive at stop i using route k as: so that the planned route k ∈ R K for passengers at stop i should result in the travel time no greater than λ times that of the alternative modes of lowest travel time. As an example, with λ = 1.5, passengers will not use the planned HHMoD service if the travel time is 50% higher than the travel time of hailing a taxi. In this regard, λ is a hyper-parameter of the route generation model and can be adjusted to account for the willingness to wait for potential passengers. Based on the above constraints, the K-MCR generation problem can be mathematically formulated as: where δ k i and s k are both indicator variables. s k = 1 if route k is selected for the K-MCR and 0 otherwise. δ k i takes the value of 1 if stop i is served by route k. Problem 2 is NP-hard. A straightforward proof is that we can reduce the longest path problem, a known NP-complete problem, to Problem 2 under K = 1 by setting λ to be an arbitrary large value. With such a λ, all travel time constraints are directly satisfied and K = 1 is a special case of the original problem. The problem becomes more difficult with active travel time constraints and K > 1. And it is impractical to solve it as an optimization problem with commercial solvers as the explicit formulation may require the enumeration of all candidate routes. In the following sections, we establish exact and heuristic solution algorithms to solve the K-MCR generation problem.

Exact algorithm
While Problem 2 represents a difficult combinatorial problem, we can still solve it efficiently by utilizing the topological properties of the network and specific hub-based nature of the problem.
First, denote R K as an optimal solution of problem 2 and the route r i ∈ R K as the ordered set of stops it traverses through, e.g. r i = {n 1 , n 2 , n 3 , ...n Li } and n 1 , n 2 , ...n Li ∈ N . Based on Assumption 3, R K can be reduced to R K where for any r i , r j ∈ R K , we have r i r j = ∅. And we term R K as the mutually disjoint K-MCR. To see this, let r i r j = S = ∅ and we may remove the set of nodes in S from either r i or r j . This step will not affect the coverage of passengers as the removed stops from one route is still served by the other route. In addition, the route after the stop removal is still feasible since the route travel time strictly decreases with the removal of an intermediate stop following the triangular inequality in Assumption 5. Based on this property, we have the subsequent proposition that allows for solving the K-MCR generation problem by sequentially identifying the highest passenger coverage route and generating the mutually disjoint MCRs. Proposition 1. Let r 1 , r 2 , ..., r K be the set of mutually disjoint routes generated as sequentially identified MCR, then R K = {r 1 , r 2 , ..., r K } is the K-MCR set that gives the highest passenger coverage.
Proof. The proposition can be proved by induction. First, we verify that the proposition holds for K = 1. Assuming that R N −1 is the (N-1)-MCR, and r N is the mutually disjoint MCR identified at step N that shares no overlapping stops with the routes in R N −1 . Let us consider that R N −1 {r N } is not the N-MCR, then there must be a route r = r N that yields higher passenger coverage than r N . But this contradicts the fact that r N is a mutually disjoint MCR. Hence R N {r N } gives the N-MCR.
On the other hand, finding the mutually disjoint route with maximum passenger coverage itself is still nontrivial. In the worst case, one may have to enumerate all possible routes in the network and select the route with the highest passenger coverage. However, the time constraint in Problem 2 introduces a promising direction for reducing the search space for finding MCR. In particular, starting from a stop i, the MCR can be obtained by tracking the potential passenger coverage of the set of stops that are reachable to and from the hub. And at stop i, we define potential passenger coverage of a stop j as the largest possible passenger demand that can be covered by extending the route from stop i to stop j. Let Γ(i) denote the operator for the maximum potential passenger coverage of a route starting from location i, we can express the MCR generation as a dynamic programming problem following: In the equation, Reach(i) represents the set of stops that are reachable from stop i, where j is said to be reachable from i if the travel time of the extended route from i to j will not exceed λ times the corresponding shortest travel time as shown in Problem 2. Instead of enumerating all routes and select the best one, we can therefore build the solution algorithm for sequentially identifying mutually disjoint MCR following equation 3. We note that the MCR generation can still be expensive with the worst-case complexity of O((e − 1)N !) as shown in Proposition 2. However, in practice, the Reach(i) may be significantly smaller in each step due to the travel time constraints and the exact MCR can still be generated efficiently with reasonable values of λ (e.g., λ <= 1.4 which represents a 40% deviation from the fastest alternate mode). We also observe the exact algorithm to generate K-MCR efficiently in our NYC case study. As a consequence, we can tackle the K-MCR problem by utilizing the hub-based structure and the travel time constraints, and the complete algorithm with details for the potential calculation and reachable set identification are summarized in Algorithm 1. The algorithm starts with all stops in the study area, finds the MCR by constructing the recursion tree based on the reachable set, and then sets the potential of the stops that are covered in the MCR to 0. This ensures that the routes generated each step are mutually disjoint MCRs. The algorithm terminates when all K routes are generated or can be terminated early if all stops have been visited. And the algorithm can be applied for passengers to and from the hubs separately, with the only difference being that the corresponding travel cost w i,j to be switched to w j,i and that the t m j to be set to minimal travel time to and from the hub respectively. Proposition 2. The worst case time complexity for identifying the MCR with N stops is O((e − 1)N !).
Proof. The complexity can be shown based on equation 3. If we denote the computation requirement for N nodes as f (N ), then we have where c is small constant cost at each step, and the equation can be generalized to Consequently we have Algorithm 1 Hub-based Route Generation Algorithm Input: G = (N, E, Q, W, λ): network of study area, K: number of routes to be generated, h: hub location, λ: deviation threshold Output: R K the set of candidate operation routes.

Heuristic algorithm
In the case when the exact algorithm may not find the solution in a timely manner, such as when passengers are insensitive to the travel time and a high λ value is required for certain real-world applications, we can further develop a heuristic solution approach to identify the near optimal K-MCR based on the following two properties for the route generation problem.
Proposition 3. The induced reachable set of a stop is a subset of the reachable set of its predecessor stop.
Proposition 4. The potential passenger coverage of a stop i is no larger than the sum of the passenger demand of the stops in its reachable set Reach(i).
Proposition 3 holds directly from the triangular inequality constraint for travel time, as in Assumption 5, which states that the set of stops that can not be reached from the predecessor of stop i is also not reachable from stop i. Proposition 4 sets the upper bound of the maximum possible potential of a stop since a HHMoD route can not serve more than the number of stops in the current reachable set. Based on these two propositions, we can relax the Γ operator in equation 3 by replacing it with a heuristic operator Γ h as the sum of the node potential in its reachable set: This step is equivalent to relaxing the M axP function in Algorithm 1 by replacing the recursion with a heuristic M axP function that calculates the sum of the potential of stops in the reachable set. In this regard, the heuristic approach represents a greedy approach in building the MCR route, and the computation time of the heuristic method is O(N 2 ). This allows us to find near-optimal K-MCR for large instances in polynomial time. We summarize the details of the heuristic M axP function in Algorithm 2.

Algorithm 2 Heuristic MaxP
end for 10: return P, r 15: end function

Route generation with connections to existing transportation systems
The effectiveness of the HHMoD can be further improved if the planned service can be aligned with the existing transportation facilities such as the subway and buses. This is especially appealing since existing transit facilities have already established a wide coverage of urban networks. In this regard, the HHMoD may avoid producing overlapping travel segments, and in many cases, the travel time using the subway can be shorter than the ground traffic. And it is likely to realize a higher passenger coverage with the same number of routes. In the following paper, we term HHMoD connect as the HHMoD service whose routes are planned by incorporating the connections to existing transit facilities, and we use HHMoD only to denote the planning of the HHMoD without such connections.
With the K-MCR generation process described in the previous sections, the framework can be directly extended to embed the connection to transit facilities by modifying the MaxP and Reach functions. And additional inputs are required to describe the costs of such connections. Let S be the set of stops of the existing transportation facilities. We denote d i,j as the connection cost between stop i of the HHMoD and the stop j ∈ S and t i,j as the known travel time for i, j ∈ S. We can then define that a transfer at stop i for the passengers who ride the HHMoD to stop k is feasible if: Equation 9 states that the summation of travel time for (1) riding HHMoD from the hub to stop i (or from stop i to the hub), (2) walking distance from i to j, (3) an inconvenience cost σ such as the extra waiting time, (4) the trip time from j to j , and (5) the walking distance from j to stop k should be no greater than the relaxed shortest travel time λ min m t m k . Such a definition is similar to that in equation 1, and subsequently, this constraint can be used to guide the search of the reachable set for the HHMoD connect problem. We note, however, that only one transfer is considered in our problem to avoid unnecessary complications. While considering multiple transfers may further the effectiveness of the route generation process, the improvement is likely to be minor as the addition of extra inconvenience costs will render most of the routes with multiple transfers infeasible.
Given the transfer costs at each of the HHMoD stops, we can formulate the dynamic programming problem for the HHMoD connect (with Γ C (·) as the max coverage operator) as: Different from equation 3, we now have both CReach(i) and HReach(i) at stop i, with CReach(i) being the set of stops that can be reached from stop i through the connections to the transportation systems. As compared to HHMoD only, we have Reach(i) ⊆ CReach(i) HReach(i) and This ensures that the K-MCR of HHMoD connect will always achieve greater or equal passenger coverage as compared to the K-MCR of HHMoD only with the same K. Nevertheless, we should also note that equation 10 does not constitute an exact solution when considering the connections as we greedily add the potential passenger coverage through connections at the first stop when such a connection is feasible. And the exact solution for HHMoD connect should be obtained by constructing the additional recursions on if each connection (i, k) should be established ∀ k ∈ CReach(i). This will result in 2 |CReach(i)| combinations of max j∈HReach(i) Γ C (j) for each step and will quickly render the problem intractable even for small λ values. As such, we limit our discussion to the greedy connection scheme as in equation 10 and the corresponding K-MCR generation algorithm for HHMoD connect is summarized in Algorithm 3. In addition, the same modification can also be extended for the heuristic algorithms of HHMoD only problem to incorporate the connections to existing transportation systems as following:

Shortest-path based route generation
In addition to the developed algorithm for MCR generation, we also introduce a shortest-path based heuristic route generation algorithm in [26] as the benchmark for comparison. The main idea is to first generate a large collection of candidate routes and then prune routes based on user-specified thresholds. We summarize the details of the benchmark heuristic in the Appendix.

Route combination
By applying the route generation algorithms for both passengers to and from the hub, we obtain R d K as the MCR for serving passengers traveling to the hub and R o K as the MCR for serving passengers traveling from the hub. The next step is to connect R d K and R o K to form round routes and therefore reduce the amount of vacant mileage for the HHMoD fleet.

Algorithm 3 Hub-based Route Generation Algorithm with Connections to Transportation Systems
Input: G = (N, E, Q, W, λ): network of study area, K: number of routes to be generated, h: hub location, D: connection costs, T : travel time in transportation systems, S: the set of stations of transportation systems, σ: inconvenience cost for transfer, λ deviation threshold Output: R K the set of candidate operation routes.
1: for i=1,2,...,K do 2: Add r i to R K 7: end for The route combination is conducted following two rules. First, the combination of routes to and from hubs with similar demand coverage should be prioritized, which avoids the unbalanced demand issue and makes the best use of the vehicle capacity. Second, the connecting distance between the end of two routes should not exceed the minimum of the trip distance of the two routes. This helps to avoid the combination of two routes of similar demand level but is distant from each other, where directly returning to the hub is the most economical solution. Based on these two criteria, we formulate a bipartite matching problem for the optimal combination of routes in R o K and R d K . In particular, we denote g r1,r2 as the gap between the demand level of route r 1 ∈ R o K and r 2 ∈ R d K . If the connecting travel time exceeds the shorter travel time of the two routes, g r1,r2 is set to an arbitrary large value to represent infeasible combination. Consequently, the optimal route combination can be obtained by solving the following optimization problem: The problem is equivalent to the minimum weight perfect bipartite matching problem and can be solved efficiently using the Hungarian algorithm. Readers may refer to [27] for implementation details.

Route scheduling under demand uncertainty
Given the set of combined candidate routes R, we next introduce a different set of notations as summarized in Table 1 to formulate the route scheduling problem. This is different from the planning of conventional public transportation services. The users of the HHMoD service are usually more time-sensitive and pay a premium for better mobility that offers shorter travel time and better riding experience, and more importantly, more reliable services. The HHMoD service is vulnerable to the loss of its market if it fails to satisfy the revealed demand that is different from the daily expected levels. In this regard, it is crucial to plan the service in the face of demand uncertainty and ensure the reliability of the scheduled service even under the worst realization of passenger arrivals at the stops. This motivates us to investigate the robust optimization for the route scheduling and identify the optimal strategy that is resilient to worst-case demand realizations.
Here, for a given time period (e.g., 30 minutes or an hour), the route scheduling problem is to determine the allocation of a fleet of B vehicles to the set of candidate routes R and the corresponding frequency of each route, termed by the headway h s , with the objective to minimize the weighted combination of the total operation cost, the waiting time cost for served passengers and the penalty for unsatisfied passenger demand. We assume the Poisson arrival of passengers at each stop and that the fleet is composed of homogeneous vehicles of the same capacity C. We consider the headway of each route is scheduled at the per-minute level, and therefore, the headway setting is modeled as an integer variable. We follow the framework in [28,29] to develop the mathematical formulation for the HHMoD frequency problem. Before accounting for the demand uncertainty, we first summarize the nominal model (deterministic demand) for the optimal route scheduling problem with individual components in the objective function can be expressed as Note that for a passenger waiting at the stop, the expected waiting time is equal to half of the vehicle route headway h s as we assume the random arrival of passengers [29]. In equation 13c, the travel distance d i by alternate mode is also included to penalize the loss of services. For instance, we may penalize longer trips more than shorter trips if the HHMoD may fail to serve these trips, which are instead served by taxis and FHVs. With the objective functions, the nominal HHMoD scheduling problem can be formally written as min y,h,κ,X,L F operation + F waiting + F loss (Nominal HHMoD Scheduling) where H describes the set of constraints that regulate the number of vehicles y s and headway setting h s for the planned HHMoD service and F characterizes the set of passengers that can be served by the scheduled services. Equation (14a) establishes the relationship between assigned vehicle number y s , vehicle headway h s and round trip travel time T s . For instance, if the the round trip travel time is T s = 36 minutes and the required time headway is h s = 12 minutes, then we need at least y s = 3 buses. Equation (14b) implies that y s is a positive integer upper-bounded by B when κ s = 1. If κ s = 0, then y s = 0. Equation (14c) shows the upper and lower bound of the time headway h s . Equation (14d) represents that the total number of assigned vehicle do not exceed to vehicle number constraint B. Besides, equations (15a) and (15b) imply that the total boarding passengers for routes from and to activity hubs should not exceed the vehicle capacity C. Equations (15c) and (15d) describe that the sum of satisfied and unsatisfied demand is equal to the total demandD F i andD T i in the city. Finally, equations (15e) and (15f) represent the served demand at stop i does not exceed the total demand of stop i. The nominal problem aims at identifying the optimal service headway h s and assigned vehicles y s for each route s ∈ R with the assignment of served and unsatisfied passengers. Nominal HHMoD Scheduling is a nonlinear and nonconvex mixed-integer programming problem, which is difficult to solve directly due to the existence of bilinear constraints y s h s and bilinear terms h s X F s,i δ s,i , h s X T s,i δ s,i in the objective function. We note that the Nominal HHMoD Scheduling considers known passenger demandD T i andD F i as a point estimation (such as expected value). Nevertheless, for the proactive planning of HHMoD service in real-time, the passenger demand may vary significantly, and the planned service based on expected passenger demand D T i andD F i is unlikely to align with the actual realization of demand profile. This may lead to a situation where there is overly supply on certain routes, with the revealed demand being significantly lower than the scheduled capacity. Meanwhile, the oversupply in certain routes may result in the shortage of supply to serve the actual demand on other routes. In this regard, if the uncertainty set of passenger demand is known, the nominal scheduling problem can be converted into a two-stage robust route scheduling problem to negate the impacts from the mismatch between scheduled supply and realized demand. Specifically, in the first stage, the master problem represents the problem of scheduling vehicle allocation and route headway with known number of passengers. As for the second stage, we have the recourse problem where the worst-caseoriented passenger demand distribution is realized in response to the master problem's route configuration. We consider the box constraint to capture the demand uncertainty as whereD i is the expected demand level and Z i is the maximum demand deviation at stop i. And we also have: To control the conservatism of the robust optimization outcomes, the budget for demand uncertainty is introduced as in [30] where the decision-maker may decide to account for the maximum deviation of passenger demand for up to Γ locations with With the above defined demand uncertainty set, the two-stage robust optimization problem (RO) can be formulated as: where F(p) follows from replacing the constraints 15c-15f by equations 16-19.

Recourse problem reformulation
As the relatively complete recourse property always holds for our recourse problem, and that the inner minimization problem is a linear programming problem if p is supplied, we can convert the max-min recourse problem RO-HHMoD Scheduling into a single maximization problem by utilizing the strong duality condition with known h * and y * from the master problem: With u, v, λ and p being the decision variables, the converted recourse problem now becomes a bilinear integer programming problem with the bilinear terms λ i p i and v i p i . As p F i and p T i are binary variables, instead of solving a nonconvex problem, we can further convert this problem into an equivalent mixed integer linear programming problem (MILP) by making use of the big-M method to replacing each p and Problem 22 can therefore be solved to global optimum using off-the-shelf commercial solvers.

Master problem reformulation
We can obtain the optimal value of p * by solving the reformulated recourse problem, which allows us to rewrite an equivalent master problem by introducing the auxiliary variable η as: To solve the master problem, the major challenge arises from the bilinear term h s X s,i in the objective function and the constraints, and the bilinear term y s h s in the constraint, where both y s and h s are integer variables. As both bilinear terms share the common elements h s as an integer variable, we introduce the unary expansion to express h s as the summation of a series of binary variables g k s as: with H = log 2 (h max + 1) − 1. We therefore have and we can further replace g k s y s with z k s and derive the boundary for z k s based on the McCormick's envelope: We note that the McCormick envelope constructed for z k s is an exact relaxation due to the binary nature of g k s , and hence equations equations (25) and equations (26) result in an exact relaxation of the binary constraint y s h s . Similarly, we can convert the bilinear terms h s X i with g k s X i and replace each g k s X i with q k s,i as: The exact relaxations in equations 26 and 27 convert the bilinear master problem into an equivalent MILP problem so that the master problem can also be solved to global optimality using the off-the-shelf MILP solvers.

Solving the two-stage robust route scheduling problem
With the reformulation of both the master and the recourse problems, we are now ready to tackle the twostage robust route scheduling problem. It is in general difficult to solve the two-stage RO problems due to its multilevel optimization structure, and there are in general two widely-adopted solution approaches in the literature: the extension of Bender's decomposition method [31] and the Column and Constraint Generation (C&CG) method [32]. In this study, we adopt the C&CG method as it has been shown to achieve superior convergence performances than the Bender's decomposition method for practical sized problems (e.g., hundreds of times faster in [33] with guaranteed algorithm convergence as both our reformulated master and recourse problems can be solved to global optimum. The C&CG algorithm solves the two-stage RO by iterating over the master and recourse problems, where optimality cuts associated with primal variables are added to the master problem until convergence.
With the master problem above, we summarize the major steps for implementing the C&CG algorithm to solve the two-stage robust route scheduling problem as following: Step 1: Initialization: Set LB = − inf, U B = inf. Set N = 0, set T = ∅.
Step 2: Find an initially feasible p 0 ∈ P, add p 0 to T. Set N = 1.
Step 5: If U B − LB ≤ , terminate and return y * and h * as the optimal solution. Otherwise update N = N + 1, add p N to T and go to Step 3.

As for
Step 2, an initially feasible p 0 can be easily selected by randomly select Γ stops and set the corresponding p i to 1. In this study, the Gurobi 9.0 optimizer is used to solve the MILP for both master and recourse problems. We note that, with the unary expansion and the relaxation of the bilinear terms, the C&CG algorithm will add 2|N |(H + 1) continuous variables and 2|R| + 8|N | + 1 constraints into the reformulated master problem. While this will quickly enlarge the size of the master problems, the additional iterations do not increase the number of integer and binary variables and fast algorithm convergence is still observed through our numerical experiments where the algorithm will usually terminate within a few iterations.

Experiments setting
We choose NYC as the study area and demonstrate the effectiveness of the proposed HHMoD framework by developing the HHMoD service at the JFK airport, which is associated with the origins and destinations of over 3.5% of total daily FHV and taxi demand in 2018. This section presents three sets of experiments: (1) the computing performances of the route generation and robust planning models, (2) the effectiveness of the route generation approaches, and (3) the resulting costs of the fleet arrangements through robust optimization. We prepare four demand scenarios during weekdays for both the HHMoD only services and HHMoD connect services: AM peak scenario from 7:00-9:00, off-peak scenario from 13:00-15:00, PM peak scenario from 17:00-19:00, and night time scenario from 21:00-23:00.
To prepare model inputs and quantify the modeling parameters, we collect information from publicly available datasets, including road information, geographical subdivisions of the city, public transit information, and demand of taxi and for-hire vehicles (FHV). The network of our study area is built from the taxi zone shapefile of NYC, which has 263 taxi zones across five major boroughs. State Island is not included due to a low demand level (less than 0.5% of total demand to hubs such as LGA, JFK, and Penn Station). Due to the lack of pick up and drop off locations, we assume the centroid of each zone being the stops for the HHMoD service and we use the combination of FHV and taxi demand to and from the JFK airport in each zone as the demand for each stop based on the 2018 trip data [34]. A more realistic scenario can be prepared by conducting spatial clustering to determine the number of stops, and the associated demand level when the pickup and drop off locations are available. We only select the top 100 highest demand stops to and from the hubs, which account for over 85% of the total demand during each time intervals. This selection results in the K-MCR problem with 100 stops each way and the robust optimization problems with the budget for demand uncertainty for up to 200 stops. We calibrate the average demand for each stop as well as the maximum deviation at each stop using the 2018 NYCTLC taxi and FHV trip data [34]. The shortest travel time between stops is calibrated by the average FHV and taxi trip time between the two areas. For HHMoD connect, we consider the planned HHMoD being connected to the NYC subway system of 36 lines and 247 stations. The walking distance is measured between each HHMoD stop and subway station and then converted to the corresponding walking time. The trip time between subway stations is obtained from the NYC subway timetable, and an inconvenience cost of 500 seconds is added to all trips that require a transfer between the HHMoD and the subway system in light of the additional waiting and delay. Finally, we summarize the modeling parameters that are required for the route generation and fleet scheduling problems in Table 2.

Route generation
We first demonstrate the computational performances for the route generation, and the results are summarized in Figure 3, where we vary the time threshold value λ and compare how the computational time scales with increasing time threshold. An inconvenience cost of 500 seconds is added to the HHMoD connect approach. And the results reported here are for the computational time required to cover all passenger demand in the study area. For both HHMoD connect and HHMoD only cases, the heuristic methods can generate the K-MCR in only a few seconds. Whereas in the exact method, the computational time for route generation scales exponentially with increasing time threshold since a higher threshold will result in more number of accessible nodes in the reachable set at each stop. As for HHMoD connect, we observe that it is much more efficient than the HHMoD only scenario, though the computational time still scales exponentially. The reason is that, by connecting to the subway system, certain candidate stops in the reachable set can be directly included in the final solution as long as they can be reached through a transfer to the existing transit system. This helps to reduce the number of nodes to be further branched on, thus saving computational time. However, the extent of the reduction that can be achieved may vary depending on the time threshold of λ, the inconvenience cost for transfers, and also the layout of the existing transit systems. In Figure 3c, we also visualize how the computational time may change with inconvenience cost ranging from 100 seconds to 1200 seconds while holding the λ = 1.3. We observe that the increase in computation time for HHMoD connect has a superlinear relationship concerning the increase in inconvenience cost until it reaches the value of around 900 seconds. Beyond this point, the transfer to the subway system is no longer viable, and the computational time is equivalent to that of the HHMoD only case. On the other hand, we report that it will take more than 30 minutes for the shortest path based heuristic to generate the set of candidate paths, and the generated path is unable to fulfill all the passenger demand in the study area (as can be seen in Figure 4). In this regard, the shortest-path-based heuristic may not apply to the route generation problem in our case. Instead, the developed exact approaches can generate high-quality solutions efficiently with reasonable λ values. In the case where a high λ is required, one can use the developed heuristic approaches to generate the K-MCR greedily in real-time. We next present the effectiveness of the route coverage algorithms for the four different times of the day and the results are summarized in Figure 4. The figure visualizes the cumulative passenger coverage rate with an increasing number of identified routes. We find that all exact and heuristic approaches in our study have superior performances compared to the shortest path based heuristic method. While the shortest path based heuristic may achieve similar passenger coverage for the first few generated routes, there is a significant gap in passenger coverage with additional routes generated. Specifically, the shortest path heuristic plateau at the passenger coverage of around 70% and is unable to reach 100% passenger coverage even with a large number of routes generated. In contrast, the exact algorithms can reach 100% passenger coverage with around 32 candidate routes, whereas the heuristic algorithm can also cover all the stops with 38 to 40 routes. And all exact and heuristic algorithms are found to deliver consistent performances with varying passenger demand profiles at different times of the day.  By varying the λ value, we can further examine the quality of the generated routes based on (1) the demand coverage level of the top candidate routes and (2) the number of routes required to cover all passenger demand. The first indicator is important for decision-makers who seek to make the best use of their limited vehicle resources and the second indicator suggests how decision-makers may maximize the service coverage with minimal cost. Here we only demonstrate the results for the routes generated for passengers departed from the hubs during AM and PM peak periods, and similar performances are observed for other scenarios and for passengers traveling to the hubs. Figures 5 (a)-(b) show the total passenger coverage of the top 5 routes with a time threshold value between 1 and 1.4. It can be observed that the exact algorithms for HHMoD only and HHMoD connect have identical performances with a small time threshold value (e.g., λ ≤ 1.15). In these cases, the additional travel time is not sufficient to compensate for the transfer and walking time to and from the subway system and none of the generated paths is connected to the NYC subway. As the time threshold increases, the top 5 routes may cover up to 70% of the passenger demand by connecting the HHMoD service to the subway system, where the corresponding coverage without such a connection is between 44% to 55% depending on the particular time of the day. The results illustrate that a sufficiently high passenger coverage level can be reached with a few routes by the proposed route generation approach, especially during day time, which only requires a reasonable relaxation of the travel time as compared to the direct service from taxis and FHVs. These results highlight the effectiveness of the HHMoD service in serving the passengers at activity hubs and also suggest notable benefits for integrating the HHMoD service with existing transportation facilities. This observation at JFK supports our initial speculation where the HHMoD can be a promising ridesharing solution to consolidates passenger to and from major activity hubs. In addition, we observe that the performances of the heuristic algorithms are 5% to 10% inferior as compared to their exact counterparts for the top 5 candidate routes and the resulting passenger coverage of the heuristic HHMoD connect approach may reach over 55%, which surpasses the coverage of the exact HHMoD only scenario. This indicates that the heuristic algorithms can also generate satisfactory candidate routes and the performances are appealing for large-scale real-world cases considering its much lower computation costs.  As for the number of required routes to cover all passenger demand, we observe that fewer than 30 routes may be sufficient with the HHMoD connect approach in all scenarios as shown in Figures ?? (c)-(d). This represents an additional 20 routes to serve the rest 30% to 40% of the passengers that are unsatisfied by the top 5 candidate routes. Without connecting to the subway system, the HHMoD only approach requires 3 to 8 more candidate routes to fulfill all passenger demand as compared to the HHMoD connect approach. And the heuristic methods would require an extra 2 to 10 routes depending on the particular times of the day. Finally, by combining the candidate routes to and from the activity hub following equation 12, we arrive at the completed round trip routes for serving the demand to and from the hub, and we visualize the outcomes of the top 5 round trip routes in Figures 6 and 7. In the figures, we distinguish the segments to and from the hubs with different lines. We also visualize the underlying layout of the NYC subway system and mark the stops that are served by connecting to the subway system in circles. As we can see in the figures, for maximum passenger coverage, the top routes generated for both HHMoD only and HHMoD connect share similar philosophy by prioritizing stops in Manhattan of high passenger demand but also visit several intermediate locations in Queens along the route which is permissible under the travel time constraint. And for all scenarios, the incoming and outgoing routes at the JFK airport are observed to be well-paired to form the round trip loops based on their trip demand and connection distance. From the results, we can also clearly tell the differences between the HHMoD only routes and HHMoD connect routes driven by the distinct philosophies of the underlying algorithms. In particular, the top 2-3 HHMoD connect routes are aligned with the NYC subway lines in Manhattan, and the passenger coverage is maximized by connecting to passengers along multiple subway lines such as lines 1-3 and lines A-F. Unlike the HHMoD only case, the HHMoD connect routes makes fewer stops in the upper west Manhattan areas with the connections to the subway system, which saves travel time from getting stuck in the heavy traffic. This allows for additional resources to operate a dedicated HHMoD route serving the Brooklyn area during AM and PM peak hours as the other routes have already achieved sufficient coverage in Manhattan. We also note that the 3rd route during night time connects to the G line at 21st Station in Queens, making it possible to serve transfer passengers whose destinations are in the west part of Brooklyn. Finally, for both HHMoD only and HHMoD connect, the direct connecting segment between the JFK airport and the LGA airport is observed to satisfy  Figure 6: Visualization of top 5 routes generated with connection to subway the passenger demand that travels between the two airport hubs during AM and PM peak hours.

Robust route scheduling
With the K-MCR finalized, we next discuss the results for the robust route scheduling problem under demand uncertainty, and here we only focus on the scheduling of the HHMoD connect routes due to its superior performances in the route generation step. As the input for the robust scheduling problems, the candidate HHMOD connect routes are generated with λ = 1.3 that cover all the candidate stops, and we consider that there are 200 available homogeneous vehicles with the capacity of 20 seats (such as medium-sized shuttle vans with rooms for luggage). We use Gurobi 9.0 package to solve the MILP for both master and recourse problems on a PC with 3.6 GHz CPU and 32GB RAM. The stopping criterion is set as the relative gap between the lower bound and the upper bound of the C&CG algorithm being smaller than 10 −4 : We first summarize the computational performances for the robust route scheduling problem from both the  Figure 8. We observe that the two-stage RO converges to the global optimum in 5-6 iterations in most cases. The convergence speed is faster for cases with very small Γ (close to 0) or large Γ (close to 200) values due to simpler problem structures that are similar to the nominal optimization problem or the single state optimization problem with box uncertainty constraints. In most of the scenarios, the robust route scheduling can be solved to global optimum within 200 seconds. The most computationally expensive scenarios are for Γ values in the range between 10 to 70, where we see a higher number of iterations and longer computational time, and a snapshot of the convergence can be found in Figure 8a for the off-peak period. We recall that each additional iteration for the C&CG algorithm will introduce a large number of variables and constraints into the master problem.
While the recourse problems are relatively easier to solve, the additional iterations and constraints lead to excessive computation time for the master problem, which accounts for almost 100% of total computation time in our experiments. On the other hand, even though certain Γ values are found to be associated with more iterations and high computational time, the C&CG converges quickly in the first 8 iterations (such as the cases with Γ = 10, 30 for PM peak scenario). Therefore, one may choose to perform early stopping and terminate with a minor optimality gap for reduced computational time with nearly identical performances under demand uncertainty. As for other Γ values, a small optimality gap can be achieved in only 3-5 iterations.
In this regard, although large-scale two-stage RO is in general difficult to solve, our numerical experiments suggest that the robust route scheduling problem can still be solved within a reasonable time frame. And this supports the applicability of the RO framework for the scheduling of urban HHMoD services over a short-time period.

Night
Operation cost Waiting and loss of service cost (b) Decompose objective function into F operation and F waiting + F loss Figure 9: Change of objective function with respect to different demand uncertainty Next, we discuss the impacts of the budget uncertainty on the performances of the scheduled HHMoD routes, and Figure 9a demonstrates how the objective function may change with increasing budget level. We observe that the objective function value resembles a concave and monotonically increasing function as the level of demand uncertainty increases. This observation is consistent across different times of the day. Based on the concave function shape, we report that majority of the impacts due to the demand uncertainty can be accounted for by modeling the demand uncertainty with a budget of 25% of the stops (Γ = 50) with a sharp increase in objective function values as compared to the nominal function value (Γ = 0). And the objective function values increase slowly for including the remaining stops in the budget of uncertainty. In addition, the changes in function values also reveal the differences in the robust route scheduling at a different time of the day. Specifically, the off-peak period is found to be associated with the highest cost, followed by the PM peak, AM peak, and night period. The reason is due to that the passenger demand level for PM peak and off-peak period at JFK is higher than that of AM peak and nigh time period. In addition, the off-peak period has the highest cost since the destinations, and origins of passengers to and from the hub are more diverse, and hence requires operating more vehicles (as shown in Figure 9b) to avoid the large penalty for the loss of services. And more conservative route schedules can be observed for the off-peak and night periods as the gap of function values between Γ = 200 and Γ = 0 is larger than the other two cases. And such conservative schedules are necessary as the demand distribution tends to have higher variations for the off-peak and night time periods.
In the end, we present the cost analysis to justify the values of the robust routing schedules by comparing the changes in objective function values under average scenarios and under heavy demand variations. We generate the same 100 random demand realizations for the average scenarios and calculate the expected cost for the robust route schedules with different Γ values. The cost-effectiveness in the average case can be measured as:Ḡ whereF Γ represents the average objective function values for the planned route schedules with budget Γ. We also measure the worst-case performances of the robust schedules by considering two demand variation scenarios: (1) the full variation case with all p i = 1 and (2) the half variation case with all p i = 0.5. And the corresponding cost-effectiveness under worst-case scenario can be computed as: Here,F Γ denotes the objective function value under demand variation. The value ofḠ Γ,0 indicates the level of conservatism of the robust solution and the value ofĜ Γ,0 measures the level of robustness of the route schedules as compared to the nominal results. Besides these three metrics, we also include the loss services rate as the percentage of unsatisfied demand over the total realized demand. We present the results on the gap under average and worst cases in Figure 10 and summarize the loss of services and the number of operated vehicles in each cases in Table 3.  (L represents the average case, L 0.5 and L 1 corresponds to the worst cases, and N denotes the number of vehicles) As we increases the budget of uncertain from 0 to 10, we observe from Table 3 that there is at least a 10% increase in the number of operated vehicles in all cases. And there may be 41 more vehicles in operation for Γ = 200 than that of Γ = 0 for the off-peak scenarios and the number is at least 26 more as in the AM scenario. Nevertheless, these significantly more number of vehicles in operation may only lead to an increase in average objective function value by 1-2% for the night time period or by up to 12.95% during the PM peak period. As a consequence, a majority of the increases in the operation cost are translated into the savings of loss of services and waiting time cost even during the normal demand scenarios. Specifically, we can verify from Table 3 that the additional operation costs contribute to the reduction of 0.55% 2.96% of unsatisfied passengers for normal conditions. And the benefits of the robust route scheduling is more evident if we look into the gains during the worst-case scenarios with both p = 0.5 and p = 1. In particular, Figure 10(d) suggests that a 2.7% increase in daily operation costs (with Γ = 60) may reduce the cost during half demand variations by 14.3% and by 28.5% with full demand variations. More importantly, this will also translate into 1.48% of unsatisfied passengers instead of 16.47% during the full demand variation cases. Therefore, these findings serve as strong support for the adoption of the robust scheduling for the planning of HHMoD services, especially considering its potential to avoid drastic loss of passengers, which is crucial to the long-term sustainable operation of HHMoD. Even for the most conservative case (Γ >= 20 for off-peak), we still observe a 6.1% reduction in objective function value with half demand variations and a 27.53% cost saving with the full demand variations during the off-peak hours. And this again can contribute to reducing the number of unsatisfied passengers from 15.29% to only 0.04% during the worst-case scenarios. One may also choose a less conservative strategy at Γ = 10 where significant savings in worst-case scenarios can still be achieved. Consequently, we report that robust route scheduling does not lead to overly conservative solutions during normal operations. At the same time, we can achieve significant savings in terms of operation cost during worst-case scenarios and reductions in the rate of unserved passengers in all scenarios. And these properties make it an ideal tool to support the scheduling of the HHMoD services over a short-time period by offering robust operation strategies against the demand uncertainties. As a final remark, the staircase patterns for the metrics in Figure 10 indicates that the performances of the robust route scheduling plateau after certain budget levels, which implies that the further increase in Γ will not lead additional cost even during normal operations. In this regard, one may directly solve the single-stage problem with box uncertainty constraints (Γ = 200) instead of the case with a medium level budget of uncertainty. This will obtain the same quality solution with minimum computational time, which could be a potential option if the robust route scheduling needs to be carried out in real-time.

Conclusion
In this study, we present the three-stage framework for optimal planning and scheduling of the HHMoD service at urban transportation hubs. The proposed framework consists of the route generation, the route combination, and the robust route schedule optimization. We develop efficient and effective route generation algorithms with connectivity to existing public transportation systems to generate the high-quality candidate route set. And a two-stage robust optimization model is developed to find the vehicle and route headway configurations that are resilient to the worst-case realization of passenger demand distributions. We conduct comprehensive numerical experiments for planning the HHMoD at the JFK airport in NYC. To best mimic the real-world settings, we use NYC taxi and FHV data to calibrate the potential passenger demand and use GoogleMap API to obtain the corresponding road traffic information. The results highlight the effectiveness of the proposed route generation algorithms for the HHMoD service, with 5 top candidates route being able to cover nearly 70% of total passenger demand and fewer than 30 routes are needed to meet all demand across the 100 stops. The results also demonstrate the cost-effectiveness of the robust route schedules during the average demand conditions and its superior performances in reducing the unsatisfied passengers during both normal and worst-case conditions. Moreover, we find that the quality of the generated routes and the performances of the robust route schedules are consistent across different times of the day.
There are several future directions to extend the scope of our study further. First, as our case study concerns the HHMoD services at the JFK airport, it will be interesting to examine the effectiveness of the solution framework at other activity hubs located in the central urban areas and test the performances in a different city. Second, while the proposed algorithms only consider the given travel time information, the exact and heuristic algorithms may be modified to incorporate the structural properties of the underlying service network (e.g. grid network or the planar graph) and design specific algorithms that can obtain optimal solutions more efficiently. Specifically, as the longest path problem is a special case of our problem, polynomial algorithms may exist on certain types of networks [35]. Finally, the robust scheduling problem can be further extended to account for the fleet with heterogeneous vehicles. This will contribute to the more efficient use of available resources and provide even less conservative operation strategies.

Initial route generation
The initial candidate routes are created by finding k-shortest routes from each stop i to hubs [36]. The process results in the candidate path set R h with k|N | candidate routes and we denote the stops that each r ∈ R h traverses through as N r . Since there may be multiple candidate routes that pass stop i, we then perform passenger demand assignment based on travel time: where t r i is the route travel time if passengers at stop i take route r and R h i is the set of heuristically generated routes that traverse from stop i. We can therefore measure the relative importance of each stop and each route as stop weight W i and route weight W r (1 r i denotes whether route r travels by stop i):

Candidate route expansion
The generated candidate routes in R h may be further combined to create new routes that extend the coverage of passenger demand. We consider expanding two routes r 1 , r 2 ∈ R h following two rules: 1) if r 1 is a sub-route of r 2 where N h r1 ⊆ N h r2 and 2) if r 1 and r 2 share common segments. The route expansion is then performed based on crossover operation over the actual link segments of the two routes following Algorithm 4. And the weight of the expanded route is also calculated by the stop weight W i .

Route pruning
The set of candidate routes may grow rapidly after the route expansion. To eliminate unnecessary routes, we conduct route pruning following two metrics as in [26]. The first metric is route circuity (total route length divided by the euclidean distance). We choose routes with length longer than l thd , and circuity smaller than C thd . The second metric is the subset relationship. If a route is a sub-route of the other route with a smaller length, then the shorter route will be discarded. The last metric is route similarity. It quantifies the spatial adjacency of two routes: S(r 1 , r 2 ) = l∈r1 min p∈r2 dist(l, p) + p∈r2 min l∈r1 dist(p, l) L r1 + L r2 (33) where L r is the length of route r and dist(p, l) is the euclidean distance between the center of two road segments p and l. To prune similar routes, we first sort the weight of the routes in R h in descending order: We then remove the following set of routes from R h based on route similarity threshold S thd : In this manner, it removes the repetitive routes with lower weight (importance).