# Optimality and Sustainability of Delayed Impulsive Harvesting

Jennifer Lawson and Elena Braverman
Dept. of Math & Stats, University of Calgary,
2500 University Dr. NW, Calgary, AB, Canada T2N1N4

(Final version)

###### Abstract

We consider a logistic differential equation subject to impulsive delayed
harvesting, where the deduction information is a function of the population
size at the time of one of the previous impulses. A close connection to the
dynamics of high-order difference equations is used to conclude that while the
inclusion of a delay in the impulsive condition does not impact the optimality
of the yield, sustainability may be highly affected and is generally delay-
dependent. Maximal and other types of yields are explored, and sharp stability
tests are obtained for the model, as well as explicit sufficient conditions.
It is also shown that persistence of the solution is not guaranteed for all
positive initial conditions, and extinction in finite time is possible, as is
illustrated in the simulations.

Keywords: Optimal harvesting; logistic equation; impulsive system; impulsive
delayed harvesting; population dynamics

AMS (MOS) subject classification: 92D25, 34A37

††journal: Communications in Nonlinear Science and Numerical Simulation

##  1 Introduction

Optimal sustainable resource management is one of the most important modern
problems, with the dynamics of harvesting models being an essential facet [1].
Harvesting can be represented in models either continuously or as only
occurring during short-time periods. Continuous harvesting can be described as
a continuous deduction term appearing in an ordinary Differential Equation
(DE) determining population dynamics. It assumes that harvesting occurs
without any interruptions, whereas impulsive harvesting corresponds to part of
the stock being removed at specific moments in time, with the duration of the
harvesting event being negligible compared to the process time. Even though
continuous harvesting may be preferable from the point of view of both
maximizing harvest and sustainability [2, 3], it is not always realistic or
easily applicable. This is why investigation of impulsive harvesting models is
important.

Impulsive DEs incorporate two parts, the DE that describes the behaviour of
the system during times of continuous dynamics, and the conditions that govern
the instantaneous change in the system at the impulsive moments. Usually, this
instantaneous change is a result of some external effect on the system, which
has a duration that is negligible compared to the overall time scale of the
process. Impulsive DEs have many practical applications such as pest control
[4], pulse vaccination strategies [5], and optimal harvesting in fisheries
[6]. For more on the theory of impulsive DEs see the monograph [7].

It is well known that including delays within a DE model of population
dynamics can lead to major changes to its behaviour, such as causing
instability, oscillations, and extinction, which are not observed in a
corresponding ordinary DE model [8]. The famous Hutchinson equation is one
such example. It can be compared to the logistic equation where the carrying
capacity is always a globally attractive equilibrium for all non-trivial
positive solutions. The inclusion of a delay in the Hutchinson equation can
cause the carrying capacity equilibrium to become unstable for certain values
of a delay. Other examples specific to harvesting models can be seen in [9,
10], where delay is incorporated into the continuous harvesting terms.

In an impulsive system, whenever control is involved, it is natural to assume
that the information available at the control implementation point is not up-
to-date, leading to delayed impulsive conditions. The incorporation of delays
in impulses goes back to the 1990s [11] with the further development of non-
instantaneous impulse theory, some of which is summarized in the monographs
[12, 13], and with recent progress reported in [14, 15, 16]. Delayed impulses
in the context of harvesting were explored in [17]. To the best of our
knowledge, there still has been no investigation of the optimal harvesting
policies and their sustainability for systems with delayed impulsive
harvesting in the literature, and we aim to fill this gap.

For logistic and other simple population models, such as Gompertz,
incorporating stochastic fluctuations or random differential equations created
an additional challenge but reflected unpredictable changes in the
environment, see [18, 19, 20] and the references therein. Impulsive harvesting
of a stochastic Gompertz model was a focus of [19]. Though logistic-type
equations are considered in most studies on optimal control, the use of other
growth rates can modify the preferred policies [21]. The close connection of
continuous models to difference equations has led to extensive study of
discrete population models [22] including harvesting [23]. The fact that
species do not naturally exist in isolation but coexist, compete or serve as
prey for others, led to extensive literature on harvesting of a single or
multiple populations in a food chain, and on optimal yields for exploited
species [24, 25, 26, 27]. Incorporating harvesting in systems of differential
or difference equations includes the case of structured populations where
selective harvesting is allowed [28, 29]. This approach is quite useful in
policy-making, for example, when only fish types within a certain size range
are eligible to take-home vs catch-and-release policies.

We consider the logistic equation with constant catch-per-unit effort
impulsive harvesting that is dependent upon delayed data. This assumes that
the information used to determine hunting or fishery quotas is based on data
for the population size or structure which was collected during one of the
previous harvesting events.

The main object of the present paper is the delayed impulsive harvesting
model, given for a fixed $k\in\mathbb{N}$ as

$\begin{cases}\frac{\displaystyle dN}{\displaystyle
dt}=rN(t)\bigg{(}1-\frac{N(t)}{K_{c}}\bigg{)},&t\neq nT,\text{
}n\in\mathbb{N}\\\ N(nT^{+})=\max\\{N(nT)-EN((n-k)T),0\\},&t=nT,\text{
}n\in\mathbb{N}\text{ }\\\
N(0^{+})=N_{0}^{+},N(0)=N_{0},...,N(-(k-1)T)=N_{-(k-1)}\\\ \end{cases}$ (1)

with prescribed initial conditions

$N_{0}^{+},N_{i}\in(0,\infty),i=-(k-1),...,0.$

In this model, the left-continuous $N(t)$ represents a size or a biomass of
the population as a function of time, $r>0$ is the intrinsic growth rate,
$K_{c}>0$ is a carrying capacity of the environment, $T>0$ is the time between
two consecutive harvesting events, $E$ is a harvesting effort and is assumed
to be $E\in(0,1)$ to avoid immediate extinction. We assume both that
restocking does not occur, and that the function $N(t)$ is left continuous. We
denote the size of the population after harvesting as $\lim\limits_{h\to
0^{+}}N(nT+h)=N(nT^{+})$, whereas the size of the population $N(t)$ before
harvesting is $N(nT)$, for the left-continuous function $N$. Note that we do
not require continuity at zero i.e. $N(0^{+})=N(0)$, nor do we assume anything
about the size of the population between $t=-iT$ and $t=-(i-1)T$. This makes
room for the possibility that the first harvesting event occurs at $t=0$, as
might be likely in a real-world setting, but does not require that this
happens. We assume $k\in\mathbb{N}$, though references to the non-delay model
with $k=0$ will also be given.

The motivation of incorporating delay in the impulsive harvesting is two-fold:

  1. 1.

Impulsive harvesting in general allows us to describe short duration
harvesting events, the effect of which makes it impossible for population
growth to remain continuously at optimal levels. It is a well developed topic,
with multiple publications appearing in the literature [2, 6, 17, 3, 30], but
not for delayed impulsive conditions. Thus, delayed impulsive models are of
theoretical interest. While the consideration of delayed impulsive systems is
not new [11], recent developments in the theory of such dynamical systems [14,
15, 16] have made its practical study feasible, as there are theoretical tools
to investigate it.

  2. 2.

Harvesting policies are dependent upon population data, and often the data
which is used is not up-to-date, leading to a delay in the harvesting term of
models. It is therefore important to be able to assess the impact of the
delay, so that managers may determine what additional modifications must be
made to the harvesting decisions. While non-delay harvesting models depict
gradual declines to extinction, they mostly failed to explain a frequently
observed phenomenon of instantaneous species disappearance due to over-
exploitation. Unlike traditional continuous or non-delayed impulsive
harvesting, delayed impulsive harvesting can describe immediate, not long-term
collapse of harvested population. Truncated models where the population or
commodity size is chosen as a maximum of the computed value or zero, are quite
common in mathematical economics and discrete dynamical systems. We found it
natural to extend this method to population dynamics in order to describe
possible extinction in finite time.

Let us also dwell on the choice of the model. Our purpose was to explore and
emphasize the effect of the harvesting delay in the impulsive form, thus we
consider the simplest autonomous logistic equations. The results will outline
the intrinsic effect of the delay in short-term harvesting.

Our goal is first, to consider the sustainability of (1) under harvesting,
which corresponds to the local asymptotic stability of a positive solution
which will be described later, and second, to explore the sustainable yield
(SY) and the maximum sustainable yield (MSY) of (1). The paper is structured
to follow this purpose. After presenting relevant definitions and auxiliary
results in Section 2, we explore stability of (1) in Section 3. All the issues
connected to SY and MSY, and relevant solutions of (1), are postponed to
Section 4. We show that, while optimality is unaffected by the magnitude of
delay, sustainability of the optimal solution is delay-dependent for $k\geq
2$. The analysis of the impact of the delay on local asymptotic stability of
the positive solution is based on the results obtained in Section 3. Finally,
Section 5 includes examples, numerical simulations, as well as discussion of
the results of the paper and possible future directions.

##  2 Preliminaries and Auxiliary Results

A solution $N^{*}(t)$ of impulsive system (1) is said to be stable if for any
$\varepsilon>0$ there exists a $\delta=\delta(\varepsilon)>0$ such that the
inequalities $|N_{j}-N^{*}(jT)|<\delta$, $j=-k+1,\dots,0$,
$|N_{0}^{+}-N^{*}(0^{+})|<\delta$ imply $|N(t)-N^{*}(t)|<\varepsilon$ for all
$t>0$. A solution of (1) is (locally) asymptotically stable if it is stable
and there exists $\eta>0$ such that
$\lim\limits_{t\to\infty}|N(t)-N^{*}(t)|=0$ for any $|N_{j}-N^{*}(jT)|<\eta$,
$j=-k+1,\dots,0$, $|N_{0}^{+}-N^{*}(0^{+})|<\eta$. A solution of (1) is
globally asymptotically stable if it is stable and
$\lim\limits_{t\to\infty}|N(t)-N^{*}(t)|=0$ for any
$N_{0}^{+}>0,N_{0},N_{-1}...,N_{-k+1}\geq 0$.

A $k$-th order difference equation has the form

$x_{n+1}=f(x_{n},x_{n-1},...,x_{n-k}),\quad n\in{\mathbb{N}},\quad n\geq k$
(2)

with the initial conditions $x_{0},...,x_{k}$.

A solution $x_{n}\equiv x^{*}$ of (2) is stable if for all $\varepsilon>0$,
there exists a $\delta>0$ such that
$\max\\{|x_{0}-x^{*}|,...,|x_{k}-x^{*}|\\}<\delta$ implies
$|x_{n}-x^{*}|<\varepsilon$ for any $n\geq k$. A solution $x^{*}$ of (2) is
(locally) asymptotically stable if it is stable and there exists $\eta>0$ such
that $\max\\{|x_{0}-x^{*}|,...,|x_{k}-x^{*}|\\}<\eta$ implies
$\lim\limits_{n\to\infty}|x_{n}-x^{*}|=0$. A solution $x^{*}$ of (2) is
globally asymptotically stable if it is stable and
$\lim\limits_{n\to\infty}|x_{n}-x^{*}|=0$ for any $x_{0},...,x_{k}$.

If harvesting is restricted to only the surplus production of a population,
then theoretically, harvesting should be able to continue indefinitely without
drastically altering the stock levels. The idea behind the Maximum Yield (MY)
is that it corresponds to an optimal solution of (1) such that the yield will
not be exceeded by any other solution, and that the optimal solution is
$T$-periodic leading to a constant yield over a time period. A Maximum
Sustainable Yield (MSY) is a MY where the optimal solution of (1) is (at least
locally) asymptotically stable.

In [30], the authors considered an MSY for (1) with $k=0$. The results are
summarized in the following.

######  Lemma 2.1.

[30] Consider (1) for $k=0$

$\begin{cases}\displaystyle\frac{dN}{dt}=rN(t)\bigg{(}1-\frac{N(t)}{K_{c}}\bigg{)},\text{
}t\neq nT,\quad n\in\mathbb{N}\\\ N(nT^{+})=\max\\{(1-E)N(nT),0\\},\text{
}t=nT,\quad n\in\mathbb{N}\\\ N(0)=N_{0}.\end{cases}$ (3)

Then the optimal harvesting effort is

$E_{opt}=1-e^{-rT/2},$ (4)

and the MSY is given by

$MSY=\frac{K_{c}(e^{rT/2}-1)}{T(e^{rT/2}+1)}$ (5)

The optimal positive periodic solution $N^{*}(t)$ of (3) corresponding to the
MSY and $E_{opt}$ satisfies

$N^{*}(nT^{+})=\frac{K_{c}}{e^{rT/2}+1}$ (6)

and is globally asymptotically stable.

In [30] analysis of a non-delayed impulsive model (3) is reduced to a
nonlinear difference equation of the first order. We also intensively exploit
connection between difference and impulsive equations. When a delay is
incorporated in impulsive condition, the difference equation becomes higher
order. We recall that for difference equations, the roots of the
characteristic equation of an associated linearized model should lie inside
the unit circle for local asymptotic stability, in contrast to differential
equations where the real parts of the roots have to be negative. Some
auxiliary results regarding difference equations are listed below.

######  Lemma 2.2 ([31]).

The roots of the characteristic equation

$p(\lambda)=\lambda^{2}-p_{0}\lambda+p_{1},\quad p_{0}>0,\quad p_{1}>0$

lie inside the unit circle if and only if $p_{0}-1<p_{1}<1$.

The result of [32, Theorem 1.1.1, Part f, P. 7] describes conditions when a
root of a quadratic equation lies on the boundary of the unit circle.

######  Lemma 2.3 ([32]).

A necessary and sufficient condition for a root of the characteristic equation

$\lambda^{2}-p_{0}\lambda+p_{1}=0$

with $p_{0},p_{1}\in{\mathbb{R}}$ to have a root satisfying $|\lambda|=1$ is
that either

$|p_{0}|=|1+p_{1}|$

or

$p_{1}=1\text{ and }|p_{0}|\leq 2.$

Lemma 2.4 is cited from [31, Theorem 5.10, P. 253].

######  Lemma 2.4 ([31]).

If $\displaystyle\sum_{i=0}^{k}|p_{i}|<1$ then the zero solution of the
difference equation

$x_{n+k+1}+p_{0}x_{n+k}+p_{1}x_{n+k-1}+...+p_{k}x_{n}=0$

is asymptotically stable.

The following result can be found in [31, Theorem 5.3, P. 249].

######  Lemma 2.5 ([31]).

Let $p_{0}>0$, $p_{k}\in{\mathbb{R}}$ be arbitrary, and $k\in{\mathbb{N}}$.
The zero solution of the equation

$x_{n+1}-p_{0}x_{n}+p_{k}x_{n-k}=0$ (7)

is asymptotically stable if and only if $|p_{0}|<(k+1)/k$ and
(i)
$\displaystyle|p_{0}|-1<p_{k}<(p_{0}^{2}+1-2|p_{0}|\cos(\theta^{*}))^{1/2}$ if
$k$ is odd
or
(ii) $|p_{k}-p_{0}|<1$ and
$\displaystyle|p_{k}|<(p_{0}^{2}+1-2|p_{0}|\cos(\theta^{*}))^{1/2}$ if $k$ is
even,
where $\theta^{*}$ is the solution of the equation

$\frac{\sin(k\theta)}{\sin((k+1)\theta)}=\frac{1}{|p_{0}|},\quad\theta\in\left(0,\frac{\pi}{k+1}\right).$
(8)

However, we do not need the general form of Lemma 2.5, since in our model
$0<p_{k}<p_{0}$. Then the left inequality in both (i) and (ii) becomes
$p_{0}<p_{k}+1$, the right inequalities coincide.

######  Corollary 2.5.1.

Let $0<p_{k}<p_{0}$. Then (7) is asymptotically stable if and only if the
following two inequalities hold:

$p_{0}<\min\left\\{p_{k}+1,\frac{k+1}{k}\right\\},\quad
p_{k}<\sqrt{p_{0}^{2}+1-2p_{0}\cos(\theta^{*})}\,,$ (9)

where $\theta^{*}$ is the solution of (8).

Lemma 2.6 is cited from [31, Theorem 5.2, P. 248].

######  Lemma 2.6 ([31]).

Let $q\in(0,2)$. The zero solution of the equation

$x_{n+1}=x_{n}-qx_{n-k}$

is asymptotically stable for $k=1$. It is asymptotically stable for $k\geq 2$
if and only if in addition

$q<2\cos\bigg{(}\frac{k\pi}{2k+1}\bigg{)}.$

##  3 Stability with Delay Impulsive Harvesting

We start with reducing the dynamics of the differential equation with delayed
impulsive harvesting to a difference equation.

######  Lemma 3.1.

The solution of (1) on the interval $t\in(nT,(n+1)T)$, $n=0,1,\dots$ is

$N(t)=\frac{K_{c}e^{r(t-nT)}N(nT^{+})}{K_{c}+N(nT^{+})(e^{r(t-nT)}-1)}\,.$
(10)

The solution of (1) with $N(nT^{+})=x_{n}$ at the points just after harvesting
at $t=nT$ satisfies the difference equation

$x_{n+1}=\max\bigg{\\{}\frac{K_{c}x_{n}e^{rT}}{K_{c}+x_{n}(e^{rT}-1)}-E\frac{K_{c}x_{n-k}e^{rT}}{K_{c}+x_{n-k}(e^{rT}-1)},0\bigg{\\}}=\max\bigg{\\{}f\left(x_{n},x_{n-k}\right),0\bigg{\\}},$
(11)

where $n\geq k$. In addition $x_{0}=N(0^{+})$, and $x_{1},...,x_{k}$ satisfy
the relation

$x_{i+1}=\max\bigg{\\{}\frac{K_{c}x_{i}e^{rT}}{K_{c}+x_{i}(e^{rT}-1)}-EN((i-k+1)T),0\bigg{\\}}$
(12)

where $N((i-k+1)T)$ are given by the initial conditions for $i=0,...,k-1$.

###### Proof.

Let $N(nT^{+})=x_{n}$ be the size of the population after a harvesting event.
For $t\in(nT,(n+1)T)$, the impulsive model is non-delayed, and the solution of
the differential equation exists, is monotone on the interval and is given by

$N(t)=\frac{K_{c}e^{r(t-nT)}N(nT^{+})}{N(nT^{+})(e^{r(t-nT)}-1)+K_{c}},\quad
t\in(nT,(n+1)T).$

The size of $N(t)$ at the end of the $n$-th time period before harvesting is

$N((n+1)T)=\frac{\displaystyle K_{c}N(nT^{+})e^{rT}}{\displaystyle
K_{c}+N(nT^{+})(e^{rT}-1)}.$

Using the value of $N(t)$ before harvesting as a function of $N(nT^{+})$,
combined with the definition of the impulsive condition in (1), leads to

$N((n+1)T^{+})=\max\left\\{\frac{\displaystyle
K_{c}N(nT^{+})e^{rT}}{\displaystyle
K_{c}+N(nT^{+})(e^{rT}-1)}-EN((n+1-k)T),0\right\\}.$

For $n=0,...,k-1$, $N((n+1-k)T)$ are given by the initial conditions, leading
to relation (12). For $n\geq k$ we once again treat $N((n+1-k)T)$ as a
function of $N((n-k)T^{+})$ and obtain

$N((n+1)T^{+})=\max\left\\{\frac{\displaystyle
K_{c}N(nT^{+})e^{rT}}{\displaystyle
K_{c}+N(nT^{+})(e^{rT}-1)}-E\frac{\displaystyle
K_{c}N((n-k)T^{+})e^{rT}}{\displaystyle
K_{c}+N((n-k)T^{+})(e^{rT}-1)},0\right\\}.$

Then, since $x_{n}=N(nT^{+})$ the difference equation (11) is obtained. ∎

After the reduction to a difference equation, let us justify that stability of
(1) can be deduced from that of (11).

######  Lemma 3.2.

The point $x^{*}$ is a locally asymptotically stable equilibrium of the
difference equation (11) if and only if the solution to (1) with
$N^{*}(nT^{+})=x^{*}$ is locally asymptotically stable.

###### Proof.

Evidently local asymptotic stability of the solution to (1) with
$N(nT^{+})=x^{*}$ implies local asymptotic stability of the equilibrium
$x^{*}$ of (11).

Further, let us assume that the solution $x^{*}$ of (11) is stable. Note that
the function

$F(x,a):=\frac{K_{c}ax}{K_{c}+x(a-1)}$ (13)

for a fixed $a=e^{rs}>1$ and any non-negative $x$, has the derivative in $x$

$\frac{\partial F\left(x,e^{rs}\right)}{\partial
x}=\frac{K_{c}^{2}e^{rs}}{(K_{c}+x(e^{rs}-1))^{2}}>0,\quad\left|\frac{\partial
F\left(x,e^{rs}\right)}{\partial x}\right|\leq e^{rs}\leq e^{rT}.$

If the equilibrium $x^{*}$ of (11) is stable, for any $\varepsilon>0$ there
exists $\delta>0$ such that, once all $|x_{j}-x^{*}|<\delta$, $j=0,\dots,k$,
we get $\displaystyle|x_{n}-x^{*}|<\varepsilon e^{-rT}$. The solution $N^{*}$
corresponding to $x^{*}$ on the interval $(nT,(n+1)T)$ is
$N^{*}(nT^{+})=x^{*}$ with

$N^{*}(t)=\frac{K_{c}e^{r(t-nT)}x^{*}}{x^{*}(e^{r(t-nT)}-1)+K_{c}}\,.$

Following (10), we note that on $(nT,(n+1)T)$ there exist $\zeta$ between
$x_{n}$ and $x^{*}$, $s\in[0,T]$ such that

$\displaystyle|N(t)-N^{*}(t)|$
$\displaystyle=\left|\frac{K_{c}e^{r(t-nT)}x_{n}}{x_{n}(e^{r(t-nT)}-1)+K_{c}}-\frac{K_{c}e^{r(t-nT)}x^{*}}{x^{*}(e^{r(t-nT)}-1)+K_{c}}\right|$
$\displaystyle=\left|\frac{\partial F(\zeta,^{rs})}{\partial
x}\right||x_{n}-x^{*}|\leq e^{rT}\varepsilon e^{-rT}=\varepsilon,$

thus $N^{*}$ is stable.

If the solution $x^{*}$ of (11) be locally asymptotically stable, and $n_{0}$
is such a number that $\displaystyle|x_{n}-x^{*}|<\varepsilon e^{-rT}$ for
$n\geq n_{0}$, as above,

$|N(t)-N^{*}(t)|\leq e^{rT}\varepsilon e^{-rT}=\varepsilon,\quad t\geq
nT,\quad n\geq n_{0},$

thus $N^{*}$ is both stable and attractive, and therefore is locally
asymptotically stable, which concludes the proof. ∎

Next, let us describe solution bounds for a harvested population.

######  Lemma 3.3.

Let $E\in(0,1)$. Then for any non-negative initial values there exists
$n_{0}\geq k$ such that the solution $x_{n}$ to (11) is in $[0,K_{c}]$ for
$n\geq n_{0}$.

###### Proof.

Let us note that if for some $n\geq k$, $x_{n}\leq K_{c}$, we also have
$x_{n+1}\leq K_{c}$. Really, from monotone increasing of $F\left(x,a\right)$
defined in (13) in $x$ and $F\left(K_{c},a\right)=K_{c}$ for any $a>1$, we get
$x_{n+1}\leq F\left(x_{n},e^{rT}\right)\leq K_{c}$. By induction, all
$x_{j}\leq K_{c}$, $j\geq n$. Also, $F\left(x,a\right)>x$ for $x\in(0,K_{c})$
and $F\left(x,a\right)<x$ for $x>K_{c}$.

Thus, we only have to exclude the case $x_{n}>K_{c}$ for any $n\geq k$. We
have $x_{n},x_{n-k}>K_{c}$, $F\left(x_{n},e^{rT}\right)<x_{n}$,
$F\left(x_{n-k},e^{rT}\right)>K_{c}$ and

$x_{n+1}=F\left(x_{n},e^{rT}\right)-EF\left(x_{n-k},e^{rT}\right)<x_{n}-EK_{c},$

therefore $x_{n+1}-x_{n}<-EK_{c}$, and after $\displaystyle
j=\left\lfloor\frac{x_{n}-K_{c}}{EK_{c}}\right\rfloor+1$ steps, where $\lfloor
y\rfloor$ is the integer part of $y$, we arrive at $x_{n+j}\leq K_{c}$, which
concludes the proof. ∎

Difference equation (11) has the trivial solution $x^{*}=0$, and when
$rT>-\ln(1-E)$ it has a single positive equilibrium

$x^{*}=\frac{((1-E)e^{rT}-1)K_{c}}{e^{rT}-1}\,.$ (14)

If $rT\leq-\ln(1-E)$ then only the trivial equilibrium exists, and as Lemma
3.4 states, all solutions of (11), and hence of (1), will inevitably go to
extinction.

######  Lemma 3.4.

If

$rT\leq-\ln(1-E),$ (15)

all solutions of (11) tend to zero.

###### Proof.

By Lemma 3.3, we can consider $x_{n}\in[0,K_{c}]$ for $n$ large enough. If for
some $n\geq k$, $x_{n}=0$, we have $x_{j}=0$ for any $j\geq n$ in (11), so we
restrict ourselves to only considering positive sequences $\\{x_{n}\\}$.
Further, let $\\{x_{n}\\}$ be an eventually monotone sequence, then it has a
limit $d$. If $d=0$, the sequence converges to zero; if $d>0$, we let
$n\to\infty$ in

$x_{n+1}=\frac{K_{c}x_{n}e^{rT}}{K_{c}+x_{n}(e^{rT}-1)}-E\frac{K_{c}x_{n-k}e^{rT}}{K_{c}+x_{n-k}(e^{rT}-1)}$

and get that $d$ is a positive equilibrium solution of (11), which is a
contradiction. Thus, we have only to consider sequences $\\{x_{n}\\}$ that are
neither eventually non-increasing nor eventually non-decreasing.

Before we handle this case, let us recall that the function $F(x,a)$ defined
in (13) is strictly increasing in both $x$ and $a$ for $K_{c}>0$ and $a>1$
(here $a=e^{rT}>1$ for $rT>0$).

Let $k=1$, since we are only considering sequences that are not non-
decreasing, this implies that there exists some $n$ such that $x_{n}<x_{n-1}$.
Then,

$\displaystyle x_{n+1}$
$\displaystyle=\frac{K_{c}x_{n}e^{rT}}{K_{c}+x_{n}(e^{rT}-1)}-E\frac{K_{c}x_{n-1}e^{rT}}{K_{c}+x_{n-1}(e^{rT}-1)}<\frac{(1-E)K_{c}x_{n}e^{rT}}{K_{c}+x_{n}(e^{rT}-1)}$
$\displaystyle\leq\frac{(1-E)K_{c}x_{n}\frac{1}{1-E}}{K_{c}+x_{n-k}(\frac{1}{1-E}-1)}=\frac{K_{c}x_{n}}{K_{c}+x_{n}\frac{E}{1-E}}<x_{n},$

and by induction we get that $\\{x_{j}\\}$ is a monotonically decreasing
sequence starting with $j=n$, and thus it converges to zero, as justified
above.

Next, consider $k\geq 2$. If there exists $n$ such that
$x_{n}=\min\\{x_{n-k},x_{n-k+1},\dots,x_{n-1},x_{n}\\}$, then since
$F(x_{n},e^{rT})\leq F(x_{n-k},e^{rT})$ and by (15), we get

$\displaystyle x_{n+1}$
$\displaystyle=\frac{K_{c}x_{n}e^{rT}}{K_{c}+x_{n}(e^{rT}-1)}-E\frac{K_{c}x_{n-k}e^{rT}}{K_{c}+x_{n-k}(e^{rT}-1)}\leq\frac{(1-E)K_{c}x_{n}e^{rT}}{K_{c}+x_{n}(e^{rT}-1)}$
$\displaystyle\leq\frac{(1-E)K_{c}x_{n}\frac{1}{1-E}}{K_{c}+x_{n-k}(\frac{1}{1-E}-1)}=\frac{K_{c}x_{n}}{K_{c}+x_{n}\frac{E}{1-E}}<x_{n}$

as above. Thus $x_{n+1}<x_{n}$ and
$x_{n+1}=\min\\{x_{n-k+1},x_{n-k+2},\dots,x_{n},x_{n+1}\\}$, which yields that
$x_{n+2}<x_{n+1}$, and by induction $\\{x_{n}\\}$ is monotonically decreasing
and thus converges to zero.

Let us, finally, consider the case when
$x_{n}>\min\\{x_{n-k},x_{n-k+1},\dots,x_{n-1}\\}$ for any $n\geq n_{0}>k$.
Then,

$m:=\liminf_{n\to\infty}x_{n}\geq\min\\{x_{n_{0}-k},...,x_{n_{0}}\\}.$ (16)

Introducing the functions $H$ and $h$ through $F$ as defined in (13) and using
inequality (15) leading to $\displaystyle e^{rT}\leq\frac{1}{1-E}$, we obtain

$H(x):=x-(1-E)F\left(x,e^{rT},\right)\geq
x-(1-E)F\left(x,\frac{1}{1-E}\right)=x-\frac{K_{c}x}{K_{c}+x\frac{E}{1-E}}=:h(x).$
(17)

Note that

$h(x)=x\left(1-\frac{K_{c}}{K_{c}+\frac{E}{1-E}x}\right)$

is monotone increasing for $x\in[0,K_{c}]$ (as a product of two non-negative
increasing functions) from $h(0)=0$ to $h(K_{c})=EK_{c}>0$. Let us choose
$\varepsilon>0$ such that $h(x)>\varepsilon$ for $x\in[\frac{m}{2},K_{c}]$.
From (17) we get $H(x)>\varepsilon$, $x\in[\frac{m}{2},K_{c}]$ as well. Define
a positive $\displaystyle\delta\in\left(0,\frac{\varepsilon}{4}\right)$
satisfying

$\delta<\min\left\\{\frac{m}{2},K_{c}-\frac{m}{2}\right\\}$

such that for any $x,y\in[0,K_{c}]$, the inequality $|x-y|\leq\delta$ implies

$\left|F\left(x,e^{rT}\right)-F\left(y,e^{rT}\right)\right|\leq\frac{\varepsilon}{2}.$

Such $\delta>0$ exists, as $F(x,e^{rT})$ defined in (13) is continuous and
thus is uniformly continuous for $x\in[0,K_{c}]$.

By the definition of $m$ in (16) for any $\delta>0$ there is $n_{1}\geq
n_{0}+k$ such that $x_{n}>m-\frac{\delta}{2}$ for any $n\geq n_{1}-k$. By
definition of $\liminf$, there is $n>n_{1}$ such that
$x_{n}<m+\frac{\delta}{2}$. Also, by the choice of $n_{1}$, we have
$x_{n-k}>m-\frac{\delta}{2}$ and thus $|x_{n}-x_{n-k}|<\delta$. Further,

$\displaystyle x_{n+1}$
$\displaystyle=F\left(x_{n},e^{rT}\right)-EF\left(x_{n-k},e^{rT}\right)$
$\displaystyle=F\left(x_{n},e^{rT},\right)-EF\left(x_{n},e^{rT}\right)+E\bigg{[}F\left(x_{n},e^{rT}\right)-F\left(x_{n-k},e^{rT}\right)\bigg{]}$
$\displaystyle=x_{n}-H(x_{n})+E\bigg{[}F\left(x_{n},e^{rT}\right)-F\left(x_{n-k},e^{rT}\right)\bigg{]}$
$\displaystyle\leq
x_{n}-h(x_{n})+E\bigg{|}F\left(e^{rT},x_{n}\right)-F\left(x_{n-k},e^{rT}\right)\bigg{|}$
$\displaystyle<x_{n}-\varepsilon+E\frac{\varepsilon}{2}<x_{n}-\varepsilon+\frac{\varepsilon}{2}=x_{n}-\frac{\varepsilon}{2}<x_{n}-2\delta\leq
m+\delta-2\delta=m-\delta<m-\frac{\delta}{2},$

which contradicts to the assumption that $x_{n}>m-\delta/2$ for any $n\geq
n_{1}-k$. Thus the scenario described by (16) is impossible, and the solution
either coincides with zero, starting at some $x_{j}$, or tends to zero. ∎

######  Theorem 3.5.

If (15) holds, all solutions of (1) tend to zero.

###### Proof.

By Lemma 3.4, if (15) holds then all solutions of the difference equation tend
to zero. This implies that $N(nT^{+})\to 0$, and by (10) the solution of (1)
$N(t)$ on $(nT,(n+1)T)$ also tends to zero. Therefore all solutions of (1)
tend to zero. ∎

By Theorem 3.5, it is sufficient to consider only the case

$E<1-e^{-rT}.$ (18)

If this condition is not satisfied, all solutions tend to zero. Everywhere
below we assume that (18) holds.

Next, let us focus on the behaviour of the difference equation when the
harvesting is delayed by a single time period $k=1$. This will lead to the
second-order difference equation, and allow us to apply necessary and
sufficient results such as Lemma 2.2 and 2.3 to obtain explicit conditions for
stability of the difference equation.

######  Lemma 3.6.

Let $k=1$. If $E\in(0,1/2]$ then there exists a positive equilibrium of (11)
which is locally asymptotically stable.

For $E\in(1/2,1)$, if

$rT<-\ln\bigg{(}\frac{(1-E)^{2}}{E}\bigg{)}$ (19)

then the positive equilibrium of (11) is locally unstable, while if

$rT>-\ln\bigg{(}\frac{(1-E)^{2}}{E}\bigg{)},$ (20)

the positive equilibrium of (11) is locally asymptotically stable.

###### Proof.

As (18) holds, there exists a unique positive equilibrium $x^{*}$.

Let $k=1$, then (11) has the form $x_{n+1}=\max\\{f(x_{n},x_{n-1}),0\\}$, and
the linearized equation around $x^{*}$ is

$u_{n+1}=p_{0}u_{n}-p_{1}u_{n-1},$

where

$\displaystyle p_{0}$ $\displaystyle=\frac{\partial f}{\partial
x_{n}}(x^{*},x^{*})=\frac{K^{2}e^{rT}}{(K+x_{n}(e^{rT}-1))^{2}}\bigg{|}_{(x^{*},x^{*})}=\frac{e^{-rT}}{(1-E)^{2}}\,,$
$\displaystyle p_{1}$ $\displaystyle=-\frac{\partial f}{\partial
x_{n-1}}(x^{*},x^{*})=E\frac{K^{2}e^{rT}}{(K+x_{n-1}(e^{rT}-1))^{2}}\bigg{|}_{(x^{*},x^{*})}=\frac{Ee^{-rT}}{(1-E)^{2}}.$

The characteristic equation of the linearized equation is
$\lambda^{2}-p_{0}\lambda+p_{1}=0$. By Lemma 2.2, for the roots of the
characteristic equation to lie inside the unit disc, we only require $p_{1}<1$
(the left inequality $p_{0}-1<p_{1}$ is automatically satisfied for all
$rT>-\ln(1-E)$), the inequality has the form

$rT>-\ln\bigg{(}\frac{(1-E)^{2}}{E}\bigg{)},$

which coincides with (20).

If $E\in(0,1/2]$ then $-\ln((1-E)^{2}/E)\leq-\ln(1-E)$, and therefore by (18),
the positive equilibrium $x^{*}$ is locally asymptotically stable, as (18)
implies (20).

If $E\in(1/2,1)$, (20) implies that the positive equilibrium $x^{*}$ exists
and is locally asymptotically stable.

When $E\in(1/2,1)$, if $-\ln(1-E)<rT\leq-\ln((1-E)^{2}/E)$, the positive
equilibrium exists and the roots of the characteristic equation satisfy
$\max\\{|\lambda_{1}|,|\lambda_{2}|\\}\geq 1$. By Lemma 2.3, a root of the
characteristic equation satisfies $|\lambda|=1$ if and only if
$rT=-\ln((1-E)^{2}/E)$ on this interval. Really, in this case $p_{1}=1$ and
$0<p_{0}=\frac{1}{E}<2$. Therefore for $-\ln(1-E)<rT<-\ln((1-E)^{2}/E)$, which
corresponds to (19), we get $\max\\{|\lambda_{1}|,|\lambda_{2}|\\}>1$ and by
linearization the equilibrium is unstable.

∎

Lemmata 3.6 and 3.2 immediately imply

######  Theorem 3.7.

Let k=1. If $rT>-\ln\bigg{(}\frac{(1-E)^{2}}{E}\bigg{)}$, then there exists a
unique positive periodic solution $N^{*}(t)$ of (1) with

$N^{*}(nT^{+})=x^{*}=\frac{((1-E)e^{rT}-1)K_{c}}{e^{rT}-1}\,,$ (21)

and this solution is locally asymptotically stable.

While Theorem 3.7 is similar in many ways to the result of [30] cited in Lemma
2.1, here we do not observe attractivity of the solution for all initial
conditions, implying that the solution is attractive but not globally
attractive. Even if the equilibrium is locally asymptotically stable, there
are sets of initial values leading to its instability, as further numerical
examples illustrate.

######  Remark 1.

Even for $k=1$ and any $x_{0}$, there is a domain of initial values $x_{1}$
guaranteeing immediate extinction. Really, as $\displaystyle F(x,e^{rT})$
defined in (13) is strictly increasing in $x$, $F(0,e^{rT})=0$, there are
values of $x_{1}<x_{0}$ such that

$F\left(x_{1},e^{rT}\right)=\frac{K_{c}x_{1}e^{rT}}{K_{c}+x_{1}(e^{rT}-1)}\leq\frac{EK_{c}x_{0}e^{rT}}{K_{c}+x_{0}(e^{rT}-1)}=EF\left(x_{0},e^{rT}\right),$

leading to $x_{2}=x_{3}=\dots=0$.

The following results extends Lemma 3.6 to all $k\in{\mathbb{N}}$, however,
for a given $k$, we have implicit, but easily verifiable conditions connecting
$rT$ with $E$.

######  Lemma 3.8.

The positive equilibrium $x^{*}$ of difference equation (11) is asymptotically
stable if and only if both inequalities hold

$\frac{e^{-rT}}{(1-E)^{2}}<\frac{k+1}{k},\quad\cos(\theta^{*})<\frac{e^{-rT}(1+E)}{2(1-E)}+\frac{(1-E)^{2}e^{rT}}{2}\,,$
(22)

where $\theta^{*}$ is a solution of

$\frac{\sin(k\theta)}{\sin((k+1)\theta)}=(1-E)^{2}e^{rT},\quad\theta\in\left(0,\frac{\pi}{k+1}\right).$
(23)

###### Proof.

Linearizing the difference equation $x_{n+1}=\max\\{f(x_{n},x_{n-k}),0\\}$
around the positive equilibrium $x^{*}$ given in (14) with

$f(x_{n},x_{n-k})=\frac{K_{c}x_{n}e^{rT}}{K_{c}+x_{n}(e^{rT}-1)}-E\frac{K_{c}x_{n-k}e^{rT}}{K_{c}+x_{n-k}(e^{rT}-1)},$

we get, similarly to the proof of Lemma 3.6,

$u_{n+1}=p_{0}u_{n}-p_{k}u_{n-k},\quad p_{0}=\frac{e^{-rT}}{(1-E)^{2}},\quad
p_{k}=\frac{Ee^{-rT}}{(1-E)^{2}}.$ (24)

By Corollary 2.5.1, the zero solution of the linearized equation is
asymptotically stable if and only if the inequalities in (9) hold, or
equivalently,

$\frac{e^{-rT}}{(1-E)^{2}}<\frac{k+1}{k},\quad\frac{Ee^{-rT}}{(1-E)^{2}}<\sqrt{\frac{e^{-2rT}}{(1-E)^{4}}+1-2\frac{e^{-rT}}{(1-E)^{2}}\cos(\theta^{*})}\,.$

Note that $p_{0}<p_{k}+1$, or $\displaystyle e^{-rT}<Ee^{-rT}+(1-E)^{2}$ is
equivalent to $\displaystyle e^{-rT}<1-E$, which is satisfied due to (18). The
first inequality is the same as in (22), while computing the squares in the
second one gives

$\frac{E^{2}e^{-2rT}}{(1-E)^{4}}<\frac{e^{-2rT}}{(1-E)^{4}}+1-\frac{2e^{-rT}}{(1-E)^{2}}\cos(\theta^{*}).$

After rearranging, the desired result is acquired. ∎

Applied to (1), Lemmata 3.8 and 3.2 immediately imply a sharp local asymptotic
stability result.

######  Theorem 3.9.

Let $E\in(0,1)$ and (18) be satisfied. Then there exists a unique positive
periodic solution $N^{*}(t)$ of (1) given by (21) which is asymptotically
stable if and only if both inequalities in (22) hold, where $\theta^{*}$ is a
solution of (23).

The following stability condition is delay-independent; moreover, it applies
if a constant $k$ in (11) is replaced with $k=k(n)$. It is also only
sufficient, meaning there may still exist values of $rT<\ln(E+1)-2\ln(1-E)$
such that the equilibrium is locally asymptotically stable.

######  Theorem 3.10.

If

$rT>\ln(E+1)-2\ln(1-E)$ (25)

then the positive periodic solution $N^{*}(t)$ of (1) exists and corresponds
to (21), and this solution is locally asymptotically stable.

###### Proof.

If (25) holds, then so does (18), and the positive equilibrium exists.
Reducing (1) to difference equation (11), we once again linearize (11) around
$x^{*}$ and obtain (24), where by (25),

$|p_{0}|+|p_{k}|=\frac{1+E}{(1-E)^{2}}e^{-rT}<1.$

By Lemma 2.4, the equilibrium $x^{*}$, and by Lemma 3.2 the positive periodic
solution $N^{*}(t)$ of (1) is locally asymptotically stable, once (25) holds.
∎

##  4 MSY with Delay Impulsive Harvesting

Next, we proceed with the analysis of a maximum yield (MY) and a maximum
sustainable yield (MSY). We recall that a yield is said to be sustainable if
it corresponds to a solution that is at least locally asymptotically stable.

######  Lemma 4.1.

The yield of (1) is associated to the solution $N^{*}(nT^{+})=x^{*}$ of (11)
and is given by

$Y(E)=\frac{K_{c}E}{(1-E)T}\bigg{(}\frac{(1-E)e^{rT}-1}{e^{rT}-1}\bigg{)}.$
(26)

This yield is a function of the harvesting effort $E$, it is an increasing
function for $E\in(0,E_{opt})$ and is decreasing for $E\in(E_{opt},1)$.

###### Proof.

For an optimal $T$-periodic solution of (1), we get
$N^{*}(nT^{+})=N^{*}((n+1)T^{+})=x^{*}$, where (21) holds. Then the associated
yield is

$Y(E)=\frac{EN^{*}(nT)}{T}=\frac{E}{1-E}\cdot\frac{N^{*}(nT^{+})}{T}=\frac{E}{1-E}\cdot\leavevmode\nobreak\
\frac{K_{c}}{T}\bigg{(}\frac{(1-E)e^{rT}-1}{e^{rT}-1}\bigg{)}.$

Its derivative in $E$,

$Y^{\prime}(E)=\frac{K_{c}}{T(e^{rT}-1)}\frac{d}{dE}\bigg{[}Ee^{rT}+1-\frac{1}{1-E}\bigg{]}=\frac{K_{c}}{T(e^{rT}-1)}\bigg{[}e^{rT}-\frac{1}{(1-E)^{2}}\bigg{]}.$

satisfies $Y^{\prime}(E)>0$ for $(1-E)^{2}>e^{-rT}$, which is equivalent to
$E\in(0,E_{opt})$, and $Y^{\prime}(E)<0$ for $E\in(E_{opt},1)$. ∎

######  Lemma 4.2.

The maximum yield (MY) for delayed impulsive harvesting model (1) with
$k\in{\mathbb{N}}$ is equal to the MY for non-delayed model (3) with the
optimal harvesting effort $E_{opt}=1-e^{-rT/2}$ and MY associated to (6).

###### Proof.

By Lemma 4.1, MY is attained at $E=E_{opt}$ and has the value of yield per
time

$MY_{delayed}=Y(E_{opt})=\frac{K_{c}}{T}\bigg{(}\frac{e^{rT/2}-1}{e^{rT/2}+1}\bigg{)}=MY_{non-
delayed}.$

In addition, the periodic solution when $E=E_{opt}$ becomes $\displaystyle
N^{*}(nT^{+})=\frac{K_{c}}{e^{rT/2}+1}$, as in (6). ∎

######  Theorem 4.3.

The solution corresponding to the MY of (1) satisfies (6) with $E=E_{opt}$. MY
is a MSY if and only if either $k=1$ or $k\geq 2$ and

$rT<-2\ln\bigg{(}1-2\cos\bigg{(}\frac{k\pi}{2k+1}\bigg{)}\bigg{)}\,.$ (27)

###### Proof.

By Lemma 3.1 the solutions to (1) satisfy (11). Linearizing the difference
equation $x_{n+1}=f(x_{n},x_{n-k})$ around the positive equilibrium $x^{*}$
given in (14), we get (24).

By Lemma 4.2, the yield is maximal whenever $E=E_{opt}=1-e^{-rT/2}$, and the
linearized difference equation becomes

$u_{n+1}=u_{n}-(1-e^{-rT/2})u_{n-k}$

with equilibrium (6). Then by Lemma 2.6, for $k=1$, as $1-e^{-rT/2}\in(0,1)$,
the solution is asymptotically stable for any $E\in(0,1)$. For $k>1$, the zero
solution of the linearized equation is locally asymptotically stable if and
only if

$1-e^{-rT/2}<2\cos\bigg{(}\frac{k\pi}{2k+1}\bigg{)}.$

The condition is equivalent to (27), leading to local asymptotic stability for
the positive equilibrium $x^{*}=K_{c}/(e^{rT/2}+1)$ of (11). Finally, by Lemma
3.2, a solution of (1) satisfying (6) is locally asymptotically stable, once
(27) is satisfied. By definition, a unique positive periodic solution
$N^{*}(t)$ with (6), for either $k=1$ or both $k\geq 2$ and $rT$ satisfying
(27), leads to MSY. ∎

Unlike the non-delay case, there is a possibility that the maximum yield is
not sustainable. To avoid extinction, the choice of harvesting efforts should
be among those leading to a sustainable yield. The set of such efforts is non-
empty, as the following statement guarantees.

######  Theorem 4.4.

Let $k\geq 2$ and

$E^{*}=\frac{2+e^{-rT}-\sqrt{e^{-rT}(e^{-rT}+8)}}{2}\,.$ (28)

Then $E^{*}<E_{opt}$, and for any $E\in(0,E^{*})$ the yield as given in (26)
is sustainable.

###### Proof.

Let us note that, first, $E^{*}$ defined in (28) is positive and, as
$4e^{-rT/2}>4e^{-rT}$, we have

$\sqrt{e^{-rT}(e^{-rT}+8)}>e^{-rT}+2e^{-rT/2}$

leading to $E^{*}<E_{opt}=1-e^{-rT/2}$.

For a fixed $E\in(0,1)$, we get a solution with $N^{*}(nT^{+})=x^{*}$,
corresponding to a yield as given in (26). As justified earlier, $Y(E)$ is an
increasing function of $E$ for $E\in(0,E_{opt})$.

By Theorem 3.10 and Lemma 3.2, the solution is locally asymptotically stable
for any $k$ if $\displaystyle e^{rT}>\frac{E+1}{(1-E)^{2}}$, which is
equivalent to $E^{2}-(2+e^{-rT})E+1-e^{-rT}>0$. The quadratic inequality is
also satisfied if $E\in(0,E^{*})$, meaning that for any $E\in(0,E^{*})$ the
yield is sustainable, which concludes the proof. ∎

######  Lemma 4.5.

If for some choice of $E_{s}\in(0,E_{opt}]$ the associated yield is
sustainable, then for any $E\in(0,E_{s}]$ the yield associated with $E$ is
also sustainable.

###### Proof.

If the yield associated with $E_{s}$ is sustainable, the associated solution
with $N^{*}(nT^{+})=x^{*}$ is locally asymptotically stable. By Lemma 3.8,
both inequalities in (22) must hold for $E_{s}$, where $\theta^{*}$ is a root
of (23).

Since for any $E\leq E_{opt}$,

$\frac{e^{-rT}}{(1-E)^{2}}\leq\frac{e^{-rT}}{(1-E_{opt})^{2}}=1<\frac{k+1}{k},$

it is clear that the first inequality in (22) is satisfied for any
$E\in(0,E_{opt}]$. Thus we can turn our attention to the second inequality,
denote its right-hand side for a fixed $rT$ as $h_{1}$,

$h_{1}(E):=\frac{(1+E)e^{-rT}}{2(1-E)}+\frac{(1-E)^{2}e^{rT}}{2}\,,$ (29)

and the left-hand side in (23) as

$h_{2}(\theta(E))=\frac{\sin(k\theta(E))}{\sin((k+1)\theta(E))},\quad\theta\in
I=(0,\pi/(k+1)).$ (30)

We have from (23),

$h_{2}^{\prime}(\theta(E))\frac{d\theta}{dE}=\frac{d\theta}{dE}\left[(1-E)^{2}e^{rT}\right]=-2(1-E)e^{rT}<0.$

Also,

$\displaystyle h_{2}^{\prime}(\theta)$
$\displaystyle=\frac{k\cos(k\theta)\sin((k+1)\theta)-(k+1)\cos((k+1)\theta)\sin(k\theta)}{\sin^{2}((k+1)\theta)}$
$\displaystyle=k\frac{\sin((k+1)\theta)\cos(k\theta)-\sin(k\theta)\cos((k+1)\theta)}{\sin^{2}((k+1)\theta)}-\frac{\cos((k+1)\theta)\sin(k\theta)}{\sin^{2}((k+1)\theta)}$
$\displaystyle=k\frac{\sin(\theta)}{\sin^{2}((k+1)\theta)}-\frac{\cos((k+1)\theta)\sin(k\theta)}{\sin^{2}((k+1)\theta)}$
$\displaystyle=\frac{\sin(\theta)}{\sin^{2}((k+1)\theta)}\bigg{(}k-\frac{\cos((k+1)\theta)\sin(k\theta)}{\sin(\theta)}\bigg{)}$

Now, $\sin(\theta)/\sin((k+1)\theta)^{2}>0$ $\forall\theta\in I$, and we will
show that $k\sin(\theta)-\cos((k+1)\theta)\sin(k\theta)>0$ leading to
$h_{2}^{\prime}(\theta)>0$. Since $\cos((k+1)\theta)\leq 1$,

$k\sin(\theta)-\cos((k+1)\theta)\sin(k\theta)>k\sin(\theta)-\sin(k\theta)=:H_{1}(\theta),$

where $H_{1}(0)=0$ and
$H_{1}^{\prime}(\theta)=k(\cos(\theta)-\cos(k\theta))>0$ since $\cos(\theta)$
is decreasing for all $\theta\in(0,\pi)$ and for $\theta\in I$, both $\theta$
and $k\theta$ are in $(0,\pi)$. Thus $H_{1}(\theta)>0$ for $\theta\in I$.

Since $h_{2}^{\prime}(\theta)>0$, the inequality $\displaystyle
h_{2}^{\prime}(\theta(E))\frac{d\theta}{dE}<0$ implies
$\displaystyle\frac{d\theta}{dE}<0$. Thus $\theta(E)$ decreases, and
$\cos(\theta(E))$ increases in $E$ as well. Further,

$h^{\prime}_{1}(E)=\frac{e^{-rT}}{(1-E)^{2}}-(1-E)e^{rT}<0$

for $E<1-e^{-2rT/3}$. By our assumption, $E\leq
E_{opt}=1-e^{-rT/2}<1-e^{-2rT/3}$, therefore $h^{\prime}_{1}(E)<0$, and
$h_{1}(E)$ decreases in E.

Since the yield associated with $E_{s}$ is sustainable, we have

$\cos(\theta(E_{s}))<h_{1}(E_{s}).$

Since $\cos(\theta(E))$ decreases, and $h_{1}(E)$ increases for decreasing
$E$, then for any $E\leq E_{s}$

$\cos(\theta(E))<h_{1}(E),$

and the second inequality in (22) is satisfied.

Since both inequalities in (22) are satisfied for $E\leq E_{s}$, the solution
associated with $E$ with $N^{*}(nT^{+})=x^{*}$ is locally asymptotically
stable, and the yield associated with $E$ is sustainable. ∎

######  Corollary 4.5.1.

Let $k=1$, then for any $E\in(0,E_{opt}]$ the yield is sustainable.

The statement follows immediately from Lemma 4.5 and Theorem 4.3, while Lemma
4.5 and Theorem 4.3 imply

######  Corollary 4.5.2.

Let $k\geq 2$, and
$rT<-2\ln\bigg{(}1-2\cos\bigg{(}\frac{k\pi}{2k+1}\bigg{)}\bigg{)}$. Then for
any $E\in(0,E_{opt}]$ the yield is sustainable.

Theorem 4.6 determines the maximum bound for a sustainable yield.

######  Theorem 4.6.

Let $k\geq 2$, and
$rT\geq-2\ln\bigg{(}1-2\cos\bigg{(}\frac{k\pi}{2k+1}\bigg{)}\bigg{)}$. Then
there exist $E^{**}\in(0,E_{opt}]$ and
$\theta^{*}\in\left(0,\frac{\pi}{k+1}\right)$ such that $(E^{**},\theta^{*})$
is a unique solution of

$\cos(\theta)=\frac{(1+E)e^{-rT}}{2(1-E)}+\frac{(1-E)^{2}e^{rT}}{2},\quad\frac{\sin(k\theta)}{\sin((k+1)\theta)}=(1-E)^{2}e^{rT},\quad\theta\in\left(0,\frac{\pi}{k+1}\right).$
(31)

For any $E\in(0,E^{**})$, the yield is sustainable, while for
$E\in[E^{**},E_{opt}]$ the yield is not sustainable.

###### Proof.

Let us prove that $E^{**}$ exists and is unique. By Theorem 4.3, when $k\geq
2$, $rT\geq-2\ln\bigg{(}1-2\cos\bigg{(}\frac{k\pi}{2k+1}\bigg{)}\bigg{)}$ and
$E=E_{opt}$, the associated solution is unstable and
$\cos(\theta(E_{opt}))\geq h_{1}(E_{opt})$. On the other hand, when $E\to
0^{+}$ we get in (22),
$\displaystyle\cos(\theta(0^{+}))=\frac{e^{-rT}}{2}+\frac{e^{rT}}{2}<\cosh(rT)=h_{1}(0^{+})$.

As $h_{2}(\theta(E))$ defined in (30) is continuous and monotone decreasing in
$E$, $E(\theta)$ is continuous and monotone for $\theta\in I=(0,\pi/(k+1))$,
and even on $[0,\pi/(k+1))$, if we define $h_{1}(0)=k/(k+1)$. The inverse
function is also monotone and continuous for $E\in[0,E_{opt}]$, as well as
$\cos(\theta(E))$. When $E\to 0^{+}$, $\cos(\theta(0^{+}))<h_{1}(0^{+})$ and
by assumption at the optimal harvesting effort $\cos(\theta(E_{opt}))\geq
h_{1}(E_{opt})$. By the continuity of $\cos(\theta(E))$ and the Intermediate
Value Theorem, there exists a solution $(\theta(E^{**}),E^{**})$ of (31) with
$E^{**}\in(0,E_{opt}]$. Moreover, by monotonicity this value is unique.

Since by Lemma 4.1 the yield $Y(E)$ is increasing for $E\in(0,E_{opt}]$, the
value of $Y(E^{**})$ is an upper bound, it does not satisfy (22) and thus is
not sustainable.

However for $E=(E^{**})^{-}=\lim_{E\to(E^{**})^{-}}E$, by the argument that as
$E$ decreases, $\cos(\theta(E))$ decreases and $h_{1}(E)$ increases,
$(E^{**})^{-}$ will satisfy the conditions in (22), and the associated yield
is sustainable. Then by Lemma 4.5, for any $E\in(0,E^{**})$ the associated
yield is sustainable. ∎

##  5 Numerical Simulations and Discussion

Let us illustrate sustainability of the optimal yield with simulations.

Figure 1: Solutions to (1) with $k=2$, $K_{c}=500$, and
$E=E_{opt}=1-e^{-rT/2}$. The figure on the left shows the optimal solution to
(1) (red), and the solution to (1) with $r=1$, $T=1$ and initial conditions
$N(-T)=120$, $N(0)=140$, $N(0^{+})=100$ (black). Since $0<rT<1.9248$ the
optimal solution is locally asymptotically stable. The figure on the right
shows the optimal solution to (1) (red), and the solution to (1) with $r=2.1$,
$T=1$ and initial conditions $N(-T)=120$, $N(0)=140$, $N(0^{+})=100$ (black).
Since $rT>1.9248$ the optimal solution is unstable.

In Fig. 1, stable and unstable solutions of (1) corresponding to the optimal
yield are compared. We assume that the delay is two impulsive periods
corresponding to $k=2$. To investigate stability, solutions were computed
until one of the following became true: the relative error of
$|N((n+1)T^{+})-N(nT^{+})|/N(nT^{+})$ was consistently less than $10^{-4}$,
$t=100T$, or the population went to extinction, i.e. $N(nT^{+})=0$. The
carrying capacity was chosen as $K_{c}=500$, and the harvesting effort was
chosen to be $E=E_{opt}=1-e^{-rT/2}$. Since $k=2>1$, Theorem 4.3 states that
the optimal solution is locally asymptotically stable if and only if
$0<rT<1.9248$. In both cases the optimal solution is plotted with a dashed red
line, while a solution with initial conditions which are different than the
optimal solution is shown via a black solid line.

In Fig. 1, left, $r=1$, $T=1$ and the initial conditions are $N(-T)=120$,
$N(0)=140$, $N(0^{+})=100$. Since $0<rT=1<1.9248$, by Theorem 4.3 the optimal
solution is locally asymptotically stable. Thus we observe a quick convergence
to the optimal solution, indicating that the population survives. In Fig. 1,
right, $r=2.1$, $T=1$ and the initial conditions are $N(-T)=120$, $N(0)=140$,
$N(0^{+})=100$. As $rT=2.1>1.9248$, by Theorem 4.3 the optimal solution is
unstable. Although our solution begins close to the optimal solution due to
the imposed initial conditions, we observe the population oscillating slightly
before going to extinction at $t=7$.

Figure 2: When $k=1$, $r=1.3747$, $T=1$, $K_{c}=307.1609$,
$E=1-e^{-(1.3747)(1)/2}=0.4971$ the positive equilibrium solution with
$N^{*}(nT^{+})\approx 102.7816$ is locally asymptotically stable. The dots in
this figure show whether or not the population survives or goes extinct with
the given initial conditions.

One thing that we wish to highlight is that local stability of the positive
periodic solution does not guarantee population survival for all possible
positive initial conditions, as Fig. 2 illustrates. In this figure, each dot
represents a solution to (1) ($k=1$) with a set of initial values $N(0)$ and
$N(0^{+})$. The parameters are $k=1$, $r=1.3747$, $T=1$, $K_{c}=307.1609$,
$E=E_{opt}=1-e^{-(1.3747)(1)/2}=0.4971$ and thus by Theorem 4.3 (see also
Theorem 3.7) the positive periodic solution with $N^{*}(nT^{+})\approx
102.7816$ is locally asymptotically stable. Fig. 2 tested the global stability
of the positive periodic solution by investigating whether or not different
combinations of initial conditions would result in survival of the population.
For each solution, both $N(0)$ and $N(0^{+})$ were chosen from a uniform
distribution of numbers in $[0,2K_{c}]=[0,614.3218]$. If after $10500$
iterations the solution $N(nT^{+})$ stayed positive (in fact, within
$[N^{*}(nT^{+})-10,N^{*}(nT^{+})+10]$), then we say that the population has
survived and the population is depicted by a green dot on Fig. 2. If at any
point within the $10500$ iterations the size of the population after
harvesting was less than $10^{-3}$ then the population was said to have gone
to extinction, and the population was depicted by a black dot.

Since the optimal positive periodic solution with $N^{*}(nT^{+})\approx
102.7816$ is locally asymptotically stable, it is natural that we observe the
population surviving for wide ranges of initial values. However, if $N(0)$ is
much larger than $N(0^{+})$ or if $N(0)$ is just generally very large, the
population does not survive and goes to extinction. This highlights the lack
of global attractivity of the optimal solution which is locally asymptotically
stable for all $rT>0$.

We proved that asymptotic stability of the difference equation which
corresponds to the optimal solution of impulsive model (1) is $k$-dependent.
In fact, we can show that as $k\to\infty$, the range of possible $rT$ values
that allow a locally asymptotically stable equilibrium shrinks.

Figure 3:  This function (5) is monotonically decreasing since
$f^{\prime}(k)<0$ for $k>1$, and $\displaystyle\lim_{k\to\infty}f(k)=0.$

Fig. 3 illustrates that as the delay $k$ increases, the range of possible $rT$
values which allow for optimal sustainable harvesting becomes smaller. In the
figure, the upper bound on sustainable $rT$ values

$f(k)=-2\ln\bigg{(}1-2\cos\bigg{(}\frac{\pi k}{2k+1}\bigg{)}\bigg{)},$

is plotted and it is easy to see that as $k$ increases, $f(k)\to 0$. This
parallels the theory of continuous delayed harvesting where a growing delay
exhibits a destabilizing effect on the model.

The results of the paper can be summarized as follows:

  1. 1.

With delayed impulsive harvesting, we cannot guarantee positivity of a
solution with non-negative and non-trivial initial conditions, moreover,
extinction in finite time is possible.

  2. 2.

The delay does not influence the maximum yield but can significantly influence
its sustainability.

  3. 3.

With sharp local stability conditions, the optimal solution associated with
the maximum yield is locally asymptotically stable in both non-delay [30] and
one-step-delay cases. For longer delays, there are bounds on $rT$ to attain
MSY: the longer the delay is, the more frequent impulses should be.

For the choice of optimal harvesting, the increasing value of the intrinsic
growth rate has negative influence on sustainability (see Figs. 1 and 3). This
is well aligned with the enrichment paradox when increased productivity can
lead to higher extinction risks.

Let us discuss time delays in the impulsive system. The real delay in
information is $kT$ corresponding to $k$ cycles with $T$-time duration of each
cycle. Assuming $r=1$ and applying the bound for $rT=T$ in (27), we evaluate
the delay $kT$. We get using L’Hôspital’s rule,

$\lim_{k\to+\infty}\bigg{[}-2k\ln\bigg{(}1-2\cos\bigg{(}\frac{\pi
k}{2k+1}\bigg{)}\bigg{)}\bigg{]}=\pi.$

Since $kf(k)$, with $f$ plotted in Fig. 3, is monotone decreasing with the
lower bound of $\pi$, the delay $kT$ not exceeding $\pi$ (generally, $\pi/r$)
will lead to MSY for any $k$. Note that $rT<\pi$, where $T$ is the delay, is a
sharp stability constant in the autonomous continuous delay case.

This research extends the classical continuous harvesting models of Clark [1]
as well as impulsive harvesting models with no delay [2, 30] which were
generalized to impulsive harvesting with delay, outlining new instability
effects. The models closest to those considered in the present paper were
explored in [15, 17, 30] and to some extent [2]. In particular, optimal
harvesting was discussed in [2, 30], with the same MSY being derived as in the
present paper, however in the absence of the delay in the impulsive part, the
MSY was unconditional. Delayed impulsive harvesting for the logistic equation
and asymptotic stability of its positive periodic solutions was explored in
[15, 17]. Compared to the previous work, the main points of similarity and
difference can be outlined as follows.

  1. 1.

We use the method of reducing an impulsive model to a difference equation,
which was extensively applied in [2, 30] and also used in [15]. However, the
paper [15] avoids high order difference equations by instead reducing an
equation with delayed impulses to a non-delay system. The main technique
applied to study impulsive models is Lyapunov-Razumikhin, leading to various
global stability results. But the system considered in this paper is non-
smooth (truncated) in the sense that the solution is assumed to be identically
zero, once $N(t)=0$ for some $t$, and we are not aware of this method applied
to such models.

  2. 2.

In both [15, 17], the impulse delay is within one-cycle period i.e. less than
or equal to $T^{-}$. If the equation is non-delayed, and we use as a reference
value the same $T$-period, this value at $t=(jT)^{-}$ can be expressed as
$\beta(r,\tau)x((jT)^{-})$, where $\beta$ is explicitly computed and dependent
on the delay $\tau$, the intrinsic growth rate $r$ and, generally, on the
growth law (logistic here). Then the situation is to some extent reduced to
non-delayed impulses when local and global stability coincide. If the delay is
equal to the impulse period $T=T^{+}$, collapse after a harvesting event is
possible when the positive periodic solution is locally stable, see Fig. 2. To
the best of our knowledge, the present paper is the first to report extinction
in finite time and discrepancy between local and global stability for
impulsively harvested equations. For applied models, this could possibly
explain population collapses. On the other hand, [15] explored subtler
properties, such as cycles and bifurcations.

Let us dwell on the extension of this research, its modifications and
generalizations. It would be advantageous to get sufficient conditions on the
initial values guaranteeing existence of a positive solution and to describe
global asymptotic stability for such initial conditions. Note that this
problem cannot be solved with modifying the per capita growth rate in the
equation. Even for the Gompertz growth model which is sustainable for any
level of harvesting, delayed impulsive harvesting can cause immediate
extinction, once the stock decays quickly and is significantly overestimated
during the harvesting event.

The approach of [2] where a certain deduction is incorporated in each
harvesting event, not contributing to the yield, can also be combined with
delayed impulsive harvesting. It is expected that the maximum yield in this
case will not change when impulses are delayed, but investigating its local or
global attractivity is still an open question. The current research can be
placed in the context of $T$-periodic or almost periodic $T(s)$.

Another minor extension will be including delays in the continuous part, for
example, switching from the logistic to the Hutchinson equation

$\frac{dN}{dt}=rN(t)\bigg{(}1-\frac{N(t-jT)}{K_{c}}\bigg{)}$

for some $j\in{\mathbb{N}}$, with the same impulsive conditions. Here again,
the maximum yield will not change, but an additional delay is expected also to
contribute to instability of the model.

Including a delay within impulsive harvesting allowed us to observe quite
realistic effects, such as population collapse even when the model parameters
predicted local stability of the unique positive periodic solution.

Certainly incorporating impulsive delayed harvesting in non-autonomous or
structured population models will lead to richer dynamics, in particular,
considering

  * 1.

non-autonomous models with variable parameters, in particular, considering
periodic and almost periodic solutions (similarly to [33]), stability,
bifurcations, as well as optimal harvesting policies and their sustainability;

  * 2.

predator-prey systems where either one of the species or both are subject to
harvesting, structured and Lotka-Volterra models which - without delay in
impulsive harvesting or an impulsive system with continuous harvesting - are
well-studied areas, see, for example, [9, 34, 35, 36, 37, 38, 39, 40, 41, 42];

  * 3.

fractional derivative models and spatially distributed populations, as well as
systems within fractional settings to describe spatial interactions.

## Acknowledgment

The authors are grateful to the anonymous reviewer whose thoughtful comments
significantly contributed to the quality of presentation. The authors also
acknowledge the support of NSERC, the grant RGPIN-2020-03934.

## References

  * [1] Clark CW. Mathematical Bioeconomics. 2nd ed. Hoboken New Jersey: Wiley-Interscience; 1990.
  * [2] Braverman E, Mamdani R. Continuous versus pulse harvesting for population models in constant and variable environment. J Math Biol 2008;57:413-34.
https://doi.org/10.1007/s00285-008-0169-z

  * [3] Xiao Y, Cheng D, Qin H. Optimal impulsive control in periodic ecosystems. Systems Control Lett 2006;55(7):558-565.
https://doi.org/10.1016/j.sysconle.2005.12.003

  * [4] Wang A, Xiao Y, Smith? R. Using non-smooth models to determine thresholds for microbial pest management. J Math Biol 2019;78(5):1389-1424.
https://doi.org/10.1007/s00285-018-1313-z

  * [5] Gao S, Chen L, Teng Z. Impulsive vaccination of an SEIRS Model with time delay and varying total population size. Bull Math Biol 2007;69:731-45.
https://doi.org/10.1007/s11538-006-9149-x

  * [6] Córdova-Lepe F, Del Valle R, Robledo G. A pulse fishery model with closures as function of the catch: Conditions for sustainability. Math Biosci 2012;239(1):169-77.
https://doi.org/10.1016/j.mbs.2012.04.004

  * [7] Bainov D, Simeonov P. Impulsive Differential Equations: Periodic Solutions and Applications. Longman Scientific; 1993.
  * [8] Smith H. An Introduction to Delay Differential Equations with Applications to the Life Sciences. New York: Springer; 2011.
  * [9] Kar TK. Selective harvesting in a prey-predator fishery with time delay. Math Comput Modelling 2003;38(3-4):449-58.
https://doi.org/10.1016/S0895-7177(03)90099-9

  * [10] Xu C, Zhang W, Li P. Periodic oscillating dynamics for a delayed Nicholson-type model with harvesting terms. Math Probl Eng Art 2021:9734342.
https://doi.org/10.1155/2021/9734342

  * [11] Akça, H, Berezansky L, Braverman E. On linear integro-differential equations with integral impulsive conditions. Z Anal Anwendungen 1996;15(3):709-27.
  * [12] Agarwal R, Hristova S, O’Regan D. Non-instantaneous Impulses in Differential Equations. Springer Cham; 2017.
  * [13] Liu X, Zhang K. Impulsive Systems on Hybrid Time Domains. Springer Cham; 2019.
  * [14] Church, KEM, Duchesne GW. Rigorous continuation of periodic solutions for impulsive delay differential equations. Appl Math Comput 2022;415:126733.
https://doi.org/10.1016/j.amc.2021.126733

  * [15] Church KEM, Lui X. Bifurcation analysis and application for impulsive systems with delayed impulses. Internat J Bifur Chaos Appl Sci Engrg 2017;27(12):1750186.
https://doi.org/10.1142/S0218127417501863

  * [16] Church KEM, Lui X. Computation of centre manifolds and some codimension-one bifurcations for impulsive delay differential equations. J Differential Equations 2019;267(6):3852-921.
https://doi.org/10.1016/j.jde.2019.04.022

  * [17] Pei Y, Chen L, Li C, Wang C. Impulsive selective harvesting in a logistic fishery model with time delay. J Biol Syst 2006;14(1):91-9.
https://doi.org/10.1142/S0218339006001738

  * [18] Cortés J-C, Delgadillo-Alemán SE, Kú-Carrillo RA, Villanueva R-J. Probabilistic analysis of a class of impulsive linear random differential equations via density functions. Appl Math Lett 2021;121:107519.
https://doi.org/10.1016/j.aml.2021.107519

  * [19] Wang Z, Liu M. Optimal impulsive harvesting strategy of a stochastic Gompertz model in periodic environments. Appl Math Lett 2022;125:107733.
https://doi.org/10.1016/j.aml.2021.107733

  * [20] Yang B, Cai Y, Wang K, Wang W. Optimal harvesting policy of logistic population model in a randomly fluctuating environment. Phys A 2019;526:120817.
  * [21] Cruz-Rivera E, Ramírez CH, Vasilieva O. Catch-to-stock dependence: the case of small pelagic fishery with bounded harvesting effort. Nat Resour Model 2019;32(1):e12193.
https://doi.org/10.1111/nrm.12193

  * [22] Brauer F, Castillo-Chavez C. Mathematical Models in Population Biology and Epidemiology. New York: Springer-Verlag; 2001.
  * [23] Grey S, Lenhart S, Hilker FM, Franco D. Optimal control of harvest timing in discrete population models. Nat Resour Model 2021;34(3):e12321.
https://doi.org/10.1111/nrm.12321

  * [24] Dawed MY, Kebedow KG. Coexistence and harvesting optimal policy in three species food chain model with general Holling type functional response. Nat Resour Model 2021;34(3):e12316.
https://doi.org/10.1111/nrm.12316

  * [25] Demir M, Lenhart S. Optimal sustainable fishery management of the Black Sea anchovy with food chain modeling framework. Nat Resour Model 2020;33(2):e12253.
https://doi.org/10.1111/nrm.12253

  * [26] Pei Y, Chen M, Liang X, Li C. Model-based on fishery management systems with selective harvest policies. Math Comput Simulation 2019;156:377-95.
https://doi.org/10.1016/j.matcom.2018.08.009

  * [27] Upadhyay RK, Tiwari SK. Ecological chaos and the choice of optimal harvesting policy. J Math Anal Appl 2017;448(2):1533-59.
https://doi.org/10.1016/j.jmaa.2016.11.054

  * [28] Jana D, Agrawal R, Upadhyay RK, Samanta GP. Ecological dynamics of age selective harvesting of fish population: Maximum sustainable yield and its control strategy. Chaos Solitons Fractals 2016;93:111-22.
https://doi.org/10.1016/j.chaos.2016.09.021

  * [29] Skonhoft A, Gong P. Maximum sustainable yield harvesting in an age-structured fishery population model. Nat Resour Model 2016;29(4):610-32.
https://doi.org/10.1111/nrm.12118

  * [30] Zhang X, Shuai Z, Wang K. Optimal impulsive harvesting policy for single population. Nonlinear Anal Real World Appl 2003;4(4):639-651.
https://doi.org/10.1016/S1468-1218(02)00084-6

  * [31] Elaydi S. An Introduction to Difference Equations. New York: Springer; 2005.
  * [32] Kulenovic MRS, Ladas G. Dynamics of Second Order Rational Difference Equations with Open Problems and Conjectures. Chapman & Hall/CRC; 2002.
  * [33] Yao Z. Existence and exponential stability of the unique positive almost periodic solution for impulsive Nicholson’s blowflies model with linear harvesting term. Appl Math Model 2015;39(23-24):7124-33.
https://doi.org/10.1016/j.apm.2015.03.002

  * [34] Jiao J, Liu Z, Li L, Nie X. Threshold dynamics of a stage-structured single population model with non-transient and transient impulsive effects. Appl Math Lett 2019;97:88-92.
https://doi.org/10.1016/j.aml.2019.05.024

  * [35] Liu J, Hu J, Yuen P. Extinction and permanence of the predator-prey system with general functional response and impulsive control. Appl Math Model 2020;88:55-67.
https://doi.org/10.1016/j.apm.2020.06.033

  * [36] Meng X, Zhang L. Evolutionary dynamics in a Lotka-Volterra competition model with impulsive periodic disturbance. Math Methods Appl Sci 2016;39(2):177-88.
https://doi.org/10.1002/mma.3467

  * [37] Quan Q, Tang W, Jiao J, Wang Y. Dynamics of a new stage-structured population model with transient and nontransient impulsive effects in a polluted environment. Adv Difference Equ 2021:518.
https://doi.org/10.1186/s13662-021-03667-4

  * [38] Sharma A, Gupta B, Dhar J, Srivastava SK, Sharma P. Stability analysis and optimal impulsive harvesting for a delayed stage-structured self dependent two compartment commercial fishery model. Int J Dyn Control 2022;10(4):1119-29.
https://doi.org/10.1007/s40435-021-00866-5

  * [39] Xiao Q, Dai B, Xu B, Bao L. Homoclinic bifurcation for a general state-dependent Kolmogorov type predator-prey model with harvesting. Nonlinear Anal Real World Appl 2015;26:263-73.
https://doi.org/10.1016/j.nonrwa.2015.05.012

  * [40] Wang H, Ou C, Dai B. The existence and stability of order-1 periodic solutions for an impulsive Kolmogorov predator-prey model with non-selective harvesting. J Appl Anal Comput 2021;11(3):1348-70.
  * [41] Wang L, Zhang H, Liu S. On the existence of almost periodic solutions of impulsive non-autonomous Lotka-Volterra predator-prey system with harvesting terms. AIMS Math 2022;7(1):925-38.
https://doi.org/10.3934/math.2022055

  * [42] Zhang M, Zhao Y, Chen L, Li Z. State feedback impulsive modeling and dynamic analysis of ecological balance in aquaculture water with nutritional utilization rate. Appl Math Comput 2020;373:125007.
https://doi.org/10.1016/j.amc.2019.125007