- Email: [email protected]
Optimal strategies to diminish a pest population via bilinear controls
Applied Mathematics Letters 40 (2015) 7–12
Contents lists available at ScienceDirect
Applied Mathematics Letters journal homepage: www.elsevier.com/locate/aml
Optimal strategies to diminish a pest population via bilinear controls Laura-Iulia Aniţa a , Sebastian Aniţa b,c,∗ , Costică Moroşanu b a
Faculty of Physics, ‘‘Al.I. Cuza’’ University of Iaşi, Iaşi 700506, Romania
b
Faculty of Mathematics, ‘‘Al.I. Cuza’’ University of Iaşi, Iaşi 700506, Romania
c
‘‘Octav Mayer’’ Institute of Mathematics, Iaşi 700506, Romania
article
info
Article history: Received 24 July 2014 Received in revised form 4 September 2014 Accepted 4 September 2014 Available online 16 September 2014 Keywords: Population dynamics Bilinear control problems Large-time behavior Age-dependent harvesting rate Gradient-type algorithm
abstract In this work we model and analyze two control strategies to diminish a pest population using traps. The action of any trap depends on its age. Both problems contain bilinear control rates. The large-time behavior of the model with time-periodic inflow is investigated. The first control strategy deals with a finite horizon problem while the second one is related to a time-periodic problem. We obtain Pontryagin’s principle for both control problems. A special attention is given to the periodic problem. Pontryagin’s principle is used to derive a conceptual gradient-type algorithm to approximate the optimal solution. Numerical tests are given. © 2014 Elsevier Ltd. All rights reserved.
1. Mathematical modeling Let h(x, t ) denote the density of a pest population (of insects) at location x ∈ Ω of the habitat (Ω ⊂ R2 ), and time t ≥ 0. Assume that a ratio c˜ of the whole pest population may diffuse and the ratio (1 − c˜ ) of it may migrate. If the diffusion rate in the absence of migration is d˜ and if the proportion of pests emigrating from x to other locations in the absence of diffusion is g˜ (x), then if we denote by d = c˜ d˜ and by g = (1 − c˜ )˜g , the dynamics of h is described by the following model:
h ∂ h = d ∆ h + r ( x ) h 1 − − g (x)h t K (x) + k(x, x′ )g (x′ )h(x′ , t )dx′ + f (x, t ) − u˜ (x, t )χω (x)h(x, t ), Ω ∂ν h = 0, h(x, jT +) = ρ(x)h(x, jT −), h(x, 0+) = h0 (x),
x ∈ Ω , t ∈ R+ \ T N
(1)
x ∈ ∂ Ω , t ∈ R+ \ T N x ∈ Ω , j ∈ N∗ .
Here, r (x) = growth rate of pests at location x; K (x) = carrying capacity of pests at location x; d = diffusion constant of pests when migration occurs; g (x) = proportion of pests emigrating from x to other locations; k(x, x′ ) = proportion of pests immigrating to x from location x′ ; T = period of one year; f (x, t ) = density of pest inflow; u˜ (x, t ) = harvesting rate (which acts only in the open subset ω ⊂ Ω ); χω = characteristic function ofω. The pest population is subject to an inflow f (x, t ) of population coming from outside Ω . We assume here that the Neumann (meaning an isolated habitat with respect
∗
Corresponding author at: Faculty of Mathematics, ‘‘Al.I. Cuza’’ University of Iaşi, Iaşi 700506, Romania. Tel.: +40 749 33 65 06. E-mail addresses: [email protected] (L.-I. Aniţa), [email protected], [email protected] (S. Aniţa), [email protected] (C. Moroşanu).
http://dx.doi.org/10.1016/j.aml.2014.09.003 0893-9659/© 2014 Elsevier Ltd. All rights reserved.
8
L.-I. Aniţa et al. / Applied Mathematics Letters 40 (2015) 7–12
to diffusion) homogeneous boundary condition holds for the pest population; ∂ν h(x, t ) denotes the normal derivative depending on location x. We also assume that h0 (x) is the initial density of pest population at location x. Recall that basic properties of the solutions to some biological reaction–diffusion systems with nonlocal terms have been recently investigated in Genieys et al. [1] and Aniţa and Capasso [2]. We assume that a certain ratio of the pest population dies during the winters. This occurs in a different proportion α(x), at different locations. In conclusion, at the beginning of the new year, the pest density is h(x, jT +) = h(x, jT −) − α(x)h(x, jT −) = ρ(x)h(x, jT −) (ρ(x) = 1 − α(x)). For an arbitrary function w, we have denoted by w(x, jT +) = limε↘0 w(x, jT + ε), w(x, jT −) = limε↘0 w(x, jT − ε), if these limits exist. We shall work under the following hypotheses: (H1) Ω ⊂ R2 is a bounded domain with a sufficiently smooth boundary ∂ Ω , d, T ∈ (0, +∞); (H2) r , K , g ∈ L∞ (Ω ), r (x) ≥ r0 , K (x) ≥ K0 , g (x) ≥ 0, a.e. x ∈ Ω , where r0 , K0 are positive constants; (H3) k ∈ L∞ (Ω × Ω ), k(x, x′ ) ≥ 0 a.e. in Ω × Ω , f ∈ L∞ (Ω × (0, +∞)), f (x, t ) ≥ 0 a.e. in Ω × (0, +∞), u˜ ∈ L∞ (ω × (0, +∞)); (H4) ρ, h0 ∈ L∞ (Ω ), c ≤ ρ(x) ≤ 1, a.e. x ∈ Ω (where c ∈ (0, 1)), 0 ≤ h0 (x), a.e. x ∈ Ω . The existence, uniqueness and nonnegativity of a solution to (1) are derived as in Capasso [3]. Our goal is to diminish the pest population by usage of traps and keeping this in balance with the cost of the intervention (control). We propose here two strategies and analyze them. Notice that unlike the use of pesticides in order to diminish the pest population, the use of traps is safe for the human population. Let us also notice that a trap attracts and captures a certain ratio of the pest population which depends on the age of the trap (the time elapsed from the moment when the trap has been situated at a certain position x ∈ ω). Actually this ratio decays exponentially. The harvesting rate at position x and moment t is obtained by cumulating the actions at moment t of all traps situated at position x. During winters, the traps loose their attracting and catching capabilities. So, if v(x, t ) is the density of traps at position x ∈ ω and moment t ≥ 0, and uv (x, t , a) is the solution to the following age-structured model:
∂t u + ∂a u + ηu = 0, u(x, t , 0) = v(x, t ), u(x, jT +, a) = 0,
x ∈ ω, t ∈ R+ \ T N, a ∈ (0, T ) x ∈ ω, t ∈ R+ x ∈ ω, j ∈ N, a ∈ (0, T )
(2)
(here η > 0 is the degradation rate of the traps), then the harvesting rate at position x and moment t is u˜ v (x, t ) =
T
uv (x, t , a)da.
(3)
0
In what follows we assume (H5) v ∈ L∞ (ω × (0, +∞)), 0 ≤ v(x, t ) ≤ M, a.e. in ω × (0, +∞), where M is a positive constant (representing the maximal density of traps). It is obvious that u˜ v ∈ L∞ (ω × (0, +∞)) and so there exists a unique and nonnegative solution hv to (1) corresponding to u˜ := u˜ v . 2. An optimal control problem with finite horizon Assume that our goal is to diminish (using traps) the negative effect of the pest population over a time period [0, L], where L = mT , m ∈ N∗ . We wish to have a diminished pest population at the end of our intervention, at a reasonable cost. Denote V = {v ∈ L∞ (ω × (0, L)); 0 ≤ v(x, t ) ≤ M a.e. in ω × (0, L)}. For a fixed v ∈ V we denote by hv the solution to (1) corresponding to u˜ := u˜ v given by (2) and (3) (where we consider v ≡ 0 in ω × (T , +∞)). Actually, the first problem to be investigated is the following one:
(P1 ) min v∈V
L 0
Ω
hv (x, t )dxdt + ζ1
Ω
hv (x, L−)dx + ζ2
L 0
ω
v(x, t )dxdt ,
where ζ1 , ζ2 are positive constants. We may take smaller or bigger values for ζ1 and ζ2 . A big ζ1 means that we are mainly interested to have a small pest population at the moment L. Taking small ζ2 means that we are interested to significantly reduce the pest population no matter the cost and a big ζ2 means that we are interested to reduce the pest population but at a small cost (with a small effort). It follows in a standard manner (see Aniţa et al. [4], Barbu [5], Lenhart and Workman [6]) that there exists at least an optimal control v ∗ for (P1 ). Our main goal in this section is to derive the first order necessary optimality conditions. Actually, the following result holds: Theorem 1 (Pontryagin’s principle). Let v ∗ be an optimal control of (P1 ) and q be a solution to the following problem:
∗ 2rhv ∗ ∂t q = −d∆q − rq + q + gq − g (x) k(x′ , x)q(x′ , t )dx′ + u˜ v χω q + 1,
x ∈ Ω , t ∈ [0, L] \ T N
∂ q = 0, ν q(x, jT −) = ρ(x)q(x, jT +), q(x, L−) = −ζ1 ,
x ∈ ∂ Ω , t ∈ [0, L] \ T N x ∈ Ω , j ∈ {1, 2, . . . , m − 1}.
K
Ω
(4)
L.-I. Aniţa et al. / Applied Mathematics Letters 40 (2015) 7–12
9
Then
v (x, t ) = ∗
0, M,
if θ (x, t , 0) < ζ2 if θ (x, t , 0) > ζ2 ,
(5)
where θ is the solution to
∗ ∂t θ + ∂a θ − ηθ = q(x, t )hv (x, t ), θ (x, t , T ) = 0, θ (x, jT −, a) = 0,
x ∈ ω, t ∈ [0, L] \ T N, a ∈ (0, T ) x ∈ ω, t ∈ (0, L) x ∈ ω, j ∈ {1, 2, . . . , m}, a ∈ (0, T ).
(6)
Problem (4) is called the dual problem. The method to construct it can be found in Aniţa et al. [4]. Existence and uniqueness of a solution to (4) follows as in Aniţa et al. [4]. Proof. Let w ∈ L∞ (ω×(0, L)) be such that for any ε > 0 sufficiently small we have v ∗ +εw ∈ V (i.e. 0 ≤ v ∗ (x, t )+εw(x, t ) ≤ M a.e. in ω × (0, L)). Since v ∗ is an optimal control for (P1 ) we deduce that ∗ ∗ hv +εw − hv
L 0≤
ε
Ω
0
dx dt + ζ1
∗ ∗ hv +εw (x, L−) − hv (x, L−)
ε
Ω
dx + ζ2
L 0
ω
w(x, t )dx dt ,
(7)
for any ε > 0 sufficiently small. The following auxiliary result will be used in what follows: Lemma. The following convergences yield ∗ ∗ hv +εw → hv ,
∗ ∗ hv +εw − hv
→ z , in C ([nT +, (n + 1)T −]; L∞ (Ω )), ε for n ∈ {0, 1, . . . , m − 1}, as ε → 0, where z is the solution to ∗ 2rhv z − gz ∂ z = d ∆ z + rz − t K ∗ ∗ + k(x, x′ )g (x′ )z (x′ , t )dx′ − u˜ w χω hv − u˜ v χω z , x ∈ Ω , t ∈ [0, L] \ T N Ω x ∈ ∂ Ω , t ∈ [0, L] \ T N ∂ν z = 0, z (x, jT +) = ρ(x)z (x, jT −), z (x, 0+) = 0, x ∈ Ω , j ∈ {1, 2, . . . , m − 1}.
(8)
The proof of the lemma follows as in Aniţa et al. [4], Barbu [5], Lenhart and Workman [6]. Passing to the limit in (7) and get
L 0≤ Ω
0
z (x, t )dx dt + ζ1
Ω
z (x, L−)dx + ζ2
L ω
0
w(x, t )dx dt .
(9)
If we multiply the equation in (4) by z and integrate over Ω × (0, L) we obtain m−1
Ω
q(x, L−)z (x, L−)dx +
j =1
[q(x, jT −)z (x, jT −)dx − q(x, jT +)z (x, jT +)]dx
q(x, 0+)z (x, 0+)dx −
− Ω
L = Ω
0
Ω
z −d∆q − rq +
L
2rh
0 v∗
K
Ω
q(x, t )∂t z (x, t )dx dt
q + gq − g (x)
Ω
v∗
k(x , x)q(x , t )dx + u˜ χω q dx dt + ′
′
′
L 0
Ω
z (x, t )dx dt .
Taking into account that q(x, L−) = −ζ1 a.e. x ∈ Ω and that z satisfies (8) we get after an easy calculation that ζ1 L−)dx +
L
Ω
0
L ω
0
z (x, t )dx dt = T
w
v ˜w ω u (x, t )h (x, t )q(x, t )dx dt , and consequently by (9) and (3) we get that ∗
L 0
v∗
u (x, t , a)h (x, t )q(x, t )da dx dt + ζ2 0
L 0
ω
L
Ω
z (x,
w(x, t )dx dt ≥ 0.
If we multiply (6) by uw and integrate over ω × (0, L) × (0, T ) we get (taking into account (2) and (6)) that v∗
L
L T 0
ω 0
uw
(x, t , a)h (x, t )q(x, t )da dx dt = − 0 ω w(x, t )θ (x, t , 0)dx dt , and consequently 0 ω w(x, t )[ζ2 − θ (x, t , 0)]dx dt ≥ 0, for any w ∈ L∞ (ω × (0, L)) such that for any ε > 0 sufficiently small we have v ∗ + εw ∈ V . We finally conclude (5). 3. A periodic optimal control problem If we consider the domain Ω to be small, then it is reasonable to assume that the rates r (x), K (x), g (x), ρ(x), k(x, x′ ) are independent with respect to x, ω ≡ Ω (we may act in the whole domain), and the inflow f satisfies f ∈ L∞ (0, +∞),
10
L.-I. Aniţa et al. / Applied Mathematics Letters 40 (2015) 7–12
f (t ) = f (t + T ), f (t ) > 0 a.e. in (0, +∞). The seasonality of the inflow is also a realistic condition. Therefore it is natural to deal with space-independent and T -periodic controls. Let
VT = {w ∈ L∞ (0, +∞); w(t ) = w(t + T ), 0 ≤ w(t ) ≤ M a.e. in (0, +∞)}. For any v ∈ VT , problem (1) has a unique solution hv which is nonnegative. Let us prove first that limt →+∞ ∥hv (·, t ) − h˜ v (t )∥L∞ (Ω ) = 0, where h˜ v is the unique T -periodic positive solution to
h − gh + k(x′ )gh(x′ , t )dx′ + f (t ) − u˜ v (t )h, ∂ h = d ∆ h + rh 1 − t K
x ∈ Ω , t ∈ (0, +∞) \ T N
∂ν h = 0, h(x, jT +) = ρ h(x, jT −),
x ∈ ∂ Ω , t ∈ (0, +∞) \ T N x ∈ Ω , j ∈ N∗ ,
Ω
(10)
which actually is space independent. Hence h˜ v is the unique positive T-periodic solution to
h′ = rh 1 − h − gh + gh k(x′ )dx′ + f (t ) − u˜ v (t )h, K Ω h(jT +) = ρ h(jT −),
t ∈ (0, +∞) \ T N
(11)
j∈N . ∗
Using the comparison result for parabolic equations with nonlocal term (see Aniţa and Capasso [2]) we obtain that the sequence n → hv1 (t +nT +) is increasing for any t ∈ [0, T ), where hv1 is the solution to (1) corresponding to h0 ≡ 0. Actually this solution is space independent. On the other hand for any sufficiently large initial space independent value h02 ≥ ∥h0 ∥L∞ (Ω ) , we get via the same comparison result that n → hv2 (t + nT +) is decreasing, where hv2 is the solution to (1) corresponding to h0 := h02 . The same comparison result implies that 0 < hv1 (t ) ≤ hv (x, t ) ≤ hv2 (t ),
a.e. in Ω × (0, +∞).
v
(12)
v
v
The monotonicity of n → h1 (t + nT +) and of n → h2 (t + nT +) for any t ∈ [0, T ), implies that h1 (·, t ) − h˜ v1 (t ) → 0, hv2 (·, t ) − h˜ v2 (t ) → 0, in L∞ (Ω ) as t → +∞, where h˜ v1 and h˜ v2 are T -periodic solutions to (11) (and (10)). We use now an h˜ v
argument from Smoller [7]: we multiply the equation in h˜ v1 by ˜1v and the equation in h˜ v2 by ˜ v1 2 and integrate their differh2 (h2 )
ence on [0, T ]. Taking into account that 0 < h˜ v1 (t ) ≤ h˜ v2 (t ), a.e. in (0, +∞), we conclude that necessarily h˜ v1 ≡ h˜ v2 . Passing to the limit in (12) we may infer that limt →+∞ ∥hv (·, t ) − h˜ v (t )∥L∞ (Ω ) = 0, where h˜ v1 ≡ h˜ v2 ≡ h˜ v is the unique positive T-periodic solution to (11). For any v ∈ VT we have that (n+1)T
Ω
nT
hv (x, t )dx dt + ζ2
(n+1)T
nT
Ω
v(t )dx dt → meas(Ω )
T
h˜ v (t )dt + ζ2
0
T
v(t )dt ,
0
as n → +∞. Here meas(Ω ) is the Lebesgue measure of Ω . Therefore, a natural approach is that of minimizing the pest population along a time interval of length T (on long term) at a small cost via a T -periodic harvesting effort v . The last convergence implies that the most natural way to mathematically formulate this problem is
(P2 ) min
T
v∈VT
h˜ v (t )dt + ζ2
0
T
v(t )dt . 0
Existence of an optimal control v ∗ for problem (P2 ) follows in a standard manner. Moreover Theorem 2. Let v ∗ be an optimal control of (P2 ) and q be the T -periodic solution to the following problem:
v∗ ∗ q′ = −rq + 2r h˜ q + gq − gq k(x′ )dx′ + u˜ v q + 1, K Ω q(jT −) = ρ q(jT +),
t ∈ (0, +∞) \ T N
(13)
j∈N . ∗
Then
v (t ) = ∗
0, M,
if θ (t , 0) < ζ2 if θ (t , 0) > ζ2 ,
(14)
where θ is the solution to
v∗ ∂t θ + ∂a θ − ηθ = q(t )h˜ (t ), θ (t , T ) = 0, θ (jT −, a) = 0,
t ∈ (0, +∞) \ T N, a ∈ (0, T ) t ∈ (0, +∞) j ∈ N∗ , a ∈ (0, T ).
(15)
It is possible to prove that problem (13) has a unique solution using the same technique as for proving that problem (11) has a unique solution. For the proof we have to follow the steps in the previous section.
L.-I. Aniţa et al. / Applied Mathematics Letters 40 (2015) 7–12
11
Remarks. For the more complex situation when r , K , g , k, f , ρ are space dependent, we may however use comparison T T results in order to find a control v ∗ which provides a low upper bound for the cost functional 0 Ω h˜ v (x, t )dx dt + ζ2 0 Ω v(x, t )dx dt. Such result is based on the last theorem. For other recent results concerning some other optimal control problems with bilinear terms (with harvesting rates) in biology and economics, but of a different type we refer to Aniţa et al. [8], Fister and Lenhart [9], Zhao et al. [10], Hritonenko and Yatsenko [11], Hritonenko et al. [12], Belyakov and Veliov [13], Behringer and Upmann [14] and He and Liu [15]. 4. A numerical algorithm for problem (P2 ) Based on the optimality conditions from Section 3 we will develop here a conceptual algorithm to approximate the optimal control of problem (P2 ) (see also Aniţa et al. [16]). We first introduce an abstract Numerical Stabilization Method (NSM) for periodic differential equations. We consider the following abstract problem: find the unique positive solution of
(F h)(t ) = f (t ), h(nT +) = ρ h(nT −), h(t ) = h(t + T ),
t ∈ [0, +∞) \ T N n∈N t ∈ [0, +∞),
(S)
where f is T -periodic and F : X → Y is an appropriate operator between the appropriate spaces X and Y . The solution of the periodic problem (S ) may be computed via the Initial Value Problem (IVP)
(F h)(t ) = f (t ), h(nT +) = ρ h(nT −), h(0) = h0 ≥ 0,
t ∈ [0, +∞) \ T N n ∈ N∗ .
(S 0 )
The corresponding algorithm is presented below. Step 0: solve (S 0 ) for t ∈ [0, T −] and denote the numerical solution h(0) ; set n := 1; Step 1: solve
(F h)(t ) = f (t ), h(0+) = ρ h(n−1) (T −).
t ∈ (0, T )
(S n )
and denote the numerical solution h(n) ; Step 2: the stopping criterion if n = m or ∥h(n) − h(n−1) ∥ < ε then stop (the T -periodic extension of h(n) is the solution) else n := n + 1; go to Step 1. In Step 2 above m is a given positive integer and the norm should be an appropriate one, while ε > 0 is a given convergence parameter. For the time interval [0, T ] we consider a discretization grid with equidistant knots 0 = t1 < t2 < · · · < (n) (n) tN = T and we approximate the corresponding values of h(n) = (hi )i , i = 1, 2, . . . , N, as hi ≈ h(n) (ti ), for i = 1, 2, . . . , N. (n)
The norm above can be a discrete one with respect to h(n) = (hi )i . When the NSM stops, the T -periodic extension of h(n) is considered to be the solution of problem (S ). The system analyzed in this paper has the stabilization property required by the NSM. Independent of the initial value h0 the solution h is stabilized after a certain number of iterations. Now we present a Projected Gradient Method (PGM) for the optimal control problem (P2), using the mathematical results from Section 3. Taking into account the control restrictions, we use Rosen’s algorithm (e.g. Aniţa et al. [4]). We also use the NSM algorithm above as a subroutine. Step 0:set j := 0, v (j) (t ) := 0, compute the solution h˜ (j) to (11) corresponding to v := v (j) using the NSM algorithm; (j) ∗ ∗ Step 1: compute the solution q(j) to (13) (for h˜ v := h˜ (j) , u˜ v := u˜ v ) ( j ) using again the NSM algorithm, where u˜ v is the solution to (2)–(3) corresponding to v := v (j) ; ∗ Step 2: compute the solution θ (j) to system (15) (for q := q(j) , h˜ v := h˜ (j) ); (j) (j) ∗ Step 3: compute w according to formula (14) (where θ := θ , v := w (j) ); Step 4: compute the new control v (j+1) Step 4.1: compute the solution λj ∈ [0, 1] of the following problem Minimize {Φ (λv (j) + (1 − λ)w (j) ) ; λ ∈ [0, 1]}, where Φ is the cost functional corresponding to the control problem (P2 ), Step 4.2: compute v (j+1) = λj v (j) + (1 − λj )w (j) ; Step 5: the stopping criterion if j + 1 = jmax or |Φ (v (j+1) ) − Φ (v (j) )| < ε or ∥v (j+1) − v (j) ∥ < ε then stop (v (j+1) is the solution) else j := j + 1; go to Step 1.
12
L.-I. Aniţa et al. / Applied Mathematics Letters 40 (2015) 7–12
Fig. 1. Left: Plot of the optimal control v ∗ over [0, mT ]; Right: Plot of Φ (v (j) ).
In Step 5 above jmax is a given positive integer and the norm should be an appropriate one, while ε > 0 is a given convergence parameter. Numerical results. Since a convex combination of two bang–bang controls is not necessarily a bang–bang control, we have replaced the convex combination λj v (j) +(1−λj )w (j) from Step 4.2 of the PGM Algorithm. We have used in our computer program a convex combination of the switching points of v (j) and of w (j) obtaining thus a system of switching points for the new bang–bang control. The computer program to implement the PGM algorithm was written in Matlab. To compute the op timal control we have considered: r = 0.5, g = r /2, Ω k(x′ )dx′ = 0.95, K = 20, f (t ) = | sin t |, η = 0.1, T = π , ρ = 0.982, m = 10, M = 0.1, ζ2 = 1. We have made experiments for different values of ε , jmax, h0 , v (0) , q. For instance for ε = 0.01, ∗ ∗ jmax = 20, h0 ≡ 20, and v (0) ≡ 0 we have obtained v ∗ , hv and h˜ v over the time interval [0, mT ]. The stabilization of hv to h˜ v occurs after 4 or 5 T -periods for any control v (the same idea of approximation is used for the dual state q, starting from an arbitrary datum qf = q(mT −) and moving backward 4 or 5 T -periods). The computer program worked quickly, requiring only few iterations for the optimal control loop (see Fig. 1). Acknowledgments 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’’. The authors gratefully acknowledge useful comments made by the referees. References [1] S. Genieys, V. Volpert, P. Auger, Pattern and waves for a model in population dynamics with nonlocal consumption of resources, Math. Modelling Nat. Phenomena 1 (2006) 65–82. [2] S. Aniţa, V. Capasso, Stabilization of a reaction–diffusion system modelling a class of spatially structured epidemic systems via feedback control, Nonlinear Anal. RWA 13 (2012) 725–735. [3] V. Capasso, Mathematical Structures of Epidemic Systems, in: Lecture Notes Biomath., vol. 97, Springer–Verlag, Heidelberg, 2008, 2nd corrected printing. [4] S. Aniţa, V. Arnăutu, V. Capasso, An introduction to optimal control problems in life sciences and economics, in: From Mathematical Models to Numerical Simulation with MATLAB, Birkhauser, Berlin, 2010. [5] V. Barbu, Analysis and Control of Nonlinear Infinite Dimensional Systems, Academic Press, San Diego, 1993. [6] S. Lenhart, J.T. Workman, Optimal Control Applied to Biological Models, Chapman & Hall/CRC, Boca Raton, 2007. [7] J. Smoller, Shock Waves and Reaction-Diffusion Equations, Springer-Verlag, Berlin, Heidelberg, 1983. [8] S. Aniţa, M. Iannelli, M.-Y. Kim, E.-J. Park, Optimal harvesting for periodic age-dependent population dynamics, SIAM J. Appl. Math. 58 (1998) 1648–1666. [9] K.R. Fister, S. Lenhart, Optimal harvesting in an age-structured predator–prey model, Appl. Math. Optim. 54 (2006) 1–15. [10] C. Zhao, P. Zhao, M. Wang, Optimal harvesting for nonlinear age-dependent population dynamics, Math. Comput. Modelling 43 (2006) 310–319. [11] N. Hritonenko, Yu. Yatsenko, Bang–bang, impulse, and sustainable harvesting in age-structured populations, J. Biol. Syst. 20 (2012) 133–154. [12] N. Hritonenko, Yu. Yatsenko, R. Goetz, A. Xabadia, Optimal harvesting in forestry: steady-state analysis and climate change impact, J. Biol. Dyn. 7 (2012) 41–58. [13] A. Belyakov, V. Veliov, Constant versus periodic fishing: age structured optimal control approach, Math. Modelling Nat. Phenomena 9 (4) (2014) 20–37. [14] S. Behringer, T. Upmann, Optimal harvesting of a spatial renewable resource, J. Econ. Dynamics Control 47 (2014) 105–120. [15] Z.-R. He, R. Liu, Theory of optimal harvesting for a nonlinear size-structured population in periodic environments, Int. J. Biomath. 7 (4) (2014) 18. [16] L.-I. Aniţa, S. Aniţa, V. Arnăutu, Optimal harvesting for periodic age-dependent population dynamics with logistic term, Appl. Math. Comput. 215 (2009) 2701–2715.
Contents lists available at ScienceDirect
Applied Mathematics Letters journal homepage: www.elsevier.com/locate/aml
Optimal strategies to diminish a pest population via bilinear controls Laura-Iulia Aniţa a , Sebastian Aniţa b,c,∗ , Costică Moroşanu b a
Faculty of Physics, ‘‘Al.I. Cuza’’ University of Iaşi, Iaşi 700506, Romania
b
Faculty of Mathematics, ‘‘Al.I. Cuza’’ University of Iaşi, Iaşi 700506, Romania
c
‘‘Octav Mayer’’ Institute of Mathematics, Iaşi 700506, Romania
article
info
Article history: Received 24 July 2014 Received in revised form 4 September 2014 Accepted 4 September 2014 Available online 16 September 2014 Keywords: Population dynamics Bilinear control problems Large-time behavior Age-dependent harvesting rate Gradient-type algorithm
abstract In this work we model and analyze two control strategies to diminish a pest population using traps. The action of any trap depends on its age. Both problems contain bilinear control rates. The large-time behavior of the model with time-periodic inflow is investigated. The first control strategy deals with a finite horizon problem while the second one is related to a time-periodic problem. We obtain Pontryagin’s principle for both control problems. A special attention is given to the periodic problem. Pontryagin’s principle is used to derive a conceptual gradient-type algorithm to approximate the optimal solution. Numerical tests are given. © 2014 Elsevier Ltd. All rights reserved.
1. Mathematical modeling Let h(x, t ) denote the density of a pest population (of insects) at location x ∈ Ω of the habitat (Ω ⊂ R2 ), and time t ≥ 0. Assume that a ratio c˜ of the whole pest population may diffuse and the ratio (1 − c˜ ) of it may migrate. If the diffusion rate in the absence of migration is d˜ and if the proportion of pests emigrating from x to other locations in the absence of diffusion is g˜ (x), then if we denote by d = c˜ d˜ and by g = (1 − c˜ )˜g , the dynamics of h is described by the following model:
h ∂ h = d ∆ h + r ( x ) h 1 − − g (x)h t K (x) + k(x, x′ )g (x′ )h(x′ , t )dx′ + f (x, t ) − u˜ (x, t )χω (x)h(x, t ), Ω ∂ν h = 0, h(x, jT +) = ρ(x)h(x, jT −), h(x, 0+) = h0 (x),
x ∈ Ω , t ∈ R+ \ T N
(1)
x ∈ ∂ Ω , t ∈ R+ \ T N x ∈ Ω , j ∈ N∗ .
Here, r (x) = growth rate of pests at location x; K (x) = carrying capacity of pests at location x; d = diffusion constant of pests when migration occurs; g (x) = proportion of pests emigrating from x to other locations; k(x, x′ ) = proportion of pests immigrating to x from location x′ ; T = period of one year; f (x, t ) = density of pest inflow; u˜ (x, t ) = harvesting rate (which acts only in the open subset ω ⊂ Ω ); χω = characteristic function ofω. The pest population is subject to an inflow f (x, t ) of population coming from outside Ω . We assume here that the Neumann (meaning an isolated habitat with respect
∗
Corresponding author at: Faculty of Mathematics, ‘‘Al.I. Cuza’’ University of Iaşi, Iaşi 700506, Romania. Tel.: +40 749 33 65 06. E-mail addresses: [email protected] (L.-I. Aniţa), [email protected], [email protected] (S. Aniţa), [email protected] (C. Moroşanu).
http://dx.doi.org/10.1016/j.aml.2014.09.003 0893-9659/© 2014 Elsevier Ltd. All rights reserved.
8
L.-I. Aniţa et al. / Applied Mathematics Letters 40 (2015) 7–12
to diffusion) homogeneous boundary condition holds for the pest population; ∂ν h(x, t ) denotes the normal derivative depending on location x. We also assume that h0 (x) is the initial density of pest population at location x. Recall that basic properties of the solutions to some biological reaction–diffusion systems with nonlocal terms have been recently investigated in Genieys et al. [1] and Aniţa and Capasso [2]. We assume that a certain ratio of the pest population dies during the winters. This occurs in a different proportion α(x), at different locations. In conclusion, at the beginning of the new year, the pest density is h(x, jT +) = h(x, jT −) − α(x)h(x, jT −) = ρ(x)h(x, jT −) (ρ(x) = 1 − α(x)). For an arbitrary function w, we have denoted by w(x, jT +) = limε↘0 w(x, jT + ε), w(x, jT −) = limε↘0 w(x, jT − ε), if these limits exist. We shall work under the following hypotheses: (H1) Ω ⊂ R2 is a bounded domain with a sufficiently smooth boundary ∂ Ω , d, T ∈ (0, +∞); (H2) r , K , g ∈ L∞ (Ω ), r (x) ≥ r0 , K (x) ≥ K0 , g (x) ≥ 0, a.e. x ∈ Ω , where r0 , K0 are positive constants; (H3) k ∈ L∞ (Ω × Ω ), k(x, x′ ) ≥ 0 a.e. in Ω × Ω , f ∈ L∞ (Ω × (0, +∞)), f (x, t ) ≥ 0 a.e. in Ω × (0, +∞), u˜ ∈ L∞ (ω × (0, +∞)); (H4) ρ, h0 ∈ L∞ (Ω ), c ≤ ρ(x) ≤ 1, a.e. x ∈ Ω (where c ∈ (0, 1)), 0 ≤ h0 (x), a.e. x ∈ Ω . The existence, uniqueness and nonnegativity of a solution to (1) are derived as in Capasso [3]. Our goal is to diminish the pest population by usage of traps and keeping this in balance with the cost of the intervention (control). We propose here two strategies and analyze them. Notice that unlike the use of pesticides in order to diminish the pest population, the use of traps is safe for the human population. Let us also notice that a trap attracts and captures a certain ratio of the pest population which depends on the age of the trap (the time elapsed from the moment when the trap has been situated at a certain position x ∈ ω). Actually this ratio decays exponentially. The harvesting rate at position x and moment t is obtained by cumulating the actions at moment t of all traps situated at position x. During winters, the traps loose their attracting and catching capabilities. So, if v(x, t ) is the density of traps at position x ∈ ω and moment t ≥ 0, and uv (x, t , a) is the solution to the following age-structured model:
∂t u + ∂a u + ηu = 0, u(x, t , 0) = v(x, t ), u(x, jT +, a) = 0,
x ∈ ω, t ∈ R+ \ T N, a ∈ (0, T ) x ∈ ω, t ∈ R+ x ∈ ω, j ∈ N, a ∈ (0, T )
(2)
(here η > 0 is the degradation rate of the traps), then the harvesting rate at position x and moment t is u˜ v (x, t ) =
T
uv (x, t , a)da.
(3)
0
In what follows we assume (H5) v ∈ L∞ (ω × (0, +∞)), 0 ≤ v(x, t ) ≤ M, a.e. in ω × (0, +∞), where M is a positive constant (representing the maximal density of traps). It is obvious that u˜ v ∈ L∞ (ω × (0, +∞)) and so there exists a unique and nonnegative solution hv to (1) corresponding to u˜ := u˜ v . 2. An optimal control problem with finite horizon Assume that our goal is to diminish (using traps) the negative effect of the pest population over a time period [0, L], where L = mT , m ∈ N∗ . We wish to have a diminished pest population at the end of our intervention, at a reasonable cost. Denote V = {v ∈ L∞ (ω × (0, L)); 0 ≤ v(x, t ) ≤ M a.e. in ω × (0, L)}. For a fixed v ∈ V we denote by hv the solution to (1) corresponding to u˜ := u˜ v given by (2) and (3) (where we consider v ≡ 0 in ω × (T , +∞)). Actually, the first problem to be investigated is the following one:
(P1 ) min v∈V
L 0
Ω
hv (x, t )dxdt + ζ1
Ω
hv (x, L−)dx + ζ2
L 0
ω
v(x, t )dxdt ,
where ζ1 , ζ2 are positive constants. We may take smaller or bigger values for ζ1 and ζ2 . A big ζ1 means that we are mainly interested to have a small pest population at the moment L. Taking small ζ2 means that we are interested to significantly reduce the pest population no matter the cost and a big ζ2 means that we are interested to reduce the pest population but at a small cost (with a small effort). It follows in a standard manner (see Aniţa et al. [4], Barbu [5], Lenhart and Workman [6]) that there exists at least an optimal control v ∗ for (P1 ). Our main goal in this section is to derive the first order necessary optimality conditions. Actually, the following result holds: Theorem 1 (Pontryagin’s principle). Let v ∗ be an optimal control of (P1 ) and q be a solution to the following problem:
∗ 2rhv ∗ ∂t q = −d∆q − rq + q + gq − g (x) k(x′ , x)q(x′ , t )dx′ + u˜ v χω q + 1,
x ∈ Ω , t ∈ [0, L] \ T N
∂ q = 0, ν q(x, jT −) = ρ(x)q(x, jT +), q(x, L−) = −ζ1 ,
x ∈ ∂ Ω , t ∈ [0, L] \ T N x ∈ Ω , j ∈ {1, 2, . . . , m − 1}.
K
Ω
(4)
L.-I. Aniţa et al. / Applied Mathematics Letters 40 (2015) 7–12
9
Then
v (x, t ) = ∗
0, M,
if θ (x, t , 0) < ζ2 if θ (x, t , 0) > ζ2 ,
(5)
where θ is the solution to
∗ ∂t θ + ∂a θ − ηθ = q(x, t )hv (x, t ), θ (x, t , T ) = 0, θ (x, jT −, a) = 0,
x ∈ ω, t ∈ [0, L] \ T N, a ∈ (0, T ) x ∈ ω, t ∈ (0, L) x ∈ ω, j ∈ {1, 2, . . . , m}, a ∈ (0, T ).
(6)
Problem (4) is called the dual problem. The method to construct it can be found in Aniţa et al. [4]. Existence and uniqueness of a solution to (4) follows as in Aniţa et al. [4]. Proof. Let w ∈ L∞ (ω×(0, L)) be such that for any ε > 0 sufficiently small we have v ∗ +εw ∈ V (i.e. 0 ≤ v ∗ (x, t )+εw(x, t ) ≤ M a.e. in ω × (0, L)). Since v ∗ is an optimal control for (P1 ) we deduce that ∗ ∗ hv +εw − hv
L 0≤
ε
Ω
0
dx dt + ζ1
∗ ∗ hv +εw (x, L−) − hv (x, L−)
ε
Ω
dx + ζ2
L 0
ω
w(x, t )dx dt ,
(7)
for any ε > 0 sufficiently small. The following auxiliary result will be used in what follows: Lemma. The following convergences yield ∗ ∗ hv +εw → hv ,
∗ ∗ hv +εw − hv
→ z , in C ([nT +, (n + 1)T −]; L∞ (Ω )), ε for n ∈ {0, 1, . . . , m − 1}, as ε → 0, where z is the solution to ∗ 2rhv z − gz ∂ z = d ∆ z + rz − t K ∗ ∗ + k(x, x′ )g (x′ )z (x′ , t )dx′ − u˜ w χω hv − u˜ v χω z , x ∈ Ω , t ∈ [0, L] \ T N Ω x ∈ ∂ Ω , t ∈ [0, L] \ T N ∂ν z = 0, z (x, jT +) = ρ(x)z (x, jT −), z (x, 0+) = 0, x ∈ Ω , j ∈ {1, 2, . . . , m − 1}.
(8)
The proof of the lemma follows as in Aniţa et al. [4], Barbu [5], Lenhart and Workman [6]. Passing to the limit in (7) and get
L 0≤ Ω
0
z (x, t )dx dt + ζ1
Ω
z (x, L−)dx + ζ2
L ω
0
w(x, t )dx dt .
(9)
If we multiply the equation in (4) by z and integrate over Ω × (0, L) we obtain m−1
Ω
q(x, L−)z (x, L−)dx +
j =1
[q(x, jT −)z (x, jT −)dx − q(x, jT +)z (x, jT +)]dx
q(x, 0+)z (x, 0+)dx −
− Ω
L = Ω
0
Ω
z −d∆q − rq +
L
2rh
0 v∗
K
Ω
q(x, t )∂t z (x, t )dx dt
q + gq − g (x)
Ω
v∗
k(x , x)q(x , t )dx + u˜ χω q dx dt + ′
′
′
L 0
Ω
z (x, t )dx dt .
Taking into account that q(x, L−) = −ζ1 a.e. x ∈ Ω and that z satisfies (8) we get after an easy calculation that ζ1 L−)dx +
L
Ω
0
L ω
0
z (x, t )dx dt = T
w
v ˜w ω u (x, t )h (x, t )q(x, t )dx dt , and consequently by (9) and (3) we get that ∗
L 0
v∗
u (x, t , a)h (x, t )q(x, t )da dx dt + ζ2 0
L 0
ω
L
Ω
z (x,
w(x, t )dx dt ≥ 0.
If we multiply (6) by uw and integrate over ω × (0, L) × (0, T ) we get (taking into account (2) and (6)) that v∗
L
L T 0
ω 0
uw
(x, t , a)h (x, t )q(x, t )da dx dt = − 0 ω w(x, t )θ (x, t , 0)dx dt , and consequently 0 ω w(x, t )[ζ2 − θ (x, t , 0)]dx dt ≥ 0, for any w ∈ L∞ (ω × (0, L)) such that for any ε > 0 sufficiently small we have v ∗ + εw ∈ V . We finally conclude (5). 3. A periodic optimal control problem If we consider the domain Ω to be small, then it is reasonable to assume that the rates r (x), K (x), g (x), ρ(x), k(x, x′ ) are independent with respect to x, ω ≡ Ω (we may act in the whole domain), and the inflow f satisfies f ∈ L∞ (0, +∞),
10
L.-I. Aniţa et al. / Applied Mathematics Letters 40 (2015) 7–12
f (t ) = f (t + T ), f (t ) > 0 a.e. in (0, +∞). The seasonality of the inflow is also a realistic condition. Therefore it is natural to deal with space-independent and T -periodic controls. Let
VT = {w ∈ L∞ (0, +∞); w(t ) = w(t + T ), 0 ≤ w(t ) ≤ M a.e. in (0, +∞)}. For any v ∈ VT , problem (1) has a unique solution hv which is nonnegative. Let us prove first that limt →+∞ ∥hv (·, t ) − h˜ v (t )∥L∞ (Ω ) = 0, where h˜ v is the unique T -periodic positive solution to
h − gh + k(x′ )gh(x′ , t )dx′ + f (t ) − u˜ v (t )h, ∂ h = d ∆ h + rh 1 − t K
x ∈ Ω , t ∈ (0, +∞) \ T N
∂ν h = 0, h(x, jT +) = ρ h(x, jT −),
x ∈ ∂ Ω , t ∈ (0, +∞) \ T N x ∈ Ω , j ∈ N∗ ,
Ω
(10)
which actually is space independent. Hence h˜ v is the unique positive T-periodic solution to
h′ = rh 1 − h − gh + gh k(x′ )dx′ + f (t ) − u˜ v (t )h, K Ω h(jT +) = ρ h(jT −),
t ∈ (0, +∞) \ T N
(11)
j∈N . ∗
Using the comparison result for parabolic equations with nonlocal term (see Aniţa and Capasso [2]) we obtain that the sequence n → hv1 (t +nT +) is increasing for any t ∈ [0, T ), where hv1 is the solution to (1) corresponding to h0 ≡ 0. Actually this solution is space independent. On the other hand for any sufficiently large initial space independent value h02 ≥ ∥h0 ∥L∞ (Ω ) , we get via the same comparison result that n → hv2 (t + nT +) is decreasing, where hv2 is the solution to (1) corresponding to h0 := h02 . The same comparison result implies that 0 < hv1 (t ) ≤ hv (x, t ) ≤ hv2 (t ),
a.e. in Ω × (0, +∞).
v
(12)
v
v
The monotonicity of n → h1 (t + nT +) and of n → h2 (t + nT +) for any t ∈ [0, T ), implies that h1 (·, t ) − h˜ v1 (t ) → 0, hv2 (·, t ) − h˜ v2 (t ) → 0, in L∞ (Ω ) as t → +∞, where h˜ v1 and h˜ v2 are T -periodic solutions to (11) (and (10)). We use now an h˜ v
argument from Smoller [7]: we multiply the equation in h˜ v1 by ˜1v and the equation in h˜ v2 by ˜ v1 2 and integrate their differh2 (h2 )
ence on [0, T ]. Taking into account that 0 < h˜ v1 (t ) ≤ h˜ v2 (t ), a.e. in (0, +∞), we conclude that necessarily h˜ v1 ≡ h˜ v2 . Passing to the limit in (12) we may infer that limt →+∞ ∥hv (·, t ) − h˜ v (t )∥L∞ (Ω ) = 0, where h˜ v1 ≡ h˜ v2 ≡ h˜ v is the unique positive T-periodic solution to (11). For any v ∈ VT we have that (n+1)T
Ω
nT
hv (x, t )dx dt + ζ2
(n+1)T
nT
Ω
v(t )dx dt → meas(Ω )
T
h˜ v (t )dt + ζ2
0
T
v(t )dt ,
0
as n → +∞. Here meas(Ω ) is the Lebesgue measure of Ω . Therefore, a natural approach is that of minimizing the pest population along a time interval of length T (on long term) at a small cost via a T -periodic harvesting effort v . The last convergence implies that the most natural way to mathematically formulate this problem is
(P2 ) min
T
v∈VT
h˜ v (t )dt + ζ2
0
T
v(t )dt . 0
Existence of an optimal control v ∗ for problem (P2 ) follows in a standard manner. Moreover Theorem 2. Let v ∗ be an optimal control of (P2 ) and q be the T -periodic solution to the following problem:
v∗ ∗ q′ = −rq + 2r h˜ q + gq − gq k(x′ )dx′ + u˜ v q + 1, K Ω q(jT −) = ρ q(jT +),
t ∈ (0, +∞) \ T N
(13)
j∈N . ∗
Then
v (t ) = ∗
0, M,
if θ (t , 0) < ζ2 if θ (t , 0) > ζ2 ,
(14)
where θ is the solution to
v∗ ∂t θ + ∂a θ − ηθ = q(t )h˜ (t ), θ (t , T ) = 0, θ (jT −, a) = 0,
t ∈ (0, +∞) \ T N, a ∈ (0, T ) t ∈ (0, +∞) j ∈ N∗ , a ∈ (0, T ).
(15)
It is possible to prove that problem (13) has a unique solution using the same technique as for proving that problem (11) has a unique solution. For the proof we have to follow the steps in the previous section.
L.-I. Aniţa et al. / Applied Mathematics Letters 40 (2015) 7–12
11
Remarks. For the more complex situation when r , K , g , k, f , ρ are space dependent, we may however use comparison T T results in order to find a control v ∗ which provides a low upper bound for the cost functional 0 Ω h˜ v (x, t )dx dt + ζ2 0 Ω v(x, t )dx dt. Such result is based on the last theorem. For other recent results concerning some other optimal control problems with bilinear terms (with harvesting rates) in biology and economics, but of a different type we refer to Aniţa et al. [8], Fister and Lenhart [9], Zhao et al. [10], Hritonenko and Yatsenko [11], Hritonenko et al. [12], Belyakov and Veliov [13], Behringer and Upmann [14] and He and Liu [15]. 4. A numerical algorithm for problem (P2 ) Based on the optimality conditions from Section 3 we will develop here a conceptual algorithm to approximate the optimal control of problem (P2 ) (see also Aniţa et al. [16]). We first introduce an abstract Numerical Stabilization Method (NSM) for periodic differential equations. We consider the following abstract problem: find the unique positive solution of
(F h)(t ) = f (t ), h(nT +) = ρ h(nT −), h(t ) = h(t + T ),
t ∈ [0, +∞) \ T N n∈N t ∈ [0, +∞),
(S)
where f is T -periodic and F : X → Y is an appropriate operator between the appropriate spaces X and Y . The solution of the periodic problem (S ) may be computed via the Initial Value Problem (IVP)
(F h)(t ) = f (t ), h(nT +) = ρ h(nT −), h(0) = h0 ≥ 0,
t ∈ [0, +∞) \ T N n ∈ N∗ .
(S 0 )
The corresponding algorithm is presented below. Step 0: solve (S 0 ) for t ∈ [0, T −] and denote the numerical solution h(0) ; set n := 1; Step 1: solve
(F h)(t ) = f (t ), h(0+) = ρ h(n−1) (T −).
t ∈ (0, T )
(S n )
and denote the numerical solution h(n) ; Step 2: the stopping criterion if n = m or ∥h(n) − h(n−1) ∥ < ε then stop (the T -periodic extension of h(n) is the solution) else n := n + 1; go to Step 1. In Step 2 above m is a given positive integer and the norm should be an appropriate one, while ε > 0 is a given convergence parameter. For the time interval [0, T ] we consider a discretization grid with equidistant knots 0 = t1 < t2 < · · · < (n) (n) tN = T and we approximate the corresponding values of h(n) = (hi )i , i = 1, 2, . . . , N, as hi ≈ h(n) (ti ), for i = 1, 2, . . . , N. (n)
The norm above can be a discrete one with respect to h(n) = (hi )i . When the NSM stops, the T -periodic extension of h(n) is considered to be the solution of problem (S ). The system analyzed in this paper has the stabilization property required by the NSM. Independent of the initial value h0 the solution h is stabilized after a certain number of iterations. Now we present a Projected Gradient Method (PGM) for the optimal control problem (P2), using the mathematical results from Section 3. Taking into account the control restrictions, we use Rosen’s algorithm (e.g. Aniţa et al. [4]). We also use the NSM algorithm above as a subroutine. Step 0:set j := 0, v (j) (t ) := 0, compute the solution h˜ (j) to (11) corresponding to v := v (j) using the NSM algorithm; (j) ∗ ∗ Step 1: compute the solution q(j) to (13) (for h˜ v := h˜ (j) , u˜ v := u˜ v ) ( j ) using again the NSM algorithm, where u˜ v is the solution to (2)–(3) corresponding to v := v (j) ; ∗ Step 2: compute the solution θ (j) to system (15) (for q := q(j) , h˜ v := h˜ (j) ); (j) (j) ∗ Step 3: compute w according to formula (14) (where θ := θ , v := w (j) ); Step 4: compute the new control v (j+1) Step 4.1: compute the solution λj ∈ [0, 1] of the following problem Minimize {Φ (λv (j) + (1 − λ)w (j) ) ; λ ∈ [0, 1]}, where Φ is the cost functional corresponding to the control problem (P2 ), Step 4.2: compute v (j+1) = λj v (j) + (1 − λj )w (j) ; Step 5: the stopping criterion if j + 1 = jmax or |Φ (v (j+1) ) − Φ (v (j) )| < ε or ∥v (j+1) − v (j) ∥ < ε then stop (v (j+1) is the solution) else j := j + 1; go to Step 1.
12
L.-I. Aniţa et al. / Applied Mathematics Letters 40 (2015) 7–12
Fig. 1. Left: Plot of the optimal control v ∗ over [0, mT ]; Right: Plot of Φ (v (j) ).
In Step 5 above jmax is a given positive integer and the norm should be an appropriate one, while ε > 0 is a given convergence parameter. Numerical results. Since a convex combination of two bang–bang controls is not necessarily a bang–bang control, we have replaced the convex combination λj v (j) +(1−λj )w (j) from Step 4.2 of the PGM Algorithm. We have used in our computer program a convex combination of the switching points of v (j) and of w (j) obtaining thus a system of switching points for the new bang–bang control. The computer program to implement the PGM algorithm was written in Matlab. To compute the op timal control we have considered: r = 0.5, g = r /2, Ω k(x′ )dx′ = 0.95, K = 20, f (t ) = | sin t |, η = 0.1, T = π , ρ = 0.982, m = 10, M = 0.1, ζ2 = 1. We have made experiments for different values of ε , jmax, h0 , v (0) , q. For instance for ε = 0.01, ∗ ∗ jmax = 20, h0 ≡ 20, and v (0) ≡ 0 we have obtained v ∗ , hv and h˜ v over the time interval [0, mT ]. The stabilization of hv to h˜ v occurs after 4 or 5 T -periods for any control v (the same idea of approximation is used for the dual state q, starting from an arbitrary datum qf = q(mT −) and moving backward 4 or 5 T -periods). The computer program worked quickly, requiring only few iterations for the optimal control loop (see Fig. 1). Acknowledgments 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’’. The authors gratefully acknowledge useful comments made by the referees. References [1] S. Genieys, V. Volpert, P. Auger, Pattern and waves for a model in population dynamics with nonlocal consumption of resources, Math. Modelling Nat. Phenomena 1 (2006) 65–82. [2] S. Aniţa, V. Capasso, Stabilization of a reaction–diffusion system modelling a class of spatially structured epidemic systems via feedback control, Nonlinear Anal. RWA 13 (2012) 725–735. [3] V. Capasso, Mathematical Structures of Epidemic Systems, in: Lecture Notes Biomath., vol. 97, Springer–Verlag, Heidelberg, 2008, 2nd corrected printing. [4] S. Aniţa, V. Arnăutu, V. Capasso, An introduction to optimal control problems in life sciences and economics, in: From Mathematical Models to Numerical Simulation with MATLAB, Birkhauser, Berlin, 2010. [5] V. Barbu, Analysis and Control of Nonlinear Infinite Dimensional Systems, Academic Press, San Diego, 1993. [6] S. Lenhart, J.T. Workman, Optimal Control Applied to Biological Models, Chapman & Hall/CRC, Boca Raton, 2007. [7] J. Smoller, Shock Waves and Reaction-Diffusion Equations, Springer-Verlag, Berlin, Heidelberg, 1983. [8] S. Aniţa, M. Iannelli, M.-Y. Kim, E.-J. Park, Optimal harvesting for periodic age-dependent population dynamics, SIAM J. Appl. Math. 58 (1998) 1648–1666. [9] K.R. Fister, S. Lenhart, Optimal harvesting in an age-structured predator–prey model, Appl. Math. Optim. 54 (2006) 1–15. [10] C. Zhao, P. Zhao, M. Wang, Optimal harvesting for nonlinear age-dependent population dynamics, Math. Comput. Modelling 43 (2006) 310–319. [11] N. Hritonenko, Yu. Yatsenko, Bang–bang, impulse, and sustainable harvesting in age-structured populations, J. Biol. Syst. 20 (2012) 133–154. [12] N. Hritonenko, Yu. Yatsenko, R. Goetz, A. Xabadia, Optimal harvesting in forestry: steady-state analysis and climate change impact, J. Biol. Dyn. 7 (2012) 41–58. [13] A. Belyakov, V. Veliov, Constant versus periodic fishing: age structured optimal control approach, Math. Modelling Nat. Phenomena 9 (4) (2014) 20–37. [14] S. Behringer, T. Upmann, Optimal harvesting of a spatial renewable resource, J. Econ. Dynamics Control 47 (2014) 105–120. [15] Z.-R. He, R. Liu, Theory of optimal harvesting for a nonlinear size-structured population in periodic environments, Int. J. Biomath. 7 (4) (2014) 18. [16] L.-I. Aniţa, S. Aniţa, V. Arnăutu, Optimal harvesting for periodic age-dependent population dynamics with logistic term, Appl. Math. Comput. 215 (2009) 2701–2715.