This paper develops a theoretical framework to investigate optimal harvesting control for stochastic delay differential systems. We first propose a novel stochastic two-predator and one-prey competitive system subject to time delays and Lévy jumps. Then we obtain sufficient conditions for persistence in mean and extinction of three species by using the stochastic qualitative analysis method. Finally, the optimal harvesting effort and the maximum of expectation of sustainable yield (ESY) are derived from Hessian matrix method and optimal harvesting theory of delay differential equations. Moreover, some numerical simulations are given to illustrate the theoretical results.

1. Introduction

Optimal control problem in the field of biological mathematics has been widely concerned by researchers. Resource exploitation always aims to maximum sustainable yield (MSY) or the profit associated with the maximum economic yield (MEY) [1]. How to obtain MSY or MEY involves the optimal harvesting control problem. Therefore, it is interesting and meaningful to investigate optimal harvesting control strategies for biological population models, especially stochastic population models. A number of researchers have investigated single-species or two-species population models [24]. However, the above two classes of models can not describe some natural phenomena completely and it is believed that models with three or more species can explain the dynamical behaviors of the population accurately [58]. Predator-prey models are arguably known as the most fundamental building blocks of the biosystems and ecosystems as the biomasses are grown out of their resource, which have been widely investigated [917]. As we all know, after a predator catches and feeds on a prey, the number of the predators will not increase at once, which needs the processes of digestion and absorption. Therefore, the time delays were considered in many differential systems [1821], especially stochastic delay models [2229].

The basic model we consider is based on the delay Lotka-Volterra model with two competitive predators and one prey. We propose the following model by assuming that the random perturbations of intrinsic growth rate are subjected to Gaussian white noise [30]:with initial valuewhere denotes the number of the prey species at time and (=1, 2) denote the predator species at time . stands for the growth rate of prey. and stand for the death rates of the two predators, respectively. is the interspecific competition rate of species , , and and are the competition rates of two predators. and are positive representing capture rates, and and are negative representing the growth rates from prey relatively. signifies the time delay, and are continuous functions on . (i=1,2,3) are mutually independent Brownian motion with defined on a complete probability space with a filtration that satisfies that it is right continuous and increasing with that contains all -null sets.

Various harvesting models have been used to investigate the optimal harvesting policies of renewable resources (e.g., Beddington and May [30], Neubert [31]). Recently, many authors explored the optimal harvesting models [3237]. According to the catch-per-unit-effort (CPUE) hypothesis, we consider that predators and are subject to exploitation with harvesting effort rates and , respectively.

However, in some cases, models just perturbed by the Gaussian white noise can not effectively and efficiently describe the circumstance when the species suffer sudden catastrophic disturbance in nature. The sudden environmental change can affect the dynamical behavior of the species significantly. Therefore, it is necessary to use the discontinuous stochastic process (e.g., Lévy jump) to model the abrupt nature phenomenon in ecosystem [3843].

To introduce the Lévy jump into the underlying stochastic model (1), we first give some facts about the Lévy jump [33]. Generally, a Lévy process can be decomposed into the sum of a linear drift, a Brownian motion , and a superposition of centered Poisson processes with different jump sizes which is the rate of arrival of the Poisson process with jump of size . According to the Lévy decomposition theorem [44], we know that where , , is a compensated Poisson process, and is a Poisson measure with characteristic measure on a measurable subset of with . The distribution of a Lévy process has the property of infinite divisibility and is characterized by its characteristic function , which is given by the following Lévy-Khintchine formula [45]: where is the indicator function of set and is the imaginary unit. The distribution of Lévy jump can be completely parameterized by .

Motivated by the above discussion, we can assume that the intrinsic growth rate and the death rates and of model (1) perturbed by the Lévy jump to signify the sudden climate change, [42, 46], and then we can obtain the following stochastic model incorporating Lévy jump:here , , , , and with initial value (2).

In this paper, we devote our main attention to obtain the optimal harvesting control strategy of system (5). To this end, we firstly investigate the dynamical behavior of the three species including persistence in mean and extinction and asymptotically stable distribution. Then we explore how the time delay, sudden environmental shock expressed by Lévy jump, and other factors affect the optimal harvesting policy and the maximum expectation of sustainable yield.

This paper is organized as follows. We discuss the persistence in mean and extinction of the three species in Section 2. Based on the conclusion of Section 2, we consider the optimal harvesting policy in Section 3. Finally, we conclude our results by numerical simulations and discussions in the last section.

2. Persistence in Mean and Extinction

For convenience, we give some notations which will be used for analyzing our main results. Let and let stand for all continuous functions from to , where is a 3-dimensional Euclidean space and is the positive cone in . Correspondingly, . Denote , , and . To make it more convenient, let , , , , , , , , and . Additionally, denote , , and , and . Accordingly, we define is the complement minor of in the deterministic , . Now we give a fundamental assumption of Lévy jump.

Assumption 1. Assume that , and there exists a constant such that

Assumption 1 implies that the intensity of Lévy jump can not be sufficiently large.

Lemma 2. For any given initial value , model (5) has a unique global positive solution on .

Proof. Firstly, we consider the following stochastic model:with initial value (2). It is easy to see that the coefficients of model (8) satisfy the local Lipschitz condition; therefore model (8) has a unique local solution on , where is the explosion time. According to Itô’s formula, we can find that is the unique positive local solution to model (5). Now let us prove . Thus, we introduce the following auxiliary model:with initial valueTaking advantage of the comparison theorem for stochastic equation [47] yields that, for ,According to Theorem 4.2 in Jiang and Shi [48], the explicit solution of model (10) iswhere , , and , [33]. It is not difficult to see that , , and are existent on ; thereby . This completes the proof of Lemma 2.

Lemma 3 (see [22]). Suppose that .
(i) If there exist three constants , , and , such thatfor all , where and are constants; then(ii) If there exist three constants , , and , such thatfor all ; then a.s.

Lemma 4. For model (10), consider the following:
(a) If , then(b) If , , and , then(c) If , , and , then(d) If , , and , then(e) If , , and , then

Proof. Firstly, we prove . Applying Itô’s formula to model (10), we can get thatIt follows from (22) that By Assumption 1, the quadratic variation of is where is Meyer’s angle bracket process. Utilizing the strong law of large numbers for local martingales gives that and . Since ,Substituting (27) into (23) gives thatwhere is small enough satisfying . Applying (i) in Lemma 3 gives Similarly, we can get thatSecondly, we prove (b). Since , applying Lemma 3 to (22) yields thatWe know thatUsing (31) in (32) gives thatMultiplying (32) and (33) by and , respectively, and, adding them, we can derive thatAn application of (31) gives thatConsequently, utilizing (35), (36), (37), and Lemma 3 yields that Similarly, we can get when , Thirdly, we prove (c). If , according to (35), (36), (37), and Lemma 3, one can obtain that Similarly, if , we can obtain The proofs of (d) and (e) are similar to that of (b) and (c), respectively, and hence are omitted.

Lemma 5. The solution of model (5) obeys

Proof. Since , , , we only need to proveAccording to the proof of Lemma 4, with the same method of [22], we can get (43). Therefore, (42) is obtained.

Theorem 6. For system (5), if , , , and , we set , , and , and then ; moreover,(i) if , then prey species and predator species go to extinction a.s.; i.e.,(ii) if , then two predators go to extinction a.s., and prey is persistent in mean; i.e.,(iii) if , then predator species goes to extinction, while prey species and predator species are persistent in mean; i.e.,(iv) if , then three species are persistent in mean; i.e.,

Proof. It is easy to prove . Applying Itô’s formula to model (5) yields thatFirstly, we prove (i). Since , and are positive, we can getNote that ; then . By Lemma 3, we getSubstituting this identity into (49), we can observe that, for sufficiently large ,where is small enough such that . By Lemma 3, we can get Similarly, we haveSecondly, we prove (ii). Since , we know Simplifying the above inequality gives that , which means . Furthermore, from , we have . According to (c) in Lemma 3 and (12), one can observe that Substituting the above identity into (48) and using Lemma 3, we can getThirdly, we prove (iii). Dividing (48), (49), and (50) by t, we can derive the following equations:Denote , as the solution of the following equations: Consequently, By (12), (43), and Lemma 5, we haveIn addition, According to Lemma 5, for arbitrarily given , there exists a such that when Multiplying (59), (60), and (61) by , , and 1, respectively, and adding them, one can observe that, for sufficiently large such that ,Using (43) in (81) yields thatSince , we can choose to be sufficiently small such that . Making use of the arbitrariness and Lemma 3 gives that Consequently, model (5) reduces to the following model:For system (70), similarly to the proof of Theorem 1 in [49], the following identities can be derived:Fourthly, we prove (iv). By (68), since , we know from the arbitrariness of and Lemma 3 thatDenote as the solution of the following equations: Then we have According to Lemma 5, for arbitrarily given , there exists a , such that Multiplying (59), (60), and (61) by , 1, and , respectively, and, adding them, we can obtain that, for ,Using (43) in (76) yields that Note that . According to the arbitrariness of and Lemma 3, we haveIt follows that, for any sufficiently small , there exist and such thatSubstituting (79) into (59) results in that, for sufficiently large ,