Multiscale heterogeneous optimal lockdown control for COVID-19 using geographic information

We study the problem of synthesizing lockdown policies—schedules of maximum capacities for different types of activity sites—to minimize the number of deceased individuals due to a pandemic within a given metropolitan statistical area (MSA) while controlling the severity of the imposed lockdown. To synthesize and evaluate lockdown policies, we develop a multiscale susceptible, infected, recovered, and deceased model that partitions a given MSA into geographic subregions, and that incorporates data on the behaviors of the populations of these subregions. This modeling approach allows for the analysis of heterogeneous lockdown policies that vary across the different types of activity sites within each subregion of the MSA. We formulate the synthesis of optimal lockdown policies as a nonconvex optimization problem and we develop an iterative algorithm that addresses this nonconvexity through sequential convex programming. We empirically demonstrate the effectiveness of the developed approach by applying it to six of the largest MSAs in the United States. The developed heterogeneous lockdown policies not only reduce the number of deceased individuals by up to 45 percent over a 100 day period in comparison with three baseline lockdown policies that are less heterogeneous, but they also impose lockdowns that are less severe.


A The Multiscale SIRD Model
Each metropolitan statistical area (MSA) is separated into geographically grouped subregions with associated subpopulations. The multiscale SIRD model of a given MSA consists of a directed graph with a node for each subregion. The spread of COVID-19 within the subpopulation of each subregion is governed by a distinct SIRD model. Following the regional model outlined by Della Rossa et al. 1 , we model the dynamics of the COVID-19 epidemic spread within the subpopulation of each subregion with the following system of equations.
D v (t) = ΓI v (t), ∀v ∈ {1, . . . , n}, R v (t) = γI v (t), ∀v ∈ {1, . . . , n}, where S v (t), I v (t), D v (t), and R v (t) are the number of susceptible, infected, recovered, and deceased within subregion v at time t. We use n to denote the number of distinct subregions within the MSA, and we use the symbols v, w, k to index the subregions. We use N v to represent the population of subregion v, and N p v (t) to represent the number of people interacting in subregion v at time t. Furthermore, we use φ vw (t) to represent the fraction of the population from subregion v visiting subregion w at time t.
The parameters in this model are as follows.
• β ∈ [0, 1] is a base infection rate of COVID-19, which is assumed to be constant across all subregions.
• ρ v ∈ [0, 1] is a subregion-specific modifier, that is used to model heterogeneity in the infection rates of different subregions.
• γ ∈ [0, 1] is the rate at which infected members of the population recover.
The summation terms in (1) and (2) capture the spread of COVID-19 between the populations of the various subregions in the MSA. Equation (1) models the rate of changeṠ v (t) of the number of people from subregion v that are susceptible at time t. The value ofṠ v (t) is given by the summation over all subregions w of all the susceptible people from subregion v who become infected while in subregion w. Mathematically, we express this statement asṠ .
Here, NumInfected w (t) represents the number of the people in subregion w at time t who are infected, and N p w represents the total number of people in subregion w at time t. The value of NumInfected w (t) can be approximated as ∑ n k=1 φ kw (t)I k (t). Putting these equations together, we arrive at (1). Meanwhile, (2) uses the same expression to define the rate of changeİ v (t) of the number of people from subregion v who are infected at time t. However, (2) also includes the term −(γ + Γ)I v (t), which captures the decrease in the number of infected people due to recovery (4), or death (3). An estimate for the value of N p w is given by (5).
Note that we assume that the members of the population return to their home subregions at the end of each timestep. More specifically, φ vw (t) represents the fraction of subregion v's population that visits an activity site in subregion w on day t. However, region v's total population N v stays constant; N v is not reduced by a factor of φ vw (t) at time t + 1.
To obtain values for the parameters β , ρ v , γ, and Γ that have biological significance and that are consistent with our network SIRD model, we use the parameter values that were presented by Della Rossa et al. 1 . The primary focus of their work was to calibrate and analyze a similar network model to ours using data from the start of the pandemic in Italy.

B Adding Activity Sites and Lockdown Control to the Multiscale SIRD Model
In the multiscale SIRD model, the variable driving inter-regional spread is φ vw (t), which is a dimensionless interaction term between subpopulations. Any attempt at controlling the inter-regional spread of the virus through limitations of population mobility, will thus depend on influencing the values of φ vw (t).
One possible explanation for people's inter-regional interactions is the frequency at which they visit various types of activity sites within the different subregions. We use m to denote the number of distinct types of activity sides included in the model. For example, grocery stores, restaurants, and fitness centers are all examples of different types of activity sites. We use the symbol i to index of the different types of activity sites. We assume that the interactions at these activity sites are the only factor influencing φ vw .
For each subregion v ∈ {1, . . . , n} and activity site type i ∈ {1, . . . , m}, we assume there is some fixed demand rate d i v , which represents the average number of people from subregion v who typically visit an activity site of type i within a day. This quantity is estimated from the SafeGraph data 2 .
During the pandemic, many members of the population may either alter their demands for certain types of activity sites or decide to switch to socially distant means of satisfying their needs; for example, switching to online shopping instead of visiting stores. The number of inter-regional interactions will be affected by the fraction of the population choosing to alter their behavior in this way. Let o i v (t) denote the number of people from subregion v that have decided not to visit activity sites of type i per day. We assume a linear relationship o i v (t) = A i v L i v (t) + B i v holds between o i v (t) and the lockdown value L i v (t) at any given time t for all subregions v and all activity site types i. Here, A i v and B i v are parameters to be learned from survey data. We also assume that each type of activity site (e.g., grocery stores or restaurants) within each subregion has some maximum capacity C i v ∈ N representing the number of people that can be served by activity sites of type i within subregion v. To influence the interactions between subpopulations, and thus the inter-regional spread of the disease, we control the lockdown of each activity site type i within each subregion v by setting the parameter L i v (t) which takes values in [0, 1] for all times t in the considered time horizon [0, T ]. This parameter affects the available capacity by scaling the value of C i v . That is, L i v (t)C i v is the total number of people that can be served by an activity site of type i in subregion v at time t.
Furthermore, we introduce the variables f i vw (t) which represent the average number of people from subregion v who visit activity site i within subregion w at time t. The value of φ vw (t) -the fraction of subregion v's population who visit any type of activity site within subregion w at time t is then given by We model the value of f i vw (t) at any given time t by comparing the relative values of , the average demand per day for activity site i in subregion v is smaller than the capacity that is being served in that subregion. In this scenario, we have there is no motivation for members of subpopulation v to travel to other subregions in order to visit an activity site of type i.
Here, W i wv is a value in the range [0, 1] that represents the likelihood that a person will choose to travel to subregion w if there is no remaining capacity for activity sites of type i within their home subregion v. These values are estimated from SafeGraph data, as described in §D.3. We may thus concisely write the value of f i vw (t) for a given time t using the equations We note that while the multiscale SIRD model of (1)-(5) is formulated as continuous-time system of ordinary differential equations, a discretization of the dynamics in time is necessary for numerical simulation. Furthermore, the synthesis of optimal lockdown policies L i v (t) in continuous time requires the solution to an infinite-dimensional optimization problem that is computationally intractable due to the aforementioned non-smooth relationship between L i v (t) and φ vw (t). We accordingly use the forward Euler method to obtain a discrete representation of the system dynamics and we model the lockdown policies L i v (t) as zero-order-hold signals (i.e. the lockdown factor is assumed to be constant for each discrete timestep). This formulation leads to the discretized form of the system dynamics given below.
Given initial values for the number of susceptible S v,0 , infected I v,0 , recovered R v,0 , and deceased D v,0 members of the subpopulations of each subregion v at time t = 0, and given a selected lockdown policy L i v (t) for each type of activity site within each subregion, we write the discretized dynamics of the system as follows.
where T s is the sampling period, which is defined to be 1 day in our formulation.

C Formulating the Heterogeneous Optimal Lockdown Problem
In the multiscale SIRD model, the variable driving inter-regional spread is φ vw (t), the dimensionless interaction term between subpopulations. Given a finite horizon T , our goal is to schedule the partial lockdown capacities of the activity sites within all subregions to minimize the total number of deaths across the entire network at time T .
The constants in this problem are: . . , m} and v, w ∈ {1, . . . , n}. The variables in this problem are S v (t), 3/9 where α > 0 is a positive economic impact parameter. The optimization problem is then defined by the above objective subject to the constraints in (6a)-(6m)

C.1 The Relaxed Optimization Problem
The constraints in (6i) and (6j) involve min and max operators over two affine functions in the variables, respectively. Including these operators in an optimization problem results in a nonsmooth problem. Obtaining locally optimal solutions for such nonsmooth problems is significantly more challenging than smooth problems, both in theory and practice. A detailed discussion on the challenges and convergence guarantees for nonsmooth optimization problems can be found in 3 . Therefore, we relax these constraints as follows, i.e., the constraint in (6i) is replaced by and the constraint in (6j) is replaced by The practical implication of these modifications is, given a subregion v and an activity site of type i, we allow more people to travel to other cities to fulfill their demand by relaxing the constraint in (6j). This relaxation may increase the inter-regional spread of the disease. On the other hand, relaxing the constraints in (6i) may result in fewer people staying in their subregion while fulfilling their demand, which may decrease the regional spread of the disease.

C.2 Computing Locally Optimal Solutions to the Optimization Problem
The optimization problem in (6a)-(6m) with the objective (7) is nonconvex as the functions in the constraints (6b)-(6c) are nonlinear in the problem variables φ , S, and I. Therefore, it is in general computationally intractable to compute a globally optimal solution to the nonconvex optimization problem 4 . We thus focus on computing a locally optimal solution to the nonconvex problem in a reasonable time.
Our method is based on sequential convex programming (SCP) approaches 5,6 . SCP methods iteratively compute a locally optimal solution to a nonconvex problem by linearizing the nonconvex constraints around a previous solution. However, when a feasible nonconvex problem is linearized around a previous solution, it may result in an infeasible convex subproblem 5 .
To ensure the feasability of the linearized problem, we add slack variables to the linearized constraints. Additionally, the linearization may not be accurate if the solutions to the linearized problem deviate significantly from the previous solution. Therefore, we utilize trust region constraints 5,7,8 to ensure that the linearized problems accurately approximate the nonconvex optimization problem.
Linearizing the nonconvex constraints. For a given v, w, k, and t, the function in the constraints (6b)-(6c) are nonlinear in the following problem variables . In our optimization approach, we fix the term N p w (t) in each iteration, which improves the efficiency and numerical stability. Therefore, we linearize the function h vwk (t) where h vwk (t;x) denotes the evaluation of the function h vwk (t) atx and ∇xh vwk (t) denotes the gradient of the function h vwk (t) at x. Therefore, we linearize the constraints (6b)-(6c) as vwk (t) , ∀v ∈ {1, . . . , n},t ∈ {0, . . . , T − 1} (11) which are convex in the problem variables.

4/9
Adding slack variables and trust region constraints. Recall that the affine constraints in (10)-(11) may be infeasible, or the linearization may be inaccurate if the solution deviates significantly from the previous computed solution. We alleviate these drawbacks by adding slack variables k vwk (t) and l vwk (t) ∈ R to (10) and (11) respectively, and trust region constraints to the φ variables as where σ is the size of the trust region to restrict the set of allowed policies in the linearized problem. Finally, we augment the cost function in (7) with the term λ |k vwk (t)| + |l vwk (t)| to ensure that we minimize the violation of the constraints, where λ is a large enough positive constant to minimize the violations in the linearized constraints.
Checking the validity of the solution and updating the trust region. After solving the relaxed problem and obtaining solutions for the interaction terms φ vw (t) for all v, w, and t, we compute the actual cost of the solution by implementing the solution for the interaction terms and running the multiscale SIRD dynamics, defined in (6a)-(6g). If the actual cost for the obtained solution for the computed interaction terms is lower than the best-obtained cost, we update our solution and linearize the nonconvex optimization problem around the updated solution. When we update the solution, we also increase the size σ of the trust region. On the other hand, if the actual cost is larger, we conclude that our approximation is not accurate and shrink the trust region without updating the solution, and re-solve the linearized optimization problem. If the size σ of the trust region is sufficiently small (indicating that our approximation to the nonconvex optimization problem is accurate) and the algorithm cannot compute a solution that obtains a lower objective value than the last obtained solution, we conclude that the last obtained solution is locally optimal.

D.1 Data Collection
Our study was conducted on six of the largest metropolitan statistical areas (MSAs) in the United States, including New York, Los Angeles, Chicago, Dallas-Fort Worth, Phoenix, and Seattle. These areas are selected considering that they are among the fifteen most populated MSAs and are distributed in different parts of the country. In each MSA, the foot traffic information of activity sites in April 2020 are collected from SafeGraph 2 . The data include information on activity site ID, activity site location, NAICS code describing the particular type of business activity associated with each activity site, the census block groups with at least four residents who visited the activity site, and the associated amount of visitors. Based on the NAICS codes, we extract six main types of activity sites, including grocery stores, restaurants, hotels, pharmacies, fitness centers, and physicians.

D.2 Data Processing
As the SafeGraph data only includes the foot traffic information pertaining to a subset of the total population, we first calibrate the data by scaling all counts by the number of devices Sa f eGraphDevices v used to conduct the data collection for each subregion v. Using the calibrated foot traffic data for each type of activity site, we calculate the demand rate d i v for each subregion v and activity site type i. We recall that d i v is the average number of people from subregion v visiting an activity site of type i per day. We compute the value of d i v by Here, f req i denotes the frequency at which a person typically visits an activity site of type i per month. This scaling is necessary as the SafeGraph dataset only reports the number of unique visitors to an activity site per month; even if a person visits an activity site type multiple times in a month, SafeGraph only reports a single visit. Meanwhile, Population v represents the total population in subregion v, and FootTra f f ic i vw is the number of people in the SafeGraph dataset from subregion v visiting an activity site of type i in subregion w in a month.
Additionally, using data showing COVID-19 cases per zip code or per county on specific dates 9-14 , we find the number of cases in each subregion on a specific date that can be used as an initial condition in (6f) and (6g).

5/9 D.3 Obtaining the Intercity Flow Matrix and Edge Weights
We also compute parameters W i vw using SafeGraph data. We recall that W i vw are included in (6j) and represent the likelihood a person from subregion v would visit an activity site of type i in subregion w, given there is not enough capacity for activity sites of type i in subregion v. We compute the values of parameters W i vw by Here, Ad j k (v) is the set of subregions w for which there exists a directed edge from subregion v to w. For any value k, representing the maximum number of outgoing edges from a subregion, we compute Ad j k (v) using the values of FootTra f f ic i vw . Specifically, we compute and we define Ad j k (v) as the set containing the subregions w with the k largest values of TotalTra f f ic vw .

E Extended Results
In this section, we demonstrate the effects of the presented multiscale heterogeneous lockdown policies on the dynamics of the multiscale SIRD model that are estimated using data for Phoenix, Seattle, New York, Chicago, Los Angeles, and Dallas-Fort Worth MSAs. In the following figures, we visualize the resulting dynamics without any imposed lockdowns in the left column, and with an optimal lockdown in the right column for each MSA. We note that in all MSAs, there is no noticeable peak for the infections after imposing the lockdown policy, and the number of cumulative infections and deceased people decreases significantly.  In the left column, we show the dynamics of the multiscale SIRD model without any lockdown policy. In the right column, we show the dynamics after applying an optimal lockdown policy. In all cases, we observe a significant decrease both in the peak infection rate and also in the number of cumulative infected and deceased individuals.