OPTIMAL CONTROL TECHNIQUES BASED ON INFECTION AGE FOR THE STUDY OF THE COVID-19 EPIDEMIC∗

We propose a model for the COVID-19 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. Received May 18, 2020. Accepted September 3, 2020.


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. [8][9][10]. 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 u(t, a) = 0, total confinement if u(t, a) = 1). We consider the following dynamics (t, a) = −δ(a)(1 − u(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)(1 − u(t, a))Z(t)y(t, a) Here, δ(a) ≥ 0 is the "transmission coefficient", Z(t) is the "global transmission factor", modelled as 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 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), (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ν(a),η(a) andγ(a) resp. Then, the coefficients take the form: Consequently, h(t, a, b) = 0 whenever b ∈ [0, n 0 ] and So, in this case, the model has four parameters for each value of a: δ(a),ν(a),η(a) andγ(a).

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: The total death toll is the weighted confinement cost, with weights c e (a) ≥ 0. The cost to be minimized is a combination of several costs. Its expression is 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 Indeed, if p M is nonzero, the optimal 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:

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 ∈ N. Let n 0 , n b , with n 0 < n b be the duration of the incubation phase and of the infection, respectively.  The global transmission factor and the size of the hospitalized population are respectively Also E k ∈ [0, 1] is a discrete counterpart to the hospital saturation estimate E(t), given by Willing to keep conservation laws, we may write the dynamics as follows: for all a = 1, . . . , n a , (3.4) The death toll per age and day, and per day, are resp.
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 k a remains nonnegative iff 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 k a 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 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 M a ∈ [0, 1], M ≥ 0, and M a ≥ 0: (3.7) The expression of the cost is with c u weighted confinement cost, D T total death toll: 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.

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 for 1 ≤ j ≤ n 0 − 1, e a,j = 1 for n 0 ≤ j ≤ n b , and ν a,j = 0 for 1 ≤ j ≤ n 0 − 1. (4.1) The above assumptions together with the initial condition h k a,1 = 0 implies that 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ν a ,η a ,γ a . We assume that we know for each age the total proportionν of transfers from z to h, that satisfieŝ so that we can estimateν Similarly, we may assume that we know for each age the proportionη a of hospitalized people that will die, when the hospital capacity is not saturated. Since h k a,j = 0, for j = 1, . . . , n 0 , we have by similar arguments As expected,ν andη take values in [0, 1]. When the total hospital capacity is exceeded we assume that it can be measured so thatγ a itself can be evaluated.

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 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 We can express c a,j as function of c a,1 : So the total and infectious population of z for each age are resp.
In the absence of confinement at initial time, estimating c a,1 with the previous relation, we deduce the value of δ from: As noted in Remark 3.1, in order that y k remains nonnegative, we need that 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" λ,ν,η, y 0 , Z 0 .
For the estimate of the γ parameters, specific methods are needed, which are beyond the scope of this paper.

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: From these data we can estimate the values ofν andη that represent the proportion of infected population in needs of hospitalization and the proportion of deaths with respect to the hospitalizations, respectively.
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 consider that 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 In order to obtain the values ofη 1 andη 2 , we assume that ratio is the same, i.e.,η 1 /η 2 = 0.0248. Normalizing the population, the global death rateη is equal tô η =η 1 y 1 +η 2 y 2 = (0.0248η 2 )0.734 +η 2 0.266. Since we do not have detailed information about the age of hospitalized people, we assume: 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 ifZ 0 is the initial confirmed cases, then we can define the (normalized) initial confirmed cases by age, as where y 0 = 6 706 3703 is the total population of France.
Since we have taken the same values ofν 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: In addition, we assume:γ a,j =η a,j , ∀a = 1, 2, 1 ≤ j ≤ n b ; C = 0.005. (5.10)

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 k 1 = u k 2 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. 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 a month 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 and the 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 that in 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.

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 The optimal control is 0. Without confinement, the death toll is greater than 10% and the peak value is very large, as expected. 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 7: p M = 1, p u = 0.0005, p D = 1, Figure 7 In the following table, we summarized the results of the tests with differentiated confinement: 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%.
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, the confinement 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 The main objective is to reduce the death toll, without considering the confinement cost. Therefore, the optimal policy assumes the maximum 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 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.

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 the current data.

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 p M = 12, p u = 0.005, p D = 1 and assume that the hospital capacity is C = 0.01. See Figure 8. In this case, we want to reduce the peak value and the death toll, taking into account the economic impact.
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.

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 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 Figure 5. Test 5: p M = 1, p u = 0.000001, p D = 0. Since hospitalizations ratesν 1 andν 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%. 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.

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ν, which at the beginning was close to 70%. Given that there are apparently many cases that are not serious, 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.
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. 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.

Conclusion
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 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 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) = z(t, a, B) + w(t, a, B) + h(t, a, B).
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 : 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: 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.

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 containment is 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: 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 /∆ to be integer. Clearly the following system is a consistent approximation of (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.