ON APPLICATION OF OPTIMAL CONTROL TO SEIR NORMALIZED MODELS: PROS AND CONS

Abstract. In this work we normalize a SEIR model that incorporates exponential natural birth and death, as well as disease-caused death. We use optimal control to control by vaccination the spread of a generic infectious disease described by a normalized model with L1 cost. We discuss the pros and cons of SEIR normalized models when compared with classical models when optimal control with L1 costs are considered. Our discussion highlights the role of the cost. Additionally, we partially validate our numerical solutions for our optimal control problem with normalized models using the Maximum Principle.

1. Introduction.Since the publication of the seminal paper [?] mathematical compartmental models are widely used to describe infectious diseases dymanics in large populations (see, for example [?], [?], [?] and [?]).It is well accepted that once an infected individual comes into contact with an una ected population, the disease will spread by contact with the infectious individuals.Compartmental models divide the population into compartments characterizing the spread of the diseases and letters are used to denote the number of individuals in each compartment.Usually, the size of the population to be studied is N , S is the number of susceptible individuals and the number of infectious individuals is I.The letter R is also used to denote the number of those who recover from the disease or, as in [?], those who are removed from the disease by death or by recovery.The nature of the disease as well the reason why the models are studied may dictate the need for di erent compartments to be included.For example, SEIR models are used for diseases where infected individuals do not become immediately infectious; they are considered to be exposed to the disease and placed in compartment E, being moved to the I compartment only after some latent period of time.
The basic reproduction number, R 0 , measures the transmission potential of a disease (see, for example, [?]).It is defined as "the expected number of secondary method, we use known software packages, and we present numerical solutions that satisfy necessary optimality conditions with high accuracy.Notably, and di ering from [?], we work with a normalized SEIR model.
The normalized SEIR model di ers from the usual SEIR model since the variables are fractions of the whole population instead of the number of individuals in each compartment.The theoretical and numerical treatment involving the latter model is usually done as if the variables are continuous and not integers; treating such variables as integers would demand the use of integer programming what is known to be very heavy computationally.When we turn to normalized models the variables are, by nature, continuous.In the literature, normalized models are common when the total population is assumed to remain constant during the time frame under study.This is not our case; here we normalize a SEIR model that incorporates exponential natural birth and death, as well as disease-caused death (similarly to what is done in [?]).As far as optimal control is concerned, normalizing such model brings out some new issues related to the choice of costs and the introduction of non standard constraints, questions we discuss here when comparing optimal control for normalized and not normalized SEIR models.
Herein, we refer to the SEIR model, where the variables S, E, I and R denote the number of individuals in each compartment as the classical SEIR model as opposed to the normalized SEIR model.
We emphasize that we do not concentrate on any particular disease.Rather, our aim is to illustrate how previously proposed optimal control formulations can be handled by this new model, when di erent scenarios are considered.Taking into account that the set of parameters for the population in [?], based on [?], are not to be found in today's world, we use di erent population's parameters closed related to some European countries.
Like other models in epidemiology, SEIR models represent only a rough approximation of reality.However, they provide new insights into the spreading of diseases and, when optimal control is applied, new insight on di erent vaccination policies.
This paper is organized in the following way.In Section 2 we introduce an optimal control problem with L 1 cost involving the classical SEIR model and its "normalized" counterpart.A brief description of our numerical methods as well as the tables with the set of parameters used in our simulation are presented in Section 3. Numerical methods give a helping hand for the discussion on the pros and cons of optimal control problems with classical and normalized models presented in Section 4. This discussion highlights the need for a criterious choice of cost.In Sections 5 and 6, we solve our optimal control problem numerically partially validating our solution using the maximum principle and the computed multipliers.Conclusions follow in the last section of this paper.2. Optimal Control Problems for SEIR models with L 1 cost.The SEIR model is a compartmental model well accepted as modelling some infectious diseases.At each instant t, S(t), E(t), I(t), and R(t) denote the number of individuals in the susceptible, exposed, infectious and recovered compartments.The total population is N (t) = S(t) + E(t) + I(t) + R(t).We assume that our population has exponential natural birth b and death d and that that all newborns are susceptible to the diseases.The disease transmission is described by the parameters f , denoting the rate at which the exposed individuals become infectious, g, denoting the rate at which infectious individuals recover, a, denoting the death rate due to the disease and c, denoting the disease incidence coe cient of horizontal transmission (vertical transmission is not considered).The rate of transmission of the disease is c I(t)  N (t) .For simplicity the parameters a, b, c, d and f are assumed to be constants.For more information about such model we refer the reader to [?], [?], [?] and references within.
Optimal control techniques for SEIR models allow the study of di erent vaccines policies; di erent policies are confronted in [?] and [?] where the minimizing cost is L 2 , and in [?], covering the case of L 1 cost.As in the aforementioned works we assume that only susceptible people is vaccinated (implying that it is possible to distinguish between exposed and susceptible individuals) and that the vaccine is e ective so that all vaccinated susceptible individuals become immune.Let u(t) represent the percentage of the susceptible individuals being vaccinated per unit of time.Taking all the above considerations into account we are led to the following dynamical system (see [?]): where S 0 , E 0 , I 0 and N 0 are nonnegative initial conditions.For some ū Ø 0, we impose the following control constraint to u: where T represents the period of time under consideration.The recovered population is related to the total population by N (t) = S(t) + E(t) + I(t) + R(t).Hence, Here, the aim of applying optimal control to SEIR models is to control the spreading of the disease with some minimum financial cost.The cost should then be a weighted sum of the society financial costs of having, at each time, I(t) infected individuals and the cost of the vaccination e ort what is mathematical translated as (7).This leads to the L 1 cost as in [?]: where X = (S, E, I, N ) and A > 0, B > 0 are weight parameters related with financial costs.
Throughout this paper we refer to the optimal control problem of minimizing J 1 (X, u) over all X oe W 1,1 ([0, T ] : R 4 ) and all piecewise continuous functions u : [0, T ] ae R satisfying (1)-( 6) and ( 5) as the classical and basic optimal control problem: Next, we associate (P ) with a normalized optimal control problem.Normalizing models are obtained considering the percentage of the total population to be 1 at each instant t.Then, defining we have s(t) + e(t) + i(t) + r(t) = 1 for all t.
(9) Notice that s(t) is the percentage of the population in compartment S(t), e(t) is the percentage of the population in compartment E(t), i(t) is the percentage of the population in compartment I(t) and r(t) be the percentage of the population in compartment R(t).The normalized counterpart of ( 1)-( 6) is then: Remarkably, the dead rate parameters does not appear in this model (a feature we discuss in Remark 1 below).It is a simple matter to see that due to (9) we can discard equation ( 13), allowing us to reduce the number of di erential equations from the normalized SEIR model ( 10)-( 13).Now we are faced with the choice of the cost for the normalized model.Taking into account that the main aim is to control or to eliminate the disease from the population under study, di erent costs is may be considered, reflecting di erent concerns.
The choice of the cost for (P ) is based on the need to define vaccine policies to minimize some social and financial burden of the disease into the population.The cost J C for (P ) captures two easily recognized social and financial costs that the infectious disease imposes on a society; the burden of having people infected among the population, translated in J C as AI(t), and the burden of vaccination to combat the spreading of the diseases, translated in J C as Bu(t).Noteworthy, the weights A and B can be easily changed to reflect di erent scenarios.
An almost straightforward translation of this reasoning to our normalized model , where x = (s, e, i) is now the state variable and fl a parameter reflecting di erent weights on i and u.Numerically speaking, this cost is clearly di erent to J C since its two terms are now both percentages while J C has a weighted sum of two distinctive amounts.Not surprising, some simple tuning of fl may lead to similar results specially in situations where the total population N does not vary much and, in particular, when the total population remains constant.However, J n demands that the cost of infected people should be calculated in terms of the percentage of people infected instead of using the number of infected people as in J C , something that depends on the economic practices.We postpone this discussion of the introduction of di erent costs to future research and we proceed now with the cost J n above for our normalized model.Putting all together we are led to the normalized basic optimal control problem: Note that the dynamics is of the form ẋ(t) = f (x(t)) + g(x(t))u, with x = (s, e, i) and appropriate functions f (x) and g(x).
Remark 1.A word of caution regarding the way the system (10)-( 13) is viewed.
We cannot interpret the dynamics between these new compartments in the same way as with the classical model.Indeed, in equation ( 10) the term +ai(t)s(t) does not mean that those who died of the disease become susceptible as if reborn in a di erent compartment.Instead, this dynamical model can be better understood using fluid analogies in the following way.Consider a system with four tanks with the same amount of water circulating between them.At each instant t the total amount of water in the system is constant but the level of water in each tank varies.This analogy is presented in Figure 1 assuming that u(t) = 0.As mentioned before, (P n ) does not depend on the parameter d.This happens because the normal dead rate is equal for any compartment and thus it does not a ect the distribution of population into compartments.On the other hand, the birth rate works as if feeding the susceptible compartment.We will discuss pros and cons of (P ) and (P n ) using numerical solutions.Before that, however, we need to briefly present some remarks on the numerical tools and to describe the set of parameters used in our simulations.
3. Numerical Methods and Data for Simulations.Optimal control problems can be solved numerically by direct or indirect methods.Here, we opt to use the direct method (for a description these two methods see, for example, [?]): first the problem is discretized and the subsequent optimization problem is then solved using software packages with large scale nonlinear continuous optimization solvers.In this work all the simulations were made with the Applied Modelling Programming Language (AMPL), developed by [?], and interfaced to the Interior-Point optimization solver IPOPT, developed by [?].Alternatively, the optimization solver WORHP (see [?]) can also be interfaced with AMPL.We refer the reader to [?] and references within for more information on software for optimal control problems.
The application of the Maximum Principle to problems in the form of (P ) (the control appearing linearly in the dynamics and the cost and with box control constraints) yields that the the solution is a concatenation of bang-bang arcs and/or singular arcs.That is the case of (P n ) and (P ).When the optimal control is bangbang su cient second order optimality conditions(SSC) can be checked numerically; SSC as described in [?] and [?] can be checked with high accuracy using the control package NUDOCCCS.Alternatively they can also be tested with AMPL interfaced with IPOPTS or with WORHP.
In all the computations we consider the time horizon to be 20 years: thus T = 20.The parameters characterizing the population and the disease are in table 1.These ).The problem (P ) is sensitive to the dimension of the population.For example, the solution changes when we merely perturb the dead rate parameter, d, or when the initial values S 0 , E 0 , I 0 and R 0 are multiplied by a given positive constant.The first situation is illustrated in the left graph of figure 2, where the optimal control for a problem (P ) with d = 0.0099 is shown together with the optimal control for a problems with d = 0.0005; although the profile of the optimal control is the same, the switching times do change.Sensitivity of (P ) to di erent initial values is illustrated in the right graph of figure 2. For a smaller population the optimal control is of the bang-singular-bang type whereas, for a larger population it is bang-bang.In both cases, the percentages of susceptible, exposed and infected initial individuals is the same.Clearly, the cost is in the core of sensitivity of (P ) with respect to the size of the initial population and also with respect to di erent death rates.
Although the two problem (P ) and (P n ) are di erent, we can get (P n ) to produce approximately the same solutions of (P ) as mentioned above.To do so, we need to choose the parameter fl in J n to be fl = A ◊ fi B , where A and B are the weights of I(t) and u(t) and fi is an average of the total population during T .An example is shown in 3. Observe that if the total population remains almost constant, we can choose fi = N 0 .Control for a population where the order of magnitude of its size is 1000 Control for a population where the order of magnitude of its size is 100 000 Left: Optimal control di erent dead rates: in red for d = 0.0099 and in blue for d = 0.0005.Right: Optimal control with di erent initial values.
In red for S0, E0, I0 and R0 as in the table 2, in blue for initial conditions S0 ◊ 100, E0 ◊ 100, I0 ◊ 100 and R0 ◊ 100 .Clearly, the reason why (P ) and the (P n ), with parameters associated to (P ), di er resides on the cost.As mentioned before, the cost of J n , introduced in section 2 as an adaptation of J C , requires a fresh approach.Indeed, (P n ) demands that the "financial" cost be expressed in terms of percentages and not numbers of individuals.
When considering problem (P n ), we are mainly concerned on how the total population is distributed into compartments.Because of the nature of the percentages we use, both the total number of the population as well as the dead rate are simply not there.The problem (P n ) may be useful to simulate the same disease acting on di erent sized populations but with similar birth rates.Observe also, that if needed, when solving (P n ) numerically we can keep track of the changes of the total the population by adding the di erential equation Ṅ (t) = (b≠d)N (t)≠ai(t)N (t) to our code.This does not cause any change on the optimal control since the solution does not depend on N .While (P n ) has the advantage of covering di erent populations in one go, it does not always subsumes or replaces (P ).For example, the situation when one seeks vaccination policies when the total number of vaccines available are bounded (as in [?]) is tailored for the classical model (P ).Mixed constraints like those introduced in [?], i.e., of the form can be mathematically translated to normalized models but they loose their meaning.However, this drawback may be overcome by considering u(t) oe [0, ū], for a suitable contant ū < 1.This technique was also used in [?], for (P ), as numerical simulations shows that the results obtained by using this strategy are similar to the solutions obtained by using the restriction described in ( 14).
5. Solution of (P n ).We now focus on the Maximum Principle for the problem (P n ).The Hamiltonian given by for appropriated f and g, x = (s, e, i) and where p(t) = (p s (t), p e (t), p i (t)) oe R 3 denotes the adjoint variable.Let (x ú , u ú ) oe W 1,OE ◊ L OE be an optimal solution to (P n ).Then the maximum principle [?] asserts the existence of a scalar ⁄ Ø 0 and an absolutely continuous function p such that the following conditions are satisfied almost everywhere: of the Hamiltonian), (iv) p(T ) = 0 (transversality condition).
Since (P n ) does not have final time constraints on the states, it is well known that the above conditions hold with ⁄ = 1.Set "(t) = H u [t] = Èp(t), g[t]Í ≠ 1.This is called the switching function and it is of help to deduce a characterization of the optimal control u ú from (i)-(iv) which we proceed to do next.
It is a simple matter to see that condition (iii) is equivalent to It follows that u ú is bang-bang in an interval I µ [0, T ], if the switching function " has a finite number of isolated zeros at which the control switches between 0 and ū.If " is zero on an interval, then u ú is singular, i.e., In terms of the data of (P n ) the adjoint condition (ii) reads Also, we have Since our computations show that singular arc s may appear, let us assume that "(t) = 0 for t in an interval I µ [0, T ] and check if we can obtain formulas for the singular controls.We work in the region1 R := ) (s, e, i) oe R 3 : s Ø 0, e Ø 0, i Ø 0 * and so we deduce that In the interior of the singular interval we have d" dt = 0 and It is important to observe that the above expression for singular controls depends on the multipliers.Since we do not establish that the multipliers are unique, we can only expect to use (22) to validate numerical findings but not to prove to optimality.In fact, to prove optimality of computed solution we need to check numerically su cient conditions.Unfortunately, there are no numerically verifiable su cient conditions for problems with singular arcs.
6. Numerical solutions for (P n ).We now present and discuss the results of our simulations for (P n ).Recall that we use the data in table 1 and 3. We treat three di erent cases: • Case 1: ū = 1 and fl = 500, • Case 2: ū = 1 and fl = 10, • Case 3: ū = 0.2 and fl = 500.
In the first two cases the computed optimal control exhibits a bang-singular-bang structure while in the last one the optimal control is bang-bang.For all the three cases we present graphs with the computed controls and trajectories.As in [?] and to o keep the exposition short, we do not present the graphs of the multipliers but we give their computed initial values, and we also present the final states, the costs and the switching times Case 1: Taking ū to be 1 depicts the situation when all the susceptible population can be vaccinated.The results of the simulations are shown in figures 4 and 5.In Fig. 4 we show numerically that the numeric optimal control and switching function, ", satisfy ( 16) while the computed singular control ( 22) coincides with the computed optimal control u ú .The optimal trajectories are presented in Fig. 5. Numerical results for Case 1: J = 327.15,t 1 = 6.4,t 2 = 12.07, s(T ) = 0.095341, e(T ) = 0.00051104, i(T ) = 0.0020380, p s (0) = ≠126.5,p e (0) = ≠2253, p i (0) = ≠3219.Case 2: The results of the simulations are shown in figures 6 and 7.In figure 6 the optimal control, u ú , and the switching function, ", are presented.As expected, u ú = 1 and u ú = 0 when, respectively, "(t) > 0 and " < 0. The optimal trajectories for the state are presented in figure 7. Numerical results for case 2: J = 10.181,t 1 = 2.71, s(T ) = 0.16598, e(T ) = 0.0033060, i(T ) = 0.0079433, p s (0) = ≠377.0,p e (0) = ≠5137, p i (0) = ≠7081.When we go from case 1 to case 2, the optimal control goes from singular to bang-bang.This is because when we decrease the value of fl, the weight of the control in the cost increases.
Case 3: While keeping fl = 500, we now consider that only 20% of the susceptible people are eligible to be vaccinated, i.e., ū = 0.2.Results of the simulations are illustrated in figures 8 and 9.In Fig. 8 the optimal vaccinate rate, u ú , and the switching function, ", are presented.As expected, u ú (t) = ū = 0.2 and u ú = 0 when, respectively, " > 0 and " < 0. The trajectories for the state are presented in Fig. 9.The profile of the optimal control is bang-bang.The control becomes equal to zero at instant t 1 = 16.67.Numerical results for case 3: J = 806.12,t 1 = 16.67,s(T ) = 0.071623, e(T ) = 0.0031999, i(T ) = 0.015275, 6.1.Approximation for optimal control.If the control is bang-bang as in case 2 and 3 second order su cient conditions may be checked numerically as described in [?]and [?].Here we refrain from engaging in such discussion to keep the exposition short.Here we compute numerically the switching times using the so called induced optimization problem as in [?].Recall that the switching times are the points t s1 at which the optimal control changes from one bound to another.
Implementing the induced optimization problem with AMPL for case 2, with fl = 10 and ū = 1 and denoting by t s1 and J s the computed switching time and cost, we get t s1 = t 1 = 2.71 and J s = 10.183 which in contrast with J = 10.181 is a good approximation.
For case 3, with fl = 500 and ū = 0.2, we have t s1 = 16.63 < t 1 = 16 ≠ 67 while we get get a match for the cost J s = J = 806.12.We now turn to case 1 where the computed control u ú is of the form of Bang-Singular-Bang.As discussed in [?], singular controls may be hard to implemented in practice.To remedy this, [?] proposes to approximate the singular control by a constant value ũ.Here we do as in [?] to calculate the the switching times defining the singular interval t s1 and t s2 , the value 0 AE ũ AE 1 and the corresponing cost for case 1 where fl = 500 and ū = 1.Using the arc parametrization method described in [?], now implemented in AMPL code, we obtain t s1 = 7, t s2 = 11.5, ũ = 0.55 with cost J s = 326.12 a value very close to the previous J = 327.15.

Conclusion.
We studied the optimal control of an epidemiological normalized SEIR model using a L 1 -type objectives.We extracted information about the optimal solution from the Maximum Principle.In particular, we determined a closed form for the singular controls.Numerical NP solvers were applied to the discretized problem enabling us to compute optimal control solutions that match the necessary conditions, in particular, the switching conditions and forms for the singular controls.Moreover, we confronted this problem with the one previously studied in [?] where the so called classical SEIR model is used.The normalized model may cover in one single problem populations of di erent size and it is defined with what may be seen as a more realistic cost.Becuase of the use of normalized model, the solution of (P n ) is dictated by the distribution of the population into the three di erent compartments and the disease characteristics.On the other hand, the normalized model has the disadvantage of not allowing the analysis of certain real situations that cannot be translated by percentages, such as when there is a limited stock of vaccine at a non-constant population.Our discussion of these two problems highlights their di erences and the need for a criterion choice of cost, a challenging and relevant subject in optimal control for biomathematics problems.8. Acknowledgement.The authors would like to thank Prof. Helmut Maurer for numerous and enlightening discussions on this topic as well as his help in writing up AMPL codes for the problems reported here.Thanks are due to the anonymous referees whose many comments and suggestions greatly improved this paper.

Figure 1 .
Figure 1.Fluid Analogy of the Normalized SEIR compartmental model.

Figure 3 .
Figure 3. Optimal control for (Pn) with fl = 500 in blue.Optimal control calculated for (P ) in blue.The parameters are described in Tables1, 3, and 2.

Figure 4 .Figure 5 .
Figure 4. Case 1: Computed optimal control u ú plotted together with the singular control computed according to (22) and with the switching function ".During the first five years "(t) > 1 and during the last eight years "(t) < 0.

Figure 8 .Figure 9 .
Figure 8. Case 3: Computed optimal control u ú plotted together with the scaled switching function ".During the first sixteen years "(t) > 0 and during the last three years "(t) < 0.

Table 1 .
Parameters for SEIR modelsEuropean countries.The disease parameters correspond to a devastating disease as one can see either by calculating the reproduction number, which is higher than 1 (see[?]) or by solving the classical SEIR system ((1)-(6)) with u(t) = 0.In the next section, and when convenient, we use di erent parameters but this will be clearly stated.The values of the initial conditions follow those in [?]: for classical SEIR models they are presented at table 2 together with the parameters A and B of the objective functional.The initial conditions for the normalized problem (P n ) are presented at table 3.
parameters do not correspond to any specific population or diseases.In fact, f , a and g are equal to those in [?], while c is adapted from that corresponding values in[?].As for b and d, these are closely related to birth rates and death rate in

Table 2 .
Initial Conditions and cost parameters for problems with classical SEIR model

Table 3 .
Initial Conditions and cost parameters for problems with classical SEIR model normalized model.