1.
Introduction
Fish are a natural resource that many human populations rely on for food. Fishery managers, ecologists, and bioeconomists are concerned about the possibility of several marine populations being overfished [1,2]. Regarding fishery management, it is important to develop harvesting strategies that both meet human demands and prevent overexploitation of fish. Methods like time-area closures, limitations on fishing seasons, and the utilization of catch quotas have been considered [1,3], but there has been a continual push towards incorporating techniques that address these issues in an optimal way [2].
Optimal control theory has been a significant mathematical tool for fishery management, and one of the key findings obtained from control theory is in how the locations of the marine reserves (regions where fishing is prohibited) can be advantageously placed. Early applications of optimal harvesting problems over metapopulations [4,5,6,7] indicate that marine reserves used for fishery management can have economic benefits, which was contrary to what prior works [8,9] suggest. As emphasized in [10,11], investigating harvesting strategies and their effects on marine populations relies on constructing spatially dependent models. Since then, optimal control theory has been applied to spatially explicit fishery models [1,3,12,13,14,15,16,17], and the beauty behind such problems is that marine reserves are not imposed a priori but can arise if the optimal harvesting strategy possesses regions where the control is zero. Investigations of such problems provide further evidence that the incorporation of no harvesting regions for prolonged periods of times can be used to improve economic quantities like the yield [3,12,14,15,16], rent [13], and discounted net profit value [1,17].
While optimal control theory is fundamental to the development of optimal fishery management strategies and in ascertaining the effects of harvesting on marine populations, some works [1,3,12,14,15,18] have raised concerns about whether or not the optimal control obtained can be realistically administered. Depending on the structure of the optimal control, one may have to resort to developing a method for finding an alternative suboptimal harvesting strategy that is implementable.
In [18], Demir and Lenhart investigated a system of ordinary differential equations (ODEs) representing a food chain model of Black Sea anchovies, and they used optimal control theory to find a seasonal harvesting strategy that maximizes the discounted net profit value. The optimal control contained a rapid decline of harvesting followed by a rapid increase in harvesting during each season. To relax this behavior, Demir and Lenhart suggested an alternative harvesting strategy that is piecewise constant where the jump discontinuity occurs in the middle of the harvesting season, and the corresponding net profit value was close to the true optimal value. In their extension of the problem to a spatially explicit model [1], they compared the optimal harvesting strategy generated with a constant harvesting strategy. In [3], Kelly, Xing, and Lenhart used a parabolic partial differential equation (PDE) to model the population density of the stock to be harvested over a multidimensional spatial and temporal domain. Logistic growth and movement through diffusion and advection were considered within the model. Their control problem aimed to find a spatiotemporal harvesting distribution that maximizes the discounted yield. Kelly et al. emphasized that the feasibility of implementing each of the optimal harvesting strategies that they found numerically is up for debate. To alleviate these concerns, they constructed approximations of the optimal control in some of their simulations by having the control be constant in space or in both space and time. Some of the associated yields were very close to the optimal yield, while others (under certain parameter conditions) deviated from the optimal yield by as much as 5%.
Another situation related to the topic of realistic implementation is in the case in which the optimal control problem is linear in the control. Harvesting problems aimed at maximizing the yield, like the ones presented in [3,12,14,15,16], are of this form. Solutions to such problems can have regions in which the control is constant with values being the bounded constraints assumed on the control, and the control is said to be bang-bang on those regions. Such regions are easy to interpret and implement. However, it is also possible that the control may possess a singular arc. In problems in which the state equations are a system of ODEs, the singular case occurs if there is a subset of the domain for which the partial derivative of the Hamiltonian with respect to the control is zero and that set has nonzero measure [19]. One major concern associated with the presence of a singular arc is the possibility of Fuller's phenomenon. Fuller's phenomenon, or chattering, is the phenomenon in which the optimal control possesses a singular arc that cannot be concatenated with the bang-bang arcs without prompting infinite oscillations over a finite region [19,20]. In such instances, one cannot realistically employ a chattering control, and an alternative strategy must be considered. What is even more challenging is that some numerical techniques are prone to generating oscillatory numerical artifacts when a singular arc is present regardless of whether or not that control truly chatters (see [37, Section 9.1] for a discussion on problems numerical solvers have with singular controls).
This leads us to a discussion of a harvesting problem that was originally presented by Neubert in [12]. On the surface, this problem looks innocent. It consists of just one state equation, which is a second order differential equation with Dirichlet boundary conditions, and it has only two parameters. In [12], Neubert was able to characterize the singular case, and on the singular region, both the control and the population density of the stock are constant. However, when he numerically solved the problem with parameters set to where a singular arc was present, wild oscillations were found near the boundary of the singular region. He attributed these oscillations to chattering, but the local and intrinsic order of the singular control are such that standard methods for determining whether or not chattering occurs, like the theorem of conjugation [20,21] and the junction theorem [22], cannot be applied. These issues potentially carried over to the following works [14,15] where the optimal control problem is extended to PDE versions of Neubert's spatially explicit model. Ding and Lenhart [14] prevented chattering from arising by adding an $ H^1_0 $ norm penalty term to the objective to minimize variation of the control. Consequently, the admissible control set was restricted to where the control is differentiable, and variational inequalities were needed in solving the problem. In all of Joshi, Herrera, and Lenhart's simulations, no singular control was found; however, they do mention in [15, Section 5] that they were unable to obtain numerical solutions to the optimal control problem when the habitat size and final time were large. They also mentioned that questions pertaining to whether or not the singular control appears in finite time and whether the presence of a singular arc produces chattering within their problem remains open.
In this paper, we revisit Neubert's problem [12] to study the following two questions: 1) How do we numerically assess whether or not a problem presents chattering in cases when we cannot analytically prove such a phenomenon? and 2) When an optimal control has regions that are difficult to implement, how can we systematically find alternative strategies that are realistically implementable and near optimal? Although the first question does not apply to all harvesting problems, fishery managers should be concerned about the second. In further analyzing Neubert's harvesting problem, we develop a strategy that seeks to answer the first question, which results in a procedure for answering the second. Our methods involve using an extended version of the SPA [23]. The SPA is a technique in which one numerically solves an optimal control problem by optimizing the objective with respect to the switches of the control rather than the control itself. For the original presentation of SPA, we direct the reader to our prior work [24] which also surveys other switch point-based methods for solving singular and bang-bang problems. In using SPA for the harvesting problem, we solve the problem with a preset number of switches, and we study the corresponding yields of the controls when the number of switches near the singular region is increased. In our numerical experiments, we obtain strong empirical evidence that the harvesting problem chatters, and we find three alternative harvesting strategies with fewer switches that are both realistic to implement and near optimal.
This paper is organized as follows. In Section 2, we present Neubert's harvesting problem [12] and summarize the analysis of the spatially explicit model and the characterization of the optimal control. Within this section, we also discuss the topic of the order of the singular arc for the harvesting problem and what this means with regard to chattering. In Section 3, we discuss the general idea of the extended version of SPA, and we discuss the assumptions needed in order to apply SPA. Implementing SPA requires a feedback form of the control and a good initial guess of the location of the switches. In Section 4, we describe how we obtained the feedback forms that we tested, and we present our numerical results that led to our conclusion that the control does chatter. We also approximate the maximum harvesting yield and the location of the singular region from the chattering control through Aitken's extrapolation. In the last section, we summarize our results and give the conclusion.
2.
Background on the harvesting problem
In this section, we present the harvesting problem which was investigated in Neubert [12]. We also discuss some key results that were discussed in Ding and Lenhart's extension [14] to the spatially explicit model that is used in the harvesting problem. We then summarize analysis pertaining to characterizing the optimal harvesting strategy and the singular case solution to an equivalent optimal control problem. The methods used in these analyses follow directly from what was presented in [12, Appendix]. We also discuss specific issues related to chattering and implementation to further motivate our interest in this problem.
In [12], Neubert uses the Schaefer harvest model (see [25] or [26]) and Fisher's equation (see [27]) to construct a spatially explicit model with harvesting. The partial differential equation presented measures the rate of change of the population density of the stock at a given location, which depends on logistic growth, movement through diffusion, and harvesting. Neubert takes the spatial domain to be a line segment of dimensionless length $ \ell $. Dirichlet boundary conditions are applied to represent the scenario of fish immediately dying if they come into contact with the boundary of the habitat. He constructs an optimal control problem over the PDE model with the objective functional measuring the harvesting yield. Through this construction, a marine reserve can only arise if the optimal control possesses a region for which the control is zero. He then simplifies the original optimal control problem by assuming that the stock is at equilibrium and by rescaling the variables to reduce the number of parameters used. The state equation then becomes the following spatially dependent second-order differential equation:
Here, the state variable $ u\colon[0, \ell]\to \mathbb{R} $ represents the population density of the fish relative to the carrying capacity, $ x $ is the spatial variable, and $ h $ is a harvesting control. The fishery management problem presented in [12] maximizes
over all harvesting strategies $ h\in U $ subject to the state equation (2.1) where
and $ h_{ \mathop{{\rm{max}}}}\geq 1 $ is the maximum harvesting rate. Note that if $ h(x) = h_{ \mathop{{\rm{max}}}} $ for all $ x\in[0, \ell] $ with $ h_{ \mathop{{\rm{max}}}}\geq 1 $, then the stock is driven to extinction.
2.1. Uniqueness and boundedness of the state and the existence of the optimal control
One observation worth noting is that for any control, $ h $, the state equation (2.1) has a trivial solution ($ u(x) = 0 $ for $ x\in[0, \ell] $). In analyzing the problem, we choose our state variable to be the nontrivial non-negative solution of the state system, where $ u $ is positive in $ (0, \ell) $. In [14], Ding and Lenhart used prior works ([28,29,30]) pertaining to investigations of an optimal control problem on a diffusive elliptic Volterra-Lotka-type equation to prove existence and uniqueness of a positive solution to an extended version of the spatially explicit model, where the extended model is defined over a multidimensional spatial domain that is bounded in $ \mathbb{R}^n $ with boundary $ C^1 $, and their results can be applied to (2.1). Additionally, they used an extension of the maximum principle applied to weak solutions [31] to show that $ 0\le u(x)\le 1 $. Moreoever, they showed the existence of an optimal control which maximizes the objective (2.2) subject to the constraint (2.3). Our problem (2.1) is a special case of their results corresponding to $ n = 1 $.
2.2. The problem in Mayer form
For simplifying discussions on applying the switch point algorithm to the harvesting problem (2.1)–(2.2), we consider a minimization problem that is equivalent to the original problem and rewrite the minimization problem into the Mayer form [32]. The equivalent minimization problem minimizes $ -J(h) $ over all admissible control strategies subject to the state equation (2.1). What makes the minimization problem equivalent to problem (2.1)–(2.2) is that the control pair $ (u^*(x), h^*(x)) $ that minimizes $ -J $ is the control pair that maximizes the yield.
We rewrite the state equations used for the equivalent minimization problem as a system of two first-order differential equations by letting $ y_1(x) = u(x) $ and $ y_2(x) = \dot{y}_1(x) $. The minimization problem is then
where $ \mathcal{U}(x) = \{{h}(x)\in \mathbb{R}: 0\le {h}(x)\le h_{ \mathop{{\rm{max}}}}\}. $ We re-express the above problem in Mayer form by introducing a new state $ y_3(x) $ where
for all $ x\in[0, \ell] $. The associated state equation is
We define the objective given in (2.4) as $ C\colon \mathbb{R}^3\to \mathbb{R} $, where $ C({ \boldsymbol y}(\ell)) = y_3(\ell) $. Then we denote $ \boldsymbol{y}(x) = [y_1(x), y_2(x), y_3(x)]^\top $ and define $ { \boldsymbol f}\colon \mathbb{R}^3\times \mathbb{R}\to \mathbb{R}^3 $ as
where $ { \boldsymbol f}({ \boldsymbol y}(x), h(x)) $ is a column vector representing the right-hand side dynamics of the state equations given in (2.4) and (2.6). The problem in Mayer form is
2.3. Characterizing the optimal control and the singular case
We summarize a characterization of the optimal control of the Mayer harvesting problem (2.8) and an analysis of the singular case solution. The procedures used here align with what is illustrated in Neubert's analysis of problem (2.1)–(2.2) (see [12, Appendix]).
We define the adjoint variables $ { \boldsymbol \lambda}(x) = [\lambda_1(x), \lambda_2(x), \lambda_3(x)]:[0, \ell]\to \mathbb{R}^3 $, which produce a row vector. The costate associated with (2.8) satisfies the linear differential equation
where $ \nabla_y { \boldsymbol f} $ denotes the Jacobian of the state dynamics (2.7) with respect to $ { \boldsymbol y} $. Since $ y_3 $ does not appear explicitly in $ { \boldsymbol f} $, $ \dot{\lambda}_3(x) = 0 $ for all $ x\in[0, \ell] $. Moreover, $ \lambda_3(x) = 1 $ for all $ x\in[0, \ell] $ by the terminal condition. By explicitly computing $ \nabla_y { \boldsymbol f} $ and using $ \lambda_3(x) = 1 $, we rewrite the remaining two adjoint equations in (2.9) as
The Hamiltonian to problem (2.8) is
which can be rewritten as
Under the assumptions of the Pontryagin minimum principle, a local optimal pair $ ({ \boldsymbol y}^*(x), h^*(x)) $ of (2.8) and the associated costate $ { \boldsymbol \lambda}^*(x) $ have the property that
for almost every $ x\in[0, \ell] $. Since $ h\in U $, it is necessary that the following holds
where $ \psi(\cdot) $ is the switching function
The following theorem describes the singular case solution and both the local and intrinsic order of the singular control.
Theorem 2.1. If the optimal solution to problem (2.8) has a singular region $ V $ (that is, $ V\subset[0, \ell] $ with nonzero measure such that $ \psi $ vanishes identically on $ V $), then the optimal control, state, and adjoints are as follows:
Moreover, the intrinsic order of the singular arc is $ q = 1 $, and the local order of the singular arc is $ p = 2 $. Additionally, the intrinsic order version of the generalized Legendre-Clebsch condition (GLC) [22, Theorem 2.4] holds trivially on $ V $, and the local order version of the strengthened GLC [22, Theorem 2.7] holds on $ V $.
Proof. Equation (2.15) can be verified by using techniques demonstrated in Neubert [12, Appendix], which involves computing multiple higher-order derivatives of $ \psi(\cdot) $ up until the fourth order, setting those derivatives and $ \psi $ to zero on region $ V $, and solving for the control, states, and costates. Regarding the derivatives of the switching function for the Mayer problem, we have the following:
The second derivative of $ \psi $ has $ h $ appear explicitly, but the coefficient of $ h $ in the second derivative of $ \psi $ is $ 2\psi $, which vanishes on $ V $. The singular control is of intrinsic order $ q = 1 $, and the intrinsic order of the GLC [22, Theorem 2.4] is trivially satisfied on $ V $ since
The fourth derivative of $ \psi $ depends explicitly on $ h $. Using the formulas (2.15) and the fact that the second derivative of $ \psi $ vanishes on $ V $, we have
Thus, the singular control is of local order $ p = 2 $, and the local order version of the strengthened GLC [22, Theorem 2.7] also holds on $ V $. □
What cannot be solved explicitly is the expressions for $ { \boldsymbol y} $, $ \lambda_1 $, and $ \lambda_2 $ on the nonsingular regions. Also, the locations of the switches can only be approximated numerically. A switch of a control variable is a location within the domain of the control for which the qualitative behavior of the control changes, and these behavioral changes are based upon the sign of the switching function. That is, switches are the $ x $-values for which $ h^* $ changes from one nonsingular arc ($ h^*(x) = h_{ \mathop{{\rm{max}}}} $) to another nonsingular arc ($ h^*(x) = 0 $) or for which $ h^* $ changes from a singular arc to a nonsingular arc (or vice versa).
2.4. Discussion on chattering
When an optimal control possesses a singular arc, one should investigate if the optimal control chatters. In practice, one uses the junction theorem (see McDanell and Powers' [33]) to determine if the optimal control with a singular arc does not chatter, and one uses the theorem of conjugation [20,21] to determine if it chatters. Neither theorems are applicable for problem (2.8) due to the local and intrinsic order of the singular arc.
As mentioned by Lewis in [22], the junction theorem relies on the order of the singular arc, and there is some ambiguity in the usage of the word "order". If by "order" in the junction theorem, one means intrinsic order, then the theorem is not applicable in the case where the GLC conditions are only trivially satisfied, which is the case for problem (2.8). Moreover, if by "order" we mean local order, then the junction theorem is no longer true. He supports the latter statement by presenting a problem for which the local order version of the theorem fails. One of the hypotheses of the theorem of conjugation is that the singular arc is of intrinsic order 2 or an arbitrary even intrinsic order. The order of the singular arc for the harvesting problem is $ q = 1 $, and, thus, we cannot apply this theorem.
There is another analytical technique for proving chattering which involves converting the Hamiltonian system to the semiconical form (see [20, Section 1.3]). In [12], Neubert uses this approach in his analysis of chattering for the harvesting problem; however, the complementary functions for ensuring nondegeneracy of the Jacobi matrix mapping are not given. Hence, an analytic proof of chattering for problem (2.8) is an open problem. We investigate numerically whether this problem chatters by applying a SPA. Using SPA fixes the number of switches being assumed on the control, so oscillations of the control are predetermined and not a product of numerical artifacts. If we observe a trend in which a control with more oscillations near the singular arc produces an increase in the yield, then we have empirical evidence in favor of chattering, and we find near optimal harvesting strategies.
3.
The switch point algorithm
We begin by discussing the general idea of the SPA which is applied to optimal control problems whose solutions are either singular, bang-bang, or concatenations of bang-bang and singular arcs. We modify the notation that is used in our original presentations of SPA [24] and the extended SPA [23] to best describe how SPA can be applied to the harvesting problem (2.8).
The switch point algorithm is a method to numerically solve an optimal control problem by reducing the problem to an optimization over the switches. In [24], formulas are derived for the derivative of the objective with respect to the switches so that gradient-based optimization techniques can be applied. The control problems investigated in [24] concern state equations having initial conditions. In [34], it is shown that SPA also applies to the harvesting problem, while in [23], this extension is given for a general class of boundary-value problems. Our discussion utilizes the setup given in [23] with a focus on the harvesting problem.
The general fixed terminal time control problem studied in [23] is of the form
where $ { \boldsymbol y}\colon[0, \ell]\to \mathbb{R}^n $ is absolutely continuous, $ \boldsymbol{h}\colon[0, \ell]\to \mathbb{R}^m $ is essentially bounded, $ C\colon \mathbb{R}^n\to \mathbb{R} $, $ { \boldsymbol f}\colon \mathbb{R}^n\times \mathbb{R}^m\to \mathbb{R}^n $, $ \mathcal{U}(x) $ is a closed and bounded set for each $ x\in[0, \ell] $, $ I $ and $ E $ are subsets of $ \{1, 2, \dots, n\} $, $ { \boldsymbol y}_I $ denotes the components of $ { \boldsymbol y} $ associated with indices $ i\in I $, and $ { \boldsymbol y}_E $ denotes the components of $ { \boldsymbol y} $ associated with indices $ j\in E $. The vectors $ \boldsymbol{b}_I $ and $ \boldsymbol{b}_E $ are given initial and terminal values for the states. In [23], it is assumed $ |I|+|E| = n $ (here $ |\cdot| $ denotes the cardinality of the set). The harvesting problem (2.8) is of form (3.1) where $ n = 3 $, $ m = 1 $, $ { \boldsymbol f} $ is given in (2.7), $ I = \{1, 3\} $, $ E = \{1\} $, $ \boldsymbol{b}_I = [0, 0]^{\top} $, and $ { \boldsymbol b}_E = 0 $.
For simplicity of discussion, let us assume that the optimal control $ \boldsymbol{h}^* $ to problem (3.1) has one switch $ s^* $. The optimal control $ \boldsymbol{h}^*(x) $ is of the following feedback form with $ s = s^* $:
for some functions $ \boldsymbol{\phi}_0 $ and $ \boldsymbol{\phi}_1 $ defined on a larger interval containing $ (0, s) $ and $ (s, \ell) $, respectively. The feedback functions used for representing the control in the harvesting problem (2.8) are constant functions where $ \phi_i\in\{0, h_{ \mathop{{\rm{max}}}}, 0.5\} $.
For SPA, problem (3.1) is solved by optimizing over the choice of $ s $. Define $ { \boldsymbol F}_0({ \boldsymbol y}, x): = { \boldsymbol f}({ \boldsymbol y}, \boldsymbol{\phi}_0({ \boldsymbol y}, x)) $, $ { \boldsymbol F}_1({ \boldsymbol y}, x): = { \boldsymbol f}({ \boldsymbol y}, \boldsymbol{\phi}_1({ \boldsymbol y}, x)) $, and
We replace problem (3.1) by
In [23], we derive a formula for the partial derivative of the cost with respect to $ s $, and this formula allows us to utilize gradient-based methods to find $ s^* $. Let $ C(s) $ denote the objective in (3.2) which depends on $ s $. Under a smoothness assumption on each $ { \boldsymbol F}_i $ and invertibility assumptions for submatrices of the related fundamental matrices, we can generate a formula for finding the partial derivative of the cost with respect to $ s $. This formula depends on the Hamiltonian where we find the difference between the value of the Hamiltonian when the control is $ \boldsymbol{\phi}_0({ \boldsymbol y}, x) $ at $ x = s $ and the value of the Hamiltonian when the control is $ \boldsymbol{\phi}_1({ \boldsymbol y}, x) $ at $ x = s $. The formula for the partial derivative of the cost with respect to the switch is
where $ H_i({ \boldsymbol y}, { \boldsymbol \lambda}, x) = { \boldsymbol \lambda} { \boldsymbol f}({ \boldsymbol y}, \boldsymbol{\phi}_i(x)) $ with $ i\in\{0, 1\} $ is the Hamiltonian, and the row vector $ { \boldsymbol \lambda}\colon[0, \ell]\to \mathbb{R}^n $ is the adjoint vector, which is the solution to the linear differential equation:
Here, $ J $ and $ K $ denote the complements of $ I $ and $ E $, respectively, and $ \nabla_KC $ is a row vector whose $ k $-th component (with $ k\in K $) is the partial derivative of $ C $ with respect to $ { \boldsymbol y}_k $. In the context of the harvesting problem (2.8), $ J = \{2\} $, $ K = \{2, 3\} $, and $ \nabla_KC = [0\; 1] $, which matches the initial and terminal conditions used in (2.9).
3.1. Assumptions for formula (3.3)
We now discuss the specific assumptions used in obtaining (3.3). The first assumption is a dynamics smoothness assumption:
Assumption 1. (Dynamics smoothness) For $ r > 0 $, define the tubes
where $ \mathcal{B}_r(\boldsymbol{c}) = \{\hat{ { \boldsymbol y}}\in \mathbb{R}^n : {\left\| {\hat{ { \boldsymbol y}}- \boldsymbol{c}} \right\|}\leq r\} $. It is assumed that on $ \mathcal{T}_j $, $ j = 0, 1 $, $ { \boldsymbol F}_j $ is continuously differentiable, while $ { \boldsymbol F}_j(\mathcal{Y}, x) $ is Lipschitz continuously differentiable in $ \mathcal{Y} $, uniformly in $ x $, with Lipschitz constant $ L $.
In showing that the Lipschitz condition in the dynamics smoothness assumption holds for problems of form (3.1) and rewritten as (3.2), one can first show that the feedback forms $ \boldsymbol{\phi}_0({ \boldsymbol y}(x), x) $ and $ \boldsymbol{\phi}_1({ \boldsymbol y}(x), x) $ are Lipschitz continuous in $ { \boldsymbol y} $ over its respective domains in $ x $. It follows that on $ \mathcal{T}_j $, $ j\in\{0, 1\} $, $ { \boldsymbol F}_j $ is Lipschitz continuously differentiable in $ \mathcal{Y} $ if both $ { \boldsymbol f} $ and $ \nabla_y { \boldsymbol f} $ are Lipschitz continuous in $ ({ \boldsymbol y}, { \boldsymbol h}) $. For the harvesting problem, it follows from (2.13) and Theorem 2.1 that the optimal feedback control $ h^*(x) $ is 0, 1/2, and $ h_{ \mathop{{\rm{max}}}} $ when $ \psi ({ \boldsymbol y}^*(x), h^*(x), { \boldsymbol \lambda}^*(x)) $ is positive, 0, and negative, respectively. Hence, in the particular case $ h_{ \mathop{{\rm{max}}}} = 1 $ that we consider in the experiments, there are three distinct $ \boldsymbol{\phi}_i $ functions, the constant functions 0, 1/2, and 1. Since the state $ { \boldsymbol y} $ lies on $ [0, 1] $, dynamic smoothness holds if $ { \boldsymbol f} $ and $ \nabla_y { \boldsymbol f} $ are Lipschitz continuous in $ ({ \boldsymbol y}, { \boldsymbol h}) $. In Appendix 5, we utilize [34] to show that this Lipschitz continuity property holds for problem (2.8).
In addition, we have two invertibility assumptions for submatrices of related fundamental matrices that are used for generating (3.3). The first fundamental matrix $ \boldsymbol{\Phi}\colon[0, \ell]\to \mathbb{R}^{n\times n} $ is associated with the state dynamics:
where $ \boldsymbol{I} $ is the $ n\times n $ identity matrix. Formula (3.3) relies on the assumption that $ \boldsymbol{\Phi}_{EJ}(\ell) $ is invertible. In the case of the harvesting problem (2.8), $ \boldsymbol{\Phi}_{EJ}(\ell) = \Phi_{12}(\ell) $ is a scalar, which is the (1, 2) entry of $ \boldsymbol{\Phi}(\ell) $. In this case, the assumption simplifies to $ \Phi_{12}(\ell)\neq0 $.
The second fundamental matrix $ \boldsymbol{\Psi} $ is associated with the costate equation (3.4):
Formula (3.3) relies on the assumption that $ \boldsymbol{\Psi}_{IK}(\ell) $ is invertible. For the harvesting problem (2.8), submatrix $ \boldsymbol{\Psi}_{IK}(\ell) $ is a $ 2\times 2 $ matrix, where the $ (i, j) $ entry of $ \boldsymbol{\Psi}(\ell) $ is an entry of $ \boldsymbol{\Psi}_{IK}(\ell) $ for $ i\in I = \{1, 3\} $ and $ j\in K = \{2, 3\} $. Since $ { \boldsymbol F} $ is independent of $ y_3 $ in the harvesting problem, the last column of $ \nabla_y { \boldsymbol F}({ \boldsymbol y}(x), x) $ is zero. Hence, the derivative of the last column of $ \boldsymbol{\Psi} $ is zero, which implies that the last column of $ \boldsymbol{\Psi} $ is the last column of the identity matrix, and the last column of $ \boldsymbol{\Psi}_{IK}(\ell) $ is $ [0 \; 1]^\top $. Consequently, $ \boldsymbol{\Psi}_{IK} $ is a lower triangular matrix which is invertible if and only if $ \Psi_{12}(\ell) \ne 0 $. In summary, for the harvesting problem, the formula (3.3) requires that the $ (1, 2) $ matrix elements of $ \boldsymbol{\Phi}(\ell) $ and $ \boldsymbol{\Psi}(\ell) $ are nonzero.
We then arrive at the following theorem.
Theorem 3.1. [23, Theorem 5.1] If Assumption 1 holds, the objective $ C $ is continuously differentiable, and both submatrices $ \boldsymbol{\Phi}_{EJ}(\ell) $ and $ \boldsymbol{\Psi}_{KI}(\ell) $ are invertible, then the formula for $ \frac{\partial C}{\partial s} $ given in (3.3) holds.
The two invertibility assumptions in Theorem 3.1 are key for both the gradient formula (3.3) and for numerical optimization of the switch point in (3.2). In order to optimize the objective over the switch point, we need to consider different choices for $ s $. A change in $ s $ to improve the objective values requires an update in the solution $ { \boldsymbol y} $ to the boundary-value problem in (3.2). This update can be accomplished using both Newton's method and the fundamental matrix (3.5). Recall that the solution $ \boldsymbol{\Phi} $ of (3.5) satisfies
where $ { \boldsymbol y} $ is the solution to the initial-value problem
Now, let us consider the initial-value problem
One way to satisfy the boundary-value problem (3.2) is to choose $ \boldsymbol{\theta} $ so that $ { \boldsymbol y}_E(\ell) = { \boldsymbol b}_E $. Numerically, $ \boldsymbol{\theta} $ can be computed using the Newton iteration
where $ { \boldsymbol y}^k $ is the solution of (3.7) associated with $ \boldsymbol{\theta} = \boldsymbol{\theta}_k $ and $ \boldsymbol{\Phi}_{EJ}^k $ is the $ EJ $ submatrix of the fundamental matrix $ \boldsymbol{\Phi} $ in (3.5) associated with $ { \boldsymbol y} = { \boldsymbol y}^k $. Thus, the same fundamental matrix that appears in Theorem 3.1 also appears in Newton's method when we update the solution of the boundary-value problem (3.2) after a change in the switch point $ s $.
The formula for the derivative of the objective with respect to $ s $, formula (3.3), also involves the solution $ \boldsymbol{\lambda} $ to the differential equation in (3.4). The solution to this differential equation can be expressed in terms of the fundamental matrix $ \boldsymbol{\Psi} $ in (3.6): $ \boldsymbol{\lambda}(\ell) = $ $ \boldsymbol{\lambda}(0) \boldsymbol{\Psi}(\ell) $. Since $ \boldsymbol{\lambda}_J(0) = \boldsymbol{0} $ by the boundary conditions given in (3.4), it follows that
Again, by the boundary conditions in (3.4), we have $ { \boldsymbol \lambda}_K(\ell) = \nabla_{K}C({ \boldsymbol y}(\ell)) $, which yields
when the submatrix $ \boldsymbol{\Psi}_{IK}(\ell) $ is invertible.
In summary, when optimizing a switch in the harvesting problem, we can at the same time check that the invertibility assumptions of Theorem 3.1 are satisfied, update $ { \boldsymbol y} $ to satisfy the boundary condition $ { \boldsymbol y}_E(\ell) = \boldsymbol{0} $, and evaluate the derivative of the objective with respect to the switch point. More precisely, after an optimization algorithm adjusts a switch point to improve the objective value, Newton's method and the fundamental matrix $ \boldsymbol{\Phi} $ can be used to iteratively update the state $ { \boldsymbol y} $ so as to satisfy the boundary condition $ { \boldsymbol y}_E(\ell) = 0 $ using algorithm (3.8). Using the updated $ { \boldsymbol y} $ together with $ { \boldsymbol \lambda}_K(\ell) = \nabla_{K}C({ \boldsymbol y}(\ell)) $ and $ \boldsymbol{\Psi}(\ell) $, we can use (3.9) to evaluate $ \boldsymbol{\lambda}_I(0) $. Finally, the known initial condition $ \boldsymbol{\lambda}(0) $ and the fundamental matrix values $ \boldsymbol{\Psi}(s) $ yield $ \boldsymbol{\lambda}(s) = \boldsymbol{\lambda}(\boldsymbol{0}) \boldsymbol{\Psi}(s) $ and the derivative of the objective with respect to the switch point from (3.3).
3.2. Generalizing formula (3.3)
As mentioned in [23], we can generalize these results to an arbitrary number of switches. Let $ 0 < s_1 < s_2 < \cdots < s_N < \ell $ be a collection of switch points and define $ s_0 = 0 $ and $ s_{N+1} = \ell $. Assuming that the boundary condition $ { \boldsymbol y}_E(\ell) = \boldsymbol{b}_E $ is satsified, the gradient of the objective is evaluated by fixing all the switch points but one, say $ s_i $, and evaluating the derivative with respect to $ s_i $. Hence, by cycling through all of the switch points, it is possible to evaluate the entire gradient vector. This gradient can be evaluated using just one integration of the state equation and the fundamental matrix $ \boldsymbol{\Psi} $. Note that these integrations need to be done interval by interval, from $ s_i $ up to $ s_{i+1} $ for $ i = 0 $ up to $ i = N $, since the control has a different feedback form $ \boldsymbol{h}(x) = \boldsymbol{\phi}_i({ \boldsymbol y}(x), x) $ on each interval. We first integrate the state equation and the fundamental matrix $ \boldsymbol{\Psi} $ across each interval until reaching $ { \boldsymbol y}(\ell) $ while storing the states, controls, and fundamental matrices at the ends of each interval $ [s_i, s_{i+1}] $. Using (3.9), we can evaluate $ \boldsymbol{\lambda}_I(0) $, and then use $ \boldsymbol{\Psi} $ to evaluate the costate at the end of each interval $ [s_i, s_{i+1}] $, which yields the partial derivative of the objective with respect to each switch point using (3.3).
Let us define
where $ { \boldsymbol s} = [s_1, \ldots, s_N] $ is an $ N $ vector whose components are the switches, while $ s_0 = 0 $ and $ s_{N+1} = \ell $. With this notation, the generalization of (3.2) is
Let $ C({ \boldsymbol s}) $ denote the objective in (3.10) parameterized by the switches $ s_i $, $ 1 \le i \le N $. In [23], the formula
is obtained under the same smoothness assumption and invertibility assumptions for submatrices of related fundamental matrices where $ H_i({ \boldsymbol y}, { \boldsymbol \lambda}, x) = { \boldsymbol \lambda} { \boldsymbol f}({ \boldsymbol y}, \boldsymbol{\phi}_i(x)) $ and $ { \boldsymbol \lambda} $ solves (3.4).
4.
Applying SPA to the harvesting problem
We use SPA to obtain some empirical evidence in determining whether or not the harvesting problem (2.8) chatters when $ h_{ \mathop{{\rm{max}}}} = 1 $ and $ \ell = 10 $, where such settings produce a singular region. The harvesting problem (2.8) is recast as (3.10) where the control is allowed $ N $ switches and the objective is optimized with respect to those switches. We utilize formula (3.11) to compute the gradient of the objective, and we use the Polyhedral Active Set Algorithm (PASA) [35,36] to solve (3.10). Let $ { \boldsymbol s}^0 = [s_1^0, s_2^0, \dots, s_N^0] $ be our starting guess for the switches, and let $ h^0 $ be the initial control with a predetermined feedback form and with switches $ { \boldsymbol s}^0 $. Then, let $ { \boldsymbol s}^k $ be the $ k $th approximation of the switches and update the switches $ { \boldsymbol s}^{k+1} $ by using PASA to solve the following subproblem:
where $ \boldsymbol{b}_I = [0, 0]^{\top} $, $ I = \{1, 3\} $, $ E = 1 $, $ \delta > 0 $ is a small number, and
Problem (4.1) is similar to problem (3.10) but with an added constraint $ s_i\in[s^k_i-\delta, s^k_i+\delta] $. This constraint is added since the updates to $ { \boldsymbol y} $ after a change in $ \boldsymbol{s} $ is accomplished using Newton's method, which is locally convergent. To prevent divergence, we need to require small changes in the switch points. In each iteration, we update the box constraints, $ [s^k_i-\delta, s^k_i+\delta] $, to where the constraints are centered about the previous iterate. We repeat this procedure until $ s_i^{k+1} $ lies in the interior of $ [s_i^k-\delta, s^k_i+\delta] $ for all $ i = 1, \dots, N $, which yields a local solution to problem (3.10).
When solving (4.1), we use the MATLAB ODE solver ode45 to solve the initial value problem:
For all uses of ode45, the maximum step size is set to $ 10^{-2} $, and both the absolute and relative error tolerances are $ 10^{-10} $.
We perform multiple experiments where we numerically solve (3.10) with $ N = 4, 6, 8 $, and 10 switches and with the control being of a particular feedback form. The control form we investigate, $ h_{N, sing} $, is a concatenation of bang-bang and singular arcs with the singular arc located in the center of the habitat. The other control form we investigate is a bang-bang control that has $ N = 6 $ switches, which we denote as $ h_{6, bang} $. Table 1 indicates the regions of maximum harvesting or no harvesting for each control being analyzed, and these regions are represented as open intervals whose endpoints are either a switch $ s_i $ or the boundary of the habitat (0 or $ \ell $). The vector of switches obtained by solving (3.10) when $ h $ is of feedback form $ h_{N, sing} $ is denoted as $ \boldsymbol{s}^* = [s_1^*, \dots, s_{N}^*] $, and $ h_{N, sing}^* $ is the resulting control. We use the same notation when it comes to $ h_{6, bang} $. Since problem (3.10) restricts the admissible control set to controls of a specific form, further tests are performed to determine if $ h^*_{N, sing} $ optimizes the yield over all possible strategies. One test that we consider is to determine if each control satisfies Pontryagin's minimum principle by studying the graph of the switching function. We also investigate the changes in the yields as we increase $ N $.
4.1. Motivation behind the feedback forms
As suggested in [24], we can use total variation (TV) regularization prior to using SPA to generate a feedback form of the optimal control. As demonstrated in our prior work [37], oscillatory numerical artifacts arising from a non-chattering singular control problem can be tuned by adding a regularization term that measures some tuning parameter $ \rho $ multiplied by the total variation of the control to the cost functional. Caponigro et al. [38] also illustrates that one can apply TV regularization to chattering control problems, in which case the regularized control obtained is a non-oscillatory quasi-optimal control. This means that the control is not optimal, but as the tuning parameter $ \rho $ approaches zero, the resulting regularized cost value approaches the true cost value. In our prior work [34], we apply TV regularization to the harvesting problem (2.4) for varying values of the tuning parameter. From our simulations, we find two control forms, namely $ h_{4, sing} $ and $ h_{6, sing} $ (see Table 1), that are potential candidates of the structure of the optimal harvesting strategy. It is worth mentioning that both regularized controls possess some numerical artifacts on the singular region in that the control did not equal 1/2 throughout the singular region but was approximately 1/2, and this had some effects on how the resulting state $ y_1 $ appeared along the singular region. Additionally, due to the discretization process of the problem, the approximation of the switches are mesh point values that were used in partitioning the habitat. With SPA, we fix the values of the control as being 0, 1/2, or $ h_{ \mathop{{\rm{max}}}} $, so the numerical artifacts arising from the TV regularization approach are no longer present, and we use the switches obtained from TV regularization as our initial guess of the switches.
We discuss later in the results section how the sign of the switching function corresponding to $ h_{4, sing}^* $ and $ h_{6, sing}^* $ along the singular region are used to determine that neither control is optimal and how the switching function can be used to construct more control forms to investigate, which led to our investigation of control forms $ h_{8, sing} $, $ h_{6, bang} $, and $ h_{10, sing} $ (see Table 1). The feedback form $ h_{8, sing} $ is a modification of $ h_{6, sing}^* $ where a switch is added to the left and right of the singular region. The added switches lead to the introduction of two marine reserves that are directly adjacent to the singular region. We use the switches from $ h_{6, sing}^* $ as our starting guess for six of the eight switches of $ h_{8, sing} $, and for the initial guess of the two remaining switches, which correspond to the boundary of its singular region, we take the boundary of $ h_{6, sing}^* $'s singular arc and perturb those boundary points by 0.05. Similarly, $ h_{10, sing} $ takes the control form $ h_{8, sing}^* $ and adds a switch to the left and right of the singular region which leads to the introduction of two regions of maximum harvesting near the singular arc. Our interest in analyzing these forms is to investigate a potential trend where increasing the switches near the singular region results in an increase in the harvesting yield. Such a trend provides some empirical evidence that the optimal harvesting strategy requires infinite switches near the singular region.
Additionally, we consider a control form that is bang-bang with 6 switches, $ h_{6, bang} $. The beauty behind the harvesting problem (2.8) is that the singular control is constant and is a value easily implementable, namely 1/2. Such a solution may be realistic if we consider half of the maximum harvesting rate ($ h_{ \mathop{{\rm{max}}}} = 1 $) based upon resources. For example, the number of boats or fishermen assigned to fish over the singular region could be half of the number of boats or fishermen that are assigned to the regions in which maximum harvesting is applied. However, not all singular control problems have a singular case solution that is easy to administer. We want to address this concern by investigating whether or not a bang-bang control could be used as an alternative. Theoretically, such a control produces a smaller yield, but we can numerically determine if the resulting harvesting yield is comparable to the yields produced by other control forms. We obtain the control form $ h_{6, bang} $ by removing the singular arc obtained from $ h_{6, sing}^* $ and replacing it with a marine reserve. The starting guess of the switches to $ h_{6, bang} $ came from taking the switches of $ h_{6, sing}^* $ and perturbing the values of those switches by $ \pm $0.04 to where the corresponding state solution $ y_1 $ is non-trivial and positive.
4.2. Results
The solution to the harvesting problem for the control geometries summarized in Table 1 was computed using the algorithm (4.1). Plots of the optimal solutions are shown in Figures 1 and 2, while Table 2 shows the numerical values for the switches and the yields*. The computing tolerance for the optimization code PASA was $ 10^{-10} $ to ensure accurate estimates for the yield associated with each control geometry. In the first column of Table 2, we also provide the choice for $ \delta $ in $ \hat{S}^{k} $ that ensured small changes in the switch points and convergence of the Newton iterations used to satisfy the boundary conditions. As the number of switch points increased, the optimization problem became more ill conditioned, and a smaller $ \delta $ was needed. As mentioned previously, the optimal solution $ h^* $ to problem (2.8) minimizes the objective $ C(h^*) = y_3(\ell) $ or equivalently maximizes the yield $ -C(h^*) $, and the control should align with Pontryagin's minimum principle. If Pontryagin's minimum principle is violated, then the computed solution is not locally optimal, and the objective value can be improved by introducing another pair of switch points. We use the signs of the resulting switching function $ \psi(x) $ along the bang-bang and singular regions to investigate whether or not the minimum principle holds. Since PASA's error tolerance is $ 10^{-10} $, we deduce a violation in the minimum principle if $ |\psi(x)| > 10^{-10} $ over the singular region.
*Simulations pertaining to these results can be found at https://github.com/srnatkins/SPA_Applied_to_Harvesting_Problem.
The resulting control $ h_{4, sing}^* $ (solid) and associated state variable $ y_1(x) $ (dash-dot) obtained from using SPA are given in Subfigure 1a and the resulting values of the switches are given in the $ h_{4, sing}^* $ row of Table 2. In Figure 1(b), a plot of the associated switching function $ \psi(x) $ is provided. Observe that $ \psi(x) $ is negative over the singular region with values on the order of $ 10^{-3} $. This allows us to conclude that the control $ h_{4, sing}^* $ is not optimal.
The behavior of the resulting state variable $ y_1(x) $ near the boundary of the singular arc plays a large role as to why $ h_{4, sing}^* $ is not optimal. Observe from Subfigure 1a that $ y_1(x)\approx0.5 $ on the singular region, but near switches $ s_2^* $ and $ s_3^* $, $ y_1(x) $ is greater than 0.5. While SPA allows us to force the control to remain constant along the singular region, we cannot force the state to remain constant. In numerically solving for the state, the solver is trying to maintain a differentiable solution of $ y_1(x) $ at switches $ s_2 $ and $ s_3 $. We support this claim by discussing how the harvesting strategy influences the curvature of the state and refer back to state equation (2.1) since $ y_1(x) $ corresponds to state variable $ u $. In Subsection 2.1, the population density of the stock is bounded above by one. Then, based upon the second order state equation (2.1), we have the following conditions:
Consequently, we have that $ y_1 $ is concave down on regions where $ h $ is zero and concave up on regions where $ h $ is $ h_{ \mathop{{\rm{max}}}} = 1 $. These are significant conditions when considering the locations where a singular and non-singular arc are concatenated. Based on the structure of $ {h}_{4, sing}^* $, state variable $ y_1 $ should switch from concave down to constant at $ x = s_2 $. Additionally, $ y_1 $ should switch from constant to concave down at $ s_3 $. In order for $ y_1 $ to remain differentiable at $ x = s_i $ for $ i\in\{2, 3\} $, $ y_1(x) $ should have a local maximum value of $ 0.5 $ at $ x = s_i $, which is not occurring here even after using SPA to improve the switches.
Figure 1(c), (d) display the control $ h_{6, sing}^* $, the state variable $ y_1(x) $, and the switching function $ \psi(x) $. The resulting switches are given in the $ h_{6, sing}^* $ row of Table 2. The table also indicates that the harvesting yield obtained when administering $ h_{6, sing}^* $ is larger than the yield corresponding to $ h_{4, sing}^* $. Observe in Figure 1(d), marine reserves are bordered by regions of maximum harvesting. In analyzing Figure 1(d), $ \psi(x) $ seems to lie on the $ x $-axis on the singular region; however, $ \psi(x) $ also appears to be lying along the $ x $-axis on other regions. Over the intervals $ (s^*_2, s^*_3) $ and $ (s^*_4, s^*_5) $, $ h_{6, {sing}}^*(x) = h_{ \mathop{{\rm{max}}}} $, so by the Pontryagin's minimum principle, the switching function should be negative. Figure 3(a) displays the switching function over intervals $ (s^*_2, s^*_3) $, $ (s^*_3, s^*_4) $, and $ (s^*_4, s^*_5) $. On $ (s^*_2, s^*_3) $ and $ (s^*_4, s^*_5) $, the switching function is below the $ x $-axis, but on the singular region, $ (s^*_3, s^*_4) $, $ \psi(x) $ is positive with its values on the order of $ 10^{-5} $. Hence, we deduce that $ h_{6, {\rm{sing}}}^* $ is not optimal. The sign of the switching function along the singular region is influential in our consideration for $ h_{8, sing} $ and $ h_{6, bang} $. Based upon Pontragin's minimum principle, $ h_{6, sing}^* $ should have had a marine reserve instead of a singular arc on $ (s_3^*, s_4^*) $. With this in mind, we investigate what happens when we take the switches and structure of $ h_{6, sing}^* $, add a marine reserve to both sides of the boundary of the singular region ($ h_{8, sing} $), and use SPA to improve those switches to see the resulting switching function is closer to zero along the singular arc. We also consider taking the switches and structure of $ h_{6, sing}^* $, but reassign the singular arc as being a marine reserve (ie $ h_{6, bang} $) and use SPA to improve those switches.
Figure 2(a) displays $ h_{8, sing}^* $ (solid) and the corresponding state $ y_1 $ (dash-dot). The values of the switches are given in Table 2. Observe also in the last column of the table that administering $ h_{8, sing}^* $ resulted in a larger yield in comparison to implementing $ h_{4, sing}^* $ and $ h_{6, sing}^* $. This illustrates a potential trend that increasing the number of switches along the singular region may result in a larger yield. In Figure 2(b), we see that the switching function $ \psi(x) $ associated with $ h_{8, sing}^* $ appears to lie on the $ x $-axis on the singular region, $ (s^*_4, s^*_5) $, but the switching function also appears to lie on the $ x $-axis on the non-singular regions $ (s^*_2, s^*_3) $, $ (s^*_3, s^*_4) $, $ (s^*_5, s^*_6), $ and $ (s^*_6, s^*_7) $. In Figure 3(b), we provide a picture of the switching function that is zoomed-in over the interval $ (s_3, s_6) $. On the regions corresponding to marine reserves, $ (s^*_3, s^*_4) $ and $ (s^*_5, s^*_6) $, the switching function lies above the $ x $-axis with the order of magnitude being $ 10^{-7} $. Regarding the singular region, the switching function is negative with order of magnitude between $ 10^{-7} $ and $ 10^{-6} $. The switching function corresponding to $ h_{8, sing}^* $ is much closer to zero over the singular region in comparison to the switching function obtained from $ h_{6, sing}^* $ and $ h_{4, sing}^* $; however, we conclude that $ h_{8, sing}^* $ is not optimal since $ \psi(x) $ does not meet the expected error tolerance.
In the $ h_{10, sing}^* $ row of Table 2, we have the resulting switches and the harvesting yield found when applying SPA. Its yield is greater than the yields corresponding to $ h_{8, sing}^* $, so this trend of the harvesting yield increasing as we increase the number of switches near the singular region continues. In Figure 2(c), the switching function corresponding to $ h_{10, sing}^* $ is given. The switching function does appear to lie along the $ x $-axis over the interval $ (s_3, s_8) $, but in looking at a zoomed-in picture of the switching function in Figure 3(c), we see that the switching function is positive over the singular region with an order of magnitude $ 10^{-8} $. Consequently, the control presented in Figure 2(d) is not locally optimal. Each time we add another pair of switch points to the control, the violation in the Pontryagin minimum principle decreases by about a factor of 0.01, while the yield improvement is about 0.001 times the improvement associated with the prior new pair of switch points.
We find the trend of an increase in switches leading to a larger harvesting yield as being empirical evidence in favor of the claim that the optimal harvesting strategy is chattering. Unfortunately, we cannot further test this trend for the case in which the control has 12, 14, or 16 switches because this pushes the limits of the optimization solver that we are using. Observe also in Figure 2(d) that the marine reserves near the singular arc (located on intervals $ (s_3, s_4) $ and $ (s_7, s_8) $) are quite small and that the regions of maximum harvesting to the left and the right of the singular region (intervals $ (s_4, s_5) $ and $ (s_6, s_7) $) are even smaller. In using the approximated switches displayed in row $ h_{10, sing}^* $ of Table 2, we have that the lengths of the marine reserves located on interval $ (s_3, s_4) $ and of the marine reserve located on interval $ (s_7, s_8) $ are both approximately 0.0775 units, and the lengths of the maximum harvesting regions to the left and right of the singular region are both approximately 0.0148 units. This behavior is very reminiscent of chattering in that the length of these bang-bang regions taper to where an infinite number of switches are present over a finite region. Although the $ h_{10, sing}^* $ control corresponds to a higher harvesting yield, these small regions near the singular arc do raise some concern when thinking about implementation of such controls. The same could be said for $ h_{8, sing}^* $.
We use SPA to improve the switches of our $ h_{6, bang} $ control. In Figures 1(e), (f), we provide a plot of the resulting $ h_{6, bang}^* $ control and of the corresponding switching function, and in Table 2, we provide the resulting switches. The sign of the switching function is positive in regions where the switching function should be negative, namely intervals $ (s_2, s_3) $ and $ (s_4, s_5) $. Also, on interval $ (s_3, s_4) $ and on parts of intervals $ (s_1, s_2) $ and $ (s_5, s_6) $, the sign of the switching function does not align with Pontryagin's minimum principle. Based upon Pontryagin's minimum principle, the $ h_{6, bang}^* $ is not locally optimal, which is to be expected since we are ignoring the singular case. In Table 2, we see that the harvesting yield that is generated by the $ h_{6, bang}^* $ control is less than the harvesting yields corresponding to all of the singular control strategies that we observed, but we do find the corresponding yield to be close in value to the other yields.
Our results plead a very strong numerical case that chattering does occur for this problem when the habitat length is set to 10 units and when $ h_{ \mathop{{\rm{max}}}} = 1 $. In this section, we use our numerical results to extrapolate the location of the singular region and the corresponding objective functional value. We take the boundaries of the singular regions of harvesting strategies $ h_{4, sing}^* $, $ h_{6, sing}^* $, $ h_{8, sing}^* $, and $ h_{10, sing}^* $ and apply Aitken's extrapolation to approximate the boundary of the singular region of the locally optimal harvesting strategy. The boundary points of the locally optimal control's singular region give us an idea of where chattering occurs for the locally optimal control since chattering is generated by an inability to concatenate a nonsingular subarc and a singular subarc without prompting infinite oscillations. In Table 3, we have that the singular region of the locally optimal control is near the interval $ (3.535, 6.465) $. Similarly, we take the resulting harvesting yields obtained from applying $ h_{4, sing}^* $, $ h_{6, sing}^* $, $ h_{8, sing}^* $, and $ h_{10, sing}^* $ and apply Aitken's extrapolation to approximate the maximum harvesting yield. In Table 3, the maximum harvesting yield for an infinite number of switches is estimated to be 1.699046237. Thus, the 6-switch yield matches the optimal yield for an infinite number of switches when both are rounded to 7 digits.
In comparing the extrapolating harvesting yield with the yields generated by the other investigated control strategies (see Table 2), we see that the harvesting strategies generated via SPA are formidable alternatives to the chattering control that produces similar yields. While $ h_{8, sing}^* $ and $ h_{10, sing}^* $ produce higher harvesting yields in comparison to the other alternative strategies, it may be difficult to administer such controls when considering the lengths of the non-singular regions bordering the singular arc. In choosing which alternative strategy to implement, we find $ h_{4, sing}^* $, $ h_{6, sing}^* $, and $ h_{6, bang}^* $ to be the best candidates. Of those three candidates, $ h_{6, sing}^* $ produces the largest yield. One reason why $ h_{6, sing}^* $ performs better than $ h_{4, sing}^* $ is that the structure of $ h_{6, sing}^* $ results in having each marine reserve be bordered by regions of maximum harvesting. Observe in Figure 1(c) that the density of the stock is at higher levels along the marine reserves $ (s_1^*, s_2^*) $ and $ (s_5^*, s_6*) $. Naturally, one would want to harvest to the left and right of these reserves as more fish will be swimming in and out of the no-take regions. Control $ h_{6, bang}^* $ also has its marine reserves surrounded by regions of maximum harvesting, but observe in Figure 1(e) that the structure of the bang-bang strategy causes the population density to oscillate between varying levels. That is, the marine reserves are useful in maintaining the fish population, but the size of the maximum harvesting regions to the left and right of the reserve impacts the marine reserve's effectiveness in increasing the population density. However, the population density corresponding to $ h_{6, sing}^* $ given in 1c stays predominately close to 0.5 near $ (s_2^*, s_5^*) $ because of the singular harvesting region being so large.
While the aim of the harvesting problem is to choose a harvesting strategy that maximizes the yield, one could consider other factors in choosing between alternative strategies. For example, one could choose an alternative strategy based upon simplicity of implementation or based upon which of the three strategies is the least costly. One could argue that either $ h_{4, sing}^* $ or $ h_{6, bang}^* $ are the simplest of the three in that $ h_{4, sing}^* $ has fewer switches but $ h_{6, bang}^* $ has no singular region. In the last column of Table 4, we compute $ \int_0^{\ell} wh(x)dx $, where $ w $ is a weight parameter, to measure the cost per unit of effort. Such an integral was used in [14], but Ding and Lenhart also considered another integral term to represent the rate at which the wages paid rises as more labor is employed. Control strategy $ h_{4, sing}^* $ is the cheapest control strategy while $ h_{6, bang}^* $ is the most costly. In Table 4, we observe that $ h_{6, bang}^* $ produces significantly larger marine reserves but also larger regions of maximum harvesting, which results in fewer fish remaining in the habitat. In taking these things into consideration, $ h_{6, bang}^* $ does extensively lower the densities throughout the ecosystem, but $ h_{6, sing}^* $ shows that replacing the central region with a reduced rate of harvesting (through using a singular arc) produces a better yield and is a cheaper strategy to implement. In comparing the number of fish that were not caught, $ h_{4, sing}^* $ allows for more fish to remain in the habitat even though this strategy resulted in a habitat with the lowest total percentage of marine reserves. One could view $ h_{4, sing}^* $ as being ecologically beneficial in that more fish are being preserved within the ecosystem; however, one could view this as a disadvantage in that the yield could be improved if the strategy was adjusted to where the population density could have been lowered (as indicated based upon the results associated with $ h_{6, sing}^* $, $ h_{8, sing}^* $, and $ h_{10, sing}^* $).
5.
Conclusions
We numerically investigated a harvesting problem to illustrate that a SPA can be used to find alternative control strategies with fewer switches in problems that potentially chatter. These alternative solutions can be used to obtain empirical evidence for determining chattering when standard methods like the theorem of conjugation or the junction theorem are inapplicable. Additionally, Aitken's extrapolation can be used to approximate the optimal value and boundary of the singular region to the chattering control, which can be a crude approximation of the locations of the chattering regions.
When selecting the "best" alternative control strategies, one may have to find a trade off between a control with fewer switches and a control that yields a value closer to the optimal value. In using the switch point algorithm to analyze the harvesting problem, we conclude that three of the control forms investigated were the best alternative strategies. The first candidate is a control with four switches and a singular arc located at the center of the habitat, which is the alternative harvesting strategy that Neubert recommended in his analysis of the harvesting problem [12]. The second control form also has a singular region at the center of the habitat but has six switches. The third control form is a bang-bang control with six switches. Of these three, the control with six switches and a singular arc produces the largest yield. In choosing an alternative strategy for a chattering control problem, one could take other features beyond simplicity of implementation into consideration, but those features depend on the context of the problem. For example, in our analysis of the harvesting problem, the alternative harvesting strategy with 4 switches was not only simple to implement but also a cheaper control that produced higher density levels of the uncaught fish within the ecosystem. These features, though, come at the expense of obtaining a lower harvesting yield. If the implementation of the singular arc is up for debate, one can also use SPA to produce an alternative control with no singular arc. For the harvesting problem that we investigated, the resulting harvesting yield corresponding to the bang-bang alternative strategy was lower in comparison to the yields produced by controls with singular regions. In comparing the cost of implementing the control, the bang-bang control was also the most expensive strategy.
A future research direction involves applying SPA to optimal control problems whose state dynamics are systems of partial differential equations. Such extensions are of great importance with regard to constructing more accurate fishery models that depend on both time and space. In our example, we follow Neubert's suggestion of reducing a PDE system to an ODE system by studying the model at a steady state [12], but if the spatial domain is multidimensional, like the one presented in Ding and Lenhart [14], such reductions lead to analysis of an optimal control problem over a PDE. One direction we are considering in applying SPA to PDE problems is to first use standard numerical techniques to obtain a spatially and temporally dependent control. If implementation of such a control is up for debate, we follow the suggestions used in [3] and [1] to construct approximations of the control by making it constant in space and/or time (reducing the dimension of the domain of the control). We then can use SPA to improve these approximations.
Use of AI tools declaration
The authors declare that they have not used Artificial Intelligence (AI) tools in the creation of this article.
Acknowledgments
We would like to thank Prof. Suzanne Lenhart of the University of Tennessee, Department of Mathematics for conversations that inspired this research. William Hager gratefully acknowledges support by the National Science Foundation under grants 1819002 and 2031213, and by Office of Naval Research under grant N00014-22-1-2397.
Conflict of interest
The authors declare there is no conflicts of interest.
Appendix
Verifying the Lipschitz assumptions
In our prior work [34, Chapter 4], we derived the same formula (3.3) for the harvesting problem (2.8). In this derivation, we took advantage of the characterization of the optimal control being piecewise constant which led to a simplification of the dynamics smoothness assumption (1) to where we only needed $ { \boldsymbol f} $, given in $ (2.7) $, to be Lipschitz continuous in $ ({ \boldsymbol y}, h) $ and $ \nabla_y { \boldsymbol f} $ to be Lipschitz continuous in $ { \boldsymbol y} $. The derivation presented in our prior work also holds if we assume $ \nabla_y { \boldsymbol f} $ to be Lipschitz continuous in $ ({ \boldsymbol y}, h) $. Below, we verify that the harvesting problem (2.8) satisfies these assumptions. Note that the derivation presented in [34, Chapter 4] also relied on assumptions corresponding to the invertiblity of multiple submatrices relating to some fundamental matrices. We also assumed that the adjoint equations (2.9) had a solution, which corresponds to a similar assumption that is presented in [23, Theorem 5.1].
Theorem A.1. The function $ { \boldsymbol f} $ defined in (2.7) is Lipschitz continuous.
Proof. Let $ ({ \boldsymbol y}, h) $ and $ (\tilde{ { \boldsymbol y}}, \tilde{h}) $ be arbitrarily chosen. We take the 1-norm of the difference between $ { \boldsymbol f}({ \boldsymbol y}, h) $ and $ { \boldsymbol f}(\tilde{ { \boldsymbol y}}, \tilde{h}) $:
Let $ L = \mathop{{\rm{max}}}\{1+|y_1(x)|+2|\tilde{y}_1(x)|+2|h(x)|\; : \; 0\leq x\leq \ell\} $, then we have the following:
We have that $ {\left\| {\cdot} \right\|}_1 $ and $ {\left\| {\cdot} \right\|}_2 $ are equivalent, where for all $ { \boldsymbol x}\in \mathbb{R}^n $,
By equivalence of norms, we then have
Therefore, $ { \boldsymbol f} $ is Lipschitz continuous with the Lipschitz constant being $ L_1 = 2L $.
□
Theorem A.2. Let $ { \boldsymbol f} $ be the function in (2.7). Then, $ \nabla_y { \boldsymbol f} $ is Lipschitz continuous with respect to $ { \boldsymbol y} $ and $ h $.
Proof. Let $ ({ \boldsymbol y}, h) $ and $ (\tilde{ { \boldsymbol y}}, \tilde{h}) $ be arbitrarily chosen. By definition of $ { \boldsymbol f} $, we have that
Let $ \boldsymbol{A} = \nabla_{y} { \boldsymbol f}({ \boldsymbol y}, h)- \nabla_{y} { \boldsymbol f}(\tilde{ { \boldsymbol y}}, \tilde{h}) $. Then, we have
We then have that
□
DownLoad:
