how aerosol particles exhaled from passengers sitting in any seat would travel through the interior of a campus bus (see Zhang et al., 2020) . There physics-based computer models of aerosol dispersion were developed under different social-distancing and bus-operating conditions and the results were validated with experiments using water vapor in a real bus. Following the results, all bus riders and drivers must wear face coverings, and the capacity of buses also needs to be reduced to half of the original capacity (i.e., approximately 78 seats on each bus) or even lower, while the length of each trip needs to be limited to 15 minutes or less, to minimize riders' possible exposure to the virus.
Shortening rides presents a challenge, as the pre-COVID campus bus system at UM uses long routes designed to minimize transfers. Some routes (e.g., the original Commute South, Commute North and Northeastern Shuttle) take 45 minutes to 1 hour to travel even without including the passenger loading/unloading time at each bus stop. Simply breaking down the original long routes into shorter ones will significantly increase the number of buses needed to cover all the original demand, especially when only half bus capacity can be utilized. Meanwhile, it is not possible to add more buses and drivers due to the high cost of purchasing a campus bus (approximately $800,000 each) and the difficulty of retaining and hiring qualified drivers during the pandemic (Dolen et al., 2020) .
In this paper, we develop a mathematical optimization model and a simulation platform to design a hub-and-spoke system that relies on short, direct routes with fewer stops. The shorter routes not only limit all the rides to 15 minutes or less, but also enable more frequent service using the same number of buses, and therefore reduce passenger wait time as well as the number of passengers on each bus. We validate the results via simulation by generating a diverse set of scenarios that take into account spatial-temporal variations of demand, changes of social-distancing requirements, and random bus broken-down (or equivalently driver no-shows). The simulation results show that our new design of the campus bus system can satisfy regular peak demand using the same number of buses but only 50% or even fewer seats in each bus, while all the riders can find a bus stop within 5 minutes of walking distance and 70% of riders still take only one transit for completing their trips. For 35% riders, the average waiting time is less than 5 minutes in all the stops and for 70% riders, the average waiting time is less than 10 minutes. Through simulation, we also identify the routes that have higher utilization on average or during peak hours, and the bus stops that have longer average waiting time. Based on these results, we develop contingency plans of adding more buses to these routes from others in case of some bus failures, to improve the system resilience.
Our work is broadly related to prior work on Transit Network Planning (TNP), including Transit Network Design (TND) (see, e.g., Yu et al., 2012; Ouyang et al., 2014) , Transit Network Timetabling (TNT) (see, e.g., Ibarra-Rojas and Rios-Solis, 2012; Fonseca et al., 2018) , Vehicle Scheduling Problem (VSP) (see, e.g., Naumann et al., 2011; Yan et al., 2012) , and Transit Network Evaluation (see, e.g., Chandrasekar et al., 2002; Zhang et al., 2019) . The methods we use for redesigning campus bus routes are also closely related to School Bus Routing Problem (SBRP) (see, e.g., Park and Kim, 2010) and hub-spoke system design (also referred to as Hub Location Problems, see, e.g., Campbell and O'Kelly, 2012) .
TND determines the lines, types of vehicles, and stop spacing to meet population's movement requirements. Yu et al. (2012) use a direct traveler density model (first proposed in Yang et al., 2007) to maximize demand density of route under resource constraints. They propose to use ant colony to solve the model and conduct numerical study with data from Dalian, China. Ouyang et al. (2014) represent passenger demand as a continuous function over a geographical space to address heterogeneous densities across a city. They propose a hierarchical system, consisting of a main grid and several local grids, to minimize a weighted sum of the total bus idle time, traveled distance, vehicle usage, passenger walking, waiting at stops, in-vehicle travel time, and penalties for transfers. Kim and Schonfeld (2014) propose a probabilistic analytical model to design and synchronize a bus transit network by integrating conventional and flexible bus services and coordinating bus arrivals at transfer terminals. We refer to Guihaire and Hao (2008) for a comprehensive review of TND problems.
TNT determines the departure and arrival times for a set of trips along the stops of the transit network in order to optimize a specific objective function. Ibarra-Rojas and Rios-Solis (2012) consider bus timetable generation for large bus networks where passenger transfers are favored.
The authors formulate an integer program to facilitate passenger transfers and avoid bus bunching in the network. They analyze the structural properties of the feasible solution space and propose a meta-heuristic method. Extending the results in Ibarra-Rojas and Rios-Solis (2012), both Fonseca et al. (2018) and Carosi et al. (2019) integrate timetable design and vehicle scheduling problems, and specifically, Fonseca et al. (2018) formulate a bi-objective mixed-integer programming model to minimize the transfers and operational cost, solved by meta-heuristic; Carosi et al. (2019) formulate the integrated problem as a multicommodity flow problem and propose a diving-type approach to obtain heuristic solutions.
For public transit systems, VSP assigns trips to each bus to meet travel demand. Forbes et al. (1994) propose a mathematical formulation for solving multi-depot scheduling problem (MDSP) and show the equivalence between MDSP and the multicommodity flow problem. The authors develop an exact algorithm for solving large-sized MDSPs by (i) solving a pure network flow problem, (ii) solving a linear relaxation of the MDSP starting from the optimal solution obtained in step (i), and then (iii) re-optimizing an integer solution via branch-and-bound. Different from the connectionbased network used by Forbes et al. (1994) , Kliewer et al. (2006) use a time-space-based network, which leads to a significant size reduction of the corresponding mathematical model and allows solving large-scale instances using commercial solvers. Recent studies on VSP (see, e.g., Naumann et al., 2011; Yan et al., 2012) Recently, several papers consider the integration of TNT and VSP, such as Wu et al. (2019) and Dou and Meng (2019) . Wu et al. (2019) propose a bi-level model for bus schedule coordination design, allowing passengers to modify routes multiple times. The upper-level problem is a mixed-integer program minimizing total system cost and the lower-level problem is a route choice model.
They develop a heuristic algorithm to solve this model and show that a more cost-effective schedule can be obtained by taking rerouting behaviour into consideration. Dou and Meng (2019) propose a mixed-integer program for a feeder bus timetabling problem to minimize transfer waiting, transfer failure cost and bus operational cost. This model can handle both transfer passengers and local passengers whose behavior is constrained by bus capacity during peak hours. They develop a hybrid artificial bee colony (ABC) algorithm to solve the problem.
To evaluate a given transit network, Chandrasekar et al. (2002) 
There has been a growing body of literature in the context of SBRP since the first paper by Newton and Thomas (1969) . However, there is no single dominant approach for SBRP. The reason is that each real-world SBRP problem has specific assumptions and constraints and thus the proposed approaches among the literature are highly problem-dependent (see, e.g., Li and Fu, 2002) . A recent survey paper by Park and Kim (2010) provides a comprehensive review of SBRP. They decompose the SBRP as sub-problems involving data preparation, bus stop selection, bus route generation, school bell time adjustment, and bus scheduling. They summarize works on SBRP and present (i) the sub-problems (see Table 1 in Park and Kim, 2010) , (ii) the various problem settings, including the assumptions, constraints, objectives and so on (see Table 4 in Park and Kim, 2010) , and (iii) mathematical models and the corresponding solution methods (see Table 5 in Park and Kim, 2010) among the literature. Mixed-integer programming models are ubiquitously adopted to formulate SBRP as well as exact or heuristic-based algorithms (see, e.g., Fügenschuh, 2009; Sun et al., 2018; Bertsimas et al., 2019; Hulagu and Ç elikoglu, 2019) . Fügenschuh (2009) considers integrated coordination of school starting time and school bus system while allowing student transshipment from a route to another. The model is efficiently solved by a branchand-cut algorithm with preprocessing techniques and valid inequalities. Bertsimas et al. (2019) develop bi-objective routing decomposition (BiRD) algorithm to design school bus route. They 
The early work on hub-spoke system design is initiated by O'kelly (1986); O'Kelly (1986), and hubspoke system design is also referred to as Hub Location Problems (HLP). From then on, hub-spoke network systems have been attracting the attention of researchers in telecommunication (see, e.g., Carello et al., 2004; Kim and O'Kelly, 2009) , cost allocation (see, e.g., Matsubayashi et al., 2005) , air transportation (see, e.g., Aykin, 1995; Bryan and O'kelly, 1999) , and good delivery networks (see, e.g., Wasner and Zäpfel, 2004; Yaman et al., 2007) . For example, Masson et al. (2017) adopt a hub-spoke system in logistic transportation. They consider a two-tiered problem where the the goods are transferred to cities from distribution centers in the first-tier problem and goods are delivered to customers by smaller vehicles in the second-tier problem. They develop a mixedinteger programming model and an adaptive neighbourhood search to solve this problem. Campbell et al. (2002) and Farahani et al. (2013) provide comprehensive review on mathematical models, solution methods, and applications of hub-spoke systems. Campbell and O'Kelly (2012) reflect on the origins of hub location research, especially in transportation, provide some commentary on the present status of the field, and discuss some promising directions for future effort.
We first describe a mathematical programming formulation for designing a hub-and-spoke system for operating campus buses at UM. We identify locations for bus stops and also the "hubs" using real travel demand data that have information of potential passengers' origins and destinations. In Figure 1 : An illustration of the hub-and-spoke system for operating UM buses in Ann Arbor campuses Figure 1 , we show an abstract hub-and-spoke system with three hub locations (i.e., the squares) in the UM North Campus, Central-South Campus and Medical Campus, respectively, where the most number of transits take place during peak hours in regular semesters. A Campus Connector route will connect these hubs with very frequent buses. Shorter routes are designed to connect other popular bus stops within each of the three campuses (i.e., sub-regions) from and to each hub. The goal is to reduce the total travel time of all the buses and to increase the bus frequency on each (shorter) route, so that we only use 50% or less capacity of each bus to meet the same amount of peak-hour travel demand.
We define notation and formulate an integer programming model for designing bus routes originating from and returning to each hub location as follows. Let I, J, K be the sets of all physical locations of bus stops, routes and stops in each route for each sub-region that contains one hub (e.g., each UM campus), respectively. We want to assign one bus stop to each stop in each route.
For notation convenience and without loss of generality, all routes contain the same number of stops, |K|. The value of |K| is set as the maximum number of stops in any route in a real bus system.
Let c j be the cost of operating route j for all j ∈ J, t i 1 ,i 2 be average travel time from stop i 1 to stop i 2 for all i 1 , i 2 ∈ I ∪ {hub}, andt stop be average loading/unloading time at each stop. We define binary variables x j , y ij , u j i,k , u j hub,k ∈ {0, 1} for all i ∈ I, j ∈ J, k ∈ K such that x j = 1 if we operate route j, y ij = 1 if we assign bus stop i to route j, u j i,k = 1 if bus stop i is the k-th stop in route j, and u j hub,k = 1 if the hub is the k-th stop in route j, respectively. A mixed-integer programming model is given below to design routes for each sub-region that contains one hub, to satisfy the travel demand from/to the hub. has not been filled. Constraints (1g) enforce a time limit T (e.g., 15 minutes given by COVID-19 guideline) for completing each route (including the riding and loading/unloading time), where |K| − |K| k=1 u j hub,k is the number of stops in route j for all j ∈ J. The objective function (1a) involves bilinear terms u j i 1 ,k u j i 2 ,k+1 , which cause the main computational difficulty. Because variables u j i 1 ,k and u j i 2 ,k+1 are binary-valued, we can provide an exact reformulation of Model (1) using McCormick envelopes (McCormick, 1976) : 
Next, we describe procedures for designing bus frequency given all the route information. Let J be the set of all new routes, and C be the number of seats on each bus. (For notation convenience, we use the same capacity C for all buses.) For each new route j ∈ J , we first extract the single-trip time T j on Google Map API, and then examine all the old routes that pass through each stop in route j, to calculate an approximate frequency provided by buses in the old system that the new route j needs to match. We denote this bus frequency as f old j (in minute) for each new route j and the corresponding transportation capability as HC old j = (60/f old j ) × C, representing the maximum number of passengers that the old bus system can ride during an hour on the new route j based on its frequency f old j . We aim at achieving the same or higher hourly transportation capability using s% of the total bus capacity C on each new route. Specifically, the number of buses we assign to new route j during each hour can be determined as
and the new bus frequency and hourly transportation capability of route j become
We provide several examples through Tables 1 and 2 to illustrate how the formulas work. 
We validate and improve the performance of the proposed route and bus frequency design using a discrete-event simulation model. The software ProModel (Harrell and Tumay, 1992 ) is used to build the simulation platform. Buses and passengers are defined as entities. Bus stops and routes are modeled as locations and path networks, respectively. We define three sub-locations-wait area (WA), load area (LA), and unload area (UA) for each stop. A summary of "Variables" and entityspecified "Attributes" of the simulation model is presented in Table 3 . Table 3 : Attributes and variables of the simulation model Variables stopID: ID associated with each stop currStop: stopID of the stop where an entity currently stays maxCap: the maximum number of passengers a bus can hold Bus-associated Attributes currCap: current capacity of the bus Route: route ID of the bus Passenger-associated Attributes origin: stopID of the origin stop that the passenger arrives at dest: stopID of the destination stop that the passenger is moving to Buses arrive at each hub according to given frequency during different hours and move cyclically along the route paths. At each stop, buses unload passengers at UA, move to LA and load passengers from WA, and then leave for UA of the next stop along the route path. A passenger entity, who either gets off a bus or newly enters the system, will arrive at UA of the stop. The entity will decide whether to leave the system, to wait at WA, or to move to a near stop for transit. A passenger moves from WA to LA once a bus with sufficient capacity arrives and if the passenger needs to get on the bus.
The logic process described above is detailed in Algorithm 1. Passengers enter the bus system by arriving at their origin stops and move with the buses to reach their destination stops. The system must know at which stop a passenger can get off and which bus the passenger can take. In addition, a passenger may take a transit at some stop if the origin and destination stops are on different routes. In this case, the system needs to select a stop for the passenger to wait. These decision-making processes are achieved by solving for the shortest path between the passenger's origin stop and destination stop through algorithms such as Dijkstra's. Functions GetOn, GetOff, and NextStop shown in Algorithm 2 use the shortest path found by Dijkstra and the current stop to determine if the passenger needs to get on a bus, get off a bus, or select a stop to wait/exit the system, respectively. An external function block written in C++ is connected to the ProModel through a DLL application extension file for performing these decision-making functions.
We instantiate the new bus route design with the following settings. We test the six routes produced by Model (1) Algorithm 2 Functions used in the simulation to decide individual passengers' actions.
function Dijkstra(origin, dest): Input: route graph, stopID origin, and stopID dest Require: stop dest is accessible to stop origin Output: the shortest path between the given stops: origin and dest end function function GetOn(Passenger, currStop): Input: Passenger's attributes (origin and dest), stopID currStop that Passenger is currently at Require: currStop is on the shortest path; Passenger is at WA Output: Route that Passenger can take to move on the shortest path Tables 11-16 indicate the wait time at that stop in most cases. The determination of actual wait time follows the following rules. A bus waits 1 minute for most stops if it needs to load/unload passengers. It waits 2 minutes at a stop that has more than 10 passengers needing to board, or at the two hubs: Central Campus Transit Center (CCTC in central campus) and Pierpont Commons (in north campus). The baseline capacity of each bus is set as 40 passengers, being around 50% of the full capacity (i.e., 78 seats). We focus on bus operations during two peak periods: 8am-10am and 12pm-2pm. The average hourly total number of passengers is an adjustable input parameter with the default total being 2625 passengers arriving at all stops per hour. Passengers' arrival follows exponential distribution and the actual hourly demand fluctuates around the given value of the parameter. The percentage of Origin-Destination pair (O-D pair) is assumed to be fixed as shown in Table 5 and Table 6 for periods 8am-10am and 12pm-2pm, respectively. The O-D pair estimation during 8am-10am is based on the analysis of number of residents in different student dorms, their majors, daily activities, and course schedules in Fall 2019. The O-D pair settings during 12pm-2pm are the reverse of the settings in the morning (i.e., 8am-10am). Additionally, 12% randomly generated passengers will arrive at one of the 44 stops with equal probability (i.e., each stop has 0.27% random passengers arriving) and get off at all the other 43 stops randomly.
We assume that passengers can make transits at the same stops or at tuples of nearby stops.
Stops, through which at least two routes pass, enable passengers to perform same-stop transit with no transit time. Some pairs or triples of stops are considered to be within walking distance, so that passengers can make a transit among them with a transit time set as 0.5 minute. List of stops that allow same-stop transit is presented in Table 7 and list of stops that allow nearby-stop transit is presented in Table 8 . 
We test the redesigned bus system's performance given three scenarios that have different bus capacities and hourly demand estimations. We also conduct tests about the system's resilience when there exists a random bus broken-down on any of the six routes. We run 40 replicates for each scenario setting and the average results of the 40 replicates are analyzed in this section. We use the statistics "average utilization rate per stop-to-stop trip" U s and "percentage stop-to-stop trip with at least 75% utilization rate" S 0.75 to measure bus operational performance. A stop-to-stop trip is defined as the buses' travel for one stop, i.e., from the previous stop to the next stop along a route. When a bus leaves a stop, we compute its utilization rate as
and then calculate U s as the accumulated utilization rate divided by the total number of stops that buses have passed by. To obtain S 0.75 , we count the percentage of stop-to-stop trips that have utilization rate no less than 75%. We present the results of bus operational performance in Table 9, where three scenarios A, B, C are tested. Note that in scenario A, we allow to use 50% of the total number of seats as the desired capacity, while the social-distancing requirement is more strict and we only allow 25% bus capacity use in both scenarios B and C. In scenario C, we have to handle 2625 passengers per hour, being around 150% of the peak-hour travel demand in regular fall terms at the UM in the past (e.g., in Year 2019). For passenger-related measurements, we focus on the time each passenger spends in the system on average and how many transits they have to take. We also measure the average wait time of each passenger spent on each bus stop. If a passenger experiences any transit, the wait time is the sum of wait time at all the transit stops. Nearby-stop transit time (0.5 minute) is also included in the wait time. In Table 9 , we present results of the passenger-related performance of the system under different scenario settings during 8am-10am. (Period 12pm-2pm has similar results and performance, which are omitted for brevity.)
From Table 9 , the maximum utilization rate of buses on any of the six routes is around 50%-60%,
indicating that on average the number of passengers on each bus would be no more than 25. This satisfies the social-distancing requirement discussed in Zhang et al. (2020) for public transit systems.
In addition, no more than 54% stop-to-stop trips load more than 75% of its desired capacity (either 40 or 20 according to the first row in Table 9 .) Although increasing hourly demand and decreasing More than half passengers wait less than 5 minutes. Only 1 passenger in 5 needs at least one transit to complete a trip.
Next, we perform stress tests by randomly breaking down one bus (i.e., reducing the total number of buses by one) on a route in each scenario setting, while keeping other parameter settings the same. Using Scenario A (C = 40, D = 2625) as the baseline, we present statistics of passengerrelated performance measures in Table 10 . The performance results in Table 10 are similar to the ones in Table 9 and Figure 4 for all three scenarios, which implies that the system is insensitive to small changes of bus availability on any route and it is resilient to capacity reduction.
A hub-and-spoke transportation system relies on short, direct routes with fewer number of stops. 