1.1 Background

Primary care is essential to equity, providing a “front door” to preventative healthcare and access to specialty care (Starfield 2011). Healthcare accessibility can be defined as the opportunities (such as public transportation) for residents at various locations to obtain healthcare services (Wang 2012). This project is an equity-focused healthcare effort to re-imagine the infrastructure for primary care in the Chicago Area.

COVID highlighted that Chicago’s healthcare infrastructure is driving inequity and harming the city’s capacity for vitality and growth. The pandemic also upended payment models and approaches to care delivery, causing new thinking and receptivity among health system leaders about future strategies. In 2022, Liu et. al. (2022) presented an analysis showing that minority‐dominated neighborhoods (Black/Hispanic) in the west and south of Chicago suffer from low transit‐based healthcare accessibility in the Chicago Area. Since 2015, the Rush College of Nursing has partnered with community-based organizations to evolve a nurse-run, community-embedded primary care service model; it currently runs 32 such programs and clinics within the Chicago metropolitan area (Moss, n.d.; Moss et al. 2022). UChicago Medicine is actively investing in community-based care capacities; Rush University has received a $10m gift and hired new leadership for its Institute for Health Equity.

Driven by persistent and deadly deficits in human health among residents of Chicago’s most under-invested neighborhoods (Musick et al. 2021; Orsi, Margellos-Anast, and Whitman 2010), our project aims to combine a cost-effective, community-tailored primary care model with a geographic optimization strategy that locates primary care centers. SoReMo fellow Alaittin Kirtisoglu, a PhD student in Applied Math, works on developing a network optimization model for facility location and writes a programming code for the actual implementation of the model. Faculties Hemanshu Kaul from Applied Math and Kim Erwin from the Institute of Design are the collaborators in the project.

1.2 Main Problem

The model is based on the well-documented Cuban model, a decentralized network of micro community-embedded clinics anchored by a single physician and nurse to serve every 1200 residents (Keck and Reed 2012). The Cuban model is widely cited as an effective public health strategy that achieved full population coverage within 20 years of launch and has dramatically transformed the health of residents since its deployment in the 1970s (Esposito et al. 2016). By comparison, the average patient panel for a US primary care provider is 2300, a number calibrated for productivity and revenue but widely cited as contributing to access and quality issues (Raffoul et al. 2016). Deriving the number and placement of clinics for Chicago requires incorporating multiple factors to produce a real-world strategy. These factors include neighborhood population density; geographic features (highways, waterways, industrial zones); and other urban features known to drive or inhibit pedestrian traffic (street patterns, public transit, shopping districts). Taking these features into consideration, we generate a network that represents the Chicago Area. The model is an optimization model designed on that network and minimizes the difference between the maximum travel times of communities from their blocks to their facilities.

1.3 Design and Implementation

The model must be integrated with data-based computational tools. We use the Census data at the block level, the Chicago neighborhood and boundary data, the Chicago primary health center location and utilization data, and the Transportation data from the CTA.

Figure 1.1: A network example.

In combinatorics, a network is a collection of nodes and edges, where each edge connects a pair of nodes. Figure 1.1 shows an example of a network consisting of 15 nodes and 22 edges. We generate a network for the Chicago area in which nodes are the Census blocks and edges connect the nodes corresponding to neighboring Census blocks. Note that Chicago has approximately 10,000 Census blocks. The network is enlarged by adding new nodes for bus and train stations obtained from the data to calculate the travel time between the nodes. We set candidate locations for new primary care centers by using Google API considering the proximity of the locations to roads, community centers, and bus and train stations. The locations of the existing and candidate primary care centers are stored in the network. A solution of the model has to answer two questions:

1. Which candidate locations should we locate new primary care centers at?
2. Which primary care center should we assign an individual to?

The second question is crucial since each facility has a maximum capacity (1200 residents). We need to assign each census block to exactly one open facility. For this reason, our model splits Chicago city into districts. A district can be seen as the service area of a facility. Since each facility is considered to have the same capacity, the districts should be almost equally populated. To do this, we utilize the political redistricting (gerrymandering) optimization models (Hess et al. 1965). Finally, we attempt to solve the model by programming implementations of large-scale computational algorithms, e.g. Tabu Search and Simulated Annealing algorithms that are proven to work on such problems efficiently (Rumpf and Kaul 2021).

1.4 Conclusion

Question 2 requires a political redistricting model. The problem turns out to be a special version of the capacitated \(p\)-median problem (Osman and Ahmadi 2007). While the general version of the \(p\)-median problem is very difficult to solve on a large scale data (Gnägi and Baumann 2021), solving this special version becomes a bigger challenge. We implemented different models and solution methods during the fellowship. Our observation is that classical approaches are not sufficient to tackle the problem. Instead of using fully deterministic solution methods, we improved a new algorithm including some probabilistic arguments. The algorithm is in the coding process, and so we have not observed the results yet.

2.1 Network Construction

Let \(V\) be the set of \(n\) nodes corresponding to \(n\) Census blocks (land parcels) in Chicago. For simplicity, we will omit to say corresponding to Census blocks. For any \(u,v \in V\), we have \((u, v) \in E\) if and only if the land parcels corresponding to \(u\) and \(v\) share a border of non-zero length. By this construction, \(G\) is a simple planar graph. The population at \(v_i \in V\) is denoted by \(p_i\), and it is assigned as the weight of the node for \(i \leq n\). The locations of candidate primary care centers are set by using Google API as we explained previously. We obtain the locations of existing centers from the Chicago data portal (C. P. C. C. H. C. Data, n.d.). Each node in the network keeps information about whether it has a candidate or existing center location within its boundary.

Let \(V_s\) be the set of nodes corresponding to stops of transportation lines. \(E_s\) defines the set of edges connecting consecutive stop nodes lying on the same line. The directed graph \(G^{\prime} = (V \cup V_s, E \cup E_s)\) is constructed for only calculating the travel times between pairs of nodes in \(G\). Precisely, calculating the minimum travel time between a pair of nodes is a route-choice problem (Carey 1994). For each pair of adjacent nodes in \(G^\prime\) such that one of the nodes must be a block node, walking time is calculated under the assumption of 4 ft/sec walking speed, and it is assigned as weight to the edge between the nodes. Travel times between adjacent stops are averaged constant values gathered from CTA transit schedule data (C. S. S. Data, n.d.). Since our concentration is not the transportation network, we safely omit the details of vehicle capacity, line frequency, and waiting time which are the variables making the route-choice problem complicated. Finally, we calculate the best routes in \(G^\prime\) in polynomial time by using an extended version of Dijkstra’s shortest path algorithm (Noto and Sato 2000).

2.2 Facility Location Problem

Discrete facility location problem has been widely studied in location science. The problem seeks to find proper locations on a discrete network to build up a sufficient number of facilities under some restrictions. The problem has many applications in industry. For instance, finding switching centers in a telephone communication network to collect messages before they are sent to their destinations is a facility location problem. The main concern is to minimize the total length of the wires in the network. Locating emergency services such as ambulances, medical services, and police stations are also facility location problems. These emergency services should be located in a city minimizing the maximum service distances of facilities. A decision maker may need to restrict the total cost, number of employees, vehicles, or the radius of the service areas according to need and problem structure. In general, the problem is classified into five categories: uncapacitated facility location, capacitated facility location, \(p\)-center problem, \(p\)-median problem, and maximal covering location problem. We study capacicated facility location and \(p\)-median problems that are shortly explained below. A detailed review of the variants and their applications can be found in (Drezner 1996).

Let \(I\) be the set of \(n\) demand points, and \(J\) be the set of candidate facilities. The problem seeks to find a subset of \(J\) such that each point in \(I\) is assigned to exactly one center in \(J\), and the sum of cost and total travel time of assigned points to their facilities is minimized. A simple model of the capacitated version is defined as follows. The decision variable for locating candidate centers is

\[\begin{equation*} y_{j} = \begin{cases} 1 & \text{ if candidate facility } j \in J \text{ is located}\\ 0 & \text{ otherwise} \end{cases} \end{equation*}\]

and the decision variable to assign demand points (census blocks) in \(I\) to centers in \(J\) is

\[\begin{equation*} x_{ij} = \begin{cases} 1 & \text{ if block } i \in I \text{ is assigned to facility } j \in J \\ 0 & \text{ otherwise } \end{cases} \end{equation*}\]

Let \(t_{ij}\) be the travel time between demand point \(i\) and center \(j\), and \(p_i\) be the demand (population) at demand point \(i\). \(P_j\) is the maximum workload of facility \(j\), and \(c_j\) is the cost of opening center \(j\). A simple facility location model can be designed as \[\begin{align} &\text{min} \quad \sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n} t_{ij}x_{ij} + \sum\limits_{j=1}^{n} c_j y_{j} \hspace{15mm} \tag{2.1} \end{align}\] subject to \[\begin{align} & \sum\limits_{j=1}^{n} x_{ij} = 1, \quad i \in I \tag{2.2}\\ & \sum\limits_{i=1}^{n} p_i x_{ij} \leq P_j, \quad j \in J \tag{2.3} \\ & x_{ij} \leq y_{j}, \quad i \in I, \hspace{2mm} j \in J \tag{2.4} \\ & x_{ij}, y_j \in \{0,1\}, \quad i \in I, \hspace{2mm} j \in J. \tag{2.5} \end{align}\]

The objective function (2.1) minimizes the sum of the total distance from demand points to their facilities and the total cost of opening facilities. Capacity constraints (2.3) ensure that the workload of facilities is not exceeded. Constraints (2.4) guarantee that each demand point \(i\) is assigned to a facility \(j\) only if facility \(j\) is open. Constraints (2.5) are integrality constraints.

Distance becomes a key factor when patients choose a healthcare service. However, the capacitated facility model does not only focus on distance. Constraints (2.3) limit the number of doctor-nurse teams to be employed at each facility. Because, a fixed number of patients is considered the workload of a doctor nurse-team for a facility in a proper model, and \(P_j\)’s are set accordingly. We designed a capacitated facility model including one more decision variable for determining the number of doctor-nurse teams at facilities. The third decision variable helps not only limit the number of teams but also determine more precise numbers. We coded up an implementation of the model. However, we did not obtain a potent result. Solving three decision variables with a local search algorithm made neighborhood searches superficial. Because the 3-dimensional solution space is too large resulting in classical search algorithms being too slow on large-scale data. To solve this issue, we designed a bilevel optimization model, which is a composition of two submodels called lower and upper-level (Colson, Marcotte, and Savard 2007; Dempe and Zemkoho 2020). The lower-level model aims to locate primary care centers, while the upper-level model assigns Census blocks to the located centers. Since an optimal solution of the lower-level model may not cause an optimal solution in the upper level, solving the bilevel optimization model was another challenging problem that might produce new issues for us. As a result, we were not convinced we could obtain a remarkable solution.

In the \(p\)-median problem, the number of facilities to be opened is pre-specified with a given \(p\). Let \(I\) be the set of \(n\) points. The problem seeks to find a set of \(p\) centers among \(n\) points such that each point is assigned to exactly one center while the total travel time of assigned points to their centers is minimized. In the capacity version of the problem, facilities have nearly equal workloads, so the total populations assigned to the centers should be almost equal. Let \(P\) be the maximum workload of a facility. The capacitated version can be formulated as a mathematical programming model as follows. The decision variable is

\[\begin{equation} x_{ij} = \begin{cases} 1 & \text{ if point } i \in I \text{ is assigned to } j \in I\\ 0 & \text{ otherwise. } \end{cases} \tag{2.6} \end{equation}\]

Facility \(j\) is opened if and only if \(x_{jj} = 1\). The general \(p\)-median model is

\[\begin{align} \text{min } &\sum\limits_{i=1}^{n} \sum\limits_{j=1}^{n} t_{ij}x_{ij} \hspace{10mm} \tag{2.7} \\ \textrm{subject to} \quad & \sum\limits_{j=1}^{n} x_{jj} = p \tag{2.8} \\ & \sum\limits_{j=1}^{n} x_{ij} = 1, \quad i \in I \tag{2.9} \\ & \sum\limits_{i=1}^{n} p_i x_{ij} \leq P, \quad j \in I \tag{2.10}\\ & x_{ij} \leq x_{jj}, \quad i,j \in I \tag{2.11}\\ & x_{ij} \in \{0,1\}, \quad i,j \in I \tag{2.12} \end{align}\]

The \(p\)-median problem can be seen as a clustering problem. The key point of the \(p\)-median model is that demand is not sensitive to the level of service, that is, the number of doctor-nurse teams in a facility, and therefore travel time is more prominent in the assignment. This encouraged us to design a \(p\)-median model. Unlike hospitals, primary care centers provide basic front-door service, so we can safely ignore the level of service and consider the centers identical. We will discuss our model and its results after providing the necessary definitions.

2.3 Facility Location for Healthcare Systems

The first discrete facility model for healthcare systems was proposed in 1964 by Hakimi (1964). For the general methods and applications in the literature, one can see the surveys (Afshari and Peng 2014; Ahmadi-Javid, Seyedi, and Syam 2017). Lin et. al. (2023) categorize the models in terms of user choice behavior as follows.

All the models we presented so far are system-optimal. A central decision-maker determines the facilities of users. Another noteworthy point for the models is that all users living in the same block are assigned to the same facility. These assignments are called all-or-nothing. We observe that some non-equilibrium models use this approach as well. On the other hand, Huff-based models (see (Gu, Wang, and McGregor 2010)) in the non-equilibrium category allocate only a portion of demand to the same facility. For this reason, there is no need to determine which demand point goes to which facility explicitly. Instead, catchment areas of blocks and facilities are considered for calculating user flows (workloads of facilities) and attractiveness of facilities. A simple formulation is presented as follows. The probability of a client at block \(i\) traveling to a facility \(j\) is defined as

\[\begin{align*} Pr(ij)= \dfrac{S_j \Big{/} t_{ij}^{\lambda}}{\sum\limits_{j=1}^{n} S_j \Big{/} t_{ij}^{\lambda}}, \end{align*}\]

where \(S_j\) is the size of facility \(j\) and \({\lambda}\) is a parameter for setting the effect of travel time. For simplicity, take \(\lambda=1\) and assume that facilities are identical, i.e., \(S_j=1\) for every facility \(j\). The probability is simplified to

\[\begin{align*} Pr(ij)= \dfrac{t_{ij}^{-1} } {\sum\limits_{j=1}^{n} t_{ij}^{-1}}. \end{align*}\]

The user flow from demand point \(i\) with demand \(p_i\) to facility \(j\) is measured by

\[\begin{align*} f_{ij}=\sum\limits_{i \in I \cap t_{ij} \leq t^\epsilon} Pr(ij) p_i. \end{align*}\]

The catchment area of facility \(j\) is the subset of demand points, \(I^\prime \subset I\), whose locations are within a travel time of \(t^\epsilon\). Let \(p_i\) be the population at facility \(j\). The regional availability of facility \(j\) is described with the facility-to-client ratio

\[\begin{align*} R_j= \dfrac{1}{\sum\limits_{i \in I \cap t_{ij} \leq t^\epsilon} p_i}. \end{align*}\]

A higher ratio indicates that fewer clients share facility \(j\) and vice versa. The catchment area of demand point \(i\) is the subset of facilities whose locations are within a travel time of \(t^\epsilon\). The accessibility to a demand point is calculated by

\[\begin{align*} A_i = \sum\limits_{j \in J \cap t_{ij} \leq t^{\epsilon}} \dfrac{R_j}{t_{ij}}y_j \end{align*}\]

where \(y_j\) is the decision variable for locating candidate facility \(j\). In this formulation, people may go to any facility within the acceptable travel time \(t^\epsilon\). More facilities within the acceptable travel time of a demand point provide clients with a higher possibility of accessing healthcare services. If the number of opening facilities is bounded in the constraints, the objective function can be defined as the maximization of the sum of \(A_i\)’s.

In real life, patients do not necessarily need to select the nearest facility for preventive health services. Huff-based or equilibrium models involve a competitive environment, unlike classical studies. Therefore, the use of these models may be more convincing for some researchers. We improved a Huff-based model, however, we did not attempt to solve it. The inequality of access to primary care centers in Chicago is not distributed homogeneously and resources in some regions are drastically inadequate (Liu, Kwan, and Kan 2022). We aim to prioritize disadvantaged areas by providing a doctor-nurse team for every 1,200 residents, following the Cuban model. Additionally, we seek to enhance access equality by reducing travel time, which requires a continued focus on these underserved regions. To ensure that underserved neighborhoods have the minimum conditions that meet their basic needs, we omit user competition and designed a version of the system optimal model \(p\)-median. We select \(p\) locations from a set of candidate locations rather than selecting those centers from the set of all vertices in the graph. We derived the value of \(p\) by dividing the population of the city of Chicago by 1200. For the solution methods of uncapacitated and capacitated p-median models, see (Resende and Werneck 2004; Osman and Ahmadi 2007; Gnägi and Baumann 2021).

Two crucial concepts are needed to be considered in the system optimal models. A district is called contiguous if its induced subgraph in the network is connected. Compactness is a powerful criterion that forces districts to have regular geometric shapes (in terms of travel times) avoiding banana or octopus-shaped districts being drawn. The natural expectation is to have contiguous and compact districts. In general, compactness is not hard to design and solve. Even without a compactness constraint, minimization of the objective function (total deviation of the maximum travel times of almost equally populated districts) forces the districts implicitly to be compact. On the contrary, it is very difficult to obtain contiguous districts, especially on large-scale data. One way of achieving contiguity is to find a connected subgraph partition of the network, which can be designed as an optimization model (Achuthan and Caccetta 1992; Gruber 2009). Note that bounding the diameters of the subgraphs makes the contiguous districts compact. In formal, for a given graph \(G\) and a number \(k\), the \(p\)-median model needs to search additionally for a partition \(V_1, V_2,...V_k\) of \(V(G)\) in such a way that induced graph \(G[V_i]\) is a connected subgraph of \(G\) for every \(i\in [k]\). The partitioning problem itself is challenging. Instead of solving it as an optimization problem, it is obtained by generating a network spanning tree as explained in the next section.

2.4 Solution Method

We attempt to solve the model by using the Simulated Annealing algorithm (SA) (Kirkpatrick, Gelatt Jr, and Vecchi 1983; Van Laarhoven et al. 1987), which works on facility location problems efficiently (Murray and Church 1996). SA needs an initial solution as a starting point for the search in the solution space. Therefore, it has to satisfy all the constraints of the model. After choosing \(p\) centers from the set of candidate locations, we generate a spanning tree of the network by using Kruskal’s algorithm. Then, we pick two arbitrary centers in the spanning tree and cut a random edge on the unique path connecting them. We repeat this process \(p-1\) times in the components with at least two centers. Ricca et. al. (2008) use the same method for a political redistricting (gerrymandering) model. Note that the political redistricting problem also seeks a connected partition with the properties of contiguity, compactness, and population equality on large-scale data. For the solution methods in gerrymandering, see (Ricca and Simeone 2008; Validi, Buchanan, and Lykhovyd 2022).

To proceed with the simulated annealing algorithm, we need to define a neighboring relationship on the points in the solution space, so the algorithm moves in the space through neighbors to find a solution with the best objective value. Boundary of a district is the set of vertices adjacent to vertices in a different district. We consider two partitions \(C\) and \(C^\prime\) as neighbors if and only if \(C^\prime\) can be obtained by changing the district of a boundary vertex. SA checks the neighbors of the current solution to find a “better” solution. Once a “better” solution is found, it updates the current solution as long as contiguity is preserved. We indicate the better solution with a quotation mark because neighbors with worse objective values might also be accepted to avoid getting stuck in local optima depending on an acceptance probability function.

There are two important issues left to handle. (i) To simplify the model, we ignored population equality. We define a second objective function to measure the population equality of districts in a solution and replace “optimal solution” with “Pareto optimal solution”, which is defined for the models with two objective values (Ricca and Simeone 2008). (ii) We chose \(p\) centers from candidate locations and ran the algorithm for those centers. In other words, we have already decided which centers the model will locate. To tackle this issue, we plan to re-run the algorithm for different combinations of \(p\) locations and find the best \(p\) centers. However, a better strategy is needed since re-running the algorithm for all possible combinations is too costly.

2.5 Future Work

The local search algorithm we describe is called the flip-based method since SA moves only one block in every iteration. Flipping only one vertex makes the algorithm too slow. In 2021, Deford et. al. (DeFord, Duchin, and Solomon 2021) introduced a Markov chain method called recombination (ReCom) for gerrymandering. A random spanning tree is generated by using Wilson’s algorithm (Wilson 1996). Cut edges are selected such that \(k\) districts are obtained with almost equal populations. To improve the initial state, ReCom selects two neighboring districts randomly, merges them, and resplits them in every state of the chain without breaking population equality. Note that merging two districts is taking the subgraph induced by their vertices. Resplitting is achieved by generating a spanning tree and cutting an edge similar to the initial state. It has two remarkable improvements for us compared to the local search algorithm. (i) ReCom moves many vertices in every state rather than doing a single flip. (ii) Population equality is provided implicitly while transitioning to the new state. We plan to revise the ReCom algorithm for facility location problems. The improvements may allow us to focus on center selection. However, we need to re-imagine the algorithm for designing the districts around centers.