Computational Methods in Systems Biology

Data-Driven Robust Control for Type 1 Diabetes 219

Q 2 and in turn, BG concentrationG. Plasma insulin levelsIdirectly increase
x 1 ,x 2 ,x 3. Uncertainty parametersMMt(active muscular mass) andO2t(target
workload in terms of oxygen consumption) increasex 1 ,x 2 ,x 3 indirectly, through
state variablesUA(mg/min) andO2m(unitless), not shown in the figure. They
characterize physical activity and describe, respectively, the glucose uptake due
to active muscular tissue and the actual percentage of maximum oxygen con-
sumption.

Initial Conditions:The initial state of the system is derived at a steady-state
BG level of 7.8 mmol/L [ 31 ], assuming no meal and exercise. We use a non-
linear equation solver (MATLAB’sfsolve) to findx(0) and the basal insulin
level ̄ιsuch thatx ̇(0) =F

(

x(0), ̄ι,u^0

)

= 0 (see Eq. 1 ), where the uncertainty
parametersu^0 are given byD^0 G=0,MM^0 =0andO2^0 = 8 (oxygen consump-
tion at rest). Following [ 13 ], we further assess the physiologic feasibility of the
initial conditions by checking that: (1) in absence of insulin, steady-state BG
is above 300 mg/dL, and (2) delivery of high-dose insulin (15 U/h) results in a
steady-state BG below 100 mg/dL.

4 Robust MPC

Since we want to optimize the BG profile against worst-case realizations of the
uncertainty parameters, at each time stept, the robust MPC computes the
insulin infusionιtas the solution of the following non-linear minimax optimiza-
tion problem:

min

ιt,...,ιt+Nc−^1

max

ut,...,ut+Np−^1

∑Np

k=1

d( ̃x(t+k)) +β·

N∑c− 1

k=0

(Διt+k)^2 (3)

subject to:ιt+k∈Dι (k=0,...,Nc−1) (4)

ιt+k= ̄ι (k=Nc,...,Np−1) (5)

ut+k∈Ut+k (k=0,...,Np−1) (6)

̃x(t)=ˆx(t) (7)

̃x ̇(t+k)=F(x ̃(t+k),ιt+k,ut+k)(k=0,...,Np−1) (8)

whereNcandNpare the control and prediction horizon (in minutes), respec-
tively; constraint ( 4 ) states that the control inputιmust belong to some set
Dιof admissible insulin infusion rates; through ( 5 ), we impose thatιis fixed to
the basal insulin rate ̄ιoutside the control horizon; ( 6 ) states that, at any time
pointt+kin the prediction horizon, uncertainty parametersut+kmust belong
to the corresponding uncertainty setsUt+k; constraint ( 7 ) and ( 8 ) restrict how
the robust MPC computes the predicted state vector ̃x: for the initial state, it
uses the estimated plant state at timet,xˆ(t), while following states are predicted
using the same plant model (see Eq. 1 ). We set control and prediction horizons
toNc= 100 min andNp= 150 min, respectively, as opposed to [ 28 ] where