Home Stochastic optimal and time-optimal control studies for additional food provided prey–predator systems involving Holling type III functional response
Article Open Access

Stochastic optimal and time-optimal control studies for additional food provided prey–predator systems involving Holling type III functional response

  • Daliparthi Bhanu Prakash and Dasu Krishna Kiran Vamsi EMAIL logo
Published/Copyright: February 28, 2023

Abstract

This article consists of a detailed and novel stochastic optimal control analysis of a coupled non-linear dynamical system. The state equations are modelled as an additional food-provided prey–predator system with Holling type III functional response for predator and intra-specific competition among predators. We first discuss the optimal control problem as a Lagrangian problem with a linear quadratic control. Second, we consider an optimal control problem in the time-optimal control setting. We initially establish the existence of optimal controls for both these problems and later characterize these optimal controls using the Stochastic maximum principle. Further numerical simulations are performed based on stochastic forward-backward sweep methods for realizing the theoretical findings. The results obtained in these optimal control problems are discussed in the context of biological conservation and pest management.

MSC 2010: 37A50; 60H10; 60J65; 60J70

1 Introduction

Millions of species coexist in this universe where the survival of each species is meticulously woven with respect to other species. The survival of one species (e.g. predator) is dependent on the existence of other species (e.g. prey). The very first mathematical models of such interactive species were given by Lotka [16] and Volterra [31]. In the last century, various mathematical models were proposed and their behaviours were studied.

One of the primary components in defining the prey–predator dynamics is the functional response. The functional response is the rate at which each predator captures prey [13]. A type III functional response is a sigmoidal response [12] that has predators foraging inefficiently at low prey densities. The Holling type III functional response is displayed by many organisms in nature [7,8,17,21]. Some studies involving additional food-provided prey–predator systems with Holling type III and type IV functional responses can be found in the study by Srinivasu et al. [28,29]. Recently, Ananth and Vamsi [13] have studied controllability of additional food systems with respect to quality and quantity of additional food as control variables for type III and type IV systems.

Also, in real world situations, some of the parameters involved in the model always fluctuate around some average value due to continuous variation in the environment. A large number of researchers introduced a stochastic environmental variation using the Brownian motion into parameters in the deterministic model to construct a stochastic model. Belabbas et al. [4] proved that a stochastic prey–predator model with a protection zone has a unique stationary distribution which is ergodic. Sengupta et al. [24] obtained stochastic permanence for a stochastic prey–predator system with Holling type III functional response. Guo and Ruan [10] dealt with a Holling type II stochastic prey–predator model with additional food that has an ergodic stationary distribution. Li and Zhao [14] studied the deterministic and stochastic dynamics of a modified Leslie-Gower prey–predator system with a simplified Holling type IV scheme. Qin et al. [20] dealt with the survival and ergodicity of a stochastic Holling type III prey–predator model with Markovian switching in an impulsive polluted environment. In recent times, optimal control theory is being applied on various stochastic models in order to achieve optimal control values, which minimize the cost. Also to our knowledge very limited research exists on stochastic optimal control theory for additional food- provided stochastic prey–predator systems involving different functional responses.

Motivated by the above discussions, in this article, we study two optimal control problems for an additional food-provided prey–predator system with Holling type III functional response for predator. We also consider the intra-specific competition among the predators to avoid their unbounded growth in the absence of target prey [28,30]. This system is a coupled non-linear dynamical system. The first optimal control problem is a Lagrange problem. This has applications for biological conservation of species where we find the optimal quality and quantity of additional food to be provided to predators to maximize the populations of predator and prey [6,9,15]. The second optimal control problem is a special kind of optimal control problem, known as time-optimal control problem, where we determine the optimal additional food to be provided to the system to reach the final state in minimum time. This has several applications to pest management [13,25,26].

The rest of the article is organized as follows. In Section 2, we formulate the stochastic model for an additional food-provided system involving Holling type III response. Section 3 deals with the discussion on corresponding linear quadratic optimal control problems with applications to biological conservation with reference to quality and quantity of additional food as stochastic control variables. Later in Section 4 we deal with the time optimal control problems for these systems with applications to both biological conservation and pest management. Numerical simulations are also done to validate the theoretical findings in Sections 3 and 4. Finally in Section 5, discussion and conclusions are given.

2 Stochastic model formulation

In this work, we consider the following deterministic prey–predator model with Holling type III functional response (i.e. c N 2 a 2 + N 2 ) and additional food for predator given by:

(1) d N ( t ) d t = r N ( t ) 1 N ( t ) K c N 2 ( t ) P ( t ) a 2 + N 2 ( t ) + α η A 2 , d P ( t ) d t = g N 2 ( t ) + η A 2 a 2 + N 2 ( t ) + α η A 2 P ( t ) m P ( t ) d P 2 ( t ) .

Here, the parameter g denotes the conversion efficiency that represents the rate at which prey biomass gets converted into predator biomass and can be obtained as a ratio of nutritive value of prey and the handling time of predator [30]. The term α denotes the ratio of the maximum growth rates of the predator when it consumes the prey and additional food, respectively. This term can be seen to be an equivalent of quality of additional food. The term η represents the ratio of the search rate of the predator for additional food and prey, respectively. The term d P 2 ( t ) accounts for the intra-specific competition among the predators to avoid their unbounded growth in the absence of target prey [28,30]. The biological meaning for all the parameters involved in system (1) is enlisted and described in Table 1. Further details regarding the derivation of the functional response, model formulation, and the parameters can be found in the study of Srinivasu et al. [30].

Table 1

Description of variables and parameters present in the systems (1), (2), and (3)

Parameter Definition Dimension
T Time Time
N Prey density Biomass
P Predator density Biomass
A Additional food Biomass
r Prey intrinsic growth rate Time 1
K Prey carrying capacity Biomass
c Rate of predation Time 1
a Half saturation value of the predators Biomass
g Conversion efficiency Time 1
m Death rate of predators in absence of prey Time 1
d Predator intra-specific competition Biomass 1 time 1
α Quality of additional food Dimensionless
ξ Quantity of additional food Biomass 2

To reduce the complexity in the analysis, we now reduce the number of parameters in model (1) by introducing the transformations N = a x , P = a y c . Then system (1) gets transformed to:

(2) d x ( t ) d t = r x ( t ) 1 x ( t ) γ x 2 ( t ) y ( t ) 1 + x 2 ( t ) + α ξ , d y ( t ) d t = g y ( t ) x 2 ( t ) + ξ 1 + x 2 ( t ) + α ξ m y ( t ) δ y 2 ( t ) ,

where γ = K a , ξ = η A a 2 , δ = d a c . Here the term ξ denotes the quantity of additional food perceptible to the predator with respect to the prey relative to the nutritional value of prey to the additional food. Hence, this can be seen to be an equivalent of quantity of additional food.

As in previous studies [4,24], we now suppose that the intrinsic growth rate of prey and the death rate of predator are mainly affected by environmental noise such that

r r + σ 1 d W 1 ( t ) , m m + σ 2 d W 2 ( t ) ,

where W i ( t ) ( i = 1 , 2 ) are the mutually independent standard Brownian motions with W i ( 0 ) = 0 and σ 1 and σ 2 are positive constants and they represent the intensities of the white noise. Hence, system (2) with the environmental noise for parameters r and m and the Holling type III predator functional response now gets modified to

(3) d x = r x ( t ) 1 x ( t ) γ x 2 ( t ) y ( t ) 1 + x 2 ( t ) + α ξ d t + σ 1 x ( t ) d W 1 ( t ) , d y = g y ( t ) x 2 ( t ) + ξ 1 + x 2 ( t ) + α ξ m y ( t ) δ y 2 ( t ) d t + σ 2 y ( t ) d W 2 ( t ) .

2.1 Existence of global positive solution

Theorem 1

For any initial value ( x 0 , y 0 ) R + 2 there exists a unique solution ( x ( t ) , y ( t ) ) of system (3) on t 0 and the solution will remain in R + 2 with probability 1.

The aforementioned theorem for the existence of solutions of (3) can be proved in similar lines to the proof in the study by Sengupta et al. [24] using the Lyapunov method.

3 Stochastic optimal control problem

In this section, we theoretically establish the existence of an optimal control for system (3) with both quality and quantity of additional food as stochastic controls, respectively. Later we use the stochastic maximum principle to find this optimal control. We then numerically simulate and depict the same with applications to biological conservation.

3.1 Quality of additional food as a stochastic optimal control

In this section, we wish to achieve biological conservation of maximizing prey and predator population for system (3) with quality of additional food as a control variable (for a fixed quantity of food) with minimum supply of the food.

To attain this we consider the following objective functional along with the state equation (3)

(4) J ( u ) = E 0 T A 1 x ( t ) A 2 y ( t ) + A 3 α 2 ( t ) 2 d t ,

where A 1 , A 2 , and A 3 are positive constants.

Here our goal is to find an optimal control α such that J ( α ) J ( α ) , α ( t ) U , where U is an admissible control set defined by U = { α ( t ) 0 α ( t ) α max t ( 0 , t f ] } , where α max R + .

Now we use the existence theorem in the study of Yong and Zhou [34] for establishing the existence of an optimal control for the stochastic control system (3). Drawing parallels with the theorem, we see that

X ( t ) = x ( t ) y ( t ) , b ( t , X ( t ) , u ( t ) ) = b 1 ( t ) b 2 ( t ) = r x ( t ) 1 x ( t ) γ x 2 ( t ) y ( t ) 1 + x 2 ( t ) + α ξ g y ( t ) x 2 ( t ) + ξ 1 + x 2 ( t ) + α ξ m y ( t ) δ y 2 ( t ) ,

σ ( t , X ( t ) , u ( t ) ) = σ 1 x ( t ) 0 0 σ 2 y ( t ) , d W ( t ) = d W 1 ( t ) d W 2 ( t ) .

We also note that

(5) b X = b 1 x b 1 y b 2 x b 2 y = r 1 2 x γ 2 x y ( 1 + α ξ ) ( 1 + x 2 + α ξ ) 2 x 2 1 + x 2 + α ξ 2 g x y ( 1 + ( α 1 ) ξ ) ( 1 + x 2 + α ξ ) 2 g x 2 + ξ 1 + x 2 + α ξ m 2 δ y ,

(6) σ 1 = σ 1 x 0 σ X 1 = x ( σ 1 x ) y ( σ 1 x ) x ( 0 ) y ( 0 ) = σ 1 0 0 0 , σ 2 = 0 σ 2 y σ X 2 = x ( 0 ) y ( 0 ) x ( σ 2 y ) y ( σ 2 y ) = 0 0 0 σ 2 .

Also we let, f ( t , X ( t ) , u ( t ) ) = A 1 x ( t ) A 2 y ( t ) + A 3 α 2 ( t ) 2 and h ( X ( T ) ) = 0 . Here, the term f ( t , X ( t ) , u ( t ) ) denotes the Lagrangian.

3.1.1 Existence of a stochastic optimal control w.r.t. quality ( α ) as control

Theorem 2

For any ( x , y ) R + 2 , if J ( u ) is finite, then the stochastic optimal control problems (3)–(4) admit an optimal control.

Proof

In order to prove the existence of an optimal control, it is enough to show only that the conditions [C1]–[C3] hold for systems (3)–(4) as the cost functional (4) at the optimal control is finite. This can be attributed to the bounded nature of the functions x ( t ) , y ( t ) , and α ( t ) .

  1. Since U = [ 0 , α max ] , U is closed and bounded. Hence, the set U is a compact metric space with a usual metric.

  2. Since the existence of positive global solutions of (3) is guaranteed in Theorem 1, the functions b , σ are Lipschitz continuous [24].

  3. From the study of Boyd et al. [5], a set C R n is convex if the line segment between any two points in C lies in C i.e. if any x 1 , x 2 C and any 0 θ 1 , we have θ x 1 + ( 1 θ ) x 2 C .

Now, we need to prove that for every ( t , x ) [ 0 , T ] × R n , the set ( b , σ σ T , f ) = { ( b 1 , b 2 , ( σ σ T ) i j , f ) α U , i , j = { 1 , 2 } } is convex in C R 7 .

Since

σ = σ 1 x ( t ) 0 0 σ 2 y ( t ) , σ T = σ 1 x ( t ) 0 0 σ 2 y ( t ) and σ σ T ( t , X ( t ) , u ( t ) ) = σ 1 2 x 2 ( t ) 0 0 σ 2 2 y 2 ( t ) ,

we see that σ σ T is independent of u, hence the set { ( σ σ T ) i j α U , i , j = { 1 , 2 } } is convex.

Now considering f = A 1 x ( t ) A 2 y ( t ) + A 3 α 2 ( t ) 2 , we let f 1 , f 2 be two such elements given by

f 1 = A 1 x ( t ) A 2 y ( t ) + A 3 α 1 2 ( t ) 2 and f 2 = A 1 x ( t ) A 2 y ( t ) + A 3 α 2 2 ( t ) 2 , α 1 , α 2 U .

Now let 0 θ 1 . Then,

θ f 1 + ( 1 θ ) f 2 = θ ( A 1 x ( t ) A 2 y ( t ) + A 3 α 1 2 ( t ) 2 ) + ( 1 θ ) A 1 x ( t ) A 2 y ( t ) + A 3 α 2 2 ( t ) 2 = A 1 x ( t ) A 2 y ( t ) + A 3 2 ( θ α 1 2 ( t ) + ( 1 θ ) α 2 2 ( t ) ) .

Since y = α 2 is a convex set, α 3 U θ α 1 2 ( t ) + ( 1 θ ) α 2 2 ( t ) = α 3 2 .

Therefore, θ f 1 + ( 1 θ ) f 2 = A 1 x ( t ) A 2 y ( t ) + A 3 α 3 2 ( t ) 2 . Hence, the set { ( f ) u U } is convex.

Now we are left to prove that { ( b 1 , b 2 ) u U } is convex.

Since b 1 = r x ( t ) 1 x ( t ) γ x 2 ( t ) y ( t ) 1 + x 2 ( t ) + α ξ , y ( t ) 1 + x 2 ( t ) + α ξ = r x 1 x γ b 1 x 2 and substituting y ( t ) 1 + x 2 ( t ) + α ξ = r x ( 1 x γ ) b 1 x 2 in b 2 = g y ( t ) x 2 ( t ) + ξ 1 + x 2 ( t ) + α ξ m y ( t ) δ y 2 ( t ) , we obtain

(7) g x 2 ( x 2 + ξ ) b 1 + b 2 = g ( x 2 + ξ ) r x 2 1 x γ m y δ y 2 .

Now we let D = g ( x 2 + ξ ) r x 2 1 x γ m y δ y 2 .

Hence, equation (7) can be seen as A b = D , where A = g ( x 2 + ξ ) x 2 1 , b = b 1 ( t ) b 2 ( t ) .

Let b 1 = b 1 1 ( t ) b 2 1 ( t ) and b 11 = b 1 11 ( t ) b 2 11 ( t ) be two distinct non-zero vectors.

Since A ( θ b 1 + ( 1 θ ) b 11 ) = θ ( A b 1 ) + ( 1 θ ) ( A b 11 ) = θ ( D ) + ( 1 θ ) D = D , we conclude that the set b = b 1 ( t ) b 2 ( t ) u U is convex.

Hence for every ( t , x ) [ 0 , T ] × R n , the set ( b , σ σ T , f ) = { ( b 1 , b 2 , ( σ σ T ) i j , f ) α U , i , j = { 1 , 2 } } is convex in C R 7 . Therefore, from Theorem 5.3 in the study of Yong and Zhou [34], the existence of optimal control is guaranteed.□

3.1.2 Characteristics of stochastic optimal control w.r.t. quality ( α ) as control

We now find the characteristics of stochastic optimal control using the stochastic maximum principle as in [34].

As the diffusion term σ in (6) is independent of the control, the solution of the second-order adjoint equations will not be helpful in calculating the optimal control values.

From the stochastic maximum principle as in the study of Yong and Zhou [34], we see that there exists a stochastic process given by

( p ( ) , q ( ) ) L F 2 ( 0 , T ; R n ) × ( L F 2 ( 0 , T ; R n ) ) m , p ( t ) = p 1 ( t ) p 2 ( t ) , q ( t ) = q 1 ( t ) q 3 ( t ) q 2 ( t ) q 4 ( t ) ,

where the vectors p ( t ) and q ( t ) satisfy the adjoint equations given by

d p 1 ( t ) d p 2 ( t ) = b 1 x b 2 x b 1 y b 2 y p 1 ( t ) p 2 ( t ) + σ X 1 ( t , X ¯ ( t ) , u ¯ ( t ) ) T q 1 ( t ) + σ X 2 ( t , X ¯ ( t ) , u ¯ ( t ) ) T q 2 ( t ) f x f y d t + q 1 ( t ) q 3 ( t ) q 2 ( t ) q 4 ( t ) d W 1 ( t ) d W 2 ( t ) p 1 ( T ) = 0 , p 2 ( T ) = 0 .

Substituting f x = A 1 , f y = A 2 and the values of b X , σ X 1 , σ X 2 from (5) to (6), we see that the adjoint equations for the optimal control α are given as:

d p 1 ( t ) = b 1 x p 1 ( t ) + b 2 x p 2 ( t ) + σ 1 q 1 + A 1 d t + q 1 ( t ) d W 1 ( t ) + q 3 ( t ) d W 2 ( t ) , d p 2 ( t ) = b 1 y p 1 ( t ) + b 2 y p 2 ( t ) + σ 2 q 4 + A 2 d t + q 2 ( t ) d W 1 ( t ) + q 4 ( t ) d W 2 ( t ) , p 1 ( T ) = 0 , p 2 ( T ) = 0

On further simplification we obtain

(8) d p 1 ( t ) = r 1 2 x γ 2 x y ( 1 + α ξ ) ( 1 + x 2 + α ξ ) 2 p 1 ( t ) + 2 g x y ( 1 + ( α 1 ) ξ ) ( 1 + x 2 + α ξ ) 2 p 2 ( t ) + σ 1 q 1 + A 1 d t + q 1 ( t ) d W 1 ( t ) + q 3 ( t ) d W 2 ( t ) , d p 2 ( t ) = x 2 1 + x 2 + α ξ p 1 ( t ) + g x 2 + ξ 1 + x 2 + α ξ m 2 δ y p 2 ( t ) + σ 2 q 4 + A 2 d t + q 2 ( t ) d W 1 ( t ) + q 4 ( t ) d W 2 ( t ) , p 1 ( T ) = 0 , p 2 ( T ) = 0 .

The solutions of the aforementioned equation (8) gives ( p 1 ( t ) , p 2 ( t ) ) , which are the co-state vectors.

Now from the study of Yong and Zhou [34], we see that the Hamiltonian for system (3) is given by:

H ( t , X , u , p , q ) = p , b + tr [ q T σ ] + A 1 x ( t ) + A 2 y ( t ) A 3 α 2 ( t ) 2 = p 1 ( t ) p 2 ( t ) b 1 ( t ) b 2 ( t ) + tr q 1 ( t ) q 2 ( t ) q 3 ( t ) q 4 ( t ) σ 1 x 0 0 σ 2 y + A 1 x ( t ) + A 2 y ( t ) A 3 α 2 ( t ) 2 = p 1 ( t ) b 1 + p 2 ( t ) b 2 + q 1 σ 1 x + q 4 σ 2 y + A 1 x ( t ) + A 2 y ( t ) A 3 α 2 ( t ) 2 = r x ( t ) 1 x ( t ) γ x 2 ( t ) y ( t ) 1 + x 2 ( t ) + α ξ p 1 ( t ) + g y ( t ) x 2 ( t ) + ξ 1 + x 2 ( t ) + α ξ m y ( t ) δ y 2 ( t ) p 2 ( t ) + q 1 σ 1 x + q 4 σ 2 y + A 1 x ( t ) + A 2 y ( t ) A 3 α 2 ( t ) 2 .

Now from the Hamiltonian maximization condition as in the study of Yong and Zhou [34], we have

( t , X ¯ ( t ) , α ¯ ( t ) ) = max α U H ( t , X ¯ ( t ) , α ( t ) ) H α = 0 α r x ( t ) 1 x ( t ) γ x 2 ( t ) y ( t ) 1 + x 2 ( t ) + α ξ p 1 + α g y ( t ) x 2 ( t ) + ξ 1 + x 2 ( t ) + α ξ m y ( t ) δ y 2 ( t ) p 2 ( t ) + α q 1 σ 1 x + q 4 σ 2 y + A 1 x ( t ) + A 2 y ( t ) A 3 α 2 ( t ) 2 = 0 ξ x 2 y p 1 ( 1 + x 2 + α ξ ) 2 ξ g y p 2 ( x 2 + ξ ) ( 1 + x 2 + α ξ ) 2 A 3 α = 0 ( A 3 ξ 2 ) α 3 + ( 2 A 3 ξ ( 1 + x 2 ) ) α 2 + ( A 3 ( 1 + x 2 ) 2 ) α + ( ξ g y p 2 ( x 2 + ξ ) ξ x 2 y p 1 ) = 0 .

Now from Descartes’ rule of signs, we see that the above cubic equation admits a positive α only if g p 2 ( x 2 + ξ ) x 2 p 1 < 0 .

On solving the aforementioned cubic equation, we see that the optimal quality control is given by

(9) α = ( 1 + x 2 c 0 ξ ) 2 3 c 0 ξ 2 ,

where c 0 = 10 ( 1 + x 2 ) 3 ξ 3 + 27 y c 1 2 A 3 ξ + 1 2 108 c 1 y ( 1 + x 2 ) 2 A 3 ξ 4 + 729 y 2 c 1 2 A 3 2 ξ 2 1 3 and c 1 = g p 2 ( x 2 + ξ ) x 2 y

3.2 Quantity of additional food as a stochastic optimal control

In this section, we wish to achieve biological conservation of maximizing prey and predator population for system (3) with quantity of additional food as a control variable (for a fixed quality of food) with minimum supply of the food.

To attain this we consider the following objective functional with the state equation (3)

(10) J ( u ) = E 0 T A 1 x ( t ) A 2 y ( t ) + A 3 ξ 2 ( t ) 2 d t ,

where A 1 , A 2 , and A 3 are positive constants.

Here our goal is to find an optimal control ξ such that J ( ξ ) J ( ξ ) , ξ ( t ) U where U is an admissible control set defined by U = { ξ ( t ) 0 ξ ( t ) ξ max t ( 0 , t f ] } where ξ max R + .

Now comparing the cost functional (10) with cost functional in the study of Yong and Zhou [34], we see that

f ( t , X ( t ) , u ( t ) ) = A 1 x ( t ) A 2 y ( t ) + A 3 ξ 2 ( t ) 2 and h ( X ( T ) ) = 0 .

3.2.1 Existence and characteristics of stochastic optimal control w.r.t. quantity ( ξ ) as control

Theorem 3

For any ( x , y ) R + 2 , if J ( u ) is finite, then the stochastic optimal control problems (3) and (10)admit an optimal control.

Proof

The above theorem for the existence of optimal control for systems (3) and (10) can be proved in similar lines to the proof of Theorem 2. The cost functional (10) at optimal control is also finite owing to the finite Lagrangian in (10).□

Since f x and f y are independent of control parameter ξ , the adjoint equations are same as (8) in the previous subsection.

Hence from the stochastic maximum principle as in the study of Yong and Zhou [34], there exist a stochastic process given by

( p ( ) , q ( ) ) L F 2 ( 0 , T ; R n ) × ( L F 2 ( 0 , T ; R n ) ) m , p ( t ) = p 1 ( t ) p 2 ( t ) , q ( t ) = q 1 ( t ) q 3 ( t ) q 2 ( t ) q 4 ( t )

(11) d p 1 ( t ) = r 1 2 x γ 2 x y ( 1 + α ξ ) ( 1 + x 2 + α ξ ) 2 p 1 ( t ) + 2 g x y ( 1 + ( α 1 ) ξ ) ( 1 + x 2 + α ξ ) 2 p 2 ( t ) + σ 1 q 1 + A 1 d t + q 1 ( t ) d W 1 ( t ) + q 3 ( t ) d W 2 ( t ) d p 2 ( t ) = x 2 1 + x 2 + α ξ p 1 ( t ) + g x 2 + ξ 1 + x 2 + α ξ m 2 δ y p 2 ( t ) + σ 2 q 4 + A 2 d t + q 2 ( t ) d W 1 ( t ) + q 4 ( t ) d W 2 ( t ) p 1 ( T ) = 0 , p 2 ( T ) = 0 .

The solutions of equation (11) give ( p 1 ( t ) , p 2 ( t ) ) , which are the co-state vectors.

Now from the study of Yong and Zhou [34], we see that the Hamiltonian for system (3) is given by:

H ( t , X , u , p , q ) = p , b + tr [ q T σ ] + A 1 x ( t ) + A 2 y ( t ) A 3 ξ 2 ( t ) 2 = p 1 ( t ) p 2 ( t ) b 1 ( t ) b 2 ( t ) + tr q 1 ( t ) q 2 ( t ) q 3 ( t ) q 4 ( t ) σ 1 x 0 0 σ 2 y + A 1 x ( t ) + A 2 y ( t ) A 3 ξ 2 ( t ) 2 = p 1 ( t ) b 1 + p 2 ( t ) b 2 + q 1 σ 1 x + q 4 σ 2 y + A 1 x ( t ) + A 2 y ( t ) A 3 ξ 2 ( t ) 2 = r x ( t ) 1 x ( t ) γ x 2 ( t ) y ( t ) 1 + x 2 ( t ) + α ξ p 1 ( t ) + g y ( t ) x 2 ( t ) + ξ 1 + x 2 ( t ) + α ξ m y ( t ) δ y 2 ( t ) p 2 ( t ) + q 1 σ 1 x + q 4 σ 2 y + A 1 x ( t ) + A 2 y ( t ) A 3 ξ 2 ( t ) 2 .

Now from the Hamiltonian maximization condition as in the study of Yong and Zhou [34], we have

( t , X ( t ) , ξ ( t ) ) = max α U H ( t , X ( t ) , ξ ( t ) ) H ξ = 0

ξ r x ( t ) 1 x ( t ) γ x 2 ( t ) y ( t ) 1 + x 2 ( t ) + α ξ p 1 + ξ g y ( t ) x 2 ( t ) + ξ 1 + x 2 ( t ) + α ξ m y ( t ) δ y 2 ( t ) p 2 ( t ) + ξ q 1 σ 1 x + q 4 σ 2 y + A 1 x ( t ) + A 2 y ( t ) A 3 ξ 2 ( t ) 2 = 0 α x 2 y p 1 ( 1 + x 2 + α ξ ) 2 + g y p 2 ( 1 + ( 1 α ) x 2 ) ( 1 + x 2 + α ξ ) 2 A 3 ξ = 0 ( A 3 α 2 ) ξ 3 + ( 2 A 3 α ( 1 + x 2 ) ) ξ 2 + ( A 3 ( 1 + x 2 ) 2 ) ξ ( x 2 y p 1 α + g y p 2 ( 1 + ( 1 α ) x 2 ) ) = 0 .

Now from Descartes’ rule of signs, we see that the above cubic equation admits a positive ξ only if x 2 y p 1 α + g y p 2 ( 1 + ( 1 α ) x 2 ) > 0 .

On solving the above cubic equation, we see that the optimal quantity control is given by

(12) ξ = ( 1 + x 2 + c 1 α ) 2 3 c 1 α 2 ,

where c 1 = 10 ( 1 + x 2 ) 3 α 3 27 c 2 2 A 3 α 2 + 1 2 108 ( 1 + x 2 ) 3 c 2 A 3 α 5 + 729 c 2 2 A 3 2 α 4 1 3 and c 2 = x 2 y p 1 α + g y p 2 ( 1 + ( 1 α ) x 2 ) .

3.3 Numerical simulations

In this section, we numerically illustrate the theoretical findings of the aforementioned sections with application to biological conservation.

Using the Taylor series expansion, the optimal control problems are simulated and plotted using the Stochastic Forward and Backward Sampling approach. The state equation (3) and the adjoint equations (8) and (14) are solved using the forward and backward processes, respectively. The forward process is simulated using the Euler-Maruyama scheme [18]. Among the various methods available to discretize the backward process, we chose an implicit scheme with a back propagation of the conditional expectations, which is of order 1/2 [35]. These methods are implemented in Python using Sympy, Numpy, and Matplotlib packages.

The sub plots present in Figures 1 and 2 give the optimal state trajectories, optimal co-state trajectories, phase diagram, and optimal control trajectories, respectively. These examples reiterate the importance of additional food as a control variable in the context of ecological conservation.

Figure 1 
                  The optimal trajectory of the optimal control problem discussed in Section 3.1 with the objective to drive the system from the initial state 
                        
                           
                           
                              
                                 (
                                 
                                    5
                                    ,
                                    1
                                 
                                 )
                              
                           
                           \left(5,1)
                        
                      to the terminal state 
                        
                           
                           
                              
                                 (
                                 
                                    7
                                    ,
                                    0.25
                                 
                                 )
                              
                           
                           \left(7,0.25)
                        
                     . The parameter values for these plots are chosen to be 
                        
                           
                           
                              r
                              =
                              1
                           
                           r=1
                        
                     , 
                        
                           
                           
                              γ
                              =
                              7
                           
                           \gamma =7
                        
                     , 
                        
                           
                           
                              g
                              =
                              0.4
                           
                           g=0.4
                        
                     , 
                        
                           
                           
                              m
                              =
                              0.37
                           
                           m=0.37
                        
                     , 
                        
                           
                           
                              δ
                              =
                              0.1
                           
                           \delta =0.1
                        
                     , 
                        
                           
                           
                              ξ
                              =
                              0.1
                           
                           \xi =0.1
                        
                     , 
                        
                           
                           
                              
                                 
                                    A
                                 
                                 
                                    1
                                 
                              
                              =
                              
                                 
                                    A
                                 
                                 
                                    2
                                 
                              
                              =
                              
                                 
                                    A
                                 
                                 
                                    3
                                 
                              
                              =
                              1
                           
                           {A}_{1}={A}_{2}={A}_{3}=1
                        
                     . The terminal value of co-state variables is 
                        
                           
                           
                              
                                 (
                                 
                                    
                                       
                                          p
                                       
                                       
                                          1
                                       
                                    
                                    
                                       (
                                       
                                          T
                                       
                                       )
                                    
                                    ,
                                    
                                       
                                          p
                                       
                                       
                                          2
                                       
                                    
                                    
                                       (
                                       
                                          T
                                       
                                       )
                                    
                                 
                                 )
                              
                              =
                              
                                 (
                                 
                                    0
                                    ,
                                    0
                                 
                                 )
                              
                           
                           \left({p}_{1}\left(T),{p}_{2}\left(T))=\left(0,0)
                        
                     . The intensities of noise are chosen as 
                        
                           
                           
                              
                                 
                                    σ
                                 
                                 
                                    1
                                 
                              
                              =
                              
                                 
                                    σ
                                 
                                 
                                    2
                                 
                              
                              =
                              0.02
                           
                           {\sigma }_{1}={\sigma }_{2}=0.02
                        
                     . In order to account the randomness, the stochastic control problem is simulated for 5,000 times and the average is plotted in these figures. This example depicts the optimal quality of additional food required to achieve biological conservation, for a fixed quantity of food. The optimal control plot shows that the predators have to be given high quality of additional food initially and then the lower quality of food in order to conserve both the species from extinction. The desired terminal state is reached in 
                        
                           
                           
                              T
                              =
                              87
                           
                           T=87
                        
                      units of time.
Figure 1

The optimal trajectory of the optimal control problem discussed in Section 3.1 with the objective to drive the system from the initial state ( 5 , 1 ) to the terminal state ( 7 , 0.25 ) . The parameter values for these plots are chosen to be r = 1 , γ = 7 , g = 0.4 , m = 0.37 , δ = 0.1 , ξ = 0.1 , A 1 = A 2 = A 3 = 1 . The terminal value of co-state variables is ( p 1 ( T ) , p 2 ( T ) ) = ( 0 , 0 ) . The intensities of noise are chosen as σ 1 = σ 2 = 0.02 . In order to account the randomness, the stochastic control problem is simulated for 5,000 times and the average is plotted in these figures. This example depicts the optimal quality of additional food required to achieve biological conservation, for a fixed quantity of food. The optimal control plot shows that the predators have to be given high quality of additional food initially and then the lower quality of food in order to conserve both the species from extinction. The desired terminal state is reached in T = 87 units of time.

Figure 2 
                  The optimal trajectory of the optimal control problem discussed in Section 3.2 with the objective to drive the system from the initial state 
                        
                           
                           
                              
                                 (
                                 
                                    5
                                    ,
                                    1
                                 
                                 )
                              
                           
                           \left(5,1)
                        
                      to the terminal state 
                        
                           
                           
                              
                                 (
                                 
                                    7
                                    ,
                                    0.25
                                 
                                 )
                              
                           
                           \left(7,0.25)
                        
                     . The parameter values for these plots are chosen to be 
                        
                           
                           
                              r
                              =
                              1
                           
                           r=1
                        
                     , 
                        
                           
                           
                              γ
                              =
                              7
                           
                           \gamma =7
                        
                     , 
                        
                           
                           
                              g
                              =
                              0.4
                           
                           g=0.4
                        
                     , 
                        
                           
                           
                              m
                              =
                              0.37
                           
                           m=0.37
                        
                     , 
                        
                           
                           
                              δ
                              =
                              0.1
                           
                           \delta =0.1
                        
                     , 
                        
                           
                           
                              α
                              =
                              2
                           
                           \alpha =2
                        
                     , 
                        
                           
                           
                              
                                 
                                    A
                                 
                                 
                                    1
                                 
                              
                              =
                              
                                 
                                    A
                                 
                                 
                                    2
                                 
                              
                              =
                              
                                 
                                    A
                                 
                                 
                                    3
                                 
                              
                              =
                              1
                           
                           {A}_{1}={A}_{2}={A}_{3}=1
                        
                     . The terminal value of co-state variables is 
                        
                           
                           
                              
                                 (
                                 
                                    
                                       
                                          p
                                       
                                       
                                          1
                                       
                                    
                                    
                                       (
                                       
                                          T
                                       
                                       )
                                    
                                    ,
                                    
                                       
                                          p
                                       
                                       
                                          2
                                       
                                    
                                    
                                       (
                                       
                                          T
                                       
                                       )
                                    
                                 
                                 )
                              
                              =
                              
                                 (
                                 
                                    0
                                    ,
                                    0
                                 
                                 )
                              
                           
                           \left({p}_{1}\left(T),{p}_{2}\left(T))=\left(0,0)
                        
                     . The intensities of noise are chosen as 
                        
                           
                           
                              
                                 
                                    σ
                                 
                                 
                                    1
                                 
                              
                              =
                              
                                 
                                    σ
                                 
                                 
                                    2
                                 
                              
                              =
                              0.02
                           
                           {\sigma }_{1}={\sigma }_{2}=0.02
                        
                     . In order to account the randomness, the stochastic control problem is simulated for 5,000 times and the average is plotted in these figures. This example depicts the optimal quantity of additional food required to achieve biological conservation, for a fixed quality of food. The optimal control plot shows that the predators have to be given high quantity of additional food initially and then the lower quantity of food in order to conserve both the species from extinction. The desired terminal state is reached in 
                        
                           
                           
                              T
                              =
                              84
                           
                           T=84
                        
                      units of time.
Figure 2

The optimal trajectory of the optimal control problem discussed in Section 3.2 with the objective to drive the system from the initial state ( 5 , 1 ) to the terminal state ( 7 , 0.25 ) . The parameter values for these plots are chosen to be r = 1 , γ = 7 , g = 0.4 , m = 0.37 , δ = 0.1 , α = 2 , A 1 = A 2 = A 3 = 1 . The terminal value of co-state variables is ( p 1 ( T ) , p 2 ( T ) ) = ( 0 , 0 ) . The intensities of noise are chosen as σ 1 = σ 2 = 0.02 . In order to account the randomness, the stochastic control problem is simulated for 5,000 times and the average is plotted in these figures. This example depicts the optimal quantity of additional food required to achieve biological conservation, for a fixed quality of food. The optimal control plot shows that the predators have to be given high quantity of additional food initially and then the lower quantity of food in order to conserve both the species from extinction. The desired terminal state is reached in T = 84 units of time.

4 Stochastic time-optimal control problem

In this section, we wish to achieve pest eradication of nearly prey-elimination stage for system (3) with quality and quantity of additional food as control variables in minimum time.

To attain this goal, we consider the following time optimal control problem with the following objective functional along with the state equations (3)

(13) J ( u ) = E 0 T 1 d t .

Here our goal is to find optimal controls α and ξ such that J ( α , ξ ) J ( α , ξ ) , ( α ( t ) , ξ ( t ) ) U , where U is an admissible control set defined by

U = { ( α ( t ) , ξ ( t ) ) 0 α ( t ) α max , 0 ξ ( t ) ξ max t ( 0 , t f ] } , where ( α max , ξ max ) ( R + ) 2 .

Now comparing the cost functional (13) with cost functional in [34], we see that f ( t , X ( t ) , u ( t ) ) = 1 and h ( X ( T ) ) = 0 .

4.1 Existence and characteristics of stochastic optimal controls w.r.t. quality ( α ) and quantity ( ξ ) as controls

Theorem 4

For any ( x , y ) R + 2 , the stochastic optimal control problem (3), (13) admits an optimal control.

Proof

The above theorem for the existence of optimal control for systems (3) and (13) can be proved in similar lines to the proof of Theorem 2. The cost functional (13) at optimal control is also finite owing to the finite Lagrangian in (13).□

Hence from the stochastic maximum principle as in the study of Yong and Zhou [34], there exist a stochastic processes given by

( p ( ) , q ( ) ) L F 2 ( 0 , T ; R n ) × ( L F 2 ( 0 , T ; R n ) ) m , p ( t ) = p 1 ( t ) p 2 ( t ) , q ( t ) = q 1 ( t ) q 3 ( t ) q 2 ( t ) q 4 ( t )

d p 1 ( t ) d p 2 ( t ) = b 1 x b 2 x b 1 y b 2 y p 1 ( t ) p 2 ( t ) + σ X 1 ( t , X ¯ ( t ) , u ¯ ( t ) ) T q 1 ( t ) + σ X 2 ( t , X ¯ ( t ) , u ¯ ( t ) ) T q 2 ( t ) f x f y d t + q 1 ( t ) q 3 ( t ) q 2 ( t ) q 4 ( t ) d W 1 ( t ) d W 2 ( t ) p 1 ( T ) = 0 , p 2 ( T ) = 0 .

Substituting f x = 0 , f y = 0 and the values of b X , σ X 1 , σ X 2 from (5), (6), we see that the adjoint equations for the optimal controls ( α , ξ ) are given as:

d p 1 ( t ) = b 1 x p 1 ( t ) + b 2 x p 2 ( t ) + σ 1 q 1 d t + q 1 ( t ) d W 1 ( t ) + q 3 ( t ) d W 2 ( t ) , d p 2 ( t ) = b 1 y p 1 ( t ) + b 2 y p 2 ( t ) + σ 2 q 4 d t + q 2 ( t ) d W 1 ( t ) + q 4 ( t ) d W 2 ( t ) , p 1 ( T ) = 0 , p 2 ( T ) = 0 .

On further simplification, we see that

(14) d p 1 ( t ) = r 1 2 x γ 2 x y ( 1 + α ξ ) ( 1 + x 2 + α ξ ) 2 p 1 ( t ) + 2 g x y ( 1 + ( α 1 ) ξ ) ( 1 + x 2 + α ξ ) 2 p 2 ( t ) + σ 1 q 1 d t + q 1 ( t ) d W 1 ( t ) + q 3 ( t ) d W 2 ( t ) , d p 2 ( t ) = x 2 1 + x 2 + α ξ p 1 ( t ) + g x 2 + ξ 1 + x 2 + α ξ m 2 δ y p 2 ( t ) + σ 2 q 4 d t + q 2 ( t ) d W 1 ( t ) + q 4 ( t ) d W 2 ( t ) , p 1 ( T ) = 0 , p 2 ( T ) = 0 .

The solutions of the above equations (8) give ( p 1 ( t ) , p 2 ( t ) ) , which are the co-state vectors.

Now from the study of Yong and Zhou [34], we see that the Hamiltonian for system (3) is given by:

H ( t , X , u , p , q ) = p , b + tr [ q T σ ] 1 = p 1 ( t ) p 2 ( t ) b 1 ( t ) b 2 ( t ) + tr q 1 ( t ) q 2 ( t ) q 3 ( t ) q 4 ( t ) σ 1 x 0 0 σ 2 y 1 = p 1 ( t ) b 1 + p 2 ( t ) b 2 + q 1 σ 1 x + q 4 σ 2 y 1 = r x ( t ) 1 x ( t ) γ x 2 ( t ) y ( t ) 1 + x 2 ( t ) + α ξ p 1 ( t ) + g y ( t ) x 2 ( t ) + ξ 1 + x 2 ( t ) + α ξ m y ( t ) δ y 2 ( t ) p 2 ( t ) + q 1 σ 1 x + q 4 σ 2 y 1 .

Now from the Hamiltonian maximization condition, we have

( t , X ( t ) , α ( t ) , ξ ( t ) ) = max ( α , ξ ) U H ( t , X ( t ) , α ( t ) , ξ ( t ) ) H α = 0 , H ξ = 0 .

Considering

H α = 0 α r x ( t ) 1 x ( t ) γ x 2 ( t ) y ( t ) 1 + x 2 ( t ) + α ξ p 1 + α g y ( t ) x 2 ( t ) + ξ 1 + x 2 ( t ) + α ξ m y ( t ) δ y 2 ( t ) p 2 ( t ) + α ( q 1 σ 1 x + q 4 σ 2 y 1 ) = 0 ξ x 2 y p 1 ( 1 + x 2 + α ξ ) 2 ξ g y p 2 ( x 2 + ξ ) ( 1 + x 2 + α ξ ) 2 = 0 ξ g y p 2 ( x 2 + ξ ) ξ x 2 y p 1 = 0 ξ = p 1 g p 2 g p 2 .

Similarly considering

H ξ = 0 ξ r x ( t ) 1 x ( t ) γ x 2 ( t ) y ( t ) 1 + x 2 ( t ) + α ξ p 1 + ξ g y ( t ) x 2 ( t ) + ξ 1 + x 2 ( t ) + α ξ m y ( t ) δ y 2 ( t ) p 2 ( t ) + α ( q 1 σ 1 x + q 4 σ 2 y 1 ) = 0 g y p 2 ( 1 + x 2 ) + α ( x 2 y p 1 g y p 2 x 2 ) = 0 α = g p 2 ( 1 + x 2 ) x 2 ( g p 2 p 1 ) .

Hence, the optimal quality and quantity variables for the stochastic time optimal control problem are given by

( α , ξ ) = g p 2 ( 1 + x 2 ) x 2 ( g p 2 p 1 ) , p 1 g p 2 g p 2 .

4.2 Numerical simulations

In this section, we numerically illustrate the theoretical findings of the above time optimal control problem with application to both biological conservation and pest management.

Using the Taylor series expansion, the time optimal control problem is simulated and plotted using the Stochastic Forward and Backward Sampling approach. The state equations (3) and the adjoint equations (8) and (14) are solved using the forward and backward processes, respectively. The forward process is simulated using the Euler-Maruyama scheme [18]. Among the various methods available to discretize the backward process, we chose an implicit scheme with a back propagation of the conditional expectations, which is of order 1/2 [35]. These methods are implemented in Python using Sympy, Numpy, and Matplotlib packages.

The subplots present in Figures 3 and 4 give the optimal state trajectories, optimal co-state trajectories, phase diagram, and optimal control trajectories, respectively. These examples reiterate the importance of additional food as control variables in the context of both ecological conservation and pest management.

Figure 3 
                  The optimal trajectory of the time-optimal control problem discussed in Section 4 with the objective to drive the system from the initial state 
                        
                           
                           
                              
                                 (
                                 
                                    5
                                    ,
                                    1
                                 
                                 )
                              
                           
                           \left(5,1)
                        
                      to the terminal state 
                        
                           
                           
                              
                                 (
                                 
                                    6.75
                                    ,
                                    0.3
                                 
                                 )
                              
                           
                           \left(6.75,0.3)
                        
                      in minimum time. The parameter values for these plots are chosen to be 
                        
                           
                           
                              r
                              =
                              1
                           
                           r=1
                        
                     , 
                        
                           
                           
                              γ
                              =
                              7
                           
                           \gamma =7
                        
                     , 
                        
                           
                           
                              g
                              =
                              0.4
                           
                           g=0.4
                        
                     , 
                        
                           
                           
                              m
                              =
                              0.37
                           
                           m=0.37
                        
                     , 
                        
                           
                           
                              δ
                              =
                              0.1
                           
                           \delta =0.1
                        
                     , 
                        
                           
                           
                              α
                              =
                              2
                           
                           \alpha =2
                        
                     , 
                        
                           
                           
                              ξ
                              =
                              0.5
                           
                           \xi =0.5
                        
                     . The initial value of co-state variables is 
                        
                           
                           
                              
                                 (
                                 
                                    
                                       
                                          p
                                       
                                       
                                          1
                                       
                                    
                                    
                                       (
                                       
                                          0
                                       
                                       )
                                    
                                    ,
                                    
                                       
                                          p
                                       
                                       
                                          2
                                       
                                    
                                    
                                       (
                                       
                                          0
                                       
                                       )
                                    
                                 
                                 )
                              
                              =
                              
                                 (
                                 
                                    0
                                    ,
                                    0
                                 
                                 )
                              
                           
                           \left({p}_{1}\left(0),{p}_{2}\left(0))=\left(0,0)
                        
                     . The intensities of noise are chosen as 
                        
                           
                           
                              
                                 
                                    σ
                                 
                                 
                                    1
                                 
                              
                              =
                              
                                 
                                    σ
                                 
                                 
                                    2
                                 
                              
                              =
                              0.02
                           
                           {\sigma }_{1}={\sigma }_{2}=0.02
                        
                     . In order to account the randomness, the stochastic time-optimal control problem is simulated for 5,000 times and the average is plotted in these figures. This example depicts the optimal quality and quantity of additional food required to achieve biological conservation. The optimal control plot shows that the predators have to be given constant quantity of additional food (1 unit) and the constant quantity of additional food (1 unit) in order to conserve both the species from extinction. The desired terminal state is reached in 
                        
                           
                           
                              T
                              =
                              63
                           
                           T=63
                        
                      units of time.
Figure 3

The optimal trajectory of the time-optimal control problem discussed in Section 4 with the objective to drive the system from the initial state ( 5 , 1 ) to the terminal state ( 6.75 , 0.3 ) in minimum time. The parameter values for these plots are chosen to be r = 1 , γ = 7 , g = 0.4 , m = 0.37 , δ = 0.1 , α = 2 , ξ = 0.5 . The initial value of co-state variables is ( p 1 ( 0 ) , p 2 ( 0 ) ) = ( 0 , 0 ) . The intensities of noise are chosen as σ 1 = σ 2 = 0.02 . In order to account the randomness, the stochastic time-optimal control problem is simulated for 5,000 times and the average is plotted in these figures. This example depicts the optimal quality and quantity of additional food required to achieve biological conservation. The optimal control plot shows that the predators have to be given constant quantity of additional food (1 unit) and the constant quantity of additional food (1 unit) in order to conserve both the species from extinction. The desired terminal state is reached in T = 63 units of time.

Figure 4 
                  The optimal trajectory of the time-optimal control problem discussed in Section 4 with the objective to drive the system from the initial state 
                        
                           
                           
                              
                                 (
                                 
                                    4.5
                                    ,
                                    5
                                 
                                 )
                              
                           
                           \left(4.5,5)
                        
                      to the terminal state 
                        
                           
                           
                              
                                 (
                                 
                                    0.26
                                    ,
                                    7.41
                                 
                                 )
                              
                           
                           \left(0.26,7.41)
                        
                      in minimum time. The parameter values for these plots are chosen to be 
                        
                           
                           
                              r
                              =
                              0.7
                           
                           r=0.7
                        
                     , 
                        
                           
                           
                              γ
                              =
                              5
                           
                           \gamma =5
                        
                     , 
                        
                           
                           
                              g
                              =
                              1.2
                           
                           g=1.2
                        
                     , 
                        
                           
                           
                              m
                              =
                              0.19
                           
                           m=0.19
                        
                     , 
                        
                           
                           
                              δ
                              =
                              0.05
                           
                           \delta =0.05
                        
                     , 
                        
                           
                           
                              α
                              =
                              1.4
                           
                           \alpha =1.4
                        
                     , 
                        
                           
                           
                              ξ
                              =
                              1.4
                           
                           \xi =1.4
                        
                     . The initial value of co-state variables is 
                        
                           
                           
                              
                                 (
                                 
                                    
                                       
                                          p
                                       
                                       
                                          1
                                       
                                    
                                    
                                       (
                                       
                                          0
                                       
                                       )
                                    
                                    ,
                                    
                                       
                                          p
                                       
                                       
                                          2
                                       
                                    
                                    
                                       (
                                       
                                          0
                                       
                                       )
                                    
                                 
                                 )
                              
                              =
                              
                                 (
                                 
                                    9
                                    ,
                                    −
                                    6
                                 
                                 )
                              
                           
                           \left({p}_{1}\left(0),{p}_{2}\left(0))=\left(9,-6)
                        
                     . The intensities of noise are chosen as 
                        
                           
                           
                              
                                 
                                    σ
                                 
                                 
                                    1
                                 
                              
                              =
                              
                                 
                                    σ
                                 
                                 
                                    2
                                 
                              
                              =
                              0.02
                           
                           {\sigma }_{1}={\sigma }_{2}=0.02
                        
                     . In order to account the randomness, the stochastic time-optimal control problem is simulated for 5,000 times and the average is plotted in these figures. This example depicts the optimal quality and optimal quantity of additional food required to achieve pest management. The optimal control plot shows that the predators have to be given a constant quality of additional food (0.67 units) and a high quantity of additional food throughout for reducing the pest population effectively. This example suggests that even with low quality of additional food and high quantity, the prey (pest) can be reduced to a threshold level at which they no longer cause significant damage to the ecosystem. The desired terminal state is reached in 
                        
                           
                           
                              T
                              =
                              100
                           
                           T=100
                        
                      units of time.
Figure 4

The optimal trajectory of the time-optimal control problem discussed in Section 4 with the objective to drive the system from the initial state ( 4.5 , 5 ) to the terminal state ( 0.26 , 7.41 ) in minimum time. The parameter values for these plots are chosen to be r = 0.7 , γ = 5 , g = 1.2 , m = 0.19 , δ = 0.05 , α = 1.4 , ξ = 1.4 . The initial value of co-state variables is ( p 1 ( 0 ) , p 2 ( 0 ) ) = ( 9 , 6 ) . The intensities of noise are chosen as σ 1 = σ 2 = 0.02 . In order to account the randomness, the stochastic time-optimal control problem is simulated for 5,000 times and the average is plotted in these figures. This example depicts the optimal quality and optimal quantity of additional food required to achieve pest management. The optimal control plot shows that the predators have to be given a constant quality of additional food (0.67 units) and a high quantity of additional food throughout for reducing the pest population effectively. This example suggests that even with low quality of additional food and high quantity, the prey (pest) can be reduced to a threshold level at which they no longer cause significant damage to the ecosystem. The desired terminal state is reached in T = 100 units of time.

5 Discussion and conclusions

The provision of additional food has proven to be very effective in conserving endangered species [11,19,22] as well as controlling invasive or harmful species [23,32,33]. A significant amount of theoretical work was also done on the impact of additional food in various ecosystems [25,27]. Both theoretically and experimentally, the technique of providing additional food for biocontrol seemed to be very effective. Nowadays researchers are also focusing on the optimal additional food to be given to the predators. Some of the findings in this direction show that the quality and quantity of additional food provided play a crucial role in these optimal studies [30]. For instance, previous studies [1,3] worked on the deterministic Holling type III predator-prey systems and found the optimal quality and quantity of additional food required to drive the system to a desired terminal state.

It is a well-known fact that the real-life systems behave more chaotic. Stochastic setting can be an ideal tool to capture such random system dynamics better than that of deterministic system dynamics. Motivated by the above discussions, in this work, we studied a stochastic predator-prey system with additional food for predator exhibiting Holling type III functional response which in a way can be considered to be a generalization of the work [1,3]. We have also incorporated the intra-specific competition among predators in the model to make the model more realistic.

In this work, we initially considered a Lagrangian optimal control problem with a cost that is linear w.r.t. state and quadratic w.r.t. control with the end goal of biological conservation of both the species. We considered two cases with the quality of additional food and the quantity of additional food as control variables, respectively. We initially established the existence of optimal controls and later characterized these optimal control values using the stochastic maximum principle. We then numerically simulated using the Forward-Backward Sampling method. These results showed that biological conservation of the system can be achieved by the optimal control strategy of providing the quality and the quantity of additional food high initially and further reducing them over the time.

Second, we studied a time-optimal control problem with the end goal of biological conservation of both the species and also to achieve the goal of pest management in minimum time. In this problem, we worked with a multi-dimensional control involving both the quality and quantity of additional food as control variables. We initially established the existence of optimal controls and later characterized these optimal control values using the stochastic maximum principle. We also plotted these solutions for two different sets of parameters, one applied to biological conservation and the other to pest management. In case of biological conservation, both the controls did not exhibit any switch over the time. In case of pest management, optimal quality of the additional food remained low and constant throughout. The optimal quantity of additional food fluctuated initially and remained high throughout.

The present stochastic optimal control studies can further be improved by incorporating alley effect and even multiple prey and predator species to make the model more realistic. Also, the model presented here does not account for time delay. As diffusion term in this system is independent of control, the results obtained will be similar to that of deterministic case. It will be very interesting to take up problems with control in the diffusion. Since the current systems are of higher orders of non-linearity, it turns out that the numerical methods play a crucial role in understanding the chaotic behaviours. Hence, the study of higher-order numerical methods such as stochastic Runge-Kutta methods will enhance the output. Also not much work has been done where multiple noise is added simultaneously. We aim to take up these studies in the future.

Acknowledgments

The authors dedicate this article to the founder chancellor of SSSIHL, Bhagawan Sri Sathya Sai Baba. The corresponding author also dedicates this article to his loving elder brother D. A. C. Prakash who still lives in his heart.

  1. Funding information: This research was supported by National Board of Higher Mathematics (NBHM), Government of India (GoI) under project grant – Time Optimal Control and Bifurcation Analysis of Coupled Nonlinear Dynamical Systems with Applications to Pest Management, Sanction number: (02011/11/2021NBHM(R.P)/R&D II/10074).

  2. Author contributions: Daliparthi Bhanu Prakash conceptualized the idea, performed the analysis, simulations and prepared draft of the manuscript. Dasu Krishna Kiran Vamsi conceptualized the flow for the entire draft of the manuscript, verified all the results and edited the draft. He is also responsible for corresponding the manuscript.

  3. Conflict of interest: The authors have no conflict of interest to disclose.

  4. Ethical approval: This research did not require ethical approval.

  5. Data availability statement: Data sharing is not applicable to this article as no data sets were generated or analysed during the current study.

References

[1] Ananth, V. S., & Vamsi, D. K. K. (2021a). Influence of quantity of additional food in achieving biological conservation and pest management in minimum-time for prey–predator systems involving Holling type III response. Heliyon, 7(8), e07699. 10.1016/j.heliyon.2021.e07699Search in Google Scholar PubMed PubMed Central

[2] Ananth, V. S., & Vamsi, D. K. K. (2021b). An optimal control study with quantity of additional food as control in prey–predator systems involving inhibitory effect. Computational and Mathematical Biophysics, 9(1), 114–145. 10.1515/cmb-2020-0121Search in Google Scholar

[3] Ananth, V. S., & Vamsi, D. K. K. (2022). Achieving minimum-time biological conservation and pest management for additional food provided predator-prey systems involving inhibitory effect: A qualitative investigation. Acta Biotheoretica, 70(1), 1–51. 10.1007/s10441-021-09430-2Search in Google Scholar PubMed

[4] Belabbas, M., Ouahab, A., & Souna, F. (2021). Rich dynamics in a stochastic predator-prey model with protection zone for the prey and multiplicative noise applied on both species. Nonlinear Dynamics, 106(3), 2761–2780. 10.1007/s11071-021-06903-4Search in Google Scholar

[5] Boyd, S., Boyd, S. P., & Vandenberghe, L. (2004). Convex optimization. Cambridge, England: Cambridge University Press. 10.1017/CBO9780511804441Search in Google Scholar

[6] El-Gohary, A., & Bukhari, F. A. (2003). Optimal control of stochastic prey–predator models. Applied Mathematics and Computation, 146(2–3), 403–415. 10.1016/S0096-3003(02)00592-1Search in Google Scholar

[7] Elkinton, J. S., Liebhold, A. M., & Muzika, R.-M. (2004). Effects of alternative prey on predation by small mammals on gypsy moth pupae. Population Ecology, 46(2), 171–178. 10.1007/s10144-004-0175-ySearch in Google Scholar

[8] Fernández-Arhex, V., & Corley, J. C. (2005). The functional response of ibalia leucospoides (hymenoptera: Ibaliidae), a parasitoid of sirex noctilio (hymenoptera: Siricidae). Biocontrol Science and Technology, 15(2), 207–212. 10.1080/09583150400016084Search in Google Scholar

[9] Gu, X., & Zhu, W. (2016). Stochastic optimal control of predator-prey ecosystem by using stochastic maximum principle. Nonlinear Dynamics, 85(2), 1177–1184. 10.1007/s11071-016-2752-ySearch in Google Scholar

[10] Guo, X., & Ruan, D. (2020). Extinction and ergodic stationary distribution of a markovian-switching prey–predator model with additional food for predator. Mathematical Modelling of Natural Phenomena, 15, 46. 10.1051/mmnp/2019055Search in Google Scholar

[11] Harwood, J. D., Sunderland, K. D., & Symondson, W. O. (2004). Prey selection by linyphiid spiders: molecular tracking of the effects of alternative prey on rates of aphid consumption in the field. Molecular Ecology, 13(11), 3549–3560. 10.1111/j.1365-294X.2004.02331.xSearch in Google Scholar PubMed

[12] Holling, C. S. (1966). The functional response of invertebrate predators to prey density. The Memoirs of the Entomological Society of Canada, 98(S48), 5–86. 10.4039/entm9848fvSearch in Google Scholar

[13] Kot, M. (2001). Elements of mathematical ecology. Cambridge, England: Cambridge University Press. 10.1017/CBO9780511608520Search in Google Scholar

[14] Li, L., & Zhao, W. (2021). Deterministic and stochastic dynamics of a modified Leslie-Gower prey–predator system with simplified Holling-type iv scheme. Mathematical Biosciences and Engineering, 18(3), 2813–2831. 10.3934/mbe.2021143Search in Google Scholar PubMed

[15] Liu, M. (2015). Optimal harvesting policy of a stochastic predator-prey model with time delay. Applied Mathematics Letters, 48, 102–108. 10.1016/j.aml.2014.10.007Search in Google Scholar

[16] Lotka, A. J. (1925). Elements of physical biology. United States: Williams & Wilkins. Search in Google Scholar

[17] Morozov, A. Y. (2010). Emergence of Holling type iii zooplankton functional response: bringing together field evidence and mathematical modelling. Journal of Theoretical Biology, 265(1), 45–54. 10.1016/j.jtbi.2010.04.016Search in Google Scholar PubMed

[18] Platen, E., & Bruti-Liberati, N. (2010). Numerical solution of stochastic differential equations with jumps in finance (Vol. 64), Germany: Springer Science & Business Media. 10.1007/978-3-642-13694-8Search in Google Scholar

[19] Putman, R., & Staines, B. W. (2004). Supplementary winter feeding of wild red deer cervus elaphus in Europe and North America: justifications, feeding practice and effectiveness. Mammal Review, 34(4), 285–306. 10.1111/j.1365-2907.2004.00044.xSearch in Google Scholar

[20] Qin, W., Zhang, H., & He, Q. (2021). Survival and ergodicity of a stochastic Holling-iii predator-prey model with Markovian switching in an impulsive polluted environment. Advances in Difference Equations, 2021(1), 1–19. 10.1186/s13662-021-03238-7Search in Google Scholar

[21] Redpath, S. M., & Thirgood, S. J. (1999). Numerical and functional responses in generalist predators: hen harriers and peregrines on Scottish grouse moors. Journal of Animal Ecology, 68(5), 879–892. 10.1046/j.1365-2656.1999.00340.xSearch in Google Scholar

[22] Redpath, S. M., Thirgood, S. J., & Leckie, F. M. (2001). Does supplementary feeding reduce predation of red grouse by hen harriers? Journal of Applied Ecology, 38(6), 1157–1168. 10.1046/j.0021-8901.2001.00683.xSearch in Google Scholar

[23] Sabelis, M. W., & Van Rijn, P. C. (2006). When does alternative food promote biological pest control? IOBC WPRS Bulletin, 29(4), 195. Search in Google Scholar

[24] Sengupta, S., Das, P., & Mukherjee, D. (2018). Stochastic non-autonomous Holling type-iii prey–predator model with predators intra-specific competition. Discrete & Continuous Dynamical Systems-B, 23(8), 3275. 10.3934/dcdsb.2018244Search in Google Scholar

[25] Srinivasu, P. D. N., & Prasad, B. S. R. V. (2010). Time optimal control of an additional food provided predator-prey system with applications to pest management and biological conservation. Journal of Mathematical Biology, 60(4), 591–613. 10.1007/s00285-009-0279-2Search in Google Scholar PubMed

[26] Srinivasu, P. D. N., & Prasad, B. S. R. V. (2011). Role of quantity of additional food to predators as a control in predator-prey systems with relevance to pest management and biological conservation. Bulletin of Mathematical Biology, 73(10), 2249–2276. 10.1007/s11538-010-9601-9Search in Google Scholar PubMed

[27] Srinivasu, P. D. N., Prasad, B. S. R. V., & Venkatesulu, M. (2007). Biological control through provision of additional food to predators: a theoretical study. Theoretical Population Biology, 72(1), 111–120. 10.1016/j.tpb.2007.03.011Search in Google Scholar PubMed

[28] Srinivasu, P. D. N., Vamsi, D. K. K., & Aditya, I. (2018a). Biological conservation of living systems by providing additional food supplements in the presence of inhibitory effect: a theoretical study using predator-prey models. Differential Equations and Dynamical Systems, 26(1), 213–246. 10.1007/s12591-016-0344-4Search in Google Scholar

[29] Srinivasu, P. D. N., Vamsi, D. K. K., & Ananth, V. S. (2018b). Additional food supplements as a tool for biological conservation of predator-prey systems involving type iii functional response: A qualitative and quantitative investigation. Journal of Theoretical Biology, 455, 303–318. Search in Google Scholar

[30] Srinivasu, P. D. N., Vamsi, D. K. K., & Ananth, V. S. (2018c). Additional food supplements as a tool for biological conservation of predator-prey systems involving type iii functional response: A qualitative and quantitative investigation. Journal of Theoretical Biology, 455, 303–318. 10.1016/j.jtbi.2018.07.019Search in Google Scholar PubMed

[31] Volterra, V. (1927). Variazioni e fluttuazioni del numero daindividui in specie animali conviventi. (Vol. 2). Societá anonima tipografica “Leonardo da Vinci”.Search in Google Scholar

[32] Wade, M. R., Zalucki, M. P., Wratten, S. D., & Robinson, K. A. (2008). Conservation biological control of arthropods using artificial food sprays: current status and future challenges. Biological Control, 45(2), 185–199. 10.1016/j.biocontrol.2007.10.024Search in Google Scholar

[33] Winkler, K., Wäckers, F. L., Stingli, A., & Van Lenteren, J. C. (2005). Plutella xylostella (diamondback moth) and its parasitoid diadegma semiclausum show different gustatory and longevity responses to a range of nectar and honeydew sugars. Entomologia Experimentalis et Applicata, 115(1), 187–192. 10.1111/j.1570-7458.2005.00254.xSearch in Google Scholar

[34] Yong, J., & Zhou, X. Y. (1999). Stochastic controls: Hamiltonian systems and HJB equations (Vol. 43). Springer Science & Business Media. Search in Google Scholar

[35] Zhang, J. (2004). A numerical scheme for BSDEs. The Annals of Applied Probability, 14(1), 459–488. 10.1214/aoap/1075828058Search in Google Scholar

Received: 2022-11-15
Revised: 2023-01-22
Accepted: 2023-01-28
Published Online: 2023-02-28

© 2023 the author(s), published by De Gruyter

This work is licensed under the Creative Commons Attribution 4.0 International License.

Articles in the same Issue

  1. Special Issue: Infectious Disease Modeling In the Era of Post COVID-19
  2. A comprehensive and detailed within-host modeling study involving crucial biomarkers and optimal drug regimen for type I Lepra reaction: A deterministic approach
  3. Application of dynamic mode decomposition and compatible window-wise dynamic mode decomposition in deciphering COVID-19 dynamics of India
  4. Role of ecotourism in conserving forest biomass: A mathematical model
  5. Impact of cross border reverse migration in Delhi–UP region of India during COVID-19 lockdown
  6. Cost-effective optimal control analysis of a COVID-19 transmission model incorporating community awareness and waning immunity
  7. Evaluating early pandemic response through length-of-stay analysis of case logs and epidemiological modeling: A case study of Singapore in early 2020
  8. Special Issue: Application of differential equations to the biological systems
  9. An eco-epidemiological model with predator switching behavior
  10. A numerical method for MHD Stokes model with applications in blood flow
  11. Dynamics of an eco-epidemic model with Allee effect in prey and disease in predator
  12. Optimal lock-down intensity: A stochastic pandemic control approach of path integral
  13. Bifurcation analysis of HIV infection model with cell-to-cell transmission and non-cytolytic cure
  14. Special Issue: Differential Equations and Control Problems - Part I
  15. Study of nanolayer on red blood cells as drug carrier in an artery with stenosis
  16. Influence of incubation delays on COVID-19 transmission in diabetic and non-diabetic populations – an endemic prevalence case
  17. Complex dynamics of a four-species food-web model: An analysis through Beddington-DeAngelis functional response in the presence of additional food
  18. A study of qualitative correlations between crucial bio-markers and the optimal drug regimen of Type I lepra reaction: A deterministic approach
  19. Regular Articles
  20. Stochastic optimal and time-optimal control studies for additional food provided prey–predator systems involving Holling type III functional response
  21. Stability analysis of an SIR model with alert class modified saturated incidence rate and Holling functional type-II treatment
  22. An SEIR model with modified saturated incidence rate and Holling type II treatment function
  23. Dynamic analysis of delayed vaccination process along with impact of retrial queues
  24. A mathematical model to study the spread of COVID-19 and its control in India
  25. Within-host models of dengue virus transmission with immune response
  26. A mathematical analysis of the impact of maternally derived immunity and double-dose vaccination on the spread and control of measles
  27. Influence of distinct social contexts of long-term care facilities on the dynamics of spread of COVID-19 under predefine epidemiological scenarios
Downloaded on 26.9.2025 from https://www.degruyterbrill.com/document/doi/10.1515/cmb-2022-0144/html
Scroll to top button