# Optimizing urban mobility for saving lives and economy during an epidemic outbreak, with application to Grenoble

Ujjwal Pratap, Carlos Canudas de Wit, Federica Garin

► To cite this version:

Ujjwal Pratap, Carlos Canudas de Wit, Federica Garin. Optimizing urban mobility for saving lives and economy during an epidemic outbreak, with application to Grenoble. IEEE Transactions on Control Systems Technology, 2025, 33 (1), pp.288-303. 10.1109/TCST.2024.3477990. hal-04711085

HAL Id: hal-04711085

https://hal.science/hal-04711085v1

Submitted on 26 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Optimizing urban mobility for saving lives and economy during an epidemic outbreak, with application to Grenoble

Ujjwal Pratap, Carlos Canudas-de-Wit and Federica Garin

***Abstract***—This paper addresses the problem of controlling human mobility in order to mitigate an epidemic in a city. We consider a discrete-time human mobility model that captures daily mobility pattern between residences and different destinations in a city and also incorporates epidemic spread at each location. For this city-wide model, we provide techniques to compute optimal mobility control policies, which tune the operating capacities of different destinations depending on their type. To obtain this kind of policies, we solve an optimization problem that takes into account the current epidemic status, and maximizes the socio-economic activity while keeping the total infections below a desired threshold. The proposed solution techniques use an outer approximation method, thanks to the monotonic nature the problem, and a receding horizon approach. We apply these techniques to the mobility network of Grenoble metropolitan area.

***Index Terms***—optimal control, urban human mobility, epidemic control, non-pharmaceutical interventions, monotonic optimization, large-scale networks, receding horizon control.

## I. INTRODUCTION

Human mobility is essential to a country or region's economy. However, it can also be responsible for the spread of an epidemic when infected people mix with the susceptible population in common places. For this reason, the government devises policies to restrict movement in the event of an epidemic. Recently, during the Covid-19 pandemic, it was seen that unrestricted mobility can lead to overburdening of the health care infrastructure and thus to a significant increase in the loss of lives. In order to mitigate the spread of epidemics and their consequences for public health, government policies tend to restrict human mobility in a conservative manner, with severe consequences for social well being and the economy [1]. Therefore, mobility restriction policies should act parsimoniously by maximizing social well being and economic activities (minimizing mobility restriction) while respecting the constraints imposed by the available health infrastructure.

Based on a rich literature on epidemics modelling (see e.g. [2], [3]), there has been a wide interest on epidemic control,

This work was supported by European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (ERC-Adg no. 694209, Scale-FreeBack, website: http://scale-freeback.eu/).

Ujjwal Pratap was with Univ. Grenoble Alpes, CNRS, Inria, Grenoble INP, GIPSA-lab, Grenoble, France and is now with Technische Universität Berlin, Germany. Email: ujjwaliitp@gmail.com.

Carlos Canudas-de-Wit and Federica Garin are with Univ. Grenoble Alpes, CNRS, Inria, Grenoble INP, GIPSA-lab, Grenoble, France. Email: carlos.canudas-de-wit@cnrs.fr, federica.garin@inria.fr.

in particular after the recent outbreak of Covid-19 motivated several studies concerning control and mitigation of epidemic spread. The epidemic control strategies discussed in these studies were either based on pharmaceutical methods such as vaccines or medicines, or non-pharmaceutical interventions (NPIs) [4] to reduce the infection and/or the contact rate. In an early stage of an epidemic outbreak, when there is a lack of vaccines and medicines, non-pharmaceutical control strategies play a crucial role. However, these NPIs can be burden on the economy and social well-being, and hence should be properly planned and optimized to cater the specific economic and social settings of the region in consideration. Several works such as [5]–[7] highlight the importance to optimally design interventions, so as to have a minimal impact. However, in these works the intervention is abstractly modelled as lowering the disease transmissibility, leaving out the practical details of policy implementation via non-pharmaceutical interventions. Many papers along these lines, such as [8]–[11], separate the effect of testing-and-isolating and mobility restrictions. Some of them, e.g., [9], [10], present interesting quantitative estimates of the economic cost associated with the NPIs; [9] also highlights the usefulness of introducing different interventions in different age groups. However, in all these papers, the mobility restrictions are represented in a very abstract form, as a change of infection rate.

The correlation between mobility and disease transmissibility has been studied by many, for example in [12], [13]. Some studies have gone further, explicitly considering mobility restrictions as control variables. Among these, there are several works that focus on the cross-region mobility. [13] analyses the effect of mobility restrictions put by government bodies during Covid-19 across Italian regions, and [14] quantifies a trade-off between allowed cross-region mobility and number of infections and gives optimal restriction policies. [15] proposes reinforced learning methods to find mobility restriction policies to and from high risk regions by using aggregated inter-regional mobility data instead of individual traces. In addition to such data-based methods, there are also models integrating cross-region mobility and epidemic spread by SIR dynamics, see e.g. [16] and [17]. Some recent works such as [18], [19] exploit cross-region models to find optimal epidemic mitigation policies.

Above-mentioned contributions are all focused on mobility between regions or states, while it is crucial to consider the everyday mobility within each region or even town, e.g. to workplaces, schools, shopping, that contributes significantly to
the epidemics spread. Sharp restrictions of such intra-region mobility correspond to the well-known lock-down policies, while milder closures should be devised, with selective mobility restrictions based on the purpose of movement. Policies restricting intra-region mobility have been considered in few papers. [20] and [21] consider partial restrictions of mobility towards different categories of destinations, such as school, work, and others. They both use a compartmental model of epidemics and they model the effect of mobility restrictions as a modification of the infection rate. The effect of mobility on the infection rate is inferred from data about the first Covid-19 waves, where reduction in mobility is obtained from google mobility reports. Uncertainty in the parameters is taken into account by a stochastic MPC approach. However, these papers only consider aggregated destinations without specific geographic locations nor a mobility model: they heavily rely on data in order to infer the correlation between mobility restrictions and the infection rate, which is only possible after some restrictions have already been put in place. Hence, these methods might be difficult to apply in the early stages of a new epidemic, before any restrictions have been put in place.

Instead, our approach relies on an integrated model of mobility and epidemics, which only requires some publicly available census and travel-survey data obtained in non-restricted mobility conditions to calibrate the mobility model. This allows to understand the effects on the epidemic evolution of mobility, possibly at a refined granularity level when including a large number of origins (e.g. municipalities or neighborhoods) and of destinations (not only divided in categories such as workplaces, schools etc, but also with several nodes per category). This approach was introduced in [22], where a dynamic mobility model is defined, which incorporates epidemic spread at each geographical location in an urban city with SIR dynamics, in continuous time. This model is then used to formulate optimal control problems for epidemic mitigation while maximizing the economic activity, with application on an academic example with a network of six nodes only. [23] extended the integrated epidemic-mobility model to large scale-urban networks, giving details on how to obtain the model parameters from publicly-available data with the example of Grenoble urban area, and also providing a discrete-time model that preserves the original continuous-time model properties (non-negativity and mass conservation) while ensuring numerical stability.

In this paper, based on the numerical version of the combined epidemic-mobility model in [23], we extend the optimal control approach in [22] to a realistic scenario of a large-scale urban city having a mobility network with thousands of nodes. We propose solution techniques to devise mobility policies which maximize the socio-economic activity while respecting the ICU-hospital constraints, in a large-scale urban area. In order to address the computational complexity and obtain scalable algorithms, we exploit a receding-horizon approach, and consider policies that are constant over such horizon. We also explore the difference between imposing a same restriction for all destinations or instead devise category-dependent rules. The latter is very useful in case some activities are given high priority, but also leads to a more complex optimization problem,

which we can address exploiting monotone optimization techniques.

We illustrate our methods and their scalability by providing effective mobility control policies for Grenoble metropolitan area, using the mobility network developed in [23].

## II. URBAN HUMAN MOBILITY AND EPIDEMIC MODEL

As in [22] and [23], consider the human mobility in an urban area between two types of locations: origins and destinations. The origins correspond to the residential areas and the destinations, on the other hand, correspond to locations that people visit daily for some purpose like work, shopping, education or leisure for example, companies, research centers, schools, restaurants, cinemas, etc. Every day a certain number of people visit these destinations during specified time intervals and then return later to their residences on the same day.

### A. Origins and destinations

Consider a daily human mobility in the timescale of hours represented by a network $G = (V_o, V_d, \mathcal{E})$ between locations of two types: origins $V_o$ and destinations $V_d$. The assumptions on the model are as follows:

* The total population of the city remains constant.

* The mobility occurs only between pairs of origins and destinations, and not among a pair of origins or a pair of destinations.

* The number of people who visit destination $j$ from origin $i$ during a day is equal to the number of people who return to $i$ from $j$ on the same day.

**Destination categories:** The destinations $V_d$ are partitioned into $q$ categories $D_1, \dots, D_q$, where each category $D_h$ could also be further partitioned into sub-categories. In this paper, as will be seen in later sections, we consider the large scale mobility network of Grenoble which was developed in [23]. This network has 60 origin nodes and 374 destination nodes partitioned in 5 categories. The categories in this case are schools, workplaces, hospitals, shopping centers and leisure.

### B. Population, operating capacity and control

Denote the total population of origin $i \in V_o$ by $P_i$ which is the total number of people residing in $i$ and let $C_i$ be the nominal capacity of the destination $i$ if $i \in V_d$. The nominal capacity of a destination is defined as the total number of people that can be present at a given time in normal times that is when there is no restriction. We introduce the possibility to describe the restrictions of capacity, such as those that were imposed by many governments during peaks of Covid-19 pandemic. Such capacity reduction is decided per category of destinations, and is described by a control coefficient $u_h(k) \in [0, 1]$, where $h$ represents the destination category $D_h$. The control $u_h(k)$ tunes the nominal capacity $C_j$ of the destinations $j$ in $D_h$ at time step $k$ to determine their allowed operating capacity which can be defined as

$$C_j^o(k) = C_j u_h(k), \quad \text{for } j \in D_h. \quad (1)$$

In (1), if for a destination $j \in D_h$, $u_h(k) = 0$, then the operating capacity $C_j^o(k) = 0$, which means that no person
is allowed to visit $j$. On the other hand, if $u_h(k) = 1$, then $C_j^o(k) = C_j$ that is the operating capacity of $j$ is equal to its nominal capacity. For the sake of simplicity and readability, here onward, we will use the notation $Z_j(k)$ for the operating capacity $C_j^o(k)$ if $j$ is a destination or the population $P_j$ if $j$ is an origin unless otherwise stated. Note that, if $j$ is an origin, $Z_j(k)$ will be constant and equal to the population $P_j$, however, if $j$ is a destination then $Z_j(k)$ will be equal to the operating capacity $C_j^o(k)$ which varies according to the control $u_h(k)$ if $j \in D_h$ as defined in (1).

### C. Discrete-time mobility model

Let $N_i(k)$ be the number of people at location $i \in V_o \cup V_d$ at time step $k$ and $\Delta t$, the time-step size, then, the discrete-time mobility model is

$$N_i(k+1) = N_i(k) + \Delta t \sum_{j \in \mathcal{N}_i} (\phi_{ji}(k) - \phi_{ij}(k)), \quad (2)$$

where $\mathcal{N}_i = \{j : (i, j) \in \mathcal{E}\}$ is the set of $i$'s neighbors and the flow $\phi_{ij}(k)$ is the flow from $i$ to $j$ defined as

$$\phi_{ij}(k) = \min(\Delta_{ij}(k), \Psi_j(k)), \quad (3)$$

where $\Delta_{ij}(k)$ is the demand of $i$ with respect to $j$ which is defined as

$$\Delta_{ij}(k) = \min \left( \delta_{ij}(k) f_{ij}(k), \frac{N_i(k) \alpha_{ij}(k)}{\Delta t} \right), \quad (4)$$

where $f_{ij}(k) = M_{ij}u_h(k)$ for $j \in D_h$ and $M_{ij}$ is the total number of people that would like to visit $j$ from $i$ daily and $\delta_{ij}(k)$ is the Demand Gating Profile (DGP), such that

$$\left\{ \begin{array}{ll} \delta_{ij}(k) > 0 & \text{if } k\Delta t \bmod 24 \in [t_{ij}, t_{ij} + \tau_{ij}) \\ \delta_{ij}(k) = 0 & \text{otherwise} \end{array} \right.$$

satisfying $\sum_{k\Delta t \in [0,24]} \delta_{ij}(k) = 1$. Here, $[t_{ij}, t_{ij} + \tau_{ij}) \subseteq [0, 24)$ is called the mobility window, defined as the time interval in which there is mobility from $i$ to $j$. In the second term of (4), $\alpha_{ij}$ is the proportion in which the number of people at location $i$ will move to its neighbors $j$. It is computed as

$$\alpha_{ij}(k) = \frac{\delta_{ij}(k) f_{ij}(k) \psi_j(k)}{\sum_{q \in \mathcal{N}_i} \delta_{iq}(k) f_{iq}(k) \psi_q(k)},$$

which ensures that $\sum_{j \in \mathcal{N}_i} \alpha_{ij}(k) = 1$.

The supply of $j$ with respect to $i$, $\Psi_{ij}(k)$ is defined as

$$\Psi_{ij}(k) = \psi_j(k) \frac{(Z_j(k) - N_j(k))\pi_{ij}(k)}{\Delta t}. \quad (5)$$

Here, $\psi_j(k)$ is the Supply Gating Function (SGF) defined as

$$\psi_j(k) = \begin{cases} 1 & \text{if } k\Delta t \bmod 24 \in [a_j, b_j) \\ 0 & \text{otherwise} \end{cases} \quad (6)$$

with $[a_j, b_j) \subseteq [0, 24)$, where $a_j$ and $b_j$ are respectively the opening and closing time of location $j$. In other words, the Supply Gating Function $\psi_j(k)$ is 1 when $j$ is open and 0 at the times $j$ is closed. This function tells us if mobility is allowed to a location or not, therefore it is relevant only in case

$j$ is a destination because origins are considered to be always open i.e., $\psi_j(k) = 1$ if $j \in V_o$. Similar to DGPs, also SGFs are defined over a day and repeated periodically, possibly with different values for different days of the week. In (5), $\pi_{ij}(k)$ is the proportion in which the remaining capacity of the location $j$ is available for its neighbors and is defined as

$$\pi_{ij}(k) = \frac{\delta_{ij}(k) f_{ij}(k)}{\sum_{q \in \mathcal{N}_j} \delta_{qj}(k) f_{qj}(k)}. \quad (7)$$

Notice that $\sum_{i \in \mathcal{N}_j} \pi_{ij}(k) = 1$. It is important to notice that the flow $\phi_{ij}(k)$ depends on the control $u(k)$ as it tunes the operating capacities of the destinations which tune the demand $\delta_{ij}(k)$ and supply $\Psi_{ij}(k)$ through the terms $f_{ij}$ and $Z_j(k)$ respectively.

Moreover, it has been shown in [23] that in the mobility model (2)-(7), the number of people $N_i(k)$ at any location $i$ is always non-negative, i.e., $N_i(k) \ge 0$, and it also remains bounded that is it remains below the population of $P_i$ or the operating capacity $C_i^o(k)$, if it is an origin or a destination respectively, i.e., $N_i(k) \le Z_i(k)$. Moreover, the total mass is also preserved: the total number of people in the system remains constant, i.e., for any time $k$, we have $\sum_{j \in V_o \cup V_d} N_j(k) = \sum_{j \in V_o \cup V_d} N_j(0)$. We will use the notation $P$ to indicate such constant number of people in the system.

### D. SIR-mobility model

One of the properties of model (2)-(7) is that it is modular. Out of many applications for the mobility model, one is to see the effect of the restrictions on mobility, done through tuning the operating capacities of destinations using controls $u_h(k)$, on epidemic spread in an urban area. In [22], a continuous-time mobility model that integrates mobility with epidemic is given. The resulting model is called the SIR-mobility model. At each location, the disease spread follows the classical SIR model, but people are also allowed to move from one location to another (irrespective of infection level) and such mobility happens according to the mobility model (2)-(7). We consider the discrete-time SIR-mobility model proposed in [23], where $S_i(k)$, $I_i(k)$ and $R_i(k)$ denote the number of susceptible, infected and recovered people at location $i$ and time $k\Delta t$. We will also use the notation $I(k)$ to denote the total number of infected people at time $k\Delta t$, i.e., $I(k) = \sum_{i \in V_o \cup V_d} I_i(k)$.

For the sake of simplicity, we denote time $k+1$ with superscript $+$ and time $k$ with no superscript. With this notation, the SIR-mobility model is given by

$$S_i^+ = S_i - \Delta t \beta_i S_i^+ \frac{I_i}{N_i} + \Delta t \sum_j \left( \phi_{ji} \frac{S_j}{N_j} - \phi_{ij} \frac{S_i}{N_i} \right) \quad (8a)$$

$$I_i^+ = I_i + \Delta t \beta_i S_i^+ \frac{I_i}{N_i} - \Delta t \gamma I_i^+ + \Delta t \sum_j \left( \phi_{ji} \frac{I_j}{N_j} - \phi_{ij} \frac{I_i}{N_i} \right) \quad (8b)$$

$$R_i^+ = R_i + \Delta t \gamma I_i^+ + \Delta t \sum_j \left( \phi_{ji} \frac{R_j}{N_j} - \phi_{ij} \frac{R_i}{N_i} \right), \quad (8c)$$

where the parameters $\beta_i(k) = \bar{\beta}_i \frac{N_i(k)}{X_i}$, $\bar{\beta}_i$ and $\gamma_i$ are the infection rate, nominal infection rate and recovery rate at $i$. Here, $X_i = P_i$ if $i \in V_o$ whereas $X_i = C_i$ if $i \in V_d$. It is
important to note that when $i \in \mathcal{V}_o$, the nominal infection rate $\bar{\beta}_i$ is very low as in the origins, people live in their individual houses and hence the infection spread remains small. The equations (8) incorporate a non-local discretization of the continuous-time SIR dynamics proposed in [23] which is inspired by [24].

They can be re-written in the following equivalent way, with explicit time-update equations which allow an easy implementation:

$$S_i^+ = \frac{1}{1 + \Delta t \beta_i \frac{I_i}{N_i}} \left[ S_i + \Delta t \sum_j \left( \phi_{ji} \frac{S_j}{N_j} - \phi_{ij} \frac{S_i}{N_i} \right) \right] \quad (9a)$$

$$I_i^+ = \frac{1}{1 + \Delta t \gamma} \left[ I_i + \Delta t \beta_i S_i^+ \frac{I_i}{N_i} + \Delta t \sum_j \left( \phi_{ji} \frac{I_j}{N_j} - \phi_{ij} \frac{I_i}{N_i} \right) \right] \quad (9b)$$

$$R_i^+ = R_i + \Delta t \gamma I_i^+ + \Delta t \sum_j \left( \phi_{ji} \frac{R_j}{N_j} - \phi_{ij} \frac{R_i}{N_i} \right). \quad (9c)$$

As shown in [23], (9) along with the mobility model (2)-(7) ensures that $S_i$, $I_i$ and $R_i$ remain non-negative and bounded, that is, $0 \le S_i(k) \le N_i(k)$, $0 \le I_i(k) \le N_i(k)$ and $0 \le R_i(k) \le N_i(k)$, and moreover $S_i(k) + I_i(k) + R_i(k) = N_i(k)$ for all $i$ and $k$.

Notice that $S_i$, $I_i$, $R_i$ at any location $i$ depend on $\mathbf{u}$ as $\phi_{ij}$, $N_i$, $N_j$ are implicit functions of the control $\mathbf{u}$. Note that, we use bold $\mathbf{u}$ to denote the vector $\mathbf{u} = [u_h]$.

### E. Extension possibilities for SIR-mobility model

The focus of this paper is on the effects of partial mobility restrictions, which are well described by the SIR-mobility model (2)-(7) and (9) presented above. It is interesting to notice that small modifications of the model can be easily done to incorporate further interventions, and in particular testing-and-isolating policies. Indeed, one could either introduce an extra compartment in addition to $S_i$, $I_i$, $R_i$, so as to separate the infected people in two classes, isolated or not, or one could use a simpler model of isolation, by setting that only a fraction of the infected people is allowed to move, i.e., to participate to a flow $\phi_{ij}$. Let's discuss the latter option in some more detail, setting for example that a fraction $\kappa_i$ of infected people $I_i$ is not allowed/not willing to move, representing the fact that people feeling significantly sick remain at home to rest and/or that some fraction of infected people is detected by testing and isolated. We might introduce such fraction term only at origins, while leaving $\kappa_i = 0$ at destinations. The fraction $\kappa_i$ could be a constant, or a control variable to be optimized, representing some more or less active testing campaign, with an associated cost. The terms $\kappa_i$ would naturally modify the mobility terms $\phi_{ij} \frac{S_i}{N_i}$, $\phi_{ij} \frac{I_i}{N_i}$, $\phi_{ij} \frac{R_i}{N_i}$ in SIR equations (8) (or equivalently (9)), where the flow $\phi_{ij}$ of people moving from node $i$ to $j$ is composed of susceptible, infected and recovered people in the same proportion as in the node $i$. In the modified equations, the proportions should rather take into account only a fraction $(1 - \kappa_i)$ of the infected people at $i$, so that

the modified mobility terms would be $\phi_{ij} \frac{S_i}{S_i+(1-\kappa_i)I_i+R_i}$, $\phi_{ij} \frac{(1-\kappa_i)I_i}{S_i+(1-\kappa_i)I_i+R_i}$, and $\phi_{ij} \frac{R_i}{S_i+(1-\kappa_i)I_i+R_i}$.

However, this is not the only modification: also the demands $\Delta_{ij}$ involved in the computation of flows $\phi_{ij}$ are modified, to represent the fact that less people are willing/allowed to move, since $\kappa_i I_i$ people surely remain in $i$. More precisely, in (4), the term $f_{ij}$ should be multiplied by a factor $\frac{S_i+(1-\kappa_i)I_i+R_i}{N_i}$ and the term $\frac{N_i \alpha_{ij}}{\Delta t}$ should be replaced by $\frac{(N_i-\kappa_i I_i) \alpha_{ij}}{\Delta t}$.

It is easy to show that with such modifications, the properties of non-negativity and boundedness of all variables and of conservation of the total mass are retained, with a similar proof as for the SIR-mobility model in [23].

Such modified SIR-mobility model could be used for similar optimization as the one presented in this paper, in the presence of a constant testing-and-isolating policy, or for jointly optimizing the mobility restrictions and the testing-and-isolating policy.

Notice that the above-mentioned modified model only accounts for the direct effects of testing-and-isolating and of mobility restrictions, while it does not incorporate rebound effects, where more testing might lead to lower adherence with social distancing rules [25].

In the remainder of this paper we will only focus on the mobility restrictions, using the SIR-mobility model without isolation term, so as to better highlight the benefits of such a control action.

### III. OPTIMAL CAPACITY CONTROL FOR EPIDEMIC MITIGATION

The SIR-mobility model (9) describes the evolution of susceptible, infected and removed population at each location in the network and also incorporates the mobility of people between different locations. It can be noticed from (9) that by applying restrictions in mobility using the control input $u_h(k)$, the infection spread in the network can also be controlled. Using this property, in [22], an optimal capacity control problem has been formulated which aims at maintaining infections below a threshold by mobility control while maximizing the economic activity defined as a function of number of people visiting the destinations. Here, we provide solution techniques for this optimal capacity control problem for a large-scale urban network and apply it to such a mobility network of Grenoble, a metropolitan city in France.

#### A. Control Action

Our main goal is to design mobility control policies which depend on different destination categories and can be updated on regular time intervals. Therefore, let us consider a time horizon of interest $T$, and $T_u$ be the duration after which a policy affecting the capacity of destinations in different categories is updated (the policy update period) and all the parameters of the model (2)-(7) and (9) be given. It is important to notice that the update period $T_u$ is a design parameters to be decided by the policy makers. It is suggested to use such $T_u$ that can be socially acceptable, taking into account the fact that extremely short values of $T_u$, of the order e.g. of a day, would create
confusion and lack of respect of the rules. However, weekly or
fortnightly updates could be realistic and socially acceptable.

Let $\mathbf{u} \in \mathcal{U}$, where $\mathcal{U}$ is the set of vectors of piecewise-
constant functions $[u_1(k), \dots, u_q(k)]$ as defined below

$$
u_h(k) = \begin{cases}
\mu_h^1 & \text{if } k\Delta t \in [0, T_u] \\
\mu_h^2 & \text{if } k\Delta t \in (T_u, 2T_u] \\
\vdots & \vdots \\
\mu_h^p & \text{if } k\Delta t \in ((p-1)T_u, T]
\end{cases} \tag{10}
$$

with $p = \lceil T/T_u \rceil$, where $[\cdot]$ is the ceiling function, with
$\mu_h^\ell \in [0, 1]$ for every $\ell \in \{1, \dots, p\}$.

Recall that in (2)-(7) and (9) there is a dependence on the
capacity control **u**, due to the fact that **u** tunes the operating
capacities, modifying the flows $\phi_{ij}$ and hence also all $N_i$, $S_i$,
$I_i$ and $R_i$.

To highlight this implicit dependence, from now on, we will
denote the number of people and number of infected people at
location *i*, at time *k*Δ*t* by *N*<sub>*i*</sub>(*k*, **u**) and *I*<sub>*i*</sub>(*k*, **u**) respectively,
when the system is operated with capacity control **u**.

B. *Objective function: Socio-economic activity*

Mobility restrictions have severe economic impact (e.g., affecting both production and demand of goods) and social impact (e.g., reducing school attendance and other cultural activities and social interactions). As suggested in [22], this can be represented by an index which counts the number of people (or more precisely the occupancy, i.e., number of people divided by maximum capacity) of all categories of destinations; some coefficients can be introduced, to give different weights to different categories, so as to describe different relevance of such activities, and allowing the policy-makers to prioritize the most crucial ones.

Hence, we consider the socio-economic activity at a given
time *k* defined as

$$
E(k, \mathbf{u}) = \sum_{j \in \mathcal{V}_d} \frac{\chi_j}{C_j} N_j(k, \mathbf{u}),
$$

where $N_j(k, \mathbf{u})$ is the number of people at destination $j$ at time $k$, $C_j$ is the nominal capacity of $j$ and $\chi_j > 0 \in \mathbb{R}$ is the weight assigned to the destination $j$ according to its socio-economic importance in the urban network.

One of the choices for the weight $\chi_j$ is

$$
\chi_j = \frac{\bar{\chi}_h}{\bar{\chi}_h^{nom}}, \quad \text{for } j \in \mathcal{D}_h \qquad (11)
$$

where $\bar{\chi}_h > 0 \in \mathbb{R}$ is a weight given to the destination category $h$ by the policymakers, representing the socio-economic importance given to this destination category, and $\bar{\chi}_h^{nom}$ is defined as

$$
\bar{\chi}_h^{nom} = \sum_{k=0}^{T/\Delta t} \sum_{i \in \mathcal{D}_h} \frac{N_i(k, \mathbf{1}_q)}{C_i},
$$

where $\mathbf{1}_q \in \mathbb{R}^q$ is the vector of all ones. In other words, $\bar{\chi}_h^{nom}$ is the economic activity done in the destinations of category $h$ when they are all open to full capacity. The effect of economic weight will be evident when devising the category-dependent control policies in Section V.

The weights $\bar{\chi}_h$ must be devised by policy-makers, to
describe their priorities. Economists can provide quantitative
estimations of the economic cost, e.g., in terms of percentage
of GDP (see [9], [10]), while social impact might be more
subtle to quantify, but nevertheless is crucially important, since
it represents effects e.g. on education (with economic impact in
future years), societal cohesion, and mental health. Notice that
our objective function depends on the number of people that
move to each category, so that the weights $\bar{\chi}_h$ represent the
importance given to each category relative to other categories.
Differently from other works such as [9], there is no need
to give a quantitative cost to human lives, since they do not
appear as a competing term in the objective function. Indeed,
we will see in Subsect. III-C that they are rather taken into
account in the constraints of the optimization problem.

C. *Optimal capacity control problem*

Our goal is to keep the socio-economic activity as high as
possible, while ensuring that the number of infected people
remains below a positive threshold $\bar{I}$ that reflects the capacity
of the health-care system, e.g. in relation with the number of
ICU beds. More precisely, our goal is to maximize the total
socio-economic activity $L(\mathbf{u})$ done in the given time interval
of interest $[0, T]$, defined as $L(\mathbf{u}) = \sum_{k=0}^{T/\Delta t} E(k, \mathbf{u})$, i.e,

$$
L(\mathbf{u}) = \sum_{k=0}^{T/\Delta t} \sum_{j \in \mathcal{V}_d} \frac{\chi_j}{C_j} N_j(k, \mathbf{u}). \quad (12)
$$

The constraint that ensures that the capacity of the health-
care system is never exceeded can be expressed as

$$
\sum_{j \in \mathcal{V}_o \cup \mathcal{V}_d} I_j(k, \mathbf{u}) \leq \bar{I} \quad \forall k \text{ s.t. } k\Delta t \in [0, T],
$$

or, equivalently, as $I_{\max}(\mathbf{u}) \le \bar{I}$, where

$$
I_{\max}(\mathbf{u}) = \max_{k=0,...,T/\Delta t} \sum_{j \in \mathcal{V}_o \cup \mathcal{V}_d} I_j(k, \mathbf{u}) \quad (13)
$$

is the maximum over the time interval $[0, T]$ of the total
number of infected people. In some cases where the interval of
interest will be $[T_0, T_f]$ instead of $[0, T]$, then both $L(\mathbf{u})$ and
$I_{\max}(\mathbf{u})$ will be computed in the time interval $[T_0, T_f]$ with $k$
in (12) and (13) ranging from $T_0/\Delta t$ to $T_f/\Delta t$.

With the above notation, the optimal control policy is
obtained by solving the following optimization problem:

$$
\max_{\mathbf{u} \in \mathcal{U}} \{L(\mathbf{u})\} \qquad (14a)
$$

$$
\text{subject to } I_{\max}(\mathbf{u}) \le \bar{I}, \qquad (14\text{b})
$$

where, $L(\mathbf{u})$ and $I_{\max}(\mathbf{u})$ are defined in (12), (13),

together with (2)-(7) and (9), and $\mathcal{U}$ is the set of vectors of piecewise-constant functions as in (10).

The problem (14) can be easily solved for smaller networks
using solvers like `fmincon` in MATLAB as shown in [22] but
for large-scale urban network like that of Grenoble in [23], the
solution becomes hard or even impossible to obtain by naive
use of off-the-shelf solvers. In this paper, at first we will start
with a simple variation of the problem (14) by taking the same
<figure><img src="image_2.png" /><figcaption>**Fig. 1:** Proportion of infected population *I*/*P* in Grenoble area when all the destinations are open up to full capacity (i.e., with no restriction in place). The dashed line depicts the threshold *I*̄/*P* representing the capacity of the health-care system.</figcaption></figure>

control for all the destinations which we call category-free policy. Then, we proceed towards the problem when we have different control policy applied to different destination based on their categories and call it category-dependent policy.

### D. Benchmark control policy

We consider as a benchmark the case where the operating capacity of all the destinations is equal to the nominal capacity that is $u_h(k) = 1$ for all $h$ and $k$. In other words, it is the case when there is no mobility restriction. The economic activity in that case is defined as the nominal economic activity denoted by $\bar{E}(k)$.

The aim of the optimization is to keep the epidemic spread bounded but to keep the economic activity as close as possible to that of the benchmark case. Therefore, for comparison between different scenarios, we define the relative economic activity $R_E(k)$ at time $k$ as

$$R_E(k) = \frac{\sum_{\ell \le k} E(\ell)}{\sum_{\ell \le k} \bar{E}(\ell)}. \qquad (15)$$

Note that $R_E(k) \in [0, 1]$, where $R_E(k) = 1$ corresponds to the benchmark scenario where there is no restriction, while $R_E(k) = 0$ corresponds to a complete lockdown.

### E. Running example

Throughout this paper, to illustrate our method and results, we consider the large-scale mobility network of Grenoble, developed in [23]. There are 60 origins and 374 aggregated destinations grouped in five destination categories, namely schools, hospitals, workplaces, shopping centers and leisure. The parameters related to mobility such as DGP, SGF, $M_{ij}$, $a_j$, $b_j$ for each destination category have been obtained from publicly available data such as travel surveys and population census, as explained in details in [23]. For the epidemic spread, the parameters have been taken from [26], [27] and curve fitting using data from [28], [29]. The epidemic parameters considered for simulations are $\gamma_i = 1/8$ per 24 hours and the nominal infection rate $\bar{\beta}_i = 0.25$ per 24 hours for $i \in \mathcal{V}_o$ and that for the destinations in categories schools, workplaces, hospitals, shopping centers, leisure are 2, 2, 4, 2, 4 per 24 hours

respectively. We take time horizon $T = 2520$ hours (15 weeks) in all the simulations and sum of total initial infections $\sum_i I_i(0) = 100$ which is distributed randomly only among the origins(residential areas). The initial conditions are

$$[N_i(0), I_i(0), S_i(0), R_i(0)] = [P_i, I_i(0), P_i - I_i(0), 0] \quad \text{for } i \in \mathcal{V}_o \qquad (16)$$

$$[N_i(0), I_i(0), S_i(0), R_i(0)] = \mathbf{0} \quad \text{for } i \in \mathcal{V}_d.$$

With this initialization, the total number of people in the system is $P = \sum_{i \in \mathcal{V}_o} P_i$.

In addition, the weights $\bar{\chi}_h$ that will be used in all the simulations in this paper are given in Table I. Note that these weights are simply meant as an example, to showcase possible effects of unbalanced weights on the socio-economic cost and control policies devised in later sections, and do not reflect actual governments' policies, nor they are meant as a recommendation.

<table><caption>**TABLE I:** Socio-economic weights for different destination categories used in simulations.</caption><thead><tr><th></th><th>Schools</th><th>Hospitals</th><th>Workplaces</th><th>Shopping centers</th><th>Leisure</th></tr></thead><tbody><tr><td>$\bar{\chi}_h$</td><td>0.1</td><td>1</td><td>0.7</td><td>0.8</td><td>0.2</td></tr><tr><td>$\bar{\chi}_h^{nom}$</td><td>$4.2 \times 10^5$</td><td>$1.5 \times 10^4$</td><td>$3.4 \times 10^5$</td><td>$9.7 \times 10^4$</td><td>$3.5 \times 10^5$</td></tr></tbody></table>

For all the plots in this paper, we will consider the weights mentioned in the table I (to be used in (11) and (12)) and the threshold on infection $\bar{I} = 0.15P$, where $P$ is the total population. This threshold has been chosen to be consistent with health-care capacity of Grenoble: the proportion of infected cases turning to ICU cases was $7.5 \times 10^{-4}$ in one of the initial waves in France [29]; with a population $P = 500 000$, the threshold 15% $P$ corresponds to 75 000 infected and hence 57 occupied ICU beds, lower than 80, the total number of ICU beds in Grenoble.

All the plots in this paper refer to the model (2)-(7) and (9) with parameters discussed in this section.

Figure 1 depicts the evolution of infection in Grenoble area with no restriction on the capacities of the destinations which is also the same as the benchmark policy. It can be seen that more than forty percent of people get infected after some weeks, which is largely above the threshold representing critical burden of the health infrastructure resulting in the loss of many lives. This highlights the need for some closure policy, to mitigate the infection level while avoiding an undesirable complete lock-down and thereby preventing socio-economic collapse.

## IV. CATEGORY-FREE POLICIES

In this section, we start with the simplest variation of the problem that restricts the optimization to the subset of control policies that have the same control for all destination categories, that is, $u_h(k) = u(k)$ for all $h$, where $u(k)$ is piece-wise constant function as in (10). This corresponds to
<figure><img src="image_3.png" /><figcaption>Fig. 2: Cost function $L(u1_q)$ and infection peak in proportion to the total population $I_{\max}(u1_q)/P$ as functions of the capacity control $u$, when a fixed control $u$ is applied over the interval $[0, T]$ (i.e., $T_u = T$), with $T = 15$ weeks.</figcaption></figure>

considering the following modification of the problem (14), where we restrict the search to a subset $\tilde{\mathcal{U}}$ of the set $\mathcal{U}$:

$$ \max_{u \in \tilde{\mathcal{U}}} \{L(u)\} \qquad (17a) $$

$$ \text{subject to } I_{\max}(u) \le \bar{I}, \qquad (17b) $$

where

$$ \tilde{\mathcal{U}} = \{u \in \mathcal{U} | u_{h_1} = u_{h_2}, \forall h_1, h_2\} $$

and $L(u)$ and $I_{\max}(u)$ are defined in (12), (13) together with (2)-(7) and (9). In Sect. IV-B, we will show how this problem can be solved efficiently, exploiting the monotonicity of the problem and applying a receding horizon approach, where a bisection search algorithm is applied over a sliding time window, whose length $T_m$ might be equal to $T$ or smaller.

### A. Fixed policy

In this subsection, we consider the optimization problem in the case in which a same fixed value of $u$ is to be applied for all categories and at all times in the interval $[0, T]$. In other words, we consider problem (17) (where a same restriction is applied to all categories) in the particular case where $T_u = T$. Indeed, when $T_u = T$, piecewise-constant functions as in (10) are actually constant functions. Recalling the notation $1_q \in \mathbb{R}^q$ for the all-ones vector, we will denote by $L(u1_q)$ the cost $L(u)$ in the case where $u$ is constant in time and is equal to the same scalar $u$ for all categories.

At first, we notice some monotonicity properties of problem (17) with $T_u = T$, and exploit them in order to find an efficient algorithm for its solution. First, we notice that $L(u1_q)$ is a non-decreasing function of $u$, as we can see from Figure 2. The intuitive explanation is that a larger $u$ allows more visitors to the destinations and hence more economic activity. For this reason, in order to maximize $L(u1_q)$, we can simply focus on maximising $u$. Then, we notice that also $I_{\max}(u1_q)$ is a non-decreasing function of $u$, except possibly for $u$ very near to zero. This property depends on the parameters of the SIR-mobility model, and corresponds to the scenario of interest, where infection rate

$\beta_i$ is significantly lower at origins than at destinations. The intuitive explanation is that larger $u$ allows more visitors to the destinations, where the infection rate is higher, and hence the number of infected grows more than when some of these persons are forbidden to reach their destinations (except possibly near $u = 0$, where destinations are so empty that the probability of infection there is small). Figure 2 shows $I_{\max}(u1_q)$ is indeed a non-decreasing function of $u$ in our example (described in Sect. III-E). Provided that $I_{\max}(0) < \bar{I}$, this property of $I_{\max}$ ensures that the feasible region is simply an interval $[0, \bar{u}] \subseteq [0, 1]$, where $\bar{u}$ is such that $I_{\max}(\bar{u}1_q) = \bar{I}$. Due to the above-discussed monotonicity of $L(u1_q)$, $\bar{u}$ is the desired optimal $u$, that maximizes $L(u1_q)$ under the constraint that $I_{\max}(\bar{u}1_q) \le \bar{I}$. Finding $\bar{u}$ can be done by bisection method, i.e., by running Algorithm 1 with $[T_0, T_f] = [0, T]$; Algorithm 1 is given on a more general interval $[T_0, T_f]$, because this will be then used within Algorithm 2 in Section IV-B.

Algorithm 1: Fixed policy in $[T_0, T_f]$

Input: Parameters for (2)-(7) and (9), initial condition $N_i(T_0), I_i(T_0), S_i(T_0), R_i(T_0) \forall i$, threshold $\bar{I}$, initial time $T_0$, final time $T_f$, maximum number of iterations $n_{\text{iter}}$ and tolerance $\epsilon_f > 0$.

Output: $u^* = u_{\min}$

1. Set $u_{\min} = 0$ and $u_{\max} = 1$.

2. **REPEAT**

3. Compute $u = (u_{\min} + u_{\max})/2$,

4. Run (2)-(7) and (9) with $u(k) = u1_q$ for $k = T_0/\Delta t, \dots, T_f/\Delta t$ and compute $I_{\max}(u1_q)$ as in (13)

5. IF $I_{\max}(u1_q) \le \bar{I}$

6. $u_{\min} = u$

7. ELSE, $u_{\max} = u$.

8. till $(\bar{I} - I_{\max}(u_{\min}1_q)) \le \epsilon_f$ or $n_{\text{iter}}$ iterations.

Before running Algorithm 1, it is necessary to verify that $I_{\max}(0) < \bar{I}$, which ensures that there is a feasible solution. Note that the obtained $u^*$ is such that $I_{\max}(u^*1_q) \le \bar{I}$ and moreover (unless the algorithm terminates because of reaching $n_{\text{iter}}$ iterations) $\bar{I} - I_{\max}(u^*1_q) < \epsilon_f$.

Under the above-mentioned monotonicity properties, this means that $u^*$ approaches the optimal solution $\bar{u}$ with any desired precision. Also notice that even in the cases where monotonicity does not hold, $u^*$ is ensured to be a feasible solution of the original problem (17).

### B. Category-free receding-horizon policy

As it was highlighted by [30], the duration and the starting time of the interventions play a crucial role. Hence, we shall go beyond the fixed-policy approach.

In this subsection, we consider the optimization problem in the case in which a same restriction is to be applied for all categories, while different restrictions can be applied over different intervals of time. In other words, we study problem (17) for general values of the update interval $T_u \le T$.
<figure><img src="image_4.png" /><figcaption>Fig. 3: Receding horizon approach with update period $T_u$ and the optimization horizon $T_m$ in the time interval $[0, T]$.</figcaption></figure>

To find a good, although possibly suboptimal, feasible solution of (17) in an efficient way, we adopt a receding horizon approach, where we find a fixed policy within an interval of optimization that is by applying Algorithm 1 over receding time intervals [$ℓT_u, ℓT_u + T_m$], where $ℓ = 0, \dots, p$, and the length of the optimization horizon is a parameter $T_m \in [T_u, T]$ to be designed. This approach has been depicted in Figure 3 and been summarised in Algorithm 2.

**Algorithm 2:** Category-free receding-horizon policy, with update period $T_u$ and optimization horizon length $T_m$

**Input:** $T, T_u, T_m$, Parameters for (2)-(7) and (9), initial condition $N_i(0), I_i(0), S_i(0), R_i(0) \forall i$, threshold $\bar{I}$, maximum number of iterations $n_{iter}$ and tolerance $\epsilon_f > 0$.

**Output:** $u^{*1}, \dots, u^{*p}$

1. Initialize $\ell = 1$. Compute $p = [T/T_u]$.

2. **REPEAT**

3. Compute $T_0 = (\ell-1)T_u$ and $T_f = (\ell-1)T_u + T_m$

4. Apply Algorithm 1 over $[T_0, T_f]$ with tolerance $\epsilon_f$ and maximum number of iterations $n_{iter}$, and denote by $u^{*\ell}$ the obtained output

5. Run (2)-(7) and (9) over the interval $[(\ell-1)T_u, \ell T_u]$, with fixed policy $u^{*\ell} \mathbf{1}_q$, so as to obtain $N_i(\ell T_u/\Delta t)$, $S_i(\ell T_u/\Delta t)$, $I_i(\ell T_u/\Delta t)$, $R_i(\ell T_u/\Delta t)$ for all $i$ (to be used as initialization for the next iteration)

6. $\ell = \ell + 1$

7. until $\ell = p$.

Notice that the piecewise constant function $u^*(k)$ defined by $u^*(k) = u^{*\ell}$ for all $k = (\ell-1)T_u/\Delta t, \dots, \ell T_u/\Delta t$ is a feasible solution of (17).

## V. CATEGORY-DEPENDENT POLICIES

During COVID-19 outbreak, government/policymakers implemented restrictions based on destination types according to their respective socio-economic importance. Therefore, in this section, we devise policies that apply different capacity controls to destinations in different categories in the optimal control problem (14). For this, we implement monotonic optimization techniques. In Subsections V-A and V-B, we introduce monotonic optimization [31], [32], particularized to our problem, while discussing some relevant implementation details. Then, in Subsection V-C, we integrate it to our receding-horizon approach.

<figure><img src="image_5.png" /><figcaption>Fig. 4: Cost function $L(\mathbf{u})$ for $T = 15$ weeks and when we freeze all variables except one. The free variable is increased in $[0, 1]$ and the others are fixed at 1. Here, $\mathbf{u} \in [0, 1]^5$ as we have $T_u = T$. The plots are obtained from two separate simulations in which a single variable is kept free namely $u_1$ and $u_3$ respectively.</figcaption></figure>

Throughout this section, for any two vectors, $\mathbf{u}, \mathbf{v} \in \mathbb{R}^n$, we use the notation $\mathbf{u} \le \mathbf{v}$ (resp. $\mathbf{u} < \mathbf{v}$) if $u_i \le v_i$ (resp. $u_i < v_i$) for all $i = 1 \dots n$.

### A. Monotonic optimization

We start with some preliminaries from [31], [32].

**Definition 1.** Non-decreasing functions: A function $f: \mathbb{R}_+^n \to \mathbb{R}$ is an non-decreasing function if $f(\mathbf{x}) \le f(\mathbf{y})$ when $0 \le \mathbf{x} \le \mathbf{y}$, where $\mathbb{R}_+^n$ denotes the set of non-negative real numbers.

**Definition 2.** Box: If $0 \le \mathbf{b}$, then the set of all $\mathbf{x} \in \mathbb{R}^n$ such that $0 \le \mathbf{x} \le \mathbf{b}$ is defined as the box $[0, \mathbf{b}]$ with the vertex $\mathbf{b}$.

**Definition 3.** Normal set: A set $G \subset \mathbb{R}_+^n$ is normal if for any element $\mathbf{x} \in G$, all other elements $\mathbf{x}'$ such that $0 \le \mathbf{x}' \le \mathbf{x}$ are also in set $G$. In other words, the set $G$ is normal if given any element $\mathbf{x} \in G$, then the box $[0, \mathbf{x}] \subset G$.

An optimization problem belongs to the class of monotonic optimization<fnref n="1" /> if it can be expressed in the following form:

$$ \max_{\mathbf{x}} F(\mathbf{x}) \qquad (18a) $$

$$ \text{subject to } \mathbf{x} \in G, \qquad (18b) $$

where, $F(\mathbf{x}) : \mathbb{R}_+^n \to \mathbb{R}$ is an non-decreasing function and $G \subset [0, \mathbf{b}] \subset \mathbb{R}_+^n$ is a compact normal set with nonempty interior.

In the optimal control problem discussed in this paper, the objective function $F(\mathbf{x})$ is given by $L(\mathbf{u})$ defined in (12) and the feasible set $G$ is defined as

$$ G = \{\mathbf{u} \in \mathcal{U} | I_{\max}(\mathbf{u}) \le \bar{I}\}. \qquad (19) $$

In the range of parameters of interest, where at least a total lockdown is enough to keep the number of infected people below $\bar{I}$, notice that $0 \in \mathbb{R}^n$ is always in $G$, therefore $G$ always

<footnote n="1">¹The general form of monotonic optimization problem from the literature [31], [32] also allows for a further constraint being a so-called co-normal set, but in our case only the normal constraint is involved, and hence we only include the normal set in this section.</footnote>
<figure><img src="image_6.png" /><figcaption>**Fig. 5:** Cost function $I_{\max}$ for $T = 15$ weeks and when we freeze all variables at feasible point $0.5625 * \mathbf{1}_5$ except one. The free variable is varied in $[0, 1]$ and the others are fixed at $0.5625 * \mathbf{1}_4$. Here, $\mathbf{u} \in [0, 1]^5$ as we have $T_u = T$. The plots are obtained from two separate simulations in which a single variable is kept free namely $u_1$ and $u_3$ respectively.</figcaption></figure>

has a non-empty interior. Also notice that for the problem (14), $G$ is enclosed in the box $[0, 1]$, where $\mathbf{1} \in \mathbb{R}^n$. In order to use monotonic optimization techniques to solve (14), we need two properties, namely the cost function $L(\mathbf{u})$ be non-decreasing and the feasible set $G = \{\mathbf{u} \in \mathcal{U} | I_{\max}(\mathbf{u}) \le \bar{I}\}$ be a normal set, where $I_{\max}(\mathbf{u})$ is computed using (2)-(7) and (9). At first, notice that the function $L(\mathbf{u})$ defined in (12) is a non-decreasing function of $\mathbf{u}$. The intuitive explanation is that given two control vectors $\mathbf{u}, \mathbf{v} \in \mathcal{U}$ such that $\mathbf{v} \ge \mathbf{u}$, implies that the number of visitors in destinations is higher with the control $\mathbf{v}$ than with the control $\mathbf{u}$. Therefore, more economic activity at the destinations which implies that $L(\mathbf{v}) \ge L(\mathbf{u})$. Figure 4 shows an example where when we freeze all variables in $\mathbf{u} \in \mathbb{R}^5$ except one variable, which we vary and we can see that indeed in this case $L(\mathbf{u})$ is non-decreasing which is according to our intuition.

In order to show that $G$ defined in (19) is a normal set, we need to show that if control $\mathbf{z} \in G$, then any $\mathbf{v} \le \mathbf{z}$ satisfies $\mathbf{v} \in G$. Consider $\mathbf{v} \le \mathbf{z}$, then the number of visitors in each destination is lower with the control $\mathbf{v}$ and higher with the control $\mathbf{z}$. Since the infection rate is higher at destinations, the number of infected people is more in case of more visitors at destinations (except possibly near $\mathbf{u} = \mathbf{0}$, where destinations are so empty that the probability of infection there is small). This certainly depends on the parameters of the SIR-mobility model, and corresponds to the scenario of interest, where infection rate $\beta_i$ is significantly lower at origins than at destinations. Therefore, if $I_{\max}(\mathbf{z}) \le \bar{I}$, then we have that also $I_{\max}(\mathbf{v}) \le \bar{I}$ for $\mathbf{v} \in [0, \mathbf{z}]$ except possibly for a very low threshold $\bar{I}$, which is outside our interest, since it would correspond to the case where even a total lockdown, i.e. $\mathbf{u} = \mathbf{0}$, would be unfeasible. It can also be seen in Figure 5, where we freeze all the variables at a feasible point and vary one variable in $[0, 0.7]$, that once a $\mathbf{z}$ is feasible, it remains feasible for all $\mathbf{u}$ such that $u_i < z_i$ for $i = 1$ while varying $u_1$ and for $i = 3$ while varying $u_3$.

## B. Polyblock outer approximation algorithm

Since the objective function (12) is non-decreasing, the optimal value lies on the boundary of the feasible set $G$ [31], [32]. However, we do not have the knowledge of the boundary of the feasible set $G$. Therefore, we aim to have an approximation of the boundary using some outer approximation methods notably Polyblock algorithm [31]. We start by defining a polyblock and related terms first.

**Definition 4.** **Polyblocks:** A set $\mathcal{P} \subset \mathbb{R}_+^n$ is called a polyblock if it is a union of a finite number of boxes $[0, \mathbf{z}]$, where $\mathbf{z} \in \mathcal{T}$ and $|\mathcal{T}| < +\infty$. The set $\mathcal{T}$ is called the vertex set of the polyblock. A vertex $\mathbf{v} \in \mathcal{T}$ is said to be proper if there is no $\mathbf{v}' \in \mathcal{T}$ such that $\mathbf{v}' \ne \mathbf{v}$ and $\mathbf{v}' \ge \mathbf{v}$. The vertices which are not proper are called improper vertices. Improper vertices can be removed from the vertex set $\mathcal{T}$ without affecting the shape of the polyblock.

**Definition 5.** **Upper boundary:** A point $\bar{x}$ of a normal closed set $G$ is called an upper boundary point of $G$ if $G \cap \{x \in \mathbb{R}_+^n | x > \bar{x}\} = \emptyset$. The set of all upper boundary points of $G$ is called its upper boundary and denoted by $\partial^+ G$.

**Definition 6.** **Projection on the upper boundary:** let $G \subset \mathbb{R}_+^n$ be a compact normal set with non-empty interior. For any point $\mathbf{z} \in \mathbb{R}_+^n \setminus G$, the line segment joining $\mathbf{0}$ to $\mathbf{z}$ meets the upper boundary $\partial^+ G$ at a unique point $\pi_G(\mathbf{z})$ defined as

$$ \pi_G(\mathbf{z}) = \lambda\mathbf{z}, \quad \lambda = \max\{\alpha > 0 | \alpha\mathbf{z} \in G\}. $$

$\pi_G(\mathbf{z})$ is called the projection of $\mathbf{z}$ on the upper boundary of $G$.

In polyblock outer approximation algorithm, the basic idea is to approximate the feasible set by constructing a sequence of polyblocks enclosing it, where at each iteration the polyblock is refined. The advantage is that the maximum of an non-decreasing function over a polyblock is attained at one of its proper vertices [31], so that is easy to find.

At each iteration $k$, we find a vertex $\mathbf{z}^k$ that maximizes $L(\mathbf{u})$ over the current polyblock approximation of $G$, i.e.,

$$ \mathbf{z}^k = \underset{\mathbf{u} \in \mathcal{T}^k}{\operatorname{argmax}} L(\mathbf{u}), $$

where $\mathcal{T}_k$ is the vertex set of the polyblock $\mathcal{P}^k$ constructed at the $k$th iteration. Therefore, $L(\mathbf{z}^k)$ is an upper bound to the optimal $L(\mathbf{u}^*)$ where

$$ \mathbf{u}^* = \underset{\mathbf{u} \in G}{\operatorname{argmax}} L(\mathbf{u}). $$

Moreover, $\mathbf{z}^k$ might not be feasible, therefore, it is projected on a point $\mathbf{x}^k$ on boundary of the feasible set $G$. Since, $\mathbf{x}^k$ is on the boundary of $G$, it is used to cut the cone $\mathcal{K}_{\mathbf{x}^k}^+$, where

$$ \mathcal{K}_{\mathbf{x}^k}^+ = \{\mathbf{y} \in \mathbb{R}_+^n | \mathbf{y} > \mathbf{x}^k\} $$

in order to refine the polyblock. Moreover, since $\mathbf{x}^k$ is a feasible point, $L(\mathbf{x}^k)$ is a lower bound to the optimal $L(\mathbf{u}^*)$. More precisely, the algorithm keeps memory of the “current best value ($CBV_k$)”, namely the largest among $L(\mathbf{x}^1), ..., L(\mathbf{x}^k)$. Clearly $CBV_k$ can be updated iteratively, as follows:

$$ CBV_k = \max\{CBV_{k-1}, L(\mathbf{x}^k)\}. $$
<figure><img src="image_7.png" /><figcaption>Fig. 6: Example of a normal set $G \subset \mathbb{R}^2$ and an illustration of the polyblock outer approximation algorithm. The red star denotes the optimal point.</figcaption></figure>

At each iteration the polyblock is refined, and hence the upper bound $L(\mathbf{z}^k)$ is decreased, while on the other hand the lower bound $CBV_k$ either remains constant or increases. The algorithm terminates when the lower and the upper bounds are within a desired precision threshold, thus certifying optimality, then the algorithm outputs $CBV_k$, together with the corresponding feasible point. Let us see in details the first iterations of the algorithm. The algorithm starts with the polyblock $\mathcal{P}^1 = [0, b]$, which contains the feasible set $G$ and corresponds to $\mathcal{T}^1 = \{b\}$, and hence $z^1 = b$. Then $\mathbf{z}^1$ is projected on $G$ along the line segment joining $0$ and $\mathbf{z}^1$ to obtain $\mathbf{x}^1$ on the boundary of $G$. Then, a smaller polyblock $\mathcal{P}^2$ enclosing $G$ is constructed based on $\mathcal{P}^1$ by replacing $\mathbf{z}^1$ with $n$ new vertices $\tilde{\mathcal{T}}_1 = \{\mathbf{v}^{1,1} \dots \mathbf{v}^{1,n}\}$. These new vertices $\mathbf{v}^{1,i}$ is computed as $\mathbf{v}^{1,i} = \mathbf{z}^1 + ((\mathbf{x})_i^1 - (\mathbf{z}^1)_i)e_i$, where $e_i \in \mathbb{R}^n$ is standard basis vector and $\mathbf{x}^1$. Thus, the vertex set of the newly constructed polyblock $\mathcal{P}^2$ is $\mathcal{T}_2 = \mathcal{T}_1 \setminus \{\mathbf{z}^1\} \cup \tilde{\mathcal{T}}_1$. Then, we choose the vertex which maximizes the objective function $L$ of the problem in (18), i.e., $\mathbf{z}^2 = \underset{z \in \mathcal{T}_2}{\operatorname{argmax}} L(\mathbf{z})$ for $\mathbf{z} \in \mathcal{T}_2$. This process is repeated successively to construct even smaller polyblocks such that

$$G \subseteq \dots \mathcal{P}_k \subseteq \dots \mathcal{P}_2 \subseteq \mathcal{P}_1.$$

Figure 6 illustrates this process for a normal set $G$ in $\mathbb{R}^2$. The algorithm terminates when $|L(\mathbf{z}^k) - L(\mathbf{x}^k)| \le \epsilon$, where $\epsilon > 0$ is the error tolerance specifying the error tolerance of the approximation. Under some mild conditions, [31], [32] show that $\mathbf{x}^k$ converges to a global optimal solution of the problem for $k \to \infty$. The Algorithm based on polyblock outer approximation algorithm has been summarised in the Algorithm 3. In particular, Step 4 searches for the vertex $\mathbf{z}^k$ which maximizes the objective function $L(\mathbf{z}^k)$ among the vertices of the polyblock $\mathcal{P}$. Step 5 computes the projection of $\mathbf{z}^k$ on the feasible set $G$ and also prevents jamming which will be discussed later. Steps 6 and 7 update the current best

value and current best feasible solution. Step 8 terminates the algorithm when $|L(\mathbf{z}^k) - L(\mathbf{x}^k)| \le \epsilon$, which implies that both the upper bound $L(\mathbf{z}^k)$ and the lower bound $L(\mathbf{x}^k))$ are within $\epsilon$ precision to the actual optimum $L(\mathbf{u}^*)$. Otherwise, Step 9 creates a smaller polyblock enclosing $G$, or, more precisely, constructs the vertex set of such polyblock (see details below). Step 10 fights against the growth of the vertex set, by removing the improper vertices (which are surely not optimal) and also removing the vertices $w$ with value $L(w)$ smaller than the current best value, as recommended in [33]. Finally, Step 13 checks if the algorithm can be terminated. The condition where the vertex set is empty corresponds to having removed all the vertices in Step 10, i.e., no vertex can have a better value than the current best one.

Now, we comment on some important steps in Algorithm 3.

* **Step 5:** The main purpose of this step is to compute the projection of $z_k$ on the upper boundary of $G$, $\pi_G(\mathbf{z}^k) = \mathbf{x}^k$ that can be done using a bisection search method summarised in Algorithm 4. Clearly, an approximation is used in bisection, therefore $\mathbf{x}^k$ ends up inside the $G$ within the approximation threshold.

Moreover, notice that the projection of $\mathbf{z}^k$, on the feasible set $G$ requires computing $L$ numerous times hence this step is computationally heavy. In our problem, it would mean running SIR-mobility model multiple times. Hence, obtaining convergence after a reasonably limited number of iterations would be desirable. However, in some cases, polyblock algorithm is known to have a very large number of iterations, due to a phenomenon known as jamming. When $\mathbf{z}^k$ is such that $z_i^k$ is close to zero for some $i$, then in order to converge from $\mathbf{z}^k$ to $(\mathbf{z}_1^k, \dots, \mathbf{z}_{i-1}^k, 0, \mathbf{z}_{i+1}^k, \dots, \mathbf{z}_n^k)$, a subsequence of such vertices are created and their difference in euclidean norm being very small. This phenomenon is termed as jamming. In particular, jamming occurs when the angle between the line from $0$ to the vertex $\mathbf{z}^k$ and one of the axes is steep. In order to prevent jamming, we used the method proposed in [34] to find the projection of $\mathbf{z}^k$ on $G$ alternatively.

Consider that $\mathbf{z}^k$ maximizes $\{L(\mathbf{z}) | \mathbf{z} \in \mathcal{T}_k\}$ at iteration $k$. Define

$$d = \min_{\{i: z_i^k > 0\}} (z_i^k), \qquad (20)$$

and define a new vertex $\mathbf{b}'$ as follows:

$$b_j' = \begin{cases} z_j^k - d & \text{when } z_j^k > 0 \\ 0 & \text{when } z_j^k = 0. \end{cases} \qquad (21)$$

In case $\mathbf{b}' \notin G$, we can use $\mathbf{b}'$ instead of $\mathbf{z}^k$, and construct the projection $\mathbf{x}^k = \pi_G(\mathbf{b}')$ instead of $\mathbf{x}^k = \pi_G(\mathbf{z}^k)$. Notice that also with this choice we obtain a point $\mathbf{x}^k$ on $\partial^+G$ in the box $[0, \mathbf{z}^k]$. Therefore, the next polyblock can be created by removing the region $(\mathbf{x}^k, \mathbf{z}^k]$. The main importance of this technique is that, when a vertex is very near to an axis, it pushes it to the axis in a single iteration, rather than slowly converging to zero.

In addition, we notice that if $d$ is large enough then $\mathbf{b}' \in G$, therefore, we set a threshold $d$ on the minimum
**Algorithm 3:** Category-dependent policy based on polyblock outer approximation method

**Input:** parameters for (2)-(7) and (9), threshold $\bar{I}$, initial time $T_0$, final time $T_f$, initial condition $N_i(T_0), I_i(T_0), S_i(T_0), R_i(T_0) \forall i$, tolerance $\epsilon, \bar{d}$, maximum number of iterations $n_{poly}$, no. of updates $p_m$ in the time interval $[T_0, T_f]$, and projection tolerance $\xi$ and projection maximum iterations $n_{it}$ (for Algorithm 4))

**Output:** $\epsilon$-optimal solution $\mathbf{u}^*$

1. Initialize the vertex set $\mathcal{T}_1 = \mathbf{1}_{p_m q}$. Set the Current Best Value($CBV_0$) = 0.

2. Set $k = 1$.

3. **REPEAT**

4. Select $\mathbf{z}^k = \operatorname{argmax}\{L(\mathbf{z})|\mathbf{z} \in \mathcal{T}_k\}$, where $L(\mathbf{z})$ is computed by running (2)-(7) and (9) for time for $k = T_0/\Delta t, \dots, T_f/\Delta t$.

5. Compute the minimum distance $d$ as in (20)

IF $d < \bar{d}$

Compute $\mathbf{b}'$ as in (21)

IF $\mathbf{b}' \in G$

Compute $\mathbf{x}^k = \pi_G(\mathbf{z}^k)$ using Algorithm 4.

ELSE

Compute $\mathbf{x}^k = \pi_G(\mathbf{b}')$ using Algorithm 4.

END IF

ELSE

Compute $\mathbf{x}^k = \pi_G(\mathbf{z}^k)$ using Algorithm 4.

END IF

6. IF $\mathbf{x}^k = \mathbf{z}^k$,

Set current best feasible solution $\bar{\mathbf{x}}^k = \mathbf{x}^k$
$CBV_k = L(\mathbf{z}^k)$

7. ELSE

IF $L(x) \ge CBV_{k-1}$
$\bar{\mathbf{x}}^k = \pi_G(\mathbf{z}^k)$ and $CBV_k = L(\pi_G(\mathbf{z}^k))$

ELSE

$\bar{\mathbf{x}}^k = \bar{\mathbf{x}}^{k-1}$ and $CBV_k = CBV_{k-1}$.

END IF

END IF

8. IF $|L(\mathbf{z}_k) - CBV_k| \le \epsilon$, GO TO Step 13 and EXIT

9. ELSE

- Create the new vertices $\mathbf{v}^j s$ from $\mathbf{z}_k$ and $x$ by following the rule:

$$ \mathbf{v}^{k,i} = \mathbf{z}^k + ((\mathbf{x})_i - (\mathbf{z}^k)_i)e_i, $$

for all $i$ such that $\mathbf{z}_i^k > \mathbf{x}_i^k$.

- $\mathcal{T}_{k+1} = \mathcal{T}_k \setminus \{\mathbf{z}^k\} \cup \mathbf{v}^{k,i}$, for all $i$, where $\mathbf{v}^{k,i}$ are the vectors computed in the previous step.

10. Remove all improper vertices from $\mathcal{T}_{k+1}$ and the vertices $w$ such that $L(\mathbf{w}) \le CBV_k + \epsilon$.

11. End IF

12. $k = k + 1$

13. Until: $\mathcal{T}_k = \emptyset$ OR $n_{poly}$ iterations.

14. $\mathbf{u}^* = \bar{\mathbf{x}}^k$

distance $d$ computed in (20) and only when $d < \bar{d}$ we check whether $\mathbf{b} \in G$, since this test is time-consuming. Furthermore, to prevent recurring jamming due to sufficiently small $x_i^k$ in the previous steps, it can be useful to round off the computed projection $\mathbf{x}^k$ if it is smaller than a set tolerance. This threshold should be chosen such that its effect on the objective function is negligible.

**Algorithm 4:** Bisection search to compute the projection $\mathbf{x} = \pi_G(\mathbf{z})$ in Step 5 of Algorithm 3

**Input:** $\mathbf{z}$, parameters for (2)-(7) and (9), threshold $\bar{I}$, initial time $T_0$, final time $T_f$, initial condition $N_i(T_0), I_i(T_0), S_i(T_0), R_i(T_0) \forall i$, maximum number of iterations $n_{it}$ and tolerance $\xi > 0$.

**Output:** $\mathbf{x} = a^*\mathbf{z}$ such that

$$ a^* = \operatorname{argmax}\{a^* > 0 | I_{\max}(a^*\mathbf{z}) \le \bar{I}\}. $$

1. Set $a_{\min} = 0$ and $a_{\max} = 1$.

2. **REPEAT**

3. Compute $\bar{a} = (a_{\min} + a_{\max})/2$,

4. Run (2)-(7) and (9) with $\mathbf{u}(k) = \bar{a}\mathbf{z}$ for $k = T_0/\Delta t, \dots, T_f/\Delta t$ and compute $I_{\max}(\bar{a}\mathbf{z})$ as in (13)

5. IF $I_{\max}(\bar{a}\mathbf{z}) \le \bar{I}$

6. $a_{\min} = \bar{a}$

7. ELSE, $a_{\max} = \bar{a}$.

8. till $(\bar{I} - I_{\max}(\bar{a}\mathbf{z})) \le \xi$ or $n_{it}$ iterations.

9. Set $a^* = a_{\min}$

10. Return $\mathbf{x} = a^*\mathbf{z}$

• **Step 9(b):** This step computes the vertex set $\mathcal{T}_{k+1}$ of the new polyblock. We know that $\mathbf{x}^k \in \partial^+G$ within the approximation threshold used in bisection search. Since $L$ is monotonic and $G$ is normal, this implies that all points in the cone $K_{\mathbf{x}^k}^+ = \{\mathbf{y} \in \mathbb{R}_+^n | \mathbf{y} > \mathbf{x}^k\}$ are unfeasible (except possibly a very small region, due to the approximation). Therefore, the new polyblock $\mathcal{P}_{k+1}$ is constructed by cutting the cone $K_{\mathbf{x}^k}^+$ out from $\mathcal{P}_k$. To create the new set of vertices $\mathcal{T}_{k+1}$, at first we compute $\{\mathbf{v}^{k,1} \dots \mathbf{v}^{k,n}\}$ as

$$ \mathbf{v}^{k,i} = \mathbf{z}^k + ((\mathbf{x})_i^k - (\mathbf{z}^k)_i)e_i, \quad (22) $$

for all $i$ such that $\mathbf{z}_i^k > \mathbf{x}_i^k$, where $\mathbf{x}^k = \pi_G(\mathbf{z}^k)$ and $e_i \in \mathbb{R}^n$ is standard basis vector. Notice that $\mathbf{v}^{k,i}$ is obtained by replacing the $i^{th}$ entry of $\mathbf{z}^k$ by the $i^{th}$ entry of $\mathbf{x}^k$. Finally, the set of vertices $\mathcal{T}_{k+1}$ of the polyblock $\mathcal{P}_{k+1}$ is computed as

$$ \mathcal{T}_{k+1} = \mathcal{T}_k \setminus \{\mathbf{z}^k\} \cup \mathbf{v}^{k,i}. $$

Although polyblock algorithm guarantees optimal solution, it is known to be slow in convergence [35] and like other non-convex optimization solver algorithms, it can be computationally demanding as the dimension increases. Therefore, in addition to enhancements, as mentioned earlier, we set a maximum number of iterations for the Algorithm 3. Even if the algorithm terminates after reaching maximum number of iterations, the obtained solution (although not optimal) is
ensured to be feasible, because of the projections on the boundary of the feasible set.

### C. Category-dependent receding-horizon policy

In this subsection, we obtain category-dependent capacity control policies, so as to apply restrictions based on the economic importance of different destination types. As in section IV-B, in this section also, the interval of interest is $[0, T]$ and the policy makers are assumed to update the policy at intervals of length $T_u < T$. For that, in principle, the problem (14) can also be solved directly using Algorithm 3 by taking $p_m = [T/T_u] = p$, and $T_0 = 0, T_f = T$, which corresponds to an optimization problem with pq variables. However, the convergence of Algorithm 3 can be very slow [35]. Therefore, to find a good, although possibly suboptimal, feasible solution of (14) in an efficient way, we adopt a receding horizon approach to obtain policies in which different control is applied on different destination categories, which can be updated on a regular interval such that we obtain a control $\mathbf{u}$ of the form (10). We compare the category-dependent policies with the category-free policies devised in Section IV.

**Algorithm 5:** Category-dependent receding-horizon policy, with update period $T_u$ and optimization horizon length $T_m$.

**Input:** $T, T_u, T_m$, Parameters for (2)-(7) and (9), initialization (16), threshold $\bar{I}$, projection tolerance $\xi$ and projection maximum iterations $n_{it}$, maximum polyblock iterations $n_{poly}$, and polyblock tolerance $\epsilon, \bar{d}$.

**Output:** $\mathbf{u}^* = \{\mu^{*1}, \dots, \mu^{*p}\}$, where $\mu^{*l} = \{\mu_h^l\}$ for all $h \in \{1, \dots, q\}$.

1. Initialize $\ell = 1$. Compute $p = [T/T_u]$.

2. **REPEAT**

3. Compute $T_0 = (\ell-1)T_u$ and $T_f = (\ell-1)T_u + T_m$

4. Apply Algorithm 3 with $p_m = 1$ over $[T_0, T_f]$, with threshold $\epsilon$, maximum number of iterations $n_{poly}$, and bisection tolerance $\xi$ and bisection maximum iterations $n_{it}$, and denote by $\mu^{*\ell}$ the obtained output.

5. Run (2)-(7) and (9) over the interval $[(\ell-1)T_u, \ell T_u]$, with control $\mu^{*\ell}$, so as to obtain $N_i(\ell T_u / \Delta t)$, $S_i(\ell T_u / \Delta t)$, $I_i(\ell T_u / \Delta t)$, $R_i(\ell T_u / \Delta t)$ for all $i$ (to be used as initialization for the next iteration).

6. $\ell = \ell + 1$

7. until $\ell = p$.

In this receding horizon approach, over each interval of optimization we choose to have a constant policy over time but allow different policies for different destination categories, which corresponds to an optimization problem with $q$ variables. In other words, we apply Algorithm 3 on the receding time intervals $[T_0, T_f]$ with the choice $p_m = 1$, where $T_0 = \ell T_u$ and $T_f = \ell T_u + T_m$ and $\ell = 0, \dots, p-1$, as described in Algorithm 5 and as depicted in Figure 3. The length of the optimization horizon is a parameter $T_m \in [T_u, T]$ to be designed.

## VI. SIMULATIONS

In this section, we present the simulations corresponding to different policies described in Sections IV-V implemented on the large scale urban mobility network of Grenoble. The parameters for the urban network are the ones described in Section III-E. Such simulations are also available in the on-line simulator GTL-Healthmob [36], which includes animations visualizing mobility and epidemics evolution on the map of Grenoble (as described in details in [37]).

### A. Category-free policy

1) **Category-free fixed policy:** Here, we show the simulation when implementing a fixed policy over the whole simulation horizon $T$ as explained in Algorithm 1. For the large-scale network of Grenoble area, Algorithm 1 with $\epsilon_f = 0.0015P$ and $n_{iter} = 20$, when applied on $[T_0, T_f] = [0, T]$ with $T = 15$ weeks, gives $u^* = 0.5625$. In Figures 7a-7c, the plots with $T_u = T = 15$ weeks correspond to the best fixed policy (since $T_u = T$), i.e., to applying the constant policy $u^* \mathbf{1}_q$ with $u^* = 0.5625$. In particular, Figure 7c depicts $I/P$, i.e., the proportion of infected population. The plot confirms that $I/P$ remains below the threshold $\bar{I}/P = 0.15$, as we expect, since by construction $u^*$ satisfies the constraint $I_{\max}(u^* \mathbf{1}_q) \le \bar{I}$.

2) **Category-free receding horizon policy:** Here, we show the simulation when implementing the category-free receding horizon policy in Grenoble as explained in Algorithm 2. For the large-scale network of Grenoble area, and with an horizon of interest $T = 15$ weeks, with $\epsilon_f = 0.0015P$ and $n_{iter} = 20$, Figure 7 illustrates the effect of different update periods $T_u$, with curves corresponding to $T_u = 1$ week, $T_u = 2$ weeks, $T_u = 3$ weeks and $T_u = T = 15$ weeks. For all curves, the optimization horizon is $T_m = 15$ weeks.

Figure 7a shows the optimal control policy, i.e., the piece-wise constant function $u^*(k)$ obtained from Algorithm 2, while Fig. 7a and Fig. 7c show the corresponding relative economic activity and proportion of infected people. We can see that all policies start with the same control $u^{*1}$, since they all start from the same initial infection status, and they all find the best constant $u$ for an optimization horizon of $T_m = 15$ weeks. They also remain equal until the infection peak (which is smaller than in the benchmark case, and below $\bar{I}$, as enforced by the constraint $I_{\max}(u\mathbf{1}_q) \le \bar{I}$). Then, after the peak, we can see the advantage of a smaller update period $T_u$, which allows a faster adaptation to the current epidemic conditions: as infections decrease, the capacity control $u$ is increased and lets the economic activity grow, and this can happen earlier for smaller values of $T_u$, leading to the best economic activity in the case of the smallest $T_u$. This suggests to use the smallest $T_u$ that can be socially acceptable, taking into account the fact that extremely short values of $T_u$, of the order e.g. of a day, would create confusion and lack of respect of the rules. As a side note, we can see the epidemic curve lasts longer in the optimized cases than in the benchmark case; this is a well-known effect of policies aiming at 'flattening the curve' so as to keep the number of infections below the rate that saturates health infrastructure.
<figure><img src="image_8.png" /><figcaption>**Fig. 7:** Category-free receding-horizon policy with different update period $T_u$: (a) Optimal control $u^*(k)$, (b) Relative economic activity $R_E(k)$, and (c) proportion of infected people ($I/P$), when Algorithm 2 is applied with $T_m = 15$ weeks and $T_u = 1, 2, 3, 15$ weeks. (c) also shows the benchmark case ($u(k)$ constantly equal to 1). The plot corresponding to $T_u = T_m = 15$ weeks is same as the fixed policy as discussed in IV-A.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 8:** Category-free receding-horizon policy with different optimization horizon $T_m$: (a) Optimal control $u^*(k)$, (b) Relative economic activity $R_E(k)$, and (c) proportion of infected people ($I/P$), when Algorithm 2 is applied with $T_u = 1$ week and $T_m = 1, 4, 7, 15$ weeks. (c) also shows the benchmark case ($u(k)$ constantly equal to 1). The plot corresponding to $T_u = T_m = 15$ weeks is same as the fixed policy as discussed in IV-A.</figcaption></figure>

Figure 8, instead, illustrates the effect of different optimization horizons $T_m$: Algorithm 2 is applied with $T_m = 1, 4, 7$ and $15$ weeks, while keeping a same small update period $T_u = 1$ week, which is chosen because a small value of $T_u$ is recommended, as discussed above.

It can be seen in Figure 8b that performance improves when $T_u$ decreases. The best solution is obtained when $T_m = 1$, which corresponds to a U-shape control (see Figure 8a): full capacity ($u = 1$) initially, when there are few infections, then strong limitations (small $u$) near the infection peak, followed by gradual reopenings, progressively augmenting $u$ while keeping the infection level near the threshold, until full reopening at the end of the infection wave. The fact that the short-sighted policy outperforms the long-sighted ones might seem somewhat counter-intuitive, but it can be clarified by recalling that Algorithm 2 looks for the best fixed policy over each receding-horizon interval. Thus, Algorithm 2 with a large $T_m$ enforces an overly-cautious small $u^{*1}$, since despite the small initial number of infections it chooses $u^{*1}$ sufficiently small to flatten the infections peak that would happen in the coming $T_m$ weeks. With small $T_m$, instead, Algorithm 2 quickly adapts the control input to current infection level,

leading to best performance, both because of a larger initial activity, and because the infection peak happens earlier, allowing for earlier partial reopenings. However, notice that this performance comes at the price of a more aggressive control policy, that uses some lower values of $u$, although for short time, as it is clearly illustrated in Figure 8a. The social acceptability of such restrictions should be taken into account in the policy design.

From these simulations in this section, we can draw the recommendation to design policies using Algorithm 2 with the smallest values of the update period that is $T_u = 1$ and of the optimization horizon $T_m = 1$, within the margin of social acceptability.

## B. Category-dependent policy

Here, we show the simulation when implementing a receding horizon policy on Grenoble as explained in Algorithm 5. Algorithm 5 is similar in spirit to Algorithm 2: both have optimization over receding intervals of length $T_m$, where for simplicity optimization only explores policies which are constant in such time interval, and re-optimize every $T_u$ weeks, while they differ in the fact that Algorithm 5 exploits the
opportunity to apply different controls to different categories. In Section IV-B, it has been seen that the approach in Algorithm 2 gives the best results when $T_u = 1$ and $T_m = 1$. For this reason, in this section we perform simulations using Algorithm 5 with $T_m = 1$ and $T_u = 1$, to obtain results with this category-dependent policy, and we compare them with the results obtained in Section IV with a category-free policy and same parameters $T_u$ and $T_m$.

For the large-scale network of Grenoble and with horizon of interest $T = 15$ weeks, Figure 9 compares the effect of category-dependent control policy, computed using Algorithm 5, with $\epsilon = 0.002$, $n_{poly} = 150$, $n_{it} = 20$, and $\bar{d} = 0.01$ for $T_m = 1$ and $T_u = 1$ with that of the category-free control policy computed using Algorithm 2 with $\epsilon_f = \bar{I}/100$ and $n_{iter} = 20$ in Section IV.

While running Algorithm 5 in our simulation, in most weeks (i.e., in most iterations of Algorithm 5) Algorithm 3 terminated by having reached a distance between lower and upper bound that is below the threshold $\epsilon$. Only in weeks 4, 5 and 6 Algorithm 3 reached the maximum number of iterations $n_{poly}$ and hence gives a suboptimal result. Recall that even in this case the obtained control policy is feasible. Despite this suboptimality, we can see in Figure 9c that the category-dependent policy obtained with Algorithm 5 significantly outperforms the category-free policy, thus confirming the interest of allowing for category-dependent controls.

Figures 9a and 9b show the piecewise-constant function $u^*(k)$ obtained from Algorithm 5, and piecewise-constant function $u^*(k)$ obtained from Algorithm 2, respectively, while Figures 9c and 9d show the corresponding relative economic activity and proportion of infected people. We can see that both the policies start with the same control $1_5$, since they are computed considering the same initial epidemics situation with few infections, and hence they both foresee that full opening is possible, without causing too many infections over the short optimization horizon $T_m = 1$. However, from the second week,

some partial closure becomes necessary in order to keep the infections below the given threshold. Here, the difference between the category-free and the category-dependent approach becomes apparent: by using a smaller control (i.e., heavier restrictions) on the categories with lower socio-economic weight $\chi_h$, the category-dependent approach allows a higher control (i.e., less restrictions) on the categories with higher socio-economic weight, whose contribution has a major impact on the socio-economic cost, while still ensuring that infections remain below the threshold. For this reason, the category-dependent approach can outperform the category-free one; indeed, Figure 9c shows that, after 15 weeks, category-dependent policy is able to retain about 92 percent of the socio-economic of the Benchmark case, in comparison to 85 percent in case of category-free policy.

From these simulations, we can draw the recommendations to design category-dependent policies using Algorithm 5, whenever the weights $\chi_h$ (which reflect socio-political priorities) are quite unbalanced.

## VII. DISCUSSION

Our method is based on a compartmental epidemic model which is an extension of a SIR model. Compartmental models are well-established and useful tools, but also have some well-known limitations, see e.g. [38]. Our model shares with standard SIR model the following limitations. First, it assumes homogeneity of infection parameters and contact probability within the population (or rather within each sub-population corresponding to a location, in our case), and hence it cannot capture the effect of possible ‘super-spreaders’, i.e., people who can cause a number of infections much higher than the average. Moreover, it considers time-invariant infection and recovery rates, disregarding seasonal effects as well as mutants. Also, it disregards births, deaths, and movements of population to and from outside the area at study: the total population size is assumed to be constant. These simplifications make this model inappropriate to predict the onset of multiple consecutive epidemic waves. The model should not be used on very long intervals of time (such as an year in the running example), but rather be used on intervals of time that might correspond to a single wave, then re-starting new computations with updated parameters.

Being SIR-based, our model also makes the assumption that recovered people are immune, while such immunity might be only partial, and even the more so when new mutants replace the original virus. In some cases it might be more appropriate to investigate different compartmental models, such as SEIRS [39], or to adjust the initialization of the number of susceptible and recovered people when starting a SIR-based simulation upon the onset of a new dominant mutant. Another well-known limitation of all compartmental models is that they are deterministic models, which only provide an estimated number of infected people, without any indication of its uncertainty. The uncertainties can be modelled either by stochastic individual-based epidemic models, or within the framework of compartmental models by letting the main parameters be uncertain. The latter approach has been taken e.g. in [20], [21], together with stochastic MPC techniques, where Monte-Carlo methods are used to ensure that the constraints are satisfied with the desired probability (when parameters appear in the constraints, as in [21]) or to approximate the expected cost (when the parameters appear in the cost, as in [20], the expectation being non-trivial to obtain due to non-linearities). Taming the complexity of a stochastic MPC formulation of our optimization problem with a fine-grained mobility model remains an open problem, beyond the scope of this paper.

## VIII. CONCLUSION

In this paper, we have designed mobility control policies to mitigate epidemics in a city-wide network.

We have designed policies for capacity control at destinations, so as to maximize the economic activity under the constraint that the number of infections remains below a desired threshold, so as to avoid saturating health-care system.

At first, we proposed algorithms to design such policies, by focusing on the case where the same control is applied to all destination categories that can be updated on a regular
<figure><img src="image_10.png" /><figcaption>**Fig. 9:** Receding-horizon policy with optimization horizon $T_m = 1$ and $T_u = 1$: (a) Category-dependent control $\mathbf{u}^*(k)$ obtained with Algorithm 5, (b) Category-free control $\mathbf{u}^*(k)$, obtained with Algorithm 2, (c) Relative economic activity $R_E(k)$, and (d) proportion of infected people ($I/P$), when Algorithm 5 and Algorithm 2 are applied with $T_u = 1$ week and $T_m = 1$ week.</figcaption></figure>

time interval. For that, a receding horizon approach which reoptimizes the cost function based on the epidemic situation over receding time windows.

Then, we address the more general optimization problem, where different controls can be applied to different destination categories and at different times. Our solution technique uses polyblock outer approximation method along with a receding horizon approach. Finally, we implemented these techniques on large-scale network of Grenoble metropolitan area to demonstrate that such minimally-invasive mobility restriction policies can be devised for a real-life large scale urban area. Our simulations highlight the importance to design category-dependent policies, whenever the weights (which reflect socio-political priorities of destinations) are quite unbalanced. From the simulations, we also obtain some design recommendations concerning the two main parameters of our algorithms: the update period and the optimization horizon should be small, so as to better capture variations in time of the epidemic evolution, but on the other hand they should not be too small, to comply with social acceptability.

## REFERENCES

[1] E. H. Dyvik, “Impact of the coronavirus pandemic on the global economy - statistics and facts,” https://www.statista.com/topics/6139/covid-19-impact-on-the-global-economy/#editorsPicks, last accessed 17 July 2024.

[2] H. W. Hethcote, “The mathematics of infectious diseases,” *SIAM Review*, vol. 42, no. 4, pp. 599–653, 2000.

[3] A. Vespignani *et al.*, “Modelling COVID-19,” *Nature Reviews Physics*, vol. 2, pp. 279–281, 2020.

[4] N. Perra, “Non-pharmaceutical interventions during the COVID-19 pandemic: A review,” *Physics Reports*, vol. 913, pp. 1–52, 2021.

[5] T. A. Perkins and G. Espana, “Optimal control of the COVID-19 pandemic with non-pharmaceutical interventions,” *Bull Math Biology*, vol. 82(9):118, 2020.

[6] J. Köhler, L. Schwenkel, A. Koch, J. Berberich, P. Pauli, and F. Allgöwer, “Robust and optimal predictive control of the COVID-19 outbreak,” *Annual Reviews in Control*, vol. 51, pp. 525–539, 2021.

[7] A. Borri, P. Palumbo, F. Papa, and C. Possieri, “Optimizing restrictions in epidemics via piecewise time-varying SIRD models: Application to the COVID-19 Italian emergency,” *European Journal of Control*, vol. 75, p. 100902, 2024.

[8] G. Giordano *et al.*, “Modelling the COVID-19 epidemic and implementation of population-wide interventions in Italy,” *Nature Medicine*, vol. 26, pp. 855–860, 2020.

[9] D. Acemoglu, V. Chernozhukov, I. Werning, and M. D. Whinston, “Optimal targeted lockdowns in a multigroup SIR model,” *American Economic Review: Insights*, vol. 3, no. 4, p. 487–502, 2021.

[10] S. T. McQuade *et al.*, “Control of COVID-19 outbreak using an extended SEIR model,” *Mathematical Models and Methods in Applied Sciences*, vol. 31, no. 12, pp. 2399–2424, 2021.

[11] A. Kasis, S. Timotheou, N. Monshizadeh, and M. Polycarpou, “Optimal intervention strategies to mitigate the COVID-19 pandemic effects,” *Scientific Reports*, vol. 12, 2022.

[12] P. Nouvellet *et al.*, “Reduction in mobility and COVID-19 transmission,” *Nature Communications*, 2021.

[13] A. Smolyak, G. Bonaccorsi, A. Flori, F. Pammolli, and S. Havlin, “Effects of mobility restrictions during COVID19 in Italy,” *Scientific Reports*, pp. 2045–2322, 2021.

[14] M. Gosgens *et al.*, “Trade-offs between mobility restrictions and transmission of SARS-CoV-2,” *J R Soc Interface*, 2021.

[15] S. Song, Z. Zong, Y. Li, X. Liu, and Y. Yu, “Reinforced epidemic control: Saving both lives and economy,” 2020. [Online]. Available: https://arxiv.org/abs/2008.01257
[16] L. Sattenspiel and K. Dietz, "A structured epidemic model incorporating geographic mobility among regions," *Mathematical Biosciences*, vol. 128, no. 1, pp. 71–91, 1995.

[17] C. Poletto, M. Tizzoni, and V. Colizza, "Human mobility and time spent at destination: impact on spatial epidemic spreading." *J Theor Biol.*, 2013.

[18] R. Carli, G. Cavone, N. Epicoco, P. Scarabaggio, and M. Dotoli, "Model predictive control to mitigate the COVID-19 outbreak in a multi-region scenario," *Annual Reviews in Control*, vol. 50, pp. 373–393, 2020.

[19] W. Hu, Y. Shi, C. Chen, and Z. Chen, "Optimal strategic pandemic control: human mobility and travel restriction," *Math Biosci Eng.*, 2021.

[20] R. Dutta, S. N. Gomes, D. Kalise, and L. Pacchiardi, "Using mobility data in the design of optimal lockdown strategies for the COVID-19 pandemic," *PLOS Computational Biology*, vol. 17, 08 2021.

[21] P. Scarabaggio, R. Carli, G. Cavone, N. Epicoco, and M. Dotoli, "Nonpharmaceutical stochastic optimal control strategies to mitigate the covid-19 spread," *IEEE Transactions on Automation Science and Engineering*, vol. 19, no. 2, pp. 560–575, 2022.

[22] M. U. Niazi, C. Canudas-de-Wit, A. Kibangou, and P.-A. Bliman, "Optimal control of urban human mobility for epidemic mitigation," *CDC 2021 - 60th IEEE Conference on Decision and Control, Dec 2021, Austin, United States*, 2021.

[23] U. Pratap, C. Canudas-de-Wit, and F. Garin, "Where, when and how people move in large-scale urban networks: the Grenoble saga," Feb. 2022, working paper or preprint. [Online]. Available: https://hal.archives-ouvertes.fr/hal-03554612

[24] A. Suryanto and I. Darti, "On the nonstandard numerical discretization of SIR epidemic model with a saturated incidence rate and vaccination," *AIMS Mathematics*, 6(1): 141–155, 2021.

[25] D. Acemoglu, A. Makhdoumi, A. Malekian, and A. Ozdaglar, "Testing, voluntary social distancing and the spread of an infection," *Operations Research*, no. 27483, July 2020.

[26] M. U. B. Niazi, A. Kibangou, C. Canudas-de-Wit, D. Nikitin, L. Tumash, and P.-A. Bliman, "Modeling and control of epidemics through testing policies," *Annual Reviews in Control*, vol. 52, pp. 554–572, 2021.

[27] S. Galmiche et al., "Exposures associated with SARS-CoV-2 infection in France: A nationwide online case-control study," *The Lancet Regional Health - Europe*, vol. 7, p. 100148, 2021.

[28] "Our world in data," https://ourworldindata.org/coronavirus/country/france, last accessed 17 April 2024.

[29] "Worldometer," https://www.worldometers.info/coronavirus/country/france/, last accessed 17 April 2024.

[30] J. L. Geertz, J. M. Greene, C. H. Sanchez-Tapia, and E. D. Sontag, "A novel COVID-19 epidemiological model with explicit susceptible and asymptomatic isolation compartments reveals unexpected consequences of timing social distancing." *Journal of Theoretical Biology*, vol. 510, p. 110539, 2021.

[31] H. Tuy, "Monotonic optimization: Problems and solution approaches," *SIAM Journal on Optimization*, vol. 11, no. 2, pp. 464–494, 2000.

[32] H. Tuy, F. Al-Khayyal, and P. T. Thach, *Monotonic Optimization: Branch and Cut Methods*. Boston, MA: Springer US, 2005, pp. 39–78.

[33] Y. J. Zhang, L. Qian, and J. Huang, "Monotonic optimization in communication and networking systems," *Foundation and trends in networking*, vol. 7, no. 1(2012), pp. 1–75, 2013.

[34] M.-S. Cheon, F. A. Al-Khayyal, and S. Ahmed, "Global optimization of monotonic programs: applications in polynomial and stochastic programming," 2005. [Online]. Available: https://smartech.gatech.edu/bitstream/handle/1853/6938/Ceon_MyunSeok_200505_phd.pdf;sequence=1

[35] H. Chour, E. A. Jorswieck, F. Bader, Y. Nasser, and O. Bazzi, "Global optimal resource allocation for efficient FD-D2D enabled cellular network," *IEEE Access*, vol. 7, pp. 59690–59707, 2019.

[36] "GTL-Healthmob simulation platform," http://gtlville.inrialpes.fr/covid-19, last accessed 17 April 2024.

[37] U. Pratap, L. Senique, and C. Canudas-de Wit, "GTL-Healthmob: Simulation platform for urban mobility and epidemic control," in *2022 - gènes journées des Démonstrateurs en Automatique*, Angers, France, Jun. 2022, pp. 1–11. [Online]. Available: https://hal.archives-ouvertes.fr/hal-03674156

[38] A. Traulsen, C. S. Gokhale, S. Shah, and H. Uecker, "The Covid-19 pandemic: Basic insights from basic mathematical models," *Nova Acta Leopoldina*, vol. 2022, no. 3, pp. 1–19, 2022.

[39] O. N. Bjørnstad, K. Shea, M. Krzywinski, and N. Altman, "The SEIRS model for infectious disease dynamics," *Nature Methods*, vol. 17, pp. 55–57, 2020.