$\displaystyle=\text{population share to which product $i=1,\cdots,n$ is
promoted at epoch $s=1,\cdots,S$}$

Again, we define a piece-wise constant function
$\overline{x}^{2}_{i}(t)=\frac{x^{2}_{is}}{\tau_{s+1}-\tau_{s}}$ for all
$\tau_{s}\leq t\leq\tau_{s+1}$.

The contagion model captures product adoption dynamics via a Bass model,
comprising susceptible users (A) and adopters (B). Product promotion increases
adoption from external sources, which, in turn, increases positive
externalities between adopters and susceptible users (Figure 2b). This model
follows the one from Lin et al. (2021), developed using data from a social
media platform. The dynamics are defined as follows:

$\displaystyle\frac{dA_{i}(t)}{dt}$
$\displaystyle=-\alpha\overline{x}^{2}_{i}(t)A_{i}(t)-\frac{\beta}{m}S_{i}(t)B_{i}(t)$
(69) $\displaystyle\frac{dB_{i}(t)}{dt}$
$\displaystyle=\alpha\overline{x}^{2}_{i}(t)A_{i}(t)+\frac{\beta}{m}S_{i}(t)B_{i}(t)$
(70)

The problem maximizes total adoption (Equation (71)), subject to a sparsity
constraint ensuring that the platform promotes at most $K$ products per period
(Constraint (72)), to a cover constraint ensuring that every user is shown a
product (Constraint (73)), to a linking constraint ensuring consistency
between the $\bm{x}^{1}$ and $\bm{x}^{2}$ variables (Constraint (74)), and to
the contagion dynamics.

$\displaystyle\max\quad$ $\displaystyle\sum_{i=1}^{n}B_{i}(T)$ (71) s.t.
$\displaystyle\sum_{i=1}^{n}x^{1}_{is}\leq K\quad\forall s=1,\cdots,S$ (72)
$\displaystyle\sum_{i=1}^{n}x^{2}_{is}=1\quad\forall s=1,\cdots,S$ (73)
$\displaystyle x^{2}_{is}\leq x^{1}_{is}\quad\forall s=1,\cdots,S,\ \forall
i=1,\cdots,n$ (74) Equations (69)–(70) (75) Initial conditions (76)
$\displaystyle\bm{S},\bm{B}\geq\bm{0}$ (77) $\displaystyle
x^{2}_{is}\in\\{0,L,2L,\cdots,1\\},\ \forall i=1,\cdots,n,\ \forall
s=1,\cdots,S$ (78) $\displaystyle\bm{x}^{1}\in\\{0,1\\}^{n\times S}$ (79)

Again, we maximize the improvements from a do-nothing baseline, with no
product promotion. We denote by $\overline{B}_{i}$ the adoption state under
the do-nothing baseline, and formulate the problem as follows:

$\displaystyle\max\quad$
$\displaystyle\sum_{i=1}^{n}(B_{i}(T)-\overline{B}_{i}(T))$ (80) s.t.
Equations (72)–(79) (81)

#### Congestion mitigation.

We consider $n$ neighborhoods, and each population corresponds to the set of
road segments in each neighborhood. The contagion model comprises six states,
shown in Figure 2c: susceptible free-flow roads (S), road work (W), roads with
accidents (A), roads with road work and accidents ($A^{\prime}$), congested
roads (I), and recovered roads (R). Details on this model are provided in
Appendix 7.2.

The centralized decision-maker can deploy $B^{1}$ emergency vehicles for
treatment purposes (namely, to respond to major accidents in order to clear
the roads faster) and $B^{2}$ emergency vehicles for prevention purposes
(namely, to respond to minor accidents and prevent the formation of
congestion). Accordingly, we define the following variables:

$\displaystyle x^{1}_{is}$ $\displaystyle=\text{number of treatment vehicles
allocated to neighborhood $i=1,\cdots,n$ at epoch $s=1,\cdots,S$}$
$\displaystyle x^{2}_{is}$ $\displaystyle=\text{number of prevention vehicles
allocated to neighborhood $i=1,\cdots,n$ at epoch $s=1,\cdots,S$}$

Again, we define piece-wise constant functions
$\overline{x}^{1}_{i}(t)=\frac{x^{1}_{is}}{\tau_{s+1}-\tau_{s}}$ and
$\overline{x}^{2}_{i}(t)=\frac{x^{2}_{is}}{\tau_{s+1}-\tau_{s}}$ for all
$\tau_{s}\leq t\leq\tau_{s+1}$.

The system in each neighborhood is governed by the following equations, where
$\phi(\cdot,\cdot)$ and $\psi(\cdot,\cdot)$ denote generic functions
reflecting the effects of the prevention and treatment interventions:

$\displaystyle\frac{dS_{i}(t)}{dt}=-\alpha^{F}S_{i}(t)(I_{i}(t)+A_{i}(t)+A^{\prime}_{i}(t))-\beta^{F}S_{i}(t)+\rho^{F}(\zeta^{F}+\psi(\overline{x}^{1}_{i}(t),\zeta^{F}))A_{i}(t)$
(82)
$\displaystyle\frac{dA_{i}(t)}{dt}=\beta^{F}S_{i}(t)-\rho^{F}(\zeta^{F}+\psi(\overline{x}^{1}_{i}(t),\zeta^{F}))A_{i}(t)-(1-\rho^{F})(\zeta^{F}+\phi(\overline{x}^{2}_{i}(t),\zeta^{F}))A_{i}(t)$
(83)
$\displaystyle\frac{dW_{i}(t)}{dt}=-\alpha^{W}W_{i}(t)(I_{i}(t)+A_{i}(t)+A^{\prime}_{i}(t))-\beta^{W}W_{i}(t)+\rho^{W}(\zeta^{W}+\psi(\overline{x}^{1}_{i}(t),\zeta^{W}))A^{\prime}_{i}(t)$
(84)
$\displaystyle\frac{dA^{\prime}_{i}(t)}{dt}=\beta^{W}W_{i}(t)-\rho^{W}(\zeta^{W}+\psi(\overline{x}^{1}_{i}(t),\zeta^{W}))A^{\prime}_{i}(t)-(1-\rho^{W})(\zeta^{W}+\phi(\overline{x}^{2}_{i}(t),\zeta^{W}))A^{\prime}_{i}(t)$
(85) $\displaystyle\frac{dR_{i}(t)}{dt}=\theta I_{i}(t)$ (86)
$\displaystyle\frac{dI_{i}(t)}{dt}=-\left(\frac{dS_{i}(t)}{dt}+\frac{dW_{i}(t)}{dt}+\frac{dA_{i}(t)}{dt}+\frac{dA^{\prime}_{i}(t)}{dt}+\frac{dR_{i}(t)}{dt}\right)$
(87)

In our experiments, we consider linear functions $\phi(\cdot,\cdot)$ and
$\psi(\cdot,\cdot)$ defined as follows, although our model can accommodate
more complex and non-linear functions.

$\displaystyle\phi(\overline{x}^{2}_{i}(t),\zeta)=\frac{n\zeta}{2D_{y}}\left(\overline{x}^{2}_{i}(t)-\frac{B^{2}}{n}\right)$
(88)
$\displaystyle\psi(\overline{x}^{1}_{i}(t),\zeta)=\frac{n\zeta}{2D_{x}}\left(\overline{x}^{1}_{i}(t)-\frac{B^{1}}{n}\right),$
(89)

The problem minimizes the costs of congestion and accidents, subject to budget
constraints and congestion dynamics governed by the contagion model:

$\displaystyle\min\quad$
$\displaystyle\sum_{i=1}^{n}\int_{0}^{T}\left(c_{I}I_{i}(t)+c_{A}A_{i}(t)+c_{A^{\prime}}A^{\prime}_{i}(t)\right)dt$
(90) s.t. $\displaystyle\sum_{i=1}^{n}x^{1}_{is}\leq B^{1}\quad\forall
s=1,\cdots,S$ (91) $\displaystyle\sum_{i=1}^{n}x^{2}_{is}\leq
B^{2}\quad\forall s=1,\cdots,S$ (92) Equations (82)–(87) (93) Initial
conditions (94)
$\displaystyle\bm{S},\bm{W},\bm{A},\bm{A}^{\prime},\bm{I},\bm{R}\geq\bm{0}$
(95) $\displaystyle\bm{x}^{1},\bm{x}^{2}\in\mathbb{Z}_{+}^{n\times T}$ (96)

Again, we maximize the improvements from a do-nothing baseline, with no
emergency vehicle. We denote by $\overline{I}_{i}$, $\overline{A}_{i}$ and
$\overline{A}^{\prime}_{i}$ the states under the do-nothing baseline, and
formulate the problem as follows:

$\displaystyle\max\quad$
$\displaystyle\sum_{i=1}^{n}\int_{0}^{T}\left(c_{I}\overline{I}_{i}(t)+c_{A}\overline{A}_{i}(t)+c_{A^{\prime}}\overline{A}^{\prime}_{i}(t)-c_{I}I_{i}(t)-c_{A}A_{i}(t)-c_{A^{\prime}}A^{\prime}_{i}(t)\right)dt$
(97) s.t. Equations (91)–(96) (98)

###  7.2 A data-driven contagion model of traffic congestion

We collected data in Singapore that tracks the average speed on each road
every five minutes, as well as the log of incidents and road works in the
city. We divide the city into five neighborhoods, and the population in each
neighborhood characterizes all the road segments. In each five-minute
interval, a road segment is labeled as congested if the average speed recorded
is less than 30% of the average speed on that road.

We considered three contagion models: a basic Susceptible-Infected-Recovered
(SIR) model inspired from Saberi et al. (2020), an SAIR model with an
“accident” state, and the full model depicted in Figure 2c where we separate
road segments with and without road work (referred to as 2-SAIR). We divide
each day into a morning period (7am to 1pm), an afternoon period (1pm to
5:30pm) and an evening period (5:30pm onward). We trained our contagion models
on those three periods separately, using historical data from March 2021. We
then tested the models on data from April 2021. The results are reported in
Figure 6 and Table 6.

(a) In sample (March 2021)

(b) Out of sample (April 01, 2021)

Figure 6: Comparison of congestion levels predicted by the three contagion models. Table 6: Root-mean-square error (RMSE) for each model on the three periods of the day. | In sample (March 2021 | Out of sample (April 2021)
---|---|---
Model | Morning | Afternoon | Evening | Morning | Afternoon | Evening
SIR | 0.0551 | 0.0546 | 0.0586 | 0.0526 | 0.1028 | 0.2196
SAIR | 0.0506 | 0.0502 | 0.0702 | 0.0434 | 0.1002 | 0.205
2-SAIR | 0.0389 | 0.0501 | 0.0533 | 0.0406 | 0.0998 | 0.1885

Note that the 2-SAIR model consistently outperforms the two benchmarks, both
in sample and out of sample—with a 23% reduction in RMSE during the morning
peak and a 14% reduction during the evening peak. Its parameters, reported in
Table 7, provide some insights on the underlying dynamics; for instance, the
accident rate and the congestion rate are both higher on roads segments with
road work than on roads segments without road work. Ultimately, the model
provides strong estimates of the congestion dynamics observed in the city, as
illustrated in Figure 6 across the training data and for one out-of-sample
day.

Table 7: Parameters of the 2-SAIR model for every region of the city. Time period | Region | $\alpha^{F}$ | $\beta^{F}$ | $\rho^{F}$ | $\zeta^{F}$ | $\alpha^{W}$ | $\beta^{W}$ | $\rho^{W}$ | $\zeta^{W}$ | $\theta$
---|---|---|---|---|---|---|---|---|---|---
Morning | West | 0.531 | 0.578 | 1.126 | 0.569 | 3.966 | 1.585 | 0.411 | 4.835 | 10.354
| Central | 0.000 | 0.048 | 0.034 | 0.076 | 1.224 | 0.875 | 1.180 | 1.544 | 2.135
| Northeast | 0.000 | 0.137 | 0.345 | 0.078 | 1.240 | 0.198 | 1.023 | 2.273 | 4.742
| East | 0.052 | 0.041 | 0.432 | 0.086 | 2.745 | 3.086 | 1.201 | 1.934 | 3.950
| North | 0.492 | 0.562 | 1.450 | 2.043 | 1.396 | 0.896 | 2.162 | 1.878 | 1.113
Afternoon | West | 0.897 | 0.180 | 1.913 | 3.180 | 2.852 | 1.321 | 1.208 | 3.444 | 2.512
| Central | 0.862 | 0.070 | 3.259 | 1.144 | 1.098 | 1.026 | 0.338 | 1.831 | 3.477
| Northeast | 0.782 | 0.024 | 5.786 | 0.283 | 1.215 | 2.360 | 0.002 | 3.555 | 1.381
| East | 0.838 | 0.115 | 2.143 | 1.260 | 0.000 | 0.288 | 0.000 | 0.137 | 1.283
| North | 1.538 | 0.196 | 1.856 | 1.656 | 0.987 | 1.085 | 0.996 | 0.930 | 1.111
Evening | West | 0.091 | 0.448 | 1.249 | 1.176 | 0.568 | 0.990 | 2.194 | 1.264 | 3.775
| Central | 0.024 | 0.331 | 1.176 | 0.793 | 0.068 | 0.733 | 1.873 | 1.335 | 3.350
| Northeast | 0.018 | 0.304 | 1.247 | 1.326 | 0.116 | 1.514 | 1.728 | 1.123 | 2.861
| East | 0.603 | 0.191 | 2.255 | 1.141 | 0.962 | 0.985 | 1.527 | 1.517 | 6.799
| North | 0.001 | 0.580 | 0.905 | 0.698 | 1.148 | 1.110 | 1.136 | 1.142 | 4.264

###  7.3 Set partitioning formulation applied to prescriptive contagion models

#### Vaccine allocation.

A plan $p$ is a list of $S$ elements $\alpha_{is}^{p}$, where
$\alpha_{is}^{p}$ represents the number of vaccines allocated to region $i$ at
epoch $s$ in plan $p$. The set partitioning formulation is given as follows:

$\displaystyle\min\quad$
$\displaystyle\sum_{i=1}^{n}\sum_{p\in\mathcal{P}_{i}}C^{p}_{i}z^{p}_{i}$ (99)
s.t $\displaystyle\sum_{i=1}^{n}\alpha_{is}^{p}z_{i}^{p}\leq
B_{s},\quad\forall s=1,\cdots,S$ (100)
$\displaystyle\sum_{p\in\mathcal{P}_{i}}z_{i}^{p}=1,\quad\forall i=1,\cdots,n$
(101) $\displaystyle z^{p}_{i}\in\\{0,1\\}\quad\forall
p\in\mathcal{P}_{i},\quad\forall i=1,\cdots,n$ (102)

Let $\lambda_{s}\in\mathbb{R}_{+}$ and $\mu_{i}\in\mathbb{R}$ denote the dual
variables associated with Constraints (100) and (101), respectively. We obtain
the following pricing problem:

$\displaystyle\min\quad$
$\displaystyle\left(c_{D}D_{i}(T)+c_{H}H_{i}(T)+c_{Q}Q_{i}(T)+c_{E}E_{i}(T)\right)+\sum_{s=1}^{S}\lambda_{s}x_{is}-\mu_{i}$
(103) s.t. $\displaystyle x_{is}\leq B_{s}\quad\forall s=1,\cdots,S$ (104)
Equations (40)–(47) (105) Initial conditions (106)
$\displaystyle\bm{S},\bm{E},\bm{I},\bm{U},\bm{H},\bm{Q},\bm{D},\bm{M}\geq\bm{0}$
(107) $\displaystyle x_{is}\in\\{0,L,2L,\cdots\\},\ \forall i=1,\cdots,n,\
\forall s=1,\cdots,S$ (108)

#### Vaccination centers.

A plan $p$ is defined as above. The set partitioning formulation is given as
follows:

$\displaystyle\min\quad$
$\displaystyle\sum_{i=1}^{n}\sum_{p\in\mathcal{P}_{i}}C^{p}_{i}z^{p}_{i}$
(109) s.t $\displaystyle\sum_{j=1}^{F}y_{j}=K$ (110)
$\displaystyle\sum_{i=i}^{n}\sum_{j=1}^{F}x_{ijs}\leq B_{s}\quad\forall
s=1,\cdots,S$ (111) $\displaystyle\sum_{i=1}^{n}x_{ijs}\leq
Cap_{js}\quad\forall s=1,\cdots,S,\quad\forall j=1,\cdots,F$ (112)
$\displaystyle\sum_{s=1}^{S}x_{ijs}\leq P_{ij}\quad\forall
i=1,\cdots,n,\quad\forall j=1,\cdots,F$ (113) $\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$ (114)
$\displaystyle\sum_{p\in\mathcal{P}_{i}}\alpha_{is}^{p}z_{i}^{p}=\sum_{j=1}^{K}x_{i,j,s}\quad\forall
s=1,\cdots,S\quad\forall i=1,\cdots,n$ (115)
$\displaystyle\sum_{p\in\mathcal{P}_{i}}z_{i}^{p}=1\quad\forall i=1,\cdots,n$
(116) $\displaystyle z^{p}_{i}\in\\{0,1\\}\quad\forall p\in\mathcal{P}_{i},\
\forall i=1,\cdots,n$ (117) $\displaystyle y_{j}\in\\{0,1\\}\quad\forall
j=1,\cdots,F$ (118)

Note that constraint (115) ensures consistency between the number of vaccines
sent to a region $i$ and the plan-based variables. Let
$\lambda_{is}\in\mathbb{R}$ and $\mu_{i}\in\mathbb{R}$ denote the dual
variables associated with Constraints (115) and (116), respectively. We obtain
the following pricing problem:

$\displaystyle\min\quad$
$\displaystyle\left(c_{D}D_{i}(T)+c_{H}H_{i}(T)+c_{Q}Q_{i}(T)+c_{E}E_{i}(T)\right)-\sum_{s=1}^{S}\lambda_{is}x_{is}-\mu_{i}$
(119) s.t. $\displaystyle x_{is}\leq B_{s}\quad\forall s=1,\cdots,S$ (120)
Equations (40)–(47) (121) Initial conditions (122)
$\displaystyle\bm{S},\bm{E},\bm{I},\bm{U},\bm{H},\bm{Q},\bm{D},\bm{M}\geq\bm{0}$
(123) $\displaystyle x_{is}\in\\{0,L,2L,\cdots\\},\ \forall i=1,\cdots,n,\
\forall s=1,\cdots,S$ (124)

#### Content promotion.

A plan $p$ is a list of $S$ tuples $(\alpha^{1p}_{is},\alpha^{2p}_{is})$ where
$\alpha^{1p}_{is}$ is equal to 1 if plan $p$ promotes product $i$ at epoch $s$
and $\alpha^{2p}_{is}$ is the share of the population to which product $i$ is
shown at epoch $s$ under plan $p$. The set partitioning formulation is given
by:

$\displaystyle\min\quad$
$\displaystyle\sum_{i=1}^{n}\sum_{p\in\mathcal{P}_{i}}C^{p}_{i}z^{p}_{i}$
(125) s.t $\displaystyle\sum_{i=1}^{n}\alpha^{1p}_{is}z_{i}^{p}\leq
K,\quad\forall s=1,\cdots,S$ (126)
$\displaystyle\sum_{i=1}^{n}\alpha^{2p}_{is}z_{i}^{p}=1,\quad\forall
s=1,\cdots,S$ (127)
$\displaystyle\sum_{p\in\mathcal{P}_{i}}z_{i}^{p}=1,\quad\forall i=1,\cdots,n$
(128) $\displaystyle z^{p}_{i}\in\\{0,1\\}\quad\forall
p\in\mathcal{P}_{i},\quad\forall i=1,\cdots,n$ (129)

Let $\lambda^{1}_{s}\in\mathbb{R}_{+}$, $\lambda^{2}_{s}\in\mathbb{R}$ and
$\mu_{i}\in\mathbb{R}$ denote the dual variables associated with Constraints
(126), (127) and (128), respectively. We obtain the following pricing problem:

$\displaystyle\max\quad$ $\displaystyle
B_{i}(T)+\sum_{s=1}^{S}\lambda^{1}_{s}x^{1}_{is}-\sum_{s=1}^{S}\lambda^{2}_{s}x^{2}_{is}-\mu_{i}$
(130) $\displaystyle x^{2}_{is}\leq x^{1}_{is}\quad\forall s=1,\cdots,S$ (131)
Equations (69)–(70) (132) $\displaystyle x^{2}_{is}\in\\{0,L,2L,\cdots,1\\},\
\forall i=1,\cdots,n,\ \forall s=1,\cdots,S$ (133)
$\displaystyle\bm{x}^{1}\in\\{0,1\\}^{n\times S}$ (134)

#### Congestion mitigation.

A plan $p$ is a list of $S$ tuples $(\alpha^{1p}_{is},\alpha^{2p}_{is})$ where
$\alpha^{1p}_{is}$ and $\alpha^{2p}_{is}$ are the number of emergency vehicles
deployed for treatment and prevention purposes, respectively, in segment $i$
at epoch $s$ under plan $p$. The set partitioning formulation is given by:

$\displaystyle\min\quad$
$\displaystyle\sum_{i=1}^{n}\sum_{p\in\mathcal{P}_{i}}C^{p}_{i}z^{p}_{i}$
(135) s.t $\displaystyle\sum_{i=1}^{n}\alpha^{1p}_{is}z_{i}^{p}\leq
B^{1},\quad\forall s=1,\cdots,S$ (136)
$\displaystyle\sum_{i=1}^{n}\alpha^{2p}_{is}z_{i}^{p}\leq B^{2},\quad\forall
s=1,\cdots,S$ (137)
$\displaystyle\sum_{p\in\mathcal{P}_{i}}z_{i}^{p}=1,\quad\forall i=1,\cdots,n$
(138) $\displaystyle z^{p}_{i}\in\\{0,1\\}\quad\forall
p\in\mathcal{P}_{i},\quad\forall i=1,\cdots,n$ (139)

Let $\lambda^{1}_{s}\in\mathbb{R}_{+}$, $\lambda^{2}_{s}\in\mathbb{R}_{+}$ and
$\mu_{i}\in\mathbb{R}$ denote the dual variables associated with Constraints
(136), (137) and (138), respectively. We obtain the following pricing problem:

$\displaystyle\min\quad$
$\displaystyle\int_{0}^{T}\left(c_{I}I_{i}(t)+c_{A}A_{i}(t)+c_{A^{\prime}}A^{\prime}_{i}(t)\right)dt+\sum_{s=1}^{S}\lambda^{1}_{s}x^{1}_{is}+\lambda^{2}_{s}x^{2}_{is}-\mu_{i}$
(140) s.t. $\displaystyle\sum_{i=1}^{n}x^{1}_{is}\leq B^{1}\quad\forall
s=1,\cdots,S$ (141) $\displaystyle\sum_{i=1}^{n}x^{2}_{is}\leq
B^{2}\quad\forall s=1,\cdots,S$ (142) Equations (82)–(87) (143)
$\displaystyle\bm{S},\bm{W},\bm{A},\bm{A}^{\prime},\bm{I},\bm{R}\geq\bm{0}$
(144) $\displaystyle\bm{x}^{1},\bm{x}^{2}\in\mathbb{Z}_{+}^{n\times T}\ $
(145)

##  8 Technical details

###  8.1 Proof of Proposition 4.1

Let $\widehat{\bm{z}},\ \widehat{\bm{y}}$ be a feasible solution to
$(\mathcal{S}\mathcal{P})$. We define:

$\displaystyle\widehat{\bm{x}}_{is}=\sum_{p\in\mathcal{P}_{i}}\bm{\alpha}_{is}^{p}z_{i}^{p}\quad\forall
i=1,\dots,n,\ \forall s=1,\dots,S$

By definition of the costs in $(\mathcal{S}\mathcal{P})$ and in
$(\mathcal{P})$, these two solutions achieve the same costs and satisfy
constraints (3) and (4). Moreover, because of constraint (16) and by
definition of $\bm{\alpha}$, we know that, for all $i=1,...,n$, there exists a
single plan $p\in\mathcal{P}_{i}$ such that $\widehat{z}_{i}^{p}=1$ and that
$\bm{\alpha}_{is}^{p}\in\mathcal{F}_{is}$ for all $s=1,\dots,S$. Therefore,
$\widehat{\bm{x}}_{is}\in\mathcal{F}_{is}$ for all $i=1,\dots,n,$ and
$s=1,\dots,S$. Finally, using constraint (15), we have:

$\displaystyle\sum_{i=1}^{n}\bm{u}^{\top}_{sji}\widehat{\bm{x}}_{is}+\bm{v}_{sj}^{\top}\widehat{\bm{y}}=\sum_{i=1}^{n}\sum_{p\in\mathcal{P}_{i}}\bm{u}_{sji}^{\top}\bm{\alpha}_{is}^{p}\widehat{z}_{i}^{p}+\bm{v}_{sj}^{\top}\widehat{\bm{y}}\geq
w_{sj}\quad\forall s=1,\dots,S,\ \forall j=1,\dots,m_{s}$

Therefore, this yields a feasible solution to $(\mathcal{P})$ with the exact
same cost.

Conversely, let $\widehat{\bm{x}},\ \widehat{\bm{y}}$ be a solution to
$(\mathcal{P})$. We define the plan $\widehat{p}_{i}$ such that:

$\bm{\alpha}_{is}^{\widehat{p}_{i}}=\widehat{\bm{x}}_{is}\quad\forall
i=1,\dots,n,\ \forall s=1,\dots,S$

Now we consider the following plan-based solution, for all $i=1,\cdots,n$:

$\displaystyle z^{p}_{i}=\begin{cases}1&\text{if }p=\widehat{p}_{i}\\\
0&\text{otherwise}\end{cases}$

By construction, constraints (16) and (17) are satisfied. Moreover:

$\displaystyle\sum_{i=1}^{n}\sum_{p\in\mathcal{P}_{i}}\bm{u}_{sji}^{\top}\bm{\alpha}_{is}^{p}z_{i}^{p}+\bm{v}_{sj}^{\top}\widehat{\bm{y}}=\sum_{i=1}^{n}\bm{u}_{sji}^{\top}\bm{\alpha}_{is}^{\widehat{p}_{i}}z_{i}^{\widehat{p}_{i}}+\bm{v}_{sj}^{\top}\widehat{\bm{y}}=\sum_{i=1}^{n}\bm{u}_{sji}^{\top}\widehat{\bm{x}}_{is}+\bm{v}_{sj}^{\top}\widehat{\bm{y}}\geq
g_{j}^{s}\quad\forall s=1,\dots,S,\ \forall j=1,\dots,m_{s}$

We define the corresponding cost parameter $C^{\widehat{p}_{i}}_{i}$ per
Equations (11)–(13). By construction, the two solutions achieve the same costs
and satisfy constraints (3) and (4). This yields a feasible solution to
$(\mathcal{S}\mathcal{P})$ with the exact same cost. Therefore, Problem
$(\mathcal{P})$ and Problem $(\mathcal{S}\mathcal{P})$ are equivalent.

###  8.2 Pseudo-code of dynamic programming algorithms

Algorithm 2 provides pseudo-code for the exact dynamic programming algorithm
comprising the forward-enumeration and backward-induction loops. The forward-
enumeration loop returns the full state space, denoted by
$\mathcal{N}^{*}_{1},\cdots,\mathcal{N}^{*}_{S+1}$ in the paper. The backward-
induction loop computes the optimal policy given a state space. Algorithm 3
provides pseudo-code for the state-clustering algorithm. It returns the
clustered state space $\mathcal{N}_{1},\cdots,\mathcal{N}_{S+1}$. The
approximate dynamic programming algorithm considered in the paper consists of
the state-clustering procedure in Algorithm 3 and of the backward-induction
loop in Algorithm 2.

Algorithm 2 Dynamic programming algorithm for pricing problem

1:procedure StateSpaceCreation($\mathcal{N}_{1}$)

2: for $s\leftarrow 1$ to $S$ do $\triangleright$ Forward enumeration

3: $\mathcal{N}_{s+1}\leftarrow\emptyset$

4: for $\bm{N}_{s}\in\mathcal{N}_{s},\ \bm{x}_{s}\in\mathcal{F}_{s}$ do

5: Transition: $\bm{M}(\tau_{s})\leftarrow\bm{N}_{s}$;
$\frac{d\bm{M}(t)}{dt}=f(\bm{M}(t),\bm{x}_{s}),\ \forall
t\in\left[\tau_{s},\tau_{s+1}\right]$;
$\mathcal{N}_{s+1}\leftarrow\mathcal{N}_{s+1}\cup\left\\{\bm{M}(\tau_{s+1})\right\\}$

6: Cost:
$c_{s}(\bm{N}_{s},\bm{x}_{s})\leftarrow\int_{\tau_{s}}^{\tau_{s+1}}g_{t}(\bm{M}(t))dt+\Gamma_{s}(\bm{x}_{s})$

7: if $s=S$ then

8: $c_{s}(\bm{N}_{s},\bm{x}_{s})\leftarrow
c_{s}(\bm{N}_{s},\bm{x}_{s})+h(\bm{N}_{s})$

9: end if

10: end for

11: end for

12: return state space $\mathcal{N}_{1},\cdots,\mathcal{N}_{S+1}$, cost
$c_{s}(\bm{N}_{s},\bm{x}_{s})$ for all $\bm{N}_{s}\in\mathcal{N}_{s},\
\bm{x}_{s}\in\mathcal{F}_{s},\ s\in\\{1,\cdots,S\\}$

13:end procedure

14:procedure DynamicProgram($\mathcal{N}_{1},\cdots,\mathcal{N}_{S},\
c_{1}(\cdot,\cdot),\cdots,c_{S}(\cdot,\cdot),\lambda_{sj}$)

15: Initialization: $c^{opt}(\bm{N}_{s})=+\infty;\
\bm{\pi}^{opt}_{s}(\bm{N}_{s})\leftarrow\emptyset$ $\forall\bm{N}_{s}$
$\forall s\in\\{1,\cdots,S\\}$

16: for $s\leftarrow S$ to $1$ do $\triangleright$ Backward induction

17: for $\bm{N}_{s}\in\mathcal{N}_{s}$ do

18: $c^{temp}\leftarrow+\infty$; $\bm{x}^{temp}\leftarrow\emptyset$

19: for $\bm{x}_{s}\in\mathcal{F}_{s}$ do

20: $\widetilde{c}_{s}(\bm{N}_{s},\bm{x}_{s})\leftarrow
c_{s}(\bm{N}_{s},\bm{x}_{s})-\sum_{j=1}^{m_{s}}\lambda_{sj}\bm{u}^{\top}_{sj}\bm{x}_{s}-\mathbf{1}\\{s=S\\}\mu$

21: If $\widetilde{c}_{s}(\bm{N}_{s},\bm{x}_{s})<c^{temp}$, update
$c^{temp}\leftarrow\widetilde{c}_{s}(\bm{N}_{s},\bm{x}_{s})$;
$\bm{x}^{temp}\leftarrow\bm{x}_{s}$

22: end for

23: $c^{opt}(\bm{N}_{s})=c^{temp};\
\bm{\pi}^{opt}_{s}(\bm{N}_{s})\leftarrow\bm{x}^{temp}$

24: end for

25: end for

26: Initialization: $\bm{N}^{opt}_{1}\leftarrow\bm{M}^{0}$, $OPT\leftarrow 0$

27: for $s\leftarrow 1$ to $S$ do $\triangleright$ Forward evaluation

28: Find optimal decision
$\bm{x}^{opt}_{s}\leftarrow\bm{\pi}^{opt}_{s}(\bm{N}^{opt}_{s})$ and update
cost $OBJ\leftarrow OBJ+c_{s}(\bm{N}^{opt}_{s},\bm{x}^{opt}_{s})$

29: Find next state: $\bm{M}(\tau_{s})\leftarrow\bm{N}^{opt}_{s}$;
$\frac{d\bm{M}(t)}{dt}=f(\bm{M}(t),\bm{x}_{s}),\ \forall
t\in\left[\tau_{s},\tau_{s+1}\right]$;
$\bm{N}^{opt}_{s+1}\leftarrow\bm{M}(\tau_{s+1})$

30: end for

31: return cost $OBJ$, decisions $\bm{x}_{s}\;\forall s\in\\{1,\cdots,S\\}$

32:end procedure

Algorithm 3 State-clustering acceleration of dynamic programming algorithm for
pricing problem

1:procedure StateClustering($\mathcal{N}_{1}$, $\varepsilon$)

2: for $s\leftarrow 1$ to $S$ do $\triangleright$ Forward pass

3: Initialization: $\mathcal{N}_{s+1}\leftarrow\emptyset$, $k\leftarrow 0$

4: for $\bm{N}_{s}\in\mathcal{N}_{s},\ \bm{x}_{s}\in\mathcal{F}_{s}$ do

5: Next state: $\bm{M}(\tau_{s})\leftarrow\bm{N}_{s}$;
$\frac{d\bm{M}(t)}{dt}=f(\bm{M}(t),\bm{x}_{s}),\ \forall
t\in\left[\tau_{s},\tau_{s+1}\right]$;
$\bm{N}_{s+1}\leftarrow\bm{M}(\tau_{s+1})$

6:
$l^{*}\leftarrow\argmin_{l\in\\{0,1,\cdots,k\\}}\max(\|\bm{N}_{s+1}-\underline{\bm{N}}(\gamma_{l})\|_{\infty},\|\bm{N}_{s+1}-\overline{\bm{N}}(\gamma_{l})\|_{\infty})$

7: if
$\max(\|\bm{N}_{s+1}-\underline{\bm{N}}(\gamma_{l^{*}})\|_{\infty},\|\bm{N}_{s+1}-\overline{\bm{N}}(\gamma_{l^{*}})\|_{\infty})\leq\varepsilon$
then

8: Update:
$\bm{N}^{\Sigma}(\gamma_{l^{*}})\leftarrow\bm{N}^{\Sigma}(\gamma_{l^{*}})+\bm{N}_{s+1}$
, $\eta(\gamma_{l^{*}})\leftarrow\eta(\gamma_{l^{*}})+1$

9: Update:
$\underline{\bm{N}}(\gamma_{l^{*}})\leftarrow\min(\underline{\bm{N}}(\gamma_{l^{*}}),\bm{N}_{s+1})$,
$\overline{\bm{N}}(\gamma_{l^{*}})\leftarrow\max(\overline{\bm{N}}(\gamma_{l^{*}}),\bm{N}_{s+1})$

10: Update:
$\mathcal{X}(\gamma_{l^{*}})\leftarrow\mathcal{X}(\gamma_{l^{*}})\cup\\{(\bm{N}_{s},\bm{x}_{s})\\}$,
$c^{\Sigma}(\gamma_{l^{*}})\leftarrow
c^{\Sigma}(\gamma_{l^{*}})+\int_{\tau_{s}}^{\tau_{s+1}}g_{t}(\bm{M}(t))dt+\mathbf{1}\\{s=S\\}(h(\bm{M}(T)))$

11: else

12: New cluster: $k\leftarrow k+1$,
$\bm{N}^{\Sigma}(\gamma_{k})\leftarrow\bm{N}_{s+1}$,
$\eta(\gamma_{k})\leftarrow 1$

13: New cluster: $\underline{\bm{N}}(\gamma_{k})\leftarrow\bm{N}_{s+1}$,
$\overline{\bm{N}}(\gamma_{k})\leftarrow\bm{N}_{s+1}$

14: New cluster:
$\mathcal{X}(\gamma_{k})\leftarrow\\{(\bm{N}_{s},\bm{x}_{s})\\}$,
$c^{\Sigma}(\gamma_{l^{*}})\leftarrow\int_{\tau_{s}}^{\tau_{s+1}}g_{t}(\bm{M}(t))dt+\mathbf{1}\\{s=S\\}(h(\bm{M}(T)))$

15: end if

16: end for

17: for $l\leftarrow 1$ to $k$ do

18:
$\mathcal{N}_{s+1}\leftarrow\mathcal{N}_{s+1}\cup\left\\{\bm{N}^{\Sigma}(\gamma_{l})/\eta(\gamma_{l})\right\\}$$\triangleright$
Cluster centroid

19: for $(\bm{N}_{s},\bm{x}_{s})\in\mathcal{X}(\gamma_{l})$ do

20: Define cost function $c_{s}(\bm{N}_{s},\bm{x}_{s})\leftarrow
c^{\Sigma}(\gamma_{l})/\eta(\gamma_{l})+\Gamma_{s}(\bm{x}_{s})$

21: end for

22: end for

23: end for

24:end procedure

###  8.3 Proof of Proposition 4.3

Let $\gamma_{l}^{(q)}$ denote the $l^{\text{th}}$ cluster after $q$ states
have been assigned to it. Let us denote by $\Phi_{l}$ the size of the
$l^{\text{th}}$ cluster upon termination of the algorithm. Thus, the cluster
at termination is denoted by $\gamma_{l}^{{\Phi_{l}}}$.

We proceed by contradiction. Assume that there exists a cluster $l$ such that:

$\|\overline{\bm{N}}(\gamma_{l}^{(\Phi_{l})})-\underline{\bm{N}}(\gamma_{l}^{(\Phi_{l})})\|_{\infty}>\varepsilon$

Note that, by construction (Step 2), each cluster gets expanded at each
iteration. Therefore, the quantity
$\|\overline{\bm{N}}(\gamma^{(i)}_{l})-\underline{\bm{N}}(\gamma^{(i)}_{l})\|_{\infty}$
is non decreasing in $i$. Thus, there must exist $i\in\\{1,\cdots,\Phi_{l}$}
such that:

$\displaystyle\|\overline{\bm{N}}(\gamma^{(i)}_{l})-\underline{\bm{N}}(\gamma^{(i)}_{l})\|_{\infty}\leq\varepsilon$
$\displaystyle\|\overline{\bm{N}}(\gamma_{l}^{(i+1)})-\underline{\bm{N}}(\gamma_{l}^{(i+1)})\|_{\infty}>\varepsilon$

By definition of the $\ell_{\infty}$ norm, we know that there exists a
coordinate $j=1,\cdots,r$ such that:

$\displaystyle|(\overline{\bm{N}}(\gamma^{(i)}_{l}))_{j}-(\underline{\bm{N}}(\gamma^{(i)}_{l}))_{j}|\leq\varepsilon$
$\displaystyle|(\overline{\bm{N}}(\gamma_{l}^{(i+1)}))_{j}-(\underline{\bm{N}}(\gamma_{l}^{(i+1)}))_{j}|>\varepsilon.$

Denote the $i+1^{\text{th}}$ sample that got clustered into $\gamma_{l}$ as
$\bm{N}_{s+1}^{*}$. By construction (Step 2), at most one of
$(\overline{\bm{N}}(\gamma^{(i)}_{l}))_{j}$ and
$(\underline{\bm{N}}(\gamma^{(i)}_{l}))_{j}$ can change after
$\bm{N}_{s+1}^{*}$ gets added to the cluster. Without loss of generality,
assume that $\bm{N}_{s+1}^{*}$ changed the minimum state, that is:

$(\bm{N}_{s+1}^{*})_{j}=(\underline{\bm{N}}(\gamma_{l}^{(i+1)}))_{j}\neq(\underline{\bm{N}}(\gamma^{(i)}_{l}))_{j}\qquad(\overline{\bm{N}}(\gamma_{l}^{(i+1)}))_{j}=(\overline{\bm{N}}(\gamma^{(i)}_{l}))_{j}$

However, this implies:

$\displaystyle\max(\|\bm{N}_{s+1}^{*}-\underline{\bm{N}}(\gamma^{(i)}_{l})\|_{\infty},\|\bm{N}_{s+1}^{*}-\overline{\bm{N}}(\gamma^{(i)}_{l})\|_{\infty})$
$\displaystyle\geq\|\bm{N}_{s+1}^{*}-\overline{\bm{N}}(\gamma^{(i)}_{l})\|_{\infty}$
$\displaystyle\geq|(\bm{N}_{s+1}^{*})_{j}-(\overline{\bm{N}}(\gamma^{(i)}_{l}))_{j}|$
$\displaystyle=|(\underline{\bm{N}}(\gamma_{l}^{(i+1)}))_{j}-(\overline{\bm{N}}(\gamma_{l}^{(i+1)}))_{j}|$
$\displaystyle>\varepsilon.$

This implies that $\bm{N}_{s+1}^{*}$ would not be accepted in the cluster
(Step 3), a contradiction. Therefore, we have:

$\|\overline{\bm{N}}(\gamma_{l}^{(\Phi_{l})})-\underline{\bm{N}}(\gamma_{l}^{(\Phi_{l})})\|_{\infty}\leq\varepsilon\;\;\forall
l$

This completes the proof.

###  8.4 Proof of Proposition 4.4

By definition of Algorithm 3, the clustered and true states coincide at
decision epoch $s=1$. By Proposition 4.3, the error at decision epoch $s=2$ is
less than $\varepsilon$. We then proceed by induction, assuming that the
result is true for some decision epoch $s$ and proving it for epoch $s+1$.

Define a true state $\bm{N}^{*}_{s+1}\in\mathcal{N}^{*}_{s+1}$. By definition
of $\mathcal{N}^{*}_{s+1}$, there exist a true state
$\bm{N}^{*}_{s}\in\mathcal{N}^{*}_{s}$ and a decision
$\bm{x}_{s}\in\mathcal{F}_{s}$ at the preceding epoch that leads to
$\bm{N}^{*}_{s+1}$, that is:

$\bm{M}^{*}(\tau_{s})\leftarrow\bm{N}^{*}_{s}\ \text{and}\
\frac{d\bm{M}^{*}(t)}{dt}=f(\bm{M}^{*}(t),\bm{x}_{s}),\ \forall
t\in\left[\tau_{s},\tau_{s+1}\right]\implies\bm{M}^{*}(\tau_{s+1})=\bm{N}^{*}_{s+1}.$

Per the induction hypothesis, there exists a clustered state
$\bm{N}_{s}\in\mathcal{N}_{s}$ such that:

$\|\bm{N}^{*}_{s}-\bm{N}_{s}\|_{\infty}\leq\varepsilon\sum_{\sigma=3}^{s+1}\exp\left(\sum_{\nu=\sigma}^{s}L_{\nu-1}(\tau_{\nu}-\tau_{\nu-1})\right)$
(146)

In order to capture the propagation of approximation error between decision
epochs $s$ and $s+1$, we make use of Gronwall’s inequality. Let us specify it
in Lemma 8.1.

###### Lemma 8.1

Let $L$ be a constant and $\bm{u}:[a,b]\mapsto\mathbb{R}^{n}$ be a real-valued
$n$-dimensional continuous function defined over an interval $[a,b]$. Assume
that $u$ is differentiable on the interior of $(a,b)$ and satisfies
$\|\bm{u}^{\prime}(t)\|_{\infty}\leq L\|\bm{u}(t)\|_{\infty}$ for all
$t\in(a,b)$. Then, for all $t\leq b$,
$\|\bm{u}(t)\|_{\infty}\leq\|\bm{u}(a)\|_{\infty}e^{L(t-a)}$.

Let us consider the state $\widehat{\bm{N}}_{s+1}$ obtained at decision epoch
$s+1$ starting from the clustered state $\bm{N}_{s}\in\mathcal{N}_{s}$ at
epoch $s$ with decision $\bm{x}_{s}\in\mathcal{F}$:

$\bm{M}(\tau_{s})\leftarrow\bm{N}_{s}\ \text{and}\
\frac{d\bm{M}(t)}{dt}=f(\bm{M}(t),\bm{x}_{s}),\ \forall
t\in\left[\tau_{s},\tau_{s+1}\right]\ \text{and}\
\widehat{\bm{N}}_{s+1}\leftarrow\bm{M}(\tau_{s+1}).$

By assumption, the transition function is $L_{s}$-Lipschitz over
$[\tau_{s},\tau_{s+1}]$. Therefore, we have:

$\|f(\bm{M}^{*}(t),\bm{x}_{s})-f(\bm{M}(t),\bm{x}_{s})\|_{\infty}\leq
L_{s}\|\bm{M}^{*}(t)-\bm{M}(t)\|_{\infty},\quad\forall
t\in[\tau_{s},\tau_{s+1}]$

Per Lemma 8.1, we obtain:

$\displaystyle\|\bm{N}^{*}_{s+1}-\widehat{\bm{N}}_{s+1}\|_{\infty}$
$\displaystyle=\|\bm{M}^{*}(\tau_{s+1})-\bm{M}(\tau_{s+1})\|_{\infty}$
$\displaystyle\leq
e^{L_{s}(\tau_{s+1}-\tau_{s})}\|\bm{M}^{*}(\tau_{s})-\bm{M}(\tau_{s})\|_{\infty}$
$\displaystyle=e^{L_{s}(\tau_{s+1}-\tau_{s})}\|\bm{N}^{*}_{s}-\bm{N}_{s}\|_{\infty}$
(147)

Combining Equations (146) and (147), we obtain:

$\displaystyle\|\bm{N}^{*}_{s+1}-\widehat{\bm{N}}_{s+1}\|_{\infty}$
$\displaystyle\leq\varepsilon\sum_{\sigma=3}^{s+1}\exp\left(\sum_{\nu=\sigma}^{s}L_{\nu-1}(\tau_{\nu}-\tau_{\nu-1})\right)\cdot
e^{L_{s}(\tau_{s+1}-\tau_{s})}$
$\displaystyle=\varepsilon\sum_{\sigma=3}^{s+1}\exp\left(\sum_{\nu=\sigma}^{s+1}L_{\nu-1}(\tau_{\nu}-\tau_{\nu-1})\right)$
(148)

Now, by construction of the state-clustering algorithm (Algorithm 3), we
assign state $\widehat{\bm{N}}_{s+1}$ to a cluster with centroid
$\bm{N}_{s+1}\in\mathcal{N}_{s+1}$ at decision epoch $s+1$. Per Proposition
4.3, we therefore know that:

$\|\bm{N}_{s+1}-\widehat{\bm{N}}_{s+1}\|_{\infty}\leq\varepsilon$ (149)

Combining Equations (148) and (149), we obtain, per the triangular inequality:

$\displaystyle\|\bm{N}^{*}_{s+1}-\bm{N}_{s+1}\|_{\infty}$
$\displaystyle\leq\|\bm{N}^{*}_{s+1}-\widehat{\bm{N}}_{s+1}\|_{\infty}+\|\widehat{\bm{N}}_{s+1}-\bm{N}_{s+1}\|_{\infty}$
$\displaystyle\leq\varepsilon\left(1+\sum_{\sigma=3}^{s+1}\exp\left(\sum_{\nu=\sigma}^{s+1}L_{\nu-1}(\tau_{\nu}-\tau_{\nu-1})\right)\right)$
$\displaystyle=\sum_{\sigma=3}^{s+2}\exp\left(\sum_{\nu=\sigma}^{s+1}L_{\nu-1}(\tau_{\nu}-\tau_{\nu-1})\right),$

where we use the convention that an empty sum is equal to zero.

This completes the induction.

###  8.5 Proof of Theorem 4.5

Consider a node of the branching tree, and let $\widehat{\bm{z}}$ denote the
solution of the linear relaxation of the set partitioning formulation along
with the corresponding branching constraints. We distinguish three cases:

  1. 1.

One of the variables $y_{l}$ is fractional for $l\in\\{1,\cdots,q\\}$. Then.
the incumbent solution is eliminated by the bi-partite branching disjunction
on $y_{l}$ (Step 2 of Algorithm 1).

  2. 2.

The variable $\bm{y}$ satisfies all integrality constraints but one of the
resource allocation constraints
$\bm{x}_{is}(\widehat{\bm{z}})\notin\mathcal{F}_{is}$ for some segment
$i\in\\{1,\cdots,n\\}$ and some decision epoch $s\in\\{1,\cdots,S\\}$. Then
there must exist a component $k$ such that
$\left\lfloor\bm{x}_{is}(\widehat{\bm{z}})\right\rfloor^{k}_{is}<\bm{x}_{is}(\widehat{\bm{z}})<\left\lceil\bm{x}_{is}(\widehat{\bm{z}})\right\rceil^{k}_{is}$.
The incumbent solution is therefore eliminated by the bi-partite branching
disjunction on $x^{k}_{is}(\bm{z})$ (Step 3 of Algorithm 1).

  3. 3.

The variable $\bm{y}$ satisfies all integrality constraints,
$\bm{x}_{is}(\widehat{\bm{z}})\in\mathcal{F}_{is}$ for all segments
$i\in\\{1,\cdots,n\\}$ and all decision epochs $s\in\\{1,\cdots,S\\}$, but
there exists a fractional plan-based variable $z^{p}_{i}\in(0,1)$. As noted in
Section 4.4, 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.$

Then, the tri-partite branching disjunction (Equation (39), Step 4 of
Algorithm 1) makes the incumbent solution infeasible. To see this, let
$\widehat{\mathcal{P}}_{i}=\\{p\in\mathcal{P}_{i}:\widehat{z}^{p}_{i}>0\\}$
denote the set of plans included in the incumbent solution for segment $i$,
partitioned into
$\widehat{\mathcal{P}}^{<}_{i}=\\{p\in\widehat{\mathcal{P}}_{i}:\alpha^{pk}_{is}<x^{k}_{is}(\widehat{\bm{z}})\\}$,
$\widehat{\mathcal{P}}^{=}_{i}=\\{p\in\widehat{\mathcal{P}}_{i}:\alpha^{pk}_{is}=x^{k}_{is}(\widehat{\bm{z}})\\}$,
and
$\widehat{\mathcal{P}}^{>}_{i}=\\{p\in\widehat{\mathcal{P}}_{i}:\alpha^{pk}_{is}>x^{k}_{is}(\widehat{\bm{z}})\\}$.
Since $\alpha^{k,p_{0}}_{is}\neq x^{k}_{is}(\widehat{\bm{z}})$, we know that
$\widehat{\mathcal{P}}^{<}_{i}\neq\emptyset$ or
$\widehat{\mathcal{P}}^{>}_{i}\neq\emptyset$. But since, in addition,
$\sum_{p\in\mathcal{P}_{i}}\alpha^{pk}_{is}\widehat{z}^{p}_{i}=x^{k}_{is}$, we
actually know that $\widehat{\mathcal{P}}^{<}_{i}\neq\emptyset$ and
$\widehat{\mathcal{P}}^{>}_{i}\neq\emptyset$. This implies that the incumbent
solution is infeasible in all three subsequent branches. Indeed:

     * –

In the left branch, we eliminate the plans in $\widehat{\mathcal{P}}^{>}_{i}$
and in $\widehat{\mathcal{P}}^{=}_{i}$, so
$\alpha^{kp}_{is}<x^{k}_{is}(\widehat{\bm{z}})$ for all $p\in\mathcal{P}$ such
that $z_{i}^{p}>0$. Therefore, there exists no linear combination of remaining
plan-based variables $\bm{z}$ such that
$x^{k}_{is}(\bm{z})=x^{k}_{is}(\widehat{\bm{z}})$.

     * –

In the right branch, the solution is infeasible for similar reasons.

     * –

In the middle branch, we eliminate the plans in
$\widehat{\mathcal{P}}^{<}_{i}$ and in $\widehat{\mathcal{P}}^{>}_{i}$, so the
branching constraints force $z_{i}^{p_{0}}=0$. Therefore, the incumbent
solution is infeasible in the $\bm{z}$ space.

This concludes that the branching scheme eliminates any infeasible incumbent
solution. Note that the tri-partite branching proceeds differently in the left
and right branches, on the one hand, and in the middle branch, on the other
hand. Indeed, the left branch and the right branch force new resource
allocation decisions. In contrast, the middle branch may result in the same
resource allocation decisions but nonetheless restores integrality of
$z_{i}^{p_{0}}$ variable and updates the cost estimates accordingly—thus
contributing to tightening the relaxation bound and to progressing toward a
feasible solution.

Finally, the branching scheme ensures that the solution is feasible upon
termination. Recall that the algorithm terminates when
$\sum_{p\in\mathcal{P}_{i}}\widehat{z}^{p}_{i}|\alpha^{pk}_{is}-x^{k}_{is}(\widehat{\bm{z}})|=0$
for all segments $i\in\\{1,\cdots,n\\}$, all decision epochs
$s\in\\{1,\cdots,S\\}$, and all components $k$. Then, we know that
$\alpha^{pk}_{is}=x^{k}_{is}(\widehat{\bm{z}})$ whenever
$\widehat{z}_{i}^{p}>0$. This ensures that a single plan $p\in\mathcal{P}_{i}$
is such that $\widehat{z}_{i}^{p}>0$ for all segments $i\in\\{1,\cdots,n\\}$,
and the solution is therefore integral. Finiteness of the algorithm follows
from the discreteness of the solution space.

##  9 Additional computational results

###  9.1 Details on the experimental setup for the vaccination centers problem

Recall that, for vaccine allocation problem, we apply our methodology to the
full United States, comprising 51 regions (50 states plus Washington, DC),
over a three-month period. The vaccination centers problem, however, is too
complex to be solved at that scale. Instead, we solve it for smaller groups of
states, defined by the CDC. Table 8 reports the composition of each of the 10
groups, along with their share of the US population, the budget of vaccines in
each group, and the resulting number of discretized decisions (assuming
pallets of 25,000 vaccines). Note that we consider here an overall budget of
2.5 million vaccines per week, which we distribution across the ten groups
proportionally to the US population.

Table 8: Group composition for the vaccination centers problem. Group | States | % US Population | Vaccines | Decisions
---|---|---|---|---
| Connecticut, Maine, Massachusetts, |  |  |
A | New Hampshire, Rhode Island, Vermont, | 11.8% | 300,000 | 12
| New York |  |  |
| Delaware, District of Columbia, Maryland, |  |  |
B | Pennsylvania, Virginia, West Virginia, | 9.8% | 250,000 | 10
| New Jersey |  |  |
C | North Carolina, South Carolina, | 12.7% | 325,000 | 13
| Georgia, Florida |  |  |
D | Kentucky, Tennessee, | 4.9% | 125,000 | 5
| Alabama, Mississippi |  |  |
E | Illinois, Indiana, Michigan, | 15.6% | 400,000 | 16
| Minnesota, Ohio, Wisconsin |  |  |
F | Arkansas, Louisiana, New Mexico, | 12.7% | 325,000 | 13
| Oklahoma, Texas |  |  |
G | Iowa, Kansas, | 3.9% | 100,000 | 4
| Missouri, Nebraska |  |  |
H | Colorado, Montana, North Dakota, | 3.9% | 100,000 | 4
| South Dakota, Utah, Wyoming |  |  |
I | Arizona, California, | 19.6% | 500,000 | 20
| Hawaii, Nevada |  |  |
J | Alaska, Idaho, | 4.9% | 125,000 | 5
| Oregon, Washington |  |  |

###  9.2 Additional results on the branch-and-price algorithm

Tables 9–12 show extended results for the four problems, complementing the
results shown in Table 2 of the paper. The observations are consistent with
those from Section 5.2, thus showing the robustness of the findings across
problem instances.

Table 9: Detailed performance assessment for the vaccine allocation problem. |  |  |  |  | CPU times (s) | Solution quality
---|---|---|---|---|---|---
$n$ | $S$ | $D$ | Method | Nodes | Init. | RMP | PP | Total | Upper bound | Solution | Gap
51 | 6 | 6 | Column Generation | 1 | 15 | 0.02 | 0.14 | 0.2 | 45,995 | 45,563 | 0.94%
|  |  | Bi-partite B&P | 14 | 15 | 0.11 | 0.69 | 5.25 | 45,821 | 45,780 | 0.09%
|  |  | Tri-partite B&P | 14 | 15 | 0.12 | 0.59 | 5.61 | 45,821 | 45,780 | 0.09%
51 | 6 | 11 | Column Generation | 1 | 27 | 0.02 | 0.55 | 0.61 | 48,240 | 48,108 | 0.27%
|  |  | Bi-partite B&P | 14 | 27 | 0.14 | 1.24 | 5.94 | 48,208 | 48,208 | 0%
|  |  | Tri-partite B&P | 14 | 27 | 0.17 | 1.5 | 7.01 | 48,208 | 48,208 | 0%
51 | 6 | 21 | Column Generation | 1 | 58 | 0.01 | 0.75 | 0.8 | 50,649 | 50,143 | 1%
|  |  | Bi-partite B&P | 9,990 | 57 | 52 | 710 | 1,249 | 50,542 | 50,437 | 0.21%
|  |  | Tri-partite B&P | 9,990 | 57 | 51 | 709 | 1,250 | 50,542 | 50,437 | 0.21%
51 | 8 | 6 | Column Generation | 1 | 21 | 0.02 | 0.22 | 0.29 | 93,694 | 92,261 | 1.53%
|  |  | Bi-partite B&P | 28 | 21 | 0.27 | 1.43 | 7.11 | 93,378 | 93,263 | 0.12%
|  |  | Tri-partite B&P | 46 | 21 | 0.42 | 2.22 | 9.62 | 93,353 | 93,263 | 0.09%
51 | 8 | 11 | Column Generation | 1 | 34 | 0.02 | 0.58 | 0.68 | 99,296 | 98,801 | 0.5%
|  |  | Bi-partite B&P | 732 | 34 | 7.8 | 52 | 137 | 99,081 | 98,969 | 0.11%
|  |  | Tri-partite B&P | 744 | 34 | 7.4 | 53 | 139 | 99,068 | 98,969 | 0.1%
51 | 8 | 21 | Column Generation | 1 | 87 | 0.02 | 2.04 | 2.13 | 103,992 | 103,668 | 0.31%
|  |  | Bi-partite B&P | 1,868 | 87 | 20 | 346 | 557 | 103,771 | 103,668 | 0.1%
|  |  | Tri-partite B&P | 1,868 | 87 | 20 | 346 | 555 | 103,771 | 103,668 | 0.1%
51 | 10 | 6 | Column Generation | 1 | 23 | 0.02 | 0.31 | 0.38 | 146,281 | 145,698 | 0.4%
|  |  | Bi-partite B&P | 16 | 23 | 0.17 | 1.39 | 6.22 | 145,947 | 145,810 | 0.09%
|  |  | Tri-partite B&P | 16 | 23 | 0.18 | 1.4 | 6.32 | 145,947 | 145,810 | 0.09%
51 | 10 | 11 | Column Generation | 1 | 43 | 0.02 | 0.81 | 0.92 | 156,464 | 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 | 1058 | 156,183 | 156,026 | 0.1%
51 | 10 | 21 | Column Generation | 1 | 118 | 0.02 | 2.55 | 2.78 | 164,927 | 163,111 | 1.1%
|  |  | Bi-partite B&P | 9,990 | 118 | 62 | 1,262 | 2,495 | 164,698 | 164,065 | 0.38%
|  |  | Tri-partite B&P | 9,990 | 118 | 62 | 1,251 | 2,486 | 164,698 | 164,065 | 0.38%
51 | 12 | 6 | Column Generation | 1 | 27 | 0.03 | 0.5 | 0.6 | 198,566 | 196,855 | 0.86%
|  |  | Bi-partite B&P | 1,224 | 27 | 17 | 63 | 238 | 197,468 | 197,272 | 0.01%
|  |  | Tri-partite B&P | 1,302 | 27 | 19 | 68 | 257 | 197,468 | 197,272 | 0.01%
51 | 12 | 11 | Column Generation | 1 | 51 | 0.03 | 0.95 | 1.04 | 212,653 | 211,412 | 0.58%
|  |  | Bi-partite B&P | 122 | 51 | 1.37 | 18 | 35 | 212,027 | 211,978 | 0.02%
|  |  | Tri-partite B&P | 122 | 51 | 1.33 | 18 | 35 | 212,027 | 211,978 | 0.02%
51 | 12 | 21 | Column Generation | 1 | 144 | 0.03 | 4.1 | 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.25%
|  |  | Tri-partite B&P | 9,990 | 144 | 103 | 1,808 | 3,535 | 225,048 | 224,495 | 0.25%
Table 10: Detailed performance assessment for vaccination centers problem.

|  |  |  |  |  |  |  | CPU times (s) | Solution quality
---|---|---|---|---|---|---|---|---|---
Group | $n$ | $S$ | $D$ | $K$ | Flex.? | Method | Nodes | Init. | RMP | PP | Total | Upper bound | Solution | Gap
A | 7 | 6 | 13 | 2 | No | Column Generation | 1 | 58 | 0.12 | 0.33 | 0.96 | 7,224 | 6,109 | 18.66%
|  |  |  |  |  | Tri-partite B&P | 36,017 | 58 | 794 | 440 | 10,813 | 7,094 | 6,665 | 6.05%
A | 7 | 6 | 13 | 2 | Yes | Column Generation | 1 | 45 | 0.19 | 0.43 | 1.11 | 8200 | 6660 | 29.91%
|  |  |  |  |  | Tri-partite B&P | 11,899 | 45 | 224 | 178 | 10,808 | 8,132 | 7,917 | 2.64%
A | 7 | 6 | 13 | 3 | No | Column Generation | 1 | 46 | 0.18 | 0.27 | 0.86 | 8,200 | 6,703 | 39.32%
|  |  |  |  |  | Tri-partite B&P | 16,586 | 45 | 383 | 241 | 10,839 | 7,867 | 7,511 | 4.53%
A | 7 | 6 | 13 | 3 | Yes | 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%
A | 7 | 8 | 13 | 2 | No | Column Generation | 1 | 79 | 0.30 | 1.34 | 2.17 | 16,562 | 12,471 | 48.54%
|  |  |  |  |  | Tri-partite B&P | 24,453 | 78 | 1,138 | 701 | 10,827 | 16,381 | 15,842 | 3.29%
A | 7 | 8 | 13 | 2 | Yes | Column Generation | 1 | 54 | 0.33 | 1.01 | 2.02 | 18,662 | 16,417 | 17.69%
|  |  |  |  |  | Tri-partite B&P | 33,034 | 53 | 979 | 754 | 10,803 | 18,628 | 18,443 | 0.99%
A | 7 | 8 | 13 | 3 | No | Column Generation | 1 | 54 | 0.41 | 1.02 | 1.87 | 18,652 | 14,619 | 27.66%
|  |  |  |  |  | Tri-partite B&P | 8,638 | 53 | 289 | 305 | 10,891 | 17,740 | 16,929 | 4.57%
A | 7 | 8 | 13 | 3 | Yes | Column Generation | 1 | 53 | 0.45 | 1.05 | 1.93 | 18,662 | 14,918 | 31.59%
|  |  |  |  |  | Tri-partite B&P | 12,287 | 53 | 343 | 375 | 10,939 | 18,274 | 17,187 | 5.94%
H | 6 | 6 | 5 | 2 | No | Column Generation | 1 | 49 | 0.03 | 0.04 | 0.55 | 786 | 652 | 16.97%
|  |  |  |  |  | Tri-partite B&P | 4,678 | 49 | 147 | 24 | 1,026 | 764 | 763 | 0.09%
H | 6 | 6 | 5 | 2 | Yes | Column Generation | 1 | 36 | 0.05 | 0.03 | 0.46 | 786 | 652 | 29.15%
|  |  |  |  |  | Tri-partite B&P | 2,919 | 36 | 75 | 11 | 604 | 764 | 763 | 0.09%
H | 6 | 6 | 5 | 3 | No | Column Generation | 1 | 42 | 0.03 | 0.02 | 0.48 | 678 | 652 | 3.75%
|  |  |  |  |  | Tri-partite B&P | 45,144 | 42 | 1,278 | 113 | 10,802 | 675 | 674 | 0.23%
H | 6 | 6 | 5 | 3 | Yes | 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%
H | 6 | 8 | 5 | 2 | No | Column Generation | 1 | 45 | 0.04 | 0.08 | 0.68 | 1,978 | 1,671 | 15.55%
|  |  |  |  |  | Tri-partite B&P | 62 | 44 | 1.43 | 0.89 | 6.99 | 1,978 | 1,978 | 0.00%
H | 6 | 8 | 5 | 2 | Yes | Column Generation | 1 | 42 | 0.08 | 0.07 | 0.58 | 1,978 | 1,753 | 49.08%
|  |  |  |  |  | Tri-partite B&P | 20,280 | 41 | 2,002 | 76 | 10,818 | 1,886 | 1,865 | 1.11%
H | 6 | 8 | 5 | 3 | No | Column Generation | 1 | 43 | 0.05 | 0.05 | 0.53 | 1791 | 1780 | 0.60%
|  |  |  |  |  | Tri-partite B&P | 14,149 | 42 | 298 | 71 | 1,835 | 1,787 | 1,787 | 0.05%
H | 6 | 8 | 5 | 3 | Yes | Column Generation | 1 | 45 | 0.03 | 0.05 | 0.54 | 1,791 | 1,780 | 0.60%
|  |  |  |  |  | Tri-partite B&P | 11,927 | 45 | 231 | 63 | 1,506 | 1,787 | 1,787 | 0.05%

  * •

“Flex.” refers to flexibility in vaccine allocations across selected
facilities. In the no-flexibility case, the capacity of each facility is equal
to the vaccine budget (shown in Table 8) divided by the number of facilities
$K$. In the flexibility case, an extra pallet is available in each facility
(Equation (112)), with the same overall budget (Equation (111)).

Table 11: Detailed performance assessment for the content promotion problem. |  |  |  |  |  | CPU times (s) | Solution quality
---|---|---|---|---|---|---|---
$n$ | $S$ | $D$ | $K$ | Method | Nodes | Init. | RMP | PP | Total | Upper bound | Solution | Gap
20 | 6 | 11 | 4 | Column Generation | 1 | 49 | 0.03 | 2.17 | 2.35 | 1.015 | 1.004 | 1.06%
|  |  |  | Bi-Partite B&P | 8 | 49 | 0.08 | 4.04 | 5.9 | 1.015 | 1.015 | 0%
|  |  |  | Tri-partite B&P | 8 | 49 | 0.09 | 4.11 | 6.1 | 1.015 | 1.015 | 0%
20 | 6 | 11 | 6 | Column Generation | 1 | 56 | 0.02 | 1.76 | 2.01 | 1.017 | 1.007 | 0.96%
|  |  |  | Bi-Partite B&P | 25,958 | 56 | 417 | 2,430 | 7,200 | 1.017 | 1.016 | 0.1%
|  |  |  | Tri-partite B&P | 26,194 | 56 | 426 | 2,407 | 7,201 | 1.017 | 1.016 | 0.1%
20 | 6 | 21 | 2 | Column Generation | 1 | 79 | 0.04 | 4.62 | 4.84 | 1.017 | 0.978 | 3.9%
|  |  |  | Bi-Partite B&P | 696 | 79 | 9.35 | 198 | 327 | 1.015 | 1.009 | 0.53%
|  |  |  | Tri-partite B&P | 142 | 79 | 1.67 | 52 | 74 | 1.013 | 1.013 | 0.09%
20 | 6 | 21 | 4 | Column Generation | 1 | 72 | 0.03 | 3.77 | 4.06 | 1.052 | 1.04 | 1.13%
|  |  |  | Bi-Partite B&P | 11,400 | 72 | 179 | 1,861 | 7,201 | 1.051 | 1.047 | 0.35%
|  |  |  | Tri-partite B&P | 5,908 | 72 | 106 | 1,158 | 3,913 | 1.049 | 1.049 | 0.06%
20 | 6 | 21 | 6 | Column Generation | 1 | 61 | 0.03 | 4.85 | 5.08 | 1.069 | 1.056 | 1.29%
|  |  |  | Bi-Partite B&P | 19,660 | 61 | 223 | 2,761 | 7,202 | 1.067 | 1.065 | 0.17%
|  |  |  | Tri-partite B&P | 19,640 | 61 | 222 | 2,766 | 7,201 | 1.067 | 1.065 | 0.17%
20 | 8 | 11 | 2 | Column Generation | 1 | 58 | 0.04 | 5.01 | 5.37 | 1.527 | 1.44 | 5.7%
|  |  |  | Bi-Partite B&P | 17,960 | 58 | 294 | 1,949 | 7,201 | 1.523 | 1.511 | 0.82%
|  |  |  | Tri-partite B&P | 17,819 | 58 | 338 | 2,057 | 7,201 | 1.52 | 1.515 | 0.31%
20 | 8 | 11 | 4 | Column Generation | 1 | 57 | 0.03 | 3.80 | 4.12 | 1.545 | 1.5 | 2.93%
|  |  |  | Bi-Partite B&P | 18,876 | 57 | 158 | 1,845 | 10,801 | 1.545 | 1.533 | 0.74%
|  |  |  | Tri-partite B&P | 19,016 | 57 | 156 | 1,783 | 10,801 | 1.545 | 1.53 | 0.74%
20 | 8 | 11 | 6 | Column Generation | 1 | 79 | 0.05 | 5.64 | 6.05 | 1.551 | 1.529 | 1.41%
|  |  |  | Bi-Partite B&P | 21,122 | 79 | 187 | 2,686 | 10,802 | 1.55 | 1.542 | 0.51%
|  |  |  | Tri-partite B&P | 23,670 | 79 | 196 | 2,588 | 10,802 | 1.55 | 1.542 | 0.51%
20 | 8 | 21 | 4 | Column Generation | 1 | 103 | 0.08 | 12 | 13 | 1.575 | 1.528 | 3.01%
|  |  |  | Bi-Partite B&P | 6,546 | 103 | 77 | 1,642 | 10,802 | 1.572 | 1.564 | 0.45%
|  |  |  | Tri-partite B&P | 6,536 | 103 | 74 | 1,685 | 10,805 | 1.572 | 1.564 | 0.45%
20 | 8 | 21 | 6 | Column Generation | 1 | 79 | 0.05 | 8.61 | 10 | 1.602 | 1.587 | 0.91%
|  |  |  | Bi-Partite B&P | 5,462 | 79 | 58 | 1,156 | 10,805 | 1.598 | 1.595 | 0.23%
|  |  |  | Tri-partite B&P | 5,252 | 79 | 59 | 1,423 | 10,806 | 1.598 | 1.595 | 0.23%
20 | 10 | 11 | 2 | Column Generation | 1 | 78 | 0.08 | 8.64 | 9.08 | 2.078 | 1.977 | 4.82%
|  |  |  | Bi-Partite B&P | 16,446 | 78 | 369 | 3,337 | 10,803 | 2.07 | 2.045 | 1.22%
|  |  |  | Tri-partite B&P | 15,688 | 78 | 347 | 3,627 | 10,803 | 2.07 | 2.045 | 1.22%
20 | 10 | 11 | 4 | Column Generation | 1 | 82 | 0.06 | 6.85 | 7.24 | 2.098 | 2.070 | 1.31%
|  |  |  | Bi-Partite B&P | 11,450 | 82 | 131 | 1,521 | 10,802 | 2.098 | 2.09 | 0.35%
|  |  |  | Tri-partite B&P | 11,304 | 82 | 130 | 1,570 | 10,802 | 2.097 | 2.09 | 0.35%
20 | 10 | 11 | 6 | Column Generation | 1 | 90 | 0.08 | 7.87 | 8.65 | 2.105 | 2.088 | 0.81%
|  |  |  | Bi-Partite B&P | 10,394 | 90 | 122 | 1,435 | 10,802 | 2.105 | 2.097 | 0.41%
|  |  |  | Tri-partite B&P | 10,536 | 90 | 120 | 1,384 | 10,802 | 2.105 | 2.097 | 0.41%
20 | 10 | 21 | 2 | Column Generation | 1 | 134 | 0.08 | 15 | 16 | 2.017 | 1.874 | 7.08%
|  |  |  | Bi-Partite B&P | 12,606 | 134 | 229 | 3,221 | 10,801 | 2.009 | 1.988 | 1.04%
|  |  |  | Tri-partite B&P | 11,864 | 134 | 229 | 3,279 | 10,801 | 2.009 | 1.988 | 1.04%
20 | 10 | 21 | 4 | Column Generation | 1 | 111 | 0.09 | 19 | 32 | 2.115 | 2.058 | 2.7%
|  |  |  | Bi-Partite B&P | 1,050 | 111 | 14 | 454 | 10,823 | 2.112 | 2.099 | 0.61%
|  |  |  | Tri-partite B&P | 1,042 | 111 | 14 | 527 | 10,827 | 2.112 | 2.099 | 0.61%
20 | 10 | 21 | 6 | Column Generation | 1 | 124 | 0.08 | 15 | 18 | 2.164 | 2.147 | 0.69%
|  |  |  | Bi-Partite B&P | 1,016 | 124 | 13 | 331 | 10,818 | 2.162 | 2.151 | 0.5%
|  |  |  | Tri-partite B&P | 1,016 | 124 | 14 | 347 | 10,821 | 2.162 | 2.151 | 0.5%
Table 12: Detailed performance assessment for the congestion mitigation problem. |  |  |  |  |  | CPU times (s) | Solution quality
---|---|---|---|---|---|---|---
$n$ | $S$ | $B^{1}$ | $B^{2}$ | Method | Nodes | Init. | RMP | PP | Total | Upper bound | Solution | Gap
5 | 2 | 4 | 2 | Column Generation | 1 | 7.54 | 0.05 | 0.02 | 0.1 | 0.50 | 0.50 | 0%
|  |  |  | Bi-partite B&P | 1 | 7.54 | 0.01 | 0.02 | 0.41 | 0.50 | 0.50 | 0%
|  |  |  | Tri-partite B&P | 1 | 7.54 | 0.01 | 0.02 | 0.41 | 0.50 | 0.50 | 0%
5 | 2 | 4 | 4 | Column Generation | 1 | 8.37 | 0.01 | 0.04 | 0.05 | 0.63 | 0.63 | 0%
|  |  |  | Bi-partite B&P | 1 | 8.37 | 0.03 | 0.04 | 0.48 | 0.63 | 0.63 | 0%
|  |  |  | Tri-partite B&P | 1 | 8.37 | 0.01 | 0.04 | 0.45 | 0.63 | 0.63 | 0%
5 | 2 | 6 | 2 | Column Generation | 1 | 8.04 | 0.01 | 0.03 | 0.06 | 0.68 | 0.67 | 0.6%
|  |  |  | Bi-partite B&P | 2 | 8.04 | 0.04 | 0.05 | 0.54 | 0.68 | 0.68 | 0.2%
|  |  |  | Tri-partite B&P | 5 | 8.04 | 0.02 | 0.07 | 0.8 | 0.68 | 0.68 | 0.02%
5 | 2 | 6 | 4 | Column Generation | 1 | 10 | 0.01 | 0.09 | 0.11 | 0.82 | 0.81 | 1.23%
|  |  |  | Bi-partite B&P | 2 | 10 | 0.02 | 0.13 | 0.56 | 0.82 | 0.81 | 0.16%
|  |  |  | Tri-partite B&P | 8 | 10s | 0.05 | 0.23 | 0.8s | 0.82 | 0.82 | 0%
5 | 2 | 8 | 2 | Column Generation | 1 | 8.85 | 0.01 | 0.26 | 0.28 | 0.87 | 0.87 | 0%
|  |  |  | Bi-partite B&P | 1 | 8.85 | 0.01 | 0.06 | 0.48 | 0.87 | 0.87 | 0%
|  |  |  | Tri-partite B&P | 1 | 8.85 | 0.01 | 0.06 | 0.51 | 0.87 | 0.87 | 0%
5 | 2 | 8 | 4 | Column Generation | 1 | 12 | 0.01 | 0.17 | 0.19 | 1.01 | 1.00 | 0.72%
|  |  |  | Bi-partite B&P | 8 | 12 | 0.06 | 0.36 | 0.95 | 1.01 | 1.01 | 0.27%
|  |  |  | Tri-partite B&P | 14 | 12 | 0.08 | 0.74 | 1.49 | 1.01 | 1.01 | 0.05%
5 | 4 | 4 | 2 | Column Generation | 1 | 14 | 0.02 | 0.41 | 0.44 | 1.34 | 1.34 | 0.31%
|  |  |  | Bi-partite B&P | 1 | 14 | 0.02 | 0.4 | 0.83 | 1.34 | 1.34 | 0.31%
|  |  |  | Tri-partite B&P | 6 | 14 | 0.05 | 0.92 | 1.74 | 1.34 | 1.34 | 0%
5 | 4 | 4 | 4 | Column Generation | 1 | 49 | 0.03 | 2.3 | 2.35 | 1.61 | 1.55 | 4.16%
|  |  |  | Bi-partite B&P | 28 | 49 | 0.2 | 8.7 | 10.8 | 1.61 | 1.59 | 1.06%
|  |  |  | Tri-partite B&P | 73 | 49 | 0.79 | 19.4 | 25 | 1.60 | 1.60 | 0.02%
5 | 4 | 6 | 2 | Column Generation | 1 | 33 | 0.02 | 1.35 | 1.39 | 1.90 | 1.90 | 0.19%
|  |  |  | Bi-partite B&P | 4 | 33 | 0.04 | 2.12 | 2.73 | 1.90 | 1.90 | 0.11%
|  |  |  | Tri-partite B&P | 7 | 33 | 0.08 | 3.38 | 4.07 | 1.90 | 1.90 | 0.04%
5 | 4 | 6 | 4 | Column Generation | 1 | 185 | 0.02 | 6.21 | 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%
5 | 4 | 8 | 2 | Column Generation | 1 | 83 | 0.03 | 4.02 | 4.1 | 2.47 | 2.44 | 1.11%
|  |  |  | Bi-partite B&P | 2 | 83 | 0.05 | 5.4 | 5.92 | 2.47 | 2.44 | 1.09%
|  |  |  | Tri-partite B&P | 63 | 83 | 0.63 | 29 | 33 | 2.46 | 2.46 | 0.07%
5 | 4 | 8 | 4 | Column Generation | 1 | 550 | 0.03 | 12 | 12 | 2.82 | 2.80 | 0.95%
|  |  |  | Bi-partite B&P | 4 | 550 | 0.049 | 18 | 18 | 2.82 | 2.80 | 0.89%
|  |  |  | Tri-partite B&P | 70 | 550 | 0.63 | 96 | 100 | 2.81 | 2.81 | 0.07%
5 | 6 | 4 | 2 | Column Generation | 1 | 24 | 0.03 | 2.02 | 2.09 | 2.05 | 1.94 | 5.32%
|  |  |  | Bi-partite B&P | 136 | 24 | 1.76 | 24 | 38 | 2.04 | 2.03 | 0.46%
|  |  |  | Tri-partite B&P | 203 | 24 | 2.9 | 32 | 56 | 2.03 | 2.03 | 0.002%
5 | 6 | 4 | 4 | Column Generation | 1 | 111 | 0.04 | 6.9 | 6.94 | 2.36 | 2.30 | 2.8%
|  |  |  | Bi-partite B&P | 724 | 111 | 20 | 328 | 656 | 2.36 | 2.33 | 1.02%
|  |  |  | Tri-partite B&P | 964 | 111 | 42 | 442 | 760 | 2.34 | 2.34 | 0.09%
5 | 6 | 6 | 2 | Column Generation | 1 | 74 | 0.04 | 5.9 | 6.07 | 2.94 | 2.80 | 4.72%
|  |  |  | Bi-partite B&P | 16 | 74 | 0.21 | 17 | 18 | 2.94 | 2.94 | 0.005%
|  |  |  | Tri-partite B&P | 16 | 74 | 1.06 | 18 | 20 | 2.94 | 2.94 | 0.005%
5 | 6 | 6 | 4 | Column Generation | 1 | 459 | 0.04 | 20 | 20 | 3.32 | 2.93 | 11.86%
|  |  |  | Bi-partite B&P | 206 | 459 | 4.1 | 350 | 389 | 3.31 | 3.31 | 0.19%
|  |  |  | Tri-partite B&P | 243 | 459 | 4.94 | 396 | 443 | 3.31 | 3.31 | 0.1%
5 | 6 | 8 | 2 | Column Generation | 1 | 209 | 0.04 | 12 | 12 | 3.93 | 3.80 | 3.39%
|  |  |  | Bi-partite B&P | 556 | 209 | 18 | 512 | 770 | 3.92 | 3.90 | 0.36%
|  |  |  | Tri-partite B&P | 757 | 209 | 31 | 676 | 1,088 | 3.91 | 3.90 | 0.1%
5 | 6 | 8 | 4 | Column Generation | 1 | 1,369 | 0.04 | 47 | 47 | 4.41 | 4.30 | 2.53%
|  |  |  | Bi-partite B&P | 56 | 1,369 | 0.8 | 273 | 278 | 4.40 | 4.40 | 0.1%
|  |  |  | Tri-partite B&P | 56 | 1,369 | 0.7 | 272 | 277 | 4.40 | 4.40 | 0.1%

###  9.3 Robustness to parameter estimation errors

Recall that we consider in this paper a deterministic prescriptive contagion
analytics model. We conduct in this appendix a sensitivity analysis to
establish the robustness of our results under parameter misspecification.
Specifically, we assess the performance of resource allocation decisions under
a perturbed contagion model in Table 13, for the vaccine allocation problem.
In these experiments, we consider all solutions obtained with each decision-
making procedure under the original DELPHI-V model (Figure 2a); we then
perturb all DELPHI-V parameters according to uniform distributions by up to
20% from their original values; and we evaluate the performance of each
solution on the perturbed model. We replicate this process 20 times and report
average performance, using the same nomenclature as in Table 4.

Results confirm the strong performance of our solution. Note that the number
of lives saved increases slightly as compared to those estimated in Table 4
with the nominal parameters, due to the non-linearities in the system
dynamics. However, the relative performance of the optimized vaccine
allocation solution is very similar to the one obtained under the nominal
parameters. Specifically, the optimized vaccine allocation solution yields
improvements of 42% (resp. 12%) with a supply of 2.5M vaccines (resp. 7M
vaccines) per week as compared to the cost-based allocation benchmark, an
improvements of 1.6% (resp. 5.3%) with a supply of 2.5M vaccines (resp. 7M
vaccines) per week as compared to the coordinate descent benchmark. These
results suggest that our solutions and our findings are highly robust against
model misspecification.

Table 13: Average death toll comparison under $20\%$ perturbation of contagion
parameters (vaccine allocation problem, full country, $D=21$, average over 20
samples).

|  | $S=6$ | $S=8$ | $S=12$
---|---|---|---|---
Budget | Method | Time (sec.) | Deaths | Time (sec.) | Deaths | Time (sec.) | Deaths
— | Do nothing | — | 576.57K | — | 607.48K | — | 657.02K
2.5M | Uniform allocation | — | -6.34K | — | -11.96K | — | -23.90K
| Cost-based allocation | — | -7.28K | — | -14.13K | — | -31.15K
| MIQO implementation | n/a | n/a | n/a | n/a | n/a | n/a
|  Discretization ($\delta=0.002$) |  1,000∗ | -5.47K | n/a | n/a | n/a | n/a
|  Discretization ($\delta=0.001$) |  1,000∗ | -5.81K | n/a | n/a | n/a | n/a
| Coordinate Descent | 5.43 | -12.90K | 4.73 | -21.75K | 6.03 | -43.49K
|  Branch-and-price ($\varepsilon=0.002$) | 426.2 | -12.98K | 1234.2 | -22.57K | 3671.4 | -44.17K
7M | Uniform allocation | — | -13.56K | — | -25.29K | — | -47.61K
| Cost-based allocation | — | -18.04K | — | -33.78K | — | -63.60K
| MIQO implementation | n/a | n/a | n/a | n/a | n/a | n/a
|  Discretization ($\delta=0.002$) |  1,000∗ | -6.07K | n/a | n/a | n/a | n/a
|  Discretization ($\delta=0.001$) |  1,000∗ | -3.18K | n/a | n/a | n/a | n/a
| Coordinate Descent | 5.46 | -21.25K | 12.7 | -38.46K | 118.5 | -67.43K
|  Branch-and-price ($\varepsilon=0.002$) | 67.5 | -22.03K | 156.8 | -39.73K | 708.7 | -70.99K

  * •

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

##  10 Benefits of $\ell_{\infty}$-based state clustering

In this appendix, we provide a detailed comparison of our state clustering
approach (referred to as $\ell_{\infty}$ clustering, for disambiguation) to a
$k$-means benchmark, detailed in Algorithm 4. Recall that the purpose of the
state clustering algorithm is to group $Q$ states
$\bm{N}_{s+1}^{1},\cdots,\bm{N}_{s+1}^{Q}$ into $k$ clusters at each epoch
$s=1,\cdots,S$. From a computational standpoint, our $\ell_{\infty}$ algorithm
scales linearly to retain tractability given the large number of states
(Algorithm 3), and so does $k$-means clustering. The key difference is that
our $\ell_{\infty}$ clustering algorithm assigns states based on their
$\ell_{\infty}$ distance to the element-wise minimum and maximum states in
each cluster; in turn, it provides guarantees on the diameter of each cluster
and on the global approximation error (Section 4.3). In comparison, the
$k$-means benchmark minimizes the average distance to cluster centroids. Note,
also, that the $k$-means benchmark requires a number of clusters as an input,
whereas our $\ell_{\infty}$ clustering algorithm assigns states to clusters
dynamically.

Algorithm 4 $k$-means state-clustering

1:procedure kMeansStateClustering($\mathcal{N}_{1}$, $k$)

2: for $s\leftarrow 1$ to $S$ do $\triangleright$ Forward pass

3: Initialization: $\mathcal{N}_{s+1}\leftarrow\emptyset$,
$\widetilde{\mathcal{N}}_{s+1}\leftarrow\emptyset$, $k\leftarrow 0$

4: for $\bm{N}_{s}\in\mathcal{N}_{s},\ \bm{x}_{s}\in\mathcal{F}_{s}$ do

5: Next state: $\bm{M}(\tau_{s})\leftarrow\bm{N}_{s}$;
$\frac{d\bm{M}(t)}{dt}=f(\bm{M}(t),\bm{x}_{s}),\ \forall
t\in\left[\tau_{s},\tau_{s+1}\right]$;
$\bm{N}_{s+1}\leftarrow\bm{M}(\tau_{s+1})$

6:
$\widetilde{\mathcal{N}}_{s+1}\leftarrow\widetilde{\mathcal{N}}_{s+1}\cup\\{\bm{N}_{s+1}\\}$

7: Compute cost:
$\overline{C}_{s}^{q}=\int_{\tau_{s}}^{\tau_{s+1}}g_{t}(\bm{M}^{q}(t))dt+\mathbf{1}\\{s=S\\}(h(\bm{M}^{q}(T)))$

8: end for

9:
$\gamma_{1},\cdots,\gamma_{k}\leftarrow\text{kMeans}(\widetilde{\mathcal{N}}_{s+1},k)$
$\triangleright$ $\gamma_{1},\cdots,\gamma_{k}$ denote the $k$-means clusters

10: Store centroids $\bm{\mu}(\gamma_{\ell})$, and binary parameters
$\xi^{q}(\gamma_{\ell})$ indicating if $\bm{N}_{s+1}^{q}$ is assigned to
cluster $\gamma_{\ell}$, $\ell=1,\cdots,k$

11: for $l\leftarrow 1$ to $k$ do

12: $\mathcal{N}_{s+1}\leftarrow\mathcal{N}_{s+1}\cup\\{\bm{\mu}(C_{\ell})\\}$

13: Compute average cost
$c_{s}^{\mu}(C_{\ell})\leftarrow\sum_{q=1}^{Q}\xi^{q}(\gamma_{l})\overline{C}^{q}_{s}$

14: end for

15: for $\bm{N}_{s}\in\mathcal{N}_{s},\ \bm{x}_{s}\in\mathcal{F}_{s}$ do

16: Define cost function
$c_{s}(\bm{N}_{s},\bm{x}_{s})\leftarrow\sum_{l=1}^{k}\xi^{q}(\gamma_{l})c_{s}^{\mu}(C_{\ell})+C_{s}(\bm{x}_{s})$

17: end for

18: end for

19:end procedure

First, we compare the $\ell_{\infty}$ clustering algorithm to the $k$-means
clustering in terms of approximation accuracy, for the vaccine allocation
problem. Table 14 reports the Median Absolute Error (MAE) and the Maximum
Absolute Error between the clustered states and the true states obtained with
exhaustive enumeration (Algorithm 2). Note that we consider here small-scale
examples for which exhaustive enumeration remains tractable. For a fair
comparison, we select values of $k$ to closely match the number of states from
$\ell_{\infty}$ clustering. As expected, $k$-means achieves a smaller median
absolute error but the $\ell_{\infty}$ algorithm achieves a smaller maximum
absolute error. For instance, in the six-period instance with the largest
number of clusters, both algorithms achieve comparable MAE but the $k$-means
benchmark results in a five-fold increase in the maximum absolute error; in
the four-period instance, the $k$-means benchmark achieves a very small MAE
but still increases the maximum absolute error by a factor 2 to 3. These
results reflect that our $\ell_{\infty}$ algorithm is precisely designed to
bound the diameter of its clusters (Proposition 4.3).

Table 14: $\ell_{\infty}$ clustering versus and $k$-means clustering (vaccine
allocation problem).

|  |  |  |  |  |  | Maximum Absolute Error ($\times 10^{-2}$) | Median Absolute Error ($\times 10^{-4}$)
---|---|---|---|---|---|---|---|---
$n$ | $S$ | $D$ | Method | Hyperparameter | $|S|$ | Time (sec.) | $s=2$ | $s=3$ | $s=4$ | $s=5$ | $s=6$ | $s=2$ | $s=3$ | $s=4$ | $s=5$ | $s=6$
51 | 4 | 6 |  $\ell_{\infty}$ Clust. | $\varepsilon=0.002$ | 2,350 | 3.33 | 0.12 | 0.19 | 0.38 | — | — | 0.00 | 0.29 | 0.27 | — | —
|  |  |  $k$-Means | $k=12$ | 2,193 | 3.41 | 0.00 | 0.51 | 1.49 | — | — | 0.00 | 0.03 | 0.02 | — | —
|  |  |  $\ell_{\infty}$ Clust. | $\varepsilon=0.005$ | 1,458 | 2.20 | 0.28 | 0.57 | 0.68 | — | — | 0.00 | 1.04 | 0.95 | — | —
|  |  |  $k$-Means | $k=8$ | 1,581 | 2.81 | 0.00 | 5.08 | 5.14 | — | — | 0.00 | 0.01 | 0.32 | — | —
|  |  |  $\ell_{\infty}$ Clust. | $\varepsilon=0.01$ | 982 | 1.57 | 0.68 | 0.87 | 0.86 | — | — | 0.01 | 3.72 | 4.02 | — | —
|  |  |  $k$-Means | $k=5$ | 1,071 | 1.95 | 4.96 | 9.22 | 14.15 | — | — | 0.00 | 0.95 | 0.97 | — | —
51 | 4 | 11 |  $\ell_{\infty}$ Clust. | $\varepsilon=0.002$ | 4,971 | 13.3 | 0.12 | 0.27 | 0.49 | — | — | 0.08 | 0.37 | 0.36 | — | —
|  |  |  $k$-Means | $k=25$ | 4,437 | 12.4 | 0.02 | 0.65 | 1.57 | — | — | 0.0 | 0.01 | 0.12 | — | —
|  |  |  $\ell_{\infty}$ Clust. | $\varepsilon=0.005$ | 2,533 | 6.74 | 0.27 | 0.68 | 0.83 | — | — | 0.29 | 1.59 | 1.42 | — | —
|  |  |  $k$-Means | $k=13$ | 2,601 | 9.34 | 0.01 | 4.47 | 4.38 | — | — | 0.03 | 0.08 | 0.24 | — | —
|  |  |  $\ell_{\infty}$ Clust. | $\varepsilon=0.01$ | 1,521 | 4.29 | 0.72 | 0.86 | 1.37 | — | — | 0.29 | 3.58 | 4.23 | — | —
|  |  |  $k$-Means | $k=8$ | 1,683 | 6.02 | 3.75 | 6.23 | 9.79 | — | — | 0.01 | 0.48 | 0.83 | — | —
51 | 4 | 21 |  $\ell_{\infty}$ Clust. | $\varepsilon=0.002$ | 10,615 | 74.6 | 0.08 | 0.27 | 0.67 | — | — | 0.08 | 0.52 | 0.41 | — | —
|  |  |  $k$-Means | $k=53$ | 9,231 | 59.3 | 0.00 | 1.14 | 1.52 | — | — | 0.00 | 0.00 | 0.01 | — | —
|  |  |  $\ell_{\infty}$ Clust. | $\varepsilon=0.005$ | 4,477 | 26.9 | 0.37 | 0.62 | 0.79 | — | — | 0.01 | 1.86 | 1.57 | — | —
|  |  |  $k$-Means | $k=22$ | 4,488 | 27.4 | 0.01 | 3.37 | 6.72 | — | — | 0.02 | 0.14 | 0.18 | — | —
|  |  |  $\ell_{\infty}$ Clust. | $\varepsilon=0.01$ | 2,473 | 14.8 | 0.82 | 0.89 | 1.47 | — | — | 0.28 | 3.78 | 4.91 | — | —
|  |  |  $k$-Means | $k=13$ | 2,703 | 18.5 | 3.82 | 6.76 | 10.47 | — | — | 0.00 | 0.51 | 0.62 | — | —
51 | 6 | 6 |  $\ell_{\infty}$ Clust. | $\varepsilon=0.002$ | 5,555 | 10.9 | 0.09 | 0.16 | 0.38 | 0.26 | 0.29 | 0.00 | 0.31 | 0.27 | 0.07 | 0.08
|  |  |  $k$-Means | $k=19$ | 5,202 | 9.11 | 0.00 | 0.32 | 1.58 | 1.66 | 1.79 | 0.00 | 0.00 | 0.02 | 0.07 | 0.08
|  |  |  $\ell_{\infty}$ Clust. | $\varepsilon=0.005$ | 2,903 | 5.51 | 0.32 | 0.63 | 0.71 | 0.84 | 0.87 | 0.00 | 1.02 | 0.88 | 0.49 | 0.16
|  |  |  $k$-Means | $k=10$ | 2,907 | 5.24 | 0.00 | 0.87 | 1.86 | 3.04 | 3.42 | 0.00 | 0.01 | 0.22 | 0.13 | 0.12
|  |  |  $\ell_{\infty}$ Clust. | $\varepsilon=0.01$ | 1,736 | 3.20 | 0.73 | 0.92 | 0.94 | 1.47 | 2.07 | 0.00 | 3.65 | 4.13 | 2.16 | 1.57
|  |  |  $k$-Means | $k=6$ | 1,887 | 3.94 | 0.00 | 9.22 | 9.05 | 17.73 | 21.47 | 0.00 | 0.24 | 0.52 | 0.47 | 0.36

Table 15 then reports the performance of the branch-and-price algorithm, with
each state clustering algorithm in the pricing problem. Again, to perform an
apples-to-apples comparison, the death toll is computed using the full
continuous-state contagion models, as opposed to either state-clustering
approximation. Note, at the outset, that both algorithms lead to similar
computational times in branch-and-price, which is expected since both scale
linearly with the number of states to cluster. However, the $\ell_{\infty}$
algorithm consistently leads to stronger solutions than the $k$-means
algorithm, resulting in hundreds of extra lives saved over a six-week horizon.
This result underscores the benefits of controlling the worst-case deviation
between the clustered states and the true states via the $\ell_{\infty}$
algorithm—as noted theoretically in Proposition 4.3 and experimentally in
Table 14. The smaller cluster diameters mitigate error propagation over
dynamical systems—as established in Proposition 4.4—and can ultimately result
in superior resource allocation decisions from the optimization methodology
developed in this paper.

Table 15: Comparison of vaccine allocation solutions obtained with
$\ell_{\infty}$ clustering vs. $k$-means clustering.

|  |  $n=51$, $S=6$, $D=6$ |  $n=51$, $S=6$, $D=11$ |  $n=51$, $S=6$, $D=21$
---|---|---|---|---
Hyperparameter | Method | $|S|$ | Time (sec.) | Deaths | $|S|$ | Time (sec.) | Deaths | $|S|$ | Time (sec.) | Deaths
— | Do nothing | — | — | 573.37K | — | — | 573.37K | — | — | 573.37K
$\varepsilon=0.01$ |  Branch-and-price ($k$-means) | 1,887 | 4.99 | -8.31K | 3,111 | 13.4 | -10.57K | 4,641 | 32.6 | -10.95K
|  Branch-and-price ($\ell_{\infty}$) | 1,736 | 4.91 | -10.37K | 2,656 | 12.5 | -10.94K | 4,328 | 26.1 | -11.10K
| Extra lives saved | — | — | $+24.8\%$ | — | — | $+3.5\%$ | — | — | $+1.4\%$
$\varepsilon=0.002$ |  Branch-and-price ($k$-means) | 5,202 | 13.11 | -10.31K | 12,342 | 62.4 | -10.70K | 29,172 | 439.1 | -11.02K
|  Branch-and-price ($\ell_{\infty}$) | 5,555 | 12.61 | -10.66K | 12,834 | 65.4 | -11.07K | 29,020 | 426.2 | -11.28K
| Extra lives saved | — | — | $+9.7\%$ | — | — | $+3.4\%$ | — | — | $+1.8\%$

  * •

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