# Network-Based Epidemic Control Through Optimal Travel and Quarantine
Management ††thanks: Research partially supported by the NSF under grants
DEB-2433726, ECCS-2317079, CCF-2200052, and IIS-1914792, by the ONR under
grant N00014-19-1-2571, by the DOE under grant DE-AC02-05CH11231, by the NIH
under grant UL54 TR004130, and by Boston University.

Mahtab Talaei, Apostolos I. Rikos, Alex Olshevsky, Laura F. White, and Ioannis
Ch. Paschalidis M. Talaei, A.I. Rikos, A. Olshevsky, and I.Ch. Paschalidis are
with the Department of Electrical and Computer Engineering, and the Division
of Systems Engineering, Boston University, Boston, MA, USA. E-mails: {mtalaei,
arikos, alexols<EMAIL_ADDRESS>I.Ch. Paschalidis is also affiliated with
the Department of Biomedical Engineering and the Faculty for Computing & Data
Sciences, Boston University, Boston, MA, USA.L.F. White is with the School of
Public Health, Boston University, Boston, MA, USA, e-mail<EMAIL_ADDRESS>results in this paper were not presented at any conference.

###### Abstract

Motivated by the swift global transmission of infectious diseases, we present
a comprehensive framework for network-based epidemic control. Our aim is to
curb epidemics using two different approaches. In the first approach, we
introduce an optimization strategy that optimally reduces travel rates. We
analyze the convergence of this strategy and show that it hinges on the
network structure to minimize infection spread. In the second approach, we
expand the classic SIR model by incorporating and optimizing quarantined
states to strategically contain the epidemic. We show that this problem
reduces to the problem of matrix balancing. We establish a link between
optimization constraints and the epidemic’s reproduction number, highlighting
the relationship between network structure and disease dynamics. We
demonstrate that applying augmented primal-dual gradient dynamics to the
optimal quarantine problem ensures exponential convergence to the KKT point.
We conclude by validating our approaches using simulation studies that
leverage public data from counties in the state of Massachusetts.

###### Index Terms:

epidemics, networked control systems, stability of nonlinear systems,
compartmental models, optimization.

##  I Introduction

As the recent COVID-19 pandemic has demonstrated, effectively controlling the
spread of infectious diseases is of paramount importance [1, 2]. When
epidemics are viewed through the lens of interconnected networks (e.g., as
cities, counties, or other societal structures), a profound challenge is to
preserve the continuity of societies and minimize socio-economic disruptions.
To address this challenge, control systems methods can play an important role
[3, 4], particularly by providing tools for designing and implementing
effective control policies that can minimize the impact of epidemics on public
health, economies, and societies in general.

Network-based epidemic spread models employ a compartmental approach.
Individuals are assigned to different compartments (e.g., susceptible,
infected, recovered) based on their disease status [3, 4]. To control disease
spreading, most works in the literature employ two types of methods: optimal
control and spectral optimization. In the optimal control literature, the main
focus is to optimize various model parameters [5, 6, 7, 8]. These parameters
are optimized with the objective of minimizing a cost function that accounts
for undesirable outcomes (e.g., number of infected). The spectral optimization
literature typically considers resource allocation problems [9, 10, 11, 12].
The aim is to contain the spread of the epidemic by minimizing the maximum
eigenvalue of the matrix that describes the model’s behavior. However, given
the complex nature of epidemic models, designing and implementing policies to
help mitigate the effects of an epidemic is computationally expensive and/or
NP hard (e.g., from the perspective of spectral optimization, the optimal node
and link removal problems are NP-complete and NP-hard, respectively [13]).
Recent works focused on designing computationally efficient epidemic control
strategies [14, 15]. Motivated by these works, we aim to provide
computationally efficient frameworks for epidemic control that optimize
allocation of available resources and minimize the impact on society.

Main Contributions. We propose two novel strategies for network-based epidemic
control. The first leverages optimal control principles. Specifically, our
model’s input parameters are optimized to minimize a cost function integrated
along the model trajectory. The second strategy draws inspiration from
spectral optimization techniques. It focuses on efficient resource allocation
to minimize the spectral abscissa of the model’s matrix, thereby achieving
specific performance guarantees. Our contributions are as follows:

  * $\bullet$

Optimizing Travel Rates: We present a novel approach to suppress the spread of
an epidemic by optimizing the travel rates between network nodes (see Section
IV). We analyze the convergence of our approach, and show that not only it
effectively curbs the transmission of infectious diseases, but also it
preserves societal and economic activities. Specifically, its operation relies
on minimizing the dominant eigenvalue of the infection spread matrix. This
underscores the centrality of the network analysis in epidemic control.

  * $\bullet$

Optimizing Quarantine Rates: We extend the SIR model to the SIQR model by
incorporating a quarantine strategy for both asymptomatic and symptomatic
individuals (see Section V). For the SIQR model, we demonstrate that the
optimization of quarantine costs can be formulated as a matrix balancing
problem, which is solvable in polynomial time. Furthermore, we elucidate the
relationship between the constraints of this problem and the basic
reproduction number of the epidemic. As an extension, we show that applying
augmented primal-dual dynamics to this problem achieves semi-global
exponential convergence to the KKT point of the quarantine optimization
problem.

  * $\bullet$

Model Validation: We validate the effectiveness of our approaches using real-
world network data (cf. Sec. VI). Namely, we calculate optimal travel and
quarantine rates designed for the counties of Massachusetts. This showcases
the practical utility of our network-driven strategies.

Paper Organization. The rest of this article is organized as follows. In Sec.
II we review necessary notation and background, and in Sec. III we outline our
problem formulation. Sec. IV details our approach for epidemic control by
optimizing travel rates. In Sec. V we introduce our extension of the SIR model
to consider quarantined individuals, and we demonstrate how epidemic control
is achievable by optimizing quarantine rates. Finally, in Sec. VI we present
numerical results obtained from simulations on networks using publicly
available data from Massachusetts counties.

##  II Notation and Preliminaries

The sets of real, integer, complex, and natural numbers are denoted by
$\mathds{R}$, $\mathds{Z}$, $\mathds{C}$, and $\mathds{N}$, respectively.
$\mathds{R}_{\geq 0}$ ($\mathds{R}_{>0}$) refers to nonnegative (positive)
real numbers, and $\mathds{R}_{\leq 0}$ ($\mathds{R}_{<0}$) refers to
nonpositive (negative) real numbers. Vectors and matrices are denoted by
lowercase and capital letters, respectively. For a matrix
$A\in\mathds{R}^{n\times n}$, the entry at row $i$ and column $j$ is denoted
by $a_{ij}$. The column vector concatenating all columns of $A$ into a single
column is $a=\text{vec}(A)$. Vectors are assumed to be column vectors unless
stated otherwise. $\mathbf{1}$ denotes the all-ones column vector, $I$ the
identity matrix, and $\mathbf{0}$ the all-zero vector or matrix. Inequalities
between vectors are elementwise. A matrix is non-negative (non-positive,
negative, positive) if all its elements are non-negative (non-positive,
negative, positive). A matrix is essentially non-negative if all off-diagonal
elements are non-negative. The dominant eigenvalue and eigenvector of a matrix
$A$ are $\lambda_{\max}(A)$ and $u_{\max}(A)$. For a vector $a$, $\|a\|$ is
its Euclidean norm, and $\mathop{\mathrm{diag}}(a)$ is a diagonal matrix with
$a$ as the diagonal elements. The element-wise multiplication of two vectors
$a$ and $b$ is $a\circ b$. The transpose, spectral radius, inverse, and Moore-
Penrose inverse of a matrix $A$ are denoted by $A^{T}$, $\rho(A)$, $A^{-1}$,
and $A^{+}$. The complex conjugate of $A$ is $\bar{A}$, and its conjugate
transpose is $A^{*}$. The spectral norm of a square matrix $A$ is
$\|A\|_{2}=\lambda_{\max}(A^{*}A)^{\frac{1}{2}}$. The derivative of a variable
$s$ with respect to time $t$ is $\dot{s}$. The gradient of a function $f$ is
$\nabla f$. A function is $L$-smooth if it is continuously differentiable with
a Lipschitz continuous gradient having Lipschitz constant $L$. A neighborhood
$\mathcal{N}(Z_{0})$ of a matrix $Z_{0}\in\mathds{C}^{n\times n}$ is a set of
matrices close to $Z_{0}$ in some metric. For $x\in\mathds{R}$,
$[x]_{+}:=\max\\{x,0\\}$.

The network is represented by a directed graph (digraph), denoted as
$\mathcal{G}=(\mathcal{V},\mathcal{E})$. In $\mathcal{G}$, the set of nodes is
expressed as $\mathcal{V}=\\{v_{1},\ldots,v_{n}\\}$, and the set of edges is
defined as
$\mathcal{E}\subseteq\\{\mathcal{V}\times\mathcal{V}\cup{(v_{i},v_{i})\mid
v_{i}\in\mathcal{V}}\\}$, where each node has a virtual self-edge. The number
of nodes and edges are denoted as $|\mathcal{V}|=n$ and $|\mathcal{E}|=m$,
respectively. A directed path of length $t$ exists from node $v_{i}$ to node
$v_{l}$ if there is a sequence of nodes $v_{i}\equiv
l_{0},l_{1},\dots,l_{t}\equiv v_{l}$ such that
$(l_{\tau},l_{\tau+1})\in\mathcal{E}$ for $\tau=0,1,\dots,t-1$. A digraph is
called “strongly connected” if there is a directed path from every node
$v_{i}$ to every other node $v_{l}$ for all pairs of nodes
$v_{i},v_{l}\in\mathcal{V}$. A square non-negative matrix
$A\in\mathds{R}_{\geq 0}^{n\times n}$ corresponds to a digraph $G(A)$ with $n$
nodes $\\{v_{1},\ldots,v_{n}\\}$ and a set edges
$\mathcal{E}=\\{(v_{i},v_{j})\mid a_{ij}>0\\}$.

######  Definition 1 (Strongly Connected Matrix).

A square non-negative matrix $A$ is said to be strongly connected if its
corresponding digraph $G(A)$ is strongly connected.

######  Definition 2 (Primitive Matrix).

A primitive matrix is a square non-negative matrix some power of which is
positive.

######  Definition 3 (Matrix Stability and Function Convergence).

A matrix $A$ is continuous time stable if all of its eigenvalues possess real
parts that are less than or equal to zero. A matrix $A$ is discrete time
stable if all of its eigenvalues have magnitudes that are less than or equal
to one. A function $y(t)$ exhibits decay at a rate of $\alpha$ starting at
$t_{0}$ if the following condition holds: $y(t)\leq y(t_{0})e^{-\alpha t}$ for
all $t\geq t_{0}$.

###  II-A Network spread model of COVID-19

To model the spread of COVID-19, we use a nuanced extension of the basic
Susceptible-Infected-Recovered (SIR) compartmental model. We separate the
infected individuals into symptomatic and asymptomatic compartments (see [16,
17, 14] and references therein). Dividing the infected population into these
subgroups allows us to better capture their epidemiological dynamics, such as
differing infectiousness and detection rates between asymptomatic and
symptomatic cases. The proposed model can be described as:

$\displaystyle\dot{s}_{i}$
$\displaystyle=-s_{i}\sum_{j=1}^{n}a_{ij}(\beta^{\mathrm{a}}x_{j}^{\mathrm{a}}+\beta^{\mathrm{s}}x_{j}^{\mathrm{s}}),$
(1a) $\displaystyle\dot{x}_{i}^{\mathrm{a}}$
$\displaystyle=s_{i}\sum_{j=1}^{n}a_{ij}(\beta^{\mathrm{a}}x_{j}^{\mathrm{a}}+\beta^{\mathrm{s}}x_{j}^{\mathrm{s}})-(\epsilon+r^{\mathrm{a}})x_{i}^{\mathrm{a}},$
(1b) $\displaystyle\dot{x}_{i}^{\mathrm{s}}$ $\displaystyle=\epsilon
x_{i}^{\mathrm{a}}-r^{\mathrm{s}}x_{i}^{\mathrm{s}},$ (1c)
$\displaystyle\dot{h}_{i}$
$\displaystyle=r^{\mathrm{a}}x_{i}^{\mathrm{a}}+r^{\mathrm{s}}x_{i}^{\mathrm{s}},$
(1d)

where $n$ is the number of locations in our network,
$s=(s_{1},\ldots,s_{n})\in\mathds{R}^{n}$ is a vector with elements $s_{i}$
equal to the proportion of the population at location $i$ being susceptible,
$x^{\mathrm{a}}$ and $x^{\mathrm{s}}$ are vectors with elements equal to the
proportion of the population at each location being asymptomatic and
symptomatic, respectively, and $h$ is the vector with elements $h_{i}$
denoting the portion of recovered population at node $i$. $a_{ij}$ is the rate
at which infections at location $j$ affect those at location $i$ (to be
computed later), $\beta^{\mathrm{a}}$ and $\beta^{\mathrm{s}}$ are the disease
transmission rates of asymptomatic and symptomatic individuals, respectively,
$r^{\mathrm{a}}$ and $r^{\mathrm{s}}$ are the recovery rates of asymptomatic
and symptomatic, respectively, and $\epsilon$ is the rate at which
asymptomatic individuals develop symptoms.

Following [14], different parameters are utilized for symptomatic and
asymptomatic individuals due to the findings of [18]. Specifically, in [18] it
is shown that asymptomatic individuals exhibit a more rapid decline in viral
load (compared to symptomatic). Consequently, not only do they experience
quicker recovery, but they are also likely to be less contagious.

######  Remark 1.

Note that we consider the model in (1) due to its versatile framework. While
(1) can be reduced to the SIR or SEIR model (by setting
$\beta^{\mathrm{s}}=r^{\mathrm{s}}=\epsilon=0$, or
$\beta^{\mathrm{a}}=r^{\mathrm{a}}=0$, respectively), its main advantage lies
in its extensibility. As detailed in the following sections, this
extensibility is one of the main research directions of this paper.

###  II-B Construction of Infection Flow and Travel Rates Matrices

Infection Flow Matrix. To construct the infection flow matrix $A$, we adopt
the approach of [17], which is particularly suitable for modeling the impact
of lockdown measures for controlling COVID-19. At each location $i$, we denote
the fixed population size by $N_{i}$. Additionally, we assume that people
travel from location $i$ to location $j$ at a rate denoted by $\tau_{ij}$
(determining $\tau_{ij}$ is addressed later in the paper). Thus, the travel
rate matrix is an $n\times n$ matrix defined as $T=[\tau_{ij}]_{i,j=n}^{n}$.
The focus on travel rates is essential because they play a pivotal role in
determining the evolution and spread of an epidemic. For instance, [19]
revealed a strong correlation between the regional progression of influenza
and people’s movements.

To construct $A$, let us track how the susceptible population moves to other
compartments over time. The flow of the susceptible population from location
$i$ to location $l$ is $s_{i}\tau_{il}$. The rate of infection at location $l$
is proportional to the fraction of total (asymptomatic and symptomatic)
infected people traveling to location $l$
($\sum_{j=1}^{n}N_{j}\tau_{jl}x_{j}^{a}$ and
$\sum_{j=1}^{n}N_{j}\tau_{jl}x_{j}^{s}$, respectively) over the total
population traveling to location $l$ ($\sum_{k=1}^{n}N_{k}\tau_{kl})$. People
at location $i$ may get infected as they travel to location $l$, which leads
to the following rate of change to the susceptible individuals at location
$i$:

$\dot{s}_{i}=-\sum_{l=1}^{n}s_{i}\tau_{il}\left(\frac{\sum_{j=1}^{n}N_{j}\tau_{jl}x_{j}^{a}}{\sum_{k=1}^{n}N_{k}\tau_{kl}}\beta^{a}+\frac{\sum_{j=1}^{n}N_{j}\tau_{jl}x_{j}^{s}}{\sum_{k=1}^{n}N_{k}\tau_{kl}}\beta^{s}\right).$
(2)

Thus, comparing (2) with (1a), we can define

$a_{ij}=\sum_{l=1}^{n}\tau_{il}\tau_{jl}\frac{N_{j}}{\sum_{k=1}^{n}N_{k}\tau_{kl}}.$
(3)

Travel Rates. To construct the travel rate matrix $T$, we rely on the Human
Mobility Flow dataset [20]. This dataset captures daily visitor flows between
Census Block Groups (CBGs) during the COVID-19 pandemic. It provides insights
into how people move between locations, measured as the daily number of CBG-
to-CBG flows for trips lasting more than $1$ minute. It also generates the
estimates for the whole population flow that reflects the general human
mobility patterns at the CBG level. In the simulations, we focus on county-
level trips within Massachusetts. We aggregate the population flow counts for
all nodes, including flows within each county and to other counties. To
generate $T$, we assume people in node $i$ spend an $t_{i}$ ratio of their day
outside of their residences and therefore,

$\tau_{ij}=t_{i}P_{f}(i,j)/\sum_{k}{P_{f}(i,k)},$

where $P_{f}$ is a population flow matrix, with $P_{f}(i,j)$ denoting the
number of daily trips (lasting in location $j$ more than $1$ minute) from
location $i$ to location $j$.

###  II-C Stability and Splitting of Matrices

We now provide the following lemmas about the stability of matrices that is
important for our subsequent development.

######  Lemma 1 (A Perron–Frobenius version from [14]).

Suppose $A$ is a matrix that is strongly connected, with nonnegative off-
diagonal elements. In this case, there exists a real eigenvalue of $A$ that is
greater in magnitude than the real part of any other eigenvalue of $A$. This
dominant eigenvalue $\lambda_{\max}(A)$ is unique (simple) and its
corresponding eigenvector $u_{\max}(A)$ is real and positive.

######  Lemma 2 ( [14]).

The following statements characterize continuous and discrete-time matrix
stability.

  * •

(Continuous-time stability). A strongly connected matrix $P$ with non-negative
off-diagonal elements is continuous-time stable if and only if there exists a
vector $d>{\mathbf{0}}$ such that $Pd\leq{\mathbf{0}}$.

  * •

(Discrete-time stability). A non-negative strongly connected matrix $B$ is
discrete-time stable if and only if there exists a vector $d>{\mathbf{0}}$
such that $Bd\leq d$.

  * •

(Connection between Continuous-time and Discrete-time stability). Suppose
$P=L-D$ where $L$ is a non-negative matrix while $D$ is a matrix with non-
positive off-diagonal elements whose inverse is elementwise non-negative.
Suppose further that both $P$ and $D^{-1}L$ are strongly connected. Then $P$
is continuous-time stable if and only if $B=D^{-1}L$ is discrete-time stable.

Lemma 1 presents a version of the Perron-Frobenius theorem for strongly
connected matrices with non-negative off-diagonal elements. It establishes the
existence of a dominant real eigenvalue, that is unique and simple, with a
positive eigenvector. Lemma 2 characterizes matrix stability for both a
continuous-time and a discrete-time setting, establishing relationships
between matrix properties and stability conditions.

###  II-D Balancing of Nonnegative Matrices

Matrix balancing is a mathematical problem that involves adjusting the
elements of a given matrix to meet certain criteria or constraints. In this
context, a square matrix $A^{n\times n}\in\mathds{R}^{n\times n}_{\geq 0}$ is
considered balanced if $A\mathbf{1}=A^{T}\mathbf{1}$. This means that for
every $i=1,...,n$, the sum of the entries in row $i$ of $A$ is equal to the
sum of the entries in column $i$. The problem of balancing matrix $A$ is
equivalent to finding a non-negative diagonal matrix $D$ such that $DAD^{-1}$
is balanced.

Early works [21] consider matrix balancing as the condition for a sequence of
matrices to convergence to a doubly stochastic limit. More recent works [22]
aimed to balance a given matrix in a distributed way by considering that each
row(/column) is the set of incoming(/outgoing) links of the corresponding node
in a given network (see [23] for a survey on matrix balancing). Here, we rely
on the approach presented in [24] in which the authors showed that matrix
balancing can be achieved in linear time. Following [24], we will summarize
the complexity by stating that it is proportional to the number of nonzero
entries of the matrix (as discussed in [14, Sec. $2.2$]).

##  III Problem Formulation

Our goal is to optimally modify the elements of the COVID-19 model (cf.Sec. 1)
in a non-uniform way to mitigate the spread of the disease while minimizing
economic costs. We focus on two problems, denoted by P1 and P2:

  * P1

In the first problem, we aim to optimally modify the travel rate matrix $T$
(e.g., by placing limited travel restrictions) in order to curb the spread of
infections (see Section IV).

  * P2

In the second problem, our objective is to extend the COVID-19 model (cf. (1))
to encompass quarantined individuals. Subsequently, we aim to optimally select
location-dependent quarantine rates in order to curb the spread of infections
in the network while simultaneously minimizing the associated economic cost
(see Section V).

##  IV Optimal Travel Rates for Epidemic Control

###  IV-A Problem Structure

The network spread model of COVID-19 in (1) can be written in matrix form as

$\begin{pmatrix}\dot{s}\\\ \dot{x^{\mathrm{a}}}\\\ \dot{x^{\mathrm{s}}}\\\
\dot{h}\\\
\end{pmatrix}=\begin{pmatrix}{\mathbf{0}}&-\beta^{\mathrm{a}}\mathop{\mathrm{diag}}{(s)}A&-\beta^{\mathrm{s}}\mathop{\mathrm{diag}}{(s)}A&{\mathbf{0}}\\\
{\mathbf{0}}&C&\beta^{\mathrm{s}}\mathop{\mathrm{diag}}{(s)}A&{\mathbf{0}}\\\
{\mathbf{0}}&\epsilon I&-r^{\mathrm{s}}I&{\mathbf{0}}\\\
{\mathbf{0}}&r^{\mathrm{a}}I&r^{\mathrm{s}}I&{\mathbf{0}}\\\
\end{pmatrix}\begin{pmatrix}s\\\ x^{\mathrm{a}}\\\ x^{\mathrm{s}}\\\ h\\\
\end{pmatrix},$ (4)

where
$C=\beta^{\mathrm{a}}\mathop{\mathrm{diag}}{(s)}A-(\epsilon+r^{\mathrm{a}})I$.
Let us consider matrix $A$ as a function of $\tau=\text{vec}(T)$ (where $T$ is
the travel matrix defined in Section II-B). For (4) we have that [11,
Assumption $1$] and [11, Assumption $2$] hold. For this reason, in (4) we can
decouple the dynamics of $\dot{x}$ from $\dot{s}$ and $\dot{h}$, and focus on
them. Therefore, for our analysis we focus on optimizing the travel rates in
matrix $M(t_{0},\tau)$ defined as

$M(t_{0},\tau)=\begin{pmatrix}C^{\prime}&\beta^{\mathrm{s}}\mathop{\mathrm{diag}}{(s(t_{0}))}A(\tau)\\\
\epsilon I&-r^{\mathrm{s}}I\\\ \end{pmatrix},$ (5)

where
$C^{\prime}=\beta^{\mathrm{a}}\mathop{\mathrm{diag}}{(s(t_{0}))}A(\tau)-(\epsilon+r^{\mathrm{a}})I$.
For a fixed $t_{0}$ we define
$f(\tau)=\lambda_{\max}\left(M(t_{0},\tau)\right)$, and we aim to solve the
following minimization problem:

$\begin{split}&{\min_{\tau}}\quad f(\tau)\\\
&\text{s.t.}\quad\|\tau-\tau_{0}\|_{1}\leq b,\\\
&\,\,\qquad\tau\geq{\mathbf{0}}\,\end{split}$ (6)

where $\tau_{0}$ is the vector of the initial travel rates and $b$ is a budget
on the amount of travel rate change (measured using an $\ell_{1}$ norm) from
their initial values. Note that $b$ is important for preventing drastic
changes over the travel rates which would result in considerable social
disruption. This constraint is essential for addressing practical
considerations related to the feasibility of changing the travel rates.

######  Remark 2.

In (6), we focus on minimizing the dominant (or maximum) eigenvalue of the
matrix $M(t_{0},\tau)$ by optimizing the travel rates $\tau$. The dominant
eigenvalue plays a crucial role in characterizing the rate of infection spread
within the COVID-19 model described in (4). By minimizing the dominant
eigenvalue, we effectively target the control of infection spread dynamics.
Additionally, by minimizing the dominant eigenvalue we can achieve a specific
decay rate of infected cases. For example, when the maximum eigenvalue is less
than a specified threshold $-\alpha$ at time $t_{0}$, the infected cases will
decay at a rate of at least $e^{-\alpha t}$ for times $t\geq t_{0}$.

###  IV-B Proposed Solution

We apply Projected Gradient Descent (PGD) by projecting onto the constraint
set, specifically the non-negative orthant of the $l_{1}$ ball. To ensure
convergence of the PGD, we employ a suitable stepsize choice using the
backtracking line search technique. This adaptive strategy iteratively adjusts
the stepsize based on sufficient decrease conditions, enhancing the robustness
and convergence of the optimization process (see [25]). Considering these
characteristics, our proposed algorithm is as follows:

$\begin{split}&\textit{Initialization:}\\\ &\text{set}\ \gamma=1\ \text{and
choose}\ \beta\in(0,1),\\\ &\textit{Iteration:}\\\ &\text{while}\
f(\tau_{k}-\gamma\nabla f(\tau_{k}))>f(\tau_{k})-\frac{\gamma}{2}\|\nabla
f(\tau_{k})\|^{2},\\\ &\ \ \ \text{set}\ \gamma\leftarrow\beta\gamma,\\\
&y_{k+1}=\tau_{k}-\gamma\nabla{f(\tau_{k}}),\\\
&\tau_{k+1}=\underset{\tau\in\Omega}{\arg\min}{\|\tau-
y_{k+1}\|}_{2}^{2},\end{split}$ (7)

where $\gamma$ is the step size adjusted via backtracking line search, $\beta$
is the scaling parameter of $\gamma$, $\tau_{k}$ is in the vector form (i.e.,
$\tau_{k}=\text{vec}(T_{k})$) of the travel rates matrix $T_{k}$ at time step
$k$, and $\Omega$ is the projection set.

Focusing on (7), the update involves a gradient descent step followed by a
projection onto the feasible set $\Omega$. Each step presents unique
challenges. Computing the gradient of the maximum eigenvalue of a matrix is a
complex task that requires specialized techniques, which we will detail in our
subsequent analysis. For the projection step, it is essential to consider
appropriate constraints and transformations to ensure the optimal travel rates
are practical and realistic.

Gradient Descent Step. Computing the gradient of a matrix maximum eigenvalue
may present challenges [26]. To address them, we leverage [27, Thm. $2$],
which we state for completeness.

######  Theorem 1 ( [27]).

Let $\lambda_{0}$ be a simple eigenvalue of a matrix
$Z_{0}\in\mathds{C}^{n\times n}$ and let $u_{0}$ be an associated eigenvector,
so that $Z_{0}u_{0}=\lambda_{0}u_{0}$. Then a (complex) function $\lambda$ and
a (complex) vector function $u$ are defined for all $Z$ in some neighborhood
$\mathcal{N}(Z_{0})\in\mathds{C}^{n\times n}$ of $Z_{0}$, such that

$\lambda(Z_{0})=\lambda_{0},\ \ u(Z_{0})=u_{0}\ ,$

and

$Zu=\lambda u,\ u_{0}^{*}u=1,\ \text{for}\ Z\in\mathcal{N}(Z_{0}).$

Moreover, the functions $\lambda$ and $u$ are $\infty$ times differentiable on
$\mathcal{N}(Z_{0})$ and the differentials at $Z_{0}$ are

$d\lambda=v_{0}^{*}(dZ)u_{0}/v_{0}^{*}u_{0}\ ,$ (8)

and

$du=(\lambda_{0}I-Z_{0})^{+}\left(I-\frac{u_{0}v_{0}^{*}}{v_{0}^{*}u_{0}}\right)(dZ)u_{0}\
,$ (9)

where $v_{0}$ is the eigenvector associated with the eigenvalue
$\bar{\lambda}_{0}$ of $Z_{0}^{*}$, so that
$Z_{0}^{*}v_{0}=\bar{\lambda}_{0}v_{0}$, $d\lambda$ is the differential of the
(complex) function $\lambda$, and $dZ$ is the differential of $Z$.

Let $Z_{0}=M(t_{0},\tau_{k})$ in Theorem 1, with $\lambda_{\max}$ and
$u_{\max}$ denoting the maximum eigenvalue and corresponding eigenvector of
$M(t_{0},\tau_{k})$, so $M(t_{0},\tau_{k})u_{\max}=\lambda_{\max}u_{\max}$.
According to Theorem 1, in a neighborhood
$\mathcal{N}(M(t_{0},\tau_{k}))\in\mathds{C}^{n\times n}$, there exists a
function $f$ such that $f(M(t_{0},\tau_{k}))=\lambda_{\max}$ and

$\frac{\partial
f(M(t_{0},\tau_{k}))}{\partial\tau_{ij}}=\frac{v^{*}_{\max}\left(\frac{\partial
M(t_{0},\tau_{k})}{\partial\tau_{ij}}\right)u_{\max}}{v^{*}_{\max}u_{\max}},$
(10)

where $v_{\max}$ is the eigenvector associated with $\bar{\lambda}_{\max}$ of
$M(t_{0},\tau_{k})^{*}$, satisfying
$M(t_{0},\tau_{k})^{*}v_{\max}=\bar{\lambda}_{\max}v_{\max}$. Since
$M(t_{0},\tau_{k})$ is real, ${M(t_{0},\tau_{k})}^{*}={M(t_{0},\tau_{k})}^{T}$
and ${M(t_{0},\tau_{k})}^{T}v_{\max}=\lambda_{\max}v_{\max}$. By the
Perron–Frobenius theorem (Lemma 1), the dominant eigenvector is real and
positive, thus $v^{*}_{\max}=v^{T}_{\max}$. Note that $v^{T}_{\max}$ and
$u_{\max}$ are the left and right eigenvectors corresponding to
$\lambda_{\max}$ of $M(t_{0},\tau_{k})$. Consequently, (10) becomes

$\frac{\partial
f(M(t_{0},\tau_{k}))}{\partial\tau_{ij}}=\frac{v^{T}_{\max}\left(\frac{\partial
M(t_{0},\tau_{k})}{\partial\tau_{ij}}\right)u_{\max}}{v^{T}_{\max}u_{\max}},$
(11)

providing the components of the gradient vector $\nabla f(M(t_{0},\tau_{k}))$.
Using (11), we perform the gradient descent step specified in (7).

Projection Step. We now focus on the projection step of (7). This step
involves minimizing $\|\tau-y_{k+1}\|_{2}^{2}$ over $\tau$ subject to the
constraint set of (6), resulting in linear constraints. By converting this
problem into a convex Quadratic Programming (QP) problem [28], the optimal
$\tau$ can be efficiently found at each step.

Convergence Analysis. For analyzing the convergence of the proposed algorithm
to a stationary point and the operation of (7), we present the following
lemmas which are necessary for the subsequent development.

Note that the spectral norm of the Hessian matrix $\nabla^{2}f(\tau)$ is
$|\lambda_{\max}(\nabla^{2}f(\tau))|$ because the Hessian matrix of $f(\tau)$
is symmetric and real.

######  Lemma 3.

The eigenvalues of the matrix $M(t_{0},\tau)$ (cf. (5)) are continuous
functions of $\tau$.

###### Proof.

The proof is a direct result of the continuity of the roots of the
characteristic polynomial of $M(t_{0},\tau)$ (see [29]). ∎

######  Lemma 4.

The function $f(\tau)=\lambda_{\max}\left(M(t_{0},\tau)\right)$ is $L$-smooth.

###### Proof.

Before proving the lemma, note that the set consisting of all $M(t_{0},\tau)$
matrices over all the possible travel rate matrices $T$ (defined in Section
II-B) is bounded and compact.

Let us consider the function $f(\tau)$ representing the dominant eigenvalue of
$M(t_{0},\tau)$. From Thm. 1, we have that $f(\tau)$ is $\infty$ times
differentiable. Therefore, every element of the Hessian of $f(\tau)$ is
continuous and differentiable. From the continuity of the Hessian and
continuity of the maximum eigenvalue (see Thm. 1) we have that the largest
element of the Hessian achieves its maximum value on the compact set
consisting of all $M(t_{0},\tau)$ over all the possible travel rate matrices
$T$. As a result, the spectral norm of the Hessian of $f(\tau)$
($\|\nabla^{2}f(x)\|_{2}$) is bounded and from [30] the maximum eigenvalue
function $f(\tau)$ has a Lipschitz continuous gradient. ∎

After showing that the dominant eigenvalue function $f(\tau)$ is $L$-smooth,
we are now ready to present the following convergence result for the PGD
operation with backtracking line search described in (7).

######  Theorem 2.

Let us consider the function
$f(\tau_{k})=\lambda_{\max}\left(M(t_{0},\tau_{k})\right)$. Executing the
operation described in (7), we have

$\lim_{k\rightarrow\infty}\nabla(f(\tau_{k}))=0.$ (12)

###### Proof.

Let $\\{\tau_{k}\\}$ be the sequence of points generated by the operation
described in (7). According to the convergence of the backtracking line search
for a differentiable and $L$-smooth function (see [25]), it follows that
$\lim_{k\rightarrow\infty}\nabla(f(\tau_{k}))=0$. ∎

##  V COVID-19 Model with Optimal Quarantine Rates

In this section, our goal is to enhance the COVID-19 model in (1) by
incorporating the concept of quarantine for infected individuals. Assigning an
economic cost to each quarantined individual, we will focus on optimization
strategies to reduce the quarantined population effectively, with the dual
objective of containing the infection spread within our network.

###  V-A Integrating Quarantining into a COVID-19 Model

In order to extend the COVID-19 model, we implement the following
modifications:

  1. 1.

Asymptomatic infected individuals $x_{i}^{\mathrm{a}}$ at each node $i$ $(a)$
exhibit automatic recovery at a rate of $r^{\mathrm{a}}$, or $(b)$ develop
symptoms at a rate of $\epsilon$, or $(c)$ are quarantined with rate
$q_{i}^{\mathrm{a}}$, with $k_{i}$ denoting the proportion of the quarantined
population at node $i$.

  2. 2.

Symptomatic cases $x_{i}^{\mathrm{s}}$ $(a)$ experience automatic recovery at
a rate of $r^{\mathrm{s}}$, or $(b)$ transition to a quarantined state at a
rate of $q_{i}^{\mathrm{s}}$.

  3. 3.

Quarantined individuals $k_{i}$ recover at a rate of $r^{q}$.

Considering the above modifications, the extended COVID-19 model, called the
Susceptible-Infected-Quarantined-Recovered (SIQR) model, is described as
follows:

$\displaystyle\dot{s}_{i}$
$\displaystyle=-s_{i}\sum_{j=1}^{n}a_{ij}(\beta^{\mathrm{a}}x_{j}^{\mathrm{a}}+\beta^{\mathrm{s}}x_{j}^{\mathrm{s}}),$
(13a) $\displaystyle\dot{x}_{i}^{\mathrm{a}}$
$\displaystyle=s_{i}\sum_{j=1}^{n}a_{ij}(\beta^{\mathrm{a}}x_{j}^{\mathrm{a}}+\beta^{\mathrm{s}}x_{j}^{\mathrm{s}})-(\epsilon+r^{\mathrm{a}}+q_{i}^{\mathrm{a}})x_{i}^{\mathrm{a}}\
,$ (13b) $\displaystyle\dot{x}_{i}^{\mathrm{s}}$ $\displaystyle=\epsilon
x_{i}^{\mathrm{a}}-(r^{\mathrm{s}}+q_{i}^{\mathrm{s}})x_{i}^{\mathrm{s}},$
(13c) $\displaystyle\dot{k}_{i}$
$\displaystyle=q_{i}^{\mathrm{a}}x_{i}^{\mathrm{a}}+q_{i}^{\mathrm{s}}x_{i}^{\mathrm{s}}-r^{q}k_{i},$
(13d) $\displaystyle\dot{h}_{i}$
$\displaystyle=r^{\mathrm{a}}x_{i}^{\mathrm{a}}+r^{\mathrm{s}}x_{i}^{\mathrm{s}}+r^{q}k_{i},$
(13e)

The proposed SIQR model uses the same notation as the model in (1) and only
differs in the incorporation of quarantining. Quarantined individuals come
from both asymptomatic and symptomatic cases and recover at a rate of $r^{q}$.

###  V-B Problem Formulation for Minimizing Quarantine Rates

Let us define the column vectors
$q^{\mathrm{a}}=(q^{\mathrm{a}}_{1},\ldots,q^{\mathrm{a}}_{n})$, and
$q^{\mathrm{s}}=(q^{\mathrm{s}}_{1},\ldots,q^{\mathrm{s}}_{n})$, and the
column vector of quarantine rates $q=(q^{\mathrm{a}},q^{\mathrm{s}})$. We
decouple the dynamics of $\dot{x}$ from the other terms in (13) for a disease-
free equilibrium and focus on analyzing the $2n\times 2n$ sub-matrix $M(t,q)$
defined as:

$M(t,q)=\begin{pmatrix}E&\beta^{\mathrm{s}}\mathop{\mathrm{diag}}{(s(t))}A\\\
\epsilon I&-(r^{\mathrm{s}}I+\mathop{\mathrm{diag}}{(q^{\mathrm{s}}}))\\\
\end{pmatrix},$ (14)

where
$E=\beta^{\mathrm{a}}\mathop{\mathrm{diag}}{(s(t))}A-(\epsilon+r^{\mathrm{a}})I-\mathop{\mathrm{diag}}{(q^{\mathrm{a}}})$.
Our objective is to optimally minimize quarantine rates, thereby curbing the
spread of infection within our network, while simultaneously minimizing the
associated economic costs. The quarantine rates $q_{i}^{s},q_{i}^{a}$ can be
assumed to be nonnegative. We note also that there is an upper bound on what
these rates could plausibly be, since the inverse of these rates equals the
amount of time a typical individual would spend before being quarantined under
a scenario where individuals do not move across classes such as symptomatic,
asymptomatic, recovered, etc. Consequently, we may assume without loss of
generality a bound of $1$ on all the quarantine rates. To achieve this, we
formulate the following optimization problem:

$\begin{split}&{\min}_{q}\quad\sum_{i=1}^{n}\frac{z_{i}^{a}}{1-q_{i}^{\mathrm{a}}}+\frac{z_{i}^{s}}{1-q_{i}^{\mathrm{s}}}\\\
&\text{s.t.}\quad\lambda_{\max}(M(t_{0},q))\leq-\alpha,\\\
&\qquad{\mathbf{0}}\leq q\leq{\mathbf{1}},\end{split}$ (15)

where $z^{a}=(z^{a}_{1},\ldots,z^{a}_{n})$ and
$z^{s}=(z^{s}_{1},\ldots,z^{s}_{n})$ represent the relative economic costs
associated with quarantining asymptomatic and symptomatic cases at each node,
respectively. The objective function in (15) is designed to minimize the
overall economic cost associated with quarantine rates. The economic cost for
each node is proportional to the effectiveness of quarantine (i.e., in (15),
increasing quarantine rates leads to higher costs) and reflects a nonlinear
increase in costs that seeks to model the increased complexity and logistical
costs associated with quarantining a large proportion of the population. The
summation across all nodes aggregates the economic costs, and the optimization
problem aims to find the quarantine rates $q$ that minimize this combined
economic cost. The constraints ensure that the maximum eigenvalue of the
associated sub-matrix $M(t_{0},q)$ is bounded by $-\alpha$. This reflects a
condition for controlling the infection spread. Therefore, our objective is to
minimize the economic costs associated with quarantining individuals while
simultaneously achieving a decay rate of $-\alpha$ for the infected cases.

###  V-C Proposed Solution

The choice of optimization method for solving problem (15) differs from our
previous approach (i.e., PGD for optimizing the travel rates) due to the
nature of the constraints involved. Specifically, unlike the linear
constraints in problem (6), the optimization problem in (15) incorporates
nonlinear constraints on quarantine rates. Therefore, in this case, we will
utilize an alternative approach.

Reduction of the Optimization Problem to Weight-Balancing. We are now ready to
show that Problem (15) can be reduced to a matrix balancing problem. To this
end, we define the following assumption, acknowledging that this assumption is
typically fulfilled by the epidemic model parameters we use. Specifically,
when fitting the models to data from Massachusetts using the existing model
parameters in the literature, it holds. In fact, the matrix balancing
algorithm provided the correct answer to problem (15).

######  Assumption 1.

Consider the matrix $B_{0}$ defined as:

$B_{0}=\begin{pmatrix}\beta^{\mathrm{a}}\operatorname{diag}(s(t_{0}))A-(\epsilon+r^{\mathrm{a}}+1)I&\beta^{\mathrm{s}}\operatorname{diag}(s(t_{0}))A\\\
\epsilon I&-(r^{\mathrm{s}}+1)I\end{pmatrix}$

and let $m$ be the largest absolute value of any diagonal element of $B_{0}$.
Then, defining
$x=\min_{j}\left(\epsilon\beta^{\mathrm{s}}s_{j}(t_{0})A_{jj}\right)$, we have

$1+\frac{x}{m^{2}}\geq m\ .$

Assumption 1 is needed to ensure the solution is non-negative. Now, we present
the main theorem of our study:

######  Theorem 3.

The minimum quarantining problem in (15) can be reduced to a matrix balancing
problem provided that Assumption 1 is satisfied, the infection flow matrix $A$
(from Sec. II-B) is strongly connected, $s(t_{0})>0$, and the problem is
feasible.

###### Proof.

See Appendix A. ∎

Inspired from [14], we reduced problem (15) to a matrix balancing problem with
a known polynomial complexity in terms of the number of unknown variables
($2n$). Note that problem (15) is potentially of higher-order complexity due
to the eigenvalue condition in the fist constraint. While the matrix balancing
approach necessitates calculating the inverse of a $2n\times 2n$ matrix
$B_{0}$ (which has a cubic complexity), our proposed solution maintains a
polynomial complexity with respect to the number of parameters ($2n$).

Connection of Constraints with the Epidemic Basic Reproduction Number. In
epidemiological analysis, significant emphasis is placed on determining the
reproduction number $R_{0}$ In our scenario $R_{0}$ is of particular
importance since it can indicate whether a disease outbreak will escalate into
an epidemic or naturally dissipate. Unfortunately, calculating $R_{0}$ for a
specific disease is a non trivial task (primarily due to the absence of
comprehensive databases containing essential data such as infection rates and
recovery rates). Motivated by this, we establish a relationship between the
constraint $\lambda_{\max}(M(t_{0}))\leq-\alpha$ and the epidemic’s
reproduction number $R_{0}$. This connection is crucial, as modifying our
constraints can prevent an epidemic from escalating to a pandemic or spreading
extensively within a population.

######  Theorem 4.

The set $R_{0}\leq 1$ is equivalent with $\lambda_{\max}\leq 0$

###### Proof.

See Appendix B. ∎

Theorem 4 implies that imposing the condition
$\lambda_{\max}(M(t_{0}))\leq-\alpha$ provides a set of solutions that
enforces a bound on the disease reproduction number $R_{0}\leq r$, for
$r\in[0,1]$.

Convergence via Augmented Primal-Dual Gradient Dynamics for Optimal Quarantine
Control. We now extend our previous approach to solving the optimal quarantine
problem by applying the augmented primal-dual gradient dynamics (Aug-PDGD).
This extension demonstrates semi-global exponential convergence to a KKT point
of the problem. Below are the relevant propositions and theorems supporting
this result.

######  Proposition 1.

Assuming the feasibility of (15), the cost function, defined as
$f(q)=\sum_{i=1}^{n}\left(\frac{z_{i}^{a}}{1-q_{i}^{\mathrm{a}}}+\frac{z_{i}^{s}}{1-q_{i}^{\mathrm{s}}}\right)$,
is continuously differentiable and strongly convex on the feasible domain
${\mathbf{0}}\leq q<{\mathbf{1}}$. Moreover, $M(t_{0},q)$ can be expressed as
the sum of an essentially non-negative matrix and a diagonal matrix with the
elements of the vector $-q$ on its diagonal. Therefore, as demonstrated in
[31], the dominant eigenvalue of the matrix $M(t_{0},q)$ is convex. Defining
$g_{1}(q)=\lambda_{\max}(M(t_{0},q))+\alpha$, $g_{(1+i)}=-q_{i}$, and
$g_{2n+1+i}=q_{i}-1$ for all $i\in\\{1,\dots,2n\\}$, $g(q)$ is continuously
differentiable and convex.

######  Theorem 5.

Given $g_{1}(q)=\lambda_{\max}(M(t_{0},q))+\alpha$, (11) and (14) imply that

$\nabla g_{1}(q)=\frac{v_{\max}\circ u_{\max}}{v_{\max}^{T}u_{\max}}.$

Recall that $v_{\max}^{T}$ and $u_{\max}$ represent the left and right
eigenvectors associated with $\lambda_{\max}(M(t_{0},q))$. Hence, $\|\nabla
g_{1}(q)\|\leq L_{g,1}$ (see Appendix C for proof), and $g_{1}(q)$ is also
$M_{g,1}$-smooth (see Lemma 4). Additionally, for all $i\geq 2$, the gradients
$\nabla g_{i}$ are Lipschitz continuous and have bounded norms, which is
evident from their definitions.

To analyze the augmented primal-dual gradient dynamics, we introduce the
augmented Lagrangian of (15) as [32]

$L_{\rho}(q,\lambda)=f(q)+\Theta_{\rho}(q,\lambda),\qquad\lambda\geq{\mathbf{0}},$
(16)

where

$\Theta_{\rho}(q,\lambda):=\sum_{i=1}^{4n+1}\frac{[\rho
g_{i}(q)+\lambda_{i}]_{+}^{2}-\lambda_{i}^{2}}{2\rho}.$

The augmented primal-dual gradient dynamics are given by

$\displaystyle\dot{q}(t)=-\nabla f(q(t))-\sum_{i=1}^{4n+1}[\rho
g_{i}(q(t))+\lambda_{i}(t)]_{+}\nabla g_{i}(q(t))$ (17a)
$\displaystyle\dot{\lambda}(t)=\sum_{i=1}^{4n+1}\frac{[\rho
g_{i}(q)+\lambda_{i}(t)]_{+}-\lambda_{i}(t)}{\rho}e_{i}.$ (17b)

######  Proposition 2 ( [32]).

Denoting $y(t)=(q(t),\lambda(t))$, $t\geq 0$ as a differentiable trajectory
satisfying (17). Then, for a KKT point of (16), $y^{*}=(q^{*},\lambda^{*})$,
we have $\|y(t)-y^{*}\|\leq\|y(0)-y^{*}\|$, for all $t\geq 0$.

######  Theorem 6.

The gradient of the cost function $f(q)$ is Lipschitz continuous over the
region traversed by the primal gradient dynamics described by (17).

###### Proof.

See Appendix D. ∎

######  Proposition 3.

Let the active set at the local optimal solution $q^{*}$ be denoted by
$\mathcal{I}:=\\{i:g_{i}(q^{*})=0\\}$. Therefore, $q^{*}$ is a regular local
minimum, i.e., $\nabla g_{i}(q^{*})$ for $i\in\mathcal{I}$ are linearly
independent.

###### Proof.

See Appendix E. ∎

######  Corollary 1.

Under Propositions 1 and 3, and Theorems 5 and 6, the application of augmented
primal-dual gradient dynamics (Aug-PDGD) to (15) achieves semi-global
exponential stability. Specifically, the KKT point $y^{*}=(q^{*},\lambda^{*})$
is a semi-globally exponentially stable equilibrium of the Aug-PDGD.

This result directly follows immediately from Theorem 6 and Proposition 3 and
the results of [32], which examines the stability of Aug-PDGD for smooth
convex optimization problems with general convex and nonlinear inequality
constraints. Consequently, the distance to the optimal solution decays
exponentially from any initial point, although the convergence rate may depend
on the initial distance to a stationary point.

It should be noted that [32] assumes the Lipschitz continuity of $\nabla
f(q)$. However, as stated in Theorem 6, what is actually required is the
Lipschitz continuity of $\nabla f(q)$ over the region invariant under the
primal gradient dynamics.

##  VI Simulation Results

In this section we present simulation results for our optimization strategies.
Specifically, for both problems P1 and P2 we focus on the scenario of a $14$
node network that represents the counties of the state of Massachusetts. We
aim to demonstrate the advantages of our approaches and emphasize their
applicability in real world scenarios.

###  VI-A Model Parameters

Massachusetts 14 Node Network. The network comprises 14 nodes, each
representing a county in Massachusetts. The population of each county, $N_{i}$
for $i\in\\{1,...,14\\}$, is based on the 2020 Census data [33]. Travel rates
$\tau$ for the infection flow matrix $A$ are derived from the Human Mobility
Flow dataset [20] using the methodology in Section II-B. We set
$t_{i}=\frac{1}{3}$ for all $i$, assuming people spend one-third of their time
outside.

Epidemic and Economic Parameters. GDP data from the Bureau of Economic
Analysis [34] is used to assess economic impact on counties. Economic loss
($z_{i}^{\mathrm{a}}$ and $z_{i}^{\mathrm{s}}$) is calculated using $g_{i}$, a
county’s GDP, as $z_{i}^{\mathrm{a}}=z_{i}^{\mathrm{s}}=g_{i}/g_{\max}$, where
$g_{\max}$ is the highest GDP among all counties.

COVID-19 spread is modeled using parameters from [17, 16]. The recovery rate
$\gamma$ is set uniformly for all disease states
($r^{a}=r^{s}=r^{q}=\gamma=0.2$). The rate of developing symptoms $\epsilon$
is 0.32. Transmission rates $\beta^{a}$ and $\beta^{s}$ follow [14] with
$\beta^{a}=\eta\beta^{s}$ and $\eta=0.6754$ [16]. $\beta^{s}$ is adjusted to
match the observed initial growth rate, and $\beta^{a}$ is calculated based on
$\eta$. Adjustments to $\beta$ ensure the infection flow matrix $A$ meets the
required growth rate.

Initial Rates. Initial susceptible rates $s(t_{0})$ are set using county-level
data on cumulative cases and deaths [35], adjusted for unreported cases by
dividing cumulative cases by 0.14 [17, 36]. The initial susceptible rate of
node $i$ is $1-\frac{c_{i}(t_{0})}{0.14N_{i}}$, where $c_{i}(t_{0})$ is the
cumulative cases for node $i$. Initial state $t_{0}$ is set to April 1, 2020.

Initial recovered cases are calculated by multiplying the total recovered
ratio of US cumulative cases on April 1, 2020 $\frac{8878}{215215}$ [14] by
each node’s cumulative infected population and adding deaths. Active cases,
both asymptomatic and symptomatic, are calculated as the remainder of
cumulative cases. We assume 14% are symptomatic and 86% are asymptomatic, with
all rates adjusted by 0.14 to account for reporting inaccuracies.

###  VI-B Massachusetts County-Level $14$ Node Network

Optimizing Travel Rates. We now validate the correctness of our optimal
solution for problem (6) over the $14$ node network that represents the state
of Massachusetts. As we mentioned in Section VI-A, for constructing the travel
rate matrix $\tau$ and the infection flow matrix $A$, we utilized the Human
Mobility Flow dataset [20].

In Fig. 1 we observe the value of
$f(\tau)=\lambda_{\max}\left(M(t_{0},\tau)\right)$ (cf. (6)) using the
aforementioned PGD algorithm. For every budget $b$, we reach the local optimal
point $\tau^{*}$ after a sufficient number of iterations. As $b$ increases,
the optimal travel rates move away from the uncontrolled values while
minimizing $f(\tau)$. So, as depicted in Fig. 1, we obtain lower
$\lambda_{\max}(M(t_{0},\tau^{*}))$ for larger values of $b$. Specifically,
for $b>20$, we have that $f(\tau^{*})$ becomes less than $0$ making the
epidemic network converge faster to a disease-free equilibrium. Note that the
degree of restrictions $b$ depends primarily on the travel rates within the
network. Therefore, the customization for a proper choice in each specific
network is essential.

Figure 1: Optimal values of $f(\tau)$ in problem (6) for different budget
parameters $b$ after changing the travel rates $\tau$ via (7).

In Fig. 2 we observe the cumulative number of infected and recovered cases for
the optimal travel rates obtained from PGD algorithm in Fig. 1 within
Massachusetts counties. As expected, the optimal travel rates associated with
larger $b$ leads to lower number of cumulative cases and a faster epidemic
control due to the stricter lockdown measures imposed.

Similarly, Fig. 3 depicts the number of active infected cases for the optimal
travel rates obtained from the PGD algorithm in Fig. 1. Notably, we achieve a
reduction rate of $-\alpha=-0.0231$, which corresponds to halving infected
cases every $30$ days, for the optimal travel rates obtained for $b>20$.

Figure 2: Number of cumulative cases (infected, quarantined, and recovered)
for different constraints over the travel rates for the state of
Massachusetts. Figure 3: Number of active (asymptomatic and symptomatic
infected) cases for different constraints over the travel rates for the state
of Massachusetts.

Optimizing Quarantine Rates. Fig. 4 and Fig. 5 compare the optimal
quarantining rates calculated by our method with several other quarantining
policies. In uniform quarantining, the quarantine rates of all locations are
chosen to be the same and such that the total economic cost equals the cost
associated with the optimal quarantine policy. In random quarantining however,
the rates are randomly chosen from a uniform distribution and the parameters
are chosen such that we have the same total economic cost as the optimal
quarantine policy. In a uniformly bounded decline policy, quarantine rates are
determined to ensure that the infection decay rates of all nodes remain within
a specified bound, resulting in the same economic cost as our optimal policy.
For the same economic costs, we observe that optimal quarantine rates minimize
the number of infected (cumulative and active) cases in the network.

##  VII Conclusions

In this paper, we presented a framework for epidemic control with two
approaches. The first focuses on strategically reducing travel rates to
contain the virus effectively, utilizing the maximum eigenvalue function and
the PGD algorithm for optimization. The second approach enhances the SIR model
with a quarantine strategy to minimize costs and decrease new infections
rapidly through node-specific rates. We propose a solution that simplifies
optimal quarantining into a weight-balancing problem and establish a link
between optimization constraints and the epidemic’s basic reproduction number.
Finally, applying (Aug-PDGD) to our optimal quarantine problem ensures
exponential stability of the solution.

Figure 4: Number of cumulative cases (infected, quarantined, and recovered)
assuming for different quarantine rates for the state of Massachusetts. For
the optimal policy, $\alpha$ is set to $0.023$, which corresponds to halving
the number of infected cases every $30$ days. Figure 5: Number of active
(asymptomatic and symptomatic infected) cases for different quarantine
policies for the state of Massachusetts. For the optimal policy, $\alpha$ is
set to $0.023$, which corresponds to halving the number of infected cases
every $30$ days.

Our numerical results highlight the efficiency of both approaches in
controlling the epidemic. Although we constructed the network within the
county-level framework, it’s noteworthy that, given data availability, these
methods can be applied to much smaller nodes (e.g., zip code or census track
levels).

##  Appendix A Proof of Theorem 3

Consider the problem formulation in (15). Setting $w={\mathbf{1}}-q$, we
equivalently have

$\begin{split}&{\min}_{w}\quad\sum_{i=1}^{n}\frac{z_{i}^{\mathrm{a}}}{w_{i}^{\mathrm{a}}}+\frac{z_{i}^{s}}{w_{i}^{\mathrm{s}}}\\\
&\text{s.t.}\quad\lambda_{\max}(A_{0})\leq 0\ ,\\\ &\qquad{\mathbf{0}}\leq
w\leq{\mathbf{1}}\ ,\end{split}$ (18)

where from (14)

$A_{0}=\begin{pmatrix}\mathop{\mathrm{diag}}{(w^{\mathrm{a}})}+E^{\prime}&\beta^{\mathrm{s}}\mathop{\mathrm{diag}}{(s(t_{0}))}A\\\
\epsilon
I&\mathop{\mathrm{diag}}{(w^{\mathrm{s}})}-(r^{\mathrm{s}}+1-\alpha)I\end{pmatrix},$

where
$E^{\prime}=\beta^{\mathrm{a}}\mathop{\mathrm{diag}}{(s(t_{0}))}A-(\epsilon+r^{\mathrm{a}}+1-\alpha)I$.

We proceed by finding the range of the decay rate $\alpha$ that ensures the
existence of a non-negative $w$ which satisfies the maximum eigenvalue
constraint in (18). Let us first assume

$\alpha<\min(r^{s}+1,\
\epsilon+r^{a}+1-\max_{i}{(\beta^{\mathrm{a}}\mathop{\mathrm{diag}}{(s(t_{0}))}A)_{ii}}),$
(19)

where $\max_{i}{(\beta^{\mathrm{a}}\mathop{\mathrm{diag}}{(s(t_{0}))}A)_{ii}}$
is the maximum diagonal element of
$\beta^{\mathrm{a}}\mathop{\mathrm{diag}}{(s(t_{0}))}A$, which is known and
fixed. If (19) does not hold, then matrix $A_{0}$ cannot be continuous time
stable (i.e., the first constraint in (18) is not satisfied for any choice of
$w$). More specifically, if (19) is not satisfied for the decay rate $\alpha$,
there exists no $d>{\mathbf{0}}$ satisfying $A_{0}d\leq{\mathbf{0}}$,
regardless of the choice of $w$ (see the first bullet of Lemma (2)).

Let us write

$\displaystyle A_{0}$
$\displaystyle=\mathop{\mathrm{diag}}{(w)}+\begin{pmatrix}E^{\prime}&\beta^{\mathrm{s}}\mathop{\mathrm{diag}}{(s(t_{0}))}A\\\
\epsilon I&-(r^{\mathrm{s}}+1-\alpha)I\end{pmatrix}$
$\displaystyle=\mathop{\mathrm{diag}}{(w)}+B_{0},$ (20)

where $w=(w^{\mathrm{a}},w^{\mathrm{s}})$ $\in\mathds{R}^{2n}$, and
$E^{\prime}=\beta^{\mathrm{a}}\mathop{\mathrm{diag}}{(s(t_{0}))}A-(\epsilon+r^{\mathrm{a}}+1-\alpha)I$.
From (20), we set $A_{0}$ to be the sum of a non-negative diagonal matrix
$\mathop{\mathrm{diag}}{(w)}$ with $B_{0}$. By the Perron-Frobenius theorem,
since adding $\mathop{\mathrm{diag}}{(w)}$ to $B_{0}$ cannot decrease the
maximum eigenvalue of $B_{0}$, in order to have a feasible problem in (18), it
is necessary for $B_{0}$ to be Hurwitz. However, although being necessary,
(19) is not a sufficient condition for $B_{0}$ to be Hurwitz. In order to make
$B_{0}$ Hurwitz, we implement the following approach. We can write

$\displaystyle B_{0}=$
$\displaystyle\begin{pmatrix}\beta^{\mathrm{a}}\mathop{\mathrm{diag}}{(s(t_{0}))}A-(\epsilon+r^{\mathrm{a}}+1)I&\beta^{\mathrm{s}}\mathop{\mathrm{diag}}{(s(t_{0}))}A\\\
\epsilon I&-(r^{\mathrm{s}}+1)I\end{pmatrix}$ $\displaystyle+\alpha
I=C_{0}+\alpha I.$

Again, to have $B_{0}$ Hurwitz, the matrix $C_{0}$ needs to be Hurwitz and
$\alpha$ must be chosen such that adding it to the diagonal elements of
$C_{0}$ does not make $C_{0}$ unstable. So, to have a feasible problem,
$\alpha$ should satisfy $\alpha<-\lambda_{\max}(C_{0})$.

Having $B_{0}$ Hurwitz, implies that $B_{0}$ is invertible and
$-B_{0}^{-1}\geq{\mathbf{0}}$ ([37, Thm. $10.3$]). Using Lemma 2, having
$A_{0}$ continuous-time stable is equivalent to
$-B_{0}^{-1}\mathop{\mathrm{diag}}{(w)}$ being discrete-time stable. Since
changing the order of the product of two matrices doesn’t affect the
eigenvalues of the product, we equivalently want the matrix
$\mathop{\mathrm{diag}}{(w)}(-B_{0}^{-1})$ to be discrete-time stable.
Applying part $(3)$ of Lemma 2, this is equivalent to
$-B_{0}^{-1}-{\mathop{\mathrm{diag}}{(w)}}^{-1}$ being continuous-time stable.

For a non-negative $w$, let us set $v_{i}={w_{i}}^{-1}$, then(18) without
considering $w\leq{\mathbf{1}}$ is equivalent to

$\begin{split}&{\min}_{v}\quad z^{T}\,v\\\ &\text{s.
t.}\quad\lambda_{max}(-B_{0}^{-1}-{\mathop{\mathrm{diag}}{(v)}})\leq 0,\\\
&v\geq{\mathbf{0}}\end{split}$ (21)

where $-B_{0}^{-1}$ is Metzler. Due to the strong connectivity of A, the two
upper blocks in $B_{0}$ are strongly connected. The two lower diagonal blocks
ensure the strong connectivity of $B_{0}$. Hence, $-B_{0}$ is irreducible as
well. The inverse of an irreducible matrix is also irreducible and we have
$-B_{0}^{-1}$ being irreducible. The objective of (21) is to find the optimal
solution $v^{*}$ such that $z^{T}v^{*}$ is minimized and the matrix
$-B_{0}^{-1}-{\mathop{\mathrm{diag}}{(v^{*})}}$ is Hurwitz. Applying, [15,
Theorem $3$], we have that (21) can be reduced to a matrix balancing problem.
More specifically, consider $d^{*}\in\mathcal{R}^{2n}_{>0}$ such that the
matrix

$(\mathop{\mathrm{diag}}{(d^{*})})^{-1}\mathop{\mathrm{diag}}{(z)}(-B_{0}^{-1})\mathop{\mathrm{diag}}{(d^{*})}$

is weight balanced. From [15, Theorem $3$] we have that (21) has a unique
optimal solution

$v^{*}=(\mathop{\mathrm{diag}}{(d^{*})})^{-1}(-B_{0}^{-1})\mathop{\mathrm{diag}}{(d^{*})}\mathbf{1}\
.$ (22)

So, from the above relation we have shown that $\eqref{revised_formulation1}$
is equivalent to a matrix balancing problem for $w\geq{\mathbf{0}}$, or
equivalently $q\leq{\mathbf{1}}$.

To complete our analysis, we also require $w\leq{\mathbf{1}}$ (or equivalently
$q\geq{\mathbf{0}}$). Let us see the solution $v^{*}$ in more details. As
discussed previously, $B_{0}$ is a Hurwitz matrix with non-negative off-
diagonal elements and negative diagonal elements. Suppose $m$ is the absolute
value of the most negative diagonal element of $B_{0}$, or

$m=\max_{i}|(B_{0})_{ii}|>0\ .$

Let us write $B_{0}=B-mI$, where $B$ is a non-negative matrix. From the Perron
Frobenius theorem we have $\rho(B)>0$ and $\rho(B)-m$ is an eigenvalue of the
Hurwitz matrix $B_{0}$.

$\displaystyle B_{0}^{-1}$
$\displaystyle=(B-mI)^{-1}=\frac{1}{m}(\frac{1}{m}B-I)^{-1}$ (23)
$\displaystyle=-\frac{1}{m}\sum_{n=0}^{\infty}\left(\frac{1}{m}B\right)^{n}=-\frac{1}{m}\sum_{n=0}^{\infty}\left(I+\frac{1}{m}B_{0}\right)^{n}.$

By substituting (23) in (22), $v^{*}$ can be derived as

$v^{*}=\frac{1}{m}(\mathop{\mathrm{diag}}{(d^{*})})^{-1}\sum_{n=0}^{\infty}\left(I+\frac{1}{m}B_{0}\right)^{n}\mathop{\mathrm{diag}}{(d^{*})}\mathbf{1}\
.$ (24)

Thus, to have $w\leq{\mathbf{1}}$, we require $v^{*}\geq{\mathbf{1}}$ or
equivalently, the sum of row elements of

$(\mathop{\mathrm{diag}}{(d^{*})})^{-1}\sum_{n=0}^{\infty}\left(I+\frac{1}{m}B_{0}\right)^{n}\mathop{\mathrm{diag}}{(d^{*})}$

should be greater or equal to $m$. Indeed, defining
$\sum_{n=0}^{\infty}\left(I+\frac{1}{m}B_{0}\right)^{n}=Q$, we need

$\sum_{j=1}^{2n}\frac{d^{*}_{j}}{d^{*}_{i}}Q_{ij}\geq m\qquad\forall
i\in\\{1,2,\dots,2n\\}\ .$ (25)

Notice that $\frac{1}{m}B_{0}$ is a matrix with non-negative off-diagonal
elements and negative diagonal elements of $-1\leq\frac{(B_{0})_{ii}}{m}<0$.
Thus, $(I+\frac{1}{m}B_{0})$ is a non-negative matrix with diagonal elements,
and

$0\leq(I+\frac{1}{m}B_{0})_{ii}<1\ .$

Hence,

$Q\geq I+[I+\frac{1}{m}B_{0}]+[(I+\frac{1}{m}B_{0})^{2}]=\bar{Q}$ (26)

Letting $(I+\frac{1}{m}B_{0})_{ii}=0$, we have

$\displaystyle(I+\frac{1}{m}B_{0})\geq\begin{pmatrix}\frac{1}{m}\tilde{A}&\frac{\beta^{\mathrm{s}}}{m}\mathop{\mathrm{diag}}{(s(t_{0}))}A\\\
\frac{\epsilon}{m}I&{\mathbf{0}}\end{pmatrix}\ ,$ (27)

where $\tilde{A}=\beta^{\mathrm{a}}\mathop{\mathrm{diag}}{(s(t_{0}))}A$ that
its diagonal elements are replaced with zero. Substituting (27) into (26), we
derive a lower bound of $\bar{Q}$ as $\tilde{Q}$:

$\displaystyle\tilde{Q}=$
$\displaystyle\begin{pmatrix}\tilde{Q}_{11}&\tilde{Q}_{12}\\\
\frac{\epsilon}{m}I+\frac{\epsilon}{m^{2}}\tilde{A}&I+\frac{\epsilon\beta^{\mathrm{s}}}{m^{2}}\mathop{\mathrm{diag}}{(s(t_{0}))}A\end{pmatrix},$
(28)

where

$\tilde{Q}_{11}=I+\frac{\tilde{A}}{m}+\frac{1}{m^{2}}(\tilde{A}^{2}+\epsilon\beta^{\mathrm{s}}\mathop{\mathrm{diag}}{(s(t_{0}))}A),$
$\tilde{Q}_{12}=\frac{\beta^{\mathrm{s}}}{m}\mathop{\mathrm{diag}}{(s(t_{0}))}A+\frac{\beta^{\mathrm{s}}}{m^{2}}\tilde{A}\mathop{\mathrm{diag}}{(s(t_{0}))}A\
.$

From the definition of $B_{0}$ in (20), $m$ can be written as

$\displaystyle m=\max\Big{(}$ $\displaystyle r^{s}+1-\alpha,$
$\displaystyle\epsilon+r^{a}+1-\alpha-\min\left(\beta^{\mathrm{a}}\operatorname{diag}(s(t_{0}))A\right)_{ii}\Big{)}\
.$

Substituting Assumption 1 into the definition of $\tilde{Q}$, we have

$\displaystyle\tilde{Q}_{ii}$
$\displaystyle=1+\frac{1}{m^{2}}((\tilde{A}^{2})_{jj}+\epsilon\beta^{\mathrm{s}}s_{i}(t_{0})A_{jj})$
$\displaystyle\geq
1+\frac{\epsilon}{m^{2}}\beta^{\mathrm{s}}s_{i}(t_{0})A_{jj}\geq
1+\frac{x}{m^{2}}$ $\displaystyle\geq m\ ,\quad\,\quad i\in\\{1,\dots,n\\}\ ,$

and

$\displaystyle\tilde{Q}_{ii}$
$\displaystyle=1+\frac{1}{m^{2}}\epsilon\beta^{\mathrm{s}}s_{j}(t_{0})A_{jj}\geq
1+\frac{x}{m^{2}}\geq m$ $\displaystyle\geq m\ ,\quad\,\quad
i\in\\{n+1,\dots,2n\\}\ .$

Finally,

$\sum_{j=1}^{2n}\frac{d^{*}_{j}}{d^{*}_{i}}Q_{ij}\geq
Q_{ii}\geq\bar{Q}_{ii}\geq\tilde{Q}_{ii}\geq m\ ,$ (29)

for all $i\in\\{1,2,\dots,2n\\}$. Hence, $v^{*}\geq{\mathbf{1}}$, and
equivalently $q^{*}\geq{\mathbf{0}}$. Subsequently, substituting the
definitions of $v$ and $w$, the optimal quarantined rates are calculated as

$q^{*}=\mathbf{1}+1./\big{(}{(\mathop{\mathrm{diag}}{(d^{*})})}^{-1}B_{0}^{-1}(\mathop{\mathrm{diag}}{(d^{*})})\mathbf{1}\big{)}\
.$

##  Appendix B Proof of Theorem 4

We will establish that initiating from $R_{0}<1$ leads us to
$\lambda_{\max}(M(t_{0}))\leq 0$. Assume

$R_{0}\leq 1\ .$

We have

$M(t_{0})=\begin{pmatrix}E&\beta^{\mathrm{s}}\mathop{\mathrm{diag}}{(s(t_{0}))}A\\\
\epsilon I&-(r^{\mathrm{s}}+\mathop{\mathrm{diag}}{(q^{\mathrm{s}})})I\\\
\end{pmatrix},$

where
$E=\beta^{\mathrm{a}}\mathop{\mathrm{diag}}{(s(t_{0}))}A-(\epsilon+r^{\mathrm{a}}+\mathop{\mathrm{diag}}{(q^{\mathrm{a}}}))I$.
Let us write $M(t_{0})$ as

$M(t_{0})=F+V,$

where

$F=\begin{pmatrix}\beta^{\mathrm{a}}\mathop{\mathrm{diag}}{(s(t_{0}))}A&\beta^{\mathrm{s}}\mathop{\mathrm{diag}}{(s(t_{0}))}A\\\
{\mathbf{0}}&{\mathbf{0}}\\\ \end{pmatrix},$

and

$V=\begin{pmatrix}-(\epsilon+r^{\mathrm{a}})I-\mathop{\mathrm{diag}}{(q^{\mathrm{a}})})&{\mathbf{0}}\\\
\epsilon I&-(r^{\mathrm{s}}+\mathop{\mathrm{diag}}{(q^{\mathrm{s}})})I\\\
\end{pmatrix}.$

Note that $F$ is nonnegative and $V$ is Metzler and Hurwitz. Using the
definition of $R_{0}$ in [11, Definition $1$] we have

$R_{0}=\rho(FV^{-1})=\rho(-FV^{-1})\leq 1.$ (30)

Utilizing the early results of Thm. 3, (30) is equivalent with

$\lambda_{\max}(F+V)\leq 0,$

(i.e., $F+V$ is Hurwitz). Hence, we have

$\lambda_{\max}(M(t_{0}))=\lambda_{\max}(F+V)\leq 0.$

Therefore, the set $\lambda_{\max}(M(t_{0}))\leq 0$ is equivalent with the set
$R_{0}\leq 1$.

## Appendix C

Given the left and right eigenvectors of $\lambda_{\max}(M(t_{0},q))$ as
$v_{\max}^{T}$ and $u_{\max}$, respectively, we argue that

$\left\|\frac{v_{\max}\circ u_{\max}}{v_{\max}^{T}u_{\max}}\right\|$

is bounded for any choice of $M(t_{0},q)$. Without loss of generality, we may
assume the left and right eigenvectors are normalized to have unit norm. In
that case, the norm of the numerator is automatically bounded.

Applying Lemma 1 to $M(t_{0},q)$ and $M(t_{0},q)^{T}$, we have that $u_{\max}$
and $v_{\max}$ are positive as are the respective dominant eigenvectors.
Hence, the dot product of the corresponding Perron-Frobenius left and right
eigenvectors is always positive. This also implies that the denominator is
bounded away from zero, for otherwise by compactness we could find a matrix
$M$ such that $v_{\max}^{T}u_{\max}=0$, contradicting positivity of these
vectors. Therefore, $\|\nabla g_{1}(q)\|$ is bounded.

##  Appendix D Proof of Theorem 6

It is obvious that the gradient of $f(q)$ is Lipschitz as long as every
$q_{i}(0)$ is bounded away from $1$. We thus need to argue that the primal-
dual dynamics keep all $q_{i}(0)$ uniformly bounded away from one.

We start by expressing the primal dynamics as given by (17a):

$\dot{q}(t)=-\nabla f(q(t))-\sum_{i=1}^{4n+1}[\rho
g_{i}(q(t))+\lambda_{i}(t)]_{+}\nabla g_{i}(q(t)).$ (31)

Here, $\lambda(t)\in\mathds{R}_{\geq 0}^{(4n+1)}$ is a vector that remains in
a compact region of space for all $t\geq 0$ by Proposition 2. Therefore, using
(5), we can state that the second part of (31) is bounded.

Note that as $q_{i}\rightarrow 1$, we have $|\nabla_{i}f(q(0))|$ approaches
$\infty$, and $\dot{q}_{i}<0$. This immediately implies there is some
$q^{\prime}<1$ such that the region $[0,q^{\prime}]^{n}$ remains invariant
under the differential equation (31). Remaining within this region implies the
Lipschitz continuity of $\nabla f(q)$.

##  Appendix E Proof of Proposition 3

Suppose $q^{*}$ is a local minimum of (15). Therefore, $q^{*}_{i}\neq 1$ for
all $i\in\\{1,2,\dots,2n\\}$. This implies the active set is

$\mathcal{I}=\\{j:g_{j}(q^{*})=0,\ j\in J\\}\ ,$

where $J=\\{1,2,\dots,2n+1\\}$.

First, assume $\mathcal{I}\neq J$. From the definition provided for $g_{1}$,
$\nabla g_{1}=\frac{v_{\max}\circ u_{\max}}{v_{\max}^{T}u_{\max}}$ is an
element-wise positive vector by the same analysis provided in Appendix D. It
is therefore evident that any $2n$ vectors of $g_{j}$, for all $j\in J$, would
have linearly independent gradients at any point $q^{*}$.

Now, consider $\mathcal{I}=J$. This implies $q^{*}=\mathbf{0}$. However, for
our model of the outbreak in (15), setting all quarantine rates equal to zero
implies having no interventions. Hence, the constraint function
$g_{1}=\lambda_{\max}(M(t_{0},q=\mathbf{0}))+\alpha$ would not be active.
Hence, $\mathcal{I}$ cannot include all $j\in J$, and this proves the non-
linearity of the active set at any local minimum $q^{*}$.

## References

  * [1] J. Center, “Covid-19 dashboard,” _Johns Hopkins University_ , 2021.
  * [2] A. Linde, M. Rotzén-Östlund, B. Zweygberg-Wirgart, S. Rubinova, and M. Brytting, “Does viral interference affect spread of influenza?” _Eurosurveillance_ , vol. 14, no. 40, p. 19354, 2009.
  * [3] C. Nowzari, V. M. Preciado, and G. J. Pappas, “Analysis and control of epidemics: A survey of spreading processes on complex networks,” _IEEE Control Systems Magazine_ , vol. 36, no. 1, pp. 26–46, 2016.
  * [4] P. E. Paré, C. L. Beck, and T. Başar, “Modeling, estimation, and analysis of epidemics over networks: An overview,” _Annual Reviews in Control_ , vol. 50, pp. 345–360, 2020.
  * [5] S. Lee, G. Chowell, and C. Castillo-Chávez, “Optimal control for pandemic influenza: The role of limited antiviral treatment and isolation,” _Journal of Theoretical Biology_ , vol. 265, no. 2, pp. 136–150, 2010.
  * [6] M. Hayhoe, F. Barreras, and V. M. Preciado, “Multitask learning and nonlinear optimal control of the covid-19 outbreak: A geometric programming approach,” _Annual Reviews in Control_ , vol. 52, pp. 495–507, 2021.
  * [7] A. Khanafer and T. Basar, “An optimal control problem over infected networks,” in _Proceedings of the International Conference of Control, Dynamic Systems, and Robotics_ , vol. 125, 2014, pp. 1–6.
  * [8] F. Liu and M. Buss, “Optimal control for heterogeneous node-based information epidemics over social networks,” _IEEE Transactions on Control of Network Systems_ , vol. 7, no. 3, pp. 1115–1126, 2020.
  * [9] A. R. Hota, J. Godbole, and P. E. Paré, “A closed-loop framework for inference, prediction, and control of SIR epidemics on networks,” _IEEE Transactions on Network Science and Engineering_ , vol. 8, no. 3, pp. 2262–2278, 2021.
  * [10] V. S. Mai, A. Battou, and K. Mills, “Distributed algorithm for suppressing epidemic spread in networks,” _IEEE Control Systems Letters_ , vol. 2, no. 3, pp. 555–560, 2018.
  * [11] K. D. Smith and F. Bullo, “Convex optimization of the basic reproduction number,” _IEEE Transactions on Automatic Control_ , vol. 68, no. 7, pp. 4398–4404, 2023.
  * [12] V. M. Preciado, M. Zargham, C. Enyioha, A. Jadbabaie, and G. J. Pappas, “Optimal resource allocation for network protection against spreading processes,” _IEEE Transactions on Control of Network Systems_ , vol. 1, no. 1, pp. 99–108, 2014.
  * [13] P. Van Mieghem, D. Stevanović, F. Kuipers, C. Li, R. Van De Bovenkamp, D. Liu, and H. Wang, “Decreasing the spectral radius of a graph by link removals,” _Physical Review E_ , vol. 84, no. 1, p. 016101, 2011.
  * [14] Q. Ma, Y.-Y. Liu, and A. Olshevsky, “Optimal lockdown for pandemic control,” _arXiv preprint arXiv:2010.12923_ , 2021.
  * [15] R. Ofir, F. Bullo, and M. Margaliot, “Minimum effort decentralized control design for contracting network systems,” _IEEE Control Systems Letters_ , vol. 6, pp. 2731–2736, 2022.
  * [16] G. Giordano, F. Blanchini, R. Bruno, P. Colaneri, A. D. Filippo, A. D. Matteo, and M. Colaneri, “Modelling the COVID-19 epidemic and implementation of population-wide interventions in Italy,” _Nature Medicine_ , vol. 26, no. 6, pp. 855–860, 2020.
  * [17] J. R. Birge, O. Candogan, and Y. Feng, “Controlling epidemic spread: Reducing economic losses with targeted closures,” _Management Science_ , vol. 68, no. 5, pp. 3175–3195, 2022.
  * [18] S. M. Kissler, J. R. Fauver, C. Mack, C. Tai, K. Y. Shiue, C. C. Kalinich, S. Jednak, I. M. Ott, C. B. Vogels, J. Wohlgemuth, J. Weisberger, J. DiFiori, D. J. Anderson, J. Mancell, D. D. Ho, N. D. Grubaugh, and Y. H. Grad, “Viral dynamics of sars-cov-2 infection and the predictive value of repeat testing,” _medRxiv_ , 2020.
  * [19] C. Viboud, O. N. Bjørnstad, D. L. Smith, L. Simonsen, M. A. Miller, and B. T. Grenfell, “Synchrony, waves, and spatial hierarchies in the spread of influenza,” _Science_ , vol. 5772, no. 312, p. 447–451, 2006.
  * [20] Y. Kang, S. Gao, Y. Liang, M. Li, J. Rao, and J. Kruse, “Multiscale dynamic human mobility flow dataset in the US during the COVID-19 epidemic,” _Scientific data_ , vol. 7, no. 1, p. 390, 2020.
  * [21] R. Sinkhorn and P. Knopp, “Concerning nonnegative matrices and doubly stochastic matrices,” _Pacific Journal of Mathematics_ , vol. 21, no. 2, pp. 343–348, 1967.
  * [22] A. I. Rikos, T. Charalambous, and C. N. Hadjicostis, “Distributed weight balancing over digraphs,” _IEEE Transactions on Control of Network Systems_ , vol. 1, no. 2, pp. 190–201, June 2014.
  * [23] M. Idel, “A review of matrix scaling and sinkhorn’s normal form for matrices and positive maps,” _arXiv preprint arXiv:1609.06349_ , 2016.
  * [24] M. B. Cohen, A. Madry, D. Tsipras, and A. Vladu, “Matrix scaling and balancing via box constrained newton’s method and interior point methods,” in _$58^{th}$ Annual Symposium on Foundations of Computer Science_, 2017, pp. 902–913.
  * [25] D. G. Luenberger, _Linear and nonlinear programming_. Klewer Academic Publishers, 2021.
  * [26] X. Chen, H. Qi, L. Qi, and K.-L. Teo, “Smooth convex approximation to the maximum eigenvalue function,” _Journal of Global Optimization_ , vol. 30, no. 2, pp. 253–270, 2004.
  * [27] J. R. Magnus, “On differentiating eigenvalues and eigenvectors,” _Econometric Theory_ , vol. 1, no. 2, p. 179–191, 1985.
  * [28] J. Nocedal and S. J. Wright, “Quadratic programming,” _Numerical optimization_ , vol. 35, pp. 448–492, 2006.
  * [29] M. Zedek, “Continuity and location of zeros of linear combinations of polynomials,” _Proceedings of the American Mathematical Society_ , vol. 16, no. 1, pp. 78–84, 1965.
  * [30] Y. Nesterov, _Introductory lectures on convex optimization: A basic course_. Springer Science & Business Media, 2003, vol. 87.
  * [31] J. E. Cohen, “Convexity of the dominant eigenvalue of an essentially nonnegative matrix,” _Proceedings of the American Mathematical Society_ , vol. 81, no. 4, pp. 657–658, 1981.
  * [32] Y. Tang, G. Qu, and N. Li, “Semi-global exponential stability of augmented primal–dual gradient dynamics for constrained convex optimization,” _Systems & Control Letters_, vol. 144, p. 104754, 2020.
  * [33] Massachusetts Legislature, “2020 Census Data by County,” https://malegislature.gov/Redistricting/MassachusettsCensusData/County.
  * [34] U.S. Bureau of Economic Analysis, “Gross Domestic Product by County, Metro, and Other Areas,” https://www.bea.gov/data/gdp/gdp-county-metro-and-other-areas.
  * [35] e. a. Mitch Smith and, “Us counties covid 19 dataset,” 2024. [Online]. Available: https://www.kaggle.com/dsv/8581321
  * [36] A. Hortaçsu, J. Liu, and T. Schwieg, “Estimating the fraction of unreported infections in epidemics with a known epicenter: An application to covid-19,” _Journal of Econometrics_ , vol. 220, no. 1, pp. 106–129, 2021.
  * [37] F. Bullo, _Lectures on Network Systems_ , 1.6 ed. Kindle Direct Publishing, 2022. [Online]. Available: https://fbullo.github.io/lns

|  Mahtab Talaei is currently pursuing her Ph.D. in the Division of Systems
Engineering at Boston University in Boston, MA. She obtained her B.Sc. and
M.Sc. degrees in Electrical Engineering from Isfahan University of Technology
(IUT) in Isfahan, Iran, in 2019 and 2022, respectively. Her research primarily
involves utilizing optimization and machine learning techniques to create
predictive models for healthcare applications.
---|---
|  Apostolos I. Rikos (M’16) is an Assistant Professor at the Artificial
Intelligence Thrust of the Information Hub, The Hong Kong University of
Science and Technology (Guangzhou), Guangzhou, China. He is also affiliated
with the Department of Computer Science and Engineering, The Hong Kong
University of Science and Technology, Clear Water Bay, Hong Kong, China. He
received his B.Sc., M.Sc., and Ph.D. degrees in Electrical Engineering from
the Department of Electrical and Computer Engineering, University of Cyprus in
2010, 2012, and 2018, respectively. In 2018, he joined the KIOS Research and
Innovation Center of Excellence in Cyprus, where he was a Research Lecturer.
He joined the Division of Decision and Control Systems at KTH Royal Institute
of Technology as a Postdoctoral Researcher in 2020 and the Department of
Electrical and Computer Engineering, Division of Systems Engineering, at
Boston University as a Postdoctoral Associate in 2023. His research interests
are in the areas of distributed optimization and learning, distributed network
control and coordination, privacy and security, and algorithmic design.
---|---
|  Alex Olshevsky (Member, IEEE) received the B.S. degrees in Applied
Mathematics and Electrical Engineering from the Georgia Institute of
Technology, Atlanta, GA, USA, both in 2004, and the M.S. and Ph.D. degrees in
electrical engineering and computer science from the Massachusetts Institute
of Technology, Cambridge, MA, USA, in 2006 and 2010, respectively. He is
currently an Associate Professor with the Department of Electrical and
Computer Engineering, Boston University, Boston, MA, USA. His research
interests include control systems, optimization, and net- work science. He was
the recipient of the National Science Foundation CAREER Award, the Air Force
Young Investigator Award, the ICS Prize from INFORMS for best paper on the
interface of operations research and computer science, and the SIAM Paper
Prize for annual paper from the SIAM Journal on Control and Optimization
chosen to be reprinted in SIAM Review.
---|---
|  Laura F. White is a Professor of Biostatistics and Associate Director of
the Population Health Data Science Program at the Boston University School of
Public Health. She received her BS in Mathematics and Statistics from Utah
State University and her PhD in Biostatistics from Harvard University. She co-
directs the graduate program in Biostatistics at BU and is an co-director of
the Data Science and Surveillance Core of the Center on Emerging Infectious
Diseases. Her research interests are in developing novel statistical
approaches to understanding infectious disease transmission dynamics and
disease burden. Her work is supported by the NIH and CDC.
---|---
|  Ioannis Ch. Paschalidis (Fellow, IEEE) received a diploma in Electrical and
Computer Engineering (ECE) from the National Technical University of Athens,
Greece, in 1991, and the MS and PhD degrees in Electrical Engineering and
Computer Science (EECS), from the Massachusetts Institute of Technology (MIT),
Cambridge, Massachusetts, in 1993 and 1996, respectively. In September 1996,
he joined Boston University where he has been ever since. He is currently a
Distinguished Professor of Engineering and the director of the Hariri
Institute for Computing. His current research interests lie in the fields of
optimization, control, stochastic systems, robust learning, and computational
medicine/biology. His has been recognized with am NSF CAREER award, several
best paper awards, and an IBM/IEEE Smarter Planet Challenge Award. He was an
invited participant at the 2002 Frontiers of Engineering Symposium organized
by the National Academy of Engineering, at the 2014 National Academies Keck
Futures Initiative Conference, and at a 2024 National Academies Symposium on
Alzheimer’s disease. He is a Fellow of IEEE, IFAC, and the Asia-Pacific AI
Association. From 2013 to 2019 he was the founding Editor-in-Chief of the IEEE
Transactions on Control of Network Systems and he is the General Co-Chair of
the 2025 IEEE Conference on Decision and Control.
---|---