Issue 
Math. Model. Nat. Phenom.
Volume 15, 2020
Coronavirus: Scientific insights and societal aspects



Article Number  48  
Number of page(s)  21  
DOI  https://doi.org/10.1051/mmnp/2020035  
Published online  14 October 2020 
Optimal control techniques based on infection age for the study of the COVID19 epidemic^{*}
^{1}
INRIASaclay and Centre de Mathématiques Appliquées, Ecole Polytechnique,
91128
Palaiseau, France.
^{2}
CIFASISCONICETUNR, Ocampo y Esmeralda, S2000EZP,
Rosario, Argentina.
^{**} Corresponding authors: frederic.bonnans@inria.fr; gianatti@cifasisconicet.gov.ar
Received:
18
May
2020
Accepted:
3
September
2020
We propose a model for the COVID19 epidemic where the population is partitioned into classes corresponding to ages (that remain constant during the epidemic). The main feature is to take into account the infection age of the infected population. This allows to better simulate the infection propagation that crucially depend on the infection age. We discuss how to estimate the coefficients from data available in the future, and introduce a confinement variable as control. The cost function is a compromise between a confinement term, the hospitalization peak and the death toll. Our numerical experiments allow to evaluate the interest of confinement varying with age classes.
Mathematics Subject Classification: 49M25 / 92D25 / 92D30 / 35Q93
Key words: COVID19 / epidemic modeling / confinement / agestructured systems / optimal control
© The authors. Published by EDP Sciences, 2020
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
The COVID19 epidemic appears for the first time in Wuhan, China at the end of 2019. It quickly spread throughout the world, and on March 12, 2020 it was declared a pandemic by the World Health Organization (WHO). So far there are no vaccines and it is not clear if there is any effective treatment to combat it. That is why the governments of the different countries began to implement containment measures in the population, trying to control the disease. Since these actions have a great impact on the daily life of the population, one of the big questions now is knowing how long to carry out these measures.
The recent COVID19 epidemic has already been modelized in e.g. [8–10]. These papers use variants and extensions of the SIR (SusceptibleInfectedRecovered) model based on ODEs (ordinary differential equations). SIR is the standard approach for the modelization of epidemics and takes its origin in [7]. Of course it is difficult to model confinement. The simplest way is to assume that it boils down to a multiplicative factor coefficient with value in [0, 1] in the formula for new infections. That is, the dynamics in the SIR setting would be $$\begin{array}{lll}\dot{S}(t)\hfill & =\hfill & \delta (1u(t))S(t)I(t)\hfill \\ \dot{I}(t)\hfill & =\hfill & \delta (1u(t))S(t)I(t)(\eta +\mu )I(t)\hfill \\ \dot{R}(t)\hfill & =\hfill & \mu I(t)\hfill \end{array}$$(1.1)
with δ > 0 infection coefficient, η > 0 death rate, and μ > 0 recovery rate. Observe that, here as in the sequel, our models do not take into account birth, as well as death for other reasons than the epidemic itself. We recover the classical SIR dynamics in the case of an uncontrolled system, i.e., when u(t) = 0 for all t. Total confinement is out of reach and therefore we may assume that we have the control constraint u(t) ≤ u^{M}, with u^{M} ∈ (0, 1). Of course it is not easy to specify a confinement cost, but we might take $$c(u):={p}_{u}{\displaystyle {\int}_{0}^{T}\phi}(u(t))\text{d}t$$(1.2)
with p_{u} ≥ 0 penalty factor, and for example φ(s) = s (as we did in our numerical tests) or φ quadratic.
Our proposal is to distinguish, among infected people, those who are hospitalized and those who are not, and to take into account the limited duration of the infection. The latter can be achieved by associating with the infected population an “infection age" variable. That also allows us to differentiate between the incubation and infectious phases of the disease. Furthermore, since the severity of this COVID19 seems to depend on the age of the infected population, it may be useful to take into account the population age, in addition to the infection age. The cost function will minimize a combination of confinement cost, value of hospitalization peak, and death toll.
Concerning the optimal control of agestructured population, see [4] (without diffusion), and [1] (with diffusion). The optimal control of age models with state constraints is discussed in [3]. About optimal control application to epidemiological models, see e.g. [11]. About the related question of models involving delays, see e.g. [13]. On optimal control approaches for COVID19, see [5]. During the revision process, we noticed the appearance of [12] that also propose an optimal control model with infection age. The main difference with respect to our model is that we consider the peak value in the objective function, which leads to a problem with state constraints. Let us also mention a mean field game approach to the COVID19 epidemic [6].
The paper is organized as follows. In Section 2, we present a continuous time (and infection age) optimal control problem. The corresponding discrete time problem is discussed in Section 3. The data analysis (how to estimate thecoefficients of the problem) as well as the initial conditions are the subject of Section 4. Finally, in Section 5, we discuss various numerical results. The paper ends with a conclusion giving a perspective of possible extensions.
Our open access code, based on the optimal control toolbox Bocop http://www.bocop.org, can be downloaded from http://www.cmap.polytechnique.fr/~bonnans/covid19/covid19.html or https://www.fceia.unr.edu.ar/~justina/covid19/.
2 Continuous time optimal control problem
2.1 An infection age structured model
We have two kinds of patients, the exposed ones, who are asymptomatic or with mild symptoms, and those who need hospitalization. Let z(t, a, b) (resp. h(t, a, b)) denote the number of exposed (resp. hospitalized) patients of age a ∈ [0, A] at time t, that have been infected for duration b ∈ [0, B]. Here A and B are positive constants. Let y(t, a) (resp. ȳ (t, a)) denote the number of nonimmunized (resp. immunized) people of age a at time t. Since in practice the age variable a does not significantly vary during the epidemic, we assume that individuals keep the same age over time. The control variable u(t, a) ∈ [0, 1] is the confinement policy (no confinement if u(t, a) = 0, total confinement if u(t, a) = 1). We considerthe following dynamics $$\{\begin{array}{l}\dot{y}(t,a)=\delta (a)(1u(t,a))Z(t)y(t,a)\hfill \\ y(0,a)={y}_{0}(a),\hfill \\ ({z}_{t}+{z}_{b})(t,a,b)=\nu (a,b)z(t,a,b)\hfill \\ z(0,a,b)={z}_{0}(a,b)>0,\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}z(t,a,0)=\delta (a)(1u(t,a))Z(t)y(t,a)\hfill \\ ({h}_{t}+{h}_{b})(t,a,b)=\nu (a,b)z(t,a,b)\left(\eta (a,b)+\gamma (a,b)E(t)\right)h(t,a,b)\hfill \\ h(0,a,b)={h}_{0}(a,b)\ge 0,\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}h(t,a,0)=0\hfill \\ \stackrel{.}{\overline{y}}(t,a)=z(t,a,B)+h(t,a,B)\hfill \\ \overline{y}(0,a)=0.\hfill \end{array}$$(2.1)
Here, δ(a) ≥ 0 is the “transmission coefficient”, Z(t) is the “global transmission factor", modelled as $$Z(t)={\displaystyle {\int}_{0}^{A}{\displaystyle {\int}_{0}^{B}e}}(a,b)z(t,a,b)\text{d}a\text{d}b,$$(2.2)
with e(a, b) ≥ 0 “transmission factor" (we assume that hospitalized patients do not contribute to the transmission of the disease). Also, ν(a, b) ≥ 0 is the hospitalization coefficient, and E(t) ∈ [0, 1] is the hospital saturation estimate, given by $$E(t):=\frac{{(H(t)C)}_{+}}{H(t)+C};\text{\hspace{1em}}H(t):={\displaystyle {\int}_{0}^{A}{\displaystyle {\int}_{0}^{B}h}}(t,a,b)\text{d}a\text{d}b,$$(2.3)
where C > 0 is the nominal capacity. We assume that only hospitalized patients die, with death rate $$d(t,a,b):=\eta (a,b)+\gamma (a,b)E(t),$$(2.4)
where η(a, b) ≥ 0 is the minimal death rate and γ(a, b) sensitivity of the death rate w.r.t. the hospital saturation estimate.
Remark 2.1.
In our numerical tests we will assume that there exists an incubation phase [0, n_{0}) during which transmission and hospitalization do not occur. Additionally, in the infectious phase for which b ∈ [n_{0}, B], we assume that e(a, b) = 1 and that ν(a, b), η(a, b) and γ(a, b) have constant values denoted by $\overline{\nu}(a)$, $\overline{\eta}(a)$ and $\overline{\gamma}(a)$ resp. Then, the coefficients take the form: $$e(a,b)={1}_{[{n}_{0},B]}(b),\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\nu (a,b)={1}_{[{n}_{0},B]}(b)\overline{\nu}(a),\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\eta (a,b)={1}_{[{n}_{0},B]}(b)\overline{\eta}(a),\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{and}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\gamma (a,b)={1}_{[{n}_{0},B]}(b)\overline{\gamma}(a).$$(2.5)
Consequently, h(t, a, b) = 0 whenever b ∈ [0, n_{0}] and $$Z(t)={\displaystyle {\int}_{0}^{A}{\displaystyle {\int}_{{n}_{0}}^{B}z}}(t,a,b)\text{d}b\text{d}a.$$(2.6)
So, in this case, the model has four parameters for each value of a: δ(a), $\overline{\nu}(a)$, $\overline{\eta}(a)$ and $\overline{\gamma}(a)$.
2.2 Optimal control problem
We next introduce the cost function to be minimized. Given the horizon T, the death toll for age a is equal to the difference between the initial and final population for this age: $$\begin{array}{ll}D(a):=\hfill & y(0,a)+\overline{y}(0,a)+{\displaystyle {\int}_{0}^{B}z}(0,a,b)+h(0,a,b)\text{d}b\hfill \\ \text{\hspace{0.17em}}\hfill & \left(y(T,a)+\overline{y}(T,a)+{\displaystyle {\int}_{0}^{B}\left(z(T,a,b)+h(T,a,b)\right)}\text{d}b\right).\hfill \end{array}$$(2.7)
The total death toll is ${D}_{T}:={\displaystyle {\int}_{0}^{A}D}(a)\text{d}a.$ Denote by $$c(u):={\displaystyle {\int}_{0}^{T}{\displaystyle {\int}_{0}^{A}{c}_{e}}}(a)u(t,a)\text{d}a\text{d}t$$(2.8)
the weighted confinement cost, with weights c_{e}(a) ≥ 0. The cost to be minimized is a combination of several costs. Its expression is $$J(M,u,{D}_{T}):={p}_{M}M+{p}_{u}c(u)+{p}_{D}{D}_{T},$$(2.9)
where the penalty coefficients p_{M}, p_{u} and p_{D} are nonnegative. Here M is an additional variable, called the peak value, subject to $$H(t)\le M,\text{\hspace{1em}}\forall t\in [0,T].$$(2.10)
Indeed, if p_{M} is nonzero, theoptimal M will be the minimal hospital capacity needed to take care in an optimal way of all patients, which justifies the peak value terminology.The ability to minimize this amount has been recognized of paramount importance.
The second term in the cost function involves the weighted confinement effort. The associated economic cost is difficult to evaluate; this question being out of the scope of the present paper, we simply take it linear with somewhat arbitrary weights, having in mind the obvious fact that the confinement of retired people has a smaller economic impact than for the rest of the population, although psychological aspects should also be taken into account.
Finally, remember that the death rate depends on the saturation of hospital capacities; so the confinement has an impact on the death toll too (assuming the coefficient γ to be nonzero and the hospital capacity to be saturated).
In addition, in our model we include confinement duration constraints for each age, more precisely: $${\int}_{0}^{T}u}(t,a)\text{d}t\le M(a),\text{\hspace{1em}}\text{fora}.\text{a}.\$\text{a}\in (\text{0},\text{A})\$.$$(2.11)
3 Discretization of the control problem
In this section, we present the discrete model associated to the previous formulation. The main objective of this section is to explain how to write the discrete problem so that an optimal control solver (based on an ODE model) can address it (see Appendix A for more details).
For the age a ∈ [0, A] we take discrete values {1, …, n_{a}}. As both time unit, time step and discretization of the infection age we choose one day.
The horizon of the problem will be $N\in \mathbb{N}$. Let n_{0}, n_{b}, with n_{0} < n_{b} be the duration of the incubation phase and of the infection, respectively.
The time index is $k\in \mathbb{N}$, starting from value 0. The infection age is j ∈{1, …, n_{b}}. The initial conditions are $$\begin{array}{cccc}{y}_{a}^{0}& =& {y}_{0}(a)& a=1,\dots ,{n}_{a},\\ {z}_{a,j}^{0}& =& {z}_{0}(a,j)& a=1,\dots ,{n}_{a},\text{\hspace{0.17em}}\text{\hspace{0.17em}}j=1,\dots ,{n}_{b},\\ {h}_{a,j}^{0}& =& {h}_{0}(a,j)& a=1,\dots ,{n}_{a},\text{\hspace{0.17em}}\text{\hspace{0.17em}}j=1,\dots ,{n}_{b}.\end{array}$$(3.1)
The global transmission factor and the size of the hospitalized population are respectively $${Z}^{k}={\sum}_{a=1}^{{n}_{a}}{\sum}_{j=1}^{{n}_{b}}{e}_{a,j}{z}_{a,j}^{k},\text{\hspace{1em}}{H}^{k}={\sum}_{a=1}^{{n}_{a}}{\sum}_{j=1}^{{n}_{b}}{h}_{a,j}^{k}.$$(3.2)
Also E^{k} ∈ [0, 1] is a discrete counterpart to the hospital saturation estimate E(t), given by $${E}^{k}:=\frac{{({H}^{k}C)}_{+}}{{H}^{k}+C}.$$(3.3)
Willing to keep conservation laws, we may write the dynamics as follows: for all a = 1, …, n_{a}, $$\begin{array}{cccc}{y}_{a}^{k+1}& =& \left(1{\delta}_{a}(1{u}_{a}^{k}){Z}^{k}\right){y}_{a}^{k}& \text{\hspace{0.17em}}\\ {z}_{a,1}^{k+1}& =& {\delta}_{a}(1{u}_{a}^{k}){Z}^{k}{y}_{a}^{k}& \text{\hspace{0.17em}}\\ {z}_{a,j}^{k+1}& =& (1{\nu}_{a,j1}){z}_{a,j1}^{k}& j=2,\dots ,{n}_{b}\\ {h}_{a,1}^{k+1}& =& 0& \text{\hspace{0.17em}}\\ {h}_{a,j}^{k+1}& =& {\nu}_{a,j1}{z}_{a,j1}^{k}+(1{\eta}_{a,j1}{\gamma}_{a,j1}{E}^{k}){h}_{a,j1}^{k}& j=2,\dots ,{n}_{b}\\ {\overline{y}}_{a}^{k+1}& =& {\overline{y}}_{a}^{k}+{z}_{a,{n}_{b}}^{k}+{h}_{a,{n}_{b}}^{k}.& \text{\hspace{0.17em}}\end{array}$$(3.4)
The death toll per age and day, and per day, are resp. $${D}_{a}^{k}:={\sum}_{j=1}^{{n}_{b}1}({\eta}_{a,j}+{\gamma}_{a,j}{E}^{k}){h}_{a,j}^{k};\text{\hspace{1em}}{D}^{k}:={\sum}_{a=1}^{{n}_{a}}{D}_{a}^{k}.$$(3.5)
The parameters of the dynamics are:
δ_{a} :  transmission rate 
ν_{a,j}:  hospitalization rate 
η_{a,j}:  minimal death rate 
γ_{a,j}:  sensitivity of the death rate w.r.t. hospital saturation. 
Remark 3.1.
The state variable ${y}_{a}^{k}$ remains nonnegative iff
$${\delta}_{a}(1{u}_{a}^{k}){Z}^{k}\le 1.$$(3.6)
In the case of normalized, (i.e., summing to 1) initial population, we have that Z^{k} ≤1. So, the previous nonnegativity condition holds whenever δ_{a} ≤ 1. In our numerical test, it may happen that δ_{a} ≥ 1, but ${y}_{a}^{k}$ always remained positive.
On the other hand, in our tests η_{a,j} + γ_{a,j} ≤ 1. Since E^{k} ≤ 1 by definition, it follows that ${h}_{a,j}^{k}$ always remains nonnegative.
We have the following constraints: maximal achievable confinement, peak value for the hospital capacity, and cumulated confinement effort. Namely, for some ${u}_{a}^{M}\in [0,1]$, M ≥ 0, and M_{a} ≥ 0: $$\begin{array}{l}0\le {u}_{a}^{k}\le {u}_{a}^{M},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}k=0,\dots ,N,\text{\hspace{0.17em}}a\in \{1,{n}_{a}\}\hfill \\ {H}^{k}\le M,\text{\hspace{1em}}k=0,\dots ,N,\hfill \\ {\sum}_{k=0}^{N1}{u}_{a}^{k}\le {M}_{a},\text{\hspace{0.17em}}a\in \{1,{n}_{a}\}.\hfill \end{array}$$(3.7)
The expression of the cost is $${J}_{d}(M,u,{D}_{T}):={p}_{M}M+{p}_{u}{c}_{u}+{p}_{D}{D}_{T},$$(3.8)
with c_{u} weighted confinement cost, D_{T} total death toll: $${c}_{u}:={\sum}_{k=0}^{N1}{\sum}_{a=1}^{{n}_{a}}{c}_{e,a}{u}_{a}^{k};\text{\hspace{1em}}{D}_{T}={\sum}_{k=0}^{N1}{D}^{k},$$(3.9)
and p_{u}, p_{M} and p_{D} are nonnegative penalty coefficients.
Due to the incubation phase, the value of control for the last n_{0} days has no influence on the hospitalizations and therefore, neither on the death toll nor on the hospital capacity constraint.
4 Estimation of coefficients and initial state
4.1 Incubation and infectious phases
As we mentioned in Remark 2.1, for our numerical tests we assume that there exists an incubation phase (days 1 to n_{0} − 1 of infection) where no transmission and hospitalizations occur. So, $${e}_{a,j}=0\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{for}\text{\hspace{0.17em}}\text{\hspace{0.17em}}1\le j\le {n}_{0}1,\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{e}_{a,j}=1\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{for}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{n}_{0}\le j\le {n}_{b},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{and}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\nu}_{a,j}=0\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{for}\text{\hspace{0.17em}}\text{\hspace{0.17em}}1\le j\le {n}_{0}1.$$(4.1)
The above assumptions together with the initial condition ${h}_{a,1}^{k}=0$ implies that $${Z}_{a}^{k}={\sum}_{a=1}^{{n}_{a}}{\sum}_{j={n}_{0}}^{{n}_{b}}{z}_{a,j}^{k},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{and}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{h}_{a,j}^{k}=0,\text{\hspace{0.17em}}\text{whenever}\text{\hspace{0.17em}}1\le j\le {n}_{0}.$$(4.2)
During the infectious phase (j varying from n_{0} + 1 to n_{b}) we may assume that the coefficients ν_{a,j}, for j ≥ n_{0} η_{a,j} and γ_{a,j}, for j ≥ n_{0} + 1, do not depend on j and therefore are equal to some positive constants ${\overline{\nu}}_{a}$, ${\overline{\eta}}_{a}$, ${\overline{\gamma}}_{a}$.
We assume that we know for each age the total proportion $\widehat{\nu}$ of transfers from z to h, that satisfies $${\widehat{\nu}}_{a}={\overline{\nu}}_{a}+{\overline{\nu}}_{a}(1{\overline{\nu}}_{a})+\cdots +{\overline{\nu}}_{a}{(1{\overline{\nu}}_{a})}^{{n}_{b}{n}_{0}1}={\overline{\nu}}_{a}{\sum}_{j=0}^{{n}_{b}{n}_{0}1}{(1{\overline{\nu}}_{a})}^{j}=1{(1{\overline{\nu}}_{a})}^{{n}_{b}{n}_{0}},$$(4.3)
so that we can estimate $${\overline{\nu}}_{a}=1{(1{\widehat{\nu}}_{a})}^{1/({n}_{b}{n}_{0})}.$$(4.4)
Similarly, we may assume that we know for each age the proportion ${\widehat{\eta}}_{a}$ of hospitalized people that will die, when the hospital capacity is not saturated. Since ${h}_{a,j}^{k}=0$, for j = 1, …, n_{0}, we have by similar arguments $${\overline{\eta}}_{a}=1{(1{\widehat{\eta}}_{a})}^{1/({n}_{b}{n}_{0}1)}.$$(4.5)
As expected, $\overline{\nu}$ and $\overline{\eta}$ take values in [0, 1].
When the total hospital capacity is exceeded we assume that it can be measured so that ${\overline{\gamma}}_{a}$ itself can beevaluated.
4.2 Initial condition
Assuming that we know a growth factor g for the infected population given a time period of p time steps, then e^{pλ} = g, so that $$\lambda =(\mathrm{log}g)/p.$$(4.6)
We assume that the initial time corresponds to an early phase of the epidemic, and therefore that the h and ȳ variables are equal to zero, and (since the order of magnitude of y is almost constant) search for the z variable a solution of the form $${z}_{a,j}^{k}={c}_{a,j}{\text{e}}^{\lambda (kj)}.$$(4.7)
We can express c_{a,j} as function of c_{a,1}: $${c}_{a,j}=\{\begin{array}{ll}{c}_{a,1}\hfill & \text{if}j\le {n}_{0}\hfill \\ {c}_{a,1}{(1{\overline{\nu}}_{a})}^{j{n}_{0}}\hfill & \text{otherwise}\text{.}\hfill \end{array}$$(4.8)
So the total and infectious population of z for each age are resp. $${\overline{Z}}_{a}^{k}:={\sum}_{j=1}^{{n}_{b}}{z}_{a,j}^{k}={c}_{a,1}{\overline{\chi}}_{a}{\text{e}}^{\lambda k};\text{\hspace{1em}}{\overline{\chi}}_{a}:=\left({\sum}_{j=1}^{{n}_{0}1}{\text{e}}^{\lambda j}+{\sum}_{j={n}_{0}}^{{n}_{b}}{(1{\overline{\nu}}_{a})}^{b{n}_{0}}{\text{e}}^{\lambda j}\right)$$(4.9) $${Z}_{a}^{k}:={\sum}_{j={n}_{0}}^{{n}_{b}}{z}_{a,j}^{k}={c}_{a,1}{\chi}_{a}{\text{e}}^{\lambda k};\text{\hspace{1em}}{\chi}_{a}:={\sum}_{j={n}_{0}}^{{n}_{b}}{(1{\overline{\nu}}_{a})}^{j{n}_{0}}{\text{e}}^{\lambda j}.$$(4.10)
In the absence of confinement at initial time, estimating c_{a,1} with the previous relation, we deduce the value of δ from: $${c}_{a,1}={z}_{a,1}^{1}={\delta}_{a}{Z}^{0}{y}_{a}^{0}.$$(4.11)
As noted in Remark 3.1, in order that y^{k} remains nonnegative, we need that $${\delta}_{a}(1{u}_{a}^{k}){Z}^{k}\le 1,\text{\hspace{1em}}\text{forall}\hspace{0.17em}k=0,\dots ,\text{N}\hspace{0.17em}.$$(4.12)
Remark 4.1
We have seen that, since the initial time corresponds to an early stage of the epidemic, we can assume that h^{0} = 0 and ȳ^{0} = 0, and additionally, we can assume that there is no saturation of the hospital capacity. We could consider that we have statistical estimates of the “initial parameters" λ, $\widehat{\nu}$, $\widehat{\eta}$, y^{0}, Z^{0}.
Then we can deduce $\overline{\nu}$ and $\overline{\eta}$ from (4.4) and (4.5), respectively. The value of χ and $\overline{\chi}$ follow. We can deduce c_{a,1} from (4.10) and the values of z^{0} follow with (4.7). Finally, we deduce δ from (4.11).
For the estimate of the γ parameters, specific methods are needed, which are beyond the scope of this paper.
5 Numerical tests
5.1 Choice of parameters and initial conditions
We consider the available data for France up to April 6, 2020. From Santé publique France and data.gouv.fr we obtain the following information:
Total confirmed cases:  74 390 
Total hospitalized cases:  54 036 
Total confirmed deaths:  8911 
From these data we can estimate the values of $\widehat{\nu}$ and $\widehat{\eta}$ that represent the proportion of infected population in needs of hospitalization and the proportion of deaths with respect to the hospitalizations, respectively. $$\widehat{\nu}=0.726,\text{\hspace{1em}}\text{and}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\widehat{\eta}=\mathrm{0.165.}$$(5.1)
Furthermore, from the abovementioned sources, we can observe that in the beginning of the epidemic, the growth factor per week is between 2.5 and 3, we can considerthat $$\lambda \in [0.13,0.15].$$(5.2)
In the French demographic institute website INED we find detailed information on deaths in hospitals, taking into account the age of the population. We deduce that the strong population y_{1} can be considered as the population up to 59 years and the weak population y_{2} as the population over 60 years. In this case, we have
Total population  %  Total confirmed deaths  %  

strong  49 213 491  73.4  416  6.5 
weak  17 850 240  26.6  6021  93.5 
total  67 063 703  100  6037  100 
Since they only consider deaths in hospitals, the number of deaths in the two previous tables is not the same. From the last table we can calculate the relative mortality rate as $${\eta}_{1}^{r}=\frac{416}{49\text{\hspace{0.17em}}213\text{\hspace{0.17em}}491}=8.45e6,\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{and}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\eta}_{2}^{r}=\frac{6021}{17\text{\hspace{0.17em}}850\text{\hspace{0.17em}}240}=3.4\text{e}4,$$(5.3)
which implies that the ratio of this rate is $$\frac{{\eta}_{1}^{r}}{{\eta}_{2}^{r}}=\mathrm{0.0248.}$$(5.4)
In order to obtain the values of ${\widehat{\eta}}_{1}$ and ${\widehat{\eta}}_{2}$, we assume that ratio is the same, i.e., ${\widehat{\eta}}_{1}/{\widehat{\eta}}_{2}=0.0248$. Normalizing the population, the global death rate $\widehat{\eta}$ is equal to $$\widehat{\eta}={\widehat{\eta}}_{1}{y}_{1}+{\widehat{\eta}}_{2}{y}_{2}=(0.0248{\widehat{\eta}}_{2})0.734+{\widehat{\eta}}_{2}\mathrm{0.266.}$$(5.5)
By (5.1) we obtain, $${\widehat{\eta}}_{1}=0.014,\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\widehat{\eta}}_{2}=\mathrm{0.58.}$$(5.6)
Since we do not have detailed information about the age of hospitalized people, we assume: $$\widehat{\nu}={\widehat{\nu}}_{2}={\widehat{\nu}}_{1}.$$(5.7)
Finally, from data.gouv.fr we know the total of confirmed cases per day, but again we do not have the information by age. For this reason, and since in our formulation, the discount term in the first dynamics equation is depending on the sum of infected population (Z^{k}), we consider that if ${\overline{Z}}^{0}$ is the initial confirmed cases, then we can define the (normalized) initial confirmed cases by age, as $${\overline{Z}}_{1}^{0}={\overline{Z}}^{0}\frac{0.734}{{y}^{0}},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\overline{Z}}_{2}^{0}={\overline{Z}}^{0}\frac{0.266}{{y}^{0}},$$(5.8)
where y^{0} = 6 706 3703 is the total population of France.
Since we have taken the same values of $\overline{\nu}$ for the two populations, and the initial infected population proportional to the total population, the values for δ are identical. The initial infected population is computed as it was explained, taking the data from data.gouv.fr.
In our numerical tests, the first day is March 15, 2020. For France this was actually the first day of confinement, with 6633 confirmed cases. We consider the following parameters:
${y}_{1}^{0}$  Initial strong population  0.734 
${y}_{2}^{0}$  Initial weak population  0.266 
${\overline{Z}}_{1}^{0}$  Initial infected strong population  7.26e − 5 
${\overline{Z}}_{2}^{0}$  Initial infected weak population  2.63e − 5 
λ  Growth coefficient given by (4.6)  0.13 
${\widehat{\nu}}_{1}$  Total proportion of transfers from z_{1} to h_{1}  0.726 
${\widehat{\nu}}_{2}$  Total proportion of transfers from z_{2} to h_{2}  0.726 
${\widehat{\eta}}_{1}$  Strong minimal death rate  0.014 
${\widehat{\eta}}_{2}$  Weak minimal death rate  0.58 
We choose an incubation duration of 6 days, and a total infection duration as 14 days. From the above table, following the methodology described in Remark 4.1, we obtain: $$\begin{array}{lll}{\delta}_{1}=1.656,\hfill & \text{\hspace{0.17em}}{\overline{\nu}}_{1}=0.149412,\hfill & \text{\hspace{0.17em}}{\overline{\eta}}_{1}=0.002012\hfill \\ {\delta}_{2}=1.656,\hfill & \text{\hspace{0.17em}}{\overline{\nu}}_{2}=0.149412,\hfill & \text{\hspace{0.17em}}{\overline{\eta}}_{2}=0.116557\hfill \end{array}$$(5.9)
In addition, we assume: $${\overline{\gamma}}_{a,j}={\overline{\eta}}_{a,j},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\forall a=1,2,\text{\hspace{0.17em}}1\le j\le {n}_{b};\text{\hspace{0.17em}}\text{\hspace{0.17em}}C=\mathrm{0.005.}$$(5.10)
5.2 Numerical results in the case of an age independent control
In these first four tests, we assume the control variables to be independent of the age a, i.e., ${u}_{1}^{k}={u}_{2}^{k}$ for all k, that the upper bound for the control is u_{M} = 0.75, and that the confinement cost, given by (3.9), satisfies c_{e,1} + c_{e,2} = 1. We present different scenarios, taking different penalty terms in the objective function. We do not consider any cumulated control constraints.
In each of the figures below, we show four graphics. The first one (upper left), is the optimal confinement, the second one (upper right) is the cumulative number of deaths for the strong and weak population, the third one (bottom left) is the number of people in need of hospitalization per day and the last one (bottom right) shows the evolution of the susceptible strong and weak population.
Test 1: p_{M} = 1, p_{u} = 1, p_{D} = 1, Figure 1.
Test 2: p_{M} = 1, p_{u} = 0.000001, p_{D} = 0, Figure 2
Test 3: p_{M} = 1e − 5, p_{u} = 0, p_{D} = 1, Figure 3
Test 4: p_{M} = 1, p_{u} = 0.0005, p_{D} = 1, Figure 4
Figure 1 Test 1: p_{M} = 1, p_{u} = 1, p_{D} = 1. The optimal control is 0. Without confinement, the death toll is greater than 10% and the peak value is very large, as expected. 
Figure 2 Test 2: p_{M} = 1, p_{u} = 0.000001, p_{D} = 0. In this case, the main objective is to reduce the peak value. We can observe the delay, between the confinement decision and the response in the hospitalized population. Indeed, when the confinement stops, the size of hospitalized population remains constant for several days. The death toll is greater than 10%. 
Figure 3 Test 3: p_{M} = 1e − 5, p_{u} = 0, p_{D} = 1. The main objective is to reduce the death toll, without considering the confinement cost. Therefore, the optimal policy assumes themaximum confinement value until the last days. Since the saturation of hospital capacity impacts in the death rate, as a consequence the peak value is reduced. 
Figure 4 Test 4: p_{M} = 1, p_{u} = 0.0005, p_{D} = 1. In this case,we want to reduce the peak value and the death toll, taking into account the economic impact. 
In the following table we show the death toll for strong, weak and total population, follow by $I(u):={\sum}_{k=0}^{N}{u}^{k}$ and the peak value M.
Test  D_{1,T}  D_{2,T}  D_{T}  I(u)  M 

1  0.0088192  0.116966  0.1257852  0  0.27665 
2  0.00816993  0.111278  0.11944793  54.49  0.0698543 
3  0.00660044  0.0906906  0.09729104  102.375  0.0693109 
4  0.00688871  0.094924  0.10181271  65.9702  0.0715698 
Discussion of the numerical results for an age independent control:
To analyze and compare the complete behavior of the epidemic when different confinement policies are applied, we assume that the maximum for this control is 0.75. Of course, this is the reason why in the tests we obtain that almost the entire population is immunized in the period considered. On the contrary, at the moment in France and in many countries around the world, containment measures are stronger and after more than amonth of confinement we are far from the end of the epidemic. If these policies continue, the time needed to reach this level of immunization will be even longer.
From the above tests, we can observe that if there is no confinement, which is the case of Test 1, then the health system will collapse andthe death toll will be very large, as expected.
It is interesting to compare Test 2 and Test 3. In Test 3, the main objective is to reduce the death toll, and we observe that in comparison with Test 2, we obtain more or less the same peak, some days after, but with a small death toll. The difference between them, is that in the first case, when the hospitalized population attains its maximum, the confinement stops, whereas in Test 3, the confinement is maximal until the end. That means that reducing the peak is not enough to reduce the final number of deaths.
Finally, as expected, in Test 4 since the control penalty is greater than those previously considered, the confinement is reduced, obtaining a greater peak value and death toll, but even though the death toll is smaller than in Test 2. We can also observe thatin Test 4 since the optimal policy is 0 during the first 20 days, the peak value is obtained several weeks before than in Test 2 and Test 3.
As we mentioned in the introduction, it is very difficult to weigh economic losses against the number of lives lost.
5.3 Numerical results with age dependent controls
Taking into account these results, it is clear that most deaths occur in the weak population, so it is natural to address the issue of considering two different levels of confinement, depending on the age of the population. We take 0.75 as upper bound for the two confinement controls. In this case, the cost function is given by (3.9). For the numerical tests, we assume that c_{e,1} = 0.734 and c_{e,2} = 0.133. With this coefficients we take into account the proportion of the population that is weak and strong and we also assume that confinement for the elderly has less impact on the economy than confinement for the young. We use the same initial conditions and parameters as in the previous tests. For Test 5 and Test 6, we do not consider any cumulated control constraints.
Test 5: p_{M} = 1, p_{u} = 0.000001, p_{D} = 0, Figure 5
Figure 5 Test 5: p_{M} = 1, p_{u} = 0.000001, p_{D} = 0. Since hospitalizations rates ${\overline{\nu}}_{1}$ and ${\overline{\nu}}_{2}$ are assumed equals, to reduce the amount of hospitalizations that occur at the same time, it makes sense that both populations are confined for almost the same period. In this case, once the peak is reached, the confinement stops. We still observe a large number of deaths, greater than 10%. 
Test 6: p_{M} = 1, p_{u} = 0.0005, p_{D} = 1, Figure 6
Figure 6 Test 6: p_{M} = 1, p_{u} = 0.0005, p_{D} = 1. As in Test 4, the aim is to balance between reducing peak value, death toll and confinement cost. We observe an important difference between optimal confinement for strong and weak populations. Of course, this occurs due to the difference between the mortality rates for both populations. 
In the last test, we take the same values of Test 6, but we include cumulated control constraints (3.7) with M_{1} = 25 and M_{2} = 45.
Test 7: p_{M} = 1, p_{u} = 0.0005, p_{D} = 1, Figure 7
Figure 7 Test 7: p_{M} = 1, p_{u} = 0.0005, p_{D} = 1, with the confinement lengths M_{1} = 25 and M_{2} = 45. The objective of this test is to show how the optimal policy changes when we add constraints in the confinement lengths. Taking the delay into account, we see that the end of confinement of the strong population corresponds to the end of the hospitalization peak. 
In the following table, we summarized the results of the tests with differentiated confinement:
Test  D_{1,T}  D_{2,T}  D_{T}  I(u_{1})  I(u_{2})  M 

5  0.00817212  0.109833  0.11800512  67.7524  73.6568  0.0693917 
6  0.00824508  0.0983062  0.10655128  35.8425  62.25  0.0729996 
7  0.00837618  0.10037  0.10874618  25  45  0.0805676 
Discussion of the numerical results for age dependent controls:
We analyze the difference between taking a confinement policy dependent or independent of the age. Test 5 is similar to Test 2, in both cases we achieve to reduce the peak value, but we obtain a very high death toll. In Test 5, there is no significant difference between confinement for the strong and weak population, and the behavior is similar to Test 2, once the peak is reached, theconfinement stops (taking into account the delay effects).
The comparison between Test 4 and Test 6 is maybe the most interesting one. With similar penalty terms, we obtain almost the same death toll and peak value. The main difference is that the strong population has a considerably shorter confinement in Test 6, while the weak population has almost the same as in Test 4. This confirms that performing differentiated confinement is perhaps the best response when seeking to minimize the economic impact of confinement, while at the same time reducing the number of deaths.
From the available data, it is evident that the majority of deaths belong to the weakest group, so it is natural to think that if the objective is to reduce the number of final deaths, the confinement for this group should be larger than for the other. When the sole objective is to reduce the peak, we do not observe this behavior because it is more important to reduce the number of infected people than the total number of deaths. Since in our model we take into account the saturation of the hospital system in the mortality rate of hospitalized people, when we minimize the number of deaths, we also obtain a reduction in the peak.
5.4 Two cases of specific interest
As we mentioned in the previous discussions, the objective of the previous tests was to analyze different scenarios that show that our numerical methods work well and that all the results obtained are consistent with our model. Now we will try to analyze cases closer to thecurrent data.
5.4.1 A quasiperiodic control
In order to show an example which may be closer to the actual available data so far in France, we consider the case when after March 15, the transmission rate δ is reduced to the half. This can be interpreted as a result of a change in the population behavior, in addition to the confinement policy itself, such as washing hands, social distancing, and the use of masks. We allow the confinement policy to take values between 0 and 0.9. We choose p_{M} = 12, p_{u} = 0.005, p_{D} = 1 and assume that the hospital capacity is C = 0.01. See Figure 8.
After several weeks of confinement, hospitalizations per day begin to decrease. As expected, since the total number of immunized people is still too small, the epidemic is far from the end and when the confinement policy is relaxed, a rebound can be observed. This quasiperiodic feature of confinement could be a good strategy not to overwhelm the health system.
Given that the containment policy is very strong, after 140 days we did not observe great changes in the number of immunized people. However, the death toll is close to 0.0011, which is equivalent to approximately 73,000 people. Again, the reason for this large number of deaths is due to the hospitalization and mortality rates considered.
Figure 8 Quasiperiodic case: p_{M} = 12, p_{u} = 0.005, p_{D} = 1, C = 0.01. 
5.4.2 Longterm behavior
In this case, we take the parameters values of the above test, in particular the transmission rate δ = 0.825 and the control takes values in [0, 0.9]. The penalty terms are p_{M} = 12, p_{u} = 0.005, p_{D} = 1, which implies that the main objective is to reduce the peak value, taking into account the death toll, but the costs associated with confinement are quite high. In this test we increase the horizon to 260 days.
We observe that some days after the peak of hospitalizations is achieved, the control is very high, and then it begins to slowly decrease to zero, which is when the epidemic is close to ending, see Figure 9. When the situation of the country does not allow to take a long and strong confinement, this may be the more accurate strategy. Of course, due to penalties considered, the death toll is still too high.
Figure 9 Longterm behavior: p_{M} = 12, p_{u} = 0.005, p_{D} = 1, C = 0.01, T = 260. 
5.5 Comments on the numerical results
A critical question is to check if data are available for evaluating the coefficients. Exchanges on the model with epidemic experts would be of course very useful.
We want to point out that all the results obtained in the tests are illustrative. It is well known that the available data on confirmed cases could be far from the actual number, due to asymptomatic cases. A higher number of infected people would imply that hospitalization and mortality rates are lower and, therefore, the death toll as well as the maximum number of hospitalizations will be lower. In fact, we took the available data on March 15, 2020, when the death toll was almost 16%. This is why,in our tests, when the epidemic is almost over, the number of deaths is so high.
In addition, the other parameter that influences the number of deaths and the peak value, is the hospitalization rate $\overline{\nu}$, which at the beginning was close to 70%. Given that there are apparently many cases that are not serious, at a later stage of the epidemic when there is more demand in hospitals, some of these cases could be treated at home, obtaining a lower rate. Furthermore, we have no information on the length of hospitalization.
6 Conclusion
6.1 Possible extensions
We want to highlight several possible extensions or modifications of our model. Of course, such changes require further information, but we may have more accurate data in the near future. For instance, if a vaccine is available we can introduce a new control v(t, a) which is the per capita rate of vaccination. The resulting modification in the dynamics would be: $$\begin{array}{lll}\dot{y}(t,a)\hfill & =\hfill & \delta (a)(1u(t,a))Z(t)y(t,a)v(t,a)y(t,a),\hfill \\ \dot{\overline{y}}(t,a)\hfill & =\hfill & z(t,a,B)+h(t,a,B)+v(t,a)y(t,a).\hfill \end{array}$$(6.1)
It is natural to think that v depends on the age of the population, and the limitation of M_{V} vaccines per day would read as $${\int}_{0}^{A}v}(t,a)y(t,a)\text{d}a\le {M}_{V},\text{\hspace{1em}}\text{for\hspace{0.17em}a}.\text{a}.\hspace{0.17em}\text{t}\in (\text{0},\text{T}).$$(6.2)
Here, we neglect any delay for the vaccination to become effective. Also, if data was available, it may be useful to split the infected population between the confirmed cases (taking advantage of a testing policy) and nonconfirmed (among then some asymptomatic cases). Keeping the variable z for the nonconfirmed cases, we could introduce a new state w(t, a, b) representing the confirmed cases. Assuming the latter to be isolated, only the z population contributes to the global transmission factor Z(t). We obtain the following dynamics: $$\{\begin{array}{l}\dot{y}(t,a)=\delta (a)(1u(t,a))Z(t)y(t,a)\hfill \\ ({z}_{t}+{z}_{b})(t,a,b)=\mu (t,a,b)z(t,a,b)\nu (a,b)z(t,a,b)\hfill \\ ({w}_{t}+{w}_{b})(t,a,b)=\mu (t,a,b)z(t,a,b)\nu (a,b)w(t,a,b)\hfill \\ ({h}_{t}+{h}_{b})(t,a,b)=\nu (a,b)(z(t,a,b)+w(t,a,b))\left(\eta (a,b)+\gamma (a,b)E(t)\right)h(t,a,b)\hfill \\ \stackrel{.}{\overline{y}}(t,a)=z(t,a,B)+w(t,a,B)+h(t,a,B).\hfill \end{array}$$(6.3)
Here μ(t, a, b) ≥ 0 is another control, which represents the detection rate. A possible associated cost for the testing policy would be, for some positive coefficients c_{1} and c_{2}: $${\int}_{0}^{T}{\displaystyle {\int}_{0}^{A}{\displaystyle {\int}_{0}^{B}\mu}}}(t,a,b)({c}_{1}z(t,a,b)+{c}_{2}y(t,a,b))\text{d}b\text{d}a\text{d}t.$$(6.4)
In addition, we can consider the hospitalization age to be a significant parameter. Denoting it by c ∈ [0, C], and assuming that h depends on (t, a, c), we get the following dynamics for the hospitalized population: $$\{\begin{array}{l}h(t,a,0)={\displaystyle {\int}_{0}^{b}\nu}(a,b)(z(t,a,b)+w(t,a,b))\text{d}b\hfill \\ ({h}_{t}+{h}_{c})(t,a,c)=\left(\eta (a,c)+\gamma (a,c)E(t)\right)h(t,a,c),\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}c\in [0,C].\hfill \end{array}$$(6.5)
It is also possible to introduce a recovery rate, that could depend on the hospitalization age, and would contribute to the dynamics of ȳ.
As we can see there are many possible adaptations of our model. On the other hand, one should avoid to excessively increase the number of states and parameters. A model with too many states is difficult to optimize; one with too many parameters makes identifiability uneasy.
6.2 Final remarks
Population models with infection age need (after infection age discretization) many states. Nevertheless, we could obtain a numerical convergence of our solvers in various cases. Sometimes we needed to solve the problem first with a small horizon (e.g. 3 months) and then to increase progressively the horizon, in order to obtain the convergence.
The main advantage of our model is the agestructure. Despite the limited accuracy of the parameters and the obtained data, in our tests we recovered the particular behavior of the epidemic with respect to the incubation phase and almost fixed duration, at least for nonsevere cases. We can see the delay in the response when a containment policy is taken, which is crucial in this type of disease where the infected population grows exponentially.
Our model is quite flexible and allows to take for instance variable coefficients during or after the incubation phase. In future studies one should of course take advantages of clinical studies in order to take into account a more accurate dynamical model, including the response to drugs, as well as the extensions just previously discussed.
It is of interest to observe how is the containment variable in the “long term” case of Section 5.4.2. The containmentis initially absent, then close to its maximal value at a time close to the saturation of the hospital capacity, then smoothly decreases and goes back to 0. This (happily) seems to be in accordance with the expected social behavior (to accept more easily strict measures at the beginning).
Besides, this study raises interesting questions related to the analysis of the optimality conditions, especially since the model involves state constraints. For this purpose the techniques of the recent paper [3] could be useful.
Appendix A General time steps
We used the opensource optimal control solver Bocop ([2], http://www.bocop.org). Since our software solves optimal control problems described by an ODE, the present problem does not enter in its format. However, note that we can rewrite the discrete system (3.4) in this way: $$\begin{array}{ccc}{y}_{a}^{k+1}{y}_{a}^{k}& =& {\delta}_{a}(1{u}_{a}^{k}){Z}^{k}{y}_{a}^{k}\\ {z}_{a,1}^{k+1}{z}_{a,1}^{k}& =& {z}_{a,1}^{k}+{\delta}_{a}(1{u}_{a}^{k}){Z}^{k}{y}_{a}^{k}\\ {z}_{a,j}^{k+1}{z}_{a,j}^{k}& =& {z}_{a,j}^{k}+(1{\nu}_{a,j1}){z}_{a,j1}^{k}\\ {h}_{a,j}^{k+1}{h}_{a,j}^{k}& =& {h}_{a,j}^{k}+{\nu}_{a,j1}{z}_{a,j1}^{k}+(1{\eta}_{a,j1}{\gamma}_{a,j1}{E}^{k}){h}_{a,j1}^{k}\\ {\overline{y}}_{a}^{k+1}{\overline{y}}_{a}^{k}& =& {z}_{a,{n}_{b}}^{k}+{h}_{a,{n}_{b}}^{k}.\end{array}$$(A.1)
Taking advantage of the fact that the time step is 1, we see that the discrete system (3.4) can be interpreted as the explicit Euler discretization scheme, for dynamics corresponding to the r.h.s. of the above array.
While a time step of one day is natural, in the case of large horizons it may be useful to conduct experiments with time steps of several days. So, consider a time step (equal to the step for the infection age) Δ > 0. We need ${{n}^{\prime}}_{b}:={n}_{b}/\Delta $ to be integer. Clearly the following system is a consistent approximation of (2.1): $$\begin{array}{ccc}\frac{{y}_{a}^{k+1}{y}_{a}^{k}}{\Delta}& =& {\delta}_{a}(1{u}_{a}^{k}){Z}^{k}{y}_{a}^{k}\\ \frac{{z}_{a,1}^{k+1}{z}_{a,1}^{k}}{\Delta}& =& \frac{{z}_{a,1}^{k}}{\Delta}+{\delta}_{a}(1{u}_{a}^{k}){Z}^{k}{y}_{a}^{k}\\ \frac{{z}_{a,j+1}^{k+1}{z}_{a,j+1}^{k}}{\Delta}& =& \frac{{z}_{a,j}^{k}{z}_{a,j+1}^{k}}{\Delta}{\nu}_{a,j1}{z}_{a,j}^{k}\\ \frac{{h}_{a,j+1}^{k+1}{h}_{a,j+1}^{k}}{\Delta}& =& \frac{{h}_{a,j}^{k}{h}_{a,j+1}^{k}}{\Delta}+{\nu}_{a,j}{z}_{a,j}^{k}({\eta}_{a,j}+{\gamma}_{a,j}{E}^{k}){h}_{a,j}^{k}\\ \frac{{\overline{y}}_{a}^{k+1}{\overline{y}}_{a}^{k}}{\Delta}& =& {z}_{a,{{n}^{\prime}}_{b}}^{k}+{h}_{a,{{n}^{\prime}}_{b}}^{k}.\end{array}$$(A.2)
In addition the above scheme is conservative, and can be interpreted as the explicit Euler discretization scheme, for dynamics corresponding to the r.h.s.
References
 B. Ainseba, S. Aniţa and M. Langlais, Optimal control for a nonlinear agestructured population dynamics model. Electron. J. Differ. Equ. 9 (2003) 28. [Google Scholar]
 J. Bonnans, D. Giorgi, V. Grélard, B. Heymann, S. Maindrault, P. Martinon, O. Tissot and J. Liu, Bocop – a collection of examples. Technical report, INRIA (2017). [Google Scholar]
 J.F. Bonnans and J. Gianatti, Optimal control of state constrained agestructured problems. SIAM J. Control Optim. 58 (2020) 2206–2235. [Google Scholar]
 M. Brokate, Pontryagin’s principle for control problems in agedependent population dynamics. J. Math. Biol. 23 (1985) 75–101. [CrossRef] [MathSciNet] [PubMed] [Google Scholar]
 R. DjidjouDemasse, Y. Michalakis, M. Choisy, M.T. Sofonea and S. Alizon, Optimal COVID19 epidemic control until vaccine deployment. Preprint MedRxiv: 20049189v2 (2020). [Google Scholar]
 R. Elie, E. Hubert and G. Turinici, Contact rate epidemic control of covid19: an equilibrium view. MMNP 15 (2020) 35. [EDP Sciences] [Google Scholar]
 W.O. Kermack, A.G. McKendrick and G.T. Walker, A contribution to the mathematical theory of epidemics. Proc. Roy. Soc. Lond. A 115 (1927) 700–721. [CrossRef] [Google Scholar]
 Q. Lin, S. Zhao, D. Gao, Y. Lou, S. Yang, S.S. Musa, M.H. Wang, Y. Cai, W. Wang, L. Yang and D. He. A conceptual model for the coronavirusdisease 2019 (covid19) outbreak in Wuhan, China with individual reaction and governmental action. Int. J. Infect. Dis. 93 (2020) 211–216. [CrossRef] [PubMed] [Google Scholar]
 Z. Liu, P. Magal, O. Seydi and G. Webb, Predicting the cumulative number of cases for the COVID19 epidemic in China from early data. Preprint MedRxiv: 20034314v1 (2020). [Google Scholar]
 Z. Liu, P. Magal, O. Seydi and G. Webb, Understanding unreported cases in the COVID19 epidemic outbreak in Wuhan, China, and the importance of major public health interventions. Biology 9 (2020) 50. [Google Scholar]
 H. Maurer and M. do Rosario de Pinho, Optimal control of epidemiological SEIR models with L^{1}objectives and controlstate constraints. Pac. J. Optim. 12 (2016) 415–436. [Google Scholar]
 Q. Richard, S. Alizon, M. Choisy, M.T. Sofonea and R. DjidjouDemasse, Agestructured nonpharmaceutical interventions for optimal control of covid19 epidemic. Preprint MedRxiv: 20138099v1 (2020). [Google Scholar]
 C. Silva, H. Maurer and D. Torres, Optimal control of a tuberculosis model with state and control delays. Math. Biosci. Eng. 14 (2017) 321–337. [CrossRef] [MathSciNet] [PubMed] [Google Scholar]
All Figures
Figure 1 Test 1: p_{M} = 1, p_{u} = 1, p_{D} = 1. The optimal control is 0. Without confinement, the death toll is greater than 10% and the peak value is very large, as expected. 

In the text 
Figure 2 Test 2: p_{M} = 1, p_{u} = 0.000001, p_{D} = 0. In this case, the main objective is to reduce the peak value. We can observe the delay, between the confinement decision and the response in the hospitalized population. Indeed, when the confinement stops, the size of hospitalized population remains constant for several days. The death toll is greater than 10%. 

In the text 
Figure 3 Test 3: p_{M} = 1e − 5, p_{u} = 0, p_{D} = 1. The main objective is to reduce the death toll, without considering the confinement cost. Therefore, the optimal policy assumes themaximum confinement value until the last days. Since the saturation of hospital capacity impacts in the death rate, as a consequence the peak value is reduced. 

In the text 
Figure 4 Test 4: p_{M} = 1, p_{u} = 0.0005, p_{D} = 1. In this case,we want to reduce the peak value and the death toll, taking into account the economic impact. 

In the text 
Figure 5 Test 5: p_{M} = 1, p_{u} = 0.000001, p_{D} = 0. Since hospitalizations rates ${\overline{\nu}}_{1}$ and ${\overline{\nu}}_{2}$ are assumed equals, to reduce the amount of hospitalizations that occur at the same time, it makes sense that both populations are confined for almost the same period. In this case, once the peak is reached, the confinement stops. We still observe a large number of deaths, greater than 10%. 

In the text 
Figure 6 Test 6: p_{M} = 1, p_{u} = 0.0005, p_{D} = 1. As in Test 4, the aim is to balance between reducing peak value, death toll and confinement cost. We observe an important difference between optimal confinement for strong and weak populations. Of course, this occurs due to the difference between the mortality rates for both populations. 

In the text 
Figure 7 Test 7: p_{M} = 1, p_{u} = 0.0005, p_{D} = 1, with the confinement lengths M_{1} = 25 and M_{2} = 45. The objective of this test is to show how the optimal policy changes when we add constraints in the confinement lengths. Taking the delay into account, we see that the end of confinement of the strong population corresponds to the end of the hospitalization peak. 

In the text 
Figure 8 Quasiperiodic case: p_{M} = 12, p_{u} = 0.005, p_{D} = 1, C = 0.01. 

In the text 
Figure 9 Longterm behavior: p_{M} = 12, p_{u} = 0.005, p_{D} = 1, C = 0.01, T = 260. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.