\documentclass[reqno]{amsart}
\usepackage{hyperref}

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2016 (2016), No. 214, pp. 1--13.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2016 Texas State University.}
\vspace{9mm}}

\begin{document}
\title[\hfilneg EJDE-2016/214\hfil 
 Optimal harvesting in diffusive population models]
{Optimal harvesting in diffusive population models with size 
random growth and \\ distributed recruitment}

\author[Q. Xie, Z.-R. He, X. Wang \hfil EJDE-2016/214\hfilneg]
{Qiangjun Xie, Ze-Rong He, Xiaohui Wang} 

\address{Qiangjun Xie (corresponding author) \newline
Institute of Operational Research and Cybernetics, 
Hangzhou Dianzi University, Zhejiang 310018, China. \newline
Department of Mathematics, 
University of Texas-Rio Grande Valley, 
Edinburg, TX 78539, USA}
\email{qjunxie@hdu.edu.cn}

\address{Ze-Rong He \newline
Institute of Operational Research and Cybernetics, 
Hangzhou Dianzi University, Zhejiang 310018, China}
\email{zrhe@hdu.edu.cn}

\address{Xiaohui Wang \newline
Department of Mathematics, 
University of Texas-Rio Grande Valley, 
Edinburg, TX 78539, USA}
\email{xiaohui.wang@utrgv.edu}

\thanks{Submitted September 10, 2015. Published August 11, 2016.}
\subjclass[2010]{92B05, 93C20, 49K20}
\keywords{Optimal harvesting; spatial diffusion; size-structured model;
\hfill\break\indent random growth; normal cone} 

\begin{abstract}
 In this article, we consider an optimal harvesting control problem for 
 a spatial diffusion population system, which incorporates
 individual's random growth of size and distributed style of recruitment. 
 The existence and uniqueness of nonnegative solutions to
 this practical model is established by means of Banach's fixed point theorem. 
 The continuous dependence of population density on the
 harvesting effort is analyzed. The optimal harvesting strategies are 
 discussed through normal cone and adjoint techniques.
 Some conditions are presented to assure that there is only one optimal policy.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section] 
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{remark}[theorem]{Remark}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{proposition}[theorem]{Proposition}
\allowdisplaybreaks

\section{Introduction} 

Structured population models provide the connection between the population 
level dynamics and individual level vital rates.
It has attracted a lot of attention from a rather diverse group of scientific 
researchers in biology and mathematics \cite{Cushing1998, Magal2008}.
Dynamic analyses  on the size-structured and age-structured population models
 are presented in \cite{ Anitabook2000,Ebenman1998}.
Optimal control and optimization analyses have also been considered extensively 
from the economical and ecological points of view
\cite{Anitabook2000, Anitabook2010, Apreutesei2014, Ding2010, He2012, 
Troltzsch2010}. 
To the optimal harvesting problems, there are quite many meaningful
results on the age-structured population systems with or without spatial 
diffusion \cite{Anitabook2000, Anitabook2010, Liu2013, Luo2007, Zhao2005} 
and the references therein.

For more realistic biological significance of modeling,  Hadeler
 \cite{Hadeler2010} proposed structured population models with diffusion 
in the size-space. The biological motivation is that the diffusion allows 
for ``stochastic noise" to be incorporated in the models, namely, the 
stochastic fluctuations around the tendency to growth. Faugeras and 
Maury \cite{Faugeras2005} established an advection-diffusion-reaction
 model of fish with length (i.e. size structure) and plane position 
(i.e. spatial structure) distributions. The diffusion-convection process 
with respect to size is also called the random growth process \cite{Chu2009}.
 For these models, the existence and asymptotic behaviors of solutions 
were shown by semigroup theories in \cite{Farkas2011} and Hopf bifurcation 
properties  with the modified Ricker type birth function were studied 
in \cite{Chu2009}.  Recently, some numerical approximate solutions by 
the method of lines were investigated in \cite{Bart.2015}. Up to present, 
it seems that very few results on the optimal harvesting control problems 
is presented for these biological models with the size random growth.

Inspired by the above results, we are concerned with the optimal harvesting 
for the following diffusive population model with the size random growth 
(Let $Q:=\Omega\times(S_0,S_1)\times(0,T)$, and 
$ \Sigma:=\partial\Omega\times(S_0,S_1)\times(0,T)$):
\begin{gather}
\begin{aligned}
\partial_tp-k\Delta p & = \partial_s(d(s)\partial_s p-g(s)p)
-\mu p-up \\
&\quad +\int_{S_0}^{S_1}\beta(x, s, t, \hat{s})p(x, \hat{s}, t)
\,{\rm d}\hat{s},\quad \text{in }Q ,
\end{aligned}  \label{e1} \\
 \frac{\partial p}{\partial s}(x, S_0, t)
=\frac{\partial p}{\partial s}(x, S_1, t)=0,\quad 
\forall (x, t)\in\Omega\times(0,T),\label{e2} \\
\frac{\partial p}{\partial n}(x, s, t)=0,\quad \text{on }\Sigma,
\label{e3} \\
  p(x, s,0)=p^0(x, s),\quad  \forall (x,s)\in\Omega\times(S_0,S_1),  \label{e4}
\end{gather}   
where $\Delta$ stands for the Laplace operator with respect to the  
spatial variable $x$, and $\Omega\subset R^N(N\leq 3)$ is a bounded 
open domain with a boundary $\partial\Omega$ smooth enough, $k>0$ 
is the spatial diffusion coefficient, $\mu:=\mu(x, s, t)$ denotes 
the death rate, and $u:=u(x, s, t)$ is the harvesting effect which can 
be controlled by the outside force. The constants $S_0$ and $S_1$ stand 
for the minimal and maximal sizes of individuals, respectively. $T>0$ is 
the finite horizon of control, independent of any initial-boundary conditions. 
$p(x, s, t)$ denotes the population density of individuals of size 
$s\in[S_0,S_1]$ at time $t\in[0,T]$ at location $x\in\Omega$. 
The homogeneous Neumann boundary conditions are introduced with respect to 
the $N$-dimension spatial variable $x$  and 1-dimension size variable $s$.

The individuals' size random growth process is described here by the term 
$\partial_s(d(s)\partial_s p-g(s)p)$ in the Eq. (1). Here, $d(s)>0$ on
 $[S_0,S_1]$ stands for the size-specific diffusion coefficient, and $g(s)$ 
is the growth modulus. Similar to \cite{Farkas2011}, we choose the non-local 
integral term in (1) as the recruitment process. The distributed recruitment 
means that individuals may be recruited into the population at different 
sizes with the rate $\beta(x, s, t, \hat{s})$. This choice is different 
from the one given in \cite{Faugeras2005}.

The aim of this article is to study the  optimal harvesting control problem
maximize
\begin{equation} 
J(u):=\int_0^T\int_{S_0}^{S_1}\int_{\Omega}[wp^uu-\frac{1}{2}\rho u^2]
\,\mathrm{d}x\,\mathrm{d}s\,\mathrm{d}t,\label{e5}
\end{equation}
subject to
$$
u\in\mathcal{U}=\{v\in L^2(Q):0\leq \zeta_1(x, s, t)\leq v(x, s, t)
\leq \zeta_2(x, s, t)\text{ a.e. in }\ Q\},
$$
where $w:=w(x, s, t)$ denotes the economic value of an individual of size 
$s\in[S_0,S_1]$ at time $t\in[0,T]$ at $x\in\Omega$. 
$\rho>0$ is a cost factor for implementing the harvesting policy $u$, 
and $p^u(x, s, t)$ is the solution of the system \eqref{e1}--\eqref{e4} 
corresponding to $u$.

We assume that the following conditions hold throughout this article:
\begin{itemize}
\item[(H1)] $g\in C^1[S_0,S_1]$, $g(S_0)>0$;

\item[(H2)] $\mu\in L_{loc}^{\infty}(\overline{\Omega}
\times[S_0,S_1)\times[0,T])$, $\mu(x, s, t)\geq 0$,  a.e. in $Q$;

\item[(H3)] $\beta(x, s, t, \hat{s})\geq 0$ a.e. in 
$\Omega\times(S_0,S_1)\times (0,T)\times(S_0,S_1)$, $
\beta\in L^{\infty}$ and let $\overline{\beta}:=\|\beta\|_{\infty}$;

\item[(H4)] $d(s)\geq d_1>0$  a.e. in $(S_0,S_1)$, $d\in L^{\infty}(S_0,S_1)$ 
and let $\overline{d}:=\|d\|_{\infty}$;

\item[(H5)] $p^0(x, s)\geq 0$  a.e. in $\Omega\times(S_0,S_1)$, 
$p_0\in L^2(\Omega\times(S_0,S_1))$;

\item[(H6)] $w(s,t,x)>0$  a.e. in $Q$ from \eqref{e5}, $w\in L^{\infty}(Q)$ and 
let $W:=\|w\|_{\infty}$.

\end{itemize}

The rest of this article is organized as follows. 
In Section 2, we deal with the existence and uniqueness of solutions 
of the state system \eqref{e1}--\eqref{e4} with the given parameters. 
Then we display the optimal strategies by feedback laws in Section 3, 
and establish the existence of optimal harvesting control and a unique 
optimal policy in Section 4. A short conclusion is given in Section 5.

\section{Existence and uniqueness of solution of the state system}

In this section, we  establish the existence and uniqueness of a positive 
weak solution to the state system \eqref{e1}--\eqref{e4}.

Let $\mathcal{Q}=\Omega \times (S_0, S_1)$ be an open subset of $R^{N+1}$. 
Then $Q=\mathcal{Q}\times (0, T)$. We regard $p(x, s, \cdot)$ as an element 
of the functional space $H:=L^2(\mathcal Q)$. For any $t\in [0,T]$
we have
\begin{equation}
\int_{S_0}^{S_1}\int_{\Omega}|p(x,s,t)|^2\,\mathrm{d}x\,\mathrm{d}s<\infty.
\label{e6}
\end{equation}
Denote by $H^1(\mathcal Q)$ the Sobolev space $W^{1,2}(\mathcal Q)$ 
endowed with the norm
\begin{equation}
\|p\|_{H^1(\mathcal Q)}=\Big(\int_{S_0}^{S_1}\int_{\Omega}
(p^2+|\nabla_x p|^2+|\partial_s p|^2)\,\mathrm{d}x\,\mathrm{d}s \Big)^{1/2}.
\label{e7}
\end{equation}
Let $H^1(\mathcal Q)^*$ denote the dual of $H^1(\mathcal Q)$.
 Then we have the chain of dense and continuous embeddings
\begin{equation}
H^1(\mathcal Q)\hookrightarrow H\hookrightarrow H^1(\mathcal Q)^*,\label{e8}
\end{equation}
and any $F\in H^1(\mathcal Q)^*$ can be continuously extended to $H$ if 
and only if there is some $f\in H$ such that
\begin{equation}
F(p)=\int_{S_0}^{S_1}\int_{\Omega}f\cdot p \,\mathrm{d}x\,\mathrm{d}s
=(f,p)_H, \quad \forall p\in H^1(\mathcal Q).\label{e9}
\end{equation}

\begin{definition} \label{def1} \rm  %W(0,T)}
 We denote by $W(0,T)$ the linear space of all $p\in L^2(0,T; H^1(\mathcal Q))$ 
which has a distributional derivative $p'\in L^2(0,T; H^1(\mathcal Q)^*)$, 
equipped with the norm
\begin{equation}
\|p\|_{W(0,T)}=\Big(\int_0^T\left(\|p\|_{H^1(\mathcal Q)}^2
+\|p'(t)\|_{H^1(\mathcal Q)^*}^2\Big)\,\mathrm{d}t\right)^{1/2}.\label{e10}
\end{equation}
\end{definition}

The space $W(0,T)=\{p\in L^2(0,T; H^1(\mathcal Q)):
 \frac{dp}{dt}\in L^2(0,T; H^1(\mathcal Q)^*)\}$ is a Hilbert space with 
the inner product
\begin{equation}
(p,q)_{W(0,T)}=\int_0^T(p,q)_{H^1(\mathcal Q)}\,\mathrm{d}t+\int_0^T(p'(t),
q'(t))_{H^1(\mathcal Q)^*}\,\mathrm{d}t.\label{e11}
\end{equation}
From \cite{Troltzsch2010}, we have
\begin{equation}
W(0,T)\hookrightarrow C([0,T]; H).\label{e12}
\end{equation} 

For the sake of convenience, we  change  the unknown function $p$ in 
the equation \eqref{e1}) by $\hat{p}=e^{-\theta t}p$ 
($\theta$ is to be determined latter). 
Then we have the following proposition.

\begin{proposition} \label{prop2} % :p_hat 
The function $p$ satisfies the state system \eqref{e1}--\eqref{e4} 
if and only if $\hat{p}$ is a solution to the  equation
\begin{equation}
\partial_t\hat{p}-k\Delta \hat{p}
= \partial_s(d(s)\partial_s \hat{p}-g(s)\hat{p})-
\mu \hat{p}-u\hat{p}-\theta \hat{p}+\int_{S_0}^{S_1}
\beta(x, s, t, \hat{s})\hat{p}(x, \hat{s}, t)\,{\rm d}\hat{s},\label{e1'}
\end{equation}
endowed with the analogical initial-boundary conditions:
\begin{gather}
  \frac{\partial \hat{p}}{\partial s}(x, S_0, t)
=\frac{\partial \hat{p}}{\partial s}(x, S_1, t)=0,\quad \forall 
(x, t)\in\Omega\times(0,T),\label{e2'} \\
\frac{\partial \hat{p}}{\partial n}(x, s, t)=0,\quad\text{on }\Sigma,\label{e3'}\\
\hat{p}(x, s,0)=p^0(x, s),\quad\forall (x,s)\in\Omega\times(S_0,S_1).\quad 
\label{e4'}
\end{gather}
\end{proposition}

Multiplying  the equation \eqref{e1'} by a function $q$ and using 
integration by parts on ${\mathcal Q}$, we arrive at the following definition.

\begin{definition}\label{def3} \rm % : a(t;p,q)
The bilinear mapping $a(t;\cdot,\cdot): H^1(\mathcal Q))
\times H^1(\mathcal Q))\to \mathbb{R}$ for $t\in[0,T]$,
is defined as
\begin{equation}
a(t;\hat{p},q)=\int_{\mathcal Q}\left(k\nabla \hat{p}
\cdot\nabla q+d(\partial_s \hat{p})(\partial_s q)
+g(\partial_s \hat{p})q+(\mu+u+\theta+g_s)\hat{p}q\right)\,{\rm d}x\,{\rm d}s.
\label{e13}
\end{equation}
\end{definition}

According to classical discussions (see, e.g. \cite{Wu2006}), we cite the 
following result and  omit the proof.

\begin{lemma}\label{lem4} %  a(t;.,.)
 For almost every $t\in(0,T)$, $a(t;\hat{p},q)$ is continuous on 
$H^1(\mathcal Q)\times H^1(\mathcal Q)$, and for $\theta$ large enough,
 $a(t;\hat{p},q)$ is coercive on $H^1(\mathcal Q)$. There exist two 
constants $M>0$ and $\delta>0$, depending on $k$, $\overline{d}$, 
$\|\mu\|_{\infty}$, $|g_s|_{\rm max}$, $d_1$, and $\theta$, such that
\begin{gather}
|a(t;\hat{p},q)|\leq M\|\hat{p}\|_{H^1(\mathcal Q)}\|q\|_{H^1(\mathcal Q)}, 
\quad \forall \hat{p}, q\in {H^1(\mathcal Q)}, \label{e14} \\
a(t;\hat{p},\hat{p})\geq \delta \|\hat{p}\|_{H^1(\mathcal Q)}^2, 
\quad \forall \hat{p}\in {H^1(\mathcal Q)}. \label{e15}
\end{gather}
\end{lemma}

Now we are ready to define the weak solutions $\hat{p}$ to 
\eqref{e1'}--\eqref{e4'}.

\begin{definition}\label{def5} \rm % weak hat p 
A function $\hat{p}\in W(0,T)$ is said to be a solution of \eqref{e1'}--\eqref{e4'} 
if the following variational equation holds for all $
q\in L^2(0,T; H^1(\mathcal Q))$:
\begin{equation}
\int_0^T \Big(\frac{d\hat{p}}{dt}, q \Big)_Hdt
+\int_0^Ta(t;\hat{p},q)dt=\int_0^T(\mathcal{I}\hat{p}, q)_Hdt,\label{e16}
\end{equation}
and
$$
\hat{p}(x,s,0)=p^0(x,s)\ \text{in}\  \Omega\times(S_0,S_1),
$$
where $\mathcal{I}\hat{p}:=\int_{S_0}^{S_1}\beta(x, s, t, \hat{s})
\hat{p}(x, \hat{s}, t)\,{\rm d}\hat{s}$.
\end{definition}

\begin{lemma}\label{lem6} % :p_hat exitence
 System \eqref{e1'}--\eqref{e4'} has a unique non-negative bounded 
solution $\hat{p}\in W(0,T)$.
\end{lemma}

\begin{proof}
 Firstly, we define an operator $\mathcal{A}$ by freezing the integral 
term $\mathcal{I}\hat{p}$, and then apply the Banach fixed-point theorem 
to $\mathcal{A}$.
So it is clear to see that the fixed point is our desired solution.

Let $\hat{p}^*$ be fixed in $W(0,T)$ and replace $(\mathcal{I}\hat{p}, q)_H$ 
by $(\mathcal{I}\hat{p}^*, q)_H$ in \eqref{e16}. 
For all $q\in L^2(0,T; H^1(\mathcal Q))$ and some appropriate $T$,  
the problem reduces to the following standard linear problem in 
the sense of distribution:
\begin{equation}
\begin{gathered}
 \Big(\frac{d\hat{p}}{dt}, q \Big)_H+a(t;\hat{p},q)
=(\mathcal{I}\hat{p}^*, q)_H,\\
 \hat{p}(x,s,0)=p^0(x,s).
\end{gathered} \label{e17}
\end{equation}
We get a unique solution $\hat{p}\in W(0,T)$ of the problem \eqref{e17} 
by the classical discussion. So this solution defines an operator 
$\mathcal{A}$ on $W(0,T)$ and $\mathcal{A}\hat{p}^*=\hat{p}$.

Taking $q=\hat{p}$ in \eqref{e17}, integrating it on $[0,t]$, 
using the coerciveness of $a(t; p,q)$ in (15) and Cauchy-Schwarz inequality, 
we have
\begin{equation}
\int_0^t\Big(\frac{1}{2}\frac{d}{dt}\|\hat{p}(\tau)\|_{H}^2
+\delta\|\hat{p}(\tau)\|_{H_1}^2\Big)d\tau
\leq \int_0^t\|\mathcal{I}\hat{p}^*(\tau)\|_H\cdot\|\hat{p}(\tau)\|_Hd\tau.
\label{e18}
\end{equation}

By using Young's inequality, for all $\alpha>0$ we obtain
\begin{equation}
\frac{1}{2}\left(\|\hat{p}(t)\|_{H}^2-\|p^0\|_{H}^2\right)
+\delta\int_0^t\|\hat{p}(\tau)\|_{H_1}^2d\tau
\leq \int_0^t\frac{1}{\alpha} \|\mathcal{I}\hat{p}^*(\tau)\|_H^2\,d\tau
+\int_0^t\alpha\|\hat{p}(\tau)\|_H^2\,d\tau.\label{e19}
\end{equation}
Choosing $\alpha=\delta$, by the norm definition of $H_1$ in \eqref{e7} 
and the assumption (H3), we derive that
\begin{equation}
\begin{aligned}
\|\hat{p}(t)\|_{H}^2
&\leq \frac{2}{\delta}\int_0^t\|\mathcal{I}\hat{p}^*(\tau)
 \|_H^2d\tau+2\|p^0\|_{H}^2\\ 
& \leq \frac{2}{\delta}\int_0^T 
\int_{\mathcal{Q}}\Big(\int_{S_0}^{S_1}\beta(x,s,\tau,
\hat{s})\hat{p}^*(x,\hat{s},\tau)d\hat{s}\Big)^2dx\,ds\,d\tau
 +2\|p^0\|_{H}^2
\\ 
& \leq \frac{2\overline{\beta}^2(S_1-S_0)^2T}{\delta}
\|\hat{p}^*\|_{H}^2+2\|p^0\|_{H}^2.
\end{aligned} \label{e20}
\end{equation}
Thus, we have
\begin{equation}
\|\hat{p}(t)\|_{L^{\infty}(0,T;H)}^2
\leq \frac{2\overline{\beta}^2(S_1-S_0)^2T}{\delta}
\|\hat{p}^*\|_{L^{\infty}(0,T;H)}^2+2\|p^0\|_{H}^2.\label{e21}
\end{equation}

Define a ball domain 
\begin{equation}
B_r:=\big\{p\in W(0,T):\|p\|_{L^{\infty}(0,T;H)}\leq r, 
r\geq \frac{\|p^0\|_H}{\sqrt{\frac{1}{2}
-\frac{\overline{\beta}^2(S_1-S_0)^2T}{\delta}}}\big\},\label{e22}
\end{equation}
where  $T<\delta/(2\overline{\beta}^2(S_1-S_0)^2)$.
Then we have $\mathcal{A}B_r\subset B_r$ by \eqref{e21}, because if 
$\|\hat{p}^*\|_{L^{\infty}(0,T;H)}\leq r$, it gives 
$\|\hat{p}\|_{L^{\infty}(0,T;H)}\leq r$ for 
$\frac{2\overline{\beta}^2(S_1-S_0)^2T}{\delta}r^2+2\|p^0\|_{H}^2\leq r^2$  
from \eqref{e22}.

 Furthermore, we claim that $\mathcal{A}$ is a strict contraction on $B_r$. In fact, 
Let $\mathcal{A}\hat{p}_i^*=\hat{p}_i,\, \hat{p}_i,\hat{p_i}^*\in B_r, i=1,2$. 
By using a similar deduction from \eqref{e17} to \eqref{e21} to 
$\hat{p}_1-\hat{p}_2$, we have
\begin{equation}
\|\hat{p}_1-\hat{p}_2|_{L^{\infty}(0,T;H)}^2
\leq \frac{2\overline{\beta}^2(S_1-S_0)^2T}{\delta}
\|\hat{p}_1^*-\hat{p}_2^*\|_{L^{\infty}(0,T;H)}^2.\label{e23}
\end{equation}
For $T<\frac{\delta}{2\overline{\beta}^2(S_1-S_0)^2}$, namely, 
$\frac{2\overline{\beta}^2(S_1-S_0)^2T}{\delta}<1$, $\mathcal{A}$ 
is a strict contraction. Banach fixed-point theorem allows us to 
conclude that there exists a unique $\hat{p}\in B_r$ such that 
$\mathcal{A}\hat{p}=\hat{p}$. This point is the desired unique bounded
 solution $\hat{p}\in W(0,T)$.

Since $T$ does not depend on $p^0$, we can apply the same procedure as 
the above on $(T,2T)$, $(2T,3T)$, \dots and so on. So we deduce that 
a solution of \eqref{e1'}--\eqref{e4'}
 can be found on the desired time interval.

We now prove the non-negativity. Let $\hat{p}_1\geq 0$ be given in $W(0,T)$
 and define the sequence $\{\hat{p}_n\}_{n\geq 1}$ by 
$\mathcal{A}\hat{p}_n=\hat{p}_{n+1}$.

Taking $\hat{p}_2^-=\max\{0,-\hat{p}_2\}$ as a test function in \eqref{e17} 
leads to
\begin{equation}
\Big(\frac{d\hat{p}_2}{dt}, \hat{p}_2^-\Big)_H+a(t;\hat{p}_2,\hat{p}_2^-)
=(\mathcal{I}\hat{p}_1, \hat{p}_2^-)_H.\label{e24}
\end{equation}

If we let $\hat{p}_2^+=\max\{0,\hat{p}_2\}$, then 
$\hat{p}_2=\hat{p}_2^+-\hat{p}_2^-$ and $\hat{p}_2^+\cdot\hat{p}_2^-=0$, 
and it gives
\begin{equation}
\frac{1}{2}\frac{d}{dt}\|\hat{p}_2^-\|_H^2
\leq \frac{1}{2}\frac{d}{dt}\|\hat{p}_2^-\|_H^2+a(t;\hat{p}_2^-,\hat{p}_2^-)
=-(\mathcal{I}\hat{p}_1, \hat{p}_2^-)_H.\label{e25}
\end{equation}
Since $\hat{p}_1\geq 0$, it leads to $\mathcal{I}\hat{p}_1\geq 0$ and 
$-(\mathcal{I}\hat{p}_1, \hat{p}_2^-)_H\leq 0$. Then we find 
$\frac{d}{dt}\|\hat{p}_2^-\|_H^2\leq 0$, and
\begin{equation}
\|\hat{p}(t)_2^-\|_H^2\leq \|\hat{p}(0)_2^-\|_H^2
=\|p^{0-}\|_H^2=0,\label{e26}
\end{equation}
which means that $\hat{p}_2\geq 0$. By induction, we can further show 
that $\hat{p}_n\geq 0$, for all $n\geq 1$. The unique solution $\hat{p}\in W(0,T)$
(the limiting  point of the sequence) is non-negative.
\end{proof}

From Lemma \ref{lem6} and Proposition \ref{prop2}, we have the following result.

\begin{theorem}\label{thm7} % p_hat exitence 
Assume that the hypotheses {\rm (H1)--(H5)} hold. 
For any $u\in\mathcal{U}$, the system \eqref{e1}--\eqref{e4} has a unique 
nonnegative solution $p^u(x,s,t)\in W(0,T)$ in $Q$ and
\begin{equation}
0\leq p^u(x,s,t)\leq M_1, \quad \text{a.e. in }  Q, \label{e25b}
\end{equation}
where $M_1>0$ is a constant independent of $p^u$ and $u$.
\end{theorem}

\section{Optimal strategies}

In this section, we derive the first-order necessary optimality conditions 
for the optimal harvesting control problem \eqref{e5}.

We present an auxiliary result for the continuous dependence of the population 
density with the harvesting effort $u$ reads as follows.

\begin{lemma}\label{lem8} % :pu1-pu2
Let $p^{u_1}, p^{u_2}$ be the solutions of \eqref{e1}--\eqref{e4} 
corresponding to the controls $u_1, u_2\in\mathcal{U}$, respectively. 
Then we have
\begin{equation}
|p^{u_1}-p^{u_2}|\leq TC_1\|u_1-u_2\|_{L^{\infty}(Q)},\label{e26b}
\end{equation}
where $C_1$ is a positive constant independent of $u_1$ and $u_2$.
\end{lemma}

\begin{proof} 
Let $y=p^{u_1}-p^{u_2}$. Then $y$ is the solution of the  system
\begin{equation}
\begin{gathered}
\begin{aligned}
\partial_ty-k\Delta y &= \partial_s(d(s)\partial_s y-g(s)y)
 -\mu y-u_1y+(u_2-u_1)p^{u_2} \\
&\quad +\int_{S_0}^{S_1}\beta(x, s, t, \hat{s})
 y(x, \hat{s}, t)\,{\rm d}\hat{s}, 
\end{aligned}\\
\frac{\partial y}{\partial s}(x, S_0, t)=\frac{\partial y}{\partial s}(x, S_1, t)=0,
\quad \forall (x, t)\in\Omega\times(0,T),\\
\frac{\partial y}{\partial n}(x, s, t)=0,\quad \text{on }\Sigma,\\
 y(x, s,0)=0,\quad \forall (x,s)\in\Omega\times(S_0,S_1).
\end{gathered} \label{e27}
\end{equation}
Multiplying the first equation of \eqref{e27} by $y$ and integrating it 
on $Q_t:=\Omega\times (S_0,S_1)\times(0,t)$, we deduce that
\begin{equation}
\begin{aligned}
&\|y(\cdot, t, \cdot)\|_H^2  \\
&\leq \int_0^t\int_{\Omega}g(S_0)y(x,S_0,\tau)^2dxd\tau-\int_{Q_t}g'(s)y^2d\sigma\\
&\quad +2\int_{Q_t}(u_2-u_1)p^{u_2}yd\sigma
 +2\int_{Q_t}\Big(\int_{S_0}^{S_1}\beta(x,s,\tau,\hat{s})
  y(x,\hat{s},\tau)d\hat{s}\Big)yd\sigma\\
&\leq |g'(s)|_{\rm max}\int_0^t\|y(\cdot, \tau, \cdot)\|_H^2d\tau
 +2\int_{Q_t}|u_1-u_2||M_1y|d\sigma\\
  &\quad +2\overline{\beta}(S_1-S_0)\int_0^t\|y(\cdot, \tau, \cdot)\|_H^2d\tau\\
  &\leq |g'(s)|_{\rm max}\int_0^t\|y(\cdot, \tau, \cdot)\|_H^2d\tau+\int_{Q_t}\left(|M_1(u_1-u_2)|^2+|y|^2\right)d\sigma\\
  &\quad +2\overline{\beta}(S_1-S_0)\int_0^t\|y(\cdot, \tau, \cdot)\|_H^2d\tau\\
  &=M_1^2\int_0^t\|u_1-u_2\|_H^2d\tau+\Big(|g'(s)|_{\rm max}
 +2\\
&\quad +2\overline{\beta}(S_1-S_0)\Big)\int_0^t\|y(\cdot, \tau, \cdot)\|_H^2d\tau.
\end{aligned}\label{e28}
\end{equation}

It follows  from Bellman's lemma that
\begin{equation}
\|y(\cdot, t, \cdot)\|_H^2\leq M_1^2e^{\left(|g'(s)|_{\rm max}
+1+2\overline{\beta}(S_1-S_0)\right)T}\|u_1-u_2\|_H^2.\label{e29}
\end{equation}
Integrating it on $(0,T)$ yields
\begin{equation}
\|y\|_{L^2(0,T;H)}^2\leq TM_1^2e^{\left(|g'(s)|_{\rm max}
+1+2\overline{\beta}(S_1-S_0)\right)T}\|u_1-u_2\|_{L^2(0,T;H)}^2.\label{e30}
\end{equation}
Thus, by the fundamental embedding inequality, we know that \eqref{e26}
 holds for some  constant $C_1>0$.
\end{proof}

To characterize the structure of the optimal controller, we need to define 
the following dual problem associated with \eqref{e1}--\eqref{e4}:
\begin{equation}
\begin{gathered}
\begin{aligned}
\partial_tq+k\Delta q
&= -\partial_s(d(s)\partial_s q)-g(s)\partial_s q+(\mu +u^*)q+wu^* \\
&\quad -\int_{S_0}^{S_1}\beta(x, \hat{s}, t, s)q(x,\hat{s}, t)\,{\rm d}\hat{s},
\end{aligned}\\
d(s)\partial_sq+g(s)q|_{s=S_0}
=d(s)\partial_sq+g(s)q|_{s=S_1}=0,\quad
 \forall (x, t)\in\Omega\times(0,T),
\\
  \frac{\partial q}{\partial n}(x, s, t)=0,\quad \text{on }\Sigma,\\
 q(x, s,T)=0,\quad \forall (x,s)\in\Omega\times(S_0,S_1).
\end{gathered}\label{e31}
\end{equation}

Under the changes $\tau=T-t$ and $\tilde{q}(x,s,\tau):=q(x,s,T-\tau)$, the above problem becomes
\begin{equation}
\begin{gathered}
\begin{aligned}
\partial_{\tau}\tilde{q}-k\Delta \tilde{q}
&= \partial_s(d(s)\partial_s \tilde{q})+g(s)\partial_s \tilde{q}
-(\mu +u^*)\tilde{q}-wu^*\\
&\quad +\int_{S_0}^{S_1}\beta(x, \hat{s}, \tau, s)
\tilde{q}(x,\hat{s}, \tau)\,{\rm d}\hat{s},
\end{aligned} \\
d(s)\partial_s\tilde{q}+g(s)\tilde{q}|_{s=S_0}
=d(s)\partial_s\tilde{q}+g(s)\tilde{q}|_{s=S_1}=0,\quad \forall 
(x, \tau)\in\Omega\times(0,T),\\
  \frac{\partial \tilde{q}}{\partial n}(x, s, \tau)=0,\quad \text{on }\Sigma,\\
 \tilde{q}(x, s,0)=0,\quad\forall (x,s)\in\Omega\times(S_0,S_1).
\end{gathered}\label{e32}
\end{equation}
Using classical results for parabolic equations associated with \eqref{e32}, 
and discussing in the same manner as that in Lemmas \ref{lem6} and \ref{lem8}, 
we can derive the following lemma.


\begin{lemma}\label{lem9} % :adjoint q
Problem \eqref{e31} has a unique solution 
$q^u\in L^2(0,T;H^1(\mathcal{Q}))$ and 
\begin{equation}
|q^u(x,s,t)|\leq M_2,\quad\text{a.e. in }Q, \label{e33}
\end{equation}
where $M_2$ is a positive constant independent of $q^u$ and $u$.

Furthermore, let $q^{u_1},\ q^{u_2}$ be the solutions of  \eqref{e31} 
corresponding to $u_1, u_2\in\mathcal{U}$, respectively. 
Then there exists a positive constant $C_2$, which is independent of 
$u_1, u_2$, such that
\begin{equation}
|q^{u_1}-q^{u_2}|\leq TC_2\|u_1-u_2\|_{L^{\infty}(Q)}.\label{e34}
\end{equation}
\end{lemma}

We now describe the structure of optimal controllers as follows.

\begin{theorem}\label{thm10} % :u* optimal control
Let $u^*(x,s,t)\in\mathcal{U}$ be an optimal control for the problem 
\eqref{e1}--\eqref{e5}, and $p^{u^*}$ and $q^{u^*}$ be the corresponding 
solutions of system \eqref{e1}--\eqref{e4} and \eqref{e31}, respectively. 
Then we have
\begin{equation}
u^*(x,s,t)=\mathcal{F}\big\{\frac{[w+q^{u^*}]p^{u^*}}{\rho}\big\}(x,s,t),
\label{e35}
\end{equation}
in which the mapping $\mathcal{F}$ is defined as
\begin{equation}
(\mathcal{F}h)(x,s,t)
=\begin{cases}
\zeta_1(x,s,t), &h(s,t,x)<\zeta_1(x,s,t),\\
h(x,s,t),&\zeta_1(x,s,t)\leq h(x,s,t)\leq\zeta_2(x,s,t),\\
\zeta_2(x,s,t),&h(x,s,t)>\zeta_2(x,s,t).
\end{cases} \label{e36}
\end{equation}
\end{theorem}

\begin{proof} 
Let $\mathcal{T}_{\mathcal{U}}(u^*)$ be the tangent cone to 
$\mathcal{U}$ at $u^*$(see \cite{Aubin1984}). 
For any $v\in\mathcal{T}_{\mathcal{U}}(u^*)$, we know that 
$u^*+\varepsilon v\in\mathcal{U}$ for the sufficient small $\varepsilon>0$. 
Since $u^*$ is optimal, it follows that
\begin{equation}
\begin{aligned}
&\int_Q\big(wu^*p^{u^*}-\frac{1}{2}\rho u^{*2}\big)\,dx\,ds\,dt \\
&\geq \int_Q \big(w(u^*+\varepsilon v)p^{u^*+\varepsilon v}
-\frac{1}{2}\rho (u^*+\varepsilon v)^2\big)\,dx\,ds\,dt,
\end{aligned} \label{e37}
\end{equation}
which implies
\begin{equation}
  \int_Q\Big(wu^*\frac{p^{u^*+\varepsilon v}-p^{u^*}}{\varepsilon}
+wvp^{u^*+\varepsilon v}-\frac{1}{2}\rho v(2u^*+\varepsilon v)\Big)
\,dx\,ds\,dt\leq 0.\label{e38}
\end{equation}

Let $z(x,s,t)$ be the solution of
\begin{equation}
\begin{gathered}
\begin{aligned}
\partial_tz-k\Delta z
&= \partial_s(d(s)\partial_s z-g(s)z)
-(\mu +u^*)z-vp^{u^*} \\
&\quad +\int_{S_0}^{S_1}\beta(x, s, t, \hat{s})
z(x, \hat{s}, t)\,{\rm d}\hat{s}, 
\end{aligned} \\
\frac{\partial z}{\partial s}(x, S_0, t)
=\frac{\partial z}{\partial s}(x, S_1, t)=0,\quad
 \forall (x, t)\in\Omega\times(0,T),\\
  \frac{\partial z}{\partial n}(x, s, t)=0,\quad \text{on }\Sigma,\\
 z(x, s,0)=0, \quad\forall (x,s)\in\Omega\times(S_0,S_1).
\end{gathered} \label{e39}
\end{equation}
The existence and uniqueness of solutions to \eqref{e39}
 follows from the theory of nonhomogeneous parabolic equations 
(see, e.g. \cite{Evans2010}). 

Let
\begin{equation}
w_\varepsilon(x,s,t)=\frac{p^{u^*+\varepsilon v}-p^{u^*}}{\varepsilon}-z(x,s,t),
\quad (x,s,t)\in Q.\label{e40}
\end{equation}
It is not hard to deduce that $w_{\varepsilon}(x, s, t)$ is the solution of
\begin{equation}
\begin{gathered}
\begin{aligned}
\partial_tw-k\Delta w &= \partial_s(d(s)\partial_s w-g(s)w)-\mu w-u^*w\\
 &\quad -v(p^{u^*+\varepsilon v}-p^{u^*})+\int_{S_0}^{S_1}
 \beta(x, s, t, \hat{s})w(x, \hat{s}, t)\,{\rm d}\hat{s},
\end{aligned} \\
\frac{\partial w}{\partial s}(x, S_0, t)
=\frac{\partial w}{\partial s}(x, S_1, t)=0,\quad
 \forall (x, t)\in\Omega\times(0,T),\\
  \frac{\partial w}{\partial n}(x, s, t)=0,\quad \text{on }\Sigma,\\
 w(x, s,0)=0,\quad \forall (x,s)\in\Omega\times(S_0,S_1).
\end{gathered}\label{e41}
\end{equation}
In what follows, we show $w_{\varepsilon}\to 0$ as $\varepsilon\to 0+$. 
By estimating \eqref{e26}, we may infer
\begin{equation}
p^{u^*+\varepsilon v}-p^{u^*}\to 0,\quad\text{in $L^2(0,T;H)$ as }
 \varepsilon\to 0^+.\label{e42}
\end{equation}
We now consider the  limit system
\begin{equation}
\begin{gathered}
\partial_tw-k\Delta w= \partial_s(d(s)\partial_s w-g(s)w)-\mu w-u^*w
+\int_{S_0}^{S_1}\beta(x, s, t, \hat{s})w(x, \hat{s}, t)\,{\rm d}\hat{s}, \\
\frac{\partial w}{\partial s}(x, S_0, t)
=\frac{\partial w}{\partial s}(x, S_1, t)=0,\quad
 \forall (x, t)\in\Omega\times(0,T),\\
  \frac{\partial w}{\partial n}(x, s, t)=0,\quad \text{on }\Sigma,\\
 w(x, s,0)=0,\quad \forall (x,s)\in\Omega\times(S_0,S_1),
\end{gathered}\label{e43}
\end{equation}
which is a homogenous linear parabolic system and has a unique solution 
$w(x,s,t)=0$ a.e. in $Q$. Hence, we have
\begin{equation}
\frac{p^{u^*+\varepsilon v}-p^{u^*}}{\varepsilon}\to z,
\quad\text{in }L^2(0,T;H)\ \mathrm{as } \varepsilon\to 0^+.\label{e44}
\end{equation}
Passing to the limit in \eqref{e38} we find
\begin{equation}
\int_Q\big(wu^*z+(wp^{u^*}-\rho u^*)v\big)\,dx\,ds\,dt\leq 0.\label{e45}
\end{equation}
Multiplying \eqref{e31} by $z(x,s,t)$ and integrating it over $Q$ 
(denoting $u$ by $u^*$), we deduce that
\begin{equation}
\int_Q wu^*z\,dx\,ds\,dt=\int_Q vp^{u^*}q^{u^*}\,dx\,ds\,dt.\label{e46}
\end{equation}
Then it follows from \eqref{e45} and \eqref{e46} that
\begin{equation}
  \int_Q \left \{\left [(w+q^{u^*})p^{u^*}-\rho u^*\right]v\right\}(s,t,x)
\,\mathrm{d}x\,\mathrm{d}t\,\mathrm{d}s\leq 0, \quad
\forall  v\in\mathcal{T}_{\mathcal{U}}(u^*).\label{e47}
\end{equation}
According to the properties of normal cone (see [22]), 
the expression in the square brackets of \eqref{e47} satisfies 
$(w+q^{u^*})p^{u^*}-\rho u^*\in\mathcal{N}_{\mathcal{U}}(u^*)$, 
the normal cone to $\mathcal{U}$ at $u^*$. Consequently the conclusion follows.
\end{proof}


\section{Existence and uniqueness of optimal solutions}

In this section, we show that there is one and only one solution for optimal 
harvesting control problem \eqref{e5}.
We need the following lemma which can be proven by the definition of normal 
cones (see, e.g. \cite{Barbu1999}).

\begin{lemma}\label{lem11} % : tangent normal cone
Suppose that $\eta(x,s,t)\in L^1(Q)$ satisfies
\begin{equation}
\int_Q[\eta(x,s,t)v(x,s,t)+\alpha|v(x,s,t)|]\,dx\,ds\,dt\geq 0, \quad 
\forall v\in\mathcal{T}_{\mathcal{U}}(u), \label{e48}
\end{equation}
where $\alpha$ is some small positive constant. 
Then there exists some $\theta\in L^{\infty}(Q)$ such that 
$|\theta|_{\infty}\leq 1$ and $u+\alpha\theta\in\mathcal{N}_{\mathcal{U}}(u)$.
\end{lemma}

The following result guarantees  the existence and uniqueness of the 
optimal strategies.

\begin{theorem}\label{thm12} % :existence of u*
Assume that {\rm (H1)--(H6)} hold. If
\begin{equation}
T(C_1(W+M_2)+C_2M_1)\rho^{-1}<1,\label{e49}
\end{equation}
where $W$ is the same as in {\rm (H6)}, and 
$M_i,C_i$ $(i=1,2)$ are given in Theorem \ref{thm7}, Lemmas \ref{lem8} 
and \ref{lem9}, then the optimal control problem \eqref{e1}--\eqref{e5} 
has a unique solution.
\end{theorem}

\begin{proof} 
Define a functional $\Phi: L^1(Q)\to (-\infty,+\infty]$ by
\begin{equation}
\Phi(u)=\begin{cases}
-J(u) =\int_Q\left(\frac{1}{2}\rho u^2-wp^uu\right)\,dx\,ds\,dt,
&\text{if } u\in\mathcal{U},\\
+\infty,&\text{if }u\notin\mathcal{U},
\end{cases} \label{e50}
\end{equation}
where $J(\cdot)$ is of the form \eqref{e5}.  
By Lemma \ref{lem8}, it is easily seen that $\Phi$ is lower semi-continuous. 
According to the Ekeland variational principle \cite{Ekeland1974}, 
for each $\varepsilon>0$ there exists $u_{\varepsilon}\in\mathcal{U}$ such that
\begin{gather}
\Phi(u_{\varepsilon})\leq \inf_{u\in\mathcal{U}}\Phi(u)+\varepsilon, \label{e51} \\
\Phi(u_{\varepsilon})\leq \inf_{u\in\mathcal{U}}\{\Phi(u)
+\sqrt{\varepsilon}|u-u_{\varepsilon}|_1\}, \label{e52}
\end{gather}
where $|\cdot|_1$ denotes the norm in $L^1(Q)$.

Note that the perturbed functional 
$\Phi_{\varepsilon}(u):=\Phi(u)+\sqrt{\varepsilon}|u-u_{\varepsilon}|_1$ 
attains its infimum at $u_{\varepsilon}$. 
By the same argument as in the previous section, we obtain the 
 condition
\begin{equation}
\int_Q\big(\rho u_{\varepsilon}-(w+q^{u_{\varepsilon}})p^{u_{\varepsilon}}
\big)v\,dx\,ds\,dt
+\sqrt{\varepsilon}\int_Q|v(x,s,t)|\,dx\,ds\,dt\geq 0, \quad \forall 
v\in \mathcal{T}_{\mathcal{U}}(u_{\varepsilon}).\label{e53}
\end{equation}
By Lemma \ref{lem11}, we see that there exists $\theta_{\varepsilon}\in L^{\infty}(Q)$, 
and $|\theta|_{\infty}\leq 1$, such that
\begin{equation}
\rho u_{\varepsilon}-(w+q^{u_{\varepsilon}})p^{u_{\varepsilon}}
-\sqrt{\varepsilon}\theta_{\varepsilon}\in\mathcal{N}_{\mathcal{U}}
(u_{\varepsilon}), \label{e54}
\end{equation}
and consequently, 
\begin{equation}
u_{\varepsilon}(x,s,t)=\mathcal{F}[(1/\rho)
\left((w+q^{u_{\varepsilon}})p^{u_{\varepsilon}}
+\sqrt{\varepsilon}\theta_{\varepsilon}\right)], \quad \text{a.e. in }
 Q.\label{e55}
\end{equation}

To show the uniqueness of the optimal controller $u$, we define 
$\mathcal{J}: \mathcal{U}\subset L^{\infty}(Q)\to \mathcal{U}$ by
\begin{equation}
(\mathcal{J}(u))(x,s,t)
=\mathcal{F}\Big(\frac{w(x,s,t)+q^u(x,s,t)}{\rho}p^u(x,s,t)\Big),\quad
 \text{a.e. in } Q. \label{e56}
\end{equation}
For  $(x,s,t)\in Q$ we have
\begin{equation}
\begin{aligned}
&|(\mathcal{J}(u))(x,s,t)-(\mathcal{J}(v))(x,s,t)| \\
& =|\mathcal{F}((1/\rho)(w+q^u)p^u)-\mathcal{F}((1/\rho)(w+q^v)p^v)|\\
& \leq\frac{W}{\rho}|p^u-p^v|+\frac{1}{\rho}|q^u||p^u-p^v|
 +\frac{1}{\rho}|p^v||q^u-q^v|.
\end{aligned}\label{e57}
\end{equation}
By the estimates \eqref{e25}, \eqref{e26}, \eqref{e33} and \eqref{e34},
 we obtain
\begin{equation}
\begin{aligned}
&\|(\mathcal{J}(u))(x,s,t)-(\mathcal{J}(v))(x,s,t)\|_{L^{\infty}(Q)} \\
&\leq T(C_1(W+M_2)+C_2M_1)\rho^{-1}\|u-v\|_{L^{\infty}(Q)}.
\end{aligned}\label{e58}
\end{equation}
So $\mathcal{J}$ is a contraction if $T(C_1(W+M_2)+C_2M_1)\rho^{-1}<1$, 
and $\mathcal{J}$ has one and only one fixed point $\overline{u}\in\mathcal{U}$. 
Theorem \ref{thm10} implies that an optimal controller $u^*$, if it exists, must
 coincide with this fixed point. Hence, the uniqueness of optimal controls
 is proved.

Next, we prove the existence of optimal controls. In fact, we only need 
to show that the fixed point $\overline{u}$ minimizes $\Phi(\cdot)$.
By \eqref{e55} and \eqref{e56}, we have
\begin{equation}
\begin{aligned}
&\|\mathcal{J}(u_{\varepsilon})-u_{\varepsilon}\|_{L^{\infty}(Q)}\\
&=\|\mathcal{F}(\rho^{-1}(w+q^u)p^u)-\mathcal{F}[\rho^{-1}
((w+q^{u_{\varepsilon}})p^{u_{\varepsilon}}
+\sqrt{\varepsilon}\theta_{\varepsilon})]\|_{L^{\infty}(Q)}\\
&\leq \rho^{-1}\sqrt{\varepsilon}.
\end{aligned}\label{e59}
\end{equation}
This leads to
\begin{equation}
\begin{aligned}
\|\overline{u}-u_{\varepsilon}\|_{L^{\infty}(Q)}
&\leq \|\mathcal{J}(\overline{u})-\mathcal{J}(u_{\varepsilon})
\|_{L^{\infty}(Q)}+\rho^{-1}\sqrt{\varepsilon}\\ 
&\leq T(C_1(W+M_2)+C_2M_1)\rho^{-1}\|\overline{u}
-u_{\varepsilon}\|_{L^{\infty}(Q)}+\rho^{-1}\sqrt{\varepsilon};
\end{aligned}\label{e60}
\end{equation}
that is,
\begin{equation}
\|\overline{u}-u_{\varepsilon}\|_{L^{\infty}(Q)}
\leq [1-T(C_1(W+M_2)+C_2M_1)\rho^{-1}]^{-1}\rho^{-1}\sqrt{\varepsilon}.\label{e61}
\end{equation}
So, we see that $u_{\varepsilon}\to\overline{u}$ in $L^{\infty}(Q)$, 
and by \eqref{e51} we have $\Phi(\overline{u})=\inf_{u\in\mathcal{U}}\Phi(u)$ 
which completes the proof. 
\end{proof}

\subsection*{Conclusions and comments}

In this article, we introduced a linear structured population model 
with spatial diffusion, size random growth (namely, diffusion in the 
size space) and the distributed recruitment.  
Application of the size diffusion is natural and significant in the 
biological phenomena \cite{Hadeler2010, Farkas2011}, since individuals 
that have the same size initially, may disperse as time progresses. 
We equipped our model with the homogeneous Neumann boundary condition 
with respect to the $N$-dimension spatial variable $x$  and 1-dimension 
size variable $s$, and presented the results of the existence and 
uniqueness of solution of the state system, which laid a sufficient 
foundation for the optimal harvesting control problems.

By applying the Ekeland variational principle \cite{Ekeland1974} and 
the properties of normal cone and adjoint  techniques \cite{Barbu1999}, 
we developed the optimal harvesting strategies and deduced the conditions 
to assure only one optimal policy. Theorem \ref{thm12} tells us that, under the 
given conditions, our optimal control problem admits one and only one solution. 
Furthermore, Theorem \ref{thm10} described
the structure, other than a specific analytical expression, of optimal strategy. 
Unfortunately, one could not derive the explicit formula for the optimal  
strategy since the strategy, the state and the costate are coupled into 
a complex system. The results at this stage may be regarded as a 
middle step to real world applications and serve as a starting point 
for numerical computations.

\subsection*{Acknowledgmetns}
This work is supported by the NSF of China under No. 11271104,
andby the  Natural Science Foundation of Zhejiang Province under 
No. LY13A010018 and No. LY16G010008.


\begin{thebibliography}{99}

\bibitem{Anitabook2000} S. Ani\c{t}a; 
\emph{Analysis and Control of Age-Dependent Population Dynamics}, 
Kluwer, Dordrecht, Netherlands, 2000.

\bibitem{Anitabook2010} S. Ani\c{t}a, V. Arna\v{u}tu, V. Capasso; 
\emph{An Introduction to Optimal Control Problems in Life Sciences 
and Economics: From Mathematical Models to Numerical Simulation with MATLAB}, 
Birkh\"auser, 2010.

\bibitem{Apreutesei2014} N. Apreutesei, G. Dimitriu, R. Strugariu; 
\emph{An optimal control problem for a two-prey and one-predator model
 with diffusion}, Computers and Mathematics with Applications, 
67 (2014), 2127-2143.

\bibitem{Aubin1984} J.-P. Aubin, I. Ekeland; 
\emph{Applied Nonlinear Analysis}, John Wiley \& Sons, New York, USA, 1984.

\bibitem{Barbu1999} V. Barbu, M. Iannelli; 
\emph{Optimal control of population dynamics}, 
Journal of optimization theory and applications, 102 (1999), 1-14.

\bibitem{Bart.2015} A. Bartlomiejczyk, H. Leszczy\'nski; 
\emph{Method of lines for physiologically structured models with diffusion}, 
Applied Numerical Mathematics 94 (2015), 140--148.

\bibitem{Chu2009} J. Chu, A. Ducrot, P. Magal, S. Ruan; 
\emph{Hopf bifurcation in a size-structured population dynamic model 
with random growth}, Journal of Differential Equations, 247 (2009), 956-1000.

\bibitem{Cushing1998} J. M. Cushing; 
\emph{An introduction to structured population dynamics}, 
SIAM, Philadelphia, 1998.

\bibitem{Ding2010} W. Ding, H. Finotti, S. Lenhart, Y. Lou, Q. Ye; 
\emph{Optimal control of growth coefficient on a steady-state population model}, 
Nonlinear Analysis: Real World Applications, 11 (2010), 688-704.

\bibitem{Ebenman1998} B. Ebenman, L. Persson; 
\emph{Size-structured populations: Ecology and Evolution}, 
Springer-Verlag, Berlin, 1988.

\bibitem{Ekeland1974} I. Ekeland; 
\emph{On the variational principle}, Journal of Mathematical Analysis and 
Applications, 47 (1974), 324-353.

\bibitem{Evans2010} L. C. Evans; 
\emph{Parital Differential Equations (2nd edition)}, 
American Mathematical Society, Rhode Island, USA, 2010.

\bibitem{Farkas2011} J. Z. Farkas, P. Hinow; 
\emph{Physiologically structured populations with diffusion and dynamic 
boundary conditions}, Mathematical Biosciences and Engineering, 
8 (2011), 503-513.

\bibitem{Faugeras2005} B. Faugeras, O. Maury; 
\emph{An advection-diffusion-reaction size-structured fish population 
dynamics model combined with a statistical parameter estimation procedure: 
application to the Indian Ocean skipjack tuna fishery}, 
Mathematical Biosciences and Engineering, 2 (2005), 719-741.

\bibitem{Hadeler2010} K. P. Hadeler; 
\emph{Structured populations with diffusion in state space}, 
Mathematical Biosciences and Engineering, 7 (2010), 37-49.

\bibitem{He2012} Z.-R. He, Y. Liu; 
\emph{An optimal birth control problem for a dynamical population model 
with size-structure},  Nonlinear Analysis: Real World Applications, 1
3 (2012), 1369-1378.

\bibitem{Liu2013} Y. Liu, X.-L. Cheng,  Z.-R. He; 
\emph{On the optimal harvesting of size-structured population dynamics}, 
Applied Mathematics B: A Journal of Chinese Universities, 28 (2013), 173-186.

\bibitem{Luo2007} Z. Luo; 
\emph{Optimal harvesting problem for an age-dependent n-dimensional food 
chain diffusion model}, Applied Mathematics and Computation, 186 (2007),
 1742-1752.

\bibitem{Magal2008} P. Magal, S. Ruan; 
\emph{Structured Population Models in Biology and Epidemiology}, 
Springer-Verlag, Berlin, 2008.

\bibitem{Troltzsch2010} F. Tr\"{o}ltzsch; 
\emph{Optimal control of partial differential equations: theory, methods, 
and applications (Translated by J\"{u}rgen Sprekels)}, 
American Mathematical Society, USA, 2010.


\bibitem{Wu2006} Z. Wu, J. Yin, C. Wang; 
\emph{Elliptic \& Parabolic Equations}, World Scientific, Singapore, 2006.


\bibitem{Zhao2005} C. Zhao, M. Wang,  P. Zhao; 
\emph{Optimal harvesting problems for age-dependent interacting species 
with diffusion}, Applied Mathematics and Computation, 163 (2005), 117-129.


\end{thebibliography}
\end{document}