Multiobjective Optimal Control of HIV Dynamics

Various aspects of the interaction of HIV with the human immune system can be modeled by a system of ordinary differential equations. This model is utilized, and a multiobjective optimal control problem MOOCP is proposed to maximize the CD4 T cells population and minimize both the viral load and drug costs. The weighted sum method is used, and continuous Pareto optimal solutions are derived by solving the corresponding optimality system. Moreover, a model predictive control MPC strategy is applied, with the final goal of implementing Pareto optimal structured treatment interruptions STI protocol. In particular, by using a fuzzy approach, the MOOCP is converted to a single-objective optimization problem to derive a Pareto optimal solution which among other Pareto optimal solutions has the best satisfaction performance. Then, by using an embedding method, the problem is transferred into a modified problem in an appropriate space in which the existence of solution is guaranteed by compactness of the space. The metamorphosed problem is approximated by a linear programming LP model, and a piecewise constant solution which shows the desired combinations of reverse transcriptase inhibitor RTI and protease inhibitor PI drug efficacies is achieved.


Introduction
Human Immunodeficiency Virus HIV infects CD4 T-cells, which are an important part of the human immune system, and other target cells.The infected cells produce a large number of viruses.Medical treatments for HIV have greatly improved during the last two decades.Highly active antiretroviral therapy HAART allows for the effective suppression of HIV-infected individuals and prolongs the time before the onset of Acquired Immune Deficiency Syndrome AIDS for years or even decades and increases life expectancy and quality for the patient but eradication of HIV infection does not seem possible with currently available antiretroviral drugs.This is due primarily to the establishment of a pool of latently infected T-cells and sites within the body where drugs may not achieve effective levels 1-3 .HAART contains two major types of anti-HIV drugs, reverse transcriptase inhibitors RTI and protease inhibitors PI .Reverse transcriptase inhibitors prevent HIV from infecting Mathematical Problems in Engineering cells by blocking the integration of the HIV viral code into the host cell genome.Protease inhibitors prevent infected cells from replication of infectious virus particles and can reduce and maintain viral load below the limit of detection in many patients.Moreover, treatment with either type of drug can also increase the CD4 T-cell count that are target cells for HIV.
Many of the host-pathogen interaction mechanisms during HIV infection and progression to AIDS are still unknown.Mathematical modeling of HIV infection is of interest to the medical community as no adequate animal models exist to test efficacy of drug regimes.These models can test different assumptions and provide new insights into questions that are difficult to answer by clinical or experimental studies.A number of mathematical models have been formulated to describe various aspects of the interaction of HIV with healthy cells 4 .The basic model of HIV infection is presented by Perelson et al. 5 that contains three state variables: healthy CD4 T-cells, infected CD4 T-cells, and concentration of free virus.Another model is presented in 6 that although maintaining a simple structure, it offers important theoretical insights into immune control of the virus based on treatment strategies.Furthermore, this modified model is developed to describe the natural evolution of HIV infection, as qualitatively described in several clinical studies 7 .
The problem of designing drug administration in HIV infected patients using mathematical models can be considered as multi-objective optimal control problems.For example, these objectives may include maximizing the level of healthy CD4 T-cells and minimizing the cost of treatment 8-15 , maximizing immune response and minimizing both the cost of treatment and viral load 16-21 , maximizing both the level of healthy CD4 Tcells and immune response and minimizing the cost of treatment 22 , maximizing the level of healthy CD4 T-cells while minimizing both the side effects and drug resistance 23 , and maximizing survival time of patient subject to drug cost 24 .
When a multi-objective problem is treated, each objective conflicts with one another and, unlike a single-objective optimization, the solution to this problem is not a single point, but a family of solutions known as the Pareto-optimal set.Among these solutions, the designer should find the best compromise taking into proper account the attributes and characteristics of the handled problem.
There exists a wide variety of methods that can be used to compute Pareto optimal points.A widely used technique consists of reducing the multi-objective problem to a singleobjective one by means of "scalarization" procedure.The weighted sum WS method is a commonly used scalarization technique which consists of assigning each objective function a weight coefficient and then optimizing the function obtained by summing up all the objective functions scaled by their weight coefficients.Nevertheless, it has intrinsic drawbacks.The WS method is highly scale dependent, an equal distribution of weights does not necessarily lead to an even spread along the Pareto front, and that points in a nonconvex part of the Pareto front cannot be obtained 25 .Recent scalar multiple objective optimization techniques such as Normal Boundary Intersection NBI 26 and Normalized Normal Constraint NNC 27 have been found to mitigate the disadvantages of the WS method.Recently, these methods have been successfully combined with direct optimal control approaches for the efficient solution of multi-objective optimal control problems.For example, in 28 , a successful application of NBI and NNC for the multiple objective optimal control of bio chemical processes has been reported, and in 29 several scalarization techniques for multi-objective optimization, for example, WS, NNC, and NBI have been integrated with fast deterministic direct optimal control approaches.
In the considered control approaches, the amount of medications can be either continuous or on-off type.This treatment is also known as structured treatment interruption STI .STI has received considerable attentions as it might reduce the risk of HIV mutating to strains which are resistant to current medication regimens.STI approach might also reduce possible long-term toxicity of the drugs.A concise summary of clinical STI studies, including protocols and results, is presented in 35 .In this paper, we consider a mathematical model of HIV dynamics that includes the effect of antiretroviral therapy, and we perform analysis of optimal control regarding maximizing the CD4 T-cell counts and minimizing both the viral load and cost of drugs.
The paper is organized as follows: in Section 2, the underlying HIV mathematical model is presented.Our formulation of the MOOCP is described in Section 3. In Section 4, the weighted sum method is used, and the MOOCP is converted to a single-objective optimal control problem.First, we assume that the cost of the drug regime varies quadratically with the amount of drug administered, and then we characterize the continuous Pareto optimal solutions using Pontryagin's Maximum Principle.Secondly, it is shown that the presence of linear cost functionals leads to STI-type Pareto optimal solutions, and an MPC-based technique is proposed to find this kind of solutions.In Section 5, a fuzzy approach is utilized and the MOOCP is converted to a single-objective optimization problem, with the final goal of finding a Pareto optimal solution which has the best satisfaction performance among other Pareto optimal solutions.Moreover, the converted problem is approximated by an LP model applying a measure-theoretical approach, and a piecewise constant solution is achieved.Some numerical experiments are provided in Section 6.The last section is the conclusion.

Presentation of a Working Model
In this paper, the pathogenesis of HIV is modeled with a system of ordinary differential equations ODEs described in 7 .This model can be viewed as an extension of basic HIV Model of Perelson et al. 5 Most of the terms in the model have straightforward interpretations as following.
The first equation represents the dynamics of the concentration of healthy CD4 T-cells x .The healthy CD4 T-cells are produced from a source, such as the thymus, at a constant rate s, and die at a rate dx.The cells are infected by the virus at a rate rxv.The second equation describes the dynamics of the concentration of infected CD4 T-cells y .The infected CD4 T-cells result from the infection of healthy CD4 T-cells and die at a rate ay and are killed by The effect of PI drugs is modeled by reducing the proliferation rate of viruses from infected cells, while the effect of RTI drugs is modeled by reducing the infection rate, and in this way, blocking the infection of CD4 T-cells by free virus.Hence, in this model the RTI drugs have an effect on virulence because their main role is halting cellular infection and preventing virus production by reducing the production rate from infected CD4 T cells.The model has several parameters that must be assigned for numerical simulations.The descriptions, numerical values, and units of the parameters are summarized in Table 1.These descriptions and values are taken from 7 .We note that 2.1 -2.6 with these parameters model dynamics of fast progressive patients FPP .

Multiobjective Optimal Control Formulation
A problem arising from the use of most chemotherapy is the multiple and sometimes harmful side effects, as well as the ineffectiveness of treatment after a certain time due to the capability of the virus to mutate and become resistant to the treatment.Although we do not intend to model effects of resistance or side effects, we impose a condition called a limited treatment window, that monitors the global effects of these phenomena.The treatment lasts for a given period from time t 0 to t f , say.In clinical practice, antiretroviral therapy is initiated at t 0 , the time at which CD4 T-cell counts reach 350 cells/μL.We would like to maximize levels of healthy CD4 T-cells, minimize level of viral load, and also we want to minimize the systemic costs of treatments.Thus, the following objective functionals are to be maximized simultaneously Therefore, the optimal drug regimen problem can be represented as Maximize u∈U {I 1 u , I 2 u , I 3 u , I 4 u }.

3.3
In general, there does not exist a pair of control functions u P , u R ∈ U that renders the maximum value to each functional I i , i 1, 2, 3, 4, simultaneously, and one uses the concept of Pareto optimality in the sense of following definition.
Definition 3.1.A pair u * u * P , u * R ∈ U is said to be an Pareto optimal solution of the problem 3.3 if, and only if, there exist no u P , u R ∈ U such that I i u * ≤ I i u for all i ∈ {1, 2, 3, 4}, and I i u * < I i u for some i ∈ {1, 2, 3, 4}.
As seen from Definition 3.1, in general there exist an infinite number of Pareto optimal solutions.Actually, we should select one control function among the set of Pareto optimal solutions.We are going to find a Pareto optimal solution of problem 3.3 , combining the techniques and methods from the classic multi-objective optimization and control theory.

The Weighted Sum Method
The weighted sum method scalarizes set of objectives into a single-objective by multiplying each objective with user supplied weights.The weight of an objective is usually chosen in proportion to the objective's relative importance in the problem.After, a composite objective functional I u can be formed by summing the weighted objectives and the MOOCP given Proof.Assume u * ∈ U is the optimal solution of problem 4.1 .If u * is not a Pareto optimal solution, then there is u ∈ U such that I i u * ≤ I i u for all i ∈ {1, 2, 3, 4}, and I j u * < I j u for some j ∈ {1, 2, 3, 4}.Since w i > 0 for all i 1, 2, 3, 4, thus, I u * < I u , and this contradicts the optimality of u * .

Continuous Solutions
Here, we followed 8, 32 in assuming that systemic costs of the PI and RTI drugs treatment are proportional to u 2 P t and u 2 R t at time t, respectively, m 2 .We note that the existence of the solution can be obtained using a result by Fleming and Rishel 36 .That is rather straightforward to show that the right-hand sides of 2.1 -2.6 are bounded by a linear function of the state and control variables and that the integrand of the composed objective functional I u is convex on U and is bounded.We now proceed to compute candidates for a Pareto optimal solution.To this end, we apply the Pontryagin's Maximum Principle 37 and begin by defining the Lagrangian which is the Hamiltonian augmented with penalty terms for the constraints to be where ω ij t ≥ 0 are the penalty multipliers satisfying

4.3
Here, u * P , u * R is the optimal control pair yet to be found.Thus, the Maximum Principle gives the existence of adjoint variables satisfying where λ i t f 0, i 1, . . ., 6, are the transversality conditions.The Lagrangian is maximized with respect to u P and u R at the optimal pair u * P , u * R , so the partial derivatives of the Lagrangian with respect to u P and u R at u Combining all the three cases in compact form gives Using similar arguments, we also obtain the following expression for the second optimal control function We point out that the optimality system consists of the state system 2.1 -2.6 with the initial conditions, the adjoint or co-state system 4.4 with the terminal conditions, together with the expressions 4.12 and 4.13 for the control functions.The results of implementing the optimal control policy obtained by this procedure are shown in Section 6.

STI Solution
From a therapeutic point of view, it may be unsafe to administrate drugs at a dose less than maximum because virus mutations may occur see, e.g., 21 and references therein for a more exhaustive discussion on this point .Therefore, standard HAART protocols require persistent drugs uptake at maximum value.However, a number of clinical and theoretical studies attempted STI protocols in which periods of therapy at maximum dosages are alternated with periods of treatment suspension 21, 38 .The reasons for these attempts can be found in several side effects of HAART, such as serious hepatic damages and the high cost of the therapy, but also in the evidence that appropriate suspension periods may enhance the immune response of the patient.
In this section, we follow 24 in assuming that the cost of the drug regime varies linearly with the amount of drug administered.Thus, I 3 and I 4 are considered as follow:

4.14
We now proceed to compute candidates for a Pareto STI solution.To this end, we apply the Pontryagin's Maximum Principle again and begin by defining the Hamiltonian to be: and the maximum principle requires that in which the optimal controls and corresponding states and costates satisfying 2.1 -2.6 and 4.4 are indicated by a star superscript * .From 4.16 and by performing straightforward calculations it is concluded that the optimal pair u * P , u * R should maximize the following expression: Therefore,

4.18
where S 1 w 4 λ * 6 t and S 2 w 3 kλ * 5 t y * t are the switching functions which determine the type of optimal control pair u * P , u * R .Now we consider an opportunity for the control pair u * P , u * R to contain singular arcs.Let us assume that the switching function S 1 is zero on the singular interval J 1 t 1 , t 2 .To find a singular control u R,sin g , we use the fact that S 1 0, and Ṡ1 0, on J 1 ; hence, λ * 6 t 0. Therefore, from the last two equations 4.4 , it is concluded that where,

4.21
Differentiating both sides of 2.1 and then substituting 4.20 into the resulting expression, lead to where Obviously, K1 can be calculated using the Chain rule.Now, we find the value of singular control u P,sin g .With respect to this fact that d S 2 /dt 0, 0, 1, 2, and from second last equation 4.4 and 2.2 it is easy to see that during the corresponding singular interval, say

4.24
Moreover, differentiating both sides of 2.2 and then substituting 4.24 into the resulting expression yield

Mathematical Problems in Engineering 11
Note that the sets J 1 and J 2 are considered such that u P,sin g , u R,sin g ∈ U.Moreover, for the control pair u * P , u * R to be optimal, the Kelley's condition 39 must be satisfied where γ and ζ are known as the order of singularity.
The above discussion shows that using linear cost functionals, while ignoring singular arcs, leads to an STI solution.Unfortunately, solving the corresponding optimality system is not analytically possible and the gradient method 40 fails to converge; therefore, we resort to model predictive control-MPC-based approach.MPC is a control technology widely used in many areas, especially in the process industries 41 , for systems with a large number of controlled and manipulated variables, which interact significantly.In general, feedback control technologies, and MPC in particular, have started to gain significant attention in the biomedical area 12, 19-21 .A thorough overview of the history of MPC and its various incarnations can be found in 42 .In order to give a formal description of the proposed MPC algorithm, some preliminary definitions are necessary.Let Δ be the sampling time, we denote with X i the state of the system at time t i iΔ, that is,

4.28
Next, let U i u P t i u R t i u P i u R i denote the vector of PI and RTI drugs taken in the time interval iΔ, i 1 Δ , and rewrite the system model 2.1 -2.6 in the integrated form X i 1 F X i , U i , where F • is a vector-value function obtained by numerical integration of 2.1 -2.6 over the sampling Δ, assuming a constant drug uptake during the sampling time.With these definitions, it is now possible to state the following finite-horizon optimal control problem FHOCP to be solved at each discrete time k: in which N is an integer control horizon , and k be a maximizing control input sequence.The first input in the optimal sequence, that is, U * k , is injected into the system, and the entire optimization is repeated at subsequent control intervals.

Piecewise Constant Solution Using Fuzzy Aggregation and Embedding Method
In this section, as in continuous solutions, we set m 2. With respect to the nature of controls u P and u R , that indicate the PI and RTI drug efficacies as a function of time, piecewise constant 1 control is more practical from the clinical point of view.In this section, we are going to propose such kind of solutions for problem 3.3 .
Setting, ξ x, y, w, z, v, r and u u P , u R , the system of differential equations 2.1 -2.6 can be represented in a generalized form as ξ t g t, ξ t , u t In the absence of importance on the objectives and without knowledge of the possible level of attainment for those objectives, finding a Pareto optimal solution that the best satisfies the decision maker can be viewed as a fuzzy problem.In this situation, for each of the objective functionals, I i , i 1, 2, 3, 4, of problem 3.3 we associate a linear membership function μ i I i u as where m i or M i denotes the value of the objective functional such that the degree of membership function is 0 or 1, respectively, and it is depicted in Figure 1.Following the Proof.If u * is not a Pareto optimal solution for the problem 3.3 , then there exists u ∈ U such that I i u * ≤ I i u for all i ∈ {1, 2, 3, 4}, and I j u * < I j u for some j ∈ {1, 2, 3, 4}.Obviously, λ * min i 1,...,4 {μ i I i u * } ≤ min i 1,...,4 {μ i I i u } λ and 4 i 1 I i u * < 4 i 1 I i u .Therefore, λ * ρ 4 i 1 I i u * < λ ρ 4 i 1 I i u .This contradicts the assumption that λ * , u * is the optimal solution for problem 5.4 -5.6 .
Using the measure theory for solving optimal control problems based on the idea of Young 44 , which was applied for the first time by Wilson and Rubio 45 , has been theoretically established by Rubio in 46 .In the rest of this section, we follow their approach and we transfer the problem 5.4 -5.6 into a modified problem in an appropriate space that the existence of the optimal solution λ * , u * is guaranteed.

Transformation to Functional Space
We assume that state variable ξ • and control input u • get their values in the compact sets Definition 5.2.A pair p ξ, u is said to be admissible if the following conditions hold.
i The vector function ξ • is absolutely continuous and belongs to A for all t ∈ J.
ii The function u • takes its values in the set U and is Lebesgue measurable on J.
iii p satisfies in 5.5 and 5.6 on J 0 , that is, the interior set of J.

Mathematical Problems in Engineering
We assume that the set of all admissible pairs is nonempty and denote it by W. Let p be an admissible pair and B an open ball in R 6 containing J × A, and C B the space of all realvalued continuous differentiable functions on it.Let ϕ ∈ C B and define ϕ g as follows: for all ϕ ∈ C B .Let D J 0 be the space of infinitely differentiable all real-valued functions with compact support in J 0 .Define ψ n t, ξ t , u t ξ n t ψ t g n t, ξ t , u t ψ t , n 1, . . ., 6, ∀ψ ∈ D J 0 .5.9 Assume p ξ, u be an admissible pair.Since the function ψ • has compact support in J 0 , so, ψ t 0 ψ t f 0. Thus, for n 1, . . ., 6, and for all ψ ∈ D J 0 , by integrating both sides of 5.9 and using integration by parts, we have where

Transformation to Measure Space
Let M Ω denote the space of all positive Radon measures on Ω.By the Riesz representation theorem 46 , there exists a unique positive Radon measure μ on Ω such that So, we may change the space of optimization problem to measure space.In other words, the optimization problem in functional space 5.13 -5.17 can be replaced by the following new problem in measure space: Obviously, λ takes its values in the compact set 0, 1 ⊆ R. We shall consider the maximization of 5.21 over the set Q of all pairs μ, λ ∈ M Ω × 0, 1 satisfying 5.20 -5.23 .The main advantages of considering this measure theoretic form of the problem is the existence of optimal pair μ * , λ * ∈ Q, where this point can be studied in a straightforward manner without having to impose conditions such as convexity which may be artificial.Define function σ : Q → R, as σ μ, λ λ μ f 0 .The following theorem guarantees the existence of an optimal solution.Theorem 5.4.Function σ attains its maximum μ * , λ * on the set Q.
Proof.The set Q can be written as where

5.25
We will show in the next section that 5.22 and 5.23 are special cases of 5.21 .Therefore, the set Q 2 can be written as Q 2 Π × 0, 1 , where

5.26
It is well known that the set {μ ∈ M Ω : μ 1 t f − t 0 } is compact in weak * -topology 46 .Furthermore, for ϕ ∈ C B , the set Π as intersection of inverse image of the closed singleton sets {Δϕ} under the continuous functions μ → μ ϕ g is also closed.Thus, Π is a closed subset of a compact set.This proves the compactness of the set Π. Hence, the compactness of 0, 1 implies the compactness of Q 2 with the product topology.Besides, for i 1, 2, 3, 4, the set Q 1 as intersection of the inverse image of the closed sets m i , ∞ under the continuous functions μ, λ → μ f i − λ M i − m i is also closed.Therefore, the set Q as a closed subset of the compact set Q 2 is compact also.Since the function σ mapping the compact set Q on the real line is continuous so it has a maximum on the compact set Q.

Approximation
Since M Ω is an infinite-dimensional space, the problem 5.19 -5.23 is an infinitedimensional linear programming problem and we are mainly interested in approximating it.First, the maximization of σ is considered not over the set Q, but over a subset of it denoted by requiring that only a finite number of constraints 5.21 -5.23 be satisfied.

Proposition 5.6. Consider the linear program consisting of maximizing the function σ over the set
Δϕ j , j 1, . . ., K.

5.28
We have μ ∈ R and functional f → μ f is linear.Therefore, μ ϕ

5.29
Since the right-hand side of the above inequality tends to zero as k → ∞, while left-hand side is independent of k, therefore, μ ϕ g Δϕ.Thus R ⊆ Q and η * ≥ ζ which implies ζ η * .
Proposition 5.7.An optimal pair μ * , λ * in the set Q K at which the function σ attains its maximum can be found where the measure μ * has the form where α * j ≥ 0, z * j ∈ Ω, and δ z is unitary atomic measure with the support being the singleton set {z * j }, characterized by δ z F F z , z ∈ Ω.
Proof.Assume that μ * , λ * is an optimal solution of maximizing the function σ over the set Q K .Denote the set of all measures μ ∈ M Ω satisfying in 5.27 for a fixed λ λ * , by Q λ * .The set Q λ * is a weak * -compact subset of M Ω 46 .Therefore, from Theorem A.5 of 46 , Q λ * has at least one extreme point in the form of 5.30 .
Therefore, our attention is restricted to finding a measure in the form of 5.30 which maximizes the function σ and satisfies in 5.20 and K number of the constraints 5.21 -5.23 .Thus, by choosing the functions ϕ i , i 1, 2, . . ., K 1 , ψ k , k 1, . . ., K 2 , and ϑ s , s 1, . . ., S, the infinite-dimensional problem 5.19 -5.23 is approximated by the following finite-dimensional nonlinear programming NLP problem:

5.31
where K K 1 6K 2 S. Clearly, 5.31 is an NLP problem with 2 K 4 1 unknowns λ, α j , z j , j 1, . . ., K 4. One is interested in LP problem.The following proposition enables us to approximate the NLP problem 5.31 by a finite-dimensional LP problem.Proposition 5.8.Let Ω N {y 1 , y 2 , . . ., y N } be a countable dense subset of Ω, where N is asufficiently large number.Given ε > 0, a measure v ∈ M Ω can be found such that
For constructing a suitable set Ω N , J is divided to S subintervals as follows: where Δ t f − t 0 /S.Moreover, the intervals A i i 1, . . ., 6 and U j j 1, 2 are divided, respectively, into n i and s j subintervals.So, the set Ω is divided into N Sn 1 n 2 n 3 n 4 n 5 n 6 s 1 s 2 cells.One point is chosen from each cell.In this way, we will have a grid of points, which are numbered sequentially as y j t j , ξ 1 j , . . ., ξ 6 j , u P j , u R j , j 1, . . ., N. The desired numbers m i and M i in definition of the membership functions μ i I i u , i 1, 2, 3, 4 can be set

5.36
Remark 5.9.With respect to this fact that the maximum level of healthy cells is attained with full HAART, we have Similarly, the other numbers m i and M i can be calculated without explicitly solving auxiliary optimal control problems.The values of these variables are shown in Table 3.
Therefore, according to 5.33 the NLP problem 5.31 is converted to the following LP problem: α j ϑ s y j a ϑ s , s 1, . . ., S.

5.38
In addition, we choose some functions with compact support in the following form 46 : Finally, the following functions are considered that are dependent on t only: where J s , s 1, . . ., S, are given by 5.35 .These functions are used to construct the approximate piecewise constant controls 46 .Of course, we need only to construct the control function u • , since ξ • can be obtained by solving the ODEs 2.1 -2.6 .By using simplex method, a nonzero optimal solution α * of the LP problem 5.37 can be found where k cannot exceed the number of constraints, that is, k ≤ K 4. Setting α * i 0 t 0 , a piecewise constant control function u • approximating the optimal solution is constructed based on these nonzero coefficients as follows 46 : , j 1, 2, . . ., k, 5.41 where u P i j and u R i j are, respectively, the 8th and the 9th components of grid point y i j .

Numerical Results
In this section we show the response of the HIV model to the continuous, STI and piecewise constant controls presented in Sections 4 and 5.

Continuous Solutions
As we discussed in Section 4.1, in the case of continuous solutions the optimality systems, is a two-point boundary value problem because of the forward-in-time nature of the state system 2.1 -2.6 with the initial conditions and the backward-in-time nature of the adjoint system 4.4 with the terminal conditions.We consider a gradient method for the solution of the optimal control problem.The algorithm proceeds as follows.
i Choose an initial guess of control.
ii Solve the state system forward in time with initial conditions using an initial guess of control.
iii Solve the adjoint system backward in time with terminal conditions.
iv Update the controls in each iteration using the optimality conditions 4.12 and 4.13 .
v Continue the iterations until convergence is achieved.
For further discussion of the gradient method, we refer the interested reader to 40 .We use the parameters for solving the optimality system from Table 1.Antiretroviral therapy is initiated at t 0 , the time at which the CD4 T-cell count falls below 350 cells/μL, and treatment was simulated for 750 days.Furthermore, control variables u P and u R are bounded by a 1 0, b 1 0.7, a 2 0, and b 2 9 × 10 −10 , and the following initial condition is used see 7 : x 0 10 3 cells μL −1 , v 0 104 copies mL −1 , y 0 0 cells μL −1 , w 0 10 −3 cells mL −1 , z 0 10 −7 cells μL −1 , r 0 2 × 10 −7 mL copies −1 day −1 .

6.1
We ran simulations with five different values of weight factor w. The corresponding control functions are presented in Figure 2. As we increase w 3 and w 4 , thereby increasing the cost of therapy, the control functions are decreasing.In Figure 8   population is much larger than the magnitude of the cost of the drugs treatment in the objective functional I u , these differences in magnitudes are balanced by this choice of w.Moreover, the performance of the gradient method with this choice of w is shown in Figure 3.The number of required iterations to achieve the convergence is 307.Note that in Figure 3, • denotes the infinity norm, and x i denotes the state x t in the ith iteration.

MPC-Based STI Solutions
We implement therapy protocols based on the MPC algorithm 4.29 described in Section 4.
All computations are carried out by MATLAB.We choose a sampling time of five days, Δ 5. Consequently, we do not worry about creating an explicit discretization of our differential equation; we simply use a numerical simulator to approximate our discretization.Moreover,   we chose horizon length N 8.In particular, a finite horizon and a finite control space mean that we have, for each horizon, a finite number of possible control sequences.We implement MPC by using the Simulated Annealing SA metaheuristic combined with the Steepest Descent Explorer 47 for searching this space.
We ran simulations with three different values of weight factor w. The corresponding control functions are plotted in Figures 4, 5, and 6.In Figure 8, we present only the response of the HIV model to the STI solution corresponding to weights w 1 1, w 2 0.1, w 3 1800, and w 4 5 × 10 5 for brevity.
6.2 Furthermore, we set K 2 2 and S 6.By using controllability on the dynamical control system, the considered ranges for states and controls and the corresponding number of partitions in the construction of the set Ω N are summarized in Table 2.Note that the selected values from the set U 1 in the construction of y j s are 0, 0.4, and 0.7.These values indicate off, moderate, and strong PI-therapy, respectively.Similarly, the corresponding values for RTI control are 0, 5 × 10 −10 , and 9 × 10 −10 7 .The desired values for m i s and M i s are summarized in Table 3. Implementing the corresponding LP model 5.37 , we achieve λ * 0.3902, which is close to the obtained number from Table 3, that is, λ min i 1,...,4 {μ i I i u } 0.3867.Figure 7 shows the resulting piecewise constant control pair.We solved the HIV model 2.1 -2.6 with no treatment u P u R ≡ 0 and with fully efficacious treatment u P ≡ 0.7, u R ≡ 9 × 10 −10 .The numerical results for these two cases as well as the response of the systems to piecewise constant control, STI control, and continuous control are shown in Figure 8 for comparison.Except for STI solution, which is derived by considering linear cost functionals, the values of the objective functionals as well as corresponding values of the membership functions for other cases are given in Table 3.
From Figures 8 a and 8 b , we see drop in the number of CD4 T cells, and a rise in viral load following the initial infection until about the third month.After this time, CD4 T cells start recovering and virus starts decreasing due to the immune response, but can never eradicate virus completely.Then CD4 T cells level decreases and viral load increases due to destruction of immune system in absence of treatment.d interestingly, a decrease in CTLs occur in response to therapy can be observed.The extent of the decrease is directly correlated with the increase in treatment effectiveness which is consistent with experimental findings 48 .From Figure 8 a , we conclude that in the presence of treatment, the level of CD4 T cells decreases very slowly as compared with untreated patient.Moreover, this figure shows that the average CD4 T cells concentration with piecewise constant control is higher than continuous control.Interestingly, from Figure 8 d we see that the intensity of immune response in the presence of continuous control is more than full treated patient until about the 42th month, and the intensity of immune response with piecewise constant control is higher than other controls after about the 31th month.

Conclusion
In this paper, a six-dimensional HIV model is considered, and a multi-objective optimal control problem is proposed to design therapy protocols to treat the HIV infection.The model permits drug therapy RTIs and PIs as controllers.We derived continuous treatment strategies by solving the corresponding optimality systems with the gradient method.In addition, structured treatment interruption STI protocols are found by using ideas from model predictive control.The results of simulations show that the proposed MPC-based method can effectively be applied to find STI-type solutions.In particular, a novel idea based on fuzzy aggregation and measure theory is proposed to find a compromise Pareto optimal solution, and numerical results confirmed the effectiveness of this approach.The method is not iterative and it does not need any initial guess of the solution and it can be utilized to solve multi-objective optimal control problems.

MathematicalTheorem 4 . 1 .
Problems in Engineering in 3.3 is converted to a single-objective optimization problem as follows: The solution of the weighted sum problem 4.1 is Pareto optimal if the weighting coefficients are positive, that is, w i > 0 for all i 1, 2, 3, 4.

Figure 1 :
Figure 1: Linear membership function for the ith fuzzy objective.

Figure 3 :
Figure 3: Convergence of the gradient method.
, we plot only the response of the HIV model to the obtained control function corresponding to weights w 1 1, w 2 0.1, w 3 21000, and w 4 10 21 for brevity.Since the magnitude of CD4 T cell, and virus 0

Figure 7 :
Figure 7: The piecewise constant control pair u * P and u * R .

Figure 8 :
Figure 8: Dynamic behavior of the state variables versus time in the presence of continues controls the dashed line , piecewise constant controls the solid line , STI controls the dotted line , with no treatment u P u R ≡ 0 the × shape , and with fully efficacious treatment u P ≡ 0.7, u R ≡ 9 × 10 −10 the dashed dotted line .
Figures 8 b and 8 c show a clear correlation between the CTLe and virus population.As the virus increases upon initial infection, CTLe increases in order to decrease the virus.Once this is accomplished, virus decreases.Then virus grows back slowly, and this triggers an increase in the CTLe population.CTLe further increases in an attempt to keep the virus at constant levels but lose the battle because of virus-induced impairment of CD4 T cell function, in absence of treatment.Memory CTL responses depend on the presence of CD4 T cell help.Figures 8 a and 8 b show that, in presence of treatment, the virus is controlled to very low levels and CD4 T cells are maintained above the critical levels.Therefore, immune response expands for relatively long time successfully.Furthermore, these Figures indicate inverse coloration between viral load and CD4 T cells level.From Figures 8 c and

Table 1 :
Parameter values for the HIV model.
∈ Ω, where Ω J × A × U.The function ϕ g is in the space C Ω , the set of all continuous functions on the compact set Ω. Since p ξ, u is an admissible pair, we have Ω is the space of all functions in C Ω that depend only on time and a ϑ is the integral of ϑ on J. Equations 5.8 , 5.10 , and 5.11 are weak forms of 5.6 .Now, we consider the following positive linear functional on C Ω Γ p : F −→ J F t, x t , u t dt, ∀F ∈ C Ω .5.12 Proposition 5.3.Transformation p → Γ p of admissible pairs in W into the linear mappings Γ p defined in 5.12 is an injection.

Table 2 :
Considered ranges for states and controls and corresponding number of partitions.

Table 3 :
The values of the objective functionals and corresponding values of the membership functions.I 1 u /μ 1 I 1 u I 2 u /μ 2 I 2 u I 3 u /μ 3 I 3 u I 4 u /μ 4 I 4 u