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

© The authors. Published by EDP Sciences, 2020

Licence Creative CommonsThis 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 COVID-19 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 COVID-19 epidemic has already been modelized in e.g. [810]. These papers use variants and extensions of the SIR (Susceptible-Infected-Recovered) 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 S ˙ (t) = δ(1u(t))S(t)I(t) I ˙ (t) = δ(1u(t))S(t)I(t)(η+μ)I(t) R ˙ (t) = μI(t) \begin{equation*} \begin{array}{lll} \dot S(t) & = & - \delta (1-u(t)) S(t) I(t) \vspace{1mm} \\ \dot I(t) & = & \delta (1-u(t)) S(t) I(t) -(\eta+\mu)I(t) \vspace{1mm} \\ \dot R(t) & =& \mu I(t) \end{array}\end{equation*}(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) ≤ uM, with uM ∈ (0, 1). Of course it is not easy to specify a confinement cost, but we might take c(u):= p u 0 T φ (u(t))dt \begin{equation*} c(u) := p_u \int_0^T \varphi (u(t)) \textrm{d} t \end{equation*}(1.2)

with pu ≥ 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 COVID-19 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 age-structured 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 COVID-19, 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 COVID-19 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/covid-19/covid-19.html or https://www.fceia.unr.edu.ar/~justina/covid-19/.

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 non-immunized (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 { y ˙ (t,a)=δ(a)(1u(t,a))Z(t)y(t,a) y(0,a)= y 0 (a), ( z t + z b )(t,a,b)=ν(a,b)z(t,a,b) z(0,a,b)= z 0 (a,b)>0,z(t,a,0)=δ(a)(1u(t,a))Z(t)y(t,a) ( h t + h b )(t,a,b)=ν(a,b)z(t,a,b)( η(a,b)+γ(a,b)E(t) )h(t,a,b) h(0,a,b)= h 0 (a,b)0,h(t,a,0)=0 y ¯ . (t,a)=z(t,a,B)+h(t,a,B) y ¯ (0,a)=0. \begin{equation*}\left\{\begin{array}{l} \dot{y}(t,a)=-\delta(a)(1-u(t,a))Z(t)y(t,a)\\[6pt] y(0,a)=y_0(a),\\[6pt] (z_t+z_b)(t,a,b)=-\nu(a,b)z(t,a,b)\\[6pt] z(0,a,b)=z_0(a,b)>0,\;\;\;\;z(t,a,0)=\delta(a)(1-u(t,a))Z(t)y(t,a)\\[6pt] (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)\\[6pt] h(0,a,b)=h_0(a,b)\geq 0, \;\;\;\;h(t,a,0)=0\\[6pt] \dot{{\bar{y}}}(t,a)=z(t,a,B)+h(t,a,B)\\[6pt] {\bar{y}}(0,a)=0.\\[6pt] \end{array} \right. \end{equation*}(2.1)

Here, δ(a) ≥ 0 is the “transmission coefficient”, Z(t) is the “global transmission factor", modelled as Z(t)= 0 A 0 B e (a,b)z(t,a,b)dadb, \begin{equation*} Z(t)=\int_0^A\int_0^B e(a,b) z(t,a,b)\textrm{d} a\textrm{d} b, \end{equation*}(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):= (H(t)C) + H(t)+C ;H(t):= 0 A 0 B h (t,a,b)dadb, \begin{equation*} E(t) := \frac{(H(t)-C)_+}{H(t)+C}; \quad H(t) := \int_0^A\int_0^B h(t,a,b)\textrm{d} a\textrm{d} b, \end{equation*}(2.3)

where C > 0 is the nominal capacity. We assume that only hospitalized patients die, with death rate d(t,a,b):=η(a,b)+γ(a,b)E(t), \begin{equation*} d(t,a,b) := \eta(a,b) + \gamma(a,b) E(t), \end{equation*}(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, n0) during which transmission and hospitalization do not occur. Additionally, in the infectious phase for which b ∈ [n0, B], we assume that e(a, b) = 1 and that ν(a, b), η(a, b) and γ(a, b) have constant values denoted by ν ¯ (a) ${\bar{\nu}}(a)$, η ¯ (a) ${\bar{\eta}}(a)$ and γ ¯ (a) ${\bar{\gamma}}(a)$ resp. Then, the coefficients take the form: e(a,b)= 1 [ n 0 ,B] (b),ν(a,b)= 1 [ n 0 ,B] (b) ν ¯ (a),η(a,b)= 1 [ n 0 ,B] (b) η ¯ (a),andγ(a,b)= 1 [ n 0 ,B] (b) γ ¯ (a). \begin{equation*} e(a,b)= \mathbb{1}_{[n_0,B]}(b), \;\;\; \nu(a,b)=\mathbb{1}_{[n_0,B]}(b){\bar{\nu}}(a), \;\;\; \eta(a,b)=\mathbb{1}_{[n_0,B]}(b){\bar{\eta}}(a),\;\;\text{and}\;\; \gamma(a,b)=\mathbb{1}_{[n_0,B]}(b){\bar{\gamma}}(a). \end{equation*}(2.5)

Consequently, h(t, a, b) = 0 whenever b ∈ [0, n0] and Z(t)= 0 A n 0 B z (t,a,b)dbda. \begin{equation*} Z(t)=\int_0^A\int_{n_0}^B z(t,a,b)\textrm{d} b\textrm{d} a. \end{equation*}(2.6)

So, in this case, the model has four parameters for each value of a: δ(a), ν ¯ (a) ${\bar{\nu}}(a)$, η ¯ (a) ${\bar{\eta}}(a)$ and γ ¯ (a) ${\bar{\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: D(a):= y(0,a)+ y ¯ (0,a)+ 0 B z (0,a,b)+h(0,a,b)db ( y(T,a)+ y ¯ (T,a)+ 0 B ( z(T,a,b)+h(T,a,b) ) db ). \begin{equation*} \begin{array}{ll} D(a) :=& y(0,a)+{\bar{y}}(0,a) +\int_0^Bz(0,a,b)+h(0,a,b)\textrm{d} b\\[6pt] \;&- \left( y(T,a)+{\bar{y}}(T,a)+\int_0^B\left(z(T,a,b)+h(T,a,b)\right)\textrm{d} b\right). \end{array} \end{equation*}(2.7)

The total death toll is D T := 0 A D (a)da. $D_T := \int_0^A D(a) \textrm{d} a.$ Denote by c(u):= 0 T 0 A c e (a)u(t,a)dadt \begin{equation*} c(u) := \int_0^T\int_0^A c_e(a) u(t,a)\textrm{d} a \textrm{d} t \end{equation*}(2.8)

the weighted confinement cost, with weights ce(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 , \begin{equation*} J(M,u,D_T) := p_M M + p_u c(u) + p_D D_T, \end{equation*}(2.9)

where the penalty coefficients pM, pu and pD are nonnegative. Here M is an additional variable, called the peak value, subject to H(t)M,t[0,T]. \begin{equation*} H(t)\leq M, \quad\forall t\in[0,T]. \end{equation*}(2.10)

Indeed, if pM 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: 0 T u (t,a)dtM(a),fora.a.$a(0,A)$. \begin{equation*}\int_0^{T} u(t,a) \textrm{d} t \leq M(a), \quad \text{for a.a. $a\in (0,A)$.} \end{equation*}(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, …, na}. As both time unit, time step and discretization of the infection age we choose one day.

The horizon of the problem will be N $N\in\mathbb{N}$. Let n0, nb, with n0 < nb be the duration of the incubation phase and of the infection, respectively.

The time index is k $k\in{\mathbb{N}}$, starting from value 0. The infection age is j ∈{1, …, nb}. The initial conditions are y a 0 = y 0 (a) a=1,, n a , z a,j 0 = z 0 (a,j) a=1,, n a ,j=1,, n b , h a,j 0 = h 0 (a,j) a=1,, n a ,j=1,, n b . \begin{equation*} \begin{array}{rcll} y^0_{a}&=&y_0(a) & a=1,\dots, n_a,\\[6pt] z^0_{a,{j}}&=&z_0(a,{j})& a=1,\dots,n_a, \;\; {j} =1,\dots,n_b,\\[6pt] h^0_{a,{j}}&=&h_0(a,{j}) & a=1,\dots, n_a, \;\; {j} =1,\dots,n_b.\\[6pt] \end{array} \end{equation*}(3.1)

The global transmission factor and the size of the hospitalized population are respectively Z k = a=1 n a j=1 n b e a,j z a,j k , H k = a=1 n a j=1 n b h a,j k . \begin{equation*} Z^k= \sum_{a=1}^{n_a}\sum_{{j} =1}^{n_b}e_{a,{j}}z^k_{a,{j}},\quad H^k=\sum_{a=1}^{n_a}\sum_{{j} =1}^{n_b}h^k_{a,{j}}. \end{equation*}(3.2)

Also Ek ∈ [0, 1] is a discrete counterpart to the hospital saturation estimate E(t), given by E k := ( H k C) + H k +C . \begin{equation*} E^k := \frac{(H^k-C)_&#x002B;}{H^k&#x002B;C}. \end{equation*}(3.3)

Willing to keep conservation laws, we may write the dynamics as follows: for all a = 1, …, na, y a k+1 = ( 1 δ a (1 u a k ) Z k ) y a k z a,1 k+1 = δ a (1 u a k ) Z k y a k z a,j k+1 = (1 ν a,j1 ) z a,j1 k j=2,, n b h a,1 k+1 = 0 h a,j k+1 = ν a,j1 z a,j1 k +(1 η a,j1 γ a,j1 E k ) h a,j1 k j=2,, n b y ¯ a k+1 = y ¯ a k + z a, n b k + h a, n b k . \begin{equation*} \begin{array}{rcll} y^{k&#x002B;1}_{a}&=& \left(1 - \delta_a (1-u^k_a) Z^k\right)y^k_{a} & \;\\[6pt] z^{k&#x002B;1}_{a,1}&=& \delta_a (1-u^k_a) Z^k y^k_a&\;\\[6pt] z^{k&#x002B;1}_{a,{j}}&=&(1-\nu_{a,{j}-1})z^k_{a,{j}-1}& {j} =2,\dots, n_b\\[6pt] h^{k&#x002B;1}_{a,1}&=&0&\;\\[6pt] h^{k&#x002B;1}_{a,{j}}&=& \nu_{a,{j}-1}z^k_{a,{j}-1} &#x002B; (1-\eta_{a,{j}-1} - \gamma_{a,{j}-1} E^k)h^k_{a,{j}-1}& {j} =2, \dots, n_b\\[6pt] {\bar{y}}^{k&#x002B;1}_{a}&=&{\bar{y}}^k_{a}&#x002B;z^k_{a,n_b}&#x002B;h^k_{a,n_b}. & \; \end{array}\end{equation*}(3.4)

The death toll per age and day, and per day, are resp. D a k := j=1 n b 1 ( η a,j + γ a,j E k ) h a,j k ; D k := a=1 n a D a k . \begin{equation*} D^k_a :=\sum_{{j} =1}^{n_b-1}(\eta_{a,{j}} &#x002B; \gamma_{a,{j}} E^k)h^k_{a,{j}}; \quad D^k := \sum_{a=1}^{n_a} D^k_a. \end{equation*}(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 $y^k_{a}$ remains nonnegative iff

δ a (1 u a k ) Z k 1. \begin{equation*} \delta_a (1-u^k_a) Z^k\leq 1.\end{equation*}(3.6)

In the case of normalized, (i.e., summing to 1) initial population, we have that Zk ≤1. So, the previous nonnegativity condition holds whenever δa ≤ 1. In our numerical test, it may happen that δa ≥ 1, but y a k $y^k_{a}$ always remained positive.

On the other hand, in our tests ηa,j + γa,j ≤ 1. Since Ek ≤ 1 by definition, it follows that h a,j k $h^k_{a,{j}}$ 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 [0,1] $u^M_a \in [0,1]$, M ≥ 0, and Ma ≥ 0: 0 u a k u a M ,k=0,,N,a{1, n a } H k M,k=0,,N, k=0 N1 u a k M a ,a{1, n a }. \begin{equation*} \begin{array}{l} 0\leq u^k_a\leq u^M_a,\;\;\;k=0,\dots, N, \; a \in \{1,n_a\} \\[6pt] H^k\leq M, \quad k=0,\dots, N,\\[6pt] \sum_{k=0}^{N-1} u^k_a \leq M_a, \; a \in \{1,n_a\}. \end{array}\end{equation*}(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 , \begin{equation*} J_d(M,u,D_T) := p_M M&#x002B;p_u c_u &#x002B; p_D D_T, \end{equation*}(3.8)

with cu weighted confinement cost, DT total death toll: c u := k=0 N1 a=1 n a c e,a u a k ; D T = k=0 N1 D k , \begin{equation*} c_u := \sum_{k=0}^{N-1} \sum_{a=1}^{n_a}c_{e,a}u_a^k;\quad D_T = \sum_{k=0}^{N-1} D^k, \end{equation*}(3.9)

and pu, pM and pD are nonnegative penalty coefficients.

Due to the incubation phase, the value of control for the last n0 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 n0 − 1 of infection) where no transmission and hospitalizations occur. So, e a,j =0for1j n 0 1, e a,j =1for n 0 j n b ,and ν a,j =0for1j n 0 1. \begin{equation*} e_{a, {j}}= 0 \;\;\text{for}\;\;1\leq {j} \leq n_0-1, \;\;\; e_{a, {j}}=1 \;\;\text{for}\;\; n_0\leq {j}\leq n_b, \;\;\;\text{and} \;\; \nu_{a,{j}}=0 \;\;\text{for}\;\; 1\leq {j}\leq n_0-1.\end{equation*}(4.1)

The above assumptions together with the initial condition h a,1 k =0 $h^k_{a,1}=0$ implies that Z a k = a=1 n a j= n 0 n b z a,j k ,and h a,j k =0, whenever 1j n 0 . \begin{equation*} Z^k_a=\sum_{a=1}^{n_a}\sum_{{j}=n_0}^{n_b}z^k_{a,{j}}, \;\;\text{and}\;\;h^k_{a,{j}}=0,\; \text{ whenever }\; 1 \leq {j}\leq n_0. \end{equation*}(4.2)

During the infectious phase (j varying from n0 + 1 to nb) we may assume that the coefficients νa,j, for jn0 ηa,j and γa,j, for jn0 + 1, do not depend on j and therefore are equal to some positive constants ν ¯ a ${\bar{\nu}}_a$, η ¯ a ${\bar{\eta}}_a$, γ ¯ a ${\bar{\gamma}}_a$.

We assume that we know for each age the total proportion ν ^ ${\hat{\nu}}$ of transfers from z to h, that satisfies ν ^ a = ν ¯ a + ν ¯ a (1 ν ¯ a )++ ν ¯ a (1 ν ¯ a ) n b n 0 1 = ν ¯ a j=0 n b n 0 1 (1 ν ¯ a ) j =1 (1 ν ¯ a ) n b n 0 , \begin{equation*} {\hat{\nu}}_a = {\bar{\nu}}_a &#x002B; {\bar{\nu}}_a (1-{\bar{\nu}}_a) &#x002B; \cdots &#x002B; {\bar{\nu}}_a (1-{\bar{\nu}}_a)^{n_b-n_0-1} = {\bar{\nu}}_a \sum_{{j} =0}^{n_b-n_0-1} (1-{\bar{\nu}}_a)^{j} = 1 - (1-{\bar{\nu}}_a)^{n_b-n_0}, \end{equation*}(4.3)

so that we can estimate ν ¯ a =1 (1 ν ^ a ) 1/( n b n 0 ) . \begin{equation*}{\bar{\nu}}_a = 1 - (1-{\hat{\nu}}_a)^{1/(n_b-n_0)}. \end{equation*}(4.4)

Similarly, we may assume that we know for each age the proportion η ^ a ${\hat{\eta}}_a$ of hospitalized people that will die, when the hospital capacity is not saturated. Since h a,j k =0 $h^k_{a,{j}}=0$, for j = 1, …, n0, we have by similar arguments η ¯ a =1 (1 η ^ a ) 1/( n b n 0 1) . \begin{equation*}{\bar{\eta}}_a = 1 - (1-{\hat{\eta}}_a)^{1/(n_b-n_0-1)}. \end{equation*}(4.5)

As expected, ν ¯ ${\bar{\nu}}$ and η ¯ ${\bar{\eta}}$ take values in [0, 1].

When the total hospital capacity is exceeded we assume that it can be measured so that γ ¯ a ${\bar{\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 = g, so that λ=(logg)/p. \begin{equation*} \lambda = (\log g)/p.\end{equation*}(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 e λ(kj) . \begin{equation*}z^k_{a,{j}} = c_{a,{j}} \textrm{e}^{\lambda (k-{j})}. \end{equation*}(4.7)

We can express ca,j as function of ca,1: c a,j ={ c a,1 if j n 0 c a,1 (1 ν ¯ a ) j n 0 otherwise. \begin{equation*} c_{a,{j}} = \left\{ \begin{array}{lll} c_{a,1} & \text{if } {j}\leq n_0 \vspace{1mm} \\ c_{a,1}(1-{\bar{\nu}}_a)^{{j}-n_0} & \text{otherwise.} \end{array}\right.\end{equation*}(4.8)

So the total and infectious population of z for each age are resp. Z ¯ a k := j=1 n b z a,j k = c a,1 χ ¯ a e λk ; χ ¯ a :=( j=1 n 0 1 e λj + j= n 0 n b (1 ν ¯ a ) b n 0 e λj ) \begin{equation*}{\bar{Z}}^k_a := \sum_{{j} =1}^{n_b} z^k_{a,{j}} = c_{a,1} {\bar{\chi}}_a \textrm{e}^{\lambda k}; \quad {\bar{\chi}}_a :=\left(\sum_{{j} =1}^{n_0-1} \textrm{e}^{-\lambda {j} } &#x002B; \sum_{{j} =n_0}^{n_b} (1-{\bar{\nu}}_a)^{b-n_0} \textrm{e}^{-\lambda {j} } \right) \end{equation*}(4.9) Z a k := j= n 0 n b z a,j k = c a,1 χ a e λk ; χ a := j= n 0 n b (1 ν ¯ a ) j n 0 e λj . \begin{equation*}Z^k_a := \sum_ {{j} =n_0}^{n_b} z^k_{a,{j}} = c_{a,1} \chi_a \textrm{e}^{\lambda k} ; \quad \chi_a := \sum_{{j} =n_0}^{n_b} (1-{\bar{\nu}}_a)^{{j}-n_0} \textrm{e}^{-\lambda {j} }. \end{equation*}(4.10)

In the absence of confinement at initial time, estimating ca,1 with the previous relation, we deduce the value of δ from: c a,1 = z a,1 1 = δ a Z 0 y a 0 . \begin{equation*}c_{a,1} = z^1_{a,1} = \delta_a Z^0 y^0_a. \end{equation*}(4.11)

As noted in Remark 3.1, in order that yk remains nonnegative, we need that δ a (1 u a k ) Z k 1,forallk=0,,N. \begin{equation*} \delta_a (1 - u^k_a) Z^k \leq 1, \quad \text{for all $k=0,\ldots,N$}.\end{equation*}(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 h0 = 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" λ, ν ^ ${\hat{\nu}}$, η ^ ${\hat{\eta}}$, y0, Z0.

Then we can deduce ν ¯ ${\bar{\nu}}$ and η ¯ ${\bar{\eta}}$ from (4.4) and (4.5), respectively. The value of χ and χ ¯ ${\bar{\chi}}$ follow. We can deduce ca,1 from (4.10) and the values of z0 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 ν ^ ${\hat{\nu}}$ and η ^ ${\hat{\eta}}$ that represent the proportion of infected population in needs of hospitalization and the proportion of deaths with respect to the hospitalizations, respectively. ν ^ =0.726,and η ^ =0.165. \begin{equation*} {\hat{\nu}} = 0.726, \quad\text{and}\;\;\; {\hat{\eta}}= 0.165.\end{equation*}(5.1)

Furthermore, from the above-mentioned 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 λ[0.13,0.15]. \begin{equation*} \lambda \in [0.13,0.15]. \end{equation*}(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 y1 can be considered as the population up to 59 years and the weak population y2 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 η 1 r = 416 49213491 =8.45e6,and η 2 r = 6021 17850240 =3.4e4, \begin{equation*} \eta^r_1= \frac{416}{49\,213\,491}=8.45 e-6, \;\;\;\text{and}\;\;\;\eta^r_2=\frac{6021}{17\,850\,240}= 3.4\textrm{e}-4, \end{equation*}(5.3)

which implies that the ratio of this rate is η 1 r η 2 r =0.0248. \begin{equation*} \frac{\eta^r_1}{\eta^r_2} = 0.0248. \end{equation*}(5.4)

In order to obtain the values of η ^ 1 ${\hat{\eta}}_1$ and η ^ 2 ${\hat{\eta}}_2$, we assume that ratio is the same, i.e., η ^ 1 / η ^ 2 =0.0248 ${\hat{\eta}}_1/{\hat{\eta}}_2 = 0.0248 $. Normalizing the population, the global death rate η ^ ${\hat{\eta}}$ is equal to η ^ = η ^ 1 y 1 + η ^ 2 y 2 =(0.0248 η ^ 2 )0.734+ η ^ 2 0.266. \begin{equation*} {\hat{\eta}} = {\hat{\eta}}_1 y_1 &#x002B; {\hat{\eta}}_2 y_2 = (0.0248 {\hat{\eta}}_2) 0.734 &#x002B;{\hat{\eta}}_2 0.266. \end{equation*}(5.5)

By (5.1) we obtain, η ^ 1 =0.014, η ^ 2 =0.58. \begin{equation*} {\hat{\eta}}_1= 0.014, \;\;\;\; {\hat{\eta}}_2= 0.58. \end{equation*}(5.6)

Since we do not have detailed information about the age of hospitalized people, we assume: ν ^ = ν ^ 2 = ν ^ 1 . \begin{equation*} {\hat{\nu}}={\hat{\nu}}_2={\hat{\nu}}_1. \end{equation*}(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 (Zk), we consider that if Z ¯ 0 ${\bar{Z}}^0$ is the initial confirmed cases, then we can define the (normalized) initial confirmed cases by age, as Z ¯ 1 0 = Z ¯ 0 0.734 y 0 , Z ¯ 2 0 = Z ¯ 0 0.266 y 0 , \begin{equation*} {\bar{Z}}^0_1= {\bar{Z}}^0 \frac{0.734}{y^0}, \;\;\;\;\;\; {\bar{Z}}^0_2= {\bar{Z}}^0 \frac{0.266}{y^0}, \end{equation*}(5.8)

where y0 = 6 706 3703 is the total population of France.

Since we have taken the same values of ν ¯ ${\bar{\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 $y^0_1$ Initial strong population 0.734

y 2 0 $y^0_2$ Initial weak population 0.266

Z ¯ 1 0 ${\bar{Z}}^0_1$ Initial infected strong population 7.26e − 5

Z ¯ 2 0 ${\bar{Z}}^0_2$ Initial infected weak population 2.63e − 5

λ Growth coefficient given by (4.6) 0.13

ν ^ 1 ${\hat{\nu}}_1$ Total proportion of transfers from z1 to h1 0.726

ν ^ 2 ${\hat{\nu}}_2$ Total proportion of transfers from z2 to h2 0.726

η ^ 1 ${\hat{\eta}}_1$ Strong minimal death rate 0.014

η ^ 2 ${\hat{\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: δ 1 =1.656, ν ¯ 1 =0.149412, η ¯ 1 =0.002012 δ 2 =1.656, ν ¯ 2 =0.149412, η ¯ 2 =0.116557 \begin{equation*} \begin{array}{lll} \delta_1=1.656, &\; {\bar{\nu}}_1=0.149412, &\; {\bar{\eta}}_1=0.002012\\[6pt] \delta_2=1.656, &\; {\bar{\nu}}_2=0.149412, &\; {\bar{\eta}}_2= 0.116557 \end{array} \end{equation*}(5.9)

In addition, we assume: γ ¯ a,j = η ¯ a,j ,a=1,2,1j n b ;C=0.005. \begin{equation*} {\bar{\gamma}}_{a,{j}}={\bar{\eta}}_{a,{j}},\;\;\;\;\forall a=1,2, \; 1\leq {j}\leq n_b; \;\; C=0.005. \end{equation*}(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 $u^k_1=u^k_2$ for all k, that the upper bound for the control is uM = 0.75, and that the confinement cost, given by (3.9), satisfies ce,1 + ce,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: pM = 1, pu = 1, pD = 1, Figure 1.

Test 2: pM = 1, pu = 0.000001, pD = 0, Figure 2

Test 3: pM = 1e − 5, pu = 0, pD = 1, Figure 3

Test 4: pM = 1, pu = 0.0005, pD = 1, Figure 4

thumbnail Figure 1

Test 1: pM = 1, pu = 1, pD = 1. The optimal control is 0. Without confinement, the death toll is greater than 10% and the peak value is very large, as expected.

thumbnail Figure 2

Test 2: pM = 1, pu = 0.000001, pD = 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%.

thumbnail Figure 3

Test 3: pM = 1e − 5, pu = 0, pD = 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.

thumbnail Figure 4

Test 4: pM = 1, pu = 0.0005, pD = 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):= k=0 N u k $I(u):=\sum_{k=0}^{N}u^k$ and the peak value M.

Test D1,T D2,T DT 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 ce,1 = 0.734 and ce,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: pM = 1, pu = 0.000001, pD = 0, Figure 5

thumbnail Figure 5

Test 5: pM = 1, pu = 0.000001, pD = 0. Since hospitalizations rates ν ¯ 1 ${\bar{\nu}}_1$ and ν ¯ 2 ${\bar{\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: pM = 1, pu = 0.0005, pD = 1, Figure 6

thumbnail Figure 6

Test 6: pM = 1, pu = 0.0005, pD = 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 M1 = 25 and M2 = 45.

Test 7: pM = 1, pu = 0.0005, pD = 1, Figure 7

thumbnail Figure 7

Test 7: pM = 1, pu = 0.0005, pD = 1, with the confinement lengths M1 = 25 and M2 = 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 D1,T D2,T DT I(u1) I(u2) 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 quasi-periodic 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 pM = 12, pu = 0.005, pD = 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 quasi-periodic 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.

thumbnail Figure 8

Quasi-periodic case: pM = 12, pu = 0.005, pD = 1, C = 0.01.

5.4.2 Long-term 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 pM = 12, pu = 0.005, pD = 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.

thumbnail Figure 9

Long-term behavior: pM = 12, pu = 0.005, pD = 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 ν ¯ ${\bar{\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: y ˙ (t,a) = δ(a)(1u(t,a))Z(t)y(t,a)v(t,a)y(t,a), y ¯ ˙ (t,a) = z(t,a,B)+h(t,a,B)+v(t,a)y(t,a). \begin{equation*} \begin{array}{lll} \dot{y}(t,a) &=& -\delta(a)(1-u(t,a))Z(t)y(t,a)-v(t,a)y(t,a), \vspace{1mm}\\ \dot{\bar{y}}(t,a) &=& z(t,a,B)&#x002B;h(t,a,B) &#x002B; v(t,a) y(t,a). \end{array} \end{equation*}(6.1)

It is natural to think that v depends on the age of the population, and the limitation of MV vaccines per day would read as 0 A v (t,a)y(t,a)da M V ,for a.a.t(0,T). \begin{equation*} \int_0^A v(t,a) y(t,a) \textrm{d} a \leq M_V, \quad \text{for a.a. $t\in (0,T)$.} \end{equation*}(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 non-confirmed (among then some asymptomatic cases). Keeping the variable z for the non-confirmed 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: { y ˙ (t,a)=δ(a)(1u(t,a))Z(t)y(t,a) ( z t + z b )(t,a,b)=μ(t,a,b)z(t,a,b)ν(a,b)z(t,a,b) ( w t + w b )(t,a,b)=μ(t,a,b)z(t,a,b)ν(a,b)w(t,a,b) ( h t + h b )(t,a,b)=ν(a,b)(z(t,a,b)+w(t,a,b))( η(a,b)+γ(a,b)E(t) )h(t,a,b) y ¯ . (t,a)=z(t,a,B)+w(t,a,B)+h(t,a,B). \begin{equation*} \left\{\begin{array}{l} \dot{y}(t,a)=-\delta(a)(1-u(t,a))Z(t)y(t,a)\\[6pt] (z_t&#x002B;z_b)(t,a,b)=-\mu(t,a,b)z(t,a,b)-\nu(a,b)z(t,a,b)\\[6pt] (w_t&#x002B;w_b)(t,a,b)=\mu(t,a,b)z(t,a,b)-\nu(a,b)w(t,a,b)\\[6pt] (h_t&#x002B;h_b)(t,a,b)=\nu(a,b)(z(t,a,b)&#x002B;w(t,a,b))-\left(\eta(a,b) &#x002B; \gamma(a,b) E(t)\right) h(t,a,b)\\[6pt] \dot{{\bar{y}}}(t,a)=z(t,a,B)&#x002B;w(t,a,B)&#x002B;h(t,a,B). \end{array} \right. \end{equation*}(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 c1 and c2: 0 T 0 A 0 B μ (t,a,b)( c 1 z(t,a,b)+ c 2 y(t,a,b))dbdadt. \begin{equation*} \int_0^T \int_0^A \int_0^B \mu(t,a,b) ( c_1 z(t,a,b)&#x002B; c_2 y(t,a,b)) \textrm{d} b \textrm{d} a \textrm{d} t. \end{equation*}(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: { h(t,a,0)= 0 b ν (a,b)(z(t,a,b)+w(t,a,b))db ( h t + h c )(t,a,c)=( η(a,c)+γ(a,c)E(t) )h(t,a,c),c[0,C]. \begin{equation*} \left\{ \begin{array}{l} h(t,a,0)=\int_0^b \nu(a,b)(z(t,a,b)&#x002B;w(t,a,b)) \textrm{d} b\\[6pt] (h_t&#x002B;h_c)(t,a,c)=-\left(\eta(a,c) &#x002B; \gamma(a,c) E(t)\right) h(t,a,c), \;\;\; c\in [0,C]. \end{array} \right. \end{equation*}(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 age-structure. 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 non-severe 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 open-source 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: y a k+1 y a k = δ a (1 u a k ) Z k y a k z a,1 k+1 z a,1 k = z a,1 k + δ a (1 u a k ) Z k y a k z a,j k+1 z a,j k = z a,j k +(1 ν a,j1 ) z a,j1 k h a,j k+1 h a,j k = h a,j k + ν a,j1 z a,j1 k +(1 η a,j1 γ a,j1 E k ) h a,j1 k y ¯ a k+1 y ¯ a k = z a, n b k + h a, n b k . \begin{equation*} \begin{array}{rcl} y^{k&#x002B;1}_{a} - y^k_{a}&=& - \delta_a (1-u^k_a) Z^ky^k_{a} \\[6pt] z^{k&#x002B;1}_{a,1} - z^k_{a,1}&=& -z^k_{a,1} &#x002B; \delta_a (1-u^k_a) Z^k y^k_a\\[6pt] z^{k&#x002B;1}_{a,{j}} - z^k_{a,{j}}&=& -z^k_{a,{j}}&#x002B;(1-\nu_{a,{j}-1})z^k_{a,{j}-1}\\[6pt] h^{k&#x002B;1}_{a,{j}} - h^{k}_{a,{j}}&=& -h^{k}_{a,{j}} &#x002B; \nu_{a,{j}-1}z^k_{a,{j}-1} &#x002B; (1-\eta_{a,{j}-1} - \gamma_{a,{j}-1} E^k) h^k_{a,{j}-1}\\[6pt] {\bar{y}}^{k&#x002B;1}_{a} - {\bar{y}}^k_{a} &=& z^k_{a,n_b}&#x002B;h^k_{a,n_b}. \end{array}\end{equation*}(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 b := n b /Δ $n&#x0027;_b :=n_b/\UpDelta$ to be integer. Clearly the following system is a consistent approximation of (2.1): y a k+1 y a k Δ = δ a (1 u a k ) Z k y a k z a,1 k+1 z a,1 k Δ = z a,1 k Δ + δ a (1 u a k ) Z k y a k z a,j+1 k+1 z a,j+1 k Δ = z a,j k z a,j+1 k Δ ν a,j1 z a,j k h a,j+1 k+1 h a,j+1 k Δ = h a,j k h a,j+1 k Δ + ν a,j z a,j k ( η a,j + γ a,j E k ) h a,j k y ¯ a k+1 y ¯ a k Δ = z a, n b k + h a, n b k . \begin{equation*} \begin{array}{rcl} \displaystyle \frac{y^{k&#x002B;1}_{a} - y^k_{a}}{\UpDelta}&=& - \delta_a (1-u^k_a) Z^ky^k_{a} \\[6pt] \displaystyle \frac{z^{k&#x002B;1}_{a,1} - z^k_{a,1}}{\UpDelta}&=& \displaystyle \frac{-z^k_{a,1}}{\UpDelta} &#x002B; \delta_a (1-u^k_a) Z^k y^k_a\\[6pt] \displaystyle \frac{z^{k&#x002B;1}_{a,j&#x002B;1} - z^k_{a, j&#x002B;1}}{\UpDelta}&=& \displaystyle \frac{z^k_{a,{j}}-z^k_{a,j&#x002B;1}}{\UpDelta}-\nu_{a,j-1}z^k_{a,{j}}\\[6pt] \displaystyle \frac{h^{k&#x002B;1}_{a,j&#x002B;1} - h^{k}_{a,j&#x002B;1}}{\UpDelta}&=& \displaystyle \frac{h^k_{a,{j}}-h^{k}_{a,j&#x002B;1}}{\UpDelta} &#x002B; \nu_{a,j}z^k_{a,j} - (\eta_{a,j} &#x002B; \gamma_{a,j} E^k) h^k_{a,j}\\[6pt] \displaystyle \frac{{\bar{y}}^{k&#x002B;1}_{a} - {\bar{y}}^k_{a}}{\UpDelta} &=& z^k_{a,n&#x0027;_b}&#x002B;h^k_{a,n&#x0027;_b}. \end{array} \end{equation*}(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

  1. B. Ainseba, S. Aniţa and M. Langlais, Optimal control for a nonlinear age-structured population dynamics model. Electron. J. Differ. Equ. 9 (2003) 28. [Google Scholar]
  2. 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]
  3. J.F. Bonnans and J. Gianatti, Optimal control of state constrained age-structured problems. SIAM J. Control Optim. 58 (2020) 2206–2235. [Google Scholar]
  4. M. Brokate, Pontryagin’s principle for control problems in age-dependent population dynamics. J. Math. Biol. 23 (1985) 75–101. [CrossRef] [MathSciNet] [PubMed] [Google Scholar]
  5. R. Djidjou-Demasse, Y. Michalakis, M. Choisy, M.T. Sofonea and S. Alizon, Optimal COVID-19 epidemic control until vaccine deployment. Preprint MedRxiv: 20049189v2 (2020). [Google Scholar]
  6. R. Elie, E. Hubert and G. Turinici, Contact rate epidemic control of covid-19: an equilibrium view. MMNP 15 (2020) 35. [EDP Sciences] [Google Scholar]
  7. 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]
  8. 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 (covid-19) outbreak in Wuhan, China with individual reaction and governmental action. Int. J. Infect. Dis. 93 (2020) 211–216. [CrossRef] [PubMed] [Google Scholar]
  9. Z. Liu, P. Magal, O. Seydi and G. Webb, Predicting the cumulative number of cases for the COVID-19 epidemic in China from early data. Preprint MedRxiv: 20034314v1 (2020). [Google Scholar]
  10. Z. Liu, P. Magal, O. Seydi and G. Webb, Understanding unreported cases in the COVID-19 epidemic outbreak in Wuhan, China, and the importance of major public health interventions. Biology 9 (2020) 50. [Google Scholar]
  11. H. Maurer and M. do Rosario de Pinho, Optimal control of epidemiological SEIR models with L1-objectives and control-state constraints. Pac. J. Optim. 12 (2016) 415–436. [Google Scholar]
  12. Q. Richard, S. Alizon, M. Choisy, M.T. Sofonea and R. Djidjou-Demasse, Age-structured non-pharmaceutical interventions for optimal control of covid-19 epidemic. Preprint MedRxiv: 20138099v1 (2020). [Google Scholar]
  13. 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

thumbnail Figure 1

Test 1: pM = 1, pu = 1, pD = 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
thumbnail Figure 2

Test 2: pM = 1, pu = 0.000001, pD = 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
thumbnail Figure 3

Test 3: pM = 1e − 5, pu = 0, pD = 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
thumbnail Figure 4

Test 4: pM = 1, pu = 0.0005, pD = 1. In this case,we want to reduce the peak value and the death toll, taking into account the economic impact.

In the text
thumbnail Figure 5

Test 5: pM = 1, pu = 0.000001, pD = 0. Since hospitalizations rates ν ¯ 1 ${\bar{\nu}}_1$ and ν ¯ 2 ${\bar{\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
thumbnail Figure 6

Test 6: pM = 1, pu = 0.0005, pD = 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
thumbnail Figure 7

Test 7: pM = 1, pu = 0.0005, pD = 1, with the confinement lengths M1 = 25 and M2 = 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
thumbnail Figure 8

Quasi-periodic case: pM = 12, pu = 0.005, pD = 1, C = 0.01.

In the text
thumbnail Figure 9

Long-term behavior: pM = 12, pu = 0.005, pD = 1, C = 0.01, T = 260.

In the text

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.