Feedback Control of HIV Antiviral Therapy with Long Measurement

23 Pages·2011·0.69 MB·English

Feedback Control of HIV Antiviral Therapy with Long
Measurement Time
H. T. Banks, Taesoo Jang and Hee-Dae Kwon ∗
Center for Research in Scientific Computation
Center for Quantitative Sciences in Biomedicine
North Carolina State University
Raleigh, NC, 27695-8212
and
Department of Mathematics
Inha University
Yonghyundong, Namgu
Incheon 402-751, Republic of Korea
January 2, 2011
Abstract
InthispresentationweapplyarecedinghorizonobservertoanHIVfeedbackcontrolproblem
in order to derive optimal treatments of HIV progression and/or optimal structured treatment
interruptions(STI)inantiviraltherapythatincludedrug-freeperiodsofimmune-mediatedcon-
trol of HIV. We use a nonlinear differential equation model which has been well-validated with
clinical patient data and shown to have reliable predictive capabilities. The basic feedback
control problem utilizes a tracking formulation (desired states for viral load and immune effec-
tor calculated from a healthy steady state for the model are tracked). Here we assume only
very realistic observerables of 𝐶𝐷4+ T cell count and viral load. Moreover, we use a second
deterministic optimal tracking problem for state estimation as opposed to stochastic filtering
approaches. In both of the optimization problems (for optimal feedback control and for state
estimator) conjugate gradient algorithms are employed.
AMS Subject Classification: 92C50, 92C42, 49K15, 93C10, 93C15.
Key Words: HIV, feedback control, receding horizon control, optimality systems and computa-
tional methods, state estimator, optimal tracking.
∗[email protected],[email protected], [email protected]
1
1 Introduction
A number of researchers have used a control theoretic approach to derive optimal drug adminis-
tration schemes for HIV infections. Open loop controls have been discussed in [1, 3, 14, 16] and
feedbackcontrolhasbeenstudiedin[5,8,9,24]. Moreover, thewaytoderiveanoptimalon-offtype
of treatment, which is also known as a structured treatment interruption (STI), has been presented
in[1,3]. OptimalfeedbackcontrolbasedonthestatedependentRiccati equation(SDRE)approach
for HIV infection is considered in [8]. The authors in [8] also consider the problem of designing a
state estimator to beusedwiththenonlinearfeedback controllaws sinceonly partialmeasurements
of thestate areavailable inmany practical problems. However, inthateffort themeasurementtime
step is too short, e.g., 15 seconds, to be useful in applications where measurement frequency is typ-
ically on the order of weeks to months. This difficulty is somewhat alleviated in [6] where it is
shown that filters are available which can be tuned to allow such infrequent measurements.
Our approach to the problem of obtaining an estimate of the state of nonlinear systems is
optimization-based. This approach is called a moving horizon observer, receding horizon observer,
or model predictive observer. Moving horizon observers for linear, time-invariant, discrete-time
systems are discussed in [19, 20, 21]. The authors derive sufficient conditions for the stability
of moving horizon state estimation with a linear model subject to constraints on the estimate.
The authors of the papers [4, 13, 18, 22, 25] consider moving horizon observers for nonlinear
systems. The stability of observers is verified under the assumption that there exists a solution of
the corresponding global optimization problem and that a uniform reconstructability condition is
satisfied[18]. Inoureffort, weapplyamovinghorizonobservertoanHIVfeedback control problem
in order to derive optimal treatments for HIV infections and/or optimal structured treatment
interruptions(STI)inantiviraltherapiesthatincludedrug-freeperiodsofimmune-mediatedcontrol
of HIV. The idea of using a moving horizon observer for combined drug/immune response control
of HIV infection is not new. Indeed earlier efforts include [12, 24, 26, 27]. Here, as in [6, 11, 12], we
use the nonlinear model of [2] which was was well-validated with clinical patient data and shown to
have reliable predictive capabilities. However, in [6, 12] the authors assume observations of 𝐶𝐷4+
T cell count, viral load, and immune effector (the latter is not generally available since this is
not a physical compartment!) while here we assume only the more realistic (see the data in [2])
observerables of 𝐶𝐷4+ T cell count and viral load. Moreover, the approaches to state estimation
in[11,12]involve asymptoticfilters (e.g., extendedKalman filtersandvariations) whilehereweuse
a second deterministic optimal tracking problem for state estimation. The basic feedback control
problem here utilizes a tracking formulation (desired states for viral load and immune effector
calculated from a healthy steady state for the model are tracked). In both of the optimization
problems (for optimal feedback control and for state estimator) conjugate gradient algorithms are
employed.
2
2 HIV Model
The system of ordinary differential equations (validated in [2]) describing the infection dynamics is
given by
𝑇˙ = 𝜆 −𝑑 𝑇 −(1−𝜖 )𝑘 𝑉 𝑇
1 1 1 1 1 1 𝐼 1
𝑇˙ = 𝜆 −𝑑 𝑇 −(1−𝑓𝜖 )𝑘 𝑉 𝑇
2 2 2 2 1 2 𝐼 2
𝑇˙∗ = (1−𝜖 )𝑘 𝑉 𝑇 −𝛿𝑇∗−𝑚 𝐸𝑇∗
1 1 1 𝐼 1 1 1 1
𝑇˙∗ = (1−𝑓𝜖 )𝑘 𝑉 𝑇 −𝛿𝑇∗−𝑚 𝐸𝑇∗
2 1 2 𝐼 2 2 2 2 (2.1)
𝑉˙ = (1−𝜖 )𝑁 𝛿(𝑇∗ +𝑇∗)−[𝑐+(1−𝜖 )𝜌 𝑘 𝑇 +(1−𝑓𝜖 )𝜌 𝑘 𝑇 ]𝑉
𝐼 2 𝑇 1 2 1 1 1 1 1 2 2 2 𝐼
𝑉˙ = 𝜖 𝑁 𝛿(𝑇∗ +𝑇∗)−𝑐𝑉
𝑁𝐼 2 𝑇 1 2 𝑁𝐼
𝑇∗+𝑇∗ 𝑇∗+𝑇∗
𝐸˙ = 𝜆 +𝑏 1 2 𝐸 −𝑑 1 2 𝐸−𝛿 𝐸.
𝐸 𝐸𝑇∗+𝑇∗+𝐾 𝐸𝑇∗+𝑇∗+𝐾 𝐸
1 2 𝑏 1 2 𝑑
In the model (2.1), the state variables (densities) are: 𝑇 , the uninfected 𝐶𝐷4+ T-cells; 𝑇 ,
1 2
the uninfected target cells of second kind; 𝑇∗, the infected T-cells; 𝑇∗, the infected target cells of
1 2
second kind; 𝑉 , the infectious virus; 𝑉 , the non infectious virus; and 𝐸, the immune effectors.
𝐼 𝑁𝐼
As in [2], we do not give precise biological definitions for the target cells of second kind and the im-
mune effectors. They could, for example, be related to macrophages and cytotoxic T-lymphocytes,
respectively.
In the model, we have two controllers 𝜖 and 𝜖 . The first control term 𝜖 represents the
1 2 1
efficacy of reverse transcriptase inhibitors (RTI) that block new infections. The second control
term 𝜖 represents the efficacy of protease inhibitors (PI). Thus, the productivity, 𝑁 , is reduced
2 𝑇
to (1−𝜖 )𝑁 . For a more detailed description of the variables and rationale for the model (2.1)
2 𝑇
we refer the reader to the articles [1, 2, 3, 7, 8]. Table 1 contains the values of parameters we used
in our computations.
parameter value unit parameter value unit
𝜆 10.0 cells 𝜆 31.98×10−3 cells
1 mm3⋅day 2 mm3⋅day
𝑑 0.01 1 𝑑 0.01 1
1 day 2 day
𝑘 8.0×10−4 mm3 𝑘 0.1 mm3
1 virions⋅day 2 virions⋅day
𝑚 0.01 mm3 𝑚 0.01 mm3
1 cells⋅day 2 cells⋅day
𝜌 1 virions 𝜌 1 virions
1 cells 2 cells
𝛿 0.7 1 𝑐 13.0 1
day day
𝑓 0.34 – 𝑁 100.0 virions
𝑇 cells
𝜆 1.0×10−3 cells 𝛿 0.1 1
𝐸 mm3⋅day 𝐸 day
𝑏 0.3 1 𝑑 0.25 1
𝐸 day 𝐸 day
𝐾 0.1 cells 𝐾 0.5 cells
𝑏 mm3 𝑑 mm3
Table 1: The values of the parameters in the HIV model.
In order to simplify our subsequent discussions, we introduce the notation 𝑥 and 𝑢 to denote
3
the state and control vectors, respectively. Thus we define
⎛ ⎞
𝑇
1
⎜ 𝑇 ⎟
⎜ 2 ⎟
⎜ 𝑇∗ ⎟ ( )
⎜ 1 ⎟ 𝜖
𝑥 =⎜ 𝑇∗ ⎟ and 𝑢 = 1 .
⎜ 2 ⎟ 𝜖
⎜ 𝑉 ⎟ 2
⎜ 𝐼 ⎟
⎝ ⎠
𝑉
𝑁𝐼
𝐸
With the above notation, the HIV model (2.1) can be expressed in the generic form
𝑥˙ = 𝑓(𝑥)+𝐵(𝑥)𝑢, (2.2)
where the precise representation for 𝑓(𝑥) and 𝐵(𝑥) will be described later.
For feedback control we need current knowledge on the state of the system. In our effort here
we assume that partial state observations (𝑇 +𝑇∗,𝑉 +𝑉 ) are available. This is representative
1 1 𝐼 𝑁𝐼
of the type of clinical data widely discussed in the literature (see for example, [1, 2, 7]). Hence, the
output or observation takes the form
( ) ( )
𝑧 1 0 1 0 0 0 0
𝑧 = 1 = 𝑥 = 𝐶𝑥, (2.3)
𝑧 0 0 0 0 1 1 0
2
where 𝑧 and 𝑧 represent the total CD4+ densities and the total viral load densities, respectively.
1 2
4
3 State and State Estimator
With the parameters given in Table 1 and the assumption of no medication (𝜖 = 𝜖 = 0), the
1 2
model exhibits several steady states. These are described and analyzed in [1, 2, 3] with the non
infectious virus 𝑉 state always being zero in these steady states. We are particularly interested
𝑁𝐼
in the so-called “healthy” steady state given by
cells cells cells cells
𝑇˜ = 967.839 , 𝑇˜ = 0.621 , 𝑇˜∗ =0.076 , 𝑇˜∗ = 0.006 ,
1 mm3 2 mm3 1 mm3 2 mm3
(3.1)
virions virions cells
𝑉˜ = 0.415 , 𝑉˜ = 0.0 , and 𝐸˜ = 353.108 ,
𝐼 mm3 𝑁𝐼 mm3 mm3
which was also shown to be locally asymptotically stable. This means that after a sufficiently small
perturbation from (3.1), the trajectory of the state 𝑥 returns to the stable equilibrium (3.1). It is
noted thatthis stable equilibriumexhibitsa strongimmuneresponse,low viralload and reasonably
large target cell (𝑇 ) count densities.
1
The approach to derive the feedback controllers for our nonlinear system is as follows. Let
{𝑡 } be a sequence of measurement times, not necessarily periodic, such that 𝑡 → ∞ as 𝑖 → ∞.
𝑖 𝑖∈ℕ 𝑖
Let {𝑧¯} be the corresponding data (observations) to be used in updating our controls. First, we
𝑖∈ℕ
formulate the optimal tracking control problems of finding an effective HIV therapeutic strategies
over the subinterval [𝑡 , 𝑡 ]. To this end, we define the objective tracking functional
𝑖 𝑗
∫
𝐽(𝑥,𝑢) = 1 𝑡𝑗{𝑃(𝑉 −𝑉˜ )2+𝑄(𝐸 −𝐸˜)2+(𝜖 /𝜖max)2+(𝜖 /𝜖max)2}𝑑𝑡, (3.2)
2 𝐼 𝐼 1 1 2 2
𝑡
𝑖
where 0 ≤ 𝜖 ≤ 𝜖max, 0 ≤ 𝜖 ≤ 𝜖max and 𝑉˜ ,𝐸˜ are the desired (healthy states) infectious viral load
1 1 2 2 𝐼
and immune factor, respectively, given above. We introduce the following two vectors
( ) ( )
0 𝜖max
𝑢 = and 𝑢¯ = 1 .
0 𝜖max
2
The optimal tracking control problem over the subinterval [𝑡 , 𝑡 ] is to find a dynamic multidrug
𝑖 𝑗
therapy 𝑢(𝑡) satisfying
min 𝐽(𝑥(𝑡),𝑢(𝑡)) (3.3)
𝑢≤𝑢(𝑡)≤𝑢¯
subject to the state equation given by (2.2) with initial condition 𝑥 over the subinterval [𝑡 , 𝑡 ].
𝑖 𝑖 𝑗
A standard optimality system methodology as in [3, 8] was used to computationally solve these
problems.
After finding the optimal control, we just consider the optimal control function 𝑢∗, over the
1
first subinterval [𝑡 , 𝑡 ]. So then we have the state equation and the observation which are given by
0 1
𝑥˙ = 𝑓(𝑥)+𝐵(𝑥)𝑢∗
1 (3.4)
𝑧 = 𝐶𝑥
with a given initial value 𝑥(𝑡 ) = 𝑥 over the interval [𝑡 , 𝑡 ]. Here 𝐶 is defined by (2.3).
0 0 0 1
Next, given measurements 𝑧¯ = (𝑧 ,𝑧 )𝑇 at time 𝑡 for 𝑇(𝑡 )+ 𝑇∗(𝑡 ) and 𝑉 (𝑡 )+ 𝑉 (𝑡 ),
1 1 2 1 1 1 1 𝐼 1 𝑁𝐼 1
respectively, to useas “feedback”, we would like to finda continuous optimal control, or an optimal
5
x
x
x0=x(t0) xe x1=xe(t1)
t t t t
0 1 2 3
Figure 1: 𝑥 is state and 𝑥 is state estimator
𝑒
STI, 𝑢∗ over the second subinterval [𝑡 , 𝑡 ]. However, we don’t have sufficient observation informa-
2 1 2
tion for the initial condition, 𝑥 , of the second subinterval since we just have partial measurements
1
𝑧¯ of the state at 𝑡 = 𝑡 such that 𝑧¯ ≈ 𝐶𝑥(𝑡 ) in (3.4). So we design a state estimator based on
1 1 1 1
optimization in order to obtain full knowledge for the initial condition of the second subinterval.
Our estimator 𝑥 has the form
𝑒
𝑥˙ = 𝑓(𝑥 )+𝐵(𝑥 )𝑢∗+𝑤 (3.5)
𝑒 𝑒 𝑒 1
withagiveninitialvalue𝑥 (𝑡 )= 𝑥(𝑡 )and𝑤(𝑡) = (𝑤 (𝑡), 𝑤 (𝑡), 𝑤 (𝑡), 𝑤 (𝑡), 𝑤 (𝑡), 𝑤 (𝑡), 𝑤 (𝑡))𝑇.
𝑒 0 0 1 2 3 4 5 6 7
Here 𝑥 = (𝑇ˆ , 𝑇ˆ , 𝑇ˆ∗, 𝑇ˆ∗, 𝑉ˆ , 𝑉ˆ , 𝐸ˆ)𝑇. The time-dependent functions 𝑤 (𝑡), 𝑖 = 1,2,⋅⋅⋅ ,7, are
𝑒 1 2 1 2 𝐼 𝑁𝐼 𝑖
obtained by solving a new optimal control problem
min1 ∫ 𝑡1∑7 𝑤2(𝑡)𝑑𝑡+𝑃ˆ(𝑇ˆ1(𝑡1)+𝑇ˆ1∗(𝑡1) −1)2+𝑄ˆ(𝑉ˆ𝐼(𝑡1)+𝑉ˆ𝑁𝐼(𝑡1) −1)2, (3.6)
𝑤 2 𝑡 𝑗 𝑧1 𝑧2
0 𝑗=1
where 𝑧 and 𝑧 represent the total CD4+ count densities and the total viral load densities of
1 2 ∑
a patient at 𝑡 = 𝑡 , respectively. Note that because of (3.5), the sum 7 𝑤2(𝑡) in the cost
1 𝑗=1 𝑗
functional is the same as ∣𝑥˙ (𝑡)−𝑓(𝑥 (𝑡))−𝐵(𝑥 (𝑡))𝑢∗(𝑡)∣2 which provides our state estimator 𝑥
𝑒 𝑒 𝑒 1 𝑒
on the interval [𝑡 ,𝑡 ]. In our simulations, we use
0 1
𝑧 = 𝑅ˆ(𝑇 (𝑡 )+𝑇∗(𝑡 )) and 𝑧 = 𝑅ˆ(𝑉 (𝑡 )+𝑉 (𝑡 )) (3.7)
1 1 1 1 1 2 𝐼 1 𝑁𝐼 1
asourgiven“data”. Nowweconsidertheoriginaloptimalcontrolprobleminthesecondsubinterval
[𝑡 , 𝑡 ]withtheinitialcondition,𝑥(𝑡 ) = 𝑥 (𝑡 )(Figure1). Werepeatthisprocesstofindanoptimal
1 2 1 𝑒 1
control and state estimator over the entire interval of interest.
6
4 Numerical simulations
Because the optimality system is atwo pointboundaryvalueproblem, weuseaconjugate gradient-
typeiterative methodto solve thesystem. TheConjugate Gradient Algorithm is apowerfulscheme
forsolvinglargescaleminimization problems[10,15,17,23]. Inparticular,weneedmanyiterations
tosolvethelongmeasurementtimeproblembecausewehavetwodifferentoptimalcontrolproblems
in each subinterval. In each optimal control problem, the state system with initial conditions is
solved forwardintimeby usingan initial guessof thecontrols andtheadjointsystem withterminal
conditions is solved backward in time. The controls are updated in each iteration by using the
conjugate gradient algorithm.
We now give an algorithm and provide a detailed description of our conjugate gradient method
for a given Lagrangian 𝐿, as follows:
[Initialization] 𝜖0 = initial guess of controls
( ) ∂𝐿 ( )
𝑑0 = −𝑔0 = −𝑔 𝜖0 := − 𝜖0
∂𝜖
[Choose step size] 𝛼𝑖 := 𝐴𝑟𝑔min𝐽(𝜖𝑖+𝛼𝑑𝑖) (4.8)
𝛼
[Update] 𝜖𝑖+1 = 𝜖𝑖 +𝛼𝑖𝑑𝑖
( ) ∂𝐿 ( )
𝑔𝑖+1 = 𝑔 𝜖𝑖+1 := 𝜖𝑖+1 (4.9)
∂𝜖
(𝑔𝑖+1,𝑔𝑖+1−𝑔𝑖)
𝛽𝑖 =
(𝑔𝑖,𝑔𝑖)
𝑑𝑖+1 = −𝑔𝑖+1+𝛽𝑖𝑑𝑖, (4.10)
where
∫
𝑇
(𝑔𝑖,𝑔𝑗)= 𝑔𝑖(𝑡)𝑔𝑗(𝑡)𝑑𝑡.
0
Let 𝜖𝑖 be the 𝑖th approximation to the optimal control 𝜖∗. The corresponding gradient 𝑔(𝜖𝑖) is
computed by solving the state equation forward with 𝜖 = 𝜖𝑖, solving theadjoint equation backward,
and then computing 𝑔(𝜖𝑖) from (4.9). The step size 𝛼 is determined by solving the one dimensional
minimization problem (4.8). For more information on this iterative method, we refer the interested
reader to [17].
The parameters used in solving the optimality system are summarized in Table 1. We set the
maximum efficacy of drug at 𝜖max = 0.7 and 𝜖max = 0.3 and ∣𝑤 ∣≤ 1 for 𝑖= 1,2,⋅⋅⋅ ,6. Because of
1 2 𝑖
differences in the magnitude of the cost of drug treatment, the viral load, and immune effectors in
the firstobjective functional (3.2), the weight constants 𝑃 = 1 and 𝑄 = 10 were selected to balance
the units. Moreover, we also choose 𝑃ˆ = 𝑄ˆ = 105 to place more weight in the tracking terms
in the second objective functional (3.6). In the numerical runs, we further varied the constant,
𝑅ˆ = 0.8, 0.9, 1.1, and 1.2 in the simulated measurement data (3.7). We assume that medication
starts immediately after HIV infection so the initial values for the state variables are given by the
7
acute state. For this, we simulate by introducing one virus particle and very low levels of infected
T cell densities as follows:
𝑇 (0) = 103, 𝑇 (0) = 3.198, 𝑇∗(0) = 10−7, 𝑇∗(0) = 10−7, 𝑉(0) =0.001 and 𝐸(0) = 0.01.
1 2 1 2
In all examples, we solve the first optimal control problem to obtain the control function 𝑢∗
over 60 days (i.e., 𝑡 −𝑡 = 60). We consider three data sampling frequency cases (15 days, 20 days
𝑗 𝑖
and 30 days) for subinterval in the second optimal control problem for the state estimator. We
assume these days are the measurement times to update our optimal treatment strategies 𝑢∗. The
following is a detailed description of our procedure for 15 day measurement times:
1. Choose the initial condition 𝑥 (0) = (103, 3.198, 10−7, 10−7, 0.001, 0.01);
𝑖𝑛𝑖
∙ Solve the first optimal control problem for the state, 𝑥 with the initial condition 𝑥 (0)
𝑖𝑛𝑖
over [0, 60],
∙ Take the optimal control functions 𝑢∗ over [0, 15],
1
∙ Solve the second optimal control problem for the state estimator 𝑥 with the initial
𝑒
condition 𝑥 (0) =𝑥 (0) over [0, 15];
𝑒 𝑖𝑛𝑖
2. Determine an initial condition, 𝑥 (15), for the state 𝑥 in the second subinterval using the
𝑖𝑛𝑖
state estimator 𝑥 , i.e., 𝑥 (15) = 𝑥 (15);
𝑒 𝑖𝑛𝑖 𝑒
∙ Solve the first optimal control problem for the state 𝑥 with the initial condition 𝑥 (15)
𝑖𝑛𝑖
over [15, 75],
∙ Take the optimal control functions 𝑢∗ over [15, 30],
2
∙ Solve the second optimal control problem for the state estimator 𝑥 with the initial
𝑒
condition 𝑥 (15) = 𝑥 (15) over [15, 30];
𝑒 𝑖𝑛𝑖
3. Determine an initial condition, 𝑥 (30), for the state 𝑥 in the third subinterval using the
𝑖𝑛𝑖
state estimator 𝑥 , i.e., 𝑥 (30) = 𝑥 (30);
𝑒 𝑖𝑛𝑖 𝑒
∙ Solve the first optimal control problem for the state 𝑥 with the initial condition 𝑥 (30)
𝑖𝑛𝑖
over [30, 90],
∙ Take the optimal control functions 𝑢∗ over [30, 45],
3
∙ Solve the second optimal control problem for the state estimator 𝑥 with the initial
𝑒
condition 𝑥 (30) = 𝑥 (30) over [30, 45];
𝑒 𝑖𝑛𝑖
4. Continue this process over the entire time interval of interest.
Theoptimalcontrolfunctions,optimalsolutions,andthefunctions𝑤(𝑡)for15daymeasurement
times are depicted in Figures 2-5. The results for 20 day and 30 day measurement times are shown
inFigures6-9andinFigures10-13, respectively. Wealsopresentthecontrolfunctionsandsolutions
(dashed line) with MPC (Model Predictive Control) for comparison purposes. MPC solves a finite
horizonopen-loopoptimalcontrolproblemiteratively. Thecurrentstateissampledatmeasurement
times and a control strategy minimizing a cost is computed for a short time interval. In the case
of MPC, we assumed a best case scenario that we have full knowledge of all the state variables at
each measurement time. That is, the values of the state 𝑥 at each sampling instant are given by
the final values of the state in the previous subinterval.
8
Control RT
0.6
0.4
0.2
0
0 50 100 150 200 250 300 350 400 450 500
Control PIs
0.4
0.3
0.2
0.1
0
0 50 100 150 200 250 300 350 400 450 500
103 T1 T2 1 W1 1 W2
100
0.5 0.5
10−1
0 0
10−2 −0.5 −0.5
102 10−3 −1 −1
0 100 200 300 400 500 0 100 200 300 400 500 0 100 200 300 400 500 0 100 200 300 400 500
T1* T2* W3 W4
1 1
102 100
0.5 0.5
100 10−1
0 0
10−2
10−2 −0.5 −0.5
10−3 −1 −1
0 100 200 300 400 500 0 100 200 300 400 500 0 100 200 300 400 500 0 100 200 300 400 500
Vi VNI W5 W6
1 1
102
102
0.5 0.5
100 100 0 0
10−2 10−2 −0.5 −0.5
−1 −1
0 100 200 300 400 500 0 100 200 300 400 500 0 100 200 300 400 500 0 100 200 300 400 500
E W7
1
102
0.5
100 0
−0.5
10−2 −1
0 50 100 150 200 250 300 350 400 450 500 0 50 100 150 200 250 300 350 400 450 500
Figure 2: Feedback solutions (solid line) and MPC solutions (dashed line) with 𝑅ˆ = 0.8, ∣𝑤 ∣ ≤
𝑖
1, 𝑃 = 1,𝑄 = 10,𝑃ˆ = 105,𝑄ˆ =105. 𝑡 −𝑡 = 15.
𝑖+1 𝑖
9
Control RT
0.6
0.4
0.2
0
0 50 100 150 200 250 300 350 400 450 500
Control PIs
0.4
0.3
0.2
0.1
0
0 50 100 150 200 250 300 350 400 450 500
103 T1 T2 1 W1 1 W2
100
0.5 0.5
10−1
0 0
10−2 −0.5 −0.5
102 10−3 −1 −1
0 100 200 300 400 500 0 100 200 300 400 500 0 100 200 300 400 500 0 100 200 300 400 500
T1* T2* W3 W4
1 1
102 100
0.5 0.5
100 10−1
0 0
10−2
10−2 −0.5 −0.5
10−3 −1 −1
0 100 200 300 400 500 0 100 200 300 400 500 0 100 200 300 400 500 0 100 200 300 400 500
Vi VNI W5 W6
1 1
102
102
0.5 0.5
100 100 0 0
10−2 10−2 −0.5 −0.5
−1 −1
0 100 200 300 400 500 0 100 200 300 400 500 0 100 200 300 400 500 0 100 200 300 400 500
E W7
1
102
0.5
100
0
10−2 −0.5
−1
0 50 100 150 200 250 300 350 400 450 500 0 50 100 150 200 250 300 350 400 450 500
Figure 3: Feedback solutions (solid line) and MPC solutions (dashed line) with 𝑅ˆ = 0.9, ∣𝑤 ∣ ≤
𝑖
1, 𝑃 = 1,𝑄 = 10,𝑃ˆ = 105,𝑄ˆ =105. 𝑡 −𝑡 = 15.
𝑖+1 𝑖
10