Optimal resource control in reaction diffusion advection population model

TL;DR

This paper studies optimal resource control in a reaction-diffusion-advection population model, establishing existence, uniqueness, and characterization of optimal controls, supported by numerical simulations under various boundary conditions.

Contribution

It introduces a new optimal control framework for a reaction-diffusion-advection population model with resource management, including existence and uniqueness results.

Findings

Optimal control exists and is unique.

Characterization of the optimal control is provided.

Numerical simulations demonstrate control strategies under different boundary conditions.

Abstract

This paper is to investigate the control problem of maximizing the net benefit of a single species while the cost of the resource allocation is minimized in a population model which can be described by a reaction diffusion advection equation of logistic type with spatial-temporal resource control coefficient. The existence of an optimal control is established and the uniqueness and characterization of the optimal control are investigated. Numerical simulation illustrate several cases with Dirichlet and Neumann boundary conditions.

Equations

$$\{\begin{array}{l}{u}_t-\nabla \cdot [\mu \nabla u-u\overrightarrow{h}]=u[m-f(x,t,u)],Q_T,\\ \mu \frac{\partial u}{\partial \nu }-u\overrightarrow{h}\cdot \nu =0,S_T,\\ u(\cdot ,0)=u_0\ge 0,\Omega ,\end{array}$$

$$U=\{m\in L^{\infty }(Q_T):0\le m\le M\},$$

$$J(m)={\int }_{Q_T}[u-(B{m}^2)]dxdt,$$

$$u_t=\lambda \Delta u+m(x)u-u^2,\Omega ,$$

$$u_t=\nabla \cdot [D\nabla u-\alpha u\nabla m(x)]+m(x)u,\Omega \times (0,\infty )$$

$$u_t=\nabla \cdot [D\nabla u-\alpha u\nabla m(x)]+m(x)u-c{u}^2,\Omega \times (0,\infty )$$

$$\{\begin{array}{ll}-p_t-\mu \Delta p-\vec{h}\cdot \nabla p-[m^{\ast }-2u^{\ast }]p=1,& \text{mbox}in{Q}_T\\ \frac{\partial p}{\partial \nu }=0,& \text{mbox}in\partial \Omega \times (0,T),\\ p(\cdot ,T)=0,& \text{mbox}in\Omega .\end{array}$$

$$m^{\ast }=\max \{\min \{M,\frac{u^{\ast }p}{2B}\},0\}.$$

$$\Vert u{\Vert }_{V_2(Q_T)}=\{\underset{0\le t\le T}{\sup }{\int }_{\Omega }u(x,t{)}^{2}dx+{\int }_{Q_T}|\nabla u(x,t){|}^{2}dxdt{\}}^{1/2}<\infty .$$

$${\int }_{\Omega }{u}_t\theta dx+{\int }_{\Omega }[\mu \nabla u-\overrightarrow{h}u]\cdot \nabla \theta dx={\int }_{\Omega }u[m(x,t)-u]\theta dx,\forall \theta \in H^1(\Omega ).$$

$${\int }_{\Omega }u(x,t)dx\le e^{t{C}_m}{\int }_{\Omega }{u}_0(x)dx,\forall t\in [0,T),$$

$$J(m)\le \frac{(e^{T{C}_m}-1)\Vert u_0{\Vert }_{L^{\infty }}}{C_m},\forall m\in U.$$

$$\Vert u(m){\Vert }_{V_2(Q_T)}\le C,\text{and}0\le u(m)\le C,\forall (x,t)\in Q_T.$$

$$\Vert {\phi }^{\epsilon }{\Vert }_{V_2(Q_T)},\Vert {\phi }^{\epsilon }{\Vert }_{L^{\infty }(Q_T)}\le C.$$

$${\phi }^{\epsilon }\rightharpoonup \phi \text{weakly in}{L}^2((0,T),H^1(\Omega ))\text{as}\epsilon \to 0,$$

$$\{\begin{array}{ll}{\phi }_t-\nabla \cdot [\mu \nabla \phi -\phi \vec{h}]-\phi [m-2u]& =-ul,\\ \mu \frac{\partial \phi }{\partial \nu }-\phi \vec{h}\cdot \nu & =0,\\ u(\cdot ,0)& =0.\end{array}$$

$$\{\begin{array}{ll}{u}_t^{\epsilon }-\nabla \cdot [\mu \nabla u^{\epsilon }-u^{\epsilon }\vec{h}]& =u^{\epsilon }[m+\epsilon l-u^{\epsilon }],\\ \mu \frac{\partial u^{\epsilon }}{\partial v}-u^{\epsilon }\vec{h}\cdot v& =0,\\ u(\cdot ,0)& =u_0.\end{array}$$

$$\Vert u^{\epsilon }{\Vert }_{V_2(Q_T)},\Vert u{\Vert }_{V_2(Q_T)},\Vert u^{\epsilon }{\Vert }_{L_{\infty }(Q_T)},\Vert u{\Vert }_{L_{\infty }(Q_T)}\le C<\infty ,\forall \epsilon >0,$$

$${\gamma }^{\epsilon }(x,t)=\frac{(u^{\epsilon }{)}^2(x,t)-u^2(x,t)}{u^{\epsilon }(x,t)-u(x,t)}=u^{\epsilon }(x,t)+u(x,t),\forall (x,t)\in Q_T.$$

$$u^{\epsilon }\to u\text{mbox}in{L}^2(Q_T);u_t^{\epsilon }\rightharpoonup u_t\text{mbox}in{L}^2(Q_T);u^{\epsilon }\rightharpoonup u\text{mbox}in{L}^2((0,T),H^1(\Omega )).$$

$$u^{\epsilon }\to ua.e.\text{in}{Q}_T,\text{mbox}and{\gamma }^{\epsilon }(x,t)\to 2u(x,t),a.e.\text{mbox}in{Q}_T.$$

$$\{\begin{array}{ll}{\phi }_t^{\epsilon }-\nabla \cdot [\mu \nabla {\phi }^{\epsilon }-{\phi }^{\epsilon }\vec{h}]-{\phi }^{\epsilon }[m-{\gamma }^{\epsilon }]& =u^{\epsilon }l,\\ \mu \frac{\partial {\phi }^{\epsilon }}{\partial \nu }-{\phi }^{\epsilon }\vec{h}\cdot \nu & =0,\\ {\phi }^{\epsilon }(\cdot ,0)& =0.\end{array}$$

$$\Vert {\phi }^{\epsilon }{\Vert }_{V_2(Q_T)}\le C,\Vert {\phi }^{\epsilon }{\Vert }_{L^{\infty }}\le C,\forall \epsilon .$$

$$0\le p(x,t)\le C\text{mbox}and|\nabla p(x,t)|\le C,\forall (x,t)\in Q_T.$$

$$\{\begin{array}{ll}& q_t-\mu \Delta q+b\cdot \nabla q+cq=1,\\ & \nabla q\cdot \nu =0,\\ & q(\cdot ,0)=0,\end{array}$$

$$\Vert p{\Vert }_{W_l^{2,1}(Q_T)}=\Vert q{\Vert }_{W_l^{2,1}(Q_T)}\le C=C(l),\forall l\ge 2,$$

$$\Vert p{\Vert }_{W_l^{2,1}(Q_T)}:=\Vert p{\Vert }_{L^l(Q_T)}+\Vert p_t{\Vert }_{L^l(Q_T)}+\Vert \nabla p{\Vert }_{L^l(Q_T)}+\sum _{i,j=1}^d\Vert p_{x_i{x}_j}{\Vert }_{L^l(Q_T)}.$$

$$|p(x,t)|,|\nabla p(x,t)|\le C(l),\forall (x,t)\in Q_T,$$

$$\underset{n\to \infty }{\lim }J(m^n)=\underset{m\in U}{\sup }J(m).$$

$$\Vert u^n{\Vert }_{V_2(Q_T)},\Vert u^n{\Vert }_{L^{\infty }(Q_T)}\le C<\infty ,\forall n\in \mathbb{N},$$

Full text

Key words. Optimal Control, Reaction-Diffusion-Advection, Population Models.

1 Introduction

In the current paper, we investigated the optimal resource control problem for the following nonlinear diffusion-reaction-advection equation:

[TABLE]

where $Q_{T}=\Omega\times(0,T)$ for $T>0,S_{T}=\partial\Omega\times(0,T)$. Here $\Omega$ is an open bounded domain in $d-$dimensional space $\mathbb{R}^{d},d\in\mathbb{N}$, with smooth boundary $\partial\Omega$, $m(x,t)$ is the controlled resource and $u(x,t)$ is the density. Moreover, $\mu>0$ is the diffusion rate, $\overrightarrow{h}:Q_{T}\rightarrow\mathbb{R}^{n}$ is the advection direction. The optimal resource control problem is stated as: Within the control set

[TABLE]

how to find $m^{*}\in U$ such that $J(m^{*})=\max_{m}J(m)$, where the objective functional is defined by

[TABLE]

subject to the system (1.1). The objective functional represents the net benefit, which is the size of the population less the cost of implementing the control. The coefficient $B>0$ is the parameter which balance the two parts of the objective functional.

The first term in $J(m)$, i.e. $\int_{Q_{T}}udxdt$, represents the total population over time and space, which not only serves as a good measurement for the conservation of a single species, but also plays an important role in preventing the invasion of alien species [22]. The second term $\int_{Q_{T}}Bm^{2}dxdt$ is a measurement of the cost of distributing the resource in the habitat. As a whole, $J(m)$ can be regarded as a way of determining the net benefit in the conservation of a single species.

Population movement and its distribution in response to its surrounding environment is an important issue in biology(see[1, 3, 4, 10, 11, 16, 18, 19, 20, 21, 17, 6, 7, 8, 12] and the references therein). Since the population abundance is a good measurement of conservation effort, it is more interesting to know how resource allocation affects population size of the species and it is a very challenging mathematical problem. For instance, given a fixed amount of resource, how can we determine the optimal spatial arrangement of the favorable and unfavorable parts of the habitat for species to survive? The question was first addressed by Cantrell and Cosner [3, 4] via the reaction-diffusion equation

[TABLE]

subject to Dirichlet, Robin, or Neumann boundary conditions, where $u(x,t)$ is the density of the species at location $x$ and time $t$, and constant $\lambda$ is the dispersal rate of the species and is assumed to be a positive constant. The coefficient $m(x)$ represents the intrinsic growth rate of the species and it measures the availability of the resource.

Among other things, Cantrell and Cosner [3] showed that there exists a “bang-bang” type optimal spatial arrangement of the favorable and unfavorable parts of the habitat for species to survive, i.e., the corresponding optimal control function $m(x)$ must be a step function in $\Omega$. When $\Omega$ is an interval, Cantrell and Cosner [4] showed that if the resource is so arranged that $m(x)$ is equal to some positive constant in one subinterval and is equal to some negative constant otherwise, then the optimal arrangement occurs when the subinterval with $m(x)$ positive is one of the two ends of the interval. It is further shown in [23] that any optimal control function $m$ must be of “bang-bang” type and when the domain $\Omega$ is an interval, there are exactly two optimal controls, for which the control $m(x)$ is positive at one end of the habitat provides the best opportunity for the species to survive. For high-dimensional habitats, very little is known about the exact shape of the optimal control.

Assume that the species can move up along the gradient of the density. In the field of ecology, organisms can often sense and respond to local environmental cues by moving towards favourable habitats, and these movement usually depend upon a combination of local biotic and abiotic factors such as stream, climate, food, chemical substance and predators. There are many examples involves advection can be found in the field of mathematical biology (see[5, 9, 15, 26, 27] and the references therein).

The reaction may result in movement with two features: directed advection and random diffusion [24, 25]. It is commonly believed that the population will move along the direction of increasing resources. With this regard, Belgacem and Cosner [1] investigated an reaction-diffusion-advection model in which the advection term is the gradient of the resource function. They investigated steady states solutions of the reaction-diffusion equations with linear and nonlinear logistic growth terms:

[TABLE]

and

[TABLE]

together with Dirichlet or Neumann boundary condition. They studied the benefit to the population[1], meaning the persistence of the population in the long run or the existence of a unique globally attracting positive steady state solution. It was found in [1] that the directed movement towards better resources could be beneficial to the population, while in the Dirichlet boundary condition, the movement towards better resources can be harmful if more favorable patches are closer to the boundary. In a further investigation [11], Cosner et al. studied the logistic reaction-diffusion models with the advection along the environment gradient with Neumann boundary condition, and they found that under the Neumann conditions, the movement along the resource gradient may not always be beneficial to the population. Indeed, it turns out that the convexity of the domains plays a major role in this situation. If the domain is not convex, moving up along the resource gradient could be harmful to the population.

Optimal control techniques were also used in other related work such as [14, 17, 12] to explore how different conditions such as limited resources, growth coefficient, advection movement, and harvesting can be optimized to be “beneficial” for populations. In the elliptic case, Ding et al. studied in [12] the effects of resource allocation on population size of the species. In their work, they investigated the maximizing the total population with the minimum cost for the resource of fixed amount.

In a similar framework as in [12], but different direction, Finotti et al. [14] investigated Equation (1.1) with the control on the advection direction $\overrightarrow{h}$, and they seek for the optimal advection direction $\overrightarrow{h}(x,t)$ that maximizes the total population while minimizing the “cost” due to movement.

In the current paper, we focus on the work using optimal control techniques to investigate Equation (1.1) with the resource control $m(x,t)$ with simplified $f(x,t,u)=u$ which can be extended to a more general function of logistic type. The main results of this paper can be stated as follows:

Theorem 1.1**.**

(Existence of the optimal control) Assume that $0<T<\infty,\vec{h}\in L^{\infty}(Q_{T}),u_{0}\in L^{\infty}(\Omega)\bigcap{H^{1}(\Omega)}$ and $u_{0}$ is non-negative. There exists an optimal control $m\in U$ maximizing the objective functional $J(m)$.

Theorem 1.2**.**

(The characterization of the optimal control) Given an optimal control $m^{*}$ and corresponding state $u^{*}$, there exists a solution $p$ in $L^{2}((0,T),H^{1}(\Omega))$ which satisfies $p_{t}\in L_{2}{((0,T),H^{1}(\Omega)^{*})}$ and

[TABLE]

Furthermore, $m^{*}$ is characterized by

[TABLE]

Theorem 1.3**.**

(The uniqueness of the optimal control) There exist two positive number $T_{0}$ and $B_{0}$ such that if $0<T<T_{0}$ and $B>B_{0}$, then there is a unique solution of the optimality system.

The rest of the paper is organized in the following way: In section 2, we present some preliminary lemmas to be used in the proofs of the main results. We prove the main results of the paper in section 3. The numerical simulation results are illustrated in section 4.

2 Preliminary lemmas

In this section, we present some results to guarantee the existence and Apriori estimate of a positive solution to Equation (1.1). We denote by $H^{1}(\Omega)$ the usual Sobolev space and its dual space is $H^{1}(\Omega)^{*}$. The space $V_{2}(Q_{T})$ is defined for all functions $u\in L^{2}((0,T),H^{1}(\Omega))$ such that its norm

[TABLE]

Definition 2.1**.**

The function $u\in L^{2}((0,T),H^{1}(\Omega))$ with $u_{t}\in L^{2}((0,T),H^{1}(\Omega)^{*})$ and $u(x,0)=u_{0}$ is said to be a weak solution of Equation (1.1) if and only if for a.e. $t\in(0,T)$

[TABLE]

In order to prove the existence of the optimal control, we need the following results and the proofs are similar in [14].

Lemma 2.1**.**

Assume that $\overrightarrow{h}$ and $m$ are in $L^{\infty}(Q_{T})$ and $u_{0}\geq 0.$ Then, any solution $u$ of Equation (1.1) must be non-negative on $Q_{T}$.

Which is natural that Equation (1.1) simulate the population density over $Q_{T}$, and the density should be non-negative.

Lemma 2.2**.**

Assume that $C_{m}=\sup\limits_{Q_{T}}|m|<\infty,u_{0}\in L^{\infty}(Q_{T})$ and $u_{0}\geq 0$. Then for each $m\in U$, any solution $u$ of Equation (1.1) satisfies

[TABLE]

and

[TABLE]

Lemma 2.3**.**

Let $0<T<\infty,\vec{h}\in{L^{\infty}(Q_{T})}$ and $u_{0}$ be non-negative,bounded and in $H^{1}(\Omega)$. Then, for each $m\in U$, there is a unique weak solution $u=u(m)$ of (1.1). Moreover, there is a finite constant $C>0$ depending only on $|\Omega|,\mu,M,T,\|u_{0}\|,\|h\|_{L^{\infty}},d$ such that

[TABLE]

The above lemmas show the solution of Equation (1.1) must be bounded and the bounds only depend on the bounds of $m(x,t)$ and $\overrightarrow{h}$ which will help to prove the existence of the optimal control problem.

In order to characterize the optimal control, we need to investigate the relationship about $m$, $u(m)$ and $J(m)$. For the given unique positive solution $u=u(m)$ of (1.1), the derivative of the mapping $m\rightarrow u(m)$ is called the sensitivity, and we have the following differentiability results of the mapping $m\rightarrow u(m).$

Lemma 2.4**.**

The mapping $m\in U\rightarrow u(m)$ is differentiable in the following sense: for each $m$, $l$ in $U$ such that $m+\varepsilon l\in U$ for all $\varepsilon$ sufficiently small, then there is a uniform constant $C>0$ such that $\varphi^{\varepsilon}=\frac{u(m+\varepsilon l)-u(m)}{\varepsilon}$ satisfies

[TABLE]

Moreover, there exists $\varphi=\varphi(m,l)\in L^{2}((0,T),H^{1}(\Omega))$, such that

[TABLE]

and the sensitivity $\varphi$ satisfies

[TABLE]

Proof.

First, let us denote $u^{\varepsilon}=u(m+\varepsilon l)$, $u=u(m)$ for each $m\in U$ and $m+\varepsilon l\in U$. Through that $u$ solves (1.1), we can prove that $u^{\varepsilon}$ solves

[TABLE]

Then, by Lemma 2.3, it follow from that

[TABLE]

where $C$ is a constant depending only on $\mu,d,|\Omega|,T,\|u_{0}\|_{L_{\infty}},\|h\|_{L_{\infty}(Q_{T})}$ and $M$. Define

[TABLE]

It is easy to see that $\gamma^{\varepsilon}(x,t)$ is uniformly bounded. In the meanwhile, it follows from the Theorem 1.1 and by the uniqueness of solution $u=u(m)$ of (1.1), we have

[TABLE]

Moreover, this gives

[TABLE]

Recall that $\varphi^{\varepsilon}=\frac{u^{\varepsilon}-u}{\varepsilon}$. Then subtracting (2.4) from (1.1) and dividing by $\varepsilon$, we obtain

[TABLE]

As in the proofs of Theorem 1.1, we have

[TABLE]

Then we can assume that $\varphi^{\varepsilon}\rightharpoonup\varphi$ in $L^{2}((0,T),H_{1}(\Omega))$. Using the estimates (2.5), the equation (2.6), and as in the proof Theorem 1.1, we obtain (2.3). This concludes the proof of Lemma 2.4. ∎

Lemma 2.5**.**

Let $p$ be the solution of (1.4). Then, there is a constant $C$ such that

[TABLE]

where $C$ depends on $|\Omega|,\mu,M,T,\|u_{0}\|,\|h\|_{L^{\infty}},d$.

Proof.

Let us denote $q(x,t)=p(x,T-t)$, then it follows that $q$ solves the equation

[TABLE]

where $b(x,t)=-\vec{h}(x,T-t)$ and $c(x,t)=2u^{*}(x,T-t)-m^{*}(x,T-t)$ for all $(x,t)\in Q_{T}$. From Lemma 2.3, it follows Equation (2.7) is a linear parabolic equation with bounded coefficients. By the maximum principle, it follows that $p\geq 0$. On the other hand, by the parabolic regularity theory, we have

[TABLE]

where

[TABLE]

Then it follows from the Sobolev embedding theorem that for $l>d+2$, we have

[TABLE]

which yields the desired estimates. ∎

3 The proof of the main results

In this section, we mainly provide the proof for Theorem 1.1, Theorem 1.2 and Theorem 1.3.

Proof of Theorem 1.1.

First,we show the existence of $\sup\limits_{m\in U}J(m)$. By Lemma 2.2, $J(m)\leq C$ where $C$ is a constant depending on $\|m\|_{L^{\infty}(Q_{T})}$ $\|u_{0}\|_{L^{\infty}}$ and $T$ only. It is obvious that there exist a maximizing sequence $\{m^{n}\}\in U$ such that

[TABLE]

Denote $u^{n}=u(m^{n})$, $u^{n}$ is the corresponding solutions of (1.1) and the control $m$ is $m^{n}$. Next,we need to show the estimates of $u^{n},\nabla u^{n}$ and $u_{t}^{n}$. By the theorem 2.3,we obtain that

[TABLE]

for some constant $C>0$ depending only on $\mu,d,|\Omega|,T,\|u_{0}\|_{L_{\infty}},\|h\|_{L_{\infty}(Q_{T})}$ and $M$, where $M$ is a fixed constant defined in (1.2). From(3.2), we can assume that

[TABLE]

And for each $n$ and each $\phi\in L^{2}((0,T),H^{1}(\Omega))$, the weak form of the solution $u^{n}$ is

[TABLE]

Thus, by combining (3.2), (3.4) and H$\ddot{o}$lder inequality, it shows that

[TABLE]

Therefore,

[TABLE]

where the constant $C$ depends only on $\mu,d,|\Omega|,T,\|u_{0}\|_{L_{\infty}},\|h\|_{L_{\infty}(Q_{T})}$ and $M$. Using the results in [2] and the forms (3.2),(3.3),(3.5), it follows that

[TABLE]

Due to the weakly compactness of $L^{2}(\Omega)$ and the bounded control set $U$ defined in (1.2), there exist $m^{*}\in U$ such that

[TABLE]

By the definition of weak solution, $u^{n}$ satisfies

[TABLE]

From the weak convergence in (3.3), (3.6) and (3.7), it follows

[TABLE]

According to the strong convergence of the sequence $\{u^{n}\}_{n\in N}$ and the weak convergence of the sequence $\{m^{n}\}_{n\in N}$, we conclude

[TABLE]

which leads

[TABLE]

Collecting those convergence terms in(3.6), (3.9) and (3.11), and using the equation (3.8), we arrive at

[TABLE]

which implies that $u^{*}$ is the solution of (1.1) with respect to the control $m^{*}$. That is to say $u^{*}=u(m^{*})$. On the other hand, using the strong convergence in $L^{2}(Q_{T})$ of the sequence $\{u^{n}\}_{n\in N}$, and the fact that the function $m\mapsto\int_{Q_{T}}m^{2}\mathrm{d}x\mathrm{d}t$ is weakly lower semi-continuous in $L^{2}(Q_{T})$, we also get

[TABLE]

This implies that $J(m^{*})=\sup\limits_{m\in U}J(m)$. Therefore, $m^{*}\in U$ is an optimal control and the proof of Theorem 1.1 is complete. ∎

Proof of Theorem 1.2.

Suppose $m^{*}$ is an optimal control. Let $l\in U$ such that $m^{*}+\varepsilon l\in U$ for sufficiently small $\varepsilon>0$ and denote $u^{\varepsilon}=u(m^{*}+\varepsilon l)$ be the unique solution of (1.1) when the control term is $m^{*}+\varepsilon l$.

The operator in the adjoint equation is the formal analysis “adjoint” of the operator in the sensitivity equation (2.3) at $m^{*}$. Equation (1.4) is linear in $p$ and its coefficients are measurable and bounded. By the change of variable $t\rightarrow T-t$, the existence and uniqueness of the weak solution $p$ of (1.4) follows by Galerkin’s method (see [13]).

Observe that the directional derivative of $J$ with respect to the control at $m^{*}$ in the direction of $l$ satisfies

[TABLE]

Using the weak solution formulation for the adjoint problem with test function $\varphi$, we obtain

[TABLE]

By Theorem 2.3 and Lemma 4.3 below, we know that there is a constant $C_{0}>0$ such that

[TABLE]

For each $0<\delta<1$, let $\tau_{\delta}$ be the set $\{(x,t)\in Q_{T}:|m(x,t)|\leq\delta M\}$. From (3.15), let $l=u^{*}p-2Bm^{*}$ and we obtain that

[TABLE]

Thus, $m^{*}=\frac{u^{*}p}{2B}$ on $\tau_{\delta}$. Since $\delta$ is arbitrary, we conclude that on the set where $|m^{*}|<M$,

[TABLE]

Now assume that $m^{*}=M$ on some non-empty $\tau\subseteq Q_{T}$. Thus, $m^{*}=M\leq\frac{u^{*}p}{2B}$. Therefore, we conclude that $m^{*}=\min\{M,\frac{u^{*}p}{2B}\}$. Since $0\leq m\leq M$, we can show that

[TABLE]

This completes the proof of Theorem 1.2. ∎

Proof of Theorem 1.3.

In the above, the existence of the optimal control and corresponding adjoint and states have been proved. In the following, we prove the uniqueness of the system.

Let $m^{*}_{1}$ and $m^{*}_{2}$ be two controls corresponding to solutions of the optimal system. Denote $u_{i}=u(m^{*}_{i})$, $p_{i}=p(m^{*}_{i})$ for $i=1,2$ to be the state solution and the solution of the adjoint problem (1.4). For some $\lambda>0$ which will be determined and denote $u_{i}=e^{\lambda t}w_{i}$ and $p_{i}=e^{-\lambda t}z_{i}$. If we let $W=w_{1}-w_{2},Z=z_{1}-z_{2}$, and $\widetilde{m}=m^{*}_{1}-m^{*}_{2}$, we then obtain

[TABLE]

By subtracting the equation of $u_{1}$ and $u_{2}$, we see that $W$ solves

[TABLE]

Similarly, $Z$ also solves

[TABLE]

Multiplying (3.17) by $W$ and using the integration by parts, Holders’s inequality, Young’s inequality, we get

[TABLE]

Thus,

[TABLE]

Integrating this inequality with respect to time, we also get

[TABLE]

Similar process to Equation (3.18), we obtain

[TABLE]

Thus,

[TABLE]

Integrating the above inequality with respect to time, we obtain

[TABLE]

Note that the constants $C_{1},C_{2},C_{3},C_{4}$ all depend on $|\Omega|,\mu,M,T,\|u_{0}\|,\|h\|_{L^{\infty}},d$. Then there exists a sufficiently small $T_{0}>0$ with $T\leq T_{0}$, sufficiently large positive numbers $B_{0}$ and $\lambda$ so that

[TABLE]

It follows from (3.19) and (3.20) that

[TABLE]

which leads

[TABLE]

This implies $u_{1}=u_{2},p_{1}=p_{2}$. Thus, $m^{*}_{1}=m_{2}^{*}$. ∎

4 Numerical simulation results

We have run several examples for the case when $f(x,t,u)=u$, with both time-independent and time-dependent function $\vec{h}$. In order to solve the optimality system, we use an iterative scheme with an explicit finite difference method. First, starting with an initial guess for the control function $m$ and using a forward-backward sweep method [21], we approximate the first state and adjoint. We then obtain the next approximation to the optimal control by evaluating our optimal control characterization. Continuing the above iteration process until the optimal state and optimal control converge.

Let $\mu=0.2$ and the final time $T=1$ in all simulations. Solution $u$ and the optimal control are depend on the settings of control set (1.2) and objective functional (1.3), initial values $u_{0}$, advection function $\vec{h}$, we illustrate each effect in the following with one variable changing and others fixed.

The effect of the advection function. The evolution of population is affected by the diffusion and advection and the initial solution. We compare four optimal controls and states with time-independent and time-dependent advection functions $\vec{h}$ and these simulations show the optimal control has strong relations with the initial solution and the advection. The optimal strategies are concentrating more resource at the maximum value point of initial solution which is corresponding to the favorable habitat in ecology and moving leftward along the positive advection direction. In the simulation, we set $M=10$ and $B=0.1$ and let $u_{0}$ be the function

[TABLE]

In the first example, we set $\vec{h}$ as a time-independent function

[TABLE]

whose graph is shown in Figure 1. $\vec{h}_{1}$ is a positive function. Due to the influence of advection, the species will move leftward with positive advection function. The corresponding optimal control and state functions are shown in Figure 2.

In the second example, we consider a time-independent quadratic advection function

[TABLE]

whose graph is shown in Figure 3 and the corresponding optimal control and state are shown in Figure 4.

The third example is for a time-dependent advection function

[TABLE]

whose graph is shown in Figure 5 and the corresponding graph of the optimal control and state are in Figure 6.

The fourth example is for a time-dependent advection function

[TABLE]

whose graph is shown in Figure 7 and the corresponding optimal control and state are in Figure 8.

In the following, we show the time slices of the optimal state at early, middle and final stages. For simplicity, we let $\vec{h}$ be a constant function, $\vec{h}=0.1$. Figure 9 shows that the species move leftward with positive advection function and the optimal resource increase when the density increase over time.

The effect of the initial values. The evolution of the population is affected by the initial solution. In this part, we fix other variables and only let the initial solution change. These figure shows the optimal strategy should concentrate the resource at the maximum value of the initial solution. Let the advection function be

[TABLE]

$B$ is taken to be $B=0.1$ and have sufficient resources which means $M$ is big enough. First, we let $u_{0}$ be the function

[TABLE]

We find where has more species, the more resources should be allocated and the corresponding graphs of the optimal control and state are in Figure 10.

The second initial function is

[TABLE]

whose optimal strategy is similar to the first one and the corresponding graph of the optimal control and state are given in Figure 11.

To show a clear relation between $m$ and $u$, we compare the resulting optimal control and corresponding state at time $t=0.1$ in Figure 12 with $u_{0}=\sin(2\pi x)+1$ and in Figure 13 with $u_{0}=2x^{2}-\frac{1}{2}x+1$.

The effect of the maximum resource and the cost constant. In the real life, the distribution of resource may have up-bound $M$ at one specific location which means the resource may not be sufficient. In the following, three different up bounds $M$ of the maximum resource are compared. The maximum resource $M$ affects the distribution of the densities, however, the trends in densities’ distributions are almost the same. In the simulation, we set

[TABLE]

[TABLE]

and the cost constant $B=0.1$. We set $M=1,10,20$. When $M=1$, the graph of the optimal control and state is given in Figure 14. The second case have been discussed before, the graph is same as Figure 2. The graph of the optimal control and state is given in Figure 15 for $M=20$.

We also investigate the effect of the cost constant on the optimal control. Setting $\vec{h}$, $u_{0}$ the same as in (4.1), (4.2), the resource restrictions are the same. Then by changing $B$ from $B=0.05$ to $B=1$, the graph of the optimal control and state are respectively given in Figure 16 and Figure17.

For different $B$ values, we present the resulting optimal controls at time $t=0.2$ in Figure 18 and Figure 19.

From the graphs we can observe that the changing of $B$ have little effect on the trends of the optimal control, but changes its scale.

The numerical simulation with Dirichlet boundary. The Dirichlet boundary condition corresponds to the hostile environment in ecology. In the following, we simulate the system with Dirichlet boundary condition. Let

[TABLE]

and $\mu=0.2$, $B=0.1$, and $M=10$. With the above initial condition, the center of the region has positive density distribution and has zero density near the hostile boundary. For simplicity, we first set

[TABLE]

which corresponds to no advection, the species have a stable growth around the center of habitat and the optimal control strategy is allocate resource at the center. The graph of the optimal control and state are given in Figure 20.

Secondly, we set $\vec{h}(x)=3x$, which means that the advection $\vec{h}$ is positive on the habitat. The species will move toward left side and the optimal control strategy should be moving the resource allocation along the advection direction. Because of the advection in one direction to the boundary, the species will die out in the habitat when time is long enough. The graph of the optimal control and state are given in Figure 21.

Thirdly, we set $\vec{h}(x)=-10x$. The species will move toward right side and the appearances and trends should be opposite to the positive advection case. The graph of the optimal control and state are given in Figure 22.

The advection direction affect the distribution of density and the corresponding optimal control strategy. Because of the harsh environment in the boundary and under one direction of advection, the species will die out as long as time is big enough.

5 Acknowledgement

We thank Prof. Yuan Lou for introducing this topic. Lianzhang Bao was partially supported by China Postdoctoral Science Foundation –183816.