Multiobjective approach to optimal control for a dengue transmission model

During the last decades, the global prevalence of dengue progressed dramatically. It is a disease which is now endemic in more than one hundred countries of Africa, America, Asia and the Western Pacific. This study addresses a mathematical model for the dengue disease transmission and finding the most effective ways of controlling the disease. The model is described by a system of ordinary differential equations representing human and vector dynamics. Multiobjective optimization is applied to find the optimal control strategies, considering the simultaneous minimization of infected humans and costs due to insecticide application. The obtained results show that multiobjective optimization is an effective tool for finding the optimal control. The set of trade-off solutions encompasses a whole range of optimal scenarios, providing valuable information about the dynamics of infection transmissions. The results are discussed for different values of model parameters.


Introduction
Dengue is a vector-borne disease transmitted from an infected human to a female Aedes mosquito by a bite. The mosquito, which needs regular meals of blood to feed their eggs, bites a potentially healthy human and transmits the disease, turning it into a cycle. There are four distinct, but closely related, viruses that cause dengue. The four serotypes, named DEN-1 to DEN-4, belong to the Flavivirus family, but they are antigenically distinct. Recovery from infection by one serotype provides lifelong immunity against that serotype but provides only partial and transient protection against subsequent infection by the other three viruses. There are strong evidences that a sequential infection increases the risk of developing dengue hemorrhagic fever.
The spread of dengue is attributed to the geographic expansion of the mosquitoes responsible for the disease: Aedes aegypti and Aedes albopictus. The Aedes aegypti mosquito is a tropical and subtropical species widely distributed around the world, mostly between latitudes 35 o N and 35 o S. In urban areas, Aedes mosquitoes breed on water collections in artificial containers such as cans, plastic cups, used tires, broken bottles and flower pots. Due to its high interaction with humans and its urban behavior, the Aedes aegypti mosquito is considered the major responsible for the dengue transmission. The life cycle of a mosquito has four distinct stages: egg, larva, pupa and adult. In the case of Aedes aegypti, the first three stages take place in, or near, the water, whereas the air is the medium for the adult stage [13].
It is very difficult to control or eliminate Aedes aegypti mosquitoes due to their resiliency, fast adaptation to changes in the environment and their ability to rapidly bounce back to initial numbers after disturbances resulting from natural phenomena (e.g., droughts) or human interventions (e.g., control measures). Primary prevention of dengue resides mainly in mosquito control. There are two primary methods: larval control and adult mosquito control, depending on the intended target. Larvicide treatment is done through long-lasting chemical in order to kill larvae and preferably have World Health Organization clearance for use in drinking water [3]. Adulticides is the most common measure. Its application can have a powerful impact on the abundance of adult mosquito vector. However, the efficacy is often constrained by the difficulty in achieving sufficiently high coverage of resting surfaces [4].
The present study addresses a problem of finding the most effective ways of controlling the dengue disease. To this end, a mathematical model for the dengue transmission, including a control variable presented in [18], is adopted. The main contributions are: (i) adapting multiobjective optimization methodologies for finding the optimal control in a mathematical model for the dengue transmission, (ii) analysis of the search ability of different methods on the resulting optimization problem, (iii) discussion of results for the optimal control problem, emphasizing advantages of the proposed approach when compared with results available in the literature. The aim of the work is to promote multiobjective optimization in epidemiological research community and practice.
The remainder of this paper is organized as follows. Section 2 describes a mathematical model for the dengue disease transmission. The problem of finding the optimal control is formulated in Section 3, including traditional and proposed approaches. Section 4 presents some general concepts in multiobjective optimization and methods adopted for solving the problem. Section 5 presents and discusses our results, obtained by numerical simulations, which includes the comparison of different multiobjective optimization methods, control strategies and dengue dynamics that correspond to different optimal scenarios. The study is ended with Section 6 of conclusions and some directions of future work.

The mathematical model
This section introduces a mathematical model for the dengue disease transmission. The model is based on a system of ordinary differential equations, and includes the real data of a dengue disease outbreak that occurred in the Cape Verde archipelago in 2009 [12,17].
The model consists of eight mutually-exclusive compartments representing the human and vector dynamics. It also includes a control parameter, an adulticide spray, as a measure to fight the disease.
The notation used in the mathematical model includes four epidemiological states for humans: S h (t) -number of susceptible individuals at time t; E h (t) -number of exposed individuals at time t; I h (t) -number of infected individuals at time t; R h (t) -number of resistant individuals at time t.
It is assumed that the total human population (N h ) is constant, so There are also four other state variables related to the mosquitoes (disease vectors): A m (t) -number of vectors in aquatic phase at time t; S m (t) -number of susceptible vectors at time t; E m (t) -number of exposed vectors at time t; I m (t) -number of infected vectors at time t.
Similarly, it is assumed that the total adult mosquito population (N m ) is constant, which means that N m = S m (t) + E m (t) + I m (t). The model includes a control variable, which represents the amount of insecticide that is continuously applied during a considered period, as a measure to fight the disease: The control variable c(t) is an adimensional value that is considered in relative terms, varying from 0 to 1.
In the following, for the sake of simplicity, the independent variable t is omitted when writing the dependent variables (for instance, S h is used instead of S h (t)).
The parameters necessary to completely describe the model are presented in Table 1. This set of parameters includes the real data related to the outbreak of dengue disease occurred in Cape Verde in 2009 [17].
Furthermore, in order to obtain a numerically stable problem, all the state variables are normalized as follows: Thus, the dengue epidemic is modeled by the following nonlinear time-varying state equations [18]: Since any mathematical model is an abstraction of a complex natural system, additional assumptions are made to make the model mathematically treatable. This is also the case for the above epidemiological model, which comprises the following assumptions: • the total human population (N h ) is constant; • there is no immigration of infected individuals into the human population; • the population is homogeneous, which means that every individual of a compartment is homogeneously mixed with other individuals; • the coefficient of transmission of the disease is fixed and does not vary seasonally; • both human and mosquitoes are assumed to be born susceptible, i.e., there is no natural protection; • there is no resistant phase for mosquitoes, due to their short lifetime.

Problem formulation
In this section, the problem of finding the optimal control in the mathematical model for the dengue transmission is formulated. The optimal control represents the most effective way of controlling the disease that can be adopted by authorities in response to its outbreak. First, a traditional approach based on optimal control theory is presented. Then, the proposed approach, based on multiobjective optimization, is introduced.
3.1. Optimal control theory. A problem of finding a control law for a given system is commonly formulated and solved using optimal control theory. Our aim is to find the optimal value c * of the control c, such that the associated state Consider the state system of ordinary differential equations (1) and the set of admissible control functions given by: The objective functional is defined by [16]: where γ D and γ S are positive constants representing the costs weights of infected individuals and spraying campaigns, respectively. The objective functional given by (3) is a function of state and control variables, representing cumulative costs due to infected population and prevention measures. The optimal control problem consists in determining , associated to an admissible control c * (·) ∈ Ω on the time interval [0, T ], satisfying (1), the initial conditions (2), and minimizing the cost functional (3), i.e., The optimal control can be derived using Pontryagin's maximum principle [15].

3.2.
Multiobjective approach. An approach based on optimal control theory allows to obtain a single optimal solution. The obtained solution represents some decision maker's perspective on controlling the disease. This is reflected by the constants γ D and γ S that must be provided in advance. Though the choice of proper values for the constants is not an easy task for all the cases. A single optimal solution provides limited information for the decision maker, while leaving a large range of alternatives unexplored. To overcome these drawbacks, the present study decomposes the cost functional (3) into two separate objectives and seeks the optimal control, minimizing simultaneously the costs due to the infected population and the costs associated with insecticide. The resulting optimization problem is defined as follows: (1) and (2) where T is a given period of time, f 1 and f 2 represent the cost incurred in the form of infected population and the cost of applying insecticide for the period T , respectively.

Multiobjective optimization
This section presents general concepts in multiobjective optimization and discusses some popular methods for multiobjective optimization based on scalarization. The outline of an algorithm for finding optimal solutions is presented. For a method that gives a direct solution to the full Pareto set, via solutions of a certain PDE, we refer the reader to [6,7,8]. For other frameworks see [19] and references therein.
4.1. General definitions. Without loss of generality, a multiobjective optimization problem (MOP) with m objectives and n decision variables can be formulated as follows: where x is the decision vector, Ω ⊆ R n is the feasible decision space, and f (x) is the objective vector defined in the attainable objective space Θ ⊆ R m . In multiobjective optimization, the Pareto dominance relation is usually used to define the concepts of optimality. For two solutions x and y from Ω, a solution x is said to dominate a solution y (denoted by x ≺ y) if for all i ∈ {1, . . . , m} f i (x) ≤ f i (y) and there exists at least one objective such that f j (x) < f j (y).
A solution x * ∈ Ω is Pareto optimal if and only if: In the presence of multiple conflicting objectives, there is a set of optimal solutions, known as the Pareto optimal set. For MOP (5), the Pareto optimal set (or Pareto set for short) is defined as: For MOP (5) and the Pareto set PS, the Pareto optimal front (or Pareto front for short) is defined as: Since it is often not possible to obtain the whole Pareto set, solving (5) is usually understood as approximating the Pareto set by obtaining a set of solutions that are as close as possible to the Pareto set and as diverse as possible.

Scalarization methods.
In the following, four scalarization methods for multiobjective optimization, able to deal with convex and nonconvex Pareto fronts, are discussed.
4.2.1. The ǫ-Constraint method. In the ǫ-constraint method [5], one of the objective functions is selected to be minimized, whereas all the other functions are converted into constrains by setting an upper bound to each of them. It can be defined as: minimize: where the lth objective is minimized and the parameter ǫ i represents an upper bound of the value of f i .

4.2.2.
Chebyshev's method. Chebyshev's method [1] minimizes the weighted metric L p , with p = ∞, which measures the distance from any solution to some reference point. It can be formulated as: (6) minimize: where z * = (z * 1 , . . . , z * m ) T is a reference point and w = (w 1 , . . . , w m ) T is a weight vector such that m j=1 w j = 1.

4.2.3.
The goal attainment method. The problem shown in (6) can be reformulated as [11]: T is a reference point, w = (w 1 , . . . , w m ) T is a weight vector ( m j=1 w j = 1) and x ∈ Ω, α ∈ R + are variables. This method is often referred to as the goal attainment method [11] or the Pascoletti-Serafini scalarization [14].

4.2.4.
The normal constraint method. The normal constraint (NC) method [9,10] minimizes one of the objective functions and uses an inequality constraint reduction of the feasible objective space. It can be formulated as: minimize: To cope with differently scaled objectives, the NC method normalizes the objective vector as: (7) f where z ideal i and z nadir i are the ith components of the ideal and nadir points, respectively.

4.
3. An algorithm for the Pareto set approximation. The above discussed methods convert a multiobjective optimization problem into a single-objective problem depending on some parameters. Solving the corresponding problem for different parameter settings allows to approximate multiple Pareto optimal solutions. The outline of an approach to approximate the Pareto set of problem (4) is shown in Algorithm 1. for initial point x (0) , find x * that minimizes f β (x, β);

5:
A ← A ∪ {x * }; 6: x (0) ← x * ; 7: end for 8: output: A; In line 1 of Algorithm 1, an initial point is generated as the null vector, x (0) = 0, and a set of parameters, B, used for scalarization, is initialized. For the ǫ-constraint method, β ∈ B corresponds to the value of ǫ. In the case of the Chebyshev and goal attainment methods, β is a tuple of the form {w, z * }, whereas β is a tuple of the form {w, z ideal , z nadir } for the normal constraint method. In lines 3-7, for each corresponding scalarizing function, f β (x, β), a minimizer, x * , is found using a single-objective optimizer and added to the approximation set, A. The minimizer of the previous scalarizing function becomes an initial point for solving the subsequent problem, to ensure efficiency of the approach.

Numerical experiments
This section provides results of the numerical experiments conducted for the dengue transmission model. The comparison of different scalarization methods is performed. The analysis of obtained solutions and the dengue dynamics are presented.

Experimental setup.
The fourth-order Runge-Kutta method is used to discretize the control and state variables of the system (1). The period [0, 84] of 84 days is discretized using the equally spaced time intervals of 0.25 (6 hours). This results in an optimization problem with x ∈ [0, 1] 337 , where x denotes the discrete control in (1). The integrals defining the objective functionals in (4) are calculated using the trapezoidal rule.
For approximating the Pareto set of (4), scalarization is performed by defining and solving 100 scalarizing problems corresponding to the above discussed methods. As a single-objective optimizer (line 4 in Algorithm 1), the MATLAB R function  Table 2. Hypervolume values for trade-off solutions obtained using different scalarization methods (the higher the better). The best value is marked bold.
fmincon with a sequential quadratic programming algorithm is used, setting the maximum number of function evaluations to 20, 000.

Different scalarization methods.
For comparison, the outcomes obtained by different scalarization methods are assessed using the hypervolume [21]. The hypervolume can be defined as the Lebesgue measure, Λ, of the union of hypercuboids in the objective space: where A = {a 1 , . . . , a |A| } denotes a set of nondominated solutions and r = (r 1 , . . . , r m ) T is a reference point. The hypervolume calculates a portion of the objective space dominated by A. It can measure both convergence and diversity. The higher the value of HV , the better the quality of A. For calculating the hypervolume, the objectives are normalized using (7) and [1,1] is used as a reference point. Table 2 presents the results for the four scalarization methods with respect to the hypervolume. The best result is obtained by the normal constraint method followed by the ǫ-constraint method. Both methods minimize f 1 , when f 2 is used to define the constraint for the objective space reduction. Since the NC method uses evenly distributed points on the ideal plane, the final points on the Pareto front approximation are likely to be more evenly distributed than using the usual ǫ-constraint method. This allows to achieve the highest hypervolume. The Chebyshev method provides a slightly worse result than the goal attainment method. This can be because the scalarizing function used by the Chebyshev method is nondifferentiable, which usually poses additional difficulties for optimization. Figure 1 shows the Pareto front approximations obtained by the methods. The plots confirm the previous observations based on the hypervolume values. The Chebyshev and goal attainment methods provide visually similar results. The ǫconstraint method attempts to find evenly distributed points with respect to the uniform division of f 2 , whereas the NC method seeks a uniform distribution of points, according to points on the hyperplane passing through the corner points of the Pareto front. Since the NC method gives clearly better results for (4), solutions obtained by this method will be employed for the analysis of the dengue transmission in what follows. Figures 2 and 3 present the trade-off solutions obtained for different values of transition probabilities β hm and β mh , respectively. As it was seen during discussion of different scalarization methods, the trade-off curve for (4) exhibits a smooth relation between the insecticide cost, f 2 , and the infected human population, f 1 , in the range 0 ≤ f 2 ≤ 5. This suggests that implementing    optimal control strategies in this range can produce the most significant reduction in infected individuals, having the highest ratio between the desirable effect and the cost. However, starting from some further point, reducing the number of infected humans can be achieved through exponential increase in spendings for insecticide. Even a small decrease corresponds to a high increase in expenses for insecticide. Scenarios represented by this part of the trade-off curve can be unacceptable from    To provide a better visualization of the most interesting parts of the trade-off curves in Figures 2 and 3, the plots are separated into two parts and presented in different figures, according to the range of f 2 . From Figures 2a and 3a, it can be seen that the higher the corresponding transition probability, the larger number of infected individuals. Though this difference reduces as more control is applied, vanishing for the extreme scenario with c(·) ≡ 1. Varying β hm and β mh produces apparently similar effects on the shape of the Pareto front. From Figures 2a and 3a, it can be observed that for higher values of β hm and β mh the part of the trade-off curve in 0 ≤ f 2 ≤ 5 becomes increasingly nonconvex. Solutions in this region can be unattainable for methods that face difficulties in dealing with nonconvexities in the Pareto front. Figure 4 illustrates dengue epidemics, corresponding to extreme scenarios (c(·) ≡ 0 and c(·) ≡ 1). It can be seen that for higher values of transition probabilities the peak in infected population occurs early, having larger values (Figures 4a and 4b).

Dengue dynamics.
On the other hand, applying the maximum control, c(·) ≡ 1, the dynamics of i h (t) remain unchanged for different values of β hm and β mh (Figures 4c and 4d).
To discuss intermediate scenarios, the concept of the knee solution of the Pareto front is adopted. According to [2], it can be defined as follows. Given a Pareto front with the normalized objectives, a boundary line L(p 1 , p 2 ) is constructed though two extreme points p 1 and p 2 . For any point z on the boundary line L(p 1 , p 2 ), a point p z on the Pareto front along the normal n of the boundary line is identified. The Pareto optimal point p z * with the maximum distance from its corresponding boundary point z * along the normal direction is defined as the knee point. This concept is illustrated in Figure 6. The knee solution of the Pareto front is attractive, since it is often considered as the optima in objective trade-offs. Figure 5 shows the dynamics of the state and control variables for different values of the transition probabilities, corresponding to the knee solution. It can be seen that more controls are required for higher values of β hm and β mh during the entire period of study, which is consistent with expectations. However, a higher value of transition probability does not necessarily mean a higher value of infected human population for the knee solution. This is because the change in the shape of the Pareto front caused by varying a value of transition probability.
From the Karush-Kuhn-Tucker condition, it can be induced that under smoothness assumptions the Pareto set of a continuous multiobjective optimization problem defines a piecewise continuous (m − 1)-dimensional manifold. The Pareto set of a continuous biobjective optimization problem is a piecewise continuous curve in R n [20]. The optimal control for (4) can be defined as a surface in i h (t) × t. Solving (4) gives a discrete representation of this surface. Thus, Figure 7 shows a discrete representation of a surface defining the Pareto set for (4). For a given value of f 1 , slicing the surface gives the optimal control trajectory. From the figure, it can bee seen how this dynamic changes from c(·) ≡ 0 to c(·) ≡ 1. A peak is observed in an early period of T , which grows in accordance with the decrease in the number of infected humans. Similarly, Figure 8 shows a discrete representation of a surface defining dynamics of infected humans across the Pareto set. A negative correlation between the optimal control and infected population can be observed. As the amount of control decreases, the number of infected humans increases. For higher values of the control, the peak in i h (t) is smaller and occurs later. When the control is reduced, the peak of disease outbreak grows up taking place earlier.     From the above discussion, it can be seen that the results for the optimal control problem, presented in Section 3, that are obtained using a multiobjective optimization approach provide comprehensive insights about optimal strategies for dealing with the dengue epidemic and dynamics resulting from implementing those strategies. The ability of trade-off solutions to reflect the underlying nature of the problem constitutes a major advantage of multiobjective optimization, motivating its practical use in the process of planning intervention measures by health authorities.

Conclusions
Due to difficulties in treating the dengue disease, controlling and preventing its outbreaks is essential for keeping people healthy, especially in regions where the threat of dengue is high. This study discussed a mathematical model for the dengue disease transmission from the optimal control point of view. Multiobjective optimization approach is suggested for finding the optimal control strategies to deal with an outbreak of the dengue epidemic. A biobjective optimization problem is formulated, involving minimization of expenses due to the infected population and costs of applying insecticide. The approach avoids the use of a priori information provided by the decision maker in the form of weight coefficients and allows to reflect the intrinsic nature of the problem.
The problem is numerically solved by discretizing the control variable and using scalarization methods for approximating multiple Pareto optimal solutions. The obtained results reveal different perspectives on applying insecticide: a low number of infected humans can be achieved spending larger financial resources, whereas low spendings for prevention campaigns result in significant portions of the population affected by the disease. Different trade-offs between the objectives are represented by the obtained solutions. Varying transmission probabilities introduces changes in relation between the objectives and their trade-offs. Because the problem is convex, the Pareto set defining the optimal control remains unchanged for different values of transmission probabilities. The results of the study suggest advantages of multiobjective optimization for finding the optimal control. Once the Pareto set is approximated, the final decision on the control strategy can be made taking into consideration available financial resources and goals of public health care.
As future work, it is intended to extend the model for the dengue transmission including different control variables that represent distinct measures for controlling the disease. Investigating impacts of different values of model parameters is also the subject of future work.