# Model Predictive Control of Spreading Processes
via Sparse Resource Allocation

Ruigang Wang, Armaghan Zafar and Ian R. Manchester This work was supported by
the Air Force Office of Scientific Research Grant FA2386-19-1-4076 and the NSW
Defence Innovation Network.The authors are with the Australian Centre for
Field Robotics, and Sydney Institute for Robotics and Intelligent Systems, The
University of Sydney, Sydney, NSW 2006, Australia (e-mail: ruigang.wang,
armaghan.zafar. ian.manchester@sydney.edu.au).

###### Abstract

In this paper, we propose a model predictive control (MPC) method for real-
time intervention of spreading processes, such as epidemics and wildfire, over
large-scale networks. The goal is to allocate budgeted resources each time
step to minimize the risk of an undetected outbreak, i.e. the product of the
probability of an outbreak and the impact of that outbreak. By using dynamic
programming relaxation, the MPC controller is reformulated as a convex
optimization problem, in particular an exponential cone programming. We also
provide sufficient conditions for the closed-loop risks to asymptotically
decrease and a method to estimate the upper bound of when the risk will
monotonically decrease. Numerical results are provided for a wildfire example.

##  I Introduction

Modeling, analysis and control of epidemics and wildfires have been studied
intensely due to their significant social and economic impacts, see the recent
surveys [1, 2, 3, 4]. These events are usually modeled as spreading process in
which an initial localized outbreak spreads rapidly to neighboring nodes in a
network, and causes harm or risk of harm to the nodes it touches. The task of
intervention is to reduce the risks of an outbreak by allocating resources
(e.g. vaccines, social distance) to nodes/links.

There are several challenges for resource allocation of large-scale spreading
processes. Firstly, the available resources at any given time are often
limited. Secondly, some types of resources (e.g. water bombing attacks on a
wildfire) can be difficult or impossible to widely distribute across the
network. Moreover, the large-scale nature of the spreading processes, e.g.
global transport network for epidemics or large geographic areas for
bushfires, imposes computational challenges in terms of real-time computation
of responses. Therefore, it is important to develop scalable resource
allocation methods resulting in sparse solutions, subject to meaningful budget
constraints.

A natural tool for developing such intervention strategies is optimal control
with the objective of minimizing the risk of an outbreak, which is defined as
the expected accumulated impact of the event over the time horizon. Typical
models for spreading processes are the Susceptible-Infected-Susceptible (SIS)
model and the Susceptible-Infected-Removed (SIR) model [5]. These models are
stochastic which can further increase the complexity of optimal control
problem. By using mean-field approximation [6], one can obtain approximate
models in the form of nonlinear deterministic ordinary differential equations
(ODEs). It is proven [7] that solutions of these deterministic ODE systems are
upper bounds of the expected value of solutions of the corresponding
stochastic models. Thus, they are usually the object of study in control
design.

Various methods have been proposed for resource allocation subject to budget
constraints. Convex optimization techniques such as semidefinite programming
(SDP) [8] and geometric programming (GP) [9, 10] have been applied. In [11], a
novel logarithmic resource model was proposed for encouraging sparsity
solutions for the GP-based approaches. Further explorations in surveillance
study of spreading processes can be found in [12]. However, all these methods
offer static intervention strategies only. That is, the resource allocation is
determined by a single-step optimization problem based on the current state
instead of allocated over time.

Dynamic resource allocation, in particular using Model Predictive Control
(MPC), has been studied recently for epidemic intervention. In [13], nonlinear
MPC is applied to reduce the epidemic spread. The works [14, 15] investigate
MPC approaches to control study of the COVID-19 outbreak in Germany. In [16],
robust economic MPC with rigorous convergence guarantees was developed for the
containment of stochastic epidemic processes. However, these approaches do not
provide sparse resource allocation.

In [17], a multi-stage sparse resource allocation method was developed by
extending the single-stage approach [11, 12]. Different from the MPC
approaches [13, 14, 16, 15], this method allocates resources by minimizing a
risk upper bound, which is a linear value (cost-to-go) function constructed by
dynamic programming. The resulting optimal control problem (OCP) can be
reformulated as a convex optimization problem, in particular an exponential
cone program (ECP), which can be efficiently solved by the state-of-art tools
(e.g. MOSEK [18]).

In this paper we, we build on [17] by providing a theoretical analysis of the
closed-loop behavior when the this multi-stage resource allocation approach is
applied recursively, in an MPC framework. In particular, we prove recursive
feasibility and provide sufficient conditions for the closed-loop risk bound
to asymptotically decrease. Moreover, we formulate a convex optimization
problem to estimate the time at which the risk will begin to monotonically
decrease. We also provide results of the numerical experiments performed by
implementing the proposed MPC approach on a wildfire example.

The remainder of this paper is structured as follows. The spreading process
model and control setup are presented in Section II. Section III gives the
proposed approach and some theoretical results with proofs found in Section
IV. Finally, an illustrative example of a wildfire is presented in Section V.

Notation. Throughout the paper, $\geq$ is an element-wise operator, e.g., the
matrices $A,B$ satisfying $A\geq B$ implies $A_{ij}\geq B_{ij}$ for all
$(i,j)$. And $A\gneqq B$ means $A_{ij}\geq B_{ij}$ for all $(i,j)$ and there
exist some pairs $(i,j)$ such that $A_{ij}>B_{ij}$. We call $A$ a nonnegative
matrix if $A\geq 0$.

##  II Problem Formulation

###  II-A Spreading Process over Networks

We consider an SIS spreading process on a graph
$\mathcal{G}(\mathcal{V},\mathcal{E})$ with node set
$\mathcal{V}=\\{1,2,\ldots,n\\}$ and edge set $\mathcal{E}=\\{(i,j):i\neq
j\in\mathcal{V}\\}$, where each node $i$ has a state $X_{i}(t)\in\\{0,1\\}$
indicating that it is infectious ($X_{i}(t)=1$) or susceptible ($X_{i}(t)=0$).
In this process, an infected node $i$ recovers with probability
$\delta_{i}(t)\Delta t$ to $X_{i}(t+\Delta t)=0$ and an susceptible node $i$
can be infected (i.e. $X_{i}(t+\Delta t)=1$) by its neighbor node $j$ with
probability $\beta_{ij}(t)\Delta t$, where $\Delta t$ is a small time
interval. By using mean-field approximation [6], we obtain $n$ coupled
deterministic differential equations:

$\dot{x}_{i}(t)=(1-x_{i}(t))\sum_{j=1}^{n}\beta_{ij}x_{j}(t)-\delta_{i}x_{i}(t),$
(1)

where $x_{i}(t):=P(X_{i}(t)=1)$ is the probability of a node $i$ being
infected at time $t$.

By applying Euler forward method to (1), we obtain a discrete-time model in
compact form of

$x(k+1)=\bigl{[}A-B(x(k))\bigr{]}x(k),$ (2)

where $x(k)=\bigl{[}x_{1}(kh),\ldots,x_{n}(kh)\bigr{]}^{\top}$ is the state at
time step $k\in\mathbb{N}$, and the matrices $A,B$ are given by

$\begin{split}A_{ij}&=\begin{cases}1-h\delta_{i},\;i=j,\\\
h\beta_{ij},\;(i,j)\in\mathcal{E},\end{cases}B_{ij}(x)=hx_{i}\beta_{ij}.\end{split}$
(3)

We assume that the time interval $h>0$ satisfies $h\delta_{j}\leq 1$ and
$h\sum_{i=1}^{n}\beta_{ij}<1$ for all $j\in\mathcal{V}$. Note that (2) is
well-defined, i.e., $x(0)\in[0,1]^{n}\Rightarrow x(k)\in[0,1]^{n}$ for all
$k$.

###  II-B Risk Model

We define the risk of an outbreak $x$ as the impact over the time horizon,
i.e.,

$R(x)=\sum_{\ell=0}^{\infty}\alpha^{k}Cx(k),$ (4)

where $x(k)$ is state of (2) with initial condition $x(0)=x$. Here
$C=\bigl{[}c_{1},\ldots,c_{n}\bigr{]}>0$ is a row vector defining the cost
associate with each node. The discount factor $\alpha\in(0,1]$ is used to
trade-off between near-term cost and long-term cost.

###  II-C Intervention Model

The goal of intervention is to reduce the risk by allocating certain amount of
resources at each time step $k$ to reduce the spreading rate $\beta_{ij}$ on
particular edges or increase recover rate $\delta_{i}$ on particular nodes. We
assume that changes on $\beta_{ij}$ and $\delta_{i}$ are bounded and persist
for all time, i.e.,

$\begin{split}&\beta_{ij}(-1)\geq\cdots\geq\beta_{ij}(k-1)\geq\beta_{ij}(k)\geq\cdots\geq\underline{\beta}_{ij}>0,\\\
&\overline{\Delta}_{i}>\overline{\delta}_{i}\geq\cdots\geq\delta_{i}(k)\geq\delta_{i}(k-1)\geq\cdots\geq\delta_{i}(-1)>0,\end{split}$

for all $k\in\mathbb{N}$, where $\beta_{ij}(-1)=\overline{\beta}_{ij}$ and
$\delta_{i}(-1)=\underline{\delta}_{i}$ are the unmodified rates. It is easy
to verify that

$\overline{A}\geq A(k-1)\geq A(k)\geq\underline{A},\quad\forall
k\in\mathbb{N},$ (5)

where the matrices $\overline{A},\underline{A}$ are given by (3) with rates
$(\overline{\beta}_{ij},\underline{\delta}_{i})$ and
$(\underline{\beta}_{ij},\overline{\delta}_{i})$, respectively. To make
$\underline{A}$ nonnegative we assume $h\overline{\Delta}_{i}\leq 1$ for all
$i\in\mathcal{V}$.

###  II-D Resource Model

For any two matrices $A^{1}\geq A^{2}$, we defined the resource required for
the transition from $A^{1}$ to $A^{2}$ as follows

$\begin{split}\Gamma(A^{1},A^{2})=\sum_{ij}W_{ij}U_{ij}(A^{1},A^{2}),\end{split}$
(6)

where $W_{ij}$ are the weightings on resource costed by reducing $\beta_{ij}$
and increasing $\delta_{i}$, respectively. And $U$ is defined as

$U_{ij}(A^{1},A^{2})=\begin{cases}\log\bigl{(}\beta_{ij}^{1}/\beta_{ij}^{2}\bigr{)},\;(i,j)\in\mathcal{E},\\\
\log\bigl{(}(\overline{\Delta}_{i}-\delta_{i}^{1})/(\overline{\Delta}_{i}-\delta_{i}^{2})\big{)},\;i=j.\end{cases}$
(7)

Given the matrix $A^{1}$ and resource allocation allocation $U$, one can
compute the matrix $A^{2}$ by

$A^{2}=f(A^{1},U),$ (8)

where $A^{2}$ can be computed via (3) with $\beta_{ij}^{2},\delta_{i}^{2}$
given by

$\begin{split}\beta_{ij}^{2}=\beta_{ij}^{1}e^{-U_{ij}},\quad\delta_{i}^{2}=\bigl{(}1-e^{-U_{ii}}\bigr{)}\overline{\Delta}+e^{-U_{ii}}\delta_{i}^{1}.\end{split}$

The above logarithmic resource model was proposed in [11] to encourage sparse
resource allocation. One implication of the above model is that certain
proportional decrease for $\beta_{ij}$ always has the same cost. Note that it
is impossible for $\beta_{ij}$ to be reduced to 0 since that would take
infinite resources.

###  II-E Problem Statement

In this paper we are interested in resource-constrained intervention of the
spreading process (2) using model predictive control. Specifically, the
deliverable resource at each time step is bounded by $\overline{\Gamma}>0$.
The goal of intervention is to minimize the risk $R(x(k))$ given the outbreak
detection $x(k)$ at time step $k$, which can be formulated as an optimal
control problem (OCP):

$\begin{split}\min_{\beta,\delta}\quad&R(x(k))\\\
\mathrm{s.t.}\quad&x(\ell+1|k)=[A(\ell|k)-B(x,\ell|k)]x(\ell|k),\\\
&x(0|k)=x(k),\\\
&\beta_{ij}({\ell-1|k})\geq\beta_{ij}({\ell|k})\geq\underline{\beta}_{ij},\\\
&\delta_{i}(\ell-1|k)\leq\delta_{i}(\ell|k)\leq\overline{\delta}_{i},\\\
&\Gamma(\ell|k)\leq\overline{\Gamma},\;\ell=0,\ldots,L-1\end{split}$ (9)

where $L$ is the time horizon and $\Gamma(\ell|k)$ is short for
$\Gamma(A(\ell-1|k),A(\ell|k))$. Here $\beta_{ij}(-1|k)=\beta_{ij}(k-1)$ and
$\delta_{i}(-1|k)=\delta_{i}(k-1)$ refer to the rates at the previous time
step. For evaluating $R(x(k))$ we use $\beta_{ij}(\ell|k)=\beta_{ij}(L-1|k)$
and $\delta_{i}(\ell|k)=\delta_{i}(L-1|k)$ for $\ell\geq L$.

Note that OCP (9) is a non-convex optimization problem due to the complex
dynamics (2), which may be intractable for real-time intervention of large-
scale spreading processes.

##  III Dynamic Programming Relaxation

In this section we present a convex relaxation for OCP (9). The key idea is to
minimize an particular upper bound of the risk (4), which can be constructed
via a dynamic programming formulation. We will show that under some mild
assumptions the risk bound will decrease asymptotically after certain steps.

###  III-A Risk Bounds

It is well-known that positive linear systems with non-negative linear costs
admit linear value (cost-to-go) functions [19, 20]. Similarly, we will use a
dynamic-programming-alike approach to construct the risk bounds of nonlinear
positive system (2).

######  Proposition 1.

Suppose that the rate sequences $\beta_{ij}$ and $\delta_{i}$ satisfy (5). If
there exists a sequence of nonnegative row vectors $0\leq
p(\ell|k)\in\mathbb{R}^{1\times n}$ such that

$\begin{split}p(\ell|k)&\geq C+\alpha p(\ell+1|k)A(\ell|k),\;\forall
0\leq\ell<L,\\\ p(L|k)&=p(L-1|k),\end{split}$ (10)

then we have $R(x(k))\leq\bar{R}(x(k)):=p(0|k)x(k)$.

The above proposition appears in [17] and its continuous-time version can be
found in [21]. The underlying idea is that
$\bar{R}(\hat{x},\ell|k):=p(\ell|k)\hat{x}(\ell|k)$ provides a risk upper
bound of a linear spreading process

$\hat{x}(\ell+1|k)=A(\ell|k)\hat{x}(\ell|k),\quad\hat{x}(0|k)=x(k),$ (11)

where the solution $\hat{x}$ bounds the state trajectory $x$ of nonlinear
system (2) with the same initial condition, i.e.
$\hat{x}(0|k)=x(k)\Rightarrow\hat{x}(\ell|k)\geq x(k+\ell)$ for all $\ell$.
Thus, $\bar{R}(x(k)):=\bar{R}(\hat{x},0|k)$ also serves as a risk bound of
(2). Since the row vector $p$ provides importance measurement of the outbreak
at node $i$, we call it a _priority vector_.

###  III-B Relaxed OCP

We now present a relaxed version of OCP (9) by minimizing the risk bounds over
priority vectors, i.e.,

$\begin{split}\min_{p,\beta,\delta}\quad&\bar{R}_{\epsilon}(p,x):=p(0|k)x(k)+\epsilon\sum_{i}p_{i}(0|k)\\\
\mathrm{s.t.}\quad&p(\ell|k)\geq 0,\;p(\ell|k)\geq C+\alpha
p(\ell+1|k)A(\ell|k),\\\
&\beta_{ij}({\ell-1|k})\geq\beta_{ij}({\ell|k})\geq\underline{\beta}_{ij},\\\
&\delta_{i}(\ell-1|k)\leq\delta_{i}(\ell|k)\leq\overline{\delta}_{i},\\\
&\Gamma(\ell|k)\leq\overline{\Gamma},\;\ell=0,\ldots,L-1\end{split}$ (12)

where $\epsilon>0$ is a small constant. We will give a convex reformulation of
the above problem in Section III-C.

Letting $(p^{*},\beta^{*},\delta^{*})$ be a local optimal solution of (12), we
use $A^{*}$ to denote the matrix obtained by substituting
$\beta^{*},\delta^{*}$ into (3). Since (12) can be reformulated as a convex
problem, $(p^{*},\beta^{*},\delta^{*})$ is also a global minimum. We describe
the properties of optimal solutions via the following propositions. In
particular, Proposition 2 shows that the optimal priority vectors satisfy the
Bellman’s equation. Proposition 3 indicates that both $A^{*}$ and $p^{*}$ are
monotonically decreasing before hitting the lower bound. Proposition 4 implies
that the proposed controller (12) is an “all-in” policy, i.e., distributing
all available resources at each time step.

######  Proposition 2.

The optimal solution of (12) satisfies

$p^{*}(\ell|k)=C+\alpha p^{*}(\ell+1|k)A^{*}(\ell|k)$ (13)

for $\ell=0,\ldots,L-1$, where $p^{*}(L|k)=p^{*}(L-1|k)$.

######  Proposition 3.

If $A(k-1)\gneqq\underline{A}$, then we have

$A^{*}(\ell-1|k)\gneqq A^{*}(\ell|k),\quad p^{*}(\ell-1|k)\gneqq
p^{*}(\ell|k)$ (14)

for $\ell=0,\ldots,\underline{L}-1$, where $\underline{L}\leq L$ is the
largest integer such that $A^{*}(\underline{L}-1|k)\gneqq\underline{A}$.

######  Proposition 4.

If $\Gamma(A(k-1),\underline{A})>\overline{\Gamma}$, then we have

$\Gamma(A(k-1),A^{*}(0|k))=\overline{\Gamma}.$ (15)

We now investigate the closed-loop behavior under the MPC controller (12),
i.e., implementing $\beta^{*}(0|k),\delta^{*}(0|k)$ at time step $k$ and
resolving (12) when $x(k+1)$ is available. We first make the following
assumptions.

######  Assumption 5.

The discount rate $\alpha$ satisfies $\alpha\rho(\overline{A})<1$ where
$\rho(\cdot)$ denotes the matrix spectrum radius.

######  Remark 6.

Assumption 5 implies that the risk (4) is well-defined (i.e. finite) for the
linear spreading process (11), despite that the state $\hat{x}(\ell|k)$ might
be unbounded.

######  Assumption 7.

The following set is nonempty:

$\mathcal{A}:=\\{A\mid\underline{A}\leq
A\leq\overline{A},\;\sum_{i}c_{i}\beta_{ij}<c_{j}\delta_{j},\;\forall
j\in\mathcal{V}\\}.$ (16)

######  Remark 8.

Note that $\mathcal{A}$ is convex and satisfies $A(k)\in\mathcal{A}\Rightarrow
A(k+1)\in\mathcal{A}$ due to (5). Later on we will prove that for any
$A\in\mathcal{A}$ the following inequality holds:

$C-(1-\alpha)p>0,$ (17)

where $p=C(I-\alpha A)^{-1}$, which will be used to show the asymptotic
stability of closed-loop system.

######  Theorem 9.

If Assumption 5 holds, then the proposed OCP (12) is recursively feasible.
Furthermore, if Assumption 7 holds, then $\bar{R}(x(k))>\bar{R}(x(k+1))$ for
all $x(k)\gneqq 0$ and $k\geq K$, where
$K=\left\lceil\frac{\Gamma_{M}}{\overline{\Gamma}}\right\rceil$ with
$\Gamma_{M}=\liminf_{A\in\mathcal{A}}\Gamma(\overline{A},A)$.

######  Remark 10.

Since $\Gamma(\overline{A},A)$ is a convex function and $\mathcal{A}$ is a
convex bounded set, then $\Gamma_{M}$ is well-defined and serves as an
estimate of the maximum resource required for making the risk bound
asymptotically decreasing. We will show that $\Gamma_{M}$ can also be computed
via exponential cone programming.

###  III-C Exponential Cone Programming

In this section we show that the proposed dynamic resource allocation problem
(12) can be reformulated as a convex program, in particular exponential cone
program (ECP). We refer the readers to [22] for details on exponential cone
programming.

First, we introduce new decision variables as follows

$y_{i}^{\ell}=\log p_{i}(\ell),\quad U_{ij}^{\ell}=U(A(\ell-1),A(\ell))$ (18)

with $\ell=0,\ldots,L-1$, where the time step $k$ is omitted for simplicity.
The following convex programming is an equivalent formulation of (12):

$\displaystyle\min_{r,y,U}\;$ $\displaystyle r$ $\displaystyle\mathrm{s.t.}\;$
$\displaystyle\log\left(\sum_{i}\exp\bigl{(}\log(x_{i}+\epsilon)+y_{i}^{0}-r\bigr{)}\right)\leq
0$ (19a)
$\displaystyle\log\bigl{(}\sum_{i}\exp\bigl{(}y_{i}^{\ell+1}-y_{j}^{\ell}+\log
a_{ij}-\sum_{s=0}^{\ell}U_{ij}^{s}\bigr{)}$
$\displaystyle\quad+\exp\bigl{(}y_{j}^{\ell+1}-y_{j}^{\ell}+\log
b_{j}-\sum_{s=0}^{\ell}U_{jj}^{s}\bigr{)}$ $\displaystyle\quad+\exp(\log
c_{j}-y_{j}^{\ell})$
$\displaystyle\quad+\exp\left(y_{j}^{\ell+1}-y_{j}^{\ell}+\log
d_{j}\right)\bigr{)}\leq 0,$ (19b) $\displaystyle\quad\quad
j=1,\ldots,n,\;\ell=0,\ldots,L-1$ $\displaystyle U^{\ell}\geq
0,\;\sum_{\ell=0}^{L-1}U^{l}\leq U(A(-1),\underline{A}),$ (19c)
$\displaystyle\sum_{ij}W_{ij}U_{ij}^{\ell}\leq\overline{\Gamma},\;\ell=0,\ldots,L-1$

where $a_{ij}=\alpha h\beta_{ij}(k-1)$, $b_{j}=\alpha
h(\overline{\Delta}_{j}-\delta_{j}(k-1))$,
$d_{j}=\alpha(1-h\overline{\Delta}_{j})$, and $y^{L}=y^{L-1}$.

######  Remark 11.

Constraint (19a) implies $\bar{R}_{\epsilon}(p,x)\leq e^{r}$ while (19b) is
equivalent to the Bellman inequality (10), see [17] for details. The
$\ell_{1}$-norm constraint in (19c), which is based on our logarithmic
resource model (6), can encourage sparse resource allocation. One can apply
the reweighted $\ell_{1}$ optimization approach [23] to further minimize the
number of nodes and edges with non-zero resource allocated, which is
beneficial when resources cannot be distributed widely.

We also give an ECP for computing $\Gamma_{M}$ in Thm. 9:

$\displaystyle\min_{U}\;$ $\displaystyle\sum_{ij}W_{ij}U_{ij}$
$\displaystyle\mathrm{s.t.}\;$ $\displaystyle 0\leq U\leq
U(\overline{A},\underline{A})$ (20a)
$\displaystyle\log\bigl{(}{\textstyle\sum}_{i}\exp(\log(c_{i}\overline{\beta}_{ij}/(c_{j}\overline{\Delta}_{j}))-U_{ij})$
$\displaystyle\;+\exp(\log((\overline{\Delta}_{j}-\underline{\delta}_{j})/\overline{\Delta}_{j})-U_{jj})+\epsilon_{2}\bigr{)}\leq
0$ (20b)

where $\epsilon_{2}>0$ is a small constant and $\Gamma_{M}$ is the optimal
objective value. The set $\mathcal{A}$ is transformed into Constraints (20a) -
(20b) since $\eqref{eq:Gamma-M-1}\Leftrightarrow\underline{A}\leq
A=f(\overline{A},U)\leq\overline{A}$ and

$\begin{split}\eqref{eq:Gamma-M-2}&\Rightarrow\sum_{i}\frac{c_{i}\overline{\beta}_{ij}e^{-U_{ij}}}{c_{j}\overline{\Delta}_{j}}+\frac{c_{j}(\overline{\Delta}_{j}-\underline{\delta}_{j})e^{-U_{jj}}}{c_{j}\overline{\Delta}_{j}}<1\\\
&\Rightarrow\sum_{i}\frac{c_{i}\beta_{ij}}{c_{j}\overline{\Delta}_{j}}+\frac{c_{j}(\overline{\Delta}_{j}-\delta_{j})}{c_{j}\overline{\Delta}_{j}}<1\\\
&\Rightarrow\sum_{i}c_{i}\beta_{ij}<c_{j}\delta_{j}.\end{split}$

##  IV Proofs

For the sake of simplicity, we use $L_{e}=L-1$ to denote the last step in the
time horizon and omit the dependency of time step $k$ if it does not cause any
confusion.

###  IV-A Proof of Proposition 1

We first show that $\bar{R}(\hat{x},\ell)=p(\ell)\hat{x}(\ell)$ with
$p(\ell)=p(L_{e})$ for all $\ell>L_{e}$ is a risk upper bound of the linear
spreading process (11). Since $\hat{x}(\ell)\geq 0$, Inequality (10) implies

$\bar{R}(\hat{x},\ell)\geq Cx(\ell)+\alpha\bar{R}(\hat{x},\ell+1),\;\forall
0\leq\ell\leq L_{e}.$ (21)

For the case where $\ell>L_{e}$ we have

$\begin{split}\bar{R}(\hat{x},\ell)&=p(L_{e})\hat{x}(\ell)\geq(C+\alpha
p(L_{e}+1)A(L_{e}))\hat{x}(\ell)\\\ &\geq C\hat{x}(\ell)+\alpha
p(L_{e}+1)A(\ell)\hat{x}(\ell)\\\ &=C\hat{x}(\ell)+\alpha
p(\ell+1)A(\ell)\hat{x}(\ell)\\\
&=Cx(\ell)+\alpha\bar{R}(\hat{x},\ell+1)\end{split}$ (22)

where the last inequality follows from (5). By telescoping sum we obtain
$\bar{R}(\hat{x},\ell)\geq\sum_{t=0}^{T-1}\alpha^{t}x(t+\ell)+\alpha^{T}\bar{R}(\hat{x},\ell+T)$
for any $T\in\mathbb{N}$. This further implies $\bar{R}(\hat{x},\ell)\geq
R(\hat{x},\ell)$ as $T\rightarrow\infty$, where
$R(\hat{x},\ell)=\sum_{t=0}^{\infty}\alpha^{t}C\hat{x}(t+\ell)$.

Secondly, we show that $R(\hat{x},\ell)$ also bounds the risk of the nonlinear
system (2). By setting $\hat{x}(0)=x(0)=x$ we can have the following
induction: if $\hat{x}(\ell)\geq x(\ell)$ then

$\begin{split}\hat{x}(\ell+1)&=A(\ell)\hat{x}(\ell)\geq A(\ell)x(\ell)\\\
&\geq[A(\ell)-B(x,\ell)]x(\ell)=x(\ell+1)\end{split}$

where inequalities follows by nonnegative $A,B$. Now we can show that
$R(x)\leq R(\hat{x},0)\leq\bar{R}(\hat{x},0)=\bar{R}(x)$.

###  IV-B Proof of Proposition 2

Suppose that $0\leq s\leq L_{e}$ is the smallest integer such that (13) does
not hold, i.e.,

$p^{*}(s)\gneqq C+\alpha p^{*}(s+1)A^{*}(s).$ (23)

If $s=0$, we construct a new sequence $\tilde{p}$ with
$\tilde{p}(\ell)=p^{*}(\ell)$ for $\ell>0$ and $\tilde{p}(0)=C+\alpha
p^{*}(1)A^{*}(0)\lneqq p^{*}(0)$. We have

$\begin{split}\bar{R}_{\epsilon}(p^{*},x)-\bar{R}_{\epsilon}(\tilde{p},x)&=\sum_{i}(p_{i}^{*}(0)-\tilde{p}_{i}(0))(x_{i}+\epsilon)\\\
&\geq\epsilon\sum_{i}(p_{i}^{*}(0)-\tilde{p}_{i}(0))>0,\end{split}$ (24)

which contradicts with the optimality of $p^{*}$.

For the case of $s\geq 1$, we construct a sequence $\tilde{p}$ with
$\tilde{p}(\ell)=p^{*}(\ell)$ for $s<\ell\leq L_{e}$ and

$\tilde{p}(\ell)=C+\alpha\tilde{p}(\ell+1)A^{*}(\ell),\quad\ell=s,\ldots,0.$
(25)

By comparing $p^{*}$ and $\tilde{p}$ we have

$\begin{split}p^{*}(\ell)-\tilde{p}(\ell)=\alpha[p^{*}(\ell+1)-\tilde{p}(\ell+1)]A(\ell),\;\forall
0\leq\ell<s.\end{split}$

Recursively applying the above relationship yields

$p^{*}(0)-\tilde{p}(0)=\alpha^{s}[p^{*}(s)-\tilde{p}(s)]\prod_{j=1}^{s}A(s-j).$
(26)

By the construction we have $p^{*}(s)\gneqq\tilde{p}(s)$. Moreover, since
$A(\ell)\geq 0$ has positive diagonal, we can conclude that
$p^{*}(0)\gneqq\tilde{p}(0)$. The rest is same as the case of $s=0$.

###  IV-C Proof of Proposition 3

For the case where $\underline{L}<L$, we can treat (12) as an OCP with horizon
of $\underline{L}$ as no resource can be distributed at time step
$\ell\geq\underline{L}$. Therefore, we only give the proof for
$\underline{L}=L$.

We first prove that $A^{*}(\ell-1)\gneqq A^{*}(\ell)$ for all $0\leq\ell\leq
L_{e}$ by contradiction. Let $0\leq s\leq L_{e}$ be the largest integer such
that $A^{*}(s-1)=A^{*}(s)$, i.e., there is no resource allocated at time step
$s$. We will construct feasible sequences $\tilde{A}$ and $\tilde{p}$ which
yield a smaller cost than $\bar{R}_{\epsilon}(p^{*},x)$.

We construct $\tilde{A}$ by allocating certain amount of additional resources
to certain link/node at the time step $s$. Since
$A^{*}(L_{e})\gneqq\underline{A}$, without loss of generality we assume that
$\beta_{ij}^{*}(L_{e})>\underline{\beta}_{ij}$ for some link
$(i,j)\in\mathcal{E}$. We now define a new sequence

$\tilde{\beta}_{ij}(\ell)=\begin{cases}\beta_{ij}^{*}(\ell),&0\leq\ell<s\\\
\eta\beta_{ij}^{*}(\ell),&s\leq\ell\leq L_{e}\\\ \end{cases}$ (27)

where
$\eta=\max\left(e^{-\overline{\Gamma}},\frac{\underline{\beta}_{ij}}{\beta_{ij}^{*}(L_{e})}\right)$.
Note that $\tilde{\beta}_{ij}$ is decreasing and satisfies
$\tilde{\beta}_{ij}(L_{e})\geq\underline{\beta}_{ij}$. The resource allocation
does not change except the time step $s$, which satisfies

$\tilde{\Gamma}(s)=w_{ij}\log\frac{\beta_{ij}^{*}(s-1)}{\eta\beta_{ij}^{*}(s)}=-w_{ij}\log\eta\leq\overline{\Gamma}.$

Thus, $\tilde{A}$ is a feasible solution.

We now construct a feasible $\tilde{p}$ based on $\tilde{A}$ as follows

$\tilde{p}(\ell)=C+\alpha\tilde{p}(\ell+1)\tilde{A}(\ell),\quad\ell=L_{e}-1,\ldots,0$
(28)

where $\tilde{p}(L_{e})=C\bigl{(}I-\alpha\tilde{A}(L_{e})\bigr{)}^{-1}$. From
(13) and (28) we can obtain

$\begin{split}p^{*}(&L_{e})-\tilde{p}(L_{e})\\\ &\quad=C(I-\alpha
A^{*}(L_{e}))^{-1}-C(I-\alpha\tilde{A}(L_{e}))^{-1}\\\ &\quad=\alpha
C(I-\alpha A^{*}(L_{e}))^{-1}[A^{*}(L_{e})-\tilde{A}(L_{e})]\\\
&\quad\quad\quad(I-\alpha\tilde{A}(L_{e}))^{-1}\gneqq 0.\end{split}$

Furthermore, by induction we have

$\begin{split}p^{*}(\ell)-\tilde{p}(\ell)=&\alpha[p^{*}(\ell+1)A^{*}(\ell)-\tilde{p}(\ell+1)\tilde{A}(\ell)]\\\
=&\alpha p^{*}(\ell+1)(A^{*}(\ell)-\tilde{A}(\ell))+\\\
&\alpha(p^{*}(\ell+1)-\tilde{p}(\ell+1))\tilde{A}(\ell)\\\
&\gneqq\alpha(p^{*}(\ell+1)-\tilde{p}(\ell+1))\tilde{A}(\ell)\gneqq
0\end{split}$ (29)

where $s\leq\ell\leq L_{e}$. For $0\leq\ell<s$, from (13), (27) and (28) we
can obtain

$p^{*}(\ell)-\tilde{p}(\ell)=\alpha(p^{*}(\ell+1)-\tilde{p}(\ell+1))A^{*}(\ell).$
(30)

Recursively applying (29) and (30) yields

$\begin{split}p^{*}&(0)-\tilde{p}(0)=\\\
&\alpha^{L_{e}}(p^{*}(L_{e})-\tilde{p}(L_{e}))\times\prod_{j=0}^{s-1}A^{*}(j)\times\prod_{j=s}^{L_{e}}\tilde{A}(j)\gneqq
0,\end{split}$

which leads to a contradiction due to (24).

We prove the second part of this proposition, i.e., $p^{*}(\ell-1)\gneqq
p^{*}(\ell)$ for $0\leq\ell\leq L_{e}$. From (13) we have

$\begin{split}p^{*}(L_{e}&-1)-p^{*}(L_{e})\\\
=&\alpha[p^{*}(L_{e})A^{*}(L_{e}-1)-p^{*}(L_{e}+1)A^{*}(L_{e})]\\\ =&\alpha
p^{*}(L_{e})[A^{*}(L_{e}-1)-A^{*}(L_{e})]\gneqq 0\end{split}$

where the second equality is due to $p^{*}(L_{e})=p^{*}(L_{e}+1)$. Now we can
obtain $p^{*}(\ell-1)\gneqq p^{*}(\ell)$ recursively since

$\begin{split}p^{*}(\ell-1)-p^{*}(\ell)&=\alpha[p^{*}(\ell)A^{*}(\ell-1)-p^{*}(\ell+1)A^{*}(\ell)]\\\
&=\alpha[p^{*}(\ell)(A^{*}(\ell-1)-A^{*}(\ell))+\\\
&\quad\quad(p^{*}(\ell)-p^{*}(\ell+1))A^{*}(\ell)]\\\
&\gneqq\alpha(p^{*}(\ell)-p^{*}(\ell+1))A^{*}(\ell),\;\forall\ell<L_{e}.\end{split}$

###  IV-D Proof of Proposition 4

First, the resource allocated at time step $\ell$ is

$\Gamma^{*}(\ell)=\Gamma(A^{*}(\ell-1),A^{*}(\ell))=\sum_{ij}W_{ij}U_{ij}^{*}(\ell),$
(31)

where $W$ is the resource weighting matrix. Supposing that
$\gamma=\overline{\Gamma}-\Gamma^{*}(0)>0$, we aim to construct a new solution
that leads to a lower risk bound. Let $\underline{L}\in[0,L]$ be the largest
integer such that $A^{*}(\underline{L}-1)\gneqq\underline{A}$. Note that such
$\underline{L}$ exists; otherwise we have
$\Gamma(A^{*}(-1),A^{*}(0))=\Gamma(A^{*}(-1),\underline{A})>\overline{\Gamma}$.
Similar to the proof of Proposition 3, we only consider the case where
$\underline{L}=L$.

First, we construct a nonnegative matrix $\widehat{U}$ with
$\widehat{U}_{ij}=0$ for $i\neq j$ and $(i,j)\notin\mathcal{E}$ such that

$f(A^{*}(L-1),\widehat{U})\geq\overline{A}\quad\text{and}\quad\sum_{ij}W_{ij}\widehat{U}_{ij}\leq\gamma^{\prime}$

where $\gamma^{\prime}=\min\\{\Gamma(A^{*}(L_{e}),\underline{A}),\gamma\\}$.
Then we consider a new sequence $\widetilde{U}$ with
$\widetilde{U}(0)=U^{*}(0)+\widehat{U}$ and $\widetilde{U}(\ell)=U^{*}(\ell)$
for $\ell=1,\ldots,L_{e}$. Then, a sequence $\tilde{A}$ can be constructed via

$\tilde{A}(\ell)=f(\tilde{A}(\ell-1),\widetilde{U}(\ell))=f(A^{*}(\ell),\widehat{U})$
(32)

where $\tilde{A}(-1)=A^{*}(-1)$. It is easy to verify that $\tilde{A}$ is a
feasible solution since

$\begin{split}\tilde{\Gamma}(0)&=\Gamma^{*}(0)+\sum_{ij}W_{ij}\widehat{U}_{ij}=\overline{\Gamma}-\gamma+\gamma^{\prime}\leq\overline{\Gamma},\\\
\tilde{\Gamma}(\ell)&=\Gamma^{*}(\ell)\leq\overline{\Gamma},\quad\ell=1,\ldots,L_{e}.\end{split}$

By construction (32) we have $A^{*}(\ell)\gneqq\tilde{A}(\ell)$. Thus, similar
technique in the proof of Proposition 3 can be used to obtain a contradiction.

###  IV-E Proof of Theorem 9

#### Recursive feasibility

Assumption 5 implies that $I-\alpha A(-1)$ has a nonnegative inverse [20, p.
113]. Since $A$ is non-increasing, we have $\rho(A(k-1))\leq\rho(A(k))$ ([20,
p. 15]) and thus $(I-\alpha A(k))^{-1}\gneqq 0$ for all $k$.

Suppose that OCP (12) is feasible at the time step $k$. We construct a
candidate sequence for the time step $(k+1)$ as follows

$\begin{split}A(\ell|k+1)&=\begin{cases}A^{*}(\ell+1|k),&0\leq\ell\leq L-2\\\
A^{*}(L-1|k),&\ell=L-1\end{cases}\\\
p(\ell|k+1)&=\begin{cases}p^{*}(\ell+1|k),&0\leq\ell\leq L-2\\\
p^{*}(L-1|k),&\ell=L-1\end{cases}\end{split}$ (33)

i.e., one time-step shift of the previous optimal solution. Note that $p$ is a
nonnegative sequence and the sequence $A$ satisfies the non-increasing
constraint (5) and resource budget $\Gamma(\ell|k+1)\leq\overline{\Gamma}$.
Thus, $(p,A)$ is a feasible solution.

#### Stability

We first prove that Assumption 7 implies (17) since

$\begin{split}\eqref{eq:CA}&\Rightarrow\alpha
h(c_{j}\delta_{j}-{\textstyle\sum}_{i}c_{i}\beta_{ij})>0\\\ &\Rightarrow\alpha
C(I-A)>0\\\ &\Rightarrow C(I-\alpha A)>(1-\alpha)C\\\ &\Rightarrow
C-(1-\alpha)C(I-\alpha A)^{-1}>0\Rightarrow\eqref{eq:stability}.\end{split}$
(34)

Then, we show $A^{*}(0|k)\in\mathcal{A}$ for $k\geq K$. From Proposition 4 we
have

$\Gamma(\overline{A},A^{*}(0|k))=(k+1)\overline{\Gamma}>\liminf_{A\in\mathcal{A}}\Gamma(\overline{A},A).$

This also implies that Condition (17) holds for $p^{*}(0|k)$ with $k\geq K$.
Finally, the risk bounds are asymptotically decreasing for $k\geq K$ and
$x(k)\gneqq 0$ since

$\begin{split}\bar{R}(x(k))&-\bar{R}(x(k+1))\\\
&=p^{*}(0|k)x(k)-p^{*}(0|k+1)x(k+1)\\\ &\geq
p^{*}(0|k)x(k)-p^{*}(1|k)x(k+1)\\\
&\geq(p^{*}(0|k)-p^{*}(1|k)A^{*}(0|k))x(k)\\\
&=\frac{1}{\alpha}[C-(1-\alpha)p^{*}(0|k)]x(k)>0.\end{split}$ (35)

##  V An Illustrative Example

Let us consider an example of a wildfire which can spread on a fictional
landscape with different area types, as shown in Fig. 1. We represent this
landscape as a graph $\mathcal{G}(\mathcal{V},\mathcal{E})$ with a total
number of nodes $n=1000$. We consider that the fire can spread in vertical,
horizontal, and diagonal direction i.e., the edge set $\mathcal{E}$ connects
each node to a maximum of 8 neighbouring nodes. We also consider the precise
location of a fire outbreak as shown in Fig. 1. For different vegetation type,
the rates of spreading are determined according to the wildfire models in [4]
and [24]. The spreading rate of an edge is given by
$\beta=\beta_{b}\beta_{veg}\beta_{w}$, where $\beta_{b}=0.5$ is the baseline,
and $\beta_{veg}=0.1,~{}1,~{}1.4$ for desert, grassland and eucalyptus forest,
respectively. Water is considered unburnable, hence those edges are removed.
Following [24]$,\beta_{w}$ is calculated for a northeasterly wind with a wind
speed of $4m/s$. The spreading rate is adjusted, following [4], for a spread
between the diagonally connected nodes. For all nodes $i\in\mathcal{V}$, the
recovery rate is set as $\delta=0.5$.

Figure 1: Fictional landscape with different type of areas and the location of
fire outbreak.

We implement our proposed MPC approach to minimize the risk associated with
the outbreak by allocating budgeted resources at each time step. We set the
weightings $w_{ij}~{}=~{}1,~{}\forall~{}(i,j)\in\mathcal{E}$ which indicates
an equal cost to apply resources to any edge. We assign a higher cost of
$c_{j}=1$ for all the city nodes, which could reflect either a higher risk to
human life or a higher economic cost of fire reaching the city. The cost for
the rest of the nodes is set as $c_{j}=0.001$. The discount rate is set as
$\alpha~{}=~{}1/(0.05+\rho(A))$ such that Assumption 5 holds.

The performance of our MPC method with $L=1$ and different resource bounds
$\overline{\Gamma}=10,~{}20,$ and $30$, is shown in Fig. 2. It shows the
evolution of risk $R(x(k))$ given in (4) as well as the risk bound
$\bar{R}(\hat{x},0|k)$ (see Proposition 1), as the number of time steps
increases. It can be seen that for all three cases, our MPC method is able to
reduce the risk associated with the outbreak as time increases. However, for a
smaller resource budget of $\overline{\Gamma}=10$, the risk tends to increase
first which reflects the unavailability of sufficient resources to suppress
the risk at initial time steps. On the other hand, with
$\overline{\Gamma}=30$, the risk starts to decrease almost instantly. Note
that in all three case, the value of risk $R(x(k))$ is always less than the
risk bound $\bar{R}(\hat{x},0|k)$. It also shows the estimated value of $K$
such that for $k\geq K$, the risk bound decreases asymptotically (see Theorem
9 and Remark 10). Note that, for instance $\overline{\Gamma}=10$, the
$K-$estimate is $836$ but the value of risk and the risk bound start to
decrease asymptotically after $k=400$.

Figure 2: Evolution of risk $R(x(k))$, the upper bound on risk $\bar{R}(p,x)$
and the value of $K-$estimate for different values of resource bound
$\overline{\Gamma}$.

Fig. 3 illustrates the effect of increasing $L$ on the evolution of the risk
$R(x(k))$ and the upper bound on risk $\bar{R}(p,x)$ while we set the resource
budget as $\overline{\Gamma}=40$ and we consider that the fire outbreak
affects $25\%$ of nodes at time $k=0$. It indicates that a larger value of $L$
performs better by resulting in a smaller value of risk at the initial time
steps.

Figure 3: Effect of increasing $L$ on risk $R(x(k))$ and the risk bound
$\bar{R}(p,x)$, with $\overline{\Gamma}=40$ and $25\%$ of nodes being infected
at time $k=0$.

## References

  * [1] 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.
  * [2] 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.
  * [3] L. Zino and M. Cao, “Analysis, prediction, and control of epidemics: A survey from scalar to dynamic network models,” _IEEE Circuits and Systems Magazine_ , vol. 21, no. 4, pp. 4–23, 2021.
  * [4] I. Karafyllidis and A. Thanailakis, “A model for predicting forest fire spreading using cellular automata,” _Ecological Modelling_ , vol. 99, no. 1, pp. 87–97, 1997.
  * [5] N. T. Bailey _et al._ , _The mathematical theory of infectious diseases and its applications_. Charles Griffin, 1975.
  * [6] P. Van Mieghem and R. van de Bovenkamp, “Accuracy criterion for the mean-field approximation in susceptible-infected-susceptible epidemics on networks,” _Physical Review E_ , vol. 91, no. 3, p. 032812, 2015.
  * [7] P. Van Mieghem, J. Omic, and R. Kooij, “Virus spread in networks,” _IEEE/ACM Transactions On Networking_ , vol. 17, no. 1, pp. 1–14, 2008.
  * [8] V. M. Preciado, M. Zargham, C. Enyioha, A. Jadbabaie, and G. Pappas, “Optimal vaccine allocation to control epidemic outbreaks in arbitrary networks,” in _52nd IEEE conference on decision and control_. IEEE, 2013, pp. 7486–7491.
  * [9] 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.
  * [10] C. Nowzari, V. M. Preciado, and G. J. Pappas, “Optimal resource allocation for control of networked epidemic models,” _IEEE Transactions on Control of Network Systems_ , vol. 4, no. 2, pp. 159–169, 2015.
  * [11] V. L. Somers and I. R. Manchester, “Sparse resource allocation for control of spreading processes via convex optimization,” _IEEE Control Systems Letters_ , vol. 5, no. 2, pp. 547–552, 2020.
  * [12] ——, “Minimizing the risk of spreading processes via surveillance schedules and sparse control,” _IEEE Transactions on Control of Network Systems_ , 2022\.
  * [13] F. Sélley, A. Besenyei, I. Z. Kiss, and P. L. Simon, “Dynamic control of modern, network-based epidemic models,” _SIAM Journal on applied dynamical systems_ , vol. 14, no. 1, pp. 168–187, 2015.
  * [14] J. Köhler, C. Enyioha, and F. Allgöwer, “Dynamic resource allocation to control epidemic outbreaks a model predictive control approach,” in _2018 Annual American Control Conference (ACC)_. IEEE, 2018, pp. 1546–1551.
  * [15] J. Köhler, L. Schwenkel, A. Koch, J. Berberich, P. Pauli, and F. Allgöwer, “Robust and optimal predictive control of the covid-19 outbreak,” _Annual Reviews in Control_ , vol. 51, pp. 525–539, 2021.
  * [16] N. J. Watkins, C. Nowzari, and G. J. Pappas, “Robust economic model predictive control of continuous-time epidemic processes,” _IEEE Transactions on Automatic Control_ , vol. 65, no. 3, pp. 1116–1131, 2019.
  * [17] V. L. Somers and I. R. Manchester, “Multi-stage sparse resource allocation for control of spreading processes over networks,” in _2022 American Control Conference (ACC)_. IEEE, 2022, pp. 3632–3639.
  * [18] M. ApS, _The MOSEK optimization toolbox for MATLAB_ , 2022.
  * [19] A. Rantzer, “Scalable control of positive systems,” _European Journal of Control_ , vol. 24, pp. 72–80, 2015.
  * [20] A. Berman and R. J. Plemmons, _Nonnegative matrices in the mathematical sciences_. SIAM, 1994.
  * [21] V. L. Somers and I. R. Manchester, “Priority maps for surveillance and intervention of wildfires and other spreading processes,” in _2019 International Conference on Robotics and Automation (ICRA)_. IEEE, 2019, pp. 739–745.
  * [22] S. Boyd, S.-J. Kim, L. Vandenberghe, and A. Hassibi, “A tutorial on geometric programming,” _Optimization and engineering_ , vol. 8, no. 1, pp. 67–127, 2007.
  * [23] E. J. Candes, M. B. Wakin, and S. P. Boyd, “Enhancing sparsity by reweighted $\ell_{1}$ minimization,” _Journal of Fourier analysis and applications_ , vol. 14, no. 5, pp. 877–905, 2008.
  * [24] A. Alexandridis, D. Vakalis, C. I. Siettos, and G. V. Bafas, “A cellular automata model for forest fire spread prediction: The case of the wildfire that swept through spetses island in 1990,” _Applied Mathematics and Computation_ , vol. 204, no. 1, pp. 191–201, 2008.