1 Introduction
The rapid growth of MobilityonDemand (MoD) services such as Uber and Lyft is disrupting urban mobility. Over the last few years, this disruption has led to much interest in studying the design, management and impacts of these systems. Several analytical and simulationbased studies have focused on the efficient design of the MoD system – allowing passengers to share rides and managing the fleet to serve travel demand with the minimum fleet size, passenger’s waiting time, and vehicle miles traveled. In one of the first studies to develop a practical framework for managing a largescale MoD fleet, ma2013t
devised a heuristicbased taxi dispatching strategy and fare management system that could handle operations of large urban scaled fleets. To incorporate ridepooling (capacity 2),
santi2014quantifying developed the notion of a shareability graph, and showed that almost all the taxi demand in Manhattan, New York (around 150 million annual trips as of 2011) could be served by pairing up requests (2 requests per taxi). Pushing this result further, alonso2017demand presented a computationally efficient anytime optimal algorithm that can enable realtime high capacity ridepooling (up to 10 riders per vehicle) at the scale of Manhattan. Results show that 98% of the taxi demand in Manhattan can be served by 3000 vehicles of capacity 4 instead of the current fleet of 13586 active taxis.With recent advancements in automation technologies and government regulations enabling this technology, autonomous MoD (AMoD) systems^{3}^{3}3Note that the design framework of MoD and AMoD is similar if behavioral aspects and driver economics are omitted. are also gaining research interest (levin2016multiclass, spieser2016shared). For example, an agentbased simulation of AMoD services are conducted for both Austin, Texas (fagnant2015operations) and Melbourne, Australia (dia2017autonomous). Both studies concluded that AMoD services can satisfy travel demand with around onetenth of privatelyowned vehicles. A few simulationbased studies also subsequently quantified the impact of MoD and AMoD services on the environment (fagnant2014travel), congestion (fiedler2017impact), and parking (zhang2015exploring).
In areas where public transit serves a large proportion of demand, the interaction between MoD services and transit should not be overlooked. Considering MoD services as a complementary mode of transit, a few studies have designed a bimodal MoD system to solve the first and last mile problem (vakayil2017integrating, ma2017demand, shen2017embedding, moorthy2017shared). ma2017demand developed a realtime dispatching policy for MOD services in cooperation with existing mass transit service and tested it on a real network between Luxembourg City and its Frenchside crossborder area. shen2017embedding studied the integration of MoD service with Singapore transit using agentbased simulation. The authors concluded that the integrated system is financially viable and has potential to make transit attractive by reducing outofvehicle time.
However, all these studies assumed a fixed and exogenous demand for MoD services. Most models assume a waiting time threshold with the demand that cannot be picked up within that threshold being considered unsatisfied or dropout demand. The inherent assumption here is that passengers (from a certain demand set) that can be served by MoD will necessarily choose it. In practice, this may not be true: MoD services coexist with other travel modes (e.g., transit and walking). The demand for MoD is not only a function of its attributes (e.g., price, travel time), but also depends on characteristics of the competing travel modes as well based on individual preferences. Therefore, even if a passenger can be served by MoD services, they may choose other travel modes for a number of reasons.
Moreover, optimizing the operational performance of AMoD or MoD systems under an exogenous demand cannot answer policy questions related to the extended impact of these services on the urban transportation system as a whole. Such questions include the impact of MoD services on transit ridership (araldoimplementation), vehicle ownership (hao2017analysis), and demand for parking (zhang2017parking). Only a handful of studies have considered an endogenous demand model where MoD or AMoD systems compete with other travel modes based on service quality. horl2016simulation and araldoimplementation integrated such models into the MATSim and SimMobility platforms which are multimodal agentbased activity simulators. However, both simulators suffer from a common problem of relying on the scenariobased supplydemand analysis of virtual cities, not optimizing the supplyside parameters, and using a synthetic mode choice model.
A few recent studies have considered supplydemand interactions in a multimodal transportation system in the presence of ondemand mobility service. These studies have done so in the context of agentbased stochastic user equilibrium with a daytoday adjustment process (djavadian2017agenta). In particular, djavadian2017agentb models the dynamic adjustment process of an MoD operator. However, unlike our designfocused approach, these studies emphasize evaluating the sensitivity of different MoD operating policies as opposed to determining the optimal supplyside response.
The contribution of this study is twofold. The first contribution is to develop a unified framework which integrates mode choice models with a stateoftheart system for modeling realtime ondemand mobility services with varying passenger capacities (based on alonso2017demand). We illustrate the proposed framework by developing a multiscale MoD fleet management system to serve the taxi demand in Manhattan where passengers are utility maximizers and choose a travel mode from a set of four mutually exclusive alternatives: ridehailing service (capacity 1), ridepooling service (capacity 4), microtransit service (capacity 10) and public transit (e.g. subway). The first three travel modes are MoD services assumed to be run by a single MoD provider. Note that the underlying mode choice model is calibrated with stated preference survey data from New York City.
The second contribution of this study is to implement and integrate a Bayesian Optimization (BO) based solver to find the optimal supplyside MoD parameters (e.g., fleet size, fare, etc). Historically, two types of methods have been used to determine the supplyside parameters of MoD systems: i) Heuristic methods to simultaneously optimize supplyside parameters and the vehicle routing problem (VRP) (liu2009effective, repoussis2010solving). Due to the complexity of VRP, these methods are inefficient for optimizing largescale MoD operations; ii) Predefine candidate combinations of parameters (or scenarios), and then conduct a grid search to find the one that gives the optimal objective function value in the simulation (li2010optimizing, brownell2014driverless, marczuk2015autonomous, azevedo2016microsimulation, boesch2016autonomous). The performance of this approach highly depends on the quality of the predefined set of candidate solutions. In addition, grid search struggles in finding even nearoptimal solutions in a high dimensional space in reasonably short time. In this article, we decouple the optimization of supplyside parameters and the MoD vehicle assignment problem. We consider the transportation simulation system as a blackbox function, and use BO as a sequential search strategy to optimize the supplyside parameters.
In summary, the proposed framework would not only help in better understanding the influence of MoD services on urban mobility and answer systemlevel policy questions, but would also illustrate a method to tune the supplyside parameters that influence the demand.
We note that our current study considers the case of a single MoD operator and does not address competition between providers, which is a limitation of the current model. Realistically, several MoD and taxi services coexist and their behaviors are dependent on competitive factors. Competition between these services and the resulting market equilibrium (wang2016pricing, qian2017taxi, heilker2018duopoly) and their impact on operator efficiency (sejourne2018price) have been studied independently. Integrating these competitive models in the proposed framework is a direction for future research.
The rest of this article is structured as follows. In Section 2, we formulate the mathematical problem and propose the unified framework. In particular, we describe two essential components: the mode choice model and the MoD simulation system. In Section 3, we introduce a BObased approach to find the optimal supplyside parameters. Section 4 consists of numerical experiments on the Manhattan network, which are focused on the numerical convergence of the proposed framework, efficiency of the BO approach introduced in this work, and practical insights that the framework provides. Finally, Section 5 provides closing remarks and discusses possible directions for future research.
2 Framework: Modeling MobilityonDemand Systems within a Multimodal Transportation System
We consider a general transportation system that includes public transit and three classes of MoD services with varying passenger capacities (i.e., maximum occupancy). The public transit option is considered to be a complimentary travel mode to the MoD services, but note that the framework is not limited to this specific case^{4}^{4}4In the framework, the interplay between all travel modes is considered in the mode choice model. After the demand for each travel mode is generated, the simulator for each mode is run independently. Therefore, we can use any set of candidate travel modes for any purposes. . The goal is to develop a simulation optimization framework to optimize both i) the systemlevel objective functions (e.g., MoD system operator’s profit or consumer surplus) when the demand for MoD services is not exogenous and rather depends on their and other competing travel modes’ level of service (e.g., waiting time, travel time, and travel cost), and ii) the transportation system operations given a certain set of supply and demand side parameters. The proposed simulation and optimization framework is shown in the Figure 1.
The proposed framework can be disentangled into inner and outer loops. The outer loop iteratively optimizes the supplyside parameters using BO and terminates when a stopping criterion of BO is satisfied. For a given set of supplyside parameters, the objective function (e.g., operator’s profit) of the outer loop is evaluated at the equilibrium of the inner loop. The inner loop iteratively performs the following steps: i) evaluate mode choice probabilities, ii) simulate the demand for and operations of public transit and MoD services, iii) update modespecific attributes (e.g., waiting time) as the inputs to the choice model in the next iteration, and iv) repeat until an equilibrium is reached, i.e. when the average difference in mode shares of two consecutive iterations is lower than a predefined threshold (see Section
4.1). In other words, if each iteration of the inner loop is interpreted as a “day”, for a given set of supplyside parameters, passengers learn from the experienced historical levelofservice of travel modes so that they can make more informed mode choices on a given day. The inner loop terminates when passengers saturate their learning and make nearly consistent mode choices on consecutive days.This section discusses the different components of the inner loop such as the mode choice model, and operations of public transit and MoD services. Section 3 discusses a BObased approach to find the optimal supplyside parameters (outer loop).
2.1 Mode Choice Model and Data
To understand preferences of New Yorkers for MoD services, we conducted a statedpreference (SP) study. In an online survey, the respondents were asked about sociodemographics, travel characteristics, and various other opinions. The survey also had a discrete choice experiment (DCE) in which each respondent was asked to choose the best and the worst travel mode from a set of three choices: Uber (without ridesharing), UberPool^{5}^{5}5UberPool represents MoD services with a passenger capacity of more than one. (with ridesharing), and their current travel mode (the one used most often on their most frequent trips). Table 1 shows the attribute levels of the DCE design. In the case of monthly payment of trip and parking costs, a per trip cost was computed by dividing monthly cost with trip frequency. The per mile cost of a private car ($.45) was obtained by summing the insurance, maintenance, and fuel cost (AAA News Room, 2016). In the DCE design, we made sure to constrain invehicle travel time (IVTT) and per mile cost of Uber to be less and more, respectively, than those of UberPool. Respondents were provided textual and visual information about all alternatives at the starting of DCE.
To obtain priors for pivotefficient DCE designs^{6}^{6}6In pivotefficient designs, attribute levels shown to the respondents are pivoted from reference alternatives for each respondent. In this study, the travel mode used on the most frequent trips was considered as the reference alternative., we first conducted a pretest. A Defficient design with zero priors containing 4 blocks (6 choice situations per block) was generated with the Ngene software (choicemetrics2014ngene) for the pilot study. The attribute levels of Table 1 were used in the design. The online pilot survey was created using the webbased Qualtrics platform, but the survey was distributed among a continuous panel provided by Survey Sampling International (SSI, a professional survey firm) in February 2017. Those who drive for any MoD service or are younger than 18 years were considered ineligible to participate in the survey. Those who completed the survey in less than 10 minutes or provided conflicting responses (e.g., reported more children than household size) were discarded. After eliminating such responses, 298 (out of 397) completed responses were used as valid pretest observations for further discrete choice analysis.
We used prior parameter estimates from the pilot study to create a pivotefficient design with 6 blocks (7 choice situations per block). All the attributes and attribute levels remained the same (Table 1). Table 2 shows an example of the choice situation of the final mode choice experiment.
We conducted the main study during OctoberNovember 2017. After data validation tests, preferences of 1507 (out of 1689) respondents were used in the model estimation. We estimated a multinomial logit (MNL) model and used MNL attribute valuation for prediction of mode choice probabilities, which are needed in the simulation. The closedform choice probability MNL expressions allow a seamless integration in the MoD simulation framework
^{7}^{7}7We tried nested logit specifications with different nesting structures, but in all scenarios, one of the nests’ elasticity estimates was above one. In other words, nest correlation was out of the accepted range for compatibility with random utility maximization. This implies preferences of the sample are not aligned with nested logit correlation.. Table 3 summarizes MNL parameter estimates. By dividing marginal utilities of attributes with that of trip cost, willingnesstopay (WTP) estimates are derived. By using these marginal rates of substitution, the model provides evidence that New Yorkers are willing to pay $25.9 and $18.6 to save an hour of OVTT and IVTT, respectively. The recommended hourly value of travel time for local commute by passenger car in downstate New York is $15.6 (Department of Transportation, New York State, 2012), which is close to our estimate of WTP to save an hour of IVTT. Another study of economic evaluation (US Department of Transportation, 2011) estimates that walking, waiting, and transfer time in personal travel should be valued at $19.10  $28.70 per hour. Our estimates of WTP to save an hour of OVTT falls in this range. Although WTP estimates can be transferred directly to the MoD simulation, alternative specific constants (ASCs) require recalibration (ASCs in SP studies are just manifestation of sample shares of alternatives). This process is consistent with the estimated mode choice model, but is not ideal because choice sets of SP study and the considered taxi demand are different^{8}^{8}8The recalibration of ASCs is much easier for the full travel demand model where all possible alternatives are available to travelers and their real mode shares are known.. Details about recalibration of ASCs can be found in Sections 4.2 and 4.4.2.2 Simulation of MoD and Transit Operations
In this section, we introduce the system we use to simulate the operations of the MoD and transit services.
2.2.1 Simulation of the MoD system
We consider a fleet of vehicles with varying passenger capacities to satisfy the demand for MoD services. The operation of MoD services includes two main tasks: i) match travel requests to vehicles; ii) rebalance the idle vehicles to areas with potential future demand. Our formulation follows a state of the art framework which was recently proposed by alonso2017demand.
The set of vehicles is denoted by . Each vehicle has a corresponding capacity . The set of trip requests is denoted by . Multiple travelers are allowed to share one single ride. A passenger is defined as a past request that has been picked up by a vehicle and is now enroute to its destination. For each request , the waiting time is , which is the difference between the pickup time and the request time . For each picked up request , the total travel delay is defined as , where is the dropoff time and is the earliest possible time at which the destination could be reached. is computed by following the shortest path between the origin and the destination with zero waiting time. We want to find the optimal assignment that minimizes a cost function under a set of constraints as follows:

For each request , the waiting time must be below the maximum waiting time . If one customer has waited for a time more than , the request is discarded in the simulation, and a large constant penalty will occur in the cost function.

For each picked up request , the total delay must be lower than the maximum travel delay .

For each vehicle at any time, the number of passengers in the vehicle must not be larger than the capacity

The requests for travel via a specific MoD service type (e.g. ridehailing service, ridepooling service or microtransit service) can only be served by the vehicle fleet corresponding to the same service type.
For convenience, the cost function is set to the sum of total delay and the penalty for unassigned requests in this work, but note that the objective function can be tuned for any other purposes. The total delay includes both the waiting time for all assigned requests and the invehicle delay caused by sharing with other passengers.
The MoD operation problem is essentially the same as the dynamic pickup and delivery problem (savelsbergh1995general) with time windows, which has shown to be difficult to solve exactly for a largescale demand (nanry2000solving, ropke2006adaptive). The anytime optimal algorithm in alonso2017demand decouples the problem by first computing feasible trips from a pairwise shareability graph (santi2014quantifying)
and then finding the optimal trip assignment to vehicles by solving an Integer Linear Program (ILP) of reduced dimensionality. In our study, the assignment of the trip requests is processed at a frequency of every 60 seconds. Requests that have not been picked up in an assignment round remain in the request pool until the maximum waiting time constraint is violated.
Specifically, each round of the assignment simulation includes the following steps:
(1) Compute a pairwise requestvehicle graph (RVgraph). RVgraph describes possible pairwise matching between vehicles and pickup requests. In the graph, we connect two requests and if an empty virtual vehicle at the origin of one of them can serve both of them without violating any of the constraints; Similarly, we connect a request to a vehicle if can serve under the constraints.
(2) Construct the RequestTripVehicle graph (RTVgraph) using the cliques of the RVgraph, which includes all the feasible trips and also the vehicles that can serve them. A feasible trip is defined as a set of requests that can be served by one vehicle under the constraints. In the RTVgraph, a request is connected to a trip if contains ; A trip is connected to a vehicle only if can complete without violating any constraints. The feasible trips and the edges in the graph are computed incrementally in trip length for each vehicle, which successfully reduces the search space (alonso2017demand). This step results in all the feasible assignments between the vehicles and the requests.
(3) Solve an ILP to compute the optimal assignment of vehicles to trips. Since each trip can possibly be served by multiple vehicles, the ILP needs to not only determine the assignment, but also ensure that each request and vehicle are assigned to only one trip. The ILP in alonso2017demand is as follows:
(1) 
where is the set of all feasible assignments between trips and vehicles, is the cost of vehicle serving trip , is the constant penalty for an unassigned request, which controls the tradeoff between the total delay term and the penalty term. In our experiments, is set to a large constant to make sure that the objective is to minimize the total delay given that the maximum possible number of requests are served. If is set to a relatively small value, some requests may be rejected, since the cost of serving some requests may be higher than the penalty for not serving them. For example, if the pickup point is far away from all available vehicles. The decision variables are and , where if vehicle is assigned to trip , if the request is ignored in the assignment. In the constraints, there are three sets: the set of trips that can be served by vehicle is ; the set of trips that contains request is ; the set of vehicles that can serve trip is . The constraints ensure that: i) each vehicle will only serve one trip; ii) each request is either served by one vehicle or ignored () in the assignment. In this article, we use the following formulation which reduces the number of variables by compared to the above formulation.
(2) 
where is the number of requests in trip . In this formulation, we track for each trip . A reward (negative penalty ) of serving all the requests in trip is added to the objective function. In formulation 1, a penalty is imposed on the objective function for each unassigned request. In formulation 2, we consider all the requests unassigned before the assignment, and thus a total penalty ^{9}^{9}9This term is ignored in the objective function since adding constants does not affect the optimal solution. is added to the objective function. Whenever we assign a vehicle to a trip , the system gets a reward for picking up requests in trip . Therefore, the solution to formulation 2 is equivalent to the solutions to formulation 1 since the objective function is still the sum of the cost for each assigned trip plus the penalty of all unassigned requests, but the number of variables is reduced by because we do not need to include for each request in the formulation.
(4) Rebalance the idle vehicles. After the assignment, there might still be idle vehicles and unsatisfied requests. It is reasonable to believe that more requests will be likely to appear in the neighborhood of the unsatisfied requests. Therefore, we solve a linear program as in alonso2017demand to match the idle vehicles to the location of unsatisfied requests while minimizing the rebalance cost.
2.2.2 Simulation of the Public Transit Service
We assume that the network and operational characteristics of the existing public transit system (headway, travel time, etc) are known. To obtain the level of service attributes of the public transit in the mode choice model, we perform a transit network assignment for each origindestination (OD). We first combine the road network we use in MoD service simulator and the public transit network. Then we use an allpair shortest path algorithm to find the walktransitwalk path for each OD pair. This procedure consists of the following steps:
(1) Construct the public transit network. The public transit data is obtained from the General Transit Feed Specification (GTFS) dataset published by the transit authority through Google’s GTFS project. The public transit network topology is created according to the schedules and associated geographic information in the GTFS data. The weight on each link is the scheduled travel time between two transit stations.
(2) Combine the road network and the public transit network. The weight on the link between two road network nodes is set to the walking time (). For each road network node, we add a link to connect it to each public transit network node within the walking range (e.g. 0.5 miles). As the origins and destinations are part of the nodes in the road network, the weight on the link connecting the road network node and the public transit node is set to a generalized cost, which is the sum of the walking time (), the expected waiting time (half of the headway for the public transit line) and the trip fare (converted to seconds) for the public transit line.
(3) Compute the allpairs shortest path using the combined network, and store them in a lookup table. During the simulation, we refer to the table to find the path and modespecific attributes for public transit in the mode choice model.
2.3 Updating Historical Trip Attributes and Inner Loop Stopping Criterion
To model the memory and learning process of the travelers, we store and update the historical trip attributes for each OD after each iteration of the inner loop, and use these values as inputs to the mode choice model at the next iteration. Since the demand is sparsely distributed on the thousands of network nodes (e.g., 4092 nodes in our network), we use Kmeans clustering algorithm
(macqueen1967some) to spatially cluster the nodes. If the walking range for the passenger is miles, we determine the number of clusters by satisfying . The historical trip attributes are stored and updated at the clusterlevel.Due to the maximum waiting time and the maximum delay constraints, some of the demand for a certain MoD service may not be satisfied by the system. We assume that the demand that is unsatisfied (due to these constraints) will switch modes and take public transit. Therefore, for each cluster pair , we also store the service rate for each MoD mode , which is the proportion of travel demand served by mode for cluster pair .
We compute the utility of the MoD services for the next iteration as follows (accounting for unsatisfied requests): Assume that for the cluster pair , the utility of taking MoD mode and being picked up in the simulation is , and the utility of taking transit is . Then, the utility of taking the MoD mode in the next iteration is computed according to the following equation:
(3) 
A constant penalty multiplier is used since transit was not the passenger’s first choice in the last iteration. A high penalizes the utility function more if a request is not satisfied, and further decreases the corresponding MoD mode’s share. In the numerical experiments, we use . This parameter should ideally be calibrated using a stated preference study, which we leave for the future research. and can be computed according to the historical trip attributes. After each iteration, the historical trip attributes are updated according to the following equation:
(4) 
where is the set of historical trip attributes for cluster pair for mode at iteration which consists of the invehicle travel time, waiting time and the service rate, ^{10}^{10}10In our numerical experiments, we let . is a constant coefficient that controls the balances between the historical information and the new information, is the trip attributes for cluster pair for mode which we obtain from the simulation in the current iteration by averaging the attributes for all ODs in the cluster pair.
The inner loop simulation procedure is iterated until the following stopping criterion is satisfied:
(5) 
where represents the average mode share difference among all the travel modes between iteration and , is the set of candidate travel modes, is the share of mode at iteration . If attains a value smaller than a predefined threshold of 0.01, the inner loop is terminated. In other words, the learning of travel demand saturates and system appears to achieve an equilibrium as a result of their consistent travel mode choices on consecutive ”days”.
3 System Optimization using Bayesian Optimization
In Section 2, we discussed a simulation framework where the inner loop considers the supplydemand interaction for a given set of operational parameters of MoD services. In this section, we study the problem of optimizing these parameters which include fleet size and the pricing rules for MoD services with varying passenger capacities.
The innerloop simulation can be considered as a blackbox function for which the set of supplyside parameters (i.e., decision variables) should be chosen such that the system performance metric at the equilibrium of inner loop attains its optimal value. The performance metric depends on the influence of different stakeholders. For example, whereas policymakers might be interested in optimizing the consumer surplus and vehicle miles traveled, MoD operator might be interested in just maximizing the profit.
This optimization problem is difficult to solve because of the following properties:

The objective function cannot be evaluated analytically and its derivatives are thus not easily available;

One can obtain observations of the objective function, but the function evaluation for even one set of decision variables is computationally expensive;

The function evaluation can be affected by simulation noise such as prediction of mode choice probability and other random factors (e.g. the initial location of the vehicles) in the MoD service simulator.
Whereas these characteristics make the problem intractable to solve using analytical optimization methods, Bayesian Optimization (BO) is an appropriate and powerful tool in such settings. BO is a sequential search strategy for the global optimization of an expensive blackbox function (mockus2012bayesian)
. It first emerged as a successful strategy in many machine learning applications
(bergstra2011algorithms, snoek2012practical, bergstra2013making, swersky2013multi, swersky2014raiders, yogatama2015bayesian), and has lately been employed in many other areas including robotics (lizotte2007automatic, calandra2016bayesian), sensor networks (garnett2010bayesian), environmental monitoring (marchant2012bayesian), information extraction and retrieval (wang2014bayesian, li2018bayesian), and game theory
(picheny2016bayesian).Since the objective function does not have a closedform expression in terms of the decision variables, BO treats it as a blackbox function with some prior belief. Here prior represents our belief about the space of possible objective functions. As we obtain more observations of the function, our belief about the objective function is updated by combining the prior and the likelihood of the data already acquired to get a potentially more informative posterior. The posterior distribution is then used to select the next set of decision variables for the function evaluation.
Specifically, a BO framework has two key ingredients: i) a probabilistic surrogate model for the expensive objective function ; ii) an acquisition function that is maximized to determines the next point for function evaluation (shahriari2016taking). The BO framework is demonstrated in Algorithm 1, where is the surrogate model, is the acquisition function, and , which is the set of historical function observations until iteration . Each observation is denoted as where is the set of inputs, and is the objective function value we obtain by evaluating the objective function .
We now provide a brief introduction of these two ingredients. Readers can refer to brochu2010tutorial and shahriari2016taking for more details of the method.
3.1 Surrogate Model
Since evaluating is expensive, the BO framework builds a surrogate model using the historical observations to approximate the objective function. In this work, we adopt the Gaussian Process (GP) as the surrogate model, which is the most commonlyused model in BO (mockus1994application)
. GP is a stochastic process which can be seen as a collection of random variables where any finite set of random variables follows a multivariate Gaussian distribution. Assume that the set of points we use to build the GP is denoted by
, then the GP is a fully specified by its mean function and covariance function . Intuitively, GP is a distribution over functions. Given a point, it will return the mean and the variance of a normally distributed variable
which is a prediction for .For simplicity, the mean function is usually defined to be zero function, i.e., (brochu2010tutorial, marchant2012bayesian, wang2014bayesian). Assume that we have historical observations , and we want to use GP to predict on an arbitrary point . GP will return a normally distributed variable with both the mean and variance. Let and
. As the definition of GP reveals, the joint distribution of the historical observations and
is as follows:(6) 
where , and is the covariance matrix with each entry .
Then, we can obtain the predictive distribution over as follows:
(7) 
where , .
The choice of covariance function for GP determines the smoothness properties of samples drawn from it (brochu2010tutorial). Similar to many other studies, we use the Matern kernel (matern2013spatial), which incorporates a smoothness parameter to provide flexibility in modeling functions:
(8) 
where and are the Gamma function and the Bessel function of order , respectively. When , the kernel reduces to the unsquared exponential kernel, and when , the kernel reduces to the squared exponential kernel. To provide appropriate smoothness, BO applications usually use and (rasmussen2006gaussian). In our experiments, we use . Readers can refer to hutter2013kernel and swersky2014raiders for more details about the GP as well as its covariance functions.
3.1.1 Acquisition Function
The acquisition function is used to select the next realization of the decision variable for the evaluation of the objective function. It usually considers both the mean and variance of the predictions provided by the surrogate model. Intuitively, the process can be seen as maximizing the utility of the next sampling.
A good acquisition function is the one that finds an elegant tradeoff between exploration and exploitation^{11}^{11}11Exploration represents searching in regions with high uncertainty but might have observations with high objective function values. Exploration helps in avoiding being trapped in a local optimum. Exploitation represents searching in the regions with high expected values of the objective function given the information provided by historical function evaluations. while searching for the optimum of the objective function. The BO literature proposes many acquisition functions such as the probability of improvement (PI) (kushner1964new), the expected improvement (EI) (szego1978towards), and the upper confidence bound (UCB) (srinivas2009gaussian). Similar to many other studies, we use UCB as the acquisition function.
(9) 
where
is a hyperparameter which establishes a balance between the exploration and the exploitation. In our experiments, we employ a special case of UCB which is called GPUCB
(srinivas2009gaussian, calandra2016bayesian) that casts the BO problem as a multiarmed bandit problem. In GPUCB, is automatilly updated according to the following euqation:(10) 
where is the number of past objective function evaluations, is a parameter of choice, is the dimensionality of the search space. We use the same as in srinivas2009gaussian.
3.2 Objective Function
In this article, we define the objective from the perspective of MoD service provider. Consider that all MoD operations are run by one service provider who is a profit maximizer. The profit is defined as the difference between the revenue and the operating cost. We use the ridehailing service (capacity 1) as the base mode for the fare calculation of other MoD services and compute its trip fare according to the equation for UberX service (Uber). The fare for all other MoD modes with higher passenger capacities are computed by applying a discount factor () on the fare of the ridehailing service. For a passenger that is served by the tier MoD service, the trip fare is computed as follows:
(11) 
where and are the minimum fare and base fare for the ridehailing service, and are the travel time and travel distance for , is the time rate for the ridehailing service, which is the fare per second, and is the distance rate for the ridehailing service, which is the fare per mile.
The objective function of the MoD service provider is computed as follows:
(12) 
where is the set of all passengers, is the leasing cost for a tier MoD vehicle in the experiment time span, is the salary that the operator pays to tier driver^{12}^{12}12The drive salary can be ignored in the case of AMoD systems., is the total distance traveled by the tier MoD service and is obtained from the simulator, and and are the fleet size and operating cost ($ per mile) of the tier MoD service. The set of discount factor and fleet size for each tier of MoD service are denoted as and , which are the decision variables.
Even though this work considers the objective from the perspective of a service provider, the proposed framework is more general and can be applied to designing policies, subsidies and regulatory strategies, guiding infrastructure planning and deployment, and operational management of transit services in the presence of MoD systems.
4 Numerical Experiments
In this section, we present a series of numerical experiments to show: i) the numerical convergence to an equilibrium of the mode choices within the multimodal system; ii) the calibration of the ASCs; iii) the performance of the BObased optimization algorithm; iv) relevant applications of the proposed framework. In the following experiments, convergence to an equilibrium represents that the system reaches a state where the average difference in mode shares of two consecutive iterations is lower than 1%. The system is implemented using Python , and all experiments are conducted on an Intel core I7 computer (3.4 gigahertz, 16 gigabytes RAM).
As a proxy for the real OD demands, we use the publicly available dataset of taxi trips in Manhattan, New York (donovan2015using). In our system, we serve this demand via either i) the ridehailing service (capacity 1); ii) the ridepooling service (capacity 4); iii) the microtransit service (capacity 10), and iv) public transit (e.g. subway). The first three modes are MoD services offered by one MoD service provider. We consider the pickup time for the taxi trips in the dataset to be their departure time . The network we use is the entire road network of Manhattan (4092 nodes and 9453 edges) (santi2014quantifying, alonso2017demand). The link travel time is given by the daily mean travel time, which is computed using the method in santi2014quantifying. For public transit, we use the subway network provided by the Metropolitan Transportation Authority (MTA). The following experiments are conducted with respect to the rush hour demand (8 am to 9 am) on an arbitrary day (Monday the 6th of May, 2013).
The system has 5 decision variables: the fleet size of the ridehailing service () and the fleet sizes and discount factors for the ridepooling and microtransit services (, , , and ). Note that at the initialization stage, there is no information about the ODspecific MoD attributes (e.g., waiting time and travel time). Therefore, we initialize the attributes as follows: Let be the maximum waiting time allowed in the MoD system and be the shortest path travel time for each . The travel time of the ridehailing service, ridepooling service, and microtransit service for is set to , , respectively and the waiting time is set to , , respectively. In our experiments, we set at 10 minutes.
To account for the labor cost of the drivers, in this study, we consider a driver salary of $17 per hour for all MoD services, which is close to the hourly mean wage of taxi drivers in New York state (United States Department of Labor, May 2017). The leasing cost includes insurance, depreciation, maintenance, and other registration charges, and the operating cost includes fuel and tier cost. We compute these costs for different MoD services using statistics provided by AAA News Room, 2016. The leasing cost of $11.97 per day per vehicle is obtained for the ridehailing and ridepooling fleet of Sedan cars and $19.32 per day per vehicle for the microtransit fleet of Minivans. We normalize these numbers based on the proportion of daily demand (5.94%) served during the experiment hour. The operating cost turns out to be $0.1473 per mile for all MoD services.
4.1 Convergence to an Equilibrium
In this section, we illustrate the proposed multimodal system converging to an equilibrium with respect to mode choices using two test cases. In the test cases, we are not optimizing the system, but rather using a fixed set of supplyside parameters to test for the numerical convergence of the inner loop simulation to an equilibrium. The experiments are run for 20 iterations for each set of parameters.
The supplyside parameters for the first case are as follows: , , , and . The iterative variation in the share of each travel mode is shown in Figure 3. The share of each mode is relatively stable after around 5 iterations, which shows the numerical convergence of the system to an equilibrium. Note that the mode share is likely to fluctuate a little even after reaching equilibrium. This is just a manifestation of the simulation noise of the probabilistic mode choice model. In this test case, the mode share at equilibrium is not significantly different from the mode share at the beginning because the initial waiting time and travel time are close to the ones that the MoD services provider can offer given the supplyside parameters.
We set , , , and in the second test case. We intentionally consider the lower fleet size of the ridehailing service (capacity 1) to make its level of service poorer than the one achieved with the initialized waiting time and travel time. The share of each mode is shown in Figure 3. As expected, the mode share of the ridehailing service quickly decreases in the first 5 iterations. According to initial conditions, the ridehailing service is attractive but its demand rapidly decreases as passengers learn that the service is not able to serve the demand. In addition, the mode share of the ridepooling service increases rapidly in the first 5 iterations as it offers a higher level of service than initial conditions due to a high fleet size of 1000 vehicles. The stop criterion value for the last five iterations are reported. For test case 1, the stop criterion value for the last 5 iterations are 0.0031, 0.0029, 0.0018, 0.0027 and 0.0027. For test case 2, the stop criterion value for the last 5 iterations are 0.0045, 0.0025, 0.0027, 0.0028, and 0.0018.
Note that even with a fleet size of 100 vehicles in the second test case, the ridehailing service could attain the larger share than both high capacity MoD services. There are two reasons for this to happen: i) The demand data of Manhattan contains many short trips. The passenger is likely to choose the ridehailing service over high capacity services for short trips because even the discount cannot compensate for the lower level of service in high capacity MoD; ii) As discussed earlier in Section 2, we assume that unsatisfied MoD demand (when the waiting time is larger than the threshold) shifts to the public transit and the penalty on the utility of MoD service is applied in the next iteration. It appears that even though the service rate for ridehailing service is low, the penalty is not high enough to shift the demand to other high capacity MoD services.
Additionally, the public transit share is very high (35%  40%) in both cases, even though the input data corresponds to trips that were taken by taxi in reality. This is a manifestation of not having an alternative specific constant (ASC) for transit in the mode choice models. Therefore, the ASC of the public transit needs to be calibrated to lower its mode share for this specific travel demand considered in this study, which we describe next.
4.2 Calibrating the Alternative Specific Constant of Transit.
In the absence of revealed preference data, calibration of ASCs is difficult. Since ridehailing service share characteristics of Uber, and ridepooling and microtransit services are similar to UberPool, we use their respective ASCs as our best guess in the case study. Moreover, MNL estimates of marginal utilities of alternativespecific attributes are used in the case study to ensure consistency of willingness to pay and scale.
If we had the true share of transit for the given demand, then the ASC of transit could have been directly calibrated using the method suggested by train2009discrete. Since, the travel demand considered in the study corresponds to trips that actually used the taxi as travel mode, a transit mode share of 010% in our results was endogenously set as an appropriate target. To calibrate the transit ASC of the choice model for this specific population, we use a grid search approach. Since the demand we consider was served by 13,586 active taxis, we ran the inner loop simulation at different values of transit ASC with only transit and a fleet of 13,586 ridehailing vehicles. The share for public transit for various ASC values are shown in Table 4. Transit attains a share of 6.3% at an ASC value of 3, which appears appropriate for the considered travel demand and thus is used in the remaining experiments. Figure 4 shows the results of running the test case 2 in Section 4.1 with 3 as the ASC for public transit. The shares of the ridehailing and the ridepooling services become and which were and when ASCs were not considered (see Section 4.1).
ASC  Mode share (%) 

1.0  32.9 
1.5  23.2 
2.0  15.5 
2.5  10.2 
3.0  6.3 
4.3 Algorithm Performance.
In this section, we illustrate the application of the proposed BO framework to optimize the supplyside parameters and evaluate its performance by: (1) Use a smaller test example (10% of the real travel demand) to compare BO solution with the bruteforce solution; (2) Compare the performance of BO and random search method (bergstra2012random) using realscale travel demand.
Optimizing the supplyside parameters for an MoD system is analytically intractable due to the complexity of the supplydemand interactions, since the fleet assignment problem is a function of the demand and vice versa. Therefore, it is difficult to find an optimal solution even for relatively small test cases. Instead, we use the following approach to approximate the optimal solution via a grid search of relatively high resolution and compare that solution to that of BO. First we decrease the travel demand to 10%, and change the maximum value for , and to be 800, 300 and 150 respectively. The discount factors and are set to constants (0.1 and 0.2, respectively) to further decrease the search space. In this setting, we enumerate every possible combination of the fleet size in and find the solution with the highest profit as the nearoptimal solution. On the same test case, we employ BO with three different acquisition functions (UCB, PI and EI) to understand the sensitivity of BO’s performance relative to the acquisition functions by comparing the gap between the respective BO solutions and the nearoptimal solution found by the full enumeration method.
In the full enumeration method, we evaluate the profit for 2912 parameter settings. The best solution is obtained with , , , and provides a profit of 13001. Figure 5 shows the variation in the highest profit with the number of function evaluations for BO using three different acquisition functions. It can be seen that PI performs worse than EI and UCB, while EI and UCB perform similarly well in our example. The best solution found by BO is , , , which gives a profit of 12674. Although the optimal supplyside parameters provided by BO are quite different from the nearoptimal solution of the full enumeration method, the profit gap is only about 2.5%.
In the second experiment, we use realscale travel demand and compare the performance of BO and random search method. In the random search method, we sample the decision variables from a predefined distribution (uniform distribution in our case), evaluate the objective function at those realizations, and report the realization of the decision variable corresponding to optimal objective function value as the optimal solution. For a fair comparison, we keep the same number of function evaluations for BO and the random search method.
Figure 6 shows the variation in the highest profit with the number of function evaluations in both optimization methods. The highest profit bound by BO and the random search are 145,015 and 126,308 respectively, which represents that BO’s solution increases the profit by 15% compared to the random search’s solution. This considerable gap between objective functions shows the advantage of using BO over the random search. We expect this gap to be larger in problems with more decision variables i.e., higher dimensional solution spaces.
The BO optimal solution is: , , , and . It is important to note that the discount factor is 0 for both high capacity MoD services, which seems quite unrealistic because passengers are not likely to use high capacity MoD services if they have to pay the same price as of ridehailing service. This experimental result implies that discounting the price of high capacity MoD services is not able to attract a level of demand that is high enough to compensate the reduction in profit due to the discounting. In general, the passengers have different perceptions about the discounted fare of the high capacity MoD services, but such latent factors are hard to quantify. The current mode choice model we use assumes that the tradeoff between sharing a ride and the cost savings of doing so can be modeled as a linear relationship, and does not consider the fact that travelers perceptions of this tradeoff may be more nuanced.
4.4 Scenario Analysis.
Modeling passenger’s perception of discounting is left for future research but we illustrate implications of including such latent factors in choice models using a scenariobased analysis. We hypothesize that a passenger is less likely to choose high capacity MoD services at low values of the discount factor, even if these services have the same attributegovernedutility^{13}^{13}13The attributegovernedutility implies the part of utility which depends on observed attributes of the travel mode such as travel time, waiting time and travel cost in our case. as of the ridehailing service. We represent this hypothesis by adding a function of the discount factor, which is denoted as , in the utility equation of the high capacity MoD services:
(13) 
We call the function of the discount factor that we add to the utility equation discount factor function. We optimize MoD operations under three different discount factor functions, which are shown in Table 5. The coefficients in are set such that three scenarios represent low disutility, medium disutility and high disutility for having a low discount factor of high capacity MoD services. Since the passenger is likely to expect a higher discount for the microtransit service than the ridepooling service, we use different values of coefficients in both functions to mimic these expectations. The plots of discount factor functions for the ridepooling and microtransit services under three scenarios are shown in Figures 8 and 8.
Scenario 




1  
2  
3 
Scenario  Profit  

With no discount factor function  817  1539  173  0  0  145015 
1  2826  235  526  0.12  0.24  122838 
2  2801  428  690  0.17  0.29  119154 
3  2900  924  19  0.33  0.69  111141 
The optimal supplyside parameters under three scenarios were found using the BObased approach and results are summarized in Table 6. The discount factors of all MoD services are positive after employing discount factor functions, and is (as expected) higher than . As the disutility of high capacity MoD services at low discount factors increases in scenario 2, the service provider has to offer more discount and increase the fleet size to attract the passengers towards high capacity services, and thus the profit also decreases. Note that the fleet size of microtransit (capacity 10) is only 19 in scenario 3. This result represents the situation when the disutility of using microtransit at a low discount factor is extremely high, and thus the service provider finds it hard to make profits by running this service.
Scenario 







23.2  56.1  13.7  7.0  
1  61.0  13.5  19.2  6.3  
2  60.5  14.8  18.6  6.1  
3  62.2  26.4  6.4  5.0 
Table 7 shows the share of each mode in all scenarios. In a base case scenario when discount factor function is ignored, a large portion of the demand (around 70%) is served by high capacity MoD services. However, in scenarios with a discount factor function, a large fraction of demand (around 60%) is served by the ridehailing service. Therefore, the operating cost of the service provider appears to increase in order to make the same revenue. An increase in the operating cost and a higher discount factor of high capacity MoD services result in a lower profit. Not only the share of microtransit service decreases in scenario 3, its early demand shifts to the ridepooling service instead. This shift can be attributed to the similar values of for the ridepooling service in scenario 3 and for the microtransit service in scenario 2 (see Figure 8 and Figure 8).
Figure 9 shows a parallel coordinates plot of decision variables for which the BO framework evaluates the objective function in process of searching for the optimal values in scenario 2. For visualization convenience, , and are normalized by dividing with their maximum values (3000, 2000 and 1000 respectively). The red, blue, and yellow lines represent the best, the top , and all remaining solutions, respectively. The figure shows that the BO framework explores a broad domain of the solution space while exploiting the solution space with high expected profits more intensively.
In the experiments discussed above, we change both the discount factors and the fleet sizes for each MoD service to find the best solution using BO. It is difficult to infer the effects of discount factors on the mode share and the profit of MoD service providers. Therefore, we conduct an additional set of experiments for each scenario mentioned above. In these experiments, the fleet size for each MoD service is kept constant, which is the solution reported in Table 6. We use different multipliers in to increase or decrease the discount factors reported in Table 6 for each scenario. For example, if the multiplier we used for Scenario 2 is 0.75, then . Figure 11 shows the profit for three scenarios with different discount multipliers. The profit first increases and then quickly decreases as the discount factors increase. The highest profit is given by the discount factors found by BO. Figure 11 shows the mode share for Scenario 2 with different discount factor multipliers. We should note that: (i) As the discount increases for MoD services, the mode share of public transit service keeps decreasing as it gradually loses the advantage of low cost; (ii) The mode share of ridepooling and microtransit service increases as the discount factor increases, but the mode share of microtransit service increases faster. The reason is that the base discount factor for microtransit service is larger than ridepooling service, and thus the discount factor increases more as the discount multiplier increases.
4.5 An Application: Perride Tax
In this section, we illustrate how the proposed framework can be used to quantify the impact of systemlevel policy interventions. We consider a scenario where the government plans to impose a perride tax on only the ridehailing (capacity 1) MoD service. This type of tax has already been considered in California (arstechnica.com).
Consider a government policy to charge the MoD service provider $2 on every ridehailing trip. In the subsequent analysis, we assume that this cost is absorbed by the MoD provider with no increase to the passenger fare, but similar experiments where the cost is borne completely or partially by the traveler can also be considered. Using the discount factor function in scenario 1, we run the BO framework to optimize the supplyside parameters under this policy intervention. The optimal supplyside parameters, vehicle miles traveled (VMT) of MoD services, the ratio of the passenger miles traveled (PMT) and VMT, and the mode shares under tax and notax scenarios are shown in Table 8. Imposing the tax decreases the VMT of MoD services by and decreases the profit by .
In this tax scenario, running ridehailing service is not as profitable as before. However, if the fleet size of the ridehailing service is lowered, its demand is likely to shift to other travel modes. If the MoD service provider wants to induce this shift to move towards the other MoD modes and not public transit, the service provider needs to offer a higher discount for the shared services and thus the profit may further decrease. The optimal solutions of tax and notax scenario validate this hypothesis. A decrease in VMT of MoD services in tax scenario can also be attributed to a shift of demand from ridehailing services to high capacity MoD services. In fact, the increase in the share of public transit by also contributes to lower VMT of MoD services. The tax policy increases the revenue of public transit agency by $14340 which can further be recirculated for improvement in the level of service of the local public transportation to even further increase the transit ridership. One major benefit of the proposed framework is the ability to perform such scenario analysis and gain insights into the impacts of certain policy interventions.
Profit  
No tax  2826  235  526  0.12  0.24  145015  
$2 tax  1388  1668  31  0.15  0.67  100748  
VMT 





No tax  33835.7  1.11  61.0  13.5  19.2  6.3  
$2 tax  30288.3  1.17  36.4  50.6  5.1  7.9 
5 Conclusion
This paper addresses two critical issues in the design of a MobilityonDemand (MoD) system, namely modeling MoD services under endogenous demand and optimization of supplyside parameters (e.g., fleet size and fare). We propose a unified framework to design MoD systems that not only allows interaction of MoD services with other competing travel modes, but also optimizes supplyside parameters using Bayesian optimization (BO), a global optimization technique.
The proposed framework is calibrated using taxidemand data from Manhattan County in New York City and a stated preference survey. Demand is served by three types of MoD services (ridehailing, ridepooling, and microtransit of capacities 1,4, and 10) and public transit, and user preference for each travel mode is predicted using a choice model. The mode choice model was estimated using stated preference data collected in New York City. We conducted numerical experiments to show the numerical convergence of the integrated supplydemand system and show the application of the BObased approach in optimizing the supplyside parameters. The BObased approach outperforms the traditional random search optimization by sequentially taking steps in the direction of potential function improvement. The optimization problem maximizes the profit of MoD service provider in the case study, but the framework is more general and can be used to optimize various objective functions such as consumer surplus, vehicle miles traveled, as well as combinations of these metrics.
This proposed framework provides a potential toolbox for policymakers and planning agencies to evaluate the impact of MoD services on the urban transportation system and vice versa. For instance, we assess the systemlevel impact of a policy intervention where the government imposes a tax on ridehailing (capacity 1) trips without allowing the service provider to increase the fare. This policy results in an increase in transit ridership, a decrease in vehicle miles traveled of MoD services, and a reduction in the profit of the MoD service provider due to a need of higher discount factors of high capacity MoD services. Whereas the government can evaluate these policy implications qualitatively, the proposed framework provides a method to properly quantify the consequences of the policy intervention for different stakeholders (government, passenger, MoD service operator). Other such examples include quantifying the benefits for transit agencies if they plan to collaborate with MoD services to better serve first and last mile trips.
We now highlight future research avenues to strengthen the proposed framework. The integrated supplydemand system can be viewed as a fixed point problem. Whereas we illustrate the existence of its equilibrium in numerical experiments, an analytical proof or guarantees can increase its credibility and versatility. In addition, we used MNL to predict travel mode preferences because it integrates simply and well with the MoD simulation due to its closedform choice probability expression. However, we have shown that optimization of supplyside parameters is very sensitive to the parameters of the mode choice model. More complex random parameter (mixed logit) models can be used to possibly mitigate this sensitivity concern but those models require another layer of simulation in the prediction of choice probabilities, which can lead to numerical convergence issues. Thus, developing robust mode choice models with closedform choice probability expressions that more closely mimic the passengers’ decisionmaking mechanism is an essential avenue for future travel behavior research. In particular, mode choice models for a multimodal system with MoD services are required to go beyond just the level of service attributes (e.g., waiting time, travel time, and travel cost), and various perceptionbased latent factors (e.g., notion of reliability, variation in disutility of sharing rides with varying number of passengers, perception of discount) need to be incorporated. In addition, we use a simple transit route planner model in this study to estimate transit travel times and do not consider a fullblown transit assignment model. The framework can, however, be extended to include a more advanced transit routing model (e.g. hyperpathbased routing) or a fullblown transit assignment model. Including a transit assignment model will increase the complexity of the overall system, since there will be more dependencies between the supply and the demand (e.g. transit traveltimes will also depend on demand and vice versa).
Acknowledgements
The authors are thankful to the National Science Foundation Award No. CMMI1462289 for financially supporting this research. We are also thankful to two anonymous referees for their constructive comments.