Optimal management of harvested population at the edge of extinction
OOPTIMAL MANAGEMENT OF HARVESTED POPULATION AT THE EDGE OF EXTINCTION
MICKAËL D. CHEKROUN AND HONGHU LIUAbstract. Optimal control of harvested population at the edge of extinction in an unprotected area, is considered.The underlying population dynamics is governed by a Kolmogorov-Petrovsky-Piskunov equation with a harvestingterm and space-dependent coefficients while the control consists of transporting individuals from a natural reserve.The nonlinear optimal control problem is approximated by means of a Galerkin scheme. Convergence resultabout the optimal controlled solutions and error estimates between the corresponding optimal controls, are derived.For certain parameter regimes, nearly optimal solutions are calculated from a simple logistic ordinary differentialequation (ODE) with a harvesting term, obtained as a Galerkin approximation of the original partial differentialequation (PDE) model. A critical allowable fraction α of the reserve’s population is inferred from the reducedlogistic ODE with a harvesting term. This estimate obtained from the reduced model allows us to distinguish sharplybetween survival and extinction for the full PDE itself, and thus to declare whether a control strategy leads to successor failure for the corresponding rescue operation while ensuring survival in the reserve’s population. In dynamicalterms, this result illustrates that although continuous dependence on the forcing may hold on finite-time intervals,a high sensitivity in the system’s response may occur in the asymptotic time. We believe that this work, by itsgenerality, establishes bridges interesting to explore between optimal control problems of ODEs with a harvestingterm and their PDE counterpart. Contents1. Introduction 22. Optimal control of harvested population 32.1. The KPP model with harvesting term in heterogeneous environment 32.2. Main results from [RC07] 42.3. Optimal control of harvested population at the edge of extinction 62.4. Control constraints from reserve management 83. Optimal control from Galerkin approximations 93.1. Convergence results and error estimates about the optimal control 93.2. PMP solution to the reduced optimal control problem (3.12) 114. Nearly optimal controls from low-dimensional surrogates: Numerical results 144.1. Numerical setup 144.2. Choice of the control constraints 154.3. Control from effective reduced Galerkin systems 164.4. First numerical results 174.5. Effects of the control constraints on the rescue operation 19Acknowledgments 23Appendix A. Galerkin approximations of nonlinear optimal control problems in Hilbert spaces:General convergence results 23Appendix B. Convergence and error estimates results for the optimal control of the KPPH equation 26References 29
Date : July 22, 2020. a r X i v : . [ m a t h . O C ] J u l MICKAËL D. CHEKROUN AND HONGHU LIU
1. IntroductionOverexploitation has led to the extinction of many species [BHTS04]. Traditionally, models of ordinarydifferential equations (ODEs) or difference equations have been used to estimate the maximum sustainable yieldsfrom populations and to perform quantitative analysis of harvesting policies and management strategies [GH89].Ignoring age or stage structures as well as delay mechanisms, which will not be treated by the present paper, theODE models are generally of the type d U d t = F ( U ) − Y ( U ) , (1.1)where U is the population biomass at time t , F ( U ) is the growth function, and Y ( U ) corresponds to the harvestfunction. In these models, the most commonly used growth function is logistic, with F ( U ) = U ( µ − νU ) [BM77, Sch54], where µ > is the intrinsic growth rate of the population and ν > models its susceptibilityto crowding effects.Different harvesting strategies Y ( U ) have been considered in the literature and are used in practical resourcemanagement. A very common one is the constant-yield harvesting strategy, where a constant number of indi-viduals are removed per unit of time: Y ( U ) = δ , with δ a positive constant. This harvesting function naturallyappears when a quota is set on the harvesters [RB99, SFRAS02]. Another frequently used harvesting strategyis the proportional harvesting strategy (also called constant-effort harvesting ), where a constant proportion ofthe population is removed. It leads to a harvesting function of the type Y ( U ) = δU .Reaction-diffusion equations have also been used extensively in modeling the spatiotemporal behavior of aspecies of organism [SK97, CC04, Mur01, Mur07, OL13], and harvesting effects have been analyzed within thisclass of models; see e.g. [OS02, RC07, RC10]. Within this global picture, we focus on a particular class ofmodels, namely the class of diffusive logistic equation in a heterogenous environment [SKT86, CC89, SK97,BH02, BHR05a, BHR05b].More specifically, the harvesting model we consider, described in Sec. 2.1 below, is the Kolmogorov-Petrovsky-Piskunov harvesting (KPPH) model (see (2.1) below) in which the local growth rate of thepopulation and the limit effects to crowding are spatial-dependent coefficients. The harvest function is aquasi-constant-yield harvesting term depending on a threshold below which harvesting is progressively aban-doned. This model is analyzed in [CR06, RC07] and the main results about its steady state analysis and itsdynamical behavior, are recalled in Sec. 2.2 below for the reader’s convenience. In particular, the asymptoticbehavior of the KPPH model was characterized as a function of the harvesting intensity δ ; see (2.1). It wasproved in [RC07] that if δ is smaller than a critical value δ ∗ , the population density converges to a “significant”state, which is everywhere above a small threshold. On the other hand, it was shown in [RC07] that if δ is largerthan another critical value δ , which is slightly above δ ∗ in practice, the population density eventually settlesdown to a “remnant” state, everywhere below the same small threshold. The population can be considered asextinct in this case.Within this context, we assume that given an unprotected area Ω in which the population (of let’s say somemammals) evolves according to a KPPH equation, is under threat of extinction due to some illegal hunting thatsuperimposes to some allowed harvesting, causing thereof a rise in the harvesting intensity such that δ endsto jump above δ . Given a natural reserve (cid:98) Ω in which the same species evolves, we address in this articlethe problem of saving the population under threat by releasing in a controlled fashion new individuals fromthe reserve, into the unprotected area, while avoiding to exert too much pressure on the reserve that would bedetrimental on it. Our control strategy is seriously constrained in time as the operation is assumed to startwhen a (large) fraction of the original population has been decimated and its course to extinction is, thus, wellengaged.The resulting optimal control problem differs fundamentally from other optimal control problems concernedwith harvesting models that appeared in previous studies which dealt mainly with the search for optimalharvesting strategies; see e.g. [Neu03, KS08, FLVP12, BXY14]. Instead, we aim at designing optimal rescue PTIMAL MANAGEMENT OF HARVESTED POPULATION AT THE EDGE OF EXTINCTION 3 strategies while taking into account plausible factors that may arise in real-life problems, from a reservemanagement perspective.The originality of this work does not only rely on its problem formulation, but also on its proposed solution.We indeed address the obtention of optimal solutions by means of rigorous Galerkin approximations of theunderlying KPPH models. To do so, we rely on the recent mathematical framework introduced in [CKL17]which allows for deriving convergence results and error estimates from Galerkin approximations of a broad classof nonlinear optimal control problems in infinite dimension. Of course, the idea of designing approximations ofoptimal solutions from finite-dimensional approximations is not new, and a great deal of works have addressedthis question in various contexts [Ded10, GK11, HV05, IK08, TV09] as well as based on various basis functions,possibly empirical [FHLW12, HK98, HK00, Rav00]. However, rigorous convergence analysis from finite-dimensional Galerkin approximations do not seem to have been much addressed for the optimal control ofnonlinear problems, and in that sense [CKL17] provides some useful elements.In that respect, we show that the framework of [CKL17] allows us to derive error estimates between theoptimal control of the KPPH model considered in this article and the optimal control built from any ofits Galerkin approximations. These error estimates are supplemented by convergence of the correspondingcontrolled solutions towards the optimally controlled solution of the full partial differential equation (PDE)model. These convergence and error estimates results are summarized in Theorem 3.1 below which is provedin Appendix B based on the theoretical apparatus from [CKL17], recalled in Appendix A.The finite-dimensional approximations are here simply obtained from projection onto the eigenmodes ofa natural underlying spectral problem (see (2.3) below). A standard Pontryagin Maximum Principle (PMP)approach [PBGM64, Kir12] is then used to find extremals of the finite-dimensional optimal control problems[BC03], which turn out to approximate, in our case, the optimal solutions of the full problem; see Sec. 3.2.To illustrate our theoretical framework and to favor reproducibility of the results, the numerical results ofSec. 4 are presented within the context of one-dimensional environments and for various level of fragmentationof the population’s habitat. For the parameter regime considered, we show that from a one-dimensional, ordinarydifferential equation (ODE) approximation, nearly optimal solutions can be derived. The corresponding Galerkinapproximation of the KPPH model, reducing to a simple logistic ODE with a harvesting term , allows us inturn to provide great insights about the optimal control problem of the PDE model. The main result, from anecological perspective, is indeed expressed in terms of a fraction α of the reserve that is allowed for exploitation.This fraction, related to e.g. management policy of the reserve, is caping the amount of individuals transportedfrom the reserve into the unprotected area. We show that a critical fraction α can be inferred from the reducedlogistic ODE with a harvesting term. This estimate obtained from the reduced equation is indeed shown to bevery useful as it allows us to distinguish sharply between survival and extinction for the full PDE itself, and thusto declare whether a control strategy leads to success or failure for the corresponding rescue operation whileensuring survival in the reserve’s population. In terms of optimal control, this result translates into the existenceof controls very close to each other, one of which leading to a significant survival of the population while theother leading to its extinction. In dynamical terms, this result illustrates that although continuous dependenceon the forcing may hold on finite-time intervals, a high sensitivity in the system’s response may occur in theasymptotic time. 2. Optimal control of harvested population2.1. The KPP model with harvesting term in heterogeneous environment.
The partial differential equation(PDE) and boundary conditions underlying the optimal control problem considered hereafter are described as Whose solution contains most of the PDE solution’s energy.
MICKAËL D. CHEKROUN AND HONGHU LIU follows. Let Ω be a smooth bounded and connected domain of R d ( d ≥ ). We consider ∂ t y = D ∇ y + µ ( x ) y − ν ( x ) y − δρ (cid:15) ( y ) , ( t, x ) ∈ (0 , ∞ ) × Ω , (2.1a) ∂y∂ n = 0 , ( t, x ) ∈ [0 , ∞ ) × ∂ Ω . (2.1b)Here n denotes the outward unit normal to the boundary ∂ Ω .This equation differs from the classical Fisher equation [Fis37] (also known as the Kolmogorov-Petrovsky-Piskunov (KPP) equation [KPP37]), by its spatially-dependent coefficients, µ ( x ) and ν ( x ) , as well as itsharvesting term δρ (cid:15) ( y ) . When δ = 0 in (2.1), the model reduces to the Shigesada-Kawasaki-Teramoto modeldescribed in [SKT86]; see also [BH02, BHR05a, BHR05b, SK97]. Such a problem fits with general speciesassessment and management problematics considered for instance in [BHTS04, BM77, GH89, Neu03, RC10,Sch91,Fah03,KS08]. Hereafter, Eq. (2.1a) will be referred to as the Kolmogorov-Petrovsky-Piskunov harvesting(KPPH) equation.The unknown function y = y ( t, x ) denotes the population density at time t and space position x . Thecoefficient µ represents the intrinsic growth rate of the population, which is assume to be a measurable functionof x in L ∞ (Ω) . The spatial dependence of µ is introduced to account for the possible impact of environmentalheterogeneity [RC07]. The function ν ( x ) > (also in L ∞ (Ω) ) represents the susceptibility to crowding effectsand is interpreted as an intraspecific competition term. Regions with higher values of µ ( x ) and lower values of ν ( x ) are qualified as being more favorable , while, on the other hand, regions with lower µ ( x ) and higher ν ( x ) values are considered as being less favorable or, equivalently, more hostile .The harvesting function ρ ε satisfies ρ (cid:15) ∈ C ( R ) , ρ (cid:48) (cid:15) ≥ , ρ (cid:15) ( s ) = 0 , ∀ s ≤ , and ρ (cid:15) ( s ) = 1 , ∀ s ≥ (cid:15), (2.2)where (cid:15) is a nonnegative parameter, taken to be sufficiently small in a sense made precise in Sec. 2.2.The term δρ (cid:15) ( y ) with δ ≥ , corresponds to a quasi-constant-yield harvesting term. Indeed, for such aharvesting function, the yield is constant in time whenever y ≥ (cid:15) , while it depends on the population densitywhen y < (cid:15) . Note that the function ρ (cid:15) ensures the nonnegativity of the solutions to (2.1); see [RC07].From a biological viewpoint, ε corresponds to a threshold below which harvesting is progressively abandoned.Considering constant-yield harvesting functions without this threshold value would be unrealistic since it wouldeventually lead to harvest on zero-populations.2.2. Main results from [RC07].
The problem (2.1) has been analyzed in [RC07]. Using sub- and supersolutionmethods and the characterization of the first eigenvalue of the linearized elliptic operator at the trivial solution,the authors obtained existence and nonexistence results as well as results on the number of stationary solutions;see also [CR06].The asymptotic behavior of the evolution equation was in particular characterized as a function of theharvesting intensity δ . In [RC07] it was proved that if δ is smaller than a critical value δ ∗ , the population densityconverges to a “significant” state, which is everywhere above a small threshold (depending on (cid:15) ) while if δ islarger than δ (another threshold that is bigger than but close to δ ∗ ), the population density y ( t, x ) converges toa “remnant” state, everywhere below this small threshold. Theorem 2.1 below summarizes the main asymptoticand existence results from [RC07].To formulate these results we recall some tools used in [RC07] and in the formulation of our optimal controlproblem below. For this purpose, we first consider the eigenvalue problem associated with the linearization of(2.1) at the trivial solution: − D ∇ φ − µ ( x ) φ = λφ, x ∈ Ω , (2.3a) ∂φ∂ n = 0 , x ∈ ∂ Ω . (2.3b) PTIMAL MANAGEMENT OF HARVESTED POPULATION AT THE EDGE OF EXTINCTION 5
Let λ denote the first eigenvalue and φ its corresponding positive eigenfunction, unique when normalized. Inparticular φ satisfies φ ( x ) > ∀ x ∈ Ω , and (cid:107) φ (cid:107) ∞ = 1 . (2.4)The above normalization is possible since the first eigenfunction has a fixed sign in Ω as consequence of theKrein-Rutman theorem; see e.g. [Ama76]. We introduce next φ = min x ∈ Ω φ , (2.5)Note that φ > since φ does not vanish on the boundary ∂ Ω ; see [Ama76]. Note also that φ ≤ thanks to(2.4).Recall from [RC07, Definition 2.5] that, a stationary solution of (2.1), p δ , is called a significant solution if min x ∈ Ω p δ ≥ (cid:15)φ . (2.6)On the other hand, a stationary solution p δ is called remnant if max x ∈ Ω p δ < (cid:15)φ . (2.7)With these tools in hand, it was derived in [RC07] the following analytic formulas allowing for estimatingthe critical harvesting intensity δ ∗ leading to decline of the population towards a remnant steady state (i.e. closeto extinction): δ = λ φ ν (1 + φ ) , δ = λ ν , (2.8)see [RC07, (2.15)] with α = 1 since h ( x ) therein is identically equal to 1 here. Finally, ν = max x ∈ Ω ν ( x ) and ν = min x ∈ Ω ν ( x ) .We are now in position to summarize the main theoretical results from [RC07] into the following theorem. Theorem 2.1. Steady states [RC07]. If λ < , then there exists a threshold δ ∗ ≥ such that ( i ) if δ ≤ δ ∗ , there exists at least one positive stationary significant solution of (2.1) , whereas ( ii ) if δ > δ ∗ , there is no positive stationary significant solution of (2.1) .Moreover, ( iii ) if λ < and δ ≤ δ , there exists a positive stationary significant solution p δ of (2.1) such that p δ ≥ − λ φ ν (1 + φ ) . (2.9) ( iv ) if λ < and δ > δ , the only possible positive bounded stationary solutions of (2.1) are remnant. Asymptotic behavior [RC07].
Assume that the initial datum y (0 , x ) for Eq. (2.1) is taken to be p , theunique nonnegative steady state of Eq. (2.1) when δ = 0 [BHR05a]. Then, the solution y ( t, x ) is non-increasingin t and its asymptotic behavior is characterized as follows: ( v ) Population resilience: If δ ≤ δ ∗ , the solution y ( t, · ) converges to p δ uniformly in Ω as t goes to infinite,where p δ in (2.9) is the maximal stationary significant solution of (2.1) . ( vi ) Population extinction: If δ > δ , the solution y ( t, · ) converges uniformly in Ω to a remnant solutionof (2.1) . For a proof of the above results, see [RC07, Theorem 2.6] for ( i ) and ( ii ) , [RC07, Theorem 2.10] for ( iii ) and ( iv ) , and [RC07, Theorem 2.11] for ( v ) and ( vi ) . In particular, ( iii ) and ( iv ) show that δ and δ definedin (2.8) provide easily computable bounds for the maximum allowable harvesting intensity δ ∗ . These boundshave been shown to be quite sharp for a broad range of habitat configurations [RC07, Figure 4].Numerical simulations from [RC07] based on the analytical estimates of δ ∗ , strongly supported that environ-mental fragmentation of the habitat has a significant impact on the maximum sustainable yield associated with MICKAËL D. CHEKROUN AND HONGHU LIU δ ∗ . Essentially, the more fragmented is the habitat, the more the population is susceptible to decline towards aremnant state under harvesting pressure.2.3.
Optimal control of harvested population at the edge of extinction.
We describe hereafter the optimalcontrol problem we aim at solving for (2.1). The motivation is as follows. Assume that the population whoseevolution is governed by (2.1) is under threat of extinction due to e.g. some illegal practices of harvesting.According to Theorem 2.1 this situation is encountered when e.g. δ > δ . Indeed, without additional externalintervention, the population will eventually settle down to a remnant state and thus will be close to extinction.Our goal is to prevent such a situation by exerting a certain control on the population by releasing new individualsfrom the same species coming e.g. from a natural reserve, into the unprotected area Ω , once the population inthat area drops below a certain warning threshold P c .The natural reserve is assumed to be limited in resources, and thus the rescue operation is itself underconstraints. We assume that the reserve contains the same species that we are aiming at saving in the unprotectedarea and that the reserve’s population is at an equilibrium, i.e. that it has reached a steady state. Because noharvesting is exerted in the reserve, we assume that the population dynamics in the reserve is governed by (2.1)with δ = 0 , namely the KPP equation in an heterogeneous environment [SKT86]. In the general case, thedomain (cid:98) Ω of the reserve is different from that of the unprotected area, although it can share similar features suchas a similar area and proportion of the population’s habitat (in which µ > ). More precise considerations aboutthe reserve will be formulated in the numerical section. For the moment, we specify our goal in general termsuseful to frame our optimal control problem. Thus, we are aiming at transporting a fraction of the populationfrom the reserve (cid:98) Ω to the unprotected area Ω . The question is to determine which fraction to use in order tosave the population at extinction threat in Ω while not only avoiding to remove all the population in the reservebut also not causing extinction within this same reserve. Warning time t = τ . We assume also that the population in Ω before harvesting was exerted ( δ = 0 ),had reached the stationary state p and that the total population size can be monitored in Ω . Based on thismonitoring, we assume that the authorities declare the population at threat when the following warning threshold P c is reached P c = β (cid:90) Ω y (0 , x ) d x, < β < , (2.10)where y (0 , x ) = p is the initial density of the population at time t = 0 .Recall that we assume δ > δ in (2.1) for the unprotected area Ω . As mentioned above, a situation for which δ would be greater than δ may arise due for instance to illegal hunting practiced without respect of some quotaand a lack of reliable monitoring by the authorities.In such a case, Theorem 2.1-( vi ) guarantees that y ( t, x ) converges uniformly towards a remnant solution to(2.1) and because y ( t, x ) is continuous in time and space, we conclude that there exists a time instant, t = τ ,such that (cid:82) Ω y ( τ, x ) d x = P c . Note that we want to avoid to reach extinction so P c and thus β in (2.10) must bechosen such that P c > (cid:15) | Ω | φ , (2.11)where | Ω | denotes the d -dimensional volume of Ω . In other words one wishes to trigger an alert signal whenthe situation corresponds still to a significant population size according to (2.6). The time instant τ serves usthus to trigger an alert signal from which one then starts to release new individuals from the natural reserve intothe unprotected area. It will be called the warning time instant. Extinction time t = T and optimal control problem. The goal consists then of restoring, through anappropriate management plan exploiting the reserve, the population to a safe dynamics away from extinction.For that purpose and based on the insights gained from Theorem 2.1, we aim at driving the population distributionto be as close as possible to a significant target population, chosen here to be the significant steady state p δ (cid:48) of(2.1) for some δ (cid:48) ≤ δ where δ is defined in (2.8). Recall that δ is a critical threshold that can be estimated PTIMAL MANAGEMENT OF HARVESTED POPULATION AT THE EDGE OF EXTINCTION 7 in practice and that guarantees that the KPPH equation possesses a significant steady state (thus away fromextinction) for any δ ≤ δ .We are limited by time in our action. We want indeed to save the population from extinction that will occurat time T , corresponding to first time instant at which the maximum of the solution equals (cid:15)/φ if nothing isdone . Denoting by u ( t, x ) the number of new individuals brought from the reserve at time t and position x ,we consider then the cost functional J ( y, u ) = 12 (cid:90) Tτ (cid:12)(cid:12) y ( t ) − p δ (cid:48) (cid:12)(cid:12) L (Ω) d t + κ (cid:90) Tτ (cid:12)(cid:12) u ( t ) (cid:12)(cid:12) L (Ω) d t. (2.12)Here y ( t ) denotes the solution in L (Ω) of the KPPH equation forced by u ( t, x ) and such that y ( t ) = y at t = τ , and κ > . The first term in the right-hand side (RHS) of (2.12) is to enforce closeness to p δ (cid:48) , whilethe second term is an energy-type term related to the effort of bringing the individuals from the reserve into theunprotected area.More precisely, we are aiming at addressing the following type of optimal control problems min u ∈ L ( τ,T ; L (Ω)) J ( y, u ) (2.13a)where y ( t, x ) is the solution of ∂ t y = D ∇ y + µ ( x ) y − ν ( x ) y − δρ (cid:15) ( y ) + u ( t, x ) , ( t, x ) ∈ [ τ, T ] × Ω , (2.13b) ∂y∂ n = 0 , ( t, x ) ∈ [ τ, T ] × ∂ Ω , (2.13c)with y ( t, x ) = y ( x ) at t = τ . (2.13d)As explained below, we restrict to controls lying within a subset of L ( τ, T ; L (Ω)) , by introducing constraintson the control. To describe these constraints we enter into more specificities about the model setting. Throughoutthis article, we focus on a particular case of growth rate, µ , defined to be µ ( x ) = mχ Λ , m > , (2.14)where χ Λ denotes the characteristic function of a (possibly disconnected) subdomain Λ of Ω : χ Λ ( x ) = (cid:40) , if x ∈ Λ , , otherwise . (2.15)Such a spatial dependence of the coefficient µ emphasizes that the population reproduces only when theindividuals are within the subdomain Λ , while they may spread outside of Λ due to the diffusion term D ∇ y .The choice of the control u in (2.13) is designed to act on this subdomain to enhance chances of naturalgrowth. For numerical applications we will consider piecewise constant controls in space while allowing forfluctuations in time. For that purpose, the domain Λ is decomposed into mutually disjoints subdomains suchthat Λ = Λ ∪ Λ ∪ · · · ∪ Λ K . (2.16)Then, by introducing ϕ j = 1 (cid:112) | Λ j | χ Λ j , ≤ j ≤ K, (2.17)where | Λ j | denotes the d -dimensional volume of Λ j , we consider the following set of admissible controls: U ad = (cid:26) ( t, x ) (cid:55)→ K (cid:88) j =1 Γ j ( t ) ϕ j ( x ) : Γ j ∈ L ([ τ, T ] , R ) , ≤ Γ j ( t ) ≤ C j for a.e. t ∈ [ τ, T ] (cid:27) , (2.18) Because we assume δ > δ , we know that such a time instant exists as the population will become remnant; see Theorem 2.1-(vi). MICKAËL D. CHEKROUN AND HONGHU LIU where the C j ’s are positive constants which, as explained below, are imposed by exploitation policy of theprotected reserve. Denoting by H = L (Ω) , the optimal control problem (2.13) becomes, for this set ofadmissible controls, min J ( y, u ) subject to ( y, u ) ∈ L ( τ, T ; H ) × U ad solving (2.13b)-(2.13d). (2.19)Due to the definition of U ad , this optimal control problem is under constraints on the control . In the next sectionwe will show that this optimal control problem can be addressed efficiently by means of finite-dimensionalGalerkin approximations. For the moment we discuss how the constraints C j in (2.18) arise from reservemanagement considerations. Remark 2.1.
Note that for any given nonnegative continuous initial datum y in (2.13d) , the existence of anoptimal pair ( y ( t ; u ∗ ) , u ∗ ) for the optimal control problem (2.19) can be ensured by using classical argumentssuch as found in [Trö10, Theorem 5.7] in combination with the maximum principle for (2.13b) - (2.13d) .Indeed, (2.13b) - (2.13d) can be put into the form [Trö10, Eq. (5.8)] by taking d ( x, t, y ) = − µ ( x ) y + ν ( x ) y + δρ (cid:15) ( y ) and b ≡ therein. The assumptions collected into [Trö10, Assumption 5.6] are satisfied here except therequirement that d y ( x, t, y ) ≥ . This latter condition is used in [Trö10] to ensure the existence of solutionsto [Trö10, Eq. (5.8)] and also to guarantee the uniform boundedness of the solutions given by [Trö10, (5.10)];see [Trö10, Assumption 5.4 and Theorem 5.5] . However, for (2.13b) - (2.13d) , the existence of solutions andthe desired uniform bound also hold as a consequence of the maximum principle here; see Appendix B. Control constraints from reserve management.
We assume that the population density in the reserve isat its steady state (cid:98) p , at the warning time τ . The total population in the reserve is thus given by M tot = (cid:90) (cid:98) Ω (cid:98) p ( x ) d x. (2.20)A naive strategy, referred to as unsupervised below, consists of assuming Γ j ( t ) to be set to a constant value C j to determine. Then the total (unsupervised) population, M T , released in the unprotected area Ω from thereserve, between times t = τ and t = T , is M T = (cid:16) K (cid:88) j =1 C j (cid:113) | Λ j | (cid:17) ( T − τ ) . (2.21)In practice, it is reasonable to assume that an exploitation policy holds for the natural reserve independently ofthe rescue operation that is aimed to take place for Ω . As a result, only a fraction α ( < α < ) of the totalpopulation M tot in the reserve is allowed to be used for populating another area such as Ω .An optimal control planning u ∗ ( t, x ) solving the optimal control problem (2.19) will be declared as efficientif it leads to an efficiency ratio E = (cid:82) Tτ (cid:82) Ω u ∗ d x d tM tot < M T M tot . (2.22)Assuming that the unsupervised strategy would exploit the allowable fraction α , we have (cid:16) K (cid:88) j =1 C j (cid:113) | Λ j | (cid:17) ( T − τ ) = αM tot . (2.23)Finally, we will assume in what follows that the d -dimensional volume | Λ j | are the same for different j ’s, andthat C j = C j (cid:48) = C for j (cid:54) = j (cid:48) . This way, Eq. (2.23) gives the value of the constraints C j in U ad . By doing so wehave thus an impartial way to compare the solution u ∗ of the optimal control problem (2.19), to that obtainedfrom a naive, unsupervised strategy consisting of taking Γ j ( t ) ≡ C . We have however to keep in mind that in Note that the set U ad will be endowed with the induced topology from that of L ( τ, T ; L (Ω)) . In [Trö10, Assumption 5.4], the condition d y ( x, t, y ) ≥ is stated as d being monotonically increasing with respect to y . PTIMAL MANAGEMENT OF HARVESTED POPULATION AT THE EDGE OF EXTINCTION 9 this design operatus of the constraints, the set of admissible controls, U ad depends on α , since the C j do. As aresult, any optimal solution u ∗ to (2.19) does also depend on the allowable fraction α from the reserve.From an exploitation viewpoint, we might ask what is the best fraction α to use in order to maximize theproductivity of new individuals in the unprotected area, by the whole operation. An natural metric to answerthis question consists of calculating the population ratio P R = (cid:82) Ω y ( T, u ∗ ) d x (cid:82) Ω p δ (cid:48) d x . (2.24)The ratio P R allows us to assess which fraction of the targeted total population has been obtained (when drivenby the optimal control planning u ∗ ) when one reaches what would have been the extinction time T if no actionwould have been taken. A natural question is then to analyze how the choice of α impacts P R and the efficiencyratio E , the goal being, at this stage of the discussion, to minimize E while maximizing P R as much as possible.This later aspect will be discussed in Sec. 4 dealing with the numerical results.Also we want to make sure that the excision of population from the reserve does not lead to extinctionthere, neither. For that, we can rely here again on the theoretical understanding from the harvesting problem asrecalled in Sec. 2.2. Let us assume that the removal of individuals from the reserve follows also a harvestinglaw in (cid:98) Ω of the form (cid:80) Kj =1 C j ϕ j ( x ) ρ (cid:15) (ˆ y ) , where ˆ y denotes the population density in the reserve. If one wants u ∗ = (cid:80) Kj =1 C j ϕ j , one must have ˆ y ( t, x ) ≥ (cid:15)/φ for t in [ τ, T ] and x in the reserve’s domain (cid:98) Ω to ensure ˆ y tobe significant (see (2.6)), imposing thus a constraint on the population on the reserve.On the other hand, let ˆ δ denotes the critical threshold ensuring resilience of the population within the reservethis time; see Theorem 2.1 again. A simple comparison argument shows that as soon as K (cid:88) j =1 C j ϕ j ( x ) ≤ ˆ δ (cid:18) K (cid:88) j =1 ϕ j ( x ) (cid:19) , x ∈ (cid:98) Ω , (2.25)then a significant steady state is ensured to exist, favoring resilience of the reserve’s population.The inequality (2.25) is equivalent to For ≤ j ≤ K, C j ≤ ˆ δ , (2.26)since the (cid:98) Λ j are mutually disjoint. Thus the rescue operation for the unprotected area must comply with theconstraint (2.26), and also from what precedes, with ˆ y ( t, x ) ≥ (cid:15)/φ for t ∈ [ τ, T ] and x in the reserve’s domain (cid:98) Ω . In particular, (2.26) imposes a restriction on the population size M T excised from the reserve to satisfy M T ≤ ˆ δ ( (cid:80) Kj =1 (cid:113) | (cid:98) Λ j | )( T − τ ) . Such aspects regarding the reserve management will be also discussed inSec. 4 below. For the moment we focus in the next section on the mathematical aspects of solving the optimalcontrol problem (2.19) via Galerkin approximations, for a given set of admissible controls U ad (and thus for agiven allowable fraction α ).3. Optimal control from Galerkin approximations3.1. Convergence results and error estimates about the optimal control.
In this section, we present con-vergence results and error estimates regarding Galerkin approximations of the optimal control problem (2.19)constructed from eigenprojections. The synthesis of nearly optimal controls based on these Galerkin approxi-mations is then provided in Section 3.2, following a standard Pontryagin Maximum Principle (PMP) approach.Denote the spectral elements of the eigenvalue problem (2.3) by { ( λ j , e j ) : j ∈ N } , where the eigenvalues λ j ’s are arranged in an increasing order, and the eigenfunctions e j ’s are normalized such that (cid:107) e j (cid:107) L (Ω) = 1 . We assume that the reserve’s habitat, (cid:98) Λ , is decomposed also into K mutually disjoints subdomains each of same size than itscorresponding counterpart Λ j in Ω . Denote also the N -dimensional Galerkin approximation of the controlled state y by y N ( t, x ) = N (cid:88) i =1 ξ i ( t ) e i ( x ) . (3.1)Recall that the control u is written as u ( t, x ) = K (cid:88) i =1 Γ i ( t ) ϕ i ( x ) . (3.2)Note that y N in (3.1) depends on the initial datum and the control u driving Eq. (2.13b) (see (3.3) below).Throughout this article, the initial datum for the Galerkin approximation is taken to be Π N y , where Π N denotes the projector onto the subspace spanned by the first N eigenmodes solving (2.3). Dependence of y N on Π N y or u will be made apparent depending on the context.The Galerkin approximation of (2.13b)-(2.13c) reads then: d ξ i d t = − λ i ξ i + N (cid:88) j,k =1 B ijk ξ j ξ k − δ (cid:68) ρ (cid:15) (cid:16) N (cid:88) j =1 ξ j e j (cid:17) , e i (cid:69) + (cid:104) u ( t, · ) , e i (cid:105) , t ∈ [ τ, T ] , (3.3)where i = 1 , . . . , N , and B ijk = −(cid:104) ν ( · ) e j e k , e i (cid:105) = − (cid:90) Ω ν ( x ) e i ( x ) e j ( x ) e k ( x ) d x. (3.4)Introducing a K × N matrix M , whose elements are defined by M ij = (cid:104) ϕ i , e j (cid:105) = (cid:90) Ω ϕ i ( x ) e j ( x ) d x, (3.5)we can rewrite (3.3) as d ξ i d t = − λ i ξ i + N (cid:88) j,k =1 B ijk ξ j ξ k − δ (cid:68) ρ (cid:15) (cid:16) N (cid:88) j =1 ξ j e j (cid:17) , e i (cid:69) + K (cid:88) j =1 M ji Γ j ( t ) , t ∈ [ τ, T ] . (3.6)The cost functional J N associated with the N -dimensional Galerkin approximation (3.6) is J N ( y N , u ) = 12 (cid:90) Tτ | y N ( t ; Π N y , u ) − Π N p δ (cid:48) | L (Ω) d t + κ (cid:90) Tτ | u ( t ) | L (Ω) d t. (3.7)By introducing ξ = ( ξ , . . . , ξ N ) T , and Γ = (Γ , . . . , Γ K ) T , (3.8)we rewrite the cost functional J N given by (3.7) as: J N ( ξ , Γ ) = 12 (cid:90) Tτ | ξ ( t ) − P | d t + κ (cid:90) Tτ | Γ ( t ) | d t, (3.9)where P = ( P , . . . , P N ) T with P i = (cid:104) p δ (cid:48) , e i (cid:105) , i = 1 , . . . , N. (3.10)In connection to the set of admissible controls U ad defined by (2.18), we introduce V ad = (cid:26) t (cid:55)→ (Γ ( t ) , . . . , Γ K ( t )) T : Γ j ∈ L ([ τ, T ] , R ) , ≤ Γ j ( t ) ≤ C j for a.e. t ∈ [ τ, T ] (cid:27) . (3.11)The Galerkin approximation of the optimal control problem (2.19) is thus given by min J N ( ξ , Γ ) defined in (3.9) subject to ( ξ , Γ ) ∈ L ( τ, T ; R N ) × V ad solvingthe N -dimensional Galerkin system (3.6) with ξ ( τ ) = Π N y . (3.12) PTIMAL MANAGEMENT OF HARVESTED POPULATION AT THE EDGE OF EXTINCTION 11
We have then the following convergence results and error estimates, linking the optimal control problem (2.19)to its Galerkin approximation (3.12).
Theorem 3.1. Error estimates about the optimal control.
Let us consider the optimal control problem (2.19) along with its Galerkin approximation (3.12) . Assume the initial datum y in (2.13d) is strictly positive.Then, there exist C > such that the optimal control u ∗ for (2.19) and the optimal control u ∗ N for the reducedproblem (3.12) admit the following error estimate: (cid:90) Tτ (cid:12)(cid:12) u ∗ ( t ) − u ∗ N ( t ) (cid:12)(cid:12) L (Ω) d t ≤ C (cid:104) √ T − τ + T − τ (cid:105) (cid:18) (cid:88) j =1 (cid:18) (cid:90) T ε jN ( t, u ∗ , u ∗ N )d t (cid:19) (cid:19) , (3.13) with sup t ∈ [ τ,T ] ε jN ( t, u ∗ , u ∗ N ) N →∞ −−−−→ , for j = 1 , . (3.14) Convergence results.
Furthermore, the solution y N ( t ; Π N y , u ) of (3.6) converges uniformly to the solution y ( t ; y , u ) of (2.13b) - (2.13d) : lim N →∞ sup u ∈U ad sup t ∈ [ τ,T ] (cid:12)(cid:12) y N ( t ; Π N y , u ) − y ( t ; y , u ) (cid:12)(cid:12) L (Ω) = 0 . (3.15)In Appendix B we show that these convergence results and error estimates are the consequence of generalresults of [CKL17] dealing with the Galerkin approximations of nonlinear optimal control problems in Hilbertspaces and recalled in Appendix A, for the reader’s convenience.More precisely, Theorem 3.1 above is a consequence of Theorem A.1 and Theorem A.2 presented inAppendix A. Theorems A.1 and A.2 deal with the optimal control of a broad class of nonlinear evolutionaryequations containing the KPPH equation considered in this article. The required assumptions of Theorem A.1and Theorem A.2 are verified in Appendix B for the optimal control problems (3.12) and (2.19).In particular, the analysis shows that the constant C in (3.13) depends on the local Lipschitz constant (in L (Ω) ) of G ( y ) = (cid:82) Ω | y − p δ (cid:48) | d x in a neighborhood of the origin, and also on the growth rate of J in someappropriate norms of its arguments as given by (A.14) in Appendix A. The precise expression of ε N in (3.14)is found in (A.15).3.2. PMP solution to the reduced optimal control problem (3.12) . Thanks to Theorem 3.1, it is sufficient tosolve the reduced optimal control problem (3.12) to obtain nearly optimal solutions as soon as N is sufficientlylarge. We will see that for the KPPH model considered here, one can already obtain nearly optimal solutions with N = 1 for certain specifications of the model such as the domain size, and choice of the model’s parameters.One can now use standard techniques from finite-dimensional optimal control theory to solve the abovereduced optimal control problem (3.12); see [BH75, BC03, Kir12]. Here, we follow an indirect approach byrelying on the Pontryagin Maximum Principle (PMP); see e.g. [PBGM64, Kir12].In that respect, we introduce the following Hamiltonian associated with the reduced optimal control problem(3.12): H G ( ξ , z , Γ ) = 12 | ξ − P | + κ | Γ ( t ) | + z T f ( ξ , Γ ) , (3.16)where z = ( z , . . . , z N ) T denotes the costate associated with the state ξ , and f is the Galerkin vector fieldwhose i -th component is given by the RHS of (3.6). That is, f i ( ξ , Γ ) = − λ i ξ i + N (cid:88) j,k =1 B ijk ξ j ξ k − δ (cid:68) ρ (cid:15) (cid:16) N (cid:88) j =1 ξ j e j (cid:17) , e i (cid:69) + K (cid:88) j =1 M ji Γ j ( t ) , i = 1 , . . . , N. (3.17)Let ( ξ ∗ , Γ ∗ ) ∈ L ( τ, T ; R N ) × V ad be an optimal pair for the reduced optimal control problem (3.12), and denote by z ∗ the costate associatedwith the state ξ ∗ . It follows from the PMP that the triplet ( ξ ∗ , z ∗ , Γ ∗ ) must satisfy the following constrainedHamiltonian system (see e.g. [Kir12, Section 5.3]): d ξ ∗ d t = ∇ z H G ( ξ ∗ ( t ) , z ∗ ( t ) , Γ ∗ ( t ))d z ∗ d t = −∇ ξ H G ( ξ ∗ ( t ) , z ∗ ( t ) , Γ ∗ ( t )) , ( Hamiltonian system for ( ξ ∗ , z ∗ )) (3.18a) H G ( ξ ∗ ( t ) , z ∗ ( t ) , Γ ∗ ( t )) ≤ H G ( ξ ∗ ( t ) , z ∗ ( t ) , Γ ( t )) , ∀ Γ ∈ V ad , ( optimality condition ) (3.18b) z ∗ ( T ) = 0 , ( terminal condition ) (3.18c)where (3.18a) and (3.18b) hold for all t in ( τ, T ) . Here ∇ x stands for the gradient operator in the x -direction. Remark 3.1.
The optimality condition states that an optimal control must minimize the Hamiltonian H G . Ingeneral it is a necessary condition and not a sufficient condition. In the general case, there may be controls thatsatisfy (3.18b) but that are not optimal controls. Yet the PMP may delineate a nonempty class of candidates.A triplet ( ξ ∗ , z ∗ , Γ ∗ ) solution to (3.18) is called an extremal. Extremal solutions play an important role inoptimal control theory; see [BC03]. Sufficient conditions for extremal to provide optimal controls can be foundin [HSV95]. For the problem at hand, since the cost functional J N in (3.9) is quadratic in Γ and the dependenceon the control is linear for the control system (3.6) , it is known that Γ ∗ obtained from such an extremal isactually the unique optimal control of the optimal control problem (3.12) ; see e.g. [Kir12, Sec. 5.3]. From (3.18), we derive now, the explicit formula of the optimal control Γ ∗ based on the optimal costate z ∗ .To do so, we first remark that − ∂H G ( ξ , z , Γ ) ∂ξ i = − ( ξ i − P i ) − N (cid:88) j =1 z j ∂f j ( ξ , Γ ) ∂ξ i . (3.19)Now since ∂f j ( ξ , Γ ) ∂ξ i = − λ i δ ij + N (cid:88) k =1 ( B jik + B jki ) ξ k − δ (cid:68) ρ (cid:48) (cid:15) (cid:16) N (cid:88) k =1 ξ k e k (cid:17) e i , e j (cid:69) , (3.20)where δ ij denotes the Kronecker delta, we get − ∂H G ( ξ , z , Γ ) ∂ξ i = − ( ξ i − P i ) − N (cid:88) j =1 z j (cid:16) − λ i δ ij + N (cid:88) k =1 ( B jik + B jki ) ξ k − δ (cid:68) ρ (cid:48) (cid:15) (cid:16) N (cid:88) k =1 ξ k e k (cid:17) e i , e j (cid:69)(cid:17) . (3.21)Since (3.21) does not depend on Γ , we denote − ∂H G ( ξ , z , Γ ) /∂ξ i by a function g i ( ξ , z ) . The costate equationbecomes then d z ∗ / d t = g ( ξ ∗ , z ∗ ) . Thus the Hamiltonian system (3.18a) together with boundary condition(3.18c) becomes d ξ ∗ d t = f ( ξ ∗ , Γ ∗ ) , t ∈ ( τ, T ) , d z ∗ d t = g ( ξ ∗ , z ∗ ) , t ∈ ( τ, T ) , ξ ∗ ( τ ) = Π N y , z ∗ ( T ) = 0 , (3.22)where the components of f and g are defined in the RHSs of (3.17) and (3.21), respectively. The condition ξ ∗ ( τ ) = Π N y is obtained by projecting the initial condition y ( τ, x ) = y in (2.13d). PTIMAL MANAGEMENT OF HARVESTED POPULATION AT THE EDGE OF EXTINCTION 13
We show next how Γ ∗ depends on the costate z ∗ . To do so, we first remark that by using M defined in (3.5),and the expressions of H G and the f i (see (3.16) and (3.17)), the optimality condition (3.18b) becomes κ | Γ ∗ ( t ) | + ( z ∗ ( t )) T M T Γ ∗ ( t ) ≤ κ | Γ ( t ) | + ( z ∗ ( t )) T M T Γ ( t ) , ∀ Γ ∈ V ad , t ∈ ( τ, T ) . (3.23)Then by introducing w ∗ = M z ∗ , (3.24)we get κ | Γ ( t ) | + ( z ∗ ( t )) T M T Γ ( t ) = κ | Γ ( t ) | + ( w ∗ ( t )) T Γ ( t )= κ K (cid:88) j =1 (cid:16) Γ j ( t ) + 1 κ w ∗ j ( t ) (cid:17) − κ K (cid:88) j =1 w ∗ j ( t ) . (3.25)As a consequence, the optimality condition (3.23) becomes K (cid:88) j =1 (cid:16) Γ ∗ j ( t ) + 1 κ w ∗ j ( t ) (cid:17) ≤ K (cid:88) j =1 (cid:16) Γ j ( t ) + 1 κ w ∗ j ( t ) (cid:17) , ∀ Γ ∈ V ad , t ∈ ( τ, T ) . (3.26)Recalling the control constraints ≤ Γ j ≤ C j resulting from (3.11), one obtains that (3.26) holds if and onlyif, for all j = 1 , . . . , K , Γ ∗ j ( t ) = h j ( w ∗ ( t )) def = , if − w ∗ j ( t ) < , − κ w ∗ j ( t ) , if ≤ − w ∗ j ( t ) ≤ κC j ,C j , if − w ∗ j ( t ) > κC j , ∀ t ∈ ( τ, T ) . (3.27)Now by substituting the expression of Γ ∗ thus obtained and by using (3.24), the problem (3.22) reduces to aboundary value problem (BVP) in the variables ξ ∗ and z ∗ . The synthesis of an optimal control Γ ∗ boils downthus to solving the following BVP: d ξ ∗ d t = f ( ξ ∗ , h ( M z ∗ )) , t ∈ ( τ, T ) , d z ∗ d t = g ( ξ ∗ , z ∗ ) , t ∈ ( τ, T ) , ξ ∗ ( τ ) = Π N y , z ∗ ( T ) = 0 , (3.28)where components of the function h are defined in (3.27).The optimal control u ∗ N to (3.12) in U ad (see (2.18)) is finally given by u ∗ N ( t, x ) = K (cid:88) i =1 Γ ∗ i ( t ) ϕ i ( x ) , (3.29)where Γ ∗ is given by (3.27) with w ∗ therein obtained as w ∗ = M z ∗ for z ∗ solving (3.28). Remark 3.2.
Note that, if the constraints ≤ Γ j ≤ C j on the admissible control are removed in V ad in (3.11) ,the optimality condition (3.18b) can be equivalently written as (see again [Kir12, Section 5.3]): ∇ Γ H G ( ξ ∗ , z ∗ , Γ ∗ ) = 0 . (3.30) The control law for the constrained case is given by (3.27) . In contrast, one obtains from (3.30) the followingcontrol law for the unconstrained case: Γ ∗ = − κ M z ∗ , (3.31) with M given by (3.5) . Note that in practice (3.27) reduces to (3.31) when the C j are sufficiently large.
4. Nearly optimal controls from low-dimensional surrogates: Numerical results4.1.
Numerical setup.
To illustrate our theoretical framework, and to simplify the reproducibility of the numer-ical results shown below, our numerical experiments take place in the case of one-dimensional environments.In that context, the unprotected area is a bounded, connected domain given by an interval, namely
Ω = (0 , (cid:96) ) ,with (cid:96) > . To simplify the analysis, we choose the domain (cid:98) Ω of the reserve to be also given by an interval ofsame length.We set D = 1 , (cid:96) = 1 , and ν = 0 . in the corresponding KPPH model (2.1). To account for the possibleeffects of heterogeneity of the habitat, we consider two cases of subdomain, Λ , appearing in the definition of µ given by (2.14): Case I: Λ = [1 / , , Case II:
Λ = [0 , / ∪ [3 / , . (4.1)Note that | Λ | = 0 . for both cases, but Case I corresponds to an habitat more aggregated that Case II. In eachcase we set m = 2 in (2.14).Regarding the set U ad of admissible controls (see (2.18)), we divide Λ into K = 8 segments of equal length,such that each segment is of length / .Following Sec. 2.3, we assume that for each case the population evolving in the unprotected area according to(2.1) is under extinction threat because δ = δ (1 + f ) , with f chosen to be equal to . . One aims at restoringthe population to a safe dynamics leading towards a significant steady state p δ (cid:48) appearing in the cost functional J given by (2.12) for some δ (cid:48) < δ . For that purpose we choose δ (cid:48) = δ (1 − f ) , also with f = 0 . . Recallthat δ and δ are critical harvesting thresholds defined in (2.8); see Sec. 2.2.The harvest function ρ (cid:15) in (2.2) is explicitly given here for (cid:15) = 0 . by ρ (cid:15) ( x ) = , if x ≥ (cid:15), . π ( x − . (cid:15) ) /(cid:15) ) + 0 . , if < x < (cid:15), , otherwise . (4.2)Recall that the optimal control problem (2.19) is carried out over the time interval [ τ, T ] . Following Sec. 2.3,given a solution y of the KPPH model (2.1) (emanating from p ), the warning time τ is chosen according to (cid:82) (cid:96) y ( τ, x ) d x = P c with P c defined in (2.10) for some < β < . We choose β = 1 / in all the numericalexperiments. In particular P c satisfies (2.11) for Cases I and II. In each case, the extinction time T correspondsto the first time instant at which the maximum (over (0 , (cid:96) ) ) of the solution y ( t, x ) equals (cid:15)/φ . Note that y used in the optimal control problem (2.19) is y ( τ, · ) (see (2.13d)).The goal is to drive, over the time window [ τ, T ] , the population distribution to be as close as possible to p δ (cid:48) , while minimizing the cost functional given by (2.12) for controls u lying in the admissible set U ad given by(2.18). The determination of the corresponding control constraints C j in (2.18) is discussed in Sec. 4.2 below.The rest of model’s parameters for the two choices of subdomain Λ given by (4.1) (Cases I and II) are listed inTable 1, rounded to the nearest thousandth. Note that δ and δ defined in (2.8) depend in each case on the firsteigenmode e solving (2.3) (associated with the first eigenvalue λ ) shown here in Fig. 1.In what follows, the KPPH model is solved using the Matlab solver pdepe for Ω = (0 , with δx = 10 − and δt = 10 − . The corresponding spectral problem (2.3) is solved using finite difference with δx = 2 − . Weuse a small mesh grid size here to reach good accuracy given the discontinuity of the coefficient µ . In bothparameter cases, the dominant eigenvalue λ is negative (see Table 1), corresponding to an unstable eigenmode φ . The rest of the spectrum corresponds here to stable eigenmodes. We recall, that once δ (cid:48) ≤ δ , there exists indeed a significant steady state, p δ (cid:48) , due to Theorem 2.1-( iii ). Note that e is normalized in L (Ω) whereas φ appearing in (2.8) is normalized in L ∞ (Ω) . The two are related according to φ = e / max( e ) . PTIMAL MANAGEMENT OF HARVESTED POPULATION AT THE EDGE OF EXTINCTION 15
Figure 1.
First eigenmode of the spectral problem (2.3) for
Ω = (0 , shown here for Λ = [1 / , (Case I)and Λ = [0 , / ∪ [3 / , (Case II). Table 1.
Parameter values in (2.13) λ δ δ δ (cid:48) δ τ T Case I -1.083 1.443 1.465 1.298 1.612 12.544 14.009Case II -1.021 1.301 1.302 1.171 1.432 14.023 15.5824.2.
Choice of the control constraints.
We assume that the reserve’s domain size | (cid:98) Ω | is equal to (cid:96) and that thepopulation’s habitat in the reserve (cid:98) Λ satisfies | (cid:98) Λ | = | Λ | . To simplify the analysis, we also assume that the KPPmodel’s parameters in the reserve (cid:98) Ω are the same as those of the KPPH equation in the unprotected area. Underthese working assumptions, the total population in the reserve is given by M tot = (cid:90) (cid:98) Ω (cid:98) p ( x ) d x, (4.3)where (cid:98) p is the steady state of the KPP model (under Neumann condition) in (cid:98) Ω . In our numerical experiments, M tot = 5 . for Case I and M tot = 5 . for Case II.As pointed out in Sec. 2.4, the control planning u ( t, x ) of transporting individuals from the natural reserveto the unprotected area is subordinated to the following reserve management factors:A) The choice of the population fraction α allowed to be transported from the reserve.B) The population displaced from the reserved must be such that the population distribution in the reservesatisfies ˆ y ( t, x ) ≥ (cid:15)/φ for t in [ τ, T ] and x in (cid:98) Ω .C) For ≤ j ≤ K, the constant C j must satisfy C j ≤ ˆ δ , where ˆ δ is the critical harvesting thresholdfavoring population persistence for the reserve.As explained in Sec. 2.4, the constraints C j in (2.18) are determined according to (2.23). We choose C j = C for every j , i.e. the constants are all the same across the subdomains Λ j . Then (2.23) allows us to find C givenby C = αM tot T − τ ) , (4.4) -12 -10 -8 -6 -4 -2 -20 -15 -10 -5 Figure 2.
Energy decomposition of the target population density for
Λ = [1 / , (Case I) and Λ = [0 , / ∪ [3 / , (Case II) shown in semilog-scale. The target population density (shown in the upper insets) corresponds ineach case to the significant steady states p δ (cid:48) of the KPPH equation (2.1) for δ = δ (cid:48) , with δ (cid:48) as indicated in Table 1.These steady states are highly correlated with the corresponding first eigenmodes shown in Fig. 1, explaining thatthese eigenmodes capture most of the energy contained in these steady states. because K = 8 and | Λ j | = 1 / . Note although the habitat dependence is not directly apparent in (4.4), theconstant C depends actually on Λ through the coefficient µ impacting the dynamics and thus the extinction time T and the warning time τ . Given our numerical setup, for the constant C chosen according to (4.4) then therequirement C) above is satisfied as long as < α ≤ α c with α c = 0 . for Case I and α c = 0 . for CaseII. Requirement B) will be assessed a posteriori, after the optimal control problem is solved; see Sec. 4.4.4.3. Control from effective reduced Galerkin systems.
To control the KPPH equation we approximate itby a Galerkin truncation of minimal dimension, namely we choose N = 1 . The reason is that for the modelparameters considered here, the energy in the KPPH’s solution is almost fully captured by the first eigenmode,the latter capturing more than 99.99% of the energy contained in the target population density p δ (cid:48) in each case;see Fig. 2. The insets of this figure show indeed that these target population densities are highly correlated withthe corresponding first eigenmodes shown in Fig. 1.For N = 1 , the Galerkin approximation (3.6) reduces to the following scalar logistic equation with harvestingand forcing terms: d X d t = X ( a − bX ) − δg (cid:15) ( X ) + K (cid:88) j =1 d j Γ j ( t ) , t ∈ [ τ, T ] , K = 8 , (4.5)where a = − λ > , b = − B (see (3.4)) and the d j = M j ( N = 1 ) are constants given by (3.5), thusobtained here by the inner product between the eigenmode e and the characteristic functions of the subdomains Λ j . The harvesting term g (cid:15) is given by g (cid:15) ( X ) = (cid:40) if X ≥ (cid:15) min e , (cid:68) ρ (cid:15) (cid:16) Xe (cid:17) , e (cid:69) if X < (cid:15) min e , (4.6)because the first eigenmode e is positive. The harvest function g (cid:15) satisfies the conditions in (2.2), and thus theGalerkin projection preserves the global structure of the original KPPH equation. PTIMAL MANAGEMENT OF HARVESTED POPULATION AT THE EDGE OF EXTINCTION 17
Due to the high-energy content captured by the first-mode amplitude, X ( t ) , we expect that solving the optimalcontrol (3.12) associated with (4.5) (i.e. for N = 1 ) should provide a nearly optimal solution for the originaloptimal control problem (2.19) with Ω = (0 , (cid:96) ) . The next section explores this intuition in more details.Based on the analysis of Sec. 3.2, in what follows, the optimal control u ∗ N ( t, x ) (for N = 1 ) is obtainedaccording to (3.29). The optimal control Γ ∗ j therein is obtained from (3.27), namely Γ ∗ j ( t ) = h j ( w ∗ ( t )) , where w ∗ = z ∗ M with z ∗ solving here the following BVP over [ τ, T ] (from (3.28)) d X ∗ d t = X ∗ ( a − bX ∗ ) − δg (cid:15) ( X ∗ ) + K (cid:88) j =1 d j h j ( z ∗ M ) , d z ∗ d t = − ( X ∗ − (cid:104) p δ (cid:48) , e (cid:105) ) − z ∗ ( a − bX ∗ ) − δk (cid:15) ( X ∗ ) ,X ∗ ( τ ) = (cid:104) y , e (cid:105) , z ∗ ( T ) = 0 , (4.7)with k (cid:15) ( X ∗ ) = (cid:68) ρ (cid:48) (cid:15) (cid:16) X ∗ e (cid:17) e , e (cid:69) . (4.8)Note that here M given by (3.5) is a column vector and the costate z ∗ is a scalar variable.In what follows, higher N -dimensional Galerkin truncation (with N > ) are also used to control the KPPHequation. The use of such higher-dimensional Galerkin approximations is to benchmark the one-dimensionalapproximation. The corresponding optimal control u ∗ N is, for each N , obtained by following the PMP approachas described in Sec. 3.2. Table 2. Coefficients of the 1D Galerkin approximation (4.5) a b δ d d d d d d d d Case I 1.083 0.202 1.612 0.252 0.259 0.265 0.270 0.274 0.277 0.279 0.280Case II 1.021 0.200 1.432 0.257 0.257 0.255 0.252 0.252 0.255 0.257 0.2574.4.
First numerical results.
We focus here on two particular choices of the allowable fraction α ( α = 0 . and α = 0 . ) of individuals transported from the reserve’s population to the unprotected area Ω , for each caseof habitat Λ in (4.1). Recall that the discrimination between the two cases of habitats considered in Sec. 4.1 isaimed at assessing the possible effects of fragmentation of the habitat on the rescue operation.Below, when one writes J ( y, u ∗ N ) (resp. J N ( ξ ∗ , Γ ∗ N ) ) it corresponds to the value of J defined in (2.12)(resp. J N defined in (3.9)) when the solution y (resp. ξ ∗ ) to the controlled KPPH model (2.13b)-(2.13d) (to the N -dimensional Galerkin approximation (3.6) of the controlled KPPH model) is driven by u ∗ N (resp. Γ ∗ N ).For both cases of habitat (Cases I and II) and for both α -values ( α = 0 . and α = 0 . ), we have computedthe relative error in the cost values J between the optimal solutions obtained from the 1D and the 10D Galerkinapproximations (see Table 3), namely error = | J ( y, u ∗ ) − J ( y, u ∗ ) | J ( y, u ∗ ) × . (4.9)We found that in all cases, this relative error is almost negligible, bounded by . × − % .Furthermore, the convergence of the cost value J N ( ξ ∗ , Γ ∗ N ) is achieved quickly as the first few digits in thecost value has already converged when N is increased; see Table 4. These numerical results confirm thus theintuition that due to the high-energy content carried out by the first eigenmode, a 1D Galerkin approximation issufficient to obtain nearly optimal solutions .Based on these results, we are in a comfortable position to discuss the numerical KPPH solution y ( t, x ; u ∗ N ) obtained when driven by the optimal control u ∗ N itself obtained according to (3.29) (for N = 1 ) and the BVP Table 3.
Cost value J ( y, u ∗ N ) ( α = 0 . , N = 1 ) ( α = 0 . , N = 10 ) ( α = 0 . , N = 1 ) ( α = 0 . , N = 10 ) Case I 4.168706 4.168700 3.203817 3.203754Case II 3.952086 3.952086 2.996004 2.995997Table 4.
Cost value J N ( ξ ∗ , Γ ∗ N ) N = 1 N = 2 N = 3 N = 4 N = 5 N = 10 Case I, α = 0 . α = 0 . α = 0 . α = 0 . α considered. The corresponding panels (a) (or (d)) show the optimal control u ∗ N ( t, x ) (for N = 1 )for different time instants while panels (b) (or (e)) show its time evolution, after integration over Ω = (0 , (cid:96) ) ,denoted by (cid:104) u ∗ N ( t ) (cid:105) .Recall that given a value of α , the constant C is determined from (4.4) which determines in turn the controlconstraints in the set U ad of admissible controls, as the C j therein are chosen to be equal to C ; see Sec. 4.2. Theimpact of α is as follows. The smaller the allowable fraction α is, the smaller C is and the more constrained isthe optimal control problem (2.19). This is visible on the temporal evolution of (cid:104) u ∗ N ( t ) (cid:105) which is set to C overa longer interval for α = 0 . than for α = 0 . , before decreasing to zero; compare Panels (b) and (e) of Fig. 4with those of Fig. 3. We conduct a more detailed analysis on the dependence on α in Sec. 4.5 below. For themoment, we discuss the effects of fragmentation of the habitat as the latter has been shown to play an importantrole in the population resilience to external perturbations [RC07].Table 5. The efficiency ratio E defined by (2.22) α = 0 . α = 0 . Case I 0.0912 0.1920Case II 0.0918 0.2002The efficiency ratio E defined in (2.22) serves us to compare the effects of different configurations of thehabitat. Recall that the smaller is this ratio the better it is in terms of exploitation of the reserve’s populationbut not necessarily for saving the population from extinction in the unprotected area. For that latter aspect oneshould have the population ratio P R (defined in (2.24)) as large as possible while keeping in mind to respect therequirements A)-C) formulated in Sec. 4.2 in order to preserve the reserve’s population from extinction.In terms of efficiency ratio, Table 5 indicates that, for a given allowable fraction α of the reserve, the habitatthat is more homogeneous (Case I) is slightly more advantageous than the less homogeneous one (Case II). Interms of population ratio P R defined by (2.24), it is Case II that is slightly more advantageous, for instance P R = 0 . for Case II vs P R = 0 . for Case I, when α = 0 . . This means that for Case II, . ofthe targeted population size is reached at time t = T , and . for Case I. Thus, the differences across thehabitat fragmentation are somewhat minor. To the opposite, the difference in terms of allowable fraction α of PTIMAL MANAGEMENT OF HARVESTED POPULATION AT THE EDGE OF EXTINCTION 19
Figure 3.
Panel (a) : Optimal control u ∗ N ( t, x ) (Case I, α = 0 . ) given by (3.29) with N = 1 . Here Γ ∗ N isobtained from (3.27), where w ∗ = z ∗ M with z ∗ solving here the BVP (4.7). Panel (b) : Time-dependence of (cid:104) u ∗ N ( t ) (cid:105) after space integration over (0 , (cid:96) ) . The first time instant at which (cid:104) u ∗ N ( t ) (cid:105) > corresponds to t = τ . Thenext time instant at which (cid:104) u ∗ N ( t ) (cid:105) = 0 corresponds to t = T . Panel (c) : Controlled KPPH solution y ( t, x ; u ∗ N ) to(2.13b)-(2.13d) when driven by u ∗ N . Panel (d) same as Panel (a), Panel (e) same as Panel (b), and Panel (f) same asPanel (c), for Case II. the reserve seem important as only of about (in both cases) of the targeted population size is reached attime t = T when α = 0 . . The next section points out the main factors responsible of such a marked difference.4.5. Effects of the control constraints on the rescue operation.
When α = 0 . , we observe in Panels (c) and(f) of Fig. 4 that the final profile y ( T ; u ∗ N ) of the controlled solution satisfies y ( T ; u ∗ N ) < y ( y corresponds tothe blue curve in Panels (c) and (f)), irrespectively of the degree of fragmentation of Λ . As time goes beyond t = T and the harvesting intensity δ is set back to δ (cid:48) in the KPPH equation while the control is abandoned , thepopulation eventually settles down to a remnant state. The rescue operation has thus failed.On the contrary, we observe in Panels (c) and (f) of Fig. 3 that for α = 0 . , y ( T ; u ∗ N ) > y , and one gets ineither Case I or II, closer than for α = 0 . to the target population, i.e. the significant steady state p δ (cid:48) . As timegoes beyond t = T and δ is set to δ (cid:48) , the population converges towards this significant steady state (not shown,but see below). The rescue operation is a success.What makes the difference so pronounced between α = 0 . and α = 0 . ? Remember that due to (4.4) andour protocol for choosing our set U ad of admissible controls, the allowable fraction α impacts the constraint C arising in U ad . The smaller α , the smaller C and the more constrained becomes our control planning. At thesame time it should be noted that when α is increased further towards , the constraint C becomes too large sothat this constraint is never activated and the optimal control problem (2.19) behaves as unconstrained. Thus Assuming we get rid of illegal harvesting by this time instant, while still allowing some harvesting respecting the quota δ (cid:48) < δ . It is easy to convince oneself of this fact by remarking that Γ ∗ given by (3.27) reduces to (3.31) when the C j are sufficiently large;see Remark 3.2. Figure 4.
Same as Fig. 3 but for α = 0 . . Figure 5.
Panel (a) : Population ratio P R (Case I). Panel (b):
Efficiency ratio E (Case I). The value α corresponds to the value of α above which the control bounds are not activated. The value α corresponds to thevalue of α below which the population in Ω evolves towards a remnant state when the control is abandoned and δ is set to δ (cid:48) after t = T ; see (4.11) below. there exists a critical value of alpha, say α , above which the control bounds are not activated in the set U ad ofadmissible controls. For Case I we found α = 0 . and α = 0 . for Case II.As a consequence, for α > α , only the same unconstrained optimal solution is obtained and the efficiencyratio E saturates to a constant value. Panel (b) of Fig. 5 shows the dependence of E in terms of α , and inparticular its saturation for α > α . Only Case I is shown here as E for Case II behaves almost identically.The dependence on α of the population ratio P R is shown in Panel (a) of Fig. 5. We observe that the P R - and PTIMAL MANAGEMENT OF HARVESTED POPULATION AT THE EDGE OF EXTINCTION 21
13 14 15 16 17 18 19 204.44.54.64.74.84.955.15.25.3 14 15 16 17 18 19 204.44.54.64.74.84.955.15.25.3
14 14.5 15 15.500.511.5
Figure 6.
Panel (a) : Population size (cid:104) ˆ y ( t ) (cid:105) in the reserve (Case I, α = α ). The bracket indicates integration overspace. The dynamics for the population in the natural reserve is governed by (2.1) forced by − u ∗ N , corresponding tothe amount of individual transported from the reserve to the unprotected area. The inset shows (cid:104) u ∗ N ( t ) (cid:105) over a timewindow slightly larger than [ τ, T ] . The first vertical dashed line indicates t = τ , while the second one, indicates t = T . Panel (b) : Same for Case II, α = α . E -curves are highly correlated, both increasing and saturating for α > α . Thus, “there is no free lunch,” andone cannot minimize E while maximizing P R . The choice of α is however a determining factor in the successof the rescue operation as pointed above. In all the cases, the excision of population to be transported from thereserve to the unprotected area can be absorbed by the reserve’s population. Figure 6 shows indeed that even forthe more demanding “excision pressure” on the population’s reserve (corresponding to the larger E achieved for α = α and beyond), the population in the reserve recovers it original steady state, independently on the habitat’sfragmentation.The critical role of α on the survival of the population in the unprotected after application to optimal planningcan gain great insights from the reduced system (4.5) and the dynamical properties of the KPPH equation asrecalled in Sec. 2.2. Let us assume that once the optimal control u ∗ N has been obtained from the BVP (4.7)for a given α , the harvesting intensity δ is set back to δ (cid:48) in the KPPH equation, for t > T , and the controlis abandoned. The idea is that by assuming that the illegal harvesting that was putting the population underextinction threat has been stopped starting from t = T , one aims at anticipating in terms of α (thus on theconstraint put on the control) whether the rescue operation that would take place over the time window [ τ, T ] ,would be a success or not.To do so, let us consider the following modification of (4.5) d X d t = X ( a − bX ) − δ (cid:48) g (cid:15) ( X ) , (4.10)which is nothing else than the projection of the KPPH equation (with δ = δ (cid:48) ) onto the first eigenmode e . Let X (cid:15) denotes the smallest positive steady state of (4.10). If X (0) < X (cid:15) , then X ( t ) converges to a remnant statewhereas if X (0) > X (cid:15) , it converges towards a significant steady state.Marking the dependence of u ∗ N on α , one defines α to be the smallest value of α for which the followingcondition holds true (cid:104) y ( T ; u ∗ N ( α )) , e (cid:105) > X (cid:15) . (4.11)
15 20 25 30 35 40 4500.511.522.533.5 12.5 13 13.5 1400.20.40.60.811.21.41.61.8
Figure 7.
Panel (a) : Population size (cid:104) y ( t ) (cid:105) (after space integration) in the unprotected area Ω for Case I. The firstvertical dashed line indicates t = τ , while the second one, indicates t = T . For t < τ , the dynamics is governedby (2.1) with δ set to its corresponding value given in Table 1 (Case I). The initial datum at t = 0 is taken to be p , the steady state of (2.1) for δ = 0 . The dynamics over the time window [ τ, T ] , is governed by (2.13b)–(2.13d)driven by u ∗ N ( t ) whose space integration is shown in Panel (b). For t > T , the dynamics is governed by (2.1) butfor δ = δ (cid:48) (Table 1) and with no control. Panel (b) : Space integration, (cid:104) u ∗ N ( t ) (cid:105) , of the optimal control u ∗ N obtainedfrom (3.29) (with N = 1 ) and the BVP (4.7). Each curve is shown here as α is varied according to the same valuesused for Panel (a) and for a time window slightly larger than [ τ, T ] . For Case I one finds α = 0 . using a mesh grid of size . to discretize the range [0 , of α -values. Becauseof the high-energy content carried by e , we expect that if (4.11) is satisfied, then the solution, y ( t ; u ∗ N ( α )) ,to the KPPH equation for δ = δ (cid:48) , converges to a significant state as t tends to ∞ . By showing the populationsize behavior (after space integration) (cid:104) y ( t ) (cid:105) as time evolves, Panel (a) of Fig. 7 shows that this is exactly whathappens for α ≥ α while the population converges towards a remnant steady state for α < α . One mightthus wish to get α as close as possible to α from above in order to lower down E while still ensuring success.However one has to keep in mind that the same figure reveals that the boundary between success or failure ofthe rescue operation is very narrow, as α is getting too close to α . So in practice a hard constraint C on theindividuals of the reserve should not correspond to an α -value too close to α , as C may e.g. not be exactlyrespected during the displacement operation.Panel (b) of Fig. 7 which shows the controls (applied over [ τ, T ] ) corresponding to the solutions shown inPanel (a), illustrates this statement. As shown in this panel, the control leading to a significant survival (redcurve) is indeed very close to that leading to extinction (blue curve). We are here in presence of an interestingphenomenon, namely, continuous dependence on the forcing may hold on finite-time intervals, but a highsensitivity in the system’s response may take place in the asymptotic time.We shall emphasize also that unlike what was reported for α = 0 . , one may have y ( T ; u ∗ N ) > y whilethe population still evolves towards a remnant steady state, when δ is set to δ (cid:48) after t = T and the control isabandoned. Figure 7 shows for instance that the population size in Ω on the rise from y up to t = T (due tocontrol), drops eventually to a remnant state in an asymptotic time for e.g. α = 0 . and α = 0 . .Finally, we stress that the numerical results reported are not limited to the particular numerical setupconsidered here. For instance if D is further reduced, then more modes become unstable, say p , and thedimension of an efficient reduced system must be at least p . Here N = p = 1 for D = 1 . By setting D = 0 . PTIMAL MANAGEMENT OF HARVESTED POPULATION AT THE EDGE OF EXTINCTION 23 in Case II while keeping the other parameters as in Sec. 4.1, the two dominant eigenvalues of the spectralproblem (2.3) are negative corresponding to two unstable modes e and e capturing most of the energy of thetarget population density, unlike the energy decomposition shown in Fig. 2. By choosing N = 2 , the BVP(4.7) becomes (3.28) for N = 2 . The corresponding efficient Galerkin approximation is then given by (3.3) for N = 2 from which a critical α separating survival from extinction can also be determined from an analogueof (4.11) in which X (cid:15) (resp. the projection onto e ) is replaced by the steady state of smaller norm (resp. theprojection onto the space spanned by e and e ). Such a remark about the inflation of the efficient reduceddimension holds as D approaches zero. For instance for D = 0 . in Case II, there exists four unstable modescapturing most of the energy, constraining thus the efficient reduced dimension to be N = 4 . The approachpresented here extends also to more realistic two-dimensional domains in which the fragmentation effects of thehabitat may play a more important role than for the one-dimensional domains; see [RC07].AcknowledgmentsThis work has been partially supported by the European Research Council (ERC) under the European Union’sHorizon 2020 research and innovation program (grant agreement No. 810370 (MDC)), by the Office of NavalResearch (ONR) Multidisciplinary University Research Initiative (MURI) grant N00014-20-1-2023 (MDC),and by the US National Science Foundation grant DMS-1616450 (HL).Appendix A. Galerkin approximations of nonlinear optimal control problems in Hilbert spaces:General convergence resultsIn this appendix, we summarize from [CKL17] key convergence results and error estimates obtained forGalerkin approximations of nonlinear optimal control problems in Hilbert spaces. As we will see in AppendixB below, the material of this section allows us to link precisely the optimal control problem (2.19) to its Galerkinapproximation (3.12), in terms of, both, error estimates about the optimal controls, and (strong) convergenceabout the controlled solutions as summarized in Theorem 3.1 in the Main Text.In that respect, we tailor here the conditions of applications of the abstract results of [CKL17] to the particularcase of Galerkin systems from eigenprojections as considered in this article. The interested reader is referredto [CKL17] for more general cases, and also for convergence results concerning the value functions. Applicationsto control in feedback form can be found in [CKL18].In the following, we consider finite-dimensional approximations of the following initial-value problem (IVP): d y d t = Ly + F ( y ) + C ( u ( t )) , t ∈ (0 , T ] ,y (0) = y , (A.1)where the unknown y evolves in a separable Hilbert space H , L : D ( L ) ⊂ H → H is a linear operator withdomain D ( L ) , F : H → H denotes the nonlinearity, and the initial datum y belongs to H . The time-dependentforcing u lives in a possibly different separable Hilbert space V ; the (possibly nonlinear) mapping C : V → H is assumed to be such that C (0) = 0 . Other assumptions regarding C will be made below when needed.We assume the following conditions for the linear operator L . (H0) The linear operator L : D ( L ) ⊂ H → H is the infinitesimal generator of a C -semigroup of boundedlinear operators T ( t ) on H . (H1) L is self-adjoint with compact resolvent. Recall that under assumption (H0) , the domain D ( L ) of L is dense in H and that L is a closed operator;see [Paz83, Cor. 2.5, p. 5].For the admissible set of controls, we assume that (H2) Given some q ≥ , the set of admissible controls U ad is the subset of L q (0 , T ; V ) constituted bymeasurable functions that take values in U , a bounded subset of the Hilbert space V . In other words, U ad = { f ∈ L q (0 , T ; V ) : f ( s ) ∈ U for a.e. s ∈ [0 , T ] } , q ≥ . (A.2)The set U ad will be endowed with the induced topology from that of L q (0 , T ; V ) .Let u be in U ad given by (A.2), a mild solution to (A.1) over [0 , T ] is a function y in C ([0 , T ] , H ) such that y ( t ) = T ( t ) y + (cid:90) t T ( t − s ) F ( y ( s )) d s + (cid:90) t T ( t − s ) C ( u ( s )) d s, t ∈ [0 , T ] . (A.3)In what follows we will often denote by t (cid:55)→ y ( t ; y , u ) a mild solution to (A.1).Since L is assumed to be self-adjoint with compact resolvent, it follows from spectral theory of self-adjointcompact operator that the eigenfunctions of L form an orthonormal basis of H ; see e.g. [Bré10]. We denote theeigenpairs of L by { ( β k , e k ) : k ∈ N } . For each N ≥ , let H N be the N -dimensional subspace of H spannedby the first N eigenfunctions of L : H N = span { e k : k = 1 , . . . , N } . (A.4)Denote also by Π N : H → H N the associated orthogonal projector. Note that H N ⊂ D ( L ) , ∀ N ≥ . (A.5)The corresponding Galerkin approximation of (A.1) associated with H N is then given by: d y N d t = L N y N + Π N F ( y N ) + Π N C ( u ( t )) , t ∈ [0 , T ] ,y N (0) = Π N y , y ∈ H , (A.6)where L N = Π N L Π N : H → H N . (A.7)In particular, the domain D ( L N ) of L N is H , because of (A.5).Throughout this article, a mapping f : W → W between two Banach spaces, W and W , is said to belocally Lipschitz if for any ball B r ⊂ W with radius r > centered at the origin, there exists a constant Lip( f | B r ) > such that (cid:107) f ( y ) − f ( y ) (cid:107) W ≤ Lip( f | B r ) (cid:107) y − y (cid:107) W , ∀ y , y ∈ B r . (A.8)We will make also use of the following assumptions. (H3) The mapping F : H → H is locally Lipschitz in the sense given in (A.8) . (H4) Let U ad be given by (A.2) . For each T > and ( y , u ) in H × U ad , the problem (A.1) admits a uniquemild solution y ( · ; y , u ) in C ([0 , T ] , H ) ; and for each N ≥ , its Galerkin approximation (A.6) admitsa unique solution y N ( · ; Π N y , u ) in C ([0 , T ] , H ) . Moreover, there exists a constant C = C ( T, y ) suchthat (cid:107) y ( t ; y , u ) (cid:107) H ≤ C , ∀ t ∈ [0 , T ] , u ∈ U ad , (A.9a) (cid:107) y N ( t ; Π N y , u ) (cid:107) H ≤ C , ∀ t ∈ [0 , T ] , N ∈ N , u ∈ U ad . (A.9b) Remark A.1.
Note that in applications, (H4) is typically satisfied via a priori estimates. Indeed, let u be in U ad given by (A.2) . Then uniform bounds such as in (A.9) are guaranteed if e.g. an a priori estimate of the followingtype holds for the IVPs (A.1) and (A.6) : sup t ∈ [0 ,T ] (cid:107) y ( t ; y , u ) (cid:107) H ≤ α ( (cid:107) y (cid:107) H + (cid:107) u (cid:107) L q (0 ,T ; V ) ) + β, α > , β ≥ . (A.10) See e.g. [CH98, Tem97] for such a priori bounds for nonlinear partial differential equations. Such bounds canalso be derived for nonlinear systems of delay differential equations (DDEs); see [CGLW16, CKL18]. For theKPPH problem (2.13b) - (2.13c) , the desired estimate (A.9a) is derived for nonnegative initial data by using themaximum principle together with energy estimates; see Appendix B. In particular as we will see, Theorem A.1and Theorem A.2 below, both apply to this context when y ≥ in (H4) . PTIMAL MANAGEMENT OF HARVESTED POPULATION AT THE EDGE OF EXTINCTION 25
We introduce next the cost functional, J : H × U ad → R + , associated with the IVP (A.1): J ( y , u ) = (cid:90) T [ G ( y ( s ; y , u )) + E ( u ( s ))] d s, y ∈ H , (A.11)where G : H → R + and E : V → R + are assumed to be continuous , and G is assumed to satisfy furthermorethe condition: G is locally Lipschitz in the sense of (A.8) . (A.12)The associated optimal control problem then writes min J ( y , u ) subject to ( y, u ) ∈ L (0 , T ; H ) × U ad solves (A.1) with y (0) = y ∈ H . ( P )The cost functional, J N : H N × U ad → R + , associated with the Galerkin approximation (A.6) is given by J N (Π N y , u ) = (cid:90) T [ G ( y N ( s ; Π N y , u )) + E ( u ( s ))] d s, y ∈ H , (A.13)and the corresponding optimal control problem reads: min J N (Π N y , u ) subject to ( y N , u ) ∈ L (0 , T ; H N ) × U ad solves (A.6)with y N (0) = Π N y ∈ H N . ( P N )Then, as a consequence of [CKL17, Corollary 2.14] we can deduce the following theorem about the errorestimates between the full optimal control and the optimal control from Galerkin approximations. Theorem A.1.
Assume (H0) – (H4) and (A.12) hold. Assume also that for each y in H , both ( P ) and ( P N ) admit an optimal control, denoted by u ∗ and u ∗ N , respectively. Assume furthermore that there exists σ > suchthat the following local growth condition is satisfied for the cost functional J defined in (A.11) : σ (cid:107) u ∗ − v (cid:107) qL q (0 ,T ; V ) ≤ J ( y , v ) − J ( y , u ∗ ) , (A.14) for all v in some neighborhood W ⊂ U ad of u ∗ , with U ad given by (A.2) . Assume finally that u ∗ N lies in W .Then there exists γ > such that (cid:107) u ∗ − u ∗ N (cid:107) qL q (0 ,T ; V ) ≤ σ Lip( G| B ) (cid:104) √ T + γT (cid:105) (cid:16) (cid:107) Π ⊥ N y ( · ; y , u ∗ ) (cid:107) L (0 ,T ; H ) + 2 (cid:107) Π ⊥ N y ( · ; y , u ∗ N ) (cid:107) L (0 ,T ; H ) (cid:17) , (A.15) where B denotes the ball in H centered at the origin with radius C , with C being the same as given in Assumption (H4) , and Π ⊥ N = Id H − Π N . As pointed out in [CKL17, Remark 2.13], it is not clear a priori that lim N →∞ (cid:107) Π ⊥ N y ( · ; y , u ∗ N ) (cid:107) L (0 ,T ; H ) = 0 . (A.16)The reason relies on the dependence on u ∗ N of (cid:107) Π ⊥ N y ( · ; y , u ∗ N ) (cid:107) L (0 ,T ; H ) , where u ∗ N denotes the controlsynthesized from the N -dimensional Galerkin approximation. However, for the special case of Galerkinapproximations constructed from eigenbasis, the convergence in (A.16) is guaranteed to hold under the assump-tions (H0) – (H4) of Theorem A.1 above; see [CKL17, Lemma 2.16]. As a result, when N tends to infinity, u ∗ N converges to the optimal control u ∗ .We have furthermore, the following uniform convergence result about the controlled solutions as a resultof [CKL17, Theorem 2.6]. Theorem A.2.
Assume that (H0) – (H4) hold and that the operator C : V → H is locally Lipschitz in the sensegiven in (A.8) (with C (0) = 0 ). Assume furthermore that the set U in Assumption (H2) is a compact subset of V , with q > therein.Then, for any ( y , u ) in H × U ad , the solution y N of the Galerkin approximation (A.6) converges uniformlyto the mild solution y of (A.1) in the sense that: lim N →∞ sup u ∈U ad sup t ∈ [0 ,T ] (cid:107) y N ( t ; Π N y , u ) − y ( t ; y , u ) (cid:107) H = 0 . (A.17)Note that assumptions (A0), (A3) and (A6) in [CKL17, Theorem 2.6] correspond respectively to (H0) , (H2) and (H4) assumed here, and assumption (A5) in [CKL17] corresponds to (H2) together with U being compact.Assumption (A7) required in [CKL17, Theorem 2.6] is actually here a consequence of (H0) - (H4) with (H2) assumed for q > .Finally, we emphasize that assumptions (A1) and (A2) required by [CKL17, Theorem 2.6] follow from (H1) .Indeed, as pointed out above, thanks to (H1) , the eigenfunctions of L form an orthonormal basis of H . Then,the operator L N defined by (A.7) as the eigen projection of L onto H N given by (A.4) clearly satisfies lim N →∞ (cid:107) L N φ − Lφ (cid:107) H = 0 , ∀ φ ∈ D ( L ) . (A.18)Assumption (A2) in [CKL17] is thus satisfied. Since L is self-adjoint, it is also clear that the linear flow e L N t : H N → H N generated by L N satisfies (cid:107) e L N t (cid:107) ≤ e β t , N ∈ N , t ≥ , (A.19)where β is the largest eigenvalue of L . Defining the extension T N ( t ) : H → H of e L N t to be T N ( t ) φ = e L N t Π N φ + (Id H − Π N ) φ, φ ∈ H , (A.20)Assumption (A1) in [CKL17] follows then by choosing the parameters M and ω therein to be M = 1 and ω = max { β , } .Appendix B. Convergence and error estimates results for the optimal control of the KPPHequationWe check here, in the context of the optimal control of the KPPH equation, the assumptions of Theorem A.2and Theorem A.1 about the convergence and error estimates results, respectively, allowing us in particular todeduce Theorem 3.1 of Sec. 3.1.Within this context, (2.19) and (3.12) play the role of ( P ) and ( P N ), respectively. The optimal problems(2.19) and (3.12) are posed on the time interval [ τ, T ] but a simple change of variable t (cid:48) = t − τ allows usto frame these as in Appendix A, that is over the time interval [0 , T ] . We operate this shift below to ease thepresentation.The verification of the required assumptions is organized in several steps. Step 1: Verifications of (H0), (H1), and (H3).
Let H = L (Ω) . We first put the IVP (2.13b)-(2.13d) intothe form (A.1). The corresponding operators L : D ( L ) → H , F : H → H and C : H → H are then naturallydefined as follows: Ly = D ∇ y + µ ( · ) y, y ∈ D ( L ) = H (Ω) ∩ (cid:110) y ∈ H (Ω) (cid:12)(cid:12) ∂y∂ n = 0 (cid:111) ,F ( y ) = − ν ( · ) y − δρ (cid:15) ( y ) , y ∈ H , C ( u ) = u, u ∈ H . (B.1)It is standard that Assumptions (H0) and (H1) are satisfied for the elliptic operator L defined in (B.1); seee.g. [Paz83]. It can be checked that F is locally Lipschitz as mapping from H to H , in the sense of (A.8). ThusAssumption (H3) is satisfied. It is also clear that C defined above is Lipschitz on H and satisfies furthermorethat C (0) = 0 . PTIMAL MANAGEMENT OF HARVESTED POPULATION AT THE EDGE OF EXTINCTION 27
Step 2: Verification of (H2).
Recall that the admissible set U ad is defined by (2.18). Since ϕ j = √ | Λ j | χ Λ j (cf. (2.17)), each admissible control u in U ad takes value in the following subset U of H at a given time instant: U = (cid:110) ψ ∈ H (cid:12)(cid:12) ψ ( x ) = K (cid:88) j =1 b j (cid:112) | Λ j | χ Λ j ( x ) , x ∈ Ω , ≤ b j ≤ C j (cid:111) . (B.2)Thus U ad in (2.18) can be rewritten as U ad = L (0 , T ; U ) . (B.3)Taking the space V in (H2) to be L (Ω) , it is clear that U is a bounded set in V . Assumption (H2) is thusverified with q = 2 because of (B.3). We show below that U is furthermore a compact set in V , which isrequired in Theorem A.2. Every sequence ( ψ n ) in U has a convergent subsequence in U . Indeed for each ≤ j ≤ K , the corresponding b nj has a convergent subsequence as taking value in the bounded and closedsubset [0 , C j ] of R . Because we have a finite number of such convergent subsequences, one can always find acommon extraction and thus ( ψ n ) for which convergence holds towards an element of U . Thus U is a compactsubset of V = L (Ω) . Step 3: Verification of (H4).
We first prove the uniform bound (A.9a) for the solution y to (2.13b)-(2.13d).As it will appear below, this uniform bound can be derived as soon as one can ensure that y ≥ when y ≥ .This property is a consequence of the maximum principle which holds for (2.13b)-(2.13d), a known fact but ofwhich we provide a proof of using the Stampacchia truncation method [MS68]. The arguments are standardbut are sketched below as they provide also useful insights to prove the required a priori bounds for the solutionto the Galerkin approximation of (2.13b)-(2.13d). For this purpose, one defines the negative part v − of ameasurable function v : Ω → R by v − ( x ) = min { v ( x ) , } . (B.4)A lemma due to Stampacchia ensures that if v is in H (Ω) , then v − lies also in H (Ω) and ∇ v − = χ { v< } ∇ v. (B.5)The so-called truncation method of Stampacchia allows us to obtain a powerful identity for a broad class ofinhomgeneous heat problem ∂ t y = D ∇ y + f, in (0 , T ) × Ω , (B.6a) ∂y∂ n = 0 , on (0 , T ) × ∂ Ω , (B.6b) y (0 , x ) = y ( x ) , x ∈ Ω . (B.6c)This identity is derived from the variational formulation of this problem, by using (B.5) with v = y − . It gives (cid:12)(cid:12) y − ( t ) (cid:12)(cid:12) L (Ω) + D (cid:90) t (cid:12)(cid:12) ∇ y − ( s ) (cid:12)(cid:12) L (Ω) d s = (cid:90) t (cid:90) Ω f y − d x d s + 12 (cid:12)(cid:12) y − (0) (cid:12)(cid:12) L (Ω) . (B.7)From this identity the classical maximum principle can be deduce in the sense that if f ≥ and y ≥ a.e. then y ≥ a.e. Let us introduce N ( x, y, u ) = µ ( x ) y − ν ( x ) y − δρ (cid:15) ( y ) + u ( t, x ) , (B.8)This result cannot be applied directly with f = N because f does not have the good sign. However the identity(B.7) is useful to conclude about the positivity of y when f = (cid:101) N , with (cid:101) N ( x, y, u ) = µ ( x ) y − ν ( x ) y | y | − δρ (cid:15) ( y ) + u ( t, x ) . (B.9)Indeed, first note that for f = (cid:101) N , (cid:90) t (cid:90) Ω f y − d x d s ≤ (cid:107) µ (cid:107) L ∞ (Ω) (cid:90) t (cid:12)(cid:12) y − ( s ) (cid:12)(cid:12) L (Ω) d s − (cid:90) t (cid:90) Ω ν ( x ) y | y | y − d x d s (B.10) since ρ (cid:15) ( y ) y − = 0 , because ρ (cid:15) ( s ) = 0 for all s ≤ , and uy − ≤ , because u lies in U ad and thus is nonnegative.Then, since νy | y | y − = ν | y | ( y − ) ≥ , (B.11)we deduce from (B.10) that (cid:90) t (cid:90) Ω (cid:101) N y − d x d s ≤ (cid:107) µ (cid:107) L ∞ (Ω) (cid:90) t (cid:12)(cid:12) y − ( s ) (cid:12)(cid:12) L (Ω) d s. (B.12)Assuming now that y ≥ a.e., we obtain from (B.7) that (cid:12)(cid:12) y − ( t ) (cid:12)(cid:12) L (Ω) ≤ (cid:107) µ (cid:107) L ∞ (Ω) (cid:90) t (cid:12)(cid:12) y − ( s ) (cid:12)(cid:12) L (Ω) d s. (B.13)Because y − (0) = ( y ) − = 0 , the Gronwall’s lemma in its integral form allows us to conclude that y − = 0 a.e. in (0 , T ) × Ω . Thus we have proved that y ≥ a.e. if y solves (B.6) with y ≥ and f = (cid:101) N a.e.From this we deduce that | y | = y a.e. x , t . Thus equation (B.6a) with f = (cid:101) N is the same as equation (2.13b)when y ≥ and in the end we have also found that y is a positive solution to (2.13b)-(2.13d) when y ≥ .On the other hand, the inner product with y on both sides of (2.13b) leads to
12 dd t (cid:12)(cid:12) y (cid:12)(cid:12) L (Ω) = (cid:104) Ly + F ( y ) + C ( u ( t )) , y (cid:105) H = − D (cid:12)(cid:12) ∇ y (cid:12)(cid:12) L (Ω) + (cid:90) Ω (cid:0) µ ( x ) y − ν ( x ) y (cid:1) d x − δ (cid:90) Ω ρ (cid:15) ( y ) y d x + (cid:90) Ω uy d x (B.14)Thus because y is a positive solution and ρ (cid:15) ≥ , we obtain
12 dd t (cid:12)(cid:12) y (cid:12)(cid:12) L (Ω) ≤ (cid:107) µ (cid:107) ∞ (cid:12)(cid:12) y (cid:12)(cid:12) L (Ω) + C (cid:90) Ω y d x, (B.15)with C = max (cid:110) C √ | Λ | , . . . , C K √ | Λ K | (cid:111) .By remarking that by Hölder’s inequality, we have (cid:90) Ω y d x ≤ (cid:112) | Ω | (cid:12)(cid:12) y (cid:12)(cid:12) L (Ω) ≤
12 ( | Ω | + (cid:12)(cid:12) y (cid:12)(cid:12) L (Ω) ) . Using this in (B.15), we obtain then
12 dd t (cid:12)(cid:12) y (cid:12)(cid:12) L (Ω) ≤ c (cid:12)(cid:12) y (cid:12)(cid:12) L (Ω) + c , (B.16)where c = (cid:107) µ (cid:107) ∞ + C , c = C | Ω | . (B.17)The uniform bound for y in (A.9a) follows then from Gronwall’s inequality.We turn now to the proof of the uniform bound (A.9b). First, we consider Galerkin approximation ˜ y N ofa modified version of (2.13b)-(2.13c) in which the reaction term in (2.13b) is replaced by (cid:101) N . We assume theinitial datum ˜ y to be nonnegative. The corresponding N -dimensional Galerkin approximation ˜ y N satisfies thus ∂ t ˜ y N = D ∇ ˜ y N + Π N (cid:101) N ( x, ˜ y N , u ) . (B.18)From this equation one gets,
12 dd t (cid:12)(cid:12) ˜ y N (cid:12)(cid:12) L (Ω) = − D (cid:12)(cid:12) ∇ ˜ y N (cid:12)(cid:12) L (Ω) + (cid:104) Π N (cid:101) N ( x, ˜ y N , u ) , ˜ y N (cid:105) . (B.19) PTIMAL MANAGEMENT OF HARVESTED POPULATION AT THE EDGE OF EXTINCTION 29
Since ˜ y N lies in H N , we have (cid:104) Π N (cid:101) N ( x, ˜ y N , u ) , y N (cid:105) = (cid:104) (cid:101) N ( x, ˜ y N , u ) , ˜ y N (cid:105) = (cid:90) Ω (cid:0) µ ( x )˜ y N − ν ( x ) | ˜ y N | ˜ y N (cid:1) d x − δ (cid:90) Ω ρ (cid:15) ( y N ) y N d x + (cid:90) Ω u ˜ y N d x ≤ (cid:107) µ (cid:107) ∞ (cid:90) Ω ˜ y N d x + ( δ + C ) (cid:90) Ω | ˜ y N | d x, (B.20)where C is the same as in (B.15) and we have used the fact that | ρ (cid:15) ( y N ) | is bounded above by ; cf. (2.2).We infer thus that
12 dd t (cid:12)(cid:12) ˜ y N (cid:12)(cid:12) L (Ω) ≤ (cid:107) µ (cid:107) ∞ (cid:12)(cid:12) ˜ y N (cid:12)(cid:12) L (Ω) + ( δ + C ) (cid:90) Ω | ˜ y N | d x, (B.21)leading as for y satisfying (B.15) to a uniform bound for ˜ y N , and consequently (A.9b) holds for ˜ y N . Due toTheorem A.2, ˜ y N satisfies the convergence property (A.17) towards the solution ˜ y to the modified IVP (2.13b)-(2.13d) in which (cid:101) N replaces the reaction term. Now we know that ˜ y ≥ since ˜ y ≥ , and since | ˜ y | = ˜ y in thiscase, we get that ˜ y solves actually the original IVP (2.13b)-(2.13d) which we denote now by y and ˜ y by y .Smoothing arguments ensure that the mild solution y in L (Ω) to (2.13b)-(2.13d) is actually continuous on Ω . Thus because we have, for U ad defined in (B.3) sup u ∈U ad sup t ∈ [0 ,T ] | ˜ y N ( t ; Π N y , u ) − y ( t ; y , u ) | L (Ω) −→ N →∞ , (B.22)we have ˜ y N > for N sufficiently large (when y > ) and thus (B.18) reduces to the Galerkin approximationof the original IVP (2.13b)-(2.13d).Finally, note that condition (A.12) holds for the cost functional J defined in (2.12). The checking of (A.14)is more involved but can be derived by adapting the proof of [CHW17, Theorem 5.3] to our context. The errorestimates (A.15) of Theorem A.1 also hold with q = 2 and V = L (Ω) and Theorem 3.1 of Sec. 3.1 is thusproved. References [Ama76] H. Amann, Fixed point equations and nonlinear eigenvalue problems in ordered banach spaces , SIAM review (1976),no. 4, 620–709.[BC03] B. Bonnard and M. Chyba, Singular Trajectories and Their Role in Control Theory , Mathématiques & Applications(Berlin), vol. 40, Springer, 2003.[BH75] A. E. Bryson, Jr. and Y. C. Ho,
Applied Optimal Control , Hemisphere Publishing Corp. Washington, D. C., 1975.[BH02] H. Berestycki and F. Hamel,
Front propagation in periodic excitable media , Communications on Pure and AppliedMathematics (2002), no. 8, 949–1032.[BHR05a] H. Berestycki, F. Hamel, and L. Roques, Analysis of the periodically fragmented environment model: I–Species persistence ,Journal of Mathematical Biology (2005), no. 1, 75–113.[BHR05b] , Analysis of the periodically fragmented environment model: IIâĂŤBiological invasions and pulsating travellingfronts , Journal de Mathématiques purés et appliquées (2005), no. 8, 1101–1146.[BHTS04] J.E.M. Baillie, C. Hilton-Taylor, and S.N. Stuart, A global species assessment , Tech. report, International Union forConservation of Nature (IUCN), 2004.[BM77] J.R. Beddington and R.M. May,
Harvesting natural populations in a randomly fluctuating environment , Science (1977), no. 4302, 463–465.[Bré10] H. Brézis,
Functional Analysis, Sobolev Spaces and Partial Differential Equations , Springer, 2010.[BXY14] W. A. Brock, A. Xepapadeas, and A. N. Yannacopoulos,
Optimal control in space and time and the management ofenvironmental resources , Annu. Rev. Resour. Econ. (2014), 33–68.[CC89] R. S. Cantrell and C. Cosner, Diffusive logistic equations with indefinite weights: population models in disrupted environ-ments , Proc. R. Soc. Edinb. Sect. A (1989), no. 3-4, 293–318.[CC04] ,
Spatial Ecology via Reaction-Diffusion Equations , John Wiley & Sons, 2004.[CGLW16] M. D. Chekroun, M. Ghil, H. Liu, and S. Wang,
Low-dimensional Galerkin approximations of nonlinear delay differentialequations , Disc. Cont. Dyn. Sys. A (2016), no. 8, 4133–4177. [CH98] T. Cazenave and A. Haraux, An Introduction to Semilinear Evolution Equations , Oxford Lecture Series in Mathematicsand its Applications, vol. 13, The Clarendon Press, Oxford, 1998.[CHW17] E. Casas, R. Herzog, and G. Wachsmuth,
Analysis of spatio-temporally sparse optimal control problems of semilinearparabolic equations , ESAIM: Control, Optimisation and Calculus of Variations (2017), no. 1, 263–295.[CKL17] M. D. Chekroun, A. Kröner, and H. Liu, Galerkin approximations of nonlinear optimal control problems in Hilbert spaces ,Electronic Journal of Differential Equations (2017), no. 189, 1–40.[CKL18] M. D. Chekroun, A. Kröner, and H. Liu,
Galerkin approximations for the optimal control of nonlinear delay differentialequations , Hamilton-Jacobi-Bellman Equations. Numerical Methods and Applications in Optimal Control. D. Kalise, K.Kunisch, and Z. Rao (Eds.), vol. 21, Berlin, Boston: De Gruyter, 2018, pp. 61–96.[CR06] M.D. Chekroun and L.J. Roques,
Models of population dynamics under the influence of external perturbations: mathe-matical results , Comptes Rendus Mathématique (2006), no. 5, 307–310.[Ded10] L. Dedè,
Reduced basis method and a posteriori error estimation for parametrized linear-quadratic optimal controlproblems , SIAM Journal on Scientific Computing (2010), 997–1019.[Fah03] L. Fahrig, Effects of habitat fragmentation on biodiversity , Annual review of ecology, evolution, and systematics (2003),no. 1, 487–515.[FHLW12] T. Franke, R. H. W. Hoppe, C. Linsenmann, and A. Wixforth, Projection based model reduction for optimal design of thetime-dependent Stokes system , Constrained Optimization and Optimal Control for Partial Differential Equations, Springer,2012, pp. 75–98.[Fis37] R. A. Fisher,
The wave of advance of advantageous genes , Annals of eugenics (1937), no. 4, 355–369.[FLVP12] H. Finotti, S. Lenhart, and T. Van Phan, Optimal control of advective direction in reaction-diffusion population models ,Evolution Equations & Control Theory (2012), 81–107.[GH89] W. M. Getz and R.G. Haight, Population harvesting: Demographic models of fish, forest, and animal resources , vol. 27,Princeton University Press, 1989.[GK11] M. A. Grepl and M. Kärcher,
Reduced basis a posteriori error bounds for parametrized linear-quadratic elliptic optimalcontrol problems , C. R. Acad. Sci. Paris, Ser. I (2011), no. 15, 873–877.[HK98] M. Hinze and K. Kunisch,
On suboptimal control strategies for the Navier-Stokes equations , ESAIM: Proceedings, vol. 4,1998, pp. 181–198.[HK00] ,
Three control methods for time-dependent fluid flow , Flow, Turbulence and Combustion (2000), 273–298.[HSV95] R. F. Hartl, S. P. Sethi, and R. G. Vickson, A survey of the maximum principles for optimal control problems with stateconstraints , SIAM review (1995), no. 2, 181–218.[HV05] M. Hinze and S. Volkwein, Proper orthogonal decomposition surrogate models for nonlinear dynamical systems: errorestimates and suboptimal control , Dimension Reduction of Large-Scale Systems, Lect. Notes Comput. Sci. Eng., vol. 45,Springer, Berlin, 2005, pp. 261–306.[IK08] K. Ito and K. Kunisch,
Reduced-order optimal control based on approximate inertial manifolds for nonlinear dynamicalsystems , SIAM J. Numer. Anal. (2008), no. 6, 2867–2891.[Kir12] D. E. Kirk, Optimal Control Theory: An Introduction , Dover Publications, 2012.[KPP37] A.N. Kolmogorov, I.G. Petrovsky, and N.S. Piskunov,
Etude de l’équation de la diffusion avec croissance de la quantitéde matière et son application à un problème biologique , Bulletin Université d’État à Moscou (Bjul. Moskowskogo Gos.Univ.), Série internationale A (1937), 1–26.[KS08] K. Kurata and J. Shi, Optimal spatial harvesting strategy and symmetry-breaking , Appl. Math. Optim. (2008), no. 1,89–110.[MS68] M.R.V Murthy and G. Stampacchia, Boundary value problems for some degenerate-elliptic operators , Annali di MatematicaPura ed Applicata (1968), no. 1, 1–122.[Mur01] J.D. Murray, Mathematical Biology II: Spatial Models and Biomedical Applications , Springer New York, 2001.[Mur07] ,
Mathematical Biology: I. An Introduction , vol. 17, Springer Science & Business Media, 2007.[Neu03] M.G. Neubert,
Marine reserves and optimal harvesting , Ecology Letters (2003), no. 9, 843–849.[OL13] A. Okubo and S. A. Levin, Diffusion and ecological problems: Modern perspectives , vol. 14, Springer Science & BusinessMedia, 2013.[OS02] S. Oruganti and J. Shi,
Diffusive logistic equation with constant yield harvesting. i: Steady states , Trans. Am. Math. Soc. (2002), no. 9, 3601–3619.[Paz83] A. Pazy,
Semigroups of Linear Operators and Applications to Partial Differential Equations , Applied MathematicalSciences, vol. 44, Springer-Verlag, New York, 1983.[PBGM64] L. S. Pontryagin, V. G. Boltyanskii, R. V. Gamkrelidze, and E. F. Mishchenko,
The Mathematical Theory of OptimalProcesses , Translated by D. E. Brown, A Pergamon Press Book. The Macmillan Co., New York, 1964.[Rav00] S.S. Ravindran,
A reduced-order approach for optimal control of fluids using proper orthogonal decomposition , Interna-tional Journal for Numerical Methods in Fluids (2000), no. 5, 425–448. PTIMAL MANAGEMENT OF HARVESTED POPULATION AT THE EDGE OF EXTINCTION 31 [RB99] J. G Robinson and R. E. Bodmer,
Towards wildlife management in tropical forests , The Journal of wildlife management(1999), 1–13.[RC07] L. Roques and M. D. Chekroun,
On population resilience to external perturbations , SIAM J. Appl. Math. (2007),133–153.[RC10] L. Roques and M.D. Chekroun, Does reaction-diffusion support the duality of fragmentation effect? , Ecological Complexity (2010), no. 1, 100–106.[Sch54] M. B. Schaefer, Some aspects of the dynamics of populations important to the management of the commercial marinefisheries , Inter-American Tropical Tuna Commission Bulletin (1954), no. 2, 23–56.[Sch91] , Some aspects of the dynamics of populations important to the management of the commercial marine fisheries ,Bulletin of Mathematical Biology (1991), no. 1-2, 253–279.[SFRAS02] P.A. Stephens, F. Frey-Roos, W. Arnold, and W.J. Sutherland, Sustainable exploitation of social species: a test andcomparison of models , Journal of Applied Ecology (2002), no. 4, 629–642.[SK97] N. Shigesada and K. Kawasaki, Biological Invasions: Theory and Practice , Oxford University Press, UK, 1997.[SKT86] N. Shigesada, K. Kawasaki, and E. Teramoto,
Traveling periodic waves in heterogeneous environments , TheoreticalPopulation Biology (1986), no. 1, 143–160.[Tem97] R. Temam, Infinite-Dimensional Dynamical Systems in Mechanics and Physics , 2nd ed., Applied Mathematical Sciences,vol. 68, Springer-Verlag, New York, 1997.[Trö10] F. Tröltzsch,
Optimal Control of Partial Differential Equations: Theory, Methods and Applications , Graduate Studies inMathematics, vol. 112, American Mathematical Society, 2010.[TV09] F. Tröltzsch and S. Volkwein,
POD a posteriori error estimates for linear-quadratic optimal control problems , Comput.Optim. Appl. (2009), 83–115.(MDC) Department of Earth and Planetary Sciences, Weizmann Institute, Rehovot 76100, Israel; Department ofAtmospheric and Oceanic Sciences, University of California, Los Angeles, CA 90095-1565, USA E-mail address : [email protected] (HL) Department of Mathematics, Virginia Tech, Blacksburg, Virginia 24061, USA E-mail address ::