\documentclass[reqno]{amsart} \usepackage{hyperref} \usepackage{graphicx} \AtBeginDocument{{\noindent\small \emph{Electronic Journal of Differential Equations}, Vol. 2012 (2012), No. 237, pp. 1--18.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu \newline ftp ejde.math.txstate.edu} \thanks{\copyright 2012 Texas State University - San Marcos.} \vspace{9mm}} \begin{document} \title[\hfilneg EJDE-2012/237\hfil Optimal control] {Optimal control applied to native-invasive species competition via a PDE model} \author[W. Ding, V. Hrynkiv, X. Mu \hfil EJDE-2012/237\hfilneg] {Wandi Ding, Volodymyr Hrynkiv, Xiaoyu Mu} % in alphabetical order \address{Wandi Ding \newline Department of Mathematical Sciences and Computational Science Program, Middle Tennessee State University, Murfreesboro, TN 37132, USA} \email{wandi.ding@mtsu.edu} \address{Volodymyr Hrynkiv \newline Department of Computer and Mathematical Sciences, University of Houston - Downtown, Houston, TX 77002, USA} \email{HrynkivV@uhd.edu} \address{Xiaoyu Mu \newline Department of Mathematics, University of Tennessee, Knoxville, TN 37996-1320, USA} \email{xiaoyumoon@gmail.com} \thanks{Submitted September 5, 2012. Published December 26, 2012.} \subjclass[2000]{49J20, 34K35, 92D25} \keywords{Optimal control; partial differential equations; native-invasive species; \hfill\break\indent salt cedar; cottonwood; spatial models} \begin{abstract} We consider an optimal control problem of a system of parabolic partial differential equations modelling the competition between an invasive and a native species. The motivating example is cottonwood-salt cedar competition, where the effect of disturbance in the system (such as flooding) is taken to be a control variable. Flooding being detrimental at low and high levels, and advantageous at medium levels led us to consider the quadratic growth function of the control. The objective is to maximize the native species and minimize the invasive species while minimizing the cost of implementing the control. An existence result for an optimal control is given. Numerical examples are presented to illustrate the results. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{lemma}[theorem]{Lemma} \newtheorem{definition}[theorem]{Definition} \allowdisplaybreaks \section{Introduction} \label{introduction} The problem of biological invasion as well as the interaction between the invasive and native species has drawn a great deal of attention from both biologists and mathematicians. The invasion of natural ecosystems by exotic species is an important part of global environmental change and poses a major threat to biodiversity \cite{nentwig}. In this context, understanding the underlying dynamics between all elements of the ecosystem plays an important role in eradication and control of the invasive species. In this paper we consider an optimal control problem for a system of parabolic partial differential equations (PDEs) modelling the competition between an invasive and a native species. As a motivating example, we take cottonwood - salt cedar competition, where the effect of disturbance in the system (such as flooding) is taken to be a control. Very few mathematical models have been applied to the cottonwood - salt cedar competition. We would like to point out the ordinary differential equations (ODEs) model for the cottonwood - salt cedar competition in \cite{miller} and the use of the optimal control approach there to find the optimal flooding pattern. A new control analysis result was developed there for the existence of an optimal control in the case of quadratic functions of the control appearing in both ODEs, under a restriction on the coefficients of the quadratics. A more detailed biological background on cottonwood - salt cedar competition and its management can be found in \cite{miller}. To understand the impact of habitat heterogeneity and other spatial features on the management of invasive-native species, we consider a model which includes explicit representation of space. The competition model as well as an optimal control problem considered in this paper are extensions of those in \cite{miller} in that respect that now they include space variable hence leading to a system of PDEs. PDEs are one of the standard modelling frameworks to account for spatial effects in ecological systems \cite{cc}. Optimal control of PDEs involves deriving an optimal spatial and temporal strategy, and the corresponding techniques are just beginning to be applied to natural resource problems with appropriate numerical algorithms for specific scenarios \cite{bhat, ding, fister, nt, sl}. The theoretical foundation of optimal control of models with underlying dynamics given by PDEs was developed by Lions \cite{lions} whereas to derive optimal control solutions for problems in which the state variable dynamics is governed by ODEs one needs to use Pontryagin's Maximum Principle \cite{pmp}. Even though there is no full generalization of Pontryagin's Maximum Principle to partial differential equations, the control theory and corresponding analysis for many special cases of PDEs have been developed \cite{at, barbu,fa,lt, sl2, sl3, sl4, yong}. In this paper we consider a model with diffusive and convective terms since many plants rely on animals/insects/humans and wind as their dispersal agents. A more detailed justification as well as description of models that account for different dispersal mechanisms for plants can be found in \cite{lewis, okubo1, okubo2, petrovskii, shigesada} and references therein. The paper is organized as follows: in section~\ref{model}, we introduce our system of partial differential equations and set up the optimal control problem. In section~\ref{existenceofthestatesystem} we prove the existence and uniqueness of solution to the state system. In section~\ref{existenceofanoptimalcontrol}, we present the existence of an optimal control in the case of quadratic functions of the control appearing in both PDEs. In sections \ref{section_sensitivity} and \ref{section_adjoint} we establish the adjoint system and characterize the optimal control. Section~\ref{uniqueness_of_OC} deals with the uniqueness of the optimal control. Finally, in section~\ref{numerical_results} we give some numerical examples to illustrate the results. \section{The model}\label{model} The state variables $N_1(x,t)$ and $N_2(x,t)$ represent the population densities (expressed as seedlings per unit of area) of a native and an invasive species, respectively. Define $Q = \Omega\times(0,T)$, where the domain $\Omega\in \mathbb{R}^n$ is bounded and smooth and $T$ denotes time. Let the PDE operators $\mathbf{L}_k, k=1,2$ be \begin{equation} \mathbf{L}_k N_k = \sum_{i,j=1}^n \Big( d_{ij}^k(x,t) (N_k)_{x_i} \Big)_{x_j} - \sum_{i=1}^n r_i^k(x,t) (N_k)_{x_i}; \ k=1,2. \end{equation} The general model for the population dynamics is given by \begin{equation} \label{eqnstate} \begin{gathered} (N_1)_t = \mathbf{L}_1 N_1+\Big(\theta_1 (t,u(x,t))-a_{11}N_1 \Big)N_1 -a_{12}N_1 N_2,\\ (N_2)_t = \mathbf{L}_2 N_2 + \Big(\theta_2 (t,u(x,t))-a_{22}N_2\Big)N_2 -a_{21}N_1 N_2, \quad \text{in } Q. \end{gathered} \end{equation} The diffusion coefficients are $d_{ij}^k(x,t)$, $i,j=1,\dots,n$, $k=1,2$, the convective coefficients are $r^k_i(x,t), i=1,\dots,n$, $k=1,2$, and the interaction coefficients are $a_{ij}\geq 0$, $i,j=1,2, $ indicating how species $j$ affects species $i$. The quadratic functions \begin{equation} \theta_i(t,u(x,t)) = \Big(a_i u^2(x,t) + b_i u(x,t) \Big)\chi_{\Gamma}+c_i, \quad i=1,2, \end{equation} are the intrinsic growth rates, with constants $a_i, b_i,c_i$, $(i=1,2)$ and the effect of disturbance in the system (such as flooding) is modelled by a control variable $u(x,t)$. Here $\chi_{\Gamma}$ is the characteristic function of the set $\Gamma$, which is given by $$ \Gamma = \cup_{i=1}^T \ [\sigma_i, \ \tau_i], $$ with the intervals $[\sigma_i, \tau_i]$ being in the $i$-th year. Thus, on the set $(0,T)\setminus \Gamma$, the control does not appear in the state system, which amounts to saying that the control does not apply during those time intervals. With our motivating example in mind, this means that only flooding during the spring thaw when cottonwoods release their seeds is allowed (see \cite{miller}). With flooding being detrimental at low and high levels and being advantageous at medium levels, we consider quadratic growth functions of the control. The control set is \begin{equation} U = \{ u\in L^\infty (\Omega\times\Gamma) : 0\leq u(x,t)\leq M \}, \end{equation} where the constant $M>0$. The growth rate parameters $a_i, b_i$, and $c_i$, are chosen to fit a specific scenario. In the cottonwood - salt cedar situation, the parameters are chosen so that with little or no flooding the salt cedar population $N_2(x,t)$ has a higher growth rate than the cottonwood population $N_1(x,t)$ (see \cite{miller}). The initial and boundary conditions are: \begin{equation}\label{stateboundary} \begin{gathered} N_1(x,0) = N_{10}(x), \quad N_2(x,0) = N_{20}(x), \quad x\in\Omega, \\ N_1(x,t) = N_2(x,t) = 0, \quad (x,t)\in\partial\Omega\times (0,T), \end{gathered} \end{equation} where the initial populations $N_{10}(x), N_{20}(x)$ are given. The goal is to maximize \begin{equation}\label{ob} J(u) = \int_\Omega \Big( A N_1(x,T) - B N_2(x,T) \Big) \,dx - \iint\limits_{\Omega\times \Gamma}\Big(B_1 u(x,t) + B_2 u^2(x,t)\Big)\,dx\, dt, \end{equation} where the constants $A,B\geq 0$, the cost coefficients are $B_1,B_2\geq 0$. This means we want to maximize the native population $N_1(x,t)$ and minimize the invasive population $N_2(x,t)$ at the final time $T$, while minimizing the cost of applying the control. The costs can include the actual financial impact to carry out the management plan but also may consider the negative impact on the environment or economic development. \section{Existence of the state system}\label{existenceofthestatesystem} The following assumptions are used this article: \begin{itemize} \item[(A1)] $\Omega$ is a smooth bounded domain in $\mathbb{R}^n$, \item[(A2)] $N_{10}(x), N_{20}(x)\in L^{\infty}(\Omega)$ and $N_{10}(x), N_{20}(x)\geq 0$ on $\Omega$, \item[(A3)] $r^k_i\in C^1(\bar Q)$ and $d_{ij}^k \in L^{\infty}(Q)$, $d_{ij}^k = d_{ji}^k$, for $ i,j=1,\dots,n$; $k=1,2$, \item[(A4)] $\sum_{i,j=1}^n d_{ij}^k(x,t)\xi_i\xi_j \geq \mu|\xi|^2$, for $k=1,2$, $\mu>0$, for all $(x,t)\in Q$, $\xi\in\mathbb{R}^n$. \item[(A5)] The underlying state space for system \eqref{eqnstate} is $V:=L^2(0,T;H_0^1(\Omega))$. \end{itemize} We define the following bilinear form for $v,\psi\in V$ and a.e. $0\leq t\leq T$, \begin{equation} B^k[v,\psi;t]= \int_{\Omega}\sum_{i,j=1}^n d_{ij}^k (x,t)v_{x_i}\psi_{x_j} - \sum_{i=1}^n (r_i^k(x,t)\psi)_{x_i} v \,dx, \ k=1,2. \end{equation} \begin{definition} \label{def3.1} \rm A pair of functions $u,v\in L^2(0,T;H_0^1(\Omega))$ is a weak solution of system \eqref{eqnstate} and \eqref{stateboundary} provided that: \begin{itemize} \item[(1)] $u_t, v_t\in L^2(0,T;H^{-1}(\Omega))$, \item[(2)] \begin{equation} \label{weaksoln} \begin{gathered} \int_0^T \{\langle(N_1)_t,\phi_1\rangle + B^1[N_1,\phi_1; t]\} = \int_{Q}\Big((\theta_1(t,u)-a_{11}N_1)N_1 - a_{12}N_1 N_2\Big)\phi_1 \\ \int_0^T \{\langle (N_2)_t,\phi_2\rangle + B^2[N_2,\phi_2; t]\} = \int_{Q}\Big((\theta_2(t,u)-a_{22}N_2)N_2 - a_{21}N_1 N_2\Big)\phi_2 \end{gathered} \end{equation} for all $\phi_1, \phi_2\in V$, and \item[(3)] $N_1(x,0)=N_{10}(x)$, $N_2(x,0)=N_{20}(x)$. \end{itemize} Here $\langle ,\rangle$ denotes the duality between $H_0^1(\Omega)$ and $H^{-1}(\Omega)$. Also, note that since $N_1, N_2\in C([0,T];L^2(\Omega))$ it follows that condition (3) makes sense \cite{evans}. \end{definition} \begin{theorem} \label{thm} Given $u\in U$, there exists a unique $(N_1,N_2)\in(V\times V)\cap L^{\infty}(Q)$ which is the solution of he system \eqref{eqnstate} and \eqref{stateboundary}. \end{theorem} \begin{proof} To obtain the existence of the solution to the system \eqref{eqnstate}-\eqref{stateboundary}, we construct two sequences by means of iteration. To obtain $L^{\infty}$ bounds, we need supersolutions and subsolutions for the $N_1,N_2$ iterates. Let $\bar{N_1}, \bar{N_2}$ be the solutions of the problem \begin{equation} \label{eqn-n0} \begin{gathered} (\bar{N_1})_t = \mathbf{L}_1 (\bar{N_1})+\theta_1 (t,u(x,t))\bar{N_1}, \quad \text{in } Q,\\ (\bar{N_2})_t = \mathbf{L}_2 (\bar{N_2})+\theta_2 (t,u(x,t))\bar{N_2}, \\ \bar{N_1}(x,0) = N_{10}(x), \quad \bar{N_2}(x,0) = N_{20}(x), \quad x\in\Omega, \\ \bar{N_1}(x,t) = \bar{N_2}(x,t) = 0, \quad (x,t)\in\partial\Omega\times (0,T). \end{gathered} \end{equation} By standard maximum principle arguments \cite{evans}, $\bar{N_1}, \bar{N_2}$ are bounded in $L^{\infty}(Q)$. Define $N_1^0 = \bar{N_1}$, and $N_2^0$ to be the solution of \begin{equation} \label{eqn-n20} (N_2^0)_t = \mathbf{L}_2 (N_2^0)+ (\theta_2 (t,u(x,t))-a_{22}\bar{N_2})N_2^0 -a_{21}N_1^0 N_2^0 \quad \text{in } Q, \end{equation} with $N_2^0(x,0)=N_{20}(x)$, $x\in\Omega$, and $ N_2^0(x,t) = 0$, $(x,t)\in\partial\Omega\times (0,T)$. By the weak maximum principle $0\leq N_2^0(x,t)\leq \bar{N_2}(x,t)$ on $\Omega$. For $k=1,2,3,\dots$, we define $N_1^k, N_2^k$ to be the solutions of the following problems, respectively: \begin{equation} \label{eqn-n1} \begin{gathered} (N_1^k)_t -\mathbf{L}_1 N_1^k + R N_1^k = (\theta_1 (t,u)-a_{11}N_1^{k-1})N_1^{k-1} -a_{12}N_1^{k-1} N_2^{k-1}+R N_1^{k-1}, \\ (N_2^k)_t -\mathbf{L}_2 N_2^k + R N_2^k = (\theta_2 (t,u)-a_{22}N_2^{k-1})N_2^{k-1} -a_{21}N_1^{k-1} N_2^{k-1}+R N_2^{k-1}, \end{gathered} \end{equation} with $N_1^k(x,0) = N_{10}(x)$, $N_2^k(x,0)= N_{20}(x)$, for $x\in\Omega$, and $N_1^k(x,t)=N_2^k(x,t) = 0$, on $\partial\Omega\times (0,T)$, and the constant $R$ is chosen to satisfy $R \geq 2a_{11}\bar{N_1}+a_{12}\bar{N_2}-\theta_1$ and $R \geq 2a_{22}\bar{N_2}+a_{21}\bar{N_1}-\theta_2$. Let $k=1$ in equation \eqref{eqn-n1}, subtracting equation \eqref{eqn-n1} from \eqref{eqn-n0}, we have \[ (N_1^0 - N_1^1)_t -\mathbf{L}_1 (N_1^0 -N_1^1) + R (N_1^0 - N_1^1) = a_{11}N_1^0 N_1^0 + a_{12}N_1^0 N_2^0 \geq 0. \] By comparison results, we have $N_1^0 \geq N_1^1$. Similarly, let $k=1$ in \eqref{eqn-n1}, subtracting equation \eqref{eqn-n20} and \eqref{eqn-n1}, we obtain \[ (N_2^0 - N_2^1)_t - \mathbf{L}_2(N_2^0-N_2^1)+ R (N_2^0 - N_2^1) = -a_{22}(\bar{N_2} - N_2^0) N_2^0 \leq 0. \] Since $\bar{N_2}\geq N_2^0$, then $N_2^0 \leq N_2^1$. Let $k=1$ and $k=2$ in \eqref{eqn-n1}. Then we have \begin{align*} &(N_1^1 - N_1^2)_t - \mathbf{L}_1(N_1^1-N_1^2) + R (N_1^1 - N_1^2) \\ &= \theta_1(t,u(x,t)) (N_1^0 - N_1^1) - a_{11}(N_1^0 N_1^0 - N_1^1 N_1^1) - a_{12}(N_1^0 N_2^0 - N_1^1 N_2^1) \\ &\quad+ R(N_1^0 - N_1^1)\\ &\geq \theta_1(t,u(x,t)) (N_1^0 - N_1^1) - a_{11}(N_1^0 - N_1^1)(N_1^0 + N_1^1) - a_{12}N_2^1 (N_1^0 - N_1^1)\\ &\quad + R(N_1^0 - N_1^1) \\ &\geq (N_1^0 - N_1^1)\Big( \theta_1(t,u(x,t)) - a_{11}(N_1^0 + N_1^1) -a_{12}N_2^1 +R\Big) \geq 0, \end{align*} where we used $R > 2a_{11}\bar{N_1}+a_{12}\bar{N_2}-\theta_1$, $N_1^1 \leq N_1^0$ and $N_2^1\geq N_2^0$. By comparison results, $N_1^1 \geq N_1^2$. Similarly, we obtain $N_2^1 \leq N_2^2$. By an induction argument coupled with comparison results, we can prove that there exist $N_1^k, N_2^k, \tilde{N_1}, \tilde{N_2}$, such that the following monotone pointwise convergence hold $N_1^k \searrow \tilde{N_1}$, $N_2^k \nearrow \tilde{N_2}$ in $Q$, and $ 0\leq N_1^k\leq \bar{N_1}$, $0\leq N_2^k \leq \bar{N_2}$ for $ k=1,2,\dots$. From the boundedness of $N_1^k$ and $N_2^k$ in $V$ and \eqref{eqn-n1} (see similar arguments in \cite{sl3}), we obtain that $\|(N_1)_t^k\|, \|(N_2)_t^k\|$ are bounded in $L^2(0,T;H^{-1}(\Omega))$. Hence, using weak compactness, we have $$ (N_1^k)_t \rightharpoonup \tilde{(N_1)}_t, \quad (N_2^k)_t \rightharpoonup \tilde{(N_2)}_t. $$ Since $L^2(0,T;H_0^1(\Omega))\subset\subset L^2(Q)$ \cite{simon}, we have $N_1^k \to \tilde{N_1}$, $N_2^k \to \tilde{N_2}$ strongly in $L^2(Q)$. Now we outline the proof of uniqueness. Let $(\check{N_1},\check{N_2})$ be another solution to the state system \eqref{eqnstate} with the initial and boundary conditions ${\check N_1}(x,0) = N_{10}(x)$, ${\check N_2}(x,0) = N_{20}(x), x\in\Omega$, ${\check N_1}(x,t) = {\check N_2}(x,t) = 0,(x,t)\in\partial\Omega\times (0,T)$. Let ${\check N_1}=e^{\lambda t}{\check w_1}$, ${\check N_2}=e^{\lambda t}{\check w_2}$, $N_1=e^{\lambda t}w_1$, and $N_2=e^{\lambda t}w_2$, where$\lambda>0$ is to be chosen later. After making this change of variables in the corresponding equations, we obtain \begin{align*} &\int_0^T \langle (w_1)_t,\phi_1\rangle\,dt +\int_Q\lambda w_1\phi_1\,dx\,dt \\ &+ \int_{Q}\sum_{i,j=1}^n d_{ij}^1 (x,t)(w_1)_{x_i}(\phi_1)_{x_j} - \int_Q\sum_{i=1}^n (r_i^1(x,t)\phi_1)_{x_i} w_1\\ &= \int_{Q}\Big((\theta_1(t,u)-e^{\lambda t}a_{11}w_1)w_1 - e^{\lambda t}a_{12}w_1 w_2\Big)\phi_1 , \end{align*} where $\phi_1, \phi_2\in V$. Similar equations can be written for $w_2, {\check w_1}$, and ${\check w_2}$. Define $Q_1:=\Omega\times(0,T_1)$ for $00$, taking into account \eqref{restriction}, \eqref{convexity3}, and \eqref{convexity4}, we get \begin{align*} &\lambda z^0+(1-\lambda)\zeta^0-f^0(t,w)\\ &\geq\lambda f^0(t,u)+(1-\lambda)f^0(t,v)-f^0(t,w)\\ &= B_2[\lambda u^2+(1-\lambda) v^2-w^2]+B_1[\lambda u+(1-\lambda)v-w]\\ &= \Big(B_1-\frac{B_2}{\kappa}\Big)[\lambda u+(1-\lambda)v-w]\geq 0. \end{align*} In the case when $\kappa=0$, \eqref{convexity4} becomes $w=\lambda u+(1-\lambda)v$, which leads to (due to convexity of the map $u\mapsto u^2$) \begin{align*} \lambda z^0+(1-\lambda)\zeta^0-f^0(t,w) &\geq\lambda f^0(t,u)+(1-\lambda)f^0(t,v)-f^0(t,w)\\ &= B_2\{\lambda u^2+(1-\lambda)v^2-[\lambda u+(1-\lambda)v]^2\}\geq 0. \end{align*} Hence, it follows that \[ \lambda(z^0,z)+(1-\lambda)(\zeta^0,\zeta)\in\mathcal{E}(t,N), \] which proves the convexity of the set $\mathcal{E}(t,N)$. \end{proof} We remark that the constraints on the parameters obtained in \eqref{restriction} are identical to those in \cite{miller}. \section{Sensitivity}\label{section_sensitivity} To characterize the optimal control, we need to differentiate the objective functional with respect to the control $u$. Since $N_i = N_i(u), i=1,2$ is involved in $J(u)$, we first must prove appropriate differentiability of the mapping $u\longrightarrow N_i(u), i=1,2$ whose derivative is called the \emph{sensitivity}. \begin{theorem} \label{theorem_sensitivity} The mapping $u\mapsto N(u)=(N_1(u),N_2(u))$ is differentiable in the following sense: \begin{equation} \label{psionetwo} \begin{gathered} \frac{N_1(u+\varepsilon h)-N_1(u)}{\varepsilon} \stackrel{w} \rightharpoonup \Psi_1\ {\rm in}\ L^2(0,T;H_0^{1}(\Omega)), \\ \frac{N_2(u+\varepsilon h)-N_2(u)}{\varepsilon} \stackrel{w} \rightharpoonup \Psi_2\ {\rm in}\ L^2(0,T;H_0^{1}(\Omega)), \end{gathered} \end{equation} as $\varepsilon\to 0$ for any $u\in U$ and $h\in L^\infty(Q)$ such that $u+\varepsilon h\in U$ for small $\varepsilon$. Moreover, the sensitivities, $\Psi_1\in L^2(0,T;H_0^{1}(\Omega))$ and $\Psi_2\in L^2(0,T;H_0^{1}(\Omega)))$ satisfy \begin{equation} \label{sensitivities} \begin{gathered} \begin{aligned} (\Psi_1)_t &= \mathbf{L}_1\Psi_1+ ((a_1 u^2 + b_1 u)\chi_{\Gamma}+c_1)\Psi_1 \\ &\quad + (2a_1 u +b_1)h N_1 -2a_{11} N_1\Psi_1 -a_{12}(N_1\Psi_2 + N_2\Psi_1), \end{aligned} \\ \begin{aligned} (\Psi_2)_t &= \mathbf{L}_2\Psi_2 + ((a_2 u^2 +b_2 u)\chi_{\Gamma}+c_2)\Psi_2 \\ &\quad + (2a_2 u +b_2)h N_2 -2a_{22} N_2\Psi_2 -a_{21}(N_1\Psi_2 + N_2\Psi_1). \end{aligned} \\ \Psi_1(x,0) = \Psi_2(x,0)=0, \quad x\in\Omega, \\ \Psi_1(x,t) = \Psi_2(x,t) = 0, \quad (x,t)\in\partial\Omega\times (0,T). \end{gathered} \end{equation} \end{theorem} \begin{proof} Define $N_1^\varepsilon = N_1(u+\varepsilon h)$, $N_2^\varepsilon = N_2(u+\varepsilon h)$. Let $N_1^\varepsilon=e^{\lambda t}W_1^\varepsilon$, $N_2^\varepsilon=e^{\lambda t}W_2^\varepsilon$, $N_1=e^{\lambda t}W_1$, and $N_2=e^{\lambda t}W_2$, where $\lambda>0$ is to be chosen later. After making this change of variables in \eqref{eqnstate}, we obtain \begin{align*} &\int_0^T \langle (W_1)_t,\phi_1\rangle\,dt +\int_Q\lambda W_1\phi_1 \\ &+ \int_{Q}\sum_{i,j=1}^n d_{ij}^1 (x,t)(W_1)_{x_i}(\phi_1)_{x_j} - \int_Q\sum_{i=1}^n (r_i^1(x,t)\phi_1)_{x_i} W_1\\ &= \int_{Q}\Big((\theta_1(t,u)-e^{\lambda t}a_{11}W_1)W_1 - e^{\lambda t}a_{12}W_1 W_2\Big)\phi_1, \end{align*} where $\phi_1\in V$. Similar equations can be written for $W_2$, $W_1^\varepsilon$ and $W_2^\varepsilon$. We derive our estimate on $Q_1:=\Omega\times(0,T_1)$. Then subtracting the equation for $W_1$ from the equation for $W_1^\varepsilon$ and dividing by $\varepsilon$, and taking $\phi_1=\frac{W_1^\varepsilon-W_1}{\varepsilon}$, using uniform ellipticity, regularity of the coefficients, and Cauchy inequality where necessary, we obtain the estimate \begin{align*} &\frac{1}{2}\int_{\Omega\times\{t=T_1\}} \Big|\frac{W_1^\varepsilon-W_1}{\varepsilon}\Big|^2\,dx + \lambda\int_{Q_1} \Big|\frac{W_1^\varepsilon-W_1}{\varepsilon}\Big|^2+\frac{\mu}{2}\int_{Q_1} \Big|\nabla\Big(\frac{W_1^\varepsilon-W_1}{\varepsilon}\Big)\Big|^2\\ &\leq C_1\int_{Q_1} \Big|\frac{W_1^\varepsilon-W_1}{\varepsilon}\Big|^2\,dx\,dt +C_2e^{\lambda T_1}\int_{Q_1}\Big|\frac{W_1^\varepsilon-W_1}{\varepsilon}\Big|^2|W_1+W_1^\varepsilon|\\ &\quad +C_3\int_{Q_1}|W_1^\varepsilon|\Big|\frac{W_1^\varepsilon-W_1}{\varepsilon}\Big|\,dx\,dt+ C_4e^{\lambda T_1}\int_{Q_1}|W_2|\Big|\frac{W_1^\varepsilon-W_1}{\varepsilon}\Big|^2\,dx\,dt \\ &\quad +C_5e^{\lambda T_1}\int_{Q_1}|W_1^{\varepsilon}|\Big|\frac{W_2^\varepsilon-W_2}{\varepsilon}\Big|\Big|\frac{W_1^\varepsilon-W_1}{\varepsilon}\Big|\,dx\,dt \\ &\leq C_6 e^{\lambda T_1}\int_{Q_1} \Big|\frac{W_1^\varepsilon-W_1}{\varepsilon}\Big|^2\,dx\,dt+C_7\int_{Q_1}|h|^2\,dx\,dt. \end{align*} A similar estimate can be derived for the quotient $\frac{W_2^\varepsilon-W_2}{\varepsilon}$. Combining the estimates for both quotients gives \begin{align*} &\frac{1}{2}\int_{\Omega\times\{t=T_1\}} \Big|\frac{W_1^\varepsilon-W_1}{\varepsilon}\Big|^2+\Big|\frac{W_2^\varepsilon-W_2}{\varepsilon}\Big|^2\,dx\\ & + (\lambda-C_9 e^{\lambda T_1}-C_8)\int_{Q_1} \Big|\frac{W_1^\varepsilon-W_1}{\varepsilon}\Big|^2+\Big|\frac{W_2^\varepsilon-W_2}{\varepsilon}\Big|^2\,dx\,dt\\ &+\frac{\mu}{2}\int_{Q_1} \Big|\nabla\Big(\frac{W_1^\varepsilon-W_1}{\varepsilon}\Big)\Big|^2 +\frac{\mu}{2}\int_{Q_1} \Big|\nabla\Big(\frac{W_2^\varepsilon-W_2}{\varepsilon}\Big)\Big|^2\\ &\leq C_{10}\int_{Q_1}|h|^2\,dx\,dt, \end{align*} where $C_i$ only depends on the coefficients and $L^\infty$ bounds of $W_i, W_i^\varepsilon$. Choosing $\lambda$ large enough and $T_1$ small, we get the estimate \begin{align*} \|\frac{W_i^\varepsilon-W_i}{\varepsilon}\|_{L^2(0,T_1;H_0^1(\Omega))}\leq C, \quad i=1,2, \end{align*} where $C>0$ denotes a generic constant. By stacking time intervals, we can obtain the estimate on $\Omega\times(0,T)$. This will give us the desired estimate for $\frac{N_i^\varepsilon-N_i}{\varepsilon}$. Hence, we conclude that there exist $\Psi_i\in V, i=1,2$, such that on a subsequence $\frac{N_i^\varepsilon-N_i}{\varepsilon}\rightharpoonup \Psi_i$ in $L^2(0,T;H_0^1(\Omega))$. Using the state system \eqref{eqnstate} and the above estimate it can be shown that \begin{align*} \big\|\Big(\frac{N_i^\varepsilon-N_i}{\varepsilon}\Big)_t \big\|_{L^2(0,T;H^{-1}(\Omega))}\leq C, \end{align*} which implies that $(\frac{N_i^\varepsilon-N_i}{\varepsilon})_t \stackrel{w}\rightharpoonup (\Psi_i)_t$ in $L^2(0,T;H^{-1}(\Omega))$. Using the compactness result by Simon \cite{simon}, we have $\frac{N_i^\varepsilon-N_i}{\varepsilon}\to \Psi_i$ strongly in $L^2(Q)$. Now, passing to the limit in the weak form of the PDEs satisfied by the quotients $\frac{N_i^\varepsilon-N_i}{\varepsilon}, \ i=1,2$, we obtain \eqref{sensitivities}. \end{proof} To rewrite the sensitivity system in matrix form, let \begin{equation} M= \begin{pmatrix} M_{11} & M_{12} \\ M_{21} &M_{22}, \end{pmatrix}, \end{equation} where $M_{11}=-((a_1 u^2 + b_1 u)\chi_{\Gamma}+c_1)+2a_{11} N_1 +a_{12}N_2$, $M_{12}=a_{12}N_1$, $M_{21}=a_{21}N_2$, and $M_{22}=-((a_2 u^2 + b_2 u)\chi_{\Gamma}+c_2)+2a_{22} N_2 +a_{21}N_1$. Then \begin{equation}\label{sensitivity} \mathbf{L} \begin{pmatrix} \Psi_1 \\ \Psi_2 \end{pmatrix} = \begin{pmatrix} (\Psi_1)_t -\mathbf{L}_1\Psi_1 \\ (\Psi_2)_t - \mathbf{L}_2\Psi_2 \end{pmatrix} +M \begin{pmatrix} \Psi_1 \\ \Psi_2 \end{pmatrix}. \end{equation} \section{The adjoint system}\label{section_adjoint} Now we are ready to characterize the optimal control, by deriving the optimality system through differentiating $J(u)$ with respect to $u$ at an optimal control. Define the adjoint system as \begin{equation}\label{adjoint} \mathbf{L}^* \begin{pmatrix} p \\ q \end{pmatrix} = \begin{pmatrix} A \\ -B \end{pmatrix}, \end{equation} where \begin{equation} \mathbf{L}^* \begin{pmatrix} p \\ q \end{pmatrix} =\begin{pmatrix} -p_t -\sum_{i,j=1}^n \Big( d_{ij}^1(x,t) p_{x_i} \Big)_{x_j} - \sum_{i=1}^n r_i^1(x,t) p_{x_i} \\ -q_t - \sum_{i,j=1}^n \Big( d_{ij}^2(x,t) q_{x_i} \Big )_{x_j} -\sum_{i=1}^n r_i^2(x,t) q_{x_i} \end{pmatrix} +M^{T} \begin{pmatrix} p \\ q \end{pmatrix}, \end{equation} where $M^T$ is the transpose of the matrix $M$. The terminal conditions are \begin{equation}\label{adjoint_terminal_condit} p(x,T)=A, \quad q(x,T)=-B, \quad x\in\Omega. \end{equation} So given an optimal control $u^*$ and the corresponding states $N_1^*$ and $N_2^*$, the adjoint variables $p$ and $q$ satisfy \begin{gather*}%\label{adjoint1} \mathbf{L}_1^* p = \Big( (a_1 (u^*)^2 + b_1 u^*)\chi_{\Gamma}+c_1 - 2a_{11}N_1^* -a_{12} N_2^*\Big) p - a_{21}N_2^* q, \quad\text{in } \Omega, \\ p = 0, \text{ on } \partial\Omega\times(0,T), \quad p(x,T)= A, \text{ for } x\in\Omega, \\ \mathbf{L}_2^* q = \Big( (a_2 (u^*)^2 + b_2 u^*)\chi_{\Gamma}+c_2 - 2a_{22}N_2^* -a_{21} N_1^*\Big) q - a_{12}N_1^* p, \quad\text{in } \Omega, \\ q = 0, \text{ on } \partial\Omega\times(0,T), \quad q(x,T)= -B, \text{ for } x\in\Omega. \end{gather*} \begin{theorem} \label{theorem_characterization} If there exists $\kappa\geq 0$, such that $B_1 \kappa -B_2<0$, $a_i=\kappa b_i$, $i=1,2$, then given an optimal control $u^*$ and the corresponding states $N_1^*$ and $N_2^*$, there exist solutions $p$ and $q$ to the adjoint system. Moreover, let $S=\{(x,t) : B_2 - a_1 N_1^* p -a_2 N_2^* q =0\}$ and $m(S)$ be a Lebesgue measure of $S$, then this optimal control $u^*$ is characterized by the following: \begin{itemize} \item[(1)] if $m(S)>0$, then $u^*=M$ on $S$; \item[(2)] if $m(S)=0$, then for $(x,t)\not\in S$, \begin{equation}\label{oc} u^* = \min\big\{M,\max\{0,\frac{b_1 N_1^* p + b_2 N_2^* q -B_1}{2B_2 - 2a_1 N_1^* p -2a_2 N_2^* q}\}\big\}, \end{equation} and it holds on $\Omega\times\Gamma$ a.e. \end{itemize} \end{theorem} \begin{proof} Since the adjoint system is linear, by the standard theory of parabolic equations there exists a weak solution $(p, q)$ satisfying the adjoint system. Let $N_i^\varepsilon = N_i(u^*+\varepsilon h)$, $i=1,2$, with $u^*+\varepsilon h\in U$ and $h\in L^\infty(Q)$. Then \begin{equation} \label{oc2} \begin{aligned} 0&\geq \lim_{\varepsilon\to 0^+}\frac{J(u^*+\varepsilon h) - J(u^*)}{\varepsilon} \\ &= \lim_{\varepsilon\to 0^+}\int_\Omega \Big(A\frac{N_1^\varepsilon-N_1}{\varepsilon}(x,T) -B \frac{N_2^\varepsilon-N_2}{\varepsilon}(x,T)\Big) - \iint_{\Omega\times\Gamma} (B_1 h + 2B_2 u^* h) \\ &= \int_\Omega \Big( A\Psi_1(x,T)-B\Psi_2(x,T)\Big) \,dx - \iint_{\Omega\times\Gamma}(B_1 +2B_2 u^*)h \,dx\,dt \\ &= \int_{\Omega\times(0,T)}(p, \ q)\mathbf{L} \begin{pmatrix} \Psi_1 \\ \Psi_2 \end{pmatrix} \,dx - \iint_{\Omega\times\Gamma}(B_1 +2B_2u^*)h \,dx\,dt \\ &= \iint_{\Omega\times\Gamma}(p, \ q) \begin{pmatrix} (2a_1u^*+b_1)h N_1^* \\ (2a_2 u^*+b_2)h N_2^* \end{pmatrix}\,dx\,dt -\iint_{\Omega\times\Gamma}(B_1 +2B_2u^*)h \,dx\,dt \\ &= \iint_{\Omega\times\Gamma}\Big((2a_1u^*+b_1)N_1^* p +(2a_2 u^*+b_2) N_2^* q -B_1 -2B_2 u^*\Big)h \,dx\,dt \\ &= \iint_{\Omega\times\Gamma}\Big(u^*(2 a_1 N_1^* p+2 a_2 N_2^* q-2 B_2)+ b_1 N_1^* p + b_2 N_2^* q - B_1 \Big)h \,dx\,dt. \end{aligned} \end{equation} Observe that we used \eqref{adjoint} so that \begin{equation} \iint_{\Omega\times(0,T)}(\Psi_1, \Psi_2)\mathbf{L}^* \begin{pmatrix} p \\ q \end{pmatrix} =\iint_{\Omega\times(0,T)}(p, q) \mathbf{L} \begin{pmatrix} \Psi_1 \\ \Psi_2 \end{pmatrix} \end{equation} in weak sense. Also we used the sensitivity system \eqref{sensitivity}. Define $S=\{(x,t): B_2 - a_1 N_1^* p -a_2 N_2^* q =0\}$. Let $m(S)>0$ and $h$ have support on $S$. Recalling that $a_i = \kappa b_i$, $i=1,2$, it follows from \eqref{oc2} that \begin{align*} 0 &\geq \iint_S \kappa(b_1 N_1^*p + b_2 N_2^*q - B_1)h \,dx\,dt\\ &=\iint_S (\kappa b_1 N_1^* p + \kappa b_2 N_2^* q -\kappa B_1)h \,dx\,dt \\ &= \iint_S (a_1 N_1^* p + a_2 N_2^* q - \kappa B_1)h \,dx\,dt = \iint_S (B_2 - \kappa B_1 )h \,dx\,dt, \end{align*} if $\kappa B_1 < B_2$, then $h\leq 0$, which implies $u^*=M$ on $S$. Otherwise, if $m(S)=0$, then for $(x,t)\not\in S$, on the set $00$. Thus, the original nonlinear system of parabolic PDE becomes: \begin{equation} \label{eqn1} \begin{gathered} (N_1)_t = d^1 (N_1)_{xx} - r^1 (N_1)_x + \Big(\theta_1(t, u(x,t)) -a_{11}N_1\Big)N_1 - a_{12}N_1 N_2, \\ (N_2)_t = d^2 (N_2)_{xx} - r^2 (N_2)_x + \Big(\theta_2(t,u(x,t)) -a_{22}N_2\Big)N_2 - a_{21}N_1 N_2, \end{gathered} \end{equation} where $$ \theta_i(t, u(x,t)) = \Big(a_i u(x,t)^2+b_i u(x,t)\Big)\chi_{\Gamma}+c_i, \quad i=1,2. $$ Here $d^i, r^i, a_{ij}$, $i,j=1,2$ are constants. Flooding is allowed during spring thaw when cottonwood seeds are dispersed from trees. The initial and boundary conditions are: \begin{equation} \label{bdry1} \begin{gathered} N_1(x,0) = N_{10}(x), \quad N_2(x,0) = N_{20}(x), \quad x\in[0,L], \\ N_1(0,t) = N_2(0,t) = 0, \quad N_1(L,t) = N_2(L,t) = 0, \quad t\in(0,T). \end{gathered} \end{equation} The adjoint equations are \begin{equation} \label{eqn2} \begin{gathered} -p_t = d^1 p_{xx}+r^1 p_x +\Big(\theta_1(t,u(x,t))-2a_{11}N_1-a_{12}N_2 \Big)p -a_{21}N_2 q, \\ -q_t = d^2 q_{xx}+r^2 q_x +\Big(\theta_2(t,u(x,t))-2a_{22}N_2-a_{21}N_1 \Big)q -a_{12}N_1 p, \end{gathered} \end{equation} with the terminal and boundary conditions \begin{equation} \label{bdry2} \begin{gathered} p(x,T) = A, \quad q(x,T) = -B, \quad x\in[0,L], \\ p(0,t) = q(0,t) = 0, \quad p(L,t) = q(L,t) = 0, \quad t\in(0,T). \end{gathered} \end{equation} The characterization for the optimal control is \begin{equation}\label{oc1} u^* = \min\big\{M,\max\{0,\frac{b_1 N_1^* p + b_2 N_2^* q -B_1}{2B_2 - 2a_1 N_1^* p -2a_2 N_2^* q}\}\big\}. \end{equation} Note that in our numerical simulations, we always have $m(S)=0$, where $S$ is defined in Theorem~\ref{theorem_characterization}. We will solve equations \eqref{eqn1}--\eqref{oc1}. In the cottonwood - salt cedar scenario, the quadratic growth the functions $\theta_i (t, u(x,t))$, $i=1,2$, represent reasonable qualitative behavior of the newly developed seedling communities. To model this interaction, we choose $\theta_1 < \theta_2$ when $u=0$ and \[ \theta_1 = 0.1 u^2+ u + 0.2, \quad \theta_2 = -0.1u^2 - u +0.4, \] where $\kappa = 0.1$ which satisfies the restriction in \eqref{restriction}. \begin{figure}[ht] \begin{center} \includegraphics[width=0.45\textwidth]{fig1a} \quad %./n1-noc.eps \includegraphics[width=0.45\textwidth]{fig1b} \\ % ./n2-noc.eps \text{Cottonwood } \hfil \text{Salt cedar} \end{center} \caption{Cottonwood and salt cedar without control, $L=1$, $T=3$} \label{fig1} \end{figure} \begin{figure}[ht] \begin{center} \includegraphics[width=0.45\textwidth]{fig2a} \quad % ./n1-B2-5.eps \includegraphics[width=0.45\textwidth]{fig2b} \\ %./n2-B2-5.eps \text{Cottonwood} \hfil \text{Salt cedar} \end{center} \caption{Cottonwood and salt cedar with $B_2=5$, $B_1=1$} \label{fig2} \end{figure} \begin{figure}[ht] \begin{center} \includegraphics[width=0.45\textwidth]{fig3a} \quad % ./n1-B2-10.eps} \includegraphics[width=0.45\textwidth]{fig3b} \\ %./n2-B2-10.eps \text{Cottonwood} \hfil \text{Salt cedar} \end{center} \caption{Cottonwood and salt cedar with $B_2=10$, $B_1=1$} \label{fig3} \end{figure} Since the salt cedar currently takes over the cottonwood, our initial conditions are $N_1(x,0) = x(1-x)$, $N_2(x,0)=3x(1-x)$. The upper bound for the control is $M=1$, the length of the one dimensional interval is $L=1$. We want to see the effect of the growth functions $\theta_i$, so we choose the diffusion and convection coefficients $d^k =1$, $r^k = 1$, $k=1,2$. The interaction coefficients $a_{ij}$ represent how species $j$ affects species $i$ and we choose $a_{11} = 0.005$, $a_{12} = 0.001$, $a_{21} = 0.12$, and $ a_{22} = 0.02$, since Sher et al.\ \cite{sher3} found that cottonwood density affects the population of both species (in a negative way) more than the salt cedar density. Without any flooding (control), cottonwood density increases slowly and the salt cedar density increases rapidly. See Figure \ref{fig1}. Consider the case when $A=10, B=1$. First we pick $B_1=1, B_2=5$. To correspond to the spring thaw, flooding is only allowed during the first $1/4$ of the year in order to illustrate the idea for a period of $3$ years. We can see the cottonwood density increases while the salt cedar decreases and the time when this happens corresponds to the flooding time every year. See Figure \ref{fig2}. Note that we need to choose parameters $B_1, B_2$ carefully, because of the restriction in \eqref{restriction}. We want to see how the parameter $B_2$, the quadratic cost coefficient in $J(u)$, affects the optimal flooding strategy. We fix $B_1=1$, as we increase the values of $B_2$, flooding levels are decreasing, since it is more expensive to apply the flooding. See Figures \ref{fig2} and \ref{fig3}. We can see in Figure \ref{fig3} cottonwood density is still increasing and salt cedar is decreasing but not as much as in Figure \ref{fig2} when it is much cheaper to apply the control. Figure \ref{fig4} gives the corresponding optimal flooding strategies for $B_1=1$ and $B_2=5, 10$, we can examine that as $B_2$ increases, optimal control strategy is decreasing. \begin{figure}[ht] \begin{center} \includegraphics[width=0.45\textwidth]{fig4a} \quad % ./u-B2-5.eps \includegraphics[width=0.45\textwidth]{fig4b} %./u-B2-10.eps \text{Control: $B_2=5$}\hfil \text{Control: $B_2=10$} \end{center} \caption{Control - flood: $B_2 = 5, 10$, fixed $B_1=1$} \label{fig4} \end{figure} \begin{figure}[ht] \begin{center} \includegraphics[width=0.45\textwidth]{fig5} % ./u-B1-5-B2-5.eps \end{center} \caption{Control - flood: $B_2=5$, $B_1=5$} \label{fig5} \end{figure} We also want to study the effect of the parameter $B_1$ - the linear cost coefficient in $J(u)$. We fix $B_2=5$ and compare $B_1=1$ with $B_1=5$ to observe as $B_1 $ increases, the flooding is decreasing. See Figures \ref{fig4}(a) and \ref{fig5} for the comparison of the optimal flooding. In summary, the control (flooding) is focused near the middle of the region due to the Dirichlet boundary conditions we imposed on our problem. When increasing the cost coefficients $B_1$ or $B_2$, the flooding level is decreasing. Compared with Figure 1 without any control, Figures \ref{fig2} and \ref{fig3} gave good illustrations that flooding can indeed control the invasive species $N_2$ and help the growth of the native species $N_1$. \subsection*{Conclusions}%\label{conclusions} This project shows that the optimal control theory can be an appropriate tool for designing the intervention strategy of the invasive - native species interaction. We proved existence of the optimal control when the control is quadratic in the growth function in the PDE system under certain conditions on the coefficients. We gave numerical examples for different parameter values that can help natural resource managers to apply the most appropriate and cost-effective control methods to the invasive - native species scenario. \subsection*{Acknowledgements} W. Ding was supported by Faculty Research and Creative Activities (FRCAC) at Middle Tennessee State University. \begin{thebibliography}{00} \bibitem{at} N. U. Ahmed, K. L. Teo; \emph{Optimal control of distributed parameter systems}, North Holland, Amsterdam, 1981. \bibitem{barbu} V. Barbu; \emph{Analysis and control of nonlinear infinite dimensional systems}, Academic Press, New York, 1993. \bibitem{bhat} M. G. Bhat, K. R. Fister, S. Lenhart; \emph{An optimal control model for surface runoff contamination of a large river basin}, Natural Resource Modeling, 2 (1999) 175-195. \bibitem{cc} R. S. Cantrell, C. Cosner; \emph{Spatial ecology via reaction-diffusion equations}, Wiley, New Jersey, 2003. \bibitem{ding} W. Ding, S. Lenhart; \emph{Optimal Harvesting of a Spatially Explicit Fishery Model}, Natural Resource Modeling, 22:2 (2009) 173-211. \bibitem{evans} L. C. Evans; \emph{Partial Differential Equations}, American Mathematical Society, Providence, 1998. \bibitem{fa} H. O. Fattorini; \emph{Infinite dimensional optimization and control theory}, Cambridge, New York, 1999. \bibitem{fister} R. Fister; \emph{Optimal control of harvesting in a predator-pray parabolic system}, Houston J. of Math, 3 (1997) 341-355. \bibitem{miller} D. Kern, S. Lenhart, R. Miller, J. Yong; \emph{Optimal control applied to native-invasive population dynamics}, Journal of Biological Dynamics, 1:4 (2007) 413-426. \bibitem{lt} I. Lasiecka, R. Triggiani; \emph{Control theory for partial differential equations}, Vol. I and II, Cambridge Univ. Press, New York, 2000. \bibitem{sl} S. Lenhart, M. Bhat; \emph{Application of distributed parameter control model in wildlife damage management}, Mathematical Models and Methods in Applied Sciences, 2 (1993) 423-439. \bibitem{sl2} S. Lenhart, H. Thieme; \emph{Editorial note on special issue: Modeling and control of natural resource}s, Natural Resource Modeling, 18 (2005) 3:235-236. \bibitem{sl3} S. Lenhart, M. Liang, V. Protopopescu; \emph{Optimal control of boundary habitat hostility for interacting species}, Math. Methods Appl. Sci., 22:13 (1999) 1061-1077. \bibitem{sl4} S. Lenhart, J. Workman; \emph{Optimal Control Applied to Biological Models}, Chapman Hall/CRC, Boca Raton, 2007. \bibitem{lewis} M. A. Lewis, J. Bulllock; \emph{Mathematical models for plant dispersal}, Demography, (2003) 21-25. \bibitem{yong} X. Li, J. Yong; \emph{Optimal control theory for infinite dimensional systems}, Birkhauser, 1995. \bibitem{lions} J. L. Lions; \emph{Optimal control of systems governed by partial differential equations}, Springer-Verlag, New York, 1971. \bibitem{nt} P. Neittaanmaki, D. Tiba; \emph{Optimal control of nonlinear parabolic systems: theory, algorithms and applications}, Marcel Dekker, New York, 1994. \bibitem{nentwig} W. Nentwig (ed.); \emph{Biological Invasions}, Ecological Studies, 193 (2007), Springer. \bibitem{okubo1} A. Okubo, S. Levin; \emph{A theoretical framework for data analysis of wind dispersal of seeds and pollen}, Ecology, 70:2 (1989) 329-338. \bibitem{okubo2} A. Okubo, S. Levin; \emph{Diffusion and Ecological Problems: Modern Perspectives}, Springer, Berlin, 2001. \bibitem{petrovskii} S. V. Petrovskii, B. Li; \emph{Exactly solvable models of biological invasion}, Chapman \& Hall/CRC, 2006 \bibitem{pmp} L. S. Pontryagin, V. S. Boltyanskii, R. V. Gamkrelidze, E. F. Mishchenko; \emph{The Mathematical Theory of Optimal Processes}, Wiley-Interscience, New York, 1962. \bibitem{sher3} A. A. Sher, D. L. Marshall, J. P. Taylor; \emph{Establishment patterns of native Populus and Salix in the presence of invasive Tamarix}, Ecological Applications, 12 (2002) 760-772. \bibitem{shigesada} N. Shigesada, K. Kawasaki; \emph{Biological Invasions: Theory and Practice}, Oxford University Press, Oxford, 1997. \bibitem{simon} J. Simon; \emph{Compact sets in the space $L^p((0,T);B)$}, Ann. Mat. Pura. Appl. {CXLVI}, (1987) 65-96. \end{thebibliography} \end{document}