Optimal harvesting policy for biological resources with uncertain heterogeneity for application in fisheries management

Conventional harvesting problems for natural resources often assume physiological homogeneity of the body length/weight among individuals. However, such assumptions generally are not valid in real-world problems, where heterogeneity plays an essential role in the planning of biological resource harvesting. Furthermore, it is difficult to observe heterogeneity directly from the available data. This paper presents a novel optimal control framework for the cost-efficient harvesting of biological resources for application in fisheries management. The heterogeneity is incorporated into the resource dynamics, which is the population dynamics in this case, through a probability density that can be distorted from the reality. Subsequently, the distortion, which is the model uncertainty, is penalized through a divergence, leading to a non-standard dynamic differential game wherein the Hamilton-Jacobi-Bellman-Isaacs (HJBI) equation has a unique nonlinear partial differential term. Here, the existence and uniqueness results of the HJBI equation are presented along with an explicit monotone finite difference method. Finally, the proposed optimal control is applied to a harvesting problem with recreationally, economically, and ecologically important fish species using collected field data.


Introduction 1.1 Research background
Biological resources, such as fishes, crops, and livestock, are essential for the sustainability of human lives [1][2][3][4].The aim of biological resource management is to maximize the utility while minimizing the harvesting cost [5,6] as well as the disutility that is triggered by resource exhaustion or extinction [7,8].
Thus, optimal control theory provides normative descriptions to deal with the general resource management problem (e.g., [9,10]).
Many classical optimal control problems for biological resource harvesting do not account for physiological heterogeneity, such as the timber diameter, body length, and body weight, among individuals, possibly to ensure modeling simplicity.However, such heterogenous qualities or individual differences are not negligible in real-world problems.In the simplest case, individual differences emerge as a probability distribution with respect to a physiological variable such as the body weight [11] and the distance between circuli in a scale [12].Ghwila and Willms [13] modeled physiologically structured population dynamics and the stability thereof under prescribed harvesting policies.
The emergence of heterogeneity in biological resource dynamics mathematically involves another dimension that is often infinite-dimensional, which increases the model complexity, although this aspect can be resolved analytically if the problem at hand admits a tractable structure.Examples include, but are not limited to, the graphon linear-quadratic game whose resolution is reduced to the Riccati equation [14], and the moment approximation closure [15].However, most problems that are encountered in applications do not admit such a useful mathematical structure.Furthermore, the engineering applications of the aforementioned models have not been sufficiently explored in the literature.Therefore, a mathematical model for biological resource management is required that accounts for heterogeneity while allowing the associated control problem to be solved efficiently, without resorting to complex numerical methods.
However, our focus is on the physiological heterogeneity among individuals in a habitat, as considered in some lumped fisheries management problems in a water body such as an aquaculture system or a river reach.Moreover, social heterogeneity [20] is an important concept for explaining and controlling specific population dynamics such as disease spreading.Ensemble control can also potentially handle heterogenous populations when the total number of individuals is not large [21].
The modeling of biological resource dynamics, regardless of whether they are homogenous or heterogenous, in applications faces the limitation of data availability [22,23].This data limitation creates a bottleneck in decision-making in resource management [24,25].The effects of model uncertainty can be evaluated using a Monte-Carlo simulation, as has been done in engineering research [26,27]; however, this approach is not easily implementable in control problems that are subject to model uncertainty even under a homogenous case, because it is inefficient for solving the control problem numerically a vast number of times.
Multiplier robust control [28] is a differential game formalism of optimal control problems that is subject to model uncertainty, and has mainly been applied to stochastic control problems in which the underlying probability measure is distorted owing to insufficient knowledge of the decision-makers [29][30][31].Barnett et al. [32] parameterized the model uncertainty using prior probability weights, which can be considered as a form of uncertain heterogeneity in harmony with the control problem.A theoretical as well as practical advantage of the multiplier robust approach is that the model uncertainty can be represented by suitably perturbing the objective function and system dynamics of the targeted control problem within the context of the major optimality principle, such as the dynamic programming and maximum principle.Nevertheless, to the best of the authors' knowledge, this approach has not been incorporated into biological resource management with physiological heterogeneity.

Research objective and contribution
The objectives of this study are twofold: the formulation of a finite-horizon optimal control problem for biological resource physiological heterogeneity, and the application thereof to a real case.To achieve the first objective, we consider a harvesting problem of a population in a habitat in which the individuals are physiologically heterogenous, with an emphasis on fisheries management.The proposed model couples the population dynamics coupled with a continuum of growth curves representing the physiological heterogeneity, leading to a unique hybrid system that has not been considered.More specifically, the heterogeneity is represented using a probabilistic parameterization of the logistic-type growth model, where the growth rate and/or maximum body weight are randomly distributed among the individuals (e.g., [33,34]).The uncertainty here concerns the heterogeneity: the probability distribution of the physiological heterogeneity.Therefore, in the proposed model there is no uncertainty if there is no heterogeneity, which simplifies the theoretical formulation but in a nontrivial way as discussed below.
Note that the approaches in the literature [20,21] are based on microscopic individual-based models, whereas ours is an aggregated one that deals with the population as a whole.The two approaches are therefore fundamentally different with each other.The approach based on the incomplete information [48] assumes that a model parameter is possibly better estimated as the time is elapsed due to a statistical filter, while ours assume that this mechanism is absent.Instead, we assume in this study that the decision-maker can manage the fishery resource considering model uncertainties through the uncertainty aversion term in the objective function.
The novelty of this study lies in the combination of multiplier robust control with this randomized formalism in a control problem to account for possible misspecifications of the probability measure when generating the randomness.The degree of randomness is represented by the Kullback-Leibler divergence, which is a widely used statistical divergence (e.g., [35]).We demonstrate that solving this control problem formally can be reduced to determining a proper solution to the Hamilton-Jacobi-Bellman-Isaacs (HJBI) equation.In this equation, the consideration of the model uncertainty emerges as a unique nonlinear term that potentially raises issues in both theory and computation.We demonstrate that, despite its complicated form, the HJBI equation admits a Lipschitz continuous Hamiltonian, to which the comparison argument to prove the unique solvability in a viscosity sense (e.g., [36,37]) and a convergent finite difference discretization (e.g., [38,39]) can be applied with modifications.
Our HJBI equation is related to problems that arise in machine learning, especially in reinforcement learning, where the decision-maker should select randomized controls (e.g., [40]).
Exploratory control [41], maximum entropy optimal control [42], and feedback relaxed control [43] lead to specific nonlinear and nonlocal optimality equations corresponding to the randomized control.Our HJBI equation differs because the control problem is aimed at maximizing some utility under uncertain randomness in the population dynamics; hence, the system dynamics, rather than the control variable, is randomized.In this sense, the linear-quadratic regulator of the uncertain system matrix [44] is closer to ours, although this methodology depends on the maximum principle rather than dynamic programming.
Nevertheless, our and the aforementioned optimality equations share some similarities, such as the viscosity property, which facilitate the understanding of all their formalisms.It should be noted that, to the best of our knowledge, our control problem has not yet been studied, even with the absence of model uncertainty.Furthermore, our control problem is not expected to satisfy the so-called Isaacs condition (i.e., the order of maximization and minimization can be changed in a differential game) owing to the lack of convexity of the Hamiltonian.
We also present a simple finite difference method for discretizing the HJBI equation to provide demonstrative examples using real data.Finally, we apply the proposed control problem to a harvesting problem of the inland fishery resource Ayu Plecoglossus altivelis altivelis, which is one of the most recreationally, economically, and ecologically important fish species in Japan [45,46].Although the authors have studied harvesting problems of this fish, their models have assumed only homogeneous populations [47,48].The proposed mathematical framework that accounts for physiological heterogeneity extends these methodologies.The population dynamics of the fish is estimated from the unique collected data of the body weights of the individuals since 2020.The computational results demonstrate the optimal harvesting policy of P. altivelis depending on the control objective as well as the potential model uncertainty.We demonstrate that the uncertainty acting on the heterogeneity gives the worst-case underestimation of the probability density of the body weights of the individual fishes in the population according to the uncertainty aversion of the decision-maker.Different probability densities having different skewness are examined to see their impacts on the harvesting policy and the worst-case growth curve.This work covers the formulation, analysis, and application of a novel optimal control problem.
The remainder of this paper is organized as follows: The mathematical framework that is employed in this study is introduced in Section 2. The mathematical analysis with a focus on the HJBI equation is presented in Section 3, along with the finite difference method and its theoretical analysis results.Section 4 outlines the model application to real data, focusing on the harvesting problem of P. altivelis.A summary and future perspectives of our research are presented in Section 5. Proofs of the propositions and lemmas are provided in Appendix A. Few additional computational results are presented in Appendix B.

Mathematical model
We consider a harvesting problem of a biological resource in a closed habitat in a finite horizon during which spawning (i.e., a population increase) does not occur, which is effectively considered as the terminal utility in the objective function.

Population dynamics
This study focuses on fisheries management, and the utility is evaluated through the biomass of aggregated individuals.The resource dynamics in this paper contains two continuous-time variables: the population (i.e., total number of individuals) and a probabilistically distributed body weight.The population ( ) 0 t t N  evolves according to the differential equation ( ) where t is the time,

 )  )
: 0, 0, R + → + is the mortality rate that is bounded and positive, increasing, and Lipschitz continuous: c  represents the harvesting rate that belongs to the admissible set ( )   0 is measurable and 0 ( 0) with a constant 0 c  , and , tn  is the stopping time such that   , inf 0 and We set 0 The population is approximated to be a real variable rather than an integer, which can be reasonable if it is significantly larger than ( ) , as considered in our application.The population-dependent mortality rate covers both constant and state-dependent cases, with the latter able to consider the density-dependent mortality if necessary (e.g., [49,50]).The constant c that serves as the upper bound of the control variable c represents a technological constraint.It also ensures the unique existence of the solution to (1), as indicated below.The set is non-empty as the constant control 0 c  belongs to it.The differential equation ( 1) is solved for 0, 0 n t   as follows: which implicitly defines the unique continuous solution to (1) because the right-hand side of (4) is well defined and R is bounded and Lipschitz continuous (e.g., Theorem 5.5 in Chapter 3, Section 5 of [51]; see also conditions (A0), (A1), and (A3) of this work).The population N is set to 0 after the hitting time 0,n  .In this manner, t N is a continuous decreasing function in  ) 0, + , where the decreasing property is immediate from the non-positivity of the right-hand side of (1).

Growth dynamics
Another component of the resource dynamics is the continuum of growth equations, which is simply referred to as the growth dynamics.We assume the widely-used logistic model (e.g., [34]): where

 
0,1 .Furthermore, the initial condition 0 x  is sufficiently small such that ( ) , the right-hand side of (5) exists as a uniformly continuous function of time 0 t  .We assume that there exists a probability density function ( ) pu such that the averaged weight in the population is defined as Subsequently, the total biomass of the population at time 0 t  is given by tt NX .Similarly, the unit-time harvested biomass is given by tt cX .
Remark 1 Employing the logistic growth curve is for analytical simplicity.Other growth curves, such as the Gompertz and generalized logistic curves (e.g., [34]), can be used instead of (5).Theoretically, ( ) t Xu is required to be uniformly continuous and strictly increasing with respect to 0 t  as well as strictly positive and strictly bounded from above for all 0 t  and  + with constants 12 ,0 XX .A key point in using the S-shaped curve is its profile having an inflection point, with which the model can be reasonably fitted against data as demonstrated in this study.Another key point, which will be critical in applications, is its boundedness such that the growth curve eventually saturates as the time is elapsed.Growth curves that are simpler than the S-shaped one, such as an exponentially increasing curve, will be easier to identify and operate, whereas its unboundedness will overestimate the growth curve and the biomass.Similarly, one will use a variety of distributions p .In Section 4, we argue that the simplest distribution will be the uniform, one and beta distributions can be used for unimodal alternatives; multi-model ones can also be used if necessary, but is not the case for the application to P. altivelis.

Model uncertainty
The model uncertainty is represented by a Radon-Nikodym derivative as a positive and measurable continuous-time field , which is referred to as the uncertainty field as well in this paper, such that the integrability condition is satisfied: and the Kullback-Leibler divergence (e.g., [35]) exists: The uncertainty field can be considered as an uncertainty that is induced by incomplete knowledge of the growth dynamics.For later use, we set the admissible set of  as follows: ( ) ( ) ( )   0 ( 0) is measurable, non-negative, and satisfies ( 7)-( 8) Note that the function ln 1 y y y −+ ( 0 y  ) is strictly convex, non-negative, and globally minimized at 1 y = with the minimum value 0. Here, we understand 0 ln 0 0 = .The set is non-empty as the constant field 1   belongs to it.
According to the formulation above, the heterogeneity is modulated by the uncertainty.This implies that the uncertainty not only changes the profile of p but also its moments, which further affect the biomass.This means that we focus on the uncertainty in the growth curve that eventually affects the biomass of the target fish and hence its resource management.

Control objective
We consider a finite-horizon control problem during a fixed period    , , , : d d where In the first term of (10), 0   is the weighting factor of the harvesting utility and ( )   is the power index of the utility.The concavity of this term with respect to s c combined with the linear harvesting cost results in a balanced optimal harvesting rate, as demonstrated later.Moreover, it contains a tractable case 1/ 2

 =
, in which the Hamiltonian of the HJBI equation is determined analytically.In the second term, 0   is the weighting factor of the harvesting cost, where the cost is proportional to the harvesting rate, which assumes that the labor that is required to harvest the fish increases linearly with respect to the population to be harvested because harvesting a larger population requires a larger effort.
This linearity is another key to obtain the optimal control in a closed form.In the third term, = is a bounded, increasing, and Lipschitz continuous function that represents the terminal utility, which would effectively suppress the resource extinction [47].In Section 4, we computationally demonstrate that the use of a simple form of h realizes a positive terminal population that potentially contributes to the spawning, and hence, maintains the life cycle of the biological resource.Finally, in the last term, 0   is a weighting factor that serves as an uncertainty aversion parameter such that a larger (resp., smaller)  more weakly (resp., more strongly) allows for the existence of the model uncertainty, and hence, supposes a model closer to (resp., farther from) the benchmark model with 1 11) is a distorted counterpart of ( 6) according to the uncertainty field  .
We set the domain the objective (10), we define the value function : This value function is understood as a pessimistic maximum value of the control objective (10) because the right-hand side is a minimum of a maximum.According to (10),  should satisfy the following boundary and terminal conditions: ( ) and respectively.The boundary and terminal conditions are compatible at ( ) ( ) In the following, we focus on the aforementioned case 1/ 2 as it is the most tractable case among ( ) , where the HJBI equation becomes the most tractable while its characteristics will not be critically lost.
The value function  is non-negative and bounded.Indeed, its non-negativity, which is a lower bound, follows by considering J with constant controls in combination with the non-negativity of , h : The other bound follows from (

Formulation
The dynamic programming principle heuristically leads to the HJBI equation.That is, the optimality equation that governs the value function ( ) is expressed in domain  as follows: ( This HJBI equation is subject to boundary condition (13) and terminal condition (14).We can more compactly rewrite (17) by explicitly calculating the "inf sup" part, while it should also be noted that the so-called Isaacs condition (i.e., the order of "inf" and "sup" can be exchanged in (17)) is not expected, as discussed in the following sub-section.

Isaacs condition and heuristic optimal controls
The discussion in this subsection applies to the case 1/ 2   as well.The Isaacs condition for (17) is not expected; that is, the equality for any ( ) is not expected in general.The primary reason for the failure of the Isaacs condition is the left-hand side of (18).With respect to  , the first term of the functional is concave, whereas the second term is convex, which suggests that the terms inside "  " of the left-hand side of ( 18) are not necessarily convex particularly for small 0   .This implies that the inner maximization on the left-hand side of ( 18) is not always achieved.In contrast, for each fixed ( ) ( ) ( ) is concave and is maximized by In the next step, it is necessary to minimize which turns out to be a strictly convex optimization problem under a suitable assumption.For example, this condition is satisfied with a sufficiently large c (the technological limitation does not restrict the optimal harvesting policy) and an increasing property of ( ) with respect to 0 n  (the value function, namely the net utility, is larger for a larger population).We present the following Assumption 1 ( c is sufficiently large) and Lemma 1 in preparation for the more detailed discussion in Section 3. Assumption 1 is assumed throughout this paper to simplify the problem at hand.

Lemma 1
For 0 tT  , 0 z  , and ( ) ( ) ( ) From Lemma 1, for each 0 z  , we obtain The second line of ( 25) is a convex optimization problem (e.g., Section 2 of [52]), and the infimum is achieved by the minimizer Hence, the right-hand side of ( 25) becomes Subsequently, we obtain the candidate for the optimal control ( ) and In this context, the determination of the optimal controls is reduced to the resolution of the HJBI equation.
Based on (28), the pessimistic estimate of the averaged body weight is obtained as follows: Formulations ( 28)- (30) imply that the partial derivative n   plays a crucial role in determining the optimal harvesting policy as well as the model uncertainty.Phenomenologically, a larger n   leads to a smaller uncertainty aversion, and hence, a more optimistic result.

Remark 2
The loss of the Isaacs condition discussed above is considered due to the concavity of its first term in (19).Using a convex alternative will be theoretically possible, but a problem in our context is that the optimal control may not be obtained explicitly as in (21), which makes the model be complicated.It implies that if one uses a convex alternative utility, then the Isaacs condition can be satisfied while the model complexity would significantly increase.Moreover, in such a case, the computational algorithm used in this study will have to be significantly customized, which is beyond its scope.By contrast, we consider that using an alternative distribution p (other than the uniform and beta ones examined later) does not affect the Isaacs condition because this condition is with respect to  .

Remark 3
Eq. ( 26) demonstrates that, for any probability density p of the heterogeneity, the worst-case distortion of the heterogeneity is proportional to the exponential given in Eq. ( 5).As , this exponential is monotone and more specifically decreasing in u .The ratio of individuals possibly having smaller (resp., larger) body weight is therefore underestimated in the parametric way. Figure 7 in Section 4 and Figure B2 in Appendix B visualize quantitative magnitudes of the worst-case distortion on the growth curve.

Value function
We analyze the monotonicity and continuity properties of the value function  .We first present a proposition stating that  is increasing in the second argument.This implies that we should focus on solutions to the HJBI equation complying with 0 n    in some sense.

Proposition 1
For each Remark 4 One may proceed to a verification argument of viscosity solutions, which potentially yields an existence proof of the viscosity solutions to the associated HJBI equation.We do not take this direction in this study owing of the possible lack of the Isaacs condition in our HJBI equation, which makes the problem more difficult than one that complies with this condition.Instead, we construct a continuous viscosity solution to our "modified" HJBI equation with a truncated Hamiltonian, as presented later, based on the finite difference method.Owing to Proposition 1, this truncation can be removed and the unique viscosity solution to the modified equation is also a viscosity solution to the HJBI equation.Later, we numerically imply that the value function is continuous but not Lipschitz continuous.Possibility of formulating an analogous problem satisfying the Isaac's condition is an open issue.

HJBI equation
According to Section 2.4.2, formally, if , the HJBI equation ( 17) is rewritten as We wish to discuss the uniqueness of the (viscosity) solutions to the HJBI equation, which is critically dependent on the regularity of the Hamiltonian H .We present the following Lemma 2, which states that H is Lipschitz continuous with respect to 0 z  and uniformly continuous with respect to   0, tT  .

Lemma 2
The Hamiltonian H in ( 27) is uniformly continuous with respect to for each 0 z  .Furthermore, it satisfies the Lipschitz continuity , is uniformly continuous with respect to 0 tT  for each 0 z  .Moreover, H is decreasing in the second argument.
Remark 5 Lemma 2 states that our Hamiltonian admits certain regularity properties despite its apparent complexity.Another key finding is that the existence of the harvesting cost, namely the positivity of  , is essential for its regularity.Indeed, the proof of Lemma 2 suggests that the Lipschitz constant of ( 33) is sharp and becomes arbitrary large as  approaches 0 .This is in contrast to [53], in which a costless formulation was considered, and a simpler and more tractable objective function was used.
The uniqueness is proven through the comparison argument owing to Lemma 2. In combination with the increasing property of the value function  (Proposition 1), if there exists a unique viscosity solution to the modified HJBI equation ( ) that is subject to the same boundary and terminal conditions as the original one, with this solution also solves the original HJBI equation (32)  ( z  ) has the Lipschitz constant 1.The modified Hamiltonian is defined over z  and covers a wider range than the original one (27).Suppose that the modified equation ( 34) admits a unique viscosity solution.Under this assumption, if a viscosity solution to the HJBI equation ( 32) that satisfies in some sense can be found, this solution will be the unique viscosity solution to the modified equation (34) as well.This implies that the solution to the modified equation (34 which suggests that the modified equation ( 34) serves as the optimality equation of our control problem as well.This type of use of a modified equation was employed in [54] to deal with a viscosity solution with a certain monotonicity property.
We define the viscosity solutions to the HJBI equation (32).The collection of all continuous functions that are defined in a domain D is denoted as

( )
CD .Similarly, the collection of all upper-semicontinuous (resp., lower-semicontinuous) functions that are defined in a domain D is denoted as ).The linear growth speed of the viscosity solutions can be justified owing to the boundedness (16).

Let
( )  , which grows at most linearly with respect to the second argument.We refer to  as a viscosity super-solution of (34) if, for any ( ) ( ) ( ) , and  , which grows at most linearly with respect to the second argument.We refer to  as a viscosity sub-solution of (34) if, for any ( ) , and

(c) Viscosity solution
A continuous function :   → is a viscosity solution if it is a viscosity super-solution in the sense of (a) and a viscosity sub-solution in the sense of (b).
The viscosity solutions to the original HJBI equation are also defined.

Definition 2
The viscosity super-solutions, viscosity sub-solutions, and viscosity solutions to the HJBI equation (32) are defined by following Definition 1, where H is formally replaced with H .
At this point, we state the uniqueness result of the modified HJBI equation.

Proposition 2
The modified HJBI equation ( 34) admits at most one viscosity solution in the sense of Definition 1.
We close this sub-section by analyzing the no-uncertainty limit 0  →+ of the HJBI equation.For this purpose, for ( )   , 0, t z T  , we set A straightforward calculation shows that 0 locally uniformly in ( )   , 0, t z T  , which leads to the following stability result (e.g., Lemma 3.2 of [41]; Section 6 of [36]).In Proposition 3, the viscosity solutions to the HJBI equation in which H is formally replaced with 0 H are defined analogously to those of Definition 1.This proposition states that the proposed HJBI equation that accounts for the model uncertainty is consistent with that without the uncertainty.

Proposition 3
is a sequence of viscosity solutions in the sense of Definition 1 for each fixed 0 ) is continuous on  , it is a viscosity super-solution (resp., sub-solution) to the equation ( ) subject to the boundary and terminal conditions (13) and (14).

Formulation
We present a monotone finite difference method that is fully explicit in time for the modified HJBI equation (34).The domain  is first discretized into a computational grid where ,0 tn    are the grid sizes in time t and population n , respectively.The degree of freedom of the grid is controlled by , IJ .We assume that I t T = as well as J n M = with a constant 0 M  .In practice, the parameter M can be selected to be the maximum population at the initial time 0 t = , which can be estimated from an ecological survey.
The discretized solution  to the HJBI equation (32) at the grid point ( )  .Our finite difference method directly specifies the boundary condition as and the terminal condition as ( ) Subsequently, the HJBI equation itself is discretized at each , P ij ( 0 ( ) ,0 which can be rearranged as (recall that we are dealing with an equation to be solved backward in time) , , , In our numerical computation in Section 4, the integration with respect to u that is involved in H is carried out by a midpoint rule with uniformly distributed at 150 points ( ( ) ), which was preliminary found to be sufficient for our application.Using ( 41), (42), and (44), we can obtain the numerical solution   ) from iI = to 0 i = without resorting to solving any matrix inversion.However, it is necessary to take a sufficiently small t  for stability (See, next subsection).
Remark 6 Using (44), we do not require information outside the computational domain; that is, we do not need to specify any artificial boundary conditions along jJ = .This is owing to the mathematical structure of the discretized and original HJBI equations that the information propagates from 0 n = to larger n values.

Analysis of numerical method
We first prove the following stability result, and state that the non-negativity of the numerical solution is satisfied for a sufficiently small t  .The proposed finite difference method can be stabilized, whereas the computational cost increases linearly with respect to 1/ t  owing to the fully explicit nature.

Proposition 4
The non-negativity Furthermore, we obtain an upper bound of the numerical solution that is uniform in time.

Proposition 5
From (45), it follows that Propositions 4 and 5 demonstrate that the numerical solutions generated by the proposed finite difference method are uniformly stable for a sufficiently small time increment that satisfies (45).Furthermore, an increasing property of the numerical solutions with respect to j is satisfied in the HJBI equations as a byproduct of the monotonicity.The proof of the proposition uses the argument of Lemma 5.2, 3(a) of [55].This is a discrete counterpart of Proposition 1.

Proposition 6
Under ( 45), the right-hand side of (44) increases with respect to both Remark 7 According to Propositions 4-6, scheme ( 44) is essentially the same as ( ) , which is the naïve finite difference discretization of the original HJBI equation (32).

Remark 8
One may expect that the numerical solutions that are generated by the finite difference method converge to a viscosity solution in the sense of Definition 1 locally uniformly in ( ) ( ) is selected for a sufficiently small constant 0   that at least complies with (45).This follows from the classical convergence argument [38] because the scheme is monotone, stable, and clearly consistent, the Hamiltonian H is Lipschitz continuous, and both the boundary and terminal conditions are Lipschitz continuous as well as bounded.See also Theorem 3.15 and Section 4.1 of [56].
Therefore, the numerical solutions yield the existence of a viscosity solution; hence, they converge to the unique viscosity solution by Proposition 2.More specifically, we prepare the monotone bilinear interpolation   of the numerical solutions: for each ( ) ( ) ( )

Remark 9
Relating to Remark 8, we find that the convergence that preserves the increasing property in n appears to be far trivial.Such a proof can be found in the proof of Theorem 5.5 of [55].However, their proof fails in our case owing to the lack of no increasing or decreasing property with respect to t .Indeed, the numerical solutions that are obtained in Section 4 suggest that the solution to the (modified) HJBI equation is non-monotone with respect to t .

Study site
For the model application, we used growth data of P. altivelis that were collected in the Hii River, which is one of the largest rivers in the San-in area, Shimane Prefecture, Japan, with the help of the Hii River Fishery Cooperative (HRFC) which authorizes inland fishery resources in the mid-to-up-stream reaches of this river.The authors had been communicating with the HRFC since 2015, and previously investigated the environment and fisheries of the Hii River (e.g., see [57,58] and the references therein).
The fish P. altivelis is a major inland fishery resource that contributes to the regional environment, ecosystem, and culture and recreation.Tomozuri (fishing with decoys [59]) and Toami (casting a net [60]) are major fishing methods for catching the fish.
P. altivelis has a one-year life history such that they spawn in a river during autumn, and the spawn fishes flow down along the river to the downstream sea or a reservoir, overwinter in the sea or a reservoir, and migrate towards the river in the spring to mature during the summer (e.g., [61]).The harvesting period of P. altivelis in the Hii River usually starts on July 1 ( 61 t = (days) below) and ends in October to November ( 181 T = (days) below), where May 1 was selected as the reference point of 0 t = (days) following previous studies (e.g., [57]), which corresponds to the growing season of the fish.Therefore, the computational time domain is ( ) 61,181 (days) with the length of 120 (days).The population dynamics of P. altivelis in the Hii River has not been clarified to date.In particular, no estimate has been made for the overall river system.However, it has recently been estimated that ( ) juveniles of P. altivelis were found in the upstream main branch of the Sakura-Orochi reservoir, which creates the largest dam, namely the Obara Dam, in the Hii River [62].Thus, we focused on the application of the proposed model to the river based on the assumption that the maximum population at the beginning of July was ( )

Computational conditions
Sampling surveys of the body weights of P. altivelis individuals are carried out by a union member of the HRFC from Summer to Autumn each year, which we used to estimate the model parameters of X .We used the most recent data that were collected during 2021 and 2022, as illustrated in Figure 1.Note that different individuals were caught at different observation times.
The fitting method of the uncertain logistic model ( 5) is based on the maximum likelihood method of Dorini et al. [33], whose conclusion demonstrated that the uniform distribution should be assumed with a minimum data requirement provided that the upper bound ( ) and lower bound ( ) are given.Within this framework, we assumed that only K is uncertain, leading to the parameters that were estimated in the uncertain logistic model being x , r , and , KK.For each year, we attempted to find the best values of ( ) , , , x r K K to minimize the common least-squares error loss function between the observed model ( obs,i X ) and fitted model ( theor,i X ): ( ) where obs N is the total number of observations (40 in 2021 and 59 in 2022), obs,i X is the observed weight at the i th observation time, and theor,i X is the theoretical mean body weight at the i th observation time.The minimizer of the loss function (48) was determined using a trial-and-error approach such that we selected the best quadruple ( ) , , , x r K K to minimize (48) among the parameter values of 0.020-0.050r = (1/day): increments of 0.001 (1/day), 5-15 x = (g): increments of 1 (g), 1-301 K = (g): increments of 1 (g), and 1-301 K = (g): increments of 1 (g) with the constraint 1 KK + (g). Figure 2 compares the collected data and the fitted model of the body weights of individuals for each year.Table 1 summarizes the fitted parameter values for each year.The fitted results tracked the observed data reasonably.The parameter values in Table 1 suggest that between 2021 and 2022, the fitted model of P. altivelis in 2021 had larger individual differences than that of 2022, with a higher growth speed.Furthermore, the individuals were estimated to be smaller than those in 2022 near May 1; that is, around the season of upstream migration.The reason for these differences can be attributed to different hydrological and meteorological conditions as well as the associated ecosystem dynamics, which are difficult to clarify based only on the available data and are beyond the scope of this paper.Nevertheless, the fitted results imply that the uncertain logistic model could effectively capture the growth dynamics of the fish for different years.
The other parameter values were similar for the two years.The population dynamics (1) was normalized to 4 10 individuals considering the discussion in the previous sub-section.Therefore, in the computation, the individuals were counted in units of 4 10 .Under this normalization, the maximum domain size M in the n direction was set as 10, corresponding to (with the sustainability concern); these values had been preliminarily explored so that realistic sample paths of the controlled resource dynamics could be obtained as well as we can visually demonstrate their influences on the optimal harvesting policy as well as the worst-case uncertainty.Particularly, we have heuristically chosen the parameter value of  so that the uncertainty affects the computational results while it will not lead to a too small estimate of the body weight (e.g., see Figure 7).This parameter is user-specific and quantifies how much the decision-maker concerns the uncertainty.Its estimation will need another study about decision-making mechanisms inside humans, which is beyond the scope of this study.The computational resolutions were 120 000 I = and 1 000 J = , which were found to be sufficient for the computation below. for the years 2021 and 2022, respectively, with the sustainability concern.The computed value functions appeared to be continuous and smooth inside the domain in both years, regardless of the existence of the sustainability concern, thereby supporting the continuity assumption of the mathematical analysis.The computed  in 2022 was larger than that in 2021 in most parts of the computational domain, suggesting that the smaller individual differences led to a larger net utility under the assumed condition.
Figures 5 and 6 depict the optimal harvesting policy and several controlled paths of the populations corresponding to Figures 3 and 4, respectively.The controlled paths were computed by numerically integrating the population dynamics (1) backwards in time based on the optimal harvesting policy that was obtained via the HJBI equation.This approach is justified in our case because the population dynamics does not contain stochastic noise.The set of terminal values of the backward paths is the same in each panel of Figures 5 and 6.A comparison of Figures 5 and 6 reveals that accounting for the sustainability of the fish population critically affected the controlled population dynamics.According to Figure 5, for each year, the backward computed paths indicate that, given a terminal value, the population dynamics to achieve a positive terminal value required initial populations that were significantly larger than 5 10 individuals.In contrast, Figure 6 suggests that populations of the order of ( ) archived the positive terminal values with the same order.Indeed, the optimal harvesting rates exhibited one order of difference between the cases with and without the sustainability concern.Although the required total number of populations to sustain the life cycles of the fish is dependent on the specific relationship between the body size and total number of eggs, as discussed for other fish species (e.g., [64,65]), the proposed mathematical model suggests that considering the terminal utility is key to the sustainable fisheries management of P. altivelis in the Hii River.
We also analyzed the growth dynamics under the uncertainty, which is a key element in our mathematical model.Figure 7 presents paths of the distorted mean body weights (6) that were subject to different levels of uncertainty aversion for 2021 and 2022 with the sustainability concern.Note that if a classical deterministic logistic model is used, we formally only obtain XX = , and hence, there will be no uncertainty effects in the control problem.In this study, the distorted X for each year was almost the same among the various sample paths with different terminal conditions, and the difference was at most several grams, so we only present the results of the growth dynamics corresponding to the terminal population of 4 10 individuals (i.e., 1 n = in the computation).The effects of the uncertainty aversion in terms of the weighting coefficient  between the two years is clear, which demonstrates that critically different mean growth curves were generated for the same value of  , particularly for the smaller  representing the larger uncertainty aversion.The larger potential individual differences in the model for 2021 were more susceptible to the model uncertainty, and an overly pessimistic prediction yielded an unrealistically small mean body weight (red curve in Figure 7(a)).One means of preventing model uncertainty to avoid such a pessimistic result is to collect the information at the beginning of the growth period, particularly before the harvesting season, which is the time 0 to 60 in our setting.Our model suggests that although such a survey would require high monetary and labor costs because the growth conditions are different among different years as demonstrated in this paper, its outcome can persist during the harvesting period to provide a less pessimistic prediction of the mean growth of the individual fishes.The computational also results suggest that the mean growth under the uncertainty approached that without the uncertainty, which is consistent with the theoretical results for the value function in Proposition 3.
We finally analyze the case with a non-uniform p as our mathematical framework theoretically covers such a case.We focus on the beta distribution with parameters ,1 ab ; that is, ( ) ( ) with a normalization constant 0 p .The beta distribution is right-skewed, non-skewed, and left-skewed for ab  , ab = , and ab  , respectively (Figure 8).The data in the year 2021 is used here because of the larger uncertainty than 2022, with which the influences if the skewness can be more clearly visible.The right-skewed case seems to be the most qualitatively reasonable because the observed body weight data of P. altivelis in the past years were right-skewed (e.g., [11,48,57]).We consider all the three cases by specifying ( ) ( ) , The parameter values other than those in p are the same with the computational cases presented above, and we set the uncertainty aversion as 0.01

 =
(1/day). Figure 9 shows the computed optimal harvesting policy and several controlled paths of populations in 2021 with sustainability concern for the right-skewed, non-skewed, and left-skewed cases.The skewness of the body weight distribution affects the optimal harvesting policy as well as the controlled population dynamics so that the smaller harvesting rate is optimal for a more right-skewed case.This is due to the modelling assumption that harvesting a larger individual is more profitable, meaning that too fast resource exploitation is less efficient for such a case.Figure 10 presents computed paths of distorted mean body weights with sustainability concern, demonstrating that, even under the uncertainty, the left-skewed case predicts the larger mean body weight than that without uncertainty that is due to the distributional shape of p concentrating more on the right-half of 01 u  as shown in Figure 8. Accordingly, the non-skewed and right-skewed cases yield more pessimistic results than the left-skewed case.A comparison between Figure 7(a) and Figure 10 suggests that using the beta distribution of the non-skewed case yields a more pessimistic result than the uniform case although the beta and uniform distributions share the same mean.The computational results thus suggest that the mean growth with the beta distribution, provided it is non-or right-skewed, is more pessimistic than the uniform distribution.Although the investigations of more realistic p remains a future topic, its influences were qualitatively suggested by our numerical analysis.

Conclusions
This study proposed an optimal control problem for harvesting biological resources with physiological heterogeneity under model uncertainty.The Radon-Nikodym derivative was effectively employed to evaluate the uncertainty of the heterogeneity, and the resulting HJBI equation was successfully obtained in a tractable form.The mathematical analysis results characterized the shape and viscosity properties of the HJBI equation.A fully explicit and monotone finite difference method was presented and implemented in a real problem using the collected data.The optimal harvesting policy, which is dependent on the control objective and degree of uncertainty aversion, was successfully computed.
The proposed model can be applied to any other biological resources having physiological heterogeneity, such as a fish species other than P. altivelis [66].We consider that each fish species can be characterized by the growth curve and population decrease/increase because they are used to estimate the biomass.Once the dynamics of the two quantities of the target fish have been identified, the proposed control setting can be applied to address its resource management problem.The HJBI equation for the management problem will then be obtained.Adding a noise process such as Brownian motions to the population dynamics is theoretically possible, while the difficulty lies in the estimation of the volatility (i.e., noise intensity) from field data as it would need time series of growth paths of individuals, which are difficult to track in practice.Nevertheless, one may be able to estimate the volatility based on the assumption that the growth rate r should be larger than growth declines due to stochastic fluctuations in some sense.For example, one may assume the increasing property of the second moment of the body weight (e.g., fishes will continue to grow in some sense).In our context, if we assume the no-uncertainty case where the body weight is a time-dependent quantity Several topics will be addressed in our future works.Resource harvesting problems in an infinite horizon [67] and that under an economic equilibrium [68] can be investigated based on the proposed approach by setting the terminal time to be sufficiently large.The size-dependent mortality [69] was not considered in this study owing to the increased complexity of the HJBI equation, such that the built-in optimization problem, namely the inf part, possibly becomes intractable and/or ill posed; this will be studied in the future.A more sophisticated numerical method, such as the fitted finite volume method [70], will be necessary to consider the noise-driven resource dynamics.The extension of the proposed mathematical approach to a multi-site problem [71] will also be an interesting future direction, wherein the efficient numerical computation of the associated optimality equation could be a potential challenge.Population dynamics with a delay sometimes arises in applications, and this can also be handled using dynamic programming, although the problem dimension becomes infinite [72] in such cases.Agent-based models as alternative representations of the heterogeneity are also interesting subjects that should consider high-dimensional dynamical systems.For example, it will be interesting to extend the proposed model so that the social interactions [20,21] can be accounted for.A model reduction will be necessary to analyze such advanced problems to make the problem computationally feasible.Incorporating a statistical filter [48] to improve the uncertainty estimate will be another option to extend the proposed model.We focused on a problem in a natural environment, while the proposed approach also applies to that in an artificial system such as an aquacultural system [73].
, which is impossible unless   = + : a contradiction.Therefore, we must have (51).As a byproduct, the continuity of ( 52) indicates that we must obtain The left-hand side of ( 60) is not larger than In the same manner, we have which is the desired inequality (31).

Proof of Lemma 2
The uniform continuity of H follows from the uniform continuity of , which is further combined with the lower bound which proves (33).The uniform continuity of H z   follows from the last line of (66).The decreasing property of H for 0 z  is immediate from its functional form.

Proof of Proposition 2
This proof is an adaptation of the classical method of doubling variables using Lemma 2. Indeed, we require the following inequality for any viscosity super-solution  and any sub-solution  : The inequality is clearly satisfied along the boundaries ( )   0, 0 T  and    ) 0, T  + .Subsequently, we prove the inequality (68) at each point in  using the auxiliary function (e.g., Proof of Theorem 7.5 of [54]).

Proof of Proposition 3
The proof is owing to the fact that 0 HH  → as 0  →+ locally uniformly in   The proposition then follows from the definition of viscosity solutions (e.g., Lemma 3.2 of [41]; Section 6 of [36]).

Proof of Proposition 4
We use a method of induction.From (41), we obtain which is satisfied if ( ) ( ) Consequently, we arrive at the desired result (45) by induction.

Proof of Proposition 5
We again use a method of induction.Using (41), we obtain ( ) ( ) Assume that we have    We further obtain We obtain the desired result directly from (76) and (77); namely, we should set  larger than the right-hand side of (76).

Proof of Proposition 6
The first statement is immediately verified by the partial differentiations ( where we used

satisfies the bounds in Propositions 4 and 5 .
According to the Ascoli-Arzelà theorem and Propositions 4 and 5, along with the selection of the discretization parameters ( tn   =  ), there exists a subsequence of

5 10( 1 /
individuals.The mortality R was set to 0.01 (1/day) for simplicity considering previous identification results[63].The parameters in the objective function were set to 0day), and ( ) 0 hn = (no sustainability concern) or ( )  1.5min , h n n M =

Figure 1 .
Figure 1.Collected data of body weights of individuals in 2021 (red) and 2022 (blue).

Figure 2 .
Figure 2. Comparison of collected data and fitted model of body weights of individuals in (a) 2021 and (b) 2022.In each plot, the circles denote observations, the black curve indicates the theoretical mean, the red curve denotes the theoretical mean ＋ standard deviation, and the blue curve indicates the theoretical mean − standard deviation.

4. 3
Figures 3(a) and 3(b) depict the computed value functions  for the years 2021 and 2022, respectively, without the sustainability concern.Similarly, Figures 4(a) and 4(b) show the computed value functions

Figure 5 .
Figure 5. Computed optimal harvesting policy and several controlled paths of populations for (a) 2021 and (b) 2022 without sustainability concern.All computed controlled paths did not cross but were very close to one another.

Figure 6 .
Figure 6.Computed optimal harvesting policy and several controlled paths of populations for (a) 2021 and (b) 2022 with sustainability concern.All computed controlled paths did not cross.

Figure 7 .
Figure 7. Computed paths of distorted mean body weights subject to different levels of uncertainty aversion for 2021 (a) without sustainability concern and (b) with sustainability concern.The values of theuncertainty aversion  were + (black, no uncertainty), 0.1 (blue), 0.01 (green), and 0.001 (red).

2 r
Kr are constants) for simplicity and also assume the stochastic logistic model with the fluctuation term of the form d the 1-D standard Brownian motion t B (e.g.,[48]), then the assumption stated above suggests that  should satisfy 2  by linearizing the nonlinearity of the logistic growth.This gives a rough upper bound of  .
. Using(44) and Lemma 2, we obtain ( ) the right derivative in the inequalities above.The second statement concerning the increasing property follows from the proof of Lemma 5.2, 3(a) of[55].□ of individuals with a distributed growth rate ( )

Table 1 .
Fitted parameter values of the uncertain logistic model in each year.
59)owing to the increasing property of h .Consequently, for each   , and fixing ,