Alexandre Jacquillat, Michael Lingzhi Li, Martin Ramé, Kai Wang Sloan School
of Management, Massachusetts Institute of Technology, Cambridge, MA
Harvard Business School, Harvard University, Cambridge, MA
Operations Research Center, Massachusetts Institute of Technology, Cambridge,
MA
School of Vehicle and Mobility, Tsinghua University, Beijing, China
Jacquillat, Li, Ramé and Wang

Branch-and-price for prescriptive contagion analytics

Branch-and-price for prescriptive contagion analytics

Predictive contagion models are ubiquitous in epidemiology, social sciences,
engineering, and management. This paper formulates a prescriptive contagion
analytics model where a decision-maker allocates shared resources across
multiple segments of a population, each governed by continuous-time dynamics.
We define four real-world problems under this umbrella: vaccine distribution,
vaccination centers deployment, content promotion, and congestion mitigation.
These problems feature a large-scale mixed-integer non-convex optimization
structure with constraints governed by ordinary differential equations,
combining the challenges of discrete optimization, non-linear optimization,
and continuous-time system dynamics. This paper develops a branch-and-price
methodology for prescriptive contagion analytics based on: (i) a set
partitioning reformulation; (ii) a column generation decomposition; (iii) a
state-clustering algorithm for discrete-decision continuous-state dynamic
programming; and (iv) a tri-partite branching scheme to circumvent non-
linearities. Extensive experiments show that the algorithm scales to very
large and otherwise-intractable instances, outperforming state-of-the-art
benchmarks. Our methodology provides practical benefits in contagion systems;
in particular, it can increase the effectiveness of a vaccination campaign by
an estimated 12-70%, resulting in 7,000 to 12,000 extra saved lives over a
three-month horizon mirroring the COVID-19 pandemic. We provide an open-source
implementation of the methodology in an online repository to enable
replication.

contagion analytics, column generation, branch and price, dynamic programming,
COVID-19

##  1 Introduction

Epidemiological models have played a central role throughout the COVID-19
pandemic. In the United States for instance, the Center for Disease Control
maintained an ensemble forecast based mainly on compartmental contagion models
(Ray et al. 2020). At their core, these models rely on susceptible-infected
(SI) dynamics, which express the number of new infections proportionally to
infected and susceptible individuals. These dynamics have been embedded into
more complex models to capture immunization upon recovery (SIR models), time
lags from exposure to infection (SEIR models), immunization from vaccinations,
asymptomatic cases, quarantine, hospitalization, mortality, etc. Collectively,
these models have been instrumental to guide governmental and societal
response to the COVID-19 pandemic (see Adam 2020, Hsiang et al. 2020, Dehning
et al. 2020, Walker et al. 2020, Li et al. 2022, Bennouna et al. 2022, among
many others).

Contagion models have a long history that predates the COVID-19 pandemic.
Starting from Kermack and McKendrick (1927), SI models have been used to model
infectious epidemics, such as influenza (Casagrandi et al. 2006) and Ebola
(Berge et al. 2017). In marketing, the seminal model of product adoption from
Bass (1969) relies on similar dynamics to capture network externalities
between existing “infected” users and potential “susceptible” adopters. Recent
applications of this model span multi-generational products (Li et al. 2013),
box office revenues (Chance et al. 2008), online content diffusion (Susarla et
al. 2012), etc. Altogether, SI models—and broader dynamical systems—are
ubiquitous in natural sciences, social sciences, engineering and management,
with applications to drug addiction (Behrens et al. 2000), urban congestion
(Saberi et al. 2020), online rumors (Shah and Zaman 2016), climate change
mitigation (Sterman et al. 2012), firm performance (Rahmandad et al. 2018),
employee compensation (Rahmandad and Ton 2020), etc.

Motivated by the success of predictive contagion models, this paper tackles a
prescriptive contagion analytics problem to optimize spatial-temporal resource
allocation over dynamical systems. Specifically, a centralized decision-maker
allocates shared resources across multiple segments of a population, each
governed by continuous-time dynamics. This paper focuses on a deterministic
setting; we show the robustness of the solution under model misspecification
but leave prescriptive contagion analytics problems under uncertainty beyond
the scope of this work. To demonstrate the broad applicability of the model,
we define four problems under its umbrella. The first one involves
distributing a vaccine stockpile to combat an epidemic, inspired by the
COVID-19 pandemic. The second one adds a discrete facility location structure
to deploy mass vaccination centers (Bertsimas et al. 2022). The third one
optimizes online content promotion to maximize product adoption (Lin et al.
2021). The last one deploys emergency vehicles to mitigate urban congestion,
based on a new contagion model of traffic congestion developed in this paper
using real-world data from Singapore.

Prescriptive dynamical systems—in particular, prescriptive contagion
models—involve large-scale mixed-integer non-convex optimization with
constraints governed by ordinary differential equations (ODE). As such, they
are challenging to even formulate, let alone to solve to optimality. The
spatial-temporal resource allocation component, by itself, involves a mixed-
integer optimization structure, but the problem is further complicated by
complexities of dynamical systems:

  1. 1.

Continuous time dynamics: dynamical systems are governed by ODEs of the form
$\frac{d\bm{M}_{i}(t)}{dt}=f_{i}(\bm{M}_{i}(t),\bm{x}_{i}(t))$, where
$\bm{M}_{i}(t)$ denotes a multi-dimensional state variable in segment
$i=1,\cdots,n$; $\bm{x}_{i}(t)$ denotes a multi-dimensional control variable;
and $f_{i}(\cdot,\cdot)$ refers to a continuous-time transition function. An
easy workaround involves time discretization to approximate the dynamics by
$\bm{M}_{i}(t+\Delta
t)-\bm{M}_{i}(t)=f_{i}(\bm{M}_{i}(t),\bm{x}_{i}(t))\cdot\Delta t$. This
approach, however, can lead to extensive computational requirements if the
time increment $\Delta t$ is too small, or to large approximation errors if it
is too large—or both. This challenge is particularly salient in contagion
systems, which are highly non-linear and therefore sensitive to even small
perturbations in initial conditions. A well-known incarnation of these non-
linearities lies in a disease spreading in a population if the basic
reproduction number $R_{0}$ satisfies $R_{0}>1$ but not if $R_{0}<1$.

  2. 2.

Non-linear interactions: contagion systems are driven by non-convex bilinear
interactions of the form $\frac{dS_{i}(t)}{dt}=-\alpha S_{i}(t)I_{i}(t)$,
where $S_{i}(t)$ and $I_{i}(t)$ refer to the susceptible and infected
populations. More broadly, dynamical systems may involve a non-convex
transition function $f_{i}(\cdot,\bm{x}_{i}(t))$ for a given control variable
$\bm{x}_{i}(t)$. So, even with time discretization, the optimization problem
would exhibit a mixed-integer non-convex structure. Accordingly, prescriptive
contagion models have typically been solved via approximations and heuristics
(see Section 2), but no exact method has been devised for resource allocation
in non-linear dynamical systems.

  3. 3.

Other non-linearities: system dynamics are endogenous to resource allocation
decisions. This dependency can be linear: in epidemiology, for example,
immunizations can be safely assumed to be proportional to vaccination rates.
But interventions can also induce a non-linear transition function
$f_{i}(\bm{M}_{i}(t),\cdot)$ for a state variable $\bm{M}_{i}(t)$ (see, e.g.,
Behrens et al. 2000, in drug addiction). Similarly, dynamical systems may
involve non-convex cost functions. Again, these non-linearities create
significant complexities in the resulting optimization problem.

#### Contributions and outline.

This paper develops a branch-and-price methodology for prescriptive contagion
analytics. This approach separates the coupled discrete resource allocation
decisions in a master problem formulated via mixed-integer linear
optimization, and the segment-specific system dynamics in a pricing problem
formulated via continuous-state dynamic programming. To our knowledge, it
provides the first exact methodology to solve spatial-temporal resource
allocation problems in contagion systems—and, more generally, in non-linear
dynamical systems.

Specifically, we define a spatial-temporal resource allocation model with
endogenous continuous-time non-linear dynamics (Section 3). We consider a
finite-horizon setting with multiple segments of a population, each
constituting a dynamical system governed by non-linear ODEs. The decision-
maker optimizes the discrete allocation of a shared resource in each period,
which impacts the system dynamics in each segment. This model allows a variety
of objective functions and constraints in each segment—even non-linear
ones—along with coupling polyhedral constraints across segments (e.g., a
shared budget). We apply this formulation to define our four problems: vaccine
distribution, deployment of vaccination centers, content promotion, and
congestion mitigation. In particular, a byproduct of this research is a new
data-driven contagion model of traffic congestion that yields significant
improvements in predictive performance against state-of-the-art benchmarks.

Next, this paper develops a branch-and-price algorithm to solve the mixed-
integer non-linear optimization problems with ODE constraints (Section 4). The
algorithm relies on four components:

  1. 1.

A set partitioning reformulation. This formulation uses composite variables to
select a resource allocation plan in each segment for the full planning
horizon, as opposed to optimizing natural resource allocation decisions for
each period. By pre-processing the system dynamics into plan-based variables,
this formulation eliminates the three complexities of prescriptive contagion
analytics: continuous time dynamics, non-linear interactions, and other non-
linearities. However, it comes at the cost of an exponential number of plan-
based variables.

  2. 2.

A scalable column generation scheme to solve its linear relaxation. A master
problem solves a coupled resource allocation problem based on a subset of
plan-based variables, via mixed-integer linear optimization. A pricing problem
adds new plans of negative reduced cost or proves that none exists. Thanks to
the problem’s structure, the pricing problem can be decomposed into
independent segment-specific dynamic programming models. Yet, the system
dynamics lead to a continuous state space—a notorious challenge in dynamic
programming.

  3. 3.

A state-clustering algorithm for discrete-decision continuous-state dynamic
programming. Forward-enumeration backward-induction algorithms can solve the
pricing problem but remain intractable in even small instances. Thus, we
propose a linear-time clustering algorithm that exploits the concentration of
states in dynamical systems without using the value function or the policy
function (which varies from iteration to iteration in column generation). The
algorithm provides guarantees on the $\ell_{\infty}$-diameter of each cluster,
which, we prove, controls the global approximation error. The reduction in the
size of the state space considerably enhances the scalability of the pricing
problem at limited costs in terms of approximation errors.

  4. 4.

A novel tri-partite branching scheme to circumvent the non-linearities of the
system. We embed the column generation procedure into a branch-and-price
structure to restore the integrality of resource allocation decisions. We
adopt the typical approach of branching on natural resource allocation
variables, as opposed to composite plan-based variables. Because of the non-
linear dynamics, however, integral resource allocation decisions may not map
into equivalent integral plan-based variables. We therefore develop a tri-
partite branching scheme that retains a natural branching structure, while
guaranteeing finite convergence and optimality.

Finally, this paper demonstrates the scalability of our methodology to
otherwise-intractable prescriptive contagion analytics problems (Section 5).
For resource allocation problems (vaccine allocation, content promotion, and
congestion mitigation), the algorithm returns provably optimal or near-optimal
solutions in manageable computational times, significantly outperforming
state-of-the-art benchmarks. For instance, our algorithm solves vaccine
allocation instances involving 21 decisions in each of 51 regions and 12
weeks, with $\mathcal{O}(21^{612})$ possible decisions. The vaccinations
centers problem features an even more challenging facility location structure
with linking constraints; still, the branch-and-price methodology generates
high-quality solutions and strong bounds in manageable computational times.
From a technical standpoint, the methodology significantly outperforms state-
of-the-art optimization benchmarks,. Specifically, off-the-shelf
implementation based on mixed-integer quadratic optimization or
discretization-based linearization do not scale to even small instances;
moreover, our methodology provides significant benefits as compared to a
tailored coordinate descent heuristic used to circumvent the bilinearities of
the problem. From a practical standpoint, the methodology can have a
significant impact on the management of contagion-based systems. For instance,
the optimized solution can increase the effectiveness of a vaccination
campaign by 12-70% compared to epidemiological benchmarks, resulting in 7,000
to 12,000 extra saved lives over a three-month horizon mirroring the COVID-19
pandemic. These benefits are found to be consistent across contagion analytics
problems and robust to model misspecification. Ultimately, our prescriptive
contagion analytics approach can deliver significant practical benefits in a
variety of domains, by fine-tuning resource allocations based on spatial-
temporal system dynamics.

##  2 Literature review

#### Prescriptive contagion analytics.

The prevalence of contagion models has motivated prescriptive methods to
optimize interventions in contagion-based systems. In epidemiology for
instance, Goldman and Lightwood (2002) showed the benefits of social planning
interventions against user equilibrium behaviors. Rowthorn et al. (2009)
found, using a two-region SIS model, that medical interventions should target
regions with fewer infections, whereas Ndeffo Mbah and Gilligan (2011)
prescribed, using a SIRS model, to prioritize regions with more infections.
Yamin and Gavious (2013) designed incentivization schemes for influenza
vaccinations by embedding SI-based dynamics into an equilibrium model of
vaccination behaviors. In the context of COVID-19, several models traded off
health impacts and economic costs to design differentiated lockdowns across
age tranches (Acemoglu et al. 2021), and to optimize the timing, duration and
intensity of lockdowns (Caulkins et al. 2020, Alvarez et al. 2021, Balderrama
et al. 2022). Similarly, the Bass model has been used to guide operational and
marketing decisions surrounding product innovation, such as sales planning and
time-to-market (Ho et al. 2002); pricing, production and inventory (Shen et
al. 2014); multi-product pricing (Li 2020); free access and premium
subscriptions (Mai and Hu 2022); dynamic pricing (Cosguner and Seetharaman
2022, Zhang et al. 2022b, Agrawal et al. 2021); etc. The Bass model has also
been used to design drug prevention and treatment policies (Behrens et al.
2000), advertising campaigns (Krishnan and Jain 2006), crowdfunding campaigns
(Zhang et al. 2022a), etc. Methodologically, this literature relies on dynamic
programming and control to characterize optimal interventions in a single
dynamical system (or two systems).

In contrast, our paper considers spatial-temporal resource allocation
decisions with coupling constraints across multiple contagion systems. A
seemingly related problem is the ventilator sharing problem from Mehrotra et
al. (2020) and Bertsimas et al. (2021). However, ventilator availability does
not impact the dynamics of the pandemic, so this problem could be formulated
via mixed-integer linear optimization by decoupling upstream epidemiological
predictions from downstream ventilator allocation decisions. In sharp
contrast, vaccinations impact contagion dynamics, which require to integrate
epidemiological dynamics into prescriptive resource allocation models.

Thus, our problem involves discrete optimization with endogenous contagion
dynamics, combining the mixed-integer optimization difficulties of resource
allocation and the non-convex dynamics of contagion models. Existing methods
for this class of problems rely on heuristics and approximations. In product
adoption, Alban et al. (2022) used a heuristic to optimize the deployment of
mobile healthcare units; and Lin et al. (2021) designed an approximation
algorithm with a $1-1/e$ guarantee. In epidemiology, Long et al. (2018) used
myopic linear optimization and approximate dynamic programming to allocate
Ebola treatment units. Bertsimas et al. (2022) optimized where to open
COVID-19 mass vaccination facilities and how to allocate vaccines in the
United States, using a coordinate descent heuristic. Fu et al. (2021) solved a
robust vaccine allocation problem with uncertainty in epidemiological
forecasts, using linear approximations based on discretized numbers of
infections or McCormick reformulations (Fu et al. 2021). Our paper tackles a
deterministic prescriptive contagion analytics problem, and contributes an
exact optimization approach that does not rely on finite difference
approximations of continuous-time ODE dynamics and that does not involve
approximations and heuristics to handle the non-convexities of the problem.

#### Mixed-integer non-linear optimization.

Our problem exhibits a mixed-integer non-linear optimization (MINLO) structure
with ODE constraints,. It therefore remains highly challenging even with
finite difference approximations of ODEs (see Burer and Letchford 2012, for a
detailed review of MINLO). Bilinear SI-based models fall into mixed-integer
quadratic optimization. Methods include the reformulation-linearization
technique (Sherali and Adams 2013), semi-definite relaxations (Fujie and
Kojima 1997, Anstreicher 2009), disjunctive programming (Saxena et al. 2010,
2011), perspective reformulations (Günlük and Linderoth 2010, Anstreicher and
Burer 2021), augmented Lagrangian duality (Gu et al. 2020), etc. General-
purpose MINLO methods include spatial branch-and-bound (Lee and Grossmann
2001), branch-and-reduce (Ryoo and Sahinidis 1996, Tawarmalani and Sahinidis
2004), and $\alpha$-branch-and-bound (Androulakis et al. 1995). Recent
versions of Gurobi Optimization (2020) tackle mixed-integer quadratic
optimization problems by combining RLT and spatial branch-and-bound, which we
use as a benchmark in Section 5.

The ODE constraints also link to the mixed-integer optimal control literature.
When the feasible region can be enumerated, partial outer approximation and
sum-up rounding can converge to the optimal solution as time discretization
becomes infinitely granular (Sager et al. 2012, Hante and Sager 2013, Manns
and Kirches 2020). In practice, however, very granular time discretization may
be required, and optimality is not guaranteed in the presence of control costs
or combinatorial constraints. Jung et al. (2015) proposed a continuous optimal
control problem and a combinatorial integral approximation problem to restore
integer feasibility; Göttlich et al. (2021) separated a mixed-integer optimal
control problem from time-coupling combinatorial constraints; and Bestehorn et
al. (2021) restored integral feasibility via a shortest path formulation. Our
paper departs from this literature by considering large-scale problems with
spatial coupling across multiple dynamical systems, and by proposing a new
branch-and-price methodology.

Since its introduction by Barnhart et al. (1998), branch-and-price has been
widely applied to mixed-integer linear optimization. In non-linear
optimization, Andreas et al. (2008) solved a reliable $h$-paths problem with
non-linearities stemming from failure probabilities. Nowak et al. (2018)
proposed an inner-and-outer-approximation of MINLO problems by combining
column generation with non-linear optimization. Allman and Zhang (2021)
synthesized a generic branch-and-price algorithm for non-convex mixed-integer
optimization. Our paper proposes a tailored branch-and-price decomposition for
spatial-temporal resource allocation problems over continuous-time dynamical
systems—including, notably, a novel tri-partite branching disjunction to
handle non-linearities.

Finally, one of the main bottlenecks of our methodology lies in the discrete-
decision continuous-state dynamic programming structure of the pricing problem
(see, e.g., Bertsekas 2015, Powell 2022, for reviews of approximate dynamic
programming and reinforcement learning). Continuous-state problems are
typically handled via state discretization. Bertsekas (1975) showed that
static state discretization converges to the optimum as the grid becomes
increasingly granular. To manage the growth in the number of states, adaptive
discretization methods build a finer grid where the cost function changes
rapidly (Grüne and Semmler 2004, Borraz-Sánchez and Haugland 2011). Several
studies combined adaptive discretization with reinforcement learning to
aggregate states with similar cost-to-go functions (see, e.g. Bertsekas et al.
1988, Pyeatt et al. 2001, Lolos et al. 2017, Sinclair et al. 2022). Bennouna
et al. (2021) used transition data to learn a discrete partition of a
continuous state space. In our problem, however, these methods are not readily
applicable because the cost function changes at each column generation
iteration. Instead, our state-clustering approach exploits stability in the
state space—as opposed to stability in the cost-to-go function.

##  3 Prescriptive contagion analytics: definition and formulation

###  3.1 Model formulation

We consider a general problem of spatial-temporal resource allocation in
dynamical systems, depicted in Figure 1. The “spatial” component refers to
resource allocation across $n$ segments of a population (e.g., across regions
in vaccine allocation, across products in content promotion). The temporal
component refers to $S$ decision epochs throughout a planning horizon of
length $T$. We denote by $\tau_{s}$ the time stamp of epoch $s=1,\cdots,S$,
with $\tau_{1}=0$ and $\tau_{S+1}=T$.

Resourceallocation$\bm{x}_{11}$$\bm{x}_{12}$$\bm{x}_{13}$$\cdots$$\bm{x}_{1,S-1}$$\bm{x}_{1,S}$Segment
$1$System
dynamics$\cdots$$\bm{x}_{n1}$$\bm{x}_{n2}$$\bm{x}_{n3}$$\cdots$$\bm{x}_{n,S-1}$$\bm{x}_{n,S}$Segment
$n$System dynamics Figure 1: Schematic representation of the spatial-temporal
resource allocation problem in dynamical systems.

We characterize centralized resource allocation decisions via variables
$\bm{x}_{is}\in\mathcal{F}_{is}\subseteq\mathbb{R}^{d_{is}}$ for each segment
$i=1,\cdots,n$ and each epoch $s=1,\cdots,S$. These variables are defined as
$d_{is}$-dimensional vectors to allow for the allocation of multiple resources
(e.g., treatment vs. prevention vehicles in congestion mitigation). We also
introduce a variable $\bm{y}\in\mathbb{Z}^{q}\times\mathbb{R}^{r}$ to capture
other decisions (e.g., facility location in our vaccination centers problem).
Let $\Delta_{s}=\sum_{i=1}^{n}d_{is}$; let
$\bm{X}_{s}\in\mathbb{R}^{\Delta_{s}}$ denote the concatenated resource
allocation variable at epoch $s=1,\cdots,S$. Resource allocation in segment
$i=1,\cdots,n$ at epoch $s=1,\cdots,S$ comes at a cost
$\Gamma_{is}(\bm{x}_{is})$, and other decisions come at a linear cost
$\bm{d}^{\top}\bm{y}$. We make no restriction on the feasible regions
$\mathcal{F}_{is}$ and the cost functions $\Gamma_{is}(\cdot)$ governing
segment-specific resource allocations, thus allowing non-linear decision
problems. Importantly, we focus on discrete resource allocation decisions,
where each region $\mathcal{F}_{is}$ has finite cardinality
$D_{is}=|\mathcal{F}_{is}|<\infty$. In many cases, resources are naturally
restricted to discrete quantities (e.g., emergency vehicles, in our congestion
mitigation problem). Even when the decision space is infinite, operational
constraints often lead to a discretized set of possible decisions. For
example, Moderna (2020) shipped COVID-19 vaccines in pallets of around 20,000
vaccines. Still, the resource allocation problem remains high-dimensional due
to coupled spatial-temporal resource allocations—if $D_{is}=\mathcal{O}(D)$
for all $i$, $s$, the decision space grows in $\mathcal{O}(D^{n\times S})$,
which quickly becomes very large.

Besides segment-specific constraints captured in $\mathcal{F}_{is}$, global
constraints are characterized by a polyhedral set
$\left\\{(\bm{X},\bm{y})\in\mathbb{R}^{\Delta_{s}+q+r}:\bm{U}_{s}\bm{X}+\bm{V}_{s}\bm{y}\geq\bm{w}_{s}\right\\}$,
where $\bm{U}_{s}\in\mathbb{R}^{m_{s}\times\Delta_{s}}$,
$\bm{V}_{s}\in\mathbb{R}^{m_{s}\times(q+r)}$ and
$\bm{w}_{s}\in\mathbb{R}^{m_{s}}$. In other words, resource allocation
decisions are subject to the following linear constraints, where
$\bm{u}_{sji}\in\mathbb{R}^{d_{is}}$ denotes the vector corresponding to the
$j^{\text{th}}$ constraint and the $i^{\text{th}}$ segment in $\bm{U}_{s}$,
and $\bm{v}_{sj}\in\mathbb{R}^{q+r}$ denotes the vector corresponding to the
$j^{\text{th}}$ constraint in $\bm{V}_{s}$:

$\left(\sum_{i=1}^{n}\bm{u}_{sji}^{\top}\bm{x}_{is}\right)+\bm{v}_{sj}^{\top}\bm{y}\geq
w_{sj},\ \forall s=1,\cdots,S,\ \forall j=1,\cdots,m_{s}$

Next, each population segment $i=1,\cdots,n$ constitutes a dynamical system,
governed by a continuous-time state variable
$\bm{M}_{i}(t)\in\mathbb{R}^{r_{i}}$. We assume that the dynamics are
independent across the $n$ segments. In the epidemiological context, for
instance, this captures intra-region interactions (e.g., within each state or
each country) but not inter-region interactions (see, e.g., Hsiang et al.
2020, Walker et al. 2020, Li et al. 2022, Bennouna et al. 2022, for similar
assumptions). Specifically, in each segment $i=1,\cdots,n$, the state variable
$\bm{M}_{i}(t)$ is determined by an initial condition $\bm{M}^{0}_{i}$, and
varies according to the following system of ODEs, where $f_{i}(\cdot,\cdot)$
denotes the transition function:

$\frac{d\bm{M}_{i}(t)}{dt}=f_{i}(\bm{M}_{i}(t),\bm{x}_{is}),\ \forall
s=1,\cdots,S,\ \forall t\in\left[\tau_{s},\tau_{s+1}\right].$

The system is associated with a continuous-time cost function
$g_{it}(\bm{M}_{i}(t))$ and a terminal cost function $h_{i}(\bm{M}_{i}(T))$.
We impose no restriction on the transition function $f_{i}(\cdot,\cdot)$ and
the cost functions $g_{it}(\cdot)$ and $h_{i}(\cdot)$ so, by design, our model
encompasses non-convex dynamical systems.

The spatial-temporal resource allocation problem, referred to as Problem
$(\mathcal{P})$, minimizes total costs subject to resource allocation
constraints and the system’s dynamics. It is written as follows:

$\displaystyle(\mathcal{P})\quad\min\quad$
$\displaystyle\sum_{i=1}^{n}\left(\int_{0}^{T}g_{it}(\bm{M}_{i}(t))dt+h_{i}(\bm{M}_{i}(T))\right)+\sum_{i=1}^{n}\sum_{s=1}^{S}\Gamma_{is}(\bm{x}_{is})+\bm{d}^{\top}\bm{y}$
(1) s.t.
$\displaystyle\left(\sum_{i=1}^{n}\bm{u}_{sji}^{\top}\bm{x}_{is}\right)+\bm{v}_{sj}^{\top}\bm{y}\geq
w_{sj},\ \forall s=1,\cdots,S,\ \forall j=1,\cdots,m_{s}$ (2)
$\displaystyle\frac{d\bm{M}_{i}(t)}{dt}=f_{i}(\bm{M}_{i}(t),\bm{x}_{is}),\
\forall i=1,\cdots,n,\ \forall s=1,\cdots,S,\ \forall
t\in\left[\tau_{s},\tau_{s+1}\right]$ (3)
$\displaystyle\bm{M}_{i}(\tau_{1})=\bm{M}^{0}_{i},\ \forall i=1,\cdots,n$ (4)
$\displaystyle\bm{x}_{is}\in\mathcal{F}_{is},\ \forall i=1,\cdots,n,\ \forall
s=1,\cdots,S$ (5) $\displaystyle\bm{y}\in\mathbb{Z}^{q}\times\mathbb{R}^{r}$
(6)

This formulation, however, is intractable due to continuous time dynamics
(Equation (3)), non-convex system dynamics (functions $f_{i}(\cdot,\cdot)$,
$g_{it}(\cdot)$ and $h_{i}(\cdot)$), and non-convex segment-wise decisions
(functions $\Gamma_{is}(\cdot)$, regions $\mathcal{F}_{is}$). One workaround
would be to model $(\mathcal{P})$ via dynamic programming, by defining a state
variable that encompasses all variables
$\bm{M}_{1}(\tau_{s}),\cdots,\bm{M}_{n}(\tau_{s})$ as well as all required
information to enforce resource allocation constraints (Equation (2)). This
approach, however, would scale in $\mathcal{O}(D^{n\times S})$ if
$D_{is}=\mathcal{O}(D)$ for all $i$, $s$, quickly growing into the curse of
dimensionality. Instead, we propose in Section 4 a branch-and-price
methodology that separates coupled resource allocation decisions across
segments and continuous-time dynamics in each segment. We still use dynamic
programming in the pricing problem, but we leverage segment-specific state
variables $\bm{M}_{i}(t)\in\mathbb{R}^{r_{i}}$ rather than a higher-
dimensional state variable of the form
$(\bm{M}_{1}(\tau_{s}),\cdots,\bm{M}_{n}(\tau_{s}))\in\mathbb{R}^{r_{1}}\times\cdots\times\mathbb{R}^{r_{n}}$.
Thus, the pricing problem scales in $\mathcal{O}(D^{S})$, instead of
$\mathcal{O}(D^{n\times S})$.

Finally, let us underscore that, following Powell (2022) we define the state
variable $\bm{M}_{i}(t)$ as a necessary and sufficient function of history to
compute the cost function, the constraints, and the transition function in
segment $i=1,\cdots,n$. For example, an inter-temporal budget constraint of
the form $\sum_{s=1}^{S}\bm{x}_{is}\leq\bm{\gamma}_{i}$ would require state
augmentation to capture the budget used up to epoch $s$. In contrast, our
model can accommodate budget constraints across segments of the form
$\sum_{i=1}^{n}\bm{x}_{is}\leq\bm{\beta}_{s}$. Although the problem is not
directly modeled as a dynamic program, this structure enables segment-wise
dynamic programming decomposition in our algorithm.

###  3.2 Prescriptive contagion analytics model, and applications

The special case of contagion systems is derived from the general formulation
by defining, in each segment $i=1,\cdots,n$, a susceptible state
$S_{i}(t)\in\mathbb{R}$, an infected state $I_{i}(t)\in\mathbb{R}$, and other
states represented via a vector $\bm{R}_{i}(t)\in\mathbb{R}^{V}$. We capture
bilinear interactions between susceptible and infected populations as follows,
with infection rate $\alpha$ and transitions $f^{S}_{i}(\cdot)$,
$f^{I}_{i}(\cdot)$, $f^{R}_{i}(\cdot)$:

$\displaystyle\frac{dS_{i}(t)}{dt}=-\alpha
S_{i}(t)I_{i}(t)+f^{S}_{i}(S_{i}(t),I_{i}(t),\bm{R}_{i}(t))$ (7)
$\displaystyle\frac{dI_{i}(t)}{dt}=+\alpha
S_{i}(t)I_{i}(t)+f^{I}_{i}(S_{i}(t),I_{i}(t),\bm{R}_{i}(t))$ (8)
$\displaystyle\frac{dR_{iv}(t)}{dt}=f^{R}_{i}(S_{i}(t),I_{i}(t),\bm{R}_{i}(t)),\quad\forall
v=1,\cdots,V$ (9)

Similarly, we define functions $g^{SIR}_{it}(\cdot)$ and $h^{SIR}_{i}(\cdot)$
to characterize the costs of the system dynamics in each segment. We then
formulate the prescriptive contagion model by minimizing the following
objective function, subject to the resource allocation constraints (Equation
(2) and Equations (5)–(6)) and the contagion dynamics (initial conditions, and
Equations (7)–(9)):

$\displaystyle\sum_{i=1}^{n}\left(\int_{0}^{T}g^{SIR}_{it}(S_{i}(t),I_{i}(t),\bm{R}_{i}(t))dt+h^{SIR}_{i}(S_{i}(T),I_{i}(T),\bm{R}_{i}(T))\right)+\sum_{i=1}^{n}\sum_{s=1}^{S}\Gamma_{is}(\bm{x}_{is})+\bm{d}^{\top}\bm{y}$
(10)

$S$$E$$I$$U$$H$$Q$$D$$S’$$E’$$I’$$M$$\beta x$ (a) Vaccine allocation and
vaccination centers (Bertsimas et al. 2022).

$A$$B$$\alpha Ax^{2}+\beta BA$ (b) Content promotion (Lin et al. 2021).

$S$$W$$A$$I$$A^{\prime}$$R$$\beta^{F}S$$\rho^{F}(\zeta^{F}+\psi(x^{1}))A$$\alpha^{F}S(S+A+A^{\prime})$$\alpha^{W}W(S+A+A^{\prime})$$\beta^{W}W$$\rho^{W}(\zeta^{W}+\psi(x^{1}))A^{\prime}$$(1-\rho^{S})(\zeta^{F}+\varphi(x^{2}))$$(1-\rho^{W})(\zeta^{W}+\varphi(x^{2}))$$\theta
I$ (c) Urban congestion (this paper, Appendix 7.2).

Figure 2: Contagion models for the four problems. Thick red lines indicate
transitions that are endogenous to resource allocation. Dotted lines indicate
bilinear susceptible-infected interactions.

We define our four problems from this generic formulation, all inspired from
real-world problems and data-driven contagion models. We defer details and
their full formulation to Appendix 7.1.

##### Vaccine allocation:

The problem optimizes the number of vaccines ($x$) received in 51 regions (50
US states plus Washington, DC) to maximize the number of saved lives. We use
the DELPHI-V model of COVID-19 from Li et al. (2022) and Bertsimas et al.
(2022), with 7 non-vaccinated states (susceptible (S), exposed (E), infected
(I), undetected (U), hospitalizations (H), quarantine (Q), death (D)) and 4
vaccinated states (susceptible (S’), exposed (E’), infected (I’), immunized
(M)). Our experiments consider allocating a weekly stockpile of 2.5 to 7
million vaccines across 51 states and 12 weeks. We discretize decisions into
up to 21 decisions per epoch, with $\mathcal{O}\left(21^{612}\right)$ possible
decisions. This setup mirrors the situation from February to April 2021 in the
United States.

##### Vaccination centers.

The problem extends the vaccine allocation problem by jointly optimizing the
location of mass vaccination centers and the subsequent distribution of
vaccines within a geographical radius (Bertsimas et al. 2022). Due to the
complexity of the problem, we break down the United States into 10 groups
defined by the Centers for Disease Control, with 4-7 states and 20 candidate
facilities in each (Appendix 9.1). Each group receives a budget of vaccines
proportionally to its share of the US population, discretized into
25,000-vaccine pallets, plus some flexibility buffer. We set a budget of 20
and 30 facilities across the country, again mirroring the situation in 2021.
With 20 (resp. 30) facilities, each group needs to select 2 (resp. 3)
facilities, and residents can travel up to 150 miles (resp. 100 miles),
possibly across states, to access a vaccination center.

##### Content promotion.

This problem optimizes product promotion to maximize adoption, subject to a
sparsity constraint to promote up to $K$ products per period, a cover
constraint, a linking constraint, and contagion dynamics (Lin et al. 2021).
The decision variables includes which products should be promoted ($x^{1}$
variables) and the share of the population to which each product shall be
promoted ($x^{2}$ variables). Each segment corresponds to a product, and
adoption follows a Bass model with susceptible users (A) and adopters (B). In
our experiments, we consider up to 20 products, 10 decision epochs, 21
decisions per epoch, and 2-6 products promoted per epoch.

##### Congestion mitigation.

This problem deploys $B^{1}$ prevention vehicles and $B^{2}$ treatment
vehicles per hour to mitigate congestion costs. Prevention resources ($x^{1}$)
respond to minor accidents and prevent congestion, whereas treatment resources
($x^{2}$) respond to major accidents and ease recovery. Contagion models
reflect the propagation of congestion from “infected” roads to free-flow
“susceptible” roads (Saberi et al. 2020). We leverage new data sources from
Singapore on traffic speeds, road work, and traffic accidents to propose a
richer six-state model, (Figure 2c): susceptible roads (S), road work (W),
accidents (A), road work and accidents ($A^{\prime}$), congested (I), and
recovered (R). We calibrate this model using a three-peak structure for
morning, afternoon, and evening traffic. A byproduct of this paper is a new
contagion model of urban congestion, which improves predictive performance
against state-of-the-art benchmarks (Appendix 7.2). In our experiments, we
consider a setting with 5 neighborhoods, up to 8 hours and up to 12 vehicles
to allocate, developed in collaboration with city officials to capture real-
world situations in Singapore.

##### Discussion.

The vaccine allocation problem has the largest state space. The vaccination
centers problem adds a discrete structure with a switching variable $\bm{y}$
(facility location) and linking constraints governing $\bm{x}_{is}$ (vaccine
allocation). The other two problems include multi-dimensional variables
$\bm{x}_{is}$, and the congestion mitigation problem involves a novel
contagion model. Together, they cover a range of applications and optimization
structures in prescriptive contagion analytics.

##  4 The method: branch-and-price for prescriptive contagion analytics

We propose a branch-and-price algorithm to solve Problem $(\mathcal{P})$,
using a set partitioning reformulation (Section 4.1), a column generation
procedure (Section 4.2), a state-clustering algorithm for the continuous-state
pricing problem (Section 4.3), and a tri-partite branching scheme to handle
non-convexities (Section 4.4). Section 4.5 summarizes the algorithm and
establishes its exactness.

###  4.1 Set partitioning reformulation

The reformulation optimizes over composite plan-based variables in each
segment (which characterize resource allocations over the full planning
horizon), as opposed to natural variables (which characterize resource
allocation in each period). Specifically, we define the set of feasible plans
in each segment $i=1,\cdots,n$ as the combination of feasible resource
allocation decisions, as follows:

$\mathcal{P}_{i}=\mathcal{F}_{i1}\times\cdots\times\mathcal{F}_{iS},\ \forall
i=1,\cdots,n$

In each segment $i=1,\cdots,n$, a plan $p\in\mathcal{P}_{i}$ defines a
sequence of resource allocations denoted by
$(\bm{\alpha}^{p}_{i1},\cdots,\bm{\alpha}^{p}_{iS})\in\mathcal{F}_{i1}\times\cdots\times\mathcal{F}_{iS}$,
as well as a cost parameter $C^{p}_{i}$ governed by the system’s dynamics:

$\displaystyle C^{p}_{i}=$
$\displaystyle\quad\int_{0}^{T}g_{it}(\bm{M}_{i}(t))dt+h_{i}(\bm{M}_{i}(T))+\sum_{s=1}^{S}\Gamma_{is}(\bm{x}_{is})$
(11) s.t.
$\displaystyle\quad\frac{d\bm{M}_{i}(t)}{dt}=f_{i}(\bm{M}_{i}(t),\bm{\alpha}^{p}_{is}),\
\forall s=1,\cdots,S,\ \forall t\in\left[\tau_{s},\tau_{s+1}\right]$ (12)
$\displaystyle\quad\bm{M}_{i}(\tau_{1})=\bm{M}^{0}_{i}$ (13)

We define the following plan-based decision variables:

$\displaystyle z^{p}_{i}=\begin{cases}1&\text{if plan $p\in\mathcal{P}_{i}$ is
selected for segment $i=1,\cdots,n$,}\\\ 0&\text{otherwise.}\end{cases}$

The set partitioning formulation, referred to as $(\mathcal{S}\mathcal{P})$,
minimizes total costs (Equation (14)), while enforcing the coupling
constraints (Equation (15)) and ensuring that one plan is selected in each
segment (Equation (16)). For completeness, we reformulate our four problems in
Appendix 7.3.

$\displaystyle(\mathcal{S}\mathcal{P})\qquad\min\quad$
$\displaystyle\sum_{i=1}^{n}\sum_{p\in\mathcal{P}_{i}}C^{p}_{i}z^{p}_{i}+\bm{d}^{\top}\bm{y}$
(14) s.t.
$\displaystyle\left(\sum_{i=1}^{n}\sum_{p\in\mathcal{P}_{i}}\bm{u}_{sji}^{\top}\bm{\alpha}^{p}_{is}z^{p}_{i}\right)+\bm{v}_{sj}^{\top}\bm{y}\geq
w_{sj},\ \forall s=1,\cdots,S,\ \forall j=1,\cdots,m_{s}$ (15)
$\displaystyle\sum_{p\in\mathcal{P}_{i}}z^{p}_{i}=1\quad\forall i=1,\cdots,n$
(16) $\displaystyle z^{p}_{i}\in\\{0,1\\}\quad\forall i=1,\cdots,n,\ \forall
p\in\mathcal{P}_{i}$ (17)
$\displaystyle\bm{y}\in\mathbb{Z}^{q}\times\mathbb{R}^{r}$ (18)

###### Proposition 4.1

The set partitioning formulation $(\mathcal{S}\mathcal{P})$ is equivalent to
Problem $(\mathcal{P})$.

This reformulation eliminates structural complexities of Problem
$(\mathcal{P})$ by pre-processing the continuous-time dynamics and the non-
linear functions into plan-based variables and the corresponding input
parameters. The $(\mathcal{S}\mathcal{P})$ formulation therefore exhibits a
mixed-integer linear optimization structure, but involves an exponential
number of composite plan-based variables. In an instance with $\mathcal{O}(D)$
decisions in each segment and at each epoch, the number of plans scales in
$\mathcal{O}(D^{n\times S})$. It is therefore intractable to even enumerate
the full set of plans, let alone to estimate the cost parameters by running a
dynamical model for each one (Equations (11)–(13)).

###  4.2 Solving the linear optimization relaxation via column generation

To address this challenge, we generate plans iteratively via column
generation. A restricted master problem (RMP) solves the linear relaxation of
$(\mathcal{S}\mathcal{P})$ over subsets of plans
$\mathcal{P}^{0}_{i}\subseteq\mathcal{P}_{i}$:

$\displaystyle(RMP)\qquad\min\quad$
$\displaystyle\sum_{i=1}^{n}\sum_{p\in\mathcal{P}^{0}_{i}}C^{p}_{i}z^{p}_{i}+\bm{d}^{\top}\bm{y}$
(19) s.t.
$\displaystyle\left(\sum_{i=1}^{n}\sum_{p\in\mathcal{P}^{0}_{i}}\bm{u}_{sji}^{\top}\bm{\alpha}^{p}_{is}z^{p}_{i}\right)+\bm{v}_{sj}^{\top}\bm{y}\geq
w_{sj},\ \forall s=1,\cdots,S,\ \forall j=1,\cdots,m_{s}$ (20)
$\displaystyle\sum_{p\in\mathcal{P}^{0}_{i}}z^{p}_{i}=1\quad\forall
i=1,\cdots,n$ (21) $\displaystyle z^{p}_{i}\geq 0\quad\forall i=1,\cdots,n,\
\forall p\in\mathcal{P}^{0}_{i}$ (22) $\displaystyle\bm{y}\in\mathbb{R}^{q+r}$
(23)

Let $\lambda_{sj}\in\mathbb{R}_{+}$ be the dual variable of the coupling
constraint (Equation (20)) and $\mu_{i}\in\mathbb{R}$ the dual variable of the
set partitioning constraint (Equation (21)). The pricing problem exploits
segment-wise decomposition to generate new plan-based variables or a
certificate that none exists. Specifically, it seeks the plan in
$\mathcal{P}_{i}$ with the minimal reduced cost, for each segment
$i=1,\cdots,n$:

$\displaystyle(PP_{i})\qquad\min\quad$
$\displaystyle\left(\int_{0}^{T}g_{it}(\bm{M}_{i}(t))dt+h_{i}(\bm{M}_{i}(T))+\sum_{s=1}^{S}\Gamma_{is}(\bm{x}_{is})-\sum_{s=1}^{S}\sum_{j=1}^{m_{s}}\lambda_{sj}\bm{u}_{sji}^{\top}\bm{x}_{is}-\mu_{i}\right)$
(24) s.t. $\displaystyle\bm{x}_{is}\in\mathcal{F}_{is},\ \forall s=1,\cdots,S$
(25)
$\displaystyle\frac{d\bm{M}_{i}(t)}{dt}=f_{i}(\bm{M}_{i}(t),\bm{x}_{is}),\
\forall s=1,\cdots,S,\ \forall t\in\left[\tau_{s},\tau_{s+1}\right]$ (26)
$\displaystyle\bm{M}_{i}(\tau_{1})=\bm{M}^{0}_{i}$ (27)

The pricing problem exhibits a discrete non-linear optimization structure with
continuous-time dynamics. However, we can leverage temporal decomposition to
formulate it via dynamic programming. By definition, the continuous-time state
variable $\bm{M}_{i}(t)$ encapsulates all necessary system history to compute
the cost function, the constraints in $\mathcal{F}_{is}$, and the transition
function $f_{i}(\cdot,\cdot)$, so the system follows Markovian dynamics. At
each epoch $s=1,\cdots,S$, the dynamic programming state variable is denoted
by $\bm{N}^{i}_{s}$, equal to $\bm{M}_{i}(\tau_{s})$. The decision is
$\bm{x}_{is}\in\mathcal{F}_{is}$, and the transition function is governed by
the continuous-time dynamics between $\tau_{s}$ and $\tau_{s+1}$, with initial
condition determined by the state variable. The problem involves a per-period
cost characterizing the continuous-time cost $g_{it}(\cdot)$ accrued between
$\tau_{s}$ and $\tau_{s+1}$, the cost of resource allocation
$\Gamma_{is}(\cdot)$, and the dual price of the coupling constraints. In
addition, the problem involves a terminal cost characterizing the end-state
cost $h_{i}(\cdot)$ and the dual price of the set partitioning constraints.
The Bellman equation is given as follows, using $J_{s}(\cdot)$ to define the
cost-to-go function at epoch $s=1,\cdots,S$:

$\displaystyle J_{s}\left(\bm{N}^{i}_{s}\right)$
$\displaystyle=\min_{\bm{x}_{is}\in\mathcal{F}_{is}}\left\\{\int_{\tau_{s}}^{\tau_{s+1}}g_{it}(\bm{M}_{i}(t))dt+\Gamma_{is}(\bm{x}_{is})-\sum_{j=1}^{m_{s}}\lambda_{sj}\bm{u}_{sji}^{\top}\bm{x}_{is}+J_{s+1}\left(\bm{N}^{i}_{s+1}\right)\right\\}$
(28) where $\bm{M}_{i}(\tau_{s})\leftarrow\bm{N}^{i}_{s}$;
$\frac{d\bm{M}_{i}(t)}{dt}=f_{i}(\bm{M}_{i}(t),\bm{x}_{is}),\ \forall
t\in\left[\tau_{s},\tau_{s+1}\right]$; and
$\bm{N}^{i}_{s+1}\leftarrow\bm{M}_{i}(\tau_{s+1})$ $\displaystyle
J_{S+1}\left(\bm{N}^{i}_{S+1}\right)$
$\displaystyle=h_{i}\left(\bm{N}^{i}_{S+1}\right)-\mu_{i}$ (29)

The column generation procedure solves the linear relaxation of the
$\mathcal{S}\mathcal{P}$ formulation to optimality. Starting with initial sets
$\mathcal{P}^{0}_{i}$ of plan-based variables, the RMP provides a feasible
primal solution at each iteration, along with the dual variables
$\lambda_{sj}$ and $\mu_{i}$. We then solve $(PP_{i})$ for each segment
$i=1,\cdots,n$; if its optimal value is negative, we expand the set
$\mathcal{P}^{0}_{i}$ with the new solution and proceed. Otherwise, the
incumbent RMP solution is optimal for the $(\mathcal{S}\mathcal{P})$
relaxation.

This column generation scheme separates the two complexities of the problem:
the RMP handles the coupled resource allocation decisions via mixed-integer
optimization (Equations (19)–(22)), and the pricing problem handles the
continuous-time system dynamics via dynamic programming (Equations (28)–(29)).
However, the pricing problem exhibits a continuous-state dynamic programming
structure—a notoriously challenging class of problems. This motivates our
state-clustering dynamic programming algorithm to solve it efficiently at each
column generation iteration.

###  4.3 A state-clustering dynamic programming algorithm for the pricing
problem

Since the pricing problem is separable across segments $i=1,\cdots,n$, we
ignore the index $i$ in this section; for instance, $\bm{N}_{s}$ refers to the
state, and $\bm{x}_{s}\in\mathcal{F}_{s}$ to the decision with
$D_{s}=|\mathcal{F}_{s}|$.

#### Exact algorithm.

Any finite-horizon discrete-decision continuous-state dynamic programming
model can be solved via forward enumeration and backward induction, as
follows:

  1. 1.

Forward enumeration: From the initial state $\bm{M}^{0}$, evaluate all
possible decisions $\bm{x}_{1},\cdots,\bm{x}_{S}$, and store all possible
states at each epoch $s=1,\cdots,S+1$ in an exhaustive set
$\mathcal{N}^{*}_{s}$.

  2. 2.

Backward induction: Starting from $s=S$, derive the optimal policy
$\bm{\pi}_{s}^{*}(\cdot)$ at epoch $s=S,\cdots,1$ and update the cost-to-go
function $J_{s}(\cdot)$, as follows:

$\displaystyle\bm{\pi}_{s}^{*}\left(\bm{N}_{s}\right)$
$\displaystyle=\argmin_{\bm{x}_{s}\in\mathcal{F}_{s}}\left\\{\int_{\tau_{s}}^{\tau_{s+1}}g_{t}(\bm{M}(t))dt+\Gamma_{s}(\bm{x}_{s})-\sum_{j=1}^{m_{s}}\lambda_{sj}\bm{u}_{sj}^{\top}\bm{x}_{s}+J_{s+1}\left(\bm{N}_{s+1}\right)\right\\}$
(30) $\displaystyle J_{s}\left(\bm{N}_{s}\right)$
$\displaystyle=\int_{\tau_{s}}^{\tau_{s+1}}g_{t}(\bm{M}(t))dt+\Gamma_{s}\left(\bm{\pi}_{s}^{*}\left(\bm{N}_{s}\right)\right)-\sum_{j=1}^{m_{s}}\lambda_{sj}\bm{u}_{sj}^{\top}\bm{\pi}_{s}^{*}\left(\bm{N}_{s}\right)+J_{s+1}\left(\bm{N}_{s+1}\right)$
(31) where $\bm{M}(\tau_{s})\leftarrow\bm{N}_{s}$;
$\frac{d\bm{M}(t)}{dt}=f(\bm{M}(t),\bm{\pi}_{s}^{*}\left(\bm{N}_{s}\right)),\
\forall t\in\left[\tau_{s},\tau_{s+1}\right]$; and
$\bm{N}_{s+1}\leftarrow\bm{M}(\tau_{s+1})$

This approach is detailed in Algorithm 2 in Appendix 8.2. Per the Bellman
principle of optimality, the optimal path is optimal at every decision point
(Proposition 4.2).

###### Proposition 4.2

Algorithm 2 returns a feasible and optimal solution to Equations (28)–(29).

However, Algorithm 2 is highly computationally expensive due to exhaustive
state enumeration in the forward-enumeration and backward-induction loops.
Complexity grows in $\mathcal{O}\left(\prod_{s=1}^{S}D_{s}\right)$, or
$\mathcal{O}(D^{S})$ if $D_{s}=\mathcal{O}(D),\ \forall s$, which can quickly
become very large as the planning horizon and the decision space increase.
Even for a simple eight-period problem with two decisions that can each take
four values, $\prod_{s=1}^{S}D_{s}=16^{8}\sim 1$ billion, which severely
hinders exact dynamic programming.

As noted earlier, most approximate dynamic programming and reinforcement
learning methods depend on learning approximations of the cost-to-go function
$J_{s}(\bm{N}_{s})$ and/or the optimal policy function
$\bm{\pi}_{s}^{*}(\bm{N}_{s})$. In our problem, however, both functions change
from one column generation iteration to the next, due to the updates in the
dual values $\lambda_{sj}$ from the RMP. Thus, standard approaches would need
be to re-learned at every iteration, resulting in significant inefficiencies.

Instead, we propose a state-clustering algorithm that exploits the
concentration of states, but not the cost function. This approach proceeds by
aggregation to reduce the state space during the forward-enumeration step of
dynamic programming. Note that forward enumeration does not depend on the cost
function, so the state clustering algorithm needs to be applied only once at
the beginning of the column generation algorithm. Then, we leverage the
backward-induction algorithm with the clustered state space to speed up the
pricing problem at each column generation iteration.

#### A state-clustering acceleration.

Many dynamical systems exhibit a natural concentration of states. In vaccine
allocation, for example, there are only mild differences between a no-vaccine
baseline and an allocation of just a few vaccines; similarly, two allocations
that merely swap decisions between epochs $s$ and $s+1$ may yield similar
future outcomes. We leverage this observation to approximate the full set of
states $\mathcal{N}_{s}^{*}$ by a clustered state space $\mathcal{N}_{s}$ at
each epoch $s=1,\cdots,S$.

Specifically, at epoch $s=1,\cdots,S$, we enumerate all pairs of states in
$\mathcal{N}_{s}$ and all decisions in $\mathcal{F}_{s}$; for each one, we
compute the next state at period $s+1$ and the corresponding cost accrued
between times $\tau_{s}$ and $\tau_{s+1}$. This yields
$Q=|\mathcal{N}_{s}|\times|\mathcal{F}_{s}|$ state-decision pairs
$(\bm{N}_{s}^{1},\bm{x}_{s}^{1}),\cdots,(\bm{N}_{s}^{Q},\bm{x}_{s}^{Q})$ with
$Q$ subsequent states $\bm{N}_{s+1}^{1},\cdots,\bm{N}_{s+1}^{Q}$ and
corresponding cost parameters $\chi_{s}^{1},\cdots,\chi_{s}^{Q}$. These are
defined as follows:

$\displaystyle\bm{M}^{q}(\tau_{s})\leftarrow\bm{N}_{s}^{q},\quad\frac{d\bm{M}^{q}(t)}{dt}=f(\bm{M}^{q}(t),\bm{x}_{s}^{q}),\quad\bm{N}^{q}_{s+1}\leftarrow\bm{M}^{q}(\tau_{s+1}),\quad\forall
q=1,\cdots,Q$
$\displaystyle\chi_{s}^{q}=\int_{\tau_{s}}^{\tau_{s+1}}g_{t}(\bm{M}^{q}(t))dt+\mathbf{1}\\{s=S\\}(h(\bm{M}^{q}(T))),\quad\forall
q=1,\cdots,Q$

We group states $\bm{N}_{s+1}^{1},\cdots,\bm{N}_{s+1}^{Q}$ into $\Omega_{s+1}$
clusters, and treat the centroids as representative states in
$\mathcal{N}_{s+1}$ at epoch $s+1$. We then apply the backward-induction
algorithm in the clustered state space. By design, this procedure can reduce
the state space significantly as long as
$\Omega_{s+1}\ll|\mathcal{N}_{s}|\times|\mathcal{F}_{s}|$. Obviously, the
resulting dynamic program is an approximation of the full dynamic program, but
the optimality loss can be small if the full state space is indeed
concentrated.

The next question lies in the design of the clustering algorithm. Algorithms
based on global similarity measures have at least quadratic complexity to
compute the entire distance matrix. In our setup,
$|\mathcal{N}_{s}|\times|\mathcal{F}_{s}|$ is already very large, so a
complexity of
$\mathcal{O}\left(|\mathcal{N}_{s}|^{2}\times|\mathcal{F}_{s}|^{2}\right)$
would be intractable. Most linear-time clustering algorithms, such as
$k$-means, cannot effectively bound the diameter of each cluster, potentially
leading to significant error propagation in non-linear dynamical systems.
Instead, we propose a linear-time clustering algorithm with guarantees on
cluster diameter in an $\ell_{\infty}$-space. Unlike $k$-means, the algorithm
creates clusters dynamically without a pre-specified number of clusters. The
benefits of our approach are showed theoretically in Proposition 4.4 (namely,
the bound on cluster diameter bounds the global approximation error) and
numerically in Section 5 (namely, the smaller global approximation error
induces stronger solutions than a $k$-means benchmark).

Specifically, for each cluster $\omega\in\\{1,\cdots,\Omega_{s+1}\\}$, we
store: (i) the number of data points $\eta(\omega)$ and the sum of all states
$\bm{N}^{\Sigma}(\omega)$ (to define the centroid); (ii) the element-wise
minimum and maximum states, $\underline{\bm{N}}(\omega)$ and
$\overline{\bm{N}}(\omega)$ (to bound the cluster diameter); and (iii) the set
of state-decision pairs $\mathcal{X}(\omega)$ that lead to the cluster and the
total cost $c(\omega)$ (to define cost functions in the clustered state
space). For a diameter tolerance $\varepsilon$, we cluster the states
$\bm{N}_{s+1}^{1},\cdots,\bm{N}_{s+1}^{Q}$ as follows

  1. 1.

We assign the first state $\bm{N}_{s+1}^{1}$ to a new cluster $\omega_{1}$. We
initialize: $\eta(\omega_{1})\leftarrow 1;\
\bm{N}^{\Sigma}(\omega_{1})\leftarrow\bm{N}_{s+1}^{1};\
\underline{\bm{N}}(\omega_{1})\leftarrow\bm{N}_{s+1}^{1};\
\overline{\bm{N}}(\omega_{1})\leftarrow\bm{N}_{s+1}^{1};\
\mathcal{X}(\omega_{1})\leftarrow\\{(\bm{N}_{s}^{1},\bm{x}_{s}^{1})\\};\
c^{\Sigma}(\omega_{1})\leftarrow\chi_{s}^{1}$

  2. 2.

We assign $\bm{N}_{s+1}^{q}$ to a cluster based on its $\ell_{\infty}$
distance to the minimum and maximum states:

     * –

If
$\min_{\ell=1,\cdots,k}\max(\|\bm{N}_{s+1}^{q}-\underline{\bm{N}}(\omega_{\ell})\|_{\infty},\|\bm{N}_{s+1}^{q}-\overline{\bm{N}}(\omega_{\ell})\|_{\infty})\leq\varepsilon$,
then $\bm{N}^{q}$ is assigned to $\omega_{l^{*}}$ where $\ell^{*}$ denotes the
index of the cluster that attains the minimum:

$\displaystyle\eta(\omega_{l^{*}})\leftarrow\eta(\omega_{l^{*}})+1;\
\bm{N}^{\Sigma}(\omega_{l^{*}})\leftarrow\bm{N}^{\Sigma}(\omega_{l^{*}})+\bm{N}_{s+1}^{q};\
\underline{\bm{N}}(\omega_{l^{*}})\leftarrow\min(\underline{\bm{N}}(\omega_{l^{*}}),\bm{N}_{s+1}^{q})$
$\displaystyle\overline{\bm{N}}(\omega_{l^{*}})\leftarrow\max(\overline{\bm{N}}(\omega_{l^{*}}),\bm{N}_{s+1}^{q});\
\mathcal{X}(\omega_{l^{*}})\leftarrow\mathcal{X}(\omega_{l^{*}})\cup\\{(\bm{N}_{s}^{q},\bm{x}_{s}^{q})\\};\
c^{\Sigma}(\omega_{l^{*}})\leftarrow c^{\Sigma}(\omega_{l^{*}})+\chi_{s}^{q},$

     * –

Otherwise, $\bm{N}_{s+1}^{q}$ is assigned to a new cluster $\omega_{k+1}$:
$\eta(\omega_{k+1})\leftarrow 1;\
\bm{N}^{\Sigma}(\omega_{k+1})\leftarrow\bm{N}_{s+1}^{q};\
\underline{\bm{N}}(\omega_{k+1})\leftarrow\bm{N}_{s+1}^{q};\
\overline{\bm{N}}(\omega_{k+1})\leftarrow\bm{N}_{s+1}^{q};\
\mathcal{X}(\omega_{k+1})\leftarrow\\{(\bm{N}_{s}^{q},\bm{x}_{s}^{q})\\};\
c^{\Sigma}(\omega_{k+1})\leftarrow\chi_{s}^{q}$

  3. 3.

We retrieve each cluster’s centroid
$\bm{N}_{s+1}(\omega_{\ell})=\bm{N}^{\Sigma}(\omega_{\ell})/\eta(\omega_{\ell})$
and define the cost function in the clustered state space, as follows:
$c(\bm{N}_{s},\bm{x}_{s})=\frac{c^{\Sigma}(\omega_{\ell})}{\eta(\omega_{\ell})}+\Gamma_{s}(\bm{x}_{s})$
for all $(\bm{N}_{s},\bm{x}_{s})\in\mathcal{X}(\omega_{\ell})$.

We solve a pricing problem approximation in the clustered state space, via
backward induction:

$\displaystyle\bm{\pi}_{s}^{*}\left(\bm{N}_{s}\right)$
$\displaystyle=\argmin_{\bm{x}_{s}\in\mathcal{F}_{s}}\left\\{c(\bm{N}_{s},\bm{x}_{s})-\sum_{j=1}^{m_{s}}\lambda_{sj}\bm{u}_{sj}^{\top}\bm{x}_{s}+J_{s+1}\left(\bm{N}_{s+1}\right)\right\\}$
(32) $\displaystyle J_{s}\left(\bm{N}_{s}\right)$
$\displaystyle=c(\bm{N}_{s},\bm{\pi}_{s}^{*}\left(\bm{N}_{s}\right))-\sum_{j=1}^{m_{s}}\lambda_{sj}\bm{u}_{sj}^{\top}\bm{\pi}_{s}^{*}\left(\bm{N}_{s}\right)+J_{s+1}\left(\bm{N}_{s+1}\right)$
(33) where $\bm{M}(\tau_{s})\leftarrow\bm{N}_{s}$;
$\frac{d\bm{M}(t)}{dt}=f(\bm{M}(t),\bm{\pi}_{s}^{*}\left(\bm{N}_{s}\right)),\
\forall t\in\left[\tau_{s},\tau_{s+1}\right]$; and
$\bm{N}_{s+1}\leftarrow\bm{M}(\tau_{s+1})$

Pseudo-code is given in Algorithm 3 in Appendix 8.2. The algorithm makes a
single pass through all the states, thus terminating in linear time (as long
as the number of clusters is much smaller than the number of states). Still,
the algorithm controls the maximum pair-wise distance within each cluster,
guaranteeing an $\ell_{\infty}$-diameter within $\varepsilon$. This is
formalized in Proposition 4.3.

###### Proposition 4.3

At each decision epoch $s=1,\cdots,S+1$, each cluster
$\omega_{1},\cdots,\omega_{k}$ satisfies:

$\|\overline{\bm{N}}(\omega_{\ell})-\underline{\bm{N}}(\omega_{\ell})\|_{\infty}\leq\varepsilon.$
(34)

This guarantee allows us to bound the global approximation error
proportionally to $\varepsilon$, in Proposition 4.4. Note that the error grows
exponentially in time, which is typical in dynamical systems (Sager et al.
2012). Nonetheless, by controlling the cluster diameter, the state-clustering
algorithm can control the global error. Our results in Section 5.1 show that,
in practice, our state-clustering algorithm induces moderate errors in small
instances and can scale to large instances. Moreover, by controlling the
cluster diameter, it leads to stronger optimization solutions than $k$-means.

###### Proposition 4.4

Assume that the transition function $f(.,\bm{x}_{s})$ is $L_{s}$-Lipschitz in
the state variables in the $l_{\infty}$ metric, for $s\in\\{1,\cdots,S+1\\}$.
For each true state $\bm{N}^{*}_{s}\in\mathcal{N}^{*}_{s}$ (from Algorithm 2),
there exists a clustered state $\bm{N}_{s}\in\mathcal{N}_{s}$ (from Algorithm
3) such that:

$\|\bm{N}^{*}_{s}-\bm{N}_{s}\|_{\infty}\leq\begin{cases}0&\text{if $s=1$}\\\
\varepsilon\sum\limits_{\sigma=3}^{s+1}\exp\left(\sum\limits_{\nu=\sigma}^{s}L_{\nu-1}(\tau_{\nu}-\tau_{\nu-1})\right),&\text{for
all}\ s\geq 2\end{cases}$ (35)

###  4.4 A tri-partite branching disjunction

We embed the column generation procedure into a branch-and-bound structure to
retrieve an integral solution to the $(\mathcal{S}\mathcal{P})$ formulation.
We proceed via multi-phase branching: we first branch on the mixed-integer
variables $\bm{y}$; we then branch on the natural resource allocation
variables $\bm{x}_{is}$; and we finally restore integrality of the plan-based
variables $z^{p}_{i}$ via a tri-partite branching disjunction.

##### Branching on mixed-integer variables.

In the case where the formulation involves discrete variables $\bm{y}$, we
first branch on those variables whenever one of its components is fractional:

$\underbrace{(y_{\ell}\leq\left\lfloor\widehat{y}_{\ell}\right\rfloor)}_{\text{left
branch}}\lor\underbrace{(y_{\ell}\geq\left\lceil\widehat{y}_{\ell}\right\rceil)}_{\text{right
branch}},\ \text{where $\widehat{y}_{\ell}$ denotes the $\ell^{\text{th}}$
component of the incumbent solution}$ (36)

##### Branching on the natural variables.

To avoid building deep and one-sided trees and to maintain the structure of
the pricing problem, a typical approach in branch-and-price involves branching
on the natural variables (i.e., the resource allocation decisions
$\bm{x}_{is}$ in our case) as opposed to branching on the composite plan-based
variables (i.e., the variables $z^{p}_{i}$). Let
$\bm{x}_{is}(\bm{z})\in\mathbb{R}^{d_{is}}$ denote the resource allocation
variable in segment $i=1,\cdots,n$ and epoch $s=1,\cdots,S$ for a plan-based
solution $\bm{z}$; and let $x^{k}_{is}(\bm{z})\in\mathbb{R}$ be its
$k^{\text{th}}$ component, for $k=1,\cdots,d_{is}$. A column generation
solution $\widehat{\bm{z}}$ leads to infeasible resource allocations if
$\bm{x}_{is}(\widehat{\bm{z}})\notin\mathcal{F}_{is}$. In that case, we
exploit the discreteness of the feasible region $\mathcal{F}_{is}$ to create a
valid disjunction. Since $\mathcal{F}_{is}$ does not necessarily comprise
contiguous integers, we introduce dedicated floor and ceiling functions:
$\lfloor\bm{a}\rfloor^{k}_{is}=\max\\{\beta\leq
a^{k}:\exists\bm{b}\in\mathcal{F}_{is},b^{k}=\beta\\}$ and
$\lceil\bm{a}\rceil^{k}_{is}=\min\\{\beta\geq
a^{k}:\exists\bm{b}\in\mathcal{F}_{is}:b^{k}=\beta\\}$ (where $a^{k}$ and
$b^{k}$ are the $k^{\text{th}}$ components of $\bm{a}$ and $\bm{b}$).

Armed with these notations, we can then define the following valid branching
disjunction:

$\underbrace{\bigg{(}x^{k}_{is}(\bm{z})\leq\left\lfloor\bm{x}_{is}(\widehat{\bm{z}})\right\rfloor^{k}_{is}\bigg{)}}_{\text{left
branch}}\lor\underbrace{\bigg{(}x^{k}_{is}(\bm{z})\geq\left\lceil\bm{x}_{is}(\widehat{\bm{z}})\right\rceil^{k}_{is}\bigg{)}}_{\text{right
branch}},\ \text{with}\
\bm{x}_{is}(\widehat{\bm{z}})=\sum_{p\in\mathcal{P}_{i}}\bm{\alpha}^{p}_{is}\widehat{z}^{p}_{i}$
(37)

Out of these disjunctions, we select one that corresponds to a variable
$x^{k}_{is}(\widehat{\bm{z}})$ with the largest value of
$\min\left\\{x^{k}_{is}(\widehat{\bm{z}})-\left\lfloor\bm{x}_{is}(\widehat{\bm{z}})\right\rfloor^{k}_{is},\left\lceil\bm{x}_{is}({\widehat{\bm{z}}})\right\rceil^{k}_{is}-x^{k}_{is}(\widehat{\bm{z}})\right\\}$.
This branching strategy is an analog to the most fractional branching strategy
in integer optimization. Importantly, this disjunction preserves the structure
of the pricing problem, by merely restricting the search from the full
feasible regions $\mathcal{F}_{is}$ to their subregions defined by the
corresponding lower-bound and upper-bound constraints.

##### Tri-partite branching.

By design, the above branching schemes enforce integrality of the $\bm{y}$
variables and ensure that the plan-based variables define feasible resource
allocation decisions $\bm{x}_{is}$. However, they may not guarantee integral
plan-based variables $z^{p}_{i}$; for instance, we may obtain
$\alpha^{p_{1}}_{is}=4$, $\alpha^{p_{2}}_{is}=8$,
$z^{p_{1}}_{i}=z^{p_{2}}_{i}=0.5$, and $6\in\mathcal{F}_{is}$. In linear
problems, this solution can be brought into an equivalent feasible solution by
considering the “average plan” $p^{*}$ with $\alpha^{p^{*}}_{is}=6$. In our
problem, however, the non-convex system dynamics break this equivalence
because $C^{p^{*}}_{i}\neq 0.5C^{p_{1}}_{i}+0.5C^{p_{2}}_{i}$ in general.
Thus, the plan $p^{*}$ is no longer guaranteed to form an optimal solution to
Problem $(\mathcal{P})$.

A direct way to enforce the integrality of the plan-based variables would be
to create disjunctions of the form $(z^{p}_{i}=0)\lor(z^{p}_{i}=1)$. However,
as noted earlier, this disjunction can lead to weak and imbalanced tree
structures. Moreover, it breaks the structure of the pricing problem—notably,
it is difficult to seek the “second-best” solution in the branch associated
with the $z^{p}_{i}=0$ disjunction.

Instead, we devise a novel tri-partite branching disjunction to handle the
non-linearities. Consider a node with a solution $\widehat{\bm{z}}$. Assume
that $\bm{x}_{is}(\widehat{\bm{z}})\in\mathcal{F}_{is}$ for all $i,s$ but that
there exists a fractional plan-based variable $z^{p}_{i}\in(0,1)$. Let
$\alpha^{pk}_{is}(\bm{z})\in\mathbb{R}$ denote the $k^{\text{th}}$ component
of $\bm{\alpha}^{p}_{is}$ for $k=1,\cdots d_{is}$. There must exist a segment
$i$, an epoch $s$, a component $k$ and a plan $p_{0}$ such that:

$\bm{x}_{is}(\widehat{\bm{z}})=\sum_{p\in\mathcal{P}_{i}}\bm{\alpha}^{p}_{is}\widehat{z}^{p}_{i}\in\mathcal{F}_{is}\
\text{and}\ \alpha^{k,p_{0}}_{is}\neq x^{k}_{is}(\widehat{\bm{z}})\
\text{and}\ \widehat{z}^{p_{0}}_{i}>0.$ (38)

Intuitively, $\bm{x}_{is}(\widehat{\bm{z}})$ defines a “promising” resource
allocation decision for segment $i$ at epoch $s$. Accordingly, we create a
corresponding branch in the tree. To retain a mutually exclusive and
collectively exhaustive tree, we create two additional branches with lower and
higher resource allocations. Formally, our tri-partite branching disjunction
is defined via the following disjunction:

$\underbrace{\bigg{(}x^{k}_{is}(\bm{z})<x^{k}_{is}(\widehat{\bm{z}})\bigg{)}}_{\text{left
branch}}\lor\underbrace{\bigg{(}x^{k}_{is}(\bm{z})=x^{k}_{is}(\widehat{\bm{z}})\bigg{)}}_{\text{middle
branch}}\lor\underbrace{\bigg{(}x^{k}_{is}(\bm{z})>x^{k}_{is}(\widehat{\bm{z}})\bigg{)}}_{\text{right
branch}},\ \text{with}\
\bm{x}_{is}(\widehat{\bm{z}})=\sum_{p\in\mathcal{P}_{i}}\bm{\alpha}^{p}_{is}\widehat{z}^{p}_{i}$
(39)

A final question lies in variable selection to determine which
$x^{k}_{is}(\widehat{\bm{z}})$ to branch upon. We extend the most fractional
principle, by selecting the variable with the largest weighted difference
between the plan-based allocations and the resource allocations from the
column generation solution. Formally, we select $x^{k}_{is}(\widehat{\bm{z}})$
such that the $(i,k,s)$ tuple maximizes
$\sum_{p\in\mathcal{P}_{i}}\widehat{z}^{p}_{i}|\alpha^{pk}_{is}-x^{k}_{is}(\widehat{\bm{z}})|$.

Again, this branching disjunction does not impact the structure of the pricing
problem, which can easily accommodate lower-bounding and upper-bounding
constraints on the decision variables. As we shall establish in Theorem 4.5,
this procedure eliminates any infeasible solution satisfying the conditions in
Equation (38) and yields integer plan-based variables $z_{i}^{p}$ upon
termination.

###  4.5 Summary of the branch-and-price algorithm

Algorithm 1 summarizes our solution approach. In each node, the algorithm
solves the linear relaxation of $\mathcal{S}\mathcal{P}$ via column generation
(Step 1), iterating between the RMP (Step 1.1) and the dynamic programming
algorithm for the pricing problem (Step 1.2), until convergence. We also apply
in Step 2 an optional upper-bounding scheme for acceleration: in our resource
allocation problems for instance, we solve a restricted master problem with
integrality constraints; in the vaccination centers problem, this becomes much
more cumbersome, so we fix facility construction variables based on the linear
relaxation solution and solve the subsequent vaccine allocation problem. The
algorithm then proceeds to bi-partite branching when the variables $\bm{y}$
are not integral or when the natural resource allocation variables $\bm{x}$
are infeasible (Steps 3 and 4), and to tri-partite branching to restore
integrality in plan-based variables (Step 5). This overall branching strategy
is illustrated in Figure 3. Step 6 corresponds to the node selection step—we
implemented breadth-first and depth-first search strategies, and found
comparable performance. The algorithm uses bounding and feasibility rules to
prune leaves (Step 2) and terminates when no remaining leaf is active (Step
6).

Solve $\mathcal{S}\mathcal{P}$ relaxation$OBJ=90,\
\bm{x}_{25}(\widehat{\bm{z}})\notin\mathcal{F}_{25}$$OBJ=92$$\bm{x}_{31}(\widehat{\bm{z}})\notin\mathcal{F}_{31}$pt$x^{k}_{25}(\bm{z})\leq\left\lfloor
x^{k}_{25}(\widehat{\bm{z}})\right\rfloor_{25}$$OBJ=91$$\bm{\alpha}^{p_{0}}_{14}\neq\bm{x}_{14}(\widehat{\bm{z}})$pt$x^{k}_{25}(\bm{z})\geq\left\lceil
x^{k}_{25}(\widehat{\bm{z}})\right\rceil_{25}$$OBJ=97$$\bm{x}_{24}(\widehat{\bm{z}})\notin\mathcal{F}_{24}$pt$x^{k}_{31}(\bm{z})\leq\left\lfloor
x^{k}_{31}(\widehat{\bm{z}})\right\rfloor_{31}$$OBJ=93$$\bm{\alpha}^{p_{0}}_{23}\neq\bm{x}_{23}(\widehat{\bm{z}})$pt$x^{k}_{31}(\bm{z})\geq\left\lceil
x^{k}_{31}(\widehat{\bm{z}})\right\rceil_{31}$$OBJ=96$$\bm{\alpha}^{p_{0}}_{33}\neq\bm{x}_{33}(\widehat{\bm{z}})$pt$x^{k}_{14}(\bm{z})<x^{k}_{14}(\widehat{\bm{z}})$$OBJ=95$$z^{p}_{i}\in\\{0,1\\},\forall
i,p$pt$x^{k}_{14}(\bm{z})=x^{k}_{14}(\widehat{\bm{z}})$$OBJ=93$$\bm{x}_{12}(\widehat{\bm{z}})\notin\mathcal{F}_{12}$pt$x^{k}_{14}(\bm{z})>x^{k}_{14}(\widehat{\bm{z}})$$OBJ=95$$\bm{x}_{13}(\widehat{\bm{z}})\notin\mathcal{F}_{24}$pt$x^{k}_{23}(\bm{z})<x^{k}_{23}(\widehat{\bm{z}})$$OBJ=94$$z^{p}_{i}\in\\{0,1\\},\forall
i,p$pt$x^{k}_{23}(\bm{z})=x^{k}_{23}(\widehat{\bm{z}})$$OBJ=95$$\bm{\alpha}^{p_{0}}_{32}\neq\bm{x}_{32}(\widehat{\bm{z}})$pt$x^{k}_{23}(\bm{z})>x^{k}_{23}(\widehat{\bm{z}})$$OBJ=98$$\bm{x}_{35}(\widehat{\bm{z}})\notin\mathcal{F}_{35}$pt$x^{k}_{12}(\bm{z})\leq\left\lceil
x^{k}_{12}(\widehat{\bm{z}})\right\rceil_{12}$$OBJ=99$$\bm{\alpha}^{p_{0}}_{21}\neq\bm{x}_{21}(\widehat{\bm{z}})$pt$x^{k}_{12}(\bm{z})\geq\left\lceil
x^{k}_{12}(\widehat{\bm{z}})\right\rceil_{12}$ Figure 3: Illustration of the
branching structure, combining bi-partite branching and tri-partite branching.
Red nodes are pruned by bound; blue nodes are pruned by feasibility; white
nodes trigger bi-partite branching; blue squares trigger tri-partite
branching; the green node indicates the optimal solution. Algorithm 1 Branch-
and-price algorithm for Problem ($\mathcal{P}$).

1:Initialization: node list $\Phi=\\{0\\}$ of the branch-and-bound tree;
parent node $cur=0$; initial sets of plans
$\mathcal{P}^{0}_{1},\cdots,\mathcal{P}^{0}_{n}$ and corresponding RMP
(19)–(22); upper bound $UB=+\infty$; incumbent solution
$\bm{z}^{*}=\emptyset$.

2:Step 1. Solve parent node $cur$ via column generation (Section 4.2); store
objective $OBJ$, solution $\widehat{\bm{z}}$, and corresponding resource
allocation decisions
$\bm{x}_{is}(\bm{z})=\sum_{p\in\mathcal{P}_{i}}\bm{\alpha}^{p}_{is}\widehat{z}^{p}_{i}$;
remove $cur$ from node list $\Phi$.

3: Step 1.1. Solve RMP relaxation with
$\mathcal{P}^{0}_{1},\cdots,\mathcal{P}^{0}_{n}$; store solution
$\widehat{\bm{z}}$, objective $OBJ$, and dual costs $\bm{\lambda},\bm{\mu}$.

4: Step 1.2. Solve $(PP_{i})$ for each segment $i=1,\cdots,n$, using the exact
dynamic programming algorithm (Algorithm 2) or the state-clustering dynamic
programming algorithm (Algorithm 3). If the objective is negative for some
segment $i$, add optimal plan to $\mathcal{P}^{0}_{i}$, and go to Step 1.1.

5:Step 2. Apply upper-bounding heuristic (e.g., restricted master problem with
integrality constraints, variables fixing); store objective $\overline{OBJ}$
and solution $\overline{z}$. If $\overline{OBJ}<UB$, update
$UB\leftarrow\overline{OBJ}$, $\bm{z}^{*}\leftarrow\overline{\bm{z}}$.

6: – If $OBJ<UB$ and $\widehat{y}_{\ell}\notin\mathbb{Z}$ for some
$\ell=1,\cdots,q$; go to Step 2.

7: – If $OBJ<UB$ and $\bm{x}_{is}(\widehat{\bm{z}})\notin\mathcal{F}_{is}$ for
some $i,s$; go to Step 3.

8: – If $OBJ<UB$ and $\bm{x}_{is}(\widehat{\bm{z}})\in\mathcal{F}_{is},\
\forall i,s$ and $z^{p}_{i}$ is fractional for some $i,p$; go to Step 4.

9: – If $OBJ<UB$ and $\bm{x}_{is}(\widehat{\bm{z}})\in\mathcal{F}_{is},\
\forall i,s$ and $\widehat{\bm{z}}$ is integral; update $UB\leftarrow OBJ$,
$\bm{z}^{*}\leftarrow\widehat{\bm{z}}$; go to Step 5.

10: – If $OBJ\geq UB$; go to Step 5.

11:Step 3: Branching on $\bm{y}$ variables. Add two children nodes to $\Phi$
(Equation (36)); go to Step 5.

12:Step 4: Bi-partite branching. Add two children nodes to $\Phi$ (Equation
(37)); go to Step 5.

13:Step 5: Tri-partite branching. Add three children nodes to $\Phi$ (Equation
(39)); go to Step 5.

14:Step 6. If $\Phi=\emptyset$, STOP; return solution $\bm{x}^{*}$ and
objective $UB$. Otherwise, choose $cur\in\Phi$; go to Step 1.

Theorem 4.5 establishes the convergence and exactness of the algorithm, as
long as the pricing problem is solved to optimality (Step 1.2). Exactness
stems from the facts that the algorithm maintains valid lower and upper
bounds. Finiteness follows from the fact that the branching disjunction
eliminates any infeasible solution and from the finiteness of the decision
space. Note that the pricing problem only needs to be solved to optimality
upon convergence; therefore, an exact acceleration involves applying the
state-clustering dynamic programming approximation in initial iterations and
then turning to the exact algorithm. In our experiments, we merely apply the
state-clustering approximation to retain tractability, but we evaluate the
performance of the optimized solution in view of the full dynamical system—not
the clustered approximation (Section 5.3).

###### Theorem 4.5

The branch-and-price algorithm with the tri-partite branching disjunction
converges in a finite number of iterations and returns an optimal solution to
the $\mathcal{S}\mathcal{P}$ formulation.

##  5 Experimental results

We evaluate our methodology on the four problems described in Section 3. All
experiments were run on Julia v1.5.2 with the JuMP package (Dunning et al.
2017) on 10-core i9-10900k CPU, with a 3-hour limit and a 0.1% tolerance. All
instances and code are made available
online.111https://github.com/martinrame24/Prescriptive-Contagion-Analytics

###  5.1 Benefits of the state-clustering dynamic programming algorithm

Our methodology hinges on the ability to solve the pricing problem
efficiently. We first compare our state-clustering dynamic programming
algorithm (Algorithm 3) to the exact dynamic programming algorithm based on
state space enumeration (Algorithm 2). Table 1 reports computational times,
broken down into initialization times to create the state space (“init.”) and
dynamic programming times (“DP’). For all instances solved by the exact
algorithm, we report the Median Average Percentage Error (MAPE) and Median
Absolute Error (MAE) between clustered and true states.

Table 1: State-clustering dynamic programming against exhaustive enumeration
(vaccine allocation problem).

|  |  |  |  |  | Time (sec.) | MAPE (%) | MAE ($\times 10^{-4}$)
---|---|---|---|---|---|---|---|---
$n$ | $S$ | $D$ | Method | Tolerance | $|S|$ | Init. | DP | $s=1$ | $s=2$ | $s=3$ | $s=4$ | $s=1$ | $s=2$ | $s=3$ | $s=4$
51 | 4 | 6 | Enumeration | — | 79,305 | 20.01 | 1.33 | — | — | — | — | — | — | — | —
|  |  | Clustering | $\varepsilon=0.002$ | 2,350 | 3.33 | 0.08 | 0.00 | 0.00 | 0.48 | 0.69 | 0.00 | 0.00 | 0.29 | 0.27
|  |  | Clustering | $\varepsilon=0.005$ | 1,458 | 2.20 | 0.04 | 0.00 | 0.00 | 1.32 | 2.08 | 0.00 | 0.00 | 1.04 | 0.95
|  |  | Clustering | $\varepsilon=0.01$ | 982 | 1.57 | 0.03 | 0.00 | 0.01 | 7.22 | 9.38 | 0.00 | 0.01 | 3.72 | 4.02
51 | 4 | 11 | Enumeration | — | 821K | 194 | 17.1 | — | — | — | — | — | — | — |
|  |  | Clustering | $\varepsilon=0.002$ | 4,971 | 13.3 | 0.33 | 0.00 | 0.08 | 0.57 | 0.86 | 0.00 | 0.08 | 0.37 | 0.36
|  |  | Clustering | $\varepsilon=0.005$ | 2,533 | 6.74 | 0.19 | 0.00 | 0.12 | 2.86 | 3.67 | 0.00 | 0.29 | 1.59 | 1.42
|  |  | Clustering | $\varepsilon=0.01$ | 1,521 | 4.29 | 0.10 | 0.00 | 0.11 | 6.27 | 10.02 | 0.00 | 0.29 | 3.58 | 4.23
51 | 4 | 21 | Enumeration | — | 10.4M | 2,950 | 234 | — | — | — | — | — | — | — | —
|  |  | Clustering | $\varepsilon=0.002$ | 10,615 | 74.6 | 1.09 | 0.00 | 0.08 | 0.87 | 1.08 | 0.00 | 0.08 | 0.52 | 0.41
|  |  | Clustering | $\varepsilon=0.005$ | 4,477 | 26.9 | 0.58 | 0.00 | 0.09 | 3.51 | 5.12 | 0.00 | 0.01 | 1.86 | 1.57
|  |  | Clustering | $\varepsilon=0.01$ | 2,473 | 14.8 | 0.29 | 0.00 | 0.22 | 7.17 | 13.64 | 0.00 | 0.28 | 3.78 | 4.91
51 | 6 | 21 | Enumeration | — | 4.59B | n/a | n/a | — | — | — | — | — | — | — | —
|  |  | Clustering | $\varepsilon=0.002$ | 29,020 | 370 | 7.41 | n/a | n/a | n/a | n/a | n/a | n/a | n/a | n/a
|  |  | Clustering | $\varepsilon=0.005$ | 9,212 | 75.0 | 1.48 | n/a | n/a | n/a | n/a | n/a | n/a | n/a | n/a
|  |  | Clustering | $\varepsilon=0.01$ | 4,328 | 31.3 | 0.55 | n/a | n/a | n/a | n/a | n/a | n/a | n/a | n/a
51 | 8 | 21 | Enumeration | — | 2.06T | n/a | n/a | — | — | — | — | — | — | — | —
|  |  | Clustering | $\varepsilon=0.002$ | 56,635 | 1,057 | 13.6 | n/a | n/a | n/a | n/a | n/a | n/a | n/a | n/a
|  |  | Clustering | $\varepsilon=0.005$ | 15,076 | 147 | 3.50 | n/a | n/a | n/a | n/a | n/a | n/a | n/a | n/a
|  |  | Clustering | $\varepsilon=0.01$ | 6,341 | 52.4 | 1.47 | n/a | n/a | n/a | n/a | n/a | n/a | n/a | n/a
51 | 10 | 21 | Enumeration | — | 893T | n/a | n/a | — | — | — | — | — | — | — | —
|  |  | Clustering | $\varepsilon=0.002$ | 89,685 | 2,586 | 32.8 | n/a | n/a | n/a | n/a | n/a | n/a | n/a | n/a
|  |  | Clustering | $\varepsilon=0.005$ | 21,308 | 238 | 5.50 | n/a | n/a | n/a | n/a | n/a | n/a | n/a | n/a
|  |  | Clustering | $\varepsilon=0.01$ | 8,353 | 77.9 | 2.13 | n/a | n/a | n/a | n/a | n/a | n/a | n/a | n/a

  * •

“n/a” means that the exact dynamic programming algorithm does not terminate
due to memory limitations.

Note, first, that the number of states grows very quickly, reflecting the
“curse of dimensionality” in dynamic programming (Powell 2022). Even by
exploiting the segment-based decomposition in the pricing problem, the state
space scales in $\mathcal{O}(n\times D^{S})$, with trillions of states in our
largest instances. Exhaustive enumeration does not scale to even small
instances, requiring minutes to converge with 4 segments and 4 decision
epochs. The state-clustering algorithm considerably reduces the state space—by
up to a factor of $10^{10}$ to $10^{11}$. Thus, it accelerates convergence by
two orders of magnitude in small instances (seconds versus minutes), and
scales to the largest instance with 51 states, 10 epochs and 21 decisions in
minutes when the exact algorithm fails to terminate.

Moreover, the state-clustering algorithm provides a high-quality approximation
of the full state space. With $\varepsilon=0.01$, the algorithm results in a
median absolute error up to 4.9$\times 10^{-4}$ and a median absolute
percentage error up to 13.6%. A tolerance of $\varepsilon=0.002$ further
reduces the MAE within 6$\times 10^{-5}$ and the MAPE within 1.1%. As
expected, the quality of the approximation deteriorates over time due to the
propagation of errors. Whereas we cannot estimate the error over longer time
horizons because of the limitations of the enumerative algorithm, this
observation could motivate dynamic implementations of our methodology, for
instance, by re-evaluating the system’s dynamics and re-optimizing resource
allocations at each epoch. Nonetheless, these results suggest that the drift
in state approximation remains moderate, underscoring the critical role of the
state-clustering algorithm to enhance the tractability of the pricing problem
at small costs in terms of accuracy.

Finally, we provide in Appendix 10 a detailed comparison of our state-
clustering algorithm (Algorithm 3) and a $k$-means benchmark. By design, the
$k$-means algorithm leads to a smaller average distance but a larger worst-
case distance within each cluster. Per Proposition 4.4, this worst-case error
can propagate over time in dynamical systems. In turn, our algorithm results
in much higher fidelity across the time horizon, with significant benefits in
the downstream optimization—achieving a 1.4% to 25% cost reduction against the
$k$-means benchmark for the vaccine allocation problem.

###  5.2 Benefits of the branch-and-price algorithm

Armed with the state-clustering algorithm, we now solve our four prescriptive
contagion analytics problems. Table 2 reports computational times, broken down
into the state-clustering initialization, the restricted master problems, and
the pricing problems. It also reports the solution and the lower bound
obtained via column generation, via bi-partite branching (i.e., Algorithm 1
without Step 5), and with the full branch-and-price algorithm with tri-partite
branching. The column generation and bi-partite branching solutions are
obtained by solving the restricted master problem with integrality constraints
upon convergence. For brevity, the table reports results for “easy”, “medium”
and “hard” instances. Full results in Appendix 9.2 show the robustness of
these findings.

Table 2: Performance of the branch-and-price algorithm.

CPU times (s) Solution quality Problem Instance Method Nodes Init. RMP PP
Total Upper bound Solution Gap Vaccine Easy Column Generation 1 21 0.02 0.22
0.29 93,695 92,261 1.53% allocation Bi-partite B&P 28 21 0.27 1.43 7.11 93,379
93,266 0.12% Tri-partite B&P 46 21 0.42 2.2 9.6 93,353 93,263 0.09% Medium
Column Generation 1 43 0.022 0.81 0.92 156,465 155,624 0.54% Bi-partite B&P
4,778 43 53 426 1,064 156,183 156,026 0.1% Tri-partite B&P 4,778 43 55 420
1,058 156,183 156,027 0.1% Hard Column Generation 1 144 0.03 4.07 4.32 225,323
222,592 1.21% Bi-partite B&P 9,990 144 101 1,765 3,475 225,048 224,495 0.24%
Tri-partite B&P 9,990 144 103 1,808 3,535 225,048 224,495 0.24% Vaccination
Easy, $K=2$ Column Generation 1 36 0.05 0.03 0.46 786 652 29.15% centers Bi-
partite B&P 2,919 36 77 11 608 764 763 0.09% Tri-partite B&P 2,919 36 75 11
604 764 763 0.09% Easy, $K=3$ Column Generation 1 43 0.02 0.03 0.57 678 652
3.75% Bi-partite B&P 56,273 43 1,272 142 10,809 675 674 0.23% Tri-partite B&P
58,496 43 1,382 137 10,801 675 674 0.23% Hard, $K=2$ Column Generation 1 45
0.19 0.43 1.11 8200 6660 29.91% Bi-partite B&P 11,203 45 216 147 10,808 8,134
7,912 2.73% Tri-partite B&P 11,899 45 224 178 10,808 8,132 7,917 2.64% Hard,
$K=3$ Column Generation 1 47 0.27 0.31 1.11 8,200 6,822 19.07% Bi-partite B&P
15,982 47 356 207 10,885 8,048 7,714 4.14% Tri-partite B&P 16,358 47 363 209
10,955 8,048 7,716 4.13% Content Easy Column Generation 1 79 0.04 4.6 4.8
1.013 0.98 3.86% promotion Bi-Partite B&P 696 79 9.3 198 327 1.013 1.009 0.53%
Tri-partite B&P 142 79 1.6 52 74 1.013 1.0125 0.09% Medium Column Generation 1
58 0.04 5.01 5.4 1.52 1.44 5.66% Bi-Partite B&P 17,960 58 294 1,949 7,201 1.52
1.51 0.83% Tri-partite B&P 17,819 58 3,380 2,0570 7,201 1.52 1.517 0.31% Hard
Column Generation 1 134 0.08 15 16 2.02 1.87 7.08% Bi-Partite B&P 12,606 134
229 3,221 10,801 2.009 1.99 1.04% Tri-partite B&P 11,864 134 229 3,279 10,801
2.009 1.99 1.04% Congestion Easy Column Generation 1 10 0.012 0.08 0.11 0.82
0.81 1.23% mitigation Bi-partite B&P 2 10 0.018 0.126 0.56 0.82 0.81 0.16%
Tri-partite B&P 8 10 0.053 0.23 0.8 0.82 0.82 0% Medium Column Generation 1
185 0.024 6.2 6.27 2.20 2.05 6.61% Bi-partite B&P 48 185 0.34 32 34 2.19 2.19
0.21% Tri-partite B&P 57 185 0.43 39 42 2.19 2.19 0.06% Hard Column Generation
1 459 0.044 20 20 3.32 2.93 11.9% Bi-partite B&P 206 459 4.14 350 389 3.31
3.31 0.2% Tri-partite B&P 243 459 4.9 396 443 3.31 3.31 0.09%

  * •

Vaccine allocation: $n=51$; $B=2.5$M; $S=8$, $D=6$ (“easy”); $S=10$, $D=11$
(“medium”); $S=12$, $D=21$ (“hard”).

  * •

Vaccination centers: Group H (“easy”) and Group A (“hard”), with $K=2$ or
$K=3$ facilities per group, $S=6$.

  * •

Content promotion: $n=20$; $K=2$ ; $S=6$, $D=21$, (“easy”); $S=8$, $D=11$,
(“medium”); $S=10$, $D=21$ (“hard”).

  * •

Congestion mitigation: $n=5$; $B^{1}=6$ and $B^{2}=4$ ($D=35$); $S=2$
(“easy”); $S=4$ (“medium”); $S=6$ (“hard”).

Table 3: Detailed results for the vaccination centers problem (full United States, 3 facilities per group) |  |  |  |  |  |  | CPU times (s) | Solution quality
---|---|---|---|---|---|---|---|---
Group | n | S | D | F | Method | Nodes | Init. | RMP | PP | Total | Upper bound | Solution | Gap
A | 7 | 6 | 12 | 3 | Column Generation | 1 | 47 | 0.27 | 0.31 | 1.11 | 8,200 | 6,822 | 19.07%
|  |  |  |  | Tri-partite B&P | 16,358 | 47 | 363 | 209 | 10,955 | 8,048 | 7,716 | 4.13%
B | 7 | 6 | 10 | 3 | Column Generation | 1 | 43 | 0.06 | 0.34 | 0.88 | 3,659 | 3,322 | 9.20%
|  |  |  |  | Tri-partite B&P | 11,965 | 43 | 329 | 134 | 10,806 | 3,569 | 3,467 | 2.87%
C | 4 | 6 | 13 | 3 | Column Generation | 1 | 38 | 0.04 | 0.11 | 0.55 | 3,880 | 3,171 | 18.28%
|  |  |  |  | Tri-partite B&P | 10,101 | 38 | 369 | 98 | 10,812 | 3,880 | 3,692 | 4.84%
D | 4 | 6 | 5 | 3 | Column Generation | 1 | 38 | 0.04 | 0.04 | 0.48 | 1,056 | 733 | 30.65%
|  |  |  |  | Tri-partite B&P | 22,600 | 38 | 539 | 45 | 10,824 | 1,052 | 1,000 | 4.95%
E | 6 | 6 | 16 | 3 | Column Generation | 1 | 55 | 0.08 | 0.51 | 1.02 | 8,331 | 6,927 | 16.86%
|  |  |  |  | Tri-partite B&P | 7,824 | 55 | 214 | 199 | 10,832 | 8,316 | 8,053 | 3.16%
F | 5 | 6 | 13 | 3 | Column Generation | 1 | 47 | 0.07 | 0.25 | 0.83 | 3,593 | 2,466 | 45.25%
|  |  |  |  | Tri-partite B&P | 16,882 | 47 | 698 | 181 | 10,822 | 3,530 | 3,328 | 5.74%
G | 4 | 6 | 4 | 3 | Column Generation | 1 | 42 | 0.02 | 0.02 | 0.50 | 1,985 | 1,679 | 15.39%
|  |  |  |  | Tri-partite B&P | 39,651 | 42 | 609 | 72 | 10,802 | 1,955 | 1,921 | 1.76%
H | 6 | 6 | 4 | 3 | Column Generation | 1 | 43 | 0.02 | 0.03 | 0.57 | 678 | 652 | 3.75%
|  |  |  |  | Tri-partite B&P | 58,496 | 43 | 1,382 | 137 | 10,801 | 675 | 674 | 0.23%
I | 4 | 6 | 20 | 3 | Column Generation | 1 | 39 | 0.07 | 0.23 | 0.68 | 3,573 | 2,587 | 34.89%
|  |  |  |  | Tri-partite B&P | 55,864 | 39 | 4,357 | 507 | 10,871 | 3,423 | 3,399 | 0.72%
J | 4 | 6 | 5 | 3 | Column Generation | 1 | 36 | 0.03 | 0.03 | 0.47 | 626 | 565 | 9.73%
|  |  |  |  | Tri-partite B&P | 25,016 | 36 | 327 | 41 | 10,803 | 624 | 565 | 9.46%
Full | 51 | 6 | — | 3 | Column Generation | — | — | — | — | — | 35,581 | 28,924 | 18.7%
|  |  |  |  | Tri-partite B&P | — | — | — | — | — | 35,072 | 33,815 | 3.58%

Note, first, that the column generation algorithm is highly scalable for the
three resource allocation problems: vaccine allocation, content promotion, and
congestion mitigation. Indeed, column generation solves the set partitioning
relaxation in seconds even in the harder instances. The restricted master
problem is virtually instantaneous but the pricing problem is more time-
consuming due to the non-linear system dynamics—reinforcing the need for an
efficient dynamic programming algorithm (Section 5.1). In fact, the state-
clustering algorithm shifts the complexity away from the dynamic programming
algorithm itself—enabling the online column generation algorithm to terminate
in seconds—toward the offline state space clustering—which can take up to
minutes.

Moreover, the column generation algorithm can derive high-quality solutions,
thanks to the tight set partitioning formulation for these problems. In the
vaccine allocation problem in particular, the column generation algorithm
yields solutions within a 1-2% optimality gap. For the content promotion and
congestion mitigation problems, column generation generally performs well but
can also leave a larger optimality gap—due to the sparsity constraint in
content promotion and to the coordination of prevention and treatment
resources in congestion mitigation. This limitation motivates the branch-and-
price algorithm to further improve solutions and tighten optimality gaps.

Turning to our main observation, the branch-and-price algorithm solves every
instance of the vaccine allocation, content promotion and congestion
mitigation problems to near-optimality. In fact, it reaches optimality within
a 0.1% tolerance in most instances, and leaves a small optimality gap
otherwise. The number of nodes generally remains moderate but can grow larger
for harder instances, leading to higher computational requirements. As a
result, the branch-and-price algorithm terminates in minutes to hours, but can
still lead to strong solution improvements from column generation. These
benefits can go up to 6-12% in the content promotion and congestion mitigation
cases. Even in vaccine allocation, where column generation leaves a small gap,
the branch-and-price solution yields a 1% improvement, amounting to an extra
2,000 lives saved over three months.

Finally, note that the standard bi-partite branching scheme is not sufficient
to guarantee convergence to an optimal solution. Notably, it leaves a 1%
optimality gap in content promotion instances and a 0.2% gap in congestion
mitigation instances. This result underscores the impact of system non-
linearities on the branch-and-price algorithm. Instead, the tri-partite
branching scheme developed in this paper can be instrumental to ensure
convergence to an optimal solution, while retaining an efficient branching
tree and an efficient pricing problem structure.

Next, recall that the vaccination centers problem features a joint facility
location and resource allocation structure, leading to a looser set-
partitioning relaxation. As a result, the column generation heuristic leaves
large optimality gaps up to 30%. Restoring integrality via branch-and-price is
much more challenging. Interestingly, the restricted master problem ends up
being more time-consuming than the pricing problem, thus highlighting the
joint complexities from coupled resource allocation and from non-linear system
dynamics. Yet, the branch-and-price algorithm leads to significant solution
improvements and gap reductions: it returns an optimal solution in the “easy”
two-facility instance, a near-optimal solution in the “hard” two-facility
instances, and moderate optimality gaps in the three-facility cases (0-4%
versus 3-30% with column generation).

To shed more light into these findings, Table 10 reports results for the all
ten groups constituting the full United States. Despite variations in the
problem’s complexity, our core observations are highly robust: the column
generation algorithm leaves consistently large optimality gaps, and the
branch-and-price algorithm significantly improves the solution and tightens
the gap. Ultimately, the branch-and-price algorithm results in an estimated
5,000 extra lives saved across the country over a 6-week period—a 17%
improvement over the column generation solution.

###  5.3 Practical impact of the methodology

Tables 4 and 5 evaluate the practical benefits of our branch-and-price
methodology, as compared to a do-nothing baseline, practical benchmarks (i.e.,
easily-implementable heuristics) and optimization benchmarks (i.e., existing
methods from the literature). To provide a fair assessment, we evaluate all
solutions with the full continuous-state contagion models, as opposed to the
state-clustering approximation used in our algorithm—thus ensuring an apples-
to-apples comparison.

We define two practical benchmarks for the vaccine allocation and the
congestion mitigation problems, as well as slight variants for the content
promotion problem due to the sparsity constraint.

  * –

uniform allocation: each segment receives the same amount of resources of at
each epoch (in content promotion, we randomize the $K$ products at each epoch,
repeated 100 times).

  * –

cost-based allocation: a domain-based benchmark in which each segment receives
a constant share of resources proportionally to the total cost under the do-
nothing baseline (in content promotion, we select the $K$ products for which
promotions have the strongest impact).

For the vaccination centers problem, we consider a demographic-based “top-K”
benchmark that selects the $K$ facilities that can serve the most people
within access distance restrictions. We test it both with uniform vaccine
allocations across the resulting $K$ facilities, and with optimized vaccine
allocations (obtained by fixing the facility variables and optimizing
subsequent resource allocation).

We also define three optimization benchmarks for the vaccine allocation and
vaccination centers problems, from the recent prescriptive contagion analytics
literature in epidemiology applications:

  * –

MIQO implementation: mixed-integer bilinear optimization implementation of
Problem $(\mathcal{P})$ in Gurobi 9.5 (Gurobi 2019), based on a time-
discretization approximation of the ODEs.

  * –

discretization to approximate Problem $(\mathcal{P})$ via mixed-integer linear
optimization, using time discretization to eliminate continuous-time dynamics
and a staircase approximation of the infected population to handle
bilinearities. This benchmark mirrors the approach from Fu et al. (2021) in a
robust optimization setting. To optimize its performance, we divide the
$[0,2\%]$ interval into sub-intervals of length $\delta$ (since the infected
population never exceeded 2%). We consider a coarse and a granular
discretization, with $\delta=0.001$ and $\delta=0.002$.

  * –

coordinate descent heuristic, which circumvents bilinearities by iterating
between optimizing vaccine allocations over discretized time increments with
fixed numbers of infections (via mixed-integer linear optimization) and re-
estimating infections (Bertsimas et al. 2022).

Table 4: Death toll comparison for the vaccine allocation problem (full
country, $D=21$).

|  | $S=6$ | $S=8$ | $S=12$
---|---|---|---|---
Budget | Method | Time (sec.) | Deaths | Time (sec.) | Deaths | Time (sec.) | Deaths
— | Do nothing | — | 573.37K | — | 604.22K | — | 649.75K
2.5M | Uniform allocation | — | -5.91K | — | -11.02K | — | -21.42K
| Cost-based allocation | — | -6.59K | — | -12.92K | — | -27.64K
| MIQO implementation | n/a | n/a | n/a | n/a | n/a | n/a
|  Discretization ($\delta=0.002$) |  1,000∗ | -5.03K | n/a | n/a | n/a | n/a
|  Discretization ($\delta=0.001$) |  1,000∗ | -5.56K | n/a | n/a | n/a | n/a
| Coordinate Descent | 5.43 | -11.36K | 4.73 | -20.02K | 6.03 | -39.01K
|  Branch-and-price ($\varepsilon=0.002$) | 426.2 | -11.28K | 1234.2 | -20.43K | 3671.4 | -39.41K
7M | Uniform allocation | — | -13.22K | — | -23.67K | — | -43.60K
| Cost-based allocation | — | -17.44K | — | -31.77K | — | -58.71K
| MIQO implementation | n/a | n/a | n/a | n/a | n/a | n/a
|  Discretization ($\delta=0.002$) |  1,000∗ | -5.88K | n/a | n/a | n/a | n/a
|  Discretization ($\delta=0.001$) |  1,000∗ | -2.88K | n/a | n/a | n/a | n/a
| Coordinate Descent | 5.46 | -20.61K | 12.7 | -36.22K | 118.5 | -62.55K
|  Branch-and-price ($\varepsilon=0.002$) | 67.5 | -21.46K | 156.8 | -37.38K | 708.7 | -65.65K

  * •

$*$ and “n/a”: no optimal and feasible solution, respectively. Bold font:
solutions within $1\%$ of the best-found solution.

Table 5: Performance comparison: vaccination centers, congestion mitigation
and content promotion problems.

Vaccination Centers Content Promotion Congestion Mitigation Method Time (sec.)
Deaths Method Time (sec.) Market share Method Time (sec.) Cost Do nothing —
573K Do nothing — 2.1 p.p. Do nothing — 18.98 Top-K, Uniform — -6.98K Uniform
(random) — +5.1 p.p. Uniform — -1.80% Top-K, Optimized — -7K Cost-based — +7.3
p.p. Cost-based — -1.79% MIQO n/a n/a — — — — — — Discretization 10,275∗ -4.3K
— — — — — — Coordinate descent 254 -7.7K — — — — — — Branch-and-price 10,978
-8.6K Branch-and-price 11,864 +10.0 p.p. Branch-and-price 416 -3.32%

  * •

$*$ and “n/a”: no optimal and feasible solution, respectively. Bold font:
solutions within $1\%$ of the best-found solution.

  * •

Vaccination centers: $S=6$, $K=3$, results for the full country across 10
groups.

  * •

Content promotion: $n=20$ $K=2$, $S=10$, $D=21$ (“hard” instance); market
share averaged across 20 products.

  * •

Congestion mitigation: $n=5$; $B^{1}=6$; $B^{2}=4$; $S=6$ (“hard” instance).

Results show that our solution provides significant benefits against all
benchmarks. Table 4 reports results for the vaccine allocation problem, with
weekly budgets of 2.5 million and 7 million vaccines—corresponding to the
supplies originally planned and actually available in the United States in
2021, respectively. Without vaccinations, the pandemic would lead to an
estimated 650,000 fatalities over three months (accounting for undetected
deaths). Under a uniform allocation strategy, the vaccination campaign can
save around 21,000 of these fatalities (or 3.3%) with the smaller vaccine
budget, and 43,000 fatalities (or 6.6%) with the larger budget. The cost-based
benchmark saves up to 30% of additional lives, by capturing epidemiological
information. But then, our solution saves 20-70% additional lives over six
weeks and 12-40% additional lives over three months, as compared to the cost-
based benchmark. In other words, optimized stockpile management can increase
the effectiveness of the vaccination campaign by a factor of 1.12 to 1.7
without increasing vaccine capacity or vaccine efficacy. In absolute terms,
these improvements represent 7,000 to 12,000 extra lives saved over three
months, demonstrating the importance of vaccine distribution and the edge of
optimization to guide resource allocation in complex contagion systems.

We obtain similar findings for the other three problems (Table 5). In the
content promotion and congestion mitigation problems, the relative gains of
the optimized solution are even more significant, estimated at 84-96% against
uniform allocation and at 37-85% against cost-based allocation. Simple
benchmarks perform comparatively worse due to the sparsity constraint in the
content promotion problem and the interdependencies between prevention and
treatment vehicles in the congestion mitigation problem—leading to stronger
benefits of optimization. In the vaccination centers problem, the optimized
solution also provides significant improvements against the simple
demographic-based top-K benchmark. Specifically, the benefit of optimization
amounts to 23-31% under uniform vaccine allocations and to 3-23% under
optimized vaccine allocations. Exactly like epidemiological proxies are
insufficient to guide tactical resource allocation, simple demographic proxies
are insufficient to guide strategic planning in dynamical and non-linear
contagion systems.

Finally, our branch-and-price methodology is instrumental to reap the benefits
of optimization, with significant impact against simpler optimization
benchmarks. First, Tables 4 and 5 show the strong limitations of simple
workarounds to circumvent the non-convex contagion dynamics in the
optimization problem. Direct implementation using mixed-integer quadratic
optimization solvers fails to even provide a feasible solution in small six-
week instances. Similarly, mixed-integer linear optimization implementation
(using a discretization-based staircase approximation of the number of
infections) performs worse than even uniform allocation in small instances and
does not return a feasible solution in larger instances. In comparison, our
algorithm achieves a Pareto improvement: faster computational times, stronger
scalability, and much higher-quality solutions. Our methodology can also
outperform the tailored coordinate descent heuristic, albeit in longer
computational times. For the vaccine allocation problem, our solution yields a
1-5% improvement, resulting in up to 3,000 extra lives saved over a three-
month horizon. In fact, the benefits of our methodology increases as we
increase the number of available vaccines, because the interaction between the
vaccine allocation decisions and the infection dynamics become stronger,
favoring our global optimization approach. In the vaccination centers problem,
the benefits can also be significant when coordinate descent leads to local
optima that induce a different set of facilities. In the $F=3$ instance, our
methodology results in an extra 900 lives saved over six weeks—a 12%
improvement.

We conclude by illustrating vaccine allocations in Figure 4 and vehicle
allocations in Figure 5. In vaccine allocation, the states where the pandemic
is more costly (e.g., those with a larger population, more infections)
generally receive more vaccines, but the relationship is not monotonic.
Moreover, vaccine allocations exhibit different temporal patterns; for
instance, New York and Illinois receive most vaccines early on; California and
Texas receive most of their vaccines later on; and Florida stands in between.
In the congestion mitigation problem, treatment vehicles are mainly sent to
dense urban centers and near the airport, to clear large accidents and ease
recovery; in contrast, prevention vehicles are primarily sent to suburban
regions to prevent the spread of congestion. In both cases, there is no one-
size-fits-all allocation strategy. Instead, the benefits of optimization stem
from fine-tuning resource allocations to mitigate near-term impact via
treatment interventions and to manage network propagation via prevention
interventions.

Figure 4: Spatial-temporal vaccine allocation over a 12 weeks. States are
ordered from top to bottom in decreasing order of the total cost under do-
nothing—a proxy for the prevalence of the pandemic.

(a) Treatment vehicles, 8AM

(b) Treatment vehicles, 9AM

(c) Treatment vehicles, 10AM

(d) Prevention vehicles, 8AM

(e) Prevention vehicles, 8AM

(f) Prevention vehicles, 8AM

Figure 5: Allocation of treatment vehicles (top) and prevention vehicles
(bottom) in Singapore over three hours.

In summary, our branch-and-price methodology generates consistently high-
quality solutions in large-scale practical instances of our prescriptive
contagion analytics problems, with significant benefits as compared to all
practical and optimization benchmarks. In Appendix 9.3, we report additional
results showing the robustness of these findings to parameter estimation
errors.

##  6 Conclusion

Predictive contagion systems have shown considerable success in epidemiology
and other domains of science, engineering and management. This paper developed
a prescriptive algorithm to optimize spatial-temporal resource allocation
decisions in systems governed by contagion dynamics—and in more general
dynamical systems. By combining the difficulties of combinatorial optimization
and those of dynamical systems, however, this class of problems exhibits a
large-scale and complex mixed-integer non-convex optimization structure with
continuous-time ODE constraints. In response, this paper has developed a
branch-and-price algorithm, using (i) a set partitioning reformulation to
eliminate non-linearities; (ii) column generation to separate combinatorial
optimization in a master problem from non-linear dynamics in a pricing
problem; (iii) a novel state-clustering algorithm for discrete-decision
continuous-state dynamic programming to solve the pricing problem; and (iv) a
novel tri-partite branching scheme on natural variables to circumvent the non-
linearities.

We implemented the methodology on four prescriptive analytics contagion
problems: vaccine allocation, mass vaccination centers deployment, content
promotion, and congestion mitigation. The branch-and-price algorithm scales to
realistic and otherwise-intractable instances, with up to 21 decisions, 51
regions and 12 decision epochs in vaccine allocation, hence
$\mathcal{O}(21^{612})$ possible decisions overall. From a practical
standpoint, the methodology significantly outperforms easily-implementable
benchmarks based on demographic or epidemiological proxies, as well as state-
of-the-art optimization algorithms. In the vaccine allocation example, the
methodology can improve the effectiveness of a vaccination campaign by 12-70%,
resulting in 7,000 to 12,000 extra lives saves in a situation mirroring the
midst of the COVID-19 pandemic in the United States. These results are robust
across problem instances and robust to parameter estimation errors.
Ultimately, our prescriptive contagion analytics methodology can deliver
significant benefits across contagion-based domains, by fine-tuning resource
allocations based on spatial-temporal system dynamics.

The promising results reported in this paper also motivate extensions of our
methodology to tackle prescriptive contagion analytics problems with different
decision-making structures (e.g., non-discrete resource allocation decisions,
dynamic resource allocations) and more complex system dynamics (e.g., inter-
population mixing). Perhaps the main limitation of this paper is the focus on
deterministic resource allocation problems; although our robustness tests in
Appendix 9.3 showed the benefits of our solutions under model
misspecification, our optimization methodology could be combined in future
work with the robust epidemiological approach from Fu et al. (2021). Yet, this
paper provides methodological foundations to support complex resource
allocation decisions in a broad class of contagion-based systems, and of
dynamical systems more generally.

## References

  * Acemoglu et al. (2021) Acemoglu D, Chernozhukov V, Werning I, Whinston MD (2021) Optimal targeted lockdowns in a multigroup sir model. _American Economic Review: Insights_ 3(4):487–502.
  * Adam (2020) Adam D (2020) Special report: The simulations driving the world’s response to covid-19. _Nature_ 580(7802):316–319.
  * Agrawal et al. (2021) Agrawal S, Yin S, Zeevi A (2021) Dynamic pricing and learning under the bass model. _Proceedings of the 22nd ACM Conference on Economics and Computation_ , 2–3.
  * Alban et al. (2022) Alban A, Blaettchen P, de Vries H, Van Wassenhove LN (2022) Resource allocation with sigmoidal demands: Mobile healthcare units and service adoption. _Manufacturing & Service Operations Management_ 24(6):2944–2961.
  * Allman and Zhang (2021) Allman A, Zhang Q (2021) Branch-and-price for a class of nonconvex mixed-integer nonlinear programs. _Journal of Global Optimization_ 81(4):861–880.
  * Alvarez et al. (2021) Alvarez F, Argente D, Lippi F (2021) A simple planning problem for covid-19 lock-down, testing, and tracing. _American Economic Review: Insights_ 3(3):367–82.
  * Andreas et al. (2008) Andreas AK, Smith JC, Küçükyavuz S (2008) Branch-and-price-and-cut algorithms for solving the reliable h-paths problem. _Journal of Global Optimization_ 42(4):443–466.
  * Androulakis et al. (1995) Androulakis IP, Maranas CD, Floudas CA (1995) $\alpha$bb: A global optimization method for general constrained nonconvex problems. _Journal of Global Optimization_ 7:337–363.
  * Anstreicher (2009) Anstreicher KM (2009) Semidefinite programming versus the reformulation-linearization technique for nonconvex quadratically constrained quadratic programming. _Journal of Global Optimization_ 43:471–484.
  * Anstreicher and Burer (2021) Anstreicher KM, Burer S (2021) Quadratic optimization with switching variables: the convex hull for n= 2. _Mathematical Programming_ 188(2):421–441.
  * Balderrama et al. (2022) Balderrama R, Peressutti J, Pinasco JP, Vazquez F, Vega CSdl (2022) Optimal control for a sir epidemic model with limited quarantine. _Scientific Reports_ 12(1):1–26.
  * Barnhart et al. (1998) Barnhart C, Johnson EL, Nemhauser GL, Savelsbergh MW, Vance PH (1998) Branch-and-price: Column generation for solving huge integer programs. _Operations research_ 46(3):316–329.
  * Bass (1969) Bass FM (1969) A new product growth for model consumer durables. _Management science_ 15(5):215–227.
  * Behrens et al. (2000) Behrens DA, Caulkins JP, Tragler G, Feichtinger G (2000) Optimal control of drug epidemics: prevent and treat—but not at the same time? _Management Science_ 46(3):333–347.
  * Bennouna et al. (2022) Bennouna A, Joseph J, Nze-Ndong D, Perakis G, Singhvi D, Lami OS, Spantidakis Y, Thayaparan L, Tsiourvas A (2022) Covid-19: Prediction, prevalence, and the operations of vaccine allocation. _Manufacturing & Service Operations Management_ .
  * Bennouna et al. (2021) Bennouna MA, Pachamanova D, Perakis G, Skali Lami O (2021) Learning the minimal representation of a dynamic system from transition data. _Available at SSRN 3785547_ .
  * Berge et al. (2017) Berge T, Lubuma JS, Moremedi G, Morris N, Kondera-Shava R (2017) A simple mathematical model for ebola in africa. _Journal of biological dynamics_ 11(1):42–74.
  * Bertsekas (1975) Bertsekas D (1975) Convergence of discretization procedures in dynamic programming. _IEEE Transactions on Automatic Control_ 20(3):415–419.
  * Bertsekas (2015) Bertsekas DP (2015) Dynamic programming and optimal control 4th edition, volume ii. _Athena Scientific_ .
  * Bertsekas et al. (1988) Bertsekas DP, Castanon DA, et al. (1988) Adaptive aggregation methods for infinite horizon dynamic programming .
  * Bertsimas et al. (2021) Bertsimas D, Boussioux L, Cory-Wright R, Delarue A, Digalakis V, Jacquillat A, Kitane DL, Lukin G, Li M, Mingardi L, et al. (2021) From predictions to prescriptions: A data-driven response to covid-19. _Health care management science_ 24(2):253–272.
  * Bertsimas et al. (2022) Bertsimas D, Digalakis Jr V, Jacquillat A, Li ML, Previero A (2022) Where to locate covid-19 mass vaccination facilities? _Naval Research Logistics (NRL)_ 69(2):179–200.
  * Bestehorn et al. (2021) Bestehorn F, Hansknecht C, Kirches C, Manns P (2021) Mixed-integer optimal control problems with switching costs: a shortest path approach. _Mathematical Programming_ 188(2):621–652.
  * Borraz-Sánchez and Haugland (2011) Borraz-Sánchez C, Haugland D (2011) Minimizing fuel cost in gas transmission networks by dynamic programming and adaptive discretization. _Computers & Industrial Engineering_ 61(2):364–372.
  * Burer and Letchford (2012) Burer S, Letchford AN (2012) Non-convex mixed-integer nonlinear programming: A survey. _Surveys in Operations Research and Management Science_ 17(2):97–106.
  * Casagrandi et al. (2006) Casagrandi R, Bolzoni L, Levin SA, Andreasen V (2006) The sirc model and influenza a. _Mathematical biosciences_ 200(2):152–169.
  * Caulkins et al. (2020) Caulkins J, Grass D, Feichtinger G, Hartl R, Kort PM, Prskawetz A, Seidl A, Wrzaczek S (2020) How long should the covid-19 lockdown continue? _Plos one_ 15(12):e0243413.
  * Chance et al. (2008) Chance DM, Hillebrand E, Hilliard JE (2008) Pricing an option on revenue from an innovation: An application to movie box office revenue. _Management Science_ 54(5):1015–1028.
  * Cosguner and Seetharaman (2022) Cosguner K, Seetharaman P (2022) Dynamic pricing for new products using a utility-based generalization of the bass diffusion model. _Management Science_ 68(3):1904–1922.
  * Dehning et al. (2020) Dehning J, Zierenberg J, Spitzner FP, Wibral M, Neto JP, Wilczek M, Priesemann V (2020) Inferring change points in the spread of covid-19 reveals the effectiveness of interventions. _Science_ 369(6500):eabb9789.
  * Dunning et al. (2017) Dunning I, Huchette J, Lubin M (2017) JuMP: A modeling language for mathematical optimization. _SIAM review_ 59(2):295–320.
  * Fu et al. (2021) Fu C, Sim M, Zhou M (2021) Robust epidemiological analytics. _Available at SSRN 3869521_ .
  * Fujie and Kojima (1997) Fujie T, Kojima M (1997) Semidefinite programming relaxation for nonconvex quadratic programs. _Journal of Global optimization_ 10:367–380.
  * Goldman and Lightwood (2002) Goldman S, Lightwood J (2002) Cost optimization in the sis model of infectious disease with treatment. _Topics in Economic Analysis & Policy_ 2(1):1007.
  * Göttlich et al. (2021) Göttlich S, Hante FM, Potschka A, Schewe L (2021) Penalty alternating direction methods for mixed-integer optimal control with combinatorial constraints. _Mathematical Programming_ 188(2):599–619.
  * Grüne and Semmler (2004) Grüne L, Semmler W (2004) Using dynamic programming with adaptive grid scheme for optimal control problems in economics. _Journal of Economic Dynamics and Control_ 28(12):2427–2456.
  * Gu et al. (2020) Gu X, Ahmed S, Dey SS (2020) Exact augmented lagrangian duality for mixed integer quadratic programming. _SIAM Journal on Optimization_ 30(1):781–797.
  * Günlük and Linderoth (2010) Günlük O, Linderoth J (2010) Perspective reformulations of mixed integer nonlinear programs with indicator variables. _Mathematical programming_ 124:183–205.
  * Gurobi (2019) Gurobi (2019) Non-convex quadratic optimization. https://www.gurobi.com/wp-content/uploads/2020-01-14$\\_$Non-Convex-Quadratic-Optimization-in-Gurobi-9.0-Webinar.pdf?x41151.
  * Gurobi Optimization (2020) Gurobi Optimization (2020) Gurobi optimizer reference manual. URL http://www.gurobi.com.
  * Hante and Sager (2013) Hante FM, Sager S (2013) Relaxation methods for mixed-integer optimal control of partial differential equations. _Computational Optimization and Applications_ 55:197–225.
  * Ho et al. (2002) Ho TH, Savin S, Terwiesch C (2002) Managing demand and sales dynamics in new product diffusion under supply constraint. _Management science_ 48(2):187–206.
  * Hsiang et al. (2020) Hsiang S, Allen D, Annan-Phan S, Bell K, Bolliger I, Chong T, Druckenmiller H, Huang LY, Hultgren A, Krasovich E, et al. (2020) The effect of large-scale anti-contagion policies on the covid-19 pandemic. _Nature_ 584(7820):262–267.
  * Jung et al. (2015) Jung MN, Reinelt G, Sager S (2015) The lagrangian relaxation for the combinatorial integral approximation problem. _Optimization Methods and Software_ 30(1):54–80.
  * Kermack and McKendrick (1927) Kermack WO, McKendrick AG (1927) A contribution to the mathematical theory of epidemics. _Proceedings of the royal society of london. Series A, Containing papers of a mathematical and physical character_ 115(772):700–721.
  * Krishnan and Jain (2006) Krishnan TV, Jain DC (2006) Optimal dynamic advertising policy for new products. _Management Science_ 52(12):1957–1969.
  * Lee and Grossmann (2001) Lee S, Grossmann IE (2001) A global optimization algorithm for nonconvex generalized disjunctive programming and applications to process systems. _Computers & Chemical Engineering_ 25(11-12):1675–1697.
  * Li (2020) Li H (2020) Optimal pricing under diffusion-choice models. _Operations Research_ 68(1):115–133.
  * Li et al. (2013) Li H, Armbruster D, Kempf KG (2013) A population-growth model for multiple generations of technology products. _Manufacturing & Service Operations Management_ 15(3):343–360.
  * Li et al. (2022) Li ML, Bouardi HT, Lami OS, Trikalinos TA, Trichakis N, Bertsimas D (2022) Forecasting covid-19 and analyzing the effect of government interventions. _Operations Research_ .
  * Lin et al. (2021) Lin Y, Wang M, Shen ZJM, Zhang H, Zhang RP (2021) Content promotion for online content platforms with network diffusion effect. _Available at SSRN 3863104_ .
  * Lolos et al. (2017) Lolos K, Konstantinou I, Kantere V, Koziris N (2017) Adaptive state space partitioning of markov decision processes for elastic resource management. _2017 IEEE 33rd International Conference on Data Engineering (ICDE)_ , 191–194 (IEEE).
  * Long et al. (2018) Long EF, Nohdurft E, Spinler S (2018) Spatial resource allocation for emerging epidemics: A comparison of greedy, myopic, and dynamic policies. _Manufacturing & Service Operations Management_ 20(2):181–198.
  * Mai and Hu (2022) Mai Y, Hu B (2022) Optimizing free-to-play multiplayer games with premium subscription. _Management Science_ .
  * Manns and Kirches (2020) Manns P, Kirches C (2020) Improved regularity assumptions for partial outer convexification of mixed-integer pde-constrained optimization problems. _ESAIM: Control, Optimisation and Calculus of Variations_ 26:32.
  * Mehrotra et al. (2020) Mehrotra S, Rahimian H, Barah M, Luo F, Schantz K (2020) A model of supply-chain decisions for resource sharing with an application to ventilator allocation to combat covid-19. _Naval Research Logistics (NRL)_ 67(5):303–320.
  * Moderna (2020) Moderna (2020) Storage & handling: Moderna covid-19 vaccine (eua). https://eua.modernatx.com/covid19vaccine-eua/providers/storage-handling.
  * Ndeffo Mbah and Gilligan (2011) Ndeffo Mbah ML, Gilligan CA (2011) Resource allocation for epidemic control in metapopulations. _PLoS one_ 6(9):e24577.
  * Nowak et al. (2018) Nowak I, Breitfeld N, Hendrix EM, Njacheun-Njanzoua G (2018) Decomposition-based inner-and outer-refinement algorithms for global optimization. _Journal of Global Optimization_ 72(2):305–321.
  * Powell (2022) Powell WB (2022) _Reinforcement Learning and Stochastic Optimization: A Unified Framework for Sequential Decisions_ (John Wiley & Sons).
  * Pyeatt et al. (2001) Pyeatt LD, Howe AE, et al. (2001) Decision tree function approximation in reinforcement learning. _Proceedings of the third international symposium on adaptive systems: evolutionary computation and probabilistic graphical models_ , volume 2, 70–77 (Cuba).
  * Rahmandad et al. (2018) Rahmandad H, Henderson R, Repenning NP (2018) Making the numbers?“short termism” and the puzzle of only occasional disaster. _Management Science_ 64(3):1328–1347.
  * Rahmandad and Ton (2020) Rahmandad H, Ton Z (2020) If higher pay is profitable, why is it so rare? modeling competing strategies in mass market services. _Organization Science_ 31(5):1053–1071.
  * Ray et al. (2020) Ray EL, Wattanachit N, Niemi J, Kanji AH, House K, Cramer EY, Bracher J, Zheng A, Yamana TK, Xiong X, et al. (2020) Ensemble forecasts of coronavirus disease 2019 (covid-19) in the us. _MedRXiv_ .
  * Rowthorn et al. (2009) Rowthorn RE, Laxminarayan R, Gilligan CA (2009) Optimal control of epidemics in metapopulations. _Journal of the Royal Society Interface_ 6(41):1135–1144.
  * Ryoo and Sahinidis (1996) Ryoo HS, Sahinidis NV (1996) A branch-and-reduce approach to global optimization. _Journal of global optimization_ 8:107–138.
  * Saberi et al. (2020) Saberi M, Hamedmoghadam H, Ashfaq M, Hosseini SA, Gu Z, Shafiei S, Nair DJ, Dixit V, Gardner L, Waller ST, et al. (2020) A simple contagion process describes spreading of traffic jams in urban networks. _Nature communications_ 11(1):1–9.
  * Sager et al. (2012) Sager S, Bock HG, Diehl M (2012) The integer approximation error in mixed-integer optimal control. _Mathematical programming_ 133(1-2):1–23.
  * Saxena et al. (2010) Saxena A, Bonami P, Lee J (2010) Convex relaxations of non-convex mixed integer quadratically constrained programs: extended formulations. _Mathematical programming_ 124(1-2):383–411.
  * Saxena et al. (2011) Saxena A, Bonami P, Lee J (2011) Convex relaxations of non-convex mixed integer quadratically constrained programs: projected formulations. _Mathematical programming_ 130:359–413.
  * Shah and Zaman (2016) Shah D, Zaman T (2016) Finding rumor sources on random trees. _Operations research_ 64(3):736–755.
  * Shen et al. (2014) Shen W, Duenyas I, Kapuscinski R (2014) Optimal pricing, production, and inventory for new product diffusion under supply constraints. _Manufacturing & Service Operations Management_ 16(1):28–45.
  * Sherali and Adams (2013) Sherali HD, Adams WP (2013) _A reformulation-linearization technique for solving discrete and continuous nonconvex problems_ , volume 31 (Springer Science & Business Media).
  * Sinclair et al. (2022) Sinclair SR, Banerjee S, Yu CL (2022) Adaptive discretization in online reinforcement learning. _Operations Research_ .
  * Sterman et al. (2012) Sterman J, Fiddaman T, Franck TR, Jones A, McCauley S, Rice P, Sawin E, Siegel L (2012) Climate interactive: the c-roads climate policy model. _System Dynamics Review_ 28(3):295–305.
  * Susarla et al. (2012) Susarla A, Oh JH, Tan Y (2012) Social networks and the diffusion of user-generated content: Evidence from youtube. _Information systems research_ 23(1):23–41.
  * Tawarmalani and Sahinidis (2004) Tawarmalani M, Sahinidis NV (2004) Global optimization of mixed-integer nonlinear programs: A theoretical and computational study. _Mathematical programming_ 99(3):563–591.
  * Walker et al. (2020) Walker PG, Whittaker C, Watson OJ, Baguelin M, Winskill P, Hamlet A, Djafaara BA, Cucunubá Z, Olivera Mesa D, Green W, et al. (2020) The impact of covid-19 and strategies for mitigation and suppression in low-and middle-income countries. _Science_ 369(6502):413–422.
  * Yamin and Gavious (2013) Yamin D, Gavious A (2013) Incentives’ effect in influenza vaccination policy. _Management Science_ 59(12):2667–2686.
  * Zhang et al. (2022a) Zhang J, Savin S, Veeraraghavan S (2022a) Revenue management in crowdfunding. _Manufacturing & Service Operations Management_ .
  * Zhang et al. (2022b) Zhang M, Ahn HS, Uichanco J (2022b) Data-driven pricing for a new product. _Operations Research_ 70(2):847–866.

##  7 Details on the prescriptive contagion models

In this appendix, we first formulate the four prescriptive contagion analytics
problems as subcases of the general problem presented in Section 3. We then
provide details on the contagion model for the congestion mitigation problem,
which we developed and calibrated using real-world data from Singapore.
Finally, we apply the set partitioning formulation proposed in Section 4.1 to
each problem.

###  7.1 Formulation of prescriptive contagion models

#### Vaccine allocation.

A centralized decision-maker distributes $B_{s}$ vaccine doses each week. Each
decision epoch $s=1,\cdots,S$ represents a week. Each segment represents one
of 51 US regions (50 states plus Washington, DC). The dynamics in each region
are characterized by a contagion model, shown in Figure 2a. This model
comprises seven states among the non-vaccinated population: susceptible (S),
exposed (E), infected (I), undetected (U), hospitalizations (H), quarantine
(Q), and death (D). In addition, the model includes four vaccinated states:
susceptible (S’), exposed (E’), infected (I’), and immunized (M). The core
model, called DELPHI, was fitted from historical data at the state level in
the United States; it was incorporated into the ensemble forecast from the
U.S. Center for Disease Control, with strong predictive performance across the
various waves of the pandemic (Li et al. 2022). It was then augmented into the
DELPHI-V model to reflect the impacts of vaccinations by Bertsimas et al.
(2022); we use this model in this paper, with the exception that we restrict
ourselves to a single age group for tractability purposes.

The decision variable $x_{is}$ denotes the number of vaccines allocated to
region $i=1,\cdots,n$ at epoch $s=1,\cdots,S$. Because vaccines are sent by
pallets, this amount must be a multiple of the number of vaccines per pallet,
denoted by $L$. For convenience, let us denote by $\overline{x}_{i}(t)$ a
piece-wise constant function such that
$\overline{x}_{i}(t)=\frac{x_{is}}{\tau_{s+1}-\tau_{s}}$ for all $\tau_{s}\leq
t\leq\tau_{s+1}$. Vaccinations induce a transfer of the susceptible population
to the alternative susceptible state, proportionally to an effectiveness
parameter $\beta$. This model assumes conservatively that vaccinations do not
prevent infections but immunize patients, thus preventing downstream
mortality.

The system dynamics are then governed by the following ordinary differential
equations:

$\displaystyle\frac{dS_{i}(t)}{dt}$
$\displaystyle=-\alpha\gamma(t)\left(S_{i}(t)-\beta\overline{x}_{i}(t)\right)I_{i}(t)-\beta\overline{x}_{i}(t)$
(40) $\displaystyle\frac{dE_{i}(t)}{dt}$
$\displaystyle=\alpha\gamma(t)\left(S_{i}(t)-\beta\overline{x}_{i}(t)\right)I_{i}(t)-r^{I}E_{i}(t)$
(41) $\displaystyle\frac{dI_{i}(t)}{dt}$
$\displaystyle=r^{I}E_{i}(t)-r^{d}I_{i}(t)$ (42)
$\displaystyle\frac{dU_{i}(t)}{dt}$
$\displaystyle=r_{k}^{U}(t)I_{i}(t)-r^{D}U_{i}(t)$ (43)
$\displaystyle\frac{dH_{i}(t)}{dt}$
$\displaystyle=r_{k}^{H}(t)I_{i}(t)-r^{D}H_{i}(t)$ (44)
$\displaystyle\frac{dQ_{i}(t)}{dt}$
$\displaystyle=r_{k}^{Q}(t)I_{i}(t)-r^{D}Q_{i}(t)$ (45)
$\displaystyle\frac{dD_{i}(t)}{dt}$
$\displaystyle=r^{D}\left(U_{i}(t)+H_{i}(t)+Q_{i}(t)\right)$ (46)
$\displaystyle\frac{dM_{i}(t)}{dt}$ $\displaystyle=\beta\overline{x}_{i}(t)$
(47)

The vaccine allocation problem is defined as follows. Equation (48) minimizes
the number of deaths (an absorbing state) along with the terminal number of
hospitalized, quarantined and exposed people. Constraint (49) applies the
budget of vaccines in each period. The remaining constraints define the
contagion dynamics and apply restrictions on the domain of definition of the
variables.

$\displaystyle\min\quad$
$\displaystyle\sum_{i=1}^{n}\left(c_{D}D_{i}(T)+c_{H}H_{i}(T)+c_{Q}Q_{i}(T)+c_{E}E_{i}(T)\right)$
(48) s.t. $\displaystyle\sum_{i=1}^{n}x_{is}\leq B_{s}\quad\forall
s=1,\cdots,S$ (49) Equations (40)–(47) (50) Initial conditions (51)
$\displaystyle\bm{S},\bm{E},\bm{I},\bm{U},\bm{H},\bm{Q},\bm{D},\bm{M}\geq\bm{0}$
(52) $\displaystyle x_{is}\in\\{0,L,2L,\cdots\\},\ \forall i=1,\cdots,n,\
\forall s=1,\cdots,S$ (53)

In fact, we maximize the number of lives saved by the vaccination campaign,
rather than minimizing the death toll. Indeed, system dynamics may lead to
“unavoidable” costs due to the initial conditions or to limited mitigation
resources. To capture the effects of the vaccination campaign, we optimize the
savings from a do-nothing baseline, with no vaccines. Obviously, the two
formulations are equivalent, but maximizing the cost savings provides better
estimates of the optimality gaps. Specifically, we denote by
$\overline{D}_{i}$, $\overline{H}_{i}$, $\overline{Q}_{i}$ and
$\overline{E}_{i}$ the states under the do-nothing baseline, and formulate the
problem as follows:

$\displaystyle\max\quad$
$\displaystyle\sum_{i=1}^{n}\left(c_{D}\overline{D}_{i}(T)+c_{H}\overline{H}_{i}(T)+c_{Q}\overline{Q}_{i}(T)+c_{E}\overline{E}_{i}(T)-c_{D}D_{i}(T)-c_{H}H_{i}(T)-c_{Q}Q_{i}(T)-c_{E}E_{i}(T)\right)$
(54) s.t. Equations (49)–(53) (55)

#### Vaccination centers.

We consider the same system dynamics as in the vaccine allocation problem, but
optimize the deployment of mass vaccination centers. Specifically, we select a
subset of $K$ vaccination centers out of a set of $F$ candidates. At each
epoch $s=1,\cdots,S$, a budget of $B_{s}$ vaccines needs to be allocated
across all regions. Each facility $j=1,\cdots,F$ has access to a maximum
capacity of $Cap_{js}$ vaccines at decision epoch $s=1,\cdots,S$—this capacity
is not necessarily equal to $B_{s}/K$ to create flexibility in vaccine
allocations. We impose a threshold on the maximum distance that a resident can
travel to reach a vaccination center. Let $P_{ij}$ denote the number of
residents in region $i$ than can be served by vaccination center $j$. For
example, if vaccination center $j$ is located in New York City and $i$ refers
to New Jersey, then $P_{ij}$ captures the number of New Jersey residents
living within the distance threshold from New York City.

We define the following variables:

$\displaystyle y_{j}$ $\displaystyle=\begin{cases}1&\text{if vaccination
center $j=1,\cdots,F$ is constructed}\\\ 0&\text{otherwise}\end{cases}$
$\displaystyle x_{ijs}$ $\displaystyle=\text{number of vaccines sent to region
$i=1,\cdots,n$ from vaccination center $j=1,\cdots,F$ at epoch
$s=1,\cdots,S$}$

The problem minimizes the number of deaths (Equation (56)). Constraint (57)
ensures that $K$ vaccination centers are selected. Constraint (58) applies the
total budget of vaccines. Constraint (59) applies the vaccine capacity in each
center, and constraint (60) ensures that vaccine allocation decisions are
consistent with the coverage of each vaccination center. Finally, constraint
(61) couples facility location decisions and subsequent vaccine allocation
decisions. This problem mirrors the problem studied by Bertsimas et al. (2022)
in the midst of the COVID-19 pandemic in the United States, again with the
restriction to a single age group.

$\displaystyle\min\quad$
$\displaystyle\sum_{i=1}^{n}\left(c_{D}D_{i}(T)+c_{H}H_{i}(T)+c_{Q}Q_{i}(T)+c_{E}E_{i}(T)\right)$
(56) s.t. $\displaystyle\sum_{j=1}^{F}y_{j}=K$ (57)
$\displaystyle\sum_{i=i}^{n}\sum_{j=1}^{F}x_{ijs}\leq B_{s}\quad\forall
s=1,\cdots,S$ (58) $\displaystyle\sum_{i=1}^{n}x_{ijs}\leq
Cap_{js}\quad\forall s=1,\cdots,S,\quad\forall j=1,\cdots,F$ (59)
$\displaystyle\sum_{s=1}^{S}x_{ijs}\leq P_{ij}\quad\forall
i=1,\cdots,n,\quad\forall j=1,\cdots,F$ (60) $\displaystyle
x_{ijs}\leq\min(Cap_{js},P_{ij})y_{j}\quad\forall i=1,\cdots,n,\quad\forall
j=1,\cdots,F,\quad\forall s=1,\cdots,S$ (61) Equations (40)–(47) (62) Initial
conditions (63)
$\displaystyle\bm{S},\bm{E},\bm{I},\bm{U},\bm{H},\bm{Q},\bm{D},\bm{M}\geq\bm{0}$
(64) $\displaystyle x_{ijs}\in\\{0,L,2L,\cdots\\},\ \forall i=1,\cdots,n,\
\forall j=1,\cdots,F,\ \forall s=1,\cdots,S$ (65)
$\displaystyle\bm{y}\in\\{0,1\\}^{F}$ (66)

Again, we actually maximize the number of lives saved, as follows:

$\displaystyle\max\quad$
$\displaystyle\sum_{i=1}^{n}\left(c_{D}\overline{D}_{i}(T)+c_{H}\overline{H}_{i}(T)+c_{Q}\overline{Q}_{i}(T)+c_{E}\overline{E}_{i}(T)-c_{D}D_{i}(T)-c_{H}H_{i}(T)-c_{Q}Q_{i}(T)-c_{E}E_{i}(T)\right)$
(67) s.t. Equations (57)–(66) (68)

#### Content promotion.

A platform needs to promote a portfolio of $n$ products. For implementation
reasons, we assume that the platform can promote each product to a fraction of
the population by increments of a constant fraction $L$ and can only promote
up to $K$ products at a time. We define the following variables:

$\displaystyle x^{1}_{is}$ $\displaystyle=\begin{cases}1&\text{if product
$i=1,\cdots,n$ is promoted at epoch $s=1,\cdots,S$}\\\
0&\text{otherwise}\end{cases}$ $\displaystyle x^{2}_{is}$