Optimal control of non-autonomous SEIRS models with vaccination and treatment

TL;DR

This paper investigates optimal control strategies for a non-autonomous SEIRS epidemic model incorporating vaccination and treatment, providing theoretical existence results and simulation comparisons between autonomous and periodic models.

Contribution

It introduces a novel optimal control framework for non-autonomous SEIRS models with general incidence functions, including existence, uniqueness, and simulation analysis.

Findings

Existence and uniqueness of optimal control solutions established.

Simulation results compare autonomous and periodic models.

Controlled models show improved epidemic outcomes.

Abstract

We study an optimal control problem for a non-autonomous SEIRS model with incidence given by a general function of the infective, the susceptible and the total population, and with vaccination and treatment as control variables. We prove existence and uniqueness results for our problem and, for the case of mass-action incidence, we present some simulation results designed to compare an autonomous and corresponding periodic model, as well as the controlled versus uncontrolled models.

Tables

Table 1. Table 1. Values of the parameters for problem ( P ) used in Section  6 .

Name	Description	Value
$S_0$	Initial susceptible population	0.98
$E_0$	Initial exposed population	0
$I_0$	Initial infective population	0.01
$R_0$	Initial recovered population	0.01
$\mu$	natural deaths	0.05
$\epsilon$	infectivity rate	0.03
$\gamma$	rate of recovery	0.05
$\eta$	loss of immunity rate	0.041
$k_1$	weight for the number of infected	1
$k_2$	weight for treatment	0.01
$k_3$	weight for vaccination	0.01
${\tau }_{\max }$	maximum rate of treatment	0.1
$v_{\max }$	maximum rate of vaccination	0.4

Equations

$$\begin{array}{c}\mathcal{J}\left(I,\text{mathbbm}T,\text{mathbbm}V\right)={\int }_0^{t_f}\left({\kappa }_{1}I(t)+{\kappa }_2\text{mathbbm}{T}^2(t)+{\kappa }_3\text{mathbbm}{V}^2(t)\right)dt⟶\min ,\\ \{\begin{array}{l}{S}^{\prime }(t)=\Lambda \left(t\right)-\beta \left(t\right)\phi (S(t),N(t),I(t))-\mu \left(t\right)S(t)+\eta \left(t\right)R(t)-\text{mathbbm}V(t)S(t),\\ E^{\prime }(t)=\beta \left(t\right)\phi (S(t),N(t),I(t))-\left(\mu \left(t\right)+\epsilon \left(t\right)\right)E(t),\\ I^{\prime }(t)=\epsilon \left(t\right)E(t)-\left(\mu \left(t\right)+\gamma \left(t\right)\right)I(t)-\text{mathbbm}T(t)I(t),\\ R^{\prime }(t)=\gamma \left(t\right)I(t)-(\mu \left(t\right)+\eta \left(t\right))R(t)+\text{mathbbm}T(t)I(t)+\text{mathbbm}V(t)S(t),\end{array}\\ (S(0),E(0),I(0),R(0))=(S_0,E_0,I_0,R_0),\end{array}$$

$$\begin{array}{c}J(x,u)={\int }_{t_0}^{t_1}\mathcal{L}(t,x(t),u(t))dt⟶\min ,\\ \{\begin{array}{l}{x}^{\prime }(t)=f\left(t,x(t),u(t)\right),\text{a.e.}t\in [t_0,t_1],\\ x(t_0)=x_0,\end{array}\\ x(\cdot )\in AC\left([t_0,t_1];{\mathbb{R}}^n\right),u(\cdot )\in L^1([t_0,t_1];U\subset {\mathbb{R}}^m).\end{array}$$

$$N(t)=N_0{e}^{{\int }_0^{t_f}\mu (s)ds}+{\int }_0^{t_f}\Lambda (u)e^{{\int }_u^{t_f}\mu (s)ds}du\le K,$$

$$S(t),E(t),I(t),R(t)\le K.$$

$$\{(x,y,z,w)\in ({\mathbb{R}}_0^{+}{)}^4:0\le x+y+z+w\le K\}$$

$$\begin{array}{rl}\mathcal{H}& (t,(S,E,I,R),(p_1,p_2,p_3,p_4),(\text{mathbbm}T,\text{mathbbm}V))={\kappa }_{1}I+{\kappa }_2\text{mathbbm}{T}^2+{\kappa }_3\text{mathbbm}{V}^2\\ & +p_1\left[\Lambda \left(t\right)-\beta \left(t\right)\phi (S,N,I)-\mu \left(t\right)S+\eta \left(t\right)R-\text{mathbbm}VS\right]\\ & +p_2\left[\beta \left(t\right)\phi (S,N,I)-\left(\mu \left(t\right)+\epsilon \left(t\right)\right)E\right]\\ & +p_3\left[\epsilon \left(t\right)E-\left(\mu \left(t\right)+\gamma \left(t\right)\right)I-\text{mathbbm}TI\right]\\ & +p_4\left[\gamma \left(t\right)I-\mu \left(t\right)R-\eta \left(t\right)R+\text{mathbbm}TI+\text{mathbbm}VS\right].\end{array}$$

$$\{\begin{array}{l}{p}_1^{\prime }=(p_1-p_2)\beta \left(t\right)\left({\partial }_1\phi \left(S,N,I\right)+{\partial }_2\phi \left(S,N,I\right)\right)+p_1\left(\mu \left(t\right)+\text{mathbbm}V\right)-p_4\text{mathbbm}V,\\ p_2^{\prime }=(p_1-p_2)\beta \left(t\right){\partial }_2\phi \left(S,N,I\right)+p_2\left(\mu \left(t\right)+\epsilon \left(t\right)\right)-p_3\epsilon \left(t\right),\\ p_3^{\prime }=p_3\left(\mu \left(t\right)+\gamma \left(t\right)+\text{mathbbm}T\right)+\left(p_1-p_2\right)\beta \left(t\right)\left({\partial }_2\phi \left(S,N,I\right)+{\partial }_3\phi \left(S,N,I\right)\right)\\ -p_4\left(\gamma \left(t\right)+\text{mathbbm}T\right)-{\kappa }_1,\\ p_4^{\prime }=(p_1-p_2)\beta \left(t\right){\partial }_2\phi \left(S,N,I\right)+\mu \left(t\right)p_4-\eta \left(t\right)p_1+\eta \left(t\right)p_4,\end{array}$$

$$p_1(t_f)=p_2(t_f)=p_3(t_f)=p_4(t_f)=0.$$

$$\text{mathbbm}{T}^{\ast }=\min \left\{\max \left\{0,\frac{I^{\ast }(p_3-p_4)}{2k_2}\right\},{\tau }_{\max }\right\}$$

$$\text{mathbbm}{V}^{\ast }=\min \left\{\max \left\{0,\frac{S^{\ast }(p_1-p_4)}{2k_3}\right\},v_{\max }\right\}.$$

$$\frac{\partial \mathcal{H}}{\partial \text{mathbbm}{V}^{\ast }}=-p_{1}S+p_{4}S+2k_3\text{mathbbm}{V}^{\ast }=0\text{and}\frac{\partial \mathcal{H}}{\partial \text{mathbbm}{T}^{\ast }}=-p_{3}I+p_{4}I+2k_2\text{mathbbm}{T}^{\ast }=0$$

$$\text{mathbbm}{V}^{\ast }=\frac{(p_1-p_4)S}{2k_3}\text{and}\text{mathbbm}{T}^{\ast }=\frac{(p_3-p_4)I}{2k_2}.$$

$$\frac{\partial \mathcal{H}}{\partial \text{mathbbm}{V}^{\ast }}=-p_{1}S+p_{4}S+2k_3\text{mathbbm}{V}^{\ast }\le 0\iff \frac{(p_1-p_4)S}{2k_3}\ge {\nu }_{\max }.$$

$$\frac{\partial \mathcal{H}}{\partial \text{mathbbm}{T}^{\ast }}=-p_{3}I+p_{4}I+2k_2\text{mathbbm}{T}^{\ast }\le 0\iff \frac{(p_3-p_4)I}{2k_2}\ge {\tau }_{\max }.$$

$$\frac{\partial \mathcal{H}}{\partial \text{mathbbm}{V}^{\ast }}=-p_{1}S+p_{4}S+2k_3\text{mathbbm}{V}^{\ast }\ge 0\iff \frac{(p_1-p_4)S}{2k_3}\le 0.$$

$$\frac{\partial \mathcal{H}}{\partial \text{mathbbm}{T}^{\ast }}=-p_{3}I+p_{4}I+2k_2\text{mathbbm}{T}^{\ast }\ge 0\iff \frac{(p_3-p_4)I}{2k_2}\le 0.$$

$$S(t)=e^{\alpha t}s(t),E(t)=e^{\alpha t}e(t),I(t)=e^{\alpha t}i(t),R(t)=e^{\alpha t}r(t)$$

$$p_1(t)=e^{-\alpha t}{\phi }_1(t),p_2(t)=e^{-\alpha t}{\phi }_2(t),p_3(t)=e^{-\alpha t}{\phi }_3(t),p_4(t)=e^{-\alpha t}{\phi }_4(t).$$

$$N(t)=S(t)+E(t)+I(t)+R(t)=e^{\alpha t}n(t).$$

$$\begin{array}{rl}|\phi (A,B,C)-\phi & (\bar{A},\bar{B},\bar{C})|\\ & \le |\phi (A,B,C)-\phi (\bar{A},B,C)|+|\phi (\bar{A},B,C)-\phi (\bar{A},\bar{B},C)|\\ & +|\phi (\bar{A},\bar{B},C)-\phi (\bar{A},\bar{B},\bar{C})|\\ & \le M_1^u|A-\bar{A}|+M_2^u|B-\bar{B}|+M_3^u|C-\bar{C}|,\end{array}$$

$$M_i^u:=\underset{x\in \Gamma }{\sup }|{\partial }_i\phi |<+\infty ,i=1,2,3.$$

$$\alpha e^{\alpha t}s+e^{\alpha t}\dot{s}=\Lambda -\beta \phi (e^{\alpha t}s,e^{\alpha t}n,e^{\alpha t}i)-\mu e^{\alpha t}s+\eta e^{\alpha t}r-\text{mathbbm}V{e}^{\alpha t}s$$

$$\alpha s+\dot{s}=\frac{\Lambda }{e^{\alpha t}}-\frac{\beta }{e^{\alpha t}}\phi (e^{\alpha t}s,e^{\alpha t}n,e^{\alpha t}i)-\mu s+\eta r-\text{mathbbm}Vs.$$

$$\begin{array}{rl}\alpha (s-\bar{s})+(\dot{s}-\dot{\bar{s}})& =-\frac{\beta }{e^{\alpha t}}(\phi (e^{\alpha t}s,e^{\alpha t}n,e^{\alpha t}i)-\phi (e^{\alpha t}\bar{s},e^{\alpha t}\bar{n},e^{\alpha t}\bar{i}))-\mu (s-\bar{s})\\ & +\eta (r-\bar{r})-(\text{mathbbm}Vs-\overline{\text{mathbbm}V}\bar{s}).\end{array}$$

$$\begin{array}{rl}& \frac{1}{2}(s(T)-\bar{s}(T){)}^2+\alpha {\int }_0^T(s-\bar{s}{)}^{2}dt\\ & =-{\int }_0^T\frac{\beta }{e^{\alpha t}}(s-\bar{s})(\phi (e^{\alpha t}s,e^{\alpha t}n,e^{\alpha t}i)-\phi (e^{\alpha t}\bar{s},e^{\alpha t}\bar{n},e^{\alpha t}\bar{i}))dt\\ & -{\int }_0^T\mu (s-\bar{s}{)}^{2}dt+{\int }_0^T\eta (s-\bar{s})(r-\bar{r})dt-{\int }_0^T(\text{mathbbm}Vs-\overline{\text{mathbbm}V}\bar{s})(s-\bar{s})dt\end{array}$$

$$\begin{array}{rl}& \frac{1}{2}(s(T)-\bar{s}(T){)}^2+\alpha {\int }_0^T(s-\bar{s}{)}^{2}dt\\ & \le {\int }_0^T\frac{\beta }{e^{\alpha t}}|s-\bar{s}|(M_1^u|e^{\alpha t}s-e^{\alpha t}\bar{s}|+M_2^u|e^{\alpha t}n-e^{\alpha t}\bar{n}|+M_3^u|e^{\alpha t}i-e^{\alpha t}\bar{i}|)dt\\ & -{\int }_0^T\mu (s-\bar{s}{)}^{2}dt+{\int }_0^T\eta (s-\bar{s})(r-\bar{r})dt-{\int }_0^T(\text{mathbbm}Vs-\overline{\text{mathbbm}V}\bar{s})(s-\bar{s})dt\\ & ={\int }_0^T\beta M_1^u(s-\bar{s}{)}^{2}dt+{\int }_0^T\beta M_2^u|s-\bar{s}||n-\bar{n}|dt+{\int }_0^T\beta M_3^u|s-\bar{s}||i-\bar{i}|dt\\ & -{\int }_0^T\mu (s-\bar{s}{)}^{2}dt+{\int }_0^T\eta (s-\bar{s})(r-\bar{r})dt-{\int }_0^T(\text{mathbbm}Vs-\overline{\text{mathbbm}V}\bar{s})(s-\bar{s})dt\\ & \le C_1{\int }_0^T(s-\bar{s}{)}^2+(i-\bar{i}{)}^2+(e-\bar{e}{)}^2+(r-\bar{r}{)}^{2}dt+K_1{e}^{\alpha T}{\int }_0^T(\text{mathbbm}V-\overline{\text{mathbbm}V}{)}^{2}dt,\end{array}$$

$${\left(\text{mathbbm}V-\overline{\text{mathbbm}V}\right)}^2\le C_9{e}^{2\alpha T}\left[(s-\bar{s}{)}^2+({\phi }_1-{\bar{\phi }}_1{)}^2+({\phi }_4-{\bar{\phi }}_4{)}^2\right],$$

$${\left(\text{mathbbm}T-\overline{\text{mathbbm}T}\right)}^2\le C_{10}{e}^{2\alpha T}\left[(i-\bar{i}{)}^2+({\phi }_3-{\bar{\phi }}_3{)}^2+({\phi }_4-{\bar{\phi }}_4{)}^2\right],$$

$$\Psi (t)=(s(t)-\bar{s}(t){)}^2+(e(t)-\bar{e}(t){)}^2+(i(t)-\bar{i}(t){)}^2+(r(t)-\bar{r}(t){)}^2$$

$$\Phi (t)=({\phi }_1(t)-{\bar{\phi }}_1(t){)}^2+({\phi }_2(t)-{\bar{\phi }}_2(t){)}^2+({\phi }_3(t)-{\bar{\phi }}_3(t){)}^2+({\phi }_4(t)-{\bar{\phi }}_4(t){)}^2.$$

Full text

Mateus was partially supported by FCT through CMA-UBI (project UID/MAT/00212/2013), Rebelo by FCT through CMA-UBI (project UID/MAT/00212/2013), Rosa by FCT through IT (project UID/EEA/50008/2013), Silva by FCT through CMA-UBI (project UID/MAT/00212/2013), and Torres by FCT through CIDMA (project UID/MAT/04106/2013) and TOCCATA (project PTDC/EEI-AUT/2933/2014 funded by FEDER and COMPETE 2020).

∗ Corresponding author: delfim@ua.pt

Joaquim P. Mateus

Research Unit for Inland Development (UDI)

Polytechnic Institute of Guarda, 6300-559 Guarda, Portugal

Centro de Matemática e Aplicações da Universidade da Beira Interior (CMA-UBI)

Departamento de Matemática

Universidade da Beira Interior, 6201-001 Covilhã, Portugal

Paulo Rebelo

Centro de Matemática e Aplicações da Universidade da Beira Interior (CMA-UBI)

Departamento de Matemática

Universidade da Beira Interior, 6201-001 Covilhã, Portugal

Silvério Rosa

Departamento de Matemática and Instituto de Telecomunicações (IT)

Universidade da Beira Interior, 6201-001 Covilhã, Portugal

César M. Silva

Centro de Matemática e Aplicações da Universidade da Beira Interior (CMA-UBI)

Departamento de Matemática

Universidade da Beira Interior, 6201-001 Covilhã, Portugal

Delfim F. M. Torres∗

Center for Research and Development in Mathematics and Applications (CIDMA)

Department of Mathematics, University of Aveiro, 3810-193 Aveiro, Portugal

1. Introduction

Frequently, decision makers must balance the effort put in treatment and vaccination in order to give the best response to outbreaks of infectious diseases. Optimal control is an important mathematical tool that can be used to find the best strategies [16, 23, 27]. Real-world problems, under recent investigation with optimal control techniques, include Ebola [2] and Zika [6].

In [10], Gaff and Schaefer studied three distinct autonomous epidemic models, having vaccination and treatment as control variables. They established existence and, in some small interval, uniqueness of solutions for the optimal control problems. In that paper, one of the main questions under study is to know if the underlying epidemic structure has a significant impact on the obtained optimal control strategy. One of the models discussed in [10] was the SEIR, which is one of the most studied models in epidemiology. In the family of SEIRS models, it is assumed that the population is divided in four compartments: additionally to the infected class $I$, the susceptible class $S$ and the recovered class $R$, present in SIR models, an exposed class $E$ is also considered, in order to divide the infected population in the group of individuals that are infected and can infect others (the infective class) and the individuals that are infected but are not yet able to infect other individuals (the exposed or latent class). Such division of the population is particularity suitable to include several infectious diseases, e.g. measles and, assuming vertical transmission, rubella [17]. When there is no recovery, the model can be also used to describe diseases such as Chagas’ disease [29]. It is also appropriate to model hepatitis B and AIDS [17], and Ebola [23]. Although influenza can be modeled by a SEIRS model [5], due to the short latency period, it is sometimes better to use the simpler SIRS formulation [7].

In [10] the parameters of the considered models are assumed to be independent of time. This is not very realistic in many situations, in particular due to periodic seasonal fluctuations. A classical example of seasonal patterns of incidence, exhibited by some infectious diseases, is given by data on weekly measles notification in England and Wales during the period 1948–1968 [1]. Other examples occur in several childhood diseases such as mumps, chicken-pox, rubella and pertussis [19]. It is also worth to mention that some environmental and demographic effects can be non-periodic. For example, for some diseases like cholera and yellow fever, it is known that the size of the latency period may decrease with global warming [26]. In this work, we consider a controlled SEIRS model with vaccination and treatment as control variables but, unlike [10], we let the parameters of our model to be time dependent. One of our main objectives is to discuss the effect of seasonal behavior on the optimal strategy.

Another important aspect of our work is that we consider a model with a general incidence function given by some function $\varphi$ of the susceptible, the infective and the total population. This allows us to prove simultaneously results on existence and uniqueness of optimal solution for models with several different incidence functions that have already been considered in the context of SEIR/SEIRS models. In particular, our setting includes not only Michaelis–Menten incidence functions, considered for instance in [4, 11, 15, 21, 30, 33], but also incidence functions that are not bilinear, which are appropriate to include saturation effects as well as other non-linear behaviors [18, 35]. Additionally, and in contrast with [10], we assume that immunity can be partial and thus a fraction of the recovered individuals return to the susceptible class.

The paper is organized in the following way. In Section 2, we describe our non-autonomous SEIRS model, including treatment and vaccination as control variables, and we formulate the optimal control problem under investigation. Then, in Section 3, we discuss the question of existence of optimal solutions. The optimal controls are characterized in Section 4 with the help of Pontryagin’s minimum principle and, in Section 5, we present a result concerning uniqueness of the optimal control. We end with Section 6 of numerical simulations.

2. The non-autonomous SEIRS model and the optimal control problem

In practice, evolution on the number of susceptible, exposed, infective and recovered, depends on some factors that can be controlled. Two of the main factors are: treatment of infective and vaccination of susceptible. For this reason, we consider a non-autonomous SEIRS model including treatment, denoted by $\mathbbm{T}$, and vaccination, denoted by $\mathbbm{V}$, as control variables. Namely, we consider the optimal control problem

[TABLE]

where $\kappa_{1},\kappa_{2},\kappa_{3}$ and $S_{0},E_{0},I_{0},R_{0}$ are non-negative, the state variables are absolutely continuous functions, i.e., $(S(\cdot),E(\cdot),I(\cdot),R(\cdot))\in AC\left([0,t_{f}];\mathbb{R}^{4}\right)$, and the controls are Lebesgue integrable, i.e., $(\mathbbm{T}(\cdot),\mathbbm{V}(\cdot))\in L^{1}\left([0,t_{f}];[0,\tau_{\max}]\times[0,\nu_{\max}]\right)$. Moreover, the following conditions hold:

1. $C$1)

$\Lambda,\beta,\mu,\eta,\varepsilon,\gamma\in C^{1}([0,t_{f}];\mathbb{R})$ are $\omega$-periodic; 2. $C$2)

$\varphi\in C^{2}(\mathbb{R}^{3};\mathbb{R})$; 3. $C$3)

$\varphi(0,y,x)=\varphi(x,y,0)=0$ for all $x\in\mathbb{R}_{0}^{+}$ and $y\in\mathbb{R}^{+}$; 4. $C$4)

function $f:\mathcal{D}\to\mathbb{R}$ defined in $\mathcal{D}=\{(x,y,z)\in[0,t_{f}]^{3}:x+z\leq y\}$ by $f(x,y,z)=\varphi(x,y,z)/z$, if $z>0$, and $f(x,y,0)=\lim_{z\to 0}\varphi(x,y,z)/z$, is continuous and bounded.

The total population $N(t)=S(t)+E(t)+I(t)+R(t)$, $t\in[0,t_{f}]$, is not constant (see Remark 1). It can be also seen, from the equations that define our control system, that we assume treatment to be applied to infective individuals only, moving a fraction of them from the infected to the recovered compartment, and that vaccination is applied to susceptible individuals only, also moving a fraction of them to the recovered class. Moreover, note that in our problem (P), besides general non-autonomous parameters, we consider a general incidence function. As examples of such incidence functions, we mention the ones obtained by $\varphi(S,N,I)=SI/(1+\alpha I)$, considered for instance in [25, 34], $\varphi(S,N,I)=I^{p}S^{q}$, considered in [12, 14, 18], and $\varphi(S,N,I)=SI^{p}/(1+\alpha I^{q})$, considered in [13, 24]. Naturally, one needs to specify the incidence function, the other parameters and functions of the model, in order to carry out simulations and investigate particular situations. This is done in Section 6, where we compare an autonomous problem with a corresponding periodic model, as well as uncontrolled and corresponding controlled models. Before that, we prove results on existence of optimal solution (Theorem 3.2), necessary optimality conditions (Theorem 4.1), and uniqueness of solution (Theorem 5.1) for (P).

3. Existence of optimal solutions

Problem (P) is an optimal control problem in Lagrange form:

[TABLE]

In the above context, we say that a pair $(x,u)\in AC\left([t_{0},t_{1}];\mathbb{R}^{n}\right)\times L^{1}([t_{0},t_{1}];U)$ is feasible if it satisfies the Cauchy problem in (1). We denote the set of all feasible pairs by $\mathcal{F}$. Next, we recall a theorem about existence of solution for optimal control problems (1), contained in Theorem III.4.1 and Corollary III.4.1 in [9].

Theorem 3.1** (See [9]).**

For problem (1), suppose that $f$ and $\mathcal{L}$ are continuous and there exist positive constants $C_{1}$ and $C_{2}$ such that, for $t\in\mathbb{R}$, $x,x_{1},x_{2}\in\mathbb{R}^{n}$ and $u\in\mathbb{R}^{m}$, we have

1. a)

$\|f(t,x,u)\|\leq C_{1}(1+\|x\|+\|u\|)$; 2. b)

$\|f(t,x_{1},u)-f(t,x_{2},u)\|\leq C_{2}\|x_{1}-x_{2}\|(1+\|u\|)$; 3. c)

$\mathcal{F}$* is non-empty;* 4. d)

$U$* is closed;* 5. e)

there is a compact set $S$ such that $x(t_{1})\in S$ for any state variable $x$; 6. f)

$U$* is convex, $f(t,x,u)=\alpha(t,x)+\beta(t,x)u$, and $\mathcal{L}(t,x,\cdot)$ is convex on $U$;* 7. g)

$\mathcal{L}(t,x,u)\geq c_{1}|u|^{\beta}-c_{2}$, for some $c_{1}>0$ and $\beta>1$.

Then, there exist $(x^{*},u^{*})$ minimizing $J$ on $\mathcal{F}$.

Before applying Theorem 3.1 to obtain an existence result to our problem (P), we make a useful remark.

Remark 1**.**

Adding the equations in (P), we obtain that $N^{\prime}(t)=\Lambda(t)-\mu(t)N(t)$, which describes the behavior of the total population $N(t)$. Since $N_{0}=N(t_{0})=S_{0}+E_{0}+I_{0}+R_{0}$, we have

[TABLE]

where $K=N_{0}+\sup_{t\in[0,t_{f}]}\Lambda(t)/\inf_{t\in[0,t_{f}]}\mu(t)$. Thus,

[TABLE]

Consider the problem obtained by replacing function $\varphi$ in (P) by some bounded and twice continuously differentiable function $\psi$ such that $\psi(S,N,I)=\varphi(S,N,I)$ for $(S,N,I)\in[0,K]^{3}$ (and maintaining the initial condition, the cost functional and the set of admissible controls). By (2), this new problem has exactly the same solutions as problem (P).

Theorem 3.2** (Existence of solutions for the optimal control problem (P)).**

There exists an optimal control pair $\left(\mathbbm{T}^{*},\mathbbm{V}^{*}\right)$ and a corresponding quadruple $(S^{*},E^{*},I^{*},\break R^{*})$ of the initial value problem in (P) that minimizes the cost functional $\mathcal{J}$ in (P) over $L^{1}\left([0,t_{f}];[0,\tau_{\max}]\times[0,\nu_{\max}]\right)$.

Proof.

We show that the problem obtained from (P) by replacing the function $\varphi$ by some of the functions $\psi$ in Remark 1 satisfies the conditions in Theorem 3.1. By Remark 1, using conditions C1) and C3), we immediately obtain a) and b). Conditions c) and d) are immediate from the definition of $\mathcal{F}$ and since $U=[0,\tau_{\max}]\times[0,\nu_{\max}]$. By (2), we conclude that all state variables are in the compact set

[TABLE]

and condition e) follows. Since the state equations are linearly dependent on the controls, we obtain f). Finally, $\mathcal{L}$ is convex in the controls since it is quadratic. Moreover, $\mathcal{L}=k_{1}I+k_{2}\mathbbm{T}^{2}+k_{3}\mathbbm{V}^{2}\geq\min\{k_{2},k_{3}\}(\mathbbm{T}^{2}+\mathbbm{V}^{2})\geq\min\{k_{2},k_{3}\}\|(\mathbbm{T},\mathbbm{V})\|^{2}$ and we establish g) with $c_{1}=\min\{k_{2},k_{3}\}$. Thus, taking into account Remark 1, the result follows from Theorem 3.1. ∎

4. Characterization of the optimal controls

Now we address the question of how to identify the solutions predicted by Theorem 3.2. We do this with the help of the celebrated Pontryagin Maximum Principle [22]. This is possible because all control minimizers $\mathbbm{T}^{*}$ and $\mathbbm{V}^{*}$ of problem (P) are in fact essentially bounded. Indeed, Theorem 3.1 requires $U$ to be a closed set, not necessarily bounded, and it may happen, in general, that the optimal controls predicted by Theorem 3.1 are in $L^{1}$ but not in $L^{\infty}$ and do not satisfy the Pontryagin Maximum Principle [28]. However, in our case, $U$ is compact and we can conclude that the $L^{1}$ optimal controls $\left(\mathbbm{T}^{*},\mathbbm{V}^{*}\right)$ of problem (P), predicted by Theorem 3.2, are in fact in $L^{\infty}$, as required by the necessary optimality conditions [22]. Moreover, our optimal control problem (P) has only given initial conditions, with the state variables being free at the terminal time, that is, $S(t_{f})$, $E(t_{f})$, $I(t_{f})$, $R(t_{f})$ are free. This implies that abnormal minimizers [3] are not possible in our context and we can fix the cost multiplier associated with the Lagrangian $\mathcal{L}$ to be one: for our optimal control problem (P), the Hamiltonian is given by

[TABLE]

In what follows, we use the operator $\partial_{i}$ to denote the partial derivative with respect to the $i$th variable.

Theorem 4.1** (Necessary optimality conditions for the optimal control problem (P)).**

If $((S^{*},E^{*},I^{*},R^{*}),(\mathbbm{T}^{*},\mathbbm{V}^{*}))$ is a minimizer of problem (P), then there exist multipliers $(p_{1}(\cdot),p_{2}(\cdot),p_{3}(\cdot),p_{4}(\cdot))\in AC([0,t_{f}];\mathbb{R}^{4})$ such that

[TABLE]

for almost all $t\in[0,t_{f}]$, with transversality conditions

[TABLE]

Furthermore, the optimal control pair is given by

[TABLE]

and

[TABLE]

Proof.

Direct computations show that equations (3) are a consequence of the adjoint system of the Pontryagin Minimum Principle (PMP) [22]. Similarly, (4) are directly given by the transversality conditions of the PMP. It remains to characterize the controls using the minimality condition of the PMP [22]. The minimality condition on the set $\left\{t\in[0,t_{f}]:0<\mathbbm{V}^{*}(t)<\nu_{\max}\ \text{ and }\ 0<\mathbbm{T}^{*}(t)<\tau_{\max}\right\}$ is

[TABLE]

and thus, on this set,

[TABLE]

If $t\in\operatorname{int}\{t\in[t_{0},t_{1}]:\mathbbm{V}^{*}(t)=\nu_{\max}\}$, then the minimality condition is

[TABLE]

Analogously, if $t\in\operatorname{int}\{t\in[t_{0},t_{1}]:\mathbbm{T}^{*}(t)=\tau_{\max}\}$, then the minimality condition is

[TABLE]

If $t\in\operatorname{int}\{t\in[t_{0},t_{1}]:\mathbbm{V}^{*}(t)=0\}$, then the minimality condition is

[TABLE]

Analogously, if $t\in\operatorname{int}\{t\in[t_{0},t_{1}]:\mathbbm{T}^{*}(t)=0\}$, then the minimality condition is

[TABLE]

Therefore, we obtain (5) and (6). ∎

5. Uniqueness of the optimal control

In this section we show that the optimal solution of (P) is unique. The result is new even in the particular case where all the parameters of the model are time-invariant (autonomous case), extending the local uniqueness results in [8, 10], which are valid only in a small time interval, to uniqueness in a global sense, that is, uniqueness of solution of (P) along all the time interval $[0,t_{f}]$ where the optimal control problem is defined. The proof of Theorem 5.1 is nontrivial, lengthy and technical. It can be summarized as follows:

1. (1)

By contradiction, we assume that there are two distinct optimal pairs of state and co-state variables $\xi=(S,E,I,R,p_{1},p_{2},p_{3},p_{4})$ and $\xi^{*}=(S^{*},E^{*},I^{*},R^{*},$ $p_{1}^{*},p_{2}^{*},p_{3}^{*},p_{4}^{*})$, which correspond to two different optimal controls $u=\left(\mathbbm{T},\mathbbm{V}\right)$ and $u^{*}=\left(\mathbbm{T}^{*},\mathbbm{V}^{*}\right)$, in agreement with (5) and (6). 2. (2)

We show that both minimizing trajectories are in a positively invariant region $\Gamma$, which is given by the total population and is independent on the controls. 3. (3)

We make a change of variable and prove that we have a contradiction unless $\xi=\xi^{*}$ in a small time interval $[0,T]$; from (5)–(6), $u=u^{*}$ in $[0,T]$. 4. (4)

If $[0,T]$ contains the time range of the optimal control problem (P), then we are done. Otherwise, taking for initial conditions at time $T$ the values of the state trajectories at the right-end of the interval $[0,T]$, we obtain uniqueness for the interval $[T,2T]$, because the estimates that allow us to obtain $T$ are only related with the maximum value of the parameters and the bounds for the state and co-state variables on the invariant region $\Gamma$ of the new control problem, and are therefore the same as the ones already obtained for $[0,T]$ 5. (5)

Iterating the procedure, we obtain uniqueness throughout the interval $[0,t_{f}]$ after a finite number of steps.

Theorem 5.1** (Uniqueness of solution for the optimal control problem (P)).**

The solution of the optimal control problem (P) is unique.

Proof.

Let us assume that we have two optimal solutions corresponding to state trajectories and adjoint variables $(S,E,I,R)$ and $(p_{1},p_{2},p_{3},p_{4})$ and $(\bar{S},\bar{E},\bar{I},\bar{R})$ and $(\bar{p}_{1},\bar{p}_{2},\bar{p}_{3},\bar{p}_{4})$. Then, we show that the two are the same, at least in some small interval. To achieve this, we make the change of variables

[TABLE]

and

[TABLE]

Naturally, setting $n(t)=s(t)+e(t)+i(t)+r(t)$, we have

[TABLE]

By Proposition 1 in [20], we can assume that the trajectories lie in a compact and forward invariant set $\Gamma$, which can be chosen independently of the control functions $\mathbbm{T}$ and $\mathbbm{V}$. Using the differentiability assumption C2), we get

[TABLE]

where, because $\Gamma$ is compact, we have

[TABLE]

Considering the first equation in (P), we get, a.e. in $t\in[0,t_{f}]$, that

[TABLE]

and thus, for a.a. $t\in[0,t_{f}]$,

[TABLE]

Subtracting from the above equation the corresponding barred equation, we obtain

[TABLE]

Multiplying by $(s-\bar{s})$, integrating from [math] to $T$, and noting that $s(0)=\bar{s}(0)$,

[TABLE]

and, by (7), we obtain

[TABLE]

where $K_{1}$ depends on the bounds for $\bar{S}$ and $\mathbbm{V}$ and $C_{1}=\beta^{u}M_{1}^{u}+2\beta^{u}M_{2}^{u}+\beta^{u}M_{3}^{u}+\eta^{u}+2K_{1}$ (recall that $M_{i}^{u}$ is given by (8)). We now use some estimates for $(\mathbbm{V}-\bar{\mathbbm{V}})^{2}$ and $(\mathbbm{T}-\bar{\mathbbm{T}})^{2}$ that we prove later. Namely, we have

[TABLE]

where $C_{9}$ depends on $\alpha$ and on the maximum of $S$, $\bar{p}_{1}$ and $\bar{p}_{4}$ on $\Gamma$, and

[TABLE]

where $C_{10}$ depends on $\alpha$ and on the maximum of $I$, $\bar{p}_{3}$ and $\bar{p}_{4}$ on $\Gamma$ (see (24) and (25)). Define

[TABLE]

and

[TABLE]

Observe that $\Psi(t)\geq 0$ and $\Phi(t)\geq 0$ for all $t$. By (11), we obtain

[TABLE]

[TABLE]

By similar computations to the ones in (9) and (10), and using the inequality $xy\leq x^{2}+y^{2}$, we get

[TABLE]

where $C_{2}=\beta^{u}M_{1}^{u}+2\beta^{u}M_{2}^{u}+\beta^{u}M_{3}^{u}$. Recalling that

[TABLE]

where $C>0$ depends on the bounds for $\bar{x}$ and $y$, from the third equation in (P), by similar computations to the ones in (9) and (10), and using (12), we conclude that

[TABLE]

where $K_{2}$ depends on the maximum of $\bar{I}$ and $\mathbbm{T}$ on $\Gamma$ and $C_{3}=\varepsilon^{u}$. From the fourth equation in (P), by similar computations to the ones in (9) and (10), and using (11) and (12), we conclude that

[TABLE]

where $K_{3}$ and $K_{4}$ depend on the bounds for $\bar{I}$, $\bar{S}$, $\mathbbm{T}$ and $\mathbbm{V}$ and $C_{4}=\gamma^{u}$. From the first equation of (3), we get that

[TABLE]

for a.a. $t\in[0,t_{f}]$ and thus

[TABLE]

for a.a. $t\in[0,t_{f}]$. Using C2), we get

[TABLE]

[TABLE]

where, by C1) and since $\Gamma$ is compact, we have

[TABLE]

Subtracting from equation (17) the corresponding barred equation, we get

[TABLE]

Multiplying by $\varphi_{1}-\bar{\varphi}_{1}$ and integrating from [math] to $T$, we obtain

[TABLE]

Multiplying (19) by $-1$, we obtain that

[TABLE]

and thus, by (8) and (18), we conclude that

[TABLE]

Finally, we have

[TABLE]

[TABLE]

where $K_{5}$ and $K_{6}$ depend on the bounds for $\bar{p}_{1}$, $\bar{p}_{4}$ and $\mathbbm{V}$ and

[TABLE]

By (11) and (12), we obtain

[TABLE]

On the other hand, from the second equation of (3), one has $-\alpha e^{-\alpha t}\varphi_{2}+e^{-\alpha t}\dot{\varphi}_{2}=e^{-\alpha t}(\varphi_{1}-\varphi_{2})\beta{\partial_{2}\varphi}\left(e^{\alpha t}s,e^{\alpha t}n,e^{\alpha t}i\right)+e^{-\alpha t}\varphi_{2}\left(\mu+\varepsilon\right)-e^{-\alpha t}\varphi_{3}\varepsilon$. Straightforward computations show that

[TABLE]

where

[TABLE]

From the third equation of (3), we conclude that

[TABLE]

and, by similar computations as the ones done for the first equation of (3),

[TABLE]

where $K_{7}$ and $K_{8}$ depend on the bounds for $\bar{p}_{3}$, $\bar{p}_{4}$ and $\mathbbm{T}$ and

[TABLE]

From the fourth equation of (3),

[TABLE]

and, by similar computations as the ones done for the second equation of (3),

[TABLE]

where

[TABLE]

We now obtain the bounds for $\left(\mathbbm{V}-\bar{\mathbbm{V}}\right)^{2}$ and $\left(\mathbbm{T}-\bar{\mathbbm{T}}\right)^{2}$ announced in (11) and (12). We have

[TABLE]

and thus

[TABLE]

where

[TABLE]

Analogously, we obtain

[TABLE]

where

[TABLE]

Finally, we have all the bounds needed to prove our result. Adding equations (13), (14), (15), (16), (20), (21), (22) and (23), we obtain, for the sum of the left-hand sides,

[TABLE]

and thus

[TABLE]

which is equivalent to

[TABLE]

We now choose $\alpha$ so that $\alpha>\widetilde{C}+\widehat{C}$ and note that $\frac{\alpha-\widetilde{C}}{\widehat{C}}>1$. Subsequently, we choose $T$ such that

[TABLE]

Then,

[TABLE]

It follows that $\alpha-\widetilde{C}-\widehat{C}e^{3\alpha T}>0$, so that inequality (26) can hold if and only if we have $s(t)=\bar{s}(t)$, $e(t)=\bar{e}(t)$, $i(t)=\bar{i}(t)$, $r(t)=\bar{r}(t)$, $\varphi_{1}(t)=\bar{\varphi}_{1}(t)$, $\varphi_{2}(t)=\bar{\varphi}_{2}(t)$, $\varphi_{3}(t)=\bar{\varphi}_{3}(t)$ and $\varphi_{4}(t)=\bar{\varphi}_{4}(t)$ for all $t\in[0,T]$. This is equivalent to $S(t)=\bar{S}(t)$, $E(t)=\bar{E}(t)$, $I(t)=\bar{I}(t)$, $R(t)=\bar{R}(t)$, $p_{1}(t)=\bar{p}_{1}(t)$, $p_{2}(t)=\bar{p}_{2}(t)$, $p_{3}(t)=\bar{p}_{3}(t)$ and $p_{4}(t)=\bar{p}_{4}(t)$. With this, the uniqueness of the optimal control is established on the interval $[0,T]$. If $T\geq t_{f}$, then we have uniqueness on the whole interval. Otherwise, we can obtain uniqueness on $[T,2T]$ for the optimal control problem whose initial conditions on time $T$ coincide with the values of $S$, $E$, $I$ and $R$ on the end-time of the interval $[0,T]$ (note that, by the forward invariance of the set $\Gamma$, and since the constants $\widetilde{C}$ and $\widehat{C}$ in (27) depend only on the values of the several state and co-state variables on $\Gamma$, we still have the same $T$). Proceeding in the same way, we conclude, after a finite number of steps, that we have uniqueness on the interval $[0,t_{f}]$. The proof is complete. ∎

6. Numerical simulations

In what follows, the incidence into the exposed class of susceptible individuals and the birth function $\Lambda\left(t\right)$ are chosen, for illustrative purposes, to be

[TABLE]

and

[TABLE]

with $\mathrm{per}\in[0,1[$. The remaining parameter functions, $\mu\left(t\right)$, $\eta\left(t\right)$, $\varepsilon\left(t\right)$ and $\gamma\left(t\right)$, are assumed constant. Their values, as well as several other used in this section, were taken from [31] and [32] and are presented in Table 1.

Note that we are in the autonomous case when $\mathrm{per}=0$ and in a periodic situation for $0<\mathrm{per}<1$.

The solutions to the optimal control problems (P) here considered exist (Theorem 3.2), are unique (Theorem 5.1) and can be found using the optimality conditions given in Theorem 4.1. Precisely, our method to solve the problem consists to use the state equations (the four ordinary differential equations of problem (P)), the initial conditions, the adjoint equations (3) and the transversality conditions (4) with the controls given by (5) and  (6), which are substituted into the state and adjoint equations. The state equations are solved with the initial conditions of Table 1, while the adjoint system is solved backward in time, with the change of variable

[TABLE]

Then, the two systems of equations constitute an initial value problem, which is solved numerically with an explicit 4th and 5th order Runge–Kutta method through the ode45 solver of MATLAB. The procedure is briefly described in Algorithm 1.

In each plot of Figures 1 and 2, we present two sets of trajectories (distinguished by using dashed and continuous lines), in order to easily compare the uncontrolled and optimally controlled situations, the former mentioned by the suffix “-u” in the figures’ legend.

The behavior of our SEIRS model with both $\mathrm{per}=0$ (autonomous case, Figure 1) and $\mathrm{per}=0.8$ (periodic case, Figure 2), show the effectiveness of optimal control theory. Indeed, in Figures 1 and 2 we observe that, if we apply treatment and vaccination as given by Theorem 4.1 (optimally controlled case), then the number of exposed and infected individuals is significantly lower, as well as the number of susceptible individuals, while the number of recovered is significantly higher. It can be also seen that the susceptible and recovered classes have a very different behavior in the controlled and uncontrolled situations.

In Figures 3 and 4, we have the same trajectories as in Figures 1 and 2. They illustrate the effect of the periodicity of $\Lambda(t)$ (29) and $\beta(t)$ (28) in the classes $S$, $E$, $I$ and $R$ of individuals. The effect is perceptible in susceptible and exposed classes, since the periodic functions are present in these classes. From our results, we claim that the periodicity effect is “softened” in the transition between classes.

In Figure 5, we show the optimal controls: treatment $\mathbbm{T}^{*}(t)$ of infective individuals (left side) and vaccination $\mathbbm{V}^{*}(t)$ of susceptible (right side). According to the minimality conditions (5) and (6), both controls go to zero when $t\to t_{f}=25$. The periodicity effect is perceptible in $\mathbbm{V}^{*}(t)$, consequence of the fact that vaccination takes place in the susceptible class. Treatment occurs in the infective class and, as we have seen, in this class the periodicity is not so perceptible. As a consequence, periodicity is only slightly perceptible in the treatment control variable $\mathbbm{T}^{*}(t)$.

In Figures 6 to 9, we present the behavior of infected individuals $I^{*}(t)$, treatment $\mathbbm{T}^{*}(t)$ and vaccination $\mathbbm{V}^{*}(t)$ optimal controls, when we vary the parameters $\mu$, $\gamma$, $\varepsilon$ and $\eta$, respectively, maintaining, in each case, the initial values and the remaining parameters as in Table 1. In all Figures 6–9, we varied the respective parameter ($\mu$, $\gamma$, $\varepsilon$ and $\eta$) from [math] to $0.1$ in steps of length $0.01$.

Referring to Figure 6, where the variation of $\mu$ is analyzed, we can see that the effect of periodicity is more perceptible in vaccination than in treatment. In the infected class, for low values of $\mu$ (low mortality), we observe that the number of infected increases. This is justified by the difference between birth and death.

Concerning Figure 7, where we vary the rate of recovery $\gamma$, the effect of periodicity is analogous to the previous situation: the bigger the value of $\gamma$, more infected individuals recover and thus the faster the infected class decreases.

In Figure 8, one can see the effect of the variation of the infectivity rate $\varepsilon$. The effect of periodicity is also in this case analogous to the previous situations. Indeed, when we have a higher value of $\varepsilon$, we have a faster transition of exposed individuals into the infected class and this is the reason why we observe in Figure 8 that an increase on the value of $\varepsilon$ leads to an increase in the infected class.

Finally, in Figure 9, the variation of the loss of immunity rate $\eta$ is highlighted. We conclude that the periodicity effect is similar to the previous considered scenarios. However, the variation of $\eta$ is the one that less influences the behavior of the three variables $I^{*}(t)$, $\mathbbm{T}^{*}(t)$ and $\mathbbm{V}^{*}(t)$.

It is worth mentioning that, in the simulations done and in the range of parameters considered, the variation of the periodic parameter $\mathrm{per}$ has a very small effect on the obtained cost functional $\mathcal{J}$. Namely, we saw numerically that

[TABLE]

for $v_{1},v_{2}\in\{0,0.8\}$.

Acknowledgments

The authors are grateful to an anonymous Referee for valuable comments and suggestions.