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.

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.

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}|

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.

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.

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.

Refer to caption — 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$ .

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.

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.

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.

The delay does not influence the maximum yield but can significantly influence its sustainability.
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.

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.

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