{'fraction_non_alphanumeric': 0.028662346410238543, 'fraction_numerical': 0.014911852832555222, 'mean_word_length': 4.6297482837528605, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A public transit network optimization model for equitable access to social services

We present a flexible public transit network design model which optimizes a social access objective while guaranteeing that the system's costs and transit times remain within a preset margin of their current levels. The purpose of the model is to find a set of minor, immediate modifications to an existing bus network that can give more communities access to the chosen services while having a minimal impact on the current network's operator costs and user costs. Design decisions consist of reallocation of existing resources in order to adjust line frequencies and capacities. We present a hybrid tabu search/simulated annealing algorithm for the solution of this optimization-based model. As a case study we apply the model to the problem of improving equity of access to primary health care facilities in the Chicago metropolitan area. The results of the model suggest that it is possible to achieve better primary care access equity through reassignment of existing buses and implementation of express runs, while leaving overall service levels relatively unaffected.


Introduction
The efficient design and operation of a public transit network is one of the most important and challenging tasks facing any urban planner. As such, the urban transit network design problem (UTNDP) has been the subject of an immense amount of attention in transportation engineering research for many decades. These problems are generally modeled using an optimization-based framework. Although there is a large amount of variety in the types of design decisions, constraints, and objectives used to define the model, the most common objectives consist of optimizing some combination of operator costs (including fleet size, vehicle maintenance, travel distance, number of stops, profit, etc.) and user costs (including travel time, waiting time, walking distance, number of transfers, vehicle crowding, etc.).
To give a few recent examples of UTNDP studies, Bielli et al. 2002 [8] proposed a genetic algorithm for bus network design with a wide variety of performance indicators, including minimizing number of required vehicles, number of lines, number of user transfers, user waiting time, and high crowding indices, while maximizing number of served users and vehicle efficiency. Chakroborty and Dwivedi et al. 2002 [11] considered a route design problem with the objective of minimizing a weighted combination of line scores based on minimizing user travel times, maximizing fraction of users able to reach their destination in zero, one, or two transfers, and minimizing fraction of users unable to reach their destination in two or fewer transfers. Fan and Machemehl 2006 [21] used a genetic algorithm to minimize a weighted sum of user costs (consisting of total travel time on both direct and transfer paths), operator costs (consisting of vehicle operation costs), and total unsatisfied transit demand. Bagloee and Ceder 2011 [4] studied transit planning decisions with the objective of minimizing a passenger discomfort index consisting of a weighted sum of in-vehicle time, walking time, waiting time, and number of transfers. Cipriani et al. 2012 [19] considered a bus route design problem attempting to minimize a weighted combination of user costs (in-vehicle travel time, access time, waiting time, and transfer penalty) and operator costs (bus travel distance and travel time), and a penalty for unsatisfied transit demand. Cheung and Shalaby 2016 [12] considered rerouting decisions with the goal of minimizing network congestion. Buba and Lee 2018 [9] used a differential evolution algorithm to solve a route design and frequency setting problem with the objective of minimizing a weighted sum of total user travel time and unmet demand. Chu 2018 [18] considered a bus route design and timetabling problem with the objective of minimizing a weighted sum of vehicle travel time, user travel time, and unsatisfied demand.
A relatively small set of operator cost and user cost objectives related to vehicle operation cost, user travel time, and travel demand satisfaction is repeated across most UTNDP studies. This is noted in literature reviews on the topic, such as Guihaire and Hao 2008 [28] and Farahani et al. 2013 [22], the latter of which specifically points out that other important objectives like environmental costs have gone largely ignored by past studies. It is understandable that these standard operator and user cost objectives are commonly used, since they apply directly to the most immediate and important measures of system performance. Even in a study which includes alternate optimization criteria these objectives are still extremely important to take into consideration. However, such a narrow set of design goals can lead to different communities experiencing extremely disparate levels of access to important public services.
Recently there has been a large amount of interest in studying such questions of differential access to services over the boundaries of a city. Radke and Mu 2000 [41] developed a method for studying household access to social services. Pearce et al. 2006 [40] studied community access levels to a variety of resources thought to impact community health, including recreation, shopping, education, and health care. Farber et al. 2014 [23] and McDermot et al. 2017 [36] both studied the phenomenon of food deserts (areas which lack access to healthy food options) and both found that differential access to food options is largely attributable to differential public transit access as opposed to the geographic locations of facilities. Finally, a large number of recent studies have attempted to measure differential access to primary health care facilities.
Most existing research in this area focuses primarily on developing a health care accessibility metric and using that metric to identify areas with particularly poor access. For example, Delamater 2013 [20] compared various catchment area-based accessibility metrics for measuring inpatient facility access in the state of Michigan, while McGrail and Humphreys 2014 [38] applied a catchment area-based accessibility metric to measure spatial accessibility of general practitioners throughout Australia. The topic of planning specific actions in order to alleviate these health care access issues has been less well studied. As an example, Gu et al. 2010 [27] considered an optimization-based facility location problem with the goal of placing new primary care facilities to maximize a global health care accessibility metric.
The facility location modeling paradigm is a useful tool for maximizing the public good provided by newly-opened facilities, but creating new facilities is not always an option. The placement of privately-owned facilities is governed by their owners' economic considerations and not necessarily by the needs of the community, and even public facilities are major investments that can generally only be placed in very limited number and in a very limited set of locations. For these reasons we instead turn our attention to an alternative to the facility location modeling paradigm for addressing the accessibility problem.
Transit Network Design Approach to Accessibility Maximization: The primary focus of this study is to consider a problem which is in some ways the inverse of the facility location problem. Rather than changing facility locations in a fixed transit network, we consider changing a transit network containing fixed facility locations. This can improve accessibility levels in the same way that placing facilities can, since reducing travel times makes effective catchment areas larger and brings fast access to more communities. Public transit planning can also be implemented much more immediately since tactical decisions like bus frequency setting can be enacted without the need for permanent changes to the infrastructure.
The main goal of this paper is to present an optimization model-based approach to the public transit design problem that attempts to achieve a social access objective while adhering to the following set of general modeling principles: 1. The model should be flexible, allowing the planner to choose any desired social access objective and any assumptions regarding user behavior.
2. The model must produce solutions that remain at or near the system's current cost and performance level while attempting to optimize the social access objective.
3. The design decisions should consist of measures which are low-cost, immediate, and easilyimplemented. 4. We will assume that travel related to our social access goal makes up a relatively small proportion of the day-to-day public transit travel volume. In particular we will assume that the capacity of the facilities in question is a much more significant limiting factor than the public transit service capacity.
Principles 2 and 3, in particular, are important for the model to be seriously considered by an urban planner. The main goal of a transit authority is to minimize the operator and user costs resulting from the large volume of day-to-day transit network usage. It would be unreasonable to implement sweeping system changes in service of a relatively small proportion of the overall travel volume, but low-cost, minor, tactical changes that have a minimal effect on most day-to-day traffic might have a better chance of being implemented. As a case study for testing the model formulation presented in this paper, we consider the problem of improving equity of primary health care access in Chicago, Illinois.
Chicago Metropolitan Area Case Study: The study of regional primary access for Chicago residents has been a subject of investigation for the Chicago Department of Public Health. Between 2005 and 2008 a sequence of reports were released to detail their findings [5,25,44]. Specifically these reports focused on the so-called safety net providers, which consist of a combination of community health centers and free clinics that focus on serving low-income patients.
A number of factors related to patient access to these facilities were studied, including socioeconomic status, ethnicity, gender, age, and geographic location. All reports identified geographic location as one of the biggest contributing factors to access. A large part of this geospatial disparity is due to the locations of the facilities relative to the patient populations. Specifically, Chicago's West Side contains a disproportionately large number of safety net facilities relative to its population, while the North and South Sides contain disproportionately small numbers. In both of these regions estimated population need exceeds estimated facility capacity (as determined by number of available physicians).
A 2003 study by Luo and Wang [33] applied various accessibility metrics to measure primary care accessibility throughout the Chicago region. They found that the distribution of accessibility levels can vary widely across even relatively small geographic areas of the city due to particular communities having an anomalously great or small amount of access relative to the surrounding areas. Specifically, their study identified three areas within the city with significantly better access than surrounding regions (The Loop, Lincolnwood-Skokie, and Elmhurst-Oak Brook) and two with significantly worse access (South Side and Midway), which agrees with the Chicago Department of Public Health findings. It should be noted that the South Side, in particular, is home to a disproportionately large number of citizens near or below the poverty level, many of whom lack health insurance and personal transportation. Figure 1 shows the distribution of primary care community health clinics throughout the city [13]. It is clear from the above studies that the city of Chicago is in need of more equitable access to primary care facilities. We will apply our modeling framework to examine changes that could potentially be made by the Chicago Transit Authority (CTA) in an effort to address this issue. Our goal is to find relatively simple, low-cost changes that could be put into practice with minimal impact on current service levels. Specifically, we consider altering the CTA bus network by reassigning vehicles between the existing bus lines, increasing the number of vehicles servicing certain lines while decreasing the number servicing others. We also consider the possibility of introducing express runs of existing routes.
Structure and Contributions of this Paper: In Section 2 we present the formulation of the optimization-based model, developed to fulfil the four modeling principles listed above, which constitutes the main contribution of this paper. We begin in Section 2.1 with its overall structure and then discuss its major components, including its underlying network representation in Section 2.2. In Section 2.3 we describe the model's objective function, which consists of a dynamic evaluation of the accessibility levels of the least-advantaged communities based on the current solution vector. In Section 2.4 we describe a set of model constraints which allow the planner to guarantee that operator and user costs remain within a specified margin of their initial levels. The structure of this model is unique in that its constraint functions are computationally expensive to evaluate while its objective function is not, and so in Section 3 we present a hybrid tabu search/simulated annealing solution algorithm designed to minimize the number of computationally expensive evaluations needed.
The remainder of the paper is devoted to the computational trials used to test the performance of this model. In Section 4 we explain the computational trials used as test cases, including the main CTA network case study in Section 4.1 and a set of artificial network tests in Section 4.2. The results of these trials are presented in Section 5 and discussed in Section 6. Appendix A includes a more in-depth description of the solution algorithm from Section 3, Appendices B-C give detailed descriptions of the data and networks used for the computational trials explained in Section 4, and Appendix D includes the full data tables for these trials.

Model Formulation
In this section we describe our proposed model for addressing the design goals discussed above. We first show the overall structure of the model followed by detailed explanations of the individual components. Due to our choice of case study we will make specific reference to primary care access and a bus transit network, but these can easily be replaced with any other sort of objective and transit network that corresponds to the above assumptions.

Overall Model Structure
In this section we give a high-level description of the overall structure of our proposed model. In the interest of keeping with our goal of considering only low-cost, immediate, easily-implemented tactical changes to the transit network, our design decisions consist of determining the daily average number of buses to assign to each line. A line's fleet size implicitly determines its average service frequency, and so this can be considered as a type of frequency-setting problem [22].
We begin by defining some variables that will be used throughout this paper. Let L be the set of transit lines under consideration. For each line l ∈ L, let y l ∈ Z be the number of vehicles that service line l, and let y be the vector of such variables arranged in a canonical order, which constitutes the decision vector for our design problem. Let y * l be the number of vehicles initially assigned to line l, and y * be the initial fleet vector. Without loss of generality we assume that we are given constant upper and lower fleet size bounds y max l and y min l , respectively, for each line l (which can be made infinite if such bounds are not needed, or which can both be set to y * l to force y l to remain constant). We assume that each line utilizes only a single type of vehicle, since a line that uses multiple vehicle types (e.g. a bus line serviced by both rigid-bodied buses and articulated buses) can be represented by multiple parallel lines in L. Let Z be the set of all vehicle types utilized throughout the transit network. For each vehicle type z ∈ Z let L z ⊆ L be the set of lines serviced by that vehicle, and for each line l ∈ L let z l ∈ Z be the vehicle type used by line l.
Let x be the flow vector for user traffic through the transit network. User flows are defined implicitly as the output of a transit assignment model which applies assumptions about user behavior to determine how the users, in aggregate, move through the system. The representation of the flow network and the transit assignment model will be discussed in Section 2.2. Because user decisions depend on the structure of the network as determined by the design decisions y, we will denote the transit assignment function as TransitAssignment(y). Let x * be the user flow vector associated with the initial fleet vector.
In a typical public transit design model like those discussed above the objective is to minimize some combination of operator cost functions and user cost functions subject to design feasibility and budgetary constraints. We denote these two functions as OperatorCost(y, x) and UserCost(x), respectively. It may be assumed that operator costs depend directly on both our design decision vector y and the user flow vector x, but the user costs generally depend directly only on the user flows (though both depend implicitly on the design decisions). Possible choices for these functions will be discussed in Section 2.4, but in general operator costs represent costs associated with running the transit network while user costs represent how difficult it is for users to travel through the network.
Finally, let Access(y) be a function describing the social access objective. Possible choices for this objective will be discussed in Section 2.3, but in general it should represent how easy it is for users to access some service within the network, weighted to focus on the most disadvantaged users. It depends only on the network design y and not the user flows x because, in line with modeling principle 4, the user flow volume is assumed to be small and infrequent enough that it does not compete for transit capacity in the same way that day-to-day flow volume does.
The three objectives OperatorCost(y, x), UserCost(x), and Access(y) cannot, in general, all be optimized simultaneously, which makes this a multiobjective program. We take the approach of optimizing only the access objective while including the operator and user costs as constraints. This allows our model to explicitly attempt to make improvements to the social access objective while still guaranteeing a specified service level. Our proposed model takes the form of the following nonlinear program, which we will refer to as the social access maximization problem (SAMP): To explain, the overall objective (1) is to maximize the access objective. Constraints (2) and (3) give explicit upper bounds of B operator and B user to the operator and user costs, respectively. These bounds allow the transit planner to specify precisely how much of an increase is acceptable in exchange for maximizing the social access objective, which can be chosen based on their requirements. Constraints (4) implicitly define the user flow vector x as a result of the decision vector's effect on the transit assignment function. Finally constraints (5) ensure that no new vehicles of any type are added, constraints (6) enforce fleet size bounds for each line, and constraints (7) enforce integrality of fleet sizes. Note that, in general, transit assignment functions are computationally expensive to evaluate, and so determining the user flows and thus the feasibility of constraints (2) and (3) for a given solution vector y is expected to be computationally expensive. As a result it is far faster to evaluate the objective of a solution to the SAMP than its feasibility, and this unusual structure dictates the form of the solution algorithm discussed in Section 3.
As the SAMP is meant to represent a general modeling framework for handling social accessoriented public transit design problems, the exact definitions of the functions Access, OperatorCost, UserCost, and TransitAssignment should be chosen based on planner's specific goals. The remainder of this section will be spent describing some possible choices for these functions and explain the modeling decisions made for our case study of Chicago primary care access.

Transit Assignment Model and Transit Network Representation
The mathematical representation of the transit system takes the form of a network with nodes to represent locations (like transit stops, population centers, and primary care facilities) and arcs to represent traversable routes between those locations (like bus trips, train trips, and sidewalk travel), but the specific structure of this network depends on the choice of transit assignment model, which may require the inclusion of auxiliary nodes and arcs to represent different parts of the user decision-making process. In this section we describe our choice of transit assignment model and the network representation that it requires.
Assignment Model: The transit assignment model, TransitAssignment(y), applies a set of basic assumptions regarding user behavior to simulate how public transit users route themselves through the transit network. The flow produced by such a model is called the user-optimal flow, which in general leads to worse objective values than would be possible if x were directly controllable. There is an extremely wide variety of transit assignment models used throughout transportation design literature, differing primarily in their assumptions regarding user behavior and how they account for user interactions in the presence of congestion [24]. For the purposes of our study we have chosen to use a nonlinear model by Spiess and Florian [48]. This model works by using a conical congestion function to increase the costs of arcs as they become more crowded, serving to discourage flows from exceeding line capacities. The output is the user-optimal flow vector which corresponds to all users choosing a strategy for for which their expected travel time cannot be improved with a unilateral change. This particular model was chosen because it is well-established and commonly used in transit forecasting systems (such as EMME/2) as well as being relatively simple, which is desirable since the solution algorithm presented in Section 3 potentially requires many evaluations of the operator cost and user cost functions during each iteration.
Network Structure: Our network representation of the public transit system was chosen to meet the requirements of this assignment model. We use directed graph G = (V, E) to represent the public transit network. The node set V describes locations of interest within the city, including origins, destinations, and public transit stops. The arc set E represents links between these locations via different modes of transportation, including walking, bus travel, and train travel. We consider only public transit and foot travel, not automobile traffic.
We partition the node set into subsets. Let V stop ⊂ V be the set of stop nodes, which represent the physical locations of public transit stops. Let V board ⊂ V be the set of boarding nodes, which are used to distinguish between different transit lines at each stop. For each stop node i ∈ V stop there is a boarding node j for all lines l ∈ L that utilize stop i. Let V + ⊂ V be the set of origin nodes, which represent locations the locations where public transit trips begin. Finally, let V − ⊂ V be the set of destination nodes, which represent locations the locations where public transit trips end. Note that these origin-destination (OD) pairs represent the high-volume, high-frequency day-to-day trips that the assignment model considers, not the social access-related objectives.
We similarly partition the arc set. Let E line ⊂ E be the set of line arcs, which represent invehicle travel for a transit line. Let E board ⊂ E be the set of boarding arcs, which represent the decision to board a particular line at a particular stop. Let E alight ⊂ E be the set of alighting arcs, which represent the decision to alight a particular line at a particular stop. Finally, let E walk ⊂ E be the set of walking arcs, which represent travel along sidewalks or other walking paths.
For each line l ∈ L, let G l = (V l , E l ) be the subgraph of G containing only nodes and arcs associated with that line. Let any of the above partitions with a superscript of l represent the subset of that partition corresponding to line l (so for example V l stop is the set of all stop nodes traversed by line l). Note that origin nodes, destination nodes, and walking arcs have no associated lines. For every stop node i ∈ V l stop traversed by line l, there is a corresponding boarding node i l ∈ V l board . These two nodes are joined by a boarding arc (i, i l ) ∈ E l board and an alighting arc (i l , i) ∈ E l alight . If stop node j ∈ V l stop is the next stop on line l after stop i, and if its corresponding boarding node is j l ∈ V l board , then there is a line arc (i l , j l ) ∈ E l line joining the consecutive boarding nodes. Walking arcs are not line-specific and may join any pairs of stop nodes, origin nodes, and destination nodes. Figure 2 shows an example of this network structure.
In order to incorporate the social access objective into this network, we also define two more types of node. Let V + ⊂ V be the set of community nodes, which represent locations the communities of concern in the social access objective. Let V − ⊂ V be the set of facility nodes, which represent locations of the service of concern in the social access objective. Whereas the objective of a typical public transit study is to facilitate travel between origin nodes V + and destination nodes V − , the social access objective is to facilitate access of nodes V + to nodes in V − . Like origin and destination nodes, walking arcs connect community and facility nodes to the rest of the network.
Frequency-Based Attributes: To relate the structure of network G to the SAMP we use a frequency based approach, a common choice for many tactical UTNDP problems [22,28] and transit assignment models [10,24,30,35,45,46], including our chosen model. Under this paradigm all traffic demands and flows are averaged over a specific time period during which it is assumed that travel demand and traffic levels remain constant. The frequency-based approach does not take vehicle timetables into consideration, instead considering only the expected arrival frequency of each line at each stop which in turn determines a user's expected waiting time. We will use τ to represent the daily time horizon for use in the frequency-based model.
For each OD pair (s, t) ∈ V + × V − , let D st be the total day-to-day transit demand over the specified time period. The |V + | × |V − | matrix for which element (s, t) is D st is known as the OD matrix. For each arc ij ∈ E, let x ij be the total day-to-day flow through that link, and let x be the vector of such flow values (ordered in a canonical way).
The decision variables y l of the SAMP represent the average number of vehicles servicing each line l ∈ L on a day-to-day basis. We assume that these tactical decisions are minor enough not to disturb any overall traffic patterns, and so each line also has a constant expected time of τ l (in minutes) to complete a full circuit (including layovers). These decision variables affect the output of the transit assignment model in two ways: First, the number of vehicles servicing a line dictates that line's average service frequency f l (in units of arrivals/minute) which, under the assumption of exponentially distributed arrival times, is determined by the inverse proportionality f l = y l /τ l [48]. Second, the number of vehicles servicing a line dictates that line's total capacity u l (in units of total passengers over the given time horizon), which for a line with a vehicle seating capacity of b l is determined by u l = τ b l f l . Within the assignment model line frequency is used to determine expected waiting time while line capacity is used to determine penalty costs for overcrowding, both of which contribute to users' path selections. It should be noted that any line which is assigned no vehicles, y l = 0, has an effective capacity of zero and an effectively infinite waiting time, causing it to be chosen by no users.
For the purposes of the social access objective we assume that the volume of traffic is so small that it makes up an insignificant fraction of the day-to-day travel volume, and so it is ignored within the assignment model. For this reason it can be assumed not to depend on the capacities of the lines, but it may still depend on the frequencies, since they determine expected waiting times and thus expected trip length. In addition, the introduction of new lines such as the express runs to be discussed in Section 4.1 can drastically affect both day-to-day traffic and social access traffic by making new paths available.

Objective Function
There are many reasonable choices for the social access objective function Access(y), but in general it is meant to describe the ease of traveling from a community to a facility using the network resulting from design decisions y, and to do so in a way that places focus on equity of access. Our proposed objective is based on maximizing primary care accessibility metrics for the set of communities that currently has the lowest level of access. As referenced in Section 1, a large amount of research has gone into the development of accessibility metrics, particularly with regard to primary care access. For example, the two-step floating catchment area (2SFCA) metric was first described in Luo and Wang 2003 [34] and has been used in a wide variety of studies since then [32,33,[37][38][39]49].
Let A i (y) be a primary care accessibility metric of community i ∈ V + . Our specific choice of A i (y) will be described below, but in general accessibility metrics depend on travel times within the network, and so are implicitly functions of the decision vector y. Then we define our objective function as the sum of the K minimum accessibility metrics, where K is a model parameter. Note that as the solution algorithm proceeds and the solution vector y is changed, the accessibility metrics A i (y) and their relative rankings may also change. This objective refers to the K current least accessibility metrics (for the current solution y).
This objective generalizes the ideas of maximizing only the current minimum accessibility metric (which is the special case of K = 1), and of maximizing the total (or equivalently the average) of all accessibility metrics (which is the special case of K = | V + |). Attempting to optimize for only the current least-accessible community leads to potential problems if that community cannot be significantly helped through changes in the transit network, since it results in the solution process focusing all of its decisions on attempting to make minor improvements in one location while ignoring decisions that might help far more communities. On the other hand, attempting to optimize the total accessibility level of all communities could result in decisions which improve the accessibility levels of areas that already have high accessibility while worsening the accessibility levels of disadvantaged areas, as long as it improves the overall total. The computational trials for our Chicago case study, as well as for artificial test networks, displayed both of these problems for extreme values of K (see Section 5).
Our specific choice of accessibility metric A i (y) is a gravity-based metric. Gravity-based metrics take into account how many facilities are near a community, however rather than a binary count of whether or not the travel time falls below a cutoff, each facility is weighed by a negative power of the pairwise travel time, causing the supposed access level to a facility to decay over distance like a gravitational force. This was chosen for use in the SAMP objective because, unlike the catchment area-based metrics, its continuous dependence on travel times generally gives it a nonzero gradient, which helps to give direction to the the solution algorithm's local search (see Section 3).
The specific metric used in our study is a competition-based model developed in Weibull 1976 [50] for describing differential access to employment opportunities. It has also been used in the study of health care access, for example in Joseph and Bantock 1982 [29] and in Luo and Wang 2003 [33]. It is defined as where S j is the quality of facility j ∈ V − (usually taken as its patient capacity), d ij (y) is the travel time from i to j (calculated as graph distance in G, which depends on the network design decisions y), β > 0 is a gravitational decay model parameter, and F j (y) is a competition-related facility metric defined by where P k is the population of community k seeking service at facility j. F j (y) can be interpreted as a measure of how crowded facility j is, and is greater if the facility is close to many populous communities. The accessibility metric A i (y) rewards a community for being close to many facilities with high quality and low overcrowding. The metric, itself, has no direct interpretation in isolation and is meant for use in comparing communities across space or time, as will be shown in Section 5.

Operator Costs and User Costs
The operator and user cost functions represent typical objectives from standard UTNDP studies. Operator costs generally represent monetary costs incurred by the public transit authority, and may include vehicle operation costs, vehicle acquisition costs, and fare revenue. User costs generally represent time costs incurred by users of the transit system, and may include travel time (which may be weighted by mode), waiting time, number of transfers, and measures of discomfort.
Because this study involves only the reassignment of existing vehicles within the public transit network and not the purchase of new vehicles, and because operator costs consist largely of budgetary concerns, the operator cost function is unlikely to change as a result of our decisions. For this reason we will exclude the operator cost constraint (2) from our model for the remainder of this paper, however it would still be appropriate to include for a different set of design decisions (for example, including vehicle purchase or line construction). Our user cost function takes the form of the weighted sum To explain, although the form of the user cost function mentioned in Section 2.1 included only a flow vector argument x, our chosen assignment model produces both a flow vector and a total waiting time scalar ω which represents the total time spent by all users waiting to board [48]. The above user cost function is then a weighted sum of three types of travel time: total in-vehicle travel time, total walking travel time, and total waiting time, with weights of θ 1 , θ 2 , and θ 3 , respectively, which can be used to penalize particular types of travel time more harshly than others. The user cost bound constraint (3) in the SAMP prohibits UserCost(x, ω) from increasing beyond a constant upper bound B user . We use a constant bound of where ≥ 0 is a model parameter which specifies the allowed margin of increase in user costs, x * is the initial user flow vector, and ω * is the initial waiting time scalar. This bound corresponds to the choice of B user = (1 + ) UserCost(x * , ω * ) and ensures that the user cost cannot undergo a relative increase of more than above its initial level. A larger value of results in a less restrictive bound and thus a larger feasible set and potentially a better optimal value, but at the cost of potentially increasing the overall user costs. In this way allows for different choices in the tradeoffs between improving social access and maintaining day-to-day service levels.

Hybrid Tabu Search/Simulated Annealing Solution Algorithm
In this section we give a brief description of our proposed solution method for the SAMP model defined in Section 2. As mentioned above, we wish to allow for a variety of choices in the definitions of Access(y), OperatorCost(y, x), UserCost(x), and TransitAssignment(y), and so we will treat all of these functions as black boxes, assuming that each can be solved through some subroutine and making use only of its outputs for a given set of inputs.
Versions of the UTNDP are well-known as computationally difficult problems for a variety of reasons, including their large number of integer variables and the bilevel structure which arises from the need to consider user behavior as an implicit function of the design decisions. As a result they are in general nonlinear and nonconvex, and most studies on the subject resort to a heuristic or metaheuristic solution method to obtain high-quality solutions within a reasonable time frame [4,8,11,19,22]. Transit assignment models are in general computationally expensive to solve [28], and so it is reasonable to assume that TransitAssignment(y) and anything that relies on its outputs (including OperatorCost(y, x) and UserCost(x)) are significantly more expensive to evaluate than the other major components of the model. This gives the SAMP the unusual structural property of having constraints (2)-(3) that are significantly more expensive to evaluate than the objective (1). Because of this, our proposed metaheuristic approach is designed to minimize the number of required constraint function evaluations.
We propose a hybrid tabu search/simulated annealing solution algorithm for the SAMP, due to the fact that there is a natural choice of initial feasible solution (the initial fleet vector, y * ) as well as a natural definition of local moves, consisting of adding, dropping, and swapping individual vehicles between compatible lines. Tabu search (TS) and simulated annealing (SA) are widelyused combinatorial optimization metaheuristics. They are both based on local search, wherein each solution has a defined neighborhood of solutions which differ by only one local move, and the solution algorithm consists of iteratively conducting neighborhood searches and making local moves to improve the objective value. Pure local search carries the risk of becoming trapped at a local optimum while much better solutions exist, and so TS and SA both include rules for occasionally making suboptimal moves in order to explore more of the feasible set and arrive at a better solution.
SA accomplishes this by introducing randomness [7, pp. 512-514], having a chance to make suboptimal moves according to a parameter T called temperature which decreases as the algorithm runs. As T decreases suboptimal moves become less likely, and so late in the run time it is more likely to remain near local optima. On the other hand, TS accomplishes this by using memory structures to avoid repeatedly cycling through the same solutions [26], using a combination of short-term memory (STM) to avoid reversing a move too soon after it is made, and long-term memory (LTM) to remember profitable past solutions that can be returned to if the algorithm becomes trapped. The STM generally consists of a list of tabu rules which exist for a set tenure before being dropped from the tabu list. Hybrids of the two algorithms have been shown in a number of recent combinatorial optimization studies to perform better than pure TS or SA alone [1,2,31,47].
The general structure of our approach is based primarily on that of Alhroob et al. 2014 [1]. It follows the TS structure discussed above, but uses an SA acceptance probability when considering whether to make a suboptimal move. Tabu rules stored in STM prevent undoing recent additions or subtractions from each fleet, and attractive solutions stored in LTM consist of a combination of the second-best solutions from neighborhood searches which were not chosen, and of suboptimal moves that were previously denied by the SA criterion. Tabu tenures are adjusted adaptively to increase when more diversification is required and to decrease when more intensification is required, measured by using two separate nonimprovement move counters: An inner counter increments along with the tabu tenures whenever no improving non-tabu move exists, resets whenever a suboptimal move is made by passing the SA criterion, and backtracks to a random attractive solution from the LTM when it reaches a set maximum. An outer counter increments alongside the inner counter as well as whenever a an attractive solution is backtracked to, resets whenever an improving local move is made, and prompts the tabu tenures to reset when it reaches a set maximum.
In order to take advantage of the fact that the social access objective can be evaluated relatively quickly while the operator and user cost constraints cannot, our algorithm uses a two-pass approach: The first pass consists of considering every possible ADD and DROP move that satisfies the trivial design constraints (5)-(7), evaluating its objective value, considering SWAP moves made up of the most profitable ADD and DROP moves combined, and evaluating their objective values. The second pass then consists of evaluating the much more computationally expensive user flow-dependent feasibility constraints (2)-(3) of the neighbors in descending order of objective value, halting as soon as two feasible solutions are found, which due to the first pass are guaranteed to be the two local moves with the maximum objective value. At all times tabu rules are enforced to limit the neighborhood, unless a new global optimum is found.
The following is a brief summary of our proposed solution algorithm. The full pseudocode can be found in Appendix A.
1. Initialize solution vector y (0) as the current real-world fleet size vector y * , and set it as the best known solution. Initialize TS and SA memory structures, including the SA temperature T , tabu tenures t, empty STM and LTM sets, and zero nonimprovement counters. d. Among the finalized ADD, DROP, and SWAP move lists, select the two solutions with the greatest objective values. Letỹ (k,1) be the best andỹ (k,2) be the second-best.
e. Ifỹ (k,1) is an improvement over the current solution y (k) , then move to the improving neighbor by defining y (k+1) asỹ (k,1) , reset the inner nonimprovement counter and tabu tenures, and add tabu rules to STM to prevent undoing any ADD or DROP that was just conducted.
f. Otherwise, increment both nonimprovement counters and then apply do the following: i. (SA Acceptance Criterion) With probability exp{−[Access(y (k) ) − Access(ỹ (k,1) )]/T }, move to the best neighbor by defining y (k+1) asỹ (k,1) , increment the tabu tenures, reset the inner nonimprovement counter, define tabu rules in STM to prevent undoing any ADD or DROP that was just conducted, and save the second-best neighborỹ (k,2) in LTM. ii. Otherwise, save the best neighborỹ (k,1) in LTM.
g. (Evaluate Nonimprovement Counters) If the inner nonimprovement counter has met its maximum, then reset it, increment the outer nonimprovement counter, and move to a random solution from LTM. If the outer nonimprovement counter has met its maximum, reset the tabu tenures.
h. (Upkeep) Drop moves from STM whose tabu tenures have elapsed. Apply a cooling schedule to the temperature T .
3. Return the best known solution and its objective value.
The source code for our C++ implementation can be viewed online [42]. The algorithm includes a large number of tuning parameters which were adjusted using trial-and-error during the computational trials described below. For the purposes of the Chicago case study, the SA criterion used an initial temperature of 10 16 and a cooling schedule consisting of multiplying the temperature by a factor of 0.999. We allowed a maximum of 40 solutions in LTM, after which attractive solutions were randomly dropped. We evaluated 50 ADD and DROP moves during each first pass and kept 2 of each during each second pass. The initial tabu tenures were 6.0 and increased by a factor of 1.15 whenever incremented. Nonimprovement counter cutoffs of 20 (inner) and 10 (outer) were used.

Computational Trials
In this section we describe the computational trials used for evaluating the performance and results of the SAMP. In addition to the Chicago primary care access case study described in Section 1, we also consider a small-scale artificial network in order to generate additional computational results for use in sensitivity analysis. The results of these tests are given in Section 5.

CTA Design Problem
In this section we briefly explain the setup of the problem and some of the tests that were run using the model. An in-depth explanation of the data and network setup is given in Appendix B. The main design problem of interest for this study is that of making minor alterations to the CTA network in order to improve equity of access to primary care facilities. To that end an underlying network G = (V, E) was constructed to represent the CTA public transit network. The set of stop nodes V stop corresponded to the set of public transit stops within the city and were also used as the origins and destinations for day-to-day transit volume. The OD matrix was estimated from stoplevel boarding data. To reduce the size of the network we clustered the more than 11, 000 defined CTA stops by identifying nearby locations with each other, resulting in a set of |V stop | = 997 stops.
Community nodes V + were taken as population centers within the city, defined using 2010 census data clustered into population-weighted centroids for each of the 77 Chicago community areas, resulting in | V + | = 77 community nodes. Facility nodes V − were taken as primary care facilities, obtained from a list of primary care community health clinics provided by the City of Chicago, resulting in | V − | = 119 facility nodes (see Figure 1). Each facility j ∈ V − was given a constant quality metric of S j ≡ 1 in (9). Walking arcs were generated to connect locations within 0.75 miles of each other, with travel times calculated under the assumption of 4 ft/sec walking speed. The final network G included a total of 3852 nodes and 17,522 arcs.
The public transit lines were defined using data available from CTA General Transit Feed Specification (GTFS) files. The network includes 126 bus lines and 8 train lines. Schedule data were used to estimate the current number of vehicles y * l as well as the average circuit time τ l for each line l ∈ L. The unit time horizon τ was taken as the full 24-hour (1440 minute) period. For simplicity we assumed that all bus lines used a single type of bus (the New Flyer D40LF, with b l = 39 seats per bus) and that all train lines used a single type of train (a 6 car train of Bombardier Transportation 5000 series cars, with b l = 228 seats per train). Throughout our computational trials we considered only changes to the bus lines, with the train lines being kept constant. The bus fleet lower bounds y min l were taken as the minimum number of vehicles required to achieve a line frequency f l of at least one arrival per 30 minutes. No upper bounds were used.
In order to evaluate whether the introduction of express lines can positively affect primary care access levels, we also defined potential express runs of existing bus lines. An express run follows the same route as the original line but stops at only a subset of the original stops, reducing its circuit completion time τ l and thus increasing its effective frequency f l as well as some of its pairwise travel times. While it is possible to build express line design decisions into the SAMP by introducing binary decisions for which stops to include and exclude, we instead elected to generate a small set of heuristically designed express lines (see Appendix B.5). These lines were built into G alongside the current CTA lines but defined to have initial fleet sizes of y * l ≡ 0, making them unavailable for use until adding at least one vehicle. The bus reassignment decisions y then allow the planner to transfer buses from existing lines to the new express lines if doing so would improve the access objective. A total of 65 express runs of existing lines were generated, each of which included a subset of approximately 20% of the original line's stops chosen to connect areas of high access to areas of low access.
The solution algorithm discussed in Section 3 was used to solve the instances of the SAMP discussed above, running for 500 search iterations followed by an exhaustive neighborhood search of the best known solution to ensure local optimality. Throughout all tests the access function (8) was evaluated using a gravitational decay parameter of β = 1.0 and a community count of K = 8 (which is approximately 10.39% of the 77 community areas), meaning that the objective was to maximize the 8 lowest accessibility metrics. The user cost function (11) used equal mode weights of θ 1 = θ 2 = θ 3 = 1 and a maximum relative increase of = 0.01, with constraint (12) preventing total user travel time from increasing by more than 1% above initial levels. The results of these trials are discussed in Section 5.1. A few modifications of the main trial were also run, including: a version which ignored the user cost constraint (3), a version for which the only allowed moves were to transfer vehicles between lines and their express runs, and a version for which the express lines were excluded. The results of these supplemental trials are discussed in Section 5.2.

Artificial Network Tests
Part of our evaluation of the SAMP model involved conducting sensitivity analysis to determine whether its results depend strongly on the particular choice of modeling parameters β, K, and .
The full CTA problem instance is too computationally expensive to solve repeatedly under varying circumstances, and so in order to generate more data we defined a small artificial network whose SAMP instance could be solved much more quickly.
A full description of the artificial test network is given in Appendix C. It represents a 0.9 × 1.5 mile city arranged in a square grid of blocks on a 7 × 11 street lattice, and contains a total of 527 nodes (including 204 stops, 30 communities, and 12 facilities) and 2206 arcs. The 18 transit lines each runs along one of the North/South or East/West roads. A total of 200 buses are in service with each route initially having between 6 and 14 vehicles and constant fleet size bounds of y min l = 1 and y max l = 20. Transit stops are located at each road intersection with some interspersed between intersections, and are used as OD pairs for day-to-day transit. Communities are distributed somewhat uniformly throughout the city, each having a population between 1000 and 3000. Facilities are distributed unevenly throughout the city with a constant quality metric of S j ≡ 1.
Trials were set up to test various values of the model parameters β, K, and . For β we tested β ∈ {0.5, 1.0, 1.5, 2.0} while holding K = 6 and = 0.01. For K we tested K ∈ {1, 6, 12, 18, 24, 30} while holding β = 1.0 and = 0.01. For we tested ∈ {0.00, 0.01, 0.05, 0.10, ∞} (with = ∞ meaning that the user cost constraints were completely ignored) while holding β = 1.0 and K = 6. For each parameter combination we solved the SAMP using 500 solution algorithm iterations followed by a neighborhood search of the best known solution. In order to reduce potential noise in the solutions from the random aspects of the search algorithm each problem instance was solved 10 times and the solution with the best objective value was kept. The results of these trials are discussed in Section 5.3.

Computational Results
In this section we summarize the results of the computational trials discussed above. Complete data tables for all computational trials can be found in Appendix D.

CTA Network Design Results
The SAMP solution algorithm for the express design problem described in Section 4.1 was run for 500 iterations, after which the primary care accessibility metrics A i (y) were evaluated on the final network and compared to their original values. Table 1 shows statistical results for the accessibility metrics of the 77 Chicago community areas. Of the 77 communities, approximately half (39 of 77) of the accessibility metrics improved while half (38 of 77) worsened. Among the 39 improved metrics the mean relative change was +30.28551%, while among the 38 worsened metrics the mean relative change was −12.27414%. The greatest individual increase among any community area's metric occurred for Mount Greenwood, which increased from 2.08703 · 10 −5 to 4.84086 · 10 −5 for a relative difference of +131.94930%. The greatest individual decrease occurred for Near North Side, which decreased from 5.30391 · 10 −5 to 4.18952 · 10 −5 for a relative difference of −21.01062%. Figure 4 shows maps comparing the initial accessibility metric distribution to the final distribution, while Figure 5 shows the relative difference in accessibility metrics before and after the trial. See Tables 3-4 in Appendix D.1 for a complete list of metric changes over all community areas. Figure 3 shows the objective value of the express route design over the 500 search iterations. The objective increased from an initial value of 1.87302 · 10 −4 to a final value of 2.70226 · 10 −4 , for an increase of 44.2%. The first 166 iterations consisted entirely of exchanging vehicles between non-express routes. The first express route was added during iteration 167, and the second was added during iteration 186, after which many more vehicles were added to express routes. This

CTA Network Supplemental Results
The results of the supplemental tests discussed at the end of Section 4.1 are summarized below.
Express Line Exchanges: A more restrictive version of the main express line design problem was run for which vehicles could not be exchanged freely between lines, but rather only exchanged between a line and its express run. The solution algorithm was run for 200 search iterations without finding a single solution that improved on the initial objective value. It was found that every possible SWAP move from the initial solution worsened the objective value, making the initial solution (with no exchanges) locally optimal.  Unconstrained Trial: A relaxed version of the main express line design problem was run in which the user cost constraint (3) was completely ignored. Eliminating this constraint also eliminated the need to evaluate the transit assignment model, drastically reducing the computation time needed for each search iteration. This trial was initialized with the final solution of the main trial and allowed to run for an additional 1000 iterations beyond the initial trial's 500 iterations. Between iterations 500 and 1500 the objective increased from 2.70226 · 10 −4 to 2.74911 · 10 −4 , for an increase of 1.7%. The mean absolute value in relative difference of community accessibility metric was less than 2.1%. Beginning approximately with iteration 800 nearly every move involved adding a vehicle to Harlem, which ended with 413 vehicles. Vicennes/111th and 31st also ended with significantly more vehicles than the other bus lines (113 and 67, respectively), although neither total changed after approximately iteration 900.
Excluding Express Lines: In order to determine whether the sharp improvement to the main trial's objective that occurred during iteration 186 was the result of the express lines, a restricted version of the trial was run which excluded all express lines. It was initialized with the main trial's solution at iteration 165, before the first express route was added, and allowed to run until reaching 250 iterations. Between iterations 165 and 250 this alternate trial's objective increased from 2.19014 · 10 −4 to 2.63517 · 10 −4 , for a 0.5% increase. During the same period the original trial had increased to 2.20147 · 10 −4 , for a 20.3% increase.

Sensitivity Results from the Artificial Network
Here we discuss the results of the different parameter tests separately. Table 11 in Appendix D.2 shows the full solution vectors for all artificial network trials. User Cost Increase Parameter: In the case of = 0.00 the solution algorithm made no changes to the initial network due to all neighbors of the initial solution being infeasible. All other values ∈ {0.01, 0.05, 0.10, ∞} produced exactly the same accessibility metrics for all communities, and thus exactly the same objective value, which increased from an initial value of 4.77960 · 10 −4 to a final value of 7.71440 · 10 −4 for a 61.4% increase. Table 8 in Appendix D.2 shows the full list of final community metrics. The solution vectors, themselves, were mostly the same across all trials. 12 of the 18 lines received the same number of vehicles in all solutions, with a maximum difference of 3 vehicles among the lines that did differ.
Community Inclusion Parameter: The number of communities K to include in the objective function directly affects its value even if the solution remains the same. For this reason we cannot directly compare objective values for separate trials, but we can still compare the accessibility metrics for individual communities. Table 9 in Appendix D.2 shows the full set of community access relative changes across the K trials.
The case of K = 1 corresponds to optimizing the single lowest community area metric. This trial made no changes to the initial network, meaning that no explored solution was able to improve the accessibility of the lowest-metric community. Most of the larger values of K showed similar trends in which communities improved their access and which worsened, with the exception of K = 30, which includes all communities in the objective (and thus corresponds to maximizing the average accessibility of all communities). Whereas the other trials generally improved the areas with the worst initial access at the cost of areas with the best initial access, the K = 30 trial reversed this pattern and improved primarily the areas with the best initial access.
The initial artificial network contains 15 communities with below-median access and 15 with above-median access. Table 2 shows for each trial the average relative change for the initially below-median communities and the initially above-median communities. As mentioned above, K = 1 resulted in no change. K between 6 and 24 all produced similar trends, achieving a small increase for the communities below median at the cost of a decrease for the communities above median. K = 30 instead produces a small decrease for the communities below the median and a significant increase for the communities above the median.
Gravitational Decay Parameter: The gravitational decay parameter β changes the evaluations of all community accessibility metrics, and so neither the objective nor the individual metrics can be directly compared between separate trials. Relative changes can still be compared, and are shown in Table 10 in Appendix D.2. All trials produced a qualitatively similar distribution of relative accessibility changes for each community, although the magnitudes of the relative differences tended to increase as β increased. All trials produced the greatest relative increase and the greatest relative decrease for the same two communities. The rankings of the relative increases remained similar between trials, and 19 of the 30 communities showed the same sign of accessibility metric difference (increase or decrease) in all trials.

Discussion and Conclusions
The results presented in Section 5 are mixed but generally positive with respect to evaluating the effectiveness of the SAMP in meeting our modeling goals. While the main trial for express route design for the CTA network was able to significantly improve the primary care access levels of the most underserved areas of the city, the artificial network trials indicate that the model's output can be sensitive to certain parameters.
The results shown in Section 5.1 indicate that the SAMP shows promise as a model for achieving better social equity via minor alterations to public transit. The main trial was able to achieve improved access levels for the least advantaged areas of the city. While this did require the access levels of the most advantaged areas to fall, the changes described in Figure 4 and Table 1 indicate that the relative gains far outweighed the relative losses. The overall distribution of accessibility metrics within the city remained relatively balanced, with approximately half improving and half worsening while the mean and median remained roughly the same, both changing by less than 3.2%. However, the spread of accessibility levels also narrowed significantly, with the standard deviation decreasing by more than 33.6%, the maximum decreasing by more than 3.6%, and the minimum increasing by more than 45.3%. Achieving a smaller spread in differential accessibility levels, with a particular focus on improving the worst case, is exactly the goal in improving equity of access, and in that regard the SAMP succeeded in its stated goal. Moreover, although the solution algorithm was halted after 500 search iterations, the relaxed run discussed in Section 5.2 shows that relatively little improvement would be achieved if the algorithm had continued.
The form of the objective function Access(y) defined in (8) seems to have performed well. The resulting network generally improved the accessibility of not only the K = 8 lowest-access areas but the surrounding areas as well. The artificial network results shown in Section 4.2, however, indicate that the algorithm's performance is dependent on the choice of K. While most intermediate values of K seem to result in qualitatively similar solutions, a very large choice of K causes no particular focus to be given to the least advantaged areas and may result in their access worsening, while a very small choice of K may lead to over-focusing on communities that cannot be significantly helped by any changes to the public transit system while ignoring other communities that can. This effect was also shown within the CTA study in the community area O'Hare, which is located at the westernmost bound of the Far North Side. O'Hare has relatively little CTA bus route access, and so none of our design decisions y could have a significant effect on the accessibility metric of O'Hare. By the end of the main trial its metric had increased by less than +0.4% in spite of being one of the 8 lowest metrics (and thus part of the objective) throughout most of the algorithm's run time. The relaxed trial that ran for 1000 iterations after the initial trial had begun to divert a very large number of vehicles to the Harlem line, one of the few public transit connections for O'Hare, but was able to increase the objective by only +1.7%.
The express run design component introduced in Section 5.2 produced some unexpected results, and while the introduction of express runs seems to help with the accessibility objective, more significant system-wide changes are required before their effects are seen. The sharp improvement in the objective seen near iteration 200 of the main search in Figure 3 did not occur in the trial without express routes available, indicating that they were essential for that improvement. However the trial in which only swaps between a route and its express run was unable to attain any improvement over the initial network. These results indicate that the express lines are not immediately helpful for the initial network, but that they do become helpful after significant alterations have been made to the network. Moreover, the total number of vehicles on express routes was fairly small, at only 4.25% of the total bus fleet, which supports the idea that it is possible to significantly improve access levels with very little investment in new transit infrastructure.
The gravity decay metric trials described in Section 5.3 indicate that the model's performance can depend greatly on the choice of parameters, echoing the results of the K trials. For most values of β a similar set of communities showed an accessibility metric increase, indicating that all versions of the trial resulted in focus being given to a similar set of communities, however the magnitudes of the relative changes and the solution vectors were significantly different, indicating that the model is sensitive to how accessibility metrics are measured.
The tests of varying user cost increase bounds discussed in Section 5.3 illustrate an important point regarding constraint (12) of the SAMP. The case of = 0.00 resulted in all neighbors of the initial solution being infeasible, which is unsurprising since the initial fleet sizes were chosen to be locally optimal with respect to user costs. For computational reasons the choice of = 0.00 is not recommended for use. Due to the minor computational errors inherent in solving the transit assignment model [48], solutions whose user costs are very close to the upper bound UserCost(x * , ω * ) may falsely be calculated as being slightly above it, marking those solutions erroneously as infeasible. In practice should be chosen to be at least a very small positive value larger than some reasonable numerical tolerance for the user cost calculation.
For the larger values of all trials produced identical objectives and community metrics, al-though not quite identical solution vectors (see Table 11). The accessibility metric definition (9)-(10) is based on shortest path distance d ij (y) in G, and solution alterations that do not affect any community's shortest path have no effect on any community's accessibility metric. This is an undesirable property in an objective function since it produces regions of zero gradient within which there is no clear local search direction. Moreover, the fact that all values of produced nearly identical results indicates that the user cost constraint was not a significant impediment beyond the initial steps. This does not seem to be a realistic outcome for a real network, where it can be assumed that the initial network is already reasonably well optimized with respect to user cost. The modeling framework of the SAMP shows promise for use in the problem of achieving social equity through use of minor public transit infrastructure changes. In the Chicago case study it was able to significantly improve the primary care access levels of underserved communities with relatively few changes to the current network. This represents just one of many potential applications of this type of model, but there is reason to be optimistic about the potential for achieving important social goals with very little investment. Due to the incompleteness of our data set and the number of simplifying assumptions made, the results of this study should not be regarded as a serious set of recommendations for changes in the CTA network design. In order to actually make such planning decisions much more data would be required and a more sophisticated user assignment model would be recommended. It would also be valuable to build the design of express runs into the model, itself, rather than simply generating a small number of candidates using heuristics.
The model presented in this study was meant to be general in its form, and the choices of objective and constraint function should be made according to the needs of the planner. While our access objective defined in (8) was generally able to drive the types of improvements we had hoped to see for the community areas, the sensitivity analysis for K indicates that this specific form could be improved. It would be reasonable, for example, to explore an objective that includes all accessibility metrics but dynamically weighted towards the current lowest values, or one that includes a set of the lowest-access communities alongside their neighboring communities.

A Solution Algorithm Pseudocode
The pseudocode for the hybrid tabu search/simulated annealing algorithm discussed in Section 3 is included below. The source code for our C++ implementation can be viewed online [42]. For simplicity, x will be used to represent the entirety of the flow vector (although the assignment model also includes a waiting time scalar ω). The canonical vector e l ∈ {0, 1} |L| represents a vector with a 1 in coordinate l and 0 elsewhere. Sets N + , N − , and N × represent the ADD, DROP, and SWAP neighbors kept as candidates after the first pass of the neighborhood search, respectively, while N + , N − , and N × are the second-pass candidates. Bounds N max and N max are the required numbers of candidates to be gathered during each pass. Q max in and Q max out are the inner and outer nonimprovement counter bounds, respectively. BestAccess ← Access(y (0) ) // best known objective T ← T 0 // simulated annealing temperature 8: t ← t 0 // tabu tenure 9: Q in ← 0 // inner nonimprovement counter 10: Q out ← 0 // outer nonimprovement counter 11: for k = 0, 1, . . . , K do // search for a predetermined number K of local moves 12: // NEIGHBORHOOD SEARCH 13: N + ← ∅ // tentative ADD neighborhood 14: repeat // first pass (ADD) 15: l ← random element of L not yet considered during this loop 16: if y (k) + e l satisfies vehicle limit bounds (5)-(7) then // consider only design-feasible moves 17: if e l / ∈ STM or Access(y (k) + e l ) > BestAccess then 18: N + ← N + ∪ {l} // keep non-tabu candidates (unless improved best) 19: until | N + | ≥ N max or all lines in L have been considered 20: N + ← ∅ // final ADD neighborhood 21: repeat // second pass (ADD) 22: l ← unprocessed element of N + with maximum Access(y (k) + e l ) 23: x (k) ← TransitAssignment(y (k) + e l ) // calculate user flows for given move 24: if (y (k) + e l , x (k) ) satisfy operator and user cost constraints (2) until |N + | ≥ N max or all lines in N + have been processed 27: N − ← ∅ // tentative DROP neighborhood 28: repeat // first pass (DROP) 29: l ← random element of L not yet considered during this loop 30: if y (k) − e l satisfies vehicle limit bounds (5)-(7) then // consider only design-feasible moves 31: if −e l / ∈ STM or Access(y (k) − e l ) > BestAccess then 32: N − ← N − ∪ {l} // keep non-tabu candidates (unless improved best) 33: until | N − | ≥ N max or all lines in L have been considered 34: N − ← ∅ // final DROP neighborhood 35: repeat // second pass (DROP) 36: l ← unprocessed element of N − with maximum Access(y (k) − e l ) 37: x (k) ← TransitAssignment(y (k) − e l ) // calculate user flows for given move 38: if (y (k) − e l , x (k) ) satisfy operator and user cost constraints (2) if |N + | = |N − | = 0 then // handle the event of a failed search 42: if bounds N max were used in the previous neighborhood search then 43: restart neighborhood search while ignoring the bounds N max 44: else 45: Q out ← Q out + 1 // work towards diversification 46: delete from STM the tabu rule with the shortest remaining tenure 47: restart neighborhood search 48: if Access( y (k,1) ) − Access(y (k) ) > 0 then 60: add rules with tenure t to STM to prevent undoing the ADD/DROP moves of y (k+1)

64:
if Access( y (k,1) ) > BestAccess then 65: y best ← y (k,1) // update best known solution 66: BestAccess ← Access( y (k,1) ) 67: else 68: Q out ← Q out + 1 // nonimprovement iteration  apply cooling schedule to T 91: return (y best , BestAccess) After the algorithm terminates and we obtain the best known solution, we conduct an exhaustive local search starting at y best . In each iteration of this local search we consider every possible ADD, DROP, and SWAP move without considering tabu rules, always taking the best neighbor and terminating when no neighbors offer any improvement. The final result is guaranteed to be locally optimal.

B CTA Network Data Processing
This section gives a brief description of the data sources and processing used to generate the CTA computational trials described in Section 4.1. The source code used for data processing can be viewed online [42], as can the raw data used to define the network [43].

B.1 Public Transit Network
Most of the structural information to construct transit network representation G = (V, E) was obtained from the CTA General Transit Feed Specification (GTFS) files [17]. This includes a Stop Clustering: The GTFS files define over 11,000 stops in the CTA network. To reduce the size of the network a k-means process was used to partition the stops into clusters such that the total variance in position within all clusters was minimized. Each stop was then identified with the cluster with the nearest centroid, resulting in a set of 997 cluster centroids for use as stop nodes. The mean geodesic distance (WGS 84 ellipsoid) between a stop and its associated cluster centroid was 0.13333 miles, with a standard deviation of 0.08327 miles, a median of 0.12103 miles, and a maximum of 0.84410 miles.
Transit Lines: The GTFS stop time file was used to obtain the sequence of stops for each line l ∈ L, which were then mapped to the cluster centroids in the reduced network to define the stop nodes V l stop . Consecutive pairs of distinct stops defined line arcs ij ∈ E l line , with arrival time differences being used to define line arc travel times c ij . If two consecutive stops had been identified with the same cluster centroid, then rather than generating a line arc as a loop, the travel time between the stops was equally distributed among the incident line arcs in order to preserve the total line travel time τ l := ij∈E l line c ij . The total number of daily stop arrivals, divided by the number of unique stops, was used to calculate a daily average number of visits per stop v l . The initial average line frequency was then calculated as f * l := v l /τ l . This was used to estimate the initial fleet size as y * l := τ l f * l , obtained from the proportionality f l = y l /τ (see Section 2.2), taking the ceiling to ensure integrality.
The operational hours for the line were used to determine the fraction τ l f ∈ [0, 1] of the full 1440 daily time horizon τ during which the line was active. During the execution of the solution algorithm the current value of y l was used to update the frequency and capacity of each line. Line capacities were scaled by a factor of τ l f and calculated as u l := τ l f τ b l f l , the maximum number of people carried with service frequency f l over the actual operational hours of that line. Average line frequency was defined as f l := y l /τ l . The fleet size lower bound y min l for each bus line was taken as the minimum y l required to achieve a service frequency of at least one arrival per 30 minutes.
Walking Arcs: Walking arcs E walk were generated between pairs of clustered stops less than 0.75 miles from each other, measured using L 1 -distance (a.k.a. taxicab distance or Manhattan distance) under the assumption that walking paths within city blocks are generally restricted to the cardinal directions. To prevent excessive numbers of arcs from being generated by dense cliques of pairwise-adjacent stops, walking arcs were not generated between a pair of stops if a path through an intermediate stop existed. For each pair of connected stops, a bidirectional pair of walking arcs was generated with a travel time calculated from the L 1 -distance between the coordinates and an assumed walking speed of 4 ft/sec.

B.2 Origin/Destination Travel Demand
An OD travel demand matrix is required by the public transit assignment model discussed in Section 2.2, and contains the travel demand between every pair of locations within the city via public transit during the day-to-day time horizon, estimated from ridership data. Because in general only stop-level data are available, the transit stops V stop , themselves, are treated as the origins and destinations V + and V − . The only available CTA stop-level bus data was from October 2012 [15], and so stop-level train data from the same month were used [16]. Each data set included total numbers of daily boardings. This lack of time-dependent data was the main reason behind using the entire 24-hour period as the time horizon τ . Boardings at each stop were added to obtain an average total number of weekday boardings at each clustered stop. Alighting numbers were set equal to the boardings at each stop under the assumption that most daily trips are two-way.
These boarding and alighting numbers were used to estimate the OD matrix through iterative proportional fitting (IPF) [3]. The initial seed matrix was generated using a trip distribution function which assumed a gamma distribution of trip lengths [6]. Mean trip length was estimated using results from the 2007-2008 Chicago Metropolitan Agency for Planning (CMAP) Travel Tracker Survey [14] which provides mean trip lengths for bus and train trips. These means were weighted by the total numbers of bus and train boardings to obtain a mode-weighted mean of 43.8 minutes. The standard deviation was taken as 20 minutes. The seed matrix was then generated by applying the resulting gamma distribution to the pairwise trip lengths in G. The IPF process stabilized after 10 iterations with an error of less than 0.001, after which travel demands were rounded to the nearest integer.

B.3 Community Nodes
Community nodes V + represent population centers within Chicago, which were identified using 2010 census data. To reduce the size of the network, the over 800 census tracts were clustered into population-weighted centroids for each of the 77 Chicago community areas, grouping census tracts according to the community area to which they belonged. A community node was generated for each community area and connected via walking arcs to all stops within an L 1 -distance of 1 mile. Exactly one population center (O'Hare) was located more than 1 mile from the nearest transit stop, and so its walking arc cutoff distance was extended to 2 miles. The total population over all census tracts was 2,739,860.

B.4 Facility Nodes
Facility nodes V − correspond to primary care facilities within Chicago (see Figure 1). A list of primary care community health clinics and their geographic coordinates was obtained from the Chicago Data Portal [13]. These facilities all provide primary care, are open to everyone regardless of ability to pay, and have a specific mission to provide care for underserved populations. A facility node was generated for each primary care facility and connected via walking arcs to all stops within an L 1 -distance of 1 mile.

B.5 Candidate Express Route Generation
Candidate express routes for use in the design problem discussed in Section 4.1 were generated by applying a heuristic algorithm to all non-express CTA bus lines with at least 18 stops. Under the assumption that a beneficial express run would be one which makes it easier to reach a high-access area from a low-access area, a gravity metric A i (y * ) as defined in (9)-(10) was calculated for all transit stops i ∈ V stop on the line in question, treating each stop as a community location with population P i ≡ 1. The express run of the line then dropped all stops except for the 10% highest-access and the 10% lowest-access stops (rounded to the nearest integer). The line arcs of the express run had their travel times reduced by 40 seconds for each skipped stop. Express runs were assumed to have the same operation period τ l f as the original line. A total of 65 of the original 126 CTA bus lines were given an express run. These lines were used to define a modified version of network G for use in some of the computational trials.

C Artificial Network Generation
In this section we describe the procedural network generation process used to generate the small artificial test case described in Section 4.2. The network generation code and its resulting data can be viewed online [42,43]. The artificial network operated over a time horizon of τ = 60 minutes, and was based on a 7 × 11 square lattice of roads representing a 0.9 × 1.5 mile city arranged in blocks. Transit stops were located on the intersections of the lattice with either 0, 1, or 2 additional stops (chosen uniformly at random) interspersed along the road segments between intersections. A total of 18 transit lines were defined, with one for each North/South and each East/West road on the grid. Arcs were assigned travel times based on Euclidean length and a bus speed chosen uniformly at random from [0.3, 0.6] miles per minute, plus a stopping time chosen uniformly at random from [0.5, 0.75] minutes. Trips began at one terminal, moved to the opposite terminal, and returned to the first terminal, followed by a layover chosen uniformly at random from [20,30] minutes. Walking arcs connected each stop node to a maximum of 5 stop nodes falling within a cutoff L 1 -distance of 0.12 miles, with travel times assuming a 4 ft/sec walking speed.
A total of 30 community nodes V + were generated, with locations chosen uniformly at random within the city's 0.9 × 1.5 mile boundaries, and with the restriction that no two communities could be within 0.15 miles of each other in order to achieve somewhat even spacing. A population P i uniformly chosen from [1000, 3000] was assigned to each community i. A total of 12 facility nodes V − were generated by again choosing coordinates uniformly at random, but with no minimum pairwise distance threshold, resulting in some densely-covered and some sparsely-covered areas of the city. All facilities j were given a quality of S j ≡ 1. Walking arcs were generated to connect communities and facilities to stops within 0.25 miles. The final network included 527 nodes and 2206 arcs.
An OD matrix of stop-level travel demand was generated using techniques similar to those discussed in Appendix B.2. The trip distribution function was a gamma distribution with a trip length mean of 2.75 minutes and standard deviation of 1.25 minutes. Stop-level total boarding numbers were determined by dividing 30% of each community's population equally among all stops for which they were the nearest community. The seed matrix was perturbed by adding a random value chosen uniformly from [−2, 2] to each element, and then IPF was used to generate the OD matrix.
The initial fleet vector y * was generated by assigning a set total of 200 vehicles proportionally to each line according to its total travel demand, calculated as total OD demand for which that line was part of the shortest path through the network. This assignment was then followed with a local search of reassignment moves to minimize the user cost function (11) under the assumption that the initial network should be reasonably-well optimized for that goal. The resulting fleet vector assigned between 6 and 14 vehicles to each line.

D Computational Trial Data Tables
This section includes the raw data tables for the computational trials described in Section 4 and summarized in Section 5. The raw data compiled in these tables can be viewed online [43]. Tables   Tables 3-4 show the changes in all 77 Chicago community area accessibility metrics before and after running the express route design model for 500 iterations, and include the initial metric, final metric, absolute difference in metric, and relative different in metric before and after solving the model. The contents of these tables correspond to the summary statistics in Table 1 and to  Tables 5-7 show the changes in vehicle fleet sizes before and after running the express route design model for 500 iterations. Table 5 shows the heuristically generated express lines described in Section 4.1, along with the name of the line whose route they follow and the final number of vehicles. Table 6 shows the original routes which received additional vehicles while 7 shows the original routes which lost vehicles. Only nonzero changes are included in tables. Tables   Tables 8-11 show the results of the artificial network trials described in Section 5.3. All tables list communities in the same order, from least to greatest initial accessibility metric. Table 8 shows the community metrics (along with absolute and relative differences) resulting from varying the allowable increase in user cost . The cases of = 0.01, 0.05, 0.10, and ∞ all produced identical metrics and so are listed together, while = 0.00 produced no change and is not listed. Table 9 shows the relative change in community access metrics resulting from varying the number K of communities included in the objective. Note that the trial for K = 1 made no changes to the network, resulting in no relative change in accessibility metrics. Table 10 shows the relative community access changes resulting from varying gravitational decay parameters β. The trials for β = 1.5 and β = 2.0 made no changes to the network. Table 11 shows the solution vectors for all artificial network trials. The columns correspond to each line's fleet size. Lines E0-E7 correspond to the East/West lines in order from South to North, while N0-N9 correspond to the North/South lines in order from West to East.     Table 6: Positive changes to the fleet sizes of each existing line after running the express bus route design solution algorithm for 500 iterations. Only nonzero increases are shown. Table 11: Solution vectors for all artificial trial sets as well as the initial solution vector. Rows listing a parameter , K, or β correspond to the trial set in which the listed parameter was varied. The three starred ( * ) rows represent a single trial with = 0.01, K = 6, and β = 1.0.