Mathematical and numerical analyses of a stochastic impulse control model with imperfect interventions

A stochastic impulse control problem with imperfect controllability of interventions is formulated with an emphasis on applications to ecological and environmental management problems. The imperfectness comes from uncertainties with respect to the magnitude of interventions. Our model is based on a dynamic programming formalism to impulsively control a 1-D diffusion process of a geometric Brownian type. The imperfectness leads to a non-local operator different from the many conventional ones, and evokes a slightly different optimal intervention policy. We give viscosity characterizations of the Hamilton–Jacobi–Bellman Quasi-Variational Inequality (HJBQVI) governing the value function focusing on its numerical computation. Uniqueness and verification results of the HJBQVI are presented and a candidate exact solution is constructed. The HJBQVI is solved with the two different numerical methods, an ordinary differential equation (ODE) based method and a finite difference scheme, demonstrating their consistency. Furthermore, the resulting controlled dynamics are extensively analyzed focusing on a bird population management case from a statistical standpoint.


Introduction
Balancing costs and benefits is essential for successful environmental and ecological management [1,2]. Harvesting invasive species requires measuring invasion dynamics, evaluating costs of control efforts, and tracking invasion damages [3]. Cooperative lake eutrophication management considers the costs for water quality improvement and their benefits under social pressure [4]. Transboundary pollution problems are regional problems where the pollution mitigation policy is critically affected by the costs and benefits evaluation [5]. Smith et al. [6] discussed uncertainties in wildlife management focusing on problems related to waterbird.
Stochastic optimal control [7] is a central mathematical concept for modeling, analyzing, and controlling stochastically-driven system dynamics. Its application areas are widely distributed in both science [8,9] and engineering [10,11]. One effective way for solving a stochastic optimal control problem is to derive the optimality equation, often called Hamilton-Jacobi-Bellman (HJB) equation, from which we can construct the optimal control [12,13]. The optimality equations are solved analytically in simple problems [7,13,14], while they must be approximated in complex problems using numerical methods [15][16][17]. Management problems related to environment and ecology, especially those on population dynamics that are inherently stochastic [18], are the ones that the stochastic control theory can handle. Biological resource exploitation considering renewal dynamics [19], invasive and nuisance species management [20][21][22], and ecosystem and biodiversity management [23,24] are such examples.
Management problems of population dynamics through harvesting by human have been efficiently described with the impulse control formalism [7] whose control variables are the timing and amount of harvesting. Impulse control models have been applied to various problems, such as predator-prey population management [25], fishery resources management [26], aquaculture operation [27], tree harvesting [28,29], and waterbird management [30]. Deep mathematical analysis of the impulse control models of population dynamics and related models have been carried out [31,32]. Other important industrial problems, such as gene regulation [33] and satellite control [34,35], have been considered with the impulse control as well.
Models based on the impulse control have thus been attractive mathematical tools for management of population dynamics. However, to the best of the authors' knowledge, they lack the important fact that the impacts of interventions on the population dynamics are not always realized as planned by the decision-maker, but are often uncertain as seen in the field survey results [36]. Loosely speaking, the same effort of interventions does not always lead to the same impacts on the population dynamics. Korn [37] considered an imperfect impulse control problem of diffusion processes, and Helmes et al. [38] generalized it in an abstract manner for continuous-time Markov processes. Considering the uncertainty, namely imperfectness of interventions, is a natural policy in environmental and ecological management. However, such a mathematical approach has not been made so far. Impulse control models in other research areas [39][40][41] also do not handle cases where the impulsive interventions have uncertain impacts. This research background motivates us to formulate and analyze an impulse control problem of population dynamics subject to imperfectness interventions.
The objectives of this paper are formulation and analysis of an impulse control problem of population dynamics subject to imperfect interventions, focusing on its application to a recent waterbird management problem [21,30]. Our model is based on Korn [37], but considers a specific problem where a deeper mathematical and numerical analysis are possible. In addition, our contributions include derivation of a candidate exact solution and numerical computation of the stochastic control problem. A 1-D stochastic differential equation (SDE) describes the population dynamics in a habitat. The population is regulated through costly and imperfect impulsive interventions. Optimizing the interventions is achieved through minimization of a performance index considering both costs and benefits of the interventions. An application of the dynamic programming principle to the system dynamics and the performance index leads to the optimality equation as a nonlinear complementary problem called the Hamilton-Jacobi-Bellman Quasi Variational Inequality (HJBQVI) [7]. The HJBQVI has a different non-locality from the conventional ones due to the imperfectness of the interventions. This gives rise to an optimal policy having a slightly more complicated form than the conventional ones. We give viscosity characterizations of the (HJBQVI) [42,43] and that uniqueness and verification results hold true for the HJBQVI. A candidate exact solutions is constructed analytically, and a system of nonlinear equations governing its unknown coefficients is derived. Validity of the candidate exact solution is checked through comparison with numerical solutions generated by a steady counterpart of a formally first-order finite difference scheme [44], and a sharper convergence estimate against our HJBQVI is obtained. Optimally controlled population dynamics is extensively analyzed as well from a statistical standpoint. Our model is relatively simple, but potentially serves as a foundation for analysis of imperfect impulse control problems. Our contribution is rather from an engineering side but contains a new application of the impulse control.
The rest of this paper is organized as follows. In Sect. 2, we present the impulse control model and derives the HJBQVI. Mathematical analysis results on the value function, HJBQVI, and the controlled dynamics are presented in Sect. 3. In Sect. 4, the HJBQVI is solved numerically by two different methods, and we see that their results agree well with each other. Sensitivity analysis of the optimal control and the controlled dynamics is carried out in this section as well. Finally, in Sect. 5, we summarize the results of this paper and present our future perspectives. Proofs of Propositions, which are a bit technical, are placed in the Appendix.

Mathematical model
An SDE for describing single-species population dynamics with imperfect impulsive harvesting is presented, and the performance index to be minimized through choosing harvesting policies is formulated. The HJBQVI of the present optimal control problem is finally derived. Our mathematical formulation is based on Korn [37], but handles a new specific problem related to ecological and environmental management.

Stochastic differential equation
The time is denoted as t. A 1-D standard Brownian motion defined on a complete probability space is denoted as B t at time t [7]. A 1-D SDE governs an adapted càdlàg stochastic process (X t ) t≥0 representing the population. Let τ i (i = 0, 1, 2, . . . , τ 0 = 0) be the time at which the intervention is performed to harvest the population. Without significant loss of generality, we assume that the sequence {τ i } i=0,1,2,... is increasing. The amount of population harvested at the time τ i is denoted as ζ i with 0 < ζ i ≤ X τ i -(i ≥ 1) and ζ 0 = 0. The couple of sequences η = {τ i , ζ i } i=0,1,2,... , which is called policy, is optimized so that a performance index is minimized.
We are interested in an imperfect impulse control problem in which the harvesting is not always realized as planned by the decision-maker, the manager of the population, but may be perturbed stochastically. Let Z = {z i } i=0,1,2,... be a sequence of i.i.d stochastic variables valued in a compact set S = [a, b] with 0 < a < b ≤ 1 having the probability density function (PDF) q. Z = {z i } i=0,1,2,... is independent from (B t ) t≥0 . The decision-maker can know the value of z i just after the harvesting at τ i .

Remark 1
The assumption 0 < a < b ≤ 1 means that over-harvesting such that the population is harvested more than that planned by the decision-maker does not occur. In this sense, the present model is more pessimistic than the conventional impulse control models with perfect interventions.
A policy η is called admissible if it complies with the following conditions, which are natural requirements to well-pose the problem (see, Definition 2.1 of Korn [37]): The conditions of admissible policy τ 0 = 0 and ζ 0 = 0, τ i is a stopping time with respect to the filtration generated by ζ i is measurable with respect to the filtration generated by σ X τ n -, (τ n , ζ n , z n ), n < i , Pr lim Here, Pr[·] and E[·] present the probability and expectation, respectively. The collection of all admissible policies, namely policies complying with (2) through (7), is denoted as A. This set depends on X 0-, and is denoted as A X 0-when necessary. Here, X t-is the left-limit of X at time t. For each η ∈ A, the SDE governing (X t ) t≥0 is set as subject to an initial condition X 0-= x ≥ 0, where the fluctuation term σ X t dB t is in the Itô's sense. Here, μ > 0 is the intrinsic growth rate and σ > 0 is the stochastic growth rate modulating temporal fluctuation involved in the dynamics [21]. Since we are interested in controlling the population through harvesting, assume 2μ > σ 2 so that X t → +∞ a.s. as t → +∞ if there is no intervention [21]. This condition is utilized in the statistical analysis of Sect. 3.3. Due to the global Lipschitz continuity of the coefficients of the SDE (8), it admits a unique non-negative solution in a path-wise sense in each τ i ≤ t ≤ τ i+1 when τ i < τ i+1 . In addition, the SDE (8) admits a unique non-negative solution in a path-wise sense global in time since X τ i = X τ i-z i ζ i ≥ 0. The second line of (8) means that the decisionmaker choses ζ i based on X τ i-. Due to the uncertain perturbation z i not controllable by the decision-maker, the realized state transition at the time τ i is X τ i-→ X τ i-z i ζ i . We assume that the decision-maker knows the PDF q. He/she thus wants to decide a policy η ∈ A that can minimize a performance index under the imperfectness of harvesting.

Performance index
Following the previous research of controlling population dynamics [21,30], the perfor- is the expected net cost depending on the initial condition: with and where δ > 0 is the discount rate, R, r, M and m are the positive constants. In (10), K(ζ ) is the cost of harvesting with the amount of ζ > 0, where k 1 ζ with k 1 > 0 is the proportional cost and k 0 > 0 is the fixed cost incurred regardless of ζ .
The performance index J in (9) consists of the two terms: the cumulative total disutility (first term, the integral) caused by the existence of the population and the total cost of harvesting the population (second term, the sum). In the first term, -RX M s represents the ecological utility provided by the population, and rX m s represents the disutility caused by the population [21]. The negativity of the first term is due to measuring the cost to be positive in the present model. The sum -RX M s + rX m s is unimodal and convex with respect to X s ≥ 0, and is bounded from below by a negative constant as The relationship (11) means that the rate of increase of the disutility caused by the population is higher than that of the utility provided by the population. The condition 0 < M < 1 in (11) means that the ecosystem services provided by the population [45,46] drastically change between the states that the population exists (X s > 0) and that the population does not exist (X s = 0). The assumption m ≤ 2 is to guarantee a square-integrability of the term rX m s in the performance index. As in the conventional impulse control problems, we assume that δ is sufficiently large so that J is bounded [47][48][49]. To guarantee this, we assume In fact, by (13), we obtain for arbitrary η ∈ A. We also have for arbitrary η ∈ A because the function -Rx M + rx m (x ≥ 0) is convex and has a minimum value.

Value function and HJBQVI
The value function = (x) is the minimized performance index with respect to η ∈ A: where η * = {τ * i , ζ * i } i=0,1,2,... ∈ A minimizing J is referred to as the optimal policy. This is bounded from below and locally bounded from above, and is well-defined. These are because, firstly Secondly, with the null control η = η 0 such that τ i = +∞ (i ≥ 1), which is the control with no intervention, (X t ) t≥0 is a classical geometric Brownian motion. Then, we get In addition, (0) = 0 because X t = 0 (t ≥ 0) if x = 0 and only the null control η 0 is admissible in this case. The boundedness results of are sharpened later. We assume the following recursion inspired from the dynamic programming principle of the conventional impulse control problems [7]. Here, τ is any finite stopping time Page 7 of 34 adapted to a natural filtration F generated by (B t ) t≥0 and {z i } i=1,2,3,... : which formally leads to the HJBQVI governing : with the degenerate elliptic operator L and the non-local intervention operator M and respectively, for generic regular W = W (x). The boundary condition for the HJBQVI is (0) = 0. Growth speed of (viscosity) solutions to (20) should be at most polynomial by (17) and (18).
The left part in "max" operator in the HJBQVI (20) corresponds to the situation where no intervention should be performed, while the right part corresponds to the situation where the intervention should be performed immediately. A difference between the present and conventional HJBQVIs is that the present non-local operator M involves the effect of uncertainty induced by the uncertainties Z, while the conventional ones do not and have the non-local terms like MW = sup 0<ζ ≤x [K(ζ ) + W (xζ )] in which no uncertainty is considered [30,[39][40][41].

Mathematical analysis 3.1 Value function
Several key mathematical properties of the value function are analyzed, from which we show that it is a viscosity solution to the HJBQVI (20). A comparison result for a modified problem that has important implications in numerical computation is also presented.
The first result concerns local boundedness of as an extension of the discussion in the previous section.
The following proposition gives a key inequality characterizing the value function , which is effectively utilized in a proof of viscosity property of .
In this paper, we assume that the value function is continuous, which is true in the problems without uncertainties in the interventions [30,50,51]. The continuity is also justified numerically as the computational results with a finite difference scheme suggest (see, Sect. 4). The following continuity result supports the definitions of viscosity solutions to the HJBQVI (20) presented later. Actually, the candidate of exact solution, which turns out to be the value function under certain conditions, is continuous.

Proposition 3 There exists a constant C > 0 such that
In addition, (x) is continuous at x = 0.
Remark 2 An obstacle to derive the other side of the inequality is the dependence of the set A x of admissible controls on x. Notice that this problem has not been encountered in the conventional research such as Guo and Wu [48] because their admissible sets are independent from the state variables.

Viscosity property
Here, we show that the HJBQVI (20) is solvable in a weak sense. The weak solutions employed in this paper are viscosity solutions [42,43] as appropriate weak solutions to degenerate elliptic problems. We show that the value function is a continuous viscosity solution to the HJBQVI (20) under certain conditions. As noted above, hereafter we assume that the value function is continuous, which is partly supported in Proposition 3. Viscosity solutions to the HJBQVI (20) are defined as follows based on the monotonicity and continuity result of non-local terms in HJBQVIs (Proposition 3 of Azimzadeh et al. [42]) whose conditions are satisfied by our non-local operator M. In what follows the function space of upper-semi continuous (resp., lower-semi continuous) functions valued in a domain is denoted as USC( ) (resp., LSC( )). Note that the difference between a test function and a viscosity sub-solution or a super-solution needs not be globally maximized or minimized because the non-local part is defined using super-or sub-solution itself (e.g., [42,Definition 1]).
The following proposition shows the viscosity property of the value function and supports our approach of finding the optimal policy through solving the HJBQVI (20).

Proposition 4 The value function is a viscosity solution to the HJBQVI (20).
Uniqueness of solutions to our HJBQVI is an important issue as well, from which we can guarantee unique existence of viscosity solutions under certain assumptions, and also convergence of numerical solutions generated by an appropriate discretization method. An interest, especially in numerical computation, is uniqueness of the localized HJBQVI: with (0) = 0 and a sufficiently large constant g > 0 such that g > AX m max . Here, X max > 0 is a constant to truncate the unbounded domain [0, +∞) as [0, X max ]. This kind of localization has been commonly used in solving degenerate elliptic problems [49,52,53]. Here, introducing the termg is from a technical reason (see, Proof of Proposition 5) to guarantee uniqueness of the system and actually innocuous because of the estimate (23).
Following Definition 1, viscosity solutions to the localized HJBQVI are defined as follows: A difference between Definitions 1 and 2 is the boundary condition at the right-end. (See, (31) and (33). They are absent in Definition 1.) • Viscosity super-solution. A function ϕ ∈ LSC[0, X max ] ∩ C(0, X max ] with ϕ(0) ≥ 0 is a viscosity super-solution to the localized HJBQVI (28)- (29) if the following conditions are satisfied at each x > 0: for any ψ ∈ C 2 [0, X max ] such that ϕψ is locally strictly

Verification
In this sub-section, we present a candidate exact solution to the HJBQVI (20). From a management viewpoint [21,30], we conjecture the following with somex > 0: and The point x =x is a free boundary. The conditions (34) and (35) imply that the optimal policy is of the following form, which is of course admissible, and is referred to as QVIcontrol [37,54]: From (36) through (38),we infer that the following threshold-type policy with the threshold valuesx and x (x > x) is optimal for t ≥ 0 (see, Korn [37]): (A) X t-<x, then no intervention is performed. (B) If X t-≥x, then the intervention with the magnitude of ζ =xx(x) is immediately performed (and then X t-is actually reduced to X t =x -(xx(X t-))z, where z is uncertain for the decision-maker). Notice that we may have X t ≥x. Then, repeat the first sentence of (B) until we get X t <x. Notice that multiple interventions are allowed at the same time [37]. This is due to the model simplification that assumes the time scale of interventions is much smaller than that of the population dynamics. Here,x > 0 is a constant while x(x) is a function of x > 0 and is assumed to be 0 < x(x) <x (x ≥x). With the two thresholdsx and x(x), the optimal amount of the intervention ζ * becomes Figure 1 presents an image of sample paths by the above-presented threshold-type policy.
The function x(·) is the population just after the harvesting. It is related to equation (38) in the way that ζ =xx(x) is its maximizer given x ≥ 0 under the conjectures (A)-(B). An example of x(·) is plotted in Fig. 2 of Sect. 4.2.
Note that the conventional 1-D impulse control problems in an infinite horizon assume constant thresholds [44,50,[54][55][56], namelyx and x(x) are positive constants, while in ours the latter is not a constant. If both the thresholds are constant, then (39) becomes a policy with a constant x: However, later we imply that this is not always true in our model, through an exploration of exact solutions. This difference between the conventional and our models is due to the form of the non-local term in the HJBQVI (20 and For each given x ≥ 0, policy η ∈ A, and stopping time τ , suppose that φ(X τ ) is uniformly integrable and satisfies the bound (23). Then, we have In addition, if the QVI-control is admissible, then

Candidate of exact solutions
We try to construct a viscosity solution potentially complying with (34) and (35). Set Let¯ be a candidate of viscosity solutions. Firstly, we consider the case x <x, under which no intervention should be performed. General solutions to are represented as with and constants c 1 and c 2 . If c 2 = 0, then lim x→+0 c 2 x β -= ±∞ since β -< 0. Hereafter, we write β = β + . This contradicts with Proposition 1. Therefore, c 2 = 0 and thus where the notations c 1 = c is used and¯ (0) = 0.
Remark 4 By Proposition 1, we should have c ≤ 0. For x ≥x where the intervention must be performed immediately, we should havē whose right-hand side is determined based on the information in [0, x]. Therefore, theoretically,¯ for x ≥x can be constructed from the left toward the right by the equation (50) once we get¯ (x). We found that the equation (50) is not analytically solvable because of its nonlinearity and nonlocality. Nevertheless, assuming that the population is not large at some time t, as demonstrated below, it is sufficient to solve the HJBQVI (20) for x ≤x considering x =x as a free boundary: From the standpoint of the optimal policy constructed above, it is sufficient to have the point value ζ * (x) but not whole ζ * (x)(x ≥x) when the initial population x is not so large. This is because of the rules (A)-(B) of the optimal policy with which the process (X t ) t≥0 is eventually confined in [0,x]. We thus only have to determine the three constants ζ * (x), x, and c. For the sake of simplicity of presentations, writeζ = ζ * (x). We assume that ζ * is sufficiently smooth. Now, we derive three equations to find the three unknownsζ ,x, and c. Assume x ≥x. By the Value matching condition (a), we get The optimality condition then gives the implicit governing equation of ζ = ζ * , if it is an interior minimizer, as The Smooth pasting condition (b) gives Combining (53) and (54) leads to Now, substituting x =x into (51), (53), and (55) with the Optimality condition (c) gives the system of three nonlinear equations to find the three unknownsζ ,x, and c: and Remark 5 We can expect an asymptotic relationship (x) ∼ dk 1 x with a constant d > 0 for x sufficiently large. In fact, substituting (x) = dk 1 x into (53) leads to This is consistent with (x) = M (x) for x sufficiently large because of by (61). From a practical viewpoint, (61) means that the increasing rate d with respect to the population x is not smaller than k 1 because of the uncertainty of the intervention impacts. We get d = 1 if and only if the interventions are perfect (a = b = 1 and q is a Dirac's delta concentrated at z = 1). Therefore, the theoretical result implies that imperfect interventions are always less efficient than the perfect ones, and the efficiency can be measured by d. The modulation of the increasing rate is the inverse of the mean value of the uncertainty Z by (61). The coefficient d is determined by the moment of q. The result presented above is consistent with our intuition that the intervention is n (> 1) times less efficient if its impact is expected to be n times smaller than that without the uncertainty.

Empirical numerical approach
The nonlinear equations (56) through (58) seem to be not analytically solvable and their numerical approximation in general may not be an easy task. We consider a simplified case where we can, at least empirically, find the unknownsζ ,x, and c. A similar numerical approach has been successfully applied to finding the coefficients in the exact solution to an HJBQVI associated with a problem with a perfect intervention [51]. We consider the uniform distribution q(z) ≤ z ≤ b). In this case, by a straightforward calculation, we get and b a d dx 0 (xzζ )zq(z) dz Then, (57) and (58) are rewritten as and respectively. The three unknowns are explored through the following empirically-derived system containing the two ODEs and a nonlinear equation whose equilibrium satisfy (56), (65), and (66): where γ > 0 is a constant and the coefficientx γ in the right-hand side of (67) is chosen for the sake of empirical computational stability of the discretized system. Here, ρ ≥ 0 is an auxiliary time variable to parameterizeζ ,x, and c, and the ODEs (67) and (68) are integrated with respect to ρ > 0 with (69) starting from an initial guesses (x,ζ , c) τ =0 . This is carried out in a standard forward Euler manner. Set the small increment ρ > 0 of ρ and a variable at ρ = k ρ (k = 0, 1, 2, . . .) is denoted with the subscript (k). For k ≥ 0, we discretize (67) through (69) as where 0,(k) is 0 with c = c (k) . From a stability viewpoint of discretized ODEs [57], we expect that the discretized system (70)-(72) can be solved numerically with a small ρ > 0.
Remark 6 We observed that convergence of solutions to (70)-(72) depends on the initial guess, implying its sensitivity against the initial guess. At the current stage, we have neither unique solvability nor stability results on the continuous system (67)-(69) and discretized system (70)- (72). Therefore, we utilize the ODE-based method as a tool to validate the finite difference scheme because the latter is more computationally stable for our problem.
Remark 7 Construction of the system is based on an empirical approach. Yaegashi et al. [51] addressed this issue for an impulse control problem with perfect interventions. Fortunately, their ODE-based method can compute the unknown coefficients as a locally asymptotically stable equilibrium point of a system. However, such a result has not been found in the present case.

Statistics of controlled state dynamics
From a practical point of view, analyzing statistics of the optimally-controlled population dynamics is as important as finding the optimal intervention policy. Under the conjectured optimal policy in the previous sub-sections, we can focus on the dynamics in the interval [0,x]. The threshold-type optimal control established in the previous section is assumed in what follows. Firstly, we derive the statistical moments of the intervals τ i+1τ i between each successive intervention. This is a stochastic variable of interest in applications, because the frequency of interventions crucially affects cost-effectiveness of the population management. The present problem is a stochastic optimal control of a time-homogenous process in an infinite horizon, and thus the statistical properties of the optimally controlled dynamics can be studied through stationary distributions.
The first hitting time τ x of the process (X t ) t≥0 (t > 0) with X 0-= x ∈ (0,x] is defined as The Laplace's transformation method (Appendix of Chap. 9, Dixit et al. [58]) gives the following exponential formula because our underlying process is a geometric Brownian motion: with λ = 2μ σ 2 -1(> 0) and u = for each θ > 0, where E x is the conditional expectation on X 0-= x. The formula of the probabilistically distributed x follows directly from (74) as where τ x is rewritten as τ , which is now the first hitting time to the thresholdx just after an intervention whose exertion time is 0. We have and thus By (77) and (78), we get and The equation (79) (79) and (80) are calculated exactly or evaluated numerically depending on the PDF q. Higher order moments of τ can also be computed based on (76) if necessary. Another statistical index of interest is the mean harvested amount per unit time [21], which is denoted as H. This quantity can be simply estimated as where the numerator of the right-hand side is the mean unit-time harvesting. Validity of this formula is examined in the next section with a Monte-Carlo method.

Numerical computation
The HJBQVI (20) is computed with the two different numerical methods, which is the ODE-based method using (70) through (72) and a verified finite difference scheme having a formally first-order accuracy [44]. In addition, the statistical indices of the controlled state dynamics are calculated.

Parameter setting
A target species of the model application is the great cormorant Phalacrocorax carbo, which is one of the most widespread waterfowls found over the world whose population management has been serious ecological issues because of their excessively high predation to fishery resources [21,59,60]. They have indirect environmental impacts through modulating nutrient cycles around their habitat [61] through complex mechanisms. Nevertheless, exterminating the bird is not be the optimal solution because they are not alien species in many cases. Their population dynamics would present seasonal characteristics; however, from a practical modeling viewpoint, the dynamics can be effectively considered to be time-homogenous with the deterministic and stochastic growth as modelled in the SDE (8). Human interventions to suppress their population has been carried out in many ways, and the most popular one is gun-shooting [21,36].  -a (a ≤ z ≤ b) is assumed. The specified range [a, b] is reasonably consistent with the survey result that the total number of harvested bird population at one harvesting event had different orders (O(10 0 ) to O(10 2 )) among hunters [36]. The computational domain for the finite difference scheme is [0, X max ] with X max = 100, which is discretized into 1000 cells and 1001 nodes. The threshold X max is determined so that X max >x in what follows based on preliminary numerical experiments. These parameter values and the computational resolution are used unless otherwise specified.

Value function and optimal policy
The HJBQVI (20) is discretized with the finite difference scheme [44] that has been found to be monotone, stable, consistent, and thus convergent in the viscosity sense under the strong comparison principle [62]. Proposition 5 supports applicability of the present scheme to discretization of the HJBQVI (20) in the truncated domain, the localized HJBQVI (28)- (29). Through comparing the computational results between the ODE-based method and the finite difference scheme, we examine validity of the functional form of the exact solution to the HJBQVI (20) conjectured in the previous section. Notice that the two numerical methods are based on fundamentally different concepts; the ODE-based method relies on the conjecture of the solution shape and the matching conditions at x =x, while the finite difference scheme on a direct discretization of the localized HJBQVI (28)- (29). The scheme does not a priori use functional shapes of solutions in the computational domain. In the scheme,x is detected between mid-points of each successive cells. Table 1 shows a summary of comparison results between the numerical solution with the ODE-based method ( ρ = 0.00005) and those with the finite difference scheme against different computational resolution. The auxiliary temporal integration is terminated when the absolute difference of the computed a,x,ξ between each successive steps become smaller than 10 -13 . The ODE-based method computesx = 59.392,ξ = 19.262, and c = -0.0001307, which are consistent with the assumptionsx >ξ and c < 0. Here, the value function computed with the ODE-based method is referred to as the numerical value function for the sake of brevity. In the table, the error (Er) between the numerical value function and each numerical solution computed with the finite difference scheme is the maximum error observed at vertices in 0 < x <x. The maximum error in each numerical solution was observed near x = 0 at which the value function is not differentiable in the   Table 1 demonstrates that the numerical solutions with the finite difference scheme converge toward the numerical value function as the computational resolution increases. Convergence ofx andξ is not monotone, but it seems that the results between the two methods become closer as the resolution of the finite difference scheme becomes fine. Figure 2 presents the numerical value function (0 ≤ x ≤x) and the numerical solution to the finite difference scheme with 1000 cells (0 ≤ x ≤ X max ). They agree well with each other as observed in the figure. Furthermore, the numerical solution does not contain spurious oscillations, demonstrating its satisfactory ability to handle the non-linear and non-local problem. The computed solution shape is consistent with that assumed in Proposition 6. The lower-threshold x = x(x) (x ≤ x ≤ X max ), which equalsxξ * (x), is also plotted in Fig. 2, showing that x(x) clearly depends on the population x as suggested in Sect. 3 where a difference between the intervention operators of the conventional and present model has been notified. Figure 3 shows that the convergence of the numerical solutions toward the value function is close to O(h 1/2 ), where h represents the cell size. In fact, the least-squares best-fitted curve is Er = 0.0126 × h 0.497 . This convergence speed agrees with that of the other monotone numerical schemes for HJBQVIs in impulse and related control problems [62,63]. Furthermore, we examine validity of Remark 5. We get d dx = dk 1 = 0.250 theoretically, while the finite difference approximation leads to 0.001X max = 0.251, demonstrating an agreement between the theoretical and numerical results. In summary, the results with both numerical methods are in good agreement with each other. They agree well with the theoretical results as well. In what follows, parameter dependence of the optimal policy with the finite difference scheme is analyzed. Table 2 presents the parameter dependence of the thresholdsx, x = x(x), and the amount of harvestingξ =x -x. Figures 4 through 12 graphically present the parameter dependence ofx, x = x(x), andξ on the parameters μ, σ , δ, r, R, k 0 , k 1 , a, and b, respectively. The computational results in Table 2 and Figs. 4 through 12 demonstrate monotone dependence ofx, x, andξ on the model parameters. The computational results suggest that increasing the stochasticity σ or the deterministic growth rate μ leads to a larger threshold values Table 2 Parameter dependence of the thresholdsx, x = x(x), and the amount of harvestingξ =x -x. In this table, the notations "+" and "-" mean means increase and decrease ofx, x, orξ with respect to each parameter, respectively that potentially result in managing the dynamics with higher population. The increas-ing/decreasing nature ofξ with respect to these parameters are contrasting, but Figs. 4 and 5 imply relatively smaller changes ofξ compared with the two threshold values. The dependence of the optimal policy on the discount rate δ and the parameters R and r for the utility and disutility shows that the decision-maker with a longer-term perspective maintains the population at a higher level. In addition, the decision-maker who evaluates the utility (resp., disutility) larger (resp., smaller) seems to choose larger (smaller) thresholds and leads to larger (smaller) population. The parameter dependence on the fixed k 0 and proportional costs k 1 are qualitatively the same for the upper thresholdx and the harvesting amountξ , suggesting that incurring a higher cost result in less frequent but more drastic harvesting. This is due to that executing small harvesting several times are rather costly than choosing a more drastic harvesting. The dependence between k 0 and k 1 is different for the lower threshold as suggested in Table 2, implying that incurring a higher fixed cost results in choosing a smaller x. On the other hand, incurring a higher proportional cost results in choosing a larger x. This difference can be explained by the different nature of the two costs. Increasing the fixed cost induces less frequent interventions because this cost is the same at each intervention, while higher proportional cost induces smaller amount of harvesting because this cost is positively proportional to the harvested amount of the population. These results are consistent with the model without uncertainty [44], demonstrating a qualitative linkage between the present and previous models. Finally, Figs. 11 and 12 on the parameters a and b suggest that increasing the lowerbound (resp., upper-bound) of uncertainties through increasing a (resp., b) lowers the thresholdsx and x, and further the amount of harvestingξ . This is due to that the harvested amount is the multiplicative form zξ , meaning that a largerξ is necessary to efficiently suppress the population with a possibly smaller z. A practical implication of this observation is that increasing the mean a+b 2 would be a key to effectively suppress the population at a lower level with smaller harvesting effort. Collecting empirical relationship between the harvesting effort and the realized impacts of the interventions can potentially make the range [a, b] of the uncertainties smaller or its means larger, so that the decisionmaker can better manage and predict the controlled population dynamics. However, in practice, such an effort may require costly additional field surveys, like performance tests of the interventions [64] and behavioral analysis of the bird in the habitat [65]. Our computational results suggest that such preliminary surveys should be carried out if they are evaluated as less costly than taking interventions under the uncertainties. If it is not the case, a cost term on a preliminary field survey must be incorporated into the performance index, leading to an advanced optimization problem. Importance of reducing and evaluating the uncertainties has also been suggested in a research with empirical models [66]. Our computational results support their suggestion.

Statistical indices
Here, we focus on parameter dependence of statistical indices on the controlled state dynamics. Before the analysis, we check their validity through a comparison with a numerical result by a standard Monte-Carlo combined with the classical Euler-Maruyama discretization of the optimally-controlled SDE. For the nominal parameter values, the formulae (79), (80), and (81)    The parameter dependence of the optimally-controlled dynamics on the other parameters are summarized in Table 3, suggesting monotone dependence of the statistical indices of the dynamics on the parameters. Especially the dependence of E[τ ] and Std[τ ] on each parameter is qualitatively the same, suggesting that increase or decrease of the mean and standard deviation is coherent under the optimal policy.
The quantitative analysis in what follows puts a focus on parameter dependence of the optimally-controlled dynamics on the ecological parameters μ and σ . Figures 13 and 14 show the computed E[τ ], Std[τ ], and H for different values of μ and σ , respectively. The results imply a critical importance of identifying the model parameters especially when μ is relatively small and/or σ is relatively large. Because of the monotone and convex dependence as plotted in Figs. 13 and 14. Such a situation is in practice encountered when the population dynamics is highly stochastic where the ratio μσ -2 is large. The standard deviation Std[τ ] is more affected by the parameter values than the mean E[τ ], implying that the optimal timing of the harvesting timing may be highly variable when the deterministic growth rate is small. Statistical analysis focusing on the uncertainty of the interventions is also carried out by examining parameter dependence of the statistical indices on a and b. Figures 15 and 16 demonstrate that reducing the potential uncertainties of interventions through increasing of a or decreasing b can effectively reduce both the mean and standard deviation of the harvesting frequency, and further the unit-time harvesting H. This observation combined with Figs. 11 and 12 on the parameter dependence of the management thresholds leads to a key implication that reducing the range and/or mean of the uncertainties is indispensable for cost-effective and less frequent interventions to suppress the population. This observation found through the present model emphasizes crucial importance of accurately estimating the intervention impacts in environmental and ecological management.

Conclusions
An impulse control model of a 1-D diffusion process with imperfect interventions was formulated, and the associated HJBQVI was presented with its solvability and verification results. A candidate exact solution to the HJBQVI was found and a system of nonlinear equations to determine its coefficients was derived. The system was solved numerically and the coefficients were found successfully. In addition, the (localized) HJBQVI was discretized with a finite difference scheme, demonstrating that the solutions generated by the determined coefficients and that by the scheme agreed well. Application of the scheme to the HJBQVI revealed both qualitative and quantitative parameter dependence of the optimal policy and the controlled dynamics extensively. The parameter dependence was found to be monotone for the ranges of the examined parameter values. Statistical analysis on the controlled dynamics demonstrated importance of accurately identifying the ecological parameters involved in the SDE.
Our model was possibly the simplest one, and can be extended to generalized models. Population dynamics may be subject to nonlinear [67,68] and jump disturbances [69,70]. A more general form of uncertainty in the interventions can be considered. In practical problems, different types of imperfectness would be encountered. An example is the model uncertainty in the state dynamics itself, which has conventionally been dis-cussed in the framework of the differential game between the decision-maker and nature [71,72]. The modern poorly known dynamics approach can also harmonize with the present framework [73]. The filtering problem as an identification through doing problem would be encountered if the system dynamics has a parameter whose value is not known a priori [74,75]. It is not always possible to obtain a continuous information flow of the system dynamics but discrete and costly observations may be possible [76,77]. Delayed execution of control policy is also an important element in realistic problems [78]. The present model is a building block for approaching these advanced issues. Models based on singular-like controls [79] can also be considered in the context of the intervention uncertainty discussed in this paper.

Appendix: Proofs of propositions
Proof of Proposition 1 With the null control η = η 0 , we have Combining (82), (17), and (18) yields (23). The last statement directly follows from (23). (19), considering a policy immediately harvesting the population with a fixed amount of ζ such that 0 < ζ ≤ x, which is admissible and sub-optimal, yields

Proof of Proposition 2 By
Taking the infimum of the right-hand side of (83) with respect to 0 < ζ ≤ x yields (24).
Proof of Proposition 3 Given (B t ) t≥0 and Z, the process (X t ) t≥0 starting from the initial condition is denoted as (X (x) t ) t≥0 . Set x, y such that 0 ≤ y ≤ x. Choose one η ∈ A y . Notice that we have Therefore, (Y t ) t≥0 is a classical geometric Brownian motion: with the initial condition Y 0-= xy ≥ 0. We get E Y l Proof of Proposition 4 The proof here is based on the technique of Federico et al. [80].
Secondly, we show that is a viscosity super-solution. The condition at x = 0 is obviously satisfied. Fix x 0 > 0 and set a test function ψ for viscosity super-solutions such that (x 0 ) = ψ(x 0 ) and ψ < near x = x 0 . If (x 0 ) = M (x 0 ), then we have nothing to prove and thus assume (x 0 ) < M (x 0 ). Then, there is a constant ω > 0 such that Assume (98) is not true. Then, there exists a constant ε > 0 such that Lψ + Rx Mrx m ≤ -ε at x = x 0 . Because of the smoothness of ψ and continuity of and M (By an application of Proposition 3 of Azimzadeh et al. [42] to M), we have 2δ 0 ∈ (0, x 0 ) such that Set the stopping time τ = inf{t ≥ 0 | |X tx 0 | > δ 0 } with Pr[τ > 0] = 1. The third line of (99) implies that undertaking an intervention is not optimal during (0,τ ). By (19), and thus However, the inequality is possible (103) only when τ = 0 a.s., which is a contradiction. Therefore, we have (98). The property of viscosity super-solution is then proven because x 0 > 0 is arbitrary.
Proof of Proposition 5 The proof is based on that of Theorem 3.1 of Ishii [81] where the comparison principle of a HJBQVI is discussed by constructing a strict viscositysupersolution that satisfies the inequalities of viscosity sub-solutions (in our case, (30) and (31)) with "<" instead of "≤". The boundary condition is different between the two problems, where ours is found to be simpler. The degenerate elliptic part Lψ + Rx Mrx m is not problematic because its coefficients are bounded and smooth in [0, X max ]. Therefore, we can apply the fundamental lemma (Example 1 of Crandall and Ishii [82]). The non-local part ϕ -Mϕ ≥ 0 has a somewhat different form from the conventional ones as pointed out above because of the uncertainty involved in the interventions. We can directly check that the followings, which correspond to the conditions (1)-(6) in Proposition 2.3 of Ishii [81]: let u, v be real-valued functions defined on [0, X max ], then we have Set C = rδ -1 X m max + 1 and α = min{1, k 0 } > 0. Let u and v be a viscosity sub-solution and a viscosity super-solution, respectively. For each n ∈ N, set u n = (1 -1 n )u -C n , which is straightforwardly checked to be a (strict) viscosity sub-solution to and From this observation, we see that u n is a strict viscosity sub-solution to the localized HJBQVI. Constructing strict viscosity solution is necessary to handle the non-local term. The remaining part of the proof follows Ishii [81].
Proof of Proposition 6 The proof here is based on that of Theorem 1 of Wu [83]. We must modify their proof in several parts because our model does not have jump but has an uncertainty in interventions. By Proposition 1 for any η ∈ A and x ≥ 0, we get BE e -δt X M t ≤ E e -δt φ(X t ) ≤ AE e -δt X m t , t ≥ 0, x ≥ 0.