# Multiobjective Optimal Control of HIV Dynamics

Authors: Hassan Zarei; Ali Vahidian Kamyad; Sohrab Effati
Journal: Mathematical Problems in Engineering
 (2010)

---

## Abstract

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.


---

## Body

## 1. 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 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+ T-cells 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 single-objective 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.The papers [8, 13, 22, 24, 30, 31] consider only RTI medication while the papers [14, 32] consider only PIs. In [6, 19–21], all effects of a HAART medication are combined to one control variable in the model. In [7, 10–12, 15, 17, 18, 33–35], dynamical multidrug therapies based on RTIs and PIs are designed.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 Section2, 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.


## 2. 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](2.1)ẋ=s-dx-rxv,(2.2)ẏ=rxv-ay-fyz,(2.3)ẇ=cxyw-qyw-bw,(2.4)ż=qyw-hz,(2.5)v̇=k(1-uP)y-τv,(2.6)ṙ=r0-uR.
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 cytotoxic T-lymphocyte effectors CTLe (z) at a rate fyz. The population of cytotoxic T-cells is divided into precursors or CTLp (w), and effectors or CTLe (z). Equations (2.3) and (2.4) describe the dynamics of these compartments. In accordance with experimental, findings establishment of a lasting CTL response depends on CD4+ T-cell help, and that HIV impairs T-helper cell function. Thus, proliferation of the CTLp population is given by  cxyw and is proportional to both virus load (y) and the number of uninfected T-helper cells (x). CTLp differentiation into effectors occurs at a rate qyw. Finally, CTLe dies at a rate hz. Equation (2.5) describes the dynamics of the free virus particles (v). These free virus particles are produced from infected CD4+ T-cells at a rate  ky and are cleared at a rate τv. The model also contains an index of the intrinsic virulence or aggressiveness of the virus (r).  This index increases linearly in the case of an untreated HIV-infected individual, with a growth rate that depends on the constant r0. Finally, (2.6) describes the dynamic of this index. In model variables u and u denote protease inhibitors (PI) and reverse transcriptase inhibitors (RTI), respectively. 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 Table1. 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).Table 1
Parameter values for the HIV model.
ParametersValue/unitDescriptions7 cellsμL-1 day-1Production(source) rate of CD4+ T-cellsd7×10-3 day-1Death rate of uninfected CD4+ T-cell populationa0.0999 day-1Infected CD4+ T-cell death ratef2μL cells-1 day-1Rate at which infected cells are killed by CTL effectorsc5×10-6 μL2 cells-2 day-1Proliferation rate of CTL precursorsq6×10-4 μL cells-1 day-1Rate at which CTL effectors are produced from CTL precursorsb0.017 day-1Death rate of CTL precursorsh0.06 day-1Death rate of CTL effectorsk300 copies mL-1 cells-1 μL day-1Rate at which virus is produced from infected cells.τ0.2 day-1Virus natural death rater010-9 copies-1 mL day-2Aggressiveness growth rate of the virus


## 3. 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 timet0 to tf, say. In clinical practice, antiretroviral therapy is initiated at t0, 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(3.1)I1(u)=∫t0tfx(t)dt,I2(u)=-∫t0tfv(t)dt,I3(u)=-∫t0tfuPm(t)dt,I4(u)=-∫t0tfuPm(t)dt,
where x(t) and v(t) are the solution of ODEs (2.1)–(2.6) corresponding to control pair u=(uP,uR) and  m is a positive integer. Let U be the set of all measurable control pairs (uP,uR) that (3.2)a1≤uP(t)≤b1,a2≤uR(t)≤b2,∀t∈[t0,tf].
Therefore, the optimal drug regimen problem can be represented as(3.3)Maximizeu∈U{I1(u),I2(u),I3(u),I4(u)}.
In general, there does not exist a pair of control functions (uP,uR)∈U that renders the maximum value to each functional Ii,i=1,2,3,4, simultaneously, and one uses the concept of Pareto optimality in the sense of following definition.Definition 3.1.
A pairu*=(uP*,uR*)∈U is said to be an Pareto optimal solution of the problem (3.3) if, and only if, there exist no (uP,uR)∈U such that Ii(u*)≤Ii(u) for all i∈{1,2,3,4}, and Ii(u*)<Ii(u) for some i∈{1,2,3,4}.As seen from Definition3.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.


## 4. 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 functionalI(u) can be formed by summing the weighted objectives and the MOOCP given in (3.3) is converted to a single-objective optimization problem as follows:(4.1)Maximizeu∈UI(u)=∑j=14wjIj(u).Theorem 4.1.
The solution of the weighted sum problem (4.1) is Pareto optimal if the weighting coefficients are positive, that is, wi>0 for all i=1,2,3,4.Proof.
Assumeu*∈U is the optimal solution of problem (4.1). If u* is not a Pareto optimal solution, then there is u∈U such that Ii(u*)≤Ii(u) for all i∈{1,2,3,4}, and Ij(u*)<Ij(u) for some j∈{1,2,3,4}. Since wi>0 for all i=1,2,3,4, thus, I(u*)<I(u), and this contradicts the optimality of u*.


### 4.1. Continuous Solutions

Here, we followed [8, 32] in assuming that systemic costs of the PI and RTI drugs treatment are proportional to uP2(t) and uR2(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(4.2)L(x,y,w,z,v,r,λ1,…,λ6,uP,uR)=w1x-w2v-w3uP2-w4uR2+λ1(s-dx-rxv)+λ2(rxv-ay-fyz)+λ3(cxyw-qyw-bw)+λ4(qyw-hz)+λ5(k(1-uP)y-τv)+λ6(r0-uR)+ω11(uP-a1)+ω12(b1-uP)+ω21(uR-a2)+ω22(b2-uR),
where ωij(t)≥0 are the penalty multipliers satisfying(4.3)ω11(uP-a1)=ω12(b1-uP)=0atuP*,ω21(uR-a2)=ω22(b2-uR)=0atuR*.
Here, (uP*,uR*) is the optimal control pair yet to be found. Thus, the Maximum Principle gives the existence of adjoint variables satisfying (4.4)λ̇1=-∂L∂x=-w1+(d+vr)λ1-vrλ2-cywλ3,λ̇2=-∂L∂y=(a+fz)λ2-(cxw-qw)λ3-qwλ4-k(1-uP)λ5,λ̇3=-∂L∂w=(-cxy+qy+b)λ3-qyλ4,λ̇4=-∂L∂z=fyλ2+hλ4,λ̇5=-∂L∂v=w2+xrλ1-xrλ2+τλ5,λ̇6=-∂L∂r=xvλ1-xvλ2,
where λi(tf)=0,i=1,…,6, are the transversality conditions. The Lagrangian is maximized with respect to uP and uR at the optimal pair (uP*,uR*), so the partial derivatives of the Lagrangian with respect to uP and uR at  uP* and  uR* are zero. That is, ∂L/∂uP=0 and ∂L/∂uR=0.Since(4.5)∂L∂uP=-2w3uP-kyλ5+ω11-ω12=0,then we have(4.6)uP*(t)=-ky(t)λ5(t)+ω11(t)-ω12(t)2w3.To determine an explicit expression for the optimal control withoutω11,ω12, we consider the following three cases.(i)
On the set{t∣a1<uP*(t)<b1}, we set ω11(t)=ω12(t)=0. Hence, the optimal control is
(4.7)uP*(t)=-ky(t)λ5(t)2w3.(ii)
On the set{t∣uP*(t)=b1}, we set ω11(t)=0. Hence, the optimal control is
(4.8)uP*(t)=b1=-ky(t)λ5(t)-ω12(t)2w3which implies that
(4.9)-ky(t)λ5(t)2w3≥b1sinceω12(t)≥0.(iii)
On the set{t∣uP*(t)=a1}, we set ω12(t)=0. Hence, the optimal control is
(4.10)a1=u1*(t)=-ky(t)λ5(t)+ω11(t)2w3
which implies that(4.11)-ky(t)λ5(t)2w3≤a1sinceω11(t)≥0.Combining all the three cases in compact form gives(4.12)uP*(t)=max⁡⁡(a1,min⁡(b1,-ky(t)λ5(t)2w3)).
Using similar arguments, we also obtain the following expression for the second optimal control function(4.13)uR*(t)=max⁡⁡(a2,min⁡(b2,-λ6(t)2w4)).
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.


### 4.2. 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, I3 and I4 are considered as follow:(4.14)I3(u)=-∫t0tfuP(t)dt,I4(u)=-∫t0tfuR(t)dt.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:(4.15)H(x,y,w,z,v,r,λ1,…,λ6,uP,uR)=w1x-w2v-w3uP-w4uR+λ1(s-dx-rxv)+λ2(rxv-ay-fyz)+λ3(cxyw-qyw-bw)+λ4(qyw-hz)+λ5(k(1-uP)y-τv)+λ6(r0-uR)and the maximum principle requires that(4.16)H(x*,y*,w*,z*,v*,r*,λ1*,…,λ6*,uP,uR)≤H(x*,y*,w*,z*,v*,r*,λ1*,…,λ6*,uP*,uR*),∀(uP,uR)∈U
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 (uP*,uR*) should maximize the following expression:
(4.17)-w3uP-w4uR-kλ5*uPy*-λ6*uR.
Therefore,(4.18)uR*(t)={b2ifS1<0,uR,singifS1=0,a2ifS1>0,uP*(t)={b1ifS2<0,uP,singifS2=0,a1ifS2>0.
where S1=w4+λ6*(t) and S2=w3+kλ5*(t)y*(t) are the switching functions which determine the type of optimal control pair (uP*,uR*). Now we consider an opportunity for the control pair (uP*,uR*) to contain singular arcs. Let us assume that the switching function S1 is zero on the singular interval J1=[t1,t2]. To find a singular control uR,sing, we use the fact that S1=0, and Ṡ1=0, on J1; hence, λ̇6*(t)=0. Therefore, from the last two equations (4.4), it is concluded that
(4.19)λ1*(t)=λ2*(t),λ5*(t)=-w2τ+(ς+w2τ)eτ(t-t1)∀t∈J1,
where ς=λ5*(t1). Since λ̇1*(t)=λ̇2*(t), equating right-hand sides of first two equations (4.4) and then substituting (4.19) into the resulting expression yield(4.20)x*=K1(t,y*,w*,z*,λ1*,λ3*,λ4*,uP*),
where,
(4.21)K1(t,y,w,z,λ1,λ3,λ4,uP)=y+w1cwλ3+qc(1-λ4λ3)+(a+fz-d)λ1cwλ3-k(1-uP)cwλ3{-w2τ+(ς+w2τ)eτ(t-t1)}.
Differentiating both sides of (2.1) and then substituting (4.20) into the resulting expression, lead to(4.22)uR,sing(t)=K̈1(t,y*,w*,z*,λ1*,λ3*,λ4*,uP*)-K2(r*,x*,y*,v*,uP*)x*v*,t∈J1,where(4.23)K2(r,x,y,v,uP)=(s-dx-rxv)(-d-rv)-(r0-τr)xv-krxy(1-uP).
Obviously, K̈1 can be calculated using the Chain rule. Now, we find the value of singular control uP,sing. With respect to this fact that dℓS2/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 J2, (4.24)ÿ*=2w3(r*x*v*-ay*-fy*z*)2-2k(w2+x*r*λ1*-x*r*λ2*+τλ5*)(r*x*v*-ay*-fy*z*)y*2w3y*+ky*2λ5*.
Moreover, differentiating both sides of (2.2) and then substituting (4.24) into the resulting expression yield
(4.25)uP,sing(t)=K3(x*,y*,w*,z*,v*,r*,uR*)kr*x*y*-2(r*x*v*-ay*-fy*z*)r*x*(w3+kλ5*){w3(r*x*v*-ay*-fy*z*)ky*2-λ5*},t∈J2,
where(4.26)K3(x,y,w,z,v,r,uR)=(r0-uR)xv+srv-(a+d+τ)xrv-r2v2x+krxy+a2y+(2af+fh)yz-fqy2w-frxvz+f2yz2.Note that the setsJ1 and J2 are considered such that (uP,sing,uR,sing)∈U. Moreover, for the control pair (uP*,uR*) to be optimal, the Kelley’s condition [39] must be satisfied
(4.27)(-1)γ∂∂uR[d2γdt2γS1]≥0,(-1)ζ∂∂uP[d2ζdt2ζS2]≥0,
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 ti=iΔ, that is,(4.28)X(i)=[x(ti)y(ti)w(ti)z(ti)v(ti)r(ti)]=[x(i)y(i)w(i)z(i)v(i)r(i)].
Next, let U(i)=[uP(ti)uR(ti)]=[uP(i)uR(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:(4.29)MaximizeUk∑i=kN+k-1w1x(i)-w2v(i)-w3uP(i)-w4uR(i)subjectto:X(k)isgiven,X(j+1)=F(X(j),U(j)),j=k,…,N+k-1,U(j)∈{a1,b1}×{a2,b2},j=k,…,N+k-1,
in which  N  is an integer(control horizon), and Uk=[U(k)⋯U(N+k-1)] is the decision variable. Let Uk* 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.


## 4.1. Continuous Solutions

Here, we followed [8, 32] in assuming that systemic costs of the PI and RTI drugs treatment are proportional to uP2(t) and uR2(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(4.2)L(x,y,w,z,v,r,λ1,…,λ6,uP,uR)=w1x-w2v-w3uP2-w4uR2+λ1(s-dx-rxv)+λ2(rxv-ay-fyz)+λ3(cxyw-qyw-bw)+λ4(qyw-hz)+λ5(k(1-uP)y-τv)+λ6(r0-uR)+ω11(uP-a1)+ω12(b1-uP)+ω21(uR-a2)+ω22(b2-uR),
where ωij(t)≥0 are the penalty multipliers satisfying(4.3)ω11(uP-a1)=ω12(b1-uP)=0atuP*,ω21(uR-a2)=ω22(b2-uR)=0atuR*.
Here, (uP*,uR*) is the optimal control pair yet to be found. Thus, the Maximum Principle gives the existence of adjoint variables satisfying (4.4)λ̇1=-∂L∂x=-w1+(d+vr)λ1-vrλ2-cywλ3,λ̇2=-∂L∂y=(a+fz)λ2-(cxw-qw)λ3-qwλ4-k(1-uP)λ5,λ̇3=-∂L∂w=(-cxy+qy+b)λ3-qyλ4,λ̇4=-∂L∂z=fyλ2+hλ4,λ̇5=-∂L∂v=w2+xrλ1-xrλ2+τλ5,λ̇6=-∂L∂r=xvλ1-xvλ2,
where λi(tf)=0,i=1,…,6, are the transversality conditions. The Lagrangian is maximized with respect to uP and uR at the optimal pair (uP*,uR*), so the partial derivatives of the Lagrangian with respect to uP and uR at  uP* and  uR* are zero. That is, ∂L/∂uP=0 and ∂L/∂uR=0.Since(4.5)∂L∂uP=-2w3uP-kyλ5+ω11-ω12=0,then we have(4.6)uP*(t)=-ky(t)λ5(t)+ω11(t)-ω12(t)2w3.To determine an explicit expression for the optimal control withoutω11,ω12, we consider the following three cases.(i)
On the set{t∣a1<uP*(t)<b1}, we set ω11(t)=ω12(t)=0. Hence, the optimal control is
(4.7)uP*(t)=-ky(t)λ5(t)2w3.(ii)
On the set{t∣uP*(t)=b1}, we set ω11(t)=0. Hence, the optimal control is
(4.8)uP*(t)=b1=-ky(t)λ5(t)-ω12(t)2w3which implies that
(4.9)-ky(t)λ5(t)2w3≥b1sinceω12(t)≥0.(iii)
On the set{t∣uP*(t)=a1}, we set ω12(t)=0. Hence, the optimal control is
(4.10)a1=u1*(t)=-ky(t)λ5(t)+ω11(t)2w3
which implies that(4.11)-ky(t)λ5(t)2w3≤a1sinceω11(t)≥0.Combining all the three cases in compact form gives(4.12)uP*(t)=max⁡⁡(a1,min⁡(b1,-ky(t)λ5(t)2w3)).
Using similar arguments, we also obtain the following expression for the second optimal control function(4.13)uR*(t)=max⁡⁡(a2,min⁡(b2,-λ6(t)2w4)).
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.


## 4.2. 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, I3 and I4 are considered as follow:(4.14)I3(u)=-∫t0tfuP(t)dt,I4(u)=-∫t0tfuR(t)dt.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:(4.15)H(x,y,w,z,v,r,λ1,…,λ6,uP,uR)=w1x-w2v-w3uP-w4uR+λ1(s-dx-rxv)+λ2(rxv-ay-fyz)+λ3(cxyw-qyw-bw)+λ4(qyw-hz)+λ5(k(1-uP)y-τv)+λ6(r0-uR)and the maximum principle requires that(4.16)H(x*,y*,w*,z*,v*,r*,λ1*,…,λ6*,uP,uR)≤H(x*,y*,w*,z*,v*,r*,λ1*,…,λ6*,uP*,uR*),∀(uP,uR)∈U
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 (uP*,uR*) should maximize the following expression:
(4.17)-w3uP-w4uR-kλ5*uPy*-λ6*uR.
Therefore,(4.18)uR*(t)={b2ifS1<0,uR,singifS1=0,a2ifS1>0,uP*(t)={b1ifS2<0,uP,singifS2=0,a1ifS2>0.
where S1=w4+λ6*(t) and S2=w3+kλ5*(t)y*(t) are the switching functions which determine the type of optimal control pair (uP*,uR*). Now we consider an opportunity for the control pair (uP*,uR*) to contain singular arcs. Let us assume that the switching function S1 is zero on the singular interval J1=[t1,t2]. To find a singular control uR,sing, we use the fact that S1=0, and Ṡ1=0, on J1; hence, λ̇6*(t)=0. Therefore, from the last two equations (4.4), it is concluded that
(4.19)λ1*(t)=λ2*(t),λ5*(t)=-w2τ+(ς+w2τ)eτ(t-t1)∀t∈J1,
where ς=λ5*(t1). Since λ̇1*(t)=λ̇2*(t), equating right-hand sides of first two equations (4.4) and then substituting (4.19) into the resulting expression yield(4.20)x*=K1(t,y*,w*,z*,λ1*,λ3*,λ4*,uP*),
where,
(4.21)K1(t,y,w,z,λ1,λ3,λ4,uP)=y+w1cwλ3+qc(1-λ4λ3)+(a+fz-d)λ1cwλ3-k(1-uP)cwλ3{-w2τ+(ς+w2τ)eτ(t-t1)}.
Differentiating both sides of (2.1) and then substituting (4.20) into the resulting expression, lead to(4.22)uR,sing(t)=K̈1(t,y*,w*,z*,λ1*,λ3*,λ4*,uP*)-K2(r*,x*,y*,v*,uP*)x*v*,t∈J1,where(4.23)K2(r,x,y,v,uP)=(s-dx-rxv)(-d-rv)-(r0-τr)xv-krxy(1-uP).
Obviously, K̈1 can be calculated using the Chain rule. Now, we find the value of singular control uP,sing. With respect to this fact that dℓS2/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 J2, (4.24)ÿ*=2w3(r*x*v*-ay*-fy*z*)2-2k(w2+x*r*λ1*-x*r*λ2*+τλ5*)(r*x*v*-ay*-fy*z*)y*2w3y*+ky*2λ5*.
Moreover, differentiating both sides of (2.2) and then substituting (4.24) into the resulting expression yield
(4.25)uP,sing(t)=K3(x*,y*,w*,z*,v*,r*,uR*)kr*x*y*-2(r*x*v*-ay*-fy*z*)r*x*(w3+kλ5*){w3(r*x*v*-ay*-fy*z*)ky*2-λ5*},t∈J2,
where(4.26)K3(x,y,w,z,v,r,uR)=(r0-uR)xv+srv-(a+d+τ)xrv-r2v2x+krxy+a2y+(2af+fh)yz-fqy2w-frxvz+f2yz2.Note that the setsJ1 and J2 are considered such that (uP,sing,uR,sing)∈U. Moreover, for the control pair (uP*,uR*) to be optimal, the Kelley’s condition [39] must be satisfied
(4.27)(-1)γ∂∂uR[d2γdt2γS1]≥0,(-1)ζ∂∂uP[d2ζdt2ζS2]≥0,
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 ti=iΔ, that is,(4.28)X(i)=[x(ti)y(ti)w(ti)z(ti)v(ti)r(ti)]=[x(i)y(i)w(i)z(i)v(i)r(i)].
Next, let U(i)=[uP(ti)uR(ti)]=[uP(i)uR(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:(4.29)MaximizeUk∑i=kN+k-1w1x(i)-w2v(i)-w3uP(i)-w4uR(i)subjectto:X(k)isgiven,X(j+1)=F(X(j),U(j)),j=k,…,N+k-1,U(j)∈{a1,b1}×{a2,b2},j=k,…,N+k-1,
in which  N  is an integer(control horizon), and Uk=[U(k)⋯U(N+k-1)] is the decision variable. Let Uk* 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.


## 5. Piecewise Constant Solution Using Fuzzy Aggregation and Embedding Method

In this section, as in continuous solutions, we setm=2. With respect to the nature of controls uP and uR, that indicate the PI and RTI drug efficacies as a function of time, piecewise constant 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=(uP,uR), the system of differential equations (2.1)–(2.6) can be represented in a generalized form as(5.1)ξ̇(t)=g(t,ξ(t),u(t))=(s-dξ1-ξ1ξ5ξ6ξ1ξ5ξ6-aξ2-fξ2ξ4cξ1ξ2ξ3-qξ2ξ3-bξ3qξ2ξ3-hξ4k(1-uP)ξ2-τξ5r0-uR),ξ(0)=ξ0.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, Ii,i=1,2,3,4, of problem (3.3) we associate a linear membership function μi(Ii(u)) as
(5.2)μi(Ii(u))={0,Ii(u)<mi,Ii(u)-miMi-mi,mi≤Ii(u)≤Mi,1,Ii(u)>Mi,where mi or Mi 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 principle of the fuzzy decision of Bellman and Zadeh [43], the MOOCP (3.3) can be interpreted as the following problem:
(5.3)Maximizeu∈Umin⁡i=1,…,4{μi(Ii(u))}.Let ρ be a given small positive number. Introducing an auxiliary variable λ and by adding the augmented term ρ∑k=14Ii(u) to it, a solution to problem (3.3), can be obtained by solving the following problem:(5.4)max⁡λ,u∈Uλ+ρ∑i=14Ii(u)(5.5)subjecttomi≤Ii(u)-λ(Mi-mi),i=1,2,3,4(5.6)ξ̇=g(t,ξ,u),ξ(t0)=ξ0.Figure 1
Linear membership function for theith fuzzy objective.Theorem 5.1.
If(λ*,u*) is the optimal solution of the problem (5.4)–(5.6), then u*  is a Pareto optimal solution of (3.3).Proof.
Ifu* is not a Pareto optimal solution for the problem (3.3), then there exists u∈U such that Ii(u*)≤Ii(u) for all i∈{1,2,3,4}, and Ij(u*)<Ij(u) for some j∈{1,2,3,4}. Obviously, λ*=mini=1,…,4{μi(Ii(u*))}≤mini=1,…,4{μi(Ii(u))}=λ and ∑i=14Ii(u*)<∑i=14Ii(u). Therefore, λ*+ρ∑i=14Ii(u*)<λ+ρ∑i=14Ii(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.


### 5.1. Transformation to Functional Space

We assume that state variableξ(·) and control input u(·) get their values in the compact sets A=A1×⋯×A6⊂ℜ6 and U=U1×U2⊂ℜ2, respectively. Set J=[t0,tf].Definition 5.2.
A pairp=[ξ,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 functionu(·) takes its values in the set U and is Lebesgue measurable on J.(iii)
p satisfies in (5.5) and (5.6) on J0, that is, the interior set of J.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 ℜ6 containing J×A, and C′(B) the space of all real-valued continuous differentiable functions on it. Let φ∈C′(B) and define φg as follows:
(5.7)φg(t,ξ(t),u(t))=dφ(t,ξ(t))dt=∑j=16∂φ(t,ξ(t))∂ξjgj(t,ξ(t),u(t))+∂φ(t,ξ(t))∂t
for each [t,ξ(t),u(t)]∈Ω, 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
(5.8)∫t0tfφg(t,ξ(t),u(t))dt=φ(tf,ξ(tf))-φ(t0,ξ(t0))=Δφ,
for all φ∈C′(B). Let D(J0) be the space of infinitely differentiable all real-valued functions with compact support in J0.  Define
(5.9)ψn(t,ξ(t),u(t))=ξn(t)ψ′(t)+gn(t,ξ(t),u(t))ψ(t),n=1,…,6,∀ψ∈D(J0).
Assume p=[ξ,u] be an admissible pair. Since the function ψ(·) has compact support in J0, so, ψ(t0)=ψ(tf)=0. Thus, for n=1,…,6, and for all ψ∈D(J0), by integrating both sides of (5.9) and using integration by parts, we have
(5.10)∫t0tfψn(t,ξ(t),u(t))dt=∫t0tfξn(t)ψ′(t)dt+∫t0tfgn(t,ξ(t),u(t))ψ(t)dt=ξn(t)ψ(t)∣t0tf-∫t0tf(ξ̇n(t)-gn(t,ξ(t),u(t)))ψ(t)dt=0.
Also by choosing the functions which are dependent only on time, we have
(5.11)∫t0tfϑ(t,ξ(t),u(t))dt=aϑ,∀ϑ∈C1(Ω),
where C1(Ω) 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(Ω)(5.12)Γp:F⟶∫JF(t,x(t),u(t))dt,∀F∈C(Ω).Proposition 5.3.
Transformationp→Γp of admissible pairs in W into the linear mappings Γp defined in (5.12) is an injection.Proof.
See [46].Thus, the problem (5.4)–(5.6) is converted to the following optimization problem in functional space:
(5.13)maximizeλ,p∈Wλ+Γp(f0)(5.14)subjecttomi≤Γp(fi)-λ(Mi-mi)i=1,2,3,4,(5.15)Γp(φg)=Δφ,φ∈C′(B),(5.16)Γp(ψn)=0,n=1,…,6,ψ∈D(J0),(5.17)Γp(ϑ)=aϑ,ϑ∈C1(Ω),
where f0(t,ξ,u)=ρ(ξ1-ξ5-uP2-uR2), f1(t,ξ,u)=ξ1, f2(t,ξ,u)=-ξ5, f3(t,ξ,u)=-uP2, and f4(t,ξ,u)=-uR2.


### 5.2. Transformation to Measure Space

LetM+(Ω) 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
(5.18)Γp(F)=∫JF(t,ξ(t),u(t))dt=∫ΩF(t,ξ,u)dμ≡μ(F),F∈C(Ω).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:
(5.19)maximizeλ,μ∈M+(Ω)λ+μ(f0)(5.20)subjecttomi≤μ(fi)-λ(Mi-mi),i=1,2,3,4,(5.21)μ(φg)=Δφ,φ∈C′(B),(5.22)μ(ψn)=0,n=1,…,6,ψ∈D(J0),(5.23)μ(ϑ)=aϑ,ϑ∈C1(Ω).
Obviously, λ takes its values in the compact set [0,1]⊆ℜ. 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→ℜ, as σ(μ,λ)=λ+μ(f0). The following theorem guarantees the existence of an optimal solution.Theorem 5.4.
Functionσ attains its maximum (μ*,λ*) on the setQ.Proof.
The setQcan be written as
(5.24)Q=Q1∩Q2,
where
(5.25)Q1={(μ,λ)∈M+(Ω)×[0,1]:mi≤μ(fi)-λ(Mi-mi),i=1,2,3,4},Q2={(μ,λ)∈M+(Ω)×[0,1]:μsatisfiesin(5.21)–(5.23)}.
We will show in the next section that (5.22) and (5.23) are special cases of (5.21). Therefore, the set Q2 can be written as Q2=Π×[0,1], where
(5.26)Π=⋂φ∈C′(B){μ∈M+(Ω):μ(φg)=Δφ}.
It is well known that the set {μ∈M+(Ω):μ(1)=tf-t0} 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 Q2 with the product topology. Besides, for i=1,2,3,4, the set Q1 as intersection of the inverse image of the closed sets [mi,∞) under the continuous functions (μ,λ)→μ(fi)-λ(Mi-mi) is also closed. Therefore, the set Q as a closed subset of the compact set Q2 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.


### 5.3. Approximation

SinceM+(Ω) is an infinite-dimensional space, the problem (5.19)–(5.23) is an infinite-dimensional 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. Let {φi:i=1,2,…}, {ψj:j=1,2,…}, and {ϑs:s=1,2,…} be the sets of total functions, respectively, in C′(B), D(J0) and C1(Ω).Remark 5.5.
From (5.7) and (5.9), it can be seen that (5.22) and (5.23) are also achieved from (5.21) by setting φ(t,ξ(t))=ξn(t)ψ(t) and φ(t,ξ(t))=∫0tϑ(τ)dτ, respectively.Now, the first approximation will be completed by choosing finite number of total functions and using following propositions.Proposition 5.6.
Consider the linear program consisting of maximizing the functionσ over the set Q(K) of pairs (μ,λ) in M+(Ω)×[0,1] satisfying
(5.27)mi≤μ(fi)-λ(Mi-mi)i=1,2,3,4,μ(φjg)=Δφj,j=1,…,K.
Then, ηK≡maxQ(K)σ tends to η*≡maxQσ as K→∞.Proof.
We haveQ(1)⊇Q(2)⊇⋯⊇Q(M)⊇⋯⊇Q and hence, η1≥η2≥⋯≥ηM≥⋯≥η*. Therefore, {ηj} is nonincreasing and bounded sequence then converges to a number ζ such that ζ≥η*. We show that, ζ=η*. Set R≡⋂K=1∞Q(K). Then, R⊇Q and ζ≡maxRσ. It is sufficient to show that R⊆Q. Assume (μ,λ)∈R and φ∈C′(B). Since the Linear combinations of the functions {φj,j=1,2,…} are uniformly dense in C′(B), there is the sequence {φ̃k}∈span{φj,j=1,2,…} such that φ̃k tends to φ uniformly as k→∞. Hence, S1, S2, and S3 tend to zero as k→∞, where,
(5.28)S1=sup⁡|φξ(t,ξ)-φ̃kξ(t,ξ)|,S2=sup⁡|φt(t,ξ)-φ̃kt(t,ξ)|,S3=sup⁡|φ(t,ξ)-φ̃k(t,ξ)|.
We have μ∈R and functional f→μ(f) is linear. Therefore, μ(φ̃kg)=Δφ̃k and
(5.29)|μ(φg)-Δφ|=|μ(φg)-Δφ-μ(φ̃kg)+Δφ̃k|=|∫Ω{[φξ(t,ξ)-φ̃kξ(t,ξ)]g(t,ξ,u)+[φt(t,ξ)-φ̃kt(t,ξ)]}dμ-(Δφ-Δφ̃k)|≤S1∫Ω|g(t,ξ,u)|dμ+S2∫Ωdμ+2S3.
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
(5.30)μ*=∑j=1K+4αj*δ(zj*),
where αj*≥0,zj*∈Ω, and δ(z) is unitary atomic measure with the support being the singleton set {zj*}, 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,…,K1, ψk,k=1,…,K2, 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)maximizeλ,αj≥0,zj∈Ωλ+∑j=1K+4αjf0(zj)subjecttomi≤∑j=1K+4αjfi(zj)-λ(Mi-mi),i=1,2,3,4,∑j=1K+4αjφig(zj)=Δφi,i=1,…,K1,∑j=1K+4αjψkn(zj)=0,k=1,…,K2,n=1,…,6,∑j=1K+4αjϑs(zj)=aϑs,s=1,…,S,
where K=K1+6K2+S. Clearly, (5.31) is an NLP problem with 2(K+4)+1 unknowns λ, αj,zj,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={y1,y2,…,yN} be a countable dense subset of Ω, where N  is a sufficiently large number. Given ε>0, a measure v∈M+(Ω) can be found such that
(5.32)|v(fi)-μ*(fi)|≤ε,i=0,1,2,3,4,|v(φjg)-μ*(φjg)|≤ε,j=1,…,K1,|v(ψkn)-μ*(ψkn)|≤ε,k=1,…,K2,n=1,…,6,|v(ϑs)-μ*(ϑs)|≤ε,s=1,…,S,
where the measure v has the form
(5.33)v=∑j=1K+4αj*δ(yj),
and the coefficients αj*,j=1,…,K+4, are the same as optimal measure (5.30), and yj∈ΩN,j=1,…,K+4.Proof.
We rename the functionsfis, φjgs, ψkns, and ϑss sequentially as hj,j=1,2,…,K+5. Then, for j=1,…,K+5,
(5.34)|(μ*-v)hj|=|∑i=1K+4αi*[hj(zi*)-hj(yi)]|≤(∑i=1K+4αi*)max⁡i,j|hj(zi*)-hj(yi)|.
The functions hj,j=1,2,…,K+5, are continues. Therefore, maxi,j|hj(zi*)-hj(yi)| can be made less than ε/∑j=1K+4αj* by choosing yi,i=1,2,…,K+4, sufficiently near zi*.For constructing a suitable setΩN, J is divided to S subintervals as follows:
(5.35)Js=[t0+(s-1)Δ,t0+sΔ],s=1,2,…,S,
where Δ=(tf-t0)/S. Moreover, the intervals Ai(i=1,…,6) and Uj(j=1,2) are divided, respectively, into ni and sj subintervals. So, the set Ω is divided into N=Sn1n2n3n4n5n6s1s2 cells. One point is chosen from each cell. In this way, we will have a grid of points, which are numbered sequentially as yj=(tj,ξ1j,…,ξ6j,uPj,uRj), j=1,…,N.The desired numbersmi and Mi in definition of the membership functions μi(Ii(u)),i=1,2,3,4 can be set(5.36)mi=min⁡u∈U∫t0tffi(t,ξ,u)dt,Mi=max⁡u∈U∫t0tffi(t,ξ,u)dt.Remark 5.9.
With respect to this fact that the maximum level of healthy cells is attained with full HAART, we haveM1=I1(b1,b2). Similarly, the other numbers mi and Mi 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:(5.37)maximizeλ,αj≥0λ+∑j=1Nαjf0(yj)subjecttomi≤∑j=1Nαjfi(yj)-λ(Mi-mi)i=1,2,3,4,∑j=1Nαjφig(yj)=Δφi,i=1,…,K1,∑j=1Nαjψkn(yj)=0,k=1,…,K2,n=1,…,6,∑j=1Nαjϑs(yj)=aϑs,s=1,…,S.Here, we discuss suitable total functionsφi,  i=1,…,K1,ψk,k=1,…,K2, and ϑs,s=1,…,S. The functions φis can be taken to be monomials of t and the components of the vector ξ as follows:(5.38)tiξj,i∈{0,1,2,…},j∈{1,…,6}.In addition, we choose some functions with compact support in the following form [46]: (5.39)ψ2r-1(t)=sin(2πr(t-t0)tf-t0),ψ2r(t)=1-cos⁡(2πr(t-t0)tf-t0),r=1,2,….Finally, the following functions are considered that are dependent ont only:(5.40)ϑs(t)={1t∈Js,0otherwise,whereJs,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 αi1*,αi2*,…,αik*, i1<i2<⋯<ik, of the LP problem (5.37) can be found where k cannot exceed the number of constraints, that is, k≤K+4. Setting αi0*=t0, a piecewise constant control function u(·) approximating the optimal solution is constructed based on these nonzero coefficients as follows [46]:(5.41)u(t)={(uPij,uRij)t∈[∑h=0j-1αih*,∑h=0jαih*)0otherwise,j=1,2,…,k,whereuPijand uRijare, respectively, the 8th and the 9th components of grid point yij.


## 5.1. Transformation to Functional Space

We assume that state variableξ(·) and control input u(·) get their values in the compact sets A=A1×⋯×A6⊂ℜ6 and U=U1×U2⊂ℜ2, respectively. Set J=[t0,tf].Definition 5.2.
A pairp=[ξ,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 functionu(·) takes its values in the set U and is Lebesgue measurable on J.(iii)
p satisfies in (5.5) and (5.6) on J0, that is, the interior set of J.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 ℜ6 containing J×A, and C′(B) the space of all real-valued continuous differentiable functions on it. Let φ∈C′(B) and define φg as follows:
(5.7)φg(t,ξ(t),u(t))=dφ(t,ξ(t))dt=∑j=16∂φ(t,ξ(t))∂ξjgj(t,ξ(t),u(t))+∂φ(t,ξ(t))∂t
for each [t,ξ(t),u(t)]∈Ω, 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
(5.8)∫t0tfφg(t,ξ(t),u(t))dt=φ(tf,ξ(tf))-φ(t0,ξ(t0))=Δφ,
for all φ∈C′(B). Let D(J0) be the space of infinitely differentiable all real-valued functions with compact support in J0.  Define
(5.9)ψn(t,ξ(t),u(t))=ξn(t)ψ′(t)+gn(t,ξ(t),u(t))ψ(t),n=1,…,6,∀ψ∈D(J0).
Assume p=[ξ,u] be an admissible pair. Since the function ψ(·) has compact support in J0, so, ψ(t0)=ψ(tf)=0. Thus, for n=1,…,6, and for all ψ∈D(J0), by integrating both sides of (5.9) and using integration by parts, we have
(5.10)∫t0tfψn(t,ξ(t),u(t))dt=∫t0tfξn(t)ψ′(t)dt+∫t0tfgn(t,ξ(t),u(t))ψ(t)dt=ξn(t)ψ(t)∣t0tf-∫t0tf(ξ̇n(t)-gn(t,ξ(t),u(t)))ψ(t)dt=0.
Also by choosing the functions which are dependent only on time, we have
(5.11)∫t0tfϑ(t,ξ(t),u(t))dt=aϑ,∀ϑ∈C1(Ω),
where C1(Ω) 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(Ω)(5.12)Γp:F⟶∫JF(t,x(t),u(t))dt,∀F∈C(Ω).Proposition 5.3.
Transformationp→Γp of admissible pairs in W into the linear mappings Γp defined in (5.12) is an injection.Proof.
See [46].Thus, the problem (5.4)–(5.6) is converted to the following optimization problem in functional space:
(5.13)maximizeλ,p∈Wλ+Γp(f0)(5.14)subjecttomi≤Γp(fi)-λ(Mi-mi)i=1,2,3,4,(5.15)Γp(φg)=Δφ,φ∈C′(B),(5.16)Γp(ψn)=0,n=1,…,6,ψ∈D(J0),(5.17)Γp(ϑ)=aϑ,ϑ∈C1(Ω),
where f0(t,ξ,u)=ρ(ξ1-ξ5-uP2-uR2), f1(t,ξ,u)=ξ1, f2(t,ξ,u)=-ξ5, f3(t,ξ,u)=-uP2, and f4(t,ξ,u)=-uR2.


## 5.2. Transformation to Measure Space

LetM+(Ω) 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
(5.18)Γp(F)=∫JF(t,ξ(t),u(t))dt=∫ΩF(t,ξ,u)dμ≡μ(F),F∈C(Ω).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:
(5.19)maximizeλ,μ∈M+(Ω)λ+μ(f0)(5.20)subjecttomi≤μ(fi)-λ(Mi-mi),i=1,2,3,4,(5.21)μ(φg)=Δφ,φ∈C′(B),(5.22)μ(ψn)=0,n=1,…,6,ψ∈D(J0),(5.23)μ(ϑ)=aϑ,ϑ∈C1(Ω).
Obviously, λ takes its values in the compact set [0,1]⊆ℜ. 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→ℜ, as σ(μ,λ)=λ+μ(f0). The following theorem guarantees the existence of an optimal solution.Theorem 5.4.
Functionσ attains its maximum (μ*,λ*) on the setQ.Proof.
The setQcan be written as
(5.24)Q=Q1∩Q2,
where
(5.25)Q1={(μ,λ)∈M+(Ω)×[0,1]:mi≤μ(fi)-λ(Mi-mi),i=1,2,3,4},Q2={(μ,λ)∈M+(Ω)×[0,1]:μsatisfiesin(5.21)–(5.23)}.
We will show in the next section that (5.22) and (5.23) are special cases of (5.21). Therefore, the set Q2 can be written as Q2=Π×[0,1], where
(5.26)Π=⋂φ∈C′(B){μ∈M+(Ω):μ(φg)=Δφ}.
It is well known that the set {μ∈M+(Ω):μ(1)=tf-t0} 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 Q2 with the product topology. Besides, for i=1,2,3,4, the set Q1 as intersection of the inverse image of the closed sets [mi,∞) under the continuous functions (μ,λ)→μ(fi)-λ(Mi-mi) is also closed. Therefore, the set Q as a closed subset of the compact set Q2 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.


## 5.3. Approximation

SinceM+(Ω) is an infinite-dimensional space, the problem (5.19)–(5.23) is an infinite-dimensional 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. Let {φi:i=1,2,…}, {ψj:j=1,2,…}, and {ϑs:s=1,2,…} be the sets of total functions, respectively, in C′(B), D(J0) and C1(Ω).Remark 5.5.
From (5.7) and (5.9), it can be seen that (5.22) and (5.23) are also achieved from (5.21) by setting φ(t,ξ(t))=ξn(t)ψ(t) and φ(t,ξ(t))=∫0tϑ(τ)dτ, respectively.Now, the first approximation will be completed by choosing finite number of total functions and using following propositions.Proposition 5.6.
Consider the linear program consisting of maximizing the functionσ over the set Q(K) of pairs (μ,λ) in M+(Ω)×[0,1] satisfying
(5.27)mi≤μ(fi)-λ(Mi-mi)i=1,2,3,4,μ(φjg)=Δφj,j=1,…,K.
Then, ηK≡maxQ(K)σ tends to η*≡maxQσ as K→∞.Proof.
We haveQ(1)⊇Q(2)⊇⋯⊇Q(M)⊇⋯⊇Q and hence, η1≥η2≥⋯≥ηM≥⋯≥η*. Therefore, {ηj} is nonincreasing and bounded sequence then converges to a number ζ such that ζ≥η*. We show that, ζ=η*. Set R≡⋂K=1∞Q(K). Then, R⊇Q and ζ≡maxRσ. It is sufficient to show that R⊆Q. Assume (μ,λ)∈R and φ∈C′(B). Since the Linear combinations of the functions {φj,j=1,2,…} are uniformly dense in C′(B), there is the sequence {φ̃k}∈span{φj,j=1,2,…} such that φ̃k tends to φ uniformly as k→∞. Hence, S1, S2, and S3 tend to zero as k→∞, where,
(5.28)S1=sup⁡|φξ(t,ξ)-φ̃kξ(t,ξ)|,S2=sup⁡|φt(t,ξ)-φ̃kt(t,ξ)|,S3=sup⁡|φ(t,ξ)-φ̃k(t,ξ)|.
We have μ∈R and functional f→μ(f) is linear. Therefore, μ(φ̃kg)=Δφ̃k and
(5.29)|μ(φg)-Δφ|=|μ(φg)-Δφ-μ(φ̃kg)+Δφ̃k|=|∫Ω{[φξ(t,ξ)-φ̃kξ(t,ξ)]g(t,ξ,u)+[φt(t,ξ)-φ̃kt(t,ξ)]}dμ-(Δφ-Δφ̃k)|≤S1∫Ω|g(t,ξ,u)|dμ+S2∫Ωdμ+2S3.
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
(5.30)μ*=∑j=1K+4αj*δ(zj*),
where αj*≥0,zj*∈Ω, and δ(z) is unitary atomic measure with the support being the singleton set {zj*}, 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,…,K1, ψk,k=1,…,K2, 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)maximizeλ,αj≥0,zj∈Ωλ+∑j=1K+4αjf0(zj)subjecttomi≤∑j=1K+4αjfi(zj)-λ(Mi-mi),i=1,2,3,4,∑j=1K+4αjφig(zj)=Δφi,i=1,…,K1,∑j=1K+4αjψkn(zj)=0,k=1,…,K2,n=1,…,6,∑j=1K+4αjϑs(zj)=aϑs,s=1,…,S,
where K=K1+6K2+S. Clearly, (5.31) is an NLP problem with 2(K+4)+1 unknowns λ, αj,zj,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={y1,y2,…,yN} be a countable dense subset of Ω, where N  is a sufficiently large number. Given ε>0, a measure v∈M+(Ω) can be found such that
(5.32)|v(fi)-μ*(fi)|≤ε,i=0,1,2,3,4,|v(φjg)-μ*(φjg)|≤ε,j=1,…,K1,|v(ψkn)-μ*(ψkn)|≤ε,k=1,…,K2,n=1,…,6,|v(ϑs)-μ*(ϑs)|≤ε,s=1,…,S,
where the measure v has the form
(5.33)v=∑j=1K+4αj*δ(yj),
and the coefficients αj*,j=1,…,K+4, are the same as optimal measure (5.30), and yj∈ΩN,j=1,…,K+4.Proof.
We rename the functionsfis, φjgs, ψkns, and ϑss sequentially as hj,j=1,2,…,K+5. Then, for j=1,…,K+5,
(5.34)|(μ*-v)hj|=|∑i=1K+4αi*[hj(zi*)-hj(yi)]|≤(∑i=1K+4αi*)max⁡i,j|hj(zi*)-hj(yi)|.
The functions hj,j=1,2,…,K+5, are continues. Therefore, maxi,j|hj(zi*)-hj(yi)| can be made less than ε/∑j=1K+4αj* by choosing yi,i=1,2,…,K+4, sufficiently near zi*.For constructing a suitable setΩN, J is divided to S subintervals as follows:
(5.35)Js=[t0+(s-1)Δ,t0+sΔ],s=1,2,…,S,
where Δ=(tf-t0)/S. Moreover, the intervals Ai(i=1,…,6) and Uj(j=1,2) are divided, respectively, into ni and sj subintervals. So, the set Ω is divided into N=Sn1n2n3n4n5n6s1s2 cells. One point is chosen from each cell. In this way, we will have a grid of points, which are numbered sequentially as yj=(tj,ξ1j,…,ξ6j,uPj,uRj), j=1,…,N.The desired numbersmi and Mi in definition of the membership functions μi(Ii(u)),i=1,2,3,4 can be set(5.36)mi=min⁡u∈U∫t0tffi(t,ξ,u)dt,Mi=max⁡u∈U∫t0tffi(t,ξ,u)dt.Remark 5.9.
With respect to this fact that the maximum level of healthy cells is attained with full HAART, we haveM1=I1(b1,b2). Similarly, the other numbers mi and Mi 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:(5.37)maximizeλ,αj≥0λ+∑j=1Nαjf0(yj)subjecttomi≤∑j=1Nαjfi(yj)-λ(Mi-mi)i=1,2,3,4,∑j=1Nαjφig(yj)=Δφi,i=1,…,K1,∑j=1Nαjψkn(yj)=0,k=1,…,K2,n=1,…,6,∑j=1Nαjϑs(yj)=aϑs,s=1,…,S.Here, we discuss suitable total functionsφi,  i=1,…,K1,ψk,k=1,…,K2, and ϑs,s=1,…,S. The functions φis can be taken to be monomials of t and the components of the vector ξ as follows:(5.38)tiξj,i∈{0,1,2,…},j∈{1,…,6}.In addition, we choose some functions with compact support in the following form [46]: (5.39)ψ2r-1(t)=sin(2πr(t-t0)tf-t0),ψ2r(t)=1-cos⁡(2πr(t-t0)tf-t0),r=1,2,….Finally, the following functions are considered that are dependent ont only:(5.40)ϑs(t)={1t∈Js,0otherwise,whereJs,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 αi1*,αi2*,…,αik*, i1<i2<⋯<ik, of the LP problem (5.37) can be found where k cannot exceed the number of constraints, that is, k≤K+4. Setting αi0*=t0, a piecewise constant control function u(·) approximating the optimal solution is constructed based on these nonzero coefficients as follows [46]:(5.41)u(t)={(uPij,uRij)t∈[∑h=0j-1αih*,∑h=0jαih*)0otherwise,j=1,2,…,k,whereuPijand uRijare, respectively, the 8th and the 9th components of grid point yij.


## 6. Numerical Results

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


### 6.1. Continuous Solutions

As we discussed in Section4.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 t0, the time at which the CD4+ T-cell count falls below 350 cells/μL, and treatment was simulated for 750 days. Furthermore, control variables uP and uR are bounded by a1=0, b1=0.7, a2=0, and b2=9×10-10, and the following initial condition is used (see [7]):(6.1)x(0)=103cellsμL-1,v(0)=104copiesmL-1,y(0)=0cellsμL-1,w(0)=10-3cellsmL-1,z(0)=10-7cellsμL-1,r(0)=2×10-7mLcopies-1day-1.We ran simulations with five different values of weight factor w. The corresponding control functions are presented in Figure 2. As we increase w3 and w4, thereby increasing the cost of therapy, the control functions are decreasing. In Figure 8, we plot only the response of the HIV model to the obtained control function corresponding to weights w1=1, w2=0.1, w3=21000, and w4=1021 for brevity. Since the magnitude of CD4+ T cell, and virus 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.Continuous solutions withw1=1 and w2=0.1 and different weights w3 and w4, on control cost: w3=25000 and w4=2×1021 (the x shape), w3=24000 and w4=1.5×1021 (the dashed dotted line), w3=19000 and w4=1021 (the solid line), w3=21000 and w4=1021 (the dotted line), w3=18000 and w4=8×1020 (the dashed line).
(a)(b)Figure 3
Convergence of the gradient method.


### 6.2. 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 factorw. 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 w1=1, w2=0.1, w3=1800, and w4=5×105 for brevity.Figure 4
STI solution withw1=1, w2=0.1, w3=2000, and w4=105.Figure 5
STI solution withw1=1, w2=0.1, w3=1800, and w4=5×105.Figure 6
STI solution withw1=1, w2=0.1, w3=1500 and w4=5×105.


### 6.3. Piecewise Constant Solutions

In our implementation, we chooseK1=12 number of functions from the set C′(B) as
(6.2)tiξj,i∈{0,1,},j∈{1,…,6}.Furthermore, we setK2=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 U1 in the construction of yjs 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 mis and Mis 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, λ=mini=1,…,4{μi(Ii(u))}=0.3867. Figure 7 shows the resulting piecewise constant control pair. We solved the HIV model (2.1)–(2.6) with no treatment (uP=uR≡0) and with fully efficacious treatment (uP≡0.7,uR≡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.Table 2
Considered ranges for states and controls and corresponding number of partitions.
StateRangeNumber of partitionsξ1 (Healthy CD4+T cells)A1=[200,1000]n1=5ξ2 (Infected CD4+T cells)A2=[5,30]n2=3ξ3 (CTLp)A3=[0,1.6]n3=2ξ4 (CTLe)A4=[0,1.3]n4=2ξ5 (Viral load)A5=[500,35000]n5=10ξ6 (aggressiveness of the virus)A6=[0,2×10-7]n6=1uP (PI Control)U1=[0,0.7]s1=3uR (RTI Control)U2=[0.9×10-10]s2=3Table 3
The values of the objective functionals and corresponding values of the membership functions.
I1(u)/μ1(I1(u))I2(u)/μ2(I2(u))I3(u)/μ3(I3(u))I4(u)/μ4(I4(u))Continuous control2.4550×105/0.3478-1.3592×107/0.7513-37.5238/0.8979-3.4147×10-16/0.4379Piecewise constant control3.0929×105/0.5022-2.5882×107/0.4837-138.3152/0.6236-3.7261×10-16/0.3867uP=0.7,uR=9×10-10M1=5.1496×105/1M2=-2.1679×106/1m3=-367.5000/0m4=-6.0750×10-16/0uP=0,uR=0m1=1.0181×105/0m2=-4.8101×107/0M3=0/1M4=0/1Figure 7
The piecewise constant control pairuP* and uR*.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 (uP=uR≡0) (the × shape), and with fully efficacious treatment (uP≡0.7,uR≡9×10-10) (the dashed dotted line).
(a)(b)(c)(d)From Figures8(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. 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 8(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.


## 6.1. Continuous Solutions

As we discussed in Section4.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 t0, the time at which the CD4+ T-cell count falls below 350 cells/μL, and treatment was simulated for 750 days. Furthermore, control variables uP and uR are bounded by a1=0, b1=0.7, a2=0, and b2=9×10-10, and the following initial condition is used (see [7]):(6.1)x(0)=103cellsμL-1,v(0)=104copiesmL-1,y(0)=0cellsμL-1,w(0)=10-3cellsmL-1,z(0)=10-7cellsμL-1,r(0)=2×10-7mLcopies-1day-1.We ran simulations with five different values of weight factor w. The corresponding control functions are presented in Figure 2. As we increase w3 and w4, thereby increasing the cost of therapy, the control functions are decreasing. In Figure 8, we plot only the response of the HIV model to the obtained control function corresponding to weights w1=1, w2=0.1, w3=21000, and w4=1021 for brevity. Since the magnitude of CD4+ T cell, and virus 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.Continuous solutions withw1=1 and w2=0.1 and different weights w3 and w4, on control cost: w3=25000 and w4=2×1021 (the x shape), w3=24000 and w4=1.5×1021 (the dashed dotted line), w3=19000 and w4=1021 (the solid line), w3=21000 and w4=1021 (the dotted line), w3=18000 and w4=8×1020 (the dashed line).
(a)(b)Figure 3
Convergence of the gradient method.


## 6.2. 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 factorw. 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 w1=1, w2=0.1, w3=1800, and w4=5×105 for brevity.Figure 4
STI solution withw1=1, w2=0.1, w3=2000, and w4=105.Figure 5
STI solution withw1=1, w2=0.1, w3=1800, and w4=5×105.Figure 6
STI solution withw1=1, w2=0.1, w3=1500 and w4=5×105.


## 6.3. Piecewise Constant Solutions

In our implementation, we chooseK1=12 number of functions from the set C′(B) as
(6.2)tiξj,i∈{0,1,},j∈{1,…,6}.Furthermore, we setK2=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 U1 in the construction of yjs 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 mis and Mis 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, λ=mini=1,…,4{μi(Ii(u))}=0.3867. Figure 7 shows the resulting piecewise constant control pair. We solved the HIV model (2.1)–(2.6) with no treatment (uP=uR≡0) and with fully efficacious treatment (uP≡0.7,uR≡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.Table 2
Considered ranges for states and controls and corresponding number of partitions.
StateRangeNumber of partitionsξ1 (Healthy CD4+T cells)A1=[200,1000]n1=5ξ2 (Infected CD4+T cells)A2=[5,30]n2=3ξ3 (CTLp)A3=[0,1.6]n3=2ξ4 (CTLe)A4=[0,1.3]n4=2ξ5 (Viral load)A5=[500,35000]n5=10ξ6 (aggressiveness of the virus)A6=[0,2×10-7]n6=1uP (PI Control)U1=[0,0.7]s1=3uR (RTI Control)U2=[0.9×10-10]s2=3Table 3
The values of the objective functionals and corresponding values of the membership functions.
I1(u)/μ1(I1(u))I2(u)/μ2(I2(u))I3(u)/μ3(I3(u))I4(u)/μ4(I4(u))Continuous control2.4550×105/0.3478-1.3592×107/0.7513-37.5238/0.8979-3.4147×10-16/0.4379Piecewise constant control3.0929×105/0.5022-2.5882×107/0.4837-138.3152/0.6236-3.7261×10-16/0.3867uP=0.7,uR=9×10-10M1=5.1496×105/1M2=-2.1679×106/1m3=-367.5000/0m4=-6.0750×10-16/0uP=0,uR=0m1=1.0181×105/0m2=-4.8101×107/0M3=0/1M4=0/1Figure 7
The piecewise constant control pairuP* and uR*.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 (uP=uR≡0) (the × shape), and with fully efficacious treatment (uP≡0.7,uR≡9×10-10) (the dashed dotted line).
(a)(b)(c)(d)From Figures8(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. 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 8(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.


## 7. 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.



---

*Source: 568315-2011-01-24.xml*