OPTIMAL VACCINATION STRATEGIES FOR AN SEIR MODEL OF INFECTIOUS DISEASES WITH LOGISTIC GROWTH

. In this paper an improved SEIR model for an infectious disease is presented which includes logistic growth for the total population. The aim is to develop optimal vaccination strategies against the spread of a generic disease. These vaccination strategies arise from the study of optimal control problems with various kinds of constraints including mixed control-state and state constraints. After presenting the new model and implementing the op- timal control problems by means of a ﬁrst-discretize-then-optimize method, numerical results for six scenarios are discussed and compared to an analytical optimal control law based on Pontrygin’s minimum principle that allows to verify these results as approximations of candidate optimal solutions.


1.
Introduction. Infectious diseases have been feared by mankind for a long time, just remember the black death in the 14th century, killing about one third of the European population. Today humans can be vaccinated against many diseases, but still not against all of them, as can be seen from the recent outbreak of the ebola virus disease. Hence, it is hardly surprising that mathematical modelling of spreading of diseases has found much interest in the relevant literature. In this paper, we restrict ourselves to the textbooks of [15], [16], [1], and [4]. More references can be found in Biswas, Paiva and de Pinho [2], on which the present work is mainly based and the results of which we are going to extend.
As in [2], we consider a population of (N ) individuals which is divided into four disjoint groups (compartments): the susceptible population (S), the people (E), that are exposed but not carrying symptoms, the infected ones, (I), and the recovered (immunized) population (R), by either natural recovering or by vaccination. The starting point is a generic SEIR model as used in [2] which goes back to [17]. This model will be expanded with logistic population growth.
Exponential population growth is realistic for a quite young or fast increasing population, but in a highly developed country the population growth is known to stagnate more and more. In this case, the dynamical behaviour may be better represented by logistic population growth. The numerical results of this paper will show that certain unrealistic growth as shown by the results of [2] can be avoided. The basic model is assumed to have a continuous flow between the compartments according to proportionality laws as indicated in Fig. 1. In the following model the population growth is, in the first instance, modelled according to an exponential growth law, where the total population is assumed to satisfẏ while the complete SEIR model is given bẏ with proportionality parameters given in the upper block of Table 1; cp. [2], p. 11. See also Fig. 1. In this system of ordinary differential equations the variable u, denoting the vaccination rate, is a degree of freedom. Later on it will turn into the control variable. Note that the growth of the entire population N in Eq. 1 must be reduced by the disease-induced death term. A continuous flow between the compartments may not be realistic, but may be accepted for qualitative conclusions for a sufficiently large population. To indicate this, we introduce unit of capita for the real valued variables S, E, I, R, and N . Herewith, the units for all proportionality parameters, hence their meanings become apparent; see upper block of Table 1. In addition, this helps to check the equations.
The above model is exactly the one investigated in [2] and traces back to [17]. Our paper follows the organization of [2] in order to facilitate a comparison of the results between exponential and logistic growth. We firstly present the new SEIR model with logistic growth. Then a short excursus shows how analytical control laws can be obtained from Pontryagin's minimum principle. These laws are used to validate our numerical results at least approximately. Finally, numerical results for six optimal control scenarios are obtained and discussed including those with a mixed control state constraint, also called state dependend control constraint or zeroth-order state constraint, as well as a pure state constraint which turns out to initial susceptible population unit of capita 1000 E 0 initial exposed population unit of capita 100 I 0 initial infected population unit of capita 50 R 0 initial recovered population unit of capita 15 N 0 initial total population unit of capita 1165 W 0 initial vaccinated population unit of capita 0 Table 1. Values from [17], also chosen in [2] 1 be of first order, first each constraint of these constraints separately then both in combination.
In addition to [2] we investigate an optimal control problem with a discounted objective function, too.
2. The SEIR model with logistic population growth. The Belgian mathematician Verhulst published in 1838 a logistic growth law [23], which sets the growth rateṄ proportional to two terms, the population N on the one hand and the residual resources of the habitat compared with the carrying capacity K on the other hand,Ṅ with r being a generic constant in the so-called Verhulst's quadratical death term. Comparing exponential and logistic population growth, i. e., we see the following equivalences: The birth term must only appear in the susceptible population, since every one is susceptible of birth and can be infected. However, the death terms must be split up over the compartments, for exponential growth by and for logistic growth by Note that the additional parameters for the logistic growth model have been chosen according to the equivalences (8) despite the unsoundness of the birth and death rates. 2 Hence, the original SEIR model turns into the new model with the algebraic equation To complete the system of ordinary differential algebraic equations we choose appropriate (consistent) initial conditions see the lower block in Table 1.
Remark 1. In [14], a slightly modified exponential-growth model is investigated, where the term c S(t) I(t) in Eqs. (2) and (3) is replaced byc S(t) I(t)/N (t). In order not to mix up two modifications, we keep here the so-called Kermack-McKendrick model, as used in the paper [2], too.
In the Kermack-McKendrick model, the force of infection, i. e. the the probability per unit time for a susceptible to become infected, is assumed to be proportional to I, analog to the law of mass action in chemistry; see [4], p. 17. The constant c is called the transmission rate constant or incidence coefficient. Furthermore, it is assumed that the infectives have a constant probability per unit time g, to become removed. It is called the natural recovery rate g. In other words, the infectious period has an exponential distribution with parameter g, i. e. the probability to be still infectious τ Because of the subsequent continuation of the infected guests' journeys, the infection has spread out worldwide. 3 This behaviour of the spread of a disease can be compared to the motion of gas molecules in random walk.
The termc S(t) I(t)/N (t), used in [14], in contrast describes the spread of an infectious disease when the per capita number of contacts per unit of time is independent of the population size N . The force of infection here equalsc I N ; see [4], p. 22. Another interpretation is based on the conditional probabilty of contacts S N I N scaled byc N . Herec [unit of time −1 ] is a rate. Moreover, this modified model is more appropriate if the total number N of individuals varies considerably. Hence, our choice of the Kermack-McKendrick model is another concession in favour of a plain comparison of our logistic growth model with the exponential one of [2].
For the sake of comparison of the two different growth models, we therefore refrain from following the modified model of [14], moreover, since this paper deals with L 1 -optimization. 3. The optimal control model. The dynamics of the flow between the compartments are given by the above four ordinary differential equations for the state variables S, E, I and R, while u represents the only control variable. Note that the differential-algebraic system (11)-(15) is overdetermined; one of the differential equations can be dropped. The hidden differential equation for the algebraic equation (15) is automatically fulfilled by summing up equations (11)- (14). We will later come back to this.
The cost of an epidemic for an economy is assumed to consist of the costs for vaccination and for medical treatment of those who are infected. The latter term may also include the loss of benefit of an economy due to sick individuals. The vaccination costs usually amount to only a fraction of the costs for curing patients. Since this turns out to be a multi-objective performance index, we have to obey the scaling of these quantities when scalarizing the multi-objective functional; see below.
If a region with only few physicians is considered, every doctor has to medicate quite a large area and the costs may therefore depend quadratically on the vaccination rate. 4 So the cost functional can be modelled by here we investigate a period of 20 years, thus we set t 0 = 0 and T = 20 and choose the weight parameters A 1 = 0. Here we have to take into account the orders of magnitude of I and u which differ by a factor of 1 000.
Remark 2. Therefore, the term quadratic in the control plays more the role of a regularization term, to guarantee the existence of unique optimal solutions and to get solutions of higher regularity, which approximate the L 1 -optimal solutions the better the smaller the quotient A 2 /A 1 is. Another senceful term for including in the Pareto functional (17) is A 3 u(t) S(t) with A 3 meassured in units of money per unit of capita. Later we will see that this can also be taken into account by an appropriate mixed constraint such as (34).
Obviously, a box constraint to the vaccination rate must be imposed for practical reasons, say which defines the set U ad of admissible controls. Looking ahead on the subsequent paper concerning L 1 -optimization, the box constraint is a must to guarantee the existence of a solution. On subarcs with u(t) ≡ u max , we have maximum rates for the susceptible population to become immune, by vaccination and natural recovering. Furthermore, it may be interesting how many vaccines have to be used over the considered period of time. Therefore, we introduce an additional state variable We will consider values of 0 ≤ W M ≤ ∞ [unit of capita], which lead to a terminal inequality constraint. This enables either to consider a restricted amount of available vaccines or to assume to have an unlimited amount of it. As we will see, the numerical results surprisingly indicate that it may be sometimes more efficient not to use all of the available serums. This common sense contradicting result however has to be thoroughly discussed.
To complete the optimal control problem we have to mention that all other state variables are unspecified at terminal time (except with respect to a trivial nonnegativity condition, In summary, the optimal control problem reads as follows: Minimize subject toṠ Note that the state variable R does not enter the other equations. Hence, we can omit its differential equation, moreover since R can be computed a posteriori by (15). Alternatively, we also can omit the differential equation for N and substitute the algebraic equation (15) into the above system. Now N can be computed a posteriori by (15). Later, we will apply both of these alternatives. Obviously we can eliminate any one of the above differential equations and substitute the algebraic equation (15).
Remark 3. Note that the cancellation of one of the differential equations stabilizes the numerical computations. Otherwise constraint qualifications would be violated when optimizing.
4. Numerical results. In the following, we will firstly investigate five optimal control scenarios differing with respect to certain inequality constraints that are additionally imposed. These optimal control problems are solved numerically by a standard first-discretize-then-optimize (direct) method based on the modelling language Ampl of [5], providing exact derivatives via automatic differentiation, and the well-known interior point algorithm for solving large-scale nonlinear programming problems of [24], resp. [25] called IpOpt. All computations were performed using the implicit Euler method with a constant discretization stepsize of 1/(20·365) [unit of time]. The functional (21) has to be discretized suitably, i. e. compatibly with the implicit Euler method. In order to get such a quadrature formula, we transform the Lagrange functional (21) into a Mayer functional of type min u∈U ad z(T ), where z is an auxiliary state variable satisfyingż(t) = A 1 I(t) + A 2 u(t) 2 with z(0) = 0. Then the suitable quadrature formula becomes apparent and it is here the right-sided rectangular formula.
For the first four of the five following scenarios, the existence of optimal solutions can be guaranteed, as for the SEIR model with exponential growth in [2]; see Theorem 2.1 there and its application to that SEIR model. This is due to the facts, that we, too, have a Lagrangian functional with an integrand quadratically dependent on the control and a constrained set of admissible controls as well as dynamics affine in the control although nonlinear in the state variables.
In this paper we will additionally verify the discrete optimal control values approximately via an analytically determined optimal control law obtained from Pontryagin's maximum, resp. minimum principle 5 , in which we substitute the discrete state and costate variables (multipliers). This consistency condition shows how accurate the discrete solution of the finite dimensional optimization problem fulfills the optimality condition of the maximum principle. For this purpose, we compute an a posteriori error estimate defined by the mean square error MSE of all discrete data points between the discrete optimal control values and the one we obtain via Pontryagin's optimality condition as mentioned above and explained in detail below.
For recipes on how to apply of the necessary conditions of optimal control theory to real-life applications, see [18]  The overall costs are ca. 3 173.8 units, which will be the reference point for making up the balance of the gain of optimization.
All subsequent scenarios will now be optimal control problems. For this purpose we make a short aside how to apply the minimum principle.
We firstly have to define . Let x ∈ W 1,∞ (0, T ; R n ) and u ∈ L 2 (0, T ; R m ) be the state resp. control variable, f 0 (x, u) the cost functional, and f (x, u) =ẋ represent the ODE constraints. The latter functions are assumed to be sufficiently differentiable. Then the Hamiltonian function is defined as with λ 0 ∈ R and λ ∈ W 1,∞ (0, T ; R n ). In general λ 0 = 1 can be set. In few cases λ 0 = 0 occurs, the so-called anormal case. With free end conditions and given initial values for the state variable, λ 0 = 0 leads to a contradiction to the minimum principle saying that a vanishing of all multipliers λ 0 and λ can not happened simultaneously.
The core of the minimum principle says that for all t ∈ [0, T ] the optimal control u * is the global minimizer of the following finite dimensional optimization problem, pointwisely evaluated along an optimal trajectory (x * (t), λ(t)), For a precise and complete formulation of the minimum principle in the most general form see [7], Informal Theorem 4.1 and Theorem 4.2.
This control law will be used to verify the discrete numerical results x d , u d , and λ d obtained by AMPL/IpOpt via (30) by substituting x d and λ d into its righthand side, yielding u * d := u * (x d , λ d ) and computing MSE(u * d , u d ). 6 We now apply the core of the minimum principle and its associated verification to five different optimal control problems differing in the kind and number of additional inequality constraints.
In case the minimizer of the here convex optimization problem lies in the interior of the admissible set of control values, it is uniquely given by If u int (t) / ∈ (0, u max ) we have to project its values to [0, u max ]. Thus we obtain the control law u * (t) = min u max ; max 0; u int (t) .
The progress of the population, i. e., the computed candidate optimal trajectory for Scenario 1 is shown in Fig. 3.
The verification as described above yields a perfect match as also indicated by the mean square error of 3.44 · 10 −9 ; see Fig. 4.
The cost to the economy turns out to be 677.8, nearly a fifth of the cost without vaccination. At a first glance it is surprising that the vaccination is stopped before the end of the time period is reached. Usually one would keep the vaccination   rate at its maximum over the whole time period, but this would cause only slightly higher costs. We later will discuss a discounted version of the functional (21), where the same effect arises, which is typical, when functionals which should properly be formulated over an infinite horizon, are truncated by a finite time.  The most common solution may be: vaccination stays at the maximum limit until all vaccines are exhausted. This optimal control strategy produces costs to the economy amounting to 2 688.1 units.
Here, the optimal control is "almost bang-bang" i. e. has a sharp decline from the upper to the lower boundary of the admissible set; see Fig. 5. Note that due to the L 2 -functional (21) the optimal control must be continuous. Likewise all state variables are at least of class C 1 ; see Fig. 6. Note that the control u enters the differential equation for R but not for I. The verification leads here to an MSE of 1.21 · 10 −8 .
Bang-bang and even singular subarcs can only appear if an L 1 -functional is employed; see [22], i. e. if the control enters functional and right hand sides of the ODEs linearely.
In summary, it is -not surprisingly -cheaper to vaccinate the population than not immunizing them.
where V 0 > 0 [unit of capita per unit of time] has to be chosen appropriately. Since the number of susceptibles is positive this constraint is obviously equivalent to We have chosen V 0 = 125, which can be interpreted as approximately the number of immunizations that can be accomplished within eight hours, if each immunization takes not more than four minutes. In the following W M = ∞ is set, since an additional limitation of the total amount of vaccines doesn't make much sense. Again Eq. 27 is dropped.

Remark 4.
It is known that constraints of type (34) can be replaced by an additional term in the scalarized Pareto functional (42) as mentioned in Remark 2.
Herewith the multiplier µ associated with that constraint, the upper bound V 0 , and the weighting factor A 3 are related to each other and determine a certain point on the Pareto manifold.
Due to the mixed control-state constraints we have to augment the Hamiltonian as follows. Since the system is overdetermined, as mentioned above, it is sufficient to consider only four of the five ODE constraints (11) - (14) in addition to (7). Now, we choose (7), (11) -(13) besides the inequality constraint (34) for the optimal control problem describing senario 3, Here, a new multiplier µ is introduced associated with the constraint (34). Note that this multiplier must satisfy an additional necessary sign condition; see [7], Informal Theorem 4.1, or [18]. This sign condition being a part of the so-called complementarity condition, says that µ vanishes on unconstrained arcs and is nonnegative on arcs where the constraint (34), resp. (35) is active. Hence, the optimal control law is, on unconstrained arcs, identical to (33). On constrained arcs, the equation H u = 0 determines the multiplier µ while the optimal control is obtained by Due to the box constraint (18), the complete control law for scenario 3 can therefore be summarized by In Figs. 7 the switching time between the control according to (32) and (36) is marked, showing that control values in the interior of the admissible set are optimal, iff these values are less than those obtained by the competing values according to (36).
with an appropriately chosen parameter S max [unit of capita].
Following the guidlines in [18] we firstly examine the order of the state constraint and derive the candidate optimal boundary control on boundary arcs. If the state constraint (38) is active on a non-vanishing time intervall, i. e., if there holds the identity in time, a further differentiation of (39) with respect to time, while substituting the right hand sides of the ODE system (11)- (14), reveals the boundary control and therefore that the state constraint is of first order, From optimal control theory, we know that the associated costate λ S will be discontinuous at entry and/or exit points or at any point along constrained subarcs. Hence. the control will generally be discontinuous at those so-called junction points; cp. (32), (33).
Moreover, since the state constraint is of first order, effective touch points cannot occur; see [9] and [6].
Following the indirect adjoining approach of Bryson, Denham and Dreyfus [3] -see also Ref.
[18] -we augment the Hamiltonian similar to the case of mixed control-state constraints by a termμṠ withμ satisfying the same necessary sign conditions as µ for mixed control-state constraints. Here, we obtain the boundary control (40) directly fromṠ ≡ 0 whereas the multiplierμ can be obtained from the vanishing partial derivative of the modified Hamiltonian with respect to the control variable. 8 Recapitulating, we obtain the same optimal control law (37) as for the mixed control-state constraint; only replace u bound by u bound1 . The discrete approximations for the candidate optimal number of susceptibles and the associated control scenario 4 are shown in Figs. 9. 9 Numerically we have found that below a value of S max = 1 200 the manifold defined by the constraint seems to be empty, at least IpOpt was not able to compute a feasible solution. Here the vaccination rate is already at its maximum level from about 10 years on, while the number of susceptible individuals anyway rises since the population is still in its linerarly increasing phase. Hence, it is obviously not possible to tighten the constraint (38) much further.
Increasing the fineness of the discretization we see that an extremely precise resolution of the switching structure can be obtained even by a method of type First discretize then optimize as used here.
Remark 5. This accuracy could be even increased when we pass on to a switching point optimization technique; see [13]. For this postprocessing step one has to introduce the switching points as new optimization variables. These can be guessedin our example surely beyond reasonable doubt -from the First-discretize-thenoptimize approximations. For this, the control variable needs only to be discretized on subintervals where the control law depends on adjoint variables. With other words, in many subinterval the control values can be prescribed either by constants, i. e. the maximal, resp. minimal values allowed by the admissble control set U ad , or by the feedback control laws, e. g. (37) with u bound replaced by u bound1 in view of scenario 4. This approach leads to new finite dimensional optimization problems with considerably less optimization variables.
Herewith, an accuracy can be obtained which is surely beyond the accuracy of the model, but gives almost as much inside as using a First-optimize-then-discretize method based on the full minimum principle and using a multipoint boundary value solver for the resulting multi-point boundary value problem with jump conditions. 10 We will not pursue these ideas furtheron. They would be inappropriate for the more qualitative models discussed here.
The This choice results in a small time interval, approximately [7.39; 7.56], in which the first-order state constraint is active (see Fig. 11, left), and, for almost the entire remaining time, the mixed control-state constraint is active (see Fig. 10). The  . Scenario 4: Susceptible population and optimal control for S max = 1 200. The upright lines mark the switching points from unconstrained boundary arcs on state-constrained ones or vice versa. Three boundary arcs occur here; the last one is extremely short; see the zoom. Note that a touch point cannot exist here [6]. The optimal control is discontinuous at entry and/or exit points due to the jump discontinuities of λ S at these junction points; see [7]. The approximate candidate optimal control: the verification test yields 8.5 · 10 −8 , hence indicating again a perfect coincidence.
The numerical results for the susceptibles and infected are shown in Fig. 11. The mixed control-state constraint (35) is active over almost the entire time interval, only interrupted by a tiny state-constrained boundary arc associated with (38).   The overall costs amount to 2 320.5 units, which is an intermediate value of the optimal control scenarios investigated. 4.3.3. Scenario 6: A discounted version of the functional. In economical applications optimal control problems often pertain to functionals that are of discounted type, with r denoting the interest rate. Here r = 0.02 is chosen.
Remark 6. The valuation method of a discounted cash flow is used to estimate the attractiveness of an investment. Discounted cash flow analysis takes into account future cash flow projections and discounts them to a value estimate at present time. This is used to evaluate the potential for an investment. If the value obtained through discounted cash flow analysis is higher than the current cost of the investment, the investment is considered to be more advantageous.
Concerning realistic models for the spread of vaccinable diseases discounted functionals may be more appropriate. Usually functionals of type (42) arise from infinite horizon optimal control problems which are naively approximated by truncation of the infinite horizon to a finite time interval. Lykina [11] and Lykina, Pickenhain and Wagner [12] have shown that this truncation is an improper modelling since one often cannot assure the existence of solutions. This is caused by the fact that improper spaces are chosen in which the solutions are assumed to exist. To overcome this hurdle Pickenhain [19] has suggested to use a weighted Sobolev space for the state and a weighted Lebesgue space for the control variables. Moreover, these authors have developed various techniques for solving such problems correctly; see particularly [11]. Unfortunately, there is no single method of choice; it depends on the particular problem.
Let us firstly consider the approximate solution for the discounted functional with time horizon T = 20 and an unlimited amount of serums, i. e., W M = ∞. Hence, Scenario 6 is a modification of Scenario 1. The numerical results are shown in Figs. 12 and 13.
There one can hardly see any difference to Figs. 3 and 4, even not for higher (more or less realistic) values of r. Again, we see the sharp decline in the control behaviour at terminal time. If we continue to vaccinate the population on the maximum vaccination rate, the costs increases only slightly from 677.8 units to 678.1 units while the number of infected individuals decreases slightly from 1 194 to 1 184, by less than 0.1 %. However, the number of immunized individuals shows a larger difference: 5 282 (constant maximum vaccination rate) versus 4 781 individuals (optimal control on finite horizont). The effect with respect to the terminal behaviour of the approximate control may be unexpected on the first glance but is typical for optimal control problems which should be modelled on an infinite horizont but are truncated by a finite horizont; see [12], [11], [19], and [26]. This can be interpreted as devil-may-care solution.

5.
Conclusion. In the present paper a modified SEIR model for infectious diseases is developed which includes logistic growth. It is more appropriate for developed societies. The results of our paper complement the results of the twin paper by Biswas, Paiva and de Pinho [2], where exponential growth is taken into account. In particular several inequality constraints have been included, among others mixed control-state and pure state constraints. Although the numerical results have been   obtained only by a first-discretize, then-optimize method -however on the basis of an existence result -, optimality is validated by use of Pontryagin's minimum principle. This verification technique shows that the results obtained are close approximations of at least candidate optimal solutions. As conclusion concerning the numerical method used here, the investigations show that the easy-to-use direct method, particularly the combination of the modelling language AMPL providing the user with exact gradients via automatic differentiation, and the high performance large-scale NLP solver IpOpt, yields numerical results which even allow the determination of the switching structure of the problems if the switching points are not too dense. However, even in these cases a switching point optimization, if meaningful in view of the model accuracy, may yield further improvements; see [13]. This approach may also be used to improve the validation technique towards approximate sufficiency conditions for local optimality; see [13], too.
Concerning the interpretation of the results in view of the application we have seen that there is a correlation between low costs to the economy and a high vaccination rate, which is, of course, to be expected by common sense. What contradicts common sense is not to vaccine at maximal rate in the case of unlimited serums, but the optimal results are quite close to those of common practise with a 100 % rate over the entire time period exhibiting only slightly higher costs and may be caused by an improper approximation of an infinite horizont optimal control problem by a finite horizont one. For this purpose and due to economical reasons we also have investigated a discounted version of the functional proposed in [2].
Finally it should be mentioned that results have been obtained also for functionals of L 1 -type with both exponential and logistic growth leading to candidate optimal solutions with bang-bang and singular subarcs; see [22]. Such functionals may be more appropriate for realistic modelling, although the differences are minor when the weighting parameters are chosen to regularize the problem as in the present paper. The publication of L 1 -optimal solutions will follow in a subsequent paper.
Further investigations should focus on the mentioned alternative contact term for time varying populations. However, then the incidence ratec must be enlarged and must undergo a parameter study. To keep the value ofc = c as used here, would considerably diminish the influence of the contact term due to the term 1 N . Furthermore, it is important to acquire more realistic data in cooperation with edidemiologists and to investigate sensitivity studies for the parameters. In addition the SEIR systems should be rewritten in dimensionless forms which certainly would simplify the numerical treatment.
Finally, infinite horizont models could also be worth of future research, since our numerical results, particularly for the L 2 -optimization, indicate this wellknown weakness for truncated infinite horizon problems by the typical singular perturbed behaviour of solutions at final time.