The vanishing discount problem for monotone systems of Hamilton-Jacobi equations. Part 1: linear coupling

1.   Introduction

Phytoplankton is the primary producer of water, it not only provides natural live bait for fish directly or indirectly, but also is the main producer of dissolved oxygen in water body, which is the basis of water productivity. Its quantity determines the density of herbivorous zooplankton, and also determines the output of herbivorous fish. Therefore, the amount of harvested adult fish yields basically depends on the phytoplankton density. Fish is not only a beneficial food source for human beings, but also an important economic source for fishermen. So, it is of great significance to study phytoplankton-fish models in Marine ecosystems. For example, the commercially valuable phytoplankton-fish predation model proposed in ^[1] discusses the impact of different factors on the level of fishery. Recent studies ^[2,3,4,5] suggest that persistent or overfishing of specific fish species may lead to the extinction of this species. This prompted the researchers to conduct a more in-depth study of the dynamic properties of the fish and their harvesting.

The basic mathematical tool for designing impulse control system is impulse differential equation theory. Impulse control is a normal form of control based on impulse differential equations ^[6]. There are many transient changes in nature, such as pest control, glucose regulation, vaccination, etc. ^[7,8,9]. When these situations occur, they are not disordered, they develop and change with regularity. To describe them, scholars put forward the time impulse differential equation ^[10,11]. Not all the changes are regular like green tide bursts, there are also some changes that are irregular and not limited by time. Scholars from mathematics put forward the state impulse differential equation model to describe them ^[12,13]. In recent years, pulsed semi-continuous dynamical systems are widely used in population dynamics ^[14,15]. For example, Zeng et al. ^[16] generalized the Poincar$\acute{\mbox{e}}$-Bendixon theorem for ordinary differential equations, and explored the properties of order-$1$ periodic solutions of predator-prey models with state feedback control. Chen et al. ^[17] first proposed the method of applying successor functions to study the dynamical simplicity of semi-continuous dynamical systems. Jiang et al. ^[18,19] applied Poincar$\acute{\mbox{e}}$ map to obtain the existence and stability conditions of the order-$1$ periodic solution of the impulsive semi-dynamic system. Hou et al. ^[20] investigated a class of predator-prey models with shelter and nonlinear impulse feedback control using Poincar$\acute{\mbox{e}}$ map. Due to the importance of pulsed semi-continuous dynamic systems in population, impulsive differential equations have attracted more and more attention in recent years, and have been applied in various fields from population dynamics to chemical regulation systems ^{[21,22,23,24,25,26,27,28]}.

In previous studies, most existing models only consider population size, rather than population growth rate, when proposing control strategies ^[29,30]. Usually, when a population reaches a threshold, the number of organisms released during the pulse is independent of population density. A biological point of view, a fixed threshold cannot be combined with the actual situation of the biological population to determine the control strategy. In reality, there are two situations of population growth: $1)$ The population density is small, but the rate of change is high, which often occurs in the early stage of population growth; $2)$ The population density is large, but the rate of change is small. Therefore, the threshold of the model needs to consider both the population density and the rate of change, which is more realistic. In the process of biological population control, it is necessary to comprehensively consider different factors, so that the adopted control method is more in line with the development law of biological populations ^[31]. Based on the above analysis, this paper proposes and studies a comprehensive control model that comprehensively considers fish population density and its current growth rate, and deeply studies the complex dynamic properties of the model, showing the impact of dynamic thresholds on system dynamics.

The rest of this article is organized as follows: In Section 2, we present a new phytoplankton-fish capture model with a new control strategy that takes into account the density of fish growth as well as the current fish population. In Section 3, we construct the Poincar$\acute{\mbox{e}}$ map for the model, and then analyze some properties of the Poincar$\acute{\mbox{e}}$ map. Then, in Section 4, we discuss the existence of periodic solutions of order-$k (k \ge 1)$ of the model, as well as the uniqueness and global stability of periodic solutions of order-$1$. Then, some numerical simulations are performed in Section 5. Finally, our findings are discussed in Section 6.

2.   Model formulation

To formulate our model, we first introduce the following notations:

Furthermore let $\frac{au}{b+u}$ represent the number of fish, due to the distribution of phytoplankton toxicants, and all parameters are positive and $0 < \delta < 1$. Then, Li et al. proposed the following phytoplankton-fish model with the impulsive feedback control ^[32]

Now, it is known that the state pulse capture strategy in (2.1) has certain shortcomings, as it only considers the population size and does not consider the population growth rate. In addition, the control of the pulse model in literature ^[33,34] is also single. We propose a feedback control model that considers both population density and population rate, which brings difficulties to the analysis of the model. In literature ^[32,33,34], the pulse set and phase set are straight lines on the plane, while in this paper, the phase set is a curve and the determination of phase set is complicated, which also adds a lot of difficulties to the numerical simulation. More precisely, the control threshold is described by

where ${a_1}$ and ${b_1}$ are positive constants, satisfying ${a_1} + {b_1} = 1$, and they are the weight coefficients that control the rate of change in fish population and the rate of change in fish population density, respectively. In real life, we should fish reasonably when ${a_1}v + {b_1}\frac{{dv}}{{dt}} \geq H$ is met. We only assumes that fishing is carried out when ${a_1}v + {b_1}\frac{{dv}}{{dt}} = H$ is reached. Then, we reach a model in the following form:

which is known as a semi-continuous dynamic system ^[33,35]. Without the feedback control, the system (2.3) reduces to an ODE system

According to reference ^[24], we know that system (2.4) has two equilibrium points in the range of real numbers: a saddle point, $O(0, 0)$ and a stable positive equilibrium $E^{*}(u^{*}, v^{*})$, if and only if, $\left (d+a-\beta b\right)^{2}\geq 4\beta bd$ and $d+a > \beta b$. Where,

determined by the intersection of the two nullclines

Next, we construct the Poincar$\acute{\mbox{e}}$ map to investigate the dynamics of system (2.3).

[Note]:When $\frac{du}{dt} < 0$ and $\frac{dv}{dt} < 0$, the growth rate of its phytoplankton and fish populations is negative. Mathematically speaking, their numbers gradually decrease, and with the increase of time, they will become extinct without control.

3.   Poincar$\acute{\mbox{e}}$ map of system (2.3) and its properties

By using the Eq (2.2) and the second equation in (2.3), we obtain

From a biological point of view, we only consider the properties of the model in region $R_{2}^{+} = \{(u, v), u\geq 0, v\geq 0\}$. To reveal the properties of the Poincar$\acute{\mbox{e}}$ map, without loss of generality, we assume the initial population size $v_0^ +$ of the fish population satisfies

and the initial point $\left({u_0^ +, v_0^ + } \right)$ satisfies

where $M\left (v\right) = \frac{\left (\theta v-1\right)+\sqrt{\left (\theta v-1\right)^{2}+4\theta \left (\tau +v\right)}}{2\theta }$. Then, in order to define the impulse set and phase set of the model, we introduce the following sets:

where

We further assume $\Sigma_{N}$ and $\Sigma_{M}$ intersect with $L_{2}$ at $A^{+}$ and $A^{-}$, respectively, then assume $\Sigma_{N}$ and $\Sigma_{M}$ intersect with $x$ at ${B^1}\left({{u^1}, 0} \right)$ and ${B^2}\left({{u^2}, 0} \right)$, respectively. Then, depending on the locations of ${\Sigma _N}$ and $v^{*}$, namely, the positive equilibrium point, we can define the impulse set and phase set of the Poincar$\acute{\mbox{e}}$ map.

3.1.   Case Ⅰ: ${v^{{*}}} \ge \frac{H}{{{a_1} + {b_1}\beta {u^{{*}}} - {b_1}d - \frac{{a{b_1}{u^{{*}}}}}{{b + {u^{{*}}}}}}}$

In this case, ${E^{\rm{*}}}$ is above the impulse set $\Sigma_{M}$ (see Figure 1(a)). Thus, there exists a point $T$ on the ${\Sigma _N}$, at which the trajectory of (2.3) is tangent to the ${\Sigma _N}$. Denote the intersection between the trajectory and $\Sigma_{M}$ by ${Q_1}$. Then, we define the impulse set by

and the phase set by

3.2.   Case Ⅱ: ${v^{{*}}} < \frac{H}{{{a_1} + {b_1}\beta {u^{{*}}} - {b_1}d - \frac{{a{b_1}{u^{{*}}}}}{{b + {u^{{*}}}}}}}$ and ${v_{{{W^ - }}}} > \frac{H}{{{a_1} + {b_1}\beta \sigma ({u^{{*}}}) - {b_1}d - \frac{{a{b_1}\sigma ({u^{{*}}})}}{{b + \sigma ({u^{{*}}})}}}}$

In this case, ${E^{\rm{*}}}$ is below the impulse set $\Sigma_{M}$ (see Figure 1(b)). We first assume that ${\Gamma _M}$ is the trajectory passing through $M$ in the system. If there is a point $W$ on the $\Sigma_{M}$, so that ${\Gamma _W}$ is just tangent to the $\Sigma_{M}$, when ${v^{\rm{*}}} < \frac{H}{{{a_1} + {b_1}\beta {u^{\rm{*}}} - {b_1}d - \frac{{a{b_1}{u^{\rm{*}}}}}{{b + {u^{\rm{*}}}}}}}$ and ${v_{{{W^ - }}}} > \frac{H}{{{a_1} + {b_1}\beta \sigma ({u^{\rm{*}}}) - {b_1}d - \frac{{a{b_1}\sigma ({u^{\rm{*}}})}}{{b + \sigma ({u^{\rm{*}}})}}}}$, there must be a point $T$ on the $\Sigma_{N}$ so that ${\Gamma _T}$ is just tangent to the $\Sigma_{N}$ and ${\Gamma _T}$ tangent to the ${\Sigma _M}$ at the point ${Q_2}$. Then, in this case, we define the impulse set by

and the corresponding phase set by

3.3.   Case Ⅲ: ${v_{{W^ - }}} < \frac{H}{{{a_1} + {b_1}\beta \sigma ({u^{{*}}}) - {b_1}d - \frac{{a{b_1}\sigma ({u^{{*}}})}}{{b + \sigma ({u^{{*}}})}}}}$

In this case, ${E^{\rm{*}}}$ is below the impulse set $\Sigma_{M}$, as shown in Figure 1(c). Thus, there is a trajectory denoted by $\Gamma _{W^{-}}$ that is tangent to ${\Sigma _M}$ at $W$, and intersects the $\Sigma_{N}$ at point ${P_1}$ and point ${P_2}$, respectively. Thus, we define the impulse set as

and the phase set as

Based on the above defined impulse and phase sets, we are now ready to construct the Poincar$\acute{\mbox{e}}$ map. To this end, let

where $0 < u_{k}^{+} < + \infty$, and define the trajectory

passing through $A_{k}^{+}$ that will reach the $\Sigma_{M}$ at point ${A_{k + 1}}\left({{u_{k + 1}}, {v_{k + 1}}} \right)$ at time ${{t_1}}$, where

Then, there is

It means that $u_{k+1}^{+}$ is determined by $u_{k}^{+}$. Since point $A_{k+1}$ is on the impulse set, $A_{k+1}$ jump to point

where

From model (2.3), we obtain

Let $v_0^ + = \frac{H}{{{a_1} + {b_1}\beta \sigma (u_k^ +) - {b_1}d - \frac{{a{b_1}\sigma (u_k^ +)}}{{b + \sigma (u_k^ +)}}}}$, $u_{0}^{+} = S$, then we obtain that $(u_{0}^{+}, v_{0}^{+})$ is in $\Sigma_{N}$.

We define

Then, according to the model (3.2),

From Eqs (3.2) and (3.3), the Poincar$\acute{\mbox{e}}$ map expression of system (2.2) is

Then, we can show that it has the following properties in Theorem 3.1.

Theorem 3.1. If $\frac{H}{{{a_1} + {b_1}\beta {u^{\rm{*}}} - {b_1}d - \frac{{a{b_1}{u^{\rm{*}}}}}{{b + {u^{\rm{*}}}}}}} < {v^*}$, then the trajectory ${\Gamma _T}$ intersects with $\Sigma_{M}$ at point ${Q_1}({u_{{Q_1}}}, {v_{{Q_1}}})$, and we have:

${\bf{(i)}}$ The domain of $G_{m}(S)$ is $(0, {{u^1}})$, $G_{m}(S)$ is monotonically decreasing on $(0, {u_{T}}]$, and monotonically increasing on $({u_{T}}, {{u^1}})$.

${\bf{(ii)}}$ $G_{m}(S)$ is continuously differentiable on $(0, {{u^1}})$.

${\bf{(iii)}}$ When $G_{m}({u_{T}}) > {u_{T}}$, it has a unique fixed point on $({u_{T}}, +\infty)$ (see Figure 2(b)). When $G_{m}({u_{T}}) < {u_{T}}$, it has a unique fixed point on $(0, {u_{T}})$ (see Figure 2(a)). When $G_{m}({u_{T}}) = {u_{T}}$, ${u_{T}}$ is the fixed point.

Proof. (ⅰ) Note that $E^{*}(u^{*}, v^{*})$ is a linear center and $\frac{H}{{{a_1} + {b_1}\beta {u^{\rm{*}}} - {b_1}d - \frac{{a{b_1}{u^{\rm{*}}}}}{{b + {u^{\rm{*}}}}}}} < {v^*}$. Taking any point $A_{k}^{+}(u_{k}^{+}, v_{k}^{+})$ on $\Sigma_{N}$, the trajectory ${\Gamma _{A_k^ + }}$ will reach the $M_{1}$. So the domain of $G_{m}(S)$ is $(0, {{u^1}})$.

For any $u_{k_{1}}^{+}, u_{k_{2}}^{+} \in \left[{u_{T}}, {{u^1}} \right)$, and where $u_{k_{1}}^{+} < u_{k_{2}}^{+}$, it is easy to get

which implies that there is $G_{m}(u_{k_{1}}^{+}) < G_{m}(u_{k_{2}}^{+})$. Therefore, $G_{m}(S)$ is monotonically increasing on $\left[{u_{T}}, {{u^1}} \right)$.

When $u_{k_{1}}^{+}$, $u_{k_{2}}^{+} \in (0, {u_{T}})$, $u_{k_{1}}^{+} < u_{k_{2}}^{+}$. The trajectory ${\Gamma _{A_{{k_1}}^ + }}$ from point $A_{{k_1}}^ + \left({u_{{k_1}}^ +, v_{{k_1}}^ + } \right)$ and the trajectory ${\Gamma _{A_{{k_2}}^ + }}$ from point $A_{{k_2}}^ + \left({u_{{k_2}}^ +, v_{{k_2}}^ + } \right)$ will pass through $L_{2 }$ and intersect $\Sigma_{N}$ at points $A_{{k_{11}}}^ + \left({u_{{k_{11}}}^ +, v_{{k_{11}}}^ + } \right)$ and $A_{{k_{21}}}^ + \left({u_{{k_{21}}}^ +, v_{{k_{21}}}^ + } \right)$, respectively. Here $u_{k_{i1}}^{+}(i = 1, 2) \in \left[{{u_{T}}, {{u^1}} } \right)$ and $u_{k_{11}}^{+} > u_{k_{21}}^{+}$. It can be seen that $G_{m}(u_{k_{1}}^{+}) > G_{m}(u_{k_{2}}^{+})$. Therefore, $G_{m}(S)$ is monotonically decreasing on $(0, {u_{T}})$.

(ⅱ) Equation (3.2) suggests that $w(u, v)$ is continuously differentiable. Thus, $G_{m}(S)$ is continuously differentiable in the first quadrant.

(ⅲ) When $G_{m}({u_{T}}) = {u_{T}}$, then ${u_{T}}$ is the fixed point of function $G_{m}(S)$.

If $G_{m}({u_{T}}) > {u_{T}}$, $G_{m}({u_{T}})- {u_{T}} > 0$. Since $0 < \delta < 1$, there must be a value ${S^{\rm{*}}}$, such that

So there is at least one $\tilde u \in \left({{u_{T}}, {S^{\rm{*}}} } \right)$, satisfying ${G_m}(\tilde u) = \tilde u$.

When $G_{m}({u_{T}}) < {u_{T}}$, let ${G_m}({u_{T}}) = {u_1} < {u_{T}}$. We know that $G_{m}(S)$ is monotonically decreasing on $(0, {u_{T}}]$, so ${G_m}({u_1}) > {G_m}({u_{T}}) = {u_1}$, and because $G_{m}({u_{T}}) < {u_{T}}$, so there is at least one $\tilde u \in \left({u_1}, {u_{T}} \right)$, satisfies ${G_m}(\tilde u) = \tilde u$.

From the above analysis, we know that ${G_m}(S)$ has at least one fixed point. Next, we prove the uniqueness by contradiction. Assume the system has two fixed points, ${\tilde u_1}$ and ${\tilde u_2}$ respectively, so that ${G_m}({\tilde u_1}) = {\tilde u_1}$ and ${G_m}({\tilde u_2}) = {\tilde u_2}$. Let ${\tilde u_1} < {\tilde u_2}$, we define

Differentiate the last equation with respect to $t$ obtain

Let

then

so

that is

i.e.

From system (2.3):

that is

It is easy to know that $\frac{{{{\tilde u}_1}}}{{\mu \left({{{\tilde u}_1}} \right)}} < 1$ if ${a_1} + {b_1}\beta {u^*} - {b_1}d - \frac{{a{b_1}{u^*}}}{{b + {u^*}}} < {v^*}$, so

which is contradictory with

The fixed point is unique.

When $\frac{H}{{{a_1} + {b_1}\beta {u^{\rm{*}}} - {b_1}d - \frac{{a{b_1}{u^{\rm{*}}}}}{{b + {u^{\rm{*}}}}}}} > {v^*}$ and $\Gamma _{T}$ intersects with $\Sigma_{M}$ (case Ⅱ), $G_{m}(S)$ of system (2.3) has similar properties with case Ⅰ.

Theorem 3.2. Suppose $\frac{H}{{{a_1} + {b_1}\beta {u^{\rm{*}}} - {b_1}d - \frac{{a{b_1}{u^{\rm{*}}}}}{{b + {u^{\rm{*}}}}}}} > {v^*}$. Then, then trajectory $\Gamma _{W^{-}}$ is tangent to ${\Sigma _M}$ at $W$ and intersect $\Sigma_{N}$ at point ${P_1}\left({{u_{{P_1}}}, {v_{{P_1}}}} \right)$ and point ${P_2}\left({{u_{{P_2}}}, {v_{{P_2}}}} \right)$, respectively, where $u_{{P_1}} > u_{{P_2}}$. Furthermore, $G_{m}(S)$ has the following properties:

${\bf{(i)}}$ The domain of $G_{m}(S)$ is $\left({0, u_{{P_2}}} \right] \cup \left[{u_{{P_1}}}, {{u^1}} \right)$, and $G_{m}(S)$ is monotonically decreasing on $\left({0, u_{{P_2}}} \right]$ and monotonically increasing on $\left[{u_{{P_1}}, {{u^1}} } \right)$.

${\bf{(ii)}}$ When $G_{m}(u_{{P_2}})\leq u_{{P_2}}$, there is a unique fixed point on $\left({0, u_{{P_2}}} \right]$ (see Figure 3(a)). When $G_{m}(u_{{P_2}}) > u_{{P_2}}$, there is no fixed point (see Figure 3(b)).

Proof. (ⅰ) Because $E^{*}(u^{*}, v^{*})$ is the center point, if $\frac{H}{{{a_1} + {b_1}\beta {u^{\rm{*}}} - {b_1}d - \frac{{a{b_1}{u^{\rm{*}}}}}{{b + {u^{\rm{*}}}}}}} > {v^*}$ and the trajectory $\Gamma _{W^{-}}$ is tangent to ${\Sigma _M}$ at $W$ and intersect with the $\Sigma_{N}$ at point ${P_1}\left({{u_{{P_1}}}, {v_{{P_1}}}} \right)$ and point ${P_2}\left({{u_{{P_2}}}, {v_{{P_2}}}} \right)$, respectively, $u_{{P_1}} > u_{{P_2}}$. Then, take any point $A_{k}^{+}(u_{k}^{+}, v_{k}^{+})$ in $\Sigma_{N}$, if $u_k^ + \in \left({0, u_{{P_2}}} \right] \cup \left[{u_{{P_1}}, + \infty } \right)$, the trajectory ${\Gamma _{A_k^ + }}$ will reach the $M_{3}$ at point $A_{k+1}(u_{k}, v_{k})$, if $u_k^ + \in \left(u_{{P_2}}, u_{{P_1}} \right)$ the ${\Gamma _{A_k^ + }}$ has no intersection with the $M_{3}$. So the domain of $G_{m}(S)$ is $\left({0, u_{{P_2}}} \right] \cup \left[{u_{{P_1}}, {{u^1}} } \right)$.

For any $u_{k_{1}}^{+}, u_{k_{2}}^{+} \in \left[{u_{{P_1}}, {{u^1}} } \right)$ and $u_{k_{1}}^{+} < u_{k_{2}}^{+}$ easy to get

so $G_{m}(u_{k_{1}}^{+}) < G_{m}(u_{k_{2}}^{+})$. Therefore, $G_{m}(S)$ is monotonically increasing on $\left[{u_{{P_1}}, {{u^1}} } \right)$.

When $u_{k_{1}}^{+}$, $u_{k_{2}}^{+} \in \left({0, u_{{P_2}}} \right]$, $u_{k_{1}}^{+} < u_{k_{2}}^{+}$. The trajectory from $A_{{k_1}}^ + \left({u_{{k_1}}^ +, v_{{k_1}}^ + } \right)$ and $A_{{k_2}}^ + \left({u_{{k_2}}^ +, v_{{k_2}}^ + } \right)$ will pass through the $L_{2 }$ intersect with the $\Sigma_{N}$ at point $A_{{k_{11}}}^ + \left({u_{{k_{11}}}^ +, v_{{k_{11}}}^ + } \right)$ and point $A_{{k_{21}}}^ + \left({u_{{k_{21}}}^ +, v_{{k_{21}}}^ + } \right)$, respectively. Then, $u_{k_{i1}}^{+}(i = 1, 2) \in \left[{{u_{32}}, {{u^1}}} \right)$ and $u_{k_{11}}^{+} > u_{k_{21}}^{+}$. So $G_{m}(u_{k_{1}}^{+}) > G_{m}(u_{k_{2}}^{+})$, therefore, $G_{m}(S)$ is monotonically decreasing on $\left({0, u_{{P_2}}} \right]$.

(ⅱ) Consider the relationship between $G_{m}(u_{{P_2}})$ and $u_{{P_2}}$:

(a) When $G_{m}(u_{{P_2}}) \le u_{{P_2}}$ (see Figure 3(a)), we assume ${G_m}(u_{{P_2}}) = {u_1} \le u_{{P_2}}$. We know that $G_{m}(S)$ is monotonically decreasing on $(0, u_{{P_2}}]$, so ${G_m}(u_{{P_2}}) \ge {G_m}(u_{{P_2}}) = {u_1}$, and because $G_{m}(u_{{P_2}}) \le u_{{P_2}}$, so there is a point $\tilde u \in \left({{u_1}, u_{{P_2}}} \right]$ satisfies ${G_m}(\tilde u) = \tilde u$.

(b) When $G_{m}(u_{{P_2}}) > u_{{P_2}}$ (see Figure 3(b)), there is no $\tilde u \in \left({{0}, u_{{P_2}}} \right]$ to satisfies ${G_m}(\tilde u) = \tilde u$.

For any ${u_k} \in \left({u_{{P_1}}, {{u^1}} } \right)$, trajectory of point $A_{k}(u_{k}, v_{k})$ on $\Sigma_{N}$ is tangent to the $\Sigma_{M}$ at point $A_k^ + (u_k^ +, v_k^ +)$. The point $A_k^ + (u_k^ +, v_k^ +)$ will be pulsed to ${A_{k + 1}}({u_{k + 1}}, {v_{k + 1}})$. To get $u_{k + 1} < u_k^ + < {u_k}$, that is ${u_{k + 1}} \ne {u_k}$, so there is no $\tilde u \in \left({u_{{P_1}}, {{u^1}} } \right)$ satisfies ${G_m}(\tilde u) = \tilde u$.

4.   The order-$k$ $(k\ge1)$ periodic solution of the semi-continuous dynamic system (2.3) and its stability

Theorem 3.1 has proved that system (2.3) has a fixed point under certain conditions, that is, the system has an order-1 periodic solution. Below, we will give more properties of the order-1 solution, and also the existence of the order-$k$ solutions.

Theorem 4.1. The order-1 periodic solution of system (2.3) is globally asymptotically stable if $\frac{H}{{{a_1} + {b_1}\beta {u^{\rm{*}}} - {b_1}d - \frac{{a{b_1}{u^{\rm{*}}}}}{{b + {u^{\rm{*}}}}}}} < {v^*}$ and $G_{m}({u_{{Q_1}}}) > {u_{{Q_1}}}$.

Proof. When $G_{m}({u_{{Q_1}}}) > {u_{{Q_1}}}$, $G_{m}(S)$ has a fixed point $\tilde u$ on $({u_{{Q_1}}}, {{u^1}})$, that is, ${G_m}(\tilde u) = \tilde u$ from the (ⅲ) of Theorem 3.1. For any point $A_{0}^{+}(u_{0}^{+}, v_{0}^{+})$ in $\Sigma_{N}$, where $u_{0}^{+} > {u_{{Q_1}}}$, ${\Gamma _{A_0^ + }}$ will intersect $\Sigma_{M}$ at the point $p_{1}^{+}(u_{1}^{+}, v_{1}^{+})$, which is $G_{m}(u_{0}^{+}) = u_{1}^{+}$, repeating the above process,

that is,

further available,

Then, we have the following three cases:

Case 1: ${u_{{Q_1}}} < u_0^ + \le \tilde u$. Since $G_{m}({u_{{Q_1}}}) > {u_{{Q_1}}}$ and $G_{m}(S)$ is monotonously increasing on $({u_{{Q_1}}}, {\rm{ + }}\infty)$. Let ${G_m}(u_i^ +) = u_{i+1}^{+}$ yields

Following the same fashion, we obtain that

Then we can get

Case 2: $\tilde u < u_0^ + < {{u^1}}$. In this case, we have

and

By mathematical induction,

Thus,

Case 3: $0 < u_0^ + < {u_{{Q_1}}}$. Since $G_{m}({u_{{Q_1}}}) > {u_{{Q_1}}}$ and $G_{m}(S)$ is decreasing on $(0, {u_{{Q_1}}})$, we obtain that ${G_m}(u_0^ +) > {G_m}({u_{{Q_1}}}) > {u_{{Q_1}}}$ for any $u_0^ + \in (0, {u_{{Q_1}}})$. So, we can conclude that ${G_m}(u_0^ +) > {u_{{Q_1}}}$. When ${u_{{Q_1}}} < {G_m}(u_0^ +) < \tilde u$, this is the situation in Case 1; when ${G_m}(u_0^ +) > \tilde u$, this is Case 2 above.

Thus, we always have

The conclusion is proved.

Theorem 4.2. The semi-continuous dynamical system (2.3) has a stable order-1 periodic solution or an order-2 periodic solution when $\frac{H}{{{a_1} + {b_1}\beta {u^{\rm{*}}} - {b_1}d - \frac{{a{b_1}{u^{\rm{*}}}}}{{b + {u^{\rm{*}}}}}}} < {v^*}$, $G_{m}({u_{{Q_1}}}) < {u_{{Q_1}}}$ and $G_{m}^{2}({u_{{Q_1}}}) < {u_{{Q_1}}}$.

Proof. Take a point $A_{0}^{+}(u_{0}^{+}, v_{0}^{+})$ on $\Sigma_{N}$. Since $E^{*}$ is a center point, ${\Gamma _{A_0^ + }}$ will intersect $\Sigma_{M}$ at point $A_{1}(u_{1}, v_{1})$, where $u_{1} = \mu(u_{0}^{+}, v_{0}^{+})$. $p_{1}$ will reach point $p_1^ + (u_1^ +, v_1^ +)$ by an impulse, so $G_{m}(u_{0}^{+}) = u_{1}^{+}$. Repeat the above process

When ${u_{{Q_1}}} \le u_0^ + < {{u^1}}$, $G_{m}(S)$ is increasing on $[{u_{{Q_1}}}, {{u^1}})$ and from the (i) of Theorem 3.1, we know there is no fixed point on $[{u_{{Q_1}}}, {{u^1}})$. Therefore, there is a positive integer $i$ that satisfies

and

When $0 < u_0^ + < {u_{{Q_1}}}$, the ${\Gamma _{A_0^ + }}$ will pass through the $L_{2}$ and intersects the ${\Sigma _N}$ at $p_{0}^{'}(u_{0}^{'}, H{\rm{ - }}\frac{\tau }{{1{\rm{ + }}\theta H}})$, where $0 < u_{0}^{'} < {u_{{Q_1}}}$, this can be transformed into the above situation.

So, for any $u_{0}^{+} \in (0, {{u^1}})$, there always be $i$ satisfies

So, we just consider the initial point $p_{0}^{+}(u_{0}^{+}, v_{0}^{+})$, where ${u_{{Q_1}}} \leq u_{0}^{+} < G_{m}({u_{{Q_1}}})$. Since $G_{m}(S)$ is monotonically decreasing on $[G_{m}({u_{{Q_1}}}), {u_{{Q_1}}}]$, we have

Let $G_ {m}(u_{0}^{+})\neq u_{0}^{+}$ and $G_m^2(u_0^ +)\; \ne u_0^ +$. Consider the following four situations:

Case Ⅰ: ${u_{{Q_1}}} \ge u_0^ + > {G_m}^2{(u_0^ +)} > {G_m}(u_0^ +) \ge {G_m}({u_{{Q_1}}})$. According to the monotonicity of ${G_m}(S)$:

furthermore,

Thus, there is

Proved by mathematical induction:

Case Ⅱ: ${u_{{Q_1}}} \ge u_1^ + > u_0^ + > u_2^ + \ge {G_m}({u_{{Q_1}}})$. $G_{m}(S)$ is monotonically decreasing on $[G_{m}({u_{{Q_1}}}), {u_{{Q_1}}}]$. So,

and

then

so

Case Ⅲ: ${u_{{Q_1}}} \ge u_1^ + > u_2^ + > u_0^ + \ge {G_m}({u_{{Q_1}}})$.

To get

Case Ⅳ: ${u_{{Q_1}}} \ge u_2^ + > u_0^ + > u_1^ + \ge {G_m}({u_{{Q_1}}}).$

This situation is similar to Case Ⅰ, can get

For Case Ⅰ and Case Ⅲ, there is $\tilde{u} \in ({u_{{Q_1}}}, G_{m}({u_{{Q_1}}}))$, so that

That is to say, in these two cases, system (2.3) has a stable order-1 periodic solution.

For case Ⅱ and case Ⅳ, $\tilde{u_{1}} \neq \tilde{u_{2}}$ and

That is to say, system (2.3) has a stable order-2 periodic solution in these two cases.

Theorem 4.3. Let $H_{1} < {u_{{Q_1}}}$ and $G_{m}({u_{{Q_1}}}) < {u_{{Q_1}}}$. Then, the necessary and sufficient condition for the order-1 periodic solution of system (2.3) to be globally stable is $G_{m}^{2}(u^{+}) < u^{+}$ for any $u^{+} \in (0, {u_{{Q_1}}}]$.

Proof. Sufficiency:

When $G_{m}({u_{{Q_1}}}) < {u_{{Q_1}}}$, there exists $\tilde{u} \in (G_{m}({u_{{Q_1}}}), {u_{{Q_1}}})$ which satisfies $G_{m}(\tilde{u}) = \tilde{u}$.

For any $u^{+} \in (\tilde{u}, {u_{{Q_1}}})$, we make $u_{1}^{+} = G_{m}(u^{+})$ and $u_{2}^{+} = G_{m}(u_{1}^{+}) = G_{m}^{2}(u^{+})$, since $G_{m}^{2}(u^{+}) < u^{+} < {u_{{Q_1}}}$ to get $\tilde{u} > u_{1}^{+} > G_{m}({u_{{Q_1}}})$, ${u_{{Q_1}}} > u_{2}^{+} > u_{4}^{+} > \tilde{u} > u_{3}^{+} > u_{1}^{+} > G_{m}({u_{{Q_1}}})$, so, ${u_{{Q_1}}} > u_{2n}^{+} > \tilde{u} > u_{2n+1}^{+} > G_{m}({u_{{Q_1}}})$.

So

Necessity: If $G_{m}^{2}(u^{+}) < u^{+}$ is not true for any $u^{+} \in [\tilde{u}, {u_{{Q_1}}}]$, there exists a maximum $u_{0} \in [\tilde{u}, {u_{{Q_1}}}]$ satisfies ${G_m}^{2}{({u_0})} \ge {u_0}$. There $u_{1}$ and $G_{m}^{2}(u_{1}) < u_{1}$ for any $\tilde{u}-\varepsilon < u_{1} < \tilde{u}+\varepsilon$, where $\varepsilon > 0$. From the continuity of $G_{m}^{2}(u)$ on the $[{u_0}, {u_1}]$ there is at least one $\vec u \in \left[{{u_0}, {u_1}} \right]$ and ${G_m}^2(\vec u) = \vec u$. This is contradictory.

Theorem 4.4. If $G_{m}({u_{{Q_1}}}) < {u_{{Q_1}}}$ and there is $u_{m}^ +$, such that $u_{m}^ + = min\left\{ {{u^ + } : {G_m}({u^ + }) = {u_{{Q_1}}}} \right\}$, the system (2.3) has a order-3 periodic solution if $G_{m}^{2}({u_{{Q_1}}}) < u_{m}^{+}$.

Proof. If $G_{m}({u_{{Q_1}}}) < {u_{{Q_1}}}$, there is a order-1 periodic solution in $\left({{G_m}({u_{{Q_1}}}), {u_{{Q_1}}}} \right)$, i.e., $G_{m}(\tilde{u}) = \tilde{u}$, here, $\tilde{u} \in \left({{G_m}({u_{{Q_1}}}), {u_{{Q_1}}}} \right)$. Since ${G_m}(\tilde u) = \tilde u$, so there is $u_{m}^ +$ such that $u_{m}^{+} \in (0, \tilde{u})$ and $G_{m}(u_{m}^{+}) = {u_{{Q_1}}}$. Further, $G_{m}^{3}(u_{m}^{+}) = G_{m}^{2}({u_{{Q_1}}}) < u_{m}^{+}$, on the other hand, $\mathop {\lim }\limits_{x \to 0} G_m^3(x) > x$, so there is at least one value of $\vec u$, such that ${G_m}^3(\vec u) = \vec u$.

Remark 4.1. Using the similar arguments, we can show that if $G_m^{k - 1}({u_{{Q_1}}}) < u_m^ +$, where $G_{m}(u_{m}^{+}) = {u_{{Q_1}}}$, then the system (2.3) has an order-$k$ $\left({k \ge 2} \right)$ periodic solution.

[Note]: Global asymptotic stability means: A solution $u\left (t\right)$ is called asymptotically stable if it is stable and there exists a number $\varepsilon _{0} > 0$, such that, for any other solution $v\left (t\right)$ with $\left \| u\left (t_{0}\right)-v\left (t_{0}\right)\right \| < \varepsilon _{0}$, the following holds:$\lim_{t\rightarrow \infty }\left \| u\left (t\right)-v\left (t\right)\right \| = 0$ ^[35].

5.   Simulations and discussion

In this section, we carry out some numerical simulations, illustrating our theoretical findings in previous sections and then discuss these results from a biological point of view.

Figure 4(a) is the trajectory of the order-2 periodic solution of the model. Figure 4(b), (c) is the time series of phytoplankton density and fish density under the simulated order-2 periodic solution, respectively. The existence of the order-2 periodic solution of the system indicates that after two fishing in a period of time, the population returns to the initial point, so the population shows this stable periodic change.

In Figure 5(a), we numerically illustrate our mathematical findings in Theorem 3.2, corresponding to the cases in Figure 3. In sub-figure a, we see the only order-1 periodic solution (see Figure 5(a)), and in the case of Figure 3(b), there is no order-1 periodic solution (see Figure 5(b)).

Figure 6 shows that when the weight parameters ${{a_1}}$ and ${{b_1}}$ change, the pulse and phase sets change significantly. That is to say, the comprehensive fishing measures that comprehensively consider the fish population density and the current fish growth rate are compared to the fishing measures that only consider the fish population density. The model will generate more complex pulse sets and phase sets. Below, we discuss which of the two fishing measures is more biologically appropriate.

From Figure 7, we see that, no matter what the value of $r$ is, the smaller the weight of $a_{1}$, the greater the weight of $b_{1}$, the longer the single period of the order-1 periodic solution, that is, the comprehensive control measures that take into account the fish population density and the current growth rate are longer than the cycle time of the system that only considers the fish population density. That is, the comprehensive fishing measures that comprehensively consider the fish population density and the current fish growth rate are more reasonable than the fishing measures that only consider the fish population density.

It can be seen from Figure 8 that the trajectories of different initial points, eventually. This illustrates the global asymptotic stability of the order-1 periodic solution. This provides theoretical support for the application of state-dependent feedback control in ecosystem balance.

6.   Conclusions

In this paper, we propose an integrated control model that takes into account fish population densities and their current growth rates. The control form adopted is more in line with the development law of biological populations, so that fish can be harvested in time, and more economic benefits can be obtained. We delve into the complex dynamics of the model and illustrate the effect of dynamic thresholds on system dynamics. In this paper, the existence and stability of periodic solutions of order-$1$, and the existence of periodic solutions of order-$k$ ($k \ge 1$) are proved. The correctness of the theoretical results was verified by numerical simulation, indicating that the phytoplankton and fish population densities can maintain periodic oscillations through effective control strategies, that is, the population size can be controlled within a stable range.

The results of this paper are an extension of literature ^[32], and we believe that the innovation of this model is as follows:

$1)$ This paper employs dynamic thresholds determined by fish population density and current fish growth rate. Fixed thresholds are a special case of dynamic thresholds, which are extensions of fixed thresholds. For example, when $a_{1} = 1$, $b_{1} = 0$ in model (2.3), it is the case of a fixed threshold.

$2)$ The dynamic properties of the model can be better investigated using the method of Poincar$\acute{\mbox{e}}$ map. The application of control methods and research methods make the model (2.3) have more complex dynamics than the literature ^[32] as Theorems 3.1 and 3.2.

Due to the synthesis of various factors affecting the development of the population, the control method adopted in this paper is more realistic, and the obtained research results have practical guiding significance. The analysis method proposed in this paper also has reference value in analyzing impulse models with complex phase sets or impulse sets.

Acknowledgements

The paper was supported by the National Natural Science Foundation of China (No. 11872335).

Conflict of interests

Declare that there no of interest regarding the publication of this paper.