Advection control in parabolic PDE systems for competitive populations

This paper investigates the response of two competing species to a given resource using optimal control techniques. We explore the choices of directed movement through controlling the advective coefficients in a system of parabolic partial differential equations. The objective is to maximize the abundance represented by a weighted combination of the two populations while minimizing the 'risk' costs caused by the movements.

1. Introduction. Population movement and its distribution in response to its surrounding environment are important concerns in ecology. Many efforts are devoted to analyze the population dynamics under different conditions affecting the population [1,2,3,4,5,6,7,9,11]. Particularly, to capture, and explain more realistic dynamics of the populations as well as the persistence of species in the long run, partial differential equations (PDE) of reaction-diffusion type with linear and nonlinear growth terms have been used in many studies [1,2,4]. The effect of the resources on the population size was also investigated in [3,7,11] using such models. For these reaction-diffusion models, the spatial movement of a population can occur in two different ways: random diffusion and directed movements. The response of the population to the surrounding environments depends on the choice of diffusion and advection coefficients in the model [1,2,5,6,9].
In the inhomogeneous environments, it is commonly believed that the population will move along the direction of increasing resources. With that regard, Belgacem and Cosner [1] investigated an advection reaction-diffusion model in which the advection term is the gradient of the resource function. They investigated steady states solutions of the reaction-diffusion equations with linear and nonlinear logistic growth terms: together with no flux boundary condition or Dirichlet boundary condition (lethal boundary). In these PDEs, the coefficient m(x) for x ∈ Ω represents the intrinsic growth rate, and measures the availability of the resources. Moreover, u(x, t) is the population density of the species at location x ∈ Ω and time t ≥ 0. They studied the benefit to the population, [1] meaning the persistence of the population in the long run or the existence of a unique globally attracting positive steady state solution. An interesting observation in [1] was that the directed movement towards better resources could be beneficial or harmful to the population depending on the boundary conditions. In the case of no flux boundary conditions, the movement towards better resources could be beneficial to the population, while in the lethal boundary conditions, the movement towards better resources can be harmful if more favorable patches are closer to the boundary. In a further investigation [6], Cosner et al. studied the logistic reaction-diffusion models with the advection along the environmental gradient, with the no flux boundary conditions. Different from the results in [1], they found that even under the zero flux boundary conditions, the movement along the resource gradient may not always be beneficial to the population. Indeed, it turns out that the convexity of the domains plays a major role in this situation. If the domain is not convex, moving up along the resource gradient could be harmful to the population. In a similar perspective as in [1], system of competing species have also been investigated. In particular, in [5], Lou et al. studied the persistence of two populations in a heterogeneous environment using Lotka-Volterra advection reaction-diffusion equations. They investigated two different movement strategies for the two completive populations: one population only moves with random diffusion, while the other population is allowed to have both the random diffusion, and the directed movement along the resource gradient. A related work [6] shows that when the advection is strong, both competing populations can stably coexist.
In an alternate but related viewpoint, optimal control techniques were used in various work such as [7,9,11] to explore how different conditions such as limited resources, growth coefficient, advection movement, and harvesting can be optimized to be "beneficial" for populations. In particular, Ding et al. studied in [7] the effects of resource allocation on population size of the species. They investigated the problem in which the resource function m(x) is the control, and the population abundance is a measurement of conservation effort. In other words, the work [7] explores the optimum resource allocation which maximizes the total population with the minimum cost for the resources.
In a similar framework as in [7], but different direction, Phan et al. [9] considered a general model of nonlinear advection reaction-diffusion equation. Precisely, they studied the following nonlinear parabolic equation: with Ω, m(x, t) and u(x, t) as before. Moreover, µ > 0 is the diffusion rate, h : Q T → R n is the advection direction, and f : Q T × R → R is a non-negative function satisfying some natural smoothness and growth conditions. In their study, the advection direction h is the control, and they seek for the optimal advection direction h(x, t) that maximizes the total population while minimizing the "cost" due to movement. In this work, the "cost" represents the risk of moving to a new location. The numerical results in [9] indeed indicate that the optimal control h * in some cases approximates the spatial gradient of the resource function m(x, t). We will discuss this issue further in the numerical results of our paper. Following the line of the research directions in [1,5,6,7,9], we move further on this work using optimal control techniques to investigate two competing species, and explore how the advection directions are chosen to maximize the total population while minimizing the costs related to the risk of movement. Our goal is to analyze this playoff between benefits and cost using the optimal control techniques.
In Section 2 of the paper, we introduce the systems of advection reaction-diffusion of the two competitive species that we will study. The set up for the optimal control problem will be also presented in this section. Section 3 of the paper discusses the existence of state solution, and the existence of the optimal control. In section 4, we derive the optimality system. Section 5 characterizes the optimal control, and proves the uniqueness of this optimal control. In section 6, we give some numerical results for the case of one space dimension. This paper finishes with some conclusions and some possible future research directions.
2. The optimal control problem formulation. Let Ω be an open bounded domain in R d with smooth boundary ∂Ω, for d ∈ N. For a given fixed time 0 < T < ∞, we denote Q = Ω × (0, T ), and the lateral boundary of Q is S = ∂Ω × (0, T ). In this paper, we consider two competing species in the spatial domain Ω. We denote u(x, t), v(x, t) their population densities with x ∈ Ω and t ∈ [0, T ). We are interested in the following parabolic system of advection reaction-diffusion equations: with the non-flux boundary condition and the initial data condition where u 0 , v 0 are given sufficiently smooth functions, d k are given positive constants, and a k , b k are given non-negative constants for k = 1, 2. Moreover, vector functions which will be specified, m : Q → R is a given measurable function, and ν is the unit outward normal vector field on ∂Ω. Biologically, in the system (1), both populations are allowed to move along the advections h 1 , h 2 , in addition to random diffusion with rate d 1 , d 2 . The constants a k , b k are the competition rates within one population or between two populations, respectively. The function m = m(x, t) represents the resources/trap. In this paper, we always assume that m ∈ L ∞ (Q).
To set up our optimal control problem, let M 1 , M 2 be two fixed positive numbers. Viewing the advection as the controls, we define our control set We assume both populations can choose their advective directions h 1 and h 2 to maximize their abundance, and minimize the "cost" due to the risk in the movements. Therefore, for given non-negative constants A, B, C, D, we seek to find, and where J( h 1 , h 2 ) is the objective functional defined by is the solution of the system (1)-(3). Our goal is to maximize a weighted sum of the two populations over time and space while trying to minimize the total cost involved in this movement subject to the PDE system (1)- (3). In this case, the nonlinear cost is due to the "risk" in moving the two populations by advection as the movement to a new location may come with a risk.
3. Solving the optimal control problem. We first recall the definition of weak solutions.
Definition 3.1. We say that u, v is a weak solution of the system (1)- We have the following lemma based on the energy estimate. such that for a given control Proof. Taking u , v as test functions and considering the weak formulation of (1)-(3) with the standard techniques of using Hölder's inequality, Young's inequality, and Gronwall's inequality we can get the above a priori estimates, see for example [8].
In order to solve the optimal control problem (4), we first need to to show that for a given control ( h 1 , h 2 ) ∈ U there exists a unique weak state solution (u, v) = (u, v)( h 1 , h 2 ) of the system (1)-(3). We indeed have the existence and uniform roundedness of the two solution as in the following theorem.
Proof. This proof uses a standard iterative method. First, we construct supersolutions of our system. Let us denote U 1 and V 1 be the unique non-negative weak solutions of the following linear equations: Moreover, there are constants C 1 , C 2 > 0 depending on M 1 , M 2 , |Q|, d 1 , d 2 and L ∞norms of the resource function m, and of the initial conditions such that For i = 2, 3, ..., consider the following iterative system of PDE's with u 1 = U 1 and v 1 = 0 on Q.
The weak solution for (8) exists by the general results from [8], where, Choosing u 1 = U 1 and v 1 = 0, Using our supersolution V 1 satisfying L 2 (V 1 ) = mV 1 +RV 1 , mathematical induction, and maximum principle give Moreover, again by induction, the sequence {u i } is monotone decreasing, and the sequence {v i } is monotone increasing, i.e.
using the maximum principle, and that F (u, v) is increasing in u and decreasing in v, and . Now, we can show that there exists a constant K depending on R such that for all i = 1, 2, ..., and h = ( h 1 h 2 ), Therefore, by these a priori bounds for u i , v i , (u i ) t , and (v i ) t , and by the monotone property of the sequences {u i } and {v i } as in (10), there exist u, v ∈ L 2 ((0, T ), Moreover, by Simon's result [13] we have the strong convergence of the sequences {u i } and {v i } in L 2 (Q) as below: Note that u, v satisfy the same initial conditions as the corresponding sequences. We now show that u, v is a weak solution of the system (1)- (3). For test functions The convergence of the terms which are nonlinear in the states uses the strong L 2 convergences, while the other terms only need weak convergences. Then, by passing through the limit in each terms, using the boundedness in (9), the convergences in (11) and (12), we see that Hence, we have proved that u, v is a weak solution of (1)-(3). Standard results using energy estimates give uniqueness of the solution of the state system.
Theorem 3.4 (Existence of an optimal control). There exists an optimal control Proof. First of all, note that because of the uniform L ∞ -bounds on u, v, and h 1 , h 2 , sup h∈U J( h) exists. Hence, we can choose a maximizing sequence, Let u n , v n be the corresponding states, i.e.
By the a priori estimates in Lemma 3.2, and by passing through a subsequence, we can assume that there exists u * , v * ∈ L 2 ((0, T ), Then, by Simon's result [13], and by passing through a subsequence, we can assume that we have the strong convergence of the sequences {u n } and {v n } in L 2 (Ω×[0, T ]) as below: Now, we show that ( h * 1 , h * 2 ) is an optimal control and that We now find the limits of each terms in this two equalities. Using the above weak convergences (13) and the L ∞ bounds on h 1 , h 2 and m we have Also, note that we can write Then, by the Cauchy's inequality, the strong convergence of { u n } , ∇φ 1 being in L 2 , and the L ∞ bounds on {h n 1 }, we see that the first term of (16) can be estimated as By the weak convergence of { h n 1 } , ∇φ 1 in L 2 (Q) and {u * } being in L ∞ -bounds, we see that the second term in (16) also satisfies Also similarly, we can obtain Hence, by collecting all of the previous estimates, we see that u * , v * are the state solutions of the system (1)-(3) corresponding to ( h * 1 , h * 2 ). On the other hand, using lower semi-continuity with respect to the weak convergences (13) on sequence {( h n 1 , h n 2 )} we complete our result: 4. The optimality system. To characterize the optimal control we differentiate the map ( h 1 , h 2 ) → J( h 1 , h 2 ) with respect to the controls h 1 and h 2 . But, as the states depend on the controls, u = u( h 1 , h 2 ), v = v( h 1 , h 2 ) and these states are in the objective functional J, we first need to differentiate the map ( the quotients converge weakly to ψ 1 , ψ 2 in L 2 ((0, T ), H 1 (Ω)) and the quotients of time derivatives converge weakly to (ψ 1 ) t , (ψ 2 ) t in L 2 ((0, T ), (H 1 (Ω)) * ). Further, the sensitivity functions ψ 1 and ψ 2 satisfy Proof. By considering the state PDEs corresponding to the optimal control ( h 1 , h 2 ), the control ( h 1 , h 2 ) and then forming the system solves by (u −u) , (v −v) , we obtain: with boundary and initial conditions For (19), using Cauchy's inequality and the Gronwall's inequality with standard estimation techniques, we can show that for some constants K 1 . Hence, by passing through a subsequence, we can assume that there exists ψ 1 (x, t) and On the other hand, by the Simon's results [13] and passing through a subsequence, as → 0 + we have the strong convergence of the quotients given below: Then, using the above convergences, as goes to zero in (19) we get the resulting PDEs (18) satisfied by ψ 1 and ψ 2 .

5.
Characterizing the optimal control. To characterize the optimal control, we develop the adjoint system with terms from the formal L 2 adjoint system from the terms in the sensitivity system. The source terms, A and B, driving the adjoint system are from differentiating the integrand of the objective functional with the respect to the states. Our adjoint system becomes: (20) Theorem 5.1 (Characterizing an optimal control h * ). Let ( h * 1 , h * 2 ) ∈ U be an optimal control for the problem (4) with the corresponding state solutions u * , v * ∈ L 2 ((0, T ), H 1 (Ω)). There exists adjoint functions λ 1 , λ 2 ∈ L 2 ((0, T ), H 1 (Ω)) with (λ 1 ) t , (λ 2 ) t ∈ L 2 ((0, T ), (H 1 (Ω)) * ) which are weak solutions of (20) corresponding to ( h * 1 , h * 2 ) and u * , v * . Then, the characterization of this optimal control is Proof. Note that the adjoint system is linear and a weak solution exists by standard results in [8]. For an optimal control ( h * 1 , h * 2 ) ∈ U we compute the directional derivative of the function J(( h 1 , h 2 ) * ) with respect to ( h * 1 , h * 2 ) in the directions ( l 1 , l 2 ) at u * and v * . Then, as J( h * 1 , h * 2 ) is the maximum value for the J( h 1 , h 2 ), we have Using the weak solution for the adjoint systems (20) to substitute for the first two terms in the right hand side of (21) Using the weak solution for the sensitivity system (18) Choosing appropriate variations l i , we can characterize ( h * 1 , h * 2 ) for the problem (4) as for i = 1, 2, ..., d.
The state equations (1)-(3), the adjoint equations (20) with the optimal control characterization in (22) form the optimality system. Under some conditions on the cost coefficients C and D, we can show this optimality system is unique for a sufficiently small time T . Proof. Suppose (u, v), (λ 1 , λ 2 ) and (ū,v), (λ 1 ,λ 2 ) are two solutions to the state pde system (1)- (3) and (20). Hence the corresponding optimal control characterizations are and where a ≥ 0 to be chosen. Then, substituting u = e at p, v = e at q andū = e atp ,v = e atq in (1) and then subtracting the corresponding two systems of equations give Next, substituting λ 1 = e −at ξ 1 , λ 2 = e −at ξ 2 andλ 1 = e −atξ 1 ,λ 2 = e −atξ 2 in (20) and subtracting the corresponding systems of equations give Note that by (23) and (24) we have Then, by using standard regularity theory for the equations of ξ k ,ξ k with k = 1, 2, and the boundedness of p,p, q,q, ξ 1 ,ξ 1 , ξ 2 ,ξ 2 , h 1 , h 2 , m, we can derive the L ∞bounds on |∇ξ k | and |∇ξ k |, where the bounds only depends on M 1 , M 2 , the initial data, and the coefficients of the system (1) (see for example, the proof [9, Lemma 5.2]). We illustrate the types of estimates needed by showing two terms below. Using the Cauchy's inequality, with C g showing L ∞ bounds on |∇ξ k |.
Then, we see the dependence in the estimate below on the size of C with similar estimate on h 2 controls with D 2 in the denominator of the bound.
We conclude the uniqueness of the solution of the optimality system for a small time interval.
6. Numerical results. To illustrate several scenarios, the optimality system, consisting of the state and adjoint system together with control characterization, was solved numerically. We considered the problem of (1) in one spatial dimension using the forward-backward sweep method [12,10]. An explicit finite difference scheme with an upwind method on advection terms was used to solve the state and adjoint systems in this iterative method. The specific optimality system to be solved corresponds to problem (4): with control characterization (22) For our numerical simulations, we used parameter values We considered several cases with different resource functions, initial conditions and competition coefficients. For resource functions, we chose m(x) = x 5 , sin(πx/5) and m = 6x 5 as shown in Figure 1. The initial conditions are shown in Figure 2. First, to study further the results in [9], we have considered one population with h 1 control under all three different resource functions in Figure 1 and the initial conditions in Figures 2a -2b. All the figures in [9] have small initial conditions, which allowed the advection control terms to be dominated by the spatial gradient of the resource function. As we will see, diffusion and an initial condition with a lot of variability may cause the populations to become more level and not to move always toward increasing resources.
The results we obtained for the u population only without competition are given in Figures 3 -5.
For the resource function given in Figure 1a and the initial conditions given in Figures 2a , 2b, we obtained the optimal control and the population distribution in Figure 3. From Figure 3a, it can be seen that when the initial population level is smaller the population tends to move towards the higher resources with a small velocity. But, compared to Figure 3a, in Figure 3c when the initial population level is higher they tend to use the directed movement with larger magnitude than in   Figures 2a , 2b, the optimal control and the population distribution are shown in Figure 4. When the resources are located at the center of the domain and the initial population is low, in Figure 4a, we see u tends to move towards the resources at the center from both sides with a smaller velocity. But, in contrast to Figure 4a, when the initial population is higher, Figure  4c shows that u tends to move towards the two boundaries even though the better resources are located at the center of the domain. Hence, it is clear that in spite of being at the location of better resources, when the population level is higher, the directed movement moves towards low density areas away from the high density area at the center.
We ran the simulations to see how a smaller population may act in a higher resource environment. In Figure 5 we have the optimal control h 1 and the population distribution of u when the resource function is higher as in Figure 1c with the    Figure 5. One population only: Population dynamics and optimal control with a smaller IC as in Figure 2a and larger resources m = 6x/5 as in Figure 1c ; (5a) Optimal control h 1 , (5b) Population distribution of u smaller initial condition in Figure 2a. In Figure 5 simulations with higher resources, the population grows faster in a shorter time. Hence, compared to Figure 3a, in Figure 5a, when the population level is smaller at the beginning, they tend to move towards the resources, and once the population reaches a higher level they move towards low density areas in spite of the location of the better resources. Hence, from the numerical simulations in Figure 3 -5, we can see that given the chance for directed movement, a population may not always move towards the resources but their movement depends on the population level also. Note that this feature was not illustrated in the numerical results of [9].
Next, we considered the numerical simulations for our competition PDE system with the two controls (28). From the numerical results in Figure 3 -5, we can see that when the population level is higher at some location, their directed movement may also be affected by diffusion. For this reason to show the competition effect separate from the diffusion effect, we have used smaller initial conditions as in Figure  2a , 2c in our numerical simulations for the competition system (28).
Using the initial conditions given in Figure 2a, with the resource functions given in Figure 1a, 1b, the numerical simulation results are in Figures 6 -7. Here, in both figures we have taken b 1 = b 2 = 4 and the both populations were taken to be with the same initial condition.  Figure 6. Two populations with competition: Population dynamics and optimal control with a smaller IC as in Figure 2a and the resource function m = x/5 as given in Figure 1a ; (6a) Optimal control h 1 , (6b) Population distribution of u Figure 6 gives the optimal control h 1 and the population distribution of u when the resources are at one boundary as in Figure 1a and the initial populations levels are smaller and located at the center of the spatial domain as in Figure 2a. As we have taken the same parameter values and the initial conditions for the two populations both u and v populations react in the same way and optimal control h 2 is as same as the optimal control h 1 in Figure 6a and the population distribution of v is as same as the population distribution of u in Figure 6b. In Figure 6a compared to the Figure 3a, we can see that when there is a competition in spite of the smaller population abundance level both populations tend to move away from each competing population's high density areas. Yet the individuals closer to the higher resources move a little faster toward resources than the individuals moving toward the other direction.
In Figure 7 we give the optimal control h 1 and the population distribution of u over Q when the resources are located at the center of the domain as in Figure  (  Further, to see the dynamics of two competing populations, we have run the numerical simulations with the initial condition in Figure 2c with the resource functions in Figure 1a, 1b. In Figure 8 we have the two optimal controls and the population distributions of the two populations over Q when the resources are higher at one boundary as in Figure 1a and the initial populations levels are smaller and are overlapping only at the center of the spatial domain as in Figure 2c. From the results we can see the two populations tend to move towards the two opposite boundaries to move away from each other. In Figure 9 we give the two optimal controls and the population distributions of the two populations over Q when the resources are higher at the center as in Figure 1b and the initial populations levels are smaller and are overlapping only at the center of the spatial domain as in Figure 2c. Again, similar to Figure 8, we can see the two populations tend to move towards the two opposite boundaries to move away from each other and not towards the better resources. Hence, from the numerical results in Figures 6 -9 we can see that when there is a high competition (b 1 = b 2 = 4), in spite of the smaller population abundance level and the favorable resource location, they tend to move away from each other. But, this response could depend on the competition coefficients b 1 and b 2 and we need further analysis to see any relationship.
To see how each population will react when they have different competition coefficients, in Figure 10, we have taken the competition coefficients to be b 1 = 4 (a) (b) (c) (d) Figure 10. Two populations with different competition rates: Population dynamics and optimal control with a smaller IC as in Figure 2a and the resource function m = sin(πx/5) as given in Figure 1b and b 1 = 4 , b 2 = 0.5 ; (10a) Optimal control h 1 , (10b) Population distribution of u, (10c) Optimal control h 2 , (10d) Population distribution of v and b 2 = 0.5 and considered the population dynamics and the optimal controls with a smaller initial condition in Figure 2a and resources in Figure 1b. Note that here again we take initial populations to be the same for both u and v. Compared to the results in Figure 7, in Figure 10 we see that when the competition coefficients are different two populations move in different ways. From Figure 10a we can see that u population tend to move towards the two boundaries to move away from v population. But, in contrast to Figure 10a, in Figure 10c we can see that v population tends to move towards the center. Hence, from the numerical results in Figure 10 we can see that the directed movements of the competing populations also depend on the competition rates.

7.
Conclusions. Optimal control analysis for controlling the advection directions in a parabolic PDE systems for competing populations was completed to obtain existence, uniqueness and characterization of the optimal control. The characterization involving the adjoint functions shows the explicit dependence on the coefficients, A, B, C and D, in the objective functional.
With numerical simulations for one population only, we were able to show the population does not always choose the advection direction to move toward increasing resources. When the initial condition has a sufficiently high population with some variation, the movement may be chosen to move to level the population, sometimes instead of moving toward increasing resources. In the systems case, the level of the competition coefficients can also influence the choice of movement direction. This numerical work indicates some interesting relationships between the optimal advective directions and the initial conditions, diffusion effects, competition rates, and resources; further work to study these relationship is warranted.
The features of these PDEs and the corresponding numerical simulations need to be generalized for further investigation. Having variable diffusion, growth and competition rates will also have an impact on the optimal controls. The parameters A, B in the objective functional could be functions to force different choices on the populations. It would also be interesting to do some simulations in two space dimensions. Population interactions besides competition may also be considered. Controlling advection in a population coming from flows in a river would be an interesting extension. In our future work, we expect to continue this investigation with more generality.