figurec

# Optimal management of harvested population at the edge of extinction

Mickaël D. Chekroun  Department of Earth and Planetary Sciences, Weizmann
Institute, Rehovot 76100, Israel; Department of Atmospheric and Oceanic
Sciences, University of California, Los Angeles, CA 90095-1565, USA
<EMAIL_ADDRESS>and  Honghu Liu  Department of Mathematics, Virginia
Tech, Blacksburg, Virginia 24061, USA<EMAIL_ADDRESS>
(Date: July 22, 2020)

###### Abstract.

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 harvesting term
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 result about 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 differential equation
(ODE) with a harvesting term, obtained as a Galerkin approximation of the
original partial differential equation (PDE) model. A critical allowable
fraction $\underline{\alpha}$ of the reserve’s population is inferred from the
reduced logistic ODE with a harvesting term. This estimate obtained from the
reduced model allows us to distinguish sharply between survival and extinction
for the full PDE itself, and thus to declare whether a control strategy leads
to success or failure for the corresponding rescue operation while ensuring
survival in the reserve’s population. In dynamical terms, 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 its generality,
establishes bridges interesting to explore between optimal control problems of
ODEs with a harvesting term and their PDE counterpart.

###### Contents

  1. 1 Introduction
  2. 2 Optimal control of harvested population
    1. 2.1 The KPP model with harvesting term in heterogeneous environment
    2. 2.2 Main results from [RC07]
    3. 2.3 Optimal control of harvested population at the edge of extinction
    4. 2.4 Control constraints from reserve management
  3. 3 Optimal control from Galerkin approximations
    1. 3.1 Convergence results and error estimates about the optimal control
    2. 3.2 PMP solution to the reduced optimal control problem (3.12)
  4. 4 Nearly optimal controls from low-dimensional surrogates: Numerical results
    1. 4.1 Numerical setup
    2. 4.2 Choice of the control constraints
    3. 4.3 Control from effective reduced Galerkin systems
    4. 4.4 First numerical results
    5. 4.5 Effects of the control constraints on the rescue operation
  5. A Galerkin approximations of nonlinear optimal control problems in Hilbert spaces: General convergence results
  6. B Convergence and error estimates results for the optimal control of the KPPH equation

##  1\. Introduction

Overexploitation has led to the extinction of many species [BHTS04].
Traditionally, models of ordinary differential equations (ODEs) or difference
equations have been used to estimate the maximum sustainable yields from
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, the ODE
models are generally of the type

$\frac{\,\mathrm{d}U}{\,\mathrm{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 harvest function. In these models, the
most commonly used growth function is logistic, with $F(U)=U(\mu-\nu U)$
[BM77, Sch54], where $\mu>0$ is the intrinsic growth rate of the population
and $\nu>0$ models its susceptibility to crowding effects.

Different harvesting strategies $Y(U)$ have been considered in the literature
and are used in practical resource management. A very common one is the
_constant-yield harvesting_ strategy, where a constant number of individuals
are removed per unit of time: $Y(U)=\delta$, with $\delta$ a positive
constant. This harvesting function naturally appears when a quota is set on
the harvesters [RB99, SFRAS02]. Another frequently used harvesting strategy is
the _proportional harvesting_ strategy (also called _constant-effort
harvesting_), where a constant proportion of the population is removed. It
leads to a harvesting function of the type $Y(U)=\delta U$.

Reaction-diffusion equations have also been used extensively in modeling the
spatiotemporal behavior of a species of organism [SK97, CC04, Mur01, Mur07,
OL13], and harvesting effects have been analyzed within this class of models;
see e.g. [OS02, RC07, RC10]. Within this global picture, we focus on a
particular class of models, 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.1b) below) in which the local growth rate of the population and the limit
effects to crowding are spatial-dependent coefficients. The harvest function
is a quasi-constant-yield harvesting term depending on a threshold below which
harvesting is progressively abandoned. This model is analyzed in [CR06, RC07]
and the main results about its steady state analysis and its dynamical
behavior, are recalled in Sec. 2.2 below for the reader’s convenience. In
particular, the asymptotic behavior of the KPPH model was characterized as a
function of the harvesting intensity $\delta$; see (2.1b). It was proved in
[RC07] that if $\delta$ is smaller than a critical value $\delta^{*}$, 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
$\delta$ is larger than another critical value $\delta_{2}$, which is slightly
above $\delta^{*}$ in practice, the population density eventually settles down
to a “remnant” state, everywhere below the same small threshold. The
population can be considered as extinct in this case.

Within this context, we assume that given an unprotected area $\Omega$ in
which the population (of let’s say some mammals) evolves according to a KPPH
equation, is under threat of extinction due to some illegal hunting that
superimposes to some allowed harvesting, causing thereof a rise in the
harvesting intensity such that $\delta$ ends to jump above $\delta_{2}$. Given
a natural reserve $\widehat{\Omega}$ in which the same species evolves, we
address in this article the problem of saving the population under threat by
releasing in a controlled fashion new individuals from the reserve, into the
unprotected area, while avoiding to exert too much pressure on the reserve
that would be detrimental on it. Our control strategy is seriously constrained
in time as the operation is assumed to start when a (large) fraction of the
original population has been decimated and its course to extinction is, thus,
well engaged.

The resulting optimal control problem differs fundamentally from other optimal
control problems concerned with harvesting models that appeared in previous
studies which dealt mainly with the search for optimal harvesting strategies;
see e.g. [Neu03, KS08, FLVP12, BXY14]. Instead, we aim at designing optimal
rescue strategies while taking into account plausible factors that may arise
in real-life problems, from a reserve management 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 the underlying 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 class of nonlinear optimal control problems
in infinite dimension. Of course, the idea of designing approximations of
optimal solutions from finite-dimensional approximations is not new, and a
great deal of works have addressed this 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 of nonlinear 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 the optimal control of the KPPH model considered in
this article and the optimal control built from any of its Galerkin
approximations. These error estimates are supplemented by convergence of the
corresponding controlled 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 proved in
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 of a natural underlying spectral problem (see (2.3b)
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 of Sec. 4 are presented within the context of
one-dimensional environments and for various level of fragmentation of the
population’s habitat. For the parameter regime considered, we show that from a
one-dimensional, ordinary differential equation (ODE) approximation, nearly
optimal solutions can be derived. The corresponding Galerkin approximation of
the KPPH model, reducing to a simple logistic ODE with a harvesting
term111Whose solution contains most of the PDE solution’s energy., allows us
in turn to provide great insights about the optimal control problem of the PDE
model. The main result, from an ecological perspective, is indeed expressed in
terms of a fraction $\alpha$ 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 transported from the reserve into the unprotected area.
We show that a critical fraction $\underline{\alpha}$ can be inferred from the
reduced logistic ODE with a harvesting term. This estimate obtained from the
reduced equation is indeed shown to be very useful as it allows us to
distinguish sharply between survival and extinction for the full PDE itself,
and thus to declare whether a control strategy leads to success or failure for
the corresponding rescue operation while ensuring survival in the reserve’s
population. In terms of optimal control, this result translates into the
existence of controls very close to each other, one of which leading to a
significant survival of the population while the other leading to its
extinction. In dynamical terms, 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.

##  2\. Optimal control of harvested population

###  2.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 follows. Let
$\Omega$ be a smooth bounded and connected domain of $\mathbb{R}^{d}$ ($d\geq
1$). We consider

$\displaystyle\;\partial_{t}y=D\nabla^{2}y+\mu(x)y-\nu(x)y^{2}-\delta\rho_{\epsilon}(y),\quad(t,x)\in(0,\infty)\times\Omega,$
(2.1a) $\displaystyle\;\frac{\partial
y}{\partial{\boldsymbol{n}}}=0,\quad(t,x)\in[0,\infty)\times\partial\Omega.$
(2.1b)

Here $\boldsymbol{n}$ denotes the outward unit normal to the boundary
$\partial\Omega$.

This equation differs from the classical Fisher equation [Fis37] (also known
as the Kolmogorov-Petrovsky-Piskunov (KPP) equation [KPP37]), by its
spatially-dependent coefficients, $\mu(x)$ and $\nu(x)$, as well as its
harvesting term $\delta\rho_{\epsilon}(y)$. When $\delta=0$ in (2.1b), the
model reduces to the Shigesada-Kawasaki-Teramoto model described in [SKT86];
see also [BH02, BHR05a, BHR05b, SK97]. Such a problem fits with general
species assessment 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$. The coefficient $\mu$ represents the intrinsic growth rate
of the population, which is assume to be a measurable function of $x$ in
$L^{\infty}(\Omega)$. The spatial dependence of $\mu$ is introduced to account
for the possible impact of environmental heterogeneity [RC07]. The function
$\nu(x)>0$ (also in $L^{\infty}(\Omega)$) represents the susceptibility to
crowding effects and is interpreted as an intraspecific competition term.
Regions with higher values of $\mu(x)$ and lower values of $\nu(x)$ are
qualified as being more favorable, while, on the other hand, regions with
lower $\mu(x)$ and higher $\nu(x)$ values are considered as being less
favorable or, equivalently, more hostile.

The harvesting function $\rho_{\varepsilon}$ satisfies

$\rho_{\epsilon}\in C^{1}(\mathbb{R}),\quad\rho_{\epsilon}^{\prime}\geq
0,\quad\rho_{\epsilon}(s)=0,\;\;\forall s\leq
0,\quad\text{and}\quad\rho_{\epsilon}(s)=1,\;\;\forall s\geq\epsilon,$ (2.2)

where $\epsilon$ is a nonnegative parameter, taken to be sufficiently small in
a sense made precise in Sec. 2.2.

The term $\delta\rho_{\epsilon}(y)$ with $\delta\geq 0$, corresponds to a
quasi-constant-yield harvesting term. Indeed, for such a harvesting function,
the yield is constant in time whenever $y\geq\epsilon$, while it depends on
the population density when $y<\epsilon$. Note that the function
$\rho_{\epsilon}$ ensures the nonnegativity of the solutions to (2.1b); see
[RC07]. From a biological viewpoint, $\varepsilon$ corresponds to a threshold
below which harvesting is progressively abandoned. Considering constant-yield
harvesting functions without this threshold value would be unrealistic since
it would eventually lead to harvest on zero-populations.

###  2.2. Main results from [RC07]

The problem (2.1b) has been analyzed in [RC07]. Using sub- and supersolution
methods 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 the harvesting intensity $\delta$. In [RC07] it
was proved that if $\delta$ is smaller than a critical value $\delta^{*}$, the
population density converges to a “significant” state, which is everywhere
above a small threshold (depending on $\epsilon$) while if $\delta$ is larger
than $\delta_{2}$ (another threshold that is bigger than but close to
$\delta^{*}$), the population density $y(t,x)$ converges to a “remnant” state,
everywhere below this small threshold. Theorem 2.1 below summarizes the main
asymptotic and existence results from [RC07].

To formulate these results we recall some tools used in [RC07] and in the
formulation of our optimal control problem below. For this purpose, we first
consider the eigenvalue problem associated with the linearization of (2.1b) at
the trivial solution:

$\displaystyle\;-D\nabla^{2}\phi-\mu(x)\phi=\lambda\phi,\quad x\in\Omega,$
(2.3a) $\displaystyle\;\frac{\partial\phi}{\partial{\boldsymbol{n}}}=0,\quad
x\in\partial\Omega.$ (2.3b)

Let $\lambda_{1}$ denote the first eigenvalue and $\phi_{1}$ its corresponding
positive eigenfunction, unique when normalized. In particular $\phi_{1}$
satisfies

$\phi_{1}(x)>0\quad\forall\;x\in\Omega,\quad\text{and}\quad\|\phi_{1}\|_{\infty}=1.$
(2.4)

The above normalization is possible since the first eigenfunction has a fixed
sign in $\Omega$ as consequence of the Krein-Rutman theorem; see e.g. [Ama76].
We introduce next

$\underline{\phi_{1}}=\min_{x\in\Omega}\phi_{1},$ (2.5)

Note that $\underline{\phi_{1}}>0$ since $\phi_{1}$ does not vanish on the
boundary $\partial\Omega$; see [Ama76]. Note also that
$\underline{\phi_{1}}\leq 1$ thanks to (2.4).

Recall from [RC07, Definition 2.5] that, a stationary solution of (2.1b),
$p_{\delta}$, is called a significant solution if

$\min_{x\in\Omega}p_{\delta}\geq\frac{\epsilon}{\underline{\phi_{1}}}.$ (2.6)

On the other hand, a stationary solution $p_{\delta}$ is called remnant if

$\max_{x\in\Omega}p_{\delta}<\frac{\epsilon}{\underline{\phi_{1}}}.$ (2.7)

With these tools in hand, it was derived in [RC07] the following analytic
formulas allowing for estimating the critical harvesting intensity
$\delta^{\ast}$ leading to decline of the population towards a remnant steady
state (i.e. close to extinction):

$\quad\delta_{1}=\frac{\lambda_{1}^{2}\underline{\phi_{1}}}{\overline{\nu}(1+\underline{\phi_{1}})^{2}},\quad\delta_{2}=\frac{\lambda_{1}^{2}}{4\underline{\nu}},$
(2.8)

see [RC07, (2.15)] with $\alpha=1$ since $h(x)$ therein is identically equal
to 1 here. Finally, $\overline{\nu}=\max_{x\in\Omega}\nu(x)$ and
$\underline{\nu}=\min_{x\in\Omega}\nu(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 $\lambda_{1}<0$, then there exists a threshold
$\delta^{*}\geq 0$ such that

  * $(i)$

if $\delta\leq\delta^{*}$, there exists at least one positive stationary
significant solution of (2.1b), whereas

  * $(ii)$

if $\delta>\delta^{*}$, there is no positive stationary significant solution
of (2.1b).

Moreover,

  * $(iii)$

if $\lambda_{1}<0$ and $\delta\leq\delta_{1}$, there exists a positive
stationary significant solution $p_{\delta}$ of (2.1b) such that

$p_{\delta}\geq-\frac{\lambda_{1}\phi_{1}}{\overline{\nu}(1+\underline{\phi_{1}})}.$
(2.9)

  * $(iv)$

if $\lambda_{1}<0$ and $\delta>\delta_{2}$, the only possible positive bounded
stationary solutions of (2.1b) are remnant.

Asymptotic behavior [RC07]. Assume that the initial datum $y(0,x)$ for Eq.
(2.1b) is taken to be $p_{0}$, the unique nonnegative steady state of Eq.
(2.1b) when $\delta=0$ [BHR05a]. Then, the solution $y(t,x)$ is non-increasing
in $t$ and its asymptotic behavior is characterized as follows:

  * $(v)$

Population resilience: If $\delta\leq\delta^{*}$, the solution $y(t,\cdot)$
converges to $p_{\delta}$ uniformly in $\Omega$ as $t$ goes to infinite, where
$p_{\delta}$ in (2.9) is the maximal stationary significant solution of
(2.1b).

  * $(vi)$

Population extinction: If $\delta>\delta_{2}$, the solution $y(t,\cdot)$
converges uniformly in $\Omega$ to a remnant solution of (2.1b).

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 $\delta_{1}$
and $\delta_{2}$ defined in (2.8) provide easily computable bounds for the
maximum allowable harvesting intensity $\delta^{*}$. These bounds have 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
$\delta^{\ast}$, strongly supported that environmental fragmentation of the
habitat has a significant impact on the maximum sustainable yield associated
with $\delta^{\ast}.$ Essentially, the more fragmented is the habitat, the
more the population is susceptible to decline towards a remnant state under
harvesting pressure.

###  2.3. Optimal control of harvested population at the edge of extinction

We describe hereafter the optimal control problem we aim at solving for
(2.1b). The motivation is as follows. Assume that the population whose
evolution is governed by (2.1b) 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. $\delta>\delta_{2}$. Indeed, without additional external
intervention, 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 individuals
from the same species coming e.g. from a natural reserve, into the unprotected
area $\Omega$, once the population in that 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 under constraints. We assume that the reserve contains the
same species that we are aiming at saving in the unprotected area and that the
reserve’s population is at an equilibrium, i.e. that it has reached a steady
state. Because no harvesting is exerted in the reserve, we assume that the
population dynamics in the reserve is governed by (2.1b) with $\delta=0$,
namely the KPP equation in an heterogeneous environment [SKT86]. In the
general case, the domain $\widehat{\Omega}$ of the reserve is different from
that of the unprotected area, although it can share similar features such as a
similar area and proportion of the population’s habitat (in which $\mu>0$).
More precise considerations about the reserve will be formulated in the
numerical section. For the moment, we specify our goal in general terms useful
to frame our optimal control problem. Thus, we are aiming at transporting a
fraction of the population from the reserve $\widehat{\Omega}$ to the
unprotected area $\Omega$. The question is to determine which fraction to use
in order to save the population at extinction threat in $\Omega$ while not
only avoiding to remove all the population in the reserve but also not causing
extinction within this same reserve.

Warning time $t=\tau$. We assume also that the population in $\Omega$ before
harvesting was exerted ($\delta=0$), had reached the stationary state $p_{0}$
and that the total population size can be monitored in $\Omega$. Based on this
monitoring, we assume that the authorities declare the population at threat
when the following warning threshold $P_{c}$ is reached

$P_{c}=\beta\int_{\Omega}y(0,x)\,\mathrm{d}x,\quad 0<\beta<1,$ (2.10)

where $y(0,x)=p_{0}$ is the initial density of the population at time $t=0$.

Recall that we assume $\delta>\delta_{2}$ in (2.1b) for the unprotected area
$\Omega$. As mentioned above, a situation for which $\delta$ would be greater
than $\delta_{2}$ may arise due for instance to illegal hunting practiced
without respect of some quota and 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.1b) and because $y(t,x)$ is
continuous in time and space, we conclude that there exists a time instant,
$t=\tau$, such that $\int_{\Omega}y(\tau,x)\,\mathrm{d}x=P_{c}$. Note that we
want to avoid to reach extinction so $P_{c}$ and thus $\beta$ in (2.10) must
be chosen such that

$P_{c}>\frac{\epsilon|\Omega|}{\underline{\phi_{1}}},$ (2.11)

where $|\Omega|$ denotes the $d$-dimensional volume of $\Omega$. In other
words one wishes to trigger an alert signal when the situation corresponds
still to a significant population size according to (2.6). The time instant
$\tau$ serves us thus to trigger an alert signal from which one then starts to
release new individuals from the natural reserve into the 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 an appropriate 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
distribution to be as close as possible to a significant target population,
chosen here to be the significant steady state $p_{\delta^{\prime}}$ of (2.1b)
for some $\delta^{\prime}\leq\delta_{1}$ where $\delta_{1}$ is defined in
(2.8). Recall that $\delta_{1}$ is a critical threshold that can be estimated
in practice and that guarantees that the KPPH equation possesses a significant
steady state (thus away from extinction) for any $\delta\leq\delta_{1}$.

We are limited by time in our action. We want indeed to save the population
from extinction that will occur at time $T$, corresponding to first time
instant at which the maximum of the solution equals
$\epsilon/\underline{\phi_{1}}$ if nothing is done222Because we assume
$\delta>\delta_{2}$, we know that such a time instant exists as the population
will become remnant; see Theorem 2.1-(vi).. 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)=\frac{1}{2}\int_{\tau}^{T}\big{|}y(t)-p_{\delta^{\prime}}\big{|}^{2}_{L^{2}(\Omega)}\,\mathrm{d}t+\frac{\kappa}{2}\int_{\tau}^{T}\big{|}u(t)\big{|}^{2}_{L^{2}(\Omega)}\,\mathrm{d}t.$
(2.12)

Here $y(t)$ denotes the solution in $L^{2}(\Omega)$ of the KPPH equation
forced by $u(t,x)$ and such that $y(t)=y_{0}$ at $t=\tau$, and $\kappa>0$. The
first term in the right-hand side (RHS) of (2.12) is to enforce closeness to
$p_{\delta^{\prime}}$, while the second term is an energy-type term related to
the effort of bringing the individuals from the reserve into the unprotected
area.

More precisely, we are aiming at addressing the following type of optimal
control problems

$\displaystyle\;\underset{u\in L^{2}(\tau,T;L^{2}(\Omega))}{\min}\;J(y,u)$
(2.13a) where $y(t,x)$ is the solution of
$\displaystyle\;\partial_{t}y=D\nabla^{2}y+\mu(x)y-\nu(x)y^{2}-\delta\rho_{\epsilon}(y)+u(t,x),\;(t,x)\in[\tau,T]\times\Omega,$
(2.13b) $\displaystyle\;\frac{\partial
y}{\partial{\boldsymbol{n}}}=0,\;(t,x)\in[\tau,T]\times\partial\Omega,$
(2.13c) with $y(t,x)=y_{0}(x)$ at $t=\tau$. (2.13d)

As explained below, we restrict to controls lying within a subset of
$L^{2}(\tau,T;L^{2}(\Omega))$, by introducing constraints on the control. To
describe these constraints we enter into more specificities about the model
setting. Throughout this article, we focus on a particular case of growth
rate, $\mu$, defined to be

$\mu(x)=m\chi_{\Lambda},\;\;m>0,$ (2.14)

where $\chi_{\Lambda}$ denotes the characteristic function of a (possibly
disconnected) subdomain $\Lambda$ of $\Omega$:

$\chi_{\Lambda}(x)=\begin{cases}1,&\text{ if }x\in\Lambda,\\\
0,&\text{otherwise}.\end{cases}$ (2.15)

Such a spatial dependence of the coefficient $\mu$ emphasizes that the
population reproduces only when the individuals are within the subdomain
$\Lambda$, while they may spread outside of $\Lambda$ due to the diffusion
term $D\nabla^{2}y$.

The choice of the control $u$ in (2.13d) is designed to act on this subdomain
to enhance chances of natural growth. For numerical applications we will
consider piecewise constant controls in space while allowing for fluctuations
in time. For that purpose, the domain $\Lambda$ is decomposed into mutually
disjoints subdomains such that

$\Lambda=\Lambda_{1}\cup\Lambda_{2}\cup\cdots\cup\Lambda_{K}.$ (2.16)

Then, by introducing

$\varphi_{j}=\frac{1}{\sqrt{|\Lambda_{j}|}}\chi_{\Lambda_{j}},\quad 1\leq
j\leq K,$ (2.17)

where $|\Lambda_{j}|$ denotes the $d$-dimensional volume of $\Lambda_{j}$, we
consider the following set of admissible controls:

$\displaystyle\mathcal{U}_{\textrm{ad}}\\!=\\!\bigg{\\{}\\!(t,x)\mapsto\sum_{j=1}^{K}\Gamma_{j}(t)\varphi_{j}(x)\,:\,\Gamma_{j}\in
L^{2}([\tau,T],\mathbb{R}),\;0\leq\Gamma_{j}(t)\leq C_{j}\text{ for a.e.
$t\in[\tau,T]$}\\!\bigg{\\}},$ (2.18)

where the $C_{j}$’s are positive constants which, as explained below, are
imposed by exploitation policy of the protected reserve. Denoting by
$\mathcal{H}=L^{2}(\Omega)$, the optimal control problem (2.13d) becomes, for
this set of admissible controls,

$\min J(y,u)\;\text{ subject to $(y,u)\in
L^{2}(\tau,T;\mathcal{H})\times\mathcal{U}_{\textrm{ad}}$ solving
\eqref{Eq_controlled_KPP1}-\eqref{Eq_controlled_KPP3}.}$ (2.19)

Due to the definition of $\mathcal{U}_{\textrm{ad}}$, this optimal control
problem is under constraints on the control333Note that the set
$\mathcal{U}_{\textrm{ad}}$ will be endowed with the induced topology from
that of $L^{2}(\tau,T;L^{2}(\Omega))$.. In the next section we will show that
this optimal control problem can be addressed efficiently by means of finite-
dimensional Galerkin approximations. For the moment we discuss how the
constraints $C_{j}$ in (2.18) arise from reserve management considerations.

######  Remark 2.1.

Note that for any given nonnegative continuous initial datum $y_{0}$ in
(2.13d), the existence of an optimal pair $(y(t;u^{*}),u^{*})$ for the optimal
control problem (2.19) can be ensured by using classical arguments such 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)=-\mu(x)y+\nu(x)y^{2}+\delta\rho_{\epsilon}(y)$ and $b\equiv 0$
therein. The assumptions collected into [Trö10, Assumption 5.6] are satisfied
here except the requirement that $d_{y}(x,t,y)\geq 0$. This latter condition
is used in [Trö10] to ensure the existence of solutions to [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]444In [Trö10,
Assumption 5.4], the condition $d_{y}(x,t,y)\geq 0$ is stated as $d$ being
monotonically increasing with respect to $y$.. However, for (2.13b)-(2.13d),
the existence of solutions and the desired uniform bound also hold as a
consequence of the maximum principle here; see Appendix B.

###  2.4. Control constraints from reserve management

We assume that the population density in the reserve is at its steady state
$\widehat{p}_{0}$, at the warning time $\tau$. The total population in the
reserve is thus given by

$M_{tot}=\int_{\widehat{\Omega}}\widehat{p}_{0}(x)\,\mathrm{d}x.$ (2.20)

A naive strategy, referred to as unsupervised below, consists of assuming
$\Gamma_{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
$\Omega$ from the reserve, between times $t=\tau$ and $t=T$, is

$M_{T}=\Big{(}\sum_{j=1}^{K}C_{j}\sqrt{|\Lambda_{j}|}\,\Big{)}(T-\tau).$
(2.21)

In practice, it is reasonable to assume that an exploitation policy holds for
the natural reserve independently of the rescue operation that is aimed to
take place for $\Omega$. As a result, only a fraction $\alpha$ ($0<\alpha<1$)
of the total population $M_{tot}$ in the reserve is allowed to be used for
populating another area such as $\Omega$.

An optimal control planning $u^{\ast}(t,x)$ solving the optimal control
problem (2.19) will be declared as efficient if it leads to an efficiency
ratio

$E=\frac{\int_{\tau}^{T}\int_{\Omega}u^{\ast}\,\mathrm{d}x\,\mathrm{d}t}{M_{tot}}<\frac{M_{T}}{M_{tot}}.$
(2.22)

Assuming that the unsupervised strategy would exploit the allowable fraction
$\alpha$, we have

$\Big{(}\sum_{j=1}^{K}C_{j}\sqrt{|\Lambda_{j}|}\,\Big{)}(T-\tau)=\alpha{M_{tot}}.$
(2.23)

Finally, we will assume in what follows that the $d$-dimensional volume
$|\Lambda_{j}|$ are the same for different $j$’s, and that
$C_{j}=C_{j^{\prime}}=C$ for $j\neq j^{\prime}$. This way, Eq. (2.23) gives
the value of the constraints $C_{j}$ in $\mathcal{U}_{\textrm{ad}}$. By doing
so we have thus an impartial way to compare the solution $u^{\ast}$ of the
optimal control problem (2.19), to that obtained from a naive, unsupervised
strategy consisting of taking $\Gamma_{j}(t)\equiv C$. We have however to keep
in mind that in this design operatus of the constraints, the set of admissible
controls, $\mathcal{U}_{\textrm{ad}}$ depends on $\alpha$, since the $C_{j}$
do. As a result, any optimal solution $u^{\ast}$ to (2.19) does also depend on
the allowable fraction $\alpha$ from the reserve.

From an exploitation viewpoint, we might ask what is the best fraction
$\alpha$ to use in order to maximize the productivity of new individuals in
the unprotected area, by the whole operation. An natural metric to answer this
question consists of calculating the population ratio

$P_{R}=\frac{\int_{\Omega}y(T,u^{\ast})\,\mathrm{d}x}{\int_{\Omega}p_{\delta^{\prime}}\,\mathrm{d}x}.$
(2.24)

The ratio $P_{R}$ allows us to assess which fraction of the targeted total
population has been obtained (when driven by the optimal control planning
$u^{\ast}$) when one reaches what would have been the extinction time $T$ if
no action would have been taken. A natural question is then to analyze how the
choice of $\alpha$ impacts $P_{R}$ and the efficiency ratio $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 extinction there, neither. For that, we can rely here again
on the theoretical understanding from the harvesting problem as recalled in
Sec. 2.2. Let us assume that the removal of individuals from the reserve
follows also a harvesting law in $\widehat{\Omega}$ of the form
$\sum_{j=1}^{K}C_{j}\varphi_{j}(x)\rho_{\epsilon}(\hat{y})$, where $\hat{y}$
denotes the population density in the reserve. If one wants
$u^{\ast}=\sum_{j=1}^{K}C_{j}\varphi_{j}$, one must have
$\hat{y}(t,x)\geq\epsilon/\underline{\phi_{1}}$ for $t$ in $[\tau,T]$ and $x$
in the reserve’s domain555We assume that the reserve’s habitat,
$\widehat{\Lambda}$, is decomposed also into $K$ mutually disjoints subdomains
each of same size than its corresponding counterpart $\Lambda_{j}$ in
$\Omega$. $\widehat{\Omega}$ to ensure $\hat{y}$ to be significant (see
(2.6)), imposing thus a constraint on the population on the reserve.

On the other hand, let $\hat{\delta}_{1}$ denotes the critical threshold
ensuring resilience of the population within the reserve this time; see
Theorem 2.1 again. A simple comparison argument shows that as soon as

$\sum_{j=1}^{K}C_{j}\varphi_{j}(x)\leq\hat{\delta}_{1}\bigg{(}\sum_{j=1}^{K}\varphi_{j}(x)\bigg{)},\;x\in\widehat{\Omega},$
(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

$\mbox{ For }1\leq j\leq K,\;\;C_{j}\leq\hat{\delta}_{1},$ (2.26)

since the $\widehat{\Lambda}_{j}$ are mutually disjoint. Thus the rescue
operation for the unprotected area must comply with the constraint (2.26), and
also from what precedes, with $\hat{y}(t,x)\geq\epsilon/\underline{\phi_{1}}$
for $t\in[\tau,T]$ and $x$ in the reserve’s domain $\widehat{\Omega}$. In
particular, (2.26) imposes a restriction on the population size $M_{T}$
excised from the reserve to satisfy
$M_{T}\leq\hat{\delta}_{1}(\sum_{j=1}^{K}\sqrt{|\widehat{\Lambda}_{j}|})(T-\tau).$
Such aspects regarding the reserve management will be also discussed in Sec. 4
below. For the moment we focus in the next section on the mathematical aspects
of solving the optimal control problem (2.19) via Galerkin approximations, for
a given set of admissible controls $\mathcal{U}_{\textrm{ad}}$ (and thus for a
given allowable fraction $\alpha$).

##  3\. Optimal control from Galerkin approximations

###  3.1. Convergence results and error estimates about the optimal control

In this section, we present convergence 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 approximations is then provided in Section 3.2, following a standard
Pontryagin Maximum Principle (PMP) approach.

Denote the spectral elements of the eigenvalue problem (2.3b) by
$\\{(\lambda_{j},e_{j})\,:\,j\in\mathbb{N}\\}$, where the eigenvalues
$\lambda_{j}$’s are arranged in an increasing order, and the eigenfunctions
$e_{j}$’s are normalized such that $\|e_{j}\|_{L^{2}(\Omega)}=1$. Denote also
the $N$-dimensional Galerkin approximation of the controlled state $y$ by

$y_{N}(t,x)=\sum_{i=1}^{N}\xi_{i}(t)e_{i}(x).$ (3.1)

Recall that the control $u$ is written as

$u(t,x)=\sum_{i=1}^{K}\Gamma_{i}(t)\varphi_{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 $\Pi_{N}y_{0}$, where
$\Pi_{N}$ denotes the projector onto the subspace spanned by the first $N$
eigenmodes solving (2.3b). Dependence of $y_{N}$ on $\Pi_{N}y_{0}$ or $u$ will
be made apparent depending on the context.

The Galerkin approximation of (2.13b)-(2.13c) reads then:

$\frac{\,\mathrm{d}\xi_{i}}{\,\mathrm{d}t}=-\lambda_{i}\xi_{i}+\sum_{j,k=1}^{N}B_{jk}^{i}\xi_{j}\xi_{k}-\delta\Big{\langle}\rho_{\epsilon}\Big{(}\sum_{j=1}^{N}\xi_{j}e_{j}\Big{)},e_{i}\Big{\rangle}+\langle
u(t,\cdot),e_{i}\rangle,\quad t\in[\tau,T],$ (3.3)

where $i=1,\ldots,N$, and

$B_{jk}^{i}=-\langle\nu(\cdot)e_{j}e_{k},e_{i}\rangle=-\int_{\Omega}\nu(x)e_{i}(x)e_{j}(x)e_{k}(x)\,\mathrm{d}x.$
(3.4)

Introducing a $K\times N$ matrix $\mathcal{M}$, whose elements are defined by

$\mathcal{M}_{ij}=\langle\varphi_{i},e_{j}\rangle=\int_{\Omega}\varphi_{i}(x)e_{j}(x)\,\mathrm{d}x,$
(3.5)

we can rewrite (3.3) as

$\frac{\,\mathrm{d}\xi_{i}}{\,\mathrm{d}t}=-\lambda_{i}\xi_{i}+\sum_{j,k=1}^{N}B_{jk}^{i}\xi_{j}\xi_{k}-\delta\Big{\langle}\rho_{\epsilon}\Big{(}\sum_{j=1}^{N}\xi_{j}e_{j}\Big{)},e_{i}\Big{\rangle}+\sum_{j=1}^{K}\mathcal{M}_{ji}\,\Gamma_{j}(t),\quad
t\in[\tau,T].$ (3.6)

The cost functional $J_{N}$ associated with the $N$-dimensional Galerkin
approximation (3.6) is

$J_{N}(y_{N},u)=\frac{1}{2}\int_{\tau}^{T}|y_{N}(t;\Pi_{N}y_{0},u)-\Pi_{N}p_{\delta^{\prime}}|_{L^{2}(\Omega)}^{2}\,\mathrm{d}t+\frac{\kappa}{2}\int_{\tau}^{T}|u(t)|_{L^{2}(\Omega)}^{2}\,\mathrm{d}t.$
(3.7)

By introducing

$\boldsymbol{\xi}=(\xi_{1},\ldots,\xi_{N})^{\textrm{T}},\textrm{ and
}\;\boldsymbol{\Gamma}=(\Gamma_{1},\ldots,\Gamma_{K})^{\textrm{T}},$ (3.8)

we rewrite the cost functional $J_{N}$ given by (3.7) as:

$J_{N}(\boldsymbol{\xi},\boldsymbol{\Gamma})=\frac{1}{2}\int_{\tau}^{T}|\boldsymbol{\xi}(t)-\boldsymbol{P}|^{2}\,\mathrm{d}t+\frac{\kappa}{2}\int_{\tau}^{T}|\boldsymbol{\Gamma}(t)|^{2}\,\mathrm{d}t,$
(3.9)

where $\boldsymbol{P}=(P_{1},\ldots,P_{N})^{\textrm{T}}$ with

$P_{i}=\langle p_{\delta^{\prime}},e_{i}\rangle,\quad i=1,\ldots,N.$ (3.10)

In connection to the set of admissible controls $\mathcal{U}_{\textrm{ad}}$
defined by (2.18), we introduce

$\displaystyle\mathcal{V}_{\mathrm{ad}}\\!=\\!\bigg{\\{}\\!t\mapsto(\Gamma_{1}(t),\ldots,\Gamma_{K}(t))^{\textrm{T}}\,:\,\Gamma_{j}\in
L^{2}([\tau,T],\mathbb{R}),\;0\leq\Gamma_{j}(t)\leq C_{j}\text{ for a.e.
$t\in[\tau,T]$}\\!\bigg{\\}}.$ (3.11)

The Galerkin approximation of the optimal control problem (2.19) is thus given
by

$\displaystyle\min J_{N}(\boldsymbol{\xi},\boldsymbol{\Gamma})$
$\displaystyle\text{ defined in \eqref{cost_JN} subject to
}(\boldsymbol{\xi},\boldsymbol{\Gamma})\in
L^{2}(\tau,T;\mathbb{R}^{N})\times\mathcal{V}_{\mathrm{ad}}\text{ solving}$
(3.12) the $N$-dimensional Galerkin system (3.6) with
$\boldsymbol{\xi}(\tau)=\Pi_{N}y_{0}$.

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_{0}$ in (2.13d) is strictly positive.

Then, there exist $C>0$ such that the optimal control $u^{\ast}$ for (2.19)
and the optimal control $u^{\ast}_{N}$ for the reduced problem (3.12) admit
the following error estimate:

$\int_{\tau}^{T}\big{|}u^{\ast}(t)-u^{\ast}_{N}(t)\big{|}_{L^{2}(\Omega)}^{2}\,\mathrm{d}t\leq
C\left[\sqrt{T-\tau}+T-\tau\right]\bigg{(}\sum_{j=1}^{2}\bigg{(}\int_{0}^{T}\varepsilon_{N}^{j}(t,u^{\ast},u_{N}^{\ast})\mathrm{d}t\bigg{)}^{\frac{1}{2}}\bigg{)},$
(3.13)

with

$\underset{t\in[\tau,T]}{\sup}\varepsilon_{N}^{j}(t,u^{\ast},u_{N}^{\ast})\xrightarrow{N\rightarrow\infty}0,\mbox{
for }\;j=1,2.$ (3.14)

Convergence results. Furthermore, the solution $y_{N}(t;\Pi_{N}y_{0},u)$ of
(3.6) converges uniformly to the solution $y(t;y_{0},u)$ of (2.13b)-(2.13d):

$\lim_{N\rightarrow\infty}\sup_{u\in\mathcal{U}_{\textrm{ad}}}\sup_{t\in[\tau,T]}\big{|}y_{N}(t;\Pi_{N}y_{0},u)-y(t;y_{0},u)\big{|}_{L^{2}(\Omega)}=0.$
(3.15)

In Appendix B we show that these convergence results and error estimates are
the consequence of general results of [CKL17] dealing with the Galerkin
approximations of nonlinear optimal control problems in Hilbert spaces 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 in Appendix A. Theorems A.1 and A.2 deal with the optimal
control of a broad class of nonlinear evolutionary equations containing the
KPPH equation considered in this article. The required assumptions of Theorem
A.1 and 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^{2}(\Omega)$) of
$\mathcal{G}(y)=\int_{\Omega}|y-p_{\delta^{\prime}}|^{2}\,\mathrm{d}x$ in a
neighborhood of the origin, and also on the growth rate of $J$ in some
appropriate norms of its arguments as given by (A.14) in Appendix A. The
precise expression of $\varepsilon_{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 to solve the reduced optimal control
problem (3.12) to obtain nearly optimal solutions as soon as $N$ is
sufficiently large. 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 above reduced optimal control problem (3.12); see [BH75,
BC03, Kir12]. Here, we follow an indirect approach by relying 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}(\boldsymbol{\xi},\boldsymbol{z},\boldsymbol{\Gamma})=\frac{1}{2}|\boldsymbol{\xi}-\boldsymbol{P}|^{2}+\frac{\kappa}{2}|\boldsymbol{\Gamma}(t)|^{2}+\boldsymbol{z}^{\textrm{T}}\boldsymbol{f}(\boldsymbol{\xi},\boldsymbol{\Gamma}),$
(3.16)

where $\boldsymbol{z}=(z_{1},\ldots,z_{N})^{\textrm{T}}$ denotes the costate
associated with the state $\boldsymbol{\xi}$, and $\boldsymbol{f}$ is the
Galerkin vector field whose $i$-th component is given by the RHS of (3.6).
That is,

$f_{i}(\boldsymbol{\xi},\boldsymbol{\Gamma})=-\lambda_{i}\xi_{i}+\sum_{j,k=1}^{N}B_{jk}^{i}\xi_{j}\xi_{k}-\delta\Big{\langle}\rho_{\epsilon}\Big{(}\sum_{j=1}^{N}\xi_{j}e_{j}\Big{)},e_{i}\Big{\rangle}+\sum_{j=1}^{K}\mathcal{M}_{ji}\,\Gamma_{j}(t),\quad
i=1,\ldots,N.$ (3.17)

Let

$(\boldsymbol{\xi}^{*},\boldsymbol{\Gamma}^{*})\in
L^{2}(\tau,T;\mathbb{R}^{N})\times\mathcal{V}_{\mathrm{ad}}$

be an optimal pair for the reduced optimal control problem (3.12), and denote
by $\boldsymbol{z}^{*}$ the costate associated with the state
$\boldsymbol{\xi}^{*}$. It follows from the PMP that the triplet
$(\boldsymbol{\xi}^{*},\boldsymbol{z}^{*},\boldsymbol{\Gamma}^{*})$ must
satisfy the following constrained Hamiltonian system (see e.g. [Kir12, Section
5.3]):

$\displaystyle\begin{rcases}\frac{\displaystyle\,\mathrm{d}\boldsymbol{\xi}^{*}}{\displaystyle\,\mathrm{d}t}=\nabla_{\boldsymbol{z}}H_{G}(\boldsymbol{\xi}^{*}(t),\boldsymbol{z}^{*}(t),\boldsymbol{\Gamma}^{*}(t))\\\
\frac{\displaystyle\,\mathrm{d}\boldsymbol{z}^{*}}{\displaystyle\,\mathrm{d}t}=-\nabla_{\boldsymbol{\xi}}H_{G}(\boldsymbol{\xi}^{*}(t),\boldsymbol{z}^{*}(t),\boldsymbol{\Gamma}^{*}(t))\end{rcases},\hskip
17.50002pt(\text{Hamiltonian system for
$(\boldsymbol{\xi}^{*},\boldsymbol{z}^{*})$})$ (3.18a) $\displaystyle
H_{G}(\boldsymbol{\xi}^{*}(t),\boldsymbol{z}^{*}(t),\boldsymbol{\Gamma}^{*}(t))\leq
H_{G}(\boldsymbol{\xi}^{*}(t),\boldsymbol{z}^{*}(t),\boldsymbol{\Gamma}(t)),\;\,\forall\;\boldsymbol{\Gamma}\in\mathcal{V}_{\mathrm{ad}},\;\;(\text{optimality
condition})$ (3.18b) $\displaystyle\boldsymbol{z}^{*}(T)=0,\hskip
200.0003pt(\text{terminal condition})$ (3.18c)

where (3.18a) and (3.18b) hold for all $t$ in $(\tau,T)$. Here $\nabla_{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}$. In general it is a necessary condition and not a
sufficient condition. In the general case, there may be controls that satisfy
(3.18b) but that are not optimal controls. Yet the PMP may delineate a
nonempty class of candidates. A triplet
$(\boldsymbol{\xi}^{*},\boldsymbol{z}^{*},\boldsymbol{\Gamma}^{*})$ solution
to (3.18) is called an extremal. Extremal solutions play an important role in
optimal control theory; see [BC03]. Sufficient conditions for extremal to
provide optimal controls can be found in [HSV95]. For the problem at hand,
since the cost functional $J_{N}$ in (3.9) is quadratic in
$\boldsymbol{\Gamma}$ and the dependence on the control is linear for the
control system (3.6), it is known that $\boldsymbol{\Gamma}^{\ast}$ obtained
from such an extremal is actually 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
$\boldsymbol{\Gamma}^{\ast}$ based on the optimal costate
$\boldsymbol{z}^{\ast}$. To do so, we first remark that

$-\frac{\partial
H_{G}(\boldsymbol{\xi},\boldsymbol{z},\boldsymbol{\Gamma})}{\partial\xi_{i}}=-(\xi_{i}-P_{i})-\sum_{j=1}^{N}z_{j}\frac{\partial
f_{j}(\boldsymbol{\xi},\boldsymbol{\Gamma})}{\partial\xi_{i}}.$ (3.19)

Now since

$\frac{\partial
f_{j}(\boldsymbol{\xi},\boldsymbol{\Gamma})}{\partial\xi_{i}}=-\lambda_{i}\delta_{ij}+\sum_{k=1}^{N}(B_{ik}^{j}+B_{ki}^{j})\xi_{k}-\delta\Big{\langle}\rho^{\prime}_{\epsilon}\Big{(}\sum_{k=1}^{N}\xi_{k}e_{k}\Big{)}e_{i},e_{j}\Big{\rangle},$
(3.20)

where $\delta_{ij}$ denotes the Kronecker delta, we get

$\displaystyle-\frac{\partial
H_{G}(\boldsymbol{\xi},\boldsymbol{z},\boldsymbol{\Gamma})}{\partial\xi_{i}}$
$\displaystyle=-(\xi_{i}-P_{i})$ (3.21)
$\displaystyle-\sum_{j=1}^{N}z_{j}\Big{(}-\lambda_{i}\delta_{ij}+\sum_{k=1}^{N}(B_{ik}^{j}+B_{ki}^{j})\xi_{k}-\delta\Big{\langle}\rho^{\prime}_{\epsilon}\Big{(}\sum_{k=1}^{N}\xi_{k}e_{k}\Big{)}e_{i},e_{j}\Big{\rangle}\Big{)}.$

Since (3.21) does not depend on $\boldsymbol{\Gamma}$, we denote $-\partial
H_{G}(\boldsymbol{\xi},\boldsymbol{z},\boldsymbol{\Gamma})/\partial\xi_{i}$ by
a function $g_{i}(\boldsymbol{\xi},\boldsymbol{z})$. The costate equation
becomes then
$\,\mathrm{d}\boldsymbol{z}^{*}/\,\mathrm{d}t=\boldsymbol{g}(\boldsymbol{\xi}^{*},\boldsymbol{z}^{*})$.
Thus the Hamiltonian system (3.18a) together with boundary condition (3.18c)
becomes

$\displaystyle\frac{\displaystyle\,\mathrm{d}\boldsymbol{\xi}^{*}}{\displaystyle\,\mathrm{d}t}=\boldsymbol{f}(\boldsymbol{\xi}^{*},\boldsymbol{\Gamma}^{*}),\quad
t\in(\tau,T),$ (3.22)
$\displaystyle\frac{\displaystyle\,\mathrm{d}\boldsymbol{z}^{*}}{\displaystyle\,\mathrm{d}t}=\boldsymbol{g}(\boldsymbol{\xi}^{*},\boldsymbol{z}^{*}),\quad
t\in(\tau,T),$
$\displaystyle\boldsymbol{\xi}^{*}(\tau)=\Pi_{N}y_{0},\quad\boldsymbol{z}^{*}(T)=0,$

where the components of $\boldsymbol{f}$ and $\boldsymbol{g}$ are defined in
the RHSs of (3.17) and (3.21), respectively. The condition
$\boldsymbol{\xi}^{*}(\tau)=\Pi_{N}y_{0}$ is obtained by projecting the
initial condition $y(\tau,x)=y_{0}$ in (2.13d).

We show next how $\boldsymbol{\Gamma}^{*}$ depends on the costate
$\boldsymbol{z}^{*}$. To do so, we first remark that by using $\mathcal{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

$\frac{\kappa}{2}|\boldsymbol{\Gamma}^{*}(t)|^{2}+(\boldsymbol{z}^{*}(t))^{\textrm{T}}\mathcal{M}^{\textrm{T}}\boldsymbol{\Gamma}^{*}(t)\leq\frac{\kappa}{2}|\boldsymbol{\Gamma}(t)|^{2}+(\boldsymbol{z}^{*}(t))^{\textrm{T}}\mathcal{M}^{\textrm{T}}\boldsymbol{\Gamma}(t),\;\,\forall\;\boldsymbol{\Gamma}\in\mathcal{V}_{\mathrm{ad}},\;t\in(\tau,T).$
(3.23)

Then by introducing

$\boldsymbol{w}^{*}=\mathcal{M}\boldsymbol{z}^{*},$ (3.24)

we get

$\displaystyle\frac{\kappa}{2}|\boldsymbol{\Gamma}(t)|^{2}+(\boldsymbol{z}^{*}(t))^{\textrm{T}}\mathcal{M}^{\textrm{T}}\boldsymbol{\Gamma}(t)$
$\displaystyle=\frac{\kappa}{2}|\boldsymbol{\Gamma}(t)|^{2}+(\boldsymbol{w}^{*}(t))^{\textrm{T}}\boldsymbol{\Gamma}(t)$
(3.25)
$\displaystyle=\frac{\kappa}{2}\sum_{j=1}^{K}\Big{(}\Gamma_{j}(t)+\frac{1}{\kappa}w^{*}_{j}(t)\Big{)}^{2}-\frac{1}{2\kappa}\sum_{j=1}^{K}w^{*}_{j}(t)^{2}.$

As a consequence, the optimality condition (3.23) becomes

$\sum_{j=1}^{K}\Big{(}\Gamma^{*}_{j}(t)+\frac{1}{\kappa}w^{*}_{j}(t)\Big{)}^{2}\leq\sum_{j=1}^{K}\Big{(}\Gamma_{j}(t)+\frac{1}{\kappa}w^{*}_{j}(t)\Big{)}^{2},\;\,\forall\;\boldsymbol{\Gamma}\in\mathcal{V}_{\mathrm{ad}},\;t\in(\tau,T).$
(3.26)

Recalling the control constraints $0\leq\Gamma_{j}\leq C_{j}$ resulting from
(3.11), one obtains that (3.26) holds if and only if, for all $j=1,\ldots,K$,

$\Gamma^{*}_{j}(t)=h_{j}(\boldsymbol{w}^{*}(t))\stackrel{{\scriptstyle\textrm{def}}}{{=}}\begin{cases}0,&\text{
if $-w^{*}_{j}(t)<0$},\\\ -\frac{1}{\kappa}w^{*}_{j}(t),&\text{ if
$0\leq-w^{*}_{j}(t)\leq\kappa C_{j}$},\\\ C_{j},&\text{ if
$-w^{*}_{j}(t)>\kappa C_{j}$},\end{cases}\qquad\;\forall\;t\in(\tau,T).$
(3.27)

Now by substituting the expression of $\boldsymbol{\Gamma}^{*}$ thus obtained
and by using (3.24), the problem (LABEL:Eq_BVP_for_PMP) reduces to a boundary
value problem (BVP) in the variables $\boldsymbol{\xi}^{*}$ and
$\boldsymbol{z}^{*}$. The synthesis of an optimal control
$\boldsymbol{\Gamma}^{*}$ boils down thus to solving the following BVP:

$\displaystyle\frac{\displaystyle\,\mathrm{d}\boldsymbol{\xi}^{*}}{\displaystyle\,\mathrm{d}t}=\boldsymbol{f}(\boldsymbol{\xi}^{*},\boldsymbol{h}(\mathcal{M\boldsymbol{z}^{*}})),\quad
t\in(\tau,T),$ (3.28)
$\displaystyle\frac{\displaystyle\,\mathrm{d}\boldsymbol{z}^{*}}{\displaystyle\,\mathrm{d}t}=\boldsymbol{g}(\boldsymbol{\xi}^{*},\boldsymbol{z}^{*}),\quad
t\in(\tau,T),$
$\displaystyle\boldsymbol{\xi}^{*}(\tau)=\Pi_{N}y_{0},\quad\boldsymbol{z}^{*}(T)=0,$

where components of the function $\boldsymbol{h}$ are defined in (3.27).

The optimal control $u^{*}_{N}$ to (3.12) in $\mathcal{U}_{\textrm{ad}}$ (see
(2.18)) is finally given by

$u^{*}_{N}(t,x)=\sum_{i=1}^{K}\Gamma^{*}_{i}(t)\varphi_{i}(x),$ (3.29)

where $\boldsymbol{\Gamma}^{*}$ is given by (3.27) with $\boldsymbol{w}^{*}$
therein obtained as $\boldsymbol{w}^{*}=\mathcal{M}\boldsymbol{z}^{*}$ for
$\boldsymbol{z}^{*}$ solving (LABEL:Eq_BVP_for_PMP2).

######  Remark 3.2.

Note that, if the constraints $0\leq\Gamma_{j}\leq C_{j}$ on the admissible
control are removed in $\mathcal{V}_{\mathrm{ad}}$ in (3.11), the optimality
condition (3.18b) can be equivalently written as (see again [Kir12, Section
5.3]):

$\nabla_{\boldsymbol{\Gamma}}H_{G}(\boldsymbol{\xi}^{*},\boldsymbol{z}^{*},\boldsymbol{\Gamma}^{*})=0.$
(3.30)

The control law for the constrained case is given by (3.27). In contrast, one
obtains from (3.30) the following control law for the unconstrained case:

$\boldsymbol{\Gamma}^{*}=-\frac{1}{\kappa}\mathcal{M}\boldsymbol{z}^{*},$
(3.31)

with $\mathcal{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
results

###  4.1. Numerical setup

To illustrate our theoretical framework, and to simplify the reproducibility
of the numerical 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
$\Omega=(0,\ell)$, with $\ell>0$. To simplify the analysis, we choose the
domain $\widehat{\Omega}$ of the reserve to be also given by an interval of
same length.

We set $D=1$, $\ell=1$, and $\nu=0.2$ in the corresponding KPPH model (2.1b).
To account for the possible effects of heterogeneity of the habitat, we
consider two cases of subdomain, $\Lambda$, appearing in the definition of
$\mu$ given by (2.14):

Case I: $\displaystyle\qquad\Lambda=[1/2,1],$ (4.1) Case II:
$\displaystyle\qquad\Lambda=[0,1/4]\cup[3/4,1].$

Note that $|\Lambda|=0.5$ for both cases, but Case I corresponds to an habitat
more aggregated that Case II. In each case we set $m=2$ in (2.14).

Regarding the set $\mathcal{U}_{\textrm{ad}}$ of admissible controls (see
(2.18)), we divide $\Lambda$ into $K=8$ segments of equal length, such that
each segment is of length $1/16$.

Following Sec. 2.3, we assume that for each case the population evolving in
the unprotected area according to (2.1b) is under extinction threat because
$\delta=\delta_{2}(1+f)$, with $f$ chosen to be equal to $0.1$. One aims at
restoring the population to a safe dynamics leading towards a significant
steady state $p_{\delta^{\prime}}$ appearing in the cost functional $J$ given
by (2.12) for some $\delta^{\prime}<\delta_{1}$.666We recall, that once
$\delta^{\prime}\leq\delta_{1}$, there exists indeed a significant steady
state, $p_{\delta^{\prime}}$, due to Theorem 2.1-($iii$). For that purpose we
choose $\delta^{\prime}=\delta_{1}(1-f)$, also with $f=0.1$. Recall that
$\delta_{1}$ and $\delta_{2}$ are critical harvesting thresholds defined in
(2.8); see Sec. 2.2.

The harvest function $\rho_{\epsilon}$ in (2.2) is explicitly given here for
$\epsilon=0.05$ by

$\rho_{\epsilon}(x)=\begin{cases}1,&\text{if $x\geq\epsilon$},\\\
0.5\sin(\pi(x-0.5\epsilon)/\epsilon)+0.5,&\text{if $0<x<\epsilon$},\\\
0,&\text{otherwise}.\end{cases}$ (4.2)

Recall that the optimal control problem (2.19) is carried out over the time
interval $[\tau,T]$. Following Sec. 2.3, given a solution $y$ of the KPPH
model (2.1b) (emanating from $p_{0}$), the warning time $\tau$ is chosen
according to $\int_{0}^{\ell}y(\tau,x)\,\mathrm{d}x=P_{c}$ with $P_{c}$
defined in (2.10) for some $0<\beta<1$. We choose $\beta=1/4$ in all the
numerical experiments. In particular $P_{c}$ satisfies (2.11) for Cases I and
II. In each case, the extinction time $T$ corresponds to the first time
instant at which the maximum (over $(0,\ell)$) of the solution $y(t,x)$ equals
$\epsilon/\underline{\phi_{1}}$. Note that $y_{0}$ used in the optimal control
problem (2.19) is $y(\tau,\cdot)$ (see (2.13d)).

The goal is to drive, over the time window $[\tau,T]$, the population
distribution to be as close as possible to $p_{\delta^{\prime}}$, while
minimizing the cost functional given by (2.12) for controls $u$ lying in the
admissible set $\mathcal{U}_{\textrm{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 $\Lambda$ given by (4.1) (Cases I and II) are listed in Table 1,
rounded to the nearest thousandth. Note that $\delta_{1}$ and $\delta_{2}$
defined in (2.8) depend in each case on the first eigenmode $e_{1}$ solving
(2.3b) (associated with the first eigenvalue $\lambda_{1}$)777Note that
$e_{1}$ is normalized in $L^{2}(\Omega)$ whereas $\phi_{1}$ appearing in (2.8)
is normalized in $L^{\infty}(\Omega)$. The two are related according to
$\phi_{1}=e_{1}/\max(e_{1})$. shown here in Fig. 1.

Figure 1. First eigenmode of the spectral problem (2.3b) for $\Omega=(0,1)$
shown here for $\Lambda=[1/2,1]$ (Case I) and $\Lambda=[0,1/4]\cup[3/4,1]$
(Case II).

In what follows, the KPPH model is solved using the Matlab solver pdepe for
$\Omega=(0,1)$ with $\delta x=10^{-2}$ and $\delta t=10^{-2}$. The
corresponding spectral problem (2.3b) is solved using finite difference with
$\delta x=2^{-12}$. We use a small mesh grid size here to reach good accuracy
given the discontinuity of the coefficient $\mu$. In both parameter cases, the
dominant eigenvalue $\lambda_{1}$ is negative (see Table 1), corresponding to
an unstable eigenmode $\phi_{1}$. The rest of the spectrum corresponds here to
stable eigenmodes.

Table 1. Parameter values in (2.13d) | $\lambda_{1}$ | $\delta_{1}$ | $\delta_{2}$ | $\delta^{\prime}$ | $\delta$ | $\tau$ | $T$
---|---|---|---|---|---|---|---
Case I | -1.083 | 1.443 | 1.465 | 1.298 | 1.612 | 12.544 | 14.009
Case II | -1.021 | 1.301 | 1.302 | 1.171 | 1.432 | 14.023 | 15.582

###  4.2. Choice of the control constraints

We assume that the reserve’s domain size $|\widehat{\Omega}|$ is equal to
$\ell$ and that the population’s habitat in the reserve $\widehat{\Lambda}$
satisfies $|\widehat{\Lambda}|=|\Lambda|.$ To simplify the analysis, we also
assume that the KPP model’s parameters in the reserve $\widehat{\Omega}$ are
the same as those of the KPPH equation in the unprotected area. Under these
working assumptions, the total population in the reserve is given by

$M_{tot}=\int_{\widehat{\Omega}}\widehat{p_{0}}(x)\,\mathrm{d}x,$ (4.3)

where $\widehat{p_{0}}$ is the steady state of the KPP model (under Neumann
condition) in $\widehat{\Omega}$. In our numerical experiments, $M_{tot}=5.34$
for Case I and $M_{tot}=5.1$ for Case II.

As pointed out in Sec. 2.4, the control planning $u(t,x)$ of transporting
individuals from the natural reserve to the unprotected area is subordinated
to the following reserve management factors:

  * A)

The choice of the population fraction $\alpha$ allowed to be transported from
the reserve.

  * B)

The population displaced from the reserved must be such that the population
distribution in the reserve satisfies
$\hat{y}(t,x)\geq\epsilon/\underline{\phi_{1}}$ for $t$ in $[\tau,T]$ and $x$
in $\widehat{\Omega}$.

  * C)

For $1\leq j\leq K,$ the constant $C_{j}$ must satisfy
$C_{j}\leq\hat{\delta}_{1},$ where $\hat{\delta}_{1}$ is the critical
harvesting threshold favoring 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 $\Lambda_{j}$. Then (2.23) allows us to
find $C$ given by

$C=\frac{\alpha{M_{tot}}}{2(T-\tau)},$ (4.4)

because $K=8$ and $|\Lambda_{j}|=1/16$. Note although the habitat dependence
is not directly apparent in (4.4), the constant $C$ depends actually on
$\Lambda$ through the coefficient $\mu$ impacting the dynamics and thus the
extinction time $T$ and the warning time $\tau$. Given our numerical setup,
for the constant $C$ chosen according to (4.4) then the requirement C) above
is satisfied as long as $0<\alpha\leq\alpha_{c}$ with $\alpha_{c}=0.792$ for
Case I and $\alpha_{c}=0.796$ for Case II. 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 it by a Galerkin truncation of
minimal dimension, namely we choose $N=1$. The reason is that for the model
parameters 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_{\delta^{\prime}}$ in
each case; see Fig. 2. The insets of this figure show indeed that these target
population densities are highly correlated with the corresponding first
eigenmodes shown in Fig. 1.

Figure 2. Energy decomposition of the target population density for
$\Lambda=[1/2,1]$ (Case I) and $\Lambda=[0,1/4]\cup[3/4,1]$ (Case II) shown in
semilog-scale. The target population density (shown in the upper insets)
corresponds in each case to the significant steady states
$p_{\delta^{\prime}}$ of the KPPH equation (2.1b) for
$\delta=\delta^{\prime}$, with $\delta^{\prime}$ as indicated in Table 1.
These steady states are highly correlated with the corresponding first
eigenmodes shown in Fig. 1, explaining that these eigenmodes capture most of
the energy contained in these steady states.

For $N=1$, the Galerkin approximation (3.6) reduces to the following scalar
logistic equation with harvesting and forcing terms:

$\frac{\,\mathrm{d}X}{\,\mathrm{d}t}=X(a-bX)-\delta
g_{\epsilon}(X)+\sum_{j=1}^{K}d_{j}\,\Gamma_{j}(t),\quad
t\in[\tau,T],\;\;K=8,$ (4.5)

where $a=-\lambda_{1}>0$, $b=-B_{11}^{1}$ (see (3.4)) and the
$d_{j}=\mathcal{M}_{j1}$ ($N=1$) are constants given by (3.5), thus obtained
here by the inner product between the eigenmode $e_{1}$ and the characteristic
functions of the subdomains $\Lambda_{j}$. The harvesting term $g_{\epsilon}$
is given by

$g_{\epsilon}(X)=\begin{cases}1&\text{ if
}X\geq\frac{\epsilon}{\min{e_{1}}},\\\
\Big{\langle}\rho_{\epsilon}\Big{(}Xe_{1}\Big{)},e_{1}\Big{\rangle}&\text{if
$X<\frac{\epsilon}{\min{e_{1}}}$},\end{cases}$ (4.6)

because the first eigenmode $e_{1}$ is positive. The harvest function
$g_{\epsilon}$ satisfies the conditions in (2.2), and thus the Galerkin
projection preserves the global structure of the original KPPH equation.

Due to the high-energy content captured by the first-mode amplitude, $X(t)$,
we expect that solving the optimal control (3.12) associated with (4.5) (i.e.
for $N=1$) should provide a nearly optimal solution for the original optimal
control problem (2.19) with $\Omega=(0,\ell)$. 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 obtained according to (3.29). The optimal
control $\Gamma_{j}^{\ast}$ therein is obtained from (3.27), namely
$\Gamma^{*}_{j}(t)=h_{j}(\boldsymbol{w}^{*}(t))$, where
$\boldsymbol{w}^{*}=z^{\ast}\mathcal{M}$ with $z^{\ast}$ solving here the
following BVP over $[\tau,T]$ (from (LABEL:Eq_BVP_for_PMP2))

$\displaystyle\frac{\displaystyle\,\mathrm{d}X^{*}}{\displaystyle\,\mathrm{d}t}=X^{*}(a-bX^{*})-\delta
g_{\epsilon}(X^{*})+\sum_{j=1}^{K}d_{j}h_{j}(z^{\ast}\mathcal{M})\,,$ (4.7)
$\displaystyle\frac{\displaystyle\,\mathrm{d}z^{*}}{\displaystyle\,\mathrm{d}t}=-(X^{\ast}-\langle
p_{\delta^{\prime}},e_{1}\rangle)-z^{\ast}(a-2bX^{\ast})-\delta
k_{\epsilon}(X^{\ast}),$ $\displaystyle X^{*}(\tau)=\langle
y_{0},e_{1}\rangle,\quad z^{*}(T)=0,$

with

$k_{\epsilon}(X^{\ast})=\Big{\langle}\rho^{\prime}_{\epsilon}\Big{(}X^{\ast}e_{1}\Big{)}e_{1},e_{1}\Big{\rangle}.$
(4.8)

Note that here $\mathcal{M}$ given by (3.5) is a column vector and the costate
$z^{\ast}$ is a scalar variable.

In what follows, higher $N$-dimensional Galerkin truncation (with $N>1$) are
also used to control the KPPH equation. The use of such higher-dimensional
Galerkin approximations is to benchmark the one-dimensional approximation. The
corresponding optimal control $u_{N}^{\ast}$ is, for each $N$, obtained by
following the PMP approach as described in Sec. 3.2.

Table 2. Coefficients of the 1D Galerkin approximation (4.5) | $a$ | $b$ | $\delta$ | $d_{1}$ | $d_{2}$ | $d_{3}$ | $d_{4}$ | $d_{5}$ | $d_{6}$ | $d_{7}$ | $d_{8}$
---|---|---|---|---|---|---|---|---|---|---|---
Case I | 1.083 | 0.202 | 1.612 | 0.252 | 0.259 | 0.265 | 0.270 | 0.274 | 0.277 | 0.279 | 0.280
Case II | 1.021 | 0.200 | 1.432 | 0.257 | 0.257 | 0.255 | 0.252 | 0.252 | 0.255 | 0.257 | 0.257

###  4.4. First numerical results

We focus here on two particular choices of the allowable fraction $\alpha$
($\alpha=0.1$ and $\alpha=0.35$) of individuals transported from the reserve’s
population to the unprotected area $\Omega,$ for each case of habitat
$\Lambda$ in (4.1). Recall that the discrimination between the two cases of
habitats considered in Sec. 4.1 is aimed 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}(\boldsymbol{\xi}^{\ast},\boldsymbol{\Gamma}_{N}^{*})$) it corresponds
to the value of $J$ defined in (2.12) (resp. $J_{N}$ defined in (3.9)) when
the solution $y$ (resp. $\boldsymbol{\xi}^{\ast}$) 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}^{\ast}$ (resp.
$\boldsymbol{\Gamma}_{N}^{\ast}$).

For both cases of habitat (Cases I and II) and for both $\alpha$-values
($\alpha=0.1$ and $\alpha=0.35$), we have computed the relative error in the
cost values $J$ between the optimal solutions obtained from the 1D and the 10D
Galerkin approximations (see Table 3), namely

$\texttt{error}=\frac{|J(y,u_{1}^{*})-J(y,u_{10}^{*})|}{J(y,u_{10}^{*})}\times
100\%.$ (4.9)

We found that in all cases, this relative error is almost negligible, bounded
by $1.5\times 10^{-4}\%$.

Furthermore, the convergence of the cost value
$J_{N}(\boldsymbol{\xi}^{\ast},\boldsymbol{\Gamma}_{N}^{*})$ is achieved
quickly as the first few digits in the cost value has already converged when
$N$ is increased; see Table 4. These numerical results confirm thus the
intuition that due to the high-energy content carried out by the first
eigenmode, a 1D Galerkin approximation is sufficient to obtain nearly optimal
solutions.

Table 3. Cost value $J(y,u_{N}^{\ast})$ | ($\alpha=0.1$, $N=1$) | ($\alpha=0.1$, $N=10$) | ($\alpha=0.35$, $N=1$) | ($\alpha=0.35$, $N=10$)
---|---|---|---|---
Case I | 4.168706 | 4.168700 | 3.203817 | 3.203754
Case II | 3.952086 | 3.952086 | 2.996004 | 2.995997
Table 4. Cost value $J_{N}(\boldsymbol{\xi}^{\ast},\boldsymbol{\Gamma}_{N}^{\ast})$ | $N=1$ | $N=2$ | $N=3$ | $N=4$ | $N=5$ | $N=10$
---|---|---|---|---|---|---
Case I, $\alpha=0.1$ | 4.159675 | 4.168644 | 4.168642 | 4.168651 | 4.168651 | 4.168651
Case II, $\alpha=0.1$ | 3.953014 | 3.953014 | 3.953425 | 3.953425 | 3.953425 | 3.953426
Case I, $\alpha=0.35$ | 3.185567 | 3.203669 | 3.203668 | 3.203694 | 3.203693 | 3.203690
Case II, $\alpha=0.35$ | 2.996189 | 2.996193 | 2.997253 | 2.997253 | 2.997253 | 2.997254

Based on these results, we are in a comfortable position to discuss the
numerical KPPH solution $y(t,x;u_{N}^{\ast})$ obtained when driven by the
optimal control $u_{N}^{\ast}$ itself obtained according to (3.29) (for $N=1$)
and the BVP (LABEL:Eq_BVP_for_PMP3). Figures 3 and 4 show these solutions in
their respective panels (c) and (f) depending on the Case of habitat or value
of $\alpha$ considered. The corresponding panels (a) (or (d)) show the optimal
control $u_{N}^{\ast}(t,x)$ (for $N=1$) for different time instants while
panels (b) (or (e)) show its time evolution, after integration over
$\Omega=(0,\ell)$, denoted by $\langle u_{N}^{\ast}(t)\rangle$.

Figure 3. Panel (a): Optimal control $u_{N}^{\ast}(t,x)$ (Case I,
$\alpha=0.35$) given by (3.29) with $N=1$. Here
$\boldsymbol{\Gamma}_{N}^{\ast}$ is obtained from (3.27), where
$\boldsymbol{w}^{*}=z^{\ast}\mathcal{M}$ with $z^{\ast}$ solving here the BVP
(LABEL:Eq_BVP_for_PMP3). Panel (b): Time-dependence of $\langle
u_{N}^{\ast}(t)\rangle$ after space integration over $(0,\ell)$. The first
time instant at which $\langle u_{N}^{\ast}(t)\rangle>0$ corresponds to
$t=\tau$. The next time instant at which $\langle u_{N}^{\ast}(t)\rangle=0$
corresponds to $t=T$. Panel (c): Controlled KPPH solution
$y(t,x;u_{N}^{\ast})$ to (2.13b)-(2.13d) when driven by $u_{N}^{\ast}$. Panel
(d) same as Panel (a), Panel (e) same as Panel (b), and Panel (f) same as
Panel (c), for Case II.

Recall that given a value of $\alpha$, the constant $C$ is determined from
(4.4) which determines in turn the control constraints in the set
$\mathcal{U}_{\textrm{ad}}$ of admissible controls, as the $C_{j}$ therein are
chosen to be equal to $C$; see Sec. 4.2. The impact of $\alpha$ is as follows.
The smaller the allowable fraction $\alpha$ is, the smaller $C$ is and the
more constrained is the optimal control problem (2.19). This is visible on the
temporal evolution of $\langle u_{N}^{\ast}(t)\rangle$ which is set to $C$
over a longer interval for $\alpha=0.1$ than for $\alpha=0.35$, before
decreasing to zero; compare Panels (b) and (e) of Fig. 4 with those of Fig. 3.
We conduct a more detailed analysis on the dependence on $\alpha$ in Sec. 4.5
below. For the moment, we discuss the effects of fragmentation of the habitat
as the latter has been shown to play an important role in the population
resilience to external perturbations [RC07].

Figure 4. Same as Fig. 3 but for $\alpha=0.1$. Table 5. The efficiency ratio $E$ defined by (2.22) | $\alpha=0.1$ | $\alpha=0.35$
---|---|---
Case I | 0.0912 | 0.1920
Case II | 0.0918 | 0.2002

The efficiency ratio $E$ defined in (2.22) serves us to compare the effects of
different configurations of the habitat. Recall that the smaller is this ratio
the better it is in terms of exploitation of the reserve’s population but not
necessarily for saving the population from extinction in the unprotected area.
For that latter aspect one should have the population ratio $P_{R}$ (defined
in (2.24)) as large as possible while keeping in mind to respect the
requirements 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 $\alpha$ of the reserve, the habitat that is more homogeneous (Case
I) is slightly more advantageous than the less homogeneous one (Case II). In
terms of population ratio $P_{R}$ defined by (2.24), it is Case II that is
slightly more advantageous, for instance $P_{R}=0.5633$ for Case II vs
$P_{R}=0.5581$ for Case I, when $\alpha=0.35$. This means that for Case II,
$56.33\%$ of the targeted population size is reached at time $t=T$, and
$55.81\%$ for Case I. Thus, the differences across the habitat fragmentation
are somewhat minor. To the opposite, the difference in terms of allowable
fraction $\alpha$ of the reserve seem important as only of about $28\%$ (in
both cases) of the targeted population size is reached at time $t=T$ when
$\alpha=0.1$. 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 $\alpha=0.1$, we observe in Panels (c) and (f) of Fig. 4 that the final
profile $y(T;u_{N}^{\ast})$ of the controlled solution satisfies
$y(T;u_{N}^{\ast})<y_{0}$ ($y_{0}$ corresponds to the blue curve in Panels (c)
and (f)), irrespectively of the degree of fragmentation of $\Lambda$. As time
goes beyond $t=T$ and the harvesting intensity $\delta$ is set back to
$\delta^{\prime}$ in the KPPH equation while the control is
abandoned888Assuming we get rid of illegal harvesting by this time instant,
while still allowing some harvesting respecting the quota
$\delta^{\prime}<\delta_{1}$., the population 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
$\alpha=0.35$, $y(T;u_{N}^{\ast})>y_{0}$, and one gets in either Case I or II,
closer than for $\alpha=0.1$ to the target population, i.e. the significant
steady state $p_{\delta^{\prime}}$. As time goes beyond $t=T$ and $\delta$ is
set to $\delta^{\prime}$, 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 $\alpha=0.1$ and
$\alpha=0.35$? Remember that due to (4.4) and our protocol for choosing our
set $\mathcal{U}_{\textrm{ad}}$ of admissible controls, the allowable fraction
$\alpha$ impacts the constraint $C$ arising in $\mathcal{U}_{\textrm{ad}}$.
The smaller $\alpha$, the smaller $C$ and the more constrained becomes our
control planning. At the same time it should be noted that when $\alpha$ is
increased further towards $1$, the constraint $C$ becomes too large so that
this constraint is never activated and the optimal control problem (2.19)
behaves as unconstrained.999It is easy to convince oneself of this fact by
remarking that $\boldsymbol{\Gamma}^{\ast}$ given by (3.27) reduces to (3.31)
when the $C_{j}$ are sufficiently large; see Remark 3.2. Thus there exists a
critical value of alpha, say $\overline{\alpha}$, above which the control
bounds are not activated in the set $\mathcal{U}_{\textrm{ad}}$ of admissible
controls. For Case I we found $\overline{\alpha}=0.49$ and
$\overline{\alpha}=0.5$ for Case II.

As a consequence, for $\alpha>\overline{\alpha}$, only the same unconstrained
optimal solution is obtained and the efficiency ratio $E$ saturates to a
constant value. Panel (b) of Fig. 5 shows the dependence of $E$ in terms of
$\alpha$, and in particular its saturation for $\alpha>\overline{\alpha}$.
Only Case I is shown here as $E$ for Case II behaves almost identically.

Figure 5. Panel (a): Population ratio $P_{R}$ (Case I). Panel (b): Efficiency
ratio $E$ (Case I). The value $\overline{\alpha}$ corresponds to the value of
$\alpha$ above which the control bounds are not activated. The value
$\underline{\alpha}$ corresponds to the value of $\alpha$ below which the
population in $\Omega$ evolves towards a remnant state when the control is
abandoned and $\delta$ is set to $\delta^{\prime}$ after $t=T$; see (4.11)
below.

The dependence on $\alpha$ of the population ratio $P_{R}$ is shown in Panel
(a) of Fig. 5. We observe that the $P_{R}$\- and $E$-curves are highly
correlated, both increasing and saturating for $\alpha>\overline{\alpha}$.
Thus, “there is no free lunch,” and one cannot minimize $E$ while maximizing
$P_{R}$. The choice of $\alpha$ is however a determining factor in the success
of the rescue operation as pointed above. In all the cases, the excision of
population to be transported from the reserve to the unprotected area can be
absorbed by the reserve’s population. Figure 6 shows indeed that even for the
more demanding “excision pressure” on the population’s reserve (corresponding
to the larger $E$ achieved for $\alpha=\overline{\alpha}$ and beyond), the
population in the reserve recovers it original steady state, independently on
the habitat’s fragmentation.

Figure 6. Panel (a): Population size $\langle\hat{y}(t)\rangle$ in the reserve
(Case I, $\alpha=\overline{\alpha}$). The bracket indicates integration over
space. The dynamics for the population in the natural reserve is governed by
(2.1b) forced by $-u^{*}_{N}$, corresponding to the amount of individual
transported from the reserve to the unprotected area. The inset shows $\langle
u_{N}^{\ast}(t)\rangle$ over a time window slightly larger than $[\tau,T]$.
The first vertical dashed line indicates $t=\tau$, while the second one,
indicates $t=T$. Panel (b): Same for Case II, $\alpha=\overline{\alpha}$.

The critical role of $\alpha$ on the survival of the population in the
unprotected after application to optimal planning can gain great insights from
the reduced system (4.5) and the dynamical properties of the KPPH equation as
recalled in Sec. 2.2. Let us assume that once the optimal control
$u_{N}^{\ast}$ has been obtained from the BVP (LABEL:Eq_BVP_for_PMP3) for a
given $\alpha$, the harvesting intensity $\delta$ is set back to
$\delta^{\prime}$ in the KPPH equation, for $t>T$, and the control is
abandoned. The idea is that by assuming that the illegal harvesting that was
putting the population under extinction threat has been stopped starting from
$t=T$, one aims at anticipating in terms of $\alpha$ (thus on the constraint
put on the control) whether the rescue operation that would take place over
the time window $[\tau,T]$, would be a success or not.

To do so, let us consider the following modification of (4.5)

$\frac{\displaystyle\,\mathrm{d}X}{\displaystyle\,\mathrm{d}t}=X(a-bX)-\delta^{\prime}g_{\epsilon}(X),$
(4.10)

which is nothing else than the projection of the KPPH equation (with
$\delta=\delta^{\prime}$) onto the first eigenmode $e_{1}$. Let
$\underline{X_{\epsilon}}$ denotes the smallest positive steady state of
(4.10). If $X(0)<\underline{X_{\epsilon}}$, then $X(t)$ converges to a remnant
state whereas if $X(0)>\underline{X_{\epsilon}}$, it converges towards a
significant steady state.

Figure 7. Panel (a): Population size $\langle y(t)\rangle$ (after space
integration) in the unprotected area $\Omega$ for Case I. The first vertical
dashed line indicates $t=\tau$, while the second one, indicates $t=T$. For
$t<\tau$, the dynamics is governed by (2.1b) with $\delta$ set to its
corresponding value given in Table 1 (Case I). The initial datum at $t=0$ is
taken to be $p_{0}$, the steady state of (2.1b) for $\delta=0$. The dynamics
over the time window $[\tau,T]$, is governed by (2.13b)–(2.13d) driven by
$u_{N}^{\ast}(t)$ whose space integration is shown in Panel (b). For $t>T$,
the dynamics is governed by (2.1b) but for $\delta=\delta^{\prime}$ (Table 1)
and with no control. Panel (b): Space integration, $\langle
u_{N}^{\ast}(t)\rangle$, of the optimal control $u_{N}^{\ast}$ obtained from
(3.29) (with $N=1$) and the BVP (LABEL:Eq_BVP_for_PMP3). Each curve is shown
here as $\alpha$ is varied according to the same values used for Panel (a) and
for a time window slightly larger than $[\tau,T]$.

Marking the dependence of $u_{N}^{\ast}$ on $\alpha$, one defines
$\underline{\alpha}$ to be the smallest value of $\alpha$ for which the
following condition holds true

$\langle y(T;u_{N}^{\ast}(\alpha)),e_{1}\rangle>\underline{X_{\epsilon}}.$
(4.11)

For Case I one finds $\underline{\alpha}=0.25$ using a mesh grid of size
$0.01$ to discretize the range $[0,1]$ of $\alpha$-values. Because of the
high-energy content carried by $e_{1}$, we expect that if (4.11) is satisfied,
then the solution, $y(t;u_{N}^{\ast}(\alpha))$, to the KPPH equation for
$\delta=\delta^{\prime}$, converges to a significant state as $t$ tends to
$\infty$. By showing the population size behavior (after space integration)
$\langle y(t)\rangle$ as time evolves, Panel (a) of Fig. 7 shows that this is
exactly what happens for $\alpha\geq\underline{\alpha}$ while the population
converges towards a remnant steady state for $\alpha<\underline{\alpha}$. One
might thus wish to get $\alpha$ as close as possible to $\underline{\alpha}$
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 of the rescue operation is very narrow, as $\alpha$ is
getting too close to $\underline{\alpha}$. So in practice a hard constraint
$C$ on the individuals of the reserve should not correspond to an
$\alpha$-value too close to $\underline{\alpha}$, as $C$ may e.g. not be
exactly respected during the displacement operation.

Panel (b) of Fig. 7 which shows the controls (applied over $[\tau,T]$)
corresponding to the solutions shown in Panel (a), illustrates this statement.
As shown in this panel, the control leading to a significant survival (red
curve) is indeed very close to that leading to extinction (blue curve). We are
here in presence of an interesting phenomenon, namely, continuous dependence
on the forcing may hold on finite-time intervals, but a high sensitivity in
the system’s response may take place in the asymptotic time.

We shall emphasize also that unlike what was reported for $\alpha=0.35$, one
may have $y(T;u_{N}^{\ast})>y_{0}$ while the population still evolves towards
a remnant steady state, when $\delta$ is set to $\delta^{\prime}$ after $t=T$
and the control is abandoned. Figure 7 shows for instance that the population
size in $\Omega$ on the rise from $y_{0}$ up to $t=T$ (due to control), drops
eventually to a remnant state in an asymptotic time for e.g. $\alpha=0.2$ and
$\alpha=0.24$.

Finally, we stress that the numerical results reported are not limited to the
particular numerical setup considered here. For instance if $D$ is further
reduced, then more modes become unstable, say $p$, and the dimension of an
efficient reduced system must be at least $p$. Here $N=p=1$ for $D=1$. By
setting $D=0.1$ in Case II while keeping the other parameters as in Sec. 4.1,
the two dominant eigenvalues of the spectral problem (2.3b) are negative
corresponding to two unstable modes $e_{1}$ and $e_{2}$ capturing most of the
energy of the target population density, unlike the energy decomposition shown
in Fig. 2. By choosing $N=2$, the BVP (LABEL:Eq_BVP_for_PMP3) becomes
(LABEL:Eq_BVP_for_PMP2) for $N=2$. The corresponding efficient Galerkin
approximation is then given by (3.3) for $N=2$ from which a critical $\alpha$
separating survival from extinction can also be determined from an analogue of
(4.11) in which $\underline{X_{\epsilon}}$ (resp. the projection onto $e_{1}$)
is replaced by the steady state of smaller norm (resp. the projection onto the
space spanned by $e_{1}$ and $e_{2}$). Such a remark about the inflation of
the efficient reduced dimension holds as $D$ approaches zero. For instance for
$D=0.01$ in Case II, there exists four unstable modes capturing most of the
energy, constraining thus the efficient reduced dimension to be $N=4$. The
approach presented here extends also to more realistic two-dimensional domains
in which the fragmentation effects of the habitat may play a more important
role than for the one-dimensional domains; see [RC07].

## Acknowledgments

This work has been partially supported by the European Research Council (ERC)
under the European Union’s Horizon 2020 research and innovation program (grant
agreement No. 810370 (MDC)), by the Office of Naval Research (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 results

In this appendix, we summarize from [CKL17] key convergence results and error
estimates obtained for Galerkin approximations of nonlinear optimal control
problems in Hilbert spaces. As we will see in Appendix B below, the material
of this section allows us to link precisely the optimal control problem (2.19)
to its Galerkin approximation (3.12), in terms of, both, error estimates about
the optimal controls, and (strong) convergence about 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 particular case of Galerkin systems from
eigenprojections as considered in this article. The interested reader is
referred to [CKL17] for more general cases, and also for convergence results
concerning the value functions. Applications to control in feedback form can
be found in [CKL18].

In the following, we consider finite-dimensional approximations of the
following initial-value problem (IVP):

$\displaystyle\frac{\,\mathrm{d}y}{\,\mathrm{d}t}$
$\displaystyle=Ly+F(y)+\mathfrak{C}(u(t)),\quad t\in(0,T],$ (A.1)
$\displaystyle y(0)$ $\displaystyle=y_{0},$

where the unknown $y$ evolves in a separable Hilbert space $\mathcal{H}$,
$L\colon D(L)\subset\mathcal{H}\rightarrow\mathcal{H}$ is a linear operator
with domain $D(L)$, $F\colon\mathcal{H}\rightarrow\mathcal{H}$ denotes the
nonlinearity, and the initial datum $y_{0}$ belongs to $\mathcal{H}$. The
time-dependent forcing $u$ lives in a possibly different separable Hilbert
space $V$; the (possibly nonlinear) mapping $\mathfrak{C}\colon
V\rightarrow\mathcal{H}$ is assumed to be such that $\mathfrak{C}(0)=0$. Other
assumptions regarding $\mathfrak{C}$ will be made below when needed.

We assume the following conditions for the linear operator $L$.

  * (H0)

The linear operator $L:D(L)\subset\mathcal{H}\rightarrow\mathcal{H}$ is the
infinitesimal generator of a $C_{0}$-semigroup of bounded linear operators
$T(t)$ on $\mathcal{H}$.

  * (H1)

$L$ is self-adjoint with compact resolvent.

Recall that under assumption (H0), the domain $D(L)$ of $L$ is dense in
$\mathcal{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\geq 1$, the set of admissible controls
$\mathcal{U}_{\textrm{ad}}$ is the subset of $L^{q}(0,T;V)$ constituted by
measurable functions that take values in $U$, a bounded subset of the Hilbert
space $V$.

In other words,

$\mathcal{U}_{\textrm{ad}}=\\{f\in L^{q}(0,T;V)\;:\;f(s)\in U\textrm{ for a.e.
}s\in[0,T]\\},\;\;q\geq 1.$ (A.2)

The set $\mathcal{U}_{\textrm{ad}}$ will be endowed with the induced topology
from that of $L^{q}(0,T;V)$.

Let $u$ be in $\mathcal{U}_{\textrm{ad}}$ given by (A.2), a mild solution to
(A.1) over $[0,T]$ is a function $y$ in $C([0,T],\mathcal{H})$ such that

$y(t)=T(t)y_{0}+\int_{0}^{t}T(t-s)F(y(s))\,\mathrm{d}s+\int_{0}^{t}T(t-s)\mathfrak{C}(u(s))\,\mathrm{d}s,\;\;t\in[0,T].$
(A.3)

In what follows we will often denote by $t\mapsto y(t;y_{0},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-adjoint compact operator that the eigenfunctions
of $L$ form an orthonormal basis of $\mathcal{H}$; see e.g. [Bré10]. We denote
the eigenpairs of $L$ by $\\{(\beta_{k},e_{k})\,:\,k\in\mathbb{N}\\}$. For
each $N\geq 1$, let $\mathcal{H}_{N}$ be the $N$-dimensional subspace of
$\mathcal{H}$ spanned by the first $N$ eigenfunctions of $L$:

$\mathcal{H}_{N}=\mathrm{span}\\{e_{k}\;:\;k=1,\ldots,N\\}.$ (A.4)

Denote also by $\Pi_{N}:\mathcal{H}\rightarrow\mathcal{H}_{N}$ the associated
orthogonal projector. Note that

$\mathcal{H}_{N}\subset D(L),\;\forall\,N\geq 1.$ (A.5)

The corresponding Galerkin approximation of (A.1) associated with
$\mathcal{H}_{N}$ is then given by:

$\displaystyle\frac{\,\mathrm{d}y_{N}}{\,\mathrm{d}t}$
$\displaystyle=L_{N}y_{N}+\Pi_{N}F(y_{N})+\Pi_{N}\mathfrak{C}(u(t)),\;t\in[0,T],$
(A.6) $\displaystyle y_{N}(0)$
$\displaystyle=\Pi_{N}y_{0},\;\;y_{0}\in\mathcal{H},$

where

$L_{N}=\Pi_{N}L\Pi_{N}:\mathcal{H}\rightarrow\mathcal{H}_{N}.$ (A.7)

In particular, the domain $D(L_{N})$ of $L_{N}$ is $\mathcal{H}$, because of
(A.5).

Throughout this article, a mapping
$f:\mathcal{W}_{1}\rightarrow\mathcal{W}_{2}$ between two Banach spaces,
$\mathcal{W}_{1}$ and $\mathcal{W}_{2}$, is said to be locally Lipschitz if
for any ball $\mathfrak{B}_{r}\subset\mathcal{W}_{1}$ with radius $r>0$
centered at the origin, there exists a constant ${\rm
Lip}(f|{{}_{\mathfrak{B}_{r}}})>0$ such that

$\|f(y_{1})-f(y_{2})\|_{\mathcal{W}_{2}}\leq{\rm
Lip}(f|{{}_{\mathfrak{B}_{r}}})\|y_{1}-y_{2}\|_{\mathcal{W}_{1}},\qquad\text{
}\forall\>y_{1},y_{2}\in\mathfrak{B}_{r}.$ (A.8)

We will make also use of the following assumptions.

  * (H3)

The mapping $F:\mathcal{H}\rightarrow\mathcal{H}$ is locally Lipschitz in the
sense given in (A.8).

  * (H4)

Let $\mathcal{U}_{\textrm{ad}}$ be given by (A.2). For each $T>0$ and
$(y_{0},u)$ in $\mathcal{H}\times\mathcal{U}_{\textrm{ad}}$, the problem (A.1)
admits a unique mild solution $y(\cdot;y_{0},u)$ in $C([0,T],\mathcal{H})$;
and for each $N\geq 1$, its Galerkin approximation (A.6) admits a unique
solution $y_{N}(\cdot;\Pi_{N}y_{0},u)$ in $C([0,T],\mathcal{H})$. Moreover,
there exists a constant $\mathcal{C}=\mathcal{C}(T,y_{0})$ such that

$\displaystyle\|y(t;y_{0},u)\|_{\mathcal{H}}\leq\mathcal{C},$
$\displaystyle\text{ }\forall\>t\in[0,T],\;u\in\mathcal{U}_{\textrm{ad}},$
(A.9a)
$\displaystyle\|y_{N}(t;\Pi_{N}y_{0},u)\|_{\mathcal{H}}\leq\mathcal{C},$
$\displaystyle\text{
}\forall\>t\in[0,T],\;N\in\mathbb{N},\;u\in\mathcal{U}_{\textrm{ad}}.$ (A.9b)

######  Remark A.1.

Note that in applications, (H4) is typically satisfied via a priori estimates.
Indeed, let $u$ be in $\mathcal{U}_{\textrm{ad}}$ given by (A.2). Then uniform
bounds such as in (A.9) are guaranteed if e.g. an a priori estimate of the
following type holds for the IVPs (A.1) and (A.6):

$\hskip
21.52771pt\underset{t\in[0,T]}{\sup}\|y(t;y_{0},u)\|_{\mathcal{H}}\leq\alpha(\|y_{0}\|_{\mathcal{H}}+\|u\|_{L^{q}(0,T;V)})+\beta,\;\;\;\;\alpha>0,\;\;\beta\geq
0.$ (A.10)

See e.g. [CH98, Tem97] for such a priori bounds for nonlinear partial
differential equations. Such bounds can also be derived for nonlinear systems
of delay differential equations (DDEs); see [CGLW16, CKL18]. For the KPPH
problem (2.13b)-(2.13c), the desired estimate (A.9a) is derived for
nonnegative initial data by using the maximum principle together with energy
estimates; see Appendix B. In particular as we will see, Theorem A.1 and
Theorem A.2 below, both apply to this context when $y_{0}\geq 0$ in (H4).

We introduce next the cost functional,
$J\colon\mathcal{H}\times\mathcal{U}_{\textrm{ad}}\rightarrow\mathbb{R}^{+}$,
associated with the IVP (A.1):

$J(y_{0},u)=\int_{0}^{T}[\mathcal{G}(y(s;y_{0},u))+\mathcal{E}(u(s))]\,\,\mathrm{d}s,\;y_{0}\in\mathcal{H},$
(A.11)

where $\mathcal{G}:\mathcal{H}\rightarrow\mathbb{R}^{+}$ and
$\mathcal{E}:V\rightarrow\mathbb{R}^{+}$ are assumed to be continuous, and
$\mathcal{G}$ is assumed to satisfy furthermore the condition:

$\mathcal{G}\mbox{ is locally Lipschitz in the sense of
}\eqref{Local_Lip_cond}.$ (A.12)

The associated optimal control problem then writes

$\displaystyle\min\,J(y_{0},u)\hskip 4.30554pt\text{ subject to }\hskip
4.30554pt(y,u)\in L^{2}(0,T;\mathcal{H})\times\mathcal{U}_{\textrm{ad}}\text{
solves}~{}\eqref{Abstract_PDE}\text{ with }\;y(0)=y_{0}\in\mathcal{H}.$
($\mathcal{P}$)

The cost functional,
$J_{N}\colon\mathcal{H}_{N}\times\mathcal{U}_{\textrm{ad}}\rightarrow\mathbb{R}^{+}$,
associated with the Galerkin approximation (A.6) is given by

$J_{N}(\Pi_{N}y_{0},u)=\int_{0}^{T}[\mathcal{G}(y_{N}(s;\Pi_{N}y_{0},u))+\mathcal{E}(u(s))]\,\,\mathrm{d}s,\;y_{0}\in\mathcal{H},$
(A.13)

and the corresponding optimal control problem reads:

$\displaystyle\min\,J_{N}(\Pi_{N}y_{0},u)\quad\text{ subject to
}\quad(y_{N},u)\in
L^{2}(0,T;\mathcal{H}_{N})\times\mathcal{U}_{\textrm{ad}}\text{
solves}~{}\eqref{ODE_Galerkin}$ ($\mathcal{P}_{N}$) $\displaystyle\hskip
150.00023pt\text{ with }\;y_{N}(0)=\Pi_{N}y_{0}\in\mathcal{H}_{N}.$

Then, as a consequence of [CKL17, Corollary 2.14] we can deduce the following
theorem about the error estimates 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_{0}$ in
$\mathcal{H}$, both ($\mathcal{P}$) and ($\mathcal{P}_{N}$) admit an optimal
control, denoted by $u^{*}$ and $u^{*}_{N}$, respectively. Assume furthermore
that there exists $\sigma>0$ such that the following local growth condition is
satisfied for the cost functional $J$ defined in (A.11):

$\sigma\|u^{*}-v\|_{L^{q}(0,T;V)}^{q}\leq J(y_{0},v)-J(y_{0},u^{*}),$ (A.14)

for all $v$ in some neighborhood $\mathcal{W}\subset\mathcal{U}_{\textrm{ad}}$
of $u^{*}$, with $\mathcal{U}_{\textrm{ad}}$ given by (A.2). Assume finally
that $u^{*}_{N}$ lies in $\mathcal{W}$. Then there exists $\gamma>0$ such that

$\displaystyle\|u^{\ast}-u^{\ast}_{N}\|_{L^{q}(0,T;V)}^{q}$
$\displaystyle\leq\frac{1}{\sigma}{\rm
Lip}(\mathcal{G}|_{\mathfrak{B}})\left[\sqrt{T}+\gamma
T\right]\Bigl{(}\|\Pi_{N}^{\perp}y(\cdot;y_{0},u^{*})\|_{L^{2}(0,T;\mathcal{H})}$
(A.15) $\displaystyle\hskip
100.00015pt+2\|\Pi_{N}^{\perp}y(\cdot;y_{0},u^{*}_{N})\|_{L^{2}(0,T;\mathcal{H})}\Bigr{)},$

where $\mathfrak{B}$ denotes the ball in $\mathcal{H}$ centered at the origin
with radius $\mathcal{C}$, with $\mathcal{C}$ being the same as given in
Assumption (H4), and $\Pi_{N}^{\perp}=\mathrm{Id}_{\mathcal{H}}-\Pi_{N}$.

As pointed out in [CKL17, Remark 2.13], it is not clear a priori that

$\lim_{N\rightarrow\infty}\|\Pi_{N}^{\perp}y(\cdot;y_{0},u^{\ast}_{N})\|_{L^{2}(0,T;\mathcal{H})}=0.$
(A.16)

The reason relies on the dependence on $u^{\ast}_{N}$ of
$\|\Pi_{N}^{\perp}y(\cdot;y_{0},u^{\ast}_{N})\|_{L^{2}(0,T;\mathcal{H})}$,
where $u^{\ast}_{N}$ denotes the control synthesized from the $N$-dimensional
Galerkin approximation. However, for the special case of Galerkin
approximations constructed from eigenbasis, the convergence in (A.16) is
guaranteed to hold under the assumptions (H0)–(H4) of Theorem A.1 above; see
[CKL17, Lemma 2.16]. As a result, when $N$ tends to infinity, $u^{\ast}_{N}$
converges to the optimal control $u^{\ast}$.

We have furthermore, the following uniform convergence result about the
controlled solutions as a result of [CKL17, Theorem 2.6].

######  Theorem A.2.

Assume that (H0)–(H4) hold and that the operator
$\mathfrak{C}:V\rightarrow\mathcal{H}$ is locally Lipschitz in the sense given
in (A.8) (with $\mathfrak{C}(0)=0$). Assume furthermore that the set $U$ in
Assumption (H2) is a compact subset of $V$, with $q>1$ therein.

Then, for any $(y_{0},u)$ in $\mathcal{H}\times\mathcal{U}_{\textrm{ad}}$, the
solution $y_{N}$ of the Galerkin approximation (A.6) converges uniformly to
the mild solution $y$ of (A.1) in the sense that:

$\lim_{N\rightarrow\infty}\sup_{u\in\mathcal{U}_{\textrm{ad}}}\sup_{t\in[0,T]}\|y_{N}(t;\Pi_{N}y_{0},u)-y(t;y_{0},u)\|_{\mathcal{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>1$.

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 $\mathcal{H}$. Then,
the operator $L_{N}$ defined by (A.7) as the eigen projection of $L$ onto
$\mathcal{H}_{N}$ given by (A.4) clearly satisfies

$\lim_{N\rightarrow\infty}\|L_{N}\phi-L\phi\|_{\mathcal{H}}=0,\quad\forall\,\phi\in
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}:\mathcal{H}_{N}\rightarrow\mathcal{H}_{N}$ generated by $L_{N}$
satisfies

$\|e^{L_{N}t}\|\leq e^{\beta_{1}t},\quad N\in\mathbb{N},\,t\geq 0,$ (A.19)

where $\beta_{1}$ is the largest eigenvalue of $L$. Defining the extension
$T_{N}(t):\mathcal{H}\rightarrow\mathcal{H}$ of $e^{L_{N}t}$ to be

$T_{N}(t)\phi=e^{L_{N}t}\Pi_{N}\phi+(\mathrm{Id}_{\mathcal{H}}-\Pi_{N})\phi,\quad\phi\in\mathcal{H},$
(A.20)

Assumption (A1) in [CKL17] follows then by choosing the parameters $M$ and
$\omega$ therein to be $M=1$ and $\omega=\max\\{\beta_{1},0\\}$.

##  Appendix B Convergence and error estimates results for the optimal control
of the KPPH equation

We check here, in the context of the optimal control of the KPPH equation, the
assumptions of Theorem A.2 and Theorem A.1 about the convergence and error
estimates results, respectively, allowing us in particular to deduce Theorem
3.1 of Sec. 3.1.

Within this context, (2.19) and (3.12) play the role of ($\mathcal{P}$) and
($\mathcal{P}_{N}$), respectively. The optimal problems (2.19) and (3.12) are
posed on the time interval $[\tau,T]$ but a simple change of variable
$t^{\prime}=t-\tau$ allows us to frame these as in Appendix A, that is over
the time interval $[0,T]$. We operate this shift below to ease the
presentation.

The verification of the required assumptions is organized in several steps.

Step 1: Verifications of (H0), (H1), and (H3). Let
$\mathcal{H}=L^{2}(\Omega)$. We first put the IVP (2.13b)-(2.13d) into the
form (A.1). The corresponding operators $L\colon D(L)\rightarrow\mathcal{H}$,
$F\colon\mathcal{H}\rightarrow\mathcal{H}$ and
$\mathfrak{C}\colon\mathcal{H}\rightarrow\mathcal{H}$ are then naturally
defined as follows:

$\displaystyle Ly=D\nabla^{2}y+\mu(\cdot)y,$ $\displaystyle y\in
D(L)=H^{2}(\Omega)\cap\Big{\\{}y\in H^{1}(\Omega)\;\big{|}\;\frac{\partial
y}{\partial\boldsymbol{n}}=0\Big{\\}},$ (B.1) $\displaystyle
F(y)=-\nu(\cdot)y^{2}-\delta\rho_{\epsilon}(y),$ $\displaystyle
y\in\mathcal{H},$ $\displaystyle\mathfrak{C}(u)=u,$ $\displaystyle
u\in\mathcal{H}.$

It is standard that Assumptions (H0) and (H1) are satisfied for the elliptic
operator $L$ defined in (B.1); see e.g. [Paz83]. It can be checked that $F$ is
locally Lipschitz as mapping from $\mathcal{H}$ to $\mathcal{H}$, in the sense
of (A.8). Thus Assumption (H3) is satisfied. It is also clear that
$\mathfrak{C}$ defined above is Lipschitz on $\mathcal{H}$ and satisfies
furthermore that $\mathfrak{C}(0)=0$.

Step 2: Verification of (H2). Recall that the admissible set
$\mathcal{U}_{\textrm{ad}}$ is defined by (2.18). Since
$\varphi_{j}=\frac{1}{\sqrt{|\Lambda_{j}|}}\chi_{\Lambda_{j}}$ (cf. (2.17)),
each admissible control $u$ in $\mathcal{U}_{\textrm{ad}}$ takes value in the
following subset $U$ of $\mathcal{H}$ at a given time instant:

$U=\Big{\\{}\psi\in\mathcal{H}\;\big{|}\;\psi(x)=\sum_{j=1}^{K}\frac{b_{j}}{\sqrt{|\Lambda_{j}|}}\chi_{\Lambda_{j}}(x),\;x\in\Omega,\;0\leq
b_{j}\leq C_{j}\Big{\\}}.$ (B.2)

Thus $\mathcal{U}_{\textrm{ad}}$ in (2.18) can be rewritten as

$\mathcal{U}_{\textrm{ad}}=L^{2}(0,T;U).$ (B.3)

Taking the space $V$ in (H2) to be $L^{2}(\Omega)$, it is clear that $U$ is a
bounded set in $V$. Assumption (H2) is thus verified with $q=2$ because of
(B.3). We show below that $U$ is furthermore a compact set in $V$, which is
required in Theorem A.2. Every sequence $(\psi_{n})$ in $U$ has a convergent
subsequence in $U$. Indeed for each $1\leq j\leq K$, the corresponding
$b_{j}^{n}$ has a convergent subsequence as taking value in the bounded and
closed subset $[0,C_{j}]$ of $\mathbb{R}$. Because we have a finite number of
such convergent subsequences, one can always find a common extraction and thus
$(\psi_{n})$ for which convergence holds towards an element of $U$. Thus $U$
is a compact subset of $V=L^{2}(\Omega)$.

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\geq 0$ when $y_{0}\geq 0$.
This property is a consequence of the maximum principle which holds for
(2.13b)-(2.13d), a known fact but of which we provide a proof of using the
Stampacchia truncation method [MS68]. The arguments are standard but are
sketched below as they provide also useful insights to prove the required a
priori bounds for the solution to the Galerkin approximation of
(2.13b)-(2.13d). For this purpose, one defines the negative part $v^{-}$ of a
measurable function $v:\Omega\rightarrow\mathbb{R}$ by

$v^{-}(x)=\min\\{v(x),0\\}.$ (B.4)

A lemma due to Stampacchia ensures that if $v$ is in $H^{1}(\Omega)$, then
$v^{-}$ lies also in $H^{1}(\Omega)$ and

$\nabla v^{-}=\chi_{\\{v<0\\}}\nabla v.$ (B.5)

The so-called truncation method of Stampacchia allows us to obtain a powerful
identity for a broad class of inhomgeneous heat problem

$\displaystyle\;\partial_{t}y=D\nabla^{2}y+f,\quad\textrm{ in
}(0,T)\times\Omega,$ (B.6a) $\displaystyle\;\frac{\partial
y}{\partial\boldsymbol{n}}=0,\quad\textrm{ on }(0,T)\times\partial\Omega,$
(B.6b) $\displaystyle\;y(0,x)=y_{0}(x),\;x\in\Omega.$ (B.6c)

This identity is derived from the variational formulation of this problem, by
using (B.5) with $v=y^{-}$. It gives

$\frac{1}{2}\big{|}y^{-}(t)\big{|}^{2}_{L^{2}(\Omega)}+D\int_{0}^{t}\big{|}\nabla
y^{-}(s)\big{|}^{2}_{L^{2}(\Omega)}\,\mathrm{d}s=\int_{0}^{t}\int_{\Omega}fy^{-}\,\mathrm{d}x\,\mathrm{d}s+\frac{1}{2}\big{|}y^{-}(0)\big{|}^{2}_{L^{2}(\Omega)}.$
(B.7)

From this identity the classical maximum principle can be deduce in the sense
that if $f\geq 0$ and $y_{0}\geq 0$ a.e. then $y\geq 0$ a.e. Let us introduce

$\mathcal{N}(x,y,u)=\mu(x)y-\nu(x)y^{2}-\delta\rho_{\epsilon}(y)+u(t,x),$
(B.8)

This result cannot be applied directly with $f=\mathcal{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=\widetilde{\mathcal{N}}$, with

$\widetilde{\mathcal{N}}(x,y,u)=\mu(x)y-\nu(x)y|y|-\delta\rho_{\epsilon}(y)+u(t,x).$
(B.9)

Indeed, first note that for $f=\widetilde{\mathcal{N}}$,

$\int_{0}^{t}\int_{\Omega}fy^{-}\,\mathrm{d}x\,\mathrm{d}s\leq\|\mu\|_{L^{\infty}(\Omega)}\int_{0}^{t}\big{|}y^{-}(s)\big{|}^{2}_{L^{2}(\Omega)}\,\mathrm{d}s-\int_{0}^{t}\int_{\Omega}\nu(x)y|y|y^{-}\,\mathrm{d}x\,\mathrm{d}s$
(B.10)

since $\rho_{\epsilon}(y)y^{-}=0$, because $\rho_{\epsilon}(s)=0$ for all
$s\leq 0$, and $uy^{-}\leq 0$, because $u$ lies in $\mathcal{U}_{\textrm{ad}}$
and thus is nonnegative. Then, since

$\nu y|y|y^{-}=\nu|y|(y^{-})^{2}\geq 0,$ (B.11)

we deduce from (B.10) that

$\int_{0}^{t}\int_{\Omega}\widetilde{\mathcal{N}}y^{-}\,\mathrm{d}x\,\mathrm{d}s\leq\|\mu\|_{L^{\infty}(\Omega)}\int_{0}^{t}\big{|}y^{-}(s)\big{|}^{2}_{L^{2}(\Omega)}\,\mathrm{d}s.$
(B.12)

Assuming now that $y_{0}\geq 0$ a.e., we obtain from (B.7) that

$\big{|}y^{-}(t)\big{|}^{2}_{L^{2}(\Omega)}\leq
2\|\mu\|_{L^{\infty}(\Omega)}\int_{0}^{t}\big{|}y^{-}(s)\big{|}^{2}_{L^{2}(\Omega)}\,\mathrm{d}s.$
(B.13)

Because $y^{-}(0)=(y_{0})^{-}=0$, the Gronwall’s lemma in its integral form
allows us to conclude that $y^{-}=0$ a.e. in $(0,T)\times\Omega$. Thus we have
proved that $y\geq 0$ a.e. if $y$ solves (B.6c) with $y_{0}\geq 0$ and
$f=\widetilde{\mathcal{N}}$ a.e.

From this we deduce that $|y|=y$ a.e. $x$, $t$. Thus equation (B.6a) with
$f=\widetilde{\mathcal{N}}$ is the same as equation (2.13b) when $y_{0}\geq 0$
and in the end we have also found that $y$ is a positive solution to
(2.13b)-(2.13d) when $y_{0}\geq 0$.

On the other hand, the inner product with $y$ on both sides of (2.13b) leads
to

$\displaystyle\frac{1}{2}\frac{\,\mathrm{d}}{\,\mathrm{d}t}\big{|}y\big{|}^{2}_{L^{2}(\Omega)}$
$\displaystyle=\langle Ly+F(y)+\mathfrak{C}(u(t)),y\rangle_{\mathcal{H}}$
(B.14) $\displaystyle=-D\big{|}\nabla
y\big{|}^{2}_{L^{2}(\Omega)}+\int_{\Omega}\big{(}\mu(x)y^{2}-\nu(x)y^{3}\big{)}\,\mathrm{d}x-\delta\int_{\Omega}\rho_{\epsilon}(y)y\,\mathrm{d}x+\int_{\Omega}uy\,\mathrm{d}x$

Thus because $y$ is a positive solution and $\rho_{\epsilon}\geq 0$, we obtain

$\frac{1}{2}\frac{\,\mathrm{d}}{\,\mathrm{d}t}\big{|}y\big{|}^{2}_{L^{2}(\Omega)}\leq\|\mu\|_{\infty}\big{|}y\big{|}^{2}_{L^{2}(\Omega)}+C\int_{\Omega}y\,\mathrm{d}x,$
(B.15)

with
$C=\max\Big{\\{}\frac{C_{1}}{\sqrt{|\Lambda_{1}|}},\ldots,\frac{C_{K}}{\sqrt{|\Lambda_{K}|}}\Big{\\}}$.

By remarking that by Hölder’s inequality, we have

$\int_{\Omega}y\,\mathrm{d}x\leq\sqrt{|\Omega|}\big{|}y\big{|}_{L^{2}(\Omega)}\leq\frac{1}{2}(|\Omega|+\big{|}y\big{|}_{L^{2}(\Omega)}^{2}).$

Using this in (B.15), we obtain then

$\frac{1}{2}\frac{\,\mathrm{d}}{\,\mathrm{d}t}\big{|}y\big{|}^{2}_{L^{2}(\Omega)}\leq
c_{1}\big{|}y\big{|}^{2}_{L^{2}(\Omega)}+c_{2},$ (B.16)

where

$c_{1}=\|\mu\|_{\infty}+\frac{C}{2},\quad c_{2}=\frac{C}{2}|\Omega|.$ (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 $\tilde{y}_{N}$ of a modified version of
(2.13b)-(2.13c) in which the reaction term in (2.13b) is replaced by
$\widetilde{\mathcal{N}}$. We assume the initial datum $\tilde{y}_{0}$ to be
nonnegative. The corresponding $N$-dimensional Galerkin approximation
$\tilde{y}_{N}$ satisfies thus

$\partial_{t}\tilde{y}_{N}=D\nabla^{2}\tilde{y}_{N}+\Pi_{N}\widetilde{\mathcal{N}}(x,\tilde{y}_{N},u).$
(B.18)

From this equation one gets,

$\displaystyle\frac{1}{2}\frac{\,\mathrm{d}}{\,\mathrm{d}t}\big{|}\tilde{y}_{N}\big{|}^{2}_{L^{2}(\Omega)}$
$\displaystyle=-D\big{|}\nabla\tilde{y}_{N}\big{|}^{2}_{L^{2}(\Omega)}+\langle\Pi_{N}\widetilde{\mathcal{N}}(x,\tilde{y}_{N},u),\tilde{y}_{N}\rangle.$
(B.19)

Since $\tilde{y}_{N}$ lies in $\mathcal{H}_{N}$, we have

$\displaystyle\langle\Pi_{N}\widetilde{\mathcal{N}}(x,\tilde{y}_{N},u),y_{N}\rangle$
$\displaystyle=\langle\widetilde{\mathcal{N}}(x,\tilde{y}_{N},u),\tilde{y}_{N}\rangle$
(B.20)
$\displaystyle=\int_{\Omega}\big{(}\mu(x)\tilde{y}_{N}^{2}-\nu(x)|\tilde{y}_{N}|\tilde{y}_{N}^{2}\big{)}\,\mathrm{d}x-\delta\int_{\Omega}\rho_{\epsilon}(y_{N})y_{N}\,\mathrm{d}x+\int_{\Omega}u\tilde{y}_{N}\,\mathrm{d}x$
$\displaystyle\leq\|\mu\|_{\infty}\int_{\Omega}\tilde{y}_{N}^{2}\,\mathrm{d}x+(\delta+C)\int_{\Omega}|\tilde{y}_{N}|\,\mathrm{d}x,$

where $C$ is the same as in (B.15) and we have used the fact that
$|\rho_{\epsilon}(y_{N})|$ is bounded above by $1$; cf. (2.2).

We infer thus that

$\displaystyle\frac{1}{2}\frac{\,\mathrm{d}}{\,\mathrm{d}t}\big{|}\tilde{y}_{N}\big{|}^{2}_{L^{2}(\Omega)}\leq\|\mu\|_{\infty}\big{|}\tilde{y}_{N}\big{|}^{2}_{L^{2}(\Omega)}+(\delta+C)\int_{\Omega}|\tilde{y}_{N}|\,\mathrm{d}x,$
(B.21)

leading as for $y$ satisfying (B.15) to a uniform bound for $\tilde{y}_{N}$,
and consequently (A.9b) holds for $\tilde{y}_{N}$. Due to Theorem A.2,
$\tilde{y}_{N}$ satisfies the convergence property (A.17) towards the solution
$\tilde{y}$ to the modified IVP (2.13b)-(2.13d) in which
$\widetilde{\mathcal{N}}$ replaces the reaction term. Now we know that
$\tilde{y}\geq 0$ since $\tilde{y}_{0}\geq 0$, and since
$|\tilde{y}|=\tilde{y}$ in this case, we get that $\tilde{y}$ solves actually
the original IVP (2.13b)-(2.13d) which we denote now by $y$ and
$\tilde{y}_{0}$ by $y_{0}$.

Smoothing arguments ensure that the mild solution $y$ in $L^{2}(\Omega)$ to
(2.13b)-(2.13d) is actually continuous on $\Omega$. Thus because we have, for
$\mathcal{U}_{\textrm{ad}}$ defined in (B.3)

$\sup_{u\in\mathcal{U}_{\textrm{ad}}}\sup_{t\in[0,T]}|\tilde{y}_{N}(t;\Pi_{N}y_{0},u)-y(t;y_{0},u)|_{L^{2}(\Omega)}\underset{N\rightarrow\infty}{\longrightarrow}0,$
(B.22)

we have $\tilde{y}_{N}>0$ for $N$ sufficiently large (when $y_{0}>0$) and thus
(B.18) reduces to the Galerkin approximation of 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 error estimates
(A.15) of Theorem A.1 also hold with $q=2$ and $V=L^{2}(\Omega)$ and Theorem
3.1 of Sec. 3.1 is thus proved.

## References

  * [Ama76] H. Amann, _Fixed point equations and nonlinear eigenvalue problems in ordered banach spaces_ , SIAM review 18 (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 Applied Mathematics 55 (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 51 (2005), no. 1, 75–113.
  * [BHR05b] by same author, _Analysis of the periodically fragmented environment model: II—Biological invasions and pulsating travelling fronts_ , Journal de Mathématiques purés et appliquées 84 (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 for Conservation of Nature (IUCN), 2004.
  * [BM77] J.R. Beddington and R.M. May, _Harvesting natural populations in a randomly fluctuating environment_ , Science 197 (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 of environmental resources_ , Annu. Rev. Resour. Econ. 6 (2014), 33–68.
  * [CC89] R. S. Cantrell and C. Cosner, _Diffusive logistic equations with indefinite weights: population models in disrupted environments_ , Proc. R. Soc. Edinb. Sect. A 112 (1989), no. 3-4, 293–318.
  * [CC04] by same author, _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 differential equations_ , Disc. Cont. Dyn. Sys. A 36 (2016), no. 8, 4133–4177.
  * [CH98] T. Cazenave and A. Haraux, _An Introduction to Semilinear Evolution Equations_ , Oxford Lecture Series in Mathematics and 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 semilinear parabolic equations_ , ESAIM: Control, Optimisation and Calculus of Variations 23 (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 (2017), no. 189, 1–40.
  * [CKL18] M. D. Chekroun, A. Kröner, and H. Liu, _Galerkin approximations for the optimal control of nonlinear delay differential equations_ , 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: mathematical results_ , Comptes Rendus Mathématique 343 (2006), no. 5, 307–310.
  * [Ded10] L. Dedè, _Reduced basis method and a posteriori error estimation for parametrized linear-quadratic optimal control problems_ , SIAM Journal on Scientific Computing 32 (2010), 997–1019.
  * [Fah03] L. Fahrig, _Effects of habitat fragmentation on biodiversity_ , Annual review of ecology, evolution, and systematics 34 (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 the time-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 7 (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 1 (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 optimal control problems_ , C. R. Acad. Sci. Paris, Ser. I 349 (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] by same author, _Three control methods for time-dependent fluid flow_ , Flow, Turbulence and Combustion 65 (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 state constraints_ , SIAM review 37 (1995), no. 2, 181–218.
  * [HV05] M. Hinze and S. Volkwein, _Proper orthogonal decomposition surrogate models for nonlinear dynamical systems: error estimates 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 dynamical systems_ , SIAM J. Numer. Anal. 46 (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 1 (1937), 1–26.
  * [KS08] K. Kurata and J. Shi, _Optimal spatial harvesting strategy and symmetry-breaking_ , Appl. Math. Optim. 58 (2008), no. 1, 89–110.
  * [MS68] M.R.V Murthy and G. Stampacchia, _Boundary value problems for some degenerate-elliptic operators_ , Annali di Matematica Pura ed Applicata 80 (1968), no. 1, 1–122.
  * [Mur01] J.D. Murray, _Mathematical Biology II: Spatial Models and Biomedical Applications_ , Springer New York, 2001.
  * [Mur07] by same author, _Mathematical Biology: I. An Introduction_ , vol. 17, Springer Science & Business Media, 2007.
  * [Neu03] M.G. Neubert, _Marine reserves and optimal harvesting_ , Ecology Letters 6 (2003), no. 9, 843–849.
  * [OL13] A. Okubo and S. A. Levin, _Diffusion and ecological problems: Modern perspectives_ , vol. 14, Springer Science & Business Media, 2013.
  * [OS02] S. Oruganti and J. Shi, _Diffusive logistic equation with constant yield harvesting. i: Steady states_ , Trans. Am. Math. Soc. 354 (2002), no. 9, 3601–3619.
  * [Paz83] A. Pazy, _Semigroups of Linear Operators and Applications to Partial Differential Equations_ , Applied Mathematical Sciences, 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 Optimal Processes_ , 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_ , International Journal for Numerical Methods in Fluids 34 (2000), no. 5, 425–448.
  * [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. 68 (2007), 133–153.
  * [RC10] L. Roques and M.D. Chekroun, _Does reaction-diffusion support the duality of fragmentation effect?_ , Ecological Complexity 7 (2010), no. 1, 100–106.
  * [Sch54] M. B. Schaefer, _Some aspects of the dynamics of populations important to the management of the commercial marine fisheries_ , Inter-American Tropical Tuna Commission Bulletin 1 (1954), no. 2, 23–56.
  * [Sch91] by same author, _Some aspects of the dynamics of populations important to the management of the commercial marine fisheries_ , Bulletin of Mathematical Biology 53 (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 and comparison of models_ , Journal of Applied Ecology 39 (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_ , Theoretical Population Biology 30 (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 in Mathematics, 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. 44 (2009), 83–115.