Autoregression with Exogenous Input
The Zone-MPC model is based on an ARX model.17,18
The ARX model that serves as the basis for zone-MPC is chosen from a pool of ARX model candidates that differ by the number of past data points utilized (the ARX order). The different ARX models are scored by their prediction abilities over a fixed horizon, and the ARX model that scores the higher prediction R2
, Equation (2)
, over a fixed horizon is then estimated again to improve its prediction abilities.
presents the R2
index, where y RN
is a vector of collected data points,
¯ is the mean of collected data points, and ŷ
is the predicted value.
The R2 index approaches 1 for perfect model predictions and 0 for models that do not predict better than a constant mean value.
Autoregression with exogenous input equations use past input/output discrete records to generate a future discrete prediction. For example, Equation (3) describes the connection between predicted output y^ at instant k and past p and q output and input records, respectively.
In the presence of observed data y RN
, where N
is the number of collected data points, linear regression can be performed to establish regression vector
θ values by minimizing the sum squares of errors between N
data records y
and predicted values ŷ
as it is formulated in Equation (4)
Regression on past output regressors, [α1. . . αp], is defined as autoregression, and the model is referred to as an ARX model.
Determination of ARX model output order, input order, and delay is subject to data validation and complexity considerations. Autoregression with exogenous input model identification is initialized by the researcher’s decision based on reasonable orders and delays. A pool of ARX models is then generated by the different combinatory combinations of orders and delays, and each ARX model is evaluated by its prediction abilities on calibration data and on validation data. Some trade-off can be formulated between model prediction abilities and ARX model complexity, for example, the Akaike information criterion.18
However, the quality of collected data for identification will eventually govern identification.
In order to best simulate an in vivo
experiment, a data-collecting protocol (
) introduced by Dassau and collaborators19
was applied to the Food and Drug Administration (FDA)-accepted University of Virginia (UVa)\University of Padova (Padova) metabolic simulator.20
In the data-collecting protocol, no adjustments to the daily routine were prescribed in day 1. Both days 2 and 3 started at 7 amwith a 25-gram carbohydrate (CHO) breakfast with no insulin bolus, followed by a correction bolus at 9 am; at 1 pma 50-gram CHO lunch was taken together with a correction bolus. A 15-gram CHO snack without a correction bolus was given at 5 pm, and a 75-gram CHO dinner accompanied by an insulin bolus was given at 8 pm. Days 4 and 5 started at 7 amwith a 25-gram CHO breakfast accompanied by an insulin bolus, then a 50-gram lunch was consumed at 12 pm, and an insulin bolus was given at 2 pm. At 8 pma 75-gram CHO dinner accompanied by an insulin bolus was given, and at 10 pman insulin bolus was administered. On the 6th
day the response to a pure bolus from fasting conditions was tested by bolus administration at 9 am.
Figure 1. Proposed protocol that facilitates the separation of meal and insulin effects on blood glucose.19
In this work, data of meal and insulin inputs were mapped through second-order transfer functions [Equation (5)] to overcome a major identification problem: a large gain uncertainty was caused by the opposite effects of meal and insulin that are delivered frequently in close time instants. Losing the ability to distinguish between meal and insulin gains can produce poor models for control. However, mapping input data using two different second-order transfer functions () separated and spread input data, providing better terms to regress for the model. Clinical observations and pharmacokinetic/pharmacodynamic data suggest that, on average, the effect of insulin on blood glucose is observed 30 minutes after injection and that the effect of meals on blood glucose is observed after 20 minutes. This a priori knowledge has been used to design the following second-order transfer functions:
Figure 2. Insulin and meal inputs are mapped through second-order transfer functions and, as a result, are spread and separated. (A) Typical input data collected from T1DM subjects: meals and insulin are assigned as pulses over relatively close discrete time instances. (more ...) Equation (5)
describes transfer functions used to map measured insulin and meal, I
, respectively, into new states Imap
. Transfer functions in the Laplace transform21
provide the engineer with a powerful method for the analytical solution of differential equations. The Laplace transform is often interpreted as a transformation from the time domain to the frequency domain. The four constants (30, 25, 10, and 45) are in units of minutes.
Horizon Prediction Optimization
The ARX model regressors are normally set by a one step shift of output/input data, such that each prediction at time instant k is evaluated by recorded output/input data up to the time instant k – 1. However, predicted output at instant k is highly correlated to the output measurement at instant k – 1, especially when the sampling time is relatively short. This can lead to ARX models that lose their predictive capabilities.
A novel improvement was suggested to parametric identification of the ARX model as follows: Data are divided into groups of the length of ↓ (N / PH) · PH where N is the number of data points, PH is the prediction horizon, and “↓” is the round down operator. Next, a parametric optimization is carried out in an attempt to find the set of parameters that maximize the prediction in each group. Each group is recursively predicted: the first term of each group is predicted by measured data; then the first term is used to predict the second term in each group; and the predicted terms, first and second, are used to predict the third term. This procedure is repeated for all data points in each group, and the ability to predict all groups is used as a cost function for optimization, Equation (6).
Horizon prediction optimization [Equation (6)
] attempts to minimize normalized sum squares of errors between output recorded data, y
, and the recursive prediction in each data group, yrec
) is the modulus after the division of
are mapped insulin and meal data records, respectively. p
, and q2
are the order of output (glycemia), insulin, and meal, respectively, and d1
are insulin and meal delays in sampling time, respectively. αk
, and γk
are regressors of glucose, insulin, and meal, respectively.
The following five constraints were applied to the optimization to prevent models from being unstable, having nonphysical gains, or having inverse responses and so that they may be suitable for predictive control:
- For stability, roots of the following of the characteristic polynomial, zP– α1zP–1– α2zP–2... – α2, are all inside the unit circle.
- Negative insulin gain requirement,
- Positive meal gain requirement,
- No inverse response in insulin, , where is defined by Equation (7).
No inverse response in meal,
Equation (7) describes a Boolean function that results in 0 if provided with a vector that has elements on only one side of the origin (i.e., all positive or all negative). The function results in 1 if provided with a vector that has at least two elements from different sides of the origin.
Zone Model Predictive Control
The different MPC algorithms can be classified into four approaches to specify future process response:10
fixed set point, zone, reference trajectory, and funnel. Using a fixed set point for the future process response can lead to large input adjustments unless settings of the controller are changed in detriment of performance. A zone control is designed to keep the CV in a zone defined by upper and lower boundaries that are usually defined as soft constraints. Some MPC algorithms define a desired response path for CV, called reference trajectory. The reference trajectory usually describes a defined path from the current CV state to a desired set point. The reference trajectory control returns to a fixed set point control when the CV approaches the defined set point. Robust Multivariable Predictive Control Technology (Honeywell Inc., 1995) attempts to keep the CV in a defined zone; however, when the CV is out of the zone, a funnel is defined to bring the CV back into the zone. Kovatchev and colleagues15
presented the concept of control to range for diabetes that optimizes glucose control to a certain glucose control range. The control-to-range concept is focused on CV performances rather than on controller architecture.
Zone-MPC is applied when the specific set point value of a CV is of low relevancy compared to a zone that is defined by upper and lower boundaries. Moreover, in the presence of noise and model mismatch there is no practical value using a fixed set point. A general description of MPC with output objectives defined as a zone can be found in Maciejowski.22
Our zone-MPC is implemented by defining fixed upper and lower bounds (
) as soft constraints by letting optimization weights switch between zero and some final values when predicted CVs are in or out of the desired zone, respectively. Predicted residuals are generally defined as the difference between the CV that is out of the desired zone and the nearest bound.
Figure 3. Illustration of zone-MPC in the context of diabetes; zone-MPC is typically divided into three different zones. The permitted range is the control target and is defined by upper and lower bounds. For example, green dots indicate predicted glycemic values (more ...)
The core of zone-MPC lies in its cost function formulation. Zone-MPC, like any other form of MPC, predicts the future ŷ RP
output by an explicit model using past P
input/output records and future input u RM
moves that need to be optimized. However, instead of driving to a specific fixed set point, the optimization attempts to keep or move the predicted outputs into a zone that is defined by upper and lower bounds.
depicts zone-MPC applied to glycemia regulation. Fixed upper and lower glycemic bounds are predefined. Using a linear difference model, glycemic dynamics are predicted and optimization reduces future glycemic excursions from the zone under constraints and weights defined in its cost function.
The zone-MPC cost function used in the presented work is defined as follows:
where Q and R are constant optimization weights on predicted outputs and future inputs, respectively, and yrange is a superposition of all the predicted outputs states that exceed the permitted range:
collects all predicted points below the lower bound by setting all other predicted values to zero:
collects all predicted points above the upper bound by setting all other predicted values to zero:
Future CM are constrained, set by the ability of the insulin pump to deliver a maximum rate of insulin and the inability to deliver negative insulin values. The objective function [Equation (8)] neglects the control move increments component to enable fast control movements.
Second-order input transfer functions described by Equation (5) are used to generate an artificial input memory in zone-MPC schema to prevent insulin overdosing and, as a result, hypoglycemia. In order to avoid overdelivery of insulin, the evaluation of any sequential insulin delivery must take into consideration the past administered insulin against the length of the insulin action. However, a one-state linear difference model with a relative low order uses output (glycemia) as the main source of past administered input (insulin) “memory.” In the face of model mismatch, noise, or change in the subject’s insulin sensitivity, this may result in under- or overdelivery of insulin. This can be mitigated by adding two additional states for mapped insulin and meal inputs.
As discussed earlier, ARX models are identified using mapped inputs. Insulin and meal measurements and transfer functions used for identification can be embedded into one ARX model with one state (glucose) and two inputs (insulin and meal), as illustrated in Equation (12).
where G, I, and M represent glucose blood concentration, insulin administration, and meals, respectively. αi, β1 i, and β2 i are model coefficients; d1and d2are insulin and meal time delays, respectively; and p, q1, and q2are orders of glucose, insulin, and meal, respectively.
In order to generate a prediction using Equation (12), we need p past measurements of G and d1+ q1and d2+ q2past insulin and meal measurements, respectively. However, while G has infinite memory, which is the outcome of the deterministic difference equation, the inputs memory is restricted to past inputs of a restricted order. For example, if delay in insulin, d1, is equal to two sample times and the order of insulin in the model is equal to three sample times, then the model captured the insulin administration of only five previous sample times. In perfect conditions where there is no model mismatch, no noise, or change in insulin sensitivity, past inputs will be contained into the G dynamics; however, in reality, these conditions are unlikely to be achieved.
An alternative formulation of the model, which provides a full reliable input memory to the system, is described in Equation (13):
where Imapand Mmapare new states representing mapped insulin and meal values, respectively. γi and δi are new additions to the set of coefficients. These new states represent the insulin and meal after being absorbed into the blood.
Equation (13) describes an improved formulation to Equation (12) that uses mapped insulin and meal inputs (Imapand Mmap, respectively) as additional states (). The two new states are evaluated using the two past state records as well as two past input records. Keeping past new state records enables an infinite input memory that is independent of G measurements. The states are updated after each zone-MPC cycle and thereby maintain the influence of all past inputs.
Figure 4. Model prediction structure in zone-MPC includes two parts. First is second-order transfer function that is unvaried over the different subjects. Second is an individualized discrete model that predicts glycemia given mapped input data. The overall prediction (more ...)
Automatic Controller Tuning
Four main tuning parameters are available for MPC. First, the prediction horizon, P, is a fixed integer indicating the number of prediction samples that will be used by the MPC cost function. P should be on the order of the settling time (ST). Choosing P that is too small will cause the loss of a valuable dynamic response. However, choosing P that is too large will reduce the influence of dynamics over steady state. Second, the control horizon, M, is an integer indicative of the CM that are optimized. A large M will create a sluggish control response; however, choosing M = 1 will cause highly aggressive control. The other two parameters are Q and R, as described earlier. When applying zone-MPC, in addition to the four MPC tuning parameters, zone dimensions can be treated as another tuning variable.
The most dominant dynamic property found to influence zone-MPC tuning was the ST of the ARX model. Settling time is defined as the time it takes for a dynamic system to accomplish 95% of a response to a step input. The prediction horizon P is set to correspond to the individual subject model ST. However, when the ratio between P and M is changed [Equation (8)], then the relative cost between output and input deviations is changed as well.
It should be noted that glucose predictions based on ARX models are limited to approximately 3 hours. Hence, to compensate correctly for prolonged settling times, the weight upon the CM (R) is a function of the ST. In this way we also conserve the balance between weights of the output and input in Equation (8). Moreover, model fitness based on 3 hours of prediction is used as additional penalty upon
Equation (14) describes the adjusted cost weight upon the CM. The ST is defined by sampling times, and FIT3 h is assumed to have a value greater than 0 and lower than 1 and is calculated by Equation (15):
RN − (p + PH)
is a vector containing the collection of ARX model predictions over 3 hours, PH
is the number of sampling times contained in 3 hours, and p
is the order of autoregressors. N
is the total number of data values, y
is the vector of data values, and
is the mean of data values. Q
and M were used as fixed values of 1 and 5, respectively.
Nine control experiments were conducted on 10 in silico adult subjects following a three-meal scenario of 75, 75, and 50 grams of CHO at 7 am, 1 pm, and 8 pm, respectively, using the FDA-accepted UVa\Padova metabolic simulator. In all MPC experiments, the weight upon predicted outputs, Q, is set to 1, while the weight upon future CM, R, is set automatically by Equation (14).
List of Experiments
The zone-MPC model was compared to patient-direct, open-loop control in order to have an objective comparison to other researchers’ work.23,24
Nine experiments were conducted as follows:
Experiment 1: Built-in open-loop preadjusted treatment is applied with nominal meals values.
Experiment 2: Zone-MPC bounds are set between 80 and 140 mg/dl and meals are unannounced.
Experiment 3: Zone-MPC bounds are set between 100 and 120 mg/dl and meals are unannounced.
Experiment 4: Model predictive control with set point at 110 mg/dl and meals are unannounced.
Experiment 5: Zone-MPC bounds are set between 80 and 140 mg/dl and nominal meals are announced.
Experiment 6: Zone-MPC bounds are set between 100 and 120 mg/dl and nominal meals are announced.
Experiment 7: Model predictive control with set point at 110 mg/dl and nominal meals are announced.
Experiment 8: Built-in open-loop preadjusted treatment is tested with meals announced with –40% mismatch with consumed meals value.
Experiment 9: Zone-MPC bounds are set between 80 and 140 mg/dl with meals announced with –40% mismatch with consumed meals value.