A simplified stochastic optimization model for logistic dynamics with control-dependent carrying capacity

ABSTRACT

A simplified stochastic control model for optimization of logistic dynamics with the control-dependent carrying capacity, which is motivated by a recent algae population management problem in the river environment, is presented. Solving the optimization problem reduces to finding a solution to a non-local first-order differential equation called tt–Jacobi–Bellman (HJB) equation. It is shown that the HJB equation has a unique viscosity solution and that the solution can be approximated with a finite difference scheme. Asymptotic estimates of the solution and the optimal control are derived and compared with numerical solutions. Finally, parameter dependence of the optimal control is examined numerically with implications to river environmental management.

KEYWORDS:

• Logistic-type equation
• stochastic optimal control
• Hamilton–Jacobi–Bellman equation
• viscosity solution
• benthic algae

MSC CODE:

• 34H05
• 35D40
• 65M06
• 92D25
• 93E20

1. Introduction

This paper focuses on a new stochastic control model and its application to algae management in rivers, in which modelling and control of population dynamics play central roles. Not limited to this problem, mathematical modelling of population dynamics in natural environment is a heart of ecological management [52]. Recent examples include invasive species management [28], tracking predator-prey dynamics for biological conservation [50], and sustainable fisheries management [19]. The population dynamics is often inherently stochastic mainly due to unresolved nonlinearities involved in the ecological and environmental processes [31]. Stochastic process models like the stochastic differential equations (SDEs) [42] are adequate mathematical tools for modelling and control of such dynamics.

SDEs serve as simple mathematical models for a variety of population dynamics and related environmental and ecological problems. Mandal et al. [38] considered phytoplankton dynamics with a system of nonlinear SDEs for analysis of noise-dependent population extinction. Mau et al. [39] presented a multiplicative jump process model for minimum mathematical description of salt and contaminants dynamics in soil environment. Sanz and Alonso [48] discussed an approximate aggregation technique for simplifying an SDE-based heterogeneous population dynamics model. Allen [3] derived a stochastic partial differential equation that governs size- and age-dependent population dynamics in noisy environment from a discretized SDE counterpart.

Management of population dynamics can be addressed in the framework of stochastic control [21,43]. In this framework, the system to be controlled is often described with a system of SDEs. With a performance index to be minimized or maximized, owing to the dynamic programming principle, finding the optimal control reduces to solving a Hamilton–Jacobi–Bellman (HJB) equation: a degenerate elliptic or parabolic equation whose functional form is problem-dependent [21]. Both mathematical and numerical approaches have been addressed for analysis of HJB equations in many problems, ranging from financial [17], economic [1], and insurance problems [14] to biological [62], environmental [12], and ecological problems [60]. Neilan et al. [40] reviewed recent progresses of numerical schemes for solving HJB and related nonlinear differential equations.

The logistic equation having the form (1) $\frac{\mathrm{d}{X}_t}{\mathrm{d}t}=r\left(1-\frac{X_t}{K}\right)X_t$(1) is one of the most widely-used mathematical models for population growth [11]. Here, $t$ is the time, $r>0$ is the intrinsic growth rate, $K>0$ is the carrying capacity, and $X_t$ is the population at the time $t$. Starting from the basic equation (1), many variants like the equation with the uncertain parameters [18], that with stochastic intrinsic growth rate [36], and that having noise-driven carrying capacity [4], have been proposed and analysed. Our problem also starts from the basic Equation (1) as explained below.

In this paper, we present a prompt report on modelling and control of a stochastic logistic equation-based model with an emphasis on algae population management in river environment [35,56]. Dense and wide spreading of the invasive macrophytes and benthic algae, such as Egeria densa [26,51] and Cladophora glomerata [64,66], has been an urgent issue. The motivation of the present mathematical model is the reports on the non-trivial dependence of the in-river population of E. densa on the flow speed [27,30]. They reported that the stable equilibrium of the algae population is a concave unimodal function of the flow speed, which cannot be reproduced with the conventional mathematical models [22,55,63,64] based on the logistic-type model (2) $\frac{\mathrm{d}{X}_t}{\mathrm{d}t}=r\left(1-\frac{X_t}{K}\right)X_t-f\left(V\right)X_t,$(2) where $V$ is the flow speed and $f$ is the decay coefficient of the algae population as an increasing function of $V$. The model (2) having the constant parameters $K$ and $r$ is able to reproduce the population decay due to high-speed river flows [35,55], but the reported unimodal nature cannot be reproduced with the model (2) unless the flow-dependent $K$ and/or $r$ are employed. In river environment, the flow discharge, or flow speed, are the variables controlled by human through operating hydraulic structures like dams and weirs [32]. Therefore, from the viewpoint of algae population management [63,64], controlling the population dynamics in a river environment can be seen as a stochastic optimization problem of a logistic-type equation having control-dependent parameters. So far, such a mathematical approach has not been examined.

The objectives of this paper are thus formulation and analysis of a stochastic optimization problem based on the logistic-type population dynamics with a control-dependent carrying capacity. As a first step of the new mathematical modelling, the focus of this paper is on analysis of a simplified problem rather than model application to practical cases. Solving the optimization problem reduces to finding a solution to an associated HJB equation, which is a non-local first-order differential equation. It is shown that the HJB equation has a unique viscosity solution: an appropriate weak solution for this kind of equations [16]. The HJB equation is not exactly-solvable as in those arising in other ecological problems [59,60], but asymptotic estimates of its solution are found. A finite difference scheme for approximation of the solution is then presented. The model is finally applied to management of benthic algae to see dependence of the optimal control on the parameters in the population dynamics and performance index.

The rest of this paper is organized as follows. Section 2 presents the mathematical model focused on in this paper and derives the HJB equation. Its mathematical analysis is carried out in this section as well. Section 3 presents the finite difference scheme for discretization of the HJB equation and its mathematical analysis results. Section 4 is devoted to numerical computation of the HJB equation focusing on the algae population management problem. Section 5 concludes this paper and presents future perspectives of our research. The appendix contains the proofs of propositions stated in Sections 2 and 3.

2. Mathematical model

The population dynamics model and performance index are presented and the HJB equation to be solved is derived and analysed. The model is firstly formulated in a generalized manner, and a specific problem related to the benthic algae management is formulated in Section 2.5.

2.1. Problem setting

A logistic-type model for describing control-dependent population dynamics in stochastic environment is presented. The time is denoted as $t\ge 0$. The population at each time $t$ is a non-negative variable denoted as $X_t$. The population dynamics is described by an SDE driven by a multiplicative jump noise modelled with the standard Poisson process $P_t$ ($P_0=0$) with the intensity parameter $\rho >0$ and the jump size $0<z<1$, where ${\rho }^{-1}$ represents the mean waiting time between each successive jumps. For the benthic algae population management, the noise corresponds to flood events such that a large part of the population is suddenly removed from riverbed. The parameter $z$ corresponds to the magnitude of the flood events. The process $P_t$ is defined on the usual complete probability space. A natural filtration generated by $P_t$ is denoted as ${\mathcal{F}}_t$ [43]. The control variable to be optimized is a measurable variable of the time $t\ge 0$, which is denoted as $q_t$ at each $t$. As in standard problem settings [43], $q_t$ is assumed to be valued in the compact set $Q=\left[q_{min},q_{max}\right]$ with $0<q_{min}<q_{max}<+\mathrm{\infty }$, and to be adapted to ${\mathcal{F}}_t$. Consequently, we define the set of admissible controls as $\mathcal{Q}$, which is a collection of all measurable continuous-time processes $q_t$ ($t\ge 0$) adapted to the filtration ${\mathcal{F}}_t$.

2.2. Population dynamics to be controlled

The SDE governing the evolution of $X_t$ is set as (3) $\mathrm{d}{X}_t=r\left(q_{t-0}\right)\left(1-\frac{X_{t-0}}{K\left(q_{t-0}\right)}\right)X_{t-0}\mathrm{d}t-f\left(q_{t-0}\right)X_{t-0}\mathrm{d}t-X_{t-0}\mathrm{d}{P}_t\mathrm{for}t>0$(3) subject to the initial condition $0\le X_0\le K_{max}\equiv \underset{q\in Q}{max}K\left(q\right)$. Here, $K=K\left(q_t\right)$ is the carrying capacity of the population, $r=r\left(q_t\right)$ is the intrinsic growth rate, and $f=f\left(q_t\right)$ is the decay coefficient. The coefficients $K$, $r$, and $f$ are assumed to be Lipschitz continuous functions in $Q$. Assume $\underset{q\in Q}{min}K\left(q\right),\underset{q\in Q}{min}r\left(q\right),\underset{q\in Q}{min}f\left(q\right)>0$ without significant loss of generality.

The SDE (3) contains the four terms. The left-hand side is the increment of the population $X_t$ in the infinitesimal time interval, the first term in the right-hand side is the logistic growth term, the second term represents the gradual population decay, and the third term represents the sudden population decay. In the context of the benthic algae population management, the second term is due to regular river flows, while the third term is due to flood events. At each jump of $P_t$, $X_t$ suddenly decreases as (4) $X_t=\left(1-z\right)X_{t-0}.$(4) By the initial condition and the sign and form of each term of (3), it is shown that the process $X_t$ satisfies $0\le X_t\le K_{max}$ for $t\ge 0$. The range of the process $X_t$ is therefore compact under the assumption and is denoted as $D=\left[0,K_{max}\right]$.

Remark 2.1:

Solutions to the SDE (3) can have a unique solution for some controls $q_t$ that are admissible since at least the constant controls such that $q_t\equiv C\in Q$ with a constant $C\in \left[q_{min},q_{max}\right]$ are admissible. This is because of the standard discussion for the uniqueness result without the jump term following Theorem 1 of Lv et al. [37] for such a control, since the noise term does not affect uniqueness of the solutions. In addition, $X_t>0$ for $t>0$ a.s. if $X_0>0$, because of the relationship (4) at each jump and the fact that the jump process is a standard Poisson process. Therefore, we have ${\tau }_x=+\mathrm{\infty }$ a.s. for $x>0$, where ${\tau }_x=\underset{t>0}{inf}\left\{t|X_t=0,X_0=x\right\}$. Strong uniqueness of solutions to the SDE (3) should be studied from a mathematically more rigorous viewpoint, but we do not address this issue since our mathematical focus is on the controlling part (HJB equation).

Remark 2.2:

The SDE (3) can be extended to models with continuous noises [49], wider class of jump terms [33], and model ambiguity [24]. These models would be practically and theoretically of more interest but may be less analytically tractable.

Remark 2.3:

The SDE (3) can also be considered in the context of optimal control of piecewise deterministic Markov processes [Chapter 8 of 9] since it is driven by a standard Poisson process.

2.3. Performance index to be minimized

The performance index $\phi =\phi \left(x;q\right)$ is a functional of the initial condition $X_0=x$ and an admissible control $q_t$ ($t\ge 0$), which is set in this paper as (5) $\phi \left(x;q\right)=\mathrm{E}\left[{\int }_0^{+\mathrm{\infty }}{e}^{-\delta t}g\left(q_t,X_t\right)\mathrm{d}t\right],$(5) where $\mathrm{E}$ is the expectation, $\delta >0$ is the discount rate, and $g$ is a non-negative and Lipschitz continuous function with respect to the first argument in $D$ and the second argument in $Q$. Therefore, $g$ is Lipschitz continuous in the compact set $D\times Q$. The function $g$ models the cost of the management and the disutility caused by the existence of the population [58,63–65] in a general manner. On the function $g$, we assume that it is increasing with respect to the second argument, meaning that the disutility caused by the population is larger for larger population. In addition, it is assumed that $\underset{q\in Q}{min}g\left(q,0\right)=0$, meaning that no management problem arises if there is no population. The discount rate $\delta$ measures perspective of the decision-maker, the manager of the population. The decision-maker has a shorter-term perspective when $\delta$ is larger as in the conventional stochastic control problems [10].

The minimized performance index $\phi$ with respect to controls $q\in \mathcal{Q}$ is the value function $\mathrm{\Phi }=\mathrm{\Phi }\left(x\right)$ in the sense of Markov control. This quantity is defined as (6) $\mathrm{\Phi }\left(x\right)=\underset{q\in \mathcal{Q}}{inf}\phi \left(x;q\right).$(6) In this paper, we restrict ourselves to consider Markovian control processes of the form $q^{\ast }=\tilde{q}\left(X_t\right)$, $t\ge 0$ with some measurable function $\tilde{q}$ in the set of admissible controls. The minimizer that archives $\mathrm{\Phi }$ is called the optimal control $q^{\ast }=q^{\ast }\left(X_t\right)$: (7) $q^{\ast }=\underset{q\in \mathcal{Q}}{\arg min}\phi .$(7) The optimal control $q^{\ast }$ is a function of the process $X_t$ and hence a feed-back control. The goal of the present problem is to find $q^{\ast }$, which is addressed in this paper through solving an HJB equation presented in the next sub-section. Clearly, $\mathrm{\Phi }$ is non-negative and bounded because of the inequality (8) $0\le \frac{1}{\delta }\underset{D\times Q}{min}g\left(q,x\right)\le \mathrm{\Phi }\left(x\right)\le \frac{1}{\delta }\underset{D\times Q}{max}g\left(q,x\right)<+\mathrm{\infty }.$(8)

2.4. HJB equation

Introduce the interval $\hat{D}=\left(0,K_{max}\right]$ for the sake of brevity of descriptions. Application of the dynamic programming principle (Chapter 3 of Øksendal and Sulem [43]) to (6) formally leads to the governing equation of $\mathrm{\Phi }$: the HJB equation (9) $\begin{array}{rl}& \delta \mathrm{\Phi }+\rho \left(\mathrm{\Phi }-\mathrm{\Phi }\left(\left(1-z\right)x\right)\right)\\ & -\underset{q\in Q}{inf}\left\{-f\left(q\right)x\frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}x}+r\left(q\right)\left(1-\frac{x}{K\left(q\right)}\right)x\frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}x}+g\left(q,x\right)\right\}=0\mathrm{in}\hat{D}\end{array}$(9) subject to the boundary condition $\mathrm{\Phi }\left(0\right)=0$, which is due to (10) $\begin{array}{rl}\mathrm{\Phi }\left(0\right)& =\underset{q\in \mathcal{Q}}{inf}\mathrm{E}\left[{\int }_0^{+\mathrm{\infty }}{e}^{-\delta t}g\left(q_t,0\right)\mathrm{d}t\right]\\ & ={\int }_0^{+\mathrm{\infty }}{e}^{-\delta t}\mathrm{d}t\left(\underset{q\in Q}{min}g\left(q,0\right)\right)\\ & =0\end{array}$(10) and the non-negativity (8). From the HJB equation (9), with an abuse of notation, the optimal control $q^{\ast }=q^{\ast }\left(x\right)$ for each $x\in \hat{D}$ is obtained as (11) $q^{\ast }\left(x\right)=\underset{q\in Q}{\arg min}\left\{-f\left(q\right)x\frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}x}+r\left(q\right)\left(1-\frac{x}{K\left(q\right)}\right)x\frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}x}+g\left(q,x\right)\right\}$(11) which can be formally determined through the derivative $\mathrm{d}\mathrm{\Phi }/\mathrm{d}x$ that exists almost everywhere inside $D$ (See, proof of Proposition 2.1). In this sense, finding $q^{\ast }$, the quantity to be evaluated through tackling the optimal control problem, can be achieved through solving the HJB equation (9).

Unfortunately, HJB equations may not have classical solutions belonging to the function space $C\left(0,K_{max}\right)\cap C^1\left(D\right)$ in general because of their nonlinearities [21]. Therefore, a weak notion of solutions, namely the concept of viscosity solutions [16] is introduced. This definition harmonizes with uniqueness of solutions to the equation and its numerical discretization. Let $USC\left(D\right)$ and $LSC\left(D\right)$ be spaces of upper- and lower- semicontiguous functions in $D$, respectively. Set $D_0=\left[0,K_{max}\right)$. Viscosity solutions to the HJB equation (9) are defined as follows [5,16,21].

Definition 2.1:

1. A function $\mathrm{\Psi }$ such that $\mathrm{\Psi }\in C\left(D_0\right)\cap USC\left(D\right)$ and $\mathrm{\Psi }\left(0\right)\le 0$ is a viscosity sub-solution to the HJB equation (9) if for all $x_0\in \hat{D}$ and for all $\varphi \in C^1\left(D\right)$, $\varphi -\mathrm{\Psi }$ is locally minimized at $x=x_0$, $\varphi \ge \mathrm{\Psi }$ on $D$, and (12) $\begin{array}{rl}& \delta \mathrm{\Psi }+\rho \left(\mathrm{\Psi }-\mathrm{\Psi }\left(\left(1-z\right)x\right)\right)\\ & -\underset{q\in Q}{inf}\left\{-f\left(q\right)x\frac{\mathrm{d}\varphi }{\mathrm{d}x}+r\left(q\right)\left(1-\frac{x}{K\left(q\right)}\right)x\frac{\mathrm{d}\varphi }{\mathrm{d}x}+g\left(q,x\right)\right\}\le 0\mathrm{at}x=x_0\end{array}$(12)

2. A function $\mathrm{\Psi }$ such that $\mathrm{\Psi }\in C\left(D_0\right)\cap LSC\left(D\right)$ and $\mathrm{\Psi }\left(0\right)\ge 0$ is a viscosity super-solution to the HJB equation (9) if for all $x_0\in \hat{D}$ and for all $\psi \in C^1\left(D\right)$, $\psi -\mathrm{\Psi }$ is locally maximized at $x=x_0$, $\psi \le \mathrm{\Psi }$ on $D$, and (13) $\begin{array}{rl}& \delta \mathrm{\Psi }+\rho \left(\mathrm{\Psi }-\mathrm{\Psi }\left(\left(1-z\right)x\right)\right)\\ & -\underset{q\in Q}{inf}\left\{-f\left(q\right)x\frac{\mathrm{d}\psi }{\mathrm{d}x}+r\left(q\right)\left(1-\frac{x}{K\left(q\right)}\right)x\frac{\mathrm{d}\psi }{\mathrm{d}x}+g\left(q,x\right)\right\}\ge 0\mathrm{at}x=x_0\end{array}$(13)

3. A bounded and continuous function $\mathrm{\Psi }:D\to \mathbf{\text{R}}$ is a viscosity solution to the HJB equation (9) if it is a viscosity sub-solution in the sense of Definition 2.1(a) as well as a viscosity super-solution in the sense of Definition 2.1(b).

Now, we prove the following propositions, which are main mathematical results of this paper. They suggest that finding the optimal control $q^{\ast }$ through solving the HJB equation (9) is a mathematically rigorous approach for the present optimization problem.

Proposition 2.1:

The value function $\mathrm{\Phi }$ is a viscosity solution to the HJB equation (9).

Proposition 2.2:

The HJB equation (9) admits at most one viscosity solution.

Proposition 2.3:

The value function $\mathrm{\Phi }$ is the unique viscosity solution to the HJB equation (9).

Proposition 2.1 shows existence of solutions to the HJB equation (9), and their uniqueness is proved in Proposition 2.2. Proposition 2.3 is a consequence of Propositions 2.1 and 2.2. Uniqueness of viscosity solutions to HJB equations is not a trivial issue even for seemingly simple problems [6,67,69].

2.5. Specific model for the algae population control

In this sub-section, the functional forms of the coefficients $K$, $r$, $f$, and $g$ are specified for more detailed analysis. A motivation of considering control-dependent coefficients $K$ and $r$ is the recently-reported semi-quadratic formula between the stable population density $\overline{x}$ of E. densa and the flow speed $V>0$ in a Japanese river [27] (14) $\overline{x}=max\left\{C_0-C_1{\left(V-\overline{V}\right)}^2,0\right\}$(14) with positive constants $C_0$, $C_1$, and $\overline{V}$, which are identified from observations at several sampling stations in the river. In the context of the equations (2) or (3) (the control $q$ is here identified as the flow speed $V$), assuming that the jumps are not so frequent such that a stable equilibrium of the population is well-established, we obtain (15) $\overline{x}=\left(1-\frac{f\left(V\right)}{r}\right)K\mathrm{if}f\left(V\right)<r.$(15) Since it is physically reasonable to consider that the decay coefficient $f$ is increasing with respect to the flow speed $V$, consistency between (14) and (15) cannot be obtained unless we assume flow dependence of $K$ and $r$. At least, it cannot be obtained if both $K$ and $r$ are constants.

In this paper, for the sake of simplicity of analysis, we consider the case $r$ is a constant while $K=K\left(V\right)$. This relationship can be justified if $K\left(V\right)$ is positive and increasing with respect to $V$, at least for not large $V$ such that E. densa are critically removed from the riverbed [64]. The increasing nature of $K\left(V\right)$ can be considered because of more efficient nutrient supply by larger nutrient flux owing to the faster flow speed $V$. In addition, without significant loss of generality, it is assumed that $f$ is smooth and strictly increasing, so that it has the unique inverse $f^{\left(-1\right)}$.

Under the stated assumptions, equating (14) and (15) formally gives (16) $K\left(V\right)=\frac{r}{r-f\left(V\right)}max\left\{C_0-C_1{\left(V-\overline{V}\right)}^2,0\right\}.$(16) Assume $C_0>C_1{\overline{V}}^2$ so that $K\left(0\right)>0$. This condition is in fact satisfied in the data [27]. Then, (16) reduces to (17) $K\left(V\right)=\frac{r}{r-f\left(V\right)}\left(C_0-C_1{\left(V-\overline{V}\right)}^2\right)$(17) at least for $V<min\left\{f^{\left(-1\right)}\left(r\right),\overline{V}+\sqrt{C_0{C}_1^{-1}}\right\}$. Positivity of the stable equilibrium $\overline{x}$ cannot be achieved for large $V$ such that $r\le f\left(V\right)$ where the equilibrium $X_t=0$ is stable. This thresholding $V$ is identified with the larger zero of the right-hand side of (17), which is $\overline{V}+\sqrt{C_0{C}_1^{-1}}$. Therefore, we assume that $f^{\left(-1\right)}\left(r\right)$ equals $\overline{V}+\sqrt{C_0{C}_1^{-1}}$, converting (17) to (18) $K\left(V\right)=\frac{C_{1}r}{r-f\left(V\right)}\left(\left(f^{\left(-1\right)}\left(r\right)-V\right)\left(V-2\overline{V}+f^{\left(-1\right)}\left(r\right)\right)\right).$(18) We consider here the most tractable model, which is the simplest form of the decay (19) $f\left(V\right)=\alpha V$(19) with a constant $\alpha >0$. Substituting (19) into (18) yields the linear result (20) $K\left(V\right)=aV+b,$(20) where (21) $a=\frac{C_{1}r}{\alpha }>0\mathrm{and}b=\frac{C_{1}r}{\alpha }\left(\frac{r}{\alpha }-2\overline{V}\right)>0.$(21) The formula (20) is in accordance with the assumption that $K$ is increasing with respect to $V$. Assuming (20) instead of the original one (16) can avoid the problem that the carrying capacity, which should be a positive value from a biological viewpoint, vanishes. The linear formula (20) is employed in the rest of this paper. Formally setting $a=0$ in (21) recovers the conventional constant carrying capacity $K=b$.

In this paper, the control $q$ appearing in the stochastic control problem is identified with the flow speed $V$. For the coefficient $g$ in the performance index $\phi$, following the models on management of harmful species population [63–66], set (22) $g\left(q,X\right)=\frac{w}{2}(q-\overline{q}{)}^2+X^m$(22) with a weighting constant $w>0$, a target value $\overline{q}$ such that $q_{min}<\overline{q}<q_{max}$, and a constant $m>1$. The first term of (22) measures the deviation of the control $q$ from the target value $\overline{q}$, while the second term represents the disutility caused by the population. This leads to the simplest form of the performance index $\phi$. For the algae population management in a river, its environmental condition can be controlled through hydraulic structures like dams and weirs, in which the control variable is the river discharge. From this viewpoint, it is not reasonable to consider the flow speed along the river; however, under moderate hydraulic conditions where the river flows are not trans-critical, the widely-used empirical formula like the Manning’s formula [61] can be utilized to estimate a one-to-one relationship between the flow speed and river discharge. In this sense, the flow speed and river discharge are interchangeable variables, at least from the view point of a conceptual mathematical modelling like that employed in this paper.

In summary, the HJB equation with the specified coefficients is (23) $\begin{array}{rl}& \delta \mathrm{\Phi }-rx\frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}x}-x^m+\rho \left(\mathrm{\Phi }-\mathrm{\Phi }\left(\left(1-z\right)x\right)\right)\\ & -\underset{q\in Q}{inf}\left\{-\alpha qx\frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}x}-\frac{r{x}^2}{aq+b}\frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}x}+\frac{w}{2}{\left(q-\overline{q}\right)}^2\right\}=0\mathrm{in}\hat{D}=\left(0,K_{max}\right]\end{array}$(23) with $\mathrm{\Phi }\left(0\right)=0$ and $K_{max}=a{q}_{max}+b$. Notice that the assumptions on the sign and Lipschitz continuity of the coefficients are not violated in the resulting equation (23).

We briefly show that minimizers of the last term of (23) exist at each $x$ if $\mathrm{d}\mathrm{\Phi }/\mathrm{d}x$ exists in the classical sense. Set $\mathrm{\Lambda }\left(q\right)=-\alpha qx\left(\mathrm{d}\mathrm{\Phi }/\mathrm{d}x\right)-(r{x}^2/aq+b)\left(\mathrm{d}\mathrm{\Phi }/\mathrm{d}x\right)+\left(w/2\right)(q-\overline{q}{)}^2$, $q\in Q$. If $\mathrm{d}\mathrm{\Phi }/\mathrm{d}x$ exists in the classical sense, its minimum is attained in $Q$. In this case, $\mathrm{\Lambda }$ has no plateau. Notice that uniqueness of viscosity solutions (Proposition 2.2) does not necessarily lead to uniqueness of optimal controls. Therefore, the minimizers of $\mathrm{\Lambda }$ may not be unique, but there exist at most a countable number of minimizers by the functional shape of $\mathrm{\Lambda }$. If $\mathrm{d}\mathrm{\Phi }/\mathrm{d}x$ does not exist in the classical sense at a point $x$, then it would be multi-valued at this point, and the above statement may not hold true. Numerically, we infer that $\mathrm{\Phi }$ is smooth and strictly increasing.

Remark 2.4:

The reported quadratic-like relationship motivated us to parameterize the carrying capacity as a positive function of the flow speed. For a more realistic mathematical modelling, the intrinsic growth rate should also depend on the flow speed. However, identifying such a relationship would require larger amount of data, especially the data of the growth curves under different flow speed. This topic is beyond the scope of this paper, and will be addressed elsewhere through hydraulic experiments.

Although the resulting HJB equation (23) still cannot be solved analytically, its solution $\mathrm{\Phi }$ and the associated optimal control $q^{\ast }$ are asymptotically estimated for small $x$. The mathematical technique utilized here is an extension of that in Proposition 4.6 of Yoshioka and Yaegashi [63], but a more sophisticated estimation of $q^{\ast }$ is achieved.

Proposition 2.4:

Assume $\delta$ is sufficiently large such that (24) $T\equiv \delta +\rho \left(1-{\left(1-z\right)}^m\right)-mr+m\alpha \overline{q}>0$(24) and the value function $\mathrm{\Phi }$ and the optimal control $q^{\ast }$ are analytic for small $x$. Then, we have the asymptotic estimate of $\mathrm{\Phi }$ and $q^{\ast }$ as (25) $\mathrm{\Phi }\left(x\right)\sim c_0{x}^m$(25) and (26) $q^{\ast }\left(x\right)-\overline{q}\sim c_1{x}^m-c_2\left(1-c_3{x}^m\right)x^{m+1}.$(26) respectively, where (27) $c_0=T^{-1},c_1=\frac{m\alpha }{w}{c}_0,c_2=\frac{r{c}_0}{w{a}^2}{\left(\overline{q}+\frac{b}{a}\right)}^{-2},c_3=2{\left(\overline{q}+\frac{b}{a}\right)}^{-1}{c}_1{c}_2.$(27)

The asymptotic estimates in Proposition 2.4 are compared with numerical solutions later. The techniques utilized here are inspired from that in [63]. However, the obtained estimate of $q^{\ast }$ is not a monomial as in the literature but a polynomial. In fact, the monomial function having a simpler form (28) $q^{\ast }\left(x\right)-\overline{q}\sim c_1{x}^m$(28) potentially serves as an asymptotic estimate as well, but it turns out to be less accurate than (26). In this sense, the second term of (26) as a higher-order polynomial is a correction term for a more accurate estimate.

The estimates in Proposition 2.4 are useful for understanding how $\mathrm{\Phi }$ and $q^{\ast }$ are effectively scaled with the parameters and variables. One of their drawbacks is that they are effective only for sufficiently large $\delta$ such that (24) is satisfied, despite the HJB equation (23) is well-posed for $\delta >0$. Assume that the jump term is negligible in the population dynamics for the sake of simplicity. As identified in Yoshioka and Yaegashi [64], the intrinsic growth rate $r$ is $O\left({10}^{-1}\right)$ (1/day) to $O\left({10}^0\right)$ (1/day) and $m$ would be $O\left({10}^0\right)$. In addition, the decay term $\alpha \overline{q}$ has the same magnitude but is often smaller than $r$. Therefore, the present asymptotic estimates are valid for the discount rate $\delta$ not smaller than $O\left({10}^{-1}\right)$ (1/day). This implies that these asymptotic estimates can be used for problems with the decision-maker who tries to manage the river environment in daily or weekly time-scales. When the jump effect is considered, the restriction on $\delta$ can be somewhat relaxed since $\rho$ is $O\left({10}^{-1}\right)$ (1/day) [64].

Despite the results in Proposition 2.4 are tractable, they are valid only for small initial population $x$. Nevertheless, they are consistent with the parameter dependence result (Proposition 2.5). The proof of the proposition is omitted here since it is essentially the same with that in Proposition 3.4 of Yoshioka and Yaegashi [63].

Proposition 2.5:

The value function $\mathrm{\Phi }$ is increasing with respect to $-\delta$, $w$ and $x$.

The proposition gives an intuitive result that the value function $\mathrm{\Phi }$, which is an optimized cost, is larger if the deviation of the control $q$ from the target value $\overline{q}$ is evaluated larger and if the initial population is larger. The result on $-\delta$ implies that a larger discount means accumulation of the costs in a longer time-scale.

Based on the asymptotic estimates in Proposition 2.4, we can obtain more detailed parameter dependence of $\mathrm{\Phi }$ and $q^{\ast }$ for the decision-maker who tries to manage the river environment in a short time-scale. The proof of the Proposition 2.6 is not presented here since it is just partial differentiations of $\mathrm{\Phi }$ and $q^{\ast }$ with respect to the parameters and variables.

Proposition 2.6:

Assume (24), then the following (a) and (b) hold true.

1. For $x>0$, the asymptotic estimate (25) of the value function $\mathrm{\Phi }$ satisfies (29) $\frac{\mathrm{\partial \Phi }}{\mathrm{\partial }x}>0,\frac{\mathrm{\partial \Phi }}{\mathrm{\partial }\delta }<0,\frac{\mathrm{\partial \Phi }}{\mathrm{\partial }r}>0,\frac{\mathrm{\partial \Phi }}{\mathrm{\partial }\overline{q}}<0,\frac{\mathrm{\partial \Phi }}{\mathrm{\partial }\alpha }<0,\frac{\mathrm{\partial \Phi }}{\mathrm{\partial }\rho }<0,\mathrm{and}\frac{\mathrm{\partial \Phi }}{\mathrm{\partial }z}<0.$(29)

2. For small $x>0$, the asymptotic estimate (26) of the optimal control $q^{\ast }$ satisfies (30) $\frac{\mathrm{\partial }{q}^{\ast }}{\mathrm{\partial }x}>0,\frac{\mathrm{\partial }{q}^{\ast }}{\mathrm{\partial }\delta }<0,\frac{\mathrm{\partial }{q}^{\ast }}{\mathrm{\partial }r}>0,\frac{\mathrm{\partial }{q}^{\ast }}{\mathrm{\partial }\overline{q}}>0,\frac{\mathrm{\partial }{q}^{\ast }}{\mathrm{\partial }\alpha }<0,\frac{\mathrm{\partial }{q}^{\ast }}{\mathrm{\partial }\rho }<0,\mathrm{and}\frac{\mathrm{\partial }{q}^{\ast }}{\mathrm{\partial }z}<0$(30) and (31) $\frac{\mathrm{\partial }{q}^{\ast }}{\mathrm{\partial }w}>0,\frac{\mathrm{\partial }{q}^{\ast }}{\mathrm{\partial }a}>0\mathrm{and}\frac{\mathrm{\partial }{q}^{\ast }}{\mathrm{\partial }b}>0.$(31)

Proposition 2.6, despite its results are simple inequalities, provide several insights into management of river environment from a qualitative standpoint. Quantitative analysis on the parameter dependence is carried out in Section 4. Proposition 2.6(a) clearly gives the parameter dependence of the value function $\mathrm{\Phi }$. The first results (30) of Proposition 2.6(b) inherit (29) except for that on $\overline{q}$. These inequalities imply that the value function $\mathrm{\Phi }$ is smaller for the environment in which the population growth can be effectively supressed. From a management view point of river environment, increasing the magnitude and frequency of flood events, if their costs are small, may be an effective way for more cost-effective algae population management. This would especially apply to the river environment in a dam downstream where the magnitude and frequency of floods are artificially decreased from natural conditions where $\rho$ is very small [34,41]. In addition, if such an environment has been established, the population can be managed subject to small deviation $q^{\ast }-\overline{q}$. The second results (31) of Proposition 2.6(b) are on the dependence of the optimal control $q^{\ast }$, or equivalently the deviation $q^{\ast }-\overline{q}$, on the parameters not appearing in (29) and (31). The result on the weighting coefficient $w$ simply means that the deviation is smaller for the decision-maker who is more concerned with minimizing the cost due to the deviation of the control. The result on the parameters $a$ and $b$ of the carrying capacity $K$ indicates that management of the algae population with higher carrying capacity, which can be identified with the maximum population abundance, should tolerate larger magnitude of deviation $q^{\ast }-\overline{q}$. The obtained result implies that management of the invasive species with potentially widespread and higher population abundance [20,47,70], like E. densa and C. glomerata, at the early growth stage may be subject to a larger deviation of environmental condition from the targeted one.

3. Numerical scheme

The finite difference scheme for discretization of the HJB equation (23) is formulated and its solution algorithm is presented, since the equation cannot be solved analytically and the asymptotic estimates in Proposition 2.4 do not apply globally in the domain $D$. The scheme is based on the exponentially-fitted formalism for nonlinear advection-decay equations [64]. This discretization achieves stable numerical computation since it complies with the monotonicity [5], which serves as a powerful condition to guarantee uniqueness and stability of numerical solutions. The exponentially-fitted discretization is motivated by the fact that linear advection-decay equations often have solutions that are exponential or exponential-like functions. As a result, the scheme can generate nodally-exact solutions to advection-decay equations having constant coefficients. Similar formalism has been applied to nonlinear advection–diffusion equations related to stochastic optimal control problems [23,57]. In this paper, solving the resulting discretized system is carried out with a standard policy iteration algorithm [2], which has effectively been applied to numerical resolution of discretized HJB equations [15,45].

3.1. Finite difference discretization

The discretization basically follows the original one with an upwind treatment for the terms involving the derivative $\mathrm{d}\mathrm{\Phi }/\mathrm{d}x$. The domain $D$ is uniformly discretized with $I\ge 2$ vertices $x_i=i{I}^{-1}{K}_{max}$ $(1\le i\le I)$. Each quantity $\vartheta$ discretized at the vertex $x_i$ is denoted as ${\vartheta }_i$. The solution $\mathrm{\Phi }$ is approximated at each $x_i$. Let $l=I^{-1}{K}_{max}$. The discretization procedure below focuses on a vertex $x_i$ $(1\le i\le I)$ without the loss of generality since the boundary condition ${\mathrm{\Phi }}_0=0$ is directly specified.

Discretization of each term of the HJB equation (23) is carried out as follows: (32) $\delta \mathrm{\Phi }-rx\frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}x}-x^m\to \delta {\mathrm{\Phi }}_i-\frac{P_i}{e^{P_i}-1}r{x}_i\frac{{\mathrm{\Phi }}_{i+1}-{\mathrm{\Phi }}_i}{l}-x_i^m\mathrm{for}0<i<I$(32) and (33) $\delta \mathrm{\Phi }-rx\frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}x}-x^m\to \delta {\mathrm{\Phi }}_I-r{x}_I\frac{{\mathrm{\Phi }}_I-{\mathrm{\Phi }}_{I-1}}{l}-x_I^m$(33) with $P_i=\delta l{r}^{-1}{x}_i^{-1}>0$, (34) $\begin{array}{rl}& \underset{q\in Q}{inf}\left\{-\alpha qx\frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}x}-\frac{r{x}^2}{aq+b}\frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}x}+\frac{w}{2}{\left(q-\overline{q}\right)}^2\right\}\\ & \to -\alpha q_i^{\ast }{x}_i\frac{{\mathrm{\Phi }}_i-{\mathrm{\Phi }}_{i-1}}{l}-\frac{r{x}_i^2}{a{q}_i^{\ast }+b}\frac{{\mathrm{\Phi }}_i-{\mathrm{\Phi }}_{i-1}}{l}+\frac{w}{2}(q_i^{\ast }-\overline{q}{)}^2\end{array}$(34) with (35) $q_i^{\ast }=\underset{q\in Q}{\arg min}\left\{-\alpha q{x}_i\frac{{\mathrm{\Phi }}_i-{\mathrm{\Phi }}_{i-1}}{l}-\frac{r{x}_i^2}{aq+b}\frac{{\mathrm{\Phi }}_i-{\mathrm{\Phi }}_{i-1}}{l}+\frac{w}{2}{\left(q-\overline{q}\right)}^2\right\},$(35) and (36) $\rho \left(\mathrm{\Phi }-\mathrm{\Phi }\left(\left(1-z\right)x\right)\right)\to \rho \left\{{\mathrm{\Phi }}_i-\left({\omega }_i{\mathrm{\Phi }}_{j_i}+\left(1-{\omega }_i\right){\mathrm{\Phi }}_{j_i+1}\right)\right\}$(36) with (37) ${\omega }_i=\frac{x_{j_i+1}-\left(1-z\right)x_i}{x_{j_i+1}-x_{j_i}}.$(37) Here, the index $0\le j_i\le i-1$ in (36) is uniquely determined as the integer with $x_{j_i}\le \left(1-z\right)x_i<x_{j_i+1}$. Combining (32) through (36) yields the discretized (23) at $x_i$ for $0<i<I$ as (38) $\begin{array}{rl}\delta {\mathrm{\Phi }}_i& -\frac{P_i}{e^{P_i}-1}r{x}_i\frac{{\mathrm{\Phi }}_{i+1}-{\mathrm{\Phi }}_i}{l}+\rho \left\{{\mathrm{\Phi }}_i-\left({\omega }_i{\mathrm{\Phi }}_{j_i}+\left(1-{\omega }_i\right){\mathrm{\Phi }}_{j_i+1}\right)\right\}\\ & +\alpha q_i^{\ast }{x}_i\frac{{\mathrm{\Phi }}_i-{\mathrm{\Phi }}_{i-1}}{l}+\frac{r{x}_i^2}{a{q}_i^{\ast }+b}\frac{{\mathrm{\Phi }}_i-{\mathrm{\Phi }}_{i-1}}{l}-\frac{w}{2}(q_i^{\ast }-\overline{q}{)}^2-x_i^m=0,\end{array}$(38) which can be rewritten as (39) $\begin{array}{rl}& \left(\delta +\frac{P_i}{e^{P_i}-1}\frac{r{x}_i}{l}+\frac{\alpha q_i^{\ast }{x}_i}{l}+\frac{r{x}_i^2}{\left(a{q}_i^{\ast }+b\right)l}+\rho \right){\mathrm{\Phi }}_i\\ & -\left(\frac{\alpha q_i^{\ast }{x}_i}{l}+\frac{r{x}_i^2}{\left(a{q}_i^{\ast }+b\right)l}\right){\mathrm{\Phi }}_{i-1}-\frac{P_i}{e^{P_i}-1}\frac{r{x}_i}{l}{\mathrm{\Phi }}_{i+1}\\ & =\rho \left({\omega }_i{\mathrm{\Phi }}_{j_i}+\left(1-{\omega }_i\right){\mathrm{\Phi }}_{j_i+1}\right)+\frac{w}{2}(q_i^{\ast }-\overline{q}{)}^2+x_i^m.\end{array}$(39) Similarly, the discretized equation for $i=I$ is obtained as (40) $\begin{array}{rl}& \left(\delta +\frac{\alpha q_I^{\ast }{x}_I}{l}+\frac{r{x}_I^2}{\left(a{q}_I^{\ast }+b\right)l}-\frac{r{x}_I}{l}+\rho \right){\mathrm{\Phi }}_I-\left(\frac{\alpha q_I^{\ast }{x}_I}{l}+\frac{r{x}_I^2}{\left(a{q}_I^{\ast }+b\right)l}-\frac{r{x}_I}{l}\right){\mathrm{\Phi }}_{I-1}\\ & =\rho \left({\omega }_I{\mathrm{\Phi }}_{j_I}+\left(1-{\omega }_I\right){\mathrm{\Phi }}_{j_I+1}\right)+\frac{w}{2}(q_I^{\ast }-\overline{q}{)}^2+x_I^m.\end{array}$(40) Using the one-sided discretization (33) is because of the fact that the present exponentially-fitted methods cannot be applied to the discretization at $i=I$, since it has been originally developed for two-point boundary value problems. The formula (35) gives the discretized optimal control at each vertex.

Assembling all (39) and (40) supplemented with the boundary condition ${\mathrm{\Phi }}_0=0$ yields the matrix system governing the numerical solution vector $\overrightarrow{\mathrm{\Phi }}=[{\mathrm{\Phi }}_i{]}_{0\le i\le I}$: (41) $A\overrightarrow{\mathrm{\Phi }}=B\overrightarrow{\mathrm{\Phi }}+\overrightarrow{S},$(41) where the matrices $A=A\left(\overrightarrow{\mathrm{\Phi }}\right)=[A_{i,j}{]}_{0\le i,j\le I}$ and $B=[B_{i,j}{]}_{0\le i,j\le I}$ are given as follows: (42) $\begin{array}{r}{A}_{i,j}=\left\{\begin{array}{ll}{\delta }_{i,0}& \left(i=0,0\le j\le I\right)\\ -{\mu }_{i,i-1}& \left(1\le i\le I,j=i-1\right)\\ {\mu }_{i,i}& \left(1\le i\le I,j=i\right)\\ -{\mu }_{i,i+1}& \left(1\le i\le I-1,j=i+1\right)\\ 0& \left(\mathrm{Otherwise}\right)\end{array}\right.,B_{i,j}=\left\{\begin{array}{ll}0& \left(i=0\right)\\ \rho {\omega }_i& \left(1\le i\le I,j=j_i\right)\\ \rho \left(1-{\omega }_i\right)& \left(1\le i\le I,j=j_i+1\right)\\ 0& \left(\mathrm{Otherwise}\right)\end{array}\right.\end{array}$(42) with each $\mu$ uniquely determined from (39) and (40). We have (43) ${\mu }_{i,i-1}=\frac{\alpha q_i^{\ast }{x}_i}{l}+\frac{r{x}_i^2}{\left(a{q}_i^{\ast }+b\right)l}>0\mathrm{for}1\le i\le I-1$(43) and (44) ${\mu }_{I,I-1}=\frac{\alpha q_I^{\ast }{x}_I}{l}+\frac{r{x}_I^2}{\left(a{q}_I^{\ast }+b\right)l}-\frac{r{x}_I}{l}\ge \frac{\alpha q_I^{\ast }{x}_I}{l}+\frac{r{x}_I^2}{x_{I}l}-\frac{r{x}_I}{l}=\frac{\alpha q_I^{\ast }{x}_I}{l}>0$(44) since $a{q}_I^{\ast }+b\le K_{max}=x_I$. Similarly, we have (45) ${\mu }_{i,i}>0\mathrm{and}{\mu }_{i,i+1}>0.$(45) The source term vector $\overrightarrow{S}=[S_i{]}_{0\le i\le I}$ is given as (46) $S_i=\left\{\begin{array}{ll}0& \left(i=0\right),\\ \frac{w}{2}{\left(q_i^{\ast }-\overline{q}\right)}^2+x_i^m\left(\ge 0\right)& \left(1\le i\le I\right).\end{array}\right.$(46) Now, the HJB equation (23) has been fully discretized. Since (47) ${\mu }_{i,i}>{\mu }_{i,i-1}+{\mu }_{i,i+1}(0<i<I)\text{and}{\mu }_{I,I}>{\mu }_{I,I-1},$(47) the matrix $A$ has positive diagonal elements and negative non-diagonal elements, and is strictly diagonally-dominant. It is therefore a weakly diagonally-dominant Z-matrix [44] having non-negative diagonal elements. In addition, it is non-singular (invertible) due to Lemma 6.5 of Santambrogio [46] by ${\mu }_{i,i-1}>0$ $(0<i\le I)$. Therefore, it is an M-matrix (Theorem 3.2.5(i)–(ii) of Azimzadeh [5]) whose inverse is a non-negative matrix. The same applies to the matrix $A-B$. We can rewrite (41) as (48) $\overrightarrow{\mathrm{\Phi }}=A^{-1}\left(B\overrightarrow{\mathrm{\Phi }}+\overrightarrow{S}\right),$(48) which is a key equation for the policy iteration algorithm.

The next proposition shows that the present scheme can generate convergent numerical solutions.

Proposition 3.1:

The present finite difference discretization generates a unique numerical solution. In addition, numerical solutions with the finite difference discretization locally uniformly converge to the unique viscosity solution to the HJB equation (23) as $l\to +0$.

The discretized system (48) is solved with a policy iteration algorithm [2]. Set a non-negative initial guess ${\overrightarrow{\mathrm{\Phi }}}_0=[{\mathrm{\Phi }}_{0,i}{]}_{0\le i\le I}$ with ${\mathrm{\Phi }}_{0,i}\ge 0$ and the allowable error $\eta$, which is a sufficiently small positive constant. Starting from the initial guess, the algorithm runs as follows.

(Policy iteration algorithm)

• Step 1: Set $k=0$.

• Step 2: Set $\overrightarrow{q_k^{\ast }}=[q_{k,i}^{\ast }{]}_{0\le i\le I}$ via (49) $q_{k,i}^{\ast }=\underset{q\in Q}{\arg min}\left\{-\alpha q{x}_i\frac{{\mathrm{\Phi }}_{k,i}-{\mathrm{\Phi }}_{k,i-1}}{l}-\frac{r{x}_i^2}{aq+b}\frac{{\mathrm{\Phi }}_{k,i}-{\mathrm{\Phi }}_{k,i-1}}{l}+\frac{w}{2}{\left(q-\overline{q}\right)}^2\right\}.$(49)

• Step 3: Find ${\overrightarrow{\mathrm{\Phi }}}_{k+1}=[{\mathrm{\Phi }}_{k+1,i}{]}_{0\le i\le I}$ such that ${\overrightarrow{\mathrm{\Phi }}}_{k+1}=A_k^{-1}\left(B_k{\overrightarrow{\mathrm{\Phi }}}_k+{\overrightarrow{S}}_k\right)$, where $A_k$, $B_k$, and ${\overrightarrow{S}}_k$ are $A$, $B$, and $\overrightarrow{S}$ with the substitution $q_i^{\ast }=q_{k,i}^{\ast }$.

• Step 4: If $\underset{1\le i\le I}{max}|{\mathrm{\Phi }}_{k,i}-{\mathrm{\Phi }}_{k+1,i}|<\eta$, choose ${\overrightarrow{\mathrm{\Phi }}}_{k+1}$ as the approximate numerical solution and terminate the algorithm. If $\underset{1\le i\le I}{max}|{\mathrm{\Phi }}_{k,i}-{\mathrm{\Phi }}_{k+1,i}|\ge \eta$, then set $k\to k+1$ and go to Step 2.

Notice that, owing to a simple induction argument, Step 3 preserves non-negativity of ${\overrightarrow{\mathrm{\Phi }}}_{k+1}$ because of the non-negativity of $A_k^{-1}\left(B_k{\overrightarrow{\mathrm{\Phi }}}_k+{\overrightarrow{S}}_k\right)$ and ${\overrightarrow{\mathrm{\Phi }}}_0$. Inverting the matrix $A_k$, which is tri-diagonal, is carried out with the Thomas method [53].

4. Model analysis

4.1. Computational conditions

The HJB equation (23) is solved numerically with the finite difference scheme. For the sake of brevity of the analysis, we employ the following non-dimensionalization to have $K_{max}=1$, which is equivalent to considering the transformations $X_t\to (a{q}_{max}+b{)}^{-1}{X}_t$ and $aq+b\to (a{q}_{max}+b{)}^{-1}\left(aq+b\right)$ so that the domain $D$ of the HJB equation reduces to the unit interval $\left[0,1\right]$. In addition, the following transformations are carried out: $\mathrm{\Phi }\to r\mathrm{\Phi }$, $\delta \to r^{-1}\delta$, $\rho \to r^{-1}\rho$, $q\to {\overline{q}}^{-1}q$, $q_{min}\to {\overline{q}}^{-1}{q}_{min}$, $q_{max}\to {\overline{q}}^{-1}{q}_{max}$, $\alpha \to r^{-1}\alpha \overline{q}$, $w\to w{\overline{q}}^2$, so that (23) simplifies to (50) $\begin{array}{rl}& \delta \mathrm{\Phi }-x\frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}x}-x^m+\rho \left(\mathrm{\Phi }-\mathrm{\Phi }\left(\left(1-z\right)x\right)\right)\\ & -\underset{q\in Q}{inf}\left\{-\alpha qx\frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}x}-\frac{x^2}{aq+b}\frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}x}+\frac{w}{2}{\left(q-1\right)}^2\right\}=0\mathrm{in}\left(0,1\right]\end{array}$(50) which is (23) with formally setting $r=1$ and $\overline{q}=1$.

In our numerical computation, the domain $D$ is uniformly discretized with 2,001 vertices $(I=2,000)$, which has preliminary been checked to be a sufficiently high computational resolution for the analysis in this paper. $\eta$ for the policy iteration algorithm is set as ${10}^{-14}$. Table 1 summarizes the parameter values chosen for the numerical computation. Assuming $r=0.65$ (1/day) in the dimensional form [64], the specified parameter values potentially correspond to the algae population dynamics in dam downstream in a Japanese river. These values are used unless otherwise specified. Convergence of the policy iteration algorithm under this computational condition is fast, satisfying the convergence criterion at most 20 iterations for the numerical solutions plotted below.

Table 1. Parameter values used in the numerical computation.

Parameter	Value
$a$	0.3
$b$	0.4
$\alpha$	0.5
$\rho$	0.1
$z$	0.9
$\delta$	2.0
$w$	0.1
$m$	2.0
$q_{max}$	2.0
$q_{min}$	0.1

4.2. Comparison of numerical and asymptotic results

This sub-section compares numerical and asymptotic profiles of the value function $\mathrm{\Phi }$ and the optimal control $q^{\ast }$, focusing on the influence of the weighting factor $w$ as a key parameter that determines the resulting $\mathrm{\Phi }$ and $q^{\ast }$. Figure 1 compares the numerical and asymptotic $\mathrm{\Phi }$ for $w=1$, 0.1, and 0.01. Notice that the asymptotic estimate (25) does not depend on $w$. Figures 2–4 compare the numerical and asymptotic $q^{\ast }$ for $w=1$, 0.1, and 0.01, respectively. For $q^{\ast }$, both the monomial (28) and polynomial estimates (26) are examined.

Figure 1. Comparison of the numerical $(w=1,0.1,0.01)$ and asymptotic $\mathrm{\Phi }$. Black (framed filled circle): asymptotic estimate (25), Red (filled circle): numerical solution with $w=1$, Green (circle): numerical solution with $w=0.1$, and Blue (double circle): numerical solution with $w=0.01$.

Figure 2. Comparison of the numerical and asymptotic $q^{\ast }$ $(w=1)$. Here, $dq=q^{\ast }-\overline{q}$. Black (filled circle): numerical solution, Red (double circle): asymptotic estimate (28), and Blu (circle): asymptotic estimate (26).

Figure 3. Comparison of the numerical and asymptotic $q^{\ast }$ $(w=0.1)$. Here, $dq=q^{\ast }-\overline{q}$. Black (filled circle): numerical solution, Red (double circle): asymptotic estimate (28), and Blu (circle): asymptotic estimate (26).

Figure 4. Comparison of the numerical and asymptotic $q^{\ast }$ $(w=0.01)$. Here, $dq=q^{\ast }-\overline{q}$. Black (filled circle): numerical solution, Red (double circle): asymptotic estimate (28), and Blu (circle): asymptotic estimate (26).

Figure 2 for the value function $\mathrm{\Phi }$ demonstrates that the asymptotic estimate is reasonably accurate for sufficiently small $x$ and captures the convex nature of the numerical solutions. As stated above, the asymptotic estimate (25) does not depend on $w$. Figure 1 demonstrates that this property is reasonably satisfied for the numerical solutions with the different values of $w$, especially for sufficiently small $x$, indicating a qualitative consistency between the asymptotic and numerical results. For the optimal control $q^{\ast }$, the computational results in Figures 2–4 suggest that the polynomial estimate (28) is more close to the numerical results for all $w$ when $x$ is small, implying that the second term of (28) truly works as a correction term to the monomial estimate (26). As visually seen in Figures 2–4, relative error between the numerical and asymptotic estimates are smaller for larger weighting coefficient $w$, suggesting that the asymptotic estimates are more accurate for the problem with the decision-maker having more concern with minimization of the deviation of the control from the target. Overall, the computational results demonstrate that the asymptotic estimates for small $x$ in Proposition 2.4 are qualitatively correct for small $x$.

4.3. Parameter dependence of the optimal control

The HJB equation (23) is numerically solved for a variety of parameter values, so that behaviour of the optimal control $q^{\ast }$, the quantity to be found through solving the optimization problem, is explored more in detail. Figures 5–9 show the computed $q^{\ast }$ for different values of $a$ in the carrying capacity, $\alpha$ in the decay term, $\rho$ in the jump term, and $w$ and $\delta$ in the performance index.

Figure 5. Comparison of $q^{\ast }$ for $a=0.45\times l/20$ $(l=0,1,2,\dots ,20)$. The relationship $a{q}_{max}+b=1$ is maintained in all the cases plotted in this figure.

Figure 6. Comparison of $q^{\ast }$ for $\alpha =0.5\times l/20$ $(l=1,2,\dots ,20)$.

Figure 7. Comparison of $q^{\ast }$ for ${\rho }^{-1}=5.0+10.0\times l/20$ $(l=0,1,2,\dots ,20)$.

Figure 8. Comparison of $q^{\ast }$ for $w=1.0\times l/20$ $(l=1,2,3,\dots ,20)$.

Figure 9. Comparison of $q^{\ast }$ for $\delta =5.0\times l/20$ $(l=1,2,3,\dots ,20)$.

Figure 5 for $a$ indicates qualitatively different profiles of $q^{\ast }$ between small $a$ and large $a$. For the smallest $a=0$, $q^{\ast }$ is continuous and increasing with respect to the population $x$ and it saturates at a point in $\left(0,1\right)$ and equals $q^{\ast }=q_{max}$ for larger $x$. The same applies to relatively small $a\left(<0.20\right)$. However, for larger $a$, there exists a subdomain near $x=1$ such that $q^{\ast }=q_{min}$ is optimal and $q^{\ast }$ is discontinuous at the left boundary of this sub-domain. This type of discontinuous optimal control has not been found in similar stochastic control models with constant $K$ [63,65]. This sudden decrease of $q^{\ast }$ is considered to be due to the functional form $K=aq+b$ of the carrying capacity that depends on the control $q$. On the one hand, the population decays with the decay rate proportional to the control $q$. On the other hand, higher population abundance is possibly achieved for larger value of $q$ if there is no decay, which is more significant for larger $a$. The computational results suggest that the control should be minimum, even if it deviates from the target value $\overline{q}$, if the population $x$ is large and its carrying capacity is a highly increasing function of the control $q$. From a view point of management of river environment focusing on the harmful benthic algae, the computational result for large $a$ suggest that the flow speed should be close to the target value when the population $x$ is small, should be maximum if $x$ is moderate, and should be minimum if $x$ is large. For the last case, the flow speed should be set smallest as possible in the admissible range $Q$. This type of the control policy may not be feasible in practice for river environment where too low flow speed (or low flow discharge) triggers secondary issues like contractions of floodplains serving as habitats of many aquatic species [24,57]. Nevertheless, the policy can become feasible if the population is maintained small through voluntary human interventions like cleaning up of the river [62]. The discontinuities appearing in $q^{\ast }$ suggest that the model parameter should be carefully calibrated especially when the population is not small. The computational results in Figures 6–9 provide insights into the impacts of the other model parameters on the optimal control $q^{\ast }$.

Figures 6 and 7 for $\alpha$ and $\rho$ show that the magnitude of the deviation $q^{\ast }-\overline{q}$ can be made smaller if larger values of the parameters are achieved at least for not large $x$. For the benthic algae management in rivers, the parameter $\alpha$ can be made larger through increasing sediment size and sediment content in river flows [35,54], which can be possible through sediment injection into upstream of the river reach where the benthic algae are blooming [66]. Figure 7 implies that increasing the parameter $\rho$, which is an inverse of the non-dimensional mean waiting time between floods, does not critically affect $q^{\ast }$. This is considered due to the fact that the jump discontinuously decreases the population as (4), but its influence is not persistent since it is not time-continuous.

Figures 8 and 9 for $w$ and $\delta$ exhibit the impacts of the performance index on the optimal control $q^{\ast }$. The result in Figure 8 for not large $x$on the weighting coefficient $w$ is in accordance with the intuition that the deviation $q^{\ast }-\overline{q}$ is smaller for larger value of $w$. Figure 9 gives a linkage between the perspective of the decision-maker and the optimal control $q^{\ast }$. According to the computational result, the decision-maker with longer management viewpoint (smaller $\delta$) allows the deviation to become a larger $q^{\ast }-\overline{q}$ except for large $x$. Qualitatively similar results have been obtained in the asymptotic estimate in Proposition 2.5, again showing consistency between the asymptotic and numerical results. The sub-domain where $q^{\ast }=q_{min}$ is optimal is narrower for larger $\delta$. The above-mentioned results on the optimal control $q^{\ast }$ exploited from Figure 9 suggest that managing the algae population dynamics with smaller deviation $q^{\ast }-\overline{q}$ requires some shorter-term perspective, such that the decision-maker tries to rapidly manage river environment.

Finally, Figure 10 plots the computed value function $\mathrm{\Phi }$ for different values of $a$ in the carrying capacity, implying that $\mathrm{\Phi }$ is strictly increasing and smooth. This computational result indicates that there exist at most a countable number of optimal controls at each point, based on the discussion made just above Remark 2.4.

Figure 10. Comparison of $\mathrm{\Phi }$ for $a=0.45\times l/20$ $(l=0,1,2,\dots ,20)$. The relationship $a{q}_{max}+b=1$ is maintained in all the cases plotted in this figure.

5. Conclusions

This paper formulated a simplified stochastic optimization problem for logistic population dynamics with a control-dependent carrying capacity. Solving the optimization problem was ultimately reduced to finding the unique viscosity solution to the HJB equation. Solving the HJB equation was carried out with a stable finite difference scheme since the solution is not an analytical solution. Numerical solutions reasonably agreed with the asymptotic estimates. Numerical computation of the HJB equation focusing on algae population management revealed strong dependence of the optimal control on the environmental parameters. It was demonstrated that the continuity of the optimal control with respect to the state variable is critically affected by the parameter values, suggesting that they should be carefully calibrated depending on environmental conditions.

In summary, this paper formulated and explored behaviour of the minimum model for the population dynamics with an emphasis on the benthic algae population management in river environment. Although the presented results only concern the simplified model, they are useful for extended problems where controlling some population dynamics plays a central role. For example, the model can be directly incorporated into the operational model for environmental management [63,68]. The present model can also be used as a sub-model for describing predator-prey dynamics arising in inland fisheries management [58]. A more realistic model would involve the delay [7] and spatial-dependence [13,25] and can be numerically simulated with some algorithms, but their mathematical analysis would require more sophisticated techniques like functional analysis for infinite-dimensional dynamical systems. Handling seasonal algae population dynamics is another interesting and practically important problem. The resulting HJB equation for the seasonal-dependent problem is expected to be a time-dependent, parabolic counterpart of the one that focused on in this paper. The established finite difference scheme can be directly extended to discretization of the HJB equation in the advanced problem, which is currently undergoing. Hydraulic and hydrological observations for a more realistic mathematical modelling of the benthic algae population dynamics are undergoing as well.

Acknowledgements

The author thanks the two reviewers for their valuable comments and suggestions.

Disclosure statement

No potential conflict of interest was reported by the author.

ORCID

Hidekazu Yoshioka http://orcid.org/0000-0002-5293-3246

Additional information

Funding

JSPS Research Grant 17K15345 and a grant for ecological survey of a life history of the landlocked ayu Plecoglossus altivelis altivelis from the Ministry of Land, Infrastructure, Transport and Tourism of Japan support this research.

Appendix

Proofs of the propositions stated in the main text are provided in this appendix.

Proof of Proposition 2.1

The proof here is just an application of Theorem 9.8 in Øksendal and Sulem [43], since the SDE (3) admits a path-wise unique solution confined in $D$, and then the continuity of the value function $\mathrm{\Phi }$ in $D$ follows from the Lipschitz continuity of $g$ in $Q\times D$. The fact that $\mathrm{\Phi }$ is a viscosity solution follows from the argument essentially similar to that of the proof for Theorem 9.8 in Øksendal and Sulem [43] combined with the fact that $X_t\in D$ ($t>0$) if $x=X_0\in D$.

Proof of Proposition 2.2

The proof here is based on a standard technique employed with the variable-doubling of the comparison theorems [16,29], but employs a special auxiliary function to deal with the issue of not enforcing boundary conditions at $x=K_{max}$, where the HJB equation (9) is directly solved.

By the positivity of $K$, $r$, and $f$ in $Q$, there exists a constant $0\le \hat{x}\le K_{max}$ close to $K_{max}$ such that (A1) $r\left(q\right)\left(1-\frac{x}{K\left(q\right)}\right)-f\left(q\right)\le 0\mathrm{for}\hat{x}\le x\le K_{max}$(A1) for any $q\in Q$ because (A2) $r\left(q\right)\left(1-\frac{K_{max}}{K\left(q\right)}\right)-f\left(q\right)<0$(A2) for any $q\in Q$. Choose such a $\hat{x}$. Set (A3) $\psi \left(x\right)=(x-\hat{x}{)}^2{\chi }_{\left\{x\ge \hat{x}\right\}},$(A3) which is a smooth, non-negative, increasing, and convex function.

Suppose that $u$ is a viscosity sub-solution. We construct a new viscosity sub-solution from $u$. Set (A4) $u_{\gamma }\left(x\right)=u\left(x\right)-\frac{\gamma }{K_{max}-x}\psi \left(x\right)$(A4) with a constant $\gamma >0$. Then, $u_{\gamma }$ with $u_{\gamma }\left(0\right)=0$ is also a viscosity sub-solution. In fact, if $u\in C\left(D_0\right)\cap USC\left(D\right)$ then $u_{\gamma }\in C\left(D_0\right)\cap USC\left(D\right)$. By (12), for any $x_0\in \hat{D}$, we have (A5) $\begin{array}{rl}& \delta u+\rho \left(u-u\left(\left(1-z\right)x\right)\right)\\ & -\underset{q\in Q}{inf}\left\{-f\left(q\right)x\frac{\mathrm{d}\varphi }{\mathrm{d}x}+r\left(q\right)\left(1-\frac{x}{K\left(q\right)}\right)x\frac{\mathrm{d}\varphi }{\mathrm{d}x}+g\left(q,x\right)\right\}\le 0\mathrm{at}x=x_0\end{array}$(A5) for any test functions $\varphi$ satisfying the conditions in Definition 2.1(a). The relationship (A6) $\frac{\gamma }{K_{max}-x}\psi \left(x\right)\ge \frac{\gamma }{K_{max}-\left(1-z\right)x}\psi \left(\left(1-z\right)x\right),x\in D_0,$(A6) shows (A7) $\begin{array}{rl}& \delta u\left(x_0\right)+\rho \left(u\left(x_0\right)-u\left(\left(1-z\right)x_0\right)\right)\\ & \ge \delta u_{\gamma }\left(x_0\right)+\rho \left(u_{\gamma }\left(x_0\right)-u_{\gamma }\left(\left(1-z\right)x_0\right)\right)\\ & +\rho \left(\frac{\gamma }{K_{max}-x_0}\psi \left(x_0\right)-\frac{\gamma }{K_{max}-\left(1-z\right)x_0}\psi \left(\left(1-z\right)x_0\right)\right)\\ & \ge \delta u_{\gamma }\left(x_0\right)+\rho \left(u_{\gamma }\left(x_0\right)-u_{\gamma }\left(\left(1-z\right)x_0\right)\right).\end{array}$(A7) Set (A8) $W\left(x,p\right)=-\underset{q\in Q}{inf}\left\{\left[r\left(q\right)\left(1-\frac{x}{K\left(q\right)}\right)-f\left(q\right)\right]xp+g\left(q,x\right)\right\}\text{for any}p\in \mathbf{\text{R}}$(A8) which is an increasing function with respect to $p$ for each $x\ge \hat{x}$. In addition, for $x\ge \hat{x}$ and $p\le 0$, we have (A9) $W\left(x,p\right)\le 0\mathrm{and}W\left(x,p\right)\text{is increasing with respect to}p$(A9) since $g\ge 0$. The function $(\gamma /K_{max}-x)\psi \left(x\right)$ is smooth, non-negative, increasing, and (A10) $h\left(x\right)\equiv \frac{\mathrm{d}}{\mathrm{d}x}\left(\frac{\gamma }{K_{max}-x}\psi \left(x\right)\right)\ge 0\left(=0\right)\mathrm{for}x\ge \hat{x}(x<\hat{x}).$(A10) By (A10), we have (A11) $W\left(x,\frac{\mathrm{d}\varphi }{\mathrm{d}x}\right)\ge W\left(x,\frac{\mathrm{d}\varphi }{\mathrm{d}x}-h\right).$(A11) By (A7) and (A11), we further have (A12) $\delta u_{\gamma }+\rho \left(u_{\gamma }-u_{\gamma }\left(\left(1-z\right)x_0\right)\right)+W\left(x,\frac{\mathrm{d}\varphi }{\mathrm{d}x}-h\right)\le 0,$(A12) implying that $u_{\gamma }$ is a viscosity sub-solution.

Now, we claim a comparison theorem: for any viscosity sub-solution $u$ and viscosity-super solution $v$ such that $u\left(0\right)\le v\left(0\right)$, $u\le v$ in $\hat{D}$. The following part of the proof is based on a contradiction argument. For any $x,y\in D$, $\epsilon >0$, set (A13) ${\phi }_{\epsilon ,\gamma }\left(x,y\right)=u_{\gamma }\left(x\right)-v\left(y\right)-\frac{1}{2\epsilon }(x-y{)}^2$(A13) as in the standard variable-doubling technique (Theorem 3.7 of Koike [29], Theorem 8.2 of Crandall et al. [16]). Assume $\underset{D}{sup}\left(u_{\gamma }\left(x\right)-v\left(x\right)\right)={\theta }_{\gamma }>0$ at some $x=\tilde{x}\in D$. By the assumption, we have $\tilde{x}\in \hat{D}$. Without loss of generality, assume $\underset{\gamma \to +0}{liminf}{\theta }_{\gamma }>0$. Otherwise, we have nothing to prove. Choose sufficiently small $\gamma$ such that ${\theta }_{\gamma }\ge \frac{1}{2}\underset{\gamma \to +0}{liminf}{\theta }_{\gamma }>0$.

Let $\left(x_{\epsilon },y_{\epsilon }\right)$ be the maximizer of ${\phi }_{\epsilon ,\gamma }$ in $D\times D$. Then, as in the standard argument of comparison theorems, we can choose a sub-sequence of $\left(x_{\epsilon },y_{\epsilon }\right)$ such that (A14) $\underset{\epsilon \to +0}{lim}{u}_{\gamma }\left(x_{\epsilon }\right)=u_{\gamma }\left(\tilde{x}\right),\underset{\epsilon \to +0}{lim}v\left(y_{\epsilon }\right)=v\left(\tilde{x}\right),$(A14) (A15) $\underset{\epsilon \to +0}{lim}\frac{1}{\epsilon }(x_{\epsilon }-y_{\epsilon }{)}^2=0,$(A15) and (A16) $\underset{\epsilon \to +0}{lim}{x}_{\epsilon }=\underset{\epsilon \to +0}{lim}{y}_{\epsilon }=\tilde{x}.$(A16) By the assumption of $u_{\gamma }\le v$ at $x=0$ and the fact that ${\phi }_{\epsilon ,\gamma }$ diverges to $-\mathrm{\infty }$ at $x=K_{max}$, we must have $0<\tilde{x}<K_{max}$. Assuming $\epsilon$ is small, we can choose $\left(x_{\epsilon },y_{\epsilon }\right)$ such that $0<x_{\epsilon },y_{\epsilon }<K_{max}$, which is hereafter assumed.

By the above argument and Definition 2.1, the functions (A17) $w_1\left(x\right)=u_{\gamma }\left(x_{\epsilon }\right)+{\phi }_{\epsilon ,\gamma }\left(x,y_{\epsilon }\right)-{\phi }_{\epsilon ,\gamma }\left(x_{\epsilon },y_{\epsilon }\right)$(A17) and (A18) $w_2\left(y\right)=v\left(y_{\epsilon }\right)-{\phi }_{\epsilon ,\gamma }\left(x_{\epsilon },y\right)+{\phi }_{\epsilon ,\gamma }\left(x_{\epsilon },y_{\epsilon }\right)$(A18) are test functions for viscosity sub- and super- solutions $u_{\gamma }$ and $v$, respectively. Therefore, we have (A19) $\delta u_{\gamma }+\rho \left(u_{\gamma }-u_{\gamma }\left(\left(1-z\right)x\right)\right)+W\left(x,\frac{x-y_{\epsilon }}{\epsilon }\right)\le 0\mathrm{at}x=x_{\epsilon }$(A19) and (A20) $\delta v+\rho \left(v-v\left(\left(1-z\right)y\right)\right)+W\left(y,\frac{x_{\epsilon }-y}{\epsilon }\right)\ge 0\mathrm{at}y=y_{\epsilon }.$(A20) We then arrive at (A21) $\begin{array}{rl}& \delta u_{\gamma }\left(x_{\epsilon }\right)+\rho \left(u_{\gamma }\left(x_{\epsilon }\right)-u_{\gamma }\left(\left(1-z\right)x_{\epsilon }\right)\right)+W\left(x_{\epsilon },\frac{x_{\epsilon }-y_{\epsilon }}{\epsilon }\right)\\ & \le \delta v\left(y_{\epsilon }\right)+\rho \left(v\left(y_{\epsilon }\right)-v\left(\left(1-z\right)y_{\epsilon }\right)\right)+W\left(y_{\epsilon },\frac{x_{\epsilon }-y_{\epsilon }}{\epsilon }\right),\end{array}$(A21) which is rewritten as (A22) $\begin{array}{rl}& \delta \left(u_{\gamma }\left(x_{\epsilon }\right)-v\left(y_{\epsilon }\right)\right)+\rho \left(u_{\gamma }\left(x_{\epsilon }\right)-v\left(y_{\epsilon }\right)\right)\\ & \le \rho \left(u_{\gamma }\left(\left(1-z\right)x_{\epsilon }\right)-v\left(\left(1-z\right)y_{\epsilon }\right)\right)+W\left(y_{\epsilon },\frac{x_{\epsilon }-y_{\epsilon }}{\epsilon }\right)-W\left(x_{\epsilon },\frac{x_{\epsilon }-y_{\epsilon }}{\epsilon }\right).\end{array}$(A22) By the definition of $\left(x_{\epsilon },y_{\epsilon }\right)$, we obtain (A23) $\begin{array}{rl}& u_{\gamma }\left(\left(1-z\right)x_{\epsilon }\right)-v\left(\left(1-z\right)y_{\epsilon }\right)={\phi }_{\epsilon ,\gamma }\left(\left(1-z\right)x_{\epsilon },\left(1-z\right)y_{\epsilon }\right)+\frac{1}{2\epsilon }(1-z{)}^2(x_{\epsilon }-y_{\epsilon }{)}^2\\ & \le {\phi }_{\epsilon ,\gamma }\left(x_{\epsilon },y_{\epsilon }\right)+\frac{1}{2\epsilon }(1-z{)}^2(x_{\epsilon }-y_{\epsilon }{)}^2\\ & =u_{\gamma }\left(x_{\epsilon }\right)-v\left(y_{\epsilon }\right)-\left(1-{\left(1-z\right)}^2\right)\frac{1}{2\epsilon }(x_{\epsilon }-y_{\epsilon }{)}^2\\ & \le u_{\gamma }\left(x_{\epsilon }\right)-v\left(y_{\epsilon }\right).\end{array}$(A23) Substituting (A23) into (A22) yields (A24) $\delta \left(u_{\gamma }\left(x_{\epsilon }\right)-v\left(y_{\epsilon }\right)\right)=\delta {\theta }_{\gamma }\le W\left(y_{\epsilon },\frac{x_{\epsilon }-y_{\epsilon }}{\epsilon }\right)-W\left(x_{\epsilon },\frac{x_{\epsilon }-y_{\epsilon }}{\epsilon }\right).$(A24) Because $W$ is Lipschitz continuous with respect to $x\in D$ and $p\in \mathbf{\text{{R}}}$, taking the limit $\epsilon \to +0$ in (A24) yields (A25) $0\ge \delta {\theta }_{\gamma }\ge \frac{\delta }{2}\underset{\gamma \to +0}{liminf}{\theta }_{\gamma }>0.$(A25) This is a contradiction since $\delta >0$.

Uniqueness of viscosity solutions immediately follows from the above-presented comparison argument and their continuity. This is because a viscosity solution is continuous in $D$ and it is a viscosity sub-solution as well as a viscosity super-solution.

Proof of Proposition 2.4

By the assumption, for small $x$, the HJB equation (23) reduces to (A26) $\delta \mathrm{\Phi }-rx\frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}x}-x^m+\rho \left(\mathrm{\Phi }-\mathrm{\Phi }\left(\left(1-z\right)x\right)\right)-\underset{q\in Q}{inf}\left\{-\alpha qx\frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}x}+\frac{w}{2}{\left(q-\overline{q}\right)}^2\right\}\sim 0.$(A26) Assuming the asymptotic form $\mathrm{\Phi }\sim c{x}^n$ with some $c\ne 0$ and $n>0$, substituting this into (A26) yields (A27) $\left(\delta -nr+\rho {\left(1-\left(1-z\right)\right)}^n\right)c{x}^n-x^m-\underset{q\in Q}{inf}\left\{-\alpha qcn{x}^n+\frac{w}{2}{\left(q-\overline{q}\right)}^2\right\}\sim 0.$(A27) Now, we have (A28) $\underset{q\in Q}{\arg min}\left\{-\alpha qcn{x}^n+\frac{w}{2}{\left(q-\overline{q}\right)}^2\right\}=\overline{q}+c\frac{1}{w}\alpha cn{x}^n.$(A28) Substituting (A28) into (A27) yields (A29) $Tc{x}^n-\frac{c^2{n}^2{\alpha }^2}{2w}{x}^{2n}\sim x^m.$(A29) Matching the coefficients in both-hand sides of (A29) yields $n=m$ and thus $Tc{x}^m\sim x^m$, implying $c=T^{-1}=c_0$, which gives (25).

For the asymptotic estimate of $q^{\ast }$, we use (25). By the definition, we have (A30) $q^{\ast }\left(x\right)=\underset{q\in Q}{\arg min}\left\{-\alpha qx\frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}x}-\frac{r{x}^2}{aq+b}\frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}x}+\frac{w}{2}{\left(q-\overline{q}\right)}^2\right\}.$(A30) Substituting (25) into (A30) yields (A31) $q^{\ast }\left(x\right)=\underset{q\in Q}{\arg min}\left\{-qm\alpha c_0{x}^m-\frac{r{c}_0{x}^{m+1}}{aq+b}+\frac{w}{2}{\left(q-\overline{q}\right)}^2\right\}.$(A31) Then, a simple differentiation gives (A32) $-m\alpha c_0{x}^m+\frac{r{c}_0{x}^{m+1}}{{\left(a{q}^{\ast }+b\right)}^2}+w\left(q^{\ast }-\overline{q}\right)=0,$(A32) which can be rewritten as (A33) $q^{\ast }-\overline{q}=\frac{m\alpha }{w}{c}_0{x}^m-\frac{r{c}_0}{w{a}^2}{\left(q^{\ast }+\frac{b}{a}\right)}^{-2}{x}^{m+1}.$(A33) The equality (A33) shows $q^{\ast }\to 0$ as $x\to 0$. Therefore, we can use the expansion (A34) $\begin{array}{rl}{\left(q^{\ast }+\frac{b}{a}\right)}^{-2}& ={\left(q^{\ast }-\overline{q}+\overline{q}+\frac{b}{a}\right)}^{-2}\\ & \sim {\left(\overline{q}+\frac{b}{a}\right)}^{-2}\left[1-2{\left(\overline{q}+\frac{b}{a}\right)}^{-1}\left(q^{\ast }-\overline{q}\right)\right].\end{array}$(A34) Then, (A33) and (A34) gives (A35) $\begin{array}{rl}{q}^{\ast }-\overline{q}& \sim \frac{m\alpha }{w}{c}_0{x}^m-\frac{r{c}_0}{w{a}^2}{\left(\overline{q}+\frac{b}{a}\right)}^{-2}{x}^{m+1}\\ & +2\frac{r{c}_0}{w{a}^2}{\left(\overline{q}+\frac{b}{a}\right)}^{-3}\left(q^{\ast }-\overline{q}\right)x^{m+1}\end{array}$(A35) and thus (A36) $\begin{array}{rl}{q}^{\ast }-\overline{q}& \sim \frac{m\alpha }{w}{c}_0{x}^m-\frac{r{c}_0}{w{a}^2}{\left(\overline{q}+\frac{b}{a}\right)}^{-2}{x}^{m+1}+2\frac{r{c}_0}{w{a}^2}{\left(\overline{q}+\frac{b}{a}\right)}^{-3}\left(\frac{m\alpha }{w}{c}_0{x}^m\right)x^{m+1}\\ & =\frac{m\alpha }{w}{c}_0{x}^m-\frac{r{c}_0}{w{a}^2}{\left(\overline{q}+\frac{b}{a}\right)}^{-2}{x}^{m+1}\left(1-2{\left(\overline{q}+\frac{b}{a}\right)}^{-1}\frac{m\alpha }{w}{c}_0{x}^m\right),\end{array}$(A36) leading to (26).

Proof of Proposition 3.1

The proof here is essentially the same with but simpler than that in Chapters 3 and 4 of Azimzadeh [5] since the present HJB equation does not involve the nonlinear non-local term. To prove the convergence, we show the monotonicity, stability, and consistency of the scheme. Monotonicity of the scheme follows from the fact that the matrix $A-B$ is an M-matrix. Consistency of the scheme is straightforwardly proven just following Lemma 4.3.16 of Azimzadeh [5].

Stability of the scheme, namely uniform boundedness of numerical solutions independent from the discretization parameter $l$, is proven as follows. Set $\overline{\mathrm{\Phi }}=\underset{0\le i\le I}{max}|{\mathrm{\Phi }}_i|$. Assume that the maximum is achieved at some $i_0$ such that $0<i_0<I$. Then, by (41), we have (A37) $-{\mu }_{i_0,i_0-1}{\mathrm{\Phi }}_{i_0-1}+{\mu }_{i_0,i_0}{\mathrm{\Phi }}_{i_0}-{\mu }_{i_0,i_0+1}{\mathrm{\Phi }}_{i_0+1}=B_{i_0,j_{i_0}}{\mathrm{\Phi }}_{j_{i_0}}+B_{i_0,j_{i_0}+1}{\mathrm{\Phi }}_{j_{i_0}+1}+S_{i_0},$(A37) which leads to (A38) $\begin{array}{rl}-{\mu }_{i_0,i_0-1}{\mathrm{\Phi }}_{i_0-1}+{\mu }_{i_0,i_0}{\mathrm{\Phi }}_{i_0}-{\mu }_{i_0,i_0+1}{\mathrm{\Phi }}_{i_0+1}& =-{\mu }_{i_0,i_0-1}{\mathrm{\Phi }}_{i_0-1}+{\mu }_{i_0,i_0}\overline{\mathrm{\Phi }}-{\mu }_{i_0,i_0+1}{\mathrm{\Phi }}_{i_0+1}\\ & \ge \left({\mu }_{i_0,i_0}-{\mu }_{i_0,i_0-1}-{\mu }_{i_0,i_0+1}\right)\overline{\mathrm{\Phi }}\end{array}$(A38) and (A39) $\begin{array}{rl}{B}_{i_0,j_{i_0}}{\mathrm{\Phi }}_{j_{i_0}}+B_{i_0,j_{i_0}+1}{\mathrm{\Phi }}_{j_{i_0}+1}+S_{i_0}& \le \left(B_{i_0,j_{i_0}}+B_{i_0,j_{i_0}+1}\right)\overline{\mathrm{\Phi }}+S_{i_0}\\ & =\rho \overline{\mathrm{\Phi }}+S_{i_0}\\ & \le \rho \overline{\mathrm{\Phi }}+\frac{w}{2}max\{q_{max}-\overline{q},\overline{q}-q_{min}{\}}^2+1.\end{array}$(A39)

Combining (A38) and (A39) yields (A40) $\left({\mu }_{i_0,i_0}-{\mu }_{i_0,i_0-1}-{\mu }_{i_0,i_0+1}-\rho \right)\overline{\mathrm{\Phi }}\le \frac{w}{2}max\{q_{max}-\overline{q},\overline{q}-q_{min}{\}}^2+1.$(A40) Since (A41) $\begin{array}{rl}{\mu }_{i_0,i_0}-{\mu }_{i_0,i_0-1}-{\mu }_{i_0,i_0+1}-\rho & =\delta +\frac{P_i}{e^{P_i}-1}\frac{r{x}_i}{l}+\frac{\alpha q_i^{\ast }{x}_i}{l}+\frac{r{x}_i^2}{\left(a{q}_i^{\ast }+b\right)l}+\rho \\ & -\left(\frac{\alpha q_i^{\ast }{x}_i}{l}+\frac{r{x}_i^2}{\left(a{q}_i^{\ast }+b\right)l}\right)-\frac{P_i}{e^{P_i}-1}\frac{r{x}_i}{l}-\rho \\ & =\delta ,\end{array}$(A41) substituting (A41) into (A40) yields (A42) $\overline{\mathrm{\Phi }}\le \frac{1}{\delta }\left(\frac{w}{2}max{\left\{q_{max}-\overline{q},\overline{q}-q_{min}\right\}}^2+1\right),$(A42) showing the stability. If $i_0=I$, then an essentially similar argument again leads to (A42). If $i_0=0$, then the stability is trivial.

In summary, the scheme is monotone, consistent, and stable. Therefore, a straightforward application of the convergence theorem (Barles and Souganidis [8], Theorem 4.3.19 of Azimzadeh [5]) combined with Proposition 2.2 leads to the convergence of the scheme.