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 satisfy

$\dot{N} (t) = (b - d) \cdot N(t) \label{eq_exp-growth}$

(1)

while the complete SEIR model is given by

$\dot{S} (t) = b\, N(t) - d\, S(t) - c\, S(t)\, I(t) - u(t) \, S(t),$

(2)

$\dot{E} (t) = c\, S(t)\, I(t) - (e + d)\, E(t),$

(3)

$\dot{I} (t) = e\, E(t) - (g + a + d)\, I(t),$

(4)

$\dot{R} (t) = g\, I(t) - d\, R(t) + u(t) \, S(t),$

(5)

$\dot{N}(t) = (b - d) \, N(t) - a \, I(t)$

(6)

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.

¹Note that the birth and death rates are far from reality if human populations are considered.

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 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 $\dot{N}$ 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,

$\dot{N} (t) = \frac{r}{K} \cdot N(t) \cdot \left[ K - N(t) \right] \label{eq_log_growth}$

(7)

with $r$ being a generic constant in the so-called Verhulst's quadratical death term.

Comparing exponential and logistic population growth, i.e.,

$\dot{N}_{\rm log} (t) = r \cdot N(t) - \frac{r}{K} \, N(t) \cdot N(t),\\ \dot{N}_{\rm exp} (t) = b \cdot N(t) - d \cdot N(t),$

we see the following equivalences:

$b \leftrightarrow r \;\;\text{and}\;\; d \leftrightarrow \frac{r}{K} \cdot N(t). \label{eq_equiv}$

(8)

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

$d \cdot N(t) = d \cdot \left[ S(t) + E(t) + I(t) + R(t) \right],$

(9)

and for logistic growth by

$\frac{r}{K} \, N(t) \cdot N(t) = \frac{r}{K} \, N(t) \, \left[ S(t) + E(t) + I(t) + R(t) \right].$

(10)

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.²

²The unrealistic data for the birth and death rates transfer here consequently to r and K and lead to the pronounced logistic behaviour with a considerable increase of the entire population over the chosen time interval; see Figs. 2 and 3: a price to be paid for a plain comparability between the two growth models.

Hence, the original SEIR model turns into the new model

$\dot{S} (t) = r N(t) - \frac{r}{K} N(t) S(t) - c S(t) I(t) - u(t) S(t),$

(11)

$\dot{E} (t) = c\, S(t)\, I(t) - \left[e + \frac{r}{K} \, N(t) \right]\, E(t),$

(12)

$\dot{I} (t) = e\, E(t) - \left[g + a + \frac{r}{K} \, N(t)\right] \, I(t),$

(13)

$\dot{R} (t) = g\, I(t) - \frac{r}{K} \, N(t) \, R(t) + u(t) \, S(t)$

(14)

with the algebraic equation

$N(t) = S(t) + E(t) + I(t) + R(t). \label{eq_alg}$

(15)

To complete the system of ordinary differential algebraic equations we choose appropriate (consistent) initial conditions

$S(0) = S_0, E(0) = E_0, I(0) = I_0 \ {\rm and } \ R(0) = R_0;$

(16)

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 by $\tilde c\, 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 $\tau$ units of time after infection is $\exp{- g\,\tau}$.

The Kermack-McKendrick model may be adequate, if huge crowds (demonstrations, open air concerts, sport events, etc.) are infected by so-called super-spreaders. An example is the SARS infection of sixteen guests on the nineth floor of the Metropole Hotel in Kowloon, Hong Kong, by one person on February 21,2003. Because of the subsequent continuation of the infected guests' journeys, the infection has spread out worldwide.³ This behaviour of the spread of a disease can be compared to the motion of gas molecules in random walk.

³See Super-spreader, http://en.wikipedia.org/wiki/Super-spreader, Feb. 28,2017. Other pandemic outbreaks of certain viral diseases are depicted on this web page, too.

The term $\tilde c\, 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 equals $\tilde c \, \frac{I}{N}$; see [4], p. 22. Another interpretation is based on the conditional probabilty of contacts $\frac{S}{N} \, \frac{I}{N}$ scaled by $\tilde c \, N$. Here $\tilde c$ $[\mbox{ 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.⁴ So the cost functional can be modelled by

⁴Instead of this so-called $L^2$-functional, we also have obtained results for the $L^1$-version. This may be more appropriate for developed countries having a dense patient-centered care. $L^1$-optimization leads to bang-bang and singular subarcs; see [22]. Those results will be published in a subsequent paper. Note that $L^2$-functionals are under critical discussion for biological and biomedical appications; see [21].

$\int\limits_{t_0}^{T} A_1 \cdot I(t) + A_2 \, u(t)^2 \, {\rm d}t\,; \label{eq_cost}$

(17)

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.1\,[\frac{\mbox{ unit of money}}{\mbox{ unit of capita} \cdot \mbox{ unit of time}}]$ and $A_2 = 1$ [unit of money $\cdot$ unit of time]; cp. Table 1, middle block. A smaller/larger value of $A_1$ indicates a lower/higher importance of healing costs in comparison to the vaccination costs. 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

$u \in U_{\rm ad} := \{u \colon [t_0,T] \to [0,u_{\max}]\ a.\,e.\} \label{eq_U_ad}$

(18)

which defines the set $U_{\rm 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) \equiv 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 $W(t)$ [unit of capita] defined by

$\dot{W} (t) = u(t) \cdot S(t),\ W(0) = 0,\ \;\;\text{and}\;\;\ \ W(T) \le W_M.$

(19)

We will consider values of $0 \leq W_M \leq \infty$ [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 non-negativity condition,

$\left( S(T), E(T), I(T), R(T) \right) \in \mathbb{R}_{\geq 0}^4.$

(20)

In summary, the optimal control problem reads as follows:

Minimize

$\int\limits_{0}^{T} A_1 \cdot I(t) + A_2 \cdot u(t)^2 \, {\rm d}t \label{eq_cost_func}$

(21)

subject to

$\dot{S} (t) = r N (t) - \left[ \frac{r}{K} N(t) + c I(t) + u(t) \right] S(t),$

(22)

$\dot{E} (t) = c \, S (t) I(t) - \left[e + \frac{r}{K} \, N (t) \right] E (t),$

(23)

$\dot{I} (t) = e \, E (t) - \left[g + a + \frac{r}{K} \, N (t) \right] I (t),$

(24)

$\dot{R} (t) = g \, I (t) - \frac{r}{K} N (t) R (t) + u (t) S(t),$

(25)

$\dot{N} (t) = \frac{r}{K} N (t) \left[ K - N (t) \right] - a I (t),$

(26)

$\dot{W} (t) = u (t) S (t), \label{eq_W}$

(27)

$u(t) \in [0, u_{\max}] \ {\rm a.e.}$

(28)

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 \cdot 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 \in U_{\rm ad}} z(T)$, where $z$ is an auxiliary state variable satisfying $\dot z(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⁵, 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 ${\rm 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 remarks on the interesting history of the maximum principle against the background of the beginning Cold War, see [20].

For recipes on how to apply of the necessary conditions of optimal control theory to real-life applications, see [18].

---

4.1. A limitation of the entire amount of vaccines used

A limitation of the entire amount of vaccines used can be modelled by the aforementioned terminal constraint,

$0 \leq W(T) \leq W_M. \label{eq_term_ineq}$

(29)

---

4.1.1. Scenario 0: No vaccination is used

If there is no serum available, the constraint (29) causes $W(T) = W_M = 0$ and the model allows solely an uncontrolled simulation of the unlimited spreading of the disease and its provoked costs; see Fig. 2. This scenario is considered in order to have a reference value for the yield of the optimal vaccination strategies to be developed thereafter. Here the "optimal control law" is the zero-function $u^\ast (t) \equiv 0$ and IPOPT deviates from it with ${\rm MSE}(u, u^*) = 4.57 \cdot 10^{-7}$ which is near the computational accuracy of this NLP solver.

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

Definition 4.1 (The Hamiltonian). Let $x \in W^{1,\infty}(0,T;\mathbb{R}^n)$ and $u \in L^2(0,T;\mathbb{R}^m)$ be the state resp. control variable, $f_0 (x, u)$ the cost functional, and $f(x, u) = \dot{x}$ represent the ODE constraints. The latter functions are assumed to be sufficiently differentiable. Then the Hamiltonian function is defined as

$\mathcal{H} (x, u, \lambda) = \lambda_0 \cdot f_0 (x, u) + \lambda^\top \cdot f(x, u)$

with $\lambda_0 \in \mathbb{R}$ and $\lambda \in W^{1,\infty}(0,T;\mathbb{R}^n)$.

In general $\lambda_0 = 1$ can be set. In few cases $\lambda_0 = 0$ occurs, the so-called anormal case. With free end conditions and given initial values for the state variable, $\lambda_0 = 0$ leads to a contradiction to the minimum principle saying that a vanishing of all multipliers $\lambda_0$ and $\lambda$ can not happened simultaneously.

The core of the minimum principle says that for all $t \in [0,T]$ the optimal control $u^\ast$ is the global minimizer of the following finite dimensional optimization problem, pointwisely evaluated along an optimal trajectory $(x^\ast(t),\lambda(t))$,

$u^\ast = \arg \min\limits_{u \in U_{\rm ad}} \mathcal{H} (x^\ast, u, \lambda)\,. \label{eq_MinPri}$

(30)

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 $\lambda^d$ obtained by AMPL/IPOPT via (30) by substituting $x^d$ and $\lambda^d$ into its righthand side, yielding $u^{\ast \, d}:= u^\ast(x^d,\lambda^d)$ and computing ${\rm MSE}(u^{\ast \, d},u^d)$.⁶

⁶We want to point out that the discrete multipliers provided by IPOPT as by-product may have the opposite sign ---depending whether a maximum or a minimum principle is applied ---and may need to be scaled appropriately. For this, one can exploit the homogenity of the Lagrange multipliers.

Since the Hamiltonian must be constant for autonomous problems, this property can be used to check whether the implementation of the ODEs is correct.

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.

---

4.1.2. Scenario 1: No limitation of serums

The normal case for, let's say, a vaccination against flu is that there are sufficiently many vaccines available, i. e.,

$W(T) < W_M := \infty$

or in other words there is no terminal inequality constraint for the state variable $W$, and Eq. 27 can be omitted. Nevertheless, we will keep $W$ in order to simultaneously derive the optimality condition if $W(T) < W_M < \infty$ is bounded.

Taking into account the above mentioned redundancy, the Hamiltonian function can be defined for example by

$\mathcal{H} (S, E, I, R, W, u, \lambda ) = A_1 \, I(t) + A_2 \, u(t)^2 + \\ + \lambda_S \, \dot{S}(t) + \lambda_E \, \dot{E}(t) + \lambda_I \, \dot{I}(t) + \lambda_R \, \dot{R}(t) + \lambda_W \, \dot{W}(t). \nonumber$

(31)

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

$u^\ast(t) := u_{\rm int} (t) = \frac{S(t) \cdot \left[ \lambda_S - \lambda_R - \lambda_W \right]}{2}. \label{eq_law_int}$

(32)

If $u_{\rm int}(t) \notin (0, u_{\max})$ we have to project its values to $[0,u_{\max}]$. Thus we obtain the control law

$u^\ast (t) = \min \bigg\{ u_{\max}; \, \max \big\{ 0; \, u_{\rm int} (t) \big\} \bigg\}. \label{eq_law_box}$

(33)

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 \cdot 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.⁷

⁷For techniques to deal with infinite horizon optimal control problems we refer to [8], [12], [11], and [26].

Using the optimal control approximately $22\,860$ vaccines have to be used in total.

---

4.1.3. Scenario 2: Limited serums

This scenario has actually happened in $2011$ in Germany. After flu vaccines of certain pharmaceutical companies did not get the admission, other companies were not able to compensate for it. Such a situation can be modelled, e. g., by

$W(T) < W_M := 2\,500.$

The optimal control law stays the same as in the unlimited case, but of course different state and costate values have to be inserted.

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 \cdot 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. In fact, a limit of $W_M = 3\,113$ serums leads to nearly the same costs ($3\,173.7$ units) than not vaccinating. The solution to this problem shows the same qualitative behaviour as other solutions with different values $W_M$, i. e. a terminal limitation of serums, particularly exhibiting the sharp terminal decline in the control values.

---

4.2. Scenario 3: Limited vaccines at each instant of time

In order to limit the amount of vaccines at each instant of time we introduce a mixed control-state constraint,

$0 \leq u(t) \cdot S(t) \leq V_0 \ {\rm a.e.}, \label{eq_mix_constr-basic}$

(34)

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

$0 \leq u(t) \leq \frac{V_0}{S(t)} \ {\rm a.e.} \label{eq_mix_constr}$

(35)

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 = \infty$ 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 $\mu$ 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,

$\mathcal{H} (S, E, I, N, u, \lambda, \mu ) = A_1 \, I(t) + A_2 \, u(t)^2 + \lambda_S \, \dot{S}(t) + \\ + \lambda_E \, \dot{E}(t) + \lambda_I \, \dot{I}(t) + \lambda_N \, \dot{N}(t) + \mu \, \left[ V_0 - S(t) \, u(t) \right].$

Here, a new multiplier $\mu$ 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 $\mu$ vanishes on unconstrained arcs and is non-negative 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 $\mathcal{H}_u = 0$ determines the multiplier $\mu$ while the optimal control is obtained by

$u_{\rm bound}(t) = V_0 / S(t) > 0. \label{eq_law_bound_mix}$

(36)

Due to the box constraint (18), the complete control law for scenario 3 can therefore be summarized by

$u^* (t) = \min \left\{ u_{\max}; u_{\rm bound}; \max \left\{ 0; u_{\rm int} \right\} \right\}. \label{eq_law_mix}$

(37)

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).

Overall, there is quite a low final amount of infected people, $I(20) \approx 2\,516$ vs. $\approx 2\,740$ compared to the model for a total limitation of vaccines, with using only about $30$ vaccines less; see Figs. 6 and 8. On the other hand, the costs of $2\,950.1$ units are a bit higher. Again, we see the typical decline of the control values at the end here, too.

---

4.3. Augmenting the model by state constraints

In order to calm down an epidemic or pandemia it may be advantageous to limit the number of susceptibles or infected individuals by isolation independent of costs or the influence on the other compartments.

---

4.3.1. Senario 4: State constraint on the number of susceptibles

A limitation of the number of susceptibles is acchieved by the following pointwise state constraint,

$S(t) \leq S_{\max} \label{eq_state_constr}$

(38)

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,

$S(t) \equiv S_{\max}, \label{eq_bound_arc}$

(39)

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,

$u_{\rm bound_1} (t) = \frac{r \, N(t) }{S_{\max}} - \frac{r}{K} \, N(t) + c \, I(t). \label{eq_bound_state_constr}$

(40)

From optimal control theory, we know that the associated costate $\lambda_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 $\bar\mu \, \dot{S}$ with $\bar\mu$ satisfying the same necessary sign conditions as $\mu$ for mixed control-state constraints. Here, we obtain the boundary control (40) directly from $\dot{S} \equiv 0$ whereas the multiplier $\bar\mu$ can be obtained from the vanishing partial derivative of the modified Hamiltonian with respect to the control variable.⁸

⁸Here we would like to point out that, besides the core of the minimum principle, i. e. Eq. 37, the other necessary conditions associated with the minimum principle can be checked, too, such as the sign of the multipliers associated with path constraints and the signs of the jumps occuring at the junction points of boundary arcs in case of pure state constraints; see [7], Theorem 4.1 and 4.2. See also [10].

Recapitulating, we obtain the same optimal control law (37) as for the mixed control-state constraint; only replace $u_{\rm bound}$ by $u_{\rm bound_1}$. The discrete approximations for the candidate optimal number of susceptibles and the associated control scenario 4 are shown in Figs. 9.⁹

⁹Here we have rounded the output along state-constrained boundary arcs by maximal $\frac{1}{2}$ in order to avoid small oscillations which seem to be typical for direct methods.

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 guessed ---in our example surely beyond reasonable doubt ---from the First-discretize-then-optimize 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_{\rm ad}$, or by the feedback control laws, e. g. (37) with $u_{\rm bound}$ replaced by $u_{\rm bound_1}$ 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.¹⁰

¹⁰This approach pays off only if an extremely accurate solution is desired enabling the check of additional necessary conditions. However, there is a tremendous price to be paid, since the user must have sufficient experience to implement all the necessary conditions and to handle the usually sensitive boundary-value-problem solvers such as, for example, multiple shooting codes; see [18].

We will not pursue these ideas furtheron. They would be inappropriate for the more qualitative models discussed here.

The costs for scenario 4 add up to $1\,271.2$ units, and, in addition, a quite low number of infected individuals, namely $I(20) \approx 1\,519$, can be observed.

---

4.3.2. Scenario 5: Combining mixed control-state and state constraint

For our last scenario we impose both constraints (35) and (38),

$S(t) \leq S_{\max}\;\;\text{and}\;\;0 \leq S(t) \, u(t) \leq V_0 \ {\rm a.e.}$

(41)

After some numerical experiments we have chosen the parameters

$V_0 = 400\;\;\text{and}\;\;S_{\max} = 1\,700.$

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 distinction between the boundary controls associated with the mixed constraint (35) and the pure state constraint (38) is below the resolution of a direct method and can only be detected by the verification technique; see the inserted enlargement. Only at the very beginning and close to the final time the unconstrained control (32) becomes active; see Fig. 10.

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,

$\int\limits_{t_0}^{T} \exp(-r \,t) \, \left(A_1 \cdot I(t) + A_2 \, u(t)^2\right) \, {\rm d}t, \label{eq_func_disc}$

(42)

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 = \infty$. 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 rate $\tilde c$ must be enlarged and must undergo a parameter study. To keep the value of $\tilde c = c$ as used here, would considerably diminish the influence of the contact term due to the term $\frac{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.

---