Optimal control of a discrete-time plant–herbivore/pest model with bistability in fluctuating environments

1.   Introduction

Discrete-time population models have been developed to understand the complicated plant-herbivore/pest dynamics ^[1,2,3,4]. Based on extensive field and laboratory data on European Corn borer, Cavalieri and Kocak ^[5] proposed and studied a discrete-time model for European Corn Borer (ECB). Simulations of their proposed pest–pathogen interaction model determine conditions that cause chaotic dynamics of the pest population, warranting a cautious approach when dealing with biocontrol agents in consideration of their long-term effects. As an another example, the growth of Gypsy moth populations has caused a serious problem in the northeastern US ^[6,7]. Kang et al. ^[8] formulated a discrete-time plant–herbivore model compatible with the conservation of biomass production ^[9]. Their model has two parameters: (1) the growth rate of the plant population and (2) the damage inflicted by herbivores, which exhibits rich dynamics such as bistability and chaos ^[8]. Continuous and discrete approaches for a coupled host–parasitoid system were employed to investigate integrated pest management (IPM) control programs identifying many factors such as host–parasitoid ratios, starting densities, timing of parasitoid releases, dosages and timings of insecticide applications and levels of host-feeding and parasitism ^[10]. Other integrated pest management (IPM) has been proposed in the switched discrete host-parasitoid framework ^[11]. They searched effective switched IPM strategies by use of the economic threshold (ET) constraints varying important parameter sets in their model but in under the constant environment.

Optimal control theory has been successfully employed in various areas of applied mathematics such as economical, social, physical, biological, and epidemiological models ^[12]. Various countermeasures are incorporated in the infectious transmission dynamics using a standard continuum approach ^[13,14,15]. Optimal control problems in regulating pest population have been studied; the soybean caterpillar-predator system was developed to investigate the impacts of pest controls on the population dynamics of Lotka-Volterra multi-species models ^[16]. Furthermore, many researchers were interested in how to implement optimal strategies between chemical and biological controls in numerous discrete-time plant–pest models ^{[17,18,19,20]}. The effects of chemical and biological control have been studied ^[17] through a discrete-time host–parasitoid model. They investigated the impact of the inoculative release of parasitoids upon the host–parasitoid interactions with the host growth rate modeled by either a Beverton–Holt or a Ricker type nonlinearity. And they studied the optimal control of the variable release of parasitoids over generations to minimize the pest population and the cost of implementation. An another discrete time optimal control theory is employed to the mathematical modeling of pest control under two scenarios: biological control and the combination of pesticide and biological control ^[18]. The optimal combination of two different pest control methods was formulated and it was shown that releasing steriles is economically preferable when the infesting population is low, while if the population is high, pesticide has to be preferred ^[19]. Recently, Abraha et al. proposed a continuous mathematical model incorporating farming awareness in crop pest management, considering plant biomass and pest (healthy and infected) ^[20]. They formulated an optimal control problem incorporating two controls (pesticide and awareness campaign) and the effect of awareness programs on the control of pests were analysed.

There has been much less work on optimal control problems in the presence of the varying environment. An optimal control problem has been formulated in a gypsy moth model with periodic environments in ^[21], where they studied the effects of the Gypchek spray as optimal controls on gypsy moths. Specifically, ^[21] developed a discrete-time Nicholson–Bailey model for the tritrophic interaction of gypsy moths, tannin, and the nuclear polyhedrosis virus to explore optimal strategies for the cost-effective control of the gypsy moth population through the application of a biocontrol agent. An optimal harvesting policy has been proposed by using the exploitation of a single population modeled by time-dependent Logistic equation with periodic coefficients ^[22]. In the work of ^[23], both autonomous and nonautonomous population models are considered subject to either impulsive or continuous harvesting. Their work demonstrated that the impulsive strategy can be as good as the continuous one while it cannot outperform the autonomous model; and the impulsive strategy can be optimal in the nonautonomous model depending on their parameter values ^[23].

Motivated by the negative ecological impacts of forest pest gypsy moths and the related literature work of ^[8,21], we proposed optimal control problems of the discrete-time plant-herbivore model with two different constant control strategies (removal versus reduction strategies) to investigate how to reduce the pest population at a minimal cost. The corresponding optimal control problems has been addressed on three scenarios of bi-stability dynamics of of the plant-herbivore/pest model studied ^[8]. Our main goals are to:

● Explore two control strategies for a discrete-time plant-pest model.

● Optimal control outcomes when the model has varied bistability dynamics.

● Explore the dynamics of discrete-time models in fluctuating environments. And

● optimal control outcomes for the model when it has both bistability and fluctuating environments.

The remaining of the article is as follows. In Section 2, we revisited the plant-pest model in ^[8] and introduced two constant control strategies in the model. In Section 3, we formulated a discrete optimal control problem for the model in ^[8] and provided necessary conditions for the optimal solution to the proposed the optimal control problem. In Section 4. we provide the optimal control results of three scenarios of bistability dynamics and the corresponding cases under fluctuating environments. In the final section, we conclude our work.

2.   A plant-pest model and its constant control

Gypsy moth is a notorious forest pest in North Central United States and its outbreaks are almost periodic and cause significant damage to the infested forests. Motivated by this, Kang ^[8] developed a general discrete-time plant–herbivore (or plant–pest) model (2.1) that is modified from the well known Nicholson–Bailey model with the following assumptions :

● In the absence of the pest, the dynamics of plant follows a Ricker's map.

● The attack of pest occurs after the growth phase of plant.

where $x_n$ and $y_n$ are population densities of plant, pest at time $n$, respectively. And the parameter $r$ represents the intrinsic growth rate of plant, and $a$ denotes the attacking rate of pest to plant. The plant–pest model (2.1) has extremely complicated dynamics that include but are not limited to different bifurcations (e.g., transcript bifurcation, saddle node bifurcation, Neimark–Sacker bifurcation), bistability between boundary attractors and interior attractor, strange interior attractor and the crisis of the strange attractor.

The baseline model ^[8] is modified through the incorporation of constant control functions. Here, we introduce two strategies of constant control on pest population:

1. Removal of a constant ratio $0 < u < 1$ of current pest population:

2. Reduce the current pest population by a constant ratio $e^{-u}, u > 0$ of current pest population:

These two constant controls are expected to have different effectiveness for a same $u > 0$. Our ultimate goal is to study the optimal control problems, thus, it is necessary to understand the constant control on pest population. For this reason, we first investigate how the control parameter $u$ affect the pest population for two control strategies (2.2) and (2.3). It is easy to check that both (2.2) and (2.3) are positively invariant and bounded in $\mathbb{R}^2_+$, which leads to the following theorems:

Theorem 1 (Positive and bounded). Both $(2.2)$ and $(2.3)$ are positively invariant and bounded in $\mathbb{R}^2_+$. More specifically, System $(2.2)$ maps $\mathbb{R}^2_+$ to a compact set $\left[0, \frac{e^{r-1}}{r}\right]\times\left[0, \frac{e^{r-1}}{r}\right]$ after the first iteration while System $(2.3)$ maps $\mathbb{R}^2_+$ to a compact set $\left[0, \frac{e^{r-1}}{r}\right]\times\left[0, \frac{e^{r-1-u}}{r}\right]$ after the first iteration.

Proof. Since both (2.2) and (2.3) are continuous maps and have both $x$, $y$-axis as their invariant manifolds, therefore, we can easily deduce that both (2.2) and (2.3) are positively invariant based on continuity arguments. Notice that

Thus, the statement of Theorem 1 holds.

Notes: Theorem 2 shows that our control problem is well defined, and we are able to do further analysis in the following subsections.

2.1.   Mathematical analysis of constant controls on pest population

In this subsection, we focus on the local and global dynamical property of both System (2.2) and (2.3) with respect with $u > 0$. First, we look at the possible equilibria of (2.2) and (2.3) when $u > 0$. We can easily verify that both System (2.2) and (2.3) have the same boundary equilibria, i.e., $(0, 0)$ and $(1, 0)$, as the original System (2.1). However, the stability of $(1, 0)$ may be affected due to the control. The following theorem is a summary of the control $u$ on the possible equilibria and their stability:

Theorem 2 (Equilibria and their stability). Both System $(2.2)$ and $(2.3)$ have $(0, 0)$ and $(1, 0)$ as their boundary equilibria where $(0, 0)$ is always a saddle. For System $(2.2)$, $(1, 0)$ is locally asymptotically stable if

without interior equilibrium while it is unstable if $a > 1+u$. For System $(2.3)$, $(1, 0)$ is locally asymptotically stable if

without interior equilibrium while it is unstable if $a > e^u$. In addition, System $(2.2)$ has a unique interior equilibrium $\left(x^*, y^*\right)$ if $a > 1+u$ and $\left(x^*, y^*\right)$ is locally asymptotically stable if

While System $(2.3)$ has a unique interior equilibrium $\left(x^*, y^*\right)$ if $a > e^u$ and $\left(x^*, y^*\right)$ is locally asymptotically stable if

Proof. Let $(x^*, y^*)$ be an equilibrium of System (2.2) and (2.3), then the Jacobian matrix of (2.2) evaluated at $(x^*, y^*)$ takes the form of (2.4):

while the Jacobian matrix of (2.3) evaluated at $(x^*, y^*)$ takes the form of (2.5):

Both System (2.2) and (2.3) have the same boundary equilibria, i.e. $(0, 0)$ and $(1, 0)$. For System (2.2), we have its Jacobian matrix evaluated these boundary equilibria as

and

The eigenvalues of $J^1_{(0, 0)}$ are $\lambda_1 = e^r$ and $\lambda_2 = -u$ while the eigenvalues of $J^1_{(1, 0)}$ are $\lambda_1 = 1-r$ and $\lambda_2 = a-u$. Since $0 < u < 1$, thus, for System (2.2), $(0, 0)$ is always a saddle and $(1, 0)$ is locally asymptotically stable when $0 < r < 2$ and $0 < a < 1+u$. For System (2.3), its Jacobian matrix evaluated at $(0, 0)$ and $(1, 0)$ are represented as follows

and

The eigenvalues of $J^1_{(0, 0)}$ are $\lambda_1 = e^r$ and $\lambda_2 = 0$ while the eigenvalues of $J^1_{(1, 0)}$ are $\lambda_1 = 1-r$ and $\lambda_2 = ae^{-u}$. Since $0 < u < 1$, thus, for System (2.3), $(0, 0)$ is always a saddle and $(1, 0)$ is locally asymptotically stable when $0 < r < 2$ and $0 < a < e^u$.

Let $(x^*, y^*)$ be an interior equilibrium of System (2.2), then $(x^*, y^*)$ should satisfy the following two equations:

If $(x^*, y^*)$ be an interior equilibrium of System (2.3), then $(x^*, y^*)$ should satisfy the following two equations:

Thus according to the proof of Proposition 3.1 & 3.2 by Kang et al ^[8], we can conclude that System (2.2) has no interior equilibrium if $0 < r < 2 \text{ and } 0 < a < 1+u$ while it has a unique interior equilibrium if $a > 1+u.$ Similarly, we can conclude that System (2.3) has no interior equilibrium if $0 < r < 2 \text{ and } 0 < a < e^u$ while it has a unique interior equilibrium if $a > e^u.$

In the case that System (2.2) has an interior equilibrium $(x^*, y^*)$, its Jacobian matrix (2.4) takes the following form:

Thus, by the Jury Test (P57, ^[24]), we can conclude that the interior equilibrium $(x^*, y^*)$ of (2.2) is locally asymptotically stable if

Similarly, we can show that the interior equilibrium $(x^*, y^*)$ of (2.3) is locally asymptotically stable if

Notes: Based on the results, it seems to suggest that control $u$ can stabilize the dynamics.

Theorem 3 (Extinction and persistence). Plant species $x$ is always persistent in $\mathbb{R}^2_+$ for both System $(2.2)$ and $(2.3)$. Pest species $y$ is also persistent in $\mathbb{R}^2_+$ for System $(2.2)$ if

while pest species $y$ is persistent in $\mathbb{R}^2_+$ for System (2.3) if

In addition, all periodic orbits on the invariant manifold $y = 0$ of System $(2.2)$ are transversally stable if

while all periodic orbits on the invariant manifold $y = 0$ of System $(2.3)$ are transversally stable if

Proof. From Theorem 1, we know that System (2.2) and (2.3) map $\mathbb{R}^2_+$ to a compact set $\left[0, \frac{e^{r-1}}{r}\right]\times \left[0, \frac{e^{r-1}}{r}\right]$, $\left[0, \frac{e^{r-1}}{r}\right]\times \left[0, \frac{e^{r-1-u}}{r}\right]$, respectively. In addition, both System (2.2) and (2.3) maps $(0, y)$ to $(0, 0)$, i.e., omega limit set of $y$-axis is the fixed point $(0, 0)$, thus, the persistence of plant species $x$ is determined by the eigenvalue of $(0, 0)$ whose eigenvector points toward to $x$-axis. According to Theorem 2, we know that this eigenvalue is $\lambda_1 = e^{r} > 1$. Therefore, plant species is persistent in $\mathbb{R}^2_+$ for both System (2.2) and (2.3) based on Theorem 2.2 and its Corollary 2.3 by Hutson's (1984) ^[25] by using an average Lyapunov function $P(x, y) = x$.

This indicates that we can restrict the dynamics of both System (2.2) and (2.3) to a compact set $\left[b, \frac{e^{r-1}}{r}\right]\times \left[0, \frac{e^{r-1}}{r}\right]$ or $\left[b, \frac{e^{r-1}}{r}\right]\times \left[0, \frac{e^{r-1-u}}{r}\right]$.

If $0 < r < 2$, then both (2.2) and (2.3) become the well-known Ricker's map $x_{n+1} = x_{n}e^{r\left(1-x_{n}\right)}$ which has the global stability at $x = 1$ (Cull 1981 ^[26]). Thus, the omega limit set of $\left[b, \frac{e^{r-1}}{r}\right]\times \{0\}$ is $(1, 0)$. By applying Hutson's (1984) Theorem 2.2 ^[25] through the average Lyapunov function $P(x, y) = y$, we can conclude that the persistence of pest species $y$ in the joint System (2.2) and (2.3) is determined by the transversal stability of $(1, 0)$, i.e., the values of $\lambda_2 = a-u$ for (2.2) and $\lambda_2 = ae^{-u}$ for (2.3). Therefore, pest species $y$ is persistent for System (2.2) if $0 < r < 2, a > u+1$ while pest species $y$ is persistent for System (2.3) if $0 < r < 2, a > e^u$.

If plant species has a periodic orbit of period $N$, i.e., $\{x_1, x_2, \cdots, x_N\}$ at its single state (i.e., $y = 0$), then the transverse eigenvalue for this $N$-periodic orbit is $\prod_{n = 1}^N (a-u) x_n$ for System (2.2) and is $\prod_{n = 1}^N a e^{-u} x_n$ for System (2.3). As we discussed earlier, if $N = 1$, i.e., at $(1, 0)$, this eigenvalue becomes $(a-u)$ for System (2.2) and becomes $a e^{-u}$ for System (2.3), hence the equilibrium $(1, 0)$ is attractive in the y-direction for System (2.2) when $0 < a < 1+u$ and for System (2.3) when $0 < a < e^u$. It is known that for the Ricker map $x_{n+1} = x_{n}e^{r\left(1-x_{n}\right)}$ (i.e., plant species at its single state), the time average of $N$-periodic orbit is identical to the equilibrium $(1, 0)$ ^[27]. An application of the inequality between arithmetic and geometric means yields

implying that all the period orbits on the invariant manifold $y = 0$ are transversally stable for for System (2.2) when $0 < a < 1+u$ and for System (2.3) when $0 < a < e^u$.

Therefore, the statement of Theorem 3 holds.

Notes: Results of both Theorem 2 and Theorem 3 indicate that control $u$ can stabilize the dynamics and large control $u$ can eliminate the pest. We would like to point out that the detailed proof shown above is very similar to the proofs of Theorem 4.1 of Reference ^[8].

3.   Modeling a discrete optimal control problem

We formulate an optimal control problem for the discrete-time plant–herbivore system by incorporating a discrete time-dependent control, $(u_1, ..., u_N)$. As discussed in the previous section of the constant control cases, we also consider the following discrete time-dependent control on the pest population (the reduction control of System (2.3)) which reduces the current pest population by a time-dependent control $u_n$ in $e^{-u_n}, u_n > 0$ of current pest $y_n$:

Now, our goal is to minimize the pest population via the cost effective use of control measures over a finite number of generations $n = 1, \cdots, N$. Therefore, the discrete version of the objective functional $\mathcal {J}$ to be minimized is given by

subject to System (3.1) where ${\bf u} = (u_1, ..., u_N)$ is a discrete control and $0 \leq u_n\leq b$ for $n = 1, 2, ..., N$. Note that for each generation, the control variable, $u_n$ lies between 0 and the constant upper bound, $b$ with $0 \le b \le 1$, which represents the maximum level or the effectiveness of control measures such as eradication efforts, pesticide or biocontrol products. In the discrete version of the objective functional (3.2), $W_1$ is the weight constant on the pest and $W_2$ is the weight constant on the control that represents the relative cost of implementation of controls. The use of these definitions and notations lead to the problem of finding an optimal pair $({\bf u}^*, {\bf x}^*)$ such that $\mathcal {J}$ has a minimum at $({\bf u}^*, {\bf x}^*)$ using the vector notation ${\bf u} = \{u_n\}$ and ${\bf x} = \{x_n, y_n\}$ for $n = 1, 2, ..., N$. The optimization problem is solved by using a discrete version of Pontryagin's maximum principle. First, the Hamiltonian associated with the problem is defined as

where

Then, we have the following theorem, which describes the discrete version of Pontryagin's Maximum Principle.

Theorem 4. Given an optimal $u^*$ and the corresponding state solutions $x^*$ and $y^*$, there exist adjoint variables $\lambda^1_n$ and $\lambda^2_n$ such that the following holds for each generation $n$

with the transversality conditions

Again, for each generation, the optimal solution $u^*_n$ satisfies

where $\bar u_n = \{ u | K(u) = 0 \}$

Proof. Using the objective functional $\mathcal{J}(u)$ and System (3.1), we form the Hamiltonian (3.3). Pontryagin's Maximum Principle can be extended to one discrete equation and a system of discrete equations as given in Sections 23.1 and 23.2, respectively ^[12]. The necessary conditions are obtained in a similar manner as the continuous models when the the regularity condition is satisfied (the state equation and the objective functional are continuously differentiable with respect to all state and control variable as given p.196 and p.199 in ^[12], Theorem 8.3 in ^[28], ^[29,30]). Note that System (3.1) and $W_1 y_n + \frac{W_2}{2} u_n^2$ are continuously differentiable with respect to all state and control variable ($x_n$, $y_n$, $u_n$). In addition, the objective functional (3.2) is a convex function of the control variable $u_n$. Therefore, this convexity of the control function $u_n$ in the objective functional combined with the regularity conditions guaranteed the existence of optimal solutions.

Now, applying an extension of Pontryagin's Maximum Principle, we obtain the following adjoint equations for $n = 1, \cdots, N-1$ with the transversality conditions at $n = N$;

with the transversality conditions (which are imposed at the final time $n = N$).

Note that $\frac{\partial F_n}{\partial x_n} = e^{r(1-x_n)-a y_n} - r x_n e^{r(1-x_n)-a y_n}$, $\frac{\partial F_n}{\partial y_n} = - a x_n e^{r(1-x_n)-a y_n}$, $\frac{\partial G_n}{\partial x_n} = e^{r(1-x_n)} \left[1-e^{-a y_n}\right] e^{-u_n} - r x_n e^{r(1-x_n)} \left[1-e^{-a y_n}\right] e^{-u_n}$, and $\frac{\partial G_n}{\partial y_n} = a x_n e^{r(1-x_n)} e^{-a y_n} e^{-u_n}$.

Again, for each generation, the optimality condition can be written as

Lastly, using the necessary conditions for bounded controls, we get $u = 0$ when $\frac{\partial H_n}{\partial u_n} < 0$, $0 < u_n < b$ when $0 = \frac{\partial H_n}{\partial u_n}$, and $u = b$ when $0 < \frac{\partial H_n}{\partial u_n}$. Solving for $u^*$subject to the bounds, we obtain the characterization of the optimal control (3.6). The details of the optimality condition are found in p.73-75 in ^[12].

The standard scheme (a discrete version of the two point boundary method as given in the chapter 4 p.49-56 of ^[12]) is employed to find numerical solutions. First, the state System (3.1) is solved forward in time with initial conditions and an initial guess for the control. Second, the adjoint System with transversality conditions (3.5) is solved backward in time. Third, the optimality condition is updated and finally, the three steps above are iterated until convergence is achieved. Note that the optimality condition (3.6) is a nonlinear function of $u$ and it can be solved for $u_n^*$ by Newton method as given in the chapter 6 p. 135-163 of ^[31].

There are several important parameter values which affect optimal solutions greatly (control weight constants, control upper bounds and initial conditions, $(x_0, y_0)$). For the weight constants of the pest and the control, $W_1 = 1$ and $W_2 = 1$ are fixed throughout our simulations in the objective functional (3.2). Another important control parameter value is the control upper bound, $b$, which represents the maximum level of control effort or the effectiveness of control measures. In this optimal formulation, we assume that the best possible reduction by implementing control is $30\%$ or $50\%$ (i.e. $b = 0.3$ and $b = 0.5$ are used). We vary the control upper bound since numerical results are sensitive to this parameter. Throughout our simulations, two distinct initial conditions are chosen from the basin of attractors.

4.   Numerical simulation results

We first provide three bistability scenarios of the plant-pest dynamics as given in Table 1 with their corresponding attractors shown in Figure 1, their basin of attractions shown in Figure 2 and discrete time-series of plant-pest dynamics in Figure 3. Figure 1 illustrates three distinct characteristics of bistability of the phase planes of pests ($y$-axis) and plants ($x$-axis): The left panel of Figure 1 displays the stable interior equilibrium (a blue-colored circle) and a period two orbit on $x$-axis (two red-colored circles). The middle panel shows the period 21 attractor (blue colored circles) and chaotic on $x$-axis (red-colored circles). The right panel displays the period 6 attractor (blue colored circles) and chaotic on $x$-axis (red-colored circles).

Next, the blue points in Figure 2 are basin's attractions for the corresponding interior attractors. For example, all blue points in the first case of Figure 2 converges to the stable interior equilibrium shown in Case 1 of Figure 1. As mentioned earlier, the initial condition and the control upper bound are two critical factors to determine the optimal controls and the resulting corresponding dynamics of the plant–herbivore/pest system. Therefore, we explore the impacts of initial conditions and control upper bounds. Two distinct initial conditions are chosen from the basin of attractors; an interior point (circle) and a near the boundary point (square) as shown in Figure 2. Discrete time-series of the three bistability are shown in Figure 3; the right most panel shows plant-pest dynamics for the stable interior equilibrium, while the middle one for the period 21 attractor, and the leftmost one for the period 6 attractor.

In the following four subsections, we will provide results on the constant controls of the three bistable scenarios above, and will explore the dynamics of discrete-time models in fluctuating environments. And optimal control outcomes for the model when it has both bistability and fluctuating environments.

4.1.   Constant controls of the three bistability scenarios

In this subsection, we have carried out numerical simulations for two different models with constant control efforts (2.2) and (2.3). The baseline model of (2.1) has been modified by incorporating constant control efforts; first, the removal control effort in (2.2) and second, the reduction control effort in (2.3). Figures 4-5 present numerical simulations for the three bistability scenarios (left ones for the removal model (2.2) and right ones for the reduction model (2.3). As shown in Figure 4, the impacts of constant controls on plant-herbivore dynamics are similar in both models. For Case 1, the constant control value, $u$ = 0.3 can eliminate pests for both models (see red pests in the top panels of Figure 4). Furthermore, Case 1 of the stable interior equilibrium was changed to a period two orbit on x-axis in Figure 4 under both models (left and right panels).

On the other hand, for Case 2 and Case 3, the control value, $u$ = 0.3 is not enough to eliminate pests in both models (see red pests in the middle and bottom panels of Figure 4). Therefore, the control value should be increased to $u$ = 0.5 for successful pest management as shown in the top and bottom panels of Figure 5. Note that the resulting plant-pest dynamics were different in both models in the presence of the constant control (see the left and right panels). In addition, a period 21 attractor and a period 6 attractor were converged to chaotic on x-axis in the left and right panels.

4.2.   Optimal controls of the three bistability scenarios

We solve the discrete optimal control problem and identify the optimal controls and the resulting plant-pest dynamics under the three bi-stability cases shown in Table 1. The impacts of two distinct initial conditions on the optimal control outcomes for the plant–herbivore/pest dynamics are illustrated by using a control upper bound, $b = 0.3$. The results of the inner initial points are shown in the top panels while the ones of the boundary initial points are given in the bottom panels in Figures 6-8.

Figure 6 provides dynamics of the plant–herbivore/pest model without controls, dynamics with controls and optimal controls in the left, middle and right panels, respectively. The pest population (red) remains higher than the black plant population (black) without optimal control in the left panels. As seen in the right panel of Figure 6, Case 1 requires the maximum level only at the very beginning and rapid reduction of the control resulting in the elimination of the pests for both initial conditions (see optimal controls in the right panels). The red pest population went down to zero very quickly in the middle panels of both top and bottom. The results for Case 1 are not significantly different under two distinct initial conditions. This implies that the optimal controls and the corresponding plant-herbivore dynamics for the stable interior equilibrium case are not sensitive to initial points and relatively easier to control the pest population. Note that the right panels of top and bottom; only a maximum level of control is needed for a very short period.

Interestingly, the results shown in Figures 7 and 8 for bistability of Case 2, and Case 3, respectively, are different under two different initial conditions. The left panel of Figure 7 shows the period 21 attractor without optimal controls; the red pest population remains higher than the black plant population. Now, the left panel of Figure 8 shows the period 6 attractor without optimal controls, and similarly, the red pest population remains higher than the black plant population. For the inner initial point, the implementation of the maximum level of control during the entire period is not enough to reduce or eliminate the pest. Note that the pest-plant dynamics became a stable interior equilibrium in the middle panel of the top in Figures 7 and 8. On the other hand, for near the boundary initial points of the basin of attractions, Case 2, and Case 3 require the maximum level at the beginning and moderate reduction of the control resulting in the elimination of the pests. See optimal controls in the right panel and the resulting pest went down to zero in the middle of the bottom in Figure 7 and 8.

When the control upper bound is increased from 0.3 to 0.5, the pest can be manageable with controls. Figure 9 shows results using the upper bound 0.5 for both Case 2 and Case 3. The maximum level at the beginning and monotone decreasing optimal control function eliminates the pests. Note that the pest population (red) went down to zero in the middle panels and optimal controls in the right panels. Our numerical simulations confirm that the upper bound of the control is a key parameter to eliminate the pest. The implementation of optimal control with the upper bound 0.3 is not enough to reduce or eliminate the pest for Case 2 and Case 3 using the inner initial points. But, using the control upper bound is 0.5, the pests are effectively eliminated. This is consistent with the results of the constant control as given in Theorem 3 ($0 < a < e^u$), where $a = 1.5 < e^{0.5} \approx 1.6487$. Moreover, this is has been confirmed in the numerical results with the constant control efforts as shown in Figures 4 and 5. Even though the impacts of both constant control and optimal control (time-varying) are similar, note that constant control is implemented during the entire time while optimal control is required at the beginning only. This time-varying characteristics of optimal control highlights the importance of the effectiveness for pest control strategies.

Our optimal control results under the three baseline scenarios indicate that Case 1 (the stable interior equilibrium) is easier to eliminate the pest regardless of control upper bounds or initial conditions. However, under the inner initial points, the pest is harder to be eliminated than the case with the boundary initial points in Case 2 and Case 3 when the control upper bound is not high ($b$ = 0.3). Moreover, Case 2 (a period 21 attractor) required a higher control level or amount than Case 3 (a period 6 attractor). The maximum level of control effectiveness or the upper bound of control is increased to 0.5 (as long as $0 < a < e^u$, see Figure 9), then all three cases are not that significantly different under distinct initial conditions, which means that all pests can be manageable. The impacts of control upper bounds are summarized in Table 2 and below.

● When the bistability dynamics are simple (e.g., Case 1), the pest can be eliminated with lower control upper bounds (see Figure 6) regardless of initial conditions taken from boundary or center of the basin attractions of the interior attractor.

● Due to bistability, initial conditions are important for pest eliminations. For example, for Case 2 and 3 (see Figure 2.4 and 2.5), pest population is constant with the control upper bound being 0.3 for initial conditions taken from the center of the basin attractions but pest can be eliminated when initial conditions taken from the boundary of the basin attractions

● When control upper bounds are not large enough, constant control is the optimal control outcome for Case 2 & 3, and led to the simple dynamics of both plant and pest.

4.3.   Bi-stability of plant–herbivore dynamics with fluctuating environments

In this subsection, we explore the impacts of fluctuating environments on plant-herbivore dynamics. Recall the discrete-time plant-herbivore System 2.1 and the three bi-stability scenarios provided in the previous section. Fluctuating environments can be constructed by replacing constant parameters, $r$ and $a$ in System 2.1 with time-dependent parameter values $r_n$ and $a_n$ as given below:

In this paper, we focus on the time-dependent parameters $r_n$ and $a_n$ (the growth rate of the plant and the damage rate inflicted by the herbivore) being periodically alternating in double or triple periods. Specifically, the corresponding parameter values of $r_n$ and $a_n$ have either period 2 (double), or period 3 (triple). For instance, the fluctuating environment with period 2 can be constructed as alternating two taken from the three baseline cases: Case 12 with alternating Case 1 and Case 2 (the details are given below) or Case 21 with alternating Case 2 and Case 1, Case 13 with alternating Case 1 and Case 3 and the rest of combinations (a total six of them). We would like to point out that the order matters for the dynamical outcomes, i.e., the dynamics of Case ij and Case ji would be totally different.

Similarly, the fluctuating environment with period 3 can be constructed in the same fashion: we alter Case 1, Case 2, and Case 3 (denoted by Case 123, and the details are given below). We continue for Case 132, Case 213, and the rest of the combinations (a total of six). For instance, we can define a double period for $C_{12}$ as

For instance, we can define a triple period for $C_{123}$ as

We investigate dynamical outcomes for the double period and triple period in Figures 10-12. Pest-plant dynamics, the basin of attractions, and phase planes are illustrated. First, the left three panels of Figure 10 display plant-pest dynamics (red pest and blue plant) in 500 generation for $C_{12}$, $C_{13}$, and $C_{23}$. The right two panels present time series for $C_{123}$ and $C_{132}$. For all cases, the plant population is lower than the pest population. Note that the dynamics of fluctuating environments are different from the ones of three baseline scenarios. Certainly, alternating bistability parameters change the dynamics (as well as the basin of attractions and phase planes).

The basin of attractions for double and triple periods are displayed in Figure 11. As shown in the middle and bottom panels of Figure 11, the basin of attractions of double periods (or triple periods) become a little wider than the original three bistabiliy scenarios (top three panels).

Figure 12 provides bistability dynamics under periodic 2 (Case 12, 13, 23) and periodic 3 (Case 123,132). For double period cases, there are three different types of bi-stability dynamics (see top three panels); Case 12 has a period 2 attractor (blue-colored circles) and a period 6 attractor on the x-axis (red-colored circles). Case 13 shows a period 2 attractor (blue-colored circles) and chaotic on the x-axis (red-colored circles). Lastly, Case 23 displays chaotic attractor (blue-colored circles) and chaotic on the x-axis (red-colored circles). Next, for triple period cases, there are two different types of bi-stability dynamics (see bottom panels); Case 123 with a period 3 attractor (blue-colored circles) and chaotic on x-axis red-colored circles). Case 132 has a period 6 attractor (blue-colored circles) and chaotic on x-axis (red-colored circles). Similarly, alternating parameters change the dynamics of bistability as well. Characteristics of bistability for double and triple periods are summarized in Table 3 and 4, respectively.

Lastly, we illustrate some selected numerical simulations on the impacts of the order of period in fluctuating environment. As shown in Figure 13, the results of C 13 and C 31 are displayed in the left and right top panels, respectively, while the results of C 23 and C 32 are shown in the left and right panels, respectively. These selected examples confirmed that the order of period changes the resulting plant-pest dynamics. Therefore, further investigation should be carried out in future study.

4.4.   Optimal controls with fluctuating environments

In this subsection, we investigate the optimal controls of the corresponding plant-pest dynamics with fluctuating environments. We explore the impacts of initial conditions on the optimal controls and corresponding plant and herbivore dynamics (using two distinct initial conditions based on the basin of attractions see Figure 11). Furthermore, the impacts of control upper bound (two distinct values, 0.3 and 0.5) are carried out as in the subsection 4.1.

First, we begin with the results under fluctuating environments with double periods. Figure 14 shows the results of Case 12 using two distinct initial conditions and a control upper bound $b = 0.3$. Case 12 requires the maximum level at the very beginning and rapid reduction eliminates the pests under both two distinct initial conditions (see optimal control in the right panels). Less controls are used when we have compared with the baseline results (recall Figures 6-7). In a similar fashion, Case 13 (or Case 23) eliminates the pest with similar control efforts (see Figures 15-16). Case 23 required the most control efforts among Case 12, Case 13, and Case 23 (see optimal control in the right panels of Figures 16). This implies that the optimal controls and the corresponding plant-pest dynamics for the double period cases are not sensitive to initial points. The summary is given in Table 5.

Next, we present the results under fluctuating environments with triple periods. Figures 17 and 18 illustrate the results of Case 123 and Case 132 under two distinct initial conditions, respectively. In both figures, the top panels illustrate the results of the initial condition taken from the center of the basin attractions, the implementation of the maximum level of control during the entire time period is not able to eliminate the pest (when a control upper bound is 0.3). However, the bottom panels (under the initial condition from the boundary) show the maximum level at the very beginning only and the corresponding rapid reduction eliminates the pest. In Figure 19, as the upper bound of control is increased to 0.5, then all triple period cases are not significantly different under all fluctuating environments (all pests can be eliminated with the maximum level of control at the beginning time only). The implementation of control with the upper bound 0.3 is not enough to eliminate the pest for Case 123 and Case 132 when initial conditions are taken from the inner initial points. However, if the control upper bound is 0.5, the pests are effectively eliminated in all cases. This is also summarized in Table 5.

As a summary, we can conclude that both the initial conditions and the control upper bounds play key roles in the optimal control solutions and the corresponding plant-pest dynamics. Interestingly, for all double period cases, note that less control efforts were required to eliminate the pest population than the original bistability scenarios in the constant environment. This might indicate that the double period fluctuating environment may help control pest population under certain conditions by comparison to the original cases in the constant environment.

On the other hand, the results of triple periods (Case 123 and Case 132) are sensitive to two different initial conditions. When the initial condition is placed at the center of the basin of attraction, then it is not enough to eliminate the pest (or the control upper bound is not large enough) while the initial condition is near the boundary of the basin of attraction, it is easier to control. It seems to support that the control upper bounds become larger, the impact of the initial conditions become less significant. This is consistent with the results of Case 2 and Case 3 in the constant environment (Figures 7-9).

5.   Conclusions

We studied the dynamics of plant–herbivore (or pest) interaction models that incorporate the use of the constant removal or reduction controls, and the corresponding optimal control problems for regulating pest/herbivore populations under constant versus fluctuating environment. Our optimal control problems focus on the three scenarios when the original model ^[8] exhibits bistability. An optimal control problem for the reduction control has been formulated and the corresponding plant-pest dynamics are identified and evaluated under those three bistable scenarios. Furthermore, the dynamics of plant-pest system in the fluctuating environment and the related optimal control problems are also investigated, i.e., parameters $r$ and $a$ (the growth rate of the plant and the damage rate inflicted by pest) in the plant-pest model are periodically alternating in each year (e.g, double or triple periods).

Our results on optimal control in the constant environment showed that, due to bistability, both initial conditions and the control upper bound (effectiveness level) are two critical parameters for pest elimination especially when the dynamics are complicated. For all three bistability cases, an early maximum level of control is optimal and it can eliminate the pest when the initial condition is from the boundary of the basin of attractions. Case 1 (the stable interior equilibrium) is easiest to eliminate the pest regardless of control upper bounds or initial conditions. On the contrast, for Case 2 and Case 3, the pest can't be eliminated when the control upper bound is not high enough ($b$ = 0.3) and initial conditions are taken from the center of the basin attractions. Case 2 (a period 21 attractor) required a higher control level/amount than Case 3 (a period 6 attractor) to eliminate the pest. Nonetheless, when the maximum level of control effectiveness or the upper bound of control is increased to 0.5 (as long as $0 < a < e^u$), all three bistability cases are not that significantly different under distinct initial conditions, and pest can be eliminated after a finite time of controls.

Next, we carried out optimal control simulations in the fluctuating environment. For double period cases, an early intensive (maximum) control is optimal and the pest can be eliminated effectively, regardless of initial conditions. Case 23 (with Chaotic interior attractors) required the most control efforts among Case 12, Case 13, and Case 23. Note that less control efforts were required to eliminate the pest population than the original bistable scenarios under the constant environment. This suggests that double period fluctuating environment may be better for pest regulations under proper conditions. On the other hand, the results of triple periods (Case 123 and Case 132) are sensitive to two different initial conditions. When the initial condition is taken from the center of the basin of attraction, the control upper bound ($b$ = 0.3) is not enough to eliminate the pest. The initial condition is near the boundary of the basin of attraction, then an early short maximum level of control eliminated the pest. Thus, we can conclude that the control upper bounds become larger, the impact of the initial conditions becomes less significant.

In summary, we explore dynamics of a discrete-time model for plant-pest interactions with bistability and periodic environments. We study the corresponding optimal controls for varied bistability cases with/without periodic environments. The environment is periodic, the order of the period play an important role for dynamical outcomes. Optimal control could be easier if the model has (1) simple dynamics such as stable equilibrium dynamics in a constant environment or (2) chaotic dynamics with fluctuating environments. Due to Bistability, initial conditions are important for optimal controls. Regardless of with or without fluctuating environments, initial conditions taken from the boundary makes optimal control easier. When the dynamics has bistability with Period two, the pest seems to be managed easier in fluctuating environment than cases with constant environments. Our simulations suggest that the control upper bound (0.3) is enough for our model in period two cases. When the dynamics has bistability with Period three, the control upper bound (0.3) may not be enough when initial conditions are far from the boundary. In general, bistability and fluctuating environments make dynamics more complicated, and the related optimal control is more difficult. More studies are needed for discrete-time models with bistability and fluctuating environment.

Acknowledgements

This research of Y.K. is partially supported by NSF-DMS (Award Number 1716802 & 2052820); NSF-IOS/DMS (Award Number 1558127) and The James S. McDonnell Foundation 21st Century Science Initiative in Studying Complex Systems Scholar Award (UHC Scholar Award 220020472). This work of Sunmi Lee is supported by a National Research Foundation of Korea (2021R1A2B5B01002611) grant funded by the Korean government (MSIP).

Conflict of interest

The authors declare there is no conflict of interest.