Saint Petersburg State University
Department of mathematical game theory and statistical
decisions
Taynitskiy Vladislav Aleksandrovich
Master’s thesis
Optimal Program of Resistance to
Propagation of Malwares in Computer
Networks
Specialization 01.04.02
Applied Mathematics and Informatics
Master’s Program Game Theory and Operations Research
Research advisor,
C. Sc.(PhD), Associate Professor
Gubar E.A.
Saint Petersburg
2016
Contents
1 Abstract
4
2 Introduction
5
3 SIR model with continuous control 1
9
3.1
Mathematical model . . . . . . . . . . . . . . . . . . . . .
9
3.2
Objective function . . . . . . . . . . . . . . . . . . . . . .
11
3.3
Structure of the optimal control . . . . . . . . . . . . . . .
12
3.4
Numerical simulation . . . . . . . . . . . . . . . . . . . . .
18
4 SIR model with continuous control 2
21
4.1
Mathematical model . . . . . . . . . . . . . . . . . . . . .
21
4.2
Objective function . . . . . . . . . . . . . . . . . . . . . .
22
4.3
Structure of optimal control . . . . . . . . . . . . . . . . .
23
4.4
Numerical simulation . . . . . . . . . . . . . . . . . . . . .
24
4.5
Stability analysis . . . . . . . . . . . . . . . . . . . . . . .
28
5 SIR model with impulse control
34
5.1
Mathematical model . . . . . . . . . . . . . . . . . . . . .
34
5.2
Impulse control and objective function . . . . . . . . . . .
36
5.3
Structure of the optimal control . . . . . . . . . . . . . . .
39
5.4
Numerical simulation . . . . . . . . . . . . . . . . . . . . .
41
6 Conclusion
45
7 References
47
8 Appendix 1
50
2
1
Abstract
The security risks incurred by the spread of malware in computer and wireless networks can be reduced by the immunization of nodes, using security
and antivirus patches. Malware, which captures personal and corporate
confidential data, induces different damages, including costs generated by
the necessity to compensate disclosure of private information, loss of money
and social damage caused by loss of reputation. The distribution of security
and antivirus patches in a network enables the control of the proliferation of
malicious software and decreases possible losses. We formulate an optimal
control problem for the case whether two different types of malware can
circulate at the same time in a computer network and present an optimal
strategy of resistance to malwares.
4
2
Introduction
Spreading of information in computer and wireless networks has become
faster and a greater number of people use their network access for different activities, i.e. to retrieving financial information, to managing their
banking accounts to purchase goods online etc. At the same time, variety applications of networks increase the probability of security threats.
Malware is able to gain illegal access to confidential data such as bank accounts, credit cards, email and social networks passwords, collect sensitive
private information and disrupt computer functionalities. The malware
attack may lead to direct and indirect losses such as the cost of repairing software and hardware and the recovery of compromised servers. For
example, well-known successful attacks of computer viruses Kido in 2009,
MyDoom and ILOVEYOU affect millions of computers worldwide, with
approximate damage, 9.1 billion, 38 and 15 billion, respectively.
The self-propagation and replication of computer viruses are similar to those processes of biological viruses [7, 8]. In biology, a virus is an
infectious agent that replicates only inside the living cells of other organisms. Once a virus invades the host cell, it copies itself, infects the host
and leaves it. Host cells are not always killed by the action of the virus.
In fact, there are viruses that leave the infected cell alive, but use it as a
continuous media to produce generations of viruses. In a similar fashion, a
computer virus generates new copies of itself, inject itself into the code of
other programs, the system memory, and distributes its copies into a variety of communication channels. Such virus behavior can be captured by
an epidemic process which is described by the system of nonlinear differential equations as well as in classical Susceptible-Infected-Recovered (SIR)
5
model [11].
Each node in the network may be considered as susceptible while it
is not invaded by the replica of the malware transmitted from the infected
nodes. The node becomes recovered after the application of antivirus or
security patches. We assume that the protection software can effectively
protect the nodes, and they cannot be reinfected once antivirus software is
installed. The adoption of antivirus software provides a mechanism to control the propagation of malware, and hence protects the network. However,
the challenge with a wide adoption of the software is the tradeoff between
resource utilization and security risks. The scanning and monitoring process of antivirus software consumes computational resources. When the
security risks are low, it is more desirable to reduce the rate of scanning
and monitoring to achieve a better usage of computer systems.
Some viruses can generate epidemic process in population periodically, for example epidemics of influenza in urban population can reach
its peak two times during one epidemic period. Similarly, many examples
of several waves of spreading the identical malicious software in computer
and wireless networks are well-known. In paper [21], authors have analyzed
that according to a surviving probability different homogeneous groups of
viruses can preserve for long time period and provoke new attacks. As
the examples of repeated virus epidemics it can be considered the attacks
of Code Red, Code Red II and Conficker between 2001 and 2013 which
caused the damages of more than 200000 of computers worldwide [14].
Due to these reasons it is possible to use a series of impulse control actions
which can be applied in certain time moments or adhere to time interval.
For this reasons, we will consider different models describing the behaviour of control schemes of antivirus spreading. We establish an optimal
6
control framework to characterize the fundamental characteristics for network security. We use epidemic dynamics to define the virus propagation,
and a finite-horizon cost function to capture the tradeoffs between antivirus adoption and the impact of virus propagation. We aim to find an
optimal control strategy of antivirus protection that minimizes resource
consumption. In conclusion we will compare impulse and continuous types
of spreading of protection software and make conclusion about the advantages of both types of antiviruses.
Our work is related to previous studies in this area, including [3,
5, 9, 23]. In current study, the network can be attacked by a heterogeneous source of malware which captures the fact of coexistence of different
types of exploits and vulnerability of the existing computing systems. The
challenge for modeling heterogeneous malware spreading is multi-fold. We
consider a network attacked by two types of malware [9], where each node
can be infected either separately by each type of malware or by both types
simultaneously. First, we model the dynamics of propagation of the malware in the network in case if both types coexist in one host node. Second,
we formulate an optimal control problem and show the structure of the
optimal strategies, which provides the minimum of the aggregated system
costs depends on the properties of value-functions.
The work is organized as follows. Sections 3 and 4 present mathematical models of epidemics and formulates the optimal control problem in
case of continuous control. In subsections 3.3 and 4.3, using Pontryagin’s
maximum principle, we define the structure optimal control policies and
proof main results. Subsection 4.4 focuses on the stability analysis of the
uncontrolled system. Numerical examples will be presented in subsection
3.4 and 4.5. Section 5 describes the mathematical model of spreading of
7
viruses as impulse control problem. Structure of impulse control problem
and numerical simulation are presented in subsections 5.3 and 5.4. Section
6 concludes the work and main results.
8
3
3.1
SIR model with continuous control 1
Mathematical model
In this section, we study a network of N nodes, where two types of malicious
software spread with different speeds. Malicious software propagates very
fast, hence Susceptible-Infected-Recovery (SIR) model needs to be adapted
to describe the epidemics of viruses in computer networks. All nodes in the
network are divided into three groups: Susceptible (S), Infected (I) and
Recovered (R) [11]. Susceptible is a group of nodes which have not contact
with infected nodes yet and may be invaded by any forms of malicious
software. Infected nodes are already attacked by the virus and Recovered
is a group of dispatched nodes. Since two types of malware circulate in
the network, the Infected consist of the subgroup of nodes infected by the
first form of malware V1 , the subgroup of nodes infected by the second
form V2 and the group infected by both forms of viruses. We model the
epidemic process as a system of nonlinear differential equations, where nS ,
nV1 , nV2 , nR correspond to the number of susceptible, infected by different
forms of malware V1 , V2 and recovered nodes, respectively. The variable
nV12 is a number of nodes simultaneously infected by both viruses. The
total number of nodes in the network during the entire process remains
constant and equal to N , nS + nV1 + nV2 + nV12 + nR = N .
Let S(t) =
nS (t)
N ,
I1 (t) =
nV1 (t)
N ,
I2 (t) =
nV2 (t)
N ,
I12 =
nV12 (t)
N ,
R(t) =
nR (t)
N
as a fraction of the Susceptible, the Infected by virus V1 , V2 , both viruses
together and the Recovered nodes, respectively. At the beginning of the
epidemics t = 0, the majority of the nodes belong to the Susceptible group,
and a small fraction of nodes is infected by different types of malware.
Hence initial states are given by 0 < S(0) = S 0 , 0 < I1 (0) = I10 ,
9
0
0
0 < I2 (0) = I20 , 0 < I12 (0) = I12
, R(0) = R0 = 1 − S 0 − I10 − I20 − I12
.
A susceptible node becomes infected whenever it accepts or installs
malicious software, spread by the infected nodes. Malware spreads in the
network at the rate of β1 and β2 for V1 and V2 , respectively. In particular, if
the node is infected by the malware V1 , then with probability εβ2 it can be
infected by the second malware V2 , and vice versa, if the node is infected
with a virus V2 , then with probability εβ1 it may be infected by V1 . Here,
the variable ε ∈ [0, 1] is a probability that a node infected by one type
virus will be infected by another type virus. If a susceptible node contacts
with a node infected by both viruses V1 and V2 , then with probability βi it
may be invaded by only one form of malware Vi , i = 1, 2.
Usually, a majority of nodes are protected by permanent antivirus
software which is effective against known viruses. Then, we can consider
a recovery rate γ, which show the probability that susceptible nodes are
recovered by permanent antivirus software. However, periodically, the epidemics of new computer viruses appear, and the permanent antivirus software is often inefficient against new or unknown malicious software. In
this case, special patches can be applied to protect the network. We define ui (t), i = 1, . . . , 4 as a control parameter which corresponds to the
application of special antivirus patches. Parameter σ1 , σ2 and σ3 may
be interpreted as self-recovery rate in biological system. We may discuss
about the rate of self-recovery if a node is repaired without any external
resources or employees.
We modeled the epidemics of two forms of malware using a system
10
Figure 1: The scheme of the propagation of two forms of malwares.
of nonlinear differential equations:
Ṡ = −β1 S(I1 + I12 ) − β2 S(I2 + I12 ) − γSu4 ,
I˙1 = β1 S(I1 + I12 ) − εβ2 I1 (I2 + I12 ) − σ1 I1 − (1 − εβ2 )u1 I1 ,
I˙2 = β2 S(I2 + I12 ) − εβ1 I2 (I1 + I12 ) − σ2 I2 − (1 − εβ1 )u2 I2 ,
I˙12 = β1 I2 (I1 + I12 ) + εβ2 I1 (I2 + I12 ) − σ3 I12 − u3 I12 ,
(1)
Ṙ = γSu4 + (1 − εβ2 )u1 I1 + (1 − εβ1 )u2 I2 + u3 I12 + σ1 I1 +
+σ2 I2 + σ3 I12 ,
S + I1 + I2 + I12 + R = 1.
Hence initial states are
0
S(0) = S 0 , I1 (0) = I10 , I2 (0) = I20 , I12 (0) = I12
,
(2)
0
R(0) = R = 1 − S(0) − I1 (0) − I2 (0) − I12 (0).
3.2
Objective function
The objective of the system designer is to minimize the aggregated cost
on time interval [0, T ]. At any given t, the overall system costs include
11
infected costs f1 (I1 (t)), f2 (I2 (t)), f3 (I12 (t)) and protection costs hi (ui (t)),
i = 1, . . . , 4. Infected costs are includes losses caused by infected nodes;
protection costs are generated by the consumption of resources for the
application of antivirus or stationary security patches. Functions fi are
non-decreasing and twice-differentiable, such as fi (0) = 0, fi (Ii ) > 0 for
Ii > 0. Functions hi (ui (t)) are twice differentiable and increasing in ui (t)
such as hi (0) = 0, hi (ui ) > 0, when ui > 0, ui ∈ [0, 1] for i = 1, . . . , 4.
Aggregated system costs are defined as the functional:
RT
J(u1 , u2 , u3 , u4 ) = (f1 (I1 (t)) + f2 (I2 (t)) + f3 (I3 (t))
0
(3)
+h1 (u1 (t)) + h2 (u2 (t)) + h3 (u3 (t)) + h4 (u4 (t)))dt.
The optimal control problem is to minimize the functional in a time
interval [0, T ], i.e.,
min
u1 ,u2 ,u3 ,u4
3.3
J(u1 , u2 , u3 , u4 ).
(4)
Structure of the optimal control
We will use Pontryagin’s maximum principle [16] to solve the problem
described above. There exist continuous and piecewise continuously differentiable co-state functions λi that at every point t ∈ [0; T ] where ui (t),
i = 1, . . . , 4 are continuous, satisfy (8) and (9).
(u1 , . . . , u4 ) ∈ arg max H = (λ, S, I1 , I2 , I12 , R, u1 , . . . , u4 ).
u1 ,...,u4
Here, u1 , u2 , u3 and u4 are feasible controls.
12
(5)
Define functions ψi , i = 1, . . . , 4 as follows:
ψ1 = (1 − εβ2 )(λR − λI1 )I1 ,
ψ2 = (1 − εβ1 )(λR − λI2 )I2 ,
ψ3 = (λR − λI12 )I12 ,
ψ4 = γS(λR − λS ).
(6)
Hamiltonian H of the system (1):
H = −(f1 (I1 ) + f2 (I2 ) + f3 (I12 ) + h1 (u1 ) + h2 (u2 )+
h3 (u3 ) + h4 (u4 )) + β1 S(I1 + I12 )(λI1 − λS )+
β2 S(I2 + I12 )(λI2 − λS ) + εβ2 I1 (I2 + I12 )(λI12 − λI1 )+
εβ1 I2 (I1 + I12 )(λI12 − λI2 ) + (λR − λI1 )σ1 I1 +
(7)
(λR − λI2 )σ2 I2 + (λR − λI12 )σ3 I12 +
(1 − εβ2 )u1 I1 (λR − λI1 ) + (1 − εβ1 )u2 I2 (λR − λI2 )+
u3 (λR − λI12 )I12 + (λR − λS )γSu4 .
Adjoint system is
dλS
= β1 (λS − λI1 )(I1 + I12 ) + β2 (λS − λI2 )(I2 + I12 )−
dt
γu4 (λR − λS );
dλI1
df1 (I1 )
=
+ β1 S(λS − λI1 ) + εβ2 (I2 + I12 )(λI1 − λI12 )+
dt
dI1
εβ1 I2 (λI2 − λI12 ) + (λI1 − λR )σ1 I1 + u1 (1 − εβ2 )(λI1 − λR );
dλI2
df2 (I2 )
=
+ β2 S(λS − λI2 ) + εβ1 (I1 + I12 )(λI2 − λI12 )+
dt
dI2
εβ2 I1 (λI1 − λI12 ) + (λI2 − λR )σ2 I2 + u2 (1 − εβ1 )(λI2 − λR );
dλI12
df3 (I12 )
=
+ β1 S(λS − λI1 ) + β2 S(λS − λI2 )+
dt
dI12
(λI12 − λR )σ3 I12 + εβ2 I1 (λI1 − λI12 ) + εβ1 I2 (λI2 − λI12 )+
u3 (λI12 − λR );
dλR
= 0.
dt
13
(8)
together with the transversality conditions
λS (T ) = λI1 (T ) = λI2 (T ) = λI12 (T ) = λR (T ) = 0.
(9)
According to [3, 5, 16] we construct the optimal program to prevent
the spreading of malicious software in a computer network.
Proposition 1 1) If hi is strictly convex function (h00i (ui ) > 0) then there
exist time moments t0 , t1 ∈ [0, T ], 0 ≤ t0 ≤ t1 ≤ T such that:
0,
if t1 < t < T ;
ui = h0i −1 (ψi ), if t0 < t ≤ t1 ;
1,
if 0 < t < t0 .
(10)
2) If hi is concave function (h00i (ui ) < 0) then there exist time moment
t1 ∈ [0, T ] such that:
ui =
0,
if t1 < t < T ;
1,
if 0 < t < t1 .
(11)
where functions ψi , i = 1, . . . , 4 are follows:
Proof of the Proposition 1. Rewrite Hamiltonian in terms of
function ψi and we obtain:
H = −(f1 (I1 ) + f2 (I2 ) + f3 (I12 )) + β1 S(I1 + I12 )(λI1 − λS )+
β2 S(I2 + I12 )(λI2 − λS ) + εβ2 I1 (I2 + I12 )(λI12 − λI1 ) + εβ1 I2 (I1 + I12 )+
(λR − λI1 )σ1 I1 + (λR − λI2 )σ2 I2 + (λR − λI12 )σ3 I12 +
(−h1 (u1 ) + u1 ψ1 ) + (−h2 (u2 ) + u2 ψ2 ) + (−h3 (u3 ) + u3 ψ3 )+
(−h4 (u4 ) + u4 ψ4 ).
(12)
14
According to the algorithm of maximum principle we consider derivatives:
∂H
= −ḣi (ui ) + ψi = 0, i = 1, . . . , 4.
∂ui
(13)
As hi (ui ) are increasing functions and Ii ≥ 0 and S ≥ 0 then Hamiltonian reaches its maximum if ψi = ḣi (ui ) ≥ 0, i = 1, . . . , 4. It follows
if and only if the following conditions are satisfied: (λR (t) − λI1 (t)) ≥ 0,
(λR (t) − λI2 (t)) ≥ 0, (λR (t) − λI12 (t)) ≥ 0 and (λR (t) − λS (t)) ≥ 0. To
complete the proof of proposition we consider auxiliary lemma.
Lemma 1 For all 0 ≤ t ≤ T , we have (λR (t) − λI1 (t)) ≥ 0,
(λR (t) − λI2 (t)) ≥ 0, (λR (t) − λI12 (t)) ≥ 0 and (λR (t) − λS (t)) ≥ 0.
Lemma 1 is proved in the similar way to those in [3], [9] and it is
based on the following properties.
Property 1 Let v(t) be a continuous and piecewise differential function
of t. Let v(t1 ) = L and v(t) > L for all t ∈ (t1 , . . . , t0 ]. Then v̇(t+
1 ) ≥ 0.
Where v(t+
1 ) = lim v(x).
x→0
Property 2 For any convex and differentiable function y(x), which is 0
at x = 0, y 0 (x)x − y(x) ≥ 0 for all x ≥ 0.
Proof of Lemma 1.
Lets split our proof into two parts. At the first part we will consider the
case when t = T and show that derivatives of the functions (λR (t)−λI1 (t)),
(λR (t) − λI2 (t)), (λR (t) − λI12 (t)) and (λR (t) − λS (t)) are less or equal to
zero to prove that they are non-increasing at t = T . In the second part we
will use the method of proof by contradiction and show that on the whole
interval [0, T ] these functions are also non-negative.
15
I. At time moment T , we have according to (9):
(λR (T ) − λI1 (T )) = 0;
(λR (T ) − λI2 (T )) = 0;
(λR (T ) − λI12 (T )) = 0;
(λR (T ) − λS (T )) = 0.
(14)
From (8) we receive
(λ̇R (T ) − λ̇I1 (T )) = −f˙1 (I1 (T )) ≤ 0;
(λ̇R (T ) − λ̇I2 (T )) = −f˙2 (I2 (T )) ≤ 0;
(λ̇R (T ) − λ̇I12 (T )) = −f˙3 (I12 (T )) ≤ 0;
(15)
(λ̇R (T ) − λ̇S (T )) = 0.
Now we have that at time moment T all functions are equal to 0 and
their derivatives are less or equal to 0 then we can obtain that
(λR (t) − λI1 (t)), (λR (t) − λI2 (t)), (λR (t) − λI12 (t)) and (λR (t) − λS (t)) are
decreasing functions at t = T .
II. (Proof by contradiction)
Let 0 ≤ t∗ ≤ T be the last instant moment at which one of the
inequality constraints are satisfied:
(λR (t) − λI1 (t)) ≥ 0, (λR (t) − λI2 (t)) ≥ 0,
(λR (t) − λI12 (t)) ≥ 0, (λR (t) − λS (t)) ≥ 0.
First, suppose that following inequalities are hold
(λR (t∗ ) − λI1 (t∗ )) = 0,
(λR (t∗ ) − λI2 (t∗ )) ≥ 0,
(λR (t∗ ) − λI12 (t∗ )) ≥ 0,
(λR (t∗ ) − λS (t∗ )) ≥ 0,
∗
∗
∗
∗
(λS (t ) − λI1 (t )) ≥ 0,
(λS (t ) − λI2 (t )) ≥ 0
(λI1 (t∗ ) − λI12 (t∗ )) ≥ 0,
(λI2 (t∗ ) − λI12 (t∗ )) ≥ 0.
16
(16)
From (8) we have:
df1
− β1 S(λS − λI1 )−
dI1
εβ2 (I2 + I12 )(λI1 − λI12 ) − εβ1 I2 (λI2 − λI12 )−
λ̇R (t∗ ) − λ̇I1 (t∗ ) = −
(17)
(λI1 − λR )σ1 I1 − u1 (1 − εβ2 )(λI1 − λR ).
Since (λR (t)−λI1 (t)) is decreasing function on the interval [0; T ] then
according to Property 1, we consider a time moment t∗+ such as:
df1
− β1 S(λS − λI1 )−
dI1
εβ2 (I2 + I12 )(λI1 − λI12 ) − εβ1 I2 (λI2 − λI12 ).
λ̇R (t∗+ ) − λ̇I1 (t∗+ ) = −
(18)
According to assumption (16), the difference (λ̇R (t) − λ̇I1 (t)) is negative at the time moment t∗+ , then function (λR (t) − λI1 (t)) decreases. It
contradicts Property 1 and proves that our function is not increasing on
interval t ∈ [0; T ]. We get that (λR (t) − λI1 (t)) ≥ 0 for t ∈ [0; T ].
Second part of the proof is similar to presented above and presented
in Appendix 1.
The proof of Lemma 1 is completed.
Now we return to the main proposition and consider two cases which
depend on the properties of hi (ui ), i = 1, . . . , 4.
1) hi is concave. Since functions hi are concave (h00i < 0,
i = 1, . . . , 4), then (ui ψi − hi (ui )) are convex functions of ui . Hamiltonian
H is a strictly convex function according to (7) and for any t ∈ [0, T ] and
it reaches its maximum either at ui = 1 or ui = 0, i = 1, . . . , 4.
ui (t) =
0,
if ψi < hi (1);
1,
if ψi > hi (1).
17
(19)
2) hi is strictly convex. If functions hi are strictly convex (h00i > 0,
i = 1, . . . , 4) then (ui ψi − hi (ui )) and Hamiltonian H are concave functions
of ui , therefore
∂H
∂ui
= −ḣi (ui ) + ψi = 0 and
ui (t) =
0,
dh−1
i (ψi )
dui
1 ,
if ψi ≤
,
if
dhi (0)
dui
if ψi >
dhi (0)
dui ;
< ψi ≤
dhi (1)
dui ;
(20)
dhi (1)
dui .
The proof of the Proposition 1 is completed.
3.4
Numerical simulation
In this section, we use numerical simulations to support our theoretical
results. We use the following values and parameters: iteration step is
δ = 0.25, initial fractions of nodes are S(0) = 0.55, I1 (0) = 0.15,
I2 (0) = 0.2, I12 (0) = 0.1, R(0) = 0. Intensive rates of transition from
susceptible to infected and recovered are β1 = 0.6, β2 = 0.7 and γ = 0.1
respectively. Virus interaction rate is ε = 0.8. Intensive rates of transition
from infected to recovered are σ1 = 0.03, σ2 = 0.02 and σ3 = 0.01.We
suppose that epidemic duration is T = 15 time units.
Figure 2: Uncontrolled SIR model.
18
First, we consider uncontrolled case. Since we do not apply any
antivirus patches to reduce the epidemics, then 80 % of the hosts are
transferred into the group infected I12 at time T = 15, this is shown by the
solid line in the Fig. 2. 20 % of the nodes are healed of the malware by
self-recovery parameters σi in the system (dot line).
Figure 3: Controlled SIR model. Costs functions hi are convex.
Also, we illustrate the case when costs for treatment measures hi (ui )
are convex functions (Fig. 3). We use h1 (u1 ) = 0.35u21 ; h2 (u2 ) = 0.4u22 ;
h3 (u3 ) = 0.5u23 ; h4 (u4 ) = 0.2u24 .
Figure 4: Optimal control in SIR model (functions hi are convex).
19
At the last time moment T = 15 we can see that there are no infected
hosts, most hosts have transferred into recovered group and only some
hosts are in a susceptible group. This fraction of susceptible have not been
impact to the epidemic process. The structure of optimal control is shown
in the Fig. 4.
The next diagram (Fig. 5) illustrates the comparison between aggregated system costs for uncontrolled and controlled system, which are equal
to J = 775.5 and J = 112.1 monetary units respectively.
Figure 5: Comparison of aggregated system costs for controlled and uncontrolled SIR
models.
20
4
4.1
SIR model with continuous control 2
Mathematical model
As in previous model (1) we study a network of homogenous N nodes, with
the assumption that two forms of malicious software spread in the network
with different speeds. In contrast to the previous model, we assume that
in this case there is no permanent antivirus software (u4 (t) = 0). Selfrecovery parameters σi which are describing transferring probabilities from
infected to susceptible state. From that reason self-recovered node can be
infected again.
Fig. 6 represents the evolution of the propagation malware in the
network.
Figure 6: The scheme of the propagation of two forms of malware.
We have extended classical SIR model, and modeled the epidemics
of two forms of malware using a system of nonlinear differential equations
[5]:
21
Ṡ = −β1 S(I1 + I12 ) − β2 S(I2 + I12 ) + σ1 I1 + σ2 I2 ;
I˙1 = β1 S(I1 + I12 ) − εβ2 I1 (I2 + I12 ) − σ1 I1 + σ2 I12 − (1 − εβ2 )u1 I1 ;
I˙2 = β2 S(I2 + I12 ) − εβ1 I2 (I1 + I12 ) − σ2 I2 + σ1 I12 − (1 − εβ1 )u2 I2 ;
I˙12 = εβ1 I2 (I1 + I12 ) + εβ2 I1 (I2 + I12 ) − (σ1 + σ2 )I12 − u3 I12 ;
Ṙ = (1 − εβ2 )u1 I1 + (1 − εβ1 )u2 I2 + u3 I12 .
(21)
4.2
Objective function
As in subsection 3.2, we find optimal control parameters which minimize
aggregated system costs over the time interval [0, T ]. At any time moment t we define infected costs fi (Ii (t)) for infected nodes and protection
costs hi (ui (t)), i = 1, 2, 3. Functions fi are non-decreasing and twicedifferentiable, such as fi (0) = 0, fi (Ii ) > 0 for Ii (t) > 0, i = 1, 2, 3.
Functions hi (ui (t)) are twice differentiable and increasing in ui (t) such as
hi (0) = 0, hi (ui ) > 0, when ui (t) > 0, u ∈ [0, 1]. g(R) is benefit rate, it is
non-decreasing and differentiable function and g(0) = 0.
Aggregated system costs are defined by the following functional
J(u1 , u2 , u3 ) =
RT
0
(f1 (I1 (t)) + f2 (I2 (t)) + f3 (I12 (t)) + h1 (u1 (t))+
(22)
h2 (u2 (t)) + h3 (u3 (t)) − g(R))dt,
and the optimal control problem is to minimize the functional on
time interval [0, T ], i.e.,
min J(u1 , u2 , u3 ).
u1 ,u2 ,u3
22
(23)
4.3
Structure of optimal control
Based on previous work [3] and [9], the optimal program was constructed
and it has the following structure:
Proposition 2 1. If hi is strictly convex function then exists time moments t0 , t1 ∈ [0, T ], 0 ≤ t0 ≤ t1 ≤ T such as:
ui (t) =
0,
if ψi ≤
dh−1
i (ψi )
dui
,
1 ,
if
dhi (0)
dui
if ψi >
dhi (0)
dui ;
< ψi ≤
dhi (1)
dui ;
(24)
dhi (1)
dui ;
2. If hi is concave function then exists time moment t ∈ [0, T ] such as:
ui (t) =
0,
if ψi < hi (1);
1,
if ψi > hi (1).
(25)
where functions ψi are defined as follows:
ψ1 = (1 − εβ2 )(λR − λI1 )I1 ,
ψ2 = (1 − εβ1 )(λR − λI2 )I2 ,
(26)
ψ3 = (λR − λI12 )I12 .
Proof of Proposition 2
Proof of the Proposition 2 is similar as in previous section model. We use
Pontryagin’s maximum principle to find optimal control. The proof is in
the Appendix 2.
23
4.4
Numerical simulation
In this section, we present numerical examples to study the SIR model. We
use the following values and parameters: iteration step is δ = 0.15, initial
fractions of nodes are S(0) = 0.55, I1 (0) = 0.15, I2 (0) = 0.2, I12 (0) = 0.1,
R(0) = 0. Intensive rates of transition from susceptible to infected and
recovered are β1 = 0.6, β2 = 0.7. Virus interaction rate is ε = 0.8. Intensive
rates of transition from infected to recovered are σ1 = 0.03, σ2 = 0.02 and
σ3 = 0.01.We suppose that epidemic duration is T = 15 time units.
First, we consider uncontrolled case (Fig. 7). Since we do not apply
any antivirus patches to reduce the epidemics, then 90 % of the hosts are
transferred into the group infected I12 at time T = 15, this is shown by the
solid line in the Fig. 2. There are no recovered nodes like at uncontrolled
case in previous model because self-recovery heal nodes and transfer them
to susceptible state. After some time they infects again.
Figure 7: Uncontrolled SIR model.
System costs in uncontrolled case are equal to J = 302.5 monetary
units (Fig. 8).
In contrast to previous case, we apply antivirus patches to heal in24
Figure 8: System costs in uncontrolled SIR model. (J = 302.5 monetary units)
fected nodes. Simulation represents that almost all nodes are in recovered
state and there are no infected nodes (Fig. 9).
Figure 9: Controlled SIR model (hi - strictly convex).
Aggregated costs reduces to J = 130.88 monetary units (Fig. 10).
The next diagram (Fig. 11) illustrates the comparison between system coasts for uncontrolled and controlled system, which are equal to
J = 302.5 and J = 130.88 monetary units respectively.
Structure of the optimal control is presented in the Fig. 12.
25
Figure 10: System costs in controlled SIR model. (J = 130.88 monetary units)
Figure 11: Comparison of system costs with and without spreading of antivirus.
Next, we change our functions hi from strictly convex to concave.
h1 (u1 (t)) = 0.36 − (u1 (t) − 0.6)2 ;
h2 (u2 (t)) = 0.49 − (u2 (t) − 0.7)2 ;
(27)
h3 (u3 (t)) = 0.64 − (u3 (t) − 0.8)2 .
As in previous case, almost all of the nodes have moved into the
recovered subgroup (Fig. 13).
26
Figure 12: The optimal security response. Protection costs functions are strictly convex
h1 (u1 ) = 0.35u21 ; h2 (u2 ) = 0.4u22 ; h3 (u3 ) = 0.5u23 .
Figure 13: Controlled SIR model (hi - concave).
Aggregated system costs in case when hi are concave functions are
J = 193.17 monetary units (Fig. 14).
Structure of the optimal control is presented in the Fig. 15.
27
Figure 14: Aggregated costs in controlled SIR model. (J = 193.17 monetary units)
Figure 15: The optimal security response. Protection costs functions are strictly convex
h1 (u1 (t)) = 0.36−(u1 −0.6)2 ; h2 (u2 (t)) = 0.49−(u2 −0.7)2 ; h3 (u3 (t)) = 0.64−(u3 −0.8)2 .
4.5
Stability analysis
In this section, we study the stability of the equilibrium states for the
uncontrolled system, where ui = 0, i = 1, 2, 3, [18, 22]. At an equilibrium
point, all derivatives
dS
dt
=
dI1
dt
=
dI2
dt
=
28
dI12
dt
=
dR
dt
= 0. We consider system
(6) without control. Therefore, the subgroup of Recovered R is eliminated.
Ṡ = −β1 S(I1 + I12 ) − β2 S(I2 + I12 ) + σ1 I1 + σ2 I2 ;
I˙1 = β1 S(I1 + I12 ) − εβ2 I1 (I2 + I12 ) − σ1 I1 + σ2 I12 ;
I˙2 = β2 S(I2 + I12 ) − εβ1 I2 (I1 + I12 ) − σ2 I2 + σ1 I12 ;
(28)
I˙12 = εβ1 I2 (I1 + I12 ) + εβ2 I1 (I2 + I12 ) − (σ1 + σ2 )I12 ;
Thus we have three equilibrium points:
• I1 = I2 = I12 = 0;
• I1 = I12 = 0,
I2 = 1 − σ2 /β2 ;
• I2 = I12 = 0,
I1 = 1 − σ1 /β1 .
Thus, we can find a simple equation for I12 :
εI1 I2 (β1 + β2 ) = (σ1 + σ2 + u3 − εβ2 I1 − εβ1 I2 )I12
(29)
Lemma 2 The number of people infected by both virus 1 and virus 2 will
obey the following equation:
dI12
dt
= εβ1 I2 (I1 + I12 ) + εβ2 I1 (I2 + I12 ) − (σ1 + σ2 )I12 = 0;
(30)
hence
I12 = (εI1 I2 (β1 + β2 ))/(σ1 + σ2 + u3 − εβ2 I1 − εβ1 I2 )
From (28) we get:
dS
dt
= −β1 S(I1 + I12 ) − β2 S(I2 + I12 ) + σ1 I1 + σ2 I2 = 0,
(31)
hence
β1 S(I1 + I12 ) + β2 S(I2 + I12 ) = σ1 I1 + σ2 I2
Hence we have the expected three equilibrium points [5]:
29
• S = 1, I1 = I2 = I12 = 0;
• I1 = I12 = 0, I2 = 1 − σ2 /β2 , S = σ1 /β1 ;
• I2 = I12 = 0, I1 = 1 − σ1 /β1 , S = σ2 /β2 .
To provide the local stability analysis, we linearize the system (28).
Together with the condition S(t) = 1 − I1 (t) − I2 (t) − I12 (t), and we obtain
Jacobi matrices for each equilibrium point, and use the Routh-Hurwitz
criterion to determine the stability [10].
According to the criterion, we construct the following characteristic
equation
z 4 + A1 z 3 + A2 z 2 + A3 z 1 = 0
(32)
and examine that all eigenvalues of the characteristic equation (32) has
negative real part Re(z) < 0. It follows if the determinants of Hurwitz
matrices are positive. Ai , i = 1, 2, 3 are constructed according to the
criterion.
I. We will check the stability of equilibrium points using Hurwitz
criterion. Consider the first state (1, 0, 0, 0, 0). Hurwitz matrix for this
case:
0 −β1 + σ1 −β2 + σ2 −β1 − β2
0 β1 − σ1
0
β1 + σ2
0
0
β2 − σ2
β2 + σ1
0
0
0
−σ1 − σ2
30
(33)
Let’s define Ai , i = 1 . . . 4 as
A1 = −a11 − a22 − a33 − a44 ;
A2 = a11 (a22 + a33 + a44 ) + a22 (a33 + a44 ) + a33 a44 −
a23 a32 ;
A3 = a11 a23 a32 + a23 a32 a44 − a11 a22 a33 −
(34)
a11 a22 a44 − a11 a33 a44 −
a22 a33 a44 − a13 a21 − a32
A4 = a11 a22 a33 a44 + a13 a21 a31 a44 −
a11 a23 a32 a44
where aij is an element of Hurwitz matrix. We calculate functions Gi as
follows:
G1 = A1 ;
G2 = A1 A2 − A3 :
(35)
2
2
G3 = A1 A2 A3 − (A1 ) A4 − (A3 ) ;
G4 = G2 (A3 A4 ) − (A1 A4 )2 .
From calculation we receive that all Gi are positive, then our system in
first equilibrium point is stable.
II. Jacobi matrix at E1 (S, I1 , I2 , I2 , R), I1 = I12 = 0, I2 = 1 − σ2 /β2 ,
S = σ1 /β1 is
−β1 µ1
β1 µ 1
0
0
0
0
− ββ2 σ1 1
+ σ2
−εβ2 µ1
0 ( ββ21 − 1)σ1 − εβ1 µ1
0
ε(β2 + β1 )µ1
−η1 σ1
σ1 − εµ1 + σ2
η1 σ1
εβ2 µ1 − σ1 − σ2
To simplify the notation we denote η1 = ( ββ12 + 1), µ1 = (1 −
31
.
σ1
β1 ).
All determinants are positive Gi > 0, i = 1, 2, 3, 4, and then this
equilibrium state is stable.
III. Jacobi matrix at E2 (S, I1 , I2 , I2 , R), I2 = I12 = 0, I1 = 1 −
S=
σ1
β1
σ1
β1 ,
is
− ββ1 σ2 2
−β2 µ2
− σ1
0
( ββ12 − 1)σ2 − εβ2 µ2
β2 µ 2
−εβ1 µ2
0
ε(β1 + β2 )µ2
where µ2 = (1 −
σ2
β2 ),
0
− ββ1 σ2 2
− σ2
0
+ σ2
,
0 σ2 − εβ1 µ2 + σ1
0 εβ1 µ2 − σ1 − σ2
− ββ1 σ2 2
η2 = ( ββ12 + 1).
From calculation we get that all determinants are positive Gi > 0,
i = 1, . . . , 4, and hence this equilibrium state is stable.
Using direct Lyapunov’s method [13], we can check the global stability of E0 . We consider the Lyapunov function as follows:
L = S + I1 + I2 + I12 + R.
The derivative
dL
dt
dL
dt |(1)
(36)
under the system (28) is
= Ṡ + I˙1 + I˙2 + I˙12 + Ṙ = −β1 S(I1 + I12 ) − β2 S(I2 + I12 )+
σ1 I1 + σ2 I2 + β1 S(I1 + I12 ) − εβ2 I1 (I2 + I12 )−
σ1 I1 + σ2 I12 − (1 − εβ2 )u1 I1 + β2 S(I2 + I12 )−
εβ1 I2 (I1 + I12 ) − σ2 I2 + σ1 I12 −
(1 − εβ1 )u2 I2 + εβ1 I2 (I1 + I12 )+
εβ2 I1 (I2 + I12 ) − (σ1 + σ2 )I12 − u3 I12 +
(1 − εβ2 )u1 I1 + (1 − εβ1 )u2 I2 + u3 I12 ≤ 0.
32
(37)
Hence according to Lyapunov’s stability theorem we have that stationary state is stable but not asymptotically stable.
Proposition 3 Stationary states E0 (S, 0, 0, 0, R), E1 (S, I1 , I2 , I2 , R),
E2 (S, I1 , I2 , I2 , R) of the model (28) are stable.
33
5
5.1
SIR model with impulse control
Mathematical model
Based on the model (1) we describe a model which include the system of
differential equation to describe behavior of viruses and discrete system of
impulses. As in previous research [9], [20] we consider a multi-virus case in
which two forms of malicious software spread in the network with different
speeds. Viruses can generate epidemic process in population periodically,
for example epidemics of influenza can reach its peak two times during one
epidemic period. Many examples of several waves of spreading the identical
malicious software in computer and wireless networks are well-known [21].
In current work we formulate the conditions for eradication of epidemics of
malwares for different cases of protection policies and compare costs and
effectiveness of impulse actions and standard method of resistance.
Figure 16: The scheme of the propagation of the system 38.
As in previous work [20] we model the epidemics of two forms of
malware using a system of nonlinear differential equations (38).
34
Ṡ = −β1 S(I1 + I12 ) − β2 S(I2 + I12 ),
I˙1 = β1 S(I1 + I12 ) − εβ2 I1 (I2 + I12 ) − σ1 I1 ,
I˙2 = β2 S(I2 + I12 ) − εβ1 I2 (I1 + I12 ) − σ2 I2 ,
(38)
I˙12 = β1 I2 (I1 + I12 ) + εβ2 I1 (I2 + I12 ) − σ3 I12 ,
Ṙ = σ1 I1 + σ2 I2 + σ3 I12 ,
together with the condition
S(t) + I1 (t) + I2 (t) + I12 (t) + R(t) = 1, t ∈ [0; T ].
(39)
However, the analysis of the behavior of computer viruses has shown
that a small shares of infected nodes might be survived in the Internet and
if the local network has a connection with the Global Network then the
epidemics resumes [21]. Then we can say the IT-security deals with the
repeated waves of epidemics of malicious software. The iterative epidemic
process can be formulated as a combined multi-virus model with series
of impulses which control the shares of Infected nodes. As a basis for
complex model can be used the system (38) or its closed modification i.e.
Susceptible-Infected-Recovered-Susceptible (SIRS) model.
Ṡ = −β1 S(I1 + I12 ) − β2 S(I2 + I12 ) + µR,
I˙1 = β1 S(I1 + I12 ) − εβ2 I1 (I2 + I12 ) − σ1 I1 ,
I˙2 = β2 S(I2 + I12 ) − εβ1 I2 (I1 + I12 ) − σ2 I2 ,
I˙12 = β1 I2 (I1 + I12 ) + εβ2 I1 (I2 + I12 ) − σ3 I12 ,
(40)
Ṙ = σ1 I1 + σ2 I2 + σ3 I12 − µR,
S + I1 + I2 + I12 + R = 1.
Here µ is a rate at which recovered nodes becomes new susceptible nodes
in the network. This SIRS model can be helpful for describing epidemics
35
that are repeated periodically. Application of antivirus in impulse control
form allows us to reduce the amount of infected nodes.
Figure 17: The scheme of the propagation of the system 40.
5.2
Impulse control and objective function
Together with initial systems (38), (40), which describe the behavior of
viruses on the time intervals (τp−1 , τp ), we formulate a model with application of impulse control strategies. To protect the network from repeated
epidemics special patches can be applied as series of impulses at certain
time moments. By using these patches as the control impulses at time
moments τ1 , . . . , τk we receive the extended system of differential equation
to describe the case of spreading of two malwares for all time period except
the sequence of time moments τp . States of the system after time moments
τp are
36
Ṡ(τp+ ) = Ṡ(τp );
I˙1 (τp+ ) = I˙1 (τp ) − u1 (τp )I1 (τp ),
I˙2 (τp+ ) = I˙2 (τp ) − u2 (τp )I2 (τp ),
(41)
I˙12 (τp+ ) = I˙12 (τp ) − u3 (τp )I3 (τp ),
Ṙ(τp+ ) = Ṙ(τp ) + u1 (τp )I1 (τp ) + u2 (τp )I2 (τp ) + u3 (τp )I12 (τp ).
We define ui (·), i = 1, 2, 3 as a control parameter which corresponds to the application of special antivirus patches at time moments
τ1 , . . . , τk . At each time moment ui is a fraction of treated nodes. Here
u1 = (u11 , . . . , uk1 ), u2 = (u12 , . . . , uk2 ), u3 = (u13 , . . . , uk3 ) are components of
control vectors correspond to the set of time moments τ1 , . . . , τk ,
uj1 ∈ [0, uj1 ], uj2 ∈ [0, uj2 ], uj3 ∈ [0, uj3 ], where uj1 , uj2 , uj3 are maximum values
of control.
The objective function of the combined system (41) is constructed
to evaluate the aggregated costs on a time interval [0, T ] including the
costs of control impulses. Aggregated costs for continuous systems (38)
and (40) are defined as follows: at any given t 6= τp , p = 1, . . . , k, the
overall system costs include infected costs f1 (I1 (t)), f2 (I2 (t)), f3 (I12 (t)).
Functions fi are non-decreasing and twice-differentiable, such as fi (0) = 0,
fi (Ij (t)) > 0 for Ii (t) > 0 for t ∈ (τp−1 , τp ) for all j, i = 1, 2, 3. For
system (41) we define infected costs as functions hi (upi (τp+ )), p = 1, . . . , k,
where hi (upi (τp+ )) > 0, upi (τp+ ) > 0, ui ∈ [0, upi ] for i = 1, 2, 3 which are
generated by the consumption of resources for the application of antivirus
or stationary security patches. Infected costs are consist of damages caused
by viruses.
37
Aggregated system costs are defined as the functional:
J(u1 , u2 , u3 ) =
RT
f1 (I1 (t)) + f2 (I2 (t)) + f3 (I3 (t))dt+
0
k
P
h1 (u1 (τp+ )) +
p=1
k
P
h2 (u2 (τp+ )) +
p=1
k
P
(42)
h3 (u3 (τp+ )).
p=1
To illustrate the special properties of impulse treatment strategies we
examine some simple examples of SIR and SIRS models which assess the
process of propagation of malicious softwares. Pulse treatment is effective
if we succeed to keep the number of Susceptible below a critical value
which is generated by the envelope curves Lj (t) (index j corresponds to
the enumeration of viruses) or certain critical values [1].
Figure 18: Behaviour of the infected nodes and critical values for infected in corresponding
subgroups I1 , I2 . System parameters are: step δ = 0.25, initial states S(0) = 0.75,
I1 (0) = 0.1, I2 (0) = 0.15, I1 2 = 0, ε = 0.5, β1 = 0.35, β1 = 0.4, σ1 = 0.002, σ2 = 0.004,
µ = 0.4. Maximal fraction of treated nodes are: u1 = 0.2, u2 = 0.3.
Figures 18, 19 illustrates the SIRS model (40) in case of impulse
treatment spreading which preserve a number of infected of viruses V1 and
V2 below critical values.
38
Figure 19: Fractions of S and R in SIRS model.
5.3
Structure of the optimal control
We construct the envelope based on the concept of basic reproduction
number which is defined as the expected number of new infections from
one infected individual in a fully susceptible population through the entire
duration of the infectious period [1], [2]. This metric is significant as well it
helps to determine whether or not an infectious disease can spread through
a population. The replacement number is defined as the expected number
of secondary infections that one infected person would produce through
the entire duration of the infectious period [15]. Reproduction numbers for
virus V1 , V2 and both viruses simultaneously are defined respectively:
R01 =
β2
β1
, R02 = ,
σ1
σ2
(43)
R03 =
ε(β1 + β2 )
.
σ3
39
Then envelop curves we define as follows on the interval [0, T ]:
L1 (t) = R01 S(t),
(44)
L2 (t) = R02 S(t),
L3 (t) = R03 S(t).
Pulse treatment is effective to keep the the stable distribution of Susceptible, Infected and Recovered in population below envelope curves that
define critical values. From that reasons according to system (38) and
from the definition of envelops Lj (t) we have:
β1
S(t),
σ1
β2
I2 (t) < S(t),
σ2
ε(β1 + β2 )
I12 (t) <
S(t).
σ3
I1 (t) <
(45)
In current model, we define the series of time moments τp , p = 1, . . . , k
at which we apply the special patches to treat the infected nodes. As
functions Ij (t) are increasing for all j (in the case without any control
measures), we turn impulse of control on for the first time when
Ij (t) = Lj (t), for all j. Series of control impulses switch on at time moments
τp , p = 1, . . . , k according to following conditions:
0, I1 (t) < L1 (t),
u1 (τp ) =
up , I (t) ≥ L (t),
(46)
0, I2 (t) < L2 (t),
u2 (τp ) =
up , I (t) ≥ L (t),
(47)
1
1
2
2
40
1
2
0, I12 (t) < L3 (t),
u3 (τp ) =
up , I (t) ≥ L (t).
12
3
(48)
3
Here upi ∈ [0, 1], i = 1, 2, 3, p = 1, . . . , k is defined as maximal value of
applied control impulses.
5.4
Numerical simulation
In this section we present numerical simulations that corroborate our results. In the example we set following values of parameters: iteration step is
δ = 0.25, initial fractions of nodes are S(0) = 0.45, I1 (0) = 0.2, I2 (0) = 0.3,
I12 (0) = 0.05, R(0) = 0. Intensive rates of transition from susceptible to
infected and recovered are β1 = 0.35, β2 = 0.4 and γ = 0.1 respectively.
Virus interaction rate is ε = 0.5. Intensive rates of transition from infected
to recovered are σ1 = 0.05, σ2 = 0.04 and σ3 = 0.03. We suppose that
epidemic duration is T = 25 time units. Infected costs for infected nodes
are f1 (I1 (t)) = 5I1 (t), f2 (I2 (t)) = 6I2 (t) and f3 (I12 (t)) = 10I12 (t).
Figure 20: Multi-virus SIR model with continuous control.
In continuous control case, we use as protection costs following func41
tions: h1 (u1 (t)) = 1.3u1 (t)2 , h2 (u2 (t)) = 1.5u22 (t), h3 (u3 (t)) = 2u23 (t) and
h4 (u4 (t)) = u24 (t), aggregated system costs is J = 65.36 monetary units
(Fig. 22). At the time T = 25 we can see that there are no infected
hosts (Fig 20). The distribution of nodes in population shows that the
most of nodes have been transferred into Recovered subgroup and a small
shares of nodes is in a susceptible subgroup, which have no influence on
the epidemics. Optimal control for this case is shown in the Fig 21.
Figure 21: Optimal control for continuous case.
Figure 22: System costs for continuous case (J = 65.36).
In case with impulse control, protection costs are h1 (u1 (τp )) = 10u1 (τp ),
42
h2 (u2 (τp )) = 15u2 (τp ), h3 (u3 (τp )) = 18u3 (τp ). At the end of time interval
T = 25 we can see that the most of infected hosts become recovered
(Fig. 23) but in contract to continuous case aggregated costs is bigger.
This fact appears from the formulation of combined problem and means
that in continuous case we treat the infected nodes only once while in
impulse case we deal with repeated procedures.
Figure 23: Multi-virus SIR model with series of impulses.
System costs J = 252.51 monetary units are illustrated in the Fig.24.
Figure 24: Aggregated system costs for the combined system with impulse controls.
43
In the next diagram (Fig. 25) presented curves Li , i = 1 . . . 3 and
fractions of infected nodes I1 , I2 and I12 . We apply pulse treatment at time
moments τp when Ii (t) = Li (t).
Figure 25: Envelope curves Li for the proportion of infected nodes Ii .
44
6
Conclusion
In this work we have studied four modifications of multi-viruses SusceptibleInfected-Recovered and Susceptible-Infected-Recovered-Susceptible models which describe propagation of two types of malware in computer networks. These extended models take into account the coexistence of heterogeneous malware and the exposure of computer systems to multiple
vulnerabilities.
Firstly, we study case when we spread antivirus in continuous form.
We have formulated an optimal control problem to study the tradeoffs
between security risks and the control investment. By using Pontryagin’s
maximum principle, we have obtained different control policies structure
that minimize the aggregated cost. The structure of the control depends
on properties of costs functions. Numerical simulations were performed
using a specially written procedures.
In second part, we have reformulated the SIR and SIRS models under the impulse control. In contrast to the previous statement we analyze
the conditions for application the series of impulses that protect the network during periodic waves of epidemics of malwares instead of continuous
control. This case is also supported with numerical simulations.
In subsection 5.4, according to the numerical simulation for the set of
initial data and parameters, it has been shown that the aggregated system
costs in continuous case are less than in impulse treatment case. Due to
selected conditions of applying treatment, protection strategies in impulse
form starts to treat our system not immediately but after some time from
the beginning of epidemics, while protection strategies of continuous form
of spreading starts to heal infected nodes immediately at time moment
45
t = 0. As a result malicious software has time to infect susceptible nodes
in system which provokes additional system costs.
In further researches, continuous control SIR and SIRS models can
be reformulated as impulse optimal control problem. These models will
cover the studying of stability and finding of the optimal control in model
with pulse treatment.
46
References
1. Agur Z., Cojocaru L., Mazor G., Anderson R. M., Danon Y. L. Pulse
mass measles vaccination across age cohorts // Proceedings of the National Academy of Sciences of the United States of America, 1993. Vol.
90, P. 11698—11702.
2. Agur Z., Stone L., Shulgin B. Pulse Vaccination Strategy in the SIR
Epidemic Model // Bulletin of Mathematical Biology, 1998. Vol. 60, P.
1123—1148.
3. Altman E., Khouzani M., Sarkar S. Optimal control of epidemic evolution // Proceedings of INFOCOM, 2011. P. 1683—1691
4. Behncke H. Optimal control of deterministic epidemics // Optimal Control Applications and Methods, 2000. Vol. 21, P. 269—285.
5. Beutel A., Faloutsos C., Prakash B., Rosenfeld R. Interacting Viruses
in Networks: Can Both Survive? // KDD-2012, 2012. P.426—434.
6. Chitnis N. Introduction to Mathematical Epidemiology: The Basic Reproductive Number // Universitat Basel Mathematisches Institut, 2011.
7. Evans A. S., Kaslow R. A. Viral Infections of Humans: Epidemiology
and Control // New York: Springer, 1997.
8. Funke B., Tortora G., Gerard J., Christine L. Case. Microbiology: An
Introduction. San Francisco, CA: Benjamin Cummings, 2004. P. 848.
9. Gubar E., Zhu Q. Optimal Control of Influenza Epidemic Model with
Virus Mutations // 12th Biannual European Control Conference, IEEE
Control Systems Society, 2012. P. 3125—3130.
47
10. Hurwitz A. Ueber die Bedingungen, unter welchen eine Gleichung nur
Wurzeln mit negativen reellen Theilen besitzt // Mathematische Annalen, 1895. Vol. 46, P. 273—-284.
11. Kermack W., Mc Kendrick A. A contribution to themathematical theory of epidemics // Proceedings of the Royal Society, 1927. P. 700—721.
12. Khouzani M. H. R., Sarkar S., Altman E. Maximum damage malware
attack mobile wireless networks // Proc. of the 29th International Conference on Computer Communications (INFOCOM), 2010. P. 749—757.
13. Lyapunov A. M. The General Problem of the Stability of Motion. Doctoral dissertation, Univ. Kharkov, 1892. P. 250.
14. Moore D., Shannon C. Code-Red: a Case Study on the Spread and Victims of an Internet Worm // Proceedings of the 2002 ACM SICGOMM
Internet Measurement Workshop, 2002. PP. 273––284.
15. Nakul C. Introduction to Mathematical Epidemiology: The Basic Reproductive Number // Einfuhrung in die Mathematische Epidemiologie,
2011. P. 5.
16. Pontryagin L., Boltyanskii V., Gamkrelidze R., Mishchenko E. The
Mathematical Theory of Optimal Processes, Russia: Interscience, 1962.
P. 360.
17. Rohloff K. R., Başar T. Deterministic and stochastic models for the
detection of random constant scanning worms // ACM Transactions on
Modeling and Computer Simulation (ACM TOMACS), 2008. Vol. 18,
PP. 8:1—8:24.
48
18. Sharma S., Samanta G. P. Stability analysis and optimal control of an
epidemic model with vaccination // International Journal of Biomathematics, 2015. Vol. 8, No. 3, P. 28.
19. Taylor P. D., Jonker L. B. Evolutionarily stable strategies and game
dynamics // Mathematical Biosciences, 1978. Vol. 40, PP.145—156.
20. Taynitskiy V.A., Gubar E.A., Zhitkova E.M. Structure of optimal control in the model of propagation of two malicious softwares // International conference ”Stability and Control Processes” in memory of V.I.
Zubov (SCP), 2015. PP. 261-264.
21. Pastor-Satorras R., Vespignani A. Epidemic Spreading in Scale-Free
Networks // Physical Review Letters, 2001. Vol. 86, PP.3200—3203.
22. Weibull J. Evolutionary Game Theory. The M.I.T. Press, Cambridge,
MA, 1995. P. 265.
23. Wickwire K. H. A note on the optimal control of carrier-borne epidemics // Journal of Applied Probability, 1975. Vol. 12, PP. 565—568.
24. Wu Q., Small M., Liu H.. Superinfection Behaviors on Scale-Free Networks with Competing Strains // Journal of Nonlinear Science, 2013.
Vol. 23, PP. 113—127.
49
Appendix 1
Now we have to prove that (λR (t) − λS (t)) ≥ 0. Analogously we
assume:
(λR (t∗ ) − λI1 (t∗ )) ≥ 0, (λR (t∗ ) − λI2 (t∗ )) ≥ 0,
(λR (t∗ ) − λI12 (t∗ )) ≥ 0, (λR (t∗ ) − λS (t∗ )) = 0,
∗
∗
∗
∗
(49)
(λS (t ) − λI1 (t )) ≥ 0, (λS (t ) − λI2 (t )) ≥ 0
(λI1 (t∗ ) − λI12 (t∗ )) ≥ 0, (λI2 (t∗ ) − λI12 (t∗ )) ≥ 0.
From (8) we obtain:
λ̇R (t∗ ) − λ̇S (t∗ ) = −β1 (λS − λI1 )(I1 + I12 )−
(50)
−β2 (λS − λI2 )(I2 + I12 ) + γu4 (λR − λS ).
According to assumption (81), the difference λ̇R (t∗+ ) − λ̇S (t∗+ ) ≤ 0
is negative at the time moment t∗+ , then function (λR (t) − λS (t)) decrease.
It contradicts Property 1. At time moment t = T (λR (t) − λS (t)) is equal
to zero and on the whole interval it is decreasing. We proved that
(λR (t) − λS (t)) ≥ 0 for t ∈ [0; T ].
Analogously we can prove that (λR (t) − λI2 (t)) ≥ 0 and
(λR (t) − λI12 (t)) ≥ 0.
We get that on the interval t ∈ [0; T ]:
(λR (t) − λI1 (t)) ≥ 0;
(λR (t) − λI2 (t)) ≥ 0;
(λR (t) − λI12 (t)) ≥ 0;
(λR (t) − λS (t)) ≥ 0.
50
(51)
Appendix 2
Proof of Proposition 2.
We use Pontryagin’s maximum principle [4], [16] to find the optimal control
u = (u1 , u2 , u3 ) which yields the minimum solution to the functional (22)
for the problem described above. Consider the Hamiltonian
H = −(f1 (I1 ) + f2 (I2 ) + f3 (I12 ) + h1 (u1 )+
h2 (u2 ) + h3 (u3 ) − g(R)) + β1 S(I1 + I12 )(λI1 − λS )+
β2 S(I2 + I12 )(λI2 − λS ) + εβ2 I1 (I2 + I12 )(λI12 − λI1 )+
εβ1 I2 (I1 + I12 )(λI12 − λI2 ) + σ1 (λS − λI1 )I1 +
(52)
σ2 (λS − λI2 )I2 + σ1 (λI2 − λI12 )I12 +
σ2 (λI1 − λI12 )I12 + (1 − εβ2 )u1 I1 (λR − λI1 )+
(1 − εβ1 )u2 I2 (λR − λI2 ) + u3 (λR − λI12 )I12 .
Let be λ0 = 1. We construct the adjoint system as follows:
dλS
= β1 (λS − λI1 )(I1 + I12 ) + β2 (λS − λI2 )(I2 + I12 );
dt
dλI1
= f10 + β1 S(λS − λI1 ) + εβ1 I2 (λI2 − λI12 ) + σ1 (λI1 − λS )+
dt
εβ2 (I2 + I12 )(λI1 − λI12 ) + u1 (1 − εβ2 )(λI1 − λR );
dλI2
= f20 + β2 S(λS − λI2 ) + εβ2 I1 (λI1 − λI12 ) + σ2 (λI2 − λS )+
dt
εβ1 (I1 + I12 )(λI2 − λI12 ) + u2 (1 − εβ1 )(λI2 − λR );
dλI12
= f30 + β1 S(λS − λI1 ) + β2 S(λS − λI2 ) + σ1 (λI12 − λI2 )+
dt
εβ2 I1 (λI1 − λI12 ) + εβ1 I2 (λI2 − λI12 ) + σ2 (λI12 − λI1 )+
(53)
+u3 (λI12 − λR );
dλR
= −g 0 (R),
dt
together with condition R(t) = 1 − S(t) − I1 (t) − I2 (t) − I12 (t), along with
51
the transversality conditions
λS (T ) = λI1 (T ) = λI2 (T ) = λI12 (T ) = λR (T ) = 0.
(54)
According to Pontryagin’s maximum principle, there exist continuous
and piecewise continuously differentiable co-state functions λi that at every
point t ∈ [0; T ] where u1 , u2 and u3 are continuous, satisfy (53) and (54).
In addition, we have λ(t) = (λS (t), λI1 (t), λI2 (t), λI12 (t), λR (t)), and
(u1 , u2 , u3 ) ∈ arg max H(λ, S, I1 , I2 , I12 , R, u1 , u2 , u3 ).
u1 ,u2 ,u3
(55)
Here, u1 , u2 , u3 are feasible controls.
Rewrite Hamiltonian in terms of function ψ, and we obtain
H = −(f1 (I1 ) + f2 (I2 ) + f3 (I12 ) − g(R))+
β1 S(I1 + I12 )(λI1 − λS ) + β2 S(I2 + I12 )(λI2 − λS )+
εβ2 I1 (I2 + I12 )(λI12 − λI1 ) + εβ1 I2 (I1 + I12 )(λI12 − λI2 )+
(56)
σ1 (λS − λI1 )I1 + σ2 (λS − λI2 )I2 + σ1 (λI2 − λI12 )I12 +
σ2 (λI1 − λI12 )I12 + (−h1 (u1 ) + u1 ψ1 )+
(−h2 (u2 ) + u2 ψ2 ) + (−h3 (u3 ) + u3 ψ3 ).
For any admissible control u1 , u2 and u3 and according to (56) for all
t ∈ [0, T ] we arrive at
(−h1 (u1 ) + u1 ψ1 − h2 (u2 ) + u2 ψ2 − h3 (u3 ) + u3 ψ3 ) ≥
(−h1 (u1 ) + u1 ψ1 − h2 (u2 ) + u2 ψ2 − h3 (u3 ) + u3 ψ3 ).
52
(57)
Hence we obtain
(u1 (t), u2 (t), u3 (t)) ∈ arg max (−h1 (x) + xψ1 −
x,y,z∈[0,1]
(58)
−h2 (y) + yψ2 − h3 (z) + zψ3 ).
Then,
max (−h1 (u1 ) + u1 ψ1 − h2 (u2 ) + u2 ψ2 − h3 (u3 ) + u3 ψ3 )
u1 ,u2 ,u3
= max(−h1 (u1 ) + u1 ψ1 ) + max(−h2 (u2 ) + u2 ψ2 )
u1
u2
(59)
+ max(−h3 (u3 ) + u3 ψ3 ).
u3
According to the algorithm of Pontryagin’s maximum principle to
∂H
:
determine the optimal control structure, we consider derivatives
∂u
∂H
∂ui
= −ḣi (ui ) + ψi = 0,
(60)
As hi (ui ), i = 1, 2, 3 are increasing functions and Ii ≥ 0 then Hamiltonian reaches its maximum if ḣi (ui ) = ψi ≥ 0.
Let us calculate time derivatives of functions ψi :
ψ̇1 = (1 − εβ2 )[(λR − λI1 )I˙1 + (λ̇IR − λ̇I1 )I1 ] =
(1 − εβ2 )[(λR − λI1 )(β1 S(I1 + I12 ) − εβ2 I1 (I2 + I12 )−
σ1 I1 + σ2 I12 − (1 − εβ2 )u1 I1 )+
(−g 0 (R) − (f10 + β1 S(λS − λI1 )+
εβ1 I2 (λI2 − λI12 ) + εβ2 (I2 + I12 )(λI1 − λI12 )+
σ1 (λI1 − λS ) + u1 (1 − εβ2 )(λI1 − λR )))I1 ].
53
(61)
ψ̇2 = (1 − εβ1 )[(λR − λI2 )I˙2 + (λ̇IR − λ̇I2 )I2 ] =
(1 − εβ2 )[(λR − λI2 )(β2 S(I2 + I12 ) − εβ1 I2 (I1 + I12 )−
σ2 I2 + σ1 I12 − (1 − εβ1 )u2 I2 )+
(−g 0 (R) − (f20 + β2 S(λS − λI2 ) + εβ2 I1 (λI1 − λI12 )+
σ2 (λI2 − λS )+
εβ1 (I1 + I12 )(λI2 − λI12 ) + u2 (1 − εβ1 )(λI2 − λR )))I2 ].
ψ̇3 = [(λR − λI12 )I˙12 + (λ̇IR − λ̇I12 )I12 ] =
(λR − λI12 )(εβ1 I2 (I1 + I12 ) + εβ2 I1 (I2 + I12 )−
(σ1 + σ2 )I12 − u3 I12 )+
(−g 0 (R) − (f30 + β1 S(λS − λI1 ) + β2 S(λS − λI2 )+
εβ2 I1 (λI1 − λI12 ) + εβ1 I2 (λI2 − λI12 ) + σ2 (λI12 − λI1 )+
σ1 (λI12 − λI2 ) + u3 (λI12 − λR )))I12 .
After arrangement:
ψ̇1 = (1 − εβ2 )[(λR − λI1 )I˙1 + (λ̇IR − λ̇I1 )I1 ]
(1 − εβ2 )[I1 (−f10 − g 0 + β1 S(λR − λS )+
εβ2 (I2 + I12 )(λI12 − λR ) + εβ1 I2 (λI12 − λI2 )+
σ1 (λS − λR )) + I12 (λR − λI1 )(β1 S + σ2 )]
ψ̇2 = (1 − εβ1 )[(λR − λI2 )I˙2 + (λ̇IR − λ̇I2 )I2 ] =
(1 − εβ1 )[I2 (−f20 − g 0 + β2 S(λR − λS )+
εβ1 (I1 + I12 )(λI12 − λR ) + εβ2 I1 (λI12 − λI1 )+
σ2 (λS − λR )) + I12 (λR − λI2 )(β2 S + σ1 )]
ψ̇3 = [(λR − λI12 )I˙12 + (λ̇IR − λ̇I12 )I12 ] =
I12 [−f30 − g 0 + (λR − λI2 )(εβ1 I2 − σ1 )+
(λR − λI1 )(εβ2 I1 − σ2 ) + S(β1 (λI1 − λS )+
β2 (λI2 − λS ))] + εI1 I2 (β1 (λR − λI12 )+
β2 (λR − λI12 )).
54
(62)
Lemma 3 ψi (t), i = 1, 2, 3 are decreasing functions of t ∈ [0, T ).
Lemma 3 is proved using following the similar methodology to those
presented in [3], [9]. We can consider two cases:
1) hi is concave.
Since functions h1 , h2 and h3 are concave (h00i < 0,
i = 1, 2, 3), then
(ui ψi − hi (ui )), i = 1, 2, 3 are convex functions of ui . Hamiltonian H is
a strictly convex function according to (56) and for any t ∈ [0, T ] and it
reaches its maximum either at ui = 1 or ui = 0, i = 1, 2, 3
ui (t) =
0,
if ψi < hi (1);
1,
if ψi > hi (1).
(63)
2) hi is strictly convex.
If functions hi are strictly convex (h00i > 0, i = 1, 2, 3) then
(−hi (ui ) + uψi , i = 1, 2, 3) and Hamiltonian is concave function, then
( dH
dui = −ḣi (ui ) + ψi = 0, ui ∈ [0, 1], i = 1, 2, 3). Then,
ui (t) =
0,
dh−1
i (ψi )
dui
if ψi ≤
,
1 ,
if
dhi (0)
dui
if ψi >
dhi (0)
dui ;
< ψi ≤
dhi (1)
dui ;
(64)
dhi (1)
dui ,
functions ψi , h0i , ui are continuous at all t ∈ [0, T ]. In this case hi is strictly
convex and h0i is strictly increasing function then h0i (0) < h0i (1). Thus there
exists such moments t0 , t1 (0 < t0 < t1 < T ) such as conditions (64) are
satisfied if ψi is as described above.
To complete the proof of proposition we consider auxiliary lemma.
It follows from (60) that Hamiltonian reaches maximum if and only if the
55
following conditions are satisfied: (λR − λI1 ) ≥ 0, (λR − λI2 ) ≥ 0 and
(λR − λI12 ) ≥ 0.
Lemma 4 For all 0 ≤ t ≤ T , we have (λR − λI1 ) ≥ 0, (λR − λI12 ) ≥ 0 and
(λR − λI2 ) ≥ 0.
Lemma 4 is proved in the similar way to those in [3], [9] and it is
based on the following two properties described before (Property 1 and 2).
Proof of Lemma 4
I. At time T , we have (λR (T ) − λI1 (T )) = 0,
(λR (T ) − λI2 (T )) = 0 and (λR (T ) − λI12 (T )) = 0 according to (54).
Consider the derivatives
(λ̇R (T ) − λ̇I1 (T )) = −f˙1 (I1 (T )) ≤ 0,
(λ̇R (T ) − λ̇I2 (T )) = −f˙2 (I2 (T )) ≤ 0,
(65)
(λ̇R (T ) − λ̇I12 (T )) = −f˙3 (I12 (T )) ≤ 0.
Moreover we may say that λi (t) ≥ 0, because λ̇Ii (t) ≥ 0 and
(λR − λI1 ) ≥ 0, (λR − λI12 ) ≥ 0 are positive on the open interval [0, T ].
Now we have that functions ψi (T ) = 0 and also we may say that
ψi (t) are positive in t ∈ [0; T ].
II.(Proof by contradiction).
Let 0 ≤ t∗ ≤ T be the last instant moment at which one of the inequality
constraints are satisfied, i.e.:
(λR − λI1 ) ≥ 0, (λR − λI2 ) ≥ 0, (λR − λI12 ) ≥ 0
and
(λR − λI1 ) = 0 or (λR − λI2 ) = 0 or (λR − λI12 ) = 0.
56
First, suppose that following inequality satisfy
(λR (t∗ ) − λI1 (t∗ )) = 0, (λR (t∗ ) − λI2 (t∗ )) ≥ 0,
(λR (t∗ ) − λI12 (t∗ )) ≥ 0, (λS (t∗ ) − λI1 (t∗ )) ≥ 0,
(66)
(λI1 (t∗ ) − λI12 (t∗ )) ≥ 0, (λI2 (t∗ ) − λI12 (t∗ )) ≥ 0.
We have to prove that (λR (t)−λI1 (t)) are non-decreasing function on
the interval [0; T ]. According to Property 1, we consider a time moment
t∗+ :
λ̇R (t∗+ ) − λ̇I1 (t∗+ ) = −f10 − β1 S(λS − λI1 ) − εβ1 I2 (λI2 − λI12 )
(67)
−u1 (1 − εβ2 )(λI1 − λR ) − εβ2 (I2 + I12 )(λI1 − λI12 ).
As well as f1 (I1 ) is non-decreasing function and all parameters are
non-negatives hence we have
d
∗+
dt (λR (t )
− λI1 (t∗+ )) ≤ 0. This contradicts
Property 1.
The system of ODE is autonomous, i.e., Hamiltonian and the constraints on the control do not have an explicit dependency on the independent variable t
H(S(t), I1 (t), I2 (t), I12 (t), R(t), u1 (t), u2 (t),
(68)
λS (t), λI1 (t), λI2 (t), λI12 (t), λR (t)) = const.
From (52) and (54) we obtain
H = H(T ) = −(f1 (I1 (T )) + f2 (I2 (T )) + f3 (I12 (T ))+
(69)
h1 (u1 (T )) + h2 (u2 (T )) + h3 (u3 (T ))).
Since fi and hi are non-decreasing functions and according to (52),
57
we receive
λR (t) − λI1 (t) =
1
(H + f1 (I1 ) + f2 (I2 ) + f3 (I12 )+
(1 − εβ2 )u1 I1
h1 (u1 ) + h2 (u2 ) + h3 (u3 ) − β1 S(I1 + I12 )(λI1 − λS )−
β2 S(I2 + I12 )(λI2 − λS ) − εβ2 I1 (I2 + I12 )(λI12 − λI1 )−
(70)
(1 − εβ1 )(λR − λI2 )u2 I2 − u3 (λI12 − λR )−
εβ1 I2 (I1 + I12 )(λI12 − λI2 ))
Moreover fi (I) is a non-decreasing function, then
fi (Ii (T )) > fi (Ii (t)) > 0, where Ii (T ) > 0 according to general assumptions
and we have
H + f2 (I2 (t)) + f3 (I12 (t)) + h1 (u1 (t)) + h2 (u2 (t))+
h3 (u3 (t)) − β1 S(t)(I1 (t) + I12 (t))(λI1 − λS )−
β2 S(t)(I2 (t) + I12 (t))(λI2 − λS )−
(71)
εβ2 I1 (t)(I2 (t) + I12 (t))(λI12 − λI1 )−
εβ1 I2 (t)(I1 (t) + I12 (t))(λI12 − λI2 )−
ψ2 u2 (t) − ψ3 u3 (t) ≤ −f1 (I1 (T )) + ψ3 u1 (T ) ≤ 0.
This follows from assumptions for functions fi (Ii ), i = 1, 2, 3 and
hi (ui ), i = 1, 2, 3, such as Ii (T ) > 0 then fi (Ii ) > 0 and ui (t) > 0 then
hi (ui ) > 0.
58
From (54) and (71) we have
λ̇R (t∗+ ) − λ̇I1 (t∗+ ) = −f10 − β1 S(λS − λI1 ) − εβ2 (I2 + I12 )(λI1 − λI12 )−
εβ1 I2 (λI2 − λI12 ) −
1
I1 (H
+ f1 (I1 ) + f2 (I2 ) + f3 (I12 ) + h1 (u1 )+
h2 (u2 ) + h3 (u3 ) − β1 S(I1 + I12 )(λI1 − λS ) − β2 S(I2 + I12 )(λI2 − λS )−
εβ2 I1 (I2 + I12 )(λI12 − λI1 ) − εβ1 I2 (I1 + I12 )(λI12 − λI2 )−
ψ2 u2 − ψ3 u3 ) =
1
˙
I1 (−f1 I1
+ f1 ) +
1
I1 (H
+ f2 (I2 ) + f3 (I12 ) + h1 (u1 ) + h2 (u2 ) + h3 (u3 )−
β1 S(I1 + I12 )(λI1 − λS ) − β2 S(I2 + I12 )(λI2 − λS )−
εβ2 I1 (I2 + I12 )(λI12 − λI1 ) − εβ1 I2 (I1 + I12 )(λI12 − λI2 ) − ψ2 u2 −
ψ3 u3 ) − β1 S(λS − λI1 ) − εβ2 (I2 + I12 )(λI1 − λI12 ) − εβ1 I2 (λI2 − λI12 ).
(72)
Here the infected cost function f1 (I1 ) is convex increasing function
and f1 (0) = 0, I1 > 0, also we have f1 (I2 ) − f˙1 (I1 )I1 ≤ 0 by Property 2.
From (54), (1) we can show that (λ̇R (t) − λ̇I1 (t)) ≤ 0 and it contradicts
Property 1, hence part II of the lemma follows.
In this part we will prove that (λR (t) − λI2 (t)) > 0. Analogously to
the part II of the Lemma 4 we suppose that
(λR (t∗ ) − λI1 (t∗ )) ≥ 0, (λR (t∗ ) − λI2 (t∗ )) = 0,
(λR (t∗ ) − λI12 (t∗ )) ≥ 0, (λS (t∗ ) − λI1 (t∗ )) ≥ 0,
(λI1 (t∗ ) − λI12 (t∗ )) ≥ 0, (λI2 (t∗ ) − λI12 (t∗ )) ≥ 0.
59
(73)
We can see that
H + f1 (I1 ) + f3 (I12 ) + h1 (u1 (t)) + h2 (u2 (t)) + h3 (u3 )−
−β1 S(t)(I1 (t) + I12 (t))(λI1 − λS ) − β2 S(t)(I2 (t)+
+I12 (t))(λI2 − λS ) − εβ2 I1 (t)(I2 (t) + I12 (t))(λI12 − λI1 )−
(74)
−εβ1 I2 (t)(I1 (t) + I12 (t))(λI12 − λI2 ) − (1 − εβ2 )(λR −
−λI1 )u1 I1 − (λR − λI12 )u3 I12 ≤ −f2 (I2 (T ))+
+(1 − εβ1 )(λR (T ) − λI2 (T ))u2 I2 (T ) ≤ 0.
At the time moment t∗+ we have:
λ̇R (t∗+ ) − λ̇I2 (t∗+ ) = −f20 − β2 S(λS − λI2 )−
εβ1 (I1 + I12 )(λI2 − λI12 ) − εβ2 I1 (λI1 − λI12 )+
1
I2 (H
+ f1 (I1 ) + f2 (I2 ) + f3 (I12 ) + h1 (u1 )+
(75)
h2 (u2 ) + h3 (u3 ) − β1 S(I1 + I12 )(λI1 − λS )−
β2 S(I2 + I12 )(λI2 − λS ) − εβ2 I1 (I2 + I12 )(λI12 − λI1 )−
εβ1 I2 (I1 + I12 )(λI12 − λI2 ) − ψ1 u1 − ψ3 u3 =
λ̇R (t∗+ ) − λ̇I2 (t∗+ ) =
1
˙
I2 (−f2 I2
+ f2 ) +
1
I2 (H
+ f1 (I1 )+
f3 (I12 ) + h1 (u1 ) + h2 (u2 ) + h3 (u3 )−
β1 S(I1 + I12 )(λI1 − λS ) − β2 S(I2 + I12 )(λI2 − λS )−
(76)
εβ2 I1 (I2 + I12 )(λI12 − λI1 ) − εβ1 I2 (I1 + I12 )(λI12 − λI2 )−
ψ1 u1 − ψ3 u3 ) − β2 S(λS − λI2 ) − εβ1 (I1 + I12 )(λI2 − λI12 )−
εβ2 I1 (λI1 − λI12 ).
Here f2 (I2 ) is convex increasing function and f2 (0) = 0, I2 > 0, we have
f2 (I2 ) − f˙2 (I2 )I2 ≤ 0 by Property 2. From (54) and similar to (71) we also
have that λ̇R (t) − λ̇I2 (t) ≤ 0 and it contradicts Property 1, then part III
of the lemma follows.
60
IV. Similar to previous parts we have proved that (λR (t) − λI12 (t)) >
0. In this case we assume that
(λR (t∗ ) − λI1 (t∗ )) ≥ 0, (λR (t∗ ) − λI2 (t∗ )) ≥ 0,
(λR (t∗ ) − λI12 (t∗ )) = 0, (λS (t∗ ) − λI1 (t∗ )) ≥ 0,
(77)
(λI1 (t∗ ) − λI12 (t∗ )) ≥ 0, (λI2 (t∗ ) − λI12 (t∗ )) ≥ 0.
Also we can see that
H + f1 (I1 ) + f2 (I2 ) + h1 (u1 (t)) + h2 (u2 (t)) + h3 (u3 )−
−β1 S(t)(I1 (t) + I12 (t))(λI1 − λS ) − β2 S(t)(I2 (t)+
+I12 (t))(λI2 − λS ) − εβ2 I1 (t)(I2 (t) + I12 (t))(λI12 − λI1 )−
(78)
−εβ1 I2 (t)(I1 (t) + I12 (t))(λI12 − λI2 ) − (1 − εβ2 )(λR −
−λI1 )u1 I1 − (1 − εβ1 )(λR − λI2 )u2 I2 ≤ −f3 (I12 )+
+(λR (T ) − λI12 (T ))u3 I12 (T ) ≤ 0.
According to Property 1, we have to check time moment t∗+ :
λ̇R (t∗+ ) − λ̇I12 (t∗+ ) = − dIdf123 − β1 S(λS − λI1 ) − β2 S(λS − λI2 )−
εβ2 I1 (λI1 − λI12 ) − εβ1 I2 (λI2 − λI12 ) +
1
I12 (H
+ f1 (I1 ) + f2 (I2 )
+f3 (I12 ) + h1 (u1 ) + h2 (u2 ) + h3 (u3 ) − β1 S(I1 + I12 )(λI1 − λS )−
β2 S(I2 + I12 )(λI2 − λS ) − εβ2 I1 (I2 + I12 )(λI12 − λI1 )−
εβ1 I2 (I1 + I12 )(λI12 − λI2 ) − (1 − εβ2 )(λR − λI1 )u1 I1 −
(1 − εβ1 )(λR − λI2 )u2 I2 ) =
(79)
61
=
1
˙
I12 (−f3 I12
+ f3 ) +
1
I12 (H
+ f1 (I1 ) + f2 (I2 )+
h1 (u1 ) + h2 (u2 ) + h3 (u3 ) − β1 S(I1 + I12 )(λI1 − λS )−
β2 S(I2 + I12 )(λI2 − λS ) − εβ2 I1 (I2 + I12 )(λI12 − λI1 )−
(80)
εβ1 I2 (I1 + I12 )(λI12 − λI2 ) − ψ1 u1 − ψ2 u2 − β1 S(λS − λI1 )−
β2 S(λS − λI2 ) − εβ2 I1 (λI1 − λI12 ) − εβ1 I2 (λI2 − λI12 ).
From (54), (71) and Property 2 we receive that (λ̇R (t) − λ̇I12 (t)) ≤
0 that contradicts Property 1 and hence part IV follows, then the time
moment t∗ does not exist. This is completed the proof of lemma 4.
Now we have to prove that (λR (t) − λS (t)) ≥ 0. Analogously we
assume:
(λR (t∗ ) − λI1 (t∗ )) ≥ 0, (λR (t∗ ) − λI2 (t∗ )) ≥ 0,
(λR (t∗ ) − λI12 (t∗ )) ≥ 0, (λR (t∗ ) − λS (t∗ )) = 0,
∗
∗
∗
∗
(81)
(λS (t ) − λI1 (t )) ≥ 0, (λS (t ) − λI2 (t )) ≥ 0
(λI1 (t∗ ) − λI12 (t∗ )) ≥ 0, (λI2 (t∗ ) − λI12 (t∗ )) ≥ 0.
From (53) we obtain:
λ̇R (t∗ ) − λ̇S (t∗ ) = −β1 (λS − λI1 )(I1 + I12 )−
(82)
−β2 (λS − λI2 )(I2 + I12 ) + γu4 (λR − λS ).
According to assumption (81), the difference λ̇R (t∗+ ) − λ̇S (t∗+ ) ≤ 0
is negative at the time moment t∗+ , then function (λR (t) − λS (t)) decrease.
It contradicts Property 1. At time moment t = T (λR (t) − λS (t)) is equal
to zero and on the whole interval it is decreasing. We proved that (λR (t) −
λS (t)) ≥ 0 for t ∈ [0; T ].
62
We get that on the interval t ∈ [0; T ]:
(λR (t) − λI1 (t)) ≥ 0;
(λR (t) − λI2 (t)) ≥ 0;
(λR (t) − λI12 (t)) ≥ 0;
(λR (t) − λS (t)) ≥ 0.
The proof of the main proposition is completed.
63
(83)
Отзывы:
Авторизуйтесь, чтобы оставить отзыв