# The Hopf Bifurcation Analysis and Optimal Control of a Delayed SIR Epidemic Model

Authors: Abdelhadi Abta; Hassan Laarabi; Hamad Talibi Alaoui
Journal: International Journal of Analysis
 (2014)

---

## Abstract

We propose a delayed SIR model with  saturated incidence rate.
The  delay is incorporated into the model in order to model the latent period. The basic reproductive numberR
0 is obtained. Furthermore, using time delay as a bifurcation parameter, it is proven that there exists a critical value of delay for the stability of diseases prevalence. When the delay exceeds the critical value, the system loses its stability and a Hopf bifurcation occurs. The model is extended to assess
the impact of some control measures, by reformulating the model as an optimal control problem with vaccination and treatment. The existence of the optimal control  is also proved.  Finally, some numerical simulations are performed to verify the theoretical
analysis.


---

## Body

## 1. Introduction

Mathematical modelling is of considerable importance in the study of epidemiology because it may provide understanding of the underlying mechanisms which influence the spread of disease and may suggest control strategies. The first known mathematical model of epidemiology is formulated and solved by Daniel Bernoulli in 1760. The foundations of the modern mathematical epidemiology based on the compartment models were laid in the early 20th century [1]. Since the middle of the 20th century, mathematical epidemiology has grown exponentially. In particular, the SIR epidemic model is known as one of the most basic epidemic models, in which total host population is divided into three classes called susceptible S, infective I, and removed R. The basic and important research subjects for these systems are the existence of the threshold value which distinguishes whether the infectious disease will die out, the local stability of the disease-free equilibrium and the endemic equilibrium, the Hopf bifurcation, the existence of periodic solutions, optimal control, and so forth. Many models in the literature represent the dynamics of disease by systems of ordinary differential equations without time delay. In order to reflect the real dynamical behaviors of models that depend on the past history of systems, it is reasonable to incorporate time delays into the systems [2]. In fact, inclusion of delays in epidemic models makes them more realistic by allowing the description of the effects of disease latency or immunity [3, 4].In this paper, we propose the delayed SIR epidemic model governed by the following equations [5]:
(1)
d
S
d
t
=
r
(
1
-
S
(
t
)
K
)
S
(
t
)
-
β
S
(
t
)
I
(
t
)
1
+
α
2
I
(
t
)
,
d
I
d
t
=
β
S
(
t
-
τ
)
I
(
t
-
τ
)
1
+
α
2
I
(
t
-
τ
)
-
(
μ
+
α
+
γ
)
I
(
t
)
,
d
R
d
t
=
γ
I
(
t
)
-
μ
R
(
t
)
,
where S is the number of susceptible individuals, I is the number of infectious individuals, R is the number of recovered individuals, r is the specific growth rate, K is the environment capacity, β is the transmission rate, μ is the natural death of the population, α is the death rate due to disease, α
2 is the parameters that measure the inhibitory effect, γ is the recovery rate of the infectious individuals, and τ is the incubation period.This model takes into account a number of key biological assumptions.(i)
It is more reasonable to assume that the population of a given region obeys logistic growth due to crowding and limited sources. Epidemic models with logistic or generalized logistic demographic structure have been extensively studied (see, e.g., [6, 7]).
(ii)
Since nonlinearity in the incidence rates has been observed in disease transmission dynamics, it has been suggested that the standard bilinear incidence rate will be modified into a nonlinear incidence rate by many authors (see, e.g., [8, 9]). In this work we use a nonlinear incidence rate of the form β
S
I
/
(
1
+
α
2
I
).The aim is to gain some insights into the best intervention for minimizing the transmission of disease within the population and to explore the impacts of various intervention scenarios, namely, vaccination and treatment. We analyse the stability and the Hopf bifurcation of the model; then, we incorporate into the model appropriate cost functions in order to study and determine the possible impacts of these strategies on controlling the disease. We give the necessary conditions for optimal control of the disease using Pontryagin’s maximum principle, in order to determine optimal strategies for controlling the spread of the disease.The organization of the paper is as follows: we examine the existence and stability behaviour of the equilibrium solutions of this model and we are particularly interested in whether the unique endemic equilibrium can be destabilized by introducing a time delay to represent the effect of fading of vaccination and whether stable endemic cycles can arise by the Hopf bifurcation from the initially stable endemic equilibrium in Section2. In Section 3, we formulate an optimal control problem and we use Pontryagin’s maximum principle with delay given in [10] to characterize it. Our conclusions are discussed in Section 4.


## 2. Stability Analysis and the Hopf Bifurcation Occurrence

In this section, local stability of each feasible equilibrium of the model (1) is established; conditions are found under which the Hopf bifurcation occurs and periodic solutions emerge as the delay crosses some critical value. Some numerical simulations to illustrate the theoretical results are given.Before going into any detail, we simplify the model since the first two equations of (1) are independent of the third one; it suffices to consider the first two equations. Thus, we restrict our attention to the following reduced model:
(2)
d
S
d
t
=
r
(
1
-
S
(
t
)
K
)
S
(
t
)
-
β
S
(
t
)
I
(
t
)
1
+
α
2
I
(
t
)
,
d
I
d
t
=
β
S
(
t
-
τ
)
I
(
t
-
τ
)
1
+
α
2
I
(
t
-
τ
)
-
(
μ
+
α
+
γ
)
I
(
t
)
.
Now, we will find the equilibria of system (2) and investigate their dynamical features. System (2) always has two disease-free equilibriaP
0
=
(
0,0
)andP
1
=
(
K
,
0
). Further, if
(3)
R
0
=
β
K
μ
+
α
+
γ
>
1
,
system (2) admits a unique endemic equilibrium P
*
=
(
S
*
,
I
*
), given by
(4)
S
*
=
(
μ
1
+
γ
)
(
1
+
α
2
I
*
)
β
,
I
*
=
[
α
2
r
β
K
-
2
α
2
r
(
μ
1
+
γ
)
-
β
2
K
]
+
Δ
2
r
α
2
2
(
μ
1
+
γ
)
,
with
(5)
Δ
=
[
α
2
r
β
K
-
2
α
2
r
(
μ
1
+
γ
)
-
β
2
K
]
2
+
4
α
2
2
r
(
μ
1
+
γ
)
(
r
β
K
-
r
(
μ
1
+
γ
)
)
.
Now, let us start to discuss the local behavior of the disease-free equilibrium P
0. We linearize system (2) around P
0; we have
(6)
d
S
d
t
=
r
S
(
t
)
,
d
I
d
t
=
-
(
μ
1
+
γ
)
I
(
t
)
.
The characteristic equation associated with system (6) is
(7)
Δ
(
λ
,
τ
)
=
(
λ
-
r
)
(
λ
+
μ
1
+
γ
)
=
0
.Proposition 1.
The disease-free equilibrium pointP
0 is unstable.Proof.
It is clear that (7) has two roots λ
1
=
r
>
0 and λ
2
=
-
(
μ
1
+
γ
)
<
0. Whence P
0 is unstable.We turn to study the local behavior of the disease-free equilibriumP
1.Letx
=
S
-
K and y
=
I. The linearized system of (2) around P
1 takes the form
(8)
d
x
d
t
=
-
r
x
(
t
)
-
β
K
y
(
t
)
,
d
y
d
t
=
β
K
y
(
t
-
τ
)
-
(
μ
1
+
γ
)
y
(
t
)
.
The characteristic equation associated with system (8) is
(9)
Δ
(
λ
,
τ
)
=
(
λ
+
r
)
[
λ
+
(
μ
1
+
γ
)
-
β
K
e
-
λ
τ
]
=
0
.Proposition 2.
IfR
0
<
1, then the disease-free equilibrium P
1 is locally asymptotically stable. And if R
0
>
1, then the equilibrium point P
1 is unstable.Proof.
Indeed, ifτ
=
0, then (9) becomes
(10)
Δ
(
λ
,
0
)
=
(
λ
+
r
)
[
λ
+
(
μ
1
+
γ
)
(
1
-
R
0
)
]
=
0
.
It is easy to see that (10) has two roots λ
1
=
-
r
<
0 and λ
2
=
(
μ
1
+
γ
)
(
R
0
-
1
). So, if R
0
<
1, the equilibrium P
1 is asymptotically stable and if R
0
>
1, the equilibrium P
1 is unstable.
By Rouché’s Theorem [8, p. 248], it follows that if instability occurs for a particular value of the delay τ, a characteristic root of (9) must intersect the imaginary axis. Suppose that (9) has a purely imaginary root i
ω, with ω
>
0. Then, by separating real and imaginary parts in (9), we have
(11)
β
K
cos
⁡
⁡
(
ω
τ
)
=
μ
+
γ
,
β
K
sin
⁡
(
ω
τ
)
=
-
ω
.
Hence,
(12)
ω
2
=
[
β
K
+
μ
1
+
γ
]
(
μ
1
+
γ
)
(
R
0
-
1
)
.
So, if R
0
<
1, (9) has no purely imaginary root, and as P
1 is asymptotically stable for τ
=
0 it remains asymptotically stable for all τ
≥
0.
IfR
0
<
1, then the disease-free equilibrium P
1 is unstable for τ
=
0. By Kuang’s Theorem [4, p. 77], it follows that P
1 is unstable for all τ
≥
0. This concludes the proof.Next, we establish the local stability of the endemic equilibriumP
*and we determine the conditions under which the Hopf bifurcation occurs.Letx
=
S
-
S
*andy
=
I
-
I
*. The linearized system of (2) aroundP
*
=
(
S
*
,
I
*
)takes the form
(13)
d
x
d
t
=
-
r
S
*
K
x
(
t
)
-
β
S
*
(
1
+
α
2
I
*
)
2
y
(
t
)
,
d
y
d
t
=
β
I
*
1
+
α
2
I
*
x
(
t
-
τ
)
+
β
S
*
(
1
+
α
2
I
*
)
2
y
(
t
-
τ
)
-
(
μ
1
+
γ
)
y
(
t
)
,
and the characteristic equation is given by
(14)
Δ
(
λ
,
τ
)
=
λ
2
+
a
λ
+
b
λ
exp
⁡
⁡
(
-
λ
τ
)
+
d
exp
⁡
⁡
(
-
λ
τ
)
+
c
=
0
,
with
(15)
a
=
(
μ
1
+
γ
)
+
r
S
*
K
,
b
=
-
β
S
*
(
1
+
α
2
I
*
)
2
,
c
=
r
(
μ
1
+
γ
)
S
*
K
,
d
=
β
2
I
*
S
*
(
1
+
α
2
I
*
)
3
-
r
β
S
*
2
K
(
1
+
α
2
I
*
)
2
.
We begin by considering the case without delayτ
=
0. We have the following proposition.Proposition 3.
IfR
0
>
1, then the equilibrium point P
* is asymptotically stable for τ
=
0.Proof.
Indeed, ifτ
=
0, then (14) becomes
(16)
Δ
(
λ
,
0
)
=
λ
2
+
(
a
+
b
)
λ
+
(
c
+
d
)
=
0
,
with
(17)
a
+
b
=
r
S
*
K
+
(
μ
1
+
γ
)
-
β
S
*
(
1
+
α
2
I
*
)
2
,
c
+
d
=
β
2
I
*
S
*
(
1
+
α
2
I
*
)
3
+
r
S
*
K
[
(
μ
1
+
γ
)
-
β
S
*
(
1
+
α
2
I
*
)
2
]
.
Since (
μ
1
+
γ
)
-
β
S
*
/
(
1
+
α
2
I
*
)
2
=
(
μ
1
+
γ
)
[
α
2
I
*
/
(
1
+
α
2
I
*
)
]
>
0, it follows that a
+
b
>
0 and c
+
d
>
0. So, according to the Hurwitz criterion, all roots of (16) have negative real parts. Hence, P
* is asymptotically stable.In the following, we treat the case of positive delayτ
>
0. Let
(18)
(
H
1
)


:
S
*
2
+
2
(
μ
1
+
γ
)
β
S
*
-
K
(
μ
1
+
γ
)
β
>
0
,
(
H
2
)


:
S
*
2
+
2
(
μ
1
+
γ
)
β
S
*
-
K
(
μ
1
+
γ
)
β
<
0
.
We have the following two results.Theorem 4.
IfR
0
>
1 and the assumption (
H
1
) is satisfied, then the equilibrium point P
* is asymptotically stable for all τ
≥
0.Proof.
Indeed, it follows that if instability occurs for a particular value of the delayτ, a characteristic root of (14) must intersect the imaginary axis. If (14) has a purely imaginary root i
ω, with ω
>
0, then, by separating real and imaginary parts in (14), we have
(19)
d
cos
⁡
⁡
(
ω
τ
)
+
b
ω
sin
⁡
(
ω
τ
)
=
ω
2
-
c
,
b
ω
cos
⁡
⁡
(
ω
τ
)
-
d
sin
⁡
(
ω
τ
)
=
-
a
ω
.
Hence,
(20)
ω
4
+
(
a
2
-
b
2
-
2
c
)
ω
2
+
c
2
-
d
2
=
0
.
From the expressions ofa,b,c, andd, we have
(21)
c
+
d
=
r
S
*
K
[
(
μ
1
+
γ
)
-
β
S
*
(
1
+
α
2
I
*
)
2
]
+
β
2
I
*
S
*
(
1
+
α
2
I
*
)
3
,
c
-
d
=
r
β
K
(
1
+
α
2
I
*
)
[
S
*
2
+
2
(
μ
1
+
γ
)
β
S
*
-
K
(
μ
1
+
γ
)
β
]
,
a
2
-
b
2
-
2
c
=
[
(
μ
1
+
γ
)
-
β
S
*
(
1
+
α
2
I
*
)
2
]
×
[
(
μ
1
+
γ
)
+
β
S
*
(
1
+
α
2
I
*
)
2
]
+
(
r
S
*
K
)
2
.
Therefore, if R
0
>
1 and if the assumption (
H
1
) is satisfied, thenc
2
-
d
2
>
0anda
2
-
b
2
-
2
c
>
0. Hence, according to the criterion of Hurwitz, (20) has no positive root. Thus, the characteristic equation (14) admits no purely imaginary root and, asP
*is asymptotically stable forτ
=
0, it is asymptotically stable for allτ
≥
0.Theorem 5.
IfR
0
>
1 and if the assumption (
H
2
) is satisfied, then there exists τ
0
>
0 such that, for all τ
∈
[
0
,
τ
0
), the equilibrium point P
* is asymptotically stable and, for all τ
>
τ
0, the equilibrium point P
* is unstable, and when τ
=
τ
0, a Hopf bifurcation of periodic solutions of system (2) occurs at P
*. With
(22)
ω
0
2
=
1
2
(
b
2
+
2
c
-
a
2
)
+
1
2
[
(
b
2
+
2
c
-
a
2
)
2
-
4
(
c
2
-
d
2
)
]
1
/
2
,
τ
0
=
1
ω
0
arccos
⁡
(
d
-
a
b
)
ω
0
2
-
c
d
b
2
ω
0
2
+
d
2
,
wherea,b,c, anddare defined in (14).Proof.
The demonstration of stability ofP
* under the assumption (
H
2
) is similar to the previous demonstration. It suffices to remark that if R
0
>
1 and if the assumption (
H
2
) is satisfied, then we deduce that
(23)
c
+
d
=
r
S
*
K
[
(
μ
1
+
γ
)
-
β
S
*
(
1
+
α
2
I
*
)
2
]
+
β
2
I
*
S
*
(
1
+
α
2
I
*
)
3
>
0
,
c
-
d
=
r
β
K
(
1
+
α
2
I
*
)
[
S
*
2
+
2
(
μ
1
+
γ
)
β
S
*
-
K
(
μ
1
+
γ
)
β
]
<
0
.
Thus, the characteristic equation (14) has only one purely imaginary root with positive imaginary part for the values of τ given by
(24)
τ
n
=
1
ω
0
arccos
⁡
(
d
-
a
b
)
ω
0
2
-
c
d
b
2
ω
0
2
+
d
2
+
2
n
π
ω
0
.
Since P
* is asymptotically stable for τ
=
0, by Kuang’s Theorem [4, p. 77], the equilibrium P
* is asymptotically stable for τ
∈
[
0
,
τ
0
) and is unstable for all τ
>
τ
0.
Lastly, to complete the proof, it remains to verify the transversality condition(25)
[
d
Re
(
λ
)
d
τ
]
τ
=
τ
0
≠
0
.
We have
(26)
Sign
[
d
Re
(
λ
)
d
τ
]
τ
=
τ
0
=
Sign
[
w
4
-
c
2
+
d
2
]
.
Consequently,
(27)
[
d
Re
(
λ
)
d
τ
]
τ
=
τ
0
>
0
.
So the conditions of the Hopf bifurcation theorem are satisfied, which completes the proof.We conclude this section by giving some numerical simulations to illustrate the theoretical results. We show numerically that the system (2) has a family of periodic solutions.Proposition 6.
Letr
=
0.1, β
=
0.1, σ
=
0.96, γ
=
0.5, μ
1
=
0.5, and K
=
70. Then, the system (2) has an equilibrium point P
*
=
(
10.45,0.875
) which is asymptotically stable if 0
≤
τ
<
1.509 and unstable if τ
>
1.509 and when τ
=
1.509, the system has a family of periodic solutions (see Figure 1).Figure 1
The equilibrium pointP
*
=
(
10.45
,
0.875
)is asymptotically stable if τ
=
0.2 and unstable if τ
=
3 and when τ
=
1.509, the system has a family of periodic solutions.


## 3. The Optimal Control Problem

Generally, the eradication of the disease may be too costly when constant controls are considered as it requires treatment/vaccination at higher levels all the time. For eradication to be achievable in a finite time, we need to consider time-dependent controls. We use the optimal control strategies in the form of vaccination and treatment to decrease the number of both susceptible and infectious individuals and increase the total number of recovered individuals with minimum investment in disease control. This problem is formulated as an optimal control problem by introducing two controlsu
1andu
2, which represents the percentage of susceptible and infected individuals being vaccinated and treated, respectively, per unit of time. Hence, (1) becomes
(28)
d
S
d
t
=
r
(
1
-
S
(
t
)
K
)
S
(
t
)
-
β
S
(
t
)
I
(
t
)
1
+
α
2
I
(
t
)
-
u
1
(
t
)
S
(
t
)
,
d
I
d
t
=
β
S
(
t
-
τ
)
I
(
t
-
τ
)
1
+
α
2
I
(
t
-
τ
)
-
(
α
+
μ
+
γ
)
I
(
t
)
-
u
2
(
t
)
I
(
t
)
,
d
R
d
t
=
γ
I
(
t
)
-
μ
R
(
t
)
+
u
1
(
t
)
S
(
t
)
+
u
2
(
t
)
I
(
t
)
,
(29)
S
(
0
)
=
S
0
,
I
(
0
)
=
I
0
,
R
(
0
)
=
R
0
.
It is easy to show that there exists a unique solution (
S
(
t
)
,
I
(
t
)
,
R
(
t
)
) of system (28) with initial data (
S
0
,
I
0
,
R
0
)
∈
(
C
+
)
3.In addition, for biological reasons, we assume that the initial data for system (28) satisfy
(30)
S
0
(
t
)
≥
0
,
I
0
(
t
)
≥
0
,
R
0
(
t
)
≥
0
,
t
∈
[
-
τ
,
0
]
.
The problem is to minimize the objective (cost) functional given by
(31)
J
(
u
1
,
u
2
)
=
∫
0
t
f
‍
[
A
1
S
(
t
)
+
A
2
I
(
t
)
+
1
2
B
1
u
1
2
(
t
)
+
1
2
B
1
u
2
2
(
t
)
]
d
t
.
Subject to the differential equations (28), where the first two terms in the functional objective represent benefit of S
(
t
) and I
(
t
) populations that we wish to reduce, the parameters A
1 and A
2 are positive constants to keep a balance in the size of S
(
t
) and I
(
t
), respectively. We use in the second term in the functional objective the quadratic term (
1
/
2
)
B
i
u
i
2
,
i
=
1
,
2, where B
i is a positive weight parameter which is associated with the control u
i
(
t
) and the square of the control variable reflects the severity of the side effects of the vaccination and treatment.Our target is to minimize the objective functional defined in (31) by decreasing the number of infected and susceptible individuals and increasing the number of recovered individuals, by using possible minimal control variables (
u
1
(
t
)
,
u
2
(
t
)
). In other words, the control variables (
u
1
(
t
)
,
u
2
(
t
)
)
∈
U
ad represent the percentage of susceptible and infected individuals being vaccinated and treated, respectively, per unit of time and U
ad is the control set defined by
(32)
U
ad
=
{
u
i
max
⁡
<
∞
,
t
∈
[
0
,
t
f
]
u
=
(
u
1
,
u
2
)
∣
u
i
(
t
)
measurable
,
kkk
0
≤
u
i
(
t
)
≤
u
i
max
⁡
<
∞
,
t
∈
[
0
,
t
f
]
,
i
=
1,2
}
.


### 3.1. Existence of an Optimal Control

The existence of the optimal control pair can be obtained using a result by Fleming and Rishel in [11] and by Lukes in [12].Theorem 7.
There exists a control functionu
1
*
(
t
),u
2
*
(
t
) so that
(33)
J
(
u
1
*
(
t
)
,
u
2
*
(
t
)
)
=
min
⁡
(
u
1
,
u
2
)
∈
U
a
d
J
(
u
1
(
t
)
,
u
2
(
t
)
)
.Proof.
To prove the existence of an optimal control pair, it is easy to verify the following.(1)
The set of controls and corresponding state variables is nonempty.
(2)
The admissible setU
ad is convex and closed.
(3)
The right-hand side of the state system (28) is bounded by a linear function in the state and control variables.
(4)
The integrand of the objective functional is convex onU
ad.
(5)
There exists constantsω
1
,
ω
2
>
0, andρ
>
1such that the integrandL
(
S
,
I
,
u
1
,
u
2
)of the objective functional satisfiesL
(
S
,
I
,
u
1
,
u
2
)
≥
ω
2
+
ω
1
(
|
u
1
|
2
+
|
u
2
|
2
)
ρ
/
2.
The result follows directly from [11].


### 3.2. Characterization of the Optimal Control

Before characterizing the optimal control pair, we first define the Lagrangian for the optimal control problem (28)–(31) by
(34)
L
(
S
,
I
,
u
1
,
u
2
)
=
A
1
S
(
t
)
+
A
2
I
(
t
)
+
1
2
B
1
u
1
2
(
t
)
+
1
2
B
2
u
2
2
(
t
)
and the HamiltonianH for the control problem by
(35)
H
(
S
,
I
,
R
,
u
1
,
u
2
,
λ
i
,
t
)
=
L
(
S
,
I
,
u
1
,
u
2
)
+
∑
i
=
1
i
=
3
‍
λ
i
f
i
,
whereλ
i, i
=
1
,
2
,
3, are the adjoint functions to be determined suitably. Next, by applying Pontryagin’s maximum principle with delay given in [10] to the Hamiltonian H, we obtain the following theorem.Theorem 8.
Given optimal controlsu
1
*
(
t
) and u
2
*
(
t
) and solutions S
*
(
t
), I
*
(
t
), and R
*
(
t
) of the corresponding state systems (31) and (28), there exist adjoint variables λ
1, λ
2, and λ
3 that satisfy
(36)
d
λ
1
(
t
)
d
t
=
-
A
1
+
λ
1
(
t
)
(
r
(
1
-
2
S
*
K
)
-
Λ
1
-
u
1
*
)
-
λ
3
(
t
)
u
1
*
(
t
)
-
χ
[
0
,
t
f
-
τ
]
λ
2
(
t
+
τ
)
Λ
1
,
d
λ
2
(
t
)
d
t
=
-
A
2
+
λ
1
(
t
)
Λ
2
+
λ
2
(
t
)
(
μ
+
α
+
γ
+
u
2
*
)
-
λ
3
(
t
)
(
γ
+
u
2
*
)
-
χ
[
0
,
t
f
-
τ
]
λ
2
(
t
+
τ
)
Λ
2
,
d
λ
3
(
t
)
d
t
=
λ
3
(
t
)
μ
,
whereΛ
1
=
(
β
I
*
(
1
+
α
2
I
*
)
/
(
1
+
α
2
I
*
)
2
)andΛ
2
=
(
β
S
*
/
(
1
+
α
2
I
*
)
2
)
(37)
w
i
t
h
t
r
a
n
s
v
e
r
s
a
l
i
t
y
c
o
n
d
i
t
i
o
n
s
λ
i
(
t
f
)
=
0
,
i
=
1
,
2
,
3
.
Furthermore, the optimal control pairu
*
(
t
)is given by
(38)
u
1
*
(
t
)
=
max
⁡
⁡
(
min
⁡
⁡
(
(
λ
1
(
t
)
-
λ
3
(
t
)
)
S
*
(
t
)
B
1
,
u
1
max
⁡
)
,
0
)
,
u
2
*
(
t
)
=
max
⁡
⁡
(
min
⁡
⁡
(
(
λ
2
(
t
)
-
λ
3
(
t
)
)
I
*
(
t
)
B
2
,
u
2
max
⁡
)
,
0
)
.Proof.
Using Pontryagin’s maximum principle with delay in state, we obtain the adjoint equations and transversality conditions such that(39)
d
λ
1
(
t
)
d
t
=
-
∂
H
∂
S
-
χ
[
0
,
t
f
-
τ
]
∂
H
∂
S
τ
(
t
+
τ
)
,
λ
1
(
t
f
)
=
0
,
d
λ
2
(
t
)
d
t
=
-
∂
H
∂
I
-
χ
[
0
,
t
f
-
τ
]
∂
H
∂
I
τ
(
t
+
τ
)
,
λ
2
(
t
f
)
=
0
,
d
λ
3
(
t
)
d
t
=
-
∂
H
∂
R
-
χ
[
0
,
t
f
-
τ
]
∂
H
∂
R
τ
(
t
+
τ
)
,
λ
3
(
t
f
)
=
0
,
and by using the optimality conditions we find
(40)
∂
H
∂
u
1
=
B
1
u
*
(
t
)
-
λ
1
(
t
)
S
*
+
λ
3
(
t
)
S
*
=
0
,
at
u
1
=
u
1
*
(
t
)
,
∂
H
∂
u
2
=
B
2
u
*
(
t
)
-
λ
2
(
t
)
I
*
+
λ
3
(
t
)
S
*
=
0
,
at
u
2
=
u
2
*
(
t
)
,
which gives
(41)
u
1
*
(
t
)
=
(
λ
1
(
t
)
-
λ
3
(
t
)
)
S
*
(
t
)
B
1
,
u
2
*
(
t
)
=
(
λ
2
(
t
)
-
λ
3
(
t
)
)
I
*
(
t
)
B
2
.
Using the property of the control space, we obtain
(3.2)
u
1
*
(
t
)
=
0
,
if
(
λ
1
(
t
)
-
λ
3
(
t
)
)
S
*
(
t
)
B
1
≤
0
,
k
k
k
k
u
1
*
(
t
)
=
(
λ
1
(
t
)
-
λ
3
(
t
)
)
S
*
(
t
)
B
1
,
if
0
<
(
λ
1
(
t
)
-
λ
3
(
t
)
)
S
*
(
t
)
B
1
<
u
1
max
⁡
,
u
1
*
(
t
)
=
u
1
max
⁡
,
if
(
λ
1
(
t
)
-
λ
3
(
t
)
)
S
*
(
t
)
B
1
≥
u
1
max
⁡
,
u
2
*
(
t
)
=
0
,
if
(
λ
2
(
t
)
-
λ
3
(
t
)
)
I
*
(
t
)
B
2
≤
0
,
k
k
k
k
u
2
*
(
t
)
=
(
λ
1
(
t
)
-
λ
3
(
t
)
)
I
*
(
t
)
B
2
,
if
0
<
(
λ
1
(
t
)
-
λ
3
(
t
)
)
I
*
(
t
)
B
2
<
u
2
max
⁡
,
u
2
*
(
t
)
=
u
2
max
⁡
,
if
(
λ
1
(
t
)
-
λ
3
(
t
)
)
I
*
(
t
)
B
2
≥
u
2
max
⁡
.
So, the optimal control pair is characterized as (38).The optimal control pair and the state are found by solving the following optimality system, which consists of the state system (28), the adjoint system (36), boundary conditions (29) and (37), and the characterization of the optimal control pair(
u
1
*
,
u
2
*
)(38):
(43)
d
S
*
(
t
)
d
t
=
r
(
1
-
S
*
K
)
S
*
-
β
S
*
I
*
1
+
α
2
I
*
-
max
⁡
⁡
(
min
⁡
⁡
(
(
λ
1
(
t
)
-
λ
3
(
t
)
)
S
*
B
1
,
u
1
max
⁡
)
,
0
)
S
*
,
d
I
*
(
t
)
d
t
=
β
S
*
(
t
-
τ
)
I
*
(
t
-
τ
)
1
+
α
2
I
*
(
t
-
τ
)
-
(
α
+
μ
+
γ
)
I
*
-
max
⁡
⁡
(
min
⁡
⁡
(
(
λ
2
(
t
)
-
λ
3
(
t
)
)
I
*
B
2
,
u
2
max
⁡
)
,
0
)
I
*
,
d
R
*
(
t
)
d
t
=
γ
I
*
-
μ
R
*
+
max
⁡
⁡
(
min
⁡
⁡
(
(
λ
1
(
t
)
-
λ
3
(
t
)
)
S
*
B
1
,
u
1
max
⁡
)
,
0
)
S
*
+
max
⁡
(
min
⁡
⁡
(
(
λ
2
(
t
)
-
λ
3
(
t
)
)
I
*
(
t
)
B
2
,
u
2
max
⁡
)
,
k
k
k
k
k
k
k
k
k
0
min
⁡
⁡
(
(
λ
2
(
t
)
-
λ
3
(
t
)
)
I
*
(
t
)
B
2
,
u
2
max
⁡
)
)
I
*
d
λ
1
(
t
)
d
t
=
-
A
1
+
λ
1
(
t
)
(
r
(
1
-
2
S
*
K
)
-
λ
1
(
t
)
-
u
1
*
)
-
λ
3
u
1
*
-
χ
[
0
,
t
f
-
τ
]
λ
2
(
t
+
τ
)
Λ
1
,
d
λ
2
(
t
)
d
t
=
-
A
2
+
λ
1
(
t
)
Λ
2
+
λ
2
(
μ
+
α
+
γ
+
u
2
*
)
-
λ
3
(
γ
+
u
2
*
)
-
χ
[
0
,
t
f
-
τ
]
λ
2
(
t
+
τ
)
Λ
2
,
d
λ
3
(
t
)
d
t
=
λ
3
μ
withλ
1
(
t
f
)
=
0
,
λ
2
(
t
f
)
=
0,λ
3
(
t
f
)
=
0,S
(
0
)
=
S
0,I
(
0
)
=
I
0, andR
(
0
)
=
R
0, whereΛ
1
=
(
β
I
*
(
1
+
α
2
I
*
)
/
(
1
+
α
2
I
*
)
2
) and Λ
2
=
(
β
S
*
/
(
1
+
α
2
I
*
)
2
).


### 3.3. Numerical Results and Discussions

In this paragraph, we solve numerically the optimality system (43) and we present the results found. In this formulation, there exist initial conditions for the state variables and terminal conditions for the adjoint variables. That is, the optimality system is a two-point boundary value problem, with separated boundary conditions at times t
=
0 and t
=
t
f.Solving the optimality system (43) requires an iterative scheme developed by Hattaf and Yousfi [13]. This involves use of an appropriate algorithm.There exist a step sizeh
>
0 and integers (
n
,
m
)
∈
ℕ
2 with τ
=
m
h and t
f
=
n
h.For reasons of programming, we considerm knots to left of0and right oft
f, and we obtain the following partition:
(44)
Δ
=
(
t
f
t
-
m
=
-
τ
<
⋯
<
t
-
1
<
0
<
t
1
kkkk
<
⋯
<
t
n
=
t
f
<
⋯
<
t
n
+
m
)
.
Then, we have t
i
=
i
h
(
-
m
≤
i
≤
n
+
m
). Next, we define the state and adjoint variables S
(
t
),I
(
t
),R
(
t
),λ
1
(
t
),λ
2
(
t
),λ
3
(
t
), andu
1
(
t
),u
2
(
t
)in terms of nodal points S
i, I
i, R
i, λ
1
i, λ
2
i, λ
3
i, u
1
i, and u
2
i. Now, with a combination of forward and backward difference approximation, we obtain Algorithm 1.Algorithm 1
Step
1
for  i
=
-
m
,
…
,
0, do
S
i
=
S
0
,
I
i
=
I
0, R
i
=
R
0, u
1
i
=
0 and u
2
i
=
0,
end for
for  i
=
n
,
…
,
n
+
m, do
λ
1
i
=
0, λ
2
i
=
0, λ
3
i
=
0,
Step
2
S
i
+
1
=
S
i
+
h
[
r
(
1
-
S
i
K
)
S
i
-
β
S
i
I
i
1
+
α
2
I
i
-
u
1
i
S
i
],
I
i
+
1
=
I
i
+
h
[
β
S
i
-
m
I
i
-
m
1
+
α
2
I
i
-
m
]
-
(
α
+
μ
+
γ
-
u
2
i
)
I
i,
R
i
+
1 = R
i
+
h
[
γ
I
i
-
μ
R
i
+
u
1
i
S
i
+
u
2
i
I
i
],
Λ
1
i
=
[
β
I
i
(
1
+
α
2
I
i
)
(
1
+
α
2
I
i
)
2
],
Λ
2
i
=
[
β
S
i
(
1
+
α
2
I
i
)
2
],
λ
1
n
-
i
-
1
=
λ
1
n
-
i
+
h
[
-
A
1
+
λ
1
n
-
i
(
r
(
1
-
2
S
i
K
)
-
Λ
1
i
-
u
1
i
)
-
λ
3
n
-
i
u
1
i
-
χ
[
0
,
t
f
-
τ
]
(
t
n
-
i
)
λ
2
n
-
i
+
m
Λ
1
i
+
1
],
λ
2
n
-
i
-
1
=
λ
2
n
-
i
+
h
[
-
A
2
+
λ
1
n
-
i
Λ
2
i
+
1
+
λ
2
n
-
i
(
α
+
μ
+
γ
+
u
2
i
)
-
λ
3
n
-
i
(
γ
+
u
1
i
)
-
χ
[
0
,
t
f
-
τ
]
(
t
n
-
i
)
λ
2
n
-
i
+
m
Λ
2
i
+
1
],
λ
3
n
-
i
-
1
=
λ
3
n
-
i
+
h
[
μ
λ
3
n
-
i
],
Θ
1
i
+
1
=
(
λ
1
n
-
i
-
λ
3
n
-
i
)
S
i
+
1
B
1
Θ
2
i
+
1
=
(
λ
2
n
-
i
-
λ
3
n
-
i
)
I
i
+
1
B
2
u
1
i
=
max
⁡
(
min
⁡
(
Θ
1
i
+
1
,
u
1
max
⁡
)
,
0
)
u
2
i
=
max
⁡
(
min
⁡
(
Θ
2
i
+
1
,
u
2
max
⁡
)
,
0
)
Step
3
for  i
=
1
,
…
,
n,   do
write
S
*
(
t
i
)
=
S
i
,
I
*
(
t
i
)
=
I
i
,
R
*
(
t
i
)
=
R
i
,
u
1
*
(
t
i
)
=
u
1
i and u
2
*
(
t
i
)
=
u
2
i,
end forFor the simulations, we use the parameter values given in Table1.Table 1
Values of the parameters.
Parameters
Descriptions
Values
S
0
Initial susceptible population
120
I
0
Initial infected population
50
R
0
Initial recovered population
100
μ
Natural death of the population
0.01
α
Death rate due to disease
0.01
α
2
Parameter that measures the inhibitory effect
0.001
β
Transmission rate
0.0001
γ
Recovery rate
0.0004
r
Intrinsic birth rate
0.5
K
Carrying capacity
300
A
1
Weight parameter
1000
A
2
Weight parameter
1000
B
1
Weight parameter
5000
B
2
Weight parameter
10
τ
Time incubation
1In Figure2, we observe that there is a significant decrease in the number of infected individuals and susceptible individuals controlled compared with those not controlled and also an increase in the number of individuals recovered controlled.Figure 2
Evolution of different classes of individuals with and without control for time delayτ
=
1.


## 3.1. Existence of an Optimal Control

The existence of the optimal control pair can be obtained using a result by Fleming and Rishel in [11] and by Lukes in [12].Theorem 7.
There exists a control functionu
1
*
(
t
),u
2
*
(
t
) so that
(33)
J
(
u
1
*
(
t
)
,
u
2
*
(
t
)
)
=
min
⁡
(
u
1
,
u
2
)
∈
U
a
d
J
(
u
1
(
t
)
,
u
2
(
t
)
)
.Proof.
To prove the existence of an optimal control pair, it is easy to verify the following.(1)
The set of controls and corresponding state variables is nonempty.
(2)
The admissible setU
ad is convex and closed.
(3)
The right-hand side of the state system (28) is bounded by a linear function in the state and control variables.
(4)
The integrand of the objective functional is convex onU
ad.
(5)
There exists constantsω
1
,
ω
2
>
0, andρ
>
1such that the integrandL
(
S
,
I
,
u
1
,
u
2
)of the objective functional satisfiesL
(
S
,
I
,
u
1
,
u
2
)
≥
ω
2
+
ω
1
(
|
u
1
|
2
+
|
u
2
|
2
)
ρ
/
2.
The result follows directly from [11].


## 3.2. Characterization of the Optimal Control

Before characterizing the optimal control pair, we first define the Lagrangian for the optimal control problem (28)–(31) by
(34)
L
(
S
,
I
,
u
1
,
u
2
)
=
A
1
S
(
t
)
+
A
2
I
(
t
)
+
1
2
B
1
u
1
2
(
t
)
+
1
2
B
2
u
2
2
(
t
)
and the HamiltonianH for the control problem by
(35)
H
(
S
,
I
,
R
,
u
1
,
u
2
,
λ
i
,
t
)
=
L
(
S
,
I
,
u
1
,
u
2
)
+
∑
i
=
1
i
=
3
‍
λ
i
f
i
,
whereλ
i, i
=
1
,
2
,
3, are the adjoint functions to be determined suitably. Next, by applying Pontryagin’s maximum principle with delay given in [10] to the Hamiltonian H, we obtain the following theorem.Theorem 8.
Given optimal controlsu
1
*
(
t
) and u
2
*
(
t
) and solutions S
*
(
t
), I
*
(
t
), and R
*
(
t
) of the corresponding state systems (31) and (28), there exist adjoint variables λ
1, λ
2, and λ
3 that satisfy
(36)
d
λ
1
(
t
)
d
t
=
-
A
1
+
λ
1
(
t
)
(
r
(
1
-
2
S
*
K
)
-
Λ
1
-
u
1
*
)
-
λ
3
(
t
)
u
1
*
(
t
)
-
χ
[
0
,
t
f
-
τ
]
λ
2
(
t
+
τ
)
Λ
1
,
d
λ
2
(
t
)
d
t
=
-
A
2
+
λ
1
(
t
)
Λ
2
+
λ
2
(
t
)
(
μ
+
α
+
γ
+
u
2
*
)
-
λ
3
(
t
)
(
γ
+
u
2
*
)
-
χ
[
0
,
t
f
-
τ
]
λ
2
(
t
+
τ
)
Λ
2
,
d
λ
3
(
t
)
d
t
=
λ
3
(
t
)
μ
,
whereΛ
1
=
(
β
I
*
(
1
+
α
2
I
*
)
/
(
1
+
α
2
I
*
)
2
)andΛ
2
=
(
β
S
*
/
(
1
+
α
2
I
*
)
2
)
(37)
w
i
t
h
t
r
a
n
s
v
e
r
s
a
l
i
t
y
c
o
n
d
i
t
i
o
n
s
λ
i
(
t
f
)
=
0
,
i
=
1
,
2
,
3
.
Furthermore, the optimal control pairu
*
(
t
)is given by
(38)
u
1
*
(
t
)
=
max
⁡
⁡
(
min
⁡
⁡
(
(
λ
1
(
t
)
-
λ
3
(
t
)
)
S
*
(
t
)
B
1
,
u
1
max
⁡
)
,
0
)
,
u
2
*
(
t
)
=
max
⁡
⁡
(
min
⁡
⁡
(
(
λ
2
(
t
)
-
λ
3
(
t
)
)
I
*
(
t
)
B
2
,
u
2
max
⁡
)
,
0
)
.Proof.
Using Pontryagin’s maximum principle with delay in state, we obtain the adjoint equations and transversality conditions such that(39)
d
λ
1
(
t
)
d
t
=
-
∂
H
∂
S
-
χ
[
0
,
t
f
-
τ
]
∂
H
∂
S
τ
(
t
+
τ
)
,
λ
1
(
t
f
)
=
0
,
d
λ
2
(
t
)
d
t
=
-
∂
H
∂
I
-
χ
[
0
,
t
f
-
τ
]
∂
H
∂
I
τ
(
t
+
τ
)
,
λ
2
(
t
f
)
=
0
,
d
λ
3
(
t
)
d
t
=
-
∂
H
∂
R
-
χ
[
0
,
t
f
-
τ
]
∂
H
∂
R
τ
(
t
+
τ
)
,
λ
3
(
t
f
)
=
0
,
and by using the optimality conditions we find
(40)
∂
H
∂
u
1
=
B
1
u
*
(
t
)
-
λ
1
(
t
)
S
*
+
λ
3
(
t
)
S
*
=
0
,
at
u
1
=
u
1
*
(
t
)
,
∂
H
∂
u
2
=
B
2
u
*
(
t
)
-
λ
2
(
t
)
I
*
+
λ
3
(
t
)
S
*
=
0
,
at
u
2
=
u
2
*
(
t
)
,
which gives
(41)
u
1
*
(
t
)
=
(
λ
1
(
t
)
-
λ
3
(
t
)
)
S
*
(
t
)
B
1
,
u
2
*
(
t
)
=
(
λ
2
(
t
)
-
λ
3
(
t
)
)
I
*
(
t
)
B
2
.
Using the property of the control space, we obtain
(3.2)
u
1
*
(
t
)
=
0
,
if
(
λ
1
(
t
)
-
λ
3
(
t
)
)
S
*
(
t
)
B
1
≤
0
,
k
k
k
k
u
1
*
(
t
)
=
(
λ
1
(
t
)
-
λ
3
(
t
)
)
S
*
(
t
)
B
1
,
if
0
<
(
λ
1
(
t
)
-
λ
3
(
t
)
)
S
*
(
t
)
B
1
<
u
1
max
⁡
,
u
1
*
(
t
)
=
u
1
max
⁡
,
if
(
λ
1
(
t
)
-
λ
3
(
t
)
)
S
*
(
t
)
B
1
≥
u
1
max
⁡
,
u
2
*
(
t
)
=
0
,
if
(
λ
2
(
t
)
-
λ
3
(
t
)
)
I
*
(
t
)
B
2
≤
0
,
k
k
k
k
u
2
*
(
t
)
=
(
λ
1
(
t
)
-
λ
3
(
t
)
)
I
*
(
t
)
B
2
,
if
0
<
(
λ
1
(
t
)
-
λ
3
(
t
)
)
I
*
(
t
)
B
2
<
u
2
max
⁡
,
u
2
*
(
t
)
=
u
2
max
⁡
,
if
(
λ
1
(
t
)
-
λ
3
(
t
)
)
I
*
(
t
)
B
2
≥
u
2
max
⁡
.
So, the optimal control pair is characterized as (38).The optimal control pair and the state are found by solving the following optimality system, which consists of the state system (28), the adjoint system (36), boundary conditions (29) and (37), and the characterization of the optimal control pair(
u
1
*
,
u
2
*
)(38):
(43)
d
S
*
(
t
)
d
t
=
r
(
1
-
S
*
K
)
S
*
-
β
S
*
I
*
1
+
α
2
I
*
-
max
⁡
⁡
(
min
⁡
⁡
(
(
λ
1
(
t
)
-
λ
3
(
t
)
)
S
*
B
1
,
u
1
max
⁡
)
,
0
)
S
*
,
d
I
*
(
t
)
d
t
=
β
S
*
(
t
-
τ
)
I
*
(
t
-
τ
)
1
+
α
2
I
*
(
t
-
τ
)
-
(
α
+
μ
+
γ
)
I
*
-
max
⁡
⁡
(
min
⁡
⁡
(
(
λ
2
(
t
)
-
λ
3
(
t
)
)
I
*
B
2
,
u
2
max
⁡
)
,
0
)
I
*
,
d
R
*
(
t
)
d
t
=
γ
I
*
-
μ
R
*
+
max
⁡
⁡
(
min
⁡
⁡
(
(
λ
1
(
t
)
-
λ
3
(
t
)
)
S
*
B
1
,
u
1
max
⁡
)
,
0
)
S
*
+
max
⁡
(
min
⁡
⁡
(
(
λ
2
(
t
)
-
λ
3
(
t
)
)
I
*
(
t
)
B
2
,
u
2
max
⁡
)
,
k
k
k
k
k
k
k
k
k
0
min
⁡
⁡
(
(
λ
2
(
t
)
-
λ
3
(
t
)
)
I
*
(
t
)
B
2
,
u
2
max
⁡
)
)
I
*
d
λ
1
(
t
)
d
t
=
-
A
1
+
λ
1
(
t
)
(
r
(
1
-
2
S
*
K
)
-
λ
1
(
t
)
-
u
1
*
)
-
λ
3
u
1
*
-
χ
[
0
,
t
f
-
τ
]
λ
2
(
t
+
τ
)
Λ
1
,
d
λ
2
(
t
)
d
t
=
-
A
2
+
λ
1
(
t
)
Λ
2
+
λ
2
(
μ
+
α
+
γ
+
u
2
*
)
-
λ
3
(
γ
+
u
2
*
)
-
χ
[
0
,
t
f
-
τ
]
λ
2
(
t
+
τ
)
Λ
2
,
d
λ
3
(
t
)
d
t
=
λ
3
μ
withλ
1
(
t
f
)
=
0
,
λ
2
(
t
f
)
=
0,λ
3
(
t
f
)
=
0,S
(
0
)
=
S
0,I
(
0
)
=
I
0, andR
(
0
)
=
R
0, whereΛ
1
=
(
β
I
*
(
1
+
α
2
I
*
)
/
(
1
+
α
2
I
*
)
2
) and Λ
2
=
(
β
S
*
/
(
1
+
α
2
I
*
)
2
).


## 3.3. Numerical Results and Discussions

In this paragraph, we solve numerically the optimality system (43) and we present the results found. In this formulation, there exist initial conditions for the state variables and terminal conditions for the adjoint variables. That is, the optimality system is a two-point boundary value problem, with separated boundary conditions at times t
=
0 and t
=
t
f.Solving the optimality system (43) requires an iterative scheme developed by Hattaf and Yousfi [13]. This involves use of an appropriate algorithm.There exist a step sizeh
>
0 and integers (
n
,
m
)
∈
ℕ
2 with τ
=
m
h and t
f
=
n
h.For reasons of programming, we considerm knots to left of0and right oft
f, and we obtain the following partition:
(44)
Δ
=
(
t
f
t
-
m
=
-
τ
<
⋯
<
t
-
1
<
0
<
t
1
kkkk
<
⋯
<
t
n
=
t
f
<
⋯
<
t
n
+
m
)
.
Then, we have t
i
=
i
h
(
-
m
≤
i
≤
n
+
m
). Next, we define the state and adjoint variables S
(
t
),I
(
t
),R
(
t
),λ
1
(
t
),λ
2
(
t
),λ
3
(
t
), andu
1
(
t
),u
2
(
t
)in terms of nodal points S
i, I
i, R
i, λ
1
i, λ
2
i, λ
3
i, u
1
i, and u
2
i. Now, with a combination of forward and backward difference approximation, we obtain Algorithm 1.Algorithm 1
Step
1
for  i
=
-
m
,
…
,
0, do
S
i
=
S
0
,
I
i
=
I
0, R
i
=
R
0, u
1
i
=
0 and u
2
i
=
0,
end for
for  i
=
n
,
…
,
n
+
m, do
λ
1
i
=
0, λ
2
i
=
0, λ
3
i
=
0,
Step
2
S
i
+
1
=
S
i
+
h
[
r
(
1
-
S
i
K
)
S
i
-
β
S
i
I
i
1
+
α
2
I
i
-
u
1
i
S
i
],
I
i
+
1
=
I
i
+
h
[
β
S
i
-
m
I
i
-
m
1
+
α
2
I
i
-
m
]
-
(
α
+
μ
+
γ
-
u
2
i
)
I
i,
R
i
+
1 = R
i
+
h
[
γ
I
i
-
μ
R
i
+
u
1
i
S
i
+
u
2
i
I
i
],
Λ
1
i
=
[
β
I
i
(
1
+
α
2
I
i
)
(
1
+
α
2
I
i
)
2
],
Λ
2
i
=
[
β
S
i
(
1
+
α
2
I
i
)
2
],
λ
1
n
-
i
-
1
=
λ
1
n
-
i
+
h
[
-
A
1
+
λ
1
n
-
i
(
r
(
1
-
2
S
i
K
)
-
Λ
1
i
-
u
1
i
)
-
λ
3
n
-
i
u
1
i
-
χ
[
0
,
t
f
-
τ
]
(
t
n
-
i
)
λ
2
n
-
i
+
m
Λ
1
i
+
1
],
λ
2
n
-
i
-
1
=
λ
2
n
-
i
+
h
[
-
A
2
+
λ
1
n
-
i
Λ
2
i
+
1
+
λ
2
n
-
i
(
α
+
μ
+
γ
+
u
2
i
)
-
λ
3
n
-
i
(
γ
+
u
1
i
)
-
χ
[
0
,
t
f
-
τ
]
(
t
n
-
i
)
λ
2
n
-
i
+
m
Λ
2
i
+
1
],
λ
3
n
-
i
-
1
=
λ
3
n
-
i
+
h
[
μ
λ
3
n
-
i
],
Θ
1
i
+
1
=
(
λ
1
n
-
i
-
λ
3
n
-
i
)
S
i
+
1
B
1
Θ
2
i
+
1
=
(
λ
2
n
-
i
-
λ
3
n
-
i
)
I
i
+
1
B
2
u
1
i
=
max
⁡
(
min
⁡
(
Θ
1
i
+
1
,
u
1
max
⁡
)
,
0
)
u
2
i
=
max
⁡
(
min
⁡
(
Θ
2
i
+
1
,
u
2
max
⁡
)
,
0
)
Step
3
for  i
=
1
,
…
,
n,   do
write
S
*
(
t
i
)
=
S
i
,
I
*
(
t
i
)
=
I
i
,
R
*
(
t
i
)
=
R
i
,
u
1
*
(
t
i
)
=
u
1
i and u
2
*
(
t
i
)
=
u
2
i,
end forFor the simulations, we use the parameter values given in Table1.Table 1
Values of the parameters.
Parameters
Descriptions
Values
S
0
Initial susceptible population
120
I
0
Initial infected population
50
R
0
Initial recovered population
100
μ
Natural death of the population
0.01
α
Death rate due to disease
0.01
α
2
Parameter that measures the inhibitory effect
0.001
β
Transmission rate
0.0001
γ
Recovery rate
0.0004
r
Intrinsic birth rate
0.5
K
Carrying capacity
300
A
1
Weight parameter
1000
A
2
Weight parameter
1000
B
1
Weight parameter
5000
B
2
Weight parameter
10
τ
Time incubation
1In Figure2, we observe that there is a significant decrease in the number of infected individuals and susceptible individuals controlled compared with those not controlled and also an increase in the number of individuals recovered controlled.Figure 2
Evolution of different classes of individuals with and without control for time delayτ
=
1.


## 4. Conclusion

In this paper, we investigated the dynamics of a delayed SIR model with saturated incidence rate and logistic growth recruitment. At first, we obtain the existence and the stability of the equilibria by analyzing the distribution of the roots of associated characteristic equation. Using the time delay as a bifurcation parameter, we get the existence of the Hopf bifurcation whenτ crosses the critical value τ
0. Then, the model is extended to assess the impact of some control measures, by reformulating the model as an optimal control problem. Existence of the optimal control pair is established, Pontryagin’s maximum principle with delay is used to characterize these optimal controls, and the optimality system is derived. Finally, in the numerical simulation, we propose an algorithm based on the forward and backward difference approximation and we show that the optimal strategy becomes more effective when we combined the vaccination and treatment strategies together.



---

*Source: 940819-2014-05-07.xml*