Regional control in optimal harvesting of population dynamics

Nonlinear Analysis 147 (2016) 191–212

Contents lists available at ScienceDirect

Nonlinear Analysis www.elsevier.com/locate/na

Regional control in optimal harvesting of population dynamics Sebastian Anit¸a a,b , Vincenzo Capasso c,d,∗ , Ana-Maria Mo¸sneagu a a

Faculty of Mathematics, “Alexandru Ioan Cuza” University of Ia¸ si, Ia¸ si 700506, Romania “Octav Mayer” Institute of Mathematics of the Romanian Academy, Ia¸ si 700506, Romania c ADAMSS (Advanced Applied Mathematical and Statistical Sciences), Universit´ a degli Studi di Milano, 20133 Milano, Italy d Department of Mathematics, Universit´ a degli Studi di Milano, 20133 Milano, Italy b

article

info

Article history: Received 9 June 2016 Accepted 13 September 2016 Communicated by Enzo Mitidieri Keywords: Optimal harvesting Regional control Optimality conditions Population dynamics Space/age structure Iterative algorithm

abstract Here we investigate the regional control for some optimal harvesting problems related to population dynamics; namely we consider the problem of maximizing the profit for spatially structured harvesting problems with respect to both the harvesting effort and the selected subregion ω (of the whole domain Ω ) where the effort is localized. For a fixed subregion ω we state necessary optimality conditions and use them to get the structure of the optimal effort and to reformulate the maximization problem with respect to the subregion ω, where the harvesting effort is localized, in a more convenient way. We derive an iterative algorithm to increase at each iteration the profit by changing the subregion where the effort is localized. Some numerical tests are given to illustrate the effectiveness of the results for a particular optimal harvesting problem. Final comments are given as well concerning further directions to extend the results and methods presented here. © 2016 Elsevier Ltd. All rights reserved.

1. Introduction and setting of the problem Optimal harvesting problems have been intensively studied for structured population dynamics by several authors. We wish further to evidence that possible space (and time) heterogeneities of the domain Ω are taken in due account via the parameters of the relevant dynamical system. The main attention has been paid in deriving first order necessary optimality conditions and to use gradient-type algorithms to approximate an optimal effort (control). Spatially structured harvesting problems have usually been treated by considering an effort applied to the whole domain Ω of interest for the relevant population (see for example [2]). Actually it should be obvious that, on the basis of suitable information, the effort is localized in a subregion ω of Ω (see [5]). ∗ Corresponding author. E-mail addresses: [email protected] (S. Anit¸a), [email protected] (V. Capasso), [email protected] (A.-M. Mo¸sneagu).

http://dx.doi.org/10.1016/j.na.2016.09.008 0362-546X/© 2016 Elsevier Ltd. All rights reserved.

192

S. Anit¸a et al. / Nonlinear Analysis 147 (2016) 191–212

Hence the point we wish to raise in this paper regards not only the optimization with respect to the intensity of the harvesting effort, but also the optimal choice of the subregion ω where it is implemented. We wish to evidence that such problems may be treated as shape optimization problems. Hence, to begin with, we remind that from a measure theoretic point of view the geometry of a set ω can be characterized in terms of its Minkowski functionals, or intrinsic volumes. In 2D there are three such functionals and these are proportional to more commonly known quantities such as the area, the perimeter and the Euler–Poincar´e characteristic; see [26], [31, p. 30]. Intrinsic volumes have been successfully used in material science to characterize and discriminate morphology of various media; see [38,7]. In this paper we shall limit ourselves to refer to two of such functionals, the area, and the perimeter. On the other hand it is of course easily understandable, from the point of view of our optimal control problems that, in order to optimize the harvesting region, we have to act on both the area and the perimeter. A convenient tool for handling the shape of a domain ω in R2 is the so called implicit interface representation, according to which the boundary of a domain is defined as the isocontour of some function ϕ (see [37, p. 3], [19, p. 75]). Assume that in the domain Ω ⊂ R2 we may represent a subset ω by means of a function ϕ : Ω → R called its implicit function as ω = {x ∈ Ω ; ϕ(x) ≥ 0} (or ω = {x ∈ Ω ; ϕ(x) ≤ 0}), while its boundary is defined as the zero level set of ϕ: ∂ω = {x ∈ Ω ; ϕ(x) = 0}. In this way the interior of ω has the following implicit representation Int(ω) = {x ∈ Ω ; ϕ(x) > 0}. A useful tool to carry out integration over the interior of the domain ω is then the Heaviside function H : R → {0, 1}, such that  1, if z ≥ 0 H(z) = 0, if z < 0. If ϕ is the implicit function of ω, in order to integrate over ω a function f defined over the whole Ω we may  write Ω f (x)H(ϕ(x))dx. If we need to integrate over the boundary of ω we may then refer to the Dirac Delta function as the generalized derivative of the Heaviside function H. In fact, if ϕ is sufficiently smooth, the directional derivative of the Heaviside function in the normal direction at a point x ∈ ∂ω is given by   δ(x) = H ′ (ϕ(x))|∇ϕ(x)|, and, by using the usual Dirac Delta δ on R, we have δ(x) = δ(ϕ(x))|∇ϕ(x)|. Hence the boundary integral over ∂ω of a function f defined over the whole Ω may be expressed as  f (x)δ(ϕ(x))|∇ϕ(x)|dx. In the above it has been taken into account that, in terms of the implicit function Ω ∇ϕ(x) . ϕ, the unit (outward) normal to the boundary at point x ∈ ∂ω can be expressed as n(x) = |∇ϕ(x)| We wish to mention that one might rephrase the above control problem in terms of suitable measures, as done in [12,17,18]; indeed they have treated the problem of existence for a generic optimal harvesting effort, with constraints. Our approach has allowed us to take advantage of the meaning and of the specific structure of our system. We denote by y u (x, t) the population density at position x and time t of a population that is free to move in a habitat Ω ⊂ R2 and is subject to a harvesting effort u(x, t) (localized in a subregion ω, with t ∈ [0, T ]). Our main purpose is to investigate the following problem: Maximize Maximize Φω (u), ω

u∈Uω

(OH)

T  where Φω (u) = 0 ω [u(x, t)y u (x, t) − cu(x, t)]dx dt − c1 length(∂ω) − c2 area(ω), and Uω = {v ∈ L∞ (ω × (0, T )); 0 ≤ v(x, t) ≤ L a.e. in ω × (0, T )}, T, L, c1 , c2 are positive constants and c ∈ [0, +∞) (c represents T  the cost per unit of the harvesting effort, 0 ω cu(x, t)dx dt is the cost of the total harvesting effort, and c1 length(∂ω) + c2 area(ω) is the cost paid to harvest in the specific subregion ω). It means that we wish to maximize (with respect to the subregion ω and to the harvesting effort) the profit, which is the difference between the gain of the harvested population and the cost of the harvesting effort plus the cost paid to harvest in the subregion ω.

S. Anit¸a et al. / Nonlinear Analysis 147 (2016) 191–212

193

In view of what we have said above about the implicit interface representation, the problem may be rewritten in terms of the function ϕ, which identifies (apart from obvious multiplicities regarding the implicit representation), the subregion ω as the set ω = {x ∈ Ω ; ϕ(x) > 0}, Maximize Maximize Ψϕ (u), ϕ

u∈Uω

T    where Ψϕ (u) = 0 ω (uy u − cu)dx dt − c1 Ω δ(ϕ(x))|∇ϕ(x)|dx − c2 Ω H(ϕ(x))dx. From now on we will refer either to the subregion ω or to the implicit function ϕ, indifferently. We shall treat such problems in two steps: firstly, for a fixed ω (or ϕ) we derive first order optimality conditions and go further using them to get the structure of the optimal harvesting effort (control). Next we investigate the second problem by adapting some shape optimization techniques. Note that some particular harvesting problems related to structured population dynamics have been studied by several authors; see [3,2,4,9–14,20,24,25,27,29,34–36,44,45]. For optimal control related to biological models see [33], and references therein (for an updated introduction and literature one may refer to [32]). For basic results and methods in shape optimization we refer to [15,19,28,37,41,42]. Our present paper is organized as follows. Sections 2–4 are devoted to some optimal harvesting problems when the subregion ω (where the harvesting effort is localized) is given: Section 2 concerns a spatially structured population with a nonlocal logistic term; Section 3 investigates the dynamics of an age-and space-structured population with a logistic term which is local with respect to space and nonlocal with respect to age; Section 4 concerns a space-structured population model with a local logistic term. The model in Section 4 is a simpler one. We wish to emphasize here that under quite general assumptions the optimal control is unique. In Section 5 we derive an iterative algorithm to increase at each iteration the profit by changing the subregion where the effort is localized. Some numerical tests to illustrate the effectiveness of the results are given for a particular optimal harvesting problem. Final comments concerning possible further directions to extend the results and methods presented here are given in Section 6. 2. Optimal harvesting for a spatially structured population with diffusion and nonlocal logistic term The following model describes the dynamics of a single population species which is free to move in the habitat Ω ⊂ RN , N ∈ {2, 3} (Ω is a bounded domain with a sufficiently smooth boundary):     ∂ y(x, t) − d∆y(x, t) = η(x)y(x, t) − y(x, t) k(x, x′ )y(x′ , t)dx′ t    Ω −χω (x)u(x, t)y(x, t), (x, t) ∈ Q (1)   ∂ y(x, t) = 0, (x, t) ∈ Σ  ν   y(x, 0) = y0 (x), x ∈ Ω, where T ∈ (0, +∞), Q = Ω × (0, T ), Σ = ∂Ω × (0, T ). Here y(x, t) is the population density at position x ∈ Ω and moment t ∈ [0, T ]; d ∈ (0, +∞) is a diffusion coefficient, η(x) is the growth rate at position x (the difference between the fertility rate and the mortality rate). The second order logistic term  y(x, t) Ω k(x, x′ )y(x′ , t)dx′ is nonlocal; it represents an additional mortality term due to the competition for resources with other individuals all over the habitat, via a kernel k. For related biological models with diffusion and nonlocal terms we refer to [22,4]. The homogeneous Neumann boundary condition shows that there is no population flow through the boundary of the habitat; y0 (x) is the initial population density at position x. In system (1), u(x, t) represents the harvesting effort (control) that acts in the subdomain ω ⊂ Ω (χω is the characteristic function of ω). Hence the term χω (x)u(x, t)y(x, t) describes the action of the harvesting at location x ∈ ω, and time t on the evolution of the population density y(x, t).

194

S. Anit¸a et al. / Nonlinear Analysis 147 (2016) 191–212

The first question to be investigated is the following optimal harvesting problem:  T [u(x, t)y u (x, t) − cu(x, t)]dx dt, Maximize γ u∈U

0

(OH1)

ω

where U = {v ∈ L∞ (ω × (0, T )); 0 ≤ v(x, t) ≤ L a.e. in ω × (0, T )} (L, γ ∈ (0, +∞), c ∈ [0, +∞) are constants) and y u is the solution to (1). Actually it means that we are interested to maximize the profit (that is the difference between the gain from the harvested population and the cost of the harvesting effort). The constant γ represents the cost per unit for the harvested population. The cost of the harvested population may be proportional to the total harvested population, and the T  total profit has the form γ˜ 0 ω u(x, t)y u (x, t)dx dt (and corresponds to c ≡ 0), or may be proportional T  to the total harvesting effort, and the total profit has the form γ 0 ω [u(x, t)y u (x, t) − cu(x, t)]dx dt (and corresponds to the cost per unit of harvesting effort γc ∈ (0, +∞)). For the sake of simplicity we assume that γ = 1. Here are the hypotheses we are going to use throughout this section. (H1) η, y0 ∈ L∞ (Ω ), y0 (x) ≥ 0 a.e. in Ω and ∥y0 ∥L∞ (Ω) > 0; (H2) k ∈ L∞ (Ω × Ω ), k(x, x′ ) ≥ 0 a.e. in Ω × Ω . Consider the following parabolic equation with boundary and   ∂t w(x, t) − d∆w(x, t) = f (x, t), ∂ν w(x, t) = 0,  w(x, 0) = w (x), 0

initial conditions (x, t) ∈ Q (x, t) ∈ Σ x ∈ Ω.

(2)

We have to be precise about what we do mean by a solution to problem (2). We say that w ∈ C([0, T ]; L2 (Ω )) ∩ AC((0, T ]; L2 (Ω )) ∩ L2 (0, T ; H 1 (Ω )) ∩ L2loc ((0, T ]; H 2 (Ω )) is a (weak) solution to (2) if satisfies   ∂t w(x, t) − d∆w(x, t) = f (x, t), a.e. in Q ∂ν w(x, t) = 0, a.e. in Σ  w(x, 0) = w (x), a.e. in Ω . 0 It is known that for any w0 ∈ L2 (Ω ) and f ∈ L2 (Q), problem (2) has a unique (weak) solution (see [2, p. 187], and [10, p. 227]). Moreover, we have that     t  t 1 |w(x, t)|2 dx + d |∇w(x, s)|2 dx ds ≤ C |w0 (x)|2 dx + |f (x, s)|2 dx ds , 2 Ω 0 Ω Ω 0 Ω for any t ∈ (0, T ], where C ∈ (0, +∞) is a constant. Actually, Problem (1) is a particular case of the following nonlinear problem   ∂t w(x, t) − d∆w(x, t) = F (w)(x, t), (x, t) ∈ Q, ∂ν w(x, t) = 0, (x, t) ∈ Σ ,  w(x, 0) = w (x), x ∈ Ω. 0

(3)

We say that w is a (weak) solution to (3) if w ∈ C([0, T ]; L2 (Ω )) ∩ AC((0, T ]; L2 (Ω )) ∩ L2 (0, T ; H 1 (Ω )) ∩ L2loc ((0, T ]; H 2 (Ω )) and satisfies   ∂t w(x, t) − d∆w(x, t) = F (w)(x, t), a.e. in Q ∂ν w(x, t) = 0, a.e. in Σ  w(x, 0) = w (x), a.e. in Ω . 0

S. Anit¸a et al. / Nonlinear Analysis 147 (2016) 191–212

195

By Banach’s fixed point theorem (in C([0, T ]; L2 (Ω ))) we get that for any u ∈ U, there exists a unique solution y u to (1). Moreover, we may remark that, as far as comparison theorems are concerned, the first equation in (1) may be written in the form ∂t y(x, t) = d∆y(x, t) + ρ(x, t)y(x, t),

(x, t) ∈ Q,

for ρ(x, t) := η(x) − Ω k(x, x′ )y(x′ , t)dx′ − χω (x)u(x, t) ∈ L∞ (Q); we then obtain that 0 < y u (x, t), a.e. (x, t) ∈ Q, via the comparison result for parabolic equations (see [21,40,5] for parabolic equations  with a nonlocal term) and using the assumptions on y0 . On the other hand, since Ω k(x, x′ )y(x′ , t)dx′ ≥ 0 a.e. (x, t) ∈ Q, we conclude, via the same comparison result, that 

0 < y u (x, t) ≤ y˜(x, t)

a.e. (x, t) ∈ Q,

where y˜ is the solution to (1) corresponding to u ≡ 0 and k ≡ 0. Actually y˜ ∈ L∞ (Q) (˜ y (x, t) ≤ ∥y0 ∥∞ exp{t∥η∥∞ } a.e. in Q) and so y u ∈ L∞ (Q) as well. Theorem 1. There exists at least one optimal control u∗ for (OH1). Proof. Let T





Φ(u) = 0

[u(x, t)y u (x, t) − cu(x, t)]dx dt,

ω

and R = Sup Φ(u). Since for any u ∈ U, we have that −cLT ≤ Φ(u) ≤ L u∈U

T  0

ω

y˜dx dt, we may infer that

R ∈ R. It follows that there exists {un }n∈N∗ ⊂ U such that 1 < Φ(un ) ≤ R, n

R− for any n ∈ N∗ . Since, the sequence given by



fn (x, t) = η(x)y un (x, t) − y un (x, t)

k(x, x′ )y un (x′ , t)dx − χω (x)un (x, t)y un (x, t),

Ω

(x, t) ∈ Q, is bounded in L∞ (Q), and {un } is bounded in L2 (ω × (0, T )), we get that on a subsequence unq ⇀ u∗ ,

fnq ⇀ f

in L2 (Q).

The set U is convex and closed in L2 (ω × (0, T )), and so u∗ ∈ U. The last convergence and the compactness (in L2 (Ω )) of the C0 -semigroup generated by the operator A, defined by Ay = d∆y,

D(A) = {v ∈ H 2 (Ω ); ∂ν v(x) = 0 on ∂Ω },

implies that y unq → y

in C([0, T ]; L2 (Ω ))

(see [8,39,4]), where y is the solution to   ∂t y(x, t) − d∆y(x, t) = f (x, t), (x, t) ∈ Q ∂ν y(x, t) = 0, (x, t) ∈ Σ  y(x, 0) = y (x), x ∈ Ω. 0 On the other hand,  fnq ⇀ ηy − y Ω

k(·, x′ )y(x′ , t)dx′ − χω u∗ y

S. Anit¸a et al. / Nonlinear Analysis 147 (2016) 191–212

196

in L2 (Q). Hence, y is the unique solution to (1) with u := u∗ . We immediately conclude that Φ(unq ) → Φ(u∗ ), and so Φ(u∗ ) = R. As far as the first order necessary optimality conditions are concerned, we may prove the following result. Theorem 2. If u∗ is an optimal control (effort) for (OH1) and if p is the solution   ∗  ∂ p(x, t) + d∆p(x, t) = −η(x)p(x, t) + p(x, t) k(x, x′ )y u (x′ , t)dx′  t     Ω   ′ u∗ ′ ′ ′ + k(x , x)y (x , t)p(x , t)dx + χω (x)u∗ (x, t)(1 + p(x, t)), Ω     ∂ν p(x, t) = 0,    p(x, T ) = 0,

to

(x, t) ∈ Q

(4)

(x, t) ∈ Σ x ∈ Ω,

then  ∗

u (x, t) =

∗

∗

0, if y u (x, t)p(x, t) + y u (x, t) − c < 0 ∗ ∗ L, if y u (x, t)p(x, t) + y u (x, t) − c > 0.

(5)

We have to notice that (4) has a unique solution p, meaning that p ∈ C([0, T ]; L2 (Ω ))∩AC([0, T ); L2 (Ω ))∩ L (0, T ; H 1 (Ω )) ∩ L2loc ([0, T ); H 2 (Ω )) and satisfies   ∗  ∂t p(x, t) + d∆p(x, t) = −η(x)p(x, t) + p(x, t) k(x, x′ )y u (x′ , t)dx′     Ω   ′ u∗ ′ ′ ′ + k(x , x)y (x , t)p(x , t)dx + χω (x)u∗ (x, t)(1 + p(x, t)), a.e. in Q  Ω    ∂ p(x, t) = 0, a.e. in Σ   ν p(x, T ) = 0, a.e. in Ω . 2

Taking into account that by the substitution T − t := s, problem (4) may be included in the general framework (3) and that the existence and uniqueness of the solution to (4) follows via Banach’s fixed point theorem (in C([0, T ]; L2 (Ω ))). For the proof of Theorem 2 the following auxiliary result shall be used. Lemma 3. Let v be arbitrary such that u∗ + εv ∈ U for sufficiently small ε > 0 ( v ∈ L∞ (ω × (0, T )) and 0 ≤ u∗ (x, t) + εv(x, t) ≤ L a.e. in ω × (0, T )). Then yu

∗

+εv

→ yu

∗

in L∞ (Q)

(i)

and ∗ 1 u∗ +εv [y − y u ] → z in L∞ (Q), ε as ε → 0+, where z is the solution to   ∗  ∂t z(x, t) − d∆z(x, t) = η(x)z(x, t) − z(x, t) k(x, x′ )y u (x′ , t)dx′     Ω   ∗   − y u (x, t) k(x, x′ )z(x′ , t)dx′

(ii)

Ω

∗  − χω (x)[u∗ (x, t)z(x, t) + v(x, t)y u (x, t)],       ∂ν z(x, t) = 0, z(x, 0) = 0,

(x, t) ∈ Q (x, t) ∈ Σ x ∈ Ω.

(6)

It is obvious that via Banach’s fixed point theorem problem (6) has a unique solution (notice that (6) may be included in the framework of (3)).

S. Anit¸a et al. / Nonlinear Analysis 147 (2016) 191–212

Proof of Lemma 3. Let wε = y u

∗

197

∗

− y u ; wε is the solution to   ∗  ∂t w(x, t) − d∆w(x, t) = η(x)w(x, t) − w(x, t) k(x, x′ )y u (x′ , t)dx′    Ω    ∗   − y u +εv (x, t) k(x, x′ )w(x′ , t)dx′ +εv

Ω

∗  − χω (x)[u∗ (x, t)w(x, t) + εv(x, t)y u +εv (x, t)],       ∂ν w(x, t) = 0, w(x, 0) = 0,

(x, t) ∈ Q (x, t) ∈ Σ x∈Ω

(7)

(this system admits a unique solution via Banach’s fixed point theorem). We get that 0 ≤ y u (x, t) ≤ M

a.e. in Q,

where M = ess sup y˜(x, t). If we multiply the first equation in (7) by wε and integrate over Ω × (0, t) then, (x,t)∈Q

after some calculation, we get   t 1 |wε (x, t)|2 dx ≤ ∥η∥∞ |wε (x, s)|2 dx ds 2 Ω 0 Ω  t  t + M ∥k∥∞ meas(Ω ) |wε (x, s)|2 dx ds + ε∥v∥∞ M |wε (x, s)|dx ds 0 Ω 0 Ω  t ε∥v∥∞ M T meas(Ω ) ˜ , ≤M |wε (x, s)|2 dx ds + 2 0 Ω ˜ = ∥η∥∞ + M ∥k∥∞ meas(Ω ) + ∥v∥∞ M (if ε < 1). By Bellman’s lemma we get that where M 2  ˜ t}, |wε (x, t)|2 dx ≤ εLM T meas(Ω ) exp{2M Ω

for any t ∈ [0, T ], and consequently wε → 0

in C([0, T ]; L2 (Ω )),

as ε → 0+. The last convergence implies that the right hand side term in (7) may be written as ξ(x, t)w(x, t) + ζε (x, t), where ξ, ζε ∈ L∞ (Q) for any ε > 0 and ζε → 0 in L∞ (Q), as ε → 0. By a comparison result for parabolic equation we get that  t |w(x, t)| ≤ e∥ξ∥∞ (t−s) ∥ζε ∥∞ ds 0 ∞

a.e. in Q, and consequently wε → 0 in L (Q), as ε → 0+, and so (i) holds. If we denote by qε = 1ε wε − z, then qε is the solution to   ∗  ∂ q(x, t) − d∆q(x, t) = η(x)q(x, t) − q(x, t) k(x, x′ )y u (x′ , t)dx′ t     Ω    u∗ +εv ′ ′ ′   −y (x, t) k(x, x )q(x , t)dx    Ω  u∗ +εv u∗ + (y (x, t) − y (x, t)) k(x, x′ )z(x′ , t)dx′   Ω    − χω (x)[u∗ (x, t)q(x, t) + v(x, t)wε (x, t)],      ∂ q(x, t) = 0,   ν q(x, 0) = 0,

(x, t) ∈ Q (x, t) ∈ Σ x ∈ Ω.

S. Anit¸a et al. / Nonlinear Analysis 147 (2016) 191–212

198

By using again the comparison result for parabolic equations and (i) we may conclude that qε → 0 in L∞ (Q), as ε → 0+, and so (ii) holds. Proof of Theorem 2. The existence and uniqueness of a solution to (4) follows via Banach’s fixed point theorem. Let now v be as in Lemma 3; then Φ(u∗ ) ≥ Φ(u∗ + εv), which implies that  T 0

∗

[(u∗ (x, t) + εv(x, t))y u

+εv

∗

(x, t) − u∗ (x, t)y u (x, t) − cεv(x, t)]dx dt ≤ 0.

ω

Consequently, we obtain that obtain that  T

T  0

ω

∗

[u∗ wεε + vy u

+εv

− cv]dx dt ≤ 0, and passing to the limit (ε → 0+) we

∗

[u∗ (x, t)z(x, t) + v(x, t)y u (x, t) − cv(x, t)]dx dt ≤ 0.

(8)

ω

0

By multiplying the first equation in (4) by z and integrating over Q we get that T1 = T2 , where  T T1 = z(x, t)[∂t p(x, t) + d∆p(x, t)]dx dt, Ω

0

and 

T





T2 =

z(x, t)[−η(x)p(x, t) + p(x, t) Ω

0



∗

k(x, x′ )y u (x′ , t)dx′

Ω ∗

k(x′ , x)y u (x′ , t)p(x′ , t)dx′ + χω (x)u∗ (x, t)(1 + p(x, t))]dx dt.

+ Ω

By integrating by parts (with respect to t) and taking into account (6) and (4) we get that   T ∗ T1 = p(x, t)[−η(x)z(x, t) + z(x, t) k(x, x′ )y u (x′ , t)dx′ 0 Ω Ω  ∗ ∗ + y u (x, t) k(x, x′ )z(x′ , t)dx′ + χω (x)(u∗ (x, t)z(x, t) + v(x, t)y u (x, t))]dx dt. Ω

T  T  ∗ Since T1 = T2 , we may conclude that 0 ω u∗ zdx dt = 0 ω vy u pdx dt and by (8) we infer that  T ∗ ∗ v(x, t)[y u (x, t)p(x, t) + y u (x, t) − c]dx dt ≤ 0, 0

ω

∗

for any v such that u + εv ∈ U for sufficiently small ε > 0. This implies (5) and concludes the proof of Theorem 2. ∗

Remark 1. If we multiply the first equation in (4) by y u and integrate over Ω × (0, T ), we obtain after some calculation (and using (1) and (4)) that   T  ∗ ∗ Φ(u∗ ) = − y0 (x)p(x, 0)dx − k(x′ , x)y u (x′ , t)p(x′ , t)y u (x, t)dx′ dx dt. Ω

0

Ω

Ω

A special attention will be paid in Sections 4 and 5 to the case c = 0, which corresponds to the case when ∗ the cost paid for harvesting is a part of the gain for the harvested population. Since y u (x, t) > 0 a.e. in Q, we get by (5) that  0, if 1 + p(x, t) < 0, u∗ (x, t) = L, if 1 + p(x, t) > 0,

S. Anit¸a et al. / Nonlinear Analysis 147 (2016) 191–212

and that p is the solution to   ∗  ∂t p(x, t) + d∆p(x, t) = −η(x)p(x, t) + p(x, t) k(x, x′ )y u (x′ , t)dx′     Ω   ∗ + k(x′ , x)y u (x′ , t)p(x′ , t)dx′ + χω (x)L(1 + p(x, t))+ ,  Ω   ∂ν p(x, t) = 0,   p(x, T ) = 0,

199

(x, t) ∈ Q, (x, t) ∈ Σ , x ∈ Ω.

If p(x, t) ̸= −1 a.e. in ω × (0, T ), then u∗ (x, t) = L · H(1 + p(x, t)), where H is the Heaviside function. This happens if for example k ≡ 0 and η(x) ̸= 0 a.e. in Ω , or if k(x, x′ ) = k0 δx (x′ ) with k0 a positive constant (δx (x′ ) is the Dirac mass centered at x and evaluated in x′ ) and if η(x) < 2k0 y L (x, t) a.e. in ω × (0, T ) or if η(x) > 2k0 y 0 (x, t) a.e. in ω × (0, T ) (we shall give more details in Section 4 in a slightly more general case; this is a limit case and is not however a particular case since such k does not satisfy the hypothesis (H2)). Actually, we may conclude that the maximal profit (for a given ω) is   T  Φ(u∗ ) = − y0 (x)p(x, 0)dx − k(x′ , x)y(x′ , t)p(x′ , t)y(x, t)dx′ dx dt, Ω

0

Ω

Ω

where (y, p) is the solution to    ∂ y(x, t) − d∆y(x, t) = η(x)y(x, t) − y(x, t) k(x, x′ )y(x′ , t)dx′ t    Ω    − χω (x)LH(1 + p(x, t))y(x, t),      ∂ p(x, t) + d∆p(x, t) = −η(x)p(x, t) + p(x, t) k(x, x′ )y(x′ , t)dx′ t  Ω     + k(x′ , x)y(x′ , t)p(x′ , t)dx′ + χω (x)L(1 + p(x, t))+ ,    Ω    ∂ y(x, t) = ∂ν p(x, t) = 0,   ν y(x, 0) = y0 (x), p(x, T ) = 0,

(x, t) ∈ Q

(x, t) ∈ Q (x, t) ∈ Σ x ∈ Ω.

∗

Remark 2. Similar conclusion may be obtained in the general case if, in addition, we get that y u (x, t)p(x, t)+ ∗ y u (x, t) − c ̸= 0 a.e. in ω × (0, T ). 3. Optimal harvesting of an age-and space-structured population with diffusion and logistic term We shall investigate here an optimal harvesting problem for the following age-dependent population dynamics which describes the dynamics of a single population species which is free to move in an isolated domain Ω satisfying the assumptions in Section 2:  ∂t y(x, t, a) + ∂a y(x, t, a) − d∆y(x, t, a) = −µ(a)y(x, t, a)    A     − k y(x, t, a) y(x, t, a)da − χω (x)u(x, t, a)y(x, t, a),  0   0 ∂ν y(x, t, a) = 0,   A     β(a)y(x, t, a)da, y(x, t, 0) =    0  y(x, 0, a) = y0 (x, a),

(x, t, a) ∈ Q × (0, A) (x, t, a) ∈ Σ × (0, A)

(9)

(x, t) ∈ Q (x, a) ∈ Ω × (0, A).

Here A ∈ (0, +∞) is the maximal age for the population species and y(x, t, a) is the population density at position x, time t and age a; d ∈ (0, +∞) is the diffusion coefficient, µ(a) is the mortality rate and β(a) is the fertility rate for individuals of age a, k0 is a nonnegative constant.

S. Anit¸a et al. / Nonlinear Analysis 147 (2016) 191–212

200

A The second order logistic term taken here, k0 y(x, t, a) 0 y(x, t, a)da, is nonlocal with respect to age and local with respect to the space, and represents an additional mortality term due to the competition for resources. The homogeneous Neumann boundary conditions describe the fact that there is no population flow A through the boundary. The newborn population is given by the renewal law y(x, t, 0) = 0 β(a)y(x, t, a)da, which is a nonlocal boundary condition (see [2,30,43]). Finally y0 (x, a) is the initial population density at position x and age a; u(x, t, a) represents a harvesting effort (control) localized in ω. Here are the assumptions we are going to use throughout this section: (H1)′ y0 ∈ L∞ (Ω × (0, A)), y0 (x, a) > 0 a.e. in Ω × (0, A); (H2)′ β ∈ L∞ (0, A), β(a) ≥ 0 a.e. in (0, A), ∃(a0 , a1 ) ⊂ (0, A) such that β(a) > 0 a.e. in (a0 , a1 ); ′

(H3) µ ∈ L∞ loc (0, A), µ(a) ≥ 0 a.e. in (0, A). Take now U˜ = {v ∈ L∞ (ω × (0, T ) × (0, A)); 0 ≤ v(x, t, a) ≤ L a.e.} (L ∈ (0, +∞), c ∈ [0, +∞)). We say that y ∈ L2 (Q × (0, A)) which belongs to C(S; L2 (Ω )) ∩ AC(S; L2 (Ω )) ∩ L2 (S; H 1 (Ω )) ∩ 2 Lloc (S; H 2 (Ω )) for almost any characteristic line S of equation a − t = const., (t, a) ∈ (0, T ) × (0, A) and satisfies  A   Dy(x, t, a) − d∆y(x, t, a) = −µ(a)y(x, t, a) − k0 y(x, t, a) y(x, t, a)da    0    − χω (x)u(x, t, a)y(x, t, a), a.e. in Q × (0, A)   ∂ν y(x, t, a) = 0, a.e. in Σ × (0, A)  A     lim y(·, t + ε, ε) = β(a)y(·, t, a)da, in L2 (Ω ), a.e. in (0, T )   ε→0+  0    lim y(·, ε, a + ε) = y (·, a), in L2 (Ω ), a.e. in (0, A), 0 ε→0+

where Dh(·, t, a) denotes the directional derivative at (t, a) in direction (1, 1). Under the above mentioned ˜ (9) admits a unique solution y u hypotheses it follows as in Lemma 4.11 in [2, p. 111], that for any u ∈ U, (see [2,43]). This solution satisfies y u (x, t, a) > 0 a.e. in Q × (0, A). Using now the comparison result for linear age-dependent population dynamics (see [2,30,43]) we get that in addition 0 < y u (x, t, a) ≤ y˜(x, t, a)

a.e. in Q × (0, A),

where y˜ is the solution to (9) corresponding to u ≡ 0 and k0 = 0. The problem to be investigated in this section is 

A



T



Maximize ˜ u∈U

0

0

[u(x, t, a)y u (x, t, a) − cu(x, t, a)]dx dt da.

(OH2)

ω

Following the ideas in [1] (based on Mazur theorem) we can prove the following result: Theorem 4. There exists at least one optimal control u∗ for (OH2). AT  Let Φ(u) = 0 0 ω [u(x, t, a)y u (x, t, a)−cu(x, t, a)]dx dt da. The following result will clarify the question of the optimality conditions for (OH2):

S. Anit¸a et al. / Nonlinear Analysis 147 (2016) 191–212

201

Theorem 5. If u∗ is an optimal control (effort) for (OH2) and if p is the solution to ∂ p(x, t, a) + ∂ p(x, t, a) + d∆p(x, t, a) = µ(a)p(x, t, a) t a   A   ∗    + k p(x, t, a) y u (x, t, a)da 0    0   A    u∗ + k0 y (x, t, a)p(x, t, a)da − β(a)p(x, t, 0) 0   + χω (x)u∗ (x, t, a)(1 + p(x, t, a)),      ∂ν p(x, t, a) = 0,      p(x, t, A) = 0, p(x, T, a) = 0,

(10) (x, t, a) ∈ Q × (0, A) (x, t, a) ∈ Σ × (0, A) (x, t) ∈ Q (x, a) ∈ Ω × (0, A),

then  ∗ ∗ 0, if y u (x, t, a)p(x, t, a) + y u (x, t, a) − c < 0 ∗ ∗ ∗ u (x, t, a) = L, if y u (x, t, a)p(x, t, a) + y u (x, t, a) − c > 0.

(11)

Proof (Sketch). We say that p ∈ L2 (Q × (0, A)) which belongs to C(S; L2 (Ω )) ∩ AC(S; L2 (Ω )) ∩ L2 (S; H 1 (Ω )) ∩ L2loc (S; H 2 (Ω )) for almost any characteristic line a − t = const., (t, a) ∈ (0, T ) × (0, A) and satisfies  A  ∗  Dp(x, t, a) + d∆p(x, t, a) = µ(a)p(x, t, a) + k p(x, t, a) y u (x, t, a)da  0   0   A   ∗  u  + k0 y (x, t, a)p(x, t, a)da − β(a)p(x, t, 0)    0 + χω (x)u∗ (x, t, a)(1 + p(x, t, a)), a.e. in Q × (0, A)    ∂ p(x, t, a) = 0, a.e. in Σ × (0, A)  ν     lim p(·, t − ε, A − ε) = 0, in L2 (Ω ), a.e. in (0, T )   ε→0+   lim p(·, T − ε, a − ε) = 0, in L2 (Ω ), a.e. in (0, A). ε→0+

Actually, by p(·, t, 0) we understand limε→0+ p(·, t + ε, ε) in L2 (Ω ). The existence and uniqueness of the solution to (10) follows via a fixed point result (for details see [2, p. 143]). Consider now an arbitrary v ∈ L∞ (ω × (0, T ) × (0, A)) such that u∗ + εv ∈ U˜ for sufficiently small ε > 0. The following lemma is a consequence of the evaluations for age-dependent population dynamics (see [2,5]) and of the ideas used for the proof of Lemma 3. Lemma 6. For any v as in Lemma 3, we have that ∗ ∗ ∗ 1 u∗ +εv y u +εv → y u , [y − y u ] → z in L∞ (Q × (0, A)), ε as ε → 0+, where z is the solution to  ∂t z(x, t, a) + ∂a z(x, t, a) − d∆z(x, t, a) = −µ(a)z(x, t, a)   A  A     u∗ u∗  − k z(x, t, a) y (x, t, a)da − k y (x, t, a) z(x, t, a)da 0 0    0 0  ∗   − χ (x)[u∗ (x, t, a)z(x, t, a) + v(x, t, a)y u (x, t, a)], (x, t, a) ∈ Q × (0, A) ω

∂ν z(x, t, a) = 0,     A     β(a)z(x, t, a)da,  z(x, t, 0) =  0   z(x, 0, a) = 0,

(x, t, a) ∈ Σ × (0, A) (x, t) ∈ Q (x, a) ∈ Ω × (0, A).

(12)

S. Anit¸a et al. / Nonlinear Analysis 147 (2016) 191–212

202

The existence and uniqueness of the solution to (12) follows via a fixed point result (for details see [2, p. 145]). Since for any v as above and for sufficiently small ε > 0 we have that Φ(u∗ ) ≥ Φ(u∗ + εv), then after an easy calculation and using Lemma 6 we get that  A T  ∗ [u∗ (x, t, a)z(x, t, a) + v(x, t, a)y u (x, t, a) − cv(x, t, a)]dx dt da ≤ 0. (13) ω

0

0

By multiplying the first equation in (10) by z and integrating over Q × (0, A) we obtain that T1 = T2 , where  A T  z(x, t, a)[∂t p(x, t, a) + ∂a p(x, t, a) + d∆p(x, t, a)]dx dt da, T1 = Ω

0

0

and 

A

T







T2 =

A

z(x, t, a)[µ(a)p(x, t, a) + k0 p(x, t, a) 0

Ω

0 A



∗

y u (x, t, a)da

0 ∗

y u (x, t, a)p(x, t, a)da − β(a)p(x, t, 0) + χω (x)u∗ (x, t, a)(1 + p(x, t, a))]dx dt da.

+ k0 0

By integrating by parts with respect to t and to a, we get that  T T1 = − z(x, t, 0)p(x, t, 0)dx dt Ω

0 A



T





−

p(x, t, a)[∂t z(x, t, a) + ∂a z(x, t, a) − d∆z(x, t, a)]dx dt da. 0

Ω

0

By (12) we obtain that A



T





T1 = −

β(a)z(x, t, a)p(x, t, 0)dx dt da 0 A



Ω

0 T







−

p(x, t, a)[−µ(a)z(x, t, a) − k0 z(x, t, a) 0

Ω

0

A

∗

y u (x, t, a)da

0

∗

− k0 y u (x, t, a)

A



z(x, t, a)da − χω (x)(u∗ (x, t, a)z(x, t, a)

0 ∗

+ v(x, t, a)y u (x, t, a))]dx dt da. AT  AT  ∗ Since T1 = T2 , we may infer that 0 0 ω u∗ zdx dt da = 0 0 ω vy u pdx dt da, and by (13) we get that  A T  ∗ ∗ v(x, t, a)[y u (x, t, a)p(x, t, a) + y u (x, t, a) − c]dx dt da ≤ 0, 0

0

ω

for any v such that u + εv ∈ U˜ for sufficiently small ε > 0. This implies (11). ∗

∗

Remark 3. If we multiply the first equation in (10) by y u and integrate over Ω × (0, T ) × (0, A), we obtain after some calculation that  A Φ(u∗ ) = − y0 (x, a)p(x, 0, a)dx da 0



Ω A A



T



− k0 0

0

0

Ω

∗

∗

y u (x, t, a′ )p(x, t, a′ )y u (x, t, a)dx dt da′ da.

S. Anit¸a et al. / Nonlinear Analysis 147 (2016) 191–212

203

In the special case when c = 0 we obtain that  0, if 1 + p(x, t, a) < 0 ∗ u (x, t, a) = L, if 1 + p(x, t, a) > 0, and that p is the solution to ∂ p(x, t, a) + ∂ p(x, t, a) + d∆p(x, t, a) = µ(a)p(x, t, a) t a   A      + k0 p(x, t, a) y(x, t, a)da    0   A    + k0 y(x, t, a)p(x, t, a)da − β(a)p(x, t, 0) 0   + χω (x)L(1 + p(x, t, a))+ , (x, t, a) ∈ Q × (0, A)      ∂ν p(x, t, a) = 0, (x, t, a) ∈ Σ × (0, A)     p(x, t, A) = 0, (x, t) ∈ Q   p(x, T, a) = 0, (x, a) ∈ Ω × (0, A),

(14)

where y is the solution of (9) corresponding to u := u∗ . If we have in addition that p(x, t, a) ̸= −1 a.e. in ω × (0, T ) × (0, A), then u∗ (x, t, a) = L · H(1 + p(x, t, a)). Sufficient conditions on β and µ may be found such that this condition be satisfied. We may conclude that the maximal harvest (for a given ω) is given by  A Φ(u∗ ) = − y0 (x, a)p(x, 0, a)dx da 0



Ω A A



T



− k0 0

0

y(x, t, a′ )p(x, t, a′ )y(x, t, a)dx dt da′ da,

Ω

0

where (y, p) is the solution to (9) and (14), and u := L · H(1 + p). 4. Optimal harvesting for a spatially structured population with diffusion and local logistic term This section is devoted to the following optimal harvesting problem related to a spatially structured population with diffusion and local logistic term:  T Maximize [u(x, t)y u (x, t) − cu(x, t)]dx dt, (OH3) u∈U

0

ω

∞

where U = {v ∈ L (ω × (0, T )); 0 ≤ v(x, t) ≤ L a.e.} (L ∈ (0, +∞), c ∈ [0, +∞) are constants) and y u is the solution to the following evolution equation  2  ∂t y − d∆y = η(x)y − ζ(x)y − χω (x)uy, (x, t) ∈ Q (15) ∂ν y(x, t) = 0, (x, t) ∈ Σ  y(x, 0) = y (x), x ∈ Ω . 0 The second order logistic term here, ζ(x)y(x, t)2 is local, thus indicating that the competition for resources is local, depending only on the population density at the same position x. Here are the hypotheses that are going to be used in this section: (H1)′′ η, y0 ∈ L∞ (Ω ), y0 (x) ≥ 0 a.e. in Ω and ∥y0 ∥L∞ (Ω) > 0; (H2)′′ ζ ∈ L∞ (Ω ), ζ(x) ≥ 0 a.e. in Ω . Using the ideas in the previous sections it is possible to prove the following results: Theorem 7. There exists at least one optimal control u∗ for (OH3).

S. Anit¸a et al. / Nonlinear Analysis 147 (2016) 191–212

204

Theorem 8. If u∗ is an optimal control for (OH3) and if p is the solution  u∗ ∗  ∂t p + d∆p = −η(x)p + 2ζ(x)y p + χω (x)u (1 + p), ∂ν p(x, t) = 0,  p(x, T ) = 0,

to (x, t) ∈ Q (x, t) ∈ Σ x ∈ Ω,

(16)

then  ∗ ∗ 0, if y u (x, t)p(x, t) + y u (x, t) − c < 0 ∗ ∗ u (x, t) = L, if y u (x, t)p(x, t) + y u (x, t) − c > 0. ∗

(17)

Remark 4. Problem (15) is a limit version of (1), for k(x, x′ ) = ζ(x)δx (x′ ). Theorems 7 and 8 are not obtained as particular cases from Theorems 1 and 2. Their proofs follow however in a similar manner. ∗

An important particular case is when c = 0. If we multiply the first equation in (16) by y u and integrate over Q, after an easy calculation (and using (15)) we obtain   T ∗ Φ(u∗ ) = − y0 (x)p(x, 0)dx − ζ(x)y u (x, t)2 p(x, t)dx dt, Ω

where Φ(u) =

T  0

Ω

Ω

0

u(x, t)y u (x, t)dx dt, and  ∗

u (x, t) =

0, if 1 + p(x, t) < 0 L, if 1 + p(x, t) > 0

∗

(by (17) and because y u (x, t) > 0 a.e. in Q). This however implies that p  u∗ +  ∂t p + d∆p = −η(x)p + 2ζ(x)y p + χω (x)L(1 + p) , ∂ν p(x, t) = 0,  p(x, T ) = 0,

(18) is the unique solution to (x, t) ∈ Q (x, t) ∈ Σ x ∈ Ω.

Remark 5. If under these last hypotheses there exists a subset A ⊂ ω × (0, T ) of positive measure, such ∗ that p(x, t) = −1 a.e. in A, then by (16) we get that 0 = η(x) − 2ζ(x)y u (x, t) a.e. inA, which implies that ∗ η(x) a.e. in A, and consequently (by (15)): y u (x, t) = 2ζ(x)     η(x) η(x) η(x) ∗ −d∆ = − u (x, t) 2ζ(x) 2ζ(x) 2 a.e. in A. If η(x) ̸= 0 a.e. in ω × (0, T ), then 2ζ(x) u (x, t) = d∆ η(x) ∗



η(x) ζ(x)

 +

η(x) 2

a.e. in A. Hence, if η(x) ̸= 0 a.e. in ω × (0, T ), and if   2ζ(x) η(x) η(x) + d∆ >L η(x) η(x) 2 a.e. in ω, or if 2ζ(x) d∆ η(x)



η(x) η(x)

 +

η(x) <0 2

S. Anit¸a et al. / Nonlinear Analysis 147 (2016) 191–212

205

a.e. in ω, then we get that p(x, t) ̸= −1 a.e. in ω × (0, T ), and so u∗ is a bang–bang control, and u∗ (x, t) = L · H(1 + p(x, t)) a.e. in ω × (0, T ). If

η ζ

is a constant function and if η < 0 or η > 2L a.e. in ω, then we have the same conclusion:   T Φ(u∗ ) = − y0 (x)p(x, 0)dx − ζ(x)y(x, t)2 p(x, t)dx dt, Ω

Ω

0

where (y, p) is the solution to   ∂t y − d∆y = η(x)y − ζ(x)y 2 − χω (x)LH(1 + p)y, (x, t) ∈ Q   ∂ p + d∆p = −η(x)p + 2ζ(x)yp + χ (x)L(1 + p)+ , (x, t) ∈ Q t ω  ∂ y(x, t) = ∂ p(x, t) = 0, (x, t) ∈ Σ ν ν   y(x, 0) = y (x), p(x, T ) = 0, x ∈ Ω. 0 Remark 6. In the particular case when c = 0 and ζ ≡ 0 (a particular case also for the problem in Section 2: for k ≡ 0), we may conclude that  Φ(u∗ ) = − y0 (x)p(x, 0)dx, Ω

where p is the unique solution to  +  ∂t p(x, t) + d∆p(x, t) = −η(x)p(x, t) + χω (x)L(1 + p(x, t)) , (x, t) ∈ Q ∂ν p(x, t) = 0, (x, t) ∈ Σ  p(x, T ) = 0, x ∈ Ω. If in addition η(x) ̸= 0 a.e. in ω, then there is a unique optimal control u∗ , which is given by (18). 5. Localizing an optimal subregion ω where the control acts We investigate here the optimal harvesting for a particular model of population dynamics and corresponding to the case when k ≡ 0 (there is no second order logistic term) just to avoid technical difficulties which might obscure the conceptual algorithm presented here. We wish to emphasize the optimal choice of the subregion ω where the control is localized. Consider a 2D domain Ω ⊂ R2 , and let ϕ be the implicit function associated with ω, the subregion of Ω in which we wish to apply our control, (ω = {x ∈ Ω ; ϕ(x) > 0}). Our optimal control problem here includes both the magnitude of the harvesting effort u ∈ Uω = {v ∈ L∞ (ω × (0, T )); 0 ≤ v(x, t) ≤ L a.e.}, and the choice of the subregion ω: Maximize Maximize Ψϕ (u), ϕ

u∈Uω

where 

T





u

u(x, t)y (x, t)dx dt − c1

Ψϕ (u) = 0

ω

 δ(ϕ(x))|∇ϕ(x)|dx − c2

Ω

H(ϕ(x))dx, Ω

and y u is the solution to   ∂t y(x, t) − d∆y(x, t) = η(x)y(x, t) − χω (x)u(x, t)y(x, t), (x, t) ∈ Q ∂ν y(x, t) = 0, (x, t) ∈ Σ  y(x, 0) = y (x), x ∈ Ω. 0 Assume that the hypotheses in Section 2 are satisfied. Using the remarks in Sections 2 and 4 we may rewrite the problem as     Minimize y0 (x)p(x, 0)dx + c1 δ(ϕ(x))|∇ϕ(x)|dx + c2 H(ϕ(x))dx , ϕ

Ω

Ω

Ω

S. Anit¸a et al. / Nonlinear Analysis 147 (2016) 191–212

206

where p is the solution to   ∂t p + d∆p = −η(x)p + LH(ϕ(x))(1 + p)H(1 + p), (x, t) ∈ Q ∂ν p(x, t) = 0, (x, t) ∈ Σ  p(x, T ) = 0, x ∈ Ω. We shall use the level set approach to treat this last problem (see also [16,23]). We will however approximate this problem by the following one, where the Heaviside function H is substituted by its mollified    version Hε (t) = 21 1 + π2 arctan εt , and its derivative by the mollified function δε (t) = π(ε2ε+t2 ) . So, for a small but fixed ε > 0, we shall deal in what follows with the following (regional harvesting) problem: Minimize J(ϕ),

(ROH)

ϕ

where ϕ : Ω −→ R is a smooth function and    ϕ J(ϕ) = y0 (x)p (x, 0)dx + c1 δε (ϕ(x))|∇ϕ(x)|dx + c2 Hε (ϕ(x))dx Ω

Ω

Ω

and pϕ is the solution to   ∂t p + d∆p = −η(x)p + LHε (ϕ(x))(1 + p)Hε (1 + p), (x, t) ∈ Q ∂ν p(x, t) = 0, (x, t) ∈ Σ  p(x, T ) = 0, x ∈ Ω.

(19)

The following result gives the directional derivative of J: Theorem 9. For any smooth functions ϕ, ψ : Ω −→ R we have that     ∇ϕ(x) dJ(ϕ)(ψ) = δε (ϕ(x)) − c1 div + c2 |∇ϕ(x)| Ω   T  ϕ ϕ ϕ −L (1 + p (x, t))Hε (1 + p (x, t))r (x, t)dt ψ(x)dx + c1 ∂Ω

0

δε (ϕ(x)) ∂ν ϕ(x)ψ(x)dσ, |∇ϕ(x)|

where rϕ is the solution to  ϕ ϕ ϕ  ∂t r − d∆r = η(x)r − LHε (ϕ(x))Hε (1 + p )r − LHε (ϕ(x))(1 + p )δε (1 + p )r, ∂ν r(x, t) = 0,  r(x, 0) = y (x), 0

(x, t) ∈ Q (x, t) ∈ Σ x ∈ Ω.

(20)

Proof. As in the proof of Lemma 3 it is possible to prove that 1 ϕ+θψ [p − pϕ ] → q ψ θ

in C([0, T ]; L∞ (Ω )),

as θ → 0, where q ψ is the solution to the problem   ∂t q + d∆q = −η(x)q + Lδε (ϕ(x))(1 + pϕ )Hε (1 + pϕ )ψ(x)    + LHε (ϕ(x))Hε (1 + pϕ )q + LHε (ϕ(x))(1 + pϕ )δε (1 + pϕ )q, (x, t) ∈ Q  ∂ν q(x, t) = 0, (x, t) ∈ Σ   q(x, T ) = 0, x ∈ Ω.

(21)

S. Anit¸a et al. / Nonlinear Analysis 147 (2016) 191–212

207

For any ϕ we have that 1 [J(ϕ + θψ) − J(ϕ)] = θ→0 θ



y0 (x)q ψ (x, 0)dx + c1

lim

Ω



δε′ (ϕ(x))ψ(x)|∇ϕ(x)|dx

Ω

 ∇ϕ(x) · ∇ψ(x) dx + c2 δε (ϕ(x))ψ(x)dx |∇ϕ(x)| Ω Ω   ∇ϕ(x) ψ ψ(x))dx = y0 (x)q (x, 0)dx + c1 div(δε (ϕ(x)) |∇ϕ(x)| Ω Ω     ∇ϕ(x) dx + c2 − c1 δε (ϕ(x))ψ(x)dx. δε (ϕ(x))ψ(x)div |∇ϕ(x)| Ω Ω 

+ c1

δε (ϕ(x))

After applying the Gauss–Ostrogradski’s formula (see [10, p. 13]) we obtain     ∇ϕ(x) ψ ψ(x)dx dJ(ϕ)(ψ) = y0 (x)q (x, 0)dx − c1 δε (ϕ(x))div |∇ϕ(x)| Ω Ω   δε (ϕ(x)) + c2 δε (ϕ(x))ψ(x)dx + c1 ∂ν ϕ(x)ψ(x)dσ. Ω ∂Ω |∇ϕ(x)|

(22)

By multiplying the first equation in (21) by r and integrating over Q we get after an easy calculation, and using (20) as well, that   T y0 (x)q ψ (x, 0)dx = − Lrϕ (x, t)δε (ϕ(x))(1 + pϕ (x, t))Hε (1 + pϕ (x, t))ψ(x)dx dt. (23) Ω

0

Ω

Now, from (22) and (23), we get that the conclusion of Theorem 9. Remark 7. One may like to notice that from Theorem 9 we get that the gradient descent with respect to ϕ (see [23]) is      ∂θ ϕ(x, θ) = δε (ϕ(x, θ)) c1 div ∇ϕ(x, θ) − c2    |∇ϕ(x, θ)|      T (24) +L (1 + pϕ (x, t))Hε (1 + pϕ (x, t))rϕ (x, t)dt , x ∈ Ω , θ > 0   0       δε (ϕ(x, θ)) ∂ν ϕ(x, θ) = 0, x ∈ ∂Ω , θ > 0. |∇ϕ(x, θ)| (θ is an artificial time). Remark 8. If there is an optimal ϕ∗ , then we get by Theorem 9 that ϕ∗ is a steady-state of (24). Numerical implementation From Theorem 9 we derive the following conceptual iterative algorithm, a semi-implicit gradient descent one, to improve at each step the region where the harvesting effort acts in order to obtain a smaller value for J. STEP 0: set k := 0, J (0) := 106 and θ0 > 0 a small constant initialize ϕ(0) = ϕ(0) (x, 0) STEP 1: compute p(k+1) the solution of (19) corresponding to ϕ(k) (·, 0)  compute J (k+1) = y0 (x)p(k+1) (x, 0)dx Ω   +c1 Ω δε (ϕ(k) (x, 0))|∇ϕ(k) (x, 0)|dx + c2 Hε (ϕ(k) (x, 0))dx. Ω

208

S. Anit¸a et al. / Nonlinear Analysis 147 (2016) 191–212

  Step 2: if J (k+1) − J (k)  < ε1 or J (k+1) ≥ J (k) then STOP else go to Step 3. Step 3: compute r(k+1) the solution of problem (20) corresponding to ϕ(k) (·, 0) and p(k+1) . Step 4: compute ϕ(k+1) using (24) and the initial condition ϕ(k+1) (x, 0) = ϕ(k) (x, θ0 ) and a semi-implicit timestep scheme Step 5: if ∥ϕ(k+1) − ϕ(k) ∥L2 < ε2 then STOP else k := k + 1 go to Step 1 ε1 > 0 in Step 2 and ε2 > 0 in Step 5 are prescribed convergence parameters. For more information about gradient (descent) methods, see [6, §2.3]. Numerical examples In order to simplify the discretization formulae for the numerical tests we have considered Ω = (0, 1) × (0, 1). We introduce equidistant discretization nodes for both axis corresponding to Ω . Thus, the domain Ω is approximated by a grid of (N + 1) × (N + 1) equidistant nodes, namely {(xi1 , xj2 ) : xi1 = (i − 1)∆x1 , xj2 = (j − 1)∆x2 , i, j = 1, N + 1, ∆x1 = ∆x2 = 1/N } ((x1 , x2 ) ∈ Ω is the space variable). The interval [0, T ] is also discretized by (M + 1) equidistant nodes, tk = (k −1)∆t, k = 1, 2, . . . , M +1, ∆t = T /M . The parabolic system from Step 1 is approximated by a finite difference method, descending with respect to time levels. An implicit scheme is used. The resulting algebraic linear system is solved by Gaussian elimination. The parabolic system from Step 3 is approximated also using a finite difference method, but ascending with respect to time levels. Integrals from Step 1 are numerically computed using Simpson’s method corresponding to the discrete grid. In Step 4, time is discretized using the Gauss–Seidel iterative method. With this method we make use of updated values of ϕ as soon as they become available. In the following, in every figure, the white area represents the subregion ω that provides the smallest value for J. x2 +x2 1 2

1 − 2 We consider a normal initial population density y0 (x1 , x2 ) = 2π , (x1 , x2 ) ∈ Ω . Let the diffusion e coefficient be d = 1, the final time T = 1, L = 1, and the regularization parameter ε = 1. We take the space discretization step ∆x1 = ∆x2 = 0.05, and the time discretization step ∆t = 0.05 since the finite difference method used is implicit. For the convergence tests we consider ε1 = ε2 = 0.001.

Test 1. Let the natural growth rate of the population be constant, for example η(x1 , x2 ) = 3, (x1 , x2 ) ∈ Ω . We consider ϕ(0) (x1 , x2 , 0) = (−x1 + 0.25)(−x2 + 0.5) + 41 sin(x1 − 0.4) sin(0.25 − x2 ), (x1 , x2 ) ∈ Ω . We penalize the length of ∂ω by c1 = 0.6 and c1 = 0.7 and the area of ω by c2 = 1. The corresponding results are shown in Fig. 1. The final iterations for c1 = 0.6 and c1 = 0.7 are quite close. Test 2. We consider the natural growth rate of the population η(x1 , x2 ) = ex1 +x2 , (x1 , x2 ) ∈ Ω . We penalize the length of ∂ω and the area of ω by varying c1 in the set {0.1, 0.2, 0.3, 0.4, 0.5}, respectively by c2 = 0.5. The initialization of ϕ is made by ϕ(0) (x1 , x2 , 0) = sin(5πx1 ) sin(5πx2 ), (x1 , x2 ) ∈ Ω , a function that produce a initial checkerboard shape. The results are shown in Fig. 2. We notice that increasing c1 , the area of the final iteration increases and the length of the boundary of the final iteration decreases. Test 3. The results from Fig. 3 are obtained using the same input data from example 2, except now c1 = 1, and varying c2 ∈ {0.3, 0.45, 1}. We start with the same initial subregion ω as in Test 2 (see Fig. 2(a)). We notice that increasing c2 , the area of the final iteration decreases.

S. Anit¸a et al. / Nonlinear Analysis 147 (2016) 191–212

(a) Initial ω.

(b) c1 = 0.6.

209

(c) c1 = 0.7.

Fig. 1. The representation of the initial iteration (a) and the final iterations of ω for c1 ∈ {0.6, 0.7} and c2 = 1.

Let us point out that if we try to increase the value of c2 , the level set function returning by the computer program is everywhere negative. The interpretation is: the penalty applied to the area is too big. In all these examples, the convergence was obtained by the first test in Step 2. 6. Final comments The methods we have developed in the previous sections may be used even for general optimal control problems of the form  T Maximize G(u(t), y u (t))dt, (OC) u∈U

0

where U is a convex and closed subset of L2 (0, T ; X) with X a real Hilbert space (a functional space) and y u is the unique (Carath´eodory) solution to  y ′ (t) = f (u(t), y(t)), t ∈ (0, T ) (25) y(0) = y0 . We assume that (25) is the abstract form of a reaction–diffusion system where the control acts in a subset ω of the whole domain Ω (and so U depends on ω), and that for any u ∈ U, (25) has a unique Carath´eodory solution. Here f : X × Y → Y , G : X × Y → R are continuously differentiable functions (Y is another real Hilbert space); we will denote by fu , fy , Gu , Gy the Gˆateaux derivatives with respect to u and y. Under appropriate assumptions (OC) admits at least one optimal control u∗ . Assume that there exists a unique (Carath´eodory) solution p to  ∗ ∗ p′ (t) = −fy (u∗ (t), y u (t))∗ p(t) + Gy (u∗ (t), y u (t)), t ∈ (0, T ) (26) p(T ) = 0. Here B ∗ denotes the adjoint of a linear and bounded operator B. A formal calculation allows us to conclude that ∗

∗

−Gu (u∗ (t), y u (t)) + fu (u∗ (t), y u (t))∗ p(t) ∈ NU (u∗ (t)) (see [4, p.59–62]). Here NU (w) ⊂ L2 (0, T ; X) is the normal cone at U in w. ∗ ∗ If we may conclude in addition that u∗ (t) = F(y u (t), p(t)), then (y u , p) is the solution to  ′  t ∈ (0, T ) y (t) = f (F(y(t), p(t)), y(t)), ′ p (t) = −fy (F(y(t), p(t)), y(t))∗ p(t) + Gy (F(y(t), p(t)), y(t)), t ∈ (0, T )  y(0) = y , p(T ) = 0. 0

S. Anit¸a et al. / Nonlinear Analysis 147 (2016) 191–212

210

(a) Initial ω.

(b) c1 = 0.1.

(c) c1 = 0.2.

(d) c1 = 0.3.

(e) c1 = 0.4.

(f) c1 = 0.5.

Fig. 2. The representation of the initial iteration (a) and the final iterations of ω for c1 ∈ {0.1, 0.2, 0.3, 0.4, 0.5} and c2 = 0.5.

(a) c2 = 0.3.

(b) c2 = 0.45.

(c) c2 = 1.

Fig. 3. The representation of the final iterations for ω for c1 = 0.1 and c2 ∈ {0.3, 0.45, 1}.

If in addition we have that f (u, y) = A(u)y, (f is linear with respect to y), then fy (u, y) = A(u). Multiplying ∗ (26) by y u we may conclude after an easy calculation (involving (25) as well) that the maximal value of the cost function is Φ(u∗ ) = −y0 · p(0) −



T

y(t) · Gy (F(y(t), p(t)), y(t))dt, 0

S. Anit¸a et al. / Nonlinear Analysis 147 (2016) 191–212

211

where (y, p) is the solution to  ′  t ∈ (0, T ) y (t) = A(F(y(t), p(t)))y(t), ′ p (t) = −A(F(y(t), p(t)))∗ p(t) + Gy (F(y(t), p(t)), y(t)), t ∈ (0, T )  y(0) = y , p(T ) = 0. 0 Appropriate maximization problems with respect to ω (when U depends on ω) may be formulated and treated in the same manner as in Section 5. Remark 9. An extremely important problem from practical point of view in fisheries corresponds to a constant harvesting effort L localized in ω = {x ∈ Ω ; ϕ(x) > 0}. The problem to be investigated does not depend anymore on a dual state p and may be formulated as       T ϕ L H(ϕ)y dx dt − c1 δ(ϕ)|∇ϕ|dx − c2 H(ϕ)dx , Maximize ϕ

0

Ω

Ω

Ω

where y ϕ is the solution to   ∂t y(x, t) − d∆y(x, t) = η(x)y(x, t) − LH(ϕ(x))y(x, t), (x, t) ∈ Q ∂ν y(x, t) = 0, (x, t) ∈ Σ  y(x, 0) = y (x), x ∈ Ω. 0 The method presented in Section 5 may be used to treat such a problem as well. Acknowledgments The work of Sebastian Anit¸a and Ana-Maria Mo¸sneagu 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”. The authors are grateful to the Anonymous Referee for addressing additional important references on a related subject; G.M. Coclite, and his coauthors are worth of thanks for their permission to access their unpublished preprints. References [1] B. Ainseba, S. Anit¸a, M. Langlais, Optimal control for a nonlinear age-structured population dynamics model, Electron. J. Differential Equations 28 (2003) 1–9. [2] S. Anit¸a, Analysis and Control of Age-Dependent Population Dynamics, Kluwer Acad. Publ, Dordrecht, 2000. [3] L.-I. Anit¸a, S. Anit¸a, V. Arn˘ autu, Optimal harvesting for periodic age-dependent population dynamics with logistic term, Appl. Math. Comput. 215 (2009) 2701–2715. [4] S. Anit¸a, V. Arn˘ autu, V. Capasso, An Introduction to Optimal Control Problems in Life Sciences and Economics. From Mathematical Models to Numerical Simulation with MATLAB, Birkh¨ auser, Basel, 2011. [5] S. Anit¸a, V. Capasso, A stabilization strategy for a reaction–diffusion system modelling a class of spatially structured epidemic systems (think globally, act locally), Nonlinear Anal. RWA 10 (2009) 2026–2035. [6] V. Arn˘ autu, P. Neittaanm¨ aki, Optimal Control from Theory to Computer Programs, Kluwer Acad. Publ, Dordrecht, 2003. [7] C. Arns, M. Knackstedt, K. Mecke, Characterising the morphology of disordered materials, in: K. Mecke, D. Stoyan (Eds.), Morphology of Condensed Matter, Springer Berlin, Heidelberg, 2002, pp. 37–74. [8] P. Baras, Compacit´ e de l’op´ erateur f → u solution d’une ´ equation non lin´ eaire (du/dt) + Au ∋ f , C. R. Acad. Sci. Paris, S´ er. A 286 (1978) 1113–1116. [9] V. Barbu, Mathematical Methods in Optimization of Differential Systems, Kluwer Acad. Publ, Dordrecht, 1994. [10] V. Barbu, Partial Differention Equations and Boundary Value Problems, Kluwer Acad. Publ, Dordrecht, 1998. [11] A.O. Belyakov, V.M. Veliov, On optimal harvesting in age-structured populations, Research Report 2015-08, ORCOS, TU Wien, 2015. [12] A. Bressan, G.M. Coclite, W. Shen, A Multidimensional optimal-harvesting problem with measure-valued solutions, SIAM J. Control Optim. 51 (2013) 1186–1202.

212

S. Anit¸a et al. / Nonlinear Analysis 147 (2016) 191–212

[13] M. Brokate, Pontyagin’s principle for control problems in age-dependent population dynamics, J. Math. Biol. 23 (1985) 75–101. [14] M. Brokate, On a certain optimal harvesting problem with continuous age structure, in: K.-H. Hoffmann, W. Krabs (Eds.), Optimal Control of Partial Differential Equations II, Birkh¨ auser, Boston, 1987, pp. 29–42. [15] D. Bucur, G. Buttazzo, Variational Methods Shape Optimization Problems, Birkh¨ auser, Basel, 2005. [16] T.F. Chan, L.A. Vese, Active contours without edges, IEEE Trans. Image Process. 10 (2001) 266–277. [17] G.M. Coclite, M. Garavello, A time dependent optimal harvesting problem with measure valued solutions. Preprint, 2016. [18] G.M. Coclite, M. Garavello, L.V. Spinolo, Optimal strategies for a time-dependent harvesting problem. Preprint, 2016. [19] M.C. Delfour, J.-P. Zolesio, Shapes and Geometries. Metrics, Analysis, Differential Calculus and Optimization, second ed., SIAM, Philadelphia, 2011. [20] K.R. Fister, S. Lenhart, Optimal harvesting in an age-structured predator–prey model, Appl. Math. Optim. 54 (2006) 1–15. [21] A. Friedman, Partial Differential Equations of Parabolic Type, in: Dover Books on Mathematics Series, Dover Publ, New York, 2008. [22] S. Genieys, V. Volpert, P. Auger, Pattern and waves for a model in population dynamics with nonlocal consumption of resources, Math. Model. Nat. Phenom. 1 (2006) 65–82. [23] P. Getreuer, Chan-Vese Segmentation, IPOL J. Image Process. Online 2 (2012) 214–224. [24] M.E. Gurtin, L.F. Murphy, On the optimal harvesting of age-structured populations: some simple models, Math. Biosci. 55 (1981) 115–136. [25] M.E. Gurtin, L.F. Murphy, On the optimal harvesting of persistent age-structured populations, J. Math. Biol. 13 (1981) 131–148. [26] H. Hadwiger, Vorlesungen u ¨ber Inhalt, Oberfl¨ ache, und Isoperimetrie, Springer-Verlag, Berlin, 1957. [27] Z.R. He, Optimal harvesting of two competing species with age dependence, Nonlinear Anal. RWA 7 (2006) 769–788. [28] A. Henrot, M. Pierre, Variation et Optimisation de Formes, Math´ ematiques et Applications, Springer-Verlag, Berlin, 2005. [29] N. Hritonenko, Y. Yatsenko, Optimization of harvesting age in integral age-dependent model of population dynamics, Math. Biosci. 195 (2005) 154–167. [30] M. Iannelli, Mathematical Theory of Age-Structured Population Dynamics, in: Applied Mathematics Monographs - C.N.R., Giardini Editori e Stampatori, Pisa, 1995. [31] D.A. Klain, G.C. Rota, Introduction to Geometric Probability, Cambridge University Press, Cambridge, U.K, 1997. [32] S. Lenhart, Using optimal control of parabolic PDEs to investigate population questions, NIMBioS, April 9–11, 2014; https://www.fields.utoronto.ca/programs/scientific/13-14/BIOMAT/presentations/lenhartToronto3.pdf. [33] S. Lenhart, J.T. Workman, Optimal Control Applied to Biological Models, Chapman and Hall, 2007. [34] Z. Luo, Optimal harvesting problem for an age-dependent n-dimensional food chain diffusion model, Appl. Math. Comput. 186 (2007) 1742–1752. [35] Z. Luo, W.T. Li, M. Wang, Optimal harvesting control problem for linear periodic age-dependent population dynamics, Appl. Math. Comput. 151 (2004) 789–800. [36] L.F. Murphy, S.J. Smith, Optimal harvesting of an age-structured population, J. Math. Biol. 29 (1990) 77–90. [37] S. Osher, R. Fedkiw, Level Set Methods and Dynamic Implicit Surfaces, Springer, New York, 2003. [38] J. Ohser, F. M¨ ucklich, Statistical Analysis of Microstructures in Materials Science. Statistics in Practice, John Wiley & Sons, New York, 2000. [39] A. Pazy, Semigroups of Linear Operators and Applications to Partial Differential Equations, Springer-Verlag, New York, 1983. [40] M.H. Protter, H.F. Weinberger, Maximum Principles in Differential Equations, Springer-Verlag, New York, 1984. [41] J.A. Sethian, Level Set Methods and Fast Marching Methods. Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision, and Materials Science, Cambridge University Press, 1999. [42] J. Sokolowski, J.-P. Zolesio, Introduction to Shape Optimization, Springer-Verlag, Berlin, 1992. [43] G. Webb, Theory of Nonlinear Age-Dependent Population Dynamics, Marcel Dekker, New York, 1985. [44] C. Zhao, M. Wang, P. Zhao, Optimal harvesting problems for age-dependent interacting species with diffusion, Appl. Math. Comput. 163 (2005) 117–129. [45] C. Zhao, P. Zhao, M. Wang, Optimal harvesting for nonlinear age-dependent population dynamics, Math. Comput. Modelling 43 (2006) 310–319.