Four Eulerian network models are implemented to model high altitude air traffic flow. Three of the models use the framework of discrete time dynamical systems, while the fourth consists of a network of partial differential equations. The construction of these models is done using one year of air traffic data. The four models are applied to high altitude traffic for six Air Route Traffic Control Centers in the National Airspace System and surrounding airspace. Simulations are carried out for a full day of data for each of the models, to assess their predictive capabilities. The models’ predictions are compared to the recorded flight data. Several error metrics are used to characterize the relative accuracy of the models. The efficiency of the respective models is also compared in terms of computational time and memory requirements for the scenarios of interest. Control strategies are designed and implemented on similar benchmark scenarios for two of the models. They use techniques such as adjoint-based optimization, as well as mixed integer linear programming. A discussion of the four models’ structural differences explains why one model may outperform another.
1.
Introduction
The predator-prey model is one of the hotspots in biomathematics. For example, Yavuz and Sene [1] considered a fractional predator-prey model with harvesting rate, Chatterjee and Pal [2] studied a predator-prey model for the optimal control of fish harvesting through the imposition of a tax and Ghosh et al. [3] presented a three-component model consisting of one prey and two predator species using imprecise biological parameters as interval numbers and applied a functional parametric form in the proposed prey-predator system. Because of its important role in the ecosystem, the food chain model has been extensively studied [4,5,6,7,8]. Specifically, the classical four-species food chain model can be expressed as follows:
where x1(t), x2(t), x3(t) and x4(t) represent the densities of prey, primary predator, intermediate predator and top predator at time t, respectively. r1 is the growth rate of prey, r2, r3 and r4 are the death rates of primary predator, intermediate predator and top predator, respectively. aij and aji (i<j) are the capture rates and food conversion rates, respectively. aii are the intraspecific competition rates of species i. All parameters in system (1.1) are positive constants.
In ecology, biology, physics, engineering and other areas of applied sciences, continuous-time models, fractional-order models as well as discrete-time models have been widely adopted [9,10]. However, "time delays occur so often that to ignore them is to ignore reality" [11,12], and in the models of population dynamics, the delay differential equations are much more realistic [13,14,15]. We know that systems with discrete time delays and those with continuously distributed time delays do not contain each other but systems with S-type distributed time delays contain both. Introducing S-type distributed time delays into system (1.1) yields
where Dji(xi)(t)=ajixi(t)+∫0−τjixi(t+θ)dμji(θ), ∫0−τjixi(t+θ)dμji(θ) are Lebesgue-Stieltjes integrals, τji>0 are time delays, μji(θ) are nondecreasing bounded variation functions defined on [−τ,0], τ=max{τji}.
On the other hand, the deterministic system has its limitation in mathematical modeling of ecosystems since the parameters involved in the system are unable to capture the influence of environmental noises [16,17]. Introducing Gaussian white noises into the corresponding deterministic model is one common way to characterize environmental noises [18,19,20,21,22,23,24,25]. As we all know, Gaussian white noise ξ(t) is a stationary and ergodic stochastic process with ⟨ξ(t)⟩=0 and ⟨ξ(t)ξ(s)⟩=σ2δ(t−s), where σ2 is the noise intensity [26]. The readers can refer to [27,28,29,30,31,32,33,34] for more related works. In this paper, we assume that ri are affected by Gaussian white noises, i.e., r1↪r1+σ1˙W1(t), −r2↪−r2+σ2˙W2(t), −r3↪−r3+σ3˙W3(t) and −r4↪−r4+σ4˙W4(t). Then, system (1.2) becomes
where Wi(t) are mutually independent standard Wiener processes defined on a complete probability space (Ω,F,P) satisfying the usual statistical properties, namely ⟨dWi(t)⟩=0 and ⟨dWi(t)dWj(s)⟩=δijδ(t−s)dt [35].
Besides, population system may be affected by telephone noises which can cause the system to switch from one environmental regime to another [36,37,38]. So, telephone noises should be taken into consideration in system (1.3), resulting the following model:
where ρ(t) is a continuous time Markov chain with finite state space S={1,2,...,S}, which describes the telephone noises.
Moreover, the behaviour of real biological species, in different ecosystems, is affected by Lévy noises [39]. Lévy processes are characterized by stationary independent increments [40]. Assume that L(t) (t≥0) is a Lévy process, using the decomposition [41]
one can observe that the probability distribution of L(t) is infinitely divisible. The most general expression for the characteristic function of L(t) is
where sgn(k) is the sign function with
where α∈(0,2] is the stability parameter, σ is the scale parameter, σα is the noise intensity, μ∈R is the location parameter and β∈[−1,1] is the skewness parameter [39]. In addition, Lévy noises are statistically independent with zero mean. Now, let us further improve system (1.4) by considering Lévy noises. Some scholars pointed out that Lévy noises can be used to describe some sudden environmental perturbations, for instance, earthquakes and hurricanes [42,43,44,45,46,47]. In the context of an epidemic situation, random jumps could refer to sudden and significant increases in the number of cases or spread of the disease that occur unpredictably [48]. System (1.4) with Lévy noises can be expressed as follows:
where Si(t,ρ(t))=σi(ρ(t))dWi(t)+∫Zγi(μ,ρ(t))˜N(dt,dμ), N is a Poisson counting measure with characteristic measure λ on a measurable subset Z of [0,+∞), where λ(Z)<+∞ and ˜N(dt,dμ)=N(dt,dμ)−λ(dμ)dt, γj(μ,ρ(t))>−1 (μ∈Z) are bounded functions (j=1,2,3,4).
Finally, environmental pollution caused by agriculture, industries and other human activities has become a big challenge that is commonly concerned by international society. For example, with the rapid development of industrial and agricultural production, some chemical plants and other industries often periodically discharge sewage or other pollutants into rivers, soil and air [49]. These pollutants can cause direct damage to ecosystems, such as species extinction, desertification and the greenhouse effect. Hence, we extend system (1.5) into the following form:
where Δxi(t)=xi(t+)−xi(t), ΔCi0(t)=Ci0(t+)−Ci0(t) and ΔCe(t)=Ce(t+)−Ce(t). For other parameters in system (1.6), see Table 1.
To the best of our knowledge to date, results about a stochastic hybrid delay four-species food chain model with jumps have not been reported. So, in this paper we investigate the dynamics of a stochastic hybrid delay four-species food chain model with jumps in an impulsive polluted environment. The organization of this paper is as follows: In Section 2, some basic preliminaries are presented. In Section 3, the sufficient and necessary conditions for stochastic persistence in mean and extinction of each species are obtained. In Section 4, some numerical examples are provided to illustrate our main results. Finally, we conclude the paper with a brief conclusion and discussion in Section 5.
2.
Preliminaries
We have four fundamental assumptions for system (1.6).
Assumption 1. W1(t), W2(t), W3(t), W4(t), ρ(t) and N are mutually independent. ρ(t), taking values in S={1,2,...,S}, is irreducible with one unique stationary distribution π=(π1,π2,...,πS)T.
Assumption 2. rj(i)>0, ajk>0 and there exist γ∗j(i)≥γj∗(i)>−1 such that γj∗(i)≤γj(μ,i)≤γ∗j(i) (μ∈Z), ∀i∈S, j,k=1,2,3,4.
Remark 1. Assumption 2 implies that the intensities of Lévy jumps are not too big to ensure that the solution will not explode in finite time.
Assumption 3. 0<ki≤gi+mi (i=1,2,3,4), 0<b≤1−e−hγ.
Remark 2. Assumption 3 means 0≤Ci0(t)<1 and 0≤Ce(t)<1, which must be satisfied to be realistic because Ci0(t) and Ce(t) are concentrations of the toxicant (i=1,2,3,4).
Assumption 4. A22A33A44|A||Ξ|>A12A21A23A32A44|Ξ|+A23A32A34A43|A|2.
Lemma 1. [50,51] Ci0(t) involved in system (1.6) satisfies
3.
Persistence in mean and extinction
Denote
Denote Σ(2)=(Σ1,Σ2)T, Σ(3)=(Σ1,Σ2,Σ3)T, Σ=(Σ1,Σ2,Σ3,Σ4)T. Denote Aj is A with column j replaced by Σ(2) (j=1,2); Ξj is Ξ with column j replaced by Σ(3) (j=1,2,3); Δj is Δ with column j replaced by Σ (j=1,2,3,4).
Theorem 1. For any initial condition ϕ∈C([−τ,0],R4+), system (1.6) has a unique global solution (x1(t),x2(t),x3(t),x4(t))T∈R4+ on t∈R+ a.s. Moreover, for any constant p>0, there exists Ki(p)>0 such that supt∈R+E[xpi(t)]≤Ki(p) (i=1,2,3,4).
Proof. The proof is rather standard and hence is omitted (see e.g., [52]). □
Lemma 2. [53] Suppose Z(t)∈C(Ω×[0,+∞),R+) and limt→+∞o(t)t=0.
(i) If there exists constant δ0>0 such that for t≫1,
then
(ii) If there exist constants δ>0 and δ0>0 such that for t≫1,
then
Lemma 3. If |Δ4|>0, then |Δj|>0 (j=1,2,3).
Proof. Compute
Noting that Σj<0 (j=2,3,4), we obtain the desired assertion. □
First, let us consider the following auxiliary system:
Lemma 4. System (3.1) satisfies Table 2, where
Proof. Consider the following stochastic hybrid delay logistic model with Lévy jump in an impulsive polluted environment:
Thanks to Lemma 1 and Lemma 2.3 in [54], system (3.2) satisfies
By Itô's formula, we compute
where
Case(1): B1<0. Based on Eq (3.3), limt→+∞X1(t)=0 a.s. Therefore, for ∀ϵ∈(0,1) and t≫1,
Since Σ2<0, then limt→+∞X2(t)=0 a.s. Similarly, limt→+∞Xj(t)=0 a.s. (j=3,4).
Case(2): B1≥0. Then,
Consider the following auxiliary function:
Then X2(t)≤~X2(t) a.s. By Itô's formula, we get
In view of Lemma 2, we obtain:
If B1≥0, B2<0, then limt→+∞~X2(t)=0 a.s.
If B1≥0, B2≥0, then
Therefore, for arbitrary ζ>0, we have
According to system (3.4) and Eq (3.6), we obtain
Thanks to Lemma 2, we deduce:
If B1≥0, B2<0, then limt→+∞Xj(t)=0 a.s. (j=2,3,4).
If B1≥0, B2≥0, then
Case(3): B1≥0, B2≥0. Consider the following SDE:
Then X3(t)≤~X3(t) a.s. By Itô's formula, we get
In the light of Lemma 2, we obtain:
If B1≥0, B2≥0, B3<0, then limt→+∞~X3(t)=0 a.s.
If B1≥0, B2≥0, B3≥0, then
Hence, for arbitrary ζ>0,
Thanks to system (3.4) and Eq (3.7), we obtan
Based on Lemma 2, we obtain:
If B1≥0, B2≥0, B3<0, then limt→+∞Xj(t)=0 a.s. (j=3,4).
If B1≥0, B2≥0, B3≥0, then
Case(4): B1≥0, B2≥0, B3≥0. Consider the following SDE:
Then X4(t)≤~X4(t) a.s. By Itô's formula, we get
In view of Lemma 2, we obtain:
If B1≥0, B2≥0, B3≥0, B4<0, then limt→+∞~X4(t)=0 a.s.
If B1≥0, B2≥0, B3≥0, B4≥0, then
Hence, for arbitrary ζ>0,
Thanks to systems (3.4) and (3.8), we deduce
Based on Lemma 2 and the arbitrariness of ϵ, we obtain:
If B1≥0,B2≥0,B3<0,B4<0, then limt→+∞X4(t)=0a.s.
If B1≥0,B2≥0,B3≥0,B4≥0, then
The proof is complete. □
Lemma 5. For system (1.6):
(i) lim supt→+∞t−1lnxi(t)≤0 a.s. (i=1,2,3,4).
(ii) limt→+∞xi(t)=0⇒limt→+∞xj(t)=0 a.s. (1≤i<j≤4).
Proof. From Lemma 4, system (3.1) satisfies limt→+∞t−1lnXi(t)=0 a.s. (i=1,2,3,4). By the stochastic comparison theorem, we obtain the desired assertion (i). The proof of (ii) is similar to that of Lemma 4 and here is omitted. □
Theorem 2. Under Assumption 4 system (1.6) satisfies Table 3, where
Proof. Compute |Δ4|<A43|Ξ3|<A32A43|A2|<A21A32A43B1. By Eq (3.8), for any ζ>0,
By Itô's formula, we compute
Case(i): |Δ4|>0. According to system (3.9), we compute
In view of Lemma 5 (i) and Lemma 2, we get
Based on system (3.9), we compute
By Lemma 5 (i) and Eq (3.12), for ∀ϵ∈(0,1) and t≫1,
In view of Eq (3.11), we deduce
where |A|>0 and |Δ|>0. From Lemma 3, we have |Δ1|>0. Based on Lemma 2 and the arbitrariness of ϵ, we obtain
According to system (3.9), we compute
Thanks to Lemma 5 (i) and Eq (3.15), for ∀ϵ∈(0,1) and t≫1,
If
then by Lemma 2 and the arbitrariness of ϵ, we obtain
If
since lim inft→+∞t−1∫t0x3(s)ds≥0, we also obtain Eq (3.16).
According to system (3.9), Eq (3.14) and Eq (3.16), for ∀ϵ∈(0,1) and t≫1,
On the basis of Eq (3.13), Eq (3.14) and Eq (3.16), we have
From Lemma 3, we have |Δ2|>0. By Lemma 2 and the arbitrariness of ϵ, we obtain
By system (3.9), Eq (3.11) and Eq (3.18), for ∀ϵ∈(0,1) and t≫1,
In view of Eq (3.17) and Eq (3.18), we obtain
In the light of Lemma 3, we have |Δ3|>0. Thanks to Lemma 2 and the arbitrariness of ϵ, we obtain
By system (3.9) and Eq (3.20), for ∀ϵ∈(0,1) and t≫1,
Thanks to Eq (3.19), we obtain
In the light of Lemma 2 and the arbitrariness of ϵ, we obtain
In other words, we have
According to Eq (3.21) and Assumption 4, we obtain
Combining Eq (3.11) and Eq (3.22) yields
Substituting Eq (3.22) into Eq (3.14) yields
Substituting Eq (3.22) into Eq (3.16) yields
Substituting Eq (3.24) and Eq (3.25) into Eq (3.18) yields
Substituting Eq (3.26) into Eq (3.20) yields
Combining Eq (3.25) and Eq (3.27) yields
In view of system (3.9), we have
By Eq (3.26) and Eq (3.29), for ∀ϵ∈(0,1) and t≫1,
where B1−A12|Δ2||Δ|=A11|Δ1||Δ|. From Lemma 3, we have |Δ1|>0. According to Lemma 2 and the arbitrariness of ϵ, we have
Combining Eq (3.24) with Eq (3.30) yields
Substituting Eq (3.31) into system (3.9) yields
From Lemma 3, we have |Δ2|>0. By Lemma 2, we obtain
Case(ii): |Ξ3|>0>|Δ4|. Thanks to Eq (3.10), we deduce
which implies
From Lemma 5 (ii) and Eq (3.33), we obtain
In other words,
According to system (3.9), we compute
Combining Lemma 5 (i) with Lemma 2 yields
Based on system (3.9), we compute
By Lemma 5 (i) and Eq (3.37), for ∀ϵ∈(0,1) and t≫1,
On the basis of Eq (3.36), we deduce
In view of Lemma 2 and the arbitrariness of ϵ, we obtain
According to system (3.9), Eq (3.36) and Eq (3.38), for ∀ϵ∈(0,1) and t≫1,
Combining Eq (3.38) with system (3.39) yields
In the light of Lemma 2 and the arbitrariness of ϵ, we obtain
From system (3.9) and Eq (3.41), for ∀ϵ∈(0,1) and t≫1,
Thanks to Eq (3.40), we obtain
In the light of Lemma 2 and the arbitrariness of ϵ, we obtain
In other words, we have
In view of Eq (3.42) and Assumption 4, we obtain
Combining Eq (3.36) and Eq (3.43) yields
Substituting Eq (3.44) into Eq (3.38) yields
Substituting Eq (3.45) into Eq (3.41) yields
In view of system (3.9), we compute
By Lemma 5 (i) and Eq (3.47), for ∀ϵ∈(0,1) and t≫1,
From Eq (3.40), we have |Ξ2|>0 and |Ξ|>0. Based on system (3.48) and Lemma 2,
Combining Eq (3.46) with Eq (3.49) yields
Substituting Eq (3.50) into system (3.9) yields
From Eq (3.38), we have |Ξ1|>0 and |Ξ|>0. Therefore, by Lemma 2, we obtain
Case(iii): |A2|>0>|Ξ3|. Then limt→+∞x4(t)=0 a.s. Thanks to Eq (3.35), we deduce
which implies
According to Lemma 5 (ii) and Eq (3.52), we obtain
In other words, we derive
In view of Eq (3.47) and Eq (3.54), we get
Based on Lemma 5 (i) and Lemma 2, we have
In the light of Eq (3.37) and Eq (3.54), we obtain
By Lemma 5 (i) and Lemma 2, we obtain
Substituting Eq (3.54) and Eq (3.57) into system (3.9) yields
On the basis of Lemma 2 and the arbitrariness of ϵ, we have
Combining Eq (3.56) with Eq (3.58) yields
By system (3.9) and Eq (3.59), we compute
In the light of Lemma 2, we obtain
Case(iv): B1>0>|A2|. Then
According to Eq (3.55), we gain
Hence, lim supt→+∞xA211(t)xA112(t)=0. By Lemma 5 (ii) and Eq (3.62),
In other words, we derive
Substituting Eq (3.64) into system (3.9) yields
On the basis of Lemma 2 and the arbitrariness of ϵ, we obtain
Case(v): B1<0. Compute
By Eq (3.66), we have
In view of Lemma 2, we obtain limt→+∞x1(t)=0 a.s. According to Lemma 5 (ii), we get
The proof is complete. □
4.
Numerical examples
In this section we introduce some numerical examples to illustrate our main results. For simplicity, we suppose that S={1,2}. Then system (1.6) is a hybrid system of the following two subsystems:
and
with initial conditions x1(θ)=2eθ, x2(θ)=1.5eθ, x3(θ)=0.8eθ x4(θ)=0.5eθ and θ∈[−ln2,0].
Let rii=0.3, τji=ln2, μji(θ)=μjieθ, γj(μ,i)=γj(i) and λ(Z)=1, see Table 4. Denote
Then system (1.6) may be regarded as the result of regime switching between subsystems (4.1) and (4.2) with the following estimated parameters, respectively,
Compute |Δ|=0.066525, |Ξ|=0.1005 and |A|=0.195. Denote
4.1. The effects of Markovian switching on the persistence in mean and extinction
Let →r(1)=(0.9,0.5,0.3,0.2). Compute
Based on Theorem 2, all species in subsystem (4.1) are persistent in mean and
Let →r(2)=(0.6,0.3,0.2,0.1). Then B1=−0.1527<0. From Theorem 2, all species in subsystem (4.2) are extinctive.
Case1: (π1,π2)=(0.8,0.2). Compute
By Theorem 2, all species in system (1.6) are persistent in mean and
Case2: (π1,π2)=(0.6,0.4). Compute
From Theorem 2, x1(t), x2(t) and x3(t) are persistent in mean, while x4(t) is extinctive and
Case3: (π1,π2)=(0.5,0.5). Compute
Thanks to Theorem 2, x1(t) and x2(t) are persistent in mean, while x3(t) and x4(t) are extinctive and
Case4: (π1,π2)=(0.3,0.7). Compute
Based on Theorem 2, x1(t) is persistent in mean, while x2(t), x3(t) and x4(t) are extinctive and
Case5: (π1,π2)=(0.1,0.9). Compute B1=−0.0499<0. On the basis of Theorem 2, all species in system (1.6) are extinctive.
4.2. The effects of Lévy jumps on the persistence in mean and extinction
Let →r(1)=(0.7,0.5,0.3,0.2). We study the effects of Lévy jumps on the persistence in mean and extinction of system (4.1) by changing the values of γj(1) and setting the remaining parameters of the examples to be the same as those in system (4.1). Denote I4={−0.3,0.4}, α4∈I4; I3={−0.6,1.1}, α3∈I3; I2={−0.9,1.9}, α2∈I2; I1={−0.8,1.7}, α1∈I1.
4.2.1. The effects of γj(1) on the persistence in mean and extinction of system (4.1)
Case1: Let →γ(1)=(0.1,0.1,0.1,α4). Then |Δ4|<0, |Ξ3|=0.0606>0. According to Theorem 2, x1(t), x2(t) and x3(t) are persistent in mean, while x4(t) is extinctive.
Let →γ(1)=(0.1,0.1,0.1,0.1). Then |Δ4|=0.0047>0. By Theorem 2, all species in system (4.1) are persistent in mean. See Table 5.
Case2: Let →γ(1)=(0.1,0.1,α3,α4). Then |Ξ3|<0, |A2|=0.2478>0. Based on Theorem 2, x1(t) and x2(t) are persistent in mean, while x3(t) and x4(t) are extinctive. See Table 6.
Case3: Let →γ(1)=(0.1,α2,α3,α4). Then |A2|<0, B1=0.6753>0. From Theorem 2, x1(t) is persistent in mean, while x2(t), x3(t) and x4(t) are extinctive. See Table 7.
Case4: Let →γ(1)=(α1,α2,α3,α4). Then B1<0. Thanks to Theorem 2, all species are extinctive. See Table 8.
4.2.2. The effects of γ1(1) on the persistence in mean and extinction of system (4.1)
Case1: Let γ1(1)=−0.8. Then B1=−0.1294<0. According to Theorem 2, all species in system (4.1) are extinctive.
Let γ1(1)=−0.7. Then |A2|=−0.0518<0, B1=0.1760>0. By Theorem 2, x1(t) is persistent in mean, while x2(t), x3(t) and x4(t) are extinctive and
Let γ1(1)=−0.6. Then |Ξ3|=−0.0329<0, |A2|=0.0608>0. Based on Theorem 2, x1(t) and x2(t) are persistent in mean, while x3(t) and x4(t) are extinctive and
Let γ1(1)=−0.3. Then |Δ4|=−0.0023<0, |Ξ3|=0.0450>0. In view of Theorem 2, x1(t), x2(t) and x3(t) are persistent in mean, while x4(t) is extinctive and
Let γ1(1)=−0.1. Then |Δ4|=0.0046>0. From Theorem 2, all species are persistent in mean and
Case2: Let γ1(1)=0.2. Then |Δ4|=0.0029>0. Thanks to Theorem 2, all species are persistent in mean and
Let γ1(1)=0.6. Then |Δ4|=−0.0122<0, |Ξ3|=0.0230>0. On the basis of Theorem 2, x1(t), x2(t) and x3(t) are persistent in mean, while x4(t) is extinctive and
Let γ1(1)=0.9. Then |Ξ3|=−0.0155<0, |A2|=0.0957>0. By Theorem 2, x1(t) and x2(t) are persistent in mean, while x3(t) and x4(t) are extinctive and
Let γ1(1)=1.3. Then |A2|=−0.0297<0, B1=0.2129>0. From Theorem 2, x1(t) is persistent in mean, while x2(t), x3(t) and x4(t) are extinctive and
Let γ1(1)=1.7. Then B1=−0.0267<0. In view of Theorem 2, all species are extinctive.
5.
Discussion and conclusions
This paper concerns the dynamics of a stochastic hybrid delay food chain model with jumps in an impulsive polluted environment. Theorem 2 establishes sufficient and necessary conditions for persistence in mean and extinction of each species. Our results reveal that the stochastic dynamics of the system is closely correlated with both time delays and environmental noises.
Some interesting topics deserve further investigation, for instance, it is meaningful to consider the optimal harvesting problem of the stochastic hybrid delay food chain model with Lévy noises in an impulsive polluted environment. One may also propose some more realistic systems, such as considering the generalized functional response and the influences of impulsive perturbations. We will leave investigation of these problems to the future.
Use of AI tools declaration
The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.
Acknowledgments
This work is supported by National Natural Science Foundation of China (No. 11901166).
Conflict of interest
The authors declare that there are no conflicts of interest.