Optimal control strategies for the sterile mosquitoes technique
Luís Almeida, Michel Duprez, Yannick Privat, Nicolas Vauchelet
aa r X i v : . [ m a t h . A P ] N ov Optimal control strategies for the sterile mosquitoes technique
Luis Almeida ∗ Michel Duprez † Yannick Privat ‡ Nicolas Vauchelet § November 10, 2020
Abstract
Mosquitoes are responsible for the transmission of many diseases such as dengue fever, zika orchigungunya. One way to control the spread of these diseases is to use the sterile insect technique(SIT), which consists in a massive release of sterilized male mosquitoes. This strategy aims atreducing the total population over time, and has the advantage being specific to the targetedspecies, unlike the use of pesticides. In this article, we study the optimal release strategies inorder to maximize the efficiency of this technique. We consider simplified models that describe thedynamics of eggs, males, females and sterile males in order to optimize the release protocol. Wedetermine in a precise way optimal strategies, which allows us to tackle numerically the underlyingoptimization problem in a very simple way. We also present some numerical simulations toillustrate our results.
Keywords:
Sterile insect technique, population dynamics, optimal control problem, PontryaginMaximum Principle (PMP).
The sterile insect technique (SIT) consists in massively releasing sterilized males in the area whereone wishes to reduce the development of certain insects (mosquitoes in this case). Since the releasedsterile males mate with females, the number of offspring is then reduced, and the size of the insectpopulation diminishes. This strategy has been first studied by R. Bushland and E. Knipling andexperimented successfully in the early 1950’s by nearly eradicating screw-worm fly in North America.Since then, this technique has been considered for different pests and disease vectors [5, 14].Among such vectors, mosquitoes, especially Aedes mosquitoes, are responsible for the transmis-sion to humans of many diseases for which there is currently no efficient vaccine. Thus, the sterileinsect technique and the closely related incompatible insect technique are very promising tools tocontrol the spread of such diseases by reducing the size of the vector population (there are caseswhere these techniques have been successfully used to drastically reduce mosquito populations insome isolated regions, e.g. [21, 24]). ∗ Sorbonne Universit´e, CNRS, Universit´e de Paris, Inria, Laboratoire Jacques-Louis Lions UMR7598, F-75005 Paris,France ( [email protected] ) † Inria, ´equipe MIMESIS, Universit´e de Strasbourg, Icube, CNRS UMR 7357, Strasbourg, France( [email protected] ). ‡ IRMA, Universit´e de Strasbourg, CNRS UMR 7501, Inria, 7 rue Ren´e Descartes, 67084 Strasbourg, France( [email protected] ). § Laboratoire Analyse, G´eom´etrie et Applications CNRS UMR 7539, Universit´e Sorbonne Paris Nord, Villetaneuse,France ( [email protected] ).
1n order to study the efficiency of this technique and to optimize it, mathematical modeling isof great use. For instance, in [3, 12, 13], the authors propose mathematical models to study thedynamics of the mosquito population when releasing sterile males. Recently in [22], the authorspropose and analyze a differential system modeling the mosquito population dynamics. Their modelis based on experimental observations and is constructed by assuming that there is a strong Alleeeffect in the insect population dynamics. A similar model, without strong Allee effect, is investigatedin [2]. Control theory allows also to study the feasibility of controlling the population thanks to thesterile insect technique, and it has been studied in several works, see e.g. [7, 8, 4]. Using suchmathematical models, authors are able to compare the impact of different strategies in releasingsterile mosquitoes (see e.g. [10, 20] and [13, 19] where periodic impulsive releases are considered).In order to find the best possible release protocol, optimal control theory may be used. In [15],optimal control methods are applied to the rate of introduction of sterile mosquitoes. An approachdevelopped in [23] attempts to control both breeding rates and the rate of introduction of sterilemosquitoes. In [16], the influence of habitat modification is also considered. Finally, existence andnumerical simulations of the solution to an optimal control problem for the SIT has been proposedin [9].In this paper, we study how to optimize the release protocol in order to minimize certain costfunctionals, such as the number of mosquitoes. Starting with the mathematical model presentedin [22] without Allee effect, we investigate some optimal control problems and focus on obtaininga precise description of the optimal control. To do so, we consider a simplified version of themathematical model and we perform a complete study of the optimizers. In particular, in our mainresult, we describe precisely the optimal release function to minimize the number of sterile malesneeded to reach a given size of the population of mosquitoes. Our theoretical results are illustratedwith some numerical simulations. We also provide some extensions to certain related optimizationproblems in order to illustrate the robustness of our approach.The outline of this paper is as follow. In Section 2, we introduce the mathematical model wewill adopt for the sterile insect technique, and describe some useful qualitative properties related tostability issues. For the sake of readability, all the proofs will be postponed to Appendix A. Section 3is devoted to the introduction of the problems modeling the search of optimal release protocols andthe statement of the main theoretical results providing a precise description of optimal strategies. Wethen derive a simple algorithm to compute them numerically, and provide illustrating simulations.The proofs of the main theoretical results are postponed to Section 4. Finally, some comments onother possible approaches are gathered in Section 5.
The life cycle of a mosquito (male or female) consists of several stages and takes place successivelyin two distinct environments: it includes an aquatic phase (egg, larva, pupa) and an aerial phase(adult). A few days after mating, a female mosquito may lay a few dozen eggs, possibly spreadover several breeding sites. Once laid, the eggs of some species can withstand hostile environments(including adverse weather conditions) for up to several months before hatching. This characteristiccontributes to the adaptability of mosquitoes and has enabled them to colonize temperate regions.After stimulation (e.g. rainfall), the eggs hatch to give birth to larvae that develop in the water andreach the pupa state. This larval phase can last from a few days to a few weeks. Then, the insectundergoes its metamorphosis. The pupa (also called nymph ) remains in the aquatic state for 1 to 3days and then becomes an adult mosquito (or imago ): it is the emergence and the beginning of the2erial phase. The lifespan of an adult mosquito is estimated to be of a few weeks.In many species, egg laying is only possible after a blood meal, i.e. the female must bite avertebrate before each egg laying. This behavior, called hematophagy, can be exploited by infectiousagents (bacteria, viruses or parasites) to spread, alternately from a vertebrate host (humans, forwhat we are interested in here) to an arthropod host (here, the mosquito).Based on these observations, a compartmental model has been introduced in [22] to model thelife cycle of mosquitoes when releasing sterile mosquitoes. In what follows, we will both deal withthe full and a simplified version of [22]. The reason for studying such a simplified model is twofold:on the one hand, the simplified model can be considered relevant from a biological point of viewwithin certain limits. On the other hand, the study of such a “prototype” model can be consideredas a first step towards the development of robust control methodologies with a wider application.To this aim, we will denote by u ( · ) a control function standing for a sterile male release function(in other words the rate of sterile male mosquitoes release at each time) and by • M s ( t ), the sterilized adult males at time t ; • F ( t ), the adult females that have been fertilized at time t .The system we will use for describing the behavior of the mosquito population under the action ofthe control u ( · ) reads dFdt = f ( F, M s ) ,dM s dt = u − δ s M s , ( S )where f : R → R denotes the nonlinear function f ( F, M s ) = ν (1 − ν ) β E ν E F (cid:0) β E FK + ν E + δ E (cid:1)(cid:0) (1 − ν ) ν E β E F + δ M γ s M s ( β E FK + ν E + δ E ) (cid:1) − δ F F, (1)with the following parameter choices: • β E > • δ E , δ M , δ F , δ s > • ν E > • ν ∈ (0 ,
1) the probability that a pupa gives rise to a female, and (1 − ν ) is therefore theprobability to give rise to a male; • K > • γ s > MM + γ s M s .In the following section, we explain and comment on the choice of this system.3 .2 Derivation of the simplified model and presentation of the original one The choice of ( S ) as model is inspired by [22]. To explain how it has been derived, let us presentthe more involved model we have considered. Let us introduce: • E ( t ), the mosquito density in aquatic phase at time t ; • M ( t ), the adult male density at time t ; • M s ( t ), the sterilized adult male density at time t ; • F ( t ), the density of adult females that has been fertilized at time t .Then, the dynamics of the mosquito population is driven by the following dynamical system: dEdt = β E F (cid:18) − EK (cid:19) − (cid:0) ν E + δ E (cid:1) E,dMdt = (1 − ν ) ν E E − δ M M,dFdt = νν E E MM + γ s M s − δ F F,dM s dt = u − δ s M s . ( S )Regarding this latter model, the main difference with the one in [22] is the absence of an exponentialterm in the equation on F to introduce an Allee effect. This effect reflects the fact that, when thepopulation density is very low, it can be difficult to find a partner to mate. This term is importantwhen considering a small population size. Here, since we are focusing on large populations that wewant to reduce in size, we will neglect this term.Assuming that the time dynamics of the mosquitoes in aquatic phase and the adult males com-partments are fast leads to assume that the equations on E ( · ) and M ( · ) are at equilibrium. We referfor instance to [1] for additional explanations on the justification for these asymptotics. Hence, weget the following equalities E = β E F β E FK + ν E + δ E and M = (1 − ν ) ν E δ M E .
Plugging such expressions into ( S ) allows us to obtain ( S ).We conclude this paragraph by numerically comparing the full model ( S ) and the simplified one( S ) that we will aim to control. We consider the numerical values taken from [22, Table 3] andrecalled in Table 1 below. 4
20 40 6000 . . . . . · t F for system ( S ) F for system ( S ) 0 20 40 600 . . . . · t F for system ( S ) F for system ( S ) Figure 1: Comparisons of the numerical solutions F (in continous blue line for System ( S ) and indashed red line for System ( S ) or equation 4). The initial conditions correspond to the “persistence”equilibrium (see propositions 2.1 and 2.2) Left: u ( · ) = 15 000. Right : the releases occur every 10days with an intensity of 20 000 mosquitoes in one day, i.e. u ( · ) = 20 000 P k =0 [10 k, k +1] . Parameter Name Value interval Chosen value β E Effective fecundity 7.46–14.85 10 γ s Mating competitivenessof sterilizing males 0–1 1 ν E Hatching parameter 0.005–0.25 δ E Mosquitoes in aquatic phasedeath rate 0.023 - 0.046 0.03 δ F Female death rate 0.033 - 0.046 0.04 δ M Males death rate 0.077 - 0.139 0.1 δ s Infected male death rate 0.12 ν Probability of emergence 0.49Table 1: Value intervals of the parameters for systems ( S ) and ( S ) (see [22])The numerical results are shown in Fig. 1. In these simulations, the time of the experiment isassumed to be T = 70 days, and we choose two different release functions u : u ( · ) = 15 000 (left), and u ( · ) = 20 000 X k =0 [10 k, k +1] (right) . The dynamics for F (female compartment) is represented for both systems ( S ) (blue, continuousline) and ( S ) (red, dashed line). We observe that both problems are very close, which indicatesthat the dynamics of fertilized females in system ( S ) may be approximated by the one in system( S ).A mathematical element of this observation lies in the fact that the equilibria Systems ( S ) and( S ) coincide. When dealing with optimal control properties, we will also numerically observe inSection 3.3 that this simplification does not affect the optimal strategies in a strong way.5 .3 Mathematical properties of the dynamical systems This section is devoted to establishing stability properties for equilibria of Systems ( S ) and ( S ) inthe absence of control, in order to qualitatively understand their behavior whenever initial data arechosen close to equilibria. These results can be considered as preliminary tools before investigatingoptimal control properties for these problems.Recall that both systems share the same steady states. We will moreover show that they enjoythe same stability properties. For the sake of readability, all proofs are postponed to Appendix A.In what follows, we will make the following assumption, in accordance with the numerical valuesgathered in Table 1: δ s > δ M and νβ E ν E > δ F ( ν E + δ E ) ( H ) Proposition 2.1 (Stability properties for System ( S )) . Let us assume that ( H ) holds.(i) If u ( · ) = 0 , then, System ( S ) has two equilibria: • the “extinction” equilibrium ( E ∗ , M ∗ , F ∗ , M ∗ s ) = (0 , , , , which is linearly unstable; • the “persistence” equilibrium ( E ∗ , M ∗ , F ∗ , M ∗ s ) = (cid:0) E, M , F , (cid:1) , where E = K − δ F (cid:0) ν E + δ E (cid:1) β E νν E ! , M = (1 − ν ) ν E δ M E, F = νν E δ F E, (2) which is linearly asymptotically stable.(ii) If the control function u is assumed to be non-negative, then the corresponding solution ( E, M,F, M s ) to System ( S ) enjoys the following stability property: E (0) ∈ (0 , E ] M (0) ∈ (0 , M ] F (0) ∈ (0 , F ] M s (0) > ⇒ E ( t ) ∈ (0 , E ] M ( t ) ∈ (0 , M ] F ( t ) ∈ (0 , F ] M s ( t ) > for all t > . Finally, let U ∗ be defined by U ∗ := Kβ E ν (1 − ν ) ν E δ s γ s ( ν E + δ E ) δ F δ M (cid:18) − δ F ( ν E + δ E ) β E νν E (cid:19) (3) and let U denote any positive number such that U > U ∗ . If u ( · ) denotes the constant con-trol function almost everywhere equal to U for all t > , then the corresponding solution ( E ( t ) , M ( t ) , F ( t )) to System ( S ) converges to the extinction equilibrium as t → + ∞ . In the following result, we will again use the notations introduced in Prop. 2.1 above.
Proposition 2.2 (Stability properties for System ( S )) . Let us assume that ( H ) holds.(i) If u ( · ) = 0 , System ( S ) has two equilibria: • the “extinction” equilibrium ( F ∗ , M ∗ s ) = (0 , , which is unstable if δ s > δ F . • the “persistence” equilibrium ( F ∗ , M ∗ s ) = (cid:0) F , (cid:1) , is linearly asymptotically stable. ii) If the control function u is assumed to be non-negative, then the corresponding solution ( F, M s ) to System ( S ) enjoys the following stability property: (cid:26) F (0) ∈ (0 , F ] M s (0) > ⇒ (cid:26) F ( t ) ∈ (0 , F ] M s ( t ) > for all t > . If u ( · ) denotes the constant control function almost everywhere equal to U > U ∗ (defined by (3) ), then F ( t ) goes to as t → + ∞ . As a consequence of Propositions 2.1 and 2.2, it follows that if the horizon of time of control T is large enough, by releasing a sufficient amount of mosquitoes, it is possible to make the wildpopulation of mosquitoes be as small as desired at time T . In the next section, we will investigatethe issue of minimizing the number of released mosquitoes in order to reach a given size of the wildpopulation. In Section 5, we will also comment on other relevant minimization problems related tothe sterile insect technique. In what follows, we will consider an initial state corresponding to the beginning of the experiment,in which there are no sterile mosquitoes. Hence, the mosquito population is at the “persistence”equilibrium at the beginning of the experiment, meaning that one chooses F (0) = F , M s (0) = 0 , for System ( S ), to which we should add E (0) = E, M (0) =
M , for System ( S ) (4)where ( E, M , F ) are defined by (2). Our aim is to determine the optimal time distribution of thereleases, such that the number of mosquitoes used to reach a desired size of the population is assmall as possible.
Let us first model the optimization of the release procedure. Given the duration of the experiment,we aim at minimizing the amount of mosquitoes required to reduce the size of the wild populationof mosquitoes to a given value:
What should be the time distribution of optimal releases in order to reach a given value of the wildpopulation at the end of the experiment, by using as few sterilized mosquitoes as possible?
To answer this question, we first define a cost functional that will stand for the objective we aretrying to achieve. We mention that other formulations of minimization problem are also proposedin Section 5.Let
T >
U > ε > S ) and ( S ). We model the optimal control problems for the sterile mosquito releasestrategy as inf u ∈U ( S T,U,ε J ( u ) , ( P ( S ) T,U,ε )and inf u ∈U ( S T,U,ε J ( u ) , ( P ( S ) T,U,ε )7here the functional J stands for the total number of released mosquitoes during the time T , namely J ( u ) := Z T u ( t ) dt. For a given ε >
0, we introduce the sets of admissible controls for Systems ( S ) and ( S ), respectivelydenoted U ( S ) T,U,ε and U ( S ) T,U,ε and defined by U ( S ) T,U,ε := n u ∈ L ∞ (0 , T ) : 0 u U a.e. , F ( T ) ε with F solution of System ( S ) o (5)and U ( S ) T,U,ε := n u ∈ L ∞ (0 , T ) : 0 u U a.e. , F ( T ) ε with F solution of System ( S ) o . Remark 3.1 (Exact null-controllability) . It is notable that there does not exist any control u ∈ L ∞ (0 , T ; R + ) such that the corresponding solution F to System ( S ) (resp. System ( S ) ) satisfies F ( T ) = 0 since one has F ( t ) > F (0) e − δ F t for every t ∈ [0 , T ] . Let us now state the main results of this article. In what follows, we will use the notation U ∗ todenote the positive number defined by (3), and the notations J ( S i ) T,U,ε , i = 1 , P ( S ) T,U,ε ) and ( P ( S ) T,U,ε ), namely J ( S i ) T,U,ε := inf u ∈U ( S i ) T,U,ε J ( u ) , i = 1 , . (6) Theorem 3.2.
Let us assume that Condition ( H ) holds true. Let ε ∈ (0 , F ) and U > U ∗ . Thereexists a minimal time T > such that for all T > T , the set U ( S ) T,U,ε is nonempty and Problem ( P ( S ) T,U,ε ) has a solution u ∗ . One has J ( S ) T,U,ε U T and the mappings T ∈ [ T , + ∞ ) J ( S ) T,U,ε and U ∈ ( U ∗ , + ∞ ) J ( S ) T,U,ε are non-increasing. Furthermore, there exists t ∈ (0 , T ) such that u ∗ = 0 on ( t , T ) , one has F ( T ) = ε and F ∈ [ ε, F ] on (0 , T ) . It is interesting to note that there is no need to release sterile mosquito at the end of theexperiment. For the simplified system ( S ), we can obtain a much more precise characterization ofoptimal controls, which is the purpose of the next result. As a preliminary remark, recall that thefunction f is defined by (1). Theorem 3.3.
Let us assume that Condition ( H ) holds true and that δ s > δ F . Let ε ∈ (0 , F ) , U > U ∗ . There exists a minimal time T > such that for all T > T , the set U ( S ) T,U,ε is non emptyand Problem ( P ( S ) T,U,ε ) has a solution u ∗ , characterized as follows:(i) Optimal value: one has J ( S ) T,U,ε U T and the mappings T ∈ [ T , + ∞ ) J ( S ) T,U,ε and U ∈ ( U ∗ , + ∞ ) J ( S ) T,U,ε are non-increasing. ii) Optimal control: let
T > T . If
U > U ∗ is large enough, there exists ( t , t ) ∈ [0 , T ] with t t such that(a) u ∗ = 0 on (0 , t ) or u ∗ = U on (0 , t ) ;(b) u ∗ ∈ (0 , U ) on ( t , t ) with u ∗ ( t ) = (cid:18) ∂ f∂M s (cid:19) − (cid:18) ∂f∂M s ∂f∂F + δ s M s ∂ f∂M s − f ∂ f∂M s ∂F (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ( F ( t ) ,M s ( t )) ; (7) (c) u ∗ = 0 on ( t , T ) .Furthermore, there exists T ∗ ( U ) (simply denoted T ∗ when no confusion is possible) such that T ∗ > T and if T > T ∗ , then the optimizer u ∗ is unique, with < t < t < T , and such that u ∗ = 0 on (0 , t ) . Finally, the mapping T ∈ ( T ∗ , + ∞ ) J ( S ) T,U,ε is constant.(iii)
Optimal trajectory: we extend the domain of definition of F to R + by setting u ( · ) = 0 on ( T, + ∞ ) . One has F ( T ) = ε , F ′ on (0 , T ) and F ′ > on ( T, + ∞ ) . In particular, F hasa unique local minimum at T which is moreover global, equal to ε and F ′ ( T ) = 0 . Let us comment on this result. Under the assumptions of Theorem 3.3 (in particular that T and U are large enough), the infinite dimensional control problem ( P ( S ) T,U,ε ) can be reduced to a twodimensional one: more precisely, one only needs to determine the two parameters t and t .As expected the horizon of time T fixed for the control to reach a desired number of adult femalemosquitoes has a strong influence on the form of the optimal control. The longer T is, the smalleris the number of sterilized males needed is. However, there is a maximal time T ∗ above which thenumber of released mosquitoes needed to reach the desired state is stationary with respect to thehorizon of time. Remark 3.4.
Stabilization and feedback Note that our results can be interpreted in terms of stabi-lization: indeed, it follows from Theorem 3.3, and in particular from the description of the optimaltrajectory that, under the assumptions of this theorem, the differential system ddt (cid:18) FM s (cid:19) = (cid:18) f ( F, M s ) u − δ s M s (cid:19) in [0 , + ∞ ) with u = ∂f∂Ms ( F,M s ) ∂f∂F ( F,M s )+ ∂ f∂M s ( F,M s ) δ s M s − ∂ f∂Ms∂F ( F,M s ) f ( F,M s ) ∂ f∂M s ( F,M s ) , complemented with the initial conditions F (0) = F and M s (0) = 0 , satisfies lim t → + ∞ F ( t ) = 0 , in other words, the control function u ( · ) above defines a feedback control stabilizing System ( S ) . In this section, we introduce an elementary algorithm for computing the (unique) solution to Prob-lem ( P ( S ) T,U,ε ). The chosen approach mainly rests upon the precise description of the optimizer (when-ever T and U are large enough) provided in Theorem 3.3.9et us describe our approach. With the notations of Theorem 3.3, we assume that T > T ∗ sothat the optimal control is unique. We will take advantage of its particular form, more preciselythat there exist 0 < t < t < T such that u ∗ = 0 on (0 , t ), u ∗ is given by (7) on ( t , t ), u ∗ = 0 on( t , T ) and F ( T ) = ε .To this aim, let us introduce for τ ∈ [0 , T ] the auxiliary Cauchy system ddt (cid:18) FM s (cid:19) = (cid:18) f ( F, M s ) u τ − δ s M s (cid:19) in [0 , + ∞ )with u τ = ∂f∂Ms ( F,M s ) ∂f∂F ( F,M s )+ ∂ f∂M s ( F,M s ) δ s M s − ∂ f∂Ms∂F ( F,M s ) f ( F,M s ) ∂ f∂M s ( F,M s ) (0 ,τ ) , (8)complemented with the initial conditions F (0) = F and M s (0) = 0, whose solution, associated tothe control function u τ , will from now on be denoted ( F τ , M τ s ). Let T > T ∗ .The uniqueness property for Problem ( P ( S ) T,U,ε ) allows us to get the following relations betweenthe solution to Problem ( P ( S ) T,U,ε ) and the control u τ . Property 3.5.
Under the assumptions and with the notations of Theorem 3.3, let
U > U ∗ (inparticular we consider T > T ∗ ( U ) so that the conclusion (ii) holds true). There exists ε > suchthat, if ε ∈ (0 , ε ) , then(i) for every τ ∈ [0 , T ) , there exists a unique τ ( τ ) ∈ ( τ , + ∞ ) such that F τ is first strictlydecreasing on (0 , τ ( τ )) , and then strictly increasing on ( τ ( τ ) , + ∞ ) .(ii) the value function ψ : τ ∈ [0 , T ) min t ∈ [0 , ∞ ] F τ ( t ) = F τ ( τ ( τ )) is decreasing.(iii) the optimal control u ∗ solving Problem ( P ( S ) T,U,ε ) satisfies u ∗ ( t ) = (cid:26) if t ∈ (0 , T − τ ( τ )) ,u τ ( t − T + τ ( τ )) otherwise (9) a.e. on (0 , T ) , where τ denotes the unique solution on [0 , T ) to the equation ψ ( τ ) = ε . We now construct an efficient algorithm based on this result to solve Problem ( P ( S ) T,U,ε ), involvinga bisection type method.
In order to illustrate our main results (Theorems 3.2 and 3.3), we provide some numerical simulationsof the optimal control problems ( P ( S ) T,U,ε ) and ( P ( S ) T,U,ε ). We use the parameters values provided inTable 1, that come from [22, Table 1-3].As in [22], in order to get results relevant for an island of 74 ha (hectares) with an estimatedmale population of about 69 ha − , the total number of males at the beginning of the experimentis taken equal to M ∗ = 69 ×
74 = 5106. Regarding System ( S ), we assume that M s = 0 at theequilibrium, and we then deduce that E = δ M (1 − ν ) ν E M , F = νν E δ F E = νδ M δ F (1 − ν ) M lgorithm 1: Solving Problem ( P ( S ) T,U,ε ) Initialization:
Let n ∈ N ∗ , τ , min = 0 and τ , max = T } ) While i ≤ n do τ , test = ( τ , min + τ , max ) / Solve (8) on (0 , T ) for τ = τ , test and let u τ be the function given by (8). if min (0 ,T ) F τ < ε then τ , max = τ , test else τ , min = τ , test End:
Let τ n = ( τ , min + τ , max ) / t n = argmin (0 ,T ) F τ n and u n ( t ) = (cid:26) t ∈ (0 , T − t n ) ,u ( t − T + t n ) otherwise , for each t ∈ (0 , T ).and we thus evaluate the environmental capacity for eggs as K = − δ F (cid:0) ν E + δ E (cid:1) β E νν E ! − E. Regarding System ( S ), we consider the same initial quantity of females F and get that theenvironmental capacity K has the same expression as above.To compute the solutions to Problem ( P ( S ) T,U,ε ) and ( P ( S ) T,U,ε ), we use the opensource optimizationroutine
GEKKO (see [6]) which solves the optimization problem thanks to the
APOPT (AdvancedProcess OPTimizer) library, a software package for solving large-scale optimization problems (see[17]). We will also compare the numerical solution obtained when solving Problem ( P ( S ) T,U,ε ) with theone obtained by Formula (9), by applying Algorithm 1.Figures 2, 3 and 4 gather the solutions of the optimal control problems ( P ( S ) T,U,ε ), ( P ( S ) T,U,ε ) andthe function given by (9) for T = 60 , , U = 15000, ν E = 0 .
05 and ε = F / , T ) is discretized with 300 points. We do not give the solution to (9) for T = 60 in Figure 2, since Algorithm 1 cannot be applied in this situation (indeed, T is not largeenough and the involved has not a convenient form to apply on it a bisection procedure).These simulations allow us to recover the theoretical results of Section 3. Moreover, it is interest-ing to observe that, even if we do not have a proof of this fact, the optimal controls for the differentproblems look very closed. This is confirmed by Table 2, where the quantity of sterile insects iscomputed for different values of ν E .Regarding Problem ( P ( S ) T,U,ε ), since for T large enough, the optimal control vanishes on an intervalat the beginning of the experiment, the system stays at the equilibrium on this interval. Hence, with11he notations of Theorem 3.3, the optimal time T opt to control the system with a singular part isequal to T − t . As illustrated in Table 3, this optimal time T opt is not very sensitive to the choiceof optimal control problem, namely ( P ( S ) T,U,ε ), ( P ( S ) T,U,ε ) and (9) and also of the value of ν E .We finally provide in Table 4 for different discretizations of the time interval the computationtime for solving Problems ( P ( S ) T,U,ε ) and ( P ( S ) T,U,ε ) with the opensource optimization routine
GEKKO (default parameters) and of (9) thanks to Algorithm 1 ( n = 50 iterations). We have here usedan Intel CORE i5 8th Gen.
As expected, one can notice that Algorithm 1 yields a much fasterresolution than the other approaches. . . . . . · t Optimizer F to ( P ( S ) T,U,ε )Optimizer F to ( P ( S ) T,U,ε ) ε , ,
000 t
Optimizer u to ( P ( S T,U,ε )Optimizer u to ( P ( S T,U,ε ) U Figure 2: Solution of the optimal control problems ( P ( S ) T,U,ε ) and ( P ( S ) T,U,ε ) with T = 60, U = 5000, ν E = 0 .
05, and ε = F / . . . . . · t Optimizer F to ( P ( S T,U,ε )Optimizer F to ( P ( S T,U,ε )Optimizer F to (9) ε , ,
000 t
Optimizer u to ( P ( S T,U,ε )Optimizer u to ( P ( S T,U,ε )Optimizer u to (9) U Figure 3: Solution of the optimal control problems ( P ( S ) T,U,ε ), ( P ( S ) T,U,ε ) and (9) with T = 150, U = 5000, ν E = 0 .
05, and ε = F /
4. 12
50 100 150 20000 . . . . . · t Optimizer F to ( P ( S T,U,ε )Optimizer F to ( P ( S T,U,ε )Optimizer F to (9) ε , ,
000 t
Optimizer u to ( P ( S T,U,ε )Optimizer u to ( P ( S T,U,ε )Optimizer u to (9) U Figure 4: Solution of the optimal control problems ( P ( S ) T,U,ε ), ( P ( S ) T,U,ε ) and (9) with T = 200, U = 5000, ν E = 0 .
05, and ε = F / ν E J ( u ∗ ) for ( P ( S ) T,U,ε ) 1.21e5 1.34e5 1.42e5 1.45e5 J ( u ∗ ) for ( P ( S ) T,U,ε ) 1.15e5 1.30e5 1.40e5 1.43e5 J ( u ∗ ) for (9) 1.15e5 1.30e5 1.40e5 1.43e5 ν E J ( u ∗ ) for ( P ( S ) T,U,ε ) 1.47e5 1.49e5 1.50e5 1.50e5 J ( u ∗ ) for ( P ( S ) T,U,ε ) 1.46e5 1.49e5 1.49e5 1.50e5 J ( u ∗ ) for (9) 1.46e5 1.48e5 1.49e5 1.50e5Table 2: J ( u ∗ ) for ( P ( S ) T,U,ε ), ( P ( S ) T,U,ε ) and (9) with respect to ν E for ε = F / F = 11037. ν E T opt for ( P ( S ) T,U,ε ) 109 107 105 105 104 104 104 104 T opt for ( P ( S ) T,U,ε ) 106 105 104 104 104 104 104 104 T opt for (9) 106 105 104 103 103 103 103 103Table 3: Optimal time T opt of control for ( P ( S ) T,U,ε ), ( P ( S ) T,U,ε ) and (9) with respect to ν E for ε = F / F = 11037. 13 ime discr.
50 100 200 400 800 1600 3200C. t. for ( P ( S ) T,U,ε ) 2.64 8.11 2.56e1 1.86e2 3.83e4 X XC. t. for ( P ( S ) T,U,ε ) 2.05 6.00 2.94e1 1.13e2 7.37e2 3.53e4 XC. t. for (9) 3.58e-1 6.74e-1 1.44 2.73 5.45 11.7 2.24e1Table 4: Time of computation (sec.) for ( P ( S ) T,U,ε ), ( P ( S ) T,U,ε ) and (9) with respect to the time dis-cretization for ν E = 0 . ε = F / F = 11037. In the whole proof, we will denote by ( f E ( E, F ), f M ( E, M ), f F ( E, M, F, M s )) ⊤ the right-hand sideof the differential subsystem of System ( S ) satisfied by ( E, M, F ) ⊤ , namely f E ( E, F ) = β E F (cid:18) − EK (cid:19) − (cid:0) ν E + δ E (cid:1) E,f M ( E, M ) = (1 − ν ) ν E E − δ M M and f F ( E, M, F, M s ) = νν E E MM + γ s M s − δ F F. We first point out that system ( S ) enjoys a monotonicity property with respect to the control u . Lemma 4.1.
Let u , u ∈ L ∞ (0 , T ; R + ) be such that u > u . Let us assume that (4) holds.Then, the associated solutions ( E , M , F , M s ) , ( E , M , F , M s ) to System ( S ) respectivelyassociated to u and u satisfy ( E , M , F ) ( E , M , F ) , the inequality being understood compo-nent by component.Proof. According to Proposition 2.1, we have (
E, M, F ) ∈ [0 , E ] × [0 , M ] × [0 , F ]. Noting that M s i ( t ) = Z t u i ( s ) e δ s ( s − t ) ds, i = 1 , , it follows that M s > M s . Hence, the monotonicity property follows since one has ∂f E ∂F > , ∂f M ∂E > , ∂f F ∂E > , ∂f F ∂M > , ∂f F ∂M s , and therefore the so-called Kamke-M¨uller conditions (see e.g. [18]) hold true.Let us investigate the existence of an optimal control for Problem ( P ( S ) T,U,ε ). Lemma 4.2.
Let ε ∈ (0 , F ) and U > U ∗ . There exists T ( U ) > such that for all T > T ( U ) , theset U ( S ) T,U,ε is nonempty and for such a choice of T , Problem ( P ( S ) T,U,ε ) has a solution u ∗ .Moreover, one has J ( S ) T,U,ε U T (where J ( S ) T,U,ε is defined by (6) ) and J ( S ) T,U,ε is non-increasingwith respect to T > T and U > U ∗ . roof. According to Proposition 2.1, if u ( · ) = U , then F ( t ) → t goes to + ∞ . Hence, for any ε ∈ (0 , F ) and since F is Lipschitz-continuous, there exists T > F ( T ) ε meaning that T > T , u = U [ T − T ,T ] ∈ U ( S ) T,U,ε . The set U ( S ) T,U,ε is thus nonempty.Let us now investigate the existence property. For the sake of clarity, we temporarily denoteby ( E u , M u , F u , M us ) the solution of System ( S ) associated to u . Let ( u n ) n ∈ N be a minimizingsequence. According to the Banach-Alaoglu Bourbaki theorem, the set U ( S ) T,U,ε is compact for theweak- ∗ topology of L ∞ (0 , T ) and, up to a subsequence, ( u n ) n ∈ N converges to u ∈ L ∞ (0 , T, [0 , U ]).The functional J is obviously continuous for the weak- ∗ topology of L ∞ (0 , T ), so that it only remainsto prove that u ∈ U ( S ) T,U,ε , in other words that F u ( T ) ε . First, one has M u n s ( t ) = Z t e − δ s ( t − s ) u n ( s ) ds for all t > M u n s ) n ∈ N is bounded in W , + ∞ (0 , T ). Thus, it converges upto a subsequence to M us strongly in C ([0 , T ]) according to the Ascoli theorem. By using the samereasoning as above, the sequence ( E u n , M u n , F u n ) is non-negative, uniformly bounded by above andtherefore, the right-hand side of the first three equations of ( S ) is bounded in W , + ∞ (0 , T ). Usingto the Ascoli theorem, we infer that ( F u n ) n ∈ N converges to F u in C ([0 , T ]). The desired conclusionfollows by passing to the limit in the inequality F u n ( T ) ε .Finally, the monotonicity of J ( S ) T,U,ε with respect to U comes from the monotonicity of the ad-missible control set U ( S ) T,U,ε with respect to U and T for the inclusion. More precisely, if U U then U ( S ) T,U ,ε ⊂ U ( S ) T,U ,ε according to Lemma 4.1. Moreover, if T T , let u ∗ a solution of Prob-lem ( P ( S ) T,U,ε ). Then, u = u ∗ ( · − T + T )( ( T − T ,T ) ∈ U T ,U,ε (since ( E, M, F ) is stationary on[0 , T − T )), and is such that J ( u ) = J ( S ) T ,U,ε . Hence, we get by minimality that J ( S ) T ,U,ε J ( S ) T ,U,ε .Finally, since u = U [ T − T ,T ] belongs to U ( S ) T,U,ε for any U > U ∗ and T > T ( U ), one has J ( S ) T,U,ε J ( u ) = U T .In the following result, one shows that the state constraint is reached by any optimal control.
Lemma 4.3.
Let
U > U ∗ and T > T ( U ) . Let u ∗ solving Problem ( P ( S ) T,U,ε ) and ( E, M, F, M s ) bethe corresponding solution to System ( S ) . Then, one has necessarily F ( T ) = ε and F ( · ) ∈ ( ε, F ] on [0 , T ) .Proof. Let u ∗ be a solution to Problem ( P ( S ) T,U,ε ). Since F is a stationary solution and F (0) = F , onehas F F on (0 , T ) by using Proposition 2.2. it remains to prove that F is bounded below by ε . Letus assume by contradiction that the corresponding solution ( E, M, F, M s ) to System ( S ) satisfies F ( T ) < ε . Let u η := (1 − η ) u ∗ with η ∈ (0 ,
1) and denote by ( E η , M η , F η , M sη ) the correspondingsolution to System ( S ). Then, by mimicking the arguments of the proof of Lemma 4.2, one easilygets that the mapping u ∈ U ( S ) T,U,ε ( E, M, F, M s ) ∈ C ([0 , T ]) is continuous for the weak- ∗ topology of L ∞ and it follows that F η ( T ) = F ( T ) + O( η ) so that F η ( T ) < ε whenever η is smallenough. Since R T u η ( t ) dt < R T u ∗ ( t ) dt , this contradicts the minimality of u ∗ .Finally, we claim that F ( · ) > ε on [0 , T ). In the converse case, there exists T ′ < T such that F ( T ′ ) = ε . One sees easily the control u T ′ defined by u T ′ = u ∗ ( · − T + T ′ ) [ T − T ′ ,T ] belongs to U ( S ) T,U,ε and that J ( u T ′ ) < J ( u ∗ ), whence a contradiction.15et us now state the necessary optimality conditions for Problem ( P ( S ) T,U,ε ). To this aim, let usintroduce (
P, Q, R, S ) as the solution to the backward adjoint system − ddt PQRS = ∂f E ∂E ( E, F ) ∂f M ∂E ( E, M ) ∂f F ∂E ( E, M, F, M s ) 00 ∂f M ∂M ( E, M ) ∂f F ∂M ( E, M, F, M s ) 0 ∂f E ∂F ( E, F ) 0 ∂f F ∂F ( E, M, F, M s ) 00 0 ∂f F ∂M s ( E, M, F, M s ) − δ s PQRS ,P ( T ) = 0 , Q ( T ) = 0 , R ( T ) = 1 , S ( T ) = 0 . (10) Lemma 4.4.
Let u ∈ U ( S ) T,U,ε and introduce the functional G defined by G ( u ) = F ( T ) , where ( E, M, F, M s ) denotes the unique solution of System ( S ) associated to u . Then, G is differen-tiable in the sense of Fr´echet and for every admissible perturbation h , the Gˆateaux-derivative of G at u in the direction h is DG ( u ) · h = Z T h ( t ) S ( t ) dt, where ( P, Q, R, S ) solves System (10) .Proof. The Fr´echet-differentiability of G is standard and follows from the differentiability of themapping u ∈ U ( S ) T,U,ε F , where ( E, M, F, M s ) denotes the unique solution of System ( S ), itselfderiving from a standard application of the implicit functions theorem combined with variationalarguments.Moreover, the Fr´echet-derivative ( ˙ E, ˙ M , ˙ F , ˙ M s ) of ( E, M, F, M s ) at u in the direction h solvesthe linearized problem ddt ˙ E ˙ M ˙ F ˙ M s = A ˙ E ˙ M ˙ F ˙ M s + h , ˙ E (0) = 0 , ˙ M (0) = 0 , ˙ F (0) = 0 , ˙ M s (0) = 0 , where the Jacobian matrix A reads A = ∂fE∂E ( E,F ) 0 ∂fE∂F ( E,F ) 0 ∂fM∂E ( E,M ) ∂fM∂M ( E,M ) 0 0 ∂fF∂E ( E,M,F,M s ) ∂fF∂M ( E,M,F,M s ) ∂fF∂F ( E,M,F,M s ) ∂fF∂Ms ( E,M,F,M s )0 0 0 − δ s . More precisely, we call “admissible perturbation” any element of the tangent cone T u, U ( S T,U,ε to the set U ( S ) T,U,ε at u . Recall that the cone T u, U ( S T,U,ε is the set of functions h ∈ L ∞ (0 , T ) such that, for any sequence of positive realnumbers ε n decreasing to 0, there exists a sequence of functions h n ∈ L ∞ (0 , T ) converging to h as n → + ∞ , and u + ε n h n ∈ U ( S ) T,U,ε for every n ∈ N (see e.g. [11]).
16y integration by parts, it follows that Z T * h , PQRS + = Z T * ddt ˙ E ˙ M ˙ F ˙ M s , PQRS + − Z T * A ˙ E ˙ M ˙ F ˙ M s , PQRS + = − Z T * ˙ E ˙ M ˙ F ˙ M s , ddt PQRS + + * ˙ E ˙ M ˙ F ˙ M s , PQRS + T − Z T * ˙ E ˙ M ˙ F ˙ M s , A T PQRS + = F ( T ) , which leads to the desired result.The following Lemma leads to a characterization of any optimal control. We chose here toprovide some quick explanations on the derivation of optimality conditions. It would have also beenpossible to use a shorter argument through the so-called Pontryagin Maximum Principle (PMP) andwe would have obtained the same result. Lemma 4.5.
Let
U > U ∗ and T > T ( U ) . Let u ∗ denote a solution to Problem ( P ( S ) T,U,ε ) . Thereexists λ > such that a.e. on { u ∗ = 0 } , one has λS ( t ) > , a.e. on { < u ∗ < U } , one has λS ( t ) = 0 , a.e. on { u ∗ = U } , one has λS ( t ) . Proof.
Let us introduce the Lagrangian function L associated to problem ( P ( S ) T,U,ε ), defined by L : ( u, Λ) ∈ U T,U × R J ( u ) − Λ( F ( T ) − ε ) , where U T,U := { u ∈ L ∞ (0 , T ) : 0 u ( · ) U } .By standard arguments, we get the existence of a Lagrange multiplier λ > u ∗ , λ )satisfies D u L ( u ∗ , λ ) · h ≥ h belonging to the tangent cone of the set U ( S ) T,U at u ∗ . Moreover,according to Lemma 4.3, we have necessarily F ( T ) = ε .Let t ∗ be a Lebesgue density-one point of { u ∗ = 0 } . Let ( H n ) n ∈ N be a sequence of measurablesubsets containing all t ∗ and such that H n is included in { u ∗ = 0 } . Let us consider h = H n andnotice that, by construction, u ∗ + ηh belongs to U ( S ) T,U whenever η is small enough. One has L ( u ∗ + ηh, λ ) > L ( u ∗ , λ ) , whenever η is small enough. Let us divide this inequality by η , and let η go to 0. By using Lemma 4.4,we obtain Z T h ( t ) dt + λ Z T h ( t ) S ( t ) dt > , which rewrites | H n | + λ R H n S ( t ) dt > . Dividing this inequality by | H n | and letting H n shrink to { t ∗ } as n → ∞ shows that 1 + λS ( t ) > { u ∗ = 0 } whence the first point of Lemma 4.5, according17o the Lebesgue Density Theorem. The proof of the third point is similar and consists in consideringperturbations of the form u ∗ − ηh , where h denotes a positive admissible perturbation of u ∗ supportedin { u ∗ = U } . Finally, the proof of the second point follows the same lines by considering bilateralperturbations of the form u ∗ ± ηh , where h denotes an admissible perturbation of u ∗ supported in { < u ∗ < M } .Let us now prove that λ >
0. We argue by contradiction, assuming that λ = 0. Then, theswitching function 1 + λS is necessarily constant, equal to 1, and we have therefore u ∗ = 0 in [0 , T ],which leads to a contradiction since the optimal trajectory has to satisfy F ( T ) = ε .Let us prove the remaining facts stated in Theorem 3.2. Proof of Theorem 3.2.
Let ε ∈ (0 , F ). According to Lemma 4.2, there exists U ∗ such that for each U > U ∗ , there exists T such for all T > T , Problem ( P ( S ) T,U,ε ) has a solution u ∗ . Let ( E, M F, M s ) bethe optimal trajectory. According to Lemma 4.3, the constraint “ F ( T ) = ε ” is reached. Lemma 4.5implies that on { S > − /λ } , one has necessarily u = 0. Since S is continuous and S ( T ) = 0, itfollows that there exists t ∈ (0 , T ) such that u ∗ = 0 on ( t , T ). Let us first point out that System ( S ) is monotone with respect to the control u : Lemma 4.6.
Let u , u ∈ L ∞ (0 , T ) such that u > u > (resp. u > u > ). Let us assume that (4) holds. Then, the corresponding solutions F , F to System ( S ) satisfy F F on (0 , T ) (resp. F < F on (0 , T ) ).Proof. This is an immediate consequence of the fact that ∂f∂M s F ∈ (0 , F ].Let us investigate the existence of an optimal control for Problem ( P ( S ) T,U,ε ). Lemma 4.7.
Let ε ∈ (0 , F ) . For every U > U ∗ , there exists T ( U ) > such that for all T > T ( U ) ,the set U ( S ) T,U,ε is nonempty and Problem ( P ( S ) T,U,ε ) has a solution u ∗ .Moreover, one has J ( S ) T,U,ε U T (where J ( S ) T,U,ε is defined by (6) ) and J ( S ) T,U,ε is non-increasingwith respect to T > T and to U > U ∗ .Proof. We first prove that U ( S ) T,U,ε is nonempty when
U > U ∗ and T is large enough. If u ( · ) = U > U ∗ ,we have seen in Proposition 2.2 that F ( t ) → t → + ∞ . Thus, for any ε ∈ (0 , F ), there exists T >
0, such that F ( T ) = ε . Hence, for T > T , u = U [ T − T ,T ] belongs to U ( S ) T,U,ε .Proceeding as in the proof of Lemma 4.2, we obtain existence of a solution to Problem ( P ( S ) T,U,ε ) byconsidering a minimizing sequence and showing that it is in fact compact. Finally, the monotonicityand the bound on J ( S ) T,U,ε are obtained exactly as in the end of the proof of Lemma 4.2.By mimicking the proof of Lemma 4.3 for Problem ( P ( S ) T,U,ε ), we can prove that the constraint F ( T ) ε is saturated. Lemma 4.8.
Let
U > U ∗ , T > T ( U ) and u ∗ be a solution to Problem ( P ( S ) T,U,ε ) . Let ( F, M s ) be theassociated optimal trajectory, solution to System ( S ) . Then, one has F ( T ) = ε and F ( · ) ∈ ( ε, F ] on [0 , T ) .
18e can also prove that F is non-increasing on (0 , T ). Lemma 4.9.
Let
U > U ∗ , T > T ( U ) and u ∗ be a solution to Problem ( P ( S ) T,U,ε ) . Let ( F, M s ) be theassociated optimal trajectory, solution to System ( S ) . Then, for every t ∈ (0 , T ) , one has F ′ ( t ) .Proof. Let us argue by contradiction, assuming that the conclusion is not true. Then, since F is C ,there exist 0 < θ < θ < T such that F ′ > θ , θ ).According to Lemma 4.6, F is non-increasing in a neighborhood of 0. Moreover, according toLemma 4.8, F decreases to ε in a neighborhood of T . Consequently, it is not restrictive to assumemoreover that F ′ ( θ ) = F ′ ( θ ) = 0, and F ′ , θ )Notice from the expression of f given in (1) that F ′ = f ( F, M s ) δ M γ s M s > φ ( F ) = (1 − ν ) β E ν E νβ E ν E − δ F (cid:0) β E FK + ν E + δ E (cid:1) δ F (cid:0) β E FK + ν E + δ E (cid:1) F. Let us show that there exist 0 < τ < τ < T such that F ( τ ) < F ( τ ) and M s ( τ ) > M s ( τ ).Indeed, there are two possibilities. Either there exists τ ∈ ( θ , θ ) such that M s ( θ ) = M s ( τ ), thenwe take τ = θ . Or for any t ∈ ( θ , θ ), M s ( θ ) < M s ( t ). In this latter case, we take τ ∈ ( θ , θ ) suchthat M s ( θ ) < M s ( τ ) < φ ( F ( τ )) (which is always possible since F ′ > θ , θ )). Then, since F (0) = F > F ( τ ) > F ( θ ) and F is continuous, there exists ˜ τ ∈ (0 , θ ) such that F (˜ τ ) = F ( τ ).Moreover, since M s > φ ( F ) on (0 , θ ), we have M s (˜ τ ) > φ ( F (˜ τ )) = φ ( F ( τ )) > M s ( τ ) > M s ( θ ).By continuity of M s , there exists τ ∈ (˜ τ , θ ) such that M s ( τ ) = M s ( τ ). Since we have F ′ , θ ), we deduce that F ( τ ) F (˜ τ ) = F ( τ ) and F ( τ ) = F ( τ ) since φ ( F ( τ )) M s ( τ ) = M s ( τ ) < φ ( F ( τ )).Then, we take u ( t ) = , for t ∈ (0 , τ − τ ) ,u ∗ ( t − τ + τ ) , for t ∈ ( τ − τ , τ ) ,u ∗ ( t ) , for t ∈ ( τ , T ) . Using Lemma 4.12, we have J ( u ) J ( u ∗ ). Moreover, if we denote by ( F u , M us ) the solution withthis control u , we have F u ( τ ) = F ( τ ) < F ( τ ) and M us ( τ ) = M s ( τ ) > M s ( τ ). Then, on ( τ , T ),we have M us > M s and since M s f ( F, M s ) is non-increasing ( F u ) ′ = f ( F u , M us ) f ( F u , M s ). Bycomparison with the solution to F ′ = f ( F, M s ) on ( τ , T ) and since F u ( τ ) < F ( τ ), we deduce that F u ( T ) < F ( T ) = ε which contradicts Lemma 4.8.Let us introduce ( Q, R ) as the solution of the forward adjoint system − ddt (cid:18) QR (cid:19) = ∂f∂F ( F, M s ) 0 ∂f∂M s ( F, M s ) − δ s ! (cid:18) QR (cid:19) ,Q ( T ) = 1 , R ( T ) = 0 . (11)Similarly to Problem ( P ( S ) T,U,ε ), any optimal control can be characterized by using the first ordernecessary optimality conditions, in terms of a switching function of the form t λR ( t ) with λ ≥ Lemma 4.10.
Let
U > U ∗ and T > T ( U ) . Consider u ∗ a solution to Problem ( P ( S ) T,U,ε ) . Then thereexists λ > such that a.e. on { u ∗ = 0 } , one has λR ( t ) > , a.e. on { < u ∗ < U } , one has λR ( t ) = 0 , a.e. on { u ∗ = U } , one has λR ( t ) . Lemma 4.11.
Let
U > U ∗ and T > T ( U ) . There exists t ∈ (0 , T ) such that u ∗ = 0 on ( t , T ) .Proof. According to Lemma 4.10, on the set { R > − /λ } , one has necessarily u = 0. We concludeusing the fact that R is continuous and R ( T ) = 0. Lemma 4.12.
Let
U > U ∗ and T > T ( U ) . There exist positive constants C , C , C > that donot depend on U and T , such that | M s | < C , < C Q C . Proof.
First, one has for all t ∈ [0 , T ], M s ( t ) = Z t u ∗ ( s ) e δ s ( s − t ) ds Z T u ∗ ( t ) dt T U .
From (11), Q solves − dQdt = ∂f∂F ( F, M s ) Q, Q ( T ) = 1 . Note that, on [0 , T ], all quantities are bounded: one has ε F ( · ) F (0) = F according toLemma 4.9. Hence, we infer the existence of µ > − µ ∂f∂F ( F, M s ) µ . Thus,integrating the inequality − µQ − dQdt µQ with Q ( T ) = 1, the conclusion follows easily by usinga Gronwall type argument. Lemma 4.13.
Let
U > U ∗ and T > T ( U ) . If < u ∗ < U on a non empty interval ( s , s ) then u ∗ satisfies (7) on ( s , s ) .Proof. According to Lemma 4.10, one has R ′ = 0 on ( s , s ). Differentiating the equation satisfiedby R in (11) yields (cid:18) ∂ f∂M s ∂F F ′ + ∂ f∂M s M ′ s (cid:19) Q + ∂f∂M s Q ′ = 0 . We deduce that (cid:18) ∂ f∂M s ∂F f + ∂ f∂M s ( u − δ s M s ) (cid:19) Q = ∂f∂M s ∂f∂F Q. According to Lemma 4.12, one has
Q > , T ]. Hence u is given by u = (cid:18) ∂ f∂M s (cid:19) − (cid:18) ∂f∂M s ∂f∂F + ∂ f∂M s δ s M s − ∂ f∂M s ∂F f (cid:19) , where f is given by (1). We thus recover (7). Lemma 4.14.
Let
U > U ∗ and T > T ( U ) .(i) If δ s > δ F . Then on each open interval of { R > − /λ } , any local extremum for R is aminimum.(ii) If U is large enough, then, on each open interval of { R < − /λ } , any local extremum for R isa maximum. roof. (i) By differentiating the equation on R in (11), we get − R ′′ = (cid:18) ∂ f∂M s ∂F f + ∂ f∂M s ( u ∗ − δ s M s ) − ∂f∂M s ∂f∂F (cid:19) Q − δ s R ′ on (0 , T ). (12)Let I be an open interval of { R > − /λ } (whenever it exists). Let τ ∈ I be a local extremumfor R , we have R ′ ( τ ) = 0. Then, thanks to Lemma 4.10, we have u ∗ = 0 on I and hence − R ′′ ( τ ) = (cid:18) ∂ f∂M s ∂F f − ∂ f∂M s δ s M s − ∂f∂M s ∂f∂F (cid:19) ( τ ) Q ( τ ) . We have seen in Lemma 4.12 that
Q >
0. Then, the sign of − R ′′ ( τ ) is the sign of U given by U = ∂ f∂M s ∂F f − ∂ f∂M s δ s M s − ∂f∂M s ∂f∂F . Let us compute U . Notice that f is of the form f ( F, M s ) = µF Λ − δ F F where Λ := 1 F + aF + M s ( αF + βF + γ ) (13)with some positive constants µ, a, α, β, γ . Observe that ∂ Λ ∂M s = − Λ ( αF + βF + γ ) , ∂ Λ ∂F = − Λ (2 F + a + M s (2 αF + β )) , and therefore one computes ∂f∂M s = − µF Λ ( αF + βF + γ ) ,∂f∂F = 2 µF Λ − µF Λ (2 F + a + M s (2 αF + β )) − δ F ,∂ f∂M s = 2 µF Λ ( αF + βF + γ ) ∂ f∂M s ∂F = − µF Λ (2( αF + βF + γ ) + F (2 αF + β )) , +2 µF Λ (2 F + a + M s (2 αF + β ))( αF + βF + γ ) . After straightforward but tedious computations, we find ∂ f∂M s ∂F f − ∂f∂M s ∂f∂F = − µF Λ (2 µF ( αF + βF + γ ) + µF (2 αF + β ))+ 2 µ F Λ (2 F + a + M s (2 αF + β ))( αF + βF + γ )+ δ F F Λ (cid:0) µF ( αF + βF + γ ) + µF (2 αF + β ) (cid:1) − µδ F F Λ (2 F + a + M s (2 αF + β ))( αF + βF + γ )+ µF Λ ( αF + βF + γ )(2 µF Λ − δ F − µ (2 F + a + M s (2 αF + β ))Λ F )= µ F Λ (cid:16) (2 F + a + M s (2 αF + β ))( αF + βF + γ ) − αF + β Λ (cid:17) + µF Λ δ F × (cid:16) αF + 2 βF + γ Λ − F (2 F + a + M s (2 αF + β ))( αF + βF + γ ) (cid:17) . / Λ from (13), we get ∂ f∂M s ∂F f − ∂f∂M s ∂f∂F = µ F Λ (( β − αa ) F + 2 γF + aγ )+ µF Λ δ F ( − αF + ( aα − β ) F − γF − aγF ))+ µF Λ δ F M s ( αF + βF + γ )( γ − αF ) . Finally, we obtain U = µ F Λ (( β − αa ) F + 2 γF + aγ )+ µF Λ δ F ( − αF + ( aα − β ) F − γF − aγ ))+ µF Λ δ F M s ( αF + βF + γ )( γ − αF ) − δ s M s µF Λ ( αF + βF + γ ) . Let us use that 2 δ s > δ F . Noting that f µF Λ δ F , we obtain U − µ F Λ αaF + µF Λ δ F ( − αF + ( aα − β ) F − γF )) (14) − µF Λ δ F M s ( αF + βF + γ )(2 αF + βF ) . We observe that this quantity is negative whenever β > aα , which is the case since a = K ( ν E + δ E ) β E , α = δ M γ s K (1 − ν ) ν E and β = 2 δ M γ s ( ν E + τ E )(1 − ν ) β E ν E . (ii) Let I an open interval of { R < − /λ } . Let τ ∈ I be a local extremum of R . Then, we have R ′ ( τ ) = 0. According to Lemma 4.10, one has u ∗ = U on I . Hence, from (12), one gets − R ′′ = (cid:18) ∂ f∂M s ∂F f + ∂ f∂M s ( U − δ s M s ) − ∂f∂M s ∂f∂F (cid:19) Q − δ s R ′ . (15)Recall that ε F ( · ) F (0) = F according to Lemma 4.9. By using also Lemma 4.12, we getthe existence of C > T and U such that ∂ f∂M s ( F, M s ) = ν (1 − ν ) β E ν E F δ M γ s ( βEFK + ν E + δ E ) (cid:0) (1 − ν ) ν E β E F + δ M γ s M s ( βEFK + ν E + δ E ) (cid:1) > ν (1 − ν ) β E ν E F δ M γ s ( ν E + δ E ) (cid:0) (1 − ν ) ν E β E F + δ M γ s C ( βEFK + ν E + δ E ) (cid:1) > C. (16)A similar reasoning shows that all the other terms in (15) are uniformly bounded with respectto T and U . Thus, we can find U (independent on T ) large enough such that the right-handside is positive. Then, R ′′ ( τ ) <
0, which implies R admits a local maximum at τ . Lemma 4.15.
Let us make the same assumption as in Proposition 2.2. Consider the dynamicalsystem ( S ) with u = 0 , i.e. dFdt = f ( F, M s ) , dM s dt = − δ s M s , where f is given in (1) . If F ′ ( τ ) > , then for all t > τ , we have F ′ ( t ) > . roof. This is a consequence of the fact that the set E := { ( F, M s ) ∈ R + × (0 , + ∞ ) , such that f ( F, M s ) > } is stable by the aforementioned dynamical system. Indeed, on ( R + ) , f ( F, M s ) > M s φ ( F ), where the function φ is obtained by solving the implicit equation f ( F, φ ( F )) = 0.If at some time τ >
0, a trajectory crosses the part of the boundary of E defined by the implicitequation, the vector field defining the right-hand side of the differential system reads[ f ( F ( τ ) , M s ( τ )) , − δ s M s ( τ )] ⊤ = [0 , − δ s M s ( τ )] ⊤ . This field is vertically directed, and it points inward for E . The conclusion is similar on the otherparts of the boundary of E , whence the stability of this zone. Proof of Theorem 3.3.
We separately deal with the characterization and uniqueness properties min-imizers.
Step 1: characterization of optimal controls.
Notice first that the sets { R > − /λ } and { R < − /λ } are open, as inverse image of continuous functions (and thus contain an interval). Bycombining Lemma 4.14, Lemma 4.11, and the fact that R is continuous with R ( T ) = 0, we get theexistence of ( t , t ) ∈ [0 , T ] such that 0 ≤ t ≤ t < T and • R > − /λ on (0 , t ) or R < − /λ on (0 , t ); • R = − /λ on ( t , t ); • R > − /λ on ( t , T ).Combined with Lemmas 4.10 and 4.13, we deduce the form of the optimal control in Theorem 3.3.Let us now prove that for T large enough, we have t > u ∗ = 0 on (0 , t ). Assume bycontradiction that u ∗ = U on (0 , t ) or t = 0. Recall that, according to Lemma 4.7, we have J ( S ) T,U,ε = Z T u ∗ ( t ) dt U T .
However, combining the expression of u ∗ given in (7) with the estimates (14) and (16) show that u ∗ ( t ) > ˜ C > , on ( t , t ) , where ˜ C does not depend on T nor U . Without loss of generality, we can assume that ˜ C < U ∗ .Hence, one has Z t u ∗ ( t ) dt > t ˜ C (since u ∗ = U > U ∗ > ˜ C on (0 , t )). It follows that t U T / ˜ C is uniformly bounded with respectto T . On ( t , T ), one has u ∗ = 0. Let F U denote the trajectory associated to the control choice u U = U [0 ,UT / ˜ C ] . According to the monotonicity property stated in Lemma 4.6, one has F > F U in R + since u u U . Furthermore, according to Proposition 2.2, F U ( T ) converges to the steady state F as T → + ∞ . Let T U > F U ( t ) > ( ε + F ) / > ε in ( T U , + ∞ ). If T > T U , onethus have F ( T ) > ε which contradicts the fact that F ( T ) = ε (see Lemma 4.8).Hence, there exists T ∗ > T > T ∗ , we have t >
0, and u ∗ = 0 on(0 , t ). Necessarily, t > t , otherwise u ∗ = 0 a.e. on (0 , T ) which is not possible, since in this case F ( t ) = F > ε on (0 , T ). 23et us notice that if
T > T ∗ , we have F ′ ( T ) > . (17)Indeed, since T > T ∗ , we have seen that the optimal control has the form u ∗ = u ∗ ( t ,t ) for0 < t < t < T . If F ′ ( T ) <
0, then there exists T > T such that F ( T ) < ε and T − T < t .Then by taking u = u ∗ ( t − T + T,t − T + T ) , we obtain a control on (0 , T ) such that J ( u ) = J ( u ∗ ) and F ( T ) < ε . It is not possible (see Lemma 4.8). Thus, F ′ ( T ) > T > T > T ∗ , since T J T,U,ε is non-increasing, then J ( S ) T ,U,ε J ( S ) T ,U,ε . Let us show that J ( S ) T ,U,ε = J ( S ) T ,U,ε . By contradiction, assume that J ( S ) T ,U,ε < J ( S ) T ,U,ε . Let us denote by u ∗ (resp. u ∗ ),the optimal solution on (0 , T ) (resp. (0 , T )). Then, from the results above, there exist t (1)0 < t (1)1 and t (2)0 < t (2)1 such that u ∗ = u ∗ ( · − t (1)0 ) ( t (1)0 ,t (1)1 ) and u ∗ = u ∗ ( · − t (2)0 ) ( t (2)0 ,t (2)1 ) with the sameexpression of u ∗ given in (7). From Z t (1)1 t (1)0 u ∗ ( t − t (1)0 ) dt = J ( S ) T ,U,ε < J ( S ) T ,U,ε = Z t (2)1 t (2)0 u ∗ ( t − t (2)0 ) dt, we deduce that t (1)1 − t (1)0 < t (2)1 − t (2)0 . Notice that, denoting F , resp. F , the solution to ( S )with u = u ∗ , resp. u = u ∗ , we have F ( t + t (1)0 ) = F ( t + t (2)0 ) for each t ∈ (0 , t (1)1 − t (1)0 ); and fromabove remark, we have F ′ ( T ) > F ′ ( T ) >
0, and, with Lemma 4.15, F ′ ( t ) > t > T . Then, F ( t (1)1 ) = F ( t (1)1 + t (2)0 − t (1)0 ) > F ( t (2)1 ). Since t F ( t + t (1)1 ) and t F ( t + t (2)1 ) verify thesame dynamical system on (0 , T − t (1)1 ), we deduce that F ( t + t (1)1 ) > F ( t + t (2)1 ) which implies inparticular that F ( T ) > F ( T + t (2)1 − t (1)1 ) > ε . It is a contradiction with the fact that F ( T ) = ε (see Lemma 4.8).It remains to show that F reaches its minimal value at T and that F ′ ( T ) = 0. Let us extendthe definition of F to R + by setting u ∗ = 0 in ( T, + ∞ ). We already know that F is non-increasingon [0 , T ] according to Lemma 4.9 and therefore, F ′ ( T ) ≤
0, which leads to the conclusion since F ′ ( T ) > Step 2: uniqueness property.
To conclude the proof, let us prove the uniqueness of the optimalcontrol if
T > T ∗ . We will use a constructive argument which is also used to derive an algorithmfor solving numerically this problem in Section 3.2.For ( t , t ) ∈ [0 , T ] such that t t , introduce ( F ( t ,t ) , M s ( t ,t ) ) as the solution of the Cauchyproblem ddt (cid:18) FM s (cid:19) = (cid:18) f ( F, M s ) u − δ s M s (cid:19) in [0 , + ∞ )with u ( t ,t ) = ∂f∂Ms ( F,M s ) ∂f∂F ( F,M s )+ ∂ f∂M s ( F,M s ) δ s M s − ∂ f∂Ms∂F ( F,M s ) f ( F,M s ) ∂ f∂M s ( F,M s ) [ t ,t ] , complemented with the initial conditions F (0) = F and M s (0) = 0.Let us first highlight that, with the notations of Theorem 3.3, it is enough to consider the casewhere t is equal to 0. Lemma 4.16.
Under the assumptions of Theorem 3.3 and if
T > T ∗ , u ( t ,t ) solves Problem ( P ( S ) T,U,ε ) if, and only, if the function ˜ u given by ˜ u ( · ) = u ( t ,t ) (cid:12)(cid:12) ( t ,T ) ( · − t ) solves Problem ( P ( S ) T − t ,U,ε ) . P ( S ) T,U,ε ) is of the form u ( t ,t ) for some 0 ≤ t < t ≤ T . Observe that Z T u ( t ,t ) ( t ) dt = Z T − t ˜ u ( t ) dt and furthermore, denoting respectively by F u ( t ,t and F ˜ u the trajectories associated to u ( t ,t ) and ˜ u ,one has by construction F u ( t ,t ( T ) = F ˜ u ( T − t ) = ε . Since T ∈ ( T ∗ , + ∞ ) J ( S ) T,U,ε is nondecreasing,it is constant on ( T − t , T ) and it follows that ˜ u necessarily solves Problem ( P ( S ) T − t ,U,ε ). The proofof converse sense is exactly similar (and left to the reader). This ends the proof of Lemma 4.16.Let us argue by contradiction to prove the uniqueness of optimizers, assuming that Prob-lem ( P ( S ) T,U,ε ) has two solutions u ( t ,t ) and u ( t ′ ,t ′ ) . According to Lemma 4.16 above, u (0 ,t − t ) and u (0 ,t ′ − t ′ ) solve Problems ( P ( S ) T − t ,U,ε ) and ( P ( S ) T − t ′ ,U,ε ) respectively. Without loss of generality, assumethat t ′ − t ′ t − t . Since these controls are non-negative and R (0 ,T ) u (0 ,t − t ) = R (0 ,T ) u (0 ,t ′ − t ′ ) , weinfer that the function ∂f∂M s ( F, M s ) ∂f∂F ( F, M s ) + ∂ f∂M s ( F, M s ) δ s M s − ∂ f∂M s ∂F ( F, M s ) f ( F, M s ) ∂ f∂M s ( F, M s )vanishes on ( t − t , t ′ − t ′ ). If t ′ − t ′ < t − t , a standard regularity argument for Cauchy problemsyields that u ,t − t ) is real analytic on (0 , t − t ), and vanishes hence identically on (0 , T ), so that F u (0 ,t − t ( · ) = F on R + , which is impossible since F u (0 ,t − t ( T − t ) = ε . Thus t ′ − t ′ = t − t and F u (0 ,t − t = F u (0 ,t ′ − t ′ on R + . F u (0 ,t − t is C on R + and one has lim t → + ∞ F u (0 ,t − t ( t ) = F , accord-ing to Proposition 2.2 and since the persistence equilibrium is the only one being not unstable. Letus denote by t the time at which F u (0 ,t − t admits its first local minimum. Thanks to Lemma 4.15, F u (0 ,t − t decreases on (0 , t ) and increases on ( t , + ∞ ). Since u ,t − t solves Problems ( P ( S ) T − t ,U,ε )and ( P ( S ) T − t ′ ,U,ε ), one has F ′ u (0 ,t − t ( T − t ) = F ′ u (0 ,t − t ( T − t ′ ) = F ′ ( t ) = 0. We deduce that T − t = T − t ′ = t , which shows the uniqueness of t and t . Notice first that F t max is C on (0 , + ∞ ). Denote by t max ∈ (0 , + ∞ ] the minimal time t at whichit holds either u ∞ ( t ) < F ∞ ) ′ ( t ) >
0. Its existence results from the form of optimal controls(see Theorem 3.3). Let us first prove that t max = + ∞ . Assume by contradiction that t max < + ∞ .According to Proposition 2.2 and since the persistence equilibrium is the only one being not unstable,one has lim t → + ∞ F t max ( t ) = F . Let τ be the time at which F t max admits its first local minimum.According to Lemma 4.15, F t max is decreasing on (0 , τ ) and increasing on ( τ , ∞ ). Hence, there exists δ > F t max ( · ) > δ . To reach a contradiction, let us consider a particular choice of ε , suchthat ε ∈ (0 , δ ). Let T > P ( S ) T,U,ε ) is well-posed (see Theorem 3.3)and let u ( t ,t ) be its unique solution. According to Lemma 4.16, u t − t solves Problem ( P ( S ) T − t ,U,ε ).Since u t − t is positive and ( F t − t ) ′ > , t − t ), one has t − t t max . By Lemma 4.6,one has F t − t > F t max > δ on R + , which is in contradiction with the fact F t − t ( T ) = ε . Thus t max = ∞ .Item (i) can be obtain with the same reasoning as above and the fact that t max = + ∞ . Property(ii) is a consequence of Lemma 4.6 and again of the equality t max = + ∞ . Finally, the proof of the25ast claim (iii), establishing connections between Problem ( P ( S ) T,U,ε ) under its general form and thecontrol u τ follows from Lemma 4.16, used to get the uniqueness of optimal controls in the proof ofTheorem 3.3. In this section, we introduce two other possible choices of functionals to minimize and compare itto the one analyzed in the previous sections. L functional Consider the optimal control problem inf u ∈U ( S T,U,ε ˜ J ( u ) , ( ˜ P ( S ) T,U,ε )where the functional ˜ J stands for the square of the L -norm of released mosquitoes over the horizonof time T , namely ˜ J ( u ) := Z T u ( t ) dt and U ( S ) T,U,ε is defined by (5).
Theorem 5.1.
Let ε ∈ (0 , F ) . For any U > U ∗ (defined by (3) ), there exists a minimal time T ( U ) > such that for all T > T ( U ) , the set U ( S ) T,U,ε is non empty and the optimal control problem ( ˜ P ( S ) T,U,ε ) has a solution u ∗ . Moreover, for U > U ∗ and T > T ( U ) and large enough, there exists µ > and t ∈ [0 , T ) such that u ∗ = (cid:26) U on [0 , t ) − µR, on [ t , T ] where R solves the dual system (11) . Moreover u ∗ ( T ) = 0 and the mappings T ∈ [ T , + ∞ ) ˜ J ( S ) T,U,C and U ∈ ( U ∗ , + ∞ ) ˜ J ( S ) T,U,C are non-increasing.
Remark 5.2.
Similarly to the L case, we have reduced the infinite dimensional control problem tofinite dimensional one, and we therefore only need to determine µ and the value of Q and R at time to see ( F, M s , Q, R ) as the solution of a well-posed Cauchy problem. To prove Theorem 5.1, we first need the equivalent of Lemma 4.10.
Lemma 5.3.
Consider u ∗ a solution to Problem ( ˜ P ( S ) T,U,ε ) . Then there exists λ > such that a.e. on { u ∗ = 0 } , one has u ∗ ( t ) + λR ( t ) > , a.e. on { < u ∗ < U } , one has u ∗ ( t ) + λR ( t ) = 0 , a.e. on { u ∗ = U } , one has u ∗ ( t ) + λR ( t ) . The proof is similar to Lemma 4.10 and will be omitted.
Proof of Theorem 5.1.
From (11), we deduce that R ′ = − ∂f∂M s Q + δ s R > δ s R, Q > ∂f /∂M s <
0. Since one has R ( T ) = 0, we infer that R < , T ) by using a standard Gronwall argument. According to Lemma 5.3, one has { u ∗ = 0 } = ∅ .According to Lemma 4.12, since 0 ∂f /∂M s ( M s , F ) O( | F | ), one infers that R is uniformlybounded on [0 , T ], by a constant that do not depend on U and T . Then, by adapting the proof ofLemma 4.14, one shows that for U large enough, on any open interval of the open set { U + λR < } ,any local extremum for R is a maximum. It follows that { U + λR < } has at most one connectedcomponent of the form (0 , t ), which leads to the conclusion.We provide on Figure 5 the solutions of Problem ( ˜ P ( S ) T,U,ε ) with T = 200, U = 4000, ν E = 0 . ε = F /
4. We recover the theoretical result above, namely that the optimal control u ∗ is positive,continuous and u ∗ ( T ) = 0. . . . . · tOptimizer F to ( ˜ P ( S ) T,U,ε ) ε , , , ,
000 tOptimizer u to ( ˜ P ( S ) T,U,ε ) U Figure 5: Solution of the optimal control problems ( ˜ P ( S ) T,U,ε ) with T = 200, U = 4000, ν E = 0 .
05 and ε = F / Remark 5.4.
In optimal control theory, the L -norm is often preferred to the L -norm for differ-entiability issues. However, from a biological point of view, the L -norm is more relevant since itstands for the amount of individuals. Moreover, as it can be seen on Figure 2, the optimal controlfor the L -norm is sparse unlike the one for the L -norm, which is interesting from a practical pointof view. Consider the optimal control problem inf u ∈U ( S T,U,C ˆ J ( u ) , ( ˆ P ( S ) T,U,C )where the functional ˆ J stands for the total number of eggs and females (with some weights) at time T , namely ˆ J ( u ) := F ( T ) , where F solves System ( S ) associated to the control u and U ( S ) T,U,C is the set of admissible controls,chosen so that: • the rate of sterile male mosquito release is non-negative, uniformly bounded by a positiveconstant U ; 27 the total number of released sterilized males over the time interval (0 , T ) is assumed to belower than C .Hence, U ( S ) T,U,C is defined by U ( S ) T,U,C := n u ∈ L ∞ (0 , T ) : 0 u U a.e. in (0 , T ) , Z T u ( t ) dt C o . In [1], a close optimal control problem has been considered. Problem ( ˆ P ( S ) T,U,C ) can be seen as adual version of ( P ( S ) T,U,ε ) in the following sense: let u ∗ ∈ U ( S ) T,U,ε be a solution to Problem ( P ( S ) T,U,ε ) forsome given T , U and ε . Then, the control u ∗ is a solution to Problem ( ˆ P ( S ) T,U,C ) for the parameterchoice C := R T u ∗ ( t ) dt . Indeed, assume by contradiction that there exists u ∈ U ( S ) T,U,C such thatˆ J ( u ) < ˆ J ( u ∗ ) , that is F ( T ) < F ∗ ( T ) = ε , where ( F, M s ) and ( F ∗ , M ∗ s ) are the solutions to system ( S ) associatedto u and u ∗ respectively. By mimicking the argument provided in the proof of Lemma 4.3, we reacha contradiction, which shows that u ∗ is a solution to Problem ( ˆ P ( S ) T,U,C ).Respectively, let ˆ u ∗ ∈ U ( S ) T,U,C be an optimizer to Problem ( ˆ P ( S ) T,U,C ) for some given T , U and C . By using the same argument, one shows that ˆ u ∗ is an optimizer to Problem ( P ( S ) T,U,ε ) for theparameter choice ε := F ∗ ( T ). In this paper, we have determined the optimal release function which minimizes the number ofsterilized males needed when performing the sterile insect technique (SIT) to reduce the size ofa population of mosquitoes to a given value. Starting from a differential system modeling thedynamics of the mosquitoes population, we simplify it to obtain a reduced system, which is agood approximation, and for which we are able to compute precisely the optimal solution. Thesetheoretical results are illustrated thanks to some numerical simulations. Notice that once the formof the theoretical solution is known, efficient algorithms may be designed to compute quickly thenumerical solution.Obviously, when the final number of mosquitoes is fixed, there is a minimal time to performthe releases in order to reach this value. Interestingly, the number of sterilized males needed isnon-increasing with respect to the time of the experiment, meaning that the longer the duration ofthe experiment, the lower the number of sterilized males. However, there is a maximal time abovewhich the minimal number of sterilized males needed is stationary. In this case, the optimal releasefunction is given by a singular arc sandwiches between two regions where it is zero. The knowledgeof the existence of this time may be interesting for practical applications since for larger time thenumber of sterilized males keeps constant.Thanks to our results, we are able to give a precise description of the time repartition of thereleases to optimize some scenarios. Then a natural extension of this work is to add the spacialrepartition of mosquitoes. It may have a big impact of the success of the SIT. In particular, mostexperiments have been performed in isolated region to avoid re-invasion from the outside. But evenin isolated regions the question to know where to perform the releases to have the best efficiency ofthe SIT in the region is still open.
Appendix Mathematical properties of the dynamical systems
Proof of Proposition 2.1
Let us assume that u ( · ) = 0. Then, the equilibria ( E, M , F , M s ) of System ( S ) solve0 = β E F (cid:18) − EK (cid:19) − (cid:0) ν E + δ E (cid:1) E = (1 − ν ) ν E E − δ M M = νν E E MM + γ s M s − δ F F = − δ s M s . We thus infer that M s = 0 and0 = β E F (cid:18) − EK (cid:19) − (cid:0) ν E + δ E (cid:1) E = (1 − ν ) ν E E − δ M M = νν E E − δ F F .
Then (0 , , ,
0) is an equilibrium, and the only non-zero equilibrium is E = K − δ F (cid:0) ν E + δ E (cid:1) β E νν E ! , M = (1 − ν ) ν E δ M E, F = νν E δ F E, M s = 0whence (2).Let us show that (0 , , ,
0) is unstable. Using that M s ( t ) = e − δ s t M s (0) for t >
0, we deduce that M ( t ) > e − δ M t M (0) for t > H ), and it follows that for any ǫ >
0, there exists t ∗ > γ s M s ( t ) < ǫM ( t ) for all t > t ∗ . By using a standard comparison principle, we get that forall t > t ∗ ( E ( t ) , F ( t )) > ( E ( t ) , F ( t )) , the inequality being understood component by component, where, ( E , F ) solves dE dt = β E F (cid:18) − E K (cid:19) − (cid:0) ν E + δ E (cid:1) E , t > t ∗ dF dt = νν E ǫ E − δ F F (18)complemented with the initial data ( E ( t ∗ ) , F ( t ∗ )) = ( E ( t ∗ ) , F ( t ∗ )). An easy computation yieldsthat the Jacobian matrix of System (18) at (0 ,
0) is (cid:18) − ν E − δ E β Eνν E ǫ − δ F (cid:19) , whose determinant expands as δ F ( ν E + δ E ) − νβ E ν E + O( ǫ ). According to ( H ), we infer that theJacobian matrix has a positive root whenever ǫ is chosen small enough, which leads to the conclusion.Let us now investigate the stability of ( E, M , F ,
0) for System ( S ). Easy computations yield29hat the Jacobian matrix of System ( S ) at ( E, M , F ,
0) reads − β E K F − ( ν E + δ E ) 0 β E (cid:16) − EK (cid:17) − ν ) ν E − δ M νν E − δ F − γ s νν E EM − δ s = − νν E β E δ F δ F ( ν E + δ E ) νν E − ν ) ν E − δ M νν E − δ F − γ s νδ M − ν − δ s . so that the four eigenvalues are − δ s , − δ M , and the two (complex conjugate) roots of the polynomial P = X + (cid:18) νν E β E δ F + δ F (cid:19) X + ( νν E β E − δ F ( ν E + δ E )) , which have a negative real part under ( H ). It follows that ( E, M , F ,
0) is linearly asymptoticallystable.Let us now prove the second part of the proposition. We first notice that the set [0 , K ] × R is stable whenever u is non-negative. We claim that System ( S ) is monotone on this set (seeLemma 4.1). If u belongs to L ∞ (0 , T, R + ), then (0 , , ,
0) is obviously a subsolution of System ( S )whereas ( E, M , F , k u k ∞ /δ s ) is a supersolution. A comparison argument allows us to conclude that[0 , E ] × [0 , M ] × [0 , F ] × R + is stable. Furthermore, if all initial data are positive, then so are thefunctions E , M , F and M s , and we get that E ( t ) > e − ( ν E + δ E ) t E (0) , M ( t ) > e − δ M t M (0) , F ( t ) > e − δ F t F (0)for all t ≥
0, implying that these quantities cannot vanish.Finally, let us consider the case where u ( · ) = U , where U > U ∗ . In that case, the non-zeroequilibria of System ( S ) solve M s = Uδ S , M = (1 − ν ) ν E δ M E, F = νν E δ F E MM + γ s M s , and 0 = β E F (cid:18) − EK (cid:19) − ( ν E + δ E ) E plugging the three first above equalities into this latter equation, we get that E satisfies E = 0or β E ν (1 − ν ) ν E δ F δ M K E − β E ν (1 − ν ) ν E δ F δ M (cid:18) − δ F ( ν E + δ E ) β E νν E (cid:19) E + γ s ( ν E + δ E ) δ s U = 0 . One easily checks that this second order polynomial with unknown E does not have any real solutionsif U > U ∗ . In this case, the only equilibrium is the extinction equilibrium. We infer that any non-negative solution to the Cauchy problem converges to this unique steady state.30 roof of Proposition 2.2 Let us assume that u ( · ) = U , the equilibria are obtained by solving the system0 = f ( F , M s ) = U − δ s M s , which is equivalent to M s = U /δ s and F (cid:18) ν (1 − ν ) β E ν E F − δ F (cid:0) β E FK + ν E + δ E (cid:1)(cid:16) (1 − ν ) ν E β E F + δ M γ s M s (cid:0) β E FK + ν E + δ E (cid:1)(cid:17)(cid:19) = 0 . It follows that if U = 0, then there are exactly two different solutions for this latter system, namely F = 0 or F given by (2). A straightforward computation yields that if U > U ∗ , then the equationabove has no positive solution and furthermore, f ( F, M s ) < F >
0. We infer that anynon-negative solution converges to the steady state (0 , M s ) if U > U ∗ .Let us assume that u ( · ) = 0. Using that M s ( t ) = e − δ s t M s (0) , F ( t ) > e − δ F t F (0)for t > δ s > δ F , we get that, for all η >
0, there exists t ∗ > δ M γ s ( β E F ( t ) K + ν E + δ E ) M s ( t ) < η (1 − ν ) ν E β E F ( t ) for all t > t ∗ . It follows that f ( F ( t ) , M s ( t )) > ˜ f ( F ( t )) := νβ E ν E F ( t ) (cid:0) β E F ( t ) K + ν E + δ E (cid:1) (1 + η ) − δ F F ( t ) , whenever t > t ∗ . We conclude by observing that ˜ f ′ (0) > η small enough, so that we get theinstability of the equilibrium (0 , F ,
0) of System ( S ). Onecomputes ∂f∂F ( F , ν (1 − ν ) β E ν E F (cid:0) βEFK + ν E + δ E (cid:1) (1 − ν ) ν E β E F − F (1 − ν ) ν E β E (cid:0) βEFK + ν E + δ E (cid:1)(cid:0) βEFK + ν E + δ E (cid:1) (cid:0) (1 − ν ) ν E β E F (cid:1) − δ F = δ F (cid:16) δ F νβ E ν E ( ν E + δ E ) − (cid:17) , which is negative under condition ( H ). Because of the form of System ( S ), it follows that ( F ,
0) isa linearly asymptotically stable steady state. Finally, a standard comparison argument for Cauchyproblems yields that if 0 < F (0) < F and u ( · ) >
0. Then, we have 0 < F ( t ) < F for all t > References [1] L. Almeida, M. Duprez, Y. Privat, and N. Vauchelet. Mosquito population control strategiesfor fighting against arboviruses.
Mathematical Biosciences and Engineering , 16(6):6274, 2019.312] R. Anguelov, Y. Dumont, and I. V. Y. Djeumen. Sustainable vector/pest control using thepermanent sterile insect technique.
Mathematical Methods in the Applied Sciences. Online. ,pages 1–22, 2020.[3] R. Anguelov, Y. Dumont, and J. Lubuma. Mathematical modeling of sterile insect technologyfor control of anopheles mosquito.
Comput. Math. Appl. , 64(3):374–389, 2012.[4] M. S. Aronna and Y. Dumont. On nonlinear pest/vector control via the sterile insect technique:impact of residual fertility. arXiv preprint arXiv:2005.05595 , 2020.[5] H. Barclay and M. Mackauer. The sterile insect release method for pest control: a density-dependent model.
Environmental Entomology , 9(6):810–817, 1980.[6] L. Beal, D. Hill, R. Martin, and J. Hedengren. Gekko optimization suite.
Processes , 6(8):106,2018.[7] P.-A. Bliman. Feedback control principles for biological control of dengue vectors. , 2019.[8] P.-A. Bliman, D. Cardona-Salgado, Y. Dumont, and O. Vasilieva. Implementation of controlstrategies for sterile insect techniques.
Math. Biosci. , 314:43–60, 2019.[9] P. A. Bliman, D. Cardona-Salgado, Y. Dumont, and O. Vasilieva. Optimal control approachfor implementation of sterile insect techniques. arXiv preprint arXiv:1911.00034 , 2019.[10] L. Cai, S. Ai, and J. Li. Dynamics of mosquitoes populations with different strategies forreleasing sterile mosquitoes.
SIAM J. Appl. Math. , 74(6):1786–1809, 2014.[11] R. Cominetti and J.-P. Penot. Tangent sets to unilateral convex sets.
C. R. Acad. Sci. ParisS´er. I Math. , 321(12):1631–1636, 1995.[12] C. Dufourd and Y. Dumont. Impact of environmental factors on mosquito dispersal in theprospect of sterile insect technique control.
Comput. Math. Appl. , 66(9):1695–1715, 2013.[13] Y. Dumont and J. M. Tchuenche. Mathematical studies on the sterile insect technique for theChikungunya disease and aedes albopictus . J. Math. Biol. , 65(5):809–854, 2012.[14] V. A. Dyck, J. Hendrichs, and A. Robinson.
Sterile insect technique: principles and practice inarea-wide integrated pest management . Springer, 2006.[15] L. Esteva and H. M. Yang. Mathematical model to assess the control of aedes aegypti mosquitoesby the sterile insect technique.
Math. Biosci. , 198(2):132–147, 2005.[16] K. R. Fister, M. L. McCarthy, S. F. Oppenheimer, and C. Collins. Optimal control of insectsthrough sterile insect release and habitat modification.
Math. Biosci. , 244(2):201–212, 2013.[17] J. Hedengren, J. Mojica, W. Cole, and T. Edgar. Apopt: Minlp solver for differential andalgebraic systems with benchmark testing. In
Proceedings of the INFORMS National Meeting,Phoenix, AZ, USA , volume 1417, page 47, 2012.[18] M. W. Hirsch and H. Smith. Monotone dynamical systems. In
Handbook of differential equations:ordinary differential equations , volume II, pages 239–257. Elsevier B. V., Amsterdam, 2005.[19] M. Huang, X. Song, and J. Li. Modelling and analysis of impulsive releases of sterile mosquitoes.
Journal of biological dynamics , 11(1):147–171, 2017.3220] J. Li and Z. Yuan. Modelling releases of sterile mosquitoes with different strategies.
Journal ofbiological dynamics , 9(1):1–14, 2015.[21] B. Stoll, H. Bossin, H. Petit, and M. Marie, J.and Cheong Sang. Suppression of an isolatedpopulation of the mosquito vector aedes polynesiensis on the atoll of tetiaroa, french polynesia,by sustained release of wolbachia-incompatible male mosquitoes. In
Conference: ICE - XXVInternational Congress of Entomology, At Orlando, Florida, USA. , 2016.[22] M. Strugarek, H. Bossin, and Y. Dumont. On the use of the sterile insect release technique toreduce or eliminate mosquito populations.
Appl. Math. Model. , 68:443–470, 2019.[23] R. Thom´e, H. Yang, and L. Esteva. Optimal control of aedes aegypti mosquitoes by the sterileinsect technique and insecticide.
Math. Biosci. , 223(1):12–23, 2010.[24] X. Zheng et al. Incompatible and sterile insect techniques combined eliminate mosquitoes.