Stationary distribution and optimal control of a stochastic population model in a polluted environment

1.   Introduction

Along with fast development of modern industry and agriculture, there is no denying the fact that environmental pollution has become increasingly serious, including air pollution, water pollution and heavy metal pollution ^[1,2]. The existence of various poisons in the environment seriously threatens the survival of exposed populations, which has prompted many scholars to study the effects of toxins on the population and assess the risk of the population in a polluted environment. The population model in a polluted environment was first proposed by Hallam et al. in the 1980s ^[3,4]. From then on, more research results on the deterministic population model in a polluted environment have been presented (see ^[5,6,7,8]). However, in the real world, population dynamics are inevitably affected by randomness in the environment. Thus, several scholars have introduced random perturbations into the population model in a polluted environment to study dynamic behavior (see ^{[9,10,11,12,13]}). In particular, Liu and Wang ^[14] established a stochastic population in a polluted environment which takes the form

where $X_1(t)$ is the density of the population at time $t$, $X_2(t)$ denotes the concentration of toxicant in the organism (internal toxicant) at time $t$, and $X_3(t)$ stands for the concentration of toxicant in the environment (external toxicant) at time $t$. $B(t)$ is the standard Brownian motion, and $\sigma$ represents the intensity of the noise. The biological meanings of other parameters are listed in Table 1. For model (1.1), Liu and Wang obtained the acute threshold between local extinction and stochastic weakly persistent in the mean for population. However, the model did not consider the effects of randomness in the environment on the toxicant, nor did it consider the effects of population on external toxicant.

As is well known, toxicants are more or less disturbed by environmental factors, such as temperature, humidity, and seasonal climate change. Due to these random factors, the concentrations of toxicant in the organism and environment may exhibit a certain extent of random fluctuations in model parameters, such as the excretion rate and loss rate. Therefore, it is more realistic to consider the effects of randomness in the environment on the toxicant ^[15]. On the other hand, the uptake and egestion of toxicant by organisms from and to the environment can lead to obvious changes in toxicant in the environment, which is an important factor in determining the concentration of toxicant in the environment and needs to be considered ^[8]. Obviously, if these factors are further taken into account in model (1.1), we will get a novel stochastic population model in a polluted environment, which is an extension of model (1.1). For this model, the stationary distribution is an important concept.

Stationary distribution is an important dynamical property of the stochastic population model in a polluted environment, as it not only means random weak stability, but also provides a better description of persistence ^[16]. In addition, it is well known that the proliferation of pollutants will lead to huge economic costs, mainly including the reduction of crop yields and the expenditures related to pollution prevention. From an ecological and economic point of view, how to formulate the control strategy for pollutants is an important and meaningful issue. In this paper, we first show that the model has a unique ergodic stationary distribution by the Lyapunov method, then the optimal control strategy for pollutants is formulated, and the necessary conditions are derived for the existence of optimal control by using the Pontryagin maximum principle. In particular, our study has the following contributions:

$\bullet$ We propose a novel stochastic population model in a polluted environment and establish the sufficient conditions for the existence of an ergodic stationary distribution of the positive solutions.

$\bullet$ We study the optimal control problem for the toxicant and derive the necessary conditions for existence of optimal control by Pontryagin's maximum principle.

The remainder of this paper is organized as follows. In Section 2, the basic model and some necessary preliminary knowledge are introduced. In Section 3, the existence and uniqueness of a global positive solution of the model is proposed, and the sufficient conditions for the existence of an ergodic stationary distribution of the positive solution are established. In Section 4, we formulate the optimal control strategy for the toxicant and derive the first-order necessary conditions for the existence of optimal control. In Section 5, some numerical simulations are given to illustrate our main results. In Section 6, we summarize the work of this paper.

2.   Model formulation and preliminaries

2.1.   Model formulation

Similar to ^[8], we consider the changes in the concentration of external toxicant caused by absorption and excretion of the population, and then the third equation of system (1.1) becomes

where $k_1X_1X_3$ represents the absorption of the toxicant from the environment by the population, and $g_1X_1X_2$ represents the toxicant excreted by the population into the environment.

Following ^[17], we consider the effects of Brownian motion as a random factor on the net organismal excretion rate of toxicant, and the total loss rate of toxicant from the environment, that is,

Substitute (2.1) and (2.2) into system (1.1) to get the following stochastic population system in a polluted environment:

The initial conditions of system (2.3) take the following form:

where $\sigma_i(i = 1, 2, 3)$ is the intensity of the noises, and $B_i(t) (i = 1, 2, 3)$ denotes independent and standard Brownian motions. The biological significance of all parameters is shown in Table 1.

2.2.   Preliminaries

Let $(\Omega, \mathcal{F}, \{\mathcal{F}_t\}_{t\geq0}, P)$ be a complete probability space with a filtration $\{\mathcal{F}_t\}_{t\geq0}$ satisfying the usual conditions, and $B_i(t)(i = 1, 2, 3)$ are mutually independent standard Brownian motions defined on $(\Omega, \mathcal{F}, \{\mathcal{F}_t\}_{t\geq0}, P)$. Define $R_+^d = (x_1, x_2, \cdots, x_d)\in R^d:x_i > 0, i = 1, 2, \cdots, d$. Next, we introduce a lemma that gives a criterion for the existence of an ergodic stationary distribution to system (2.3).

Let $Y(t)$ be a homogeneous Markov process in $E_d$ ($d$-dimensional Euclidean space) described by the following stochastic differential equation:

The diffusion matrix is defined as follows:

Lemma 2.1. ^[18] The Markov process $Y(t)$ has a unique ergodic stationary distribution $\mu(\cdot)$ if there exists a bounded domain $D\subset E_d$ with regular boundary $\Gamma$ and the following:

(i) There is a positive number $M$ such that $\sum_{i, j = 1}^d\tilde{a}_{ij}(x)\xi_i\xi_j\geq M\lvert|\xi\rvert|^2, x\in D, \xi\in R^d$.

(ii) There exists a nonnegative $C^2$-function $V$ such that $LV$ is negative for any $E_d\backslash D$.

Then

for all $x\in E_d$, where $\lvert|\cdot|\rvert$ denotes the Euclidean norm, and $f(\cdot)$ is a function integrable with respect to measure $\mu$.

3.   Stationary distribution

Before establishing the sufficient conditions for the existence of an ergodic stationary distribution, we first prove that there exists a unique positive solution of system (2.3).

Theorem 3.1. For any initial value $(X_1^0, X_2^0, X_3^0)\in R_+^3$, there exists a unique solution $(X_1(t), X_2(t), X_3(t))$ of system (2.3) on $t\geq0$, and the solution will remain in $R_+^3$ with probability one, namely, $(X_1(t), X_2(t), X_3(t))$ $\in R_+^3$ for all $t\geq0$ almost surely (a.s.).

The proof is shown in the Appendix.

In the following, we establish sufficient conditions for the existence of a unique ergodic stationary distribution. Let

Theorem 3.2. If there exists a constant $0 < \theta < 1$ such that $\Lambda > 0$ holds, then for any initial date $(X_1^0, X_2^0, X_3^0)\in R_+^3$, system (2.3) exists a unique stationary distribution $\mu(\cdot)$, and it has the ergodic property.

Proof. To prove Theorem 1.1, only conditions (ⅰ) and (ⅱ) in Lemma 2.1 need to be verified. Now, we prove condition (ⅰ). The diffusion matrix of system (2.3) is given by

Choose $\tilde{M} = \min\limits_{(X_1, X_2, X_3)\in \bar{D}_k\subset R_+^3}\left\{\sigma_1^2X_1^2, \sigma_2^2X_2^2, \sigma_3^2X_3^2\right\}$, where $\bar{D}_k = \left[\frac{1}{k}, k\right]\times\left[\frac{1}{k}, k\right]\times\left[\frac{1}{k}, k\right]$ and $k > 1$ is a sufficiently large integer. By calculation, we have

for any $(X_1, X_2, X_3)\in\bar{D}_k$, $\xi = (\xi_1, \xi_2, \xi_3)\in R^3$. This completes the proof of condition (ⅰ) in Lemma 2.1.

Now, we will prove condition (ⅱ). We first define a $C^2$-function $\tilde{V}:R_+^3\rightarrow R$ as follows:

where $a_1$ and $b_1$ are positive constants, and their values are determined later, $\theta > 0$ is a sufficiently small constant such that $\Lambda > 0$, and $B > 0$ is a sufficiently large number that satisfies $-B(r_0+u)+B_1\leq-2$, where

Furthermore, $\tilde{V}$ is continuous and tends to $+\infty$ as $(X_1, X_2, X_3)$ approaches the boundary of $R_+^3$ and as $\lvert|(X_1, X_2, X_3)|\rvert\rightarrow +\infty$. Therefore, it must have a minimum point $(X_1^0, X_2^0, X_3^0)$ in the interior of $R_+^3$. Then, we define a $C^2$-function $\bar{V}:R_+^3\rightarrow R_+\cup\{0\}$ as

where $V_1(X_1, X_2, X_3) = B(-\ln X_1-X_3)$, $V_2(X_2) = -\ln X_2$, $V_3(X_3) = -\ln X_3$, $V_4(X_1, X_2, X_3) = \frac{1}{\theta+1}(X_1+a_1X_2+b_1X_3)^{\theta+1}-\tilde{V}(X_1^0, X_2^0, X_3^0)$.

Choosing $a_1 = \frac{r_1h}{2g_1k}$, $b_1 = \frac{r_1}{g_1}$, applying Itô's formula to $V_1(X_1, X_2, X_3)$, we have

Similarly,

and

where

Therefore, in view of (3.1)–(3.3), we obtain

Define a bounded closed set as follows:

where $0 < \epsilon < 1$ is a sufficiently small number and satisfies

where $B_i$ and $D_i > 0$, $i = 1, 2, 3,$ are positive constants, which will be given explicitly later. Choose

Clearly, $\Theta_\epsilon^c = \mathbb{R}_+^3\backslash\Theta^\epsilon = \cup_{i = 1}^6\Theta_i$. Then, we will prove that

We are going to prove it in the following six cases.

Case 1: If $(X_1, X_2, X_3)\in\Theta_1$, due to $X_1X_3 < \epsilon X_3\leq\epsilon\left(1+X_3^{\theta+1}\right)$, we have

where

Case 2: If $(X_1, X_2, X_3)\in\Theta_2$, by Eq (3.4), we have

where

Case 3: If $(X_1, X_2, X_3)\in\Theta_3$, according to Eq (3.4) and $X_1X_3 < \epsilon X_1\leq\epsilon\left(1+X_1^{\theta+1}\right)$, we have

where

Case 4: If $(X_1, X_2, X_3)\in\Theta_4$, by Eq (3.4), we have

where

Case 5: If $(X_1, X_2, X_3)\in\Theta_5$, it follows from Eq (3.4) that

where

Case 6: If $(X_1, X_2, X_3)\in\Theta_6$, from Eq (3.4), we can obtain that

where

Clearly, through the analysis of cases $1$ to $6$, we have

Therefore, condition (ⅱ) in Lemma 2.1 is satisfied. In view of Theorem 1.1, we obtain that system (2.3) admits a unique ergodic stationary distribution. This completes the proof.

Remark 1. The above analysis shows that the positive solution of system (2.3) has a unique ergodic stationary distribution. This means that the density of the population and the concentrations of internal toxicant and external toxicant will tend to a steady state under certain conditions.

4.   Optimal control problem

In this section, the optimal control problem for system (2.3) is formulated. Our goal is to reduce the concentration of internal toxicant and external toxicant, while keeping the cost to apply the control at the minimum level. We use $v(t)$ as the control variable to reduce the concentration of external toxicant (such as carrying out greening and tree planting activities, reducing the use of private cars, prohibiting littering, discharging factory wastewater after purification, etc.). Thus, the optimal control problem of system (2.3) is as follows:

According to our purpose, we construct the following objective function:

where $[0, T]$ is the entire time interval over which the control strategies are applied, and the constants $\tau_i$ $(i = 1, 2, 3)$ are positive weighted constants to make the terms of the integrand keep balance in the objective functional $J(v)$.

Remark 2. For the specific cost function, the practical significance is given. The term $\int_0^T\left(\tau_1X_2+\right.$ $\left.\tau_2X_3\right)dt$ gives the sum of the concentrations of internal toxicant and external toxicant over the time period $T$. The term $\int_0^T\left(\frac{1}{2}\tau_3v^2\right)dt$ gives the total cost of reducing the concentration of external toxicant. The term $h(X_2(T), X_3(T))$ is a function of toxicant at terminal time $T$. In particular, we assume that the cost is proportional to the degree of pollution and quadratic to the intensity of control.

The aim of the control problem is to seek an admissible control such that

where the control set $\mathcal{U}_{ad}$ is considered as

Next, we give the following theorem to illustrate the existence of the optimal control.

Theorem 4.1. There exist an optimal control $v^*\in\mathcal{U}_{ad}$ and the corresponding optimal state $X_1^*$, $X_2^*$, $X_3^*$ such that

subject to the control system (4.1).

Proof. We use the results in ^[19,20,21] to complete the proof. Note that both the control variable and state variables are nonnegative, and the control set $\mathcal{U}_{ad}$ is closed and convex. Then, the objective function (4.2) is convex with respect to the control variable $v(t)$. Furthermore, the optimal control is bounded. Therefore, the necessary condition for the existence of the optimal control $v^*$ is satisfied. This completes the proof.

By constructing the Hamiltonian function ^[20] and using Pontryagin's maximum principle ^[22], the first-order necessary conditions for the optimal control problem are given as follows.

Theorem 4.2. Let $v^*$ be the optimal control variable, and $X_1^*$, $X_2^*$ and $X_3^*$ are corresponding optimal state variables of the control system (4.1). Then, we have the following optimal control:

where $p(t) = (p_1(t), p_2(t), p_3(t))$ satisfies the following adjoint equation.

Proof. We define a Hamiltonian function [20,Section 3.3.1] $H:[0, T]\times\mathcal{U}_{ad}\times\mathbb{R}_+^3\times\mathbb{R}_+^3$ as follows:

Applying the general results in ^[20], $v^*$ is obtained by using the optimality condition $\frac{\partial H(t, v, p, q)}{\partial v} = 0$, and thus we have $v = \frac{p_3X_3}{\tau_3}$. Then, according to the condition that $0\leq v(x)\leq1$, the optimal control $v^*$ is obtained as follows:

Hence, the optimal value of the function can be obtained. This completes the proof.

5.   Numerical simulations

In this section, some numerical examples are presented to verify the theoretical results obtained above. The parameter values are chosen as follows:

The following subsection shows numerical examples of the stationary distribution.

5.1.   Numerical simulations of stationary distribution

Based on Milstein's higher-order method ^[23], the corresponding discrete equations of system (2.3) are

where $\varpi_{j, i}(j = 1, 2, 3, \; i = 1, \cdots, n)$ are independent random variables following the standard normal distribution $N(0, 1)$.

In Figure 2, we present the time series plot of the system solution and its corresponding histogram. It shows that there is a unique ergodic stationary distribution for system (2.3), and this is consistent with the result in Theorem 3.2.

Next, we show the different dynamic results of stochastic system (2.3) and the corresponding deterministic system under the same set of parameters. In Figure 3, we give the effects of different noise intensities for the stationary distribution of system (2.3) depending on time. The results show that the stochastic path fluctuates around the deterministic path, with larger fluctuations as the noise intensity becomes larger. The next subsection shows numerical examples of optimal control.

5.2.   Numerical simulations of the optimal control

The discretizations of system (4.1) and adjoint Eq (4.4) are as follows.

Then, we give the nonlinear conjugate gradient algorithm ^[24]. Assume the step size is $\Delta > 0$, and $T = n\Delta$, where $n$ is a positive integer. The time interval $[0, T]$ can be divided as

where

Next, we give a numerical example to demonstrate the effectiveness of the control strategy. The results show that the concentrations of internal and external toxicant decreased obviously after the control measures were implemented. In particular, the density of the population increased slightly as the concentration of the toxicant decreased. The corresponding simulation is shown in Figure 4. The optimal control $v(t)$ and optimal states $X_1(t)$, $X_2(t)$ and $X_3(t)$ are shown in Figures 5 and 6, respectively.

6.   Conclusions and discussion

In this paper, a stochastic population model in a polluted environment is developed and analyzed. By constructing a suitable Lyapunov function, the existence and uniqueness of a global positive solution is obtained, and then the sufficient conditions for existence of the unique ergodic stationary distribution of the positive solution are established. This means that the density of the population and the concentration of toxicant will tend to a steady state. Furthermore, we study optimal control of stochastic system (2.3). By using Pontryagin's maximum principle, the first-order necessary conditions are derived for the existence of optimal control. The results show that our control strategy can not only reduce the concentration of toxicant, but also have a positive effect on density of population.

Some interesting topics deserve further study. For example, some more realistic models can be considered, such as considering the effects of spatial heterogeneity, spatial diffusion and the impulse input of toxicant on system (2.3). In addition, in model (2.3), we only consider white noise, but one can also introduce colored noise into it and study the existence of an ergodic stationary distribution of the positive solutions to the considered model. These thoughts are interesting. We leave these to future work.

Acknowledgments

The research was supported by the Ningxia Key R & D Program Key Projects (No. 2021BEG03012), the Natural Science Foundation of China (No. 12161068) and the Natural Science Foundation of Ningxia (No. 2022AAC03265).

Conflict of interest

The authors declare there is no conflict of interest.

Appendix: The proof of Theorem 3.1.

Proof. According to the theory of stochastic differential equations ^[25] and the local Lipschitz continuity of the coefficients of the system (2.3), we know that there exists a unique local solution $(X_1(t), X_2(t), X_3(t))$ on $t\in[0, \tau_e)$ for any initial value $(X_1^0, X_2^0, X_3^0)\in R_+^3$, where $\tau_e$ is the explosion time ^[25]. To complete the proof, we only need to show that $\tau_e = +\infty$ a.s. To this end, let $k_0 > 0$ be sufficiently large such that $X_1^0$, $X_2^0$ and $X_3^0$ all lie within the interval $[\frac{1}{k_0}, k_0]$. For each integer $k\geq k_0$, define the stopping time

Obviously, $\tau_k$ increases as $k\rightarrow \infty$. Set $\tau_{\infty} = \lim\limits_{k\rightarrow +\infty}\tau_k$, and then $\tau_{\infty}\leq \tau_k$ a.s.. In the following, we only need to prove $\tau_{\infty} = +\infty$ a.s.. Therefore, let us consider a $C^2$-function $V$: $R_+^3\rightarrow R_+$ as follows:

where $a$ and $b$ are positive constants, and their values are determined later. The nonnegativity of this function can be seen from

For any $0\leq t\leq \tau_k\wedge T$, where $k\geq k_0$ and $T > 0$ are arbitrary, according to Itô's formula, we have

Choose $a = r_1/2(g+m), b = kr_1/2h(g+m)$, such that $\frac{1}{2}r_1-(g+m)a = 0, ka-hb = 0$. Note that

Then, we have

Now, if $X_1^{\frac{1}{2}}-2\geq0$ (i.e., $X_1\geq4$), then there exists a positive number $G_1$ independent of $x$ and $t$ such that

If $0 < X_1 < 4$, it is clear that there exists a positive number $G_2$ independent of $x$ and $t$ such that

In other words, we have already shown that there exists a positive number $G$ independent of $x$ and $t$ such that

The remaining proof is similar to [26,Theorem 2.1] and thus not introduced here. This completes the proof.