A SEIR MODEL FOR CONTROL OF INFECTIOUS DISEASES WITH CONSTRAINTS

Optimal control can be of help to test and compare different vaccination strategies of a certain disease. In this paper we propose the introduction of constraints involving state variables on an optimal control problem applied to a compartmental SEIR (Susceptible. Exposed, Infectious and Recovered) model. We study the solution of such problems when mixed state control constraints are used to impose upper bounds on the available vaccines at each instant of time. We also explore the possibility of imposing upper bounds on the number of susceptible individuals with and without limitations on the number of vaccines available. In the case of mere mixed constraints a numerical and analytical study is conducted while in the other two situations only numerical results are presented.


1.
Introduction. In the last decades mathematical models in epidemiology have been important tools in analysing the propagation and control of infectious diseases. The complexity of the dynamics of any disease dictates the use of simplified mathematical models to gain some insight on the spreading of diseases and to test control strategies ( [1]). One of the early models in epidemiology was introduced in 1927 (see [12]). Since then, many epidemic models have been derived; cf. [7], [20], [22] and [19] for some similar and recent models on such epidemic diseases.
Since the first application of optimal control in biomedical engineering around 1980s ( [16]), several vaccination strategies for infectious diseases have been successfully modelled as optimal control problems. The so called compartmental models have been used extensively in this regard. Here we concentrate on one of such models, the SEIR models. Such models are based on the division of the population under study into four compartments; an individual can either be susceptible (S), exposed to the disease but not yet infectious (E), infectious (I), or recovered (immune) (R). SEIR models can represent many human infectious diseases such as measles, pox, flu, dengue, etc, but in this paper we focus on a generic SEIR model.
Our starting point is an optimal control problem proposed by Neilan and Lenhart in [17] designed to determine a vaccination strategy over a period of time T so as to minimize a certain cost function. In [17] the propagation of a disease is controlled by a limited number of vaccines while minimizing a percentage of the overall number of dead by infection and a cost associated with vaccination. From the optimal control point of view this problem exhibits some features of interest. It is a Lagrange problem with control set constraints, an integrand cost quadratic with respect to the control and with dynamics affine with respect to the control although nonlinear with respect to some state variables.
Using the model proposed in [17] and keeping, for the sake of comparison, the same cost, we explore the introduction of non standard constraints. In contrast to [17], where the case of a limited supply of vaccines during the whole period is treated, we mainly concentrate on the case where the supply of vaccines is limited at each unit of time. In control terms we introduce a mixed state-control constraint, also known as state dependent control constraint. While our interest in this formulation resides on the study of mixed constrained optimal control problems, we believe that our problem depicts a possible situation. Even in developed regions there are restrictions on the vaccination effort dictated by difficulties associated with vaccines distribution and on facilities and qualified personal.
Taking into account that the spreading of the disease under study is based on the contacts between those susceptible (not yet vaccinated) and the infectious, it seems that a reasonable policy would be to keep an upper bound on the number of susceptible individuals during the vaccination period. To tackle such requirement we introduce state constraints, leading to a second problem, now with pure state constraints. But we go a step further and we propose a third problem with both mixed and state constraints. It is not difficult to conclude that a tight upper bound on the number of susceptible population can only be achieved when large supplies of vaccines are available. However, as mentioned before, the vaccination effort at each instant is naturally bounded by the resources and once more this implies adding mixed constraints. In summary, we propose three new problems; one with only mixed constraints, another one with only state constraints and finally one with both type of constraints.
In this paper we study analytically the first problem. We prove that this problem has a unique solution and, appealing to a variant of the maximum principle, we determine an analytical closed form for the optimal control expressed in terms of the multipliers and state variables. The problem is then solved numerically using the Imperial College London Optimal Control Software -ICLOCS -and the solution is validated by the theoretical findings. The analysis of the remaining two problems is much more complicated due to the presence of state constraints. It is then reasonable to study the problems numerically (assuming the numerical solution is indeed the solution) and decide then if an analytical study is worth the effort. We present in this paper the numerical results of those two last problems but we omit their analytical analysis.
A key ingredient in the analysis of the first problem is a variant of the maximum principle for problems with mixed constraints which will be the focus of future work. Since there is a vast literature on the maximum principle we do not make any introduction to the subject. Mixed constrained problems are amply discussed in the Russian literature. For more general problems with nonsmooth data we refer the reader to the recent books [23] and [3] and for a geometric approach to [21]. This paper is organized in the following way. In the next section we state a result that will play an important role in the analysis of our mixed constrained problem done in section 4. Section 3 includes a short description of the SEIR model used in this paper as well as a description of the optimal control problem proposed in [17]. Our mixed constrained optimal control problem is presented in Section 4 followed by the determination of the optimal solution. The two problems involving pure state constraints are then stated and discussed in section 5. Numerical results and comparison of our three constrained problems as well as that in [17] appear in section 6. The final section contains our conclusions.
We finish this section with a short description of the notation used throughout this paper. Notation. If g ∈ R m , the inequality g ≤ 0 is interpreted component-wise. The inner product of two vectors x and y in R m is written as x, y . Also | · | denotes the Euclidean norm or the induced matrix norm on R p×q . The closed ball centered at x with radius δ is denotedB(x, δ) and likewise for the open ball, regardless of the dimension of the underlying space. On the other hand,B and B denote the closed and open unit ball centred at the origin.
2. Auxiliary result. We state a result that will be of relevance in our setting. Consider the autonomous optimal control problem with mixed constraints and scalar control: x(0) = x 0 , x(T ) ∈ R n Here the state x takes values in R n while the control u is a scalar and U , the control set, is a subset of R. As for the functions we have l : R n × R n → R, L : R n × R → R, f, g : R n → R n and m : R n → R. Let R > 0 and set X :=B(0, R) ⊂ R n . Also V 0 ∈ R is a given constant.
Theorem 2.1. Assume that the data of (P) satisfies the following conditions: (a) The functions l, L, f , g and h are continuous in x; (b) There exists M > 0 such that (c) The set U is compact; (d) The function u → L(x, u) is convex for each x ∈ X and there exists a constant η such that x ∈ X, u ∈ U =⇒ L(x, u) ≥ η; (e) There is at least an admissible process (x, u) such that l(x(0), x(T )) + T 0 L(x(t), u(t)) dt is finite. Then (P) has a solution.
The above result can be stated in general terms covering problems more general than (P ). For example, it applies when the control is vector-valued instead of scalar-valued. We opt to state it in the present form since this will be enough for our purpose. Its proof is a simple adaptation of the proof of Theorem 23.11 in [3]. 3. The SEIR mathematical model. To model the progress of infectious diseases in a certain population SEIR models place the individuals into four different compartments relevant to the epidemic. Those are: • susceptible (S), An individual is in the S compartment if he/she is vulnerable (or susceptible) to catching the disease. Those already infected with the disease but not able to transmit it are called exposed. Infected individuals capable of spreading the disease are infectious and so are in the I compartment and those who are immune are in the R compartment.
Since immunity is not hereditary SEIR models assume that everyone is susceptible to the disease by birth. The disease is also assumed to be transmitted to the individual by horizontal incidence, i.e., a susceptible individual becomes infected when in contact with infectious individuals. This contact may be direct (touching or biting) or indirect (air cough or sneeze). The infectious population can either die or recover completely and all those recovered (vaccinated or recovered from infection) are considered immune.
In this four compartmental model, let S(t), E(t), I(t), and R(t) denote the number of individuals in the susceptible, exposed, infectious and recovered class at time t respectively. The total population at time t is represented by N (t) = S(t)+E(t)+I(t)+R(t). To describe the disease transmission in a certain population, let e be the rate at which the exposed individuals become infectious, g is the rate at which infectious individuals recover and a denotes the death rate due to the disease. Let b be the natural birth rate and d denotes the natural death rate. These parameters are assumed constant in a finite horizon of interest. The rate of transmission is described by the number of contacts between susceptible and infectious individuals. If c is the incidence coefficient of horizontal transmission, such rate is cS(t)I(t).
A schematic diagram of the disease transmission among the individuals for an SEIR-type model is shown in Figure 1. For more information about the model we refer the reader to [11], [17] and references within.
Now we turn to the problem of controlling the spread of the disease by vaccination. Assume that the vaccine is effective so that all vaccinated susceptible individuals become immune. Let u(t) represent the percentage of susceptible individuals being vaccinated per unit of time. Taking all the above considerations into account we are led to the following dynamical system: with the initial conditions Observe that u acts as the control variable of such system. If u = 0, then no vaccination is done and u = 1 indicates that all susceptible population is vaccinated. Using the above system Neilan and Lenhart in [17] propose an optimal control problem to determine the vaccination strategy over a fixed vaccination interval [0, T ]. The idea is then to determine the vaccination policy u so as to minimize the functional Since u is less than 1, the cost of vaccination is kept low. As for the term AI(t) in the cost care should be taken when choosing A. If A equals the parameter a (the rate of death by infection) then the idea is to minimize the number of dead by infection. Choosing A much smaller than a means that little importance is given to those who die as consequence of the disease. A larger A means that the burden of dead by infection is more important than the vaccination cost. Differing from what we do here, in [17] the rate of vaccination is assumed to take values in [0, 0.9] instead of [0, 1] to eliminate the case where the entire susceptible population is vaccinated.
The remarkable feature of Neilan and Lenhart model is that they assume that the supply of vaccines is limited. To handle such situation they introduce an extra variable W which denotes the number of vaccines used. Assuming that W M is the total amount of vaccines available during the whole period of time, the constraint can be represented bẏ With W (T ) = W M one assumes that all the vaccines should be used. We believe that it is more realistic to work with W (T ) ≤ W M . As we will see from the numerical point of view this change is complete harmless. Putting all together, with slight changes, the problem studied in [17] is then The differential equation for the recovered compartment (R) is not present here. This is due to the fact that the state variable R only appears in the corresponding differential equation and so it has no role in the overall system. Also, the number of recovered individuals at each instant t can be obtained from This problem has quadratic cost and affine dynamics with respect to u. Moreover the control set [0, 1] (or [0, 0.9] for that matter) is a convex and compact set. Such special structure permits the determination of an explicit analytical expression for the optimal control in terms of the state variables and the multipliers ( [17]). The fact that the control set is now [0, 1] instead of [0, 0.9] does not affect the analytical treatment of (P NL ).
4. Mixed constrained model. In some cases, if the vaccination interval [0, T ] is big enough, the overall limit on the vaccines may be not a reasonable option. The number of vaccines available at each instant may be limited. Also the capability to vaccinate at each unit of time may dictate the need for a constraint on the number of people vaccinated at each unit of time. To deal with this situation we propose a modification of (P NL ). We replace the overall limit of vaccines W (T ) ≤ W M by: where V 0 is an upper bound taking positive values. The inequality (8) is a mixed state control constraint, also known in the literature as state dependent control constraint. The constraint in [17] is an endpoint state constraint and so it should be satisfied only in the end point T while the mixed constraint (8) is to be satisfied at almost every instant of time during the whole vaccination program. We now consider the optimal control problem we get when in (P NL ) we replace W (T ) ≤ X by (8). We emphasize that for the sake of comparison with the results in [17] we keep the same cost. Recall that we allow the control rate to take values in [0, 1] instead of [0, 0.9] as in [17]. As far as the optimal control problem is concerned the differential equationẆ is redundant. In fact, for our new problem the variable W does not appear in neither the cost, any other differential equation nor the mixed constraint. So we remove it from our system.
To simplify the forthcoming analysis of our problem it is convenient to rewrite it in the following form: Here A is a weight parameter to be defined. Our differential equationẋ(t) = f (x(t)) + g(x(t))u(t) is affine in the control and it is nonlinear in the state x due to the term f 1 .
Before engaging into seeking a solution to (P mixed ) we need first to assert that it exists. In this respect Theorem 2.1 is of help. First observe that taking u = 0 we get an admissible solution of (P mixed ) with finite cost. Thus (e) of Theorem 2.1 is satisfied. It is a simple matter to see that both (a) and (c) of Theorem 2.1 also hold. Next we construct a set X ⊂ R 4 . This can be done in various ways. One such way is the following. For each constant value of the control u in [0, 1], determine the solution x : [0, T ] → R 4 of the differential equatioṅ that satisfies x(0) = x 0 . It can be proved that such solution always exists in our setting. Let then Z be the set of all values x(t), with t ∈ [0, T ], obtained in this way 1 . The continuity of f and g guarantees that the set Z is bounded and so we can find a R > 0 such that Z ⊂B(0, R). Set X =B(0, R). It is then a simple matter to see that (b) and (d) of Theorem 2.1 hold. Thus (P mixed ) has a solution.
Let (x * , u * ) denote such solution. To determine it we now turn to the Maximum Principle.
Define the pseudo Hamiltonian for (P mixed ) as It is a simple matter to see that the maximum principle for mixed constrained problem given as Theorem 7.1 in [4] (see also [10]) applies to our problem. We now seek its consequences. To simplify the notation in what follows we will write φ[t] to indicate that the function φ is evaluated at x * (t). We shall say that the mixed , and m(x * (t))u * (t) < V 0 for t in some intervals when t > t e and when t < t b . Then we say that t e is an exit point and t b is an entry point. We can also have points where the mixed constraint is active only at isolated points but we can discard such points since the mixed constraint is to be satisfied for almost every point. For related discussion for state constraints see [9].
Theorem 7.1 in [4] asserts the existence of an absolutely continuous function p, an integrable function q and a scalar λ ≥ 0 such that Furthermore, it is possible to prove the existence of a constant K 1 q such that for almost every t ∈ [0, T ] (see [4]). We emphasize that the inequality in (iv) holds over all controls u satisfying the mixed constraint m(x * (t))u ≤ V 0 . In (iii), N [0,1] (u * (t)) stands for the normal cone from convex analysis to [0, 1] at the optimal control u * (t) (see e.g. [2]) and it is defined as We claim that the conditions (i)-(vi) hold with λ = 0. Seeking a contradiction suppose that λ = 0. Then from (ii), we get This together with (9) yields |ṗ(t)| ≤ K|p(t)| for some K. Appealing to Gronwall's inequality (see [6] and [23]), we conclude that p(t) = 0 which is impossible in view of (i). Thus we must take λ > 0. Scalling the multipliers, we can further deduce that (i)-(vi) hold with λ = 1. In optimal control terms, this means that (x * , u * ) is normal.
We now explore the consequences of (i)-(vi) with λ = 1. To do so write p(t) = (p s , p e , p i , p n ) and x * (t) = (S * (t), E * (t), I * (t), N * (t)). Then equation (ii) reads −ṗ e (t) = −(e + d)p e (t) + ep i (t), and, for some µ(t), (iii) reduces to The Weierstrass condition (iv) reads for all u ∈ [0, 1] such that S * (t)u ≤ V 0 . We now seek some characterization of the optimal control. The information we gather here will be of importance when testing the numerical solution in Section 6. We assume without loss of generality that x * (t) > 0 everywhere.
When the mixed constraint is active we have and thus u * (t) ∈ ]0, 1[ is satisfied. In this case (15) (with µ(t) = 0) and (19) assert that If the mixed constraint is not active, then by (iv), we have q(t) = 0 and, appealing to (15), we conclude that Now suppose that the mixed constraint is active (or inactive) over a certain interval T ⊂ [0, T ]. Then both u * and q are continuous in T . Assume that t e ∈]0, T [ is any entry or exit point. Using the constancy of the pseudo Hamiltonian (18) we further deduce that u * is continuous on t e . By (15) it follows that q is also continuous at t e . The above conclusions hold when u * (t) > 0. If u * (t) = 0, then the mixed constraint should be inactive and consequently q(t) = 0. Thus we can conclude that Although (P mixed ) is an optimal control problem with mixed constraints, we get a closed form of the optimal control. As in [17], the convexity of the cost with respect to u, together with the fact that the mixed constraint is linear with respect to u and the dynamics is affine in u, plays an important role in this regard.
One last word about our problem. In view of (17), we could use second order sufficient conditions (see for example [13], [14] and [15]) via Riccati equations to test any computed solution. However this falls out of the scope of this paper and it will the focus of future work. 5. State constraints. After adding mixed constraints to the original problem we are now compelled to consider problems with state constraints. A first temptation is to replace the mixed constraints by a single state constraint Analysis of the dynamics of the disease (1) -(5) with u = 0, shows (see next section) that the increase of I is small after the outburst of the disease. This constraint can be hard to implement and so it is discarded. Since the spreading of the disease is given by cS(t)I(t), it is, however, reasonable to keep an upper bound on the number of susceptible individuals. The number of susceptible individuals will certainly increase given that any newborn is considered susceptible. However, after vaccination a susceptible individual becomes immune. Thus an upper bound on the number of susceptible individuals will imply the need for an increase of vaccination. The translation in mathematical terms of the upper bound on the number of susceptible individuals is the state constraint Another possibility is to add the term S(t) to our integrand cost. That would be a soft constraint (with respect to soft and hard constraints see [5]). Since the horizon is assumed big we opt for hard state constraints so as to keep the same cost allowing comparison with previous cases. With the notation of the previous section, this yields the optimal control problem (P state ) A special feature of our problem is that the state constraint is of order one (see [9] and references [M5] and [M8] within). Nevertheless, the analysis of the Maximum Principle for optimal control problems with state constraints is demanding due to the presence of measures as multipliers associated with the state constraint (see also, and for example, [21] or [23]). We then leave out any analytical analysis (this will be future work) and we shall concentrate only on numerical simulations. Clearly, lacking an analytical analysis any such results should be viewed with care since there is no guarantee that they are optimal.
Taking into account that the vaccination effort is, at each instant, bounded by the resources we also dare to consider another problem, this one with both mixed and state constraints: This problem reflects the need to find vaccination strategies that would keep the number of susceptible individuals low when the number of vaccines available at each instant is limited. As with (P state ), we omit here any analytical study of the solution of (P SM ); we only discuss the results of numerical simulations, a task conducted in the next section. 2 6. Numerical results. We now turn to the numerical treatment and we compare the solutions of the various problems under consideration.
In Table 1 we present the values of the parameters and constants used in all our simulations. Such values are as in [17]. Observe that the initial values of the state variables S 0 , E 0 , I 0 , N 0 and W 0 appear in the last five lines of Table 1. To do our simulations we use the Imperial College London Optimal Control Software -ICLOCS -version 0.1b ( [8]). ICLOCS is an optimal control interface, implemented in Matlab, for solving optimal control problems with general path and boundary constraints and free or fixed final time. ICLOCS uses the IPOPT -Interior Point OPTimizer -solver which is an open-source software package for large-scale nonlinear optimization [24]. Our choice of software is not random. It was chosen for it simplicity based on a study conducted using different solvers for a number of optimal control problems in [18].
As in [17], and since any model chosen for the real situation is only an approximation of reality, we treat the state variables S, E, I, N and W as continuous variables instead of integers. In what follows, and in most of the cases, we show their computed values rounded to 0 decimals unless we need to compare values of approximate magnitude.
Considering a time interval of 20 years (T = 20), a time-grid with 10000 nodes was created, i.e., for t ∈ [0, 20] we get ∆t = 0.002. Since our problem is solved by a direct method and, consequently, an iterative approach, we impose an acceptable convergence tolerance, at each step, of ε rel = 10 −9 . In this respect, we refer the reader to [8].
Without Control (WoC). We start by considering (P NL ) when the control is assumed to be 0 over the whole time interval [0,20]. This means that we want to see how the disease behaves in the absent of any vaccination policy. This is equivalent to solving the the system of differential equations consisting of (1), (2), (3) and (5) when the control is set to 0. For the sake of simplicity, we present only the evolution of I and N in Figure 3. Without any preventive control, the number of Unlimited supply of vaccines (UV). Let us now turn to (P NL ) considering W M = +∞, i.e., we treat the case where unlimited number of vaccines is available over the period of time. As mentioned before, in this case we analyse a problem that differs from an analogous problem in [17] in so far as u takes values in [0, 1] instead of [0, 0.9]. The results obtained in this case are shown in Figure 4 and do not differ much from those in [17] albeit our modification. It is a simple matter to see that necessary conditions for this problem hold with λ = 1 and they assert that the optimal control is (see [17]) IOPTS gives us the values of the computed multipliers. Taking into account the numerical values of S * (t) and the associated multipliers p s (t), we can see that the numerical values of the optimal control do coincide with the values of (23) (see Figure 5).
Overall limited supply of vaccines (LV). We now consider the case where W M = 2500 in (P NL ). This means we do have a limit number of 2500 vaccines to be used during the 20 years interval. Note that here we impose the end state W (20) to be less or equal to 2500, while in [17] the equality is imposed. Thus we allow the possibility of having no used vaccines in the end of 20 years. It turns out that this slight difference between our model and that in [17] does not translate into a different solution since the optimal solution dictates that the whole amount of vaccines should be used way before the end of the 20 years, as shown in Figure  6.  Limited supply of vaccines at each instant of time (MC). Now we solve problem (P mixed ). Recall that we impose the constraint S(t)u(t) ≤ 125 during the whole vaccination period. As mentioned before, there is no need to consider the differential equationẆ (t) = S(t)u(t) in our model. However, since we want to compare the results with those of case (UV) and (LV) we include such equation when running the solver. The results for the optimal vaccination schedule is presented in Figure 7 and the total number of vaccines used over the period is 2 330, a value less than 2 500, imposed in the case (LV). Analysis of the computed control shows that it takes values in the interior of the set [0, 1], approaching 0 when t tends to 20. The mixed constraint is active during most of the period having an exit point around t = 18. Figure 7. The optimal trajectories and optimal vaccination rate for (P mixed ).
The number of vaccinated people remains below 2500 at the end of the vaccination program. As expected, and in contrast to the two previous cases, the optimal control (rate of vaccination) never achieves it maximum 1. This is not due to any smoothing effect of the mixed constraint. In fact, the reduce number of vaccines at each instant prevents the use of the maximum rate of vaccination. As illustrated in Figure 10, the optimal control dictates that the mixed constraint should remain active during almost all the period of vaccination; it only becomes inactive in the very end, dropping to zero at t = 20. The lack of vaccines available is also responsible for a cost slightly higher than that obtained before (see Table 2 bellow), a situation to which the weight of the variable I in the cost is not strange.
Since the cost is inferior to that of the previous case, one might be tempted to discard such policy. However, (P mixed ) is of interest since it can depict what happens in some situations and it can alert to the need of having large supplies of vaccines available in the early stages of the outburst of the infection.
To test the validation of the numerical values we show in Figure 8 the graphs of the computed (or numerical) control and that of the control defined by (22) showing that we have a match of these two values. Figure 8. Computed or numerical control and the "analytical" control for (P mixed ).
However, we go a step further with our validation of the optimal solution, and we check if our conclusion in Section 4, for the mixed constraint multiplier, q are satisfied. To do so, we get the computed multipliers from IPOPT. Figure 9 shows the computed values of the adjoint multipliers (i.e. p s , p e , p i and p n ) calculated by IPOPT and the values of the multiplier q associated with the mixed constraint. At T = 20 the adjoint multipliers are all 0 in agreement with conclusion (vi) of the necessary conditions.
In Figure 10 we show the values of the mixed constraint S * (t)u * (t), as well as the values of what we call the analytical mixed multiplier and the computed values of q, during the whole period. The analytical mixed multiplier is given by (20) (during the whole period). As we know, the mixed multiplier q should be 0 from the moment the mixed constraint becomes inactive. And, in contrast with the analytical multiplier, this is what happens with the numerical (or computed) multiplier given by IPOPT, as illustrated in the right hand graph of Figure 10.
It is no surprise that the mixed constraint remains active most of the time. What one may ask is why it is not active during the whole period. This is due to the cost. In fact, solving our differential equations with u(t) = V 0 /S(t) we get a slightly higher cost (33.67 as compared to 33.66) without any significant reduction of I (20).
Finally, we have also run the solver using different values for V 0 . The pattern for the mixed constraint, the optimal control as well as that of the multiplier q is as when V 0 = 125. The greater the value of V 0 , the smaller is the exit time of the Figure 9. The computed adjoint multipliers and the mixed constraint multiplier for (P mixed ).
mixed constraint, as one would expect. We got a cost similar to that of case (UV) for V 0 = 700 but with a smaller I (20). Figure 10. The mixed constraint together and the (scaled) analytical (satisfying (20)) and computed mixed multipliers for (P mixed ).
Constraint on the number of susceptible individuals during the whole interval (SC). We now show the numerical simulations of (P state ). Theoretical analysis of the data obtained will be the focus of future work. Obviously without some analytical test, the numerical results we show next should be seen with some scepticism. We consider the state constraint S(t) ≤ S max with S max = 1100. Figure 11. The optimal trajectories and optimal vaccination rate for (P state ).
The solver ICLOCS has been quite quick until now. In this case, and not surprisingly, it took 89 iterations to calculate what may be the optimal solution. A remarkable feature of results is that, by constraining S, we are able to keep the values of infectious individuals low at the expenses of a very high number of the vaccinated people. This reflects into a low cost (in this case the cost is 24.7) due to the weight of I. To achieve these values a massive program of vaccination is needed. In fact 6 345 individuals are vaccinated during the whole period. As can be seen in Figure 11, the numerical optimal control is 1 in the beginning dropping to approximately 0.2 and increasing from then on. Such increase may be explained by the need to keep the number of susceptible individuals equal or below 1 100.
An important feature of this problem is the fact that the state constraint is active at the end point (i.e., S(20) = 1 100). Observe that the multiplier associated with the S variable, p s is not 0 when T = 20, as shown in Figure 12. This behaviour can be explained by the necessary conditions but, as mentioned, such study is out of the scope of this paper and it will be the focus of future work.
Constraint on the number of susceptible individuals during the whole interval with limited limited supply of vaccines at each instant (SM). The results obtained in (SC) alone justify our next case. The question is obviously if we can keep the number of susceptible individuals bounded if the number of vaccines available at each instant is limited? Thus we solve problem (P SM ) numerically. Once again, no analytical analysis is presented. The purpose of these numerical simulations is thus to check if it is worth to overtake an extensive study of such problems.
Here, as before, we define S max = 1 100. After some preliminary simulations we opt to consider V 1 = 400 (recall this is the maximum number of vaccines available at each instant) instead of 125 as in (SM). With a smaller number of vaccines ICLOCS was not able to calculate an "optimal" solution after 200 iterations. The graphs of the states and the control computed are presented in Figure 13.
The end state values (see Table 2) do not differ much from the previous case. However the vaccination effort over the whole period is slightly less with a total of 6 193 vaccinated (recall that now we allow a 400 individuals to be vaccinated in contrast with 125 in the (MC) case). The cost is less than in the case (SC): it is now 25.7 and it is found after 51 iterations.
It is of interest to plot the graph of mixed constraint S * (t)u * (t) against that of S * (t). This can be seen in Figure 14. Observe that the mixed constraint is scaled. Indeed, to help the illustration we multiply S * (t)u * (t) by 2.5. As it can be seen the mixed constraint is active at the beginning and at the end. Observe that the instant when the mixed constraint changes from active to inactive roughly coincide with the instant where the state S takes it minimum value. Another feature of interest is the fact that the mixed constraint becomes again active when the state Figure 13. The optimal trajectories and optimal vaccination rate for (P SM ). constraint turns inactive. Note also that the state constraint is active at the end point. As expect we have p s (20) = 0 while all the other adjoint multipliers are 0 when T = 20. To keep the exposition short, and because here we refrain from presenting an analysis of the numerical finding, we do not present the graphs of the multipliers for this problem.
Brief Comparison. In Table 2 we summarize our findings. In Figure 15 we compare the evolution of the infectious class in all the cases studied. Figure 15. A comparison of the evolution of infectious population over time.
In Table 2 we show the cost, number of iterations needed to compute the numerical solutions as well as the values of all the states at the end point T = 20. The costs of both the (LV) case (that in [17]) and the (SC) case are the best. However the number of vaccines used (W (20)) during the whole period for the state constrained problems (SC) and (SM) seems to be unrealistic. This happens because of the small weight of the vaccination cost in the cost. The mixed constraint (MC) case has the larger cost but also the smaller number of vaccines. This means that the price to pay for limitation on the vaccines per unit of time is an increase in the number of infectious individuals. If what we look for is a vaccination policy with the smaller cost (and if the cost chosen is adequate), then the (LV) case produces the most attractive values pointing to the need to vaccinate intensively in the early moments of the disease.
It is clear from the comparison of the values of I(20) for the different cases studied here (see Table 2) that (UV), (SC) and (SM) are effective as far as reducing the number of infectious individuals at the end of the program. This information may be of value in some situations (for example, to prevent any new outburst of the disease). Such reduction is obtained with a high number of vaccines.
Finally, for reasons of numerical comparison, we present in Table 3 the initial values of the adjoint multipliers.

7.
Conclusion. Optimal control theory in the form of maximum principle was applied to obtain optimal vaccination schedules and control strategies for an SEIR epidemic model of human infectious diseases. We first propose an extended model subject to mixed constraint to study the case where the vaccination process is bounded at each instant t. The mixed constraint is enforced during the whole process of a vaccination program. For such problem we get a closed form of the  Table 3. Initial values of the adjoint multipliers.
Cases p s (0) * 10 −3 p e (0) * 10 −3 p i (0) * 10 −3 p w (0) * 10 −3 p n (0) * 10 −3 optimal control problem. We also propose two models with state constraints but no analytical analysis is conducted. The relation between intensive vaccination and low cost is clear and expected. Our findings indicate that if a limit on vaccines is to be considered the model with best cost values is indeed that considered in [17]. The behaviour of the solutions to the problems with merely mixed constraints and that of the (LV) case indicates the need to provide vaccines in the early stages of the diseases. The pure state constrained cases (SC) and (SM) produce what may be seen as unrealistic solutions since they imply intensive vaccination. However, together with (UV), these two cases generate the smaller values for I (20), a criterion that may be of importance.
It is our belief that all these constrained problems should be studied considering different (and possibly more realistic) costs. Also, since any approach will depend on the parameters, tests of the numerical solution using second order sufficient conditions should be performed so that the sensitivity analysis is studied.