Some regional control problems for population dynamics

TL;DR

This paper studies regional control strategies for population dynamics, focusing on optimal harvesting and eradication of pest populations, using theoretical analysis and numerical methods to optimize interventions in specific regions.

Contribution

It introduces a novel iterative numerical method for optimizing regional harvesting efforts and provides new theoretical conditions for population eradication control.

Findings

Optimal harvesting effort structure derived from necessary conditions.

Iterative method effectively increases total harvest in numerical tests.

Conditions for zero-stabilizability of pest populations established.

Abstract

This paper deals with some control problems related to structured population dynamics with diffusion. Firstly, we investigate the regional control for an optimal harvesting problem (the control acts in a subregion $\omega$ of the whole domain $\Omega$). Using the necessary optimality conditions, for a fixed $\omega$, we get the structure of the harvesting effort which gives the maximum harvest; with this optimal effort we investigate the best choice of the subregion $\omega$ in order to maximize the harvest. We introduce an iterative numerical method to increase the total harvest at each iteration by changing the subregion where the effort acts. Numerical tests are used to illustrate the effectiveness of the theoretical results. We also consider the problem of eradication of an age-structured pest population dynamics with diffusion and logistic term, which is a zero-stabilization problem with constraints. We derive a necessary condition and a sufficient condition for zero-stabilizability. We formulate a related optimal control problem which takes into account the cost of intervention in the subregion $\omega$.

Equations

$$\left\{\begin{array}[]{ll}\partial_{t}y(x,t)-d\Delta y(x,t)=a(x)y(x,t)-\chi_{\omega}(x)u(x,t)y(x,t),&(x,t)\in Q_{T}\\ \partial_{\nu}y(x,t)=0,&(x,t)\in\Sigma_{T}\\ y(x,0)=y_{0}(x),&x\in\Omega,\\ \end{array}\right.$$

$$Maximize{\int }_0^T{\int }_{\omega }u(x,t)y^u(x,t)dxdt,$$

$$\left\{\begin{array}[]{ll}\partial_{t}p(x,t)+d\Delta p(x,t)=-a(x)p(x,t)\\ ~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\chi_{\omega}(x)u^{*}(x,t)(1+p(x,t)),&(x,t)\in Q_{T}\\ \partial_{\nu}p(x,t)=0,&(x,t)\in\Sigma_{T}\\ p(x,T)=0,&x\in\Omega,\\ \end{array}\right.$$

$$u^{*}(x,t)=\left\{\begin{array}[]{l l}0,&\quad 1+p(x,t)<0\\ L,&\quad 1+p(x,t)>0\end{array}\right.\quad\mbox{a.e.}\ (x,t)\in\omega\times(0,T).$$

$$H(z)=\left\{\begin{array}[]{ll}1,\quad\textnormal{if}\quad z\geq 0\\ 0,\quad\textnormal{if}\quad z<0.\\ \end{array}\right.$$

$$\underset{\phi }{\mathrm{Maximize}}\underset{u\in K_{\omega }}{\mathrm{Maximize}}\{{\int }_0^T{\int }_{\omega }u(x,t)y^u(x,t)dxdt$$

$$-\alpha length(\partial \omega )-\beta area(\omega )\},$$

$$\underset{\phi }{\mathrm{Maximize}}\underset{u\in K_{\omega }}{\mathrm{Maximize}}\{{\int }_0^T{\int }_{\omega }u(x,t)y^u(x,t)dxdt$$

$$-\alpha {\int }_{\Omega }\delta (\phi (x))|\nabla \phi (x)|dx-\beta {\int }_{\Omega }H(\phi (x))dx\}.$$

$$\underset{\phi }{\mathrm{Maximize}}\{{\int }_0^T{\int }_{\omega }{u}_{\phi }^{\ast }(x,t)y_{\phi }^{\ast }(x,t)dxdt$$

$$-\alpha {\int }_{\Omega }\delta (\phi (x))|\nabla \phi (x)|dx-\beta {\int }_{\Omega }H(\phi (x))dx\},$$

$$\left\{\begin{array}[]{ll}\partial_{t}y(x,t)-d\Delta y(x,t)=a(x)y(x,t)-H(\varphi(x))u_{\varphi}^{*}(x,t)y(x,t),&(x,t)\in Q_{T}\\ \partial_{\nu}y(x,t)=0,&(x,t)\in\Sigma_{T}\\ y(x,0)=y_{0}(x),&x\in\Omega.\\ \end{array}\right.$$

$$u^{*}_{\varphi}(x,t)=\left\{\begin{array}[]{l l}0,&\quad 1+p_{\varphi}(x,t)<0\\ L,&\quad 1+p_{\varphi}(x,t)>0\end{array}\right.$$

$${\int }_0^T{\int }_{\Omega }[{\partial }_t{y}_{\phi }^{\ast }{p}_{\phi }+y_{\phi }^{\ast }{\partial }_t{p}_{\phi }]dxdt+{\int }_0^T{\int }_{\Omega }[-d\Delta y_{\phi }^{\ast }{p}_{\phi }+d\Delta p_{\phi }{y}_{\phi }^{\ast }]dxdt+$$

$$+{\int }_0^T{\int }_{\Omega }[-a(x)y_{\phi }^{\ast }{p}_{\phi }+a(x)y_{\phi }^{\ast }{p}_{\phi }]dxdt=$$

$$={\int }_0^T{\int }_{\Omega }[-H(\phi (x))u_{\phi }^{\ast }{y}_{\phi }^{\ast }{p}_{\phi }+H(\phi (x))u_{\phi }^{\ast }{y}_{\phi }^{\ast }+H(\phi (x))u_{\phi }^{\ast }{y}_{\phi }^{\ast }{p}_{\phi }]dxdt.$$

$$-{\int }_{\Omega }{y}_0(x)p_{\phi }(x,0)dx={\int }_0^T{\int }_{\Omega }H(\phi (x))u_{\phi }^{\ast }(x,t)y_{\phi }^{\ast }(x,t)dxdt$$

$${\int }_0^T{\int }_{\omega }{u}_{\phi }^{\ast }(x,t)y_{\phi }^{\ast }(x,t)dxdt=-{\int }_{\Omega }{y}_0(x)p_{\phi }(x,0)dx$$

$$\underset{\phi }{\mathrm{Minimize}}\left\{{\int }_{\Omega }{y}_0(x)p_{\phi }(x,0)dx+\alpha {\int }_{\Omega }\delta (\phi (x))|\nabla \phi (x)|dx+\beta {\int }_{\Omega }H(\phi (x))dx\right\}.$$

$$\left\{\begin{array}[]{ll}\partial_{t}p+d\Delta p=-a(x)p+LH(\varphi(x))(1+p)H(1+p),&(x,t)\in Q_{T}\\ \partial_{\nu}p(x,t)=0,&(x,t)\in\Sigma_{T}\\ p(x,T)=0,&x\in\Omega.\end{array}\right.$$

$$\underset{\phi }{\mathrm{Minimize}}J(\phi ),$$

$$J(\phi )={\int }_{\Omega }{y}_0(x)p_{\phi }(x,0)dx+\alpha {\int }_{\Omega }{\delta }_{\epsilon }(\phi (x))|\nabla \phi (x)|dx+\beta {\int }_{\Omega }{H}_{\epsilon }(\phi (x))dx,$$

$$\left\{\begin{array}[]{ll}\partial_{t}p+d\Delta p=-a(x)p+LH_{\varepsilon}(\varphi(x))(1+p)H_{\varepsilon}(1+p),&(x,t)\in Q_{T}\\ \partial_{\nu}p(x,t)=0,&(x,t)\in\Sigma_{T}\\ p(x,T)=0,&x\in\Omega.\end{array}\right.$$

$$dJ(\phi )(\psi )={\int }_{\Omega }{\delta }_{\epsilon }(\phi (x))[-\alpha \text{div}\left(\frac{\nabla \phi (x)}{|\nabla \phi (x)|}\right)+\beta$$

$$-L{\int }_0^T(1+p_{\phi }(x,t))H_{\epsilon }(1+p_{\phi }(x,t))r_{\phi }(x,t)dt]\psi (x)dx+\alpha {\int }_{\partial \Omega }\frac{{\delta }_{\epsilon }(\phi (x))}{|\nabla \phi (x)|}{\partial }_{\nu }\phi (x)\psi (x)d\sigma ,$$

$$\left\{\begin{array}[]{ll}\partial_{t}r-d\Delta r=a(x)r-LH_{\varepsilon}(\varphi(x))H_{\varepsilon}(1+p_{\varphi})r&{}\\ \ \ \ \ \ \ \ \ \ \ \ -LH_{\varepsilon}(\varphi(x))(1+p_{\varphi})\delta_{\varepsilon}(1+p_{\varphi})r,&(x,t)\in Q_{T}\\ \partial_{\nu}r(x,t)=0,&(x,t)\in\Sigma_{T}\\ r(x,0)=y_{0}(x),&x\in\Omega.\\ \end{array}\right.$$

$$\frac{1}{\theta }[p_{\phi +\theta \psi }-p_{\phi }]\to q_{\psi }\text{mbox}\mathrm{in}C([0,T];L^{\infty }(\Omega )),$$

$$\left\{\begin{array}[]{ll}\partial_{t}q+d\Delta q=-a(x)q+L\delta_{\varepsilon}(\varphi(x))(1+p_{\varphi})H_{\varepsilon}(1+p_{\varphi}){\psi}(x)&{}\\ \ \ \ \ \ \ \ \ \ \ \ +LH_{\varepsilon}(\varphi(x))H_{\varepsilon}(1+p_{\varphi})q+LH_{\varepsilon}(\varphi(x))(1+p_{\varphi})\delta_{\varepsilon}(1+p_{\varphi})q,&(x,t)\in Q_{T}\\ \partial_{\nu}q(x,t)=0,&(x,t)\in\Sigma_{T}\\ q(x,T)=0,&x\in\Omega.\\ \end{array}\right.$$

$$\underset{\theta \to 0}{\lim }\frac{1}{\theta }[J(\phi +\theta \psi )-J(\phi )]={\int }_{\Omega }{y}_0(x)q_{\psi }(x,0)dx+\alpha {\int }_{\Omega }{\delta }_{\epsilon }^{\prime }(\phi (x))\psi (x)|\nabla \phi (x)|dx$$

$$+\alpha {\int }_{\Omega }{\delta }_{\epsilon }(\phi (x))\frac{\nabla \phi (x)\cdot \nabla \psi (x)}{|\nabla \phi (x)|}dx+\beta {\int }_{\Omega }{\delta }_{\epsilon }(\phi (x))\psi (x)dx$$

Full text

Keywords Optimal harvesting; population dynamics; diffusive models; regional control; numerical methods.

1 Introduction

An extensive literature was developed for the optimal harvesting problems of population dynamics (e.g. [1]–[4], [5], [8], [9], [12], [14]–[16], [17], [18]–[20], [22], [23]). In this paper we firstly remind an optimal harvesting problem for a spatially structured population with diffusion which has been introduced in [4]. For spatially structured harvesting problems it has usually been taken into consideration an effort that acts in the whole habitat $\Omega$ (see for example [2]). Instead here we consider the case in which the effort is localized in a suitably chosen subregion $\omega$ of $\Omega.$ In addition to the problem of finding the magnitude of of the control to act on a given subdomain $\omega,$ the most important task will be to identify an optimal subregion $\omega$, where the control acts, in order to maximize the harvest. To this aim, at first we have derived necessary optimality conditions for the situation when the support of the control is fixed; as a fall out we have obtained information concerning the structure of the optimal control. Hence we have taken into account this structure to investigate the optimal subregion $\omega$ where the control is localized, by taking into account the cost paid for harvesting in $\omega$. Here we have adapted some shape optimization methods, based on the level set method. These results have been previously presented in [4]. In this paper we consider also the problem of eradication of an age-structured pest population dynamics with diffusion and logistic term. We consider a related optimal control problem which can be again investigated by means of the level set method.

We consider the following population dynamics model with diffusion. A single population species is free to move in an isolated habitat $\Omega\subset{\mathbb{R}^{2}}$, with $\Omega$ a bounded domain with a sufficiently smooth boundary:

[TABLE]

where $Q_{T}=\Omega\times(0,T)$, $\Sigma_{T}=\partial\Omega\times(0,T)$, $T>0$, $y=y(x,t)$ is the population density at position $x\in\overline{\Omega}$ and time $t\in[0,T]$, while $y_{0}(x)$ is the initial population density. Here $a(x)$ denotes the natural growth rate of the population, and $d\in(0,+\infty)$ is the diffusion coefficient. No-flux boundary conditions are considered.

In System (1), $u(x,t)$ represents the harvesting effort (control), bounded and localized in the subdomain $\omega\subset\Omega$ ($\chi_{\omega}$ is the characteristic function of $\omega$). The term $\chi_{\omega}(x)u(x,t)y(x,t)$ represents the rate of the harvested population at position $x\in\omega$ and time $t\in[0,T]$.

The following hypotheses are considered:

(H1)

$a\in L^{\infty}(\Omega)$;

(H2)

$y_{0}\in L^{\infty}(\Omega),\quad y_{0}(x)\geqslant 0\quad\mbox{a.e.}\ x\in\Omega\,$ with $\,\|y_{0}\|_{L^{\infty}(\Omega)}>0$.

We consider a related optimal harvesting problem

[TABLE]

subject to $u\in K_{\omega}$, where $K_{\omega}=\{w\in L^{\infty}(\omega\times(0,T));\ 0\leq w(x,t)\leq L\ \ \mbox{\rm a.e. in }\omega\times(0,T)\}$. Here $L>0$ is a constant and $y^{u}$ is the solution to (1) corresponding to a harvesting effort $u\in K_{\omega}$.

The existence result of an optimal control for Problem (2) follows [3] or [5].

Theorem 1**.**

Problem (2) admits at least one optimal control.

We denote by $p$ the adjoint state, i.e. $p$ satisfies

[TABLE]

where $(u^{*},y^{u^{*}})$ is an optimal pair for (2). For the construction of the adjoint problems in optimal control theory we refer to [8]. Concerning the first order necessary optimality conditions it can be proved the following result (as in [3] and [5]):

Theorem 2**.**

If $(u^{*},y^{u^{*}})$ is an optimal pair for Problem (2) and if $p$ is the solution of Problem (3), then we have:

[TABLE]

In Section 2 we will treat the regional harvesting problem as a shape optimization problem. We remind that the geometry of a set $\omega$ can be characterized in terms of its Minkowski functionals. There are three such functionals and these are proportional to the area, the perimeter and the Euler-Poincaré characteristic. In this paper we control the shape of $\omega$ as follows: by the length of the boundary of $\omega$, and by the area of $\omega$.

We shall use the implicit interface representation to control the shape of the 2D domain $\omega$. Therefore, the boundary of a domain is defined as the isocontour of some function $\varphi$ (see [11] or [21]). By using the level set method, we introduce a level set function $\varphi:\overline{\Omega}\rightarrow\mathbb{R}$ such that $\omega=\{x\in\Omega;\ \varphi(x)>0\ \mbox{\rm a.e.}\}$ and $\partial\omega=\{x\in\Omega:\varphi(x)=0\ \mbox{\rm a.e.}\}$ (the boundary is defined as the zero level set of $\varphi$). We will then manipulate $\omega$ implicitly, through the function $\varphi$. This function $\varphi$ is assumed to take positive values inside the region delimited by the curve $\partial\omega$ and negative values outside.

If $\varphi$ is the implicit function of $\omega,$ in order to integrate over $\omega$ a function $f$ defined over the whole $\Omega$ we may write $\int_{\Omega}f(x)H(\varphi(x))dx,$ where we have used the Heaviside function $H:{\mathbb{R}}\rightarrow\{0,1\},$ such that

[TABLE]

If $\varphi$ is sufficiently smooth, the directional derivative of the Heaviside function in the normal direction at a point $x\in\partial\omega$ is given by $H^{\prime}(\varphi(x))|\nabla\varphi(x)|,$ and by using the usual Dirac Delta $\delta$ on $\mathbb{R}$, we have $\delta(\varphi(x))|\nabla\varphi(x)|$. If we need to integrate over $\partial\omega$ a function $f$ defined over the whole $\Omega$ we may write $\int_{\Omega}f(x)\delta(\varphi(x))|\nabla\varphi(x)|dx$.

We find the derivative of the optimal cost value with respect to the implicit function $\varphi$ of the subregion $\omega$. In order to improve the region where the control acts we derive a conceptual iterative algorithm based on these theoretical results. We also present the numerical implementation of this conceptual algorithm and some numerical tests. Basically, the theoretical results in Section 2 have been obtained in [4]. Here we give some additional details concerning the numerical scheme and its implementation. Further we present here some new numerical tests.

In Section 3 we treat the problem of eradication of an age-structured pest population with diffusion, which is a zero-stabilization problem with constraints. We derive a necessary condition and a sufficient condition of zero-stabilization. We consider a related optimal control problem which takes into account the cost paid by acting in the subregion $\omega$. We formulate this optimal control problem by means of the level set method. The results in this section are new.

2 An iterative method to localize an optimal subdomain $\omega$ where the control

acts

Here we intend to use the level set method in order to obtain the optimal subregion $\omega$ where the control is localized. Consider $\varphi:\overline{\Omega}\rightarrow\mathbb{R}$ the implicit function of $\omega$, the subregion of $\Omega$ where the control acts.

We rewrite the optimal control problem (2) such that will include both the magnitude of the harvesting effort $u\in K_{\omega}$, and the choice of the subdomain $\omega$ with respect to its implicit function $\varphi$:

[TABLE]

[TABLE]

where $y^{u}$ is the solution to (1) corresponding to a harvesting effort $u\in K_{\omega}$ and $\alpha,\beta$ are positive constants. $\alpha\;\text{length}(\partial\omega)+\beta\;\text{area}(\omega)$ represents the cost paid to harvest in the subregion $\omega$.

By using De Giorgi’s formula for the length (perimeter) of a set and assuming that $\varphi$ is sufficiently smooth , the optimal problem becomes

[TABLE]

[TABLE]

We have now two maximization problems: firstly, for a fixed $\varphi$ (and implicitly, $\omega$) we have to find the structure of the harvesting effort which gives the maximum harvest, as a function of $\varphi$ (or $\omega$); secondly, using this structure of the optimal control we investigate the optimal choice of the subregion $\omega$ with respect to its implicit function $\varphi$ in order to maximize the harvest.

For any arbitrary but fixed $\varphi$, we denote by $(u_{\varphi}^{*},y_{\varphi}^{*})$ an optimal pair for the harvesting problem (2). Now we have to investigate the following optimal control problem:

[TABLE]

[TABLE]

where $y_{\varphi}^{*}$ is the solution to

[TABLE]

Note that $H(\varphi)$ represents the characteristic function of $\omega$.

Assume that the hypotheses (H1, H2) are satisfied. We denote by $p_{\varphi}$ the adjoint state. From Theorem 2, the optimal control is given by

[TABLE]

a.e. $(x,t)\in\omega\times(0,T)$, where $p_{\varphi}$ is the solution to (3).

By multiplying (1) by $p_{\varphi}$ and (3) by $y_{\varphi}^{*}$, and integrating both of them on $Q_{T}$ we obtain:

[TABLE]

[TABLE]

[TABLE]

This means that

[TABLE]

and therefore

[TABLE]

Our problem of optimal harvesting becomes a problem of minimizing another functional with respect to the implicit function of $\omega$. Therefore, we may rewrite the optimal problem as

[TABLE]

By using (3) and (4) we get that $p_{\varphi}$ is the solution to

[TABLE]

We shall adapt some shape optimization techniques to treat this last harvesting problem (see also [10], [13]). As usual, we will approximate this problem by the following one, where the Heaviside function $H$ is substituted by its mollified version $H_{\varepsilon}(t)=\frac{1}{2}\left(1+\frac{2}{\pi}\arctan\left(\frac{t}{\varepsilon}\right)\right),$ and its derivative by the mollified function $\delta_{\varepsilon}(t)=\frac{\varepsilon}{\pi(\varepsilon^{2}+t^{2})}.$

Therefore, for a small but fixed $\varepsilon>0$, the harvesting problem to be investigated is:

[TABLE]

where $\varphi:\overline{\Omega}\longrightarrow\mathbb{R}$ is a smooth function,

[TABLE]

and $p_{\varphi}=p_{\varphi}(x,t)$ is the solution to

[TABLE]

In the following we derive the directional derivative of $J$ (see [4]).

Theorem 3**.**

For any smooth functions $\varphi,\psi:\overline{\Omega}\longrightarrow\mathbb{R}$ we have that

[TABLE]

[TABLE]

where $r_{\varphi}$ is the solution to

[TABLE]

Proof.

It is possible to prove that

[TABLE]

as $\theta\rightarrow 0$, where $q_{\psi}$ is the solution to the problem

[TABLE]

For any $\varphi$ we have that

[TABLE]

[TABLE]

Using Gauss-Ostrogradski’s formula we get

[TABLE]

By multiplying the first equation in (7) by $r_{\varphi}$ and integrating over $Q_{T}$, using (6) we get that

[TABLE]

Now, from (8) and (9), we get the conclusion of Theorem 3. ∎

Let us remark that the gradient descent with respect to $\varphi$ is

[TABLE]

($\theta$ is an artificial time).

2.1 Numerical implementation

From Theorem 3 we derive the following conceptual iterative algorithm, a semi-implicit gradient descent method, to improve at each step the region where the harvesting effort acts in order to obtain a smaller value for $J$.

STEP 0: set $n:=0$, $J^{(0)}:=10^{6}$ and $\theta_{0}>0$ a small constant

initialize $\varphi^{(0)}=\varphi^{(0)}(x,0)$

STEP 1: compute $p^{(n+1)}$ the solution of (5) corresponding to $\varphi^{(n)}(\cdot,0)$

compute $\displaystyle J^{(n+1)}=\int_{\Omega}y_{0}(x)p^{(n+1)}(x,0)dx$

$+\alpha\;\int_{\Omega}\delta_{\varepsilon}(\varphi^{(n)}(x,0))|\nabla\varphi^{(n)}(x,0)|dx\displaystyle+\beta\;\int_{\Omega}H_{\varepsilon}(\varphi^{(n)}(x,0))dx$.

Step 2: if $\left|J^{(n+1)}-J^{(n)}\right|<\varepsilon_{1}$ or $J^{(n+1)}\geq J^{(n)}$ then STOP

else go to Step 3.

Step 3: compute $r^{(n+1)}$ the solution of problem (6) corresponding to

$\varphi^{(n)}(\cdot,0)$ and $p^{(n+1)}$.

Step 4: compute $\varphi^{(n+1)}$ using (10) and the initial condition

$\varphi^{(n+1)}(x,0)=\varphi^{(n)}(x,\theta_{0})$ and a semi-implicit timestep scheme

Step 5: if $\|\varphi^{(n+1)}-\varphi^{(n)}\|_{L^{2}}<\varepsilon_{2}$ then STOP

else $n:=n+1$

go to Step 1

$\varepsilon_{1}>0$ in Step 2 and $\varepsilon_{2}>0$ in Step 5 are prescribed convergence parameters.

For the implementation we consider $\Omega=(0,1)\times(0,1)$ such that the sides are parallel with $Ox_{1}$ and $Ox_{2}$ axes. We introduce equidistant discretization nodes for both axes corresponding to $\Omega$. Thus, the domain $\Omega$ is approximated by a grid of $(N+1)\times(N+1)$ equidistant nodes, namely

[TABLE]

The interval $[0,T]$ is also discretized by $M+1$ equidistant nodes, $t^{k}=(k-1)\Delta t,k=1,2,...M,M+1,\Delta t=T/M$. We take $M$ and $N$ to be even. We denote by $\varphi_{ij}=\varphi(x_{1}^{i},x_{2}^{j}),i,j=\overline{1,N+1}$.

In order to approximate the solution of the parabolic system from Step 1 we use a finite difference method, an implicit one, descending with respect to time levels.

We denote by $h=\Delta x_{1}=\Delta x_{2}$, $p_{ij}^{k}=p(x_{1}^{i},x_{2}^{j},t^{k})$, $a_{ij}=a(x_{1}^{i},x_{2}^{j})$, $G_{ij}^{k}=\Delta tLH_{\varepsilon}(\varphi_{ij})(1+p_{ij}^{k+1})H_{\varepsilon}(1+p_{ij}^{k+1})$, $k=\overline{1,M+1}$. The numerical scheme is

[TABLE]

We take the diffusion coefficient $d=1$ for the implementation, and denote by $\lambda=\Delta t/h^{2}.$ For the interior nodes we get

[TABLE]

for $i,j=\overline{2,N},k=M,\dots,1$.

by using the Neumann conditions on the boundary, the numerical scheme becomes

[TABLE]

We denote by

[TABLE]

the vector formed by the values of $p$ at time level $k$ for the interior nodes. This is a vector of dimension $(N-1)^{2}$. We also use the following notations $P=p_{ij}^{k+1}$, $G=\Delta tG_{ij}^{k+1}$, $E_{1}=\Delta t(-a_{ij})$, $E_{2}=1+2\lambda+E_{1}$, $E_{3}=1+3\lambda+E_{1}$, $E_{4}=1+4\lambda+E_{1}$. This quantities must be evaluated at each time step $k=M,M-1,\dots,1$ and for all $i,j=\overline{2,N}$. The algebraic linear system to solve at each time step $k=M,M-1,\dots,1$ is of the form $Ax^{k}=B$, with the system matrix $A$ of dimension $(N-1)^{2}\times(N-1)^{2}$ and the vector of constant terms $B$ of dimension $(N-1)^{2}$. Based on (11) and using also the final condition, for each time level $k=M,M-1,\dots,1$ we generate the matrix $A$ and the vector $B$ with the following algorithm: we denote by $q$ the row index of matrix $A$; at the beginning of each time iteration we make the initializations: $q=0$, $A=0_{(N-1)^{2}\times(N-1)^{2}}$, and $B=0_{(N-1)^{2}\times 1}$. Then, for $i$ from $2$ to $N$ and for $j$ from $2$ to $N$, after the evaluation of $E_{1},E_{2},E_{3},E_{4},G,P$, we start the construction of $A$ and $B$. The index $q$ is incremented for each $i$ and $j$. Therefore,

• •

if i = 2 and j = 2 then q = q + 1; A(q,1) = $E_{2}$; A(q,2) = -$\lambda$;

A(q,N) = -$\lambda$; B(q) = P - G;

• •

if i = 2 and 2 < j < N then q = q + 1; A(q,j-2) = -$\lambda$;

A(q,j-1) = $E_{3}$; A(q,N+j-2) = -$\lambda$; A(q,j) = -$\lambda$; B(q) = P - G;

• •

if i = 2 and j = N then q = q + 1; A(q,N-2) = -$\lambda$;

A(q,N-1) = $E_{2}$; A(q,2*N-2) = -$\lambda$; B(q) = P - G;

• •

if 2 < i < N and j = 2 then q = q + 1;

A(q,(i-3)(N-1)+1) = -$\lambda$; A(q,(i-2)(N-1)+1) = $E_{3}$;

A(q,(i-1)(N-1)+1) = -$\lambda$; A(q,(i-2)(N-1)+2)= -$\lambda$; B(q) = P - G;

• •

if 2 < i < N and 2 < j< N then q = q + 1;

A(q,(i-3)(N-1)+j-1) = -$\lambda$; A(q,(i-1)(N-1)+j-1) = -$\lambda$

; A(q,(i-2)(N-1)+j-1) = $E_{4}$; A(q,(i-2)(N-1)+j-2) = -$\lambda$;

A(q,(i-2)*(N-1)+j) = -$\lambda$; B(q) = P - G;

• •

if 2 < i < N and j = N then q = q + 1; A(q,(i-2)*(N-1)) = -$\lambda$;

A(q,(i-2)(N-1)+N-2) = -$\lambda$; A(q,(i-1)(N-1)) = $E_{3}$;

A(q,i*(N-1)) = -$\lambda$; B(q) =P - G;

• •

if i = N and j = 2 then q = q + 1; A(q,(N-3)*(N-1)+1) = -$\lambda$;

A(q,(N-2)(N-1)+1) = $E_{2}$; A(q,(N-2)(N-1)+2) = -$\lambda$; B(q) = P - G;

• •

if i = N and 2 < j < N then q = q + 1;

A(q,(N-3)(N-1)+j-1) = -$\lambda$; A(q,(N-2)(N-1)+j-2) = -$\lambda$;

A(q,(N-2)(N-1)+j-1) = $E_{3}$; A(q,(N-2)(N-1)+j) = -$\lambda$; B(q) = P - G;

• •

if i = N and j = N then q = q + 1; A(q,N*(N-2)) = -$\lambda$;

A(q,(N-2)(N-1)) = -$\lambda$; A(q,(N-1)(N-1)) = $E_{2}$; B(q) = P - G;

Then, the resulting algebraic linear system is solved by Gaussian elimination. The solution obtained is a vector $D$ of dimension $(N-1)^{2}$. Therefore, we get the corresponding solution $p_{ij}$ at time step $k$ by the process:

q = 0;

for i = 2 to N

for j = 2 to N

q = q + 1; $p_{ij}^{k}$ = D(q);

By using the boundary condition, the solution is completed for $i=1$, $i=N+1$ and $j=\overline{1,N+1}$ and for $j=1$, $j=N+1$ and $i=\overline{1,N+1}$. Now we have the complete solution $p_{ij}^{k}$ and we can proceed with the time step $k-1$.

The integrals from Step 1 are numerical computed using Simpson’s method corresponding to the discrete grid. For each iteration $n=1,2,3,\dots$, we have to evaluate the first integral

[TABLE]

where

[TABLE]

In order to approximate this integral we first calculate, for all $i=\overline{1,N+1}$,

[TABLE]

and then

[TABLE]

To numerical evaluate of the second integral we must approximate $|\nabla\varphi(x_{1}^{i},x_{2}^{j})|$. In order to do this, we use central difference both in $x_{1}$ and in $x_{2}$ direction.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The parabolic system from Step 3 is approximated also using a finite difference method, but now ascending with respect to time levels. For each iteration $n=1,2,\dots$ and for each time level $k=1,2,3,\dots,M$, the matrix of the resulting algebraic system is the same as matrix A previously determinated, with $B(q)=p_{ij}^{k},q=\overline{1,(N-1)^{2}}$, and $E_{1}=\Delta t(-a_{ij}+G_{ij})$, $G_{ij}=LH_{\varepsilon}(\varphi_{ij})H_{\varepsilon}(1+p_{ij}^{k+1})+LH_{\varepsilon}(\varphi_{ij})(1+p_{ij}^{k+1})\delta_{\varepsilon}(1+p_{ij}^{k+1})$, which are evaluated for each $i,j=\overline{2,N}$. The resulting algebraic linear system is solved by Gaussian elimination. By using the boundary conditions we complete the solution of the parabolic system for each time level.

Numerical examples

We consider a normal initial population density $y_{0}(x_{1},x_{2})=\frac{1}{2\pi}e^{-\frac{x_{1}^{2}+x_{2}^{2}}{2}}$, where $(x_{1},x_{2})\in\Omega$. Let the diffusion coefficient be $d=1$, the final time $T=1$, $L=1$, and the regularization parameter $\varepsilon=1$. We take the space discretization step and the time discretization step to be equal $\Delta x_{1}=\Delta x_{2}=\Delta t=0.05$. For the convergence tests we consider $\varepsilon_{1}=\varepsilon_{2}=0.001$.

In the following figure, the white area represents the subregion $\omega$ that provides a small value for $J$.

Test 1. We take the natural growth rate of the population to be a constant, e.g. $a(x_{1},x_{2})=3$, $(x_{1},x_{2})\in\Omega$. The initialization of $\varphi$ is made by $\varphi^{(0)}(x_{1},x_{2})=0.25-\sqrt{(x_{1}-0.5)^{2}+(x_{2}-0.5)^{2}}$, $(x_{1},x_{2})\in\Omega$. We penalize the length of $\partial\omega$ by $\alpha=0.4$ and the area of $\omega$ by $\beta=0.6$. The corresponding results are shown in Figure 1.

Test 2. We use the same input data from Test 1 and the initialization of $\varphi$ with $\varphi^{(0)}(x_{1},x_{2})=sin(3\pi x_{1})sin(3\pi x_{2})$, $(x_{1},x_{2})\in\Omega$, a function that produce a initial checkerboard shape. We penalize the length of $\partial\omega$ by $\alpha=0.5$ and the area of $\omega$ by $\beta=0.5$. The results are shown in Figure 2.

3 Eradicating an age-structured pest population with diffusion

Consider here an age-structured population dynamics with diffusion and logistic term:

[TABLE]

Here $A\in(0,+\infty)$ is the maximal age for the population species and $y(x,a,t)$ is the population density at position $x$, age $a$ and time $t$; $d\in(0,+\infty)$ is the diffusion coefficient, $\mu(a)$ is the mortality rate and $\beta(a)$ is the fertility rate for individuals of age $a$; $y_{0}(x,a)$ is the initial density of population at position $x$ and age $a$. $u(x,t)$ is a harvesting effort (the control) and is localized in the subregion $\omega$; $u$ does not depend on age.

Assume that $\Omega$, $\omega$ satisfy the same assumptions as in the introduction, and that the following hypotheses are satisfied as well

• (H1’)

$\beta\in C([0,A])$, $\beta(a)\geq 0,\forall a\in[0,A]$;

• (H2’)

$\mu\in C([0,A))$, $\mu(a)\geq 0,\forall a\in[0,A]$, $\int_{0}^{A}\mu(a)da=+\infty$;

• (H3’)

$y_{0}\in L^{\infty}(\Omega\times(0,A))$, $y_{0}(x,a)\geq 0$ a.e. in $\Omega\times(0,A)$;

• (H4’)

$\mathcal{M}:[0,+\infty)\longrightarrow[0,+\infty)$ is continuously differentiable, $\mathcal{M}^{\prime}(r)>0,\forall r>0,\mathcal{M}(0)=0,\lim\limits_{r\to+\infty}\mathcal{M}(r)=+\infty$.

For any $u\in L_{loc}^{\infty}(\overline{\omega}\times[0,+\infty))$, such that $L\geq u(x,t)\geq 0$ a.e., there exists a unique solution $y^{u}$ to (12) (here $L\in(0,+\infty)$ is a constant; $L$ is the maximal affordable effort). This solution is nonnegative (see [4]).

Our goal is to eradicate this population which is considered to be a pest population.

Definition 1**.**

We say that the population is eradicable (zero-stabilizable) if for any $y_{0}$ satisfying the hypothesis (H3’) there exists $u\in L_{loc}^{\infty}(\overline{\omega}\times[0,+\infty))$, satisfying $L\geq u(x,t)\geq 0$ a.e., such that

[TABLE]

(and $y^{u}(x,a,t)\geq 0$ a.e. in $\Omega\times(0,A)\times(0,+\infty)$).

Note that this is a problem of zero-stabilization with control and state constraints.

Denote by $r^{*}\in\mathbb{R}$ the solution to the equation

[TABLE]

and $\lambda_{1}^{\omega}$ the principal eigenvalue for

[TABLE]

Theorem 4**.**

(i) If the population is eradicable then

[TABLE]

(ii) If $\lambda_{1}^{\omega}>r^{*}$ then the population is eradicable and the harvesting effort $u\equiv L$ diminishes exponentially the population.

Proof.

(i) Assume that the population is eradicable and let

[TABLE]

with $h_{0}\in C([0,A]),h_{0}(a)>0,\forall a\in[0,A]$, to be specified later, $g_{0}\in L^{\infty}(\Omega),g_{0}(x)>0$ a.e in $\Omega$.

Let $u\in L_{loc}^{\infty}(\overline{\omega}\times[0,+\infty))$, $L\geq u(x,t)\geq 0$ a.e., such that

[TABLE]

The unique solution $y^{u}$ to (12) may be written as

[TABLE]

where $h$ is the solution to

[TABLE]

and $g$ is the solution to

[TABLE]

It is known that the set of solutions for $\eqref{4.3}_{1-2}$ is a real vector space of dimension $1$ and there exists a time independent solution $\tilde{h}$ satisfying $\tilde{h}(a)>0,$ for all $a\in[0,A)$ ([2]).

If we consider $h_{0}=\tilde{h}$, then

[TABLE]

where $g$ is the solution to

[TABLE]

Here $H=\int_{0}^{A}\tilde{h}(a)da$.

The eradicability for (12) implies the nonnegative zero-stabilizability for (15). However, the nonnegative zero-stabilizability for (15) implies that

[TABLE]

This follows as in [1] by using of the comparison results for the solutions to parabolic equations.

(ii) If $\lambda_{1}^{\omega}>r^{*}$, then we consider $u(x,t)=L$ a.e. in $\omega\times(0,+\infty)$. Using the comparison result for linear age-structured population dynamics (see [2]) we get that

[TABLE]

where $\tilde{y}$ is the solution to

[TABLE]

Let $h_{0}(a)=1,\forall a\in[0,A]$, $g_{0}(x)=||y_{0}||_{\infty}$ a.e. $x\in\Omega$.

Using again the comparison result for linear age-structured population dynamics we get that

[TABLE]

where $h$ is the solution to (13) and $g$ is the solution to (14) corresponding to $\mathcal{M}\equiv 0$ and $u\equiv L$. Since $h(\cdot,t)\rightarrow\bar{h}$ in $L^{\infty}(0,A)$ as $t\rightarrow+\infty$ ([2]), and $g(\cdot,t)\rightarrow 0$ in $L^{\infty}(\Omega)$ as $t\rightarrow+\infty$ (because $\lambda_{1}^{\omega}>r^{*}$), we get by (16) and (17) that

[TABLE]

and the conclusion. ∎

Since our goal was actually to eradicate a pest population corresponding to a initial density $y_{0}$ with a harvesting effort less or equal than $L$ (tacking into account the above theorem) and since we have however to pay a certain cost to harvest in a subdomain $\omega$, we can consider the following related optimal control problem

[TABLE]

where $T>0$ is a certain moment and $y$ is the solution to (12) corresponding to $u\equiv L$.

This problem may be investigated by using the level set method described in Section 2 and rewriting it in the following form

[TABLE]

where $y_{\varphi}$ is the solution to

[TABLE]

with $\varphi$ the implicit function of $\omega$. The approach is similar to the one in Section 2. We will approximate this problem using the mollified version of the Heaviside function, $H_{\varepsilon}$, and its derivative by the mollified function $\delta_{\varepsilon}$.

Actually, if we denote by

[TABLE]

for a small but fixed $\varepsilon>0$, the harvesting problem to be investigated is

[TABLE]

where $\varphi:\overline{\Omega}\longrightarrow\mathbb{R}$ is a smooth function and $y_{\varphi}$ is the solution to

[TABLE]

By following the same lines as in Section 2 we can get the directional derivative of $\Psi$. We reach a similar conclusion as in Section 2 concerning the gradient descent with respect to $\varphi$:

[TABLE]

where $\theta$ is an artificial time and $r$ is solution to

[TABLE]

Acknowledgements

This work was supported by the CNCS-UEFISCDI (Romanian National Authority for Scientific Research) grant 68/2.09.2013, PN-II-ID-PCE-2012-4-0270: “Optimal Control and Stabilization of Nonlinear Parabolic Systems with State Constraints. Applications in Life Sciences and Economics”.