# Harvesting of a stochastic population under a mixed regular-singular control
formulation

Ky Q. Tran Department of Applied Mathematics and Statistics, State University
of New York - Korea Campus, Yeonsu-Gu, Incheon, Korea 21985,
<EMAIL_ADDRESS>The research of this author was supported by the
National Research Foundation of Korea grant funded by the Korea Government
(MIST) NRF-2021R1F1A1062361. Le T. N. Bich Department of Mathematics, Hue
University of Education, Hue University, Hue city, Vietnam,
<EMAIL_ADDRESS> George Yin Department of Mathematics, University of
Connecticut, Storrs, CT 06269, USA<EMAIL_ADDRESS>The research of this
author was supported in part by the Air Force Office of Scientific Research.

###### Abstract

This work focuses on optimal harvesting-renewing for a stochastic population.
A mixed regular-singular control formulation with a state constraint and
regime-switching is introduced. The decision makers either harvest or renew
with finite or infinite harvesting/renewing rates. The payoff functions depend
on the harvesting/renewing rates. Several properties of the value functions
are established. The limiting value function as the white noise intensity
approaches infinity is identified. The Markov chain approximation method is
used to find numerical approximation of the value function and optimal
strategies.

Key words. Harvesting problem, controlled diffusion, singular control, state
constraint, Markovian switching

##  1 Introduction

This work focuses on harvesting and renewing strategies for stochastic
ecosystems that are represented by controlled stochastic differential
equations. Mathematically, the problem we consider belongs to a class of
singular stochastic control problems, namely, harvesting-type problems. Such
problems have been studied extensively in various settings for different
domain of applications; see [1, 2, 3, 9, 10, 11, 12, 13, 19, 24]. To mention
just a few of the recent developments, we refer to [25, 26] for single species
and interacting population systems with regime switching. The paper [6]
focuses on sustainable harvesting policies under long-run average criteria,
which are further studied in [18]; see also [21] for a related work on a
predator-prey system and [15] for an ergodic two-sided singular control
formulation. In [7, 8], the authors study ecosystems that allow for the
modelling of both renewing and harvesting. Optimal exploitation problems of
renewable natural resources, which are harvesting-type problems, are studied
in [14, 27]. Related works on a general one-dimensional diffusion that is
reflected at zero can be found in [4] and references therein. Intensive
treatment of stochastic population systems and hybrid diffusions can be found
in [20, 28], while applications in various areas of singular stochastic
control can be found in [5, 23] and many references therein.

In this work, we propose a generalized harvesting model for a stochastic
population. The controller can perform either a regular control or a singular
control to harvest and renew the species. Moreover, we also consider the
control objective associated with renewing and harvesting. It should be noted
that in optimal harvesting formulations in [1, 2, 25, 7, 8], the control cost
is simply combined with the price functions. As a result, it is frequently
seen that when the manager decides to harvest (resp., renew), she should do
that with the maximum possible harvesting rates (resp. renewing rates); see
[1, 7, 8]. In [14, 27], interesting phenomena appear when the control
objective is taken into account. For instance, the maximal possible rates are
no longer optimal for harvesting and renewing in certain cases. With the
generalized formulation proposed in this work, the distinctions are even more
pronounced. In addition, a new twist, namely, a state constraint is considered
in this work. In particular, the time horizon of the control problem is
$[0,\tau]$ where $\tau$ is the first time the population process is below a
predetermined level. Another important issue of interest is the impact of a
white noise with large intensity. For Kolmogorov-type ecosystems, it is known
that very large white noise make the species extinct, and have a major impact
on harvesting actions; see [1, 20, 26]. A novelty of the paper is the
identification of the limiting value function as the white noise intensity
approaches infinity for a general stochastic population.

In contrast to the existing results, our new contributions in this paper are
as follows. (i) We formulate a harvesting problem with renewing and the
consideration of control objective, a state constraint, and regime-switching.
Both bounded and unbounded harvesting-renewing rates are allowed. (ii) We
establish the finiteness and the continuity of the value function. We show
that in common cases, it is optimal to keep the population size in a compact
set. (iii) We study the impact of a white noise with large intensity on
harvesting. (iv) Based on the Markov chain approximation method, we construct
a controlled Markov chain to approximate the given controlled population. It
enables us to approximate the value function and near-optimal strategies.

The rest of our work is organized as follows. Section 2 begins with the
problem formulation. Section 3 focuses on properties of the value function and
the impact of a white noise with large intensity. In Section 4, we construct a
controlled Markov chain to approximate the given controlled population system.
Finally, the paper is concluded with several numerical examples for
illustration together with additional remarks in the last section.

##  2 Formulation

We work with a complete filtered probability space
$(\Omega,\mathcal{F},{\mathbb{P}},\\{\mathcal{F}_{t}\\})$ with the filtration
$\\{\mathcal{F}_{t}\\}$ satisfying the usual condition (i.e., it is right-
continuous and ${\mathcal{F}}_{0}$ contains all the null sets). Suppose that
the population size $\xi(t)$ of a species at time $t$ is given by

$d\xi(t)=b\big{(}\xi(t),\Lambda(t)\big{)}dt+\sigma\big{(}\xi(t),\Lambda(t)\big{)}dw(t).$
(2.1)

In the model, $w(\cdot)$ is a one-dimensional standard Brownian motion and
$\Lambda(\cdot)$ is a continuous time Markov chain taking values in a finite
set $\mathcal{M}=\\{1,2,\dots,m_{0}\\}$ generated by
$\Gamma=(\Gamma_{\alpha\ell})_{m_{0}\times m_{0}}$. Moreover, $\Lambda(\cdot)$
and $w(\cdot)$ are independent and $\\{\mathcal{F}_{t}\\}$-adapted. Let
$\mathbb{R}_{+}=[0,\infty)$ and $\mathbb{Z}_{+}=\\{0,1,2,\dots\\}$. For a real
number $x$, we denote $x^{+}=\max\\{x,0\\}$ and $x^{-}=\max\\{-x,0\\}$. Thus,
$x=x^{+}-x^{-}$ and $|x|=x^{+}+x^{-}$.

To proceed, we introduce the generator of the process
$\big{(}\xi(t),\Lambda(t)\big{)}$. For a function
$\Phi(\cdot,\cdot):\mathbb{R}\times{\cal{M}}\mapsto\mathbb{R}$ satisfying
$\Phi(\cdot,\alpha)\in C^{2}(\mathbb{R})$ for each $\alpha\in{\cal{M}}$, we
define

$\begin{array}[]{ll}{\cal
L}\Phi(x,\alpha)&\\!\\!\\!\displaystyle=b(x,\alpha)\Phi^{\prime}(x,\alpha)+\dfrac{1}{2}\sigma^{2}(x,\alpha)\Phi^{\prime\prime}(x,\alpha),\end{array}$

where $\Phi^{\prime}(\cdot,\alpha)$ and $\Phi^{\prime\prime}(\cdot,\alpha)$
denote the first and second derivatives of $\Phi(\cdot,\alpha)$, respectively.

Next, we suppose that the controller can renew or harvesting by exercising an
impulsive control or a regular control. To take into account both cases, we
assume the dynamics of the species is given by

$X(t)=x+\int\limits_{0}^{t}b\big{(}X(s),\Lambda(s)\big{)}ds+\int\limits_{0}^{t}\sigma\big{(}X(s),\Lambda(s)\big{)}dw(s)-\int_{0}^{t}f\big{(}X(s),C(s)\big{)}ds-Y(t)+Z(t),$
(2.2)

where $x\in\mathbb{R}_{+}$, $X(t)$ is the population size at time $t\geq 0$,
$f:\mathbb{R}_{+}\times\mathcal{U}\mapsto\mathbb{R}$ is the harvesting-
renewing rate corresponding to the control $C(\cdot)$ taking values in a
nonempty compact set $\mathcal{U}$ in $\mathbb{R}$, while $Y(\cdot)$ and
$Z(\cdot)$ are impulsive controls.

Notation. For each time $t$, $X(t-)$ represents the state before harvesting
starts at time $t$, while $X(t)$ is the state immediately after. We assume the
initial population size to be $X(0-)=x$ and initial regime to be
$\Lambda(0)=\alpha$, respectively. Hence $X(0)$ may not be equal to $X(0-)$
due to an impulsive harvesting $Y(0)$ or an impulsive renewing $Z(0)$.
Throughout the paper we use the convention that $Y(0-)=Z(0-)=0$. The jump
sizes of $Y(\cdot)$ and $Z(\cdot)$ at time $t$ are denoted by $\Delta
Y(t)=Y(t)-Y(t-)$ and $\Delta Z(t)=Z(t)-Z(t-)$, respectively. Thus,

$Y(t)=\sum\limits_{0\leq s\leq t}\Delta Y(s),\quad Z(t)=\sum\limits_{0\leq
s\leq t}\Delta Z(s).$

Suppose $\lambda\in\mathbb{R}_{+}$ and denote

$\SS=\\{x\in\mathbb{R}_{+}:x\geq\lambda\\}.$

In this work, we consider the harvesting problem on the time horizon
$[0,\tau]$, where

$\tau=\inf\\{t\geq 0:X(t)\notin\SS\\}.$

The price per unit of the species is a positive constant $a_{1}$. The
harvesting and renewing control is costly. Consider a function
$g:\mathbb{R}_{+}\times{\cal{M}}\times\mathcal{U}\mapsto\mathbb{R}_{+}$ and
two positive constants $a_{2}$ and $a_{3}$. The accumulated cost of the
regular control $C(\cdot)$ is

$\int_{0}^{\tau}e^{-\delta s}g\big{(}X(s),\Lambda(s),C(s)\big{)}ds,$

and the accumulated cost of the impulsive controls $Y(\cdot)$ and $Z(\cdot)$
are

$\int\limits_{0}^{{\tau}}e^{-\delta s}a_{2}dY(s)\ \hbox{ and }\
\int\limits_{0}^{{\tau}}e^{-\delta s}a_{3}dZ(s),$

respectively, where $\delta>0$ is the discount factor. Note that we separate
cost and income here. Later we will take into account of both and set up the
payoff. Regarding the regular control $C(\cdot)$, the accumulated income of
selling the harvested amount is $\int\limits_{0}^{{\tau}}e^{-\delta
s}a_{1}f^{+}\big{(}X(s),C(s)\big{)}ds$ and the accumulated expense of the
renewed amount is $\int\limits_{0}^{{\tau}}e^{-\delta
s}a_{1}f^{-}\big{(}X(s),C(s)\big{)}ds$. For notational simplicity, we define
the price-cost function
$p:\mathbb{R}_{+}\times{\cal{M}}\times\mathcal{U}\to\mathbb{R}$ given by

$p(x,\alpha,c)=a_{1}f(x,c)-g(x,\alpha,c)\quad\text{for}\quad(x,\alpha,c)\in\mathbb{R}_{+}\times{\cal{M}}\times\mathcal{U}.$

We also define

$q=a_{1}-a_{2},\quad r=a_{1}+a_{3}.$

Thus, for a harvesting-renewing strategy $\Psi\equiv(C,Y,Z)$, we define the
payoff functional as

$J(x,\alpha,\Psi)={\mathbb{E}}_{x,\alpha}\Big{[}\int_{0}^{{\tau}}e^{-\delta
s}p\big{(}X(s),\Lambda(s),C(s)\big{)}ds+\int_{0}^{\tau}e^{-\delta
s}qdY(s)-\int_{0}^{\tau}e^{-\delta s}rdZ(s)\Big{]},$ (2.3)

where ${\mathbb{E}}_{x,\alpha}$ denotes the expectation with respect to the
probability law when $X(0-)=x$ and $\Lambda(0)=\alpha$.

Control strategy. Let $\mathcal{A}_{x,\alpha}$ denote the collection of all
admissible controls with initial condition $X(0-)=x$, $\Lambda(0)=\alpha$. A
harvesting-renewing strategy $\Psi=(C,Y,Z)$ will be in
$\mathcal{A}_{x,\alpha}$ if it satisfies the following conditions.

  * (a)

The processes $C(\cdot)$, $Y(\cdot)$, and $Z(\cdot)$ are adapted to
$\sigma\\{w(s),\Lambda(s):0\leq s\leq t\\}$; $C(\cdot)$ takes values in
$\mathcal{U}$, $Y(\cdot)$ and $Z(\cdot)$ are impulsive controls with non-
decreasing, nonnegative, and right continuous sample paths; $\Delta Y(s)\Delta
Z(s)=0$ for any $s\in\mathbb{R}_{+}$; $\Delta Y(s)=\Delta Z(s)=0$ for any
$s\geq\tau$;

  * (b)

The system (2.2) has a unique solution $X(\cdot)$ with $X(t)\in\SS$ for any
$t\in[0,\tau]$;

  * (c)

$0\leq J(x,\alpha,\Psi)<\infty$, where $J(\cdot)$ is the functional defined in
(2.3).

The problem we are interested in is to maximize the payoff functional and find
an optimal strategy $\Psi^{*}=(C^{*},Y^{*},Z^{*})\in\mathcal{A}_{x,\alpha}$
such that

$J(x,\alpha,\Psi^{*})=V(x,\alpha):=\sup\limits_{\Psi\in\mathcal{A}_{x,\alpha}}J(x,\alpha,\Psi).$
(2.4)

The function $V(\cdot)$ is called the value function.

The standing assumptions are given below.

  * (A)

    * (a)

For any $n\in\mathbb{Z}_{+}$, there exists a positive constant $K_{n}$ such
that for any $x,y\in\mathbb{R}_{+}$ with $|x|\leq n,|y|\leq n$ and any
$\alpha\in{\cal{M}}$,

$|b(x,\alpha)-b(y,\alpha)|+|\sigma(x,\alpha)-\sigma(y,\alpha)|\leq
K_{n}|x-y|.$

Moreover, for each initial condition
$(x,\alpha)\in\mathbb{R}_{+}\times{\cal{M}}$, the population system (2.1) has
a unique global solution.

    * (b)

The control set $\mathcal{U}$ is a nonempty compact set of real numbers,
$0\in\mathcal{U}$, and $a_{1},a_{2},a_{3}$ are positive constants. The
function $g(\cdot,\cdot,\cdot)$ is continuous and bounded on
$\mathbb{R}_{+}\times{\cal{M}}\times\mathcal{U}$. The function
$f(\cdot,\cdot)$ is continuous and bounded on
$\mathbb{R}_{+}\times\mathcal{U}$ and $f(\cdot,0)=0$.

######  Remark 2.1.

Here for simplicity, we require the systems of equations having a global
solution. In fact, the existence and uniqueness of global solutions for a
large class of Kolmogorov systems are guaranteed under suitable conditions;
see the recent work [22] and the references therein, and also [20] for various
sufficient conditions so that the population system (2.1) has a unique global
solution. Of course, here we are dealing with a controlled systems so some
modifications are needed. Nevertheless, in order to concentrate on our main
task, we chose to simply assume this condition.

The mixed regular-singular control formulation together with the consideration
of a state constraint, cost functions allow us to take into account various
aspects of harvesting-type problems that have not been considered to date. For
instance, depending on the costs, one can choose to apply either an impulsive
harvesting/renewing or harvest and renew through the regular control.

######  Remark 2.2.

The consideration of $\lambda$ is motivated by sustainability. To be best of
our knowledge, the available literature focuses on the case $\lambda=0$ in
which the optimal or near-optimal harvesting strategies might drive the
population process to a very low level or extinction; see [1, 7, 25]. Because
of the state constraint, one needs to treat carefully control actions when the
population size is close to $\lambda$.

The function $f:\mathbb{R}_{+}\times\mathcal{U}\mapsto\mathbb{R}$ is the
harvesting-renewing rate corresponding to the control $C(\cdot)$. In [8], the
authors consider the case $f(x,c)=c$ for any $c\in\mathcal{U}$. We can also
consider $f(x,c)=cx$ for $(x,c)\in[\lambda,\kappa]\times\mathcal{U}$, where
$\kappa$ is a positive constant. In other words, the bounded rate at which the
species can be harvested or renewed is proportional and positively related to
the current size. Such a consideration is motivated by the observations in [1,
Section 3] and [14].

##  3 Properties of the value function

This section is devoted to several properties of the value function. We begins
with a lemma that allows us to establish the finiteness and construct upper
bounds of the value function.

######  Lemma 3.1.

Let (A) be satisfied. Suppose that there exists a function
$\Phi:\SS\times{\cal{M}}\mapsto\mathbb{R}_{+}$ such that
$\Phi(\cdot,\alpha)\in C^{2}({\SS})$ for each $\alpha\in{\cal{M}}$ and
$\Phi(\cdot)$ solves the following coupled system of quasi-variational
inequalities

$\max\limits_{i}\bigg{\\{}\mathcal{G}\Phi(x,\alpha),q-\Phi^{\prime}(x,\alpha),\Phi^{\prime}(x,\alpha)-r\bigg{\\}}\leq
0\quad\text{for}\quad(x,\alpha)\in\SS\times{\cal{M}},$ (3.1)

where

$\mathcal{G}\Phi(x,\alpha)=(\mathcal{L}-\delta)\Phi(x,\alpha)+\max\limits_{c\in\mathcal{U}}\Big{[}p(x,\alpha,c)-\Phi^{\prime}(x,\alpha)f(x,c)\Big{]}.$

Recall that $\delta$ is the discount factor given in the payoff functional
(2.3). Then we have

$V(x,\alpha)\leq\Phi(x,\alpha)\quad\text{for}\quad(x,\alpha)\in\SS\times{\cal{M}}.$

###### Proof.

For a fixed $(x,\alpha)\in\SS\times{\cal{M}}$ and
$\Psi=(C,Y,Z)\in\mathcal{A}_{x,\alpha}$, let $X$ denote the corresponding
harvested process. Choose $N$ sufficiently large so that $|x|<N$. Define

$\tau_{N}=\inf\\{t\geq 0:X(t)\geq N\\},\quad T_{N}=N\wedge\tau_{N}\wedge\tau.$

Then

$\tau_{N}\to\infty\quad\text{and}\quad T_{N}\to\tau\quad\text{almost surely as
}\quad N\to\infty,$ (3.2)

where $\tau=\inf\\{t\geq 0:X(t)\notin\SS\\}.$ Then Dynkin’s formula leads to

$\begin{array}[]{ll}{\mathbb{E}}&\\!\\!\\!\displaystyle\big{[}e^{-\delta
T_{N}}\Phi\big{(}X(T_{N}),\Lambda(T_{N})\big{)}\big{]}=\Phi(x,\alpha)+{\mathbb{E}}\int_{0}^{T_{N}}e^{-\delta
s}({\cal L}-\delta)\Phi\left(X(s),\Lambda(s)\right)ds\\\
&\\!\\!\\!\displaystyle\qquad\qquad-{\mathbb{E}}\int_{0}^{T_{N}}e^{-\delta
s}\Phi^{\prime}\left(X(s),\Lambda(s)\right)f\big{(}X(s),C(s)\big{)}ds\\\
&\\!\\!\\!\displaystyle\qquad\qquad+{\mathbb{E}}\sum\limits_{0\leq s\leq
T_{N}}e^{-\delta
s}\Big{[}\Phi\left(X(s),\Lambda(s-)\right)-\Phi\left(X(s-),\Lambda(s-)\right)\Big{]}.\end{array}$
(3.3)

It follows from (3.1) that

$\begin{array}[]{ll}&\\!\\!\\!\displaystyle{\mathbb{E}}\int_{0}^{T_{N}}e^{-\delta
s}({\cal
L}-\delta)\Phi\left(X(s),\Lambda(s)\right)ds-{\mathbb{E}}\int_{0}^{T_{N}}e^{-\delta
s}\Phi^{\prime}\left(X(s),\Lambda(s)\right)f\big{(}X(s),C(s)\big{)}ds\\\
&\\!\\!\\!\displaystyle\quad\leq-{\mathbb{E}}\int_{0}^{T_{N}}e^{-\delta
s}p\left(X(s),C(s),\Lambda(s)\right)ds.\end{array}$ (3.4)

For each $s\in[0,T_{N}]$, by the mean value theorem, there exists a point
$\widetilde{X}(s)$ between $X(s)$ and $X(s-)$ such that

$\Phi\left(X(s),\Lambda(s-)\right)-\Phi\big{(}X(s-),\Lambda(s-)\big{)}=-\Delta
Y(s)\Phi^{\prime}\big{(}\widetilde{X}(s),\Lambda(s-)\big{)}+\Delta
Z(s)\Phi^{\prime}\big{(}\widetilde{X}(s),\Lambda(s-)\big{)}.$

By (3.1), we have

$\begin{array}[]{ll}\Phi\big{(}X(s),\Lambda(s-)\big{)}-\Phi\big{(}X(s-),\Lambda(s-)\big{)}\leq-q\Delta
Y(s)+r\Delta Z(s).\end{array}$ (3.5)

It follows from (3.3), (3.4), (3.5), and the nonnegativity of $\Phi(\cdot)$
that

$\begin{array}[]{ll}\Phi(x,\alpha)&\\!\\!\\!\displaystyle\geq{\mathbb{E}}\int_{0}^{T_{N}}e^{-\delta
s}p\big{(}X(s),\Lambda(s),C(s)\big{)}ds\\\
&\\!\\!\\!\displaystyle\qquad+{\mathbb{E}}\int_{0}^{T_{N}}e^{-\delta
s}qdY(s)-{\mathbb{E}}\int_{0}^{T_{N}}e^{-\delta s}rdZ(s).\end{array}$

Letting $N\to\infty$, it follows from (3.2) and the bounded convergence
theorem that $\Phi(x,\alpha)\geq J(x,\alpha,\Psi).$ Taking supremum over all
$\Psi\in\mathcal{A}_{x,\alpha}$, we obtain $\Phi(x,\alpha)\geq V(x,\alpha)$.
The conclusion follows. $\qquad\Box$

Using Lemma 3.1, we proceed to present an easily verifiable condition for the
finiteness of the value function.

######  Theorem 3.2.

Let (A) be satisfied. Moreover, suppose that there is a positive constant $K$
such that

$b(x,\alpha)\leq\delta
x+K\quad\text{for}\quad(x,\alpha)\in\SS\times{\cal{M}}.$ (3.6)

Then there exist a positive constant $M$ such that

$V(x,\alpha)\leq qx+M\quad\text{for}\quad(x,\alpha)\in\SS\times{\cal{M}}.$

###### Proof.

Define

$\Phi(x,\alpha)=qx+M\quad\text{for}\quad(x,\alpha)\in\SS\times{\cal{M}},$

where $M$ is a positive number to be specified. It can be seen that
$\Phi(\cdot)$ solves the system of inequalities (3.1) for sufficiently large
$M.$ By virtue of Lemma 3.1, $V(x,\alpha)\leq\Phi(x,\alpha)$ for any
$(x,\alpha)\in\SS\times{\cal{M}}.$ The details are omitted for brevity.
$\qquad\Box$

Next, we establish the continuity of the value function.

######  Theorem 3.3.

Let (A) be satisfied. Then the following assertions hold.

  * (a)

For any $x,y\in\SS$ and $\alpha\in{\cal{M}}$,

$V(x,\alpha)\geq q(x-y)^{+}-r(y-x)^{+}+V(y,\alpha).$ (3.7)

  * (b)

$V(\cdot)$ is Lipschitz continuous on $\SS\times{\cal{M}}$.

###### Proof.

(a) Fix $\Psi=(C,Y,Z)\in\mathcal{A}_{y,\alpha}$. Define

$\widetilde{C}(t)=C(t),\quad\widetilde{Y}(t)=Y(t)+(x-y)^{+},\quad\widetilde{Z}(t)=Z(t)+(y-x)^{+},\quad
t\geq 0,$

and $\widetilde{\Psi}=(\widetilde{C},\widetilde{Y},\widetilde{Z})$. Then
$\widetilde{\Psi}\in\mathcal{A}_{x,\alpha}$ and

$J(x,\alpha,\widetilde{\Psi})=q(x-y)^{+}-r(y-x)^{+}+J(y,\alpha,\Psi).$

Since $V(x,\alpha)\geq J(x,\alpha,\widetilde{\Psi})$, we have

$V(x,\alpha)\geq q(x-y)^{+}-r(y-x)^{+}+J(y,\alpha,\Psi),$

from which, (3.7) follows by taking supremum over
$\Psi\in\mathcal{A}_{y,\alpha}$.

(b) Similar to (3.7), we have

$V(y,\alpha)\geq q(y-x)^{+}-r(x-y)^{+}+V(x,\alpha).$ (3.8)

In view of (3.7), (3.8), for any $x,y\in\SS$ and $\alpha\in{\cal{M}}$,

$\begin{array}[]{ll}|V(x,\alpha)-V(y,\alpha)|&\\!\\!\\!\displaystyle\leq\big{(}|q|+|r|\big{)}|x-y|.\end{array}$

Thus, $V(\cdot)$ is Lipschitz continuous on $\SS\times{\cal{M}}$. $\qquad\Box$

For population models given by diffusion processes, under optimal or near-
optimal harvesting strategies, one should keep the population sizes in a
bounded set. In other words, if the initial population is too high, an
impulsive harvesting should be performed instantaneously; see [1, 7, 8]. We
proceed to provide a proof of that result for a general setting of a
controlled regime-switching diffusions with a mixed regular-singular control.

######  Theorem 3.4.

Suppose that there exists a number $U>\lambda$ such that

$q\Big{(}b(x,\alpha)-\delta(x-U)\Big{)}+\sup\limits_{c\in\mathcal{U}}\Big{(}a_{2}f(x,c)-g(x,\alpha,c)\Big{)}<0,\quad(x,\alpha,c)\in(U,\infty)\times{\cal{M}}\times\mathcal{U}.$
(3.9)

Then for each $(x,\alpha)\in(U,\infty)\times{\cal{M}}$,

$V(x,\alpha)=V(U,\alpha)+q(x-U).$

######  Remark 3.5.

Suppose $f(\cdot,c)=c$ for any $c\in\mathcal{U}$. If

$\limsup\limits_{x\to\infty}b(x,\alpha)<-\sup\limits_{c\in\mathcal{U}}(a_{2}c)/q,$

then we can take

$U=\sup\\{x:b(x,\alpha)\geq-\sup\limits_{c\in\mathcal{U}}(a_{2}c)/q\\}.$

###### Proof.

Fix some $(x,\alpha)\in(U,\infty)\times{\cal{M}}$,
$\Psi=(C,Y,Z)\in\mathcal{A}_{x,\alpha}$, and let $X$ be the corresponding
harvested process. Let $\varepsilon\in(0,1)$ be a constant and define

$\Phi(y,\alpha)=q(y-U)+\varepsilon\quad\text{for}\quad(y,\alpha)\in[U,\infty)\times{\cal{M}}.$
(3.10)

We can extend $\Phi(\cdot)$ to the entire $\SS\times{\cal{M}}$ so that
$\Phi(\cdot,\alpha)$ is a $C^{2}$ function for each $\alpha\in{\cal{M}}$, and
$\Phi(y,\alpha)>0$ for all $(y,\alpha)\in\SS\times{\cal{M}}$. By assumption
(3.9), we can check that

$(\mathcal{L}-\delta)\Phi(y,\alpha)+a_{2}f(y,c)-g(x,\alpha,c)<0,\quad(y,\alpha,c)\times[U,\infty)\times{\cal{M}}\times\mathcal{U}.$

Thus,

$\begin{array}[]{ll}(\mathcal{L}-\delta)\Phi(y,\alpha)-qf(y,c)&\\!\\!\\!\displaystyle<-\big{[}a_{2}f(y,c)-g(y,\alpha,c)\big{]}-qf(y,c)\\\
&\\!\\!\\!\displaystyle=-\big{[}(a_{2}+q)f(y,c)-g(y,\alpha,c)\big{]}\\\
&\\!\\!\\!\displaystyle=-p(y,\alpha,c),\quad(y,\alpha,c)\times[U,\infty)\times{\cal{M}}\times\mathcal{U}.\end{array}$
(3.11)

For an integer $N$ satisfying $N>U$, we define

$\tau_{N}=\inf\\{t\geq 0:X(t)\geq
N\\},\quad\widetilde{\gamma}_{U}=\inf\\{t\geq 0:X(t)\leq U\\},\quad
T_{N}=N\wedge\tau_{N}\wedge\widetilde{\gamma}_{U}.$

We have $T_{N}\to\widetilde{\gamma}_{U}$ almost surely as $N\to\infty$. Note
that $\widetilde{\gamma}_{U}\leq\tau$. By Dynkin’s formula,

$\begin{array}[]{ll}{\mathbb{E}}&\\!\\!\\!\displaystyle\big{[}e^{-\delta
T_{N}}\Phi\left(X(T_{N}),\Lambda(T_{N})\right)\big{]}-\Phi(x,\alpha)={\mathbb{E}}\int_{0}^{T_{N}}e^{-\delta
s}({\cal L}-\delta)\Phi\left(X(s),\Lambda(s)\right)ds\\\
&\\!\\!\\!\displaystyle\qquad-{\mathbb{E}}\int_{0}^{T_{N}}e^{-\delta
s}\Phi^{\prime}\left(X(s),\Lambda(s)\right)f\big{(}X(s),C(s)\big{)}ds\\\
&\\!\\!\\!\displaystyle\qquad+{\mathbb{E}}\sum\limits_{0\leq s\leq
T_{N}}e^{-\delta
s}\Big{[}\Phi\left(X(s),\Lambda(s-)\right)-\Phi\left(X(s-),\Lambda(s-)\right)\Big{]}.\end{array}$
(3.12)

For each $s\in[0,T_{N}]$, we have

$\begin{array}[]{ll}\Phi\left(X(s),\Lambda(s-)\right)-\Phi\left(X(s-),\Lambda(s-)\right)&\\!\\!\\!\displaystyle=-q\Delta
Y(s)+q\Delta Z(s)\\\ &\\!\\!\\!\displaystyle\leq-q\Delta Y(s)+r\Delta
Z(s).\end{array}$ (3.13)

We obtain from (3.11), (3.12), and (3.13) that

$\begin{array}[]{ll}{\mathbb{E}}&\\!\\!\\!\displaystyle\big{[}e^{-\delta
T_{N}}\Phi\left(X(T_{N}),\Lambda(T_{N})\right)\big{]}-\Phi(x,\alpha)\leq-{\mathbb{E}}\int_{s_{0}}^{T_{N}}e^{-\delta
s}p\big{(}X(s),\Lambda(s),C(s)\big{)}ds\\\
&\\!\\!\\!\displaystyle\qquad-{\mathbb{E}}\sum\limits_{0\leq s\leq
T_{N}}e^{-\delta s}q\Delta Y(s)+{\mathbb{E}}\sum\limits_{0\leq s\leq
T_{N}}e^{-\delta s}r\Delta Z(s).\end{array}$ (3.14)

Since $\Phi(y,\alpha)>0$ for any $(y,\alpha)\in\SS\times{\cal{M}}$, it follows
from (3.14) that

$\begin{array}[]{ll}&\\!\\!\\!\displaystyle{\mathbb{E}}\int_{0}^{T_{N}}e^{-\delta
s}p\big{(}X(s),\Lambda(s),C(s)\big{)}ds\\\
&\\!\\!\\!\displaystyle\qquad\qquad+{\mathbb{E}}\sum\limits_{0\leq s\leq
T_{N}}e^{-\delta s}q\Delta Y(s)-{\mathbb{E}}\sum\limits_{0\leq s\leq
T_{N}}e^{-\delta s}r\Delta Z(s)\\\
&\\!\\!\\!\displaystyle\qquad\leq\Phi(x,\alpha).\end{array}$

Letting $N\to\infty$, we obtain

$\begin{array}[]{ll}&\\!\\!\\!\displaystyle{\mathbb{E}}\int_{0}^{\widetilde{\gamma}_{U}}e^{-\delta
s}p\big{(}X(s),\Lambda(s),C(s)\big{)}ds\\\
&\\!\\!\\!\displaystyle\qquad\qquad+{\mathbb{E}}\sum\limits_{0\leq
s\leq\widetilde{\gamma}_{U}}e^{-\delta s}q\Delta
Y(s)-{\mathbb{E}}\sum\limits_{0\leq s\leq\widetilde{\gamma}_{U}}e^{-\delta
s}r\Delta Z(s)\\\ &\\!\\!\\!\displaystyle\qquad\leq\Phi(x,\alpha).\end{array}$
(3.15)

As a result,

$\begin{array}[]{ll}J(x,\alpha,\Psi)&\\!\\!\\!\displaystyle\leq
V(U,\alpha)+\Phi(x,\alpha)\\\
&\\!\\!\\!\displaystyle=V(U,\alpha)+q(x-U)+\varepsilon.\end{array}$

Letting $\varepsilon\to 0$ yields

$\begin{array}[]{ll}J(x,\alpha,\Psi)&\\!\\!\\!\displaystyle\leq
V(U,\alpha)+q(x-U).\end{array}$ (3.16)

Since (3.16) holds for any $\Psi\in\mathcal{A}_{x,\alpha}$, then

$V(x,\alpha)\leq V(U,\alpha)+q(x-U).$ (3.17)

On the other hand, it is obvious (by harvesting instantaneously $x-U$ at time
$t=0$) that

$V(x,\alpha)\geq V(U,\alpha)+q(x-U).$ (3.18)

In view of (3.17) and (3.18),

$V(x,\alpha)=V(U,\alpha)+q(x-U).$

The conclusion follows. $\qquad\Box$

We proceed to discuss the impact of large white noise. By using Lemma 3.1, we
construct an upper bound for the value function. Then we show that the value
function approaches the upper bound as the white noise intensity approaches
infinity. The following result is motivated by the present of multiplicative
noise.

######  Theorem 3.6.

Suppose that

$\lambda>0\quad\text{or}\quad\sup\limits_{\alpha\in{\cal{M}}}b(0,\alpha)\leq
0.$ (3.19)

Moreover, there exist positive constants $\beta\in(0,1)$, $K$, and $N$ such
that

$xb(x,\alpha)\leq K(1+|x|^{2}),\quad b(x,\alpha)q-\delta q(x-\lambda)\leq
K(x^{\beta}+1),$ (3.20)

and

$|\sigma(x,\alpha)|\geq
Nx\quad\text{for}\quad(x,\alpha)\in\SS\times{\cal{M}}.$

Then

$\lim\limits_{N\to\infty}V(x,\alpha)=q(x-\lambda),$

uniformly on $[\lambda,M]\times{\cal{M}}$ for any positive constant
$M>\lambda$.

###### Proof.

Let $M>\lambda$. For a fixed $\varepsilon\in(0,r-q)$, let $K_{\varepsilon}>0$
be sufficiently large such that

$\dfrac{(M+1)^{\beta}}{K_{\varepsilon}}\leq\varepsilon,\quad\dfrac{\beta}{K_{\varepsilon}}\leq
r-q-\varepsilon.$

Define

$\Phi(x,\alpha)=q(x-\lambda)+\dfrac{(x+1)^{\beta}}{K_{\varepsilon}},\quad(x,\alpha)\in\SS\times{\cal{M}}.$

Detailed computations lead to

$\begin{array}[]{ll}&\displaystyle\Phi^{\prime}(x,\alpha)=q+\dfrac{\beta}{K_{\varepsilon}(x+1)^{1-\beta}},\\\
&\displaystyle\Phi^{\prime\prime}(x,\alpha)=-\dfrac{\beta(1-\beta)}{K_{\varepsilon}(x+1)^{2-\beta}},\quad(x,\alpha)\in\SS\times{\cal{M}}.\end{array}$

It is clear that

$q-\Phi^{\prime}(x,\alpha)\leq 0,\quad\Phi^{\prime}(x,\alpha)-r\leq
0\quad\text{for}\quad(x,\alpha)\in\SS\times{\cal{M}}.$ (3.21)

Moreover, we have

$\begin{array}[]{ll}({\cal L}-\delta)\Phi(x,\alpha)&\\!\\!\\!\displaystyle\leq
b(x,\alpha)q+\dfrac{\beta}{K_{\varepsilon}(x+1)^{1-\beta}}b(x,\alpha)\\\
&\\!\\!\\!\displaystyle-\dfrac{\beta(1-\beta)N^{2}x^{2}}{2K_{\varepsilon}(x+1)^{2-\beta}}-\delta
q(x-\lambda)-\dfrac{\delta(x+1)^{\beta}}{K_{\varepsilon}}.\end{array}$ (3.22)

By assumption (3.19), there exist positive numbers $\lambda_{1}>\lambda$ and
$N_{1}$ such that for any $N\geq N_{1}$,

$b(x,\alpha)q+\dfrac{\beta}{K_{\varepsilon}(x+1)^{1-\beta}}b(x,\alpha)-\dfrac{\beta(1-\beta)N^{2}x^{2}}{2K_{\varepsilon}(x+1)^{2-\beta}}-\dfrac{\delta}{K_{\varepsilon}}<0,\quad(x,\alpha)\in[\lambda,\lambda_{1})\times{\cal{M}}.$

We obtain from (3.22) that

$({\cal L}-\delta)\Phi(x,\alpha)\leq
0\quad\text{for}\quad(x,\alpha)\in[\lambda,\lambda_{1})\times{\cal{M}}.$
(3.23)

For $x\geq\lambda_{1}>0$, it follows from (3.20) and (3.22) that

$\begin{array}[]{ll}({\cal L}-\delta)\Phi(x,\alpha)&\\!\\!\\!\displaystyle\leq
K(x^{\beta}+1)+\dfrac{\beta}{K_{\varepsilon}(x+1)^{1-\beta}}b(x,\alpha)-\dfrac{\beta(1-\beta)N^{2}x^{2}}{K_{\varepsilon}(x+1)^{2-\beta}}\\\
&\\!\\!\\!\displaystyle\leq\Big{[}K+\dfrac{\beta
K}{K_{\varepsilon}}\Big{(}1+\dfrac{1}{\lambda_{1}^{2}}\Big{)}-\dfrac{\beta(1-\beta)N^{2}x^{2-\beta}}{2K_{\varepsilon}(x+1)^{2-\beta}}\Big{]}x^{\beta}+K.\end{array}$

Hence there exists a positive number $N_{2}$ such that for any $N\geq N_{2}$,

$({\cal L}-\delta)\Phi(x,\alpha)\leq
0\quad\text{for}\quad(x,\alpha)\in[\lambda_{1},\infty)\times{\cal{M}}.$ (3.24)

Combining (3.21), (3.23), and (3.24) together with Lemma 3.1, we obtain

$V(x,\alpha)\leq
q(x-\lambda)+\dfrac{(x+1)^{\beta}}{K_{\varepsilon}},\quad(x,\alpha)\in\SS\times{\cal{M}}$

provided that $N\geq\max\\{N_{1},N_{2}\\}$. In particular, since
$\dfrac{(M+1)^{\beta}}{K_{\varepsilon}}\leq\varepsilon$, then for
$N\geq\max\\{N_{1},N_{2}\\}$, we have

$V(x,\alpha)\leq
q(x-\lambda)+\varepsilon\quad\text{for}\quad(x,\alpha)\in[\lambda,M]\times{\cal{M}}.$
(3.25)

Let $\Psi=(C,Y,Z)\in\mathcal{A}_{x,\alpha}$ given by

$C(t)=0,\quad Y(t)=x-\lambda,\quad Z(t)=0\quad\text{for}\quad t\geq 0.$

Then $J(x,\alpha,\Psi)=q(x-\lambda)$. It follows that $V(x,\alpha)\geq
q(x-\lambda)$. This together with (3.25) yields that

$0\leq
V(x,\alpha)-q(x-\lambda)\leq\varepsilon\quad\text{for}\quad(x,\alpha)\in[\lambda,M]\times{\cal{M}},$

which leads to the desired conclusion. $\qquad\Box$

We have a similar result regarding to the case of additive noise.

######  Theorem 3.7.

Suppose that

$\lambda>0\quad\text{or}\quad\sup\limits_{\alpha\in{\cal{M}}}b(0,\alpha)\leq
0.$ (3.26)

Moreover, there exist positive constants $K$ and $N$ such that

$xb(x,\alpha)\leq K(1+|x|^{2}),\quad b(x,\alpha)q-\delta q(x-\lambda)\leq K,$
(3.27)

and

$|\sigma(x,\alpha)|\geq N\quad\text{for}\quad(x,\alpha)\in\SS\times{\cal{M}}.$

Then

$\lim\limits_{N\to\infty}V(x,\alpha)=q(x-\lambda),$

uniformly on $[\lambda,M]\times{\cal{M}}$ for any positive constant
$M>\lambda$.

###### Proof.

Let $M>\lambda$. For a fixed $\varepsilon\in(0,r-q)$, let $K_{\varepsilon}>0$
be sufficiently large such that $\dfrac{1}{K_{\varepsilon}}\leq\varepsilon$
and then let $\beta\in(0,1)$ be such that

$\dfrac{\beta}{K_{\varepsilon}}\leq r-q-\varepsilon,\quad\beta
K<\dfrac{\delta}{K_{\varepsilon}}.$

Define

$\Phi(x,\alpha)=q(x-\lambda)+\dfrac{(x+1)^{\beta}}{K_{\varepsilon}}\quad\text{for}\quad(x,\alpha)\in\SS\times{\cal{M}}.$

Similar to the proof of Theorem 3.4, we have

$\begin{array}[]{ll}&\displaystyle\Phi^{\prime}(x,\alpha)=q+\dfrac{\beta}{K_{\varepsilon}(x+1)^{1-\beta}},\\\
&\displaystyle\Phi^{\prime\prime}(x,\alpha)=-\dfrac{\beta(1-\beta)}{K_{\varepsilon}(x+1)^{2-\beta}},\quad(x,\alpha)\in\SS\times{\cal{M}}.\end{array}$

we can check that

$q-\Phi^{\prime}(x,\alpha)\leq 0,\quad\Phi^{\prime}(x,\alpha)-r\leq
0\quad\text{for}\quad(x,\alpha)\in\SS\times{\cal{M}}.$ (3.28)

and

$\begin{array}[]{ll}({\cal L}-\delta)\Phi(x,\alpha)&\\!\\!\\!\displaystyle\leq
b(x,\alpha)q+\dfrac{\beta}{K_{\varepsilon}(x+1)^{1-\beta}}b(x,\alpha)\\\
&\\!\\!\\!\displaystyle-\dfrac{\beta(1-\beta)N^{2}}{2K_{\varepsilon}(x+1)^{2-\beta}}-\delta
q(x-\lambda)-\dfrac{\delta(x+1)^{\beta}}{K_{\varepsilon}}.\end{array}$ (3.29)

By assumption (3.26), there exist positive numbers $\lambda_{1}>\lambda$ and
$N_{1}$ such that for any $N\geq N_{1}$,

$b(x,\alpha)q+\dfrac{\beta}{K_{\varepsilon}(x+1)^{1-\beta}}b(x,\alpha)-\dfrac{\beta(1-\beta)N^{2}}{2K_{\varepsilon}(x+1)^{2-\beta}}-\dfrac{\delta}{K_{\varepsilon}}<0,\quad(x,\alpha)\in[\lambda,\lambda_{1})\times{\cal{M}}.$

We obtain from (3.29) that for $N\geq N_{1}$,

$({\cal L}-\delta)\Phi(x,\alpha)\leq
0\quad\text{for}\quad(x,\alpha)\in[\lambda,\lambda_{1})\times{\cal{M}}.$

For $x\geq\lambda_{1}>0$, it follows from (3.27) and (3.29) that

$\begin{array}[]{ll}({\cal L}-\delta)\Phi(x,\alpha)&\\!\\!\\!\displaystyle\leq
K+\dfrac{\beta}{K_{\varepsilon}(x+1)^{1-\beta}}b(x,\alpha)-\dfrac{\beta(1-\beta)N^{2}}{2K_{\varepsilon}(x+1)^{2-\beta}}-\dfrac{\delta(x+1)^{\beta}}{K_{\varepsilon}}\\\
&\\!\\!\\!\displaystyle\leq K+\beta
K\Big{(}x^{\beta}+\dfrac{1}{\lambda_{1}^{1-\beta}}\Big{)}-\dfrac{\beta(1-\beta)N^{2}}{2K_{\varepsilon}(x+1)^{2-\beta}}-\dfrac{\delta(x+1)^{\beta}}{K_{\varepsilon}}\\\
&\\!\\!\\!\displaystyle\leq K+\dfrac{\beta
K}{\lambda_{1}^{1-\beta}}-\Big{(}\delta/K_{\varepsilon}-\beta
K\Big{)}x^{\beta}-\dfrac{\beta(1-\beta)N^{2}}{2K_{\varepsilon}(x+1)^{2-\beta}}.\end{array}$

Let $\lambda_{2}>\lambda_{1}$ be such that

$K+\dfrac{\beta K}{\lambda^{1-\beta}}-\Big{(}\delta/K_{\varepsilon}-\beta
K\Big{)}x^{\beta}\leq
0\quad\text{for}\quad(x,\alpha)\in[\lambda_{2},\infty)\times{\cal{M}}.$

There exists a positive number $N_{2}$ such that for any $N\geq N_{2}$, we
have

$K+\dfrac{\beta
K}{\lambda_{1}^{1-\beta}}-\dfrac{\beta(1-\beta)N^{2}}{2K_{\varepsilon}(x+1)^{2-\beta}}\leq
0\quad\text{for}\quad(x,\alpha)\in[\lambda_{1},\lambda_{2}].$

Hence for $N\geq N_{2}$,

$({\cal L}-\delta)\Phi(x,\alpha)\leq
0\quad\text{for}\quad(x,\alpha)\in[\lambda_{1},\infty)\times{\cal{M}}.$ (3.30)

Combining (3.28) and (3.30), together with Lemma 3.1, we obtain

$V(x,\alpha)\leq
q(x-\lambda)+\dfrac{(x+1)^{\beta}}{K_{\varepsilon}}\quad\text{for}\quad(x,\alpha)\in\SS\times{\cal{M}}$

provided that $N\geq\max\\{N_{1},N_{2}\\}$. In particular, since
$\dfrac{(M+1)^{\beta}}{K_{\varepsilon}}\leq\varepsilon$, for
$N\geq\max\\{N_{1},N_{2}\\}$,

$V(x,\alpha)\leq
q(x-\lambda)+\varepsilon\quad\text{for}\quad(x,\alpha)\in[\lambda,M]\times{\cal{M}}.$
(3.31)

Let $\Psi=(C,Y,Z)\in\mathcal{A}_{x,\alpha}$ given by

$C(t)=0,\quad Y(t)=x-\lambda,\quad Z(t)=0\quad\text{for}\quad t\geq 0.$

Then $J(x,\alpha,\Psi)=q(x-\lambda)$. It follows that $V(x,\alpha)\geq
q(x-\lambda)$. This together with (3.31) gives us that

$0\leq
V(x,\alpha)-q(x-\lambda)\leq\varepsilon\quad\text{for}\quad(x,\alpha)\in[\lambda,M]\times{\cal{M}},$

which leads to the desired conclusion. $\qquad\Box$

##  4 Numerical Approximation

Formally, the associated Hamilton-Jacobi-Bellman equation of the underlying
problem is given by

$\begin{array}[]{ll}&\displaystyle\max\bigg{\\{}\mathcal{G}\Phi(x,\alpha),q-\Phi^{\prime}(x,\alpha),\Phi^{\prime}(x,\alpha)-r\bigg{\\}}=0\quad\text{for}\quad(x,\alpha)\in\SS\times{\cal{M}},\\\
&\displaystyle\Phi(x,\alpha)=0\quad\text{for}\quad(x,\alpha)\notin\SS\times{\cal{M}},\end{array}$
(4.1)

where

$\mathcal{G}\Phi(x,\alpha)=(\mathcal{L}-\delta)\Phi(x,\alpha)+\max\limits_{c\in\mathcal{U}}\Big{[}p(x,\alpha,c)-\Phi^{\prime}(x,\alpha)f(x,c)\Big{]}.$

A closed-form solution to (4.1) is virtually impossible to obtain. Thus, we
proceed with the Markov chain approximation method; see [11, 16, 17]. That is,
we construct a discrete time, finite-state, controlled Markov chain to
approximate the controlled switching diffusions. In [11], the authors applied
that method to solve a harvesting-type problem for a regime-switching
diffusion with a regular control and a singular control. In this paper, we
need to adopt and modify the approximation to fit the combination of a two-
sided singular control and regular control formulation. In view of Theorem
3.4, we only need to choose a large positive integer $U$ and compute the value
function on $\SS\cap[0,U]=[\lambda,U]$. With $x\in[\lambda,U]$, we can rewrite
(2.2) as

$\begin{array}[]{ll}&\\!\\!\\!\displaystyle
X(t)=x+\int\limits_{0}^{t}\big{(}b(X(s),\Lambda(s)\big{)}-f\big{(}X(s),C(s)\big{)}ds+\int\limits_{0}^{t}\sigma\big{(}X(s),\Lambda(s)\big{)}dw(s)-Y(t)+Z(t).\end{array}$
(4.2)

The payoff functional is

$\begin{array}[]{ll}&\displaystyle
J(x,\alpha,\Psi)={\mathbb{E}}\bigg{[}\int\limits_{0}^{\tau}e^{-\delta
s}p(X(s),\Lambda(s),C(s))ds+\int\limits_{0}^{\tau}e^{-\delta
s}qdY(s)-\int\limits_{0}^{\tau}e^{-\delta s}rdZ(s)\bigg{]}.\end{array}$ (4.3)

###  4.1 Approximating Markov chains

Let $h$ be a discretization parameter for $X(\cdot)$. Assume without loss of
generality that $U$ is a multiple of $h$. Define

$S_{h}:=[0,U]\cap\\{x\in[\lambda-h,U]:x=kh,k\in\mathbb{Z}_{+}\\}.$

Let $\\{(X^{h}_{n},\Lambda^{h}_{n}):n\in\mathbb{Z}_{+}\\}$ be a discrete-time
controlled Markov chain with state space $S_{h}\times{\cal{M}}$. For each $n$,
the increments of the chain $\Delta X^{h}_{n}=X^{h}_{n+1}-X^{h}_{n}$
approximates exactly one of the following quantities dynamically.

  * •

Diffusion step:
$\big{(}b(X(t),\Lambda(t))-f(X(t),C(t)\big{)}dt+\sigma(X(t),\Lambda(t))dw(t)$.

  * •

Impulsive harvesting step: $-dY(t)$.

  * •

Impulsive renewing step: $dZ(t)$.

We use $\pi^{h}_{n}$ to denote the type of step $n$. We set $\pi^{h}_{n}=0$ if
step $n$ is a diffusion step, $\pi^{h}_{n}=1$ if step $n$ is an impulsive
harvesting step, and $\pi^{h}_{n}=-1$ if step $n$ is an impulsive renewing
step. Each of these steps is described precisely in what follows.

(a) If $\pi^{h}_{n}=0$, a regular control $C^{h}_{n}$ is exercised so that
$\Delta X^{h}_{n}$ is to behave like an increment of $\int(b-f)dt+\sigma dw$
over a small time interval. (b) If $\pi^{h}_{n}=1$, a regular control
$C^{h}_{n}$ and an impulsive harvesting $\Delta Y^{h}_{n}=h$ are applied. (c)
If $\pi^{h}_{n}=-1$, a regular control $C^{h}_{n}$ and an impulsive renewing
$\Delta Z^{h}_{n}=h$ are applied.

Thus, the amounts of impulsive harvesting and impulsive renewing at any step
$n$ are $\Delta Y^{h}_{n}=hI_{\\{\pi^{h}_{n}=1\\}}$ and $\Delta
Z^{h}_{n}=hI_{\\{\pi^{h}_{n}=-1\\}}$, respectively. Let
$C^{h}=\\{C^{h}_{n}\\}$ be a sequence of regular controls. The space of
regular controls is given by $\mathcal{U}$. Let
$\pi^{h}=\\{\pi_{n}^{h}\\}_{n\geq 0}$ denote the sequence of control types. We
denote by $p^{h}\left((x,\alpha),(y,\ell)|\pi,c\right)$ the transition
probability from state $(x,\alpha)$ to another state $(y,\ell)$ under the
control type $\pi$ and regular control $c$. Denote
$\mathcal{F}^{h}_{n}=\sigma\\{X^{h}_{k},\Lambda^{h}_{k},\pi^{h}_{k},C^{h}_{k},k\leq
n\\}$.

The sequence $\\{(\pi^{h}_{n},C^{h}_{n})\\}$ is said to be admissible if it
satisfies the following conditions:

  * (a)

$\pi^{h}_{n},C^{h}_{n}$ are
$\sigma\\{X^{h}_{0},\dots,X^{h}_{n},\Lambda^{h}_{0},\dots,\Lambda^{h}_{n},\pi^{h}_{0},\dots,\pi^{h}_{n-1},C^{h}_{1},\dots,C^{h}_{n-1}\\}-\text{adapted},$

  * (b)

For any $(x,\alpha)\in S_{h}\times{\cal{M}}$, we have

$\begin{array}[]{ll}{\mathbb{P}}\\{X^{h}_{n+1}=x,\Lambda^{h}_{n+1}=\alpha|\mathcal{F}^{h}_{n}\\}&={\mathbb{P}}\\{X^{h}_{n+1}=x,\Lambda^{h}_{n+1}=\alpha|X^{h}_{n},\Lambda^{h}_{n},\pi^{h}_{n},C^{h}_{n}\\}\\\
&=p^{h}\big{(}\big{(}X^{h}_{n},\Lambda^{h}_{n}),(x,\alpha)|\pi^{h}_{n},C^{h}_{n}\big{)},\end{array}$

  * (c)

$\pi^{h}_{n}\in\\{0,\pm 1\\},C^{h}_{n}\in\mathcal{U}$, $X^{h}_{n}\in
S_{h},\Lambda^{h}_{n}\in{\cal{M}}$ for all $n\in\mathbb{Z}_{+}$.

The class of all admissible control sequences $(\pi^{h},C^{h})$ for initial
state $(x,\alpha)$ will be denoted by $\mathcal{A}^{h}_{x,\alpha}$.

For each $(x,\alpha,i,c)\in S_{h}\times{\cal{M}}\times\\{0,\pm
1\\}\times\mathcal{U}$, we define a family of the interpolation intervals
$\Delta t^{h}(x,\alpha,i,c)$. The values of $\Delta t^{h}(x,\alpha,i,c)$ will
be specified later. Then we define

$\begin{array}[]{ll}&\displaystyle t^{h}_{0}=0,\quad\Delta t^{h}_{k}=\Delta
t^{h}(X^{h}_{k},\Lambda^{h}_{k},\pi^{h}_{k},C^{h}_{k}),\quad
t^{h}_{n}=\sum\limits_{k=0}^{n-1}\Delta t^{h}_{k}.\\\ \end{array}$ (4.4)

For $(x,\alpha)\in S_{h}\times{\cal{M}}$ and
$(\pi^{h},C^{h})\in\mathcal{A}^{h}_{x,\alpha}$, the payoff functional for the
controlled Markov chain is defined as

$\begin{array}[]{ll}J^{h}(x,\alpha,\pi^{h},C^{h})&={\mathbb{E}}\sum\limits_{k=0}^{\eta_{h}-1}e^{-\delta
t_{k}^{h}}\bigg{\\{}p(X^{h}_{k},\Lambda^{h}_{k},C^{h}_{k})\Delta
t^{h}_{k}+q\Delta Y^{h}_{k}-r\Delta Z^{h}_{k}\bigg{\\}},\end{array}$ (4.5)

with $\eta_{h}=\inf\\{n\geq 0:X^{h}_{n}\notin\SS\\}$. The value function of
the controlled Markov chain is

$V^{h}(x,\alpha)=\sup\limits_{(\pi^{h},C^{h})\in\mathcal{A}^{h}_{x,\alpha}}J^{h}(x,\alpha,\pi^{h},C^{h}).$
(4.6)

The corresponding dynamic programming equation for the discrete approximation
is given by

$\begin{array}[]{ll}V^{h}(x,\alpha)=\max\limits_{(i,c)}\bigg{\\{}&\\!\\!\\!\displaystyle
e^{-\delta\Delta t^{h}(x,\alpha,i,c)}\sum\limits_{(y,\ell)\in
S_{h}}V^{h}(y,\ell)p^{h}\big{(}(x,\alpha),(y,\ell)|i,c\big{)}+p(x,\alpha,c)\Delta
t^{h}(x,\alpha,i,c),\\\
&\\!\\!\\!\displaystyle\qquad\qquad\big{(}V^{h}(x-h,\alpha)+qh\big{)}I_{\\{x>\lambda\\}},\big{(}V^{h}(x+h,\alpha)-rh\big{)}I_{\\{x<U\\}}\bigg{\\}},\end{array}$

for any $(x,\alpha)\in S_{h}\times{\cal{M}}.$ Note that $V^{h}(x,\alpha)=0$
for $(x,\alpha)\notin S_{h}\times{\cal{M}}$.

###  4.2 Transition probabilities

Let ${\mathbb{E}}^{h,\pi,c}_{x,\alpha,n}$,
${\mathbb{C}ov}^{h,\pi,c}_{x,\alpha,n}$ denote the conditional expectation and
covariance given by

$\\{X_{k}^{h},\Lambda_{k}^{h},\pi_{k}^{h},C^{h}_{k},k\leq
n-1,X_{n}^{h}=x,\Lambda_{n}^{h}=\alpha,\pi^{h}_{n}=\pi,C^{h}_{n}=c\\},$

respectively. Our objective is to define transition probabilities
$p^{h}\big{(}(x,\alpha),(y,\ell)|\pi,c\big{)}$ so that the controlled Markov
chain $\\{(X^{h}_{n},\Lambda^{h}_{n})\\}$ is locally consistent with respect
to $(X(\cdot),\Lambda(\cdot))$ in the sense that the following hold.

$\begin{array}[]{ll}&\displaystyle{\mathbb{E}}^{h,0,c}_{x,\alpha,n}\Delta
X_{n}^{h}={b}(x,\alpha)\Delta t^{h}(x,\alpha,0,c)+o(\Delta
t^{h}(x,\alpha,0,c)),\\\ &\displaystyle Cov^{h,0,c}_{x,\alpha,n}\Delta
X_{n}^{h}=\sigma^{2}(x,\alpha)\Delta t^{h}(x,\alpha,0,c)+o(\Delta
t^{h}(x,\alpha,0,c)),\\\ &\displaystyle\sup\limits_{n,\ \omega}|\Delta
X_{n}^{h}|\to 0\ \hbox{ as }\ h\to 0.\end{array}$ (4.7)

To this end, motivated by the procedure in [11, 17], we construct the
transition probabilities below. For $(x,\alpha)\in S_{h}\times{\cal{M}}$,
$\lambda\leq x<U$ and $c\in\mathcal{U}$, we define

$Q_{h}(x,\alpha,0,c)=\sigma^{2}(x,\alpha)+h|b(x,\alpha)-f(x,c)|-h^{2}q_{\alpha\alpha}+\zeta(h).$

We set $\zeta(h)=h$. If, in addition, $\sigma^{2}(x,\alpha)>0$ for any
$(x,\alpha)\in\SS\times{\cal{M}}$, we can simply take $\zeta(h)=0$. Then we
define

$\begin{array}[]{ll}&\displaystyle
p^{h}\left((x,\alpha),(x+h,\alpha)|0,c\right)=\dfrac{\sigma^{2}(x,\alpha)/2+\big{(}b(x,\alpha)-f(x,c)\big{)}^{+}h}{Q_{h}(x,\alpha,0,c)},\\\
&\displaystyle
p^{h}\left((x,\alpha),(x-h,\alpha)|0,c\right)=\dfrac{\sigma^{2}(x,\alpha)/2+\big{(}b(x,\alpha)-f(x,c)\big{)}^{-}h}{Q_{h}(x,\alpha,0,c)},\\\
&\displaystyle
p^{h}\left((x,\alpha),(x,\ell)|0,c\right)=\dfrac{h^{2}\Lambda_{\alpha\ell}}{Q_{h}(x,\alpha,0,c)}\text{
for }\alpha\neq\ell,\quad
p^{h}\left((x,\alpha),(x,\alpha)|0,c\right)=\dfrac{\zeta(h)}{Q_{h}(x,\alpha,0,c)},\\\
&\displaystyle\Delta
t^{h}(x,\alpha,0,c)=\dfrac{h^{2}}{Q_{h}(x,\alpha,0,c)},\\\ &\displaystyle
p^{h}\left((x,\alpha),(y,\ell)|0,c\right)=0\quad\text{for all not listed
values of $(y,\ell)\in S_{h}\times{\cal{M}}$},\end{array}$ (4.8)

At impulsive harvesting and renewing steps, we define

$\begin{array}[]{ll}&p^{h}\left((x,\alpha),(x-h,\alpha)|1,c\right)=1,\quad\Delta
t^{h}(x,\alpha,1,c)=0,\\\
&p^{h}\left((x,\alpha),(x+h,\alpha)|-1,c\right)=1,\quad\Delta
t^{h}(x,\alpha,-1,c)=0.\end{array}$ (4.9)

Thus, $p^{h}\left((x,\alpha),(y,\ell)|\pm 1,c\right)=0$ for all nonlisted
values of $(y,\ell)\in S_{h}\times{\cal{M}}$. The definition of
$p^{h}((x,\alpha),\cdot|i,c\big{)}$ for $(x,\alpha)\notin\SS\times{\cal{M}}$
is not important since in the analysis of the control problem, the chain will
be stopped at the first time $\\{X^{h}_{n}\\}$ exits $\SS$. Using the above
transition probabilities, we can check that the local consistency conditions
of $\\{(X^{h}_{n},\Lambda^{h}_{n})\\}$ in (4.7) are satisfied.

The convergence result is based on a continuous-time interpolation of the
chain and relaxed control representations. In addition, a “stretched-out” time
scale is introduced to overcome the possible non-tightness of the piecewise
constant interpolations associated with $Y(\cdot)$ and $Z(\cdot)$; see [11]
and [16]. Using weak convergence methods, we can obtain the convergence of the
value function. The main convergence result is given below. The proof is a
modification of that in [11], we omit it for brevity.

######  Theorem 4.1.

Let $V(x,\alpha)$ and $V^{h}(x,\alpha)$ be value functions defined in (2.4)
and (4.6), respectively. Then $V^{h}(x,\alpha)\to V(x,\alpha)$ as $h\to 0$.

##  5 Numerical Examples

We consider a stochastic population with harvesting and renewing given by

$dX(t)=b(X(t),\Lambda(t))dt+\sigma(X(t),\Lambda(t))dw(t)-C(t)dt-dY(t)+dZ(t),$

where

$\begin{array}[]{ll}&b(x,\alpha)=x(\alpha-1.5x),\quad\sigma(x,\alpha)=(\alpha/2)x.\end{array}$

Suppose that ${\cal{M}}=\\{1,2\\}$ and the generator $\Gamma$ of the Markov
chain $\alpha(\cdot)$ is given by

$\Gamma_{11}=-1,\quad\Gamma_{12}=1,\quad\Gamma_{21}=1.5,\quad\Gamma_{22}=-1.5.$

We also assume that

$\begin{array}[]{ll}&q=1,\quad r=2.25,\quad\delta=0.05,\\\
&p(x,\alpha,c)=\dfrac{3c}{2}-\dfrac{\alpha
c^{2}}{8(1+x)}\quad\text{for}\quad(x,\alpha,c)\in\mathbb{R}_{+}\times{\cal{M}}\times\mathcal{U},\\\
&f(x,c)=c\quad\text{for}\quad(x,c)\in\mathbb{R}_{+}\times\mathcal{U}.\\\
\end{array}$

The original price and costs are $a_{1}=3/2,a_{2}=1/2,a_{3}=0.75$ and the
regular control cost function is $g(x,\alpha,c)=\frac{\alpha c^{2}}{8(1+x)}$,
which has a quadratic form for a fixed $x$. By the state constraint, the time
horizon of the control problem is $[0,\tau]$, where $\tau=\inf\\{t\geq
0:X(t)<\lambda\\}$. The value of $\lambda$ and the control set $\mathcal{U}$
are to be set. Using the transition probabilities constructed in (4.8) and
(4.9), we carry out the computation by using value iteration and policy
iteration method. For each $(x,\alpha)\in S_{h}\times{\cal{M}}$, we use the
following notations for the $n$th iteration: $C^{h,n}(x,\alpha)$ is the
regular control, $\pi^{h,n}(x,\alpha)$ is the control type, and
$V^{h,n}(x,\alpha)$ is the value function. Initially, we take

$C^{h,0}(x,\alpha)=0,\quad\pi^{h,0}(x,\alpha)=1,\quad
V^{h,0}(x,\alpha)=q(x-\lambda)\quad\text{for}\quad(x,\alpha)\in
S_{h}\times{\cal{M}}.$

Note that $\big{(}\pi^{h,0}(x,\alpha),C^{h,0}(x,\alpha)\big{)}$ corresponds to
the control $\Psi=(C,Y,Z)\in\mathcal{A}_{x,\alpha}$ given by

$C(t)=0,\quad Y(t)=x-\lambda,\quad Z(t)=0\quad\text{for}\quad t\geq 0,$

of the controlled switching diffusion (2.2). Based on our computations of the
$n$th iteration, we define

$\begin{array}[]{ll}&\displaystyle
V^{h,n+1}_{1,c}(x,\alpha)=\big{(}V^{h,n}(x-h,\alpha)+qh\big{)}I_{\\{x>\lambda\\}},\\\
&\displaystyle
V^{h,n+1}_{-1,c}(x,\alpha)=\big{(}V^{h,n}(x+h,\alpha)-rh\big{)}I_{\\{x<U\\}},\\\
&\displaystyle V^{h,n+1}_{0,c}(x,\alpha)=e^{-\delta\Delta
t^{h}(x,0,c,\alpha)}\big{[}V^{h,n}(x+h,\alpha)p^{h}\left((x,\alpha),(x+h,\alpha)|0,c\right)\\\
&\displaystyle\qquad\qquad\qquad\qquad+V^{h,n}(x-h,\alpha)p^{h}\left((x,\alpha),(x-h,\alpha)|0,c\right)\\\
&\displaystyle\qquad\qquad\qquad\qquad+V^{h,n}(x,\alpha)p^{h}\left((x,\alpha),(x,\alpha)|0,c\right)\big{]}I_{\\{x<U\\}}+p(x,\alpha,c)\Delta
t^{h}(x,\alpha,i,c)I_{\\{x<U\\}}\end{array}$

We find an improved value $V^{h,n+1}(x,\alpha)$ and record the corresponding
controls by

$\begin{array}[]{ll}&\displaystyle\left(\pi^{h,n+1}(x,\alpha),C^{h,n+1}(x,\alpha)\right)=\hbox{argmax}\left\\{V^{h,n+1}_{i,c}(x,\alpha):i\in\\{-1,0,1\\},c\in\mathcal{U}\right\\},\\\
&\displaystyle
V^{h,n+1}(x,\alpha)=\max\left\\{V^{h,n+1}_{i,c}(x,\alpha):i\in\\{-1,0,1\\},c\in\mathcal{U}\right\\}.\end{array}$

The iterations stop as soon as the increment
$V^{h,n+1}(x,\alpha)-V^{h,n}(x,\alpha)$ reaches a pre-specified tolerance
level. We set the error tolerance to be $10^{-6}$.

(a)

Figure 1: The control type (left) and regular control (right) as functions of
$(x,\alpha)$ for $\lambda=0.2$ and $\mathcal{U}=\\{0\\}$ (Example 5.1) Figure
2: The value function as a function of $(x,\alpha)$ for $\lambda=0.2$ and
$\mathcal{U}=\\{0\\}$ (Example 5.1)

######  Example 5.1.

Let $\lambda=0.2$ and $\mathcal{U}=\\{0\\}$. In view of Theorem 3.4, we can
take $U=2$. Figure 1 shows the control type, regular control rate as functions
of population size $x$ and regime $\alpha$. The corresponding value function
is given in Figure 2. Since $\mathcal{U}=\\{0\\}$, the regular control does
not work. The controller can harvest or renew by exercising impulsive controls
only. It appears that in both regimes, if $x=\lambda=0.2$, an impulsive
renewing is performed. There is a level $U_{\alpha}$ depending on $\alpha$ so
that we should apply an impulsive harvesting whenever the population size is
above $U_{\alpha}$ and should not do harvest nor renew whenever the population
size is in $(\lambda,U_{\alpha})$. Thus, under this strategy, the population
size is in $[\lambda,U_{\alpha}]$ for $t\in(0,\infty)$.

(a)

(b)

Figure 3: The control type (left) and regular control (right) as functions of
$(x,\alpha)$ for $\lambda=0.2$ and $\mathcal{U}=\\{k\in\mathbb{Z}:-2\leq k\leq
3\\}$ (Example 5.2) Figure 4: The value function as a function of $(x,\alpha)$
for $\lambda=0.2$ and $\mathcal{U}=\\{k\in\mathbb{Z}:-2\leq k\leq 3\\}$
(Example 5.2)

######  Example 5.2.

Let $\lambda=0.2$ and $\mathcal{U}=\\{k\in\mathbb{Z}:-2\leq k\leq 3\\}$. In
view of Theorem 3.4, we can take $U=2$. Figure 3 shows the control type,
regular control rate as functions of population size $x$ and regime $\alpha$.
The corresponding value function is given in Figure 4. Compared to the
preceding example, the controller have more options to harvest or renew.
Therefore, the value function in Figure 4 is much larger than the preceding
one. Moreover, as shown in Figure 3, the larger the population is, the higher
the harvesting rate is.

(a)

(b)

Figure 5: The control type (left) and regular control (right) as functions of
$(x,\alpha)$ for $\lambda=0.2$ and $\mathcal{U}=\\{k\in\mathbb{Z}:-2\leq k\leq
3\\}$ (Example 5.3) Figure 6: The value function as a function of $(x,\alpha)$
for $\lambda=0.2$ and $\mathcal{U}=\\{k\in\mathbb{Z}:-2\leq k\leq 3\\}$
(Example 5.3)

######  Example 5.3.

Let $\lambda=0.4$ and $\mathcal{U}=\\{k\in\mathbb{Z}:-2\leq k\leq 3\\}$. In
view of Theorem 3.4, we can take $U=2$. Figure 5 shows the control type,
regular control rate as functions of population size $x$ and regime $\alpha$.
The corresponding value function is given in Figure 6. It can be seen that the
value function is smaller than that in Example 5.2. This observation fits the
fact that the one can harvest only if the population size is higher than
$\lambda=0.4$ compared to 0.2 in Example 5.2. In this and the previous
numerical experiments, we observe the same phenomena as follows: (a) if the
population size hits $\lambda$, an impulsive renewing is performed to keep the
species in $\SS=[\lambda,\infty)$, which is the harvesting-renewing domain;
(b) For each $\alpha$, there are levels $L^{(1)}_{\alpha},L^{(2)}_{\alpha}$,
and $L^{(3)}_{\alpha}$ so that we should renew with bounded rates if
$x\in(\lambda,L^{(1)}_{\alpha})$, we should not harvest nor renew if
$x\in[L^{(1)}_{\alpha},L^{(2)}_{\alpha})$, we should harvest with bounded
rates if $x\in[L^{(2)}_{\alpha},L^{(3)}_{\alpha})$, and we should perform an
impulsive harvesting if $x\in[L^{(3)}_{\alpha},U]$.

It can be seen from the numerical experiments that we should keep the
population size in the interval $[\lambda,U]$ for $t>0$. In other words, if
the initial population size $x>\lambda$, then $\tau=\infty$ almost surely.
Thus, compared with the well-known formulations for harvesting-type problems,
the proposed model offers an effective planning for the balance between
economical aspect and the sustainable purpose.

##  6 Further remarks

This paper has been devoted to the study of a generalized harvesting problem
for a stochastic population. We have established the finiteness and continuity
of the value function. Moreover, we have shown that for common population
systems, it is important to maintain the population size of species in a
bounded set. We also study the impact of large white noise on harvesting. A
numerical algorithm has been constructed by the Markov chain approximation
methods. The numerical experiments have revealed the effect of Markovian
switching and control costs associated with harvesting/renewing activities. It
has been observed that under a mixed singular-control formulation, the
decision makers have more options to harvest or renew the species, which is
much more beneficial than the known models with no renewing, or with only one
control component. Moreover, the state constraint provides an effective
strategy for sustainable purpose.

In this paper, we focus on single species. The development in this paper can
be carried over to the multidimensional case. However, in order to identify
the bounded set that the population sizes should be in, one needs to be
careful to handle the interaction between species. Although multidimensional
systems can be handled, the multidimensional structure added more difficulty
to study the impact of large white noise. The computation of multidimensional
ecosystems with many constraints and controls is another challenging task. In
addition, one can also consider the problem under other constraints, random
prices, and random cost functions.

## References

  * [1] L. H. R. Alvarez, L. A. Shepp, Optimal harvesting of stochastically fluctuating populations, J. Math. Biol.37 (1998), pp. 155–177.
  * [2] L. H. R. Alvarez, Singular stochastic control in the presence of a state-dependent yield structure, Stochastic Process. Appl. 86 (2000), pp. 323–343.
  * [3] R. Lande, S. Engen, B. Sæther, Optimal harvesting of fluctuating populations with a risk of extinction, Am. Naturalist.145(5) (1995), pp. 728–745.
  * [4] G. Ferrari, On a class of singular stochastic control problems for reflected diffusions, J. Math. Anal. Appl. 473(2) (2019), pp. 952–979.
  * [5] W. H. Fleming, H. M. Soner, Controlled Markov Processes and Viscosity Solutions, Springer-Verlag, New York, 1993.
  * [6] L.H.R. Alvarez, A.Hening, Optimal sustainable harvesting of populations in random environments, Stochastic Processes and their Applications, DOI: http://dx.doi.org/10.1016/j.spa.2019.02.008.
  * [7] A. Hening, K. Tran, T. T. Phan, G. Yin, Harvesting of interacting stochastic populations. J. Math. Biol. 79(2) (2019), pp. 533–570.
  * [8] A. Hening, K. Tran, Harvesting and seeding of stochastic populations: analysis and numerical approximation, J. Math. Biol. 81(1) (2020), pp. 65–112.
  * [9] U. Haussmann, W. Suo, Singular optimal stochastic controls. I. Existence, SIAM J. Control Optim. 33(3) (1995), pp. 916–936.
  * [10] U. Haussmann, W. Suo, Singular optimal stochastic controls. II. Dynamic programming, SIAM J. Control Optim. 33(3) (1995), pp. 937-959.
  * [11] Z. Jin, G. Yin, C. Zhu, Numerical solutions of optimal risk control and dividend optimization policies under a generalized singular control formulation, Automatica 48(8) (2012), pp.1489–1501.
  * [12] I. Karatzas, Ioannis S.E. Shreve, Connections between optimal stopping and singular stochastic control. I. Monotone follower problems, SIAM J. Control Optim. 22(6) (1984), pp. 856–877.
  * [13] I. Karatzas, Ioannis S.E. Shreve, Connections between optimal stopping and singular stochastic control. II. Reflected follower problems, SIAM J. Control Optim. 23(3) (1985), pp. 433–451.
  * [14] I. Kharroubi, T. Lim, T. Mastrolia, Regulation of renewable resource exploitation, SIAM J. Control Optim. 58(1) (2020), pp. 551–579.
  * [15] K. Kunwai, F. Xi, G. Yin, C. Zhu, On an Ergodic Two-Sided Singular Control Problem, 2021, arXiv preprint arXiv:2108.05500.
  * [16] H. J. Kushner and L. F. Martins, Numerical methods for stochastic singular control problems, SIAM J. Control Optim. 29(6) (1991), pp. 1443–1475.
  * [17] H. J. Kushner, P. G. Dupuis, Numerical methods for stochastic control problems in continuous time, Springer-Verlag, 1992.
  * [18] G. Liang, M. Zervos, Ergodic singular stochastic control motivated by the optimal sustainable exploitation of an ecosystem, 2020, arXiv preprint arXiv:2008.05576.
  * [19] E. M. Lungu, B. Øksendal, Optimal harvesting from a population in a stochastic crowded environment, Mathematical Biosciences, 145(1) (1997), pp. 47–75.
  * [20] X. Mao, C. Yuan, Stochastic differential equations with Markovian switching, Imperial College Press, 2006.
  * [21] D.H. Nguyen, G. Yin, Sustainable harvesting policies under long-run average criteria: near optimality, Appl. Math. Optim. 81(2) (2020), pp. 443–478.
  * [22] D. Nguyen, N. Nguyen, and G. Yin, Stochastic functional Kolmogorov equations I: Persistence, to appear in Stochastic Process Appl.
  * [23] H. Pham, Continuous-Time Stochastic Control and Optimization with Financial Applications, Stoch. Model. Appl. Probab. 61, Springer-Verlag, Berlin, 2009.
  * [24] R. Radner, Larry Shepp, Risk vs. profit potential: A model for corporate strategy, Journal of Economic Dynamics and Control 20 (8), 1996, pp. 1373–1393.
  * [25] Q. Song, R.H. Stockbridge, C. Zhu, On optimal harvesting problems in random environments. SIAM J. Control Optim. 49(2) (2011), pp. 859–889.
  * [26] K. Tran, G. Yin, Optimal harvesting strategies for stochastic ecosystems, IET Control Theory & Applications, 11(15) (2017), pp. 2521–2530.
  * [27] K. Tran, Optimal exploitation for hybrid systems of renewable resources under partial observation, Nonlinear Anal. Hybrid Syst. 40 (2021), 101013.
  * [28] G. Yin and C. Zhu, Hybrid Switching Diffusions: Properties and Applications, Springer, New York, 2010.
  * [29]