Contents lists available at ScienceDirect Mathematical Biosciences journal homepage: www.elsevier.com/locate/mbs An approach of the exact linearization techniques to analysis of population dynamics of the mosquito Aedes aegypti Célia A. dos Reisa, Helenice de O. Florentinob,⁎, Diego Cólonc, Suélia R. Fleury Rosad, Daniela R. Cantaneb a Departamento de Matemática, FC, UNESP, 17033-360, Bauru, SP, Brazil bDepartamento de Bioestatística, IB, UNESP, 18618-670, Botucatu, SP, Brazil c LAC - PTC, POLI, USP, 05508-900, São Paulo, SP, Brazil dUniversidade de Brasília UnB, Brazil A R T I C L E I N F O Keywords: Aedes aegypti Dynamic model Exact linearization techniques Dengue Chikungunya Zika A B S T R A C T Dengue fever, chikungunya and zika are caused by different viruses and mainly transmitted by Aedes aegypti mosquitoes. These diseases have received special attention of public health officials due to the large number of infected people in tropical and subtropical countries and the possible sequels that those diseases can cause. In severe cases, the infection can have devastating effects, affecting the central nervous system, muscles, brain and respiratory system, often resulting in death. Vaccines against these diseases are still under development and, therefore, current studies are focused on the treatment of diseases and vector (mosquito) control. This work focuses on this last topic, and presents the analysis of a mathematical model describing the population dynamics of Aedes aegypti, as well as present the design of a control law for the mosquito population (vector control) via exact linearization techniques and optimal control. This control strategy optimizes the use of resources for vector control, and focuses on the aquatic stage of the mosquito life. Theoretical and computational results are also presented. 1. Introduction The Aedes aegypti mosquito is currently considered one of the most dangerous creatures on the planet due to their ability to spread deadly diseases, their resistance, adaptability and proximity of human popu- lations, especially in tropical and subtropical countries. In addition to the yellow fever, the Aedes aegypti transmits viruses that cause dengue fever, chikungunya and zika [1]. Dengue fever is caused by an arbovirus of the flaviviridae family, which includes five virus serotypes, called DEN 1, DEN 2, DEN 3, DEN 4 and DEN 5. The most common symptoms are fever, headache, muscle and joint pain and skin rash. Sometimes it can progress into a life- threatening disease, such as dengue hemorrhagic fever. The chi- kungunya is caused by CHIKV virus, arboviruses of the alphavirus genus. The symptoms are similar to those of dengue, though a portion of the infected people with this virus can develop the chronic form and remain with the symptoms for up to a year. Zika fever, or simply zika, is an infectious disease caused by the ZIKV virus. In most cases there are no symptoms, and, when present, they are usually mild and can resemble dengue fever. It has not caused any reported deaths during the initial infection. However, mother-to-child transmission during pregnancy can cause microcephaly and other brain malformations in some babies. In the year 2015, it was observed a large number of cases in children that were born with microcephaly in Brazil, which were associated with the mother's infection by zika virus, called congenital zika syndrome. Infections in adults have been linked to Guillain–Barré syndrome (GBS), [2–4]. There is no fully developed vaccine and no specific treatment for these diseases, therefore the techniques of vector population control are widely used, though not always effectively, whether for reasons of environmental contamination of chemical controls or the high eco- nomic cost of the biological, genetic and mechanical controls. For a correct application of control, it is necessary to understand the popu- lation dynamics of the mosquito [1,5]. In this work, we propose to study the population dynamics of Aedes aegypti mosquitoes and apply vector control techniques for the aquatic population phase of the mosquitoes (that are egg, larva and pupa, with interventions for reduction of aquatic population). The mathematical model addressed here is based in [6], and we proposed to apply the technique of exact feedback linearization to linearize this nonlinear https://doi.org/10.1016/j.mbs.2017.12.001 Received 29 March 2017; Received in revised form 29 September 2017; Accepted 8 December 2017 ⁎ Corresponding author. E-mail addresses: celia@fc.unesp.br (C.A. dos Reis), helenice@ibb.unesp.br (H.d.O. Florentino), diego@lac.usp.br (D. Cólon), suelia@unb.br (S.R.F. Rosa). Mathematical Biosciences 299 (2018) 51–57 Available online 15 December 2017 0025-5564/ © 2017 Elsevier Inc. All rights reserved. T http://www.sciencedirect.com/science/journal/00255564 https://www.elsevier.com/locate/mbs https://doi.org/10.1016/j.mbs.2017.12.001 https://doi.org/10.1016/j.mbs.2017.12.001 mailto:celia@fc.unesp.br mailto:helenice@ibb.unesp.br mailto:diego@lac.usp.br mailto:suelia@unb.br https://doi.org/10.1016/j.mbs.2017.12.001 http://crossmark.crossref.org/dialog/?doi=10.1016/j.mbs.2017.12.001&domain=pdf model, and design a control law for the purpose of reducing the adult population by acting in the aquatic phase of the mosquito life. An op- timization procedure is also used in order to minimize the cost of the control. The exact feedback linearization is a procedure that can transform the dynamics of a nonlinear system in a linear dynamics through a nonlinear feedback of the states or outputs (that are functions of the states). Furthermore, it is an essential part for the development of nonlinear robust and adaptive controllers [7–9]. This nonlinear dynamic analysis methodology has been used suc- cessfully in a wide range of applications, such as tracking problems, control of robot arms and manipulators, artillery, helicopters, airplanes and satellites, moreover, has been used in medical equipment, phar- maceutical and chemical industries [7,9–11]. This work is organized as follows. Section 2 presents an introduction of the life cycle of the mosquito. Section 3 presents the mathematical model which describes the population dynamics of the mosquito po- pulation. Section 4 discusses the feedback linearization technique and presents concepts of internal and zero dynamics and normal form. In Section 5, the system is transformed to the normal form. In Section 6, the two open loop equilibrium points are determined, and two different exact feedback linearization laws are determined for those equilibrium points. Section 7 presents a optimization model in order to determine a control strategy that use a minimum quantity possible of aquatic con- trol and to keep the mosquito population as low as possible. Further- more, a VNS algorithm is proposed to solve the optimization model. In addition, we discuss some basic concepts and present the mathematical model describing the population dynamics of the mosquito. Finally, computational results are presented. 2. Population dynamics of the mosquito Aedes aegypti The life cycle of mosquito Aedes aegypti is a process that lasts 1.5 week (periods of intense sunlight) or 3 weeks (cold periods), and comprises four stages: egg, larva, pupa and adult. The life cycle varies with temperature, food availability and quantity of larvae in breeding, and can be divided into two phases: aquatic phase (egg, larva and pupa) and terrestrial phase (adult stage). In favorable environmental condi- tions, after the egg hatching, the development of mosquitoes into the adult form can take in up to 10 days. Therefore, mosquito control in the aquatic phase is a good practice to interrupt the mosquito life cycle. Fig. 1 presents an illustration of the mosquito life cycle. Males and females of the Aedes aegypti in adult phase feed on sugary substances such as nectar and sap of plants. Only the female bites hu- mans, because they need blood for development and maturation of their eggs. It is in the adult stage that occur the transmission of several viruses that cause diseases, such as yellow fever, dengue fever, zika fever and chikungunya. The larva develops within clear or dirty standing water, and the population control is difficult due to its ver- satility in the choice of breeding to lay its eggs, and the extremely eggs’ resistance (they can survive for several months until water becomes available to provides incubation). Once immersed, the eggs develop rapidly into larvae, which give rise to pupae, from which the adult emerges [5]. Aiming to assist the study and understanding the behavior of the mosquito population dynamics, a mathematical model is presented and discussed in Section 3. 2.1. Mathematical modeling Considering the life cycle of mosquitoes divided into aquatic and terrestrial phases, the objective here is to study the population control of Aedes aegypti in aquatic phase. Therefore, based on [6], we present the mathematical model in Eq. (1) describing the population dynamics of mosquitoes: ⎧ ⎨ ⎪⎪ ⎩ ⎪ ⎪ = − + + + − − = − + = − = − − ( )x γ μ x ϕ x x x x u t x rγ x β μ x x β x μ x x r γ x μ x ˙ ( ) ( ) ˙ ( ) , ˙ ˙ (1 ) ϕ C1 1 1 3 1 3 1 2 1 2 2 3 2 3 3 4 1 4 4 (1) where x1(t) is the density of individuals of the aquatic population phase (eggs, larvae and pupae), x2(t) is the density of individuals of the im- mature females population (before mating), x3(t) is the density of in- dividuals of the fertilized females population (after mating) and x4(t) the density of male individuals (natural male). Mortality rates asso- ciated with each segment of the population (water, immature females, fertilized females and males) are respectively denoted by μ1, μ2, μ3 and μ4. The oviposition rate of fertilized female is proportional to its density and also depend on the availability of food, which is given by −ϕ (1 )x C 1 being ϕ the rate of intrinsic oviposition and C means the environmental capability related to the number of nutrients, space, etc. γ is the per capita rate at which mosquitoes in aquatic phase pass to the terrestrial phase, and a ratio r develop in females and (1 - r) in males. β is the mating rate of immature females with male mosquitoes. The control variable u(t) acts as a time-varying mortality rate, that must be de- signed in order to control the population of mosquitoes. Without loss of generality, the following changes in the Model (1) are proposed: = − + = − = = − + = − A γ μ A ϕ C A rγ A β μ A r γ ( ), , , ( ), (1 ) . 1 1 2 3 4 2 5 Then the system of Eq. (1) can be written as: Fig. 1. Life cycle of Aedes aegypti: aquatic phase (larvae, pupae, eggs) and terrestrial phase (adults). C.A. dos Reis et al. Mathematical Biosciences 299 (2018) 51–57 52 ⎡ ⎣ ⎢ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ ⎥ = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ + + + − − ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ + ⎡ ⎣ ⎢ ⎢ − ⎤ ⎦ ⎥ ⎥ x x x x A x ϕx A x x A x A x βx μ x A x μ x x u t ˙ ˙ ˙ ˙ 0 0 0 ( ), 1 2 3 4 1 1 3 2 1 3 3 1 4 2 2 3 3 5 1 4 4 1 (2) i.e. = +x f x g x u˙ ( ) ( ) , (2a) where, = ⎡ ⎣ ⎢ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ ⎥ = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ + + + − − ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ = ⎡ ⎣ ⎢ ⎢ − ⎤ ⎦ ⎥ ⎥ = x x x x x f x A x ϕx A x x A x A x βx μ x A x μ x g x x u u t ˙ ˙ ˙ ˙ ˙ , ( ) , ( ) 0 0 0 , and ( ). 1 2 3 4 1 1 3 2 1 3 3 1 4 2 2 3 3 5 1 4 4 1 (2b) 3. Exact linearization techniques and zero dynamics The exact feedback linearization is a procedure that can transform the dynamics of a nonlinear system in a linear dynamics through a nonlinear feedback of the states or outputs previously chosen. The ad- vantage of this linearization, when compared with classical Taylor series linearization (or Jacobian), is that feedback linearization is global, i.e. it can be applied to all state space domain or output, with the possible exception of isolated points, while the Taylor series lineariza- tion is local, i.e., is valid only for a region around a given equilibrium point. Furthermore, while the Jacobian linearization is approximated, the feedback linearization is exact [7–9]. We also consider that an output variable y= h(x), that is a function of the states (population densities), can be measured in real-time by some way (using sensors). Thus, the mathematical model is of SISO (single input single output) type, and is presented in Eq. (3) below: = + = x f x g x u y h x ˙ ( ) ( ) ( ), (3) where x ∈ ℜn is a state vector, u ∈ ℜ is the control input, f, g: D⊆ℜn→ ℜn are smooth vector fields in D ⊆ ℜn. For exact feedback linearization, there are two main methods: input-to-state linearization and input-to-output linearization. The last one is the object of study in this paper. The aim in input-output line- arization is to transform the nonlinear dynamic system (3) in the so called normal form, that is, a linear outer part (input-output sub-model) and a nonlinear inner and non-observable part (also called internal dynamic). For this purpose, a new set of states is defined (that are the output y and its time derivatives) and the existence of a diffeo- morphism, which transforms the nonlinear system in another linear system, is supposed. The new outer states are =y z1, …, =y zr r ( ) , where r is the relative degree of the system. The internal dynamics associated with the input- output linearization corresponds to (n – r) equations of the form =ψ w z ψ˙ ( , ) and is part of the normal form, and represents the un- observable dynamics of the system. Note that it does not affect the system output or its derivatives, it is affected by the output y and its derivatives. Let us consider the subset of the state space which corre- sponds to its outputs identically zero, i.e. =y t( ) 0r( ) . This is called the zero dynamics of the system, and corresponds to a sub-manifold Mo = {X ∈ℜn: z = 0} of the state space, such that all state trajectories within it will not affect the output. In other words, the system restricted to Mo corresponds to the internal dynamics which describes the re- stricted movement of the system in the sub-manifold Mo= {X ∈ℜn: z= 0}. In addition, the system dynamics in this manifold cannot be changed by the control law, which means that we must check if the internal dynamics is stable in order to have the population control ef- fective (we cannot change its stability with the control input). In order to the system to operates in zero dynamics, i.e., to the state trajectory x remain in Mo, the initial system state x(0) must be in Mo, and fur- thermore, the output y must remain null. This dynamics is described in the normal form in Eq. (4) [7–9]. = = z ψ w ψ ˙ 0 ˙ (0, ). (4) For the use of the exact feedback linearization techniques, it is ne- cessary that the internal dynamics in (4) to be asymptotically stable. Otherwise, we will have a system that behaves linearly and stable in its observable part, that is, however, internally unstable. This comes im- mediately from the Slotine theorem [9]. 4. Input output linearization of system that describes the population dynamics of the Aedes aegypti Considering in the System (2)–(2b) the (measurable) output y as a linear combination of the states x2 and x3, this dynamic can be written as: ⎡ ⎣ ⎢ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ ⎥ = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ + + + − − ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ + ⎡ ⎣ ⎢ ⎢ − ⎤ ⎦ ⎥ ⎥ = + x x x x A x ϕx A x x A x A x βx μ x A x μ x x u t y ax bx ˙ ˙ ˙ ˙ 0 0 0 ( ), . 1 2 3 4 1 1 3 2 1 3 3 1 4 2 2 3 3 5 1 4 4 1 2 3 (5) According to [7]–[9] the dynamics presented in Eq. (5) has relative degree r= 2. Furthermore, a diffeomorphism that turns this in the normal form, has the expression: = z z ψ ψΨ ( , , , ),1 2 1 2 (6) where = =z y z y, ˙1 2 and ψj, j = 1, 2 are solutions of partial differential equation (PDE) ∇ψ.g= 0, j= 1, 2. Therefore, from (6): = + = + + − = + + z ax bx z aA x aA bβ x bμ x B x B x B x , ( ) , 1 2 3 2 3 1 4 2 3 3 1 1 2 2 3 3 (6a) where, = = + = −B aA B aA bβ B bμ, ( ) and .1 3 2 4 3 3 (7) Let us try the following solutions to the PDE ∇ψ.g= 0, j= 1, 2: = = ψ x ψ x , , 1 3 2 4 (8) which gives: = = + + +z z ψ ψ ax bx B x B x B x x xΨ ( , , , ) ( , , , ),1 2 1 2 2 3 1 1 2 2 3 3 3 4 (9) and B1, B2, B3 are defined in Eq. (7). From Eq. (9), Ψ is a global diffeomorphism if only if a≠ 0, ∀ x1, x2, x3, x4. Furthermore, the inverse diffeomorphism is, = = = − = − + + ⎛ ⎝ − ⎞ ⎠ x ψ x ψ x a z b a ψ x B aB z B z B bB a B ψ , , 1 , 1 1 , 3 1 4 2 2 1 1 1 2 1 1 1 2 1 2 3 1 (10) From Eqs. (9) and (10), the nonlinear dynamics (4) has a normal form: ⎡ ⎣⎢ ⎤ ⎦⎥ = ⎡ ⎣⎢ + ⎤ ⎦⎥ z z z z z ψ ψ z z ψ ψ u t ˙ ˙ Γ( , , , ) Y( , , , ) ( ) ,1 2 2 1 2 1 2 1 2 1 2 (11) C.A. dos Reis et al. Mathematical Biosciences 299 (2018) 51–57 53 ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ + ⎛ ⎝ + ⎞ ⎠ − + + ⎛ ⎝ − ⎞ ⎠ − ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ = ψ ψ β a z bβ a μ ψ A B aB z A B z A B bB a B ψ μ ψ y z ˙ ˙ , , 1 2 1 3 1 5 2 1 1 5 1 2 5 1 2 3 1 4 2 1 (11a) where: ⎜ ⎟ ⎜ ⎟ + = ⎛ ⎝ + − + ⎞ ⎠ + + + ⎛ ⎝ ⎛ ⎝ − ⎞ ⎠ + − + + − ⎞ ⎠ + − + − − + − + − z z ψ ψ z z ψ ψ u t a A B βB B B A B A B z A B A B B z B bB a B A B A B b a A B β B B ϕ μ B ψ A a B z a z aB bB ψ ψ a B z a z aB bB ψ u t Γ( , , , ) Y( , , , ) ( ) 1 ( ) ( ) ( ) 1 ( ) ( ) ( ( ) ) 1 ( ( ) ) ( ). 1 2 1 2 1 2 1 2 4 2 3 2 1 1 1 3 2 1 1 1 3 2 1 2 1 2 3 1 1 3 2 4 2 3 1 3 3 1 2 2 1 2 3 2 1 1 2 1 2 3 2 1 (11b) The nonlinear dynamics (11) is the observable part of the dynamics Eq. (4), i.e., the part that relates input-output. The linear dynamics Eq. (11a) is the internal dynamics of the System (5), i.e., it is the non- observable part of this dynamic. Note that the internal dynamics does not depend on the input u(t). According to [7]–[9], the knowledge of the internal dynamics is important because it allows the analysis of the zero dynamics, which can be used for asymptotic stability analysis of the nonlinear dynamics in Eq. (5). 5. Analysis of the equilibrium points and the asymptotic stability of the zero dynamics It is possible to determine the equilibrium points of the population dynamics in Eq. (5) by using standard methods, and there are two of them, namely: = =P P x x x x(0, 0, 0, 0) and ( , , , ),1 2 1 2 3 4 (12) where x x x x, , ,1 2 3 4 in (12) are calculated by: = − − = − = − = − x A βϕ A A μ A A u t x A A x x β μ x x A A A μ x ( ( )), , 3 and . 2 3 2 4 3 2 1 1 4 3 2 3 2 4 4 5 3 4 2 (13) It indicates that the population can stabilize in a specific point (that depends on the parameter values) or in the origin (extinction of the population). 5.1. Control law for stabilization of the origin P1 According to [7]–[9], the zero dynamics is obtained from 11a) using z1= z2=0 ∀t. With these conditions and from ((4), we have: ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ = ⎡ ⎣⎢ − ⎤ ⎦⎥ = ψ ψ C ψ C ψ μ ψ y ˙ ˙ , 0, 1 2 1 1 2 1 4 2 (14) Where = − ⎛ ⎝ + ⎞ ⎠ = ⎛ ⎝ − ⎞ ⎠ C bβ a μ C A B bB a B , . 1 3 2 5 1 2 3 (15) From (14) we note that P= (0, 0) is the only critical point of the zero dynamics. Moreover, its characteristic equation is given by: − − − =C λ μ λ( )( ) 0,1 4 (16) whose eigenvalues are: = − ⎛ ⎝ + ⎞ ⎠ = − λ bβ a μ λ μ , . 1 3 2 4 (17) The dynamics in Eq. (14) is asymptotically stable around the point P= (0, 0) if and only if > −b a μ β 3 , implying that the Slotine theorem can be applied under this conditions. In fact, the nonlinear control (18) transforms the open loop dynamics (5) globally asymptotically stable around P1 and is given by: = − ⎡ ⎣ ⎢ + + + + − + + + − ⎤ ⎦ ⎥u t a A aA A aA bβ A x aA bβ A bβμ x aA ϕ bμ x aA A x x v ( ) 1 ( ( ) ) (( ) ) ( ) , 1 3 1 3 4 3 1 4 4 3 2 3 3 3 2 3 1 3 (18) where = − − + + − −v x ak A x ak A bk β ak x bk bk μ x( ) ( ) ( )o o1 3 1 1 4 1 2 1 3 3 and k1, ko are conveniently chosen such that the polynomial p (λ)= λ2+ k1λ+ko has negatives roots [9]. We can note that p(λ)= λ2+ k1λ+ko has negative roots if and only if ko> 0 and: = = − = ∈ −∞ − ∪ +∞ > ∈ − < k k k k k k k k k k 1. 2 , if the discriminant Δ 4 0; 2. ( , 2 ) (2 , ) if the discriminant Δ 0; 3. ( 2 , 2 ) if the discriminant Δ 0. 1 0 1 2 0 1 0 0 1 0 0 (18-a) The proposed control law (with the proposed parameter values) turns the equilibrium point P1, which is naturally unstable, in globally stable point. 5.2. Control law for stabilization of the equilibrium point P2 In order to design a control law to stabilize the population in the equilibrium point P2= x x x x( , , , )1 2 3 4 , it is necessary to translate this point to the origin by a coordinate transformation ( =x x0 1 1 , =x x0 2 2 , =x x0 3 3 , =x x0 4 4 ), which means: = − = − = − = − y x x y x x y x x y x x , , , . 1 1 0 1 2 2 0 2 3 3 0 3 4 4 0 4 (19) The equations in (19) results in: = − + + ⎛ ⎝ − ⎞ ⎠ − = − − = − = − y B aB z B z B bB a B ψ x y a z b a ψ x y ψ x y ψ x 1 1 , 1 , , . o o o o 1 2 1 1 1 2 1 2 3 1 1 2 1 1 2 3 1 3 4 2 4 (20) The zero dynamics in this case can be expressed as: ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ = ⎡ ⎣⎢ − ⎤ ⎦⎥ = ψ ψ C ψ C ψ μ ψ y ˙ ˙ , 0. 1 2 1 1 2 1 4 2 (21) The asymptotic analysis of the dynamic stability is analogous to the C.A. dos Reis et al. Mathematical Biosciences 299 (2018) 51–57 54 analysis of the point P1. Eqs. (16) and (17) imply that (5) is asympto- tically stable around P2 if and only if > −b a μ β 3 . In addition, the control (18) makes the dynamic (4) asymptotically stable around P2. To determine the control that stabilizes the system (5) (according to Slotine theorem) firstly observe that from the translation (19), it is possible to write the dynamics (5) as: ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ + + + + + + + + + + − + − − + − ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ + ⎡ ⎣ ⎢ ⎢ ⎢ − + ⎤ ⎦ ⎥ ⎥ ⎥ = + + + y y y y A A x y ϕ A x y A y y A x ϕx A x x A y A y A x A x βy μ y βx μ x A y μ y A x μ x y x u t y ay by ax bx ˙ ˙ ˙ ˙ ( ) ( ) ( ) 0 0 0 ( ), . o o o o o o o o o o o o o o o 1 2 3 4 1 2 3 1 2 1 3 2 1 3 1 1 3 2 1 3 3 1 4 2 3 1 4 2 2 3 3 2 3 3 5 1 4 4 5 1 4 4 1 1 2 3 2 3 (22) Therefore, the nonlinear control is given in Eq. (23): = − − + + − − − − + + − − u t ak A y ak A bk β ak y bk bk μ y ak A x ak A k bβ ak x bk k bμ x ( ) ( ) ( ) ( ) ( ) . o o o o o o o 1 3 1 1 4 1 2 1 3 3 1 3 1 1 4 1 2 1 3 3 (23) Where k1 and ko are conveniently chosen such that p (λ)= λ2+ k1λ+ko has negative roots [9]. We can note that p (λ)= λ2+ k1λ+ko (with discriminant Δ= k12-4 ko) has negative roots if and only if ko> 0 and: • =k k2 ,1 0 if the discriminant is Δ=0; • ∈ −∞ − ∪ +∞k k k( , 2 ) (2 , )1 0 0 if the discriminant is Δ>0; • ∈ −k k k( 2 , 2 )1 0 0 if the discriminant is Δ<0. 6. Optimized control The control law u(t) in (18) depends of the values of the parameters a, b, k1 and k0 and there are infinite possibilities to choose values for these constants in order to have closed-loop stability. Thereore, we propose a methodology to choose the values to these constants such that the complete system has other desirable properties beyond stabi- lity, as optimality related to use the minimum quantity possible of aquatic control effort and to keep the mosquito population as low as possible. In this way, in this section we propose a multiobjective opti- mization model, presented in (24)–(28). aiming to determine a strategy that use a minimum quantity possible of aquatic control and to keep the mosquito population as low as possible: ∫ ∫⎡ ⎣⎢ = = + + + ⎤ ⎦⎥= = Minimize J udt J x x x x dt, ( ) t T t T 1 0 2 0 1 2 3 4 (24) Subject to System in Equation and Equation(5) (18), (25) >k 0,o (26) >b a μ β 3 (27) ≥ ≥ ≥ ≥ ≥ ∈u t x t x t x t x t t T( ) 0, ( ) 0, ( ) 0, ( ) 0, ( ) 0, [0, ].1 2 3 4 (28) where = = = =x x t x x t x x t x x t( ), ( ), ( ), ( )1 1 2 2 3 3 4 4 are the state variables of System (5) and =u u t a b k k( , , , , )0 1 is the control defined in (18). The functionals J1 and J2 proposed in (24) are related to the control applied in the time interval [0, T] and the total density of mosquitoes in this interval. Constraints (25) define the dynamics to be followed by the system and the control law. Constraints (26) and (27) are imposed to ensure that the proposed control stabilizes the System (5). The con- straints (28) define the non-negativity of the variables, guaranteeing the biological meaning of these variables. We define Sp={s=(a,b,k0,k1) such that k0>0, (b/a)> (-μ3/β) and k1 ∈ R}. Let SF={X=(x1(t),x2(t),x3(t), x4(t),u(t,s)) such that x1(t)> 0, x2(t)> 0, x3(t)> 0, x4(t)> 0, u(t,s)≥ 0, s∈SP, t∈[0,T]} be a set of feasible solutions to the Problem (24)–(28). Therefore, a feasible solution ∈ ⌣ X SF can be generated by taking some ∈ ⌣ s SP and solving the System (5) with u(t,s) defined as in (18). Aiming to solve the multiobjective optimization problem (24)–(28), we determine a pair (s*, X*), s*∈ Sp and X*∈ SF, which simultaneously minimizes the functionals J1 and J2 defined in (24) and are related respectively to the amount of control applied and the density of total population of mosquitoes in the predetermined time interval [0, T]. However, it is known that J1 and J2 are conflicting, i.e., the smaller the J1 value, the larger J2 will be (the less control, the more mosquito po- pulation grows, until it stabilizes at high levels). Therefore, the multi- objective optimization techniques can assist in the determination of Efficient Solutions for this problem [12]. We propose to use an adaptation of the ε–restricted method defined in Deb [12] in order to determine a efficient solution, which redefines the problem in a set of sub-problems in the form: (there is a different problem for each ε value). ∫⎡ ⎣⎢ = ⎤ ⎦⎥= Minimize J u t s dt( , ) t T 1 0 (29) Subject to + + + ≤ ≥x t x t x t x t t t( ) ( ) ( ) ( ) ɛ, ,fix1 2 3 4 (30) = ∈X x t x t x t x t u t s S( ( ), ( ), ( ), ( ), ( , )) .F1 2 3 4 (31) where ε is the required level for mosquito density from a defined time tfix (ε and tfix are fixed values). We propose the heuristic VNS-Variable neighborhood search [13] to determine an approximated solution to the Optimization Problem (29)–(31). Problem (29)–(31) can be summarized in: Let given the values for T, tfix and ε, we need to find the value D*=(s*, X*), where s*∈ Sp and X*∈ SF , such that the constraint (30) is satisfied and J1 is minimum, where: Sp={s=(a, b, k0, k1) such that k0>0, (b/a)> (-μ3/β) and k1∈ R}. SF={X=(x1(t), x2(t), x3(t), x4(t), u(t,s)) such that x1(t)> 0, x2(t)> 0, x3(t)> 0, x4(t)> 0, u(t,s)≥ 0, s ∈ SP, t ∈ [0,T]}. VNS is a stochastic algorithm in which a set of neighborhood structures Nk (k= 1,…, kmax) are defined a priori. Each iteration of the algorithm is composed of three steps: perturbation, local search and movement. At each iteration, a solution is disturbed in the current neighborhood Nk. i.e., a feasible solution D' is generated randomly. A local search procedure is applied in neighborhood Nk(D’), reaching the solution D”. The current solution is replaced by the new optimal local D” if and only if a better solution has been found, i.e., if J1(D”) < J1(D’). And then, the same search procedure is restarted from solution D” in the first neighborhood N1. Otherwise, if J1(D”) ≥ J1(D’), the algorithm moves to the next neighborhood, and attempts to im- prove the current solution. The proposed VNS Algorithm is presented in the following: Neighborhood structure considered: Let D’=(s’, X’), s’ ∈ Sp and X’ ∈ SF , such that s’=(a’, b’,k0’, k1’). A neighborhood N1(D’) is defined fixing value to b’,k0’, k1’ in s’ and by varying the value of a, i.e., N1(D’) = {(s, X) such that s= ( a, b’, k0’, k1’)∈Sp, a≠a’ and X∈ SF }. In the same way, we can define: N2(D’) = {(s, X) such that s= ( a’, b, k0’, k1’) ∈ Sp, b ≠ b’ and X ∈ SF C.A. dos Reis et al. Mathematical Biosciences 299 (2018) 51–57 55 }; N3(D’) = {(s, X) such that s= ( a’, b’, k0, k1’) ∈ Sp, k0 ≠ k0’ and X ∈ SF }; N4(D’) = {(s, X) such that s= ( a’, b’, k0’, k1) ∈ Sp, k1 ≠ k1’ and X ∈ SF }. VNS Algorithm Let kmax=4 and the neighborhood Nk(D), k=1, 2, 3, 4 as defined above, do: 1. Determine a random initial feasible solution Do=(so, Xo), so∈ Sp and Xo∈ SF; 2. K ←1 , D’ ← Do; 3. While k≤ kmax 3.1 Find the Best Neighbor D”∈Nk (D’); (perturbation and local search); 3.2 If ( J1(D”) < J1(D’) ); 3.2.1 D’ ←D”; k ←1; (movement) 3.2.2 Else 3.2.3 k ← k+ 1; 3.3 End-If; 4. End-While; 5. Do D*=D’. Best Neighbor Algorithm Let k and D’ be as previously defined, and lmax the maximum number of iterations, do: 1. Generate randomly D”∈ Nk (D’), l ←1 2. While l≤ lmax 2.1 Perturb D”,( generate randomly D” ∈ Nk (D’), (D” ≠ D”)); 2.2 If ( J1(D”) < J1(D”) ); 2.2.1 D” ← D”; 2.3 End-If; 3. l ← l + 1; 4. End-While; 5. Take back D”. The VNS algorithm determines a approximated solution to the Problem (29)–(31), determining an optimized control strategy u(t), in the time interval [0, T] that keeps the population density of mosquitoes at a required low level ε. We present in the following section a com- putational simulation for the methodology described in this work. Table 1 Parameter values. Source: ([11])*, ([14])⁎⁎, ([15])⁎⁎⁎ Parameters C* γ⁎⁎ R* β⁎⁎ ϕ * μ1⁎⁎ μ2⁎⁎⁎ μ3⁎⁎⁎ μ4* Value 3 0.0941 0.5 0.927 6.353 0.0583 0.0202 0.0319 0.06 Fig. 2. Simulation for the open loop situation (different initial conditions). 0 50 100 150 0 20 40 60 80 100 120 140 160 Time (day) u(t) Fig. 3. Optimized control strategy obtained applying the proposed VNS heuristic. C.A. dos Reis et al. Mathematical Biosciences 299 (2018) 51–57 56 7. Numerical simulation Aiming to implement and analyze the proposed methodology, some computational experiment was made. We use for these implementations the MATLAB 7.14.0.739 (R2012a) software and run on a micro- computer Dual Core i5-650 with 4GB memory and 400GB hard drive running Windows Operating System. Table 1 shows the parameters used in the simulations for a parti- cular case when the average local temperature is 25 °C, based in [11, 14] and [15]. All parameters are measured per day (d−1) except r. Fig. 2 shows the effect on the response of the dynamic of Eq. (5) without the control law (i.e. in open loop situation) for 150 days. We can observe in Fig. 2, that the chosen control using the Eq. (18) and without optimization process, is capable of stabilizing the system with different initial condition. However, the density in all segments of the population (aquatic phase, immature female, fertilized female and male) stabilize in the equilibrium point P2. Fig. 3 shows the optimized control strategy, obtained using the heuristic VNS discussed in Section 7. We used the data in Table 1, the total control application time T=150 days, ε=1 and tfix=70 days. The control was applied from the equilibrium point P2, for being a time in which the population starts the stabilization in higher levels (period where begins the critical phase of the proliferation of diseases trans- mitted by this vector). The determined values for the optimized para- meters of u(t) were a*= 0.0672, b*= 6.0197, k0*= 3.8677, k1*=−0.8401. Fig. 3 shows that we must apply the control with higher intensity in the first days, reducing it as the mosquitos populations reduces, as shown in Fig. 4, which shows the effect of the optimized control on the population dynamic of the mosquitoes. 8. Discussion The control of the mosquito population in aquatic phase is a common practice in Brazil. Healthcare agents search and remove mosquito breeding sites in homes, baud lands, businesses and industries throughout the year. Another common practice is the television and press advertisements, as well as educational campaigns of the popula- tion so that the breeding sites could be removed more efficiently. However, these practices are very expensive for healthcare agencies and often performed poorly. In this sense, the methodology presented here can be useful as a cost effective tool. The feedback linearization control technique was used in order to obtain a nonlinear control law that linearizes globally the mosquito population dynamics. The resulting control law has four free para- meters to choose by the designer. The heuristic VNS method was used to aid in the determination of optimized values for these parameters, obtaining an optimized control strategy to be applied in the aquatic phase of the mosquito population. This strategy consists in applying the control with higher intensity in the first days, afterwards reducing it as the populations reduces. This strategy reduces the cost of control as well as the mosquito population to desired levels in the defined time. The reduction of the Aedes aegypti population favors the reduction of the diseases transmitted by this vector. 9. Conclusion This work analyzed a mathematical model describing the popula- tion dynamics of the Aedes aegypti, and designed a generic control law for the mosquito population through exact linearization techniques. Based in this control law, we proposed a heuristic VNS algorithm to determine optimized controls strategies for the aquatic phase, which produced the minimum possible cost of control, making the process cheaper and lowering the mosquito population, which favors the re- duction of the diseases transmitted by this mosquito, such as dengue, chikungunya, zika and yellow fever. Other advantages of this control technique is the reduction of environmental impact, since reducing the amount of mosquitoes also reduces the need for chemical control, as for example the insecticide. The results of the computer simulations show that the presented methodology has great potential as a tool for the decision makers in the fight against the Aedes aegypti. Acknowledgements We wish to thank FAPESP (Grant No. 2009/15098-0 and 2014/ 01604-0), FUNDUNESP/PROPE/UNESP, CNPq (302454/2016-0), CAPES and PROPG UNESP for their financial support. References [1] WHO - World Health Organization (2017). Dengue and severe dengue. Available at: http://www.who.int/mediacentre/factsheets/fs117/en/. Accessed February 23, 2017. [2] N. Vasilakis, J. Cardosa, K.A. Hanley, E.C. Holmes, S.C. Weaver, Fever from the forest: prospects for the continued emergence of sylvatic dengue virus and its im- pact on public health, Nat. Rev. Microbiol. 9 (2011) 532–541. [3] M.S. Mustafa, V. Rasotgi, S. Jain, V. Gupta, Discovery of fifth serotype of dengue virus (DENV-5): A new public health dilemma in dengue control, Med. J. Armed Forces India 71 (2015) 67–70. [4] J.F.W. Chan, G.K.Y. Choi, C.C.Y. Yip, V.C.C. Cheng, K.Y. Yueng, Zika fever and congenital Zika syndrome: An unexpected emerging arboviral disease, J. Infect. 72 (2016) 507–524 http://dx.doi.org/10.1016/j.jinf.2016.02.011. [5] H.O. Florentino, B.F. Bannwart, D.R. Cantane, F.L.P. Santos, Multiobjective genetic algorithm applied to dengue control, Math. Biosci. 258 (2014) 77–84. [6] L. Esteva, H.M. Yang, Mathematical model to assess the control of Aedes aegypti mosquitoes by the sterile insect technique, Math. Biosci. 198 (2005) 132–147. [7] A. Isidori, Nonlinear Control Systems, Third ed., Springer-Verlag, Roma, 1995. [8] H.K. Khalil, Nonlinear Systems, Prentice Hall, New Jersey, 2002. [9] J. Slotine, W. Li, Applied Nonlinear Control, Prentice Hall, New Jersey, 1991. [10] C. Liqun, L Yanzhu, A modified exact linearization control for chaotic oscillators, Nonlinear Dyn. 20 (1999) 309–317. [11] R.C.A. Thomé, H.M. Yang, L. Esteva, Optimal control of Aedes aegypti mosquitoes by the sterile insect technique and insecticide, Math. Biosci. 223 (2010) 12–23. [12] K. Deb, Multi-Objective Optimization Using Evolutionary Algorithms, John Wiley & Sons, 2001. [13] P. Hansen, N. Mladenović, J.A.M Perez, Variable neighbourhood search: methods and applications, 4OR. 6 (2008) 319–360, http://dx.doi.org/10.1007/s10288-008- 0089-1 2008. [14] H.M. Yang, M.L.G. Macoris, K.C. Galvani, M.T.M. Andrighetti, D.M.V. Wanderley, Assessing the effects of temperature on the population of Aedes aegypti, vector of dengue, Epidemiol. Infect. 137 (2009) 1188–1202. [15] M.A.S. Koiller, M.O. Souza, C. Codeço, A. Iggidr, G. Sallet, Aedes, Wolbachia and Dengue, (2014) INRIA Research Report RR8462. 0 50 100 150 0 0.5 1 1.5 2 2.5 3 3.5 4 x4(t) x2(t) x3(t) x1(t) Fig. 4. Population dynamic of the mosquitoes in the aquatic phase (x1(t)), immature females (x2(t)), fertilized females (x3(t)), males (x4(t)) under the application of the control. C.A. dos Reis et al. Mathematical Biosciences 299 (2018) 51–57 57 http://www.who.int/mediacentre/factsheets/fs117/en/ http://refhub.elsevier.com/S0025-5564(17)30172-4/sbref0001 http://refhub.elsevier.com/S0025-5564(17)30172-4/sbref0001 http://refhub.elsevier.com/S0025-5564(17)30172-4/sbref0001 http://refhub.elsevier.com/S0025-5564(17)30172-4/sbref0002 http://refhub.elsevier.com/S0025-5564(17)30172-4/sbref0002 http://refhub.elsevier.com/S0025-5564(17)30172-4/sbref0002 http://dx.doi.org//10.1016/j.jinf.2016.02.011 http://refhub.elsevier.com/S0025-5564(17)30172-4/sbref0004 http://refhub.elsevier.com/S0025-5564(17)30172-4/sbref0004 http://refhub.elsevier.com/S0025-5564(17)30172-4/sbref0005 http://refhub.elsevier.com/S0025-5564(17)30172-4/sbref0005 http://refhub.elsevier.com/S0025-5564(17)30172-4/sbref0006 http://refhub.elsevier.com/S0025-5564(17)30172-4/sbref0007 http://refhub.elsevier.com/S0025-5564(17)30172-4/sbref0008 http://refhub.elsevier.com/S0025-5564(17)30172-4/sbref0009 http://refhub.elsevier.com/S0025-5564(17)30172-4/sbref0009 http://refhub.elsevier.com/S0025-5564(17)30172-4/sbref0010 http://refhub.elsevier.com/S0025-5564(17)30172-4/sbref0010 http://refhub.elsevier.com/S0025-5564(17)30172-4/sbref0011 http://refhub.elsevier.com/S0025-5564(17)30172-4/sbref0011 http://dx.doi.org/10.1007/s10288-008-0089-1 http://dx.doi.org/10.1007/s10288-008-0089-1 http://refhub.elsevier.com/S0025-5564(17)30172-4/sbref0013 http://refhub.elsevier.com/S0025-5564(17)30172-4/sbref0013 http://refhub.elsevier.com/S0025-5564(17)30172-4/sbref0013 http://refhub.elsevier.com/S0025-5564(17)30172-4/sbref0014 http://refhub.elsevier.com/S0025-5564(17)30172-4/sbref0014 An approach of the exact linearization techniques to analysis of population dynamics of the mosquito Aedes aegypti Introduction Population dynamics of the mosquito Aedes aegypti Mathematical modeling Exact linearization techniques and zero dynamics Input output linearization of system that describes the population dynamics of the Aedes aegypti Analysis of the equilibrium points and the asymptotic stability of the zero dynamics Control law for stabilization of the origin P1 Control law for stabilization of the equilibrium point P2 Optimized control Numerical simulation Discussion Conclusion Acknowledgements References