The Dubovitskii and Milyutin Methodology Applied to an Optimal Control Problem Originating in an Ecological System

: We research a control problem for an ecological model given by a reaction–diffusion system. The ecological model is given by a nonlinear parabolic PDE system of three equations modelling the interaction of three species by considering the standard Lotka-Volterra assumptions. The optimal control problem consists of the determination of a coefﬁcient such that the population density of predator decreases. We reformulate the control problem as an optimal control problem by introducing an appropriate cost function. Then, we introduce and prove three types of results. A ﬁrst contribution of the paper is the well-posedness framework of the mathematical model by considering that the interaction of the species is given by a general functional responses. Second, we study the differentiability properties of a cost function. The third result is the existence of optimal solutions, the existence of an adjoint state, and a characterization of the control function. The ﬁrst result is proved by the application of semigroup theory and the second and third result are proved by the application of Dubovitskii and Milyutin formalism.


Introduction
In recent decades, there has been an increasing interest in the mathematical modelling of several biological phenomena: pattern formation, epidemic disease transmission, blood circulation, atherosclerosis, species competency, migration of species, and so on, see for instance [1][2][3][4][5][6]. In particular, in this paper we are interested in a reaction-diffusion system with initial and boundary value conditions of the following type ∇u · ν = ∇v · ν = 0, on Σ T = (0, T) × ∂Ω, (1d) where Ω ⊂ R d (d = 1, 2, 3) is a bounded set with boundary ∂Ω; T > 0 is a fixed time denoting the total time of the process; the coefficients d 1 , d 2 , γ 1 , γ 3 , α i , b i (i = 1, 2, 3) are all positive constants; f , g : R 2 → R are some given functions modelling the functional responses or interactions; ν is the outward unit normal to ∂Ω, the boundary of Ω; and u 0 , v 0 , w 0 are some given functions from Ω to R + 0 modelling the initial conditions. The system (1) arose in the mathematical model of an ecological system under the considerations of Lotka-Volterra competence assumptions for three species, for instance u, v and w denotes the population densities of herbivorous, carnivorous and plants or pest, predator and plants, respectively. The boundary conditions (1d) are considered to model the case of an isolated environment, i.e., there is neither immigration nor emigration of species during the process. The maintenance of an adequate number of individuals of each species in the ecosystem is crucial to preserve the harmony and stability between living beings and the environment in which they inhabit. Therefore, understanding the mechanism of control of species is necessary to maintain the equilibrium of the ecosystem. For instance in the case of agriculture, the control of density for pest individuals can be done by incorporation of a predator and by using a chemical pesticide. Now, the usage of pesticides can dangerous for some species in the environment and should be optimized. In particular, considering mainly the following two assumptions: the density of eradicated pests is proportional to the density of living pests and the pesticide is uniformly distributed in the ecosystem, and to the determination of the proportional factor introducing the following control problem ∇u · ν = ∇v · ν = 0, on Σ T , (2d) where ζ from [0, T] to [0, 1] is the control function. The control problem (2) can be recast as an optimal control problem. Let us consider the notations: H = L 2 (Ω) 3 ; H the topological dual of H; and W 1,2 (0, T; H) the set of functions such that f ∈ L 2 (0, T; H) and f t ∈ L 2 (0, T; H ), which is a Banach space with the norm f W 1,2 (0,T;H) = f L 2 (H) + f t L 2 (H) . Moreover, we consider the notation E = W 1,2 (0, T; H) 3 × L 2 (0, T), E = L 2 (Q) 3 For a major details on Sobolev space we refer to [7]. We define the cost functional J : E → R as follows We notice that the cost function J is constructed in a appropriate sense, for instance in the case of pest-predator-plants ecosystem, we have that the first term minimizes the total density of pests and maximizes the total density of plants, and the second term is introduced to reduce the exposition to pesticides. Thus, we observe that (2) can be rewritten as the optimal control problem Find the control function ζ and the state variables of the ecosystem (u, v, w), such that the cost functional J given on (4), subject to (u, v, w, ζ) ∈ E solution of (2), is minimized,    (5) or equivalently, in the context of Dubovitskii and Milyutin formalism, as the generic optimization problem where the operator M : E → E is defined by if and only if We observe that the system (2) is a particular case of the system (8) when ψ i = 0 for i = 1, . . . , 6.
The main contributions of this paper are the introduction of appropriate assumption and a functional framework such that we can prove three kinds of results: (i) the existence and uniqueness of a positive solution of system (2) (see Theorem 1); (ii) the explicit calculus of descent and dual cones of J and other differentiability properties of the operators M defined on (7) and (8) (see Lemma 1); and (iii) the existence of solutions for (5), the existence of solutions for the adjoint system for system (2), and the characterization of control function (see Theorem 2).
We remark that there are several optimal control problems for reaction-diffusion equations in the recent literature, for instance . Some relevant results for semilinear parabolic differential equations are presented in [26]. In [14,15,22] the authors study the solution of inverse problems related with the reconstruction of coefficients in reactiondiffusion systems arising in epidemiology, from observations of the state solution in a fixed time, by application of some results of standard optimal control theory. A closer study of the problem (5) is the analysis introduced by [10] for the particular case of functional responses given by f (u, v) = uv and g(u, w) = uw. In a broad sense, the author of [10] assume the dependence of the state variables on a given ζ and instead of (5) she study the optimization problem where the notation (u ζ , v ζ , w ζ ) to emphasize the dependence of the state variables on ζ. The cost functional J is called the reduced cost functional [26]. Other relevant works for optimal control problems for Lotka-Volterra-like systems using the reduced form are [20,28]. Moreover, the application of Dubovitskii and Milyutin formalism to control problems arising in heat and solidification phenomena was conducted in [13,21], respectively. However, to the best of our knowledge, there is not yet in the literature an application of Dubovitskii and Milyutin to study optimal control problems in reaction-diffusion systems arising in the competition of organisms or species in an ecosystem.
The article is organized as follows. In Section 2 we introduce some additional notation; we precise the mathematical assumptions used on the paper for the physical domain, the coefficients, the functional responses, the control function and the initial conditions; and we introduce the statements of main results. In Section 3 we recall some useful results related with differential equations on Banach spaces and the Dubovitskii and Milyutin formalism. In Section 4 we introduce the proofs of main results. Finally, in Section 5 we give the conclusions and the guidance for future work.

Assumptions
Hereinafter, we consider the following assumptions: is a bounded domain with a boundary of class C 2+α with α > 0. Assumption 2. The coefficients d 1 , d 2 , α 1 , α 2 , α 3 , γ 1 , and γ 3 are strictly positive. Assumption 3. The initial condition (u 0 , v 0 , w 0 ) is belong the set U 0 defined as follows and also the functions u 0 , v 0 and w 0 are strictly positive on Ω.

Statements of Main Results
Theorem 1. Consider that the Assumptions 1-5 are satisfied, then there is at most one strictly positive global strong solution of the system (2) belong to W 1,2 (0, T; L 2 (Ω) 3 for a. a. t ∈ [0, T] and for some generic positive constant C (independent of (u, v, w) and ζ).

Lemma 1.
Consider that the hypotheses of Theorem 1 are satisfied, and the operators J and M defined on (4), (7) and (8), respectively. Then the following assertions are satisfied (a) The sets defines the descent and dual cones of J. (b) The application M is Gâteaux differentiable and the derivative of M in (u, v, w, ζ) is defined by M G (u, v, w, f )(û,v,ŵ,ζ) = (ψ 1 ,ψ 2 ,ψ 3 ,ψ 4 ,ψ 5 ,ψ 6 ) if and only if The application M is strictly differentiable and M (u, v, w, ζ) = M G (u, v, w, ζ) is a surjective operator.

Theorem 2.
Consider that the Assumptions 1-4 are satisfied. Then, the following assertions are satisfied (a) The optimization problem (6) has at least one solution. (6), then the adjoint system has at least one solution such that (p, q, r) ∈ W 1,2 (0, T; H). (c) Let (ū,v,w,ζ) ∈ D a solution of (6) and (p, q, r) a solution the adjoint system (16). Then, the relationζ is a characterization of the control functionζ on [0, T].

Remark 1.
The switching function Ω (ūp)(t, x)dx + 1 considered in (17) depends of the optimal control solution. However, we observe that it can rewritten in a more simply sense by getting precise bounds of the controlled and adjoint systems. For instance, one way to get that is detailed below. From Theorem 1, we have thatū is strictly positive and bounded on Q T . To fix ideas, let us consider the positive constantsū m andū M such that 0 < u m ≤ū(x, t) ≤ u M on Q T , i.e., we have the bounds of Then, if we were able to deduce upper and lower bounds for Ω p(t, x)dx, we would rewrite (17).

Preliminaries
For completeness of the presentation, we recall some results related with the theory of differential equations on Banach spaces, the Dubovitskii and Milyutin formalism and differential calculus on Banach spaces. Theorem 3 ([10], Theorem 2.1). Let X a Banach space; the linear operator A : D(A) ⊂ X be the infinitesimal generator of a C 0 -semigroup of contractions on X; and F : [0, T] × X → X be a measurable function in t and Lipschitz in x ∈ X, uniformly with in t ∈ [0, T]. Then, for any y 0 ∈ X, the following Cauchy problem has a unique mild solution belong C([0, T]; X). Moreover, if X is a Hilbert space, A is a selfadjoint dissipative operator on X, y 0 ∈ D(A), the mild solution is a strong solution and is belong To introduce the results of the differential calculus on Banach Spaces, we consider the normed vector spaces X and Y, U a neighborhood of x 0 ∈ X, and F : U ⊂ X → Y an application, for more details we refer to the book of Brezis [29]. Definition 1. We say that F has a derivative in the direction h ∈ X at x 0 if there exists the Y-limit Definition 2. The first variation of the application F at x 0 is the application δF( Definition 3. Let us denote by Λ ∈ L(X, Y) the space of linear continuous operators from X to Y. Consider that there exists Λ ∈ L(X, Y) satisfying the relation δF(x 0 , h) = Λh. We call to Λ the Gâteaux derivative of F at x 0 and is appropriately denoted by Definition 4. Consider that F satisfy the following relation at some neighborhood of x 0 . Then the application F is called Fréchet differentiable at x 0 and Λ, usually denoted by F (x 0 ), is called the Fréchet derivative (or only the derivative) of the application F in x 0 .

Definition 5.
The application F is called strictly differentiable at x 0 if there exists the operator Λ ∈ L(X, Y) such that for all x 1 , x 2 ∈ X satisfying the restrictions The terminology and results related with the Dubovitskii and Milyutin formalism are presented below, for major details consult [30,31]. Firstly, we consider the generic optimization problem where J from X to R is a functional and M from X to Y is an operator, with X and Y Banach spaces. Definition 6. The vector h belong the Banach space X is said a descent direction of the functional J : X → R at x 0 ∈ X if there is a neighborhood U of h and α = α(J, x 0 , h) > 0 such that the inequality J(x 0 + εh) ≤ J(x 0 ) − εα is satisfied for all ε ∈ (0, ε 0 ) and for any h ∈ U. Additionally, we said that the functional J is regularly decreasing at x 0 ∈ X if the set of all descent directions at x 0 is a convex set.

Definition 7.
Let h belong the Banach space X and Q i with int(Q i ) = ∅ the set defining the i-th inequality restriction. The vector h is said a feasible direction at x 0 ∈ X if x 0 + εh ∈ Q i for all ε ∈ (0, ε 0 ) and for any h ∈ U with U of a neighborhood of h. Additionally, the inequality restriction set Q i is called regular at x 0 ∈ X if the set of all feasible directions on Q i at x 0 is a convex set.

Definition 8.
Let h belong the Banach space X and Q i with int(Q i ) = ∅ the set defining the i-th inequality restriction. The vector h is said a tangent direction at . Additionally, if the the set of all possible tangent directions is a vectorial subspace it is said a tangent space, and also the inequality restriction Q i is called regular at x 0 if the set of all possible tangent directions for Q i at x 0 define a convex set.

Definition 9.
A set K ⊂ X, with X a real Banach space, is said a cone with vertex at zero if λx ∈ K for all λ > 0 and x ∈ K. Moreover, we call the dual cone for K to the set denoted by K * and defined as follows Theorem 4 (Dubovitskii and Milyutin theorem). Consider the optimization problem (19). Assume that J has a local minimum at x 0 ∈ Q = n+1 i=1 Q i , J is regularly decreasing at x 0 , with descent directions cone K 0 , Q i , i = 1, ..., n, are regular at x 0 , with feasible directions cone K i , and Q n+1 is regular at x 0 , with tangent directions cone K n+1 . Then, there exist n + 1 continuous linear functionals G i ∈ K * i , not all identically zero, such that Theorem 5 ([32]). (Lyuternik theorem) Let X and Z be real Banach spaces with norms · X and · Z , respectively; H : X → Z a given mapping; and S := {x ∈ X : H(x) = 0}. If H satisfy the following assumptions: h is Fréchet differentiable on a neighborhood of x, H (·) be continuous on x, and H (x) be surjective, then {x ∈ X : H (x)(x) = 0} ⊂ T(S, x).
However, we cannot directly apply the Theorem 1 since F does not satisfies the Lipschitz condition. Then, to study the existence and uniqueness for (18) for the particular case of functions and spaces given in (20), we proceed in three steps: we prove the existence and uniqueness of a positive local solution, we prove that the local solution is a global solution, and we deduce the estimates (10)-(12).
Step 1: Local solutions via a truncated problem. Let us consider the truncated Cauchy problem where N > 0 and F(t, N), with N = (N, N, N).
We observe that the Assumption 4 implies that F N satisfies the Lipschitz hypothesis required by Theorem 3. Then we can deduce that the Cauchy problem (21) has a unique strong solution belong W 1,2 ([0, T]; X).
We prove that the solution of (21) is bounded as follows. Let us consider the uncoupled Cauchy problems wherê The strong solution of (22) satisfy the relation Using the definition of M andŷ 0 , we deduce that each component of y 0 −ŷ 0 and F N (s, y N (s)) − M for any s ∈ [0, T] are negative. Then, from (23) Then, y N ∈ L ∞ (Q T ) and the upper bound is MT + u 0 L ∞ (Ω) which clearly does not depends of N.
On the other hand, we observe that y N satisfy the system Let us consider y N and y ⊥ N defined by y N,i (t, x) = sup{y N,i (t, x), 0} and y ⊥ N,i (t, x) = − inf{y N,i (t, x), 0} for i = 1, 2, 3, such that y N = y N − y ⊥ N . Testing (25a) by y ⊥ N , integrating by parts on [0, t] × Ω, and using the boundary conditions (25b), we have that Now, using the Hypotheses (A2), (A3) and (A4), we deduce the following bounds for some positive constants C i , i = 1, 2, 3. Then, by application of the integral form of Gronwall's inequality, we have that Ω |y ⊥ N,i (x, t)| 2 dx ≤ 0 for i = 1, 2, 3, i.e., y ⊥ N (x, t) = 0 on [0, t] × Ω or y N,i = y N,i ≥ 0, i = 1, 2, 3 on [0, t] × Ω. Thus, from the arbitrariness of t ∈ [0, T] and the fact that the components of y 0 are positive on Ω, we get that the components of y N are strictly positive on Q T .
To conclude the step 1, we deduce the existence and uniqueness of a positive local solution of (2) as follows. We observe that, if we select N > 2 y 0 L ∞ (Ω) , we get that Mθ + y 0 L ∞ (Ω) ≤ N/2 for some θ ∈ [0, T]. Now, from the estimate (24) we deduce that |y N,i (x, t)| ≤ N on [0, θ] × Ω for i = 1, 2, 3. Then, from the definition of F N we have that F N = F for t ∈]0, θ[ and y N is a solution of (2) on [0, θ] × Ω.
Step 2: The local solution is a global solution. To prove that the local solution y N on [0, θ] × Ω is a global solution it suffices to prove that y N is bounded on [0, θ] × Ω. Indeed, from (2c), the positivity of u, v on Q T deduced on Step 1 together with the positivity of g on R 2 + given by Assumption 4, the strictly positivity of α 3 assumed in Assumption 2, and the fact that and w 0 is strictly positive on Ω (considered on Assumption 3), we deduce that 0 < w(t, x) ≤ w 0 (x) exp(α 3 t) for (x, t) ∈ Q T . Then, w ∈ L ∞ ((0, θ) × Ω).
Step 3: Estimates (10)- (12). Squaring both sides of (2a) and applying an integration by parts on (0, t) × Ω, we get Then using the fact that u, v are bounded and positive and the Assumptions 3-5, we follow the estimate (10). The estimates (11) and (12) can be proved by similar arguments, starting from (2b) and (2c), respectively.

Proof of Lemma 1
Proof of item (a). From definition of J given in (4), we deduce that J is linear and consequently we get We notice that J is a convex and continuous operator. Then from ( [30], Theorem 7.3), we get that J is regularly decreasing for all (û,v,ŵ,ζ) ∈ E and the set is the descent cone for J. Moreover, from result of ( [30], pp. 69) we obtain that the set is the dual dual cone of K 0 .
To prove that the application M(·, ·, ·, ·) is strictly differentiable is sufficient to verify that the function (u, v, w, ζ) → M (u, v, w, ζ) is continuous or equivalently that it will suffice to prove that M (u, v, w, ζ) is bounded in E , since M (u, v, w, ζ) is linear. From (13), we get , which implies that ψ i L 2 (Q T ) for i = 1, 2, 3 are bounded by application of estimates (10)- (12) and (32)-(34). Then, by using the norm of E we have that Therefore M (u, v, w, ζ) is bounded in E and thus the proof is complete.
Proof of item (d). From Lemma 1-(b), we have that the application M is is Gâteaux differentiable and in a neighborhood of (û,v,ŵ,ζ) and by Lemma 1-(c) we have that M (û,v,ŵ,ζ) is continuous in a neighborhood of (û,v,ŵ,f ) and surjective. Consequently, the hypothesis of Theorem 5 are satisfied and the tangent cone to D at (u, v, w, ζ) is the kernel of the differential operator M (u, v, w, ζ) given on (14). Moreover, by linear algebra result we deduce that TC(D, u, v, w, ζ) is a vector space, since the kernel of a linear operator is a vector space. Now, the proof of (15) is deduced by the definition.

Proof of Theorem 2
Proof of item (a). We can develop the proof of existence of at least one solution (u, v, w, ζ) for the optimization problem (6) via the standard method of minimizing sequences. Let us consider that {(u n , v n , w n , ζ n )} n∈N a sequence in D, such that ∂ t u n − d 1 ∆u n = α 1 u n − β 1 f (u n , v n ) + γ 1 g(u n , w n ) − ζ n u n , in Q T , (35a) From application of Theorem 1 we deduce the following bounds for a. a. t ∈ [0, T] and for some generic positive constant C (independent of u n , v n , w n , ζ n and n). We observe that the sequence {(u n , v n , w n )} n∈N is uniformly bounded in W 1,2 (0, T; L 2 (Ω) 3 ) 3 and in C([0, T]; L 2 (Ω)) 3 . The compactness of {(u n , v n , w n )(t, ·)} n∈N in L 2 (Ω) for any t ∈ [0, T] is proved as follows. From (35c) we deduce that for any t ∈ [0, T], which implies that there is C > 0 such that Then by the Ascoli-Arzela Theorem there is at least one subsequence of w n (also labeled by n ) and there are w such that w n → w in L 2 (Ω) uniformly in t ∈ [0, T]. Now, from the compact embedding of H 1 (Ω) in L 2 (Ω), we deduce that {(u n , v n )(t, ·)} n∈N is compact in L 2 (Ω) 2 for any t ∈ [0, T]. Then by the Ascoli-Arzela Theorem {(u n , v n )} n∈N is compact in C([0, T]; L 2 (Ω)) and there is at least one subsequence of (u n , v n ) (also labeled by n ) and there are (u, v) Hence, we have that there is subsequence {(u n , v n , w n , ζ n )} n∈N ⊂ D (also labeled by n ) and there is (u, v, w, ζ) ∈ D, such that We notice that ζ ∈ Z and by the convergence properties of the sequence (u n , v n , w n ), we can take the limits in the system (35) and deduce that (u, v, w, ζ) satisfy a system associated withζ. Using the lower semi-continuity of cost functional we deduce that J(u, v, w,f ) ≤ lim n→∞ inf J(u n , v n , w n , f n ).
Thus, (u, v, w,f ) ∈ D is a solution of the problem (6).
Then by Theorem 4, we have follow the existence of two continuous functionals G 1 ∈ [DC(J)] * and G 2 ∈ [TC(D)] * , not both identically zero, satisfying the Euler-Lagrange equation Now, in order to prove the relation (17) we consider arbitrarily ζ ∈ Z and assume that (u , v , w , ζ ) is a solution of the following system The existence of solutions of (39) can be developed similarly to the system (37). From Lemma 1, we follow that (u , v , w , ζ ) ∈ TC(D, (ū,v,w,ζ))) and consequently G 2 (u , v , w , ζ ) = 0. Now, from (38) and some λ ≥ 0, we follow that We note that λ = 0, since if we assume that λ = 0, we have that G 1 = 0 and from (38) we deduce that G 2 = 0, which is a contradiction with the fact deduced by Theorem 4: there exist continuous functionals G 1 ∈ [DC(J)] * and G 2 ∈ [TC(D)] * , not both identically zero. Then, from (40) we get by dividing by λ and also without loss of generality we can fix λ = 1.

Conclusions and Future Work
This paper presents a the application of Dubovitskii and Milyutin formalism to study a control problem for a reaction-diffusion system arising in the competency of three species in a bounded ecosystem. The control problem is reformulated as an optimal control problem by incorporating a cost functional which maximizes the total density of beneficial species to the ecosystem and minimizes the pests and the exposition to the human intervention by incorporation of some substances to control pests. The main assumptions considered in order to deduce the results are the following: the ecosystem is modeled by a set which is bounded and has sufficiently smooth boundaries, the coefficients and initial conditions are strictly positive, the initial conditions satisfies a consistent relationship with the boundary conditions, the functions modelling the interaction of species are assumed positive and locally Lipschitz, and the admissible control set is assumed to be the square integrable functions. The main results obtained in the paper are: (i) the existence of a strictly positive global solution for controlled system which is deduced by using the theory of differential equations on Banach spaces; (ii) the construction of descent and dual cones of cost function and the differentiability of the operator defining the restriction set of the optimization problem is based on the application of differential calculus on Banach spaces; (iii) the existence of of at least one optimal control solution; (iv) the existence of solutions for the adjoint system; and (v) a characterization of control function by a function of "bang-bang" type with the switching function depending of the controlled and adjoint systems.
In our future work, we plan to extend the present research in at least three ways as follows. First, we plan to study other cost functions, for instance the incorporation of observations and/or attainable states in the L2-norm. Second, we will study the construction of some explicit bounds for the controlled and adjoint systems, such that the switching function can be determined only in terms of coefficients, initial conditions and the geometry of the ecosystem (see Remark 1). This is an important issue to be researched in order to develop numerical simulations of the control problem without solving the optimization problem. Thus, an another idea to study the switching function by follow similar ideas to those given in [33]. Third, we plan to develop a numerical study of the optimal control problem by applying a finite volume method and make some simulations of published experimental results.