{'fraction_non_alphanumeric': 0.03569088377426489, 'fraction_numerical': 0.015480383323777542, 'mean_word_length': 4.359047514539669, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 18, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Planning of integrated mobility-on-demand and urban transit networks

We envision a multimodal transportation system where Mobility-on-Demand (MoD) service is used to serve the first mile and last mile of transit trips. For this purpose, the current research formulates an optimization model for designing an integrated MoD and urban transit system. The proposed model is a mixed-integer non-linear programming model that captures the strategic behavior of passengers in a multimodal network through a passenger assignment model. It determines which transit routes to operate, the frequency of the operating routes, the fleet size of vehicles required in each transportation analysis zone to serve the demand, and the passenger flow on both road and transit networks. A Benders decomposition approach with several enhancements is proposed to solve the given optimization program. Computational experiments are presented for the Sioux Falls multimodal network. The results show a significant improvement in the congestion in the city center with the introduction and optimization of an integrated transportation system. The proposed design allocates more vehicles to the outskirt zones in the network (to serve the first mile and last mile of transit trips) and more frequency to the transit routes in the city center. The integrated system significantly improves the share of transit passengers and their level of service in comparison to the base optimized transit system. The sensitivity analysis of the bus and vehicle fleet shows that increasing the number of buses has more impact on improving the level of service of passengers compared to increasing the number of MoD vehicles. Finally, we provide managerial insights for deploying such multimodal service.


Introduction
The introduction of Mobility-on-Demand (MoD) services such as Uber, Lyft, and others as transportation alternatives has created many opportunities as well as challenges. On one hand, they provide a seamless mobility service with just a few taps on a cellphone application. On the other hand, they have increased congestion in densely populated areas due to an increase in the relocation and pickup trips made by the participating drivers in the network (Laris (2019)).
Furthermore, the transportation agencies envision the introduction of Autonomous Vehicles as a shared mobility service in the near future (Motavalli (2020)), which would lead to severe congestion in densely populated areas as predicted by various simulation studies (Levin and Boyles (2015), Fagnant et al. (2016), Levin et al. (2017)).
Public transportation, which can carry multiple passengers, is widely considered as a practical solution to the congestion problem by reducing vehicle-miles traveled (VMT) on roads (Aftabuzzaman et al. (2015)). However, due to its fixed routes and schedules, limited network coverage, and long waiting time, sometimes, it is less attractive to travelers in comparison to the auto mode.
The limited network coverage makes it difficult or sometimes impossible to access transit service in some areas. This inaccessibility problem is also known as the first mile/last mile (FMLM) problem for transit. The problem is commonly faced by travelers commuting from low-density areas where transit service is not available or less frequent because of the economic in-viability of providing such service.
A few studies have argued that the MoD service provided using autonomous vehicles would become a competitor of public transit mode (Chen and Kockelman (2016), Levin and Boyles (2015), Mo et al. (2020)), reducing its ridership, and other studies have even raised the question of whether urban mobility is possible without it (OECD (2015), Mendes et al. (2017)). On the contrary, Salazar et al. (2018) showed that the integration of the MoD system with transit system could help in achieving better results, such as a significant reduction in travel time, emissions, and costs as compared to the standalone MoD system. Through the current research, we also envisage an integrated MoD and transit system that aims to achieve the following potential benefits: 1. Providing fast and reliable mobility in low-density areas (i.e., by providing a first mile/last mile service) by means of characteristics of MoD service such as demand responsiveness, fleet repositioning, and reachability.
2. Allocation of resources from less congested areas to providing high-frequency transit service in congested areas through such integration.
3. Using existing transit infrastructure to reduce the number of vehicles needed for serving trips.
4. Reducing congestion and carbon emissions in the network, improving the mobility of travelers, and reducing the overall cost of providing transit service.
To achieve the above-mentioned benefits, we focus on the strategic planning of the transportation network that allows for intermodal trips with the first or last leg of the trips being served by the MoD service. To be specific, we try to answer questions such as which transit routes to operate when MoD vehicles are deployed to serve the FMLM connection, what should be the size of the vehicle fleet to be deployed, and what should be the frequency of operating transit routes. We attempt to answer these questions and make the following contributions through this article: 1. Propose a passenger assignment model that predicts the travel behavior of passengers in a multimodal network. This step extends the idea of a hyperpath transit assignment model proposed by Spiess and Florian (1989) to a multimodal transportation system with on-demand services.
2. Develop an optimization model to decide which transit routes to operate, frequency of operating transit routes, and MoD fleet size required to serve the FMLM of trips.

Develop a fast Benders decomposition implementation that uses efficient cutting planes to
solve the large instances of the current problem.

Conduct numerical experiments to show the efficacy of the proposed model and solution
methods and discuss the steps to implement such service in practice.

Related work
The passenger journeys that consist of auto, as well as transit mode, create a new mode of transportation known as intermodal or multimodal transportation. In intermodal transportation, different agencies are responsible for providing the service for different legs of the passenger journey involving multiple modes. On the other hand, in multimodal transportation, only one agency is responsible for providing service for the entire passenger journey traversed using multiple modes.
Recently, the studies are being focused on modeling the integration of MoD and transit service for future mobility. They can be divided into two categories: simulation-based and optimizationbased approaches. Under a set of assumptions on vehicle operations operations and dispatching strategies, the simulation-based studies simulate the passenger flow to assess the service quality of providing such mobility service (Gurumurthy et al. (2020)). By using a four-step travel demand simulation model, Levin and Boyles (2015) predicted that the transit ridership will decrease and the number of personal vehicles will sharply increase as a result of the repositioning of vehicles resulting in congestion on the network. Vakayil et al. (2017) developed a simulation model that accounts for transit frequency, transfer costs, and MoD fleet re-balancing to use MoD as the FMLM solution to the transit mode. Their results show that such an integrated system can reduce VMT in the network by up to 50%. Mendes et al. (2017) developed an event-based simulation model to compare the performance of the MoD system with the light rail system under the same demand patterns, alignment, and operating speed. They found that 150 vehicles with 12 passenger capacity would be needed to match the 39-vehicle light rail system if operated as a demand responsive system. Similar findings were also shown by the simulation model developed by Basu et al. (2018).
They showed that the introduction of MoD will act as the competitor of mass transit, however, to reduce congestion and maintain a sustainable urban transportation system, it cannot replace mass transit. Shen et al. (2018) also proposes and simulates an integrated autonomous vehicle and public transportation system based on the fixed modal split assumption. Using Singapore's organizational structure and demand characteristics, they propose to preserve high-demand bus routes while re-purposing low-demand bus routes and using shared MoD as an alternative. They found that the integrated system has the potential of serving the trips with less congestion, less passenger discomfort, and economically viable service. Wen et al. (2018) included mode choice and various vehicle capacities and hailing strategies in an agent-based model to provide insights into fleet sizing and frequency of transit routes for the integrated system. A few studies have used an optimization-based approach to developing an integrated passenger flow model. Salazar et al. (2018) developed a network flow model for intermodal service that couples the interaction between MoD and transit by maximizing social welfare. Using this model, they proposed a tolling scheme for this intermodal system that helps in reducing the travel time, costs, and emissions as compared to standalone vehicle mode. Liu et al. (2019) used Bayesian optimization to predict the mode choice of passengers in such a multimodal transportation system.
The above-cited studies show that an integrated MoD and transit system can provide an efficient mode of transportation that is sustainable, fast, eco-friendly, and economically viable. The design of such a system requires solving a multimodal transportation network design problem that can decide various aspects of MoD and transit modes. The problem of designing transit routes and their corresponding frequencies, which is commonly referred to as the Transit Network Design Problem (TNDP) or Line Planning Problem (LPP) in the literature, is itself a complex problem (Ceder and Wilson (1986), Baaj and Mahmassani (1991)). There has been a significant amount of research in modeling TNDP and developing solution algorithms for it. For a review on transit network design literature, we refer the interested reader to Guihaire and Hao (2008). Some aspects of the multimodal network design problem have been explored in a related research problem known as hub and arc location problem (Mahéo et al. (2019), Campbell et al. (2005a,b)). For example, Mahéo et al. (2019) proposed the design of a hub and shuttle public transit system in Canberra.
They formulated a mixed-integer program to design high-frequency bus routes between key-hubs, where the first mile or last mile of trips is covered by the shuttles. However, the hub and arc location problem have a major limitation of not able to capture passenger behavior in the transit network. Recently, a couple of studies have proposed models for the transit network design in the context of integrated MoD and transit system (Manser (2017), Pinto et al. (2020), Steiner and Irnich (2020)). Pinto et al. (2020) develops a bi-level optimization model to design a transit network integrated with MoD service. The upper-level optimization problem modifies the frequency of the transit routes and determines the fleet size of MoD service and the lower-level model simulates the passenger trajectories based on a simulation-based traveler assignment model. Due to the complexity of the model, they presented a heuristic approach to solving the current problem. Steiner and Irnich (2020) presents various aspects of this problem and develops a path-based mixed-integer programming model to decide which sections of the transit routes to operate and locate the transfer stops to allow for intermodal trips in the network. Due to an enormous number of possible paths in the network, they solve the current model using a branch-and-price approach.
The design of an integrated MoD and transit system is an important problem that can influence the future mobility of travelers. Recent studies have made important contributions to this complex problem but have several limitations. It includes designing networks by routing passengers on the shortest path using a multi-commodity network flow model and presenting solution algorithms not applicable to solving large-scale instances of the problem. We attempt to address these limitations and outline the motivation to pursue the current research in the following points: 1. Before designing the integrated system, we should understand how passengers would behave in an integrated system. It is common for studies to use the classic multi-commodity flow model to predict the behavior of travelers in the network design (Conforti et al. (2014), Chapter 4).
This may be true if passenger trajectories are completely influenced by the mobility provider.
However, this is certainly not applicable in the case of transit systems when passengers try to reduce the expected travel time based on wait time, travel time, and fare. Through this study, we extend the idea of hyperpath passenger assignment for a multimodal transportation system.
2. We develop a mixed-integer optimization model that incorporates the multimodal passenger assignment and evaluates various aspects of an integrated system. The optimization program is difficult to solve, and we need efficient techniques to solve this problem. For this purpose, an exact method based on the Benders Decomposition is proposed to solve the large-scale instances of the problem. The method improves the classic Benders decomposition strategy by precluding the infeasibility cuts and including new cuts, such as disaggregated cuts, multiple cuts, and clique/cover cuts.
The rest of the article is structured as follows. §3 discusses the notations and definitions used in this article. Then, we present the multimodal passenger assignment model, which is incorporated in the design model of the integrated MoD and transit system in §4. The solution algorithm to solve the design model is discussed in §5, which is followed by the results of numerical experiments conducted on Sioux Falls network in §6. Finally, conclusions and recommendations for future research are presented in §7.

Preliminaries and Background
In this section, we get familiarize ourselves with the notations and concepts to be used in this article. Let us begin by considering a multimodal transportation network characterized by a digraph G(N, A), where N denotes the set of nodes that includes road intersections N R , transit stops/stations N T , and centroids of traffic analysis zones Z 4 and A denotes the set of links. We associate every node i ∈ N in the network with exactly one zone Z(i). The set of links coming out and going into a node i ∈ N are denote by F S(i) = {(i, j) : (i, j) ∈ A} and BS(i) = {(j, i) : (j, i) ∈ A} respectively. Let d : N × N → R + be the distance function between two nodes in the network.
Depending on the mode, the links are also divided into three categories, namely transit, road, and walking links represented by A T , A R , and A W respectively. Let O ⊂ Z and D ⊂ Z be the subsets of centroids representing the origins and destinations respectively. The demand between various origin-destination pairs is represented by {d od } (o,d)∈O×D . The overall network can be divided into three sub-networks which are described below: 1. Transit network : The transit network is characterized by the subgraph G T (N T , A T ) which consists of a set of candidate transit lines/routes denoted by the set L. The terms "route" and "line" are used interchangeably throughout this article. Each line l ∈ L is composed of a set of stops N l T ⊂ N T which are connected by edges A l T ⊂ A T . The network also consists of transfer links A tr T between two nodes if the walking distance between those is less than the acceptable walking distance ζ (say 0.75mi), i.e., A tr T = {(n 1 , n 2 ) ∈ N × N : n 1 ∈ N l 1 T , n 2 ∈ N l 2 T for some l 1 , l 2 ∈ L s.t. l 1 = l 2 and d(n 1 , n 2 ) ≤ ζ}.
2. Road network : The road network is characterized by the subgraph G R (N R , A R ), where N R denotes the set of nodes and A R denotes the set of links in the road network.
3. Walking links: The walking links consists of access, egress, and mode transfer links. The access and egress links are defined as A a = {(n 1 , n 2 ) ∈ Z × (N T ∪ N R ) : d(n 1 , n 2 ) ≤ ζ} and A e = {(n 1 , n 2 ) ∈ (N T ∪ N R ) × Z : d(n 1 , n 2 ) ≤ ζ} respectively. Similarly, the mode transfer links are defined as A m = {(n 1 , n 2 ) ∈ N R × N T : d(n 1 , n 2 ) ≤ ζ} ∪ {(n 1 , n 2 ) ∈ N T × N R : d(n 1 , n 2 ) ≤ ζ}. The access and egress walking links connect the centroids of various zones with the road/transit nodes and vice-versa, whereas mode transfer links are used to transfer between nodes of various modes.

Costs
There is a subset of nodes in the network where passengers have to wait for the service. The collection of head nodes of links in the sets A a , A tr T , and A m constitutes the waiting nodes N w . Let us assume that c : A → R + and w : N w → R + denote the cost (e.g., walking time, in-vehicle time, and fare) associated with the links in A and wait time associated with the nodes in N w respectively.
The cost of links is known beforehand (and is computed by adding the travel time and possible fare multiplied by the value of time). On the other hand, the wait time depends on the availability of MoD or transit service.

Wait time computation
Unlike a personal vehicle, the MoD or transit service is not readily available, and passengers have to wait to access these services. So, it is important to quantify the expected wait time of these services, the computation of which is discussed below:

MoD service
We assume MoD operations in a network as a queuing system to compute the average wait time experienced by the passengers to access such service. The average wait time may not be justified for the planning of day-to-days operations but can be used to approximate the actual wait time experienced by the passengers for long-term strategic planning of the network, which is the focus of the current study. Therefore, we consider a stationary state of an MoD system, where the number of waiting customers C and vacant vehicles V are time-invariant. Using the Cobb-Douglas production function, the matching time between the customers and the vacant vehicles can be expressed as a function of C and V.
where, α 1 and α 2 are defined as the elasticities of the matching function and A is a parameter specific to a zone, which is a function of the market area divided by the running speed in that zone (Zha et al. (2016)). According to Little's law, the long-term average number of customers/drivers in a stationary system is equal to the long-term average arrival rate Q multiplied by the average wait time (w c /w t ) that a customer/driver spends in the system before being matched (Zha et al. (2016)).
Using (3) and assuming α 1 = α 2 ≈ 1 (Douglas (1972)), we can represent the stationary state (m c−v = Q) as below: Equation (5) shows that the average wait time of customers wait in a zone to access the MoD service is a function of the vacant number of vehicles. To achieve the desired level of service (i.e., average wait time), a transportation agency needs to provide V vehicles at any point in time.

Transit service
Let us now discuss the wait time computation to access transit service at the head node of an access or transfer link in the transit network. Let f : A T → R be the frequency of the transit line associated with various links of the transit network. Let g i (w) be the probability distribution function of the wait time for line i. According to Larson and Odoni (1981), for the passengers arriving randomly at a node, the probability density function of the wait time of line i is related to the headway or bus inter-arrival time distribution h i (h) as: To evaluate the wait time distribution, we make the following assumptions: Assumption 1. The inter-arrival time of a transit line i ∈ L follows an exponential distribution with rate f i .
Assumption 2. Passengers want to minimize the expected wait time to get to their destination.
Therefore, at any node, passengers waiting to be served by the transit service have selected a list of attractive transit lines that can help them to get to their destination.
Both assumption 1 and 2 are common in the transit assignment literature (e.g., see Desaulniers and Hickman (2007)). By using the assumption 1 and equation (6), one can evaluate the distribution function of the wait time g i (w) as: Proposition 1. (Spiess and Florian (1989), Gentile et al. (2005)) Assuming that a passenger waiting at node n ∈ N w is served by the set of attractive transit lines F S * (n) and let F = j∈F S * (n) f j .
With assumptions 1 and 2, the following holds: 1. The probability that a passenger would choose transit line i ∈ F S * (n) is given by 2. The expected wait time conditional to boarding line i ∈ F S * (n) is given by 3. The probability of wait time at node n follows an exponential distribution with rate F. Therefore, the expected wait time at stop n is given by EW n = 1 F .
Proof. See Appendix A.

Combined MoD and transit service wait time
Before discussing the computation of the expected wait time involving both modes, we need to make an assumption about the wait time distribution of MoD service by utilizing the value of the average wait time of MoD service calculated in equation (5).
Assumption 3. The wait time distribution of MoD service for passengers waiting at node n follows is the zone associated to node n and V Z(n) is the number of vehicles deployed in zone Z(n).
A passenger waiting at the head node of an access link faces the choice between MoD or transit mode. This is because the wait time of both services can vary based on the frequency provided, and a passenger will include one or both modes in their strategy to reduce the overall expected cost. This assumption simplifies the operation of MoD service as a transit service available at any stop of the network. The following proposition evaluates the expected wait time of that passenger.
Proposition 2. Given that the wait time for transit and MoD mode follow an exponential distribution with rate F and f M oD respectively and F = F + f M oD , the following holds: 1. The probabilities of taking transit and MoD are given by P M oD = f M oD F and P transit = F F respectively.
2. The expected wait time of the passenger departing from an access node n served by both MoD and transit service is given by To get more insights into the wait time computation, let us consider an example. Figure 1(a) shows an illustration of a multimodal transportation network. It consists of 2 zones, 6 nodes and 8 links as part of the road network, and 5 nodes and 10 links as part of the transit network. The transit network has 3 transit lines (color-coded) whose frequencies are shown in Figure 1  (assuming A 1 = 0.0017). The overall expected wait time at node Z 1 is 1.19 minutes which is less than 1.5 minutes by only considering transit service as part of the strategy.
We further use Proposition 1 and 2 to formulate the multimodal passenger assignment model.
For this purpose, we extend the frequency-based transit assignment model proposed by Spiess and Florian (1989) to a multimodal transportation system. Before moving forward, we must make the following assumptions: Assumption 4. (a) Ridepooling is not allowed, i.e., the MoD service serves one passenger at a time.
(b) The transit lines are assumed to have unlimited capacity.
(c) Passengers want to reduce their expected generalized travel cost consisting of travel time, wait time, and fare to get to their destination.
The ridepooling problem requires matching of customers using a specific algorithm. This is an important aspect to accurately estimate the cost of day-to-day operations. Nevertheless, ignoring ridepooling will give us an upper bound on the number of vehicles required to serve various zones.
The modeling of passenger behavior while incorporating the capacity constraints (congestion) is a difficult problem. The congestion is important to consider since it causes denied boarding, which leads to increased wait time, travel time, and discomfort. Several authors have tried to include congestion into frequency-based transit assignment models through various approaches, namely, discomfort function (Spiess and Florian (1989)), effective frequency (De Cea and Fernández (1993), Cominetti and Correa (2001), Cepeda et al. (2006), Leurent et al. (2014)), and failure-to-board probabilities (Kurauchi et al. (2003)). Despite the effort, there is no tractable closed-form of congested frequency-based transit assignment model. On the other hand, it would not be ideal to include transit vehicle capacity constraints into the assignment program (e.g., in Szeto and Jiang (2014)) because doing so may lead to unrealistic passenger behavior, which previous studies on congested frequency-based transit assignments were trying to avoid. Therefore, we use an uncapaciated assignment for the design problem. Assumption 4(c) is a common in the assignment literature. The relaxation of above assumptions are research topics in their own right, therefore, a discussion on possible ways to relax them is provided in §7. To proceed further, let us define a variable {g ik } i∈N,k∈D as below: Furthermore, let us denote v ak and W ik as the flow of passengers on link a ∈ A and waiting at node i ∈ N w resp. destined to k ∈ D. The assignment optimization program is presented below: The assignment program (10)  and multiplies it with the number of passengers waiting at that node. Note that the probability of selecting an option (MoD or transit line) is calculated in Proposition 2.

Design of an integrated MoD and transit system
In this section, we present an optimization model incorporating the assignment program proposed in previous section for the design of an integrated MoD and transit system. The optimization program is formalized as a Mixed Integer Non-linear Program (MINLP). In this model, we determine which transit routes to keep operating among the current transit routes in the city network, Decision Variables decide the optimal frequency of those operating routes, and finally, determine the fleet size of vehicles required to provide MoD service in various zones. Note that one can also include new candidate transit routes as part of the design plan. The sets, parameters, and decision variables for the optimization model are summarized in Table 1. The 0.01 vehicles in the set of possible number of vehicles to be deployed in a zone is a dummy element to represent that no vehicles are assigned in a zone, and that zone can be served by transit service only. This will help design an efficient solution algorithm for the design problem (see Proposition 4).
The design of an integrated transit and MoD system should consider both passenger and operator perspectives. The operator's perspective is to provide the service at minimum cost, and the passengers' perspective is to minimize the overall cost of travel (including travel time, wait time, and fare). Based on these perspectives, the design optimization model is presented as (11).
The objective function is the sum of the total expected travel cost and wait time experienced by the passengers in the network. The mapping B : L × Θ → N used in (11b) is defined as 2t a , which describes that the number of buses required to provide frequency f ∈ Θ for a line l ∈ L is equal to the product of the frequency and round trip travel time.
(11b) constrain the total number of buses needed to be less than or equal toB, which can be evaluated for a given budget. (11c) describes the flow conservation constraints at every node for every destination. For a given MoD and bus fleet assignment, (11d)-(11e) describe the passenger flow on each link based on the frequency of the bus route and MoD service. A frequency value can be assigned to a route if that route is kept operating as constrained by (11f). (11g) describe that exactly one of the fleet sizes can be adopted for each zone. (11h) constrain the required number of vehicles to be less than or equal toF . Finally, (11i), (11j) and (11k)-(11m) are the non-negativity constraints of the flow, wait time being free variables, and binary constraints of design variables respectively. One can also incorporate other constraints related to the budget of operating MoD and transit service but for the sake of simplicity, we do not include them here.
The optimization program (11) is a mixed-integer non-linear program (MINLP). The nonlinearity arise from the constraints (11d)-(11e). It is computationally difficult to solve this program for large instances, which can be attributed to the integer constraints (11k)-(11m) and the bilinear constraints (11d)-(11e). The bilinear constraints are particularly difficult to handle due to the nonconvex nature even if the integrality constraints of the involved variables are relaxed. Fortunately, in this case, the non-convexity arises due to the product of continuous and binary variables, which can be exactly relaxed by employing McCormick relaxations. Let t f aik = y l(a)f W ik , ∀f ∈ Θ, ∀a ∈ F S(i) : a ∈ A T , ∀i ∈ N w , ∀k ∈ D and ω ink = N Z(i)n W ik , ∀n ∈ Ω, ∀i ∈ N w ∩ N R , ∀k ∈ D.
Further, let us assume that there exists a finite upper and lower bound on the variable W ik , i.e., W ik ≤ W ik ≤ W ik . Then, t f aik and ω ink can be expressed as the set of linear constraints (12a)-(12d) and (13a)-(13d) respectively:

Solution methodology
After reformulating.0 the bilinear constraints (11d)-(11e), the resulting model is a Mixed Integer Linear Program (MILP). The program is still difficult to solve efficiently for large instances.
However, the structure of the problem allows us to use decomposition techniques such as Benders decomposition to efficiently solve it. In this section, we present the details of the Benders reformulation for this problem, along with the proposed algorithmic enhancements.

Benders Reformulation
Benders decomposition (Geoffrion (1972)  To speed up the process of Benders decomposition by avoiding feasibility cuts, we need to put some restrictions so that the subproblem (14) is always feasible.
Proposition 4. Given that 0 / ∈ Ω and G R (N R , A R ) is connected, then X SP is non-empty for any given feasible value ofx,ŷ,N .
By applying an outer linearization procedure to the inner (sub) problem of the original problem, we can restate it as (16), which is referred to as the Benders Master problem (MP).

Classic Benders decomposition implementation
The issue with the Benders reformulation is that there could be a large number of extreme points of the polyhedron associated with the feasible region of DSP, therefore, one applies an iterative process of solving two problems, namely, the relaxed master problem (RMP) and the subproblem (SP) repeatedly. The relaxed master problem is the master problem with constraints (16f) being defined only for a subset of extreme points, i.e., K ⊂ K. The overall implementation of the classic Benders Decomposition is summarized in Algorithm 1. We start by finding the feasible value of (x 0 , y 0 , N 0 ). This can be done by solving (16) without (16f) and including a constraint η ≥ 0. Then, in each iteration t, the algorithm solves RMP with the given set of extreme points and then SP with the current value of (x t , y t , N t ). Since RMP is relaxation and SP is solved for a feasible value (x t , y t , N t ), they provide a lower bound and upper bound respectively to the original problem. The subproblem also provides inequalities (optimality cuts) to strengthen the formulation of RMP in each iteration. Thus, it is guaranteed to have non-decreasing lower bounds. In our case, there are no feasibility cuts since our subproblem is always feasible (Proposition 4). The algorithm terminates when both the upper bound and lower bound are close to each other.
2: while U B − LB > do is the tolerance parameter 3: Set t = t + 1. Solve RMP (16) and obtain its optimal solution (x t , y t , N t ).

5:
Solve SP (14) for (x t , y t , N t ), obtain dual solutions and append that to K .

Enhanced Benders decomposition implementation
The classic Benders decomposition may take prohibitive computational effort to converge, thus making it difficult to solve the problem for large instances. The slow convergence can be attributed to the low strength of the optimality cuts, degeneracy in the subproblem, no guarantee of nondecreasing upper bounds in each iteration, or not formulating the problem "properly" (Saharidis and Ierapetritou (2010), Tang et al. (2013), Magnanti and Wong (1981)). To accelerate the Benders decomposition algorithm, we make use of several enhancements that are described below:

Use of multiple cuts via disaggregated cuts
For this design problem, we can further utilize the decomposable structure of the Benders subproblem (11) as it is decomposable for each destination k ∈ D. That is, we can solve several (smaller) subproblems and generate multiple optimality cuts for the master problem. The disaggregated cuts have a higher probability of finding facet-defining inequalities characterizing Π. For this purpose, we modify RMP to allow for the disaggregated cuts as (17): Note that by adding the disaggregated cuts for every destination, we can get back the optimality cuts defined in the classic Benders relaxed master problem.

Use of multiple cuts via multiple solutions
To further improve the convergence of the algorithm, Beheshti Asl and MirHassani (2019) used a strategy known as multiple cuts via multiple solutions. When solving the RMP, any commercial solver such as AIMMS or GUROBI can be asked to generate multiple solutions of an integer program (optimal as well as suboptimal) by using pool solution option. These multiple solutions can be used to generate multiple classic (16f) or disaggregated cuts (17c) to be added in next iteration of RMP. This strategy is expected to decrease the overall iterations and possibly the solution time of the algorithm.

Use of clique/cover cuts
Due to the limited availability of bus and vehicle fleet, one can use the clique/cover cuts to tighten the feasible region of the master problem.
Proposition 5. For every n ∈ Ω, if F n < |Z| then the clique inequality z∈Z N zn ≤ F n is valid for (11).

Proof. See Appendix A.
If for any n ∈ Ω, we have Proof. See Appendix A.
To generate some of the minimal cover cuts, one can use the heuristic given in Algorithm 2. In this algorithm, for each frequency f ∈ Θ, we keep the list of lines G for which the number of buses required to provide the frequency f does not exceedB. Then, any line which is not in G, along with G forms a minimal cover.
Algorithm 2 Cover cut generation heuristic 1: procedure 2: Compute the value of mapping B(l, f ) for all (l, f ) ∈ L × Θ. for l ∈ L\G do 13: append C to CC return CC Furthermore, one can use other efficient techniques to produce maximal clique or minimal cover cuts for the problem.

Other recommendations
One of the problems with the Benders subproblem (14) is that it assumes the value of W ik as a given upper bound. The value of W ik is a big-M introduced to relax the non-linearity in the original model. If the value of the big-M is not chosen properly, then one can face serious issues with the convergence of the algorithm. For example, choosing W ik < W ik can make the subproblem infeasible, and choosing W ik too high would generate weak optimality cuts, which would increase the computational time of the algorithm. One way to avoid this issue is to solve the assignment problem (10) for given design variables (x, y, N ) and compute the optimal value of W ik and use that as an upper bound. Further improvements in the Benders decomposition method can involve the use of pareto-optimal cuts proposed by Magnanti and Wong (1981). They help in avoiding the generation of multiple optimality cuts for a degenerate subproblem. We tried this strategy, however, we did not find any significant improvement in the solution time using these cuts, therefore, we do not discuss it here. Finally, when RMP is loaded with a large number of cuts we recommend removing the non-active cuts from the model by checking the dual value. There is no guarantee that they will not be generated again, but it will be faster to solve the RMP. The overall steps of the Benders implementation with possible acceleration techniques are summarized in Algorithm 3.
2: Prepare the master problem with clique and cover inequalities.
3: Assume an initial feasible value (x 0 , y 0 , N 0 ). Solve the SP (14), obtain the optimal dual solutions and append that to the set K k , ∀k ∈ D.
4: while U B − LB > do is the tolerance parameter 5: Set t = t + 1. Solve RMP (16), obtain its optimal solution (x t 0 , y t 0 , N t 0 ) and other optimal/suboptimal solutions {(x t s , y t s , N t s )} 1≤s≤l , where l is specified by the user.

Computational results
In this section, we present the computational study based on the model (11), (16), and acceleration techniques presented in §5.3. We start by describing the details of the experiment used to show the application of the proposed method. Then, we present the details of the network design results in §6.2, which is followed by the comparison of the computational performance of the solver, the classic Benders implementation, and the enhanced Benders techniques described in §6.3. Finally, we discuss the results of the sensitivity analysis on two important parameters in the model, namely, the available fleet of busesB and vehiclesF in §6.4, which is followed by the comparison of the performance of optimized existing transit system and proposed integrated system in §6.5.

Experiment details
The computational experiments are based on the Sioux Falls road and transit network. The road network has 24 nodes, whereas the static transit network has 384 stops. A walking distance of 0.5 miles is used to create walking links. An illustration of the two networks is shown in Fig

Network design results
We solve the network design problem (11) for the instance explained in §6.1. The selected transit routes with their optimal frequency are given in Table 3. Out of 12 candidate routes, 6 routes are kept operating. The transit network with active and inactive routes is shown in Figure   3. We observe that most of the routes are located in the central region of the network. All the routes have been assigned the highest frequency, i.e., 12 buses/hr, except route 8, which has been assigned a frequency of 3 buses/hr. To provide this service, 69 buses are required. The average number of vehicles deployed in each zone is given in Table 4. In the optimal allocation of vehicles, it is decided not to deploy any vehicles in 10 zones out of 24 zones. Most of the zones have been allocated 200 vehicles providing an average wait time of 3 minutes. We further observe that the vehicles are deployed in the outskirt zones of the network where transit routes are not located.   12  5  2  No  --3  Yes  12  5  4  Yes  12  5  5  No  --6  No  --7  No  --8  Yes  3  20  9 No --10 No --11 Yes 12 5 12 Yes 12 5 The total time spent in the system is equal to 11,301 passenger-hours including 8,673 passengerhours of travel time on various links, 1,881 passenger-hours of wait time spent on the transit network, and 747 passenger-hours of wait time spent on the road network. We found that more passengers take transit than MoD service. The share of passengers using the road, transit, and multimodal service are 23 %, 61 %, and 16 % respectively. The passenger flow on various links and wait time  1  200  3  13  200  3  2  100  6  14  200  3  3  200  3  15  200  3  4  200  3  16  --5  --17  --6  --18  200  3  7  200  3  19  200  3

Computational performance
In this section, we compare the computational performance of various models and implementation techniques. We consider the following approaches to compare: 1. Solving model (11) using Gurobi solver 2. Solving model (16) using Gurobi solver The computational performance of every method is shown in Table 5. available vehicles as well as buses. This is because more passengers take the transit mode as more buses are made available. The overall expected travel cost also reduces with the increase in the bus and vehicle fleet. However, the effect of an increase in the vehicle fleet size of more than 5,000 is negligible. Figure 6 shows the mode share as a function of the available bus and vehicle fleet size.
As expected, the transit and MoD share increase with the increase in the available bus and vehicle fleet size respectively. The share of multimodal service increases with the number of vehicles and buses up to 5,000 and 75 respectively but declines after that. The decline in multimodal share is because of the reduced wait time for both services, which drives passengers to use single mode rather than multiple modes. 6.5. Comparison of optimized base transit system with proposed integrated system In this section, we present a comparison of the operation of the "optimized base transit system" corresponding to the existing transit system with optimized frequencies versus the design of the integrated system evaluated in §6.2. For the optimized base case, we solve the optimization program (18) for the instance described in §6.1. The results of optimized frequencies of various routes are given in Table 6. The network provides an average wait time of 18 minutes to the passengers.  3  6  10  4  12  5  5  3  20  6  2  30  7  3  20  8  3  20  9  3  20  10  2  30  11  12  5  12  6  10 The results comparing the performance of the optimized base transit system and integrated system are provided in Table 7. The base transit system has 12 two-way bus services operated by a bus fleet of 69 buses, whereas the new integrated system has 6 two-way bus services operated by 69 buses. Along with 69 buses, the new integrated system deploys 3,000 vehicles to serve the demand.
The deployment of these extra vehicles can be costly to the transportation agencies. However, they provide several benefits. First, the optimized base transit system is not able to serve 13 % of the demand due to the non-availability of transit service in 2 zones in the network. On the other hand, the first mile and last mile of these zones are covered by vehicles in the integrated system. Second, the average in-vehicle travel time of passengers using the integrated system is only 14.43 minutes in comparison to the 21.16 minutes for passengers using the base transit system. However, the average wait time of the integrated system users is increased slightly in comparison to the base transit system. This is due to the increased number of transfers to access MoD and transit service. For implementing the proposed model in practice, we need to follow the following procedure.
First, we divide the region into zones. Second, we collect the peak hour demand data in the region.
Third, solve the proposed design model for varying fleet sizes. This step will be similar to the sensitivity analysis given in Section 6.4. This analysis will help us decide the optimal fleet size of buses and vehicles for our service. It will also provide us with the allocation of vehicles and buses for different zones and bus routes respectively. This allocation is designed for peak hours. We can reduce the vehicle operation in non-peak hours. For the bus service, scheduling needs to be performed to publish a schedule for the service.
The numerical results show that the proposed integrated MoD and transit system provides several benefits. First, it relieves congestion in the city center by allocating the transit frequency from low-population density areas to congested areas. Second, it improves the level of service and access of transit passengers, which leads to an increase in the overall share of transit trips in the city. Third, improvement in transit passenger mobility further leads to various socio-environmental benefits.

Conclusions and Future Research
Advances in mobility services have paved the way for the development of a new type of MoD service, which can help in serving the first mile and last mile of transit trips. Such a system requires rethinking the design of a transit system that allows for intermodal trips with MoD as the first or last leg of trips. We developed a mixed-integer non-linear program (MINLP) to design such system. The MINLP was relaxed to a mixed-integer linear program with the help of McCormick relaxations. To solve the resulting MILP model efficiently, we proposed the Benders decomposition method with several enhancements. These enhancements include the use of disaggregated cuts, clique/cover cuts, and multiple cuts via multiple solutions. The numerical results show that disaggregated cuts with clique/cover cuts and multiple cuts via multiple solutions are efficient techniques to solve the current problem. Furthermore, the experiment results on the Sioux Falls network show that the congestion in the city center is improved with such design as most of the passengers were found to take the transit in that region. The sensitivity analysis on bus fleet size and vehicle fleet size shows that the passenger hours spent in the system as in-vehicle time and wait time reduces with an increase in the number of available buses and vehicles. The share of multimodal service was observed to be highest for the vehicle fleet size and bus fleet size of 3,000 and 75 respectively.
We also compared the proposed integrated system with the optimized base transit system. We found that the integrated system can be costly due to the deployment of vehicles, but it reduces the passenger in-vehicle time and serves more demand than the optimized base case.
This research can be expanded in multiple directions. First, ridepooling was not allowed in the current study. Further research is needed to explore the ideas of including the matching of passengers for ridepooling, which will further reduce the size of the vehicle fleet required to provide the service. Second, better calibration of parameter A used in the wait time computation of MoD service is needed. The data from ridehailing services can be used for this purpose. Third, stochasticity from both supply-side and demand-side should be studied in the future as stochasticity can be good or bad for a transport system (Wang and Wu (2021)).
Summing over all the lines F S(n) gives us the expected wait time at the stop, i.e., where, i γ i (w) is the probability density function of the wait time at stop n.
i∈F S(n) Therefore, the expected wait time at stop n if given by EW n = 1 F

Proof of Proposition 2
The probability of taking transit is given by: Similarly, the probability of taking MoD is given by P M oD = F F . The expected wait time of the passenger departing from an access node n is given by: Proof of Proposition 3 To prove this, we need to show that (a) proj v,W X SP ⊆ X M A and (b) X M A ⊆ proj v,W X SP . Let us first start by proving (b). Let (v, W ) ∈ X M A . For all a ∈ F S(i) : a ∈ A T , ∀i ∈ N w , ∀k ∈ D, we have v ak ≤ f a W ik . Let y l(a)f = 1 if the frequency of the line associated to arc a is f ∈ Θ and 0, otherwise. Let t f aik =ŷ l(a)f W ik , then f a W ik = f ∈Θ fŷ l(a)f W ik = f ∈Θ f t f aik , which is same as (14c). Also, t f aik =ŷ l(a)f W ik can be expressed as (14d)-(14f) and (14l). Using a similar argument, we can show that for all a ∈ F S(i) : a ∈ A R , ∀i ∈ N w , ∀k ∈ D, the inequality v ak ≤ f a W ik can be expressed as (14g)-(14j) and (14m). This shows that X M A ⊆ proj v,W X SP . To prove part (a), let (v, W, ω, t) ∈ X SP , then using Fourier-Motzkin elimination, we have (v, W ) ∈ X M A (Conforti et al. (2014, Chapter 3)).

Proof of Proposition 4
The set X SP can be empty in two cases i.e., when there is no flow balance ( k∈D i∈N g ik = 0) or there does not exist a directed path from any node i ∈ N to any destination k. However, it is not possible to have any of these cases because from the definition of g ik , we have k∈D i∈N g ik = 0 and since there is always at least 0.01 vehicle assigned to all the zones and the road network is connected, there always exists a path from any node i ∈ N to any destination k ∈ D. Therefore, X SP = φ.
Proof of Proposition 5 Due to limited fleet available, one can allocate n ∈ Ω vehicles in at most F n zones. Proof of Proposition 6 The proposition follows from the definition of minimal cover that we cannot provide the number of buses required in the minimal cover.