Interval Safety Layer Coupled With an Impulsive MPC for Artificial Pancreas to Handle Intrapatient Variability

The aim of control strategies for artificial pancreas systems is to calculate the insulin doses required by a subject with type 1 diabetes to regulate blood glucose levels by reducing hyperglycemia and avoiding the induction of hypoglycemia. Several control formulations developed for this end involve a safety constraint given by the insulin on board (IOB) estimation. This constraint has the purpose of reducing hypoglycemic episodes caused by insulin stacking. However, intrapatient variability constantly changes the patient’s response to insulin, and thus, an adaptive method is required to restrict the control action according to the current situation of the subject. In this work, the control action computed by an impulsive model predictive controller is modulated with a safety layer to satisfy an adaptive IOB constraint. This constraint is established with two main steps. First, upper and lower IOB bounds are generated with an interval model that accounts for parameter uncertainty, and thus, define the possible system responses. Second, the constraint is selected according to the current value of glycemia, an estimation of the plant-model mismatch, and their corresponding first and second time derivatives to anticipate the changes of both glucose levels and physiological variations. With this strategy satisfactory results were obtained in an adult cohort where random circadian variability and sensor noise were considered. A 92% time in normoglycemia was obtained, representing an increase of time in range compared to previous MPC strategies, and a reduction of time in hypoglycemia to 0% was achieved without dangerously increasing the time in hyperglycemia.

1 Introduction

Managing type 1 diabetes (T1D) has proven to be challenging. People with T1D need exogenous insulin to regulate their blood glucose (BG) levels. The therapy required involves a risk of severe hypoglycemia, with all its consequences, if the insulin dose is too high. For this reason, the therapeutic goal is to minimize the number of hypoglycemic episodes and maximize time in the healthy glycemic range, known as normoglycemia zone. Accordingly, automated insulin delivery, also known as the artificial pancreas (AP), has emerged as the best solution to automatically modulate insulin and remove the threat posed by hypoglycemia (1, 2).

The control strategies more extended in the literature for AP systems are based on model predictive control (MPC) (1), proportional-integral-derivative (PID) (3), and fuzzy logic (4). There are several studies evaluating AP performance using different control strategies that have shown efficiency in clinical and free-living-condition trials (5–7). However, AP systems with dedicated safety schemes and/or adaptive laws are preferred over traditional control systems to reduce the risk of hypoglycemia in both fully closed-loop or hybrid developments (8–10).

Different adaptive control formulations based on MPC have been developed either to update the parameters of glycemia-insulin-carbohydrate models or to directly tune the control parameters. Strategies as those developed in (11, 12) and summarized in the review made by (13) consist of identifying the prediction model at every fixed period or when a condition is triggered like the detection of variations in glucose behavior. On the other hand, in (9, 14, 15), adaptive MPC strategies were formulated to change the penalization matrices of the MPC according to the current situation of the glycemia. Other works like (16–18) update other parameters of the strategy as the basal insulin, insulin-carbohydrate ratio, or the set-point, based on historical data and risk indexes to improve the controller performance.

In addition, the amount of active insulin that remains in the body, also known as insulin on board (IOB), has been used in several closed-loop safety schemes to prevent the insulin overdose, especially for the postprandial period (19–21). Different IOB models have been proposed to be part of open- and closed-loop controls with hypoglycemic prevention strategies (8, 22). In any case, an accurate value of the patient’s duration of insulin action (DIA) is required for a good performance.

Insulin pumps used in AP developments work with DIA values ranged from 2 to 8 hours in order to adjust their prevention strategies, as well as for the bolus calculators in hybrid systems (22, 23). A dynamic IOB constraint with estimated insulin action decay curves was first incorporated into the MPC problem by (20) and introduced in the safety module of the strategy proposed in (24). This idea was extended in (25, 26) where adaptive IOB rules were defined as a function of BG levels and the insulin delivery history was also considered to set the IOB boundary. Also (27), proposed an IOB decaying curve to set the input reference in the MPC cost function. The direct estimation of the IOB from a minimal physiological model was used in (28) where a control law was derived to avoid hypoglycemia. Additionally, a sliding-mode safety layer that runs outside the main control strategy was developed in (8) where the control action is modulated according to a surface computed as the difference between the estimation of the IOB and an IOB boundary. This strategy has been included in different works. In (29), the safety layer was coupled with a PID. In (30), the layer was coupled with a Linear Quadratic Gaussian (LQG) controller, and the IOB boundary was set as a piece-wise constant limit that changes according to a meal classification. In (6), the authors proposed a dynamic IOB boundary that depended on a factor of the open-loop IOB profile. In (10), it was shown that a single IOB boundary does not fit every situation of BG fluctuations as physiological variations in the subject change insulin requirements. Thus, a preliminary safety layer with adaptive IOB boundary was introduced.

Indeed, it has been shown that the range of insulin decay curves encompasses a variety of sources of uncertainty affecting the patient’s DIA, modifying the pharmacokinetics of the rapid-acting insulin analogues used in AP systems, whose value is usually between 3 and 5 hours (31). Common sources of uncertainty affecting how long it takes to absorb insulin, and thus an accurate estimate of IOB, are the patient’s insulin sensitivity, exercise, or heat.

Studies have revealed that insulin absorption can vary 10-30% in an individual and 20-50% between individuals (22). This can lead to overdoses of insulin when a shorter DIA value than actual is used in safety schemes based on the IOB, triggering episodes of hypoglycemia since the algorithm assumes there is less IOB than there actually is. In addition, it has been proved that selecting an inappropriate DIA setting is a common deficiency among clinicians and insulin pump users due to widespread misunderstanding of its concept (23).

Uncertainty associated to DIA setting, and therefore to the IOB estimation, has been addressed previously for AP designs through a method known as modal interval analysis (MIA) (32), where it is possible to consider a parametric uncertainty present in a dynamic model to be mathematically rewritten as an interval model (33). MIA allows to obtain a feasible simulation of the interval model providing an envelope of all the possible responses according to the percentage of parametric uncertainty established, avoiding under- and overestimation.

So far, only one AP design proposal implementing a safety scheme with MIA has been made (32) which was coupled with a PID-type controller. Although an IOB interval estimate was used, it was not implemented to generate the IOB boundaries required in the safety layer, and instead, the boundaries were chosen to be constant and to act like additional tuning parameters. Thus, there is no adaptation in this approach. In fact, when there are physiological variations towards hyperglycemia, the upper IOB limit (which is constant) could not be enough to allow the required insulin doses, and then prolonged hyperglycemia occurs.

This paper addresses the limitation found in previous works (referred to as the lack of consideration of the uncertainty in the IOB boundary selection) to obtain a control strategy that includes an adaptive safety module. To that end, the impulsive offset-free zone MPC (iZMPC) developed in (34) was used as the main controller, and concepts from the safety layer (6), and the uncertainty captured with the interval model in (32) were considered to build a new scheme. Now, the main features of the overall strategy are: (i) the safety layer uses the IOB interval model to generate upper and lower IOB boundaries that act as envelopes of the possible IOB trajectories; (ii) an auxiliary signal is generated to estimate the plant-model mismatch mainly given by physiological variations; (iii) with the information of the current glycemia, the estimated mismatch, and their corresponding trends and second derivatives, adaptive rules are developed to select the IOB boundary that is more appropriate for the current situation of the subject; (iv) with the selected IOB boundary at each time step, a gain is computed to modulate the control action of the MPC when required, and thus, to switch between aggressive and conservative control actions. To the best of the authors’ knowledge, this is the first time that the estimation of the plant-model mismatch, its derivatives, and the derivatives of the glycemia are used to adapt the safety layer.

A demanding assessment scenario with circadian variations in insulin sensitivity, the hepatic autoregulation, and endogenous glucose production, as well as comparisons with previous MPC designs, were considered in order to highlight the main improvements achieved. With the strategy in here formulated, a 92.7% of time in normoglycemia was obtained in comparison to the 83.8% resulting with the baseline iZMPC strategy and the 82.8% obtained with the preliminary strategy with safety layer developed in eViLeoGarRiv2021. In addition, the time in hypoglycemia was reduced to 0% in comparison to the 2.1% and 1.5% obtained with the previous strategies.

2 Methods

2.1 Glucose-Insulin-Carbohydrates System Dynamics

2.1.1 Impulsive System Dynamics for Glucose Control

The impulsive scheme for glucose control is selected to emulate the natural treatment of T1D. This because insulin doses are administered as small, spaced pulses rather than a continuous or a discrete input. Here, the impulsive discretization of a minimal physiological model based on five compartments to represent glucose dynamics, insulin absorption and action, and meal absorption dynamics is used. It is a linear model that satisfy global structural identifiability (35). The state-space representation of the model is given by

$\begin{array}{ll}\begin{array}{l}x(k+1)=Ax(k)+B_{u}u(k)+B_{r}r(k)+E,\hfill \\ y(k)=Cx(k),\hfill \end{array}& (1)\end{array}$

where matrices are related with their continuous counterpart considering a fixed sampling time T as $A=e^{A_{c}T}$, $B_u=e^{A_{c}T}{B}_{uc}$, $B_r={\displaystyle {\int }_0^T}{e}^{A_{c}s}ds{B}_{rc}$, $E=e^{A_{c}T}{E}_c$, with:

$\begin{array}{ll}{A}_c=\left[\begin{array}{ccccc}-p_0& -p_1& 0& p_2& 0\\ 0& -1/p_4& 1/p_4& 0& 0\\ 0& 0& -1/p_4& 0& 0\\ 0& 0& 0& -1/p_5& 1/p_5\\ 0& 0& 0& 0& -1/p_5\end{array}\right],B_{uc}=[\begin{array}{ccccc}0& 0& 1& 0& 0\end{array}]\text{'},B_{rc}={\left[\begin{array}{ccccc}0& 0& 0& 0& \frac{1}{p_5}\end{array}\right]}^{\prime },& (2)\end{array}$

E_c = [p₃ 0 0 0 0]', and C = [1 0 0 0 0]. All model variables and parameters have physiological interpretation as seen in Table 1. The model parameters have been identified for the 10-adult virtual subjects available in the commercial version of the UVA/Padova simulator (36). The identification process is detailed in (37), and the resulting parameters, personalized for each subject, can be seen in Table 2. In addition, an analysis of the properties of the impulse model can be found in the Appendix.

TABLE 1

Table 1 Description of variables and parameters of the model.

TABLE 2

Table 2 Model parameters identified from the 10-adult cohort of the UVA/Padova simulator.

2.1.2 Insulin-On-Board Model

An IOB model allows to estimate how much insulin is still to act in the body, which depends on the pharmacokinetics of the insulin analogue used and the patient’s glucose-insulin dynamics. To perform an IOB estimation, it is necessary to know the corresponding duration of insulin action (DIA) of the patient. In practice, there is no agreed-on standard for DIA. Many patients and clinicians enter inappropriately short DIA times into insulin pumps even though practical recommendations have been suggested for a proper choice (23).

In this work, the IOB model is obtained from subsystem (x_2,x₃) of model (1) using the IOB mathematical definition:

$\begin{array}{ll}IOB(k)=x_2(k)+x_3(k),& (3)\end{array}$

representing the amount of insulin in the subcutaneous and plasma compartments from previous boluses. This model is consistent with IOB models where the corresponding DIA time is implemented through K_DIA, a parameter that has been characterized to emulate typical DIA times (8). In model (3), K_DIA ≡1/p₄, and its value has been previously identified for each virtual subject considered here, see Table 2.

The Euler discrete time approximation of x₂ and x₃ of model (3), is given by

$\begin{array}{ll}\begin{array}{c}{x}_3(k+1)=Tu(k)-T{K}_{DIA}{x}_3(k)+x_3(k),\\ x_2(k+1)=T{K}_{DIA}(x_3(k)-x_2(k))+x_2(k),\end{array}& (4)\end{array}$

where u(k) is the insulin dose. This representation is used to obtain the IOB interval model considered in this work.

2.1.3 Insulin-On-Board Interval Model

One of the main challenges in AP control is the intra-patient variability that has to be taken into account. In different works as (38, 39), a variability of 30% has been observed and used to evaluate control performance. This uncertainty is represented here by an interval model in which the parameters, inputs, and/or initial states take interval values (33).

The simulation of a model with particular real values for the parameters, starting from any initial state, yields trajectories of the output variables across time. When the quantities involved in the simulation take values inside the intervals of variation, the set of trajectories determines a plane band bounded by an envelope, as depicted in Figure 1. At each time step of the simulation, the envelope, i.e., the possible maximum and minimum values of the variable, must be determined. This is a range computation problem. The function whose range must be determined is defined by the interval model of the system, and the parameter space is determined by the interval values of the parameters, the input, and the initial state. The simulation of an interval model provides intervals (ranges) that can be estimates of the envelopes. These envelopes can be obtained by numerical integration, qualitative reasoning, fuzzy logic, etc. A way to compute these estimates is by interval arithmetic. However, the exact range of a function is not always computable. Therefore, the results are often very over bounded, and if tighter results are needed, high computational efforts are required.

FIGURE 1

Figure 1 Envelopes obtained from the interval IOB model: Upper and lower bounds that define the possible system responses (shaded area).

To consider uncertainty in the IOB model (3), a modal interval analysis (MIA) is applied here to get an interval representation. MIA provides a strong theoretical background for dealing with problems involving uncertainty and logical quantifiers. In contrast to techniques such as Monte Carlo Simulation (MCS), MIA computational time is independent of the number of uncertain parameter. Thus, worst-case analysis can be performed efficiently, which is extremely important in the context of diabetes. A comparison between MCS and MIA was presented in (33). For further information on MIA, the reader can consult the appendix presented in (40).

MIA allows to obtain the whole range of possible responses of the model with uncertainty but avoiding, or minimizing under some conditions, the overestimation of interval computations. Overestimation is one of the main problems that arises from the existence of multiple instances of the same variable (multi-incident variable) with uncertainty. In model (4), uncertainty is considered through the patient’s K_DIA parameter, implementing it as $K_{DI{A}_i}=[0.7K_{DIA},1.3K_{DIA}]$.

To achieve an optimal computation, the IOB model (4) was initially analyzed to try to obtain its optimal form (the expression is rewritten in such a way that the exact range is obtained) to avoid multi-incident interval variables. This was possible for x₃ but not for x₂ of model (4), and therefore the theorem of coercion to optimality from MIA was also applied. To do this, the monotonic behaviour of x₂ was studied with respect to the uncertain variable K_DIA, leading to create a new interval state to obtain an optimal rational computation for the IOB model (3) as follows:

$\begin{array}{ll}\begin{array}{l}{x}_{3_i}(k+1)=x_{3_i}(k)(1-\Delta t{K}_{DI{A}_i})+\Delta tu(k),\hfill \\ x_{{23}_i}(k+1)=x_{{23}_i}(k)(1-\Delta t{K}_{DI{A}_i})+\Delta tu(k)+\Delta tDual(K_{DI{A}_i})x_{3_i}(k)\hfill \\ x_{2_i}(k+1)=x_{{23}_i}(k+1)-Dual(x_{3_i}(k+1)),\hfill \\ IO{B}_i(k)=x_{3_i}(k)+x_{2_i}(k),\hfill \end{array}& (5)\end{array}$

where the time step Δt, used for the interval model simulation, is defined so that constraint $\Delta t<1/x_{{23}_i}$ holds according to the theorem of coercion to optimality from MIA. $x_{3_i}$ and $x_{2_i}$ are interval representations of x₃ and x₂ in (1), respectively, and $x_{{23}_i}$ is the new interval state created to avoid under- and overestimation of the model. Note that, $x_{3_i}$,$x_{2_i}$,$x_{{23}_i}$, and IOBi are 2-dimension vectors containing lower and upper values. The dual operator is defined as Dual ([a₁, a₂]) = [a₂, a₁], according to MIA.

The output of this interval model corresponds to an upper (IOB_iu) and a lower (IOB_il) bound, that delimit an IOB envelope, see Figure 1. The envelope formed by these bounds is determined by the set of IOB trajectories involved in the simulation from all the possible responses taking values inside the interval of variation of the parameters with uncertainty (33). It is noteworthy that the output of the model is not simply the nominal response of the IOB in (3)by a factor of ±30%, in which case the envelope would be the same width all the time. On the contrary, as seen in Figure fig:IOBint, the width of the envelope, representing the range of possible instances, is also dynamic.

2.2 Control Strategy With Interval Safety Layer

2.2.1 Impulsive Model Predictive Control

In T1D treatment, the control objective is to maintain the system’s state within a safety zone X^Tar selected inside the desired range 70 mg/dl ≤ BG ≤ 140mg/dl. Bearing the above in mind, the MPC formulation to be used here is the impulsive offset-free zone MPC (denoted as iZMPC for simplicity). This formulation was developed in (34) and used in T1D context in (15, 41). The optimization problem solved by the iZMPC at every time step k is

$\begin{array}{ll}\begin{array}{l}\underset{u,x_a,u_a}{\min }{V}_N(x,X_s^{Tar},U_s^{Tar};u,x_a,u_a)\hfill \\ s.t.x(0)=\widehat{x}(k),d(0)=\widehat{d}(k),\hfill \\ x(j+1)=Ax(j)+B_{u}u(j)+B_{r}r(j)+B_{d}d(j)+E,\hfill \\ d(j+1)=d(j),\hfill \\ u(j)\in U,x(j)\in X,\hfill \\ x(H_p)=x_a,\hfill \\ y_a=C{x}_a+C_{d}d(j),\hfill \\ x_a=A{x}_a+B_u{u}_a+B_{d}d+E.\hfill \end{array}& (6)\end{array}$

with cost function stated as

$\begin{array}{ll}V(x,X_s^{Tar},U_s^{Tar};\mathit{u},u_a,x_a)={\displaystyle \sum _{j=0}^{H_p-1}\parallel x(j)-x_a{\parallel }_Q^2}+{\displaystyle \sum _{j=0}^{H_u-1}\parallel u(j)-u_a{\parallel }_R^2}+P\left(dis{t}_{X_s^{Tar}}(x_a)+dis{t}_{U_s^{Tar}}(u_a)\right).& (7)\end{array}$

From this iZMPC problem, it is to remark the following features:

1. The decision variables are u, which represents the control policy, and (x_a,u_a) ϵ (X_s, U_s), which are artificial equilibrium variables to be reached in the prediction horizon H_p, where the sets X_s, U_s denote the equilibrium sets for system (1). The control horizon is denoted by H_u.

2. The last term of the cost function ensures that the final point reached in the prediction remains within the target zone, where $dis{t}_{X_s^{Tar}}(x_a)$ denotes the distance of x_a to set $X_s^{Tar}$, and $X_s^{Tar}$, $U_s^{Tar}$are generalized equilibrium sets for the target (42).

3. The state is augmented with a disturbance d(k + 1) = d(k) to handle a constant plant-model mismatch and avoid offset error. Thus, the model is extended as x(k + 1)=Ax(k)+B_u u(k)+B_rr(k)+B_dd(k) + E, y(k) = Cx(k) + C_dd(k), with $B_d=e^{A_{c}T}{B}_{dc}$ and B_dc selected such that the augmented model be observable. In (34) it was shown that B_dc = [100 0 0 0 0]' satisfies this condition. Then, the augmented state is estimated with a state estimator to capture both the system state $\widehat{x}(k)$ and the mismatch $\widehat{d}(k)$ and it is used to initialize the optimization problem. In this work, the Kalman Filter is implemented.

4. In the MPC problem (6), the prediction model and equilibrium constraints are corrected to consider the disturbance effect, i.e., to consider d.

5. In this work, the control strategy has been designed as hybrid, in which meal announcement is provided to the controller by using the term B_rr(·) in the prediction model. In addition, sensing and actuation delays are not considered in the design of the control strategy, thus, the strategy should be able to compensate for these additional challenges despite the lack of information about it.

6. To solve the optimization problem at each time step, the quadprog solver of Matlab was used with the interior-point-convex algorithm and with tolerances TolPGG= 1×10^-5, Tolcon=1×10^-4, TolX= 1×10^-4, and Tolfun=1×10^-4 (see (43) for more details on the selection of the solver for quadratic programming in artificial pancreas context).

2.2.2 Interval Safety Layer

The main idea of the interval safety layer is to modulate the control action computed by the iZMPC to avoid hypoglycemia events. To this end, a gain γ ϵ [0 1] is calculated, such that the final command to the insulin pump u_f (k)=γ(k)u(k) satisfies an imposed constraint on the IOB. As shown in Figure 2, the layer consists of 6 blocks, from which the IOB is estimated and compared to a boundary selected to be safe in the current situation of the system.

FIGURE 2

Figure 2 Block diagram of the iZMPC strategy with the interval layer coupling.

First, the interval IOB is performed to obtain the upper (IOB_iu) and lower bound (IOB_il). In addition, a middle band (IOB_im) is obtained as the middle point between IOB_iu and IOB_il. These values are obtained with model (4) and the input u_OL corresponding to the open-loop insulin dose that should be administered to the subject, i.e., u_OL = u_basal for fasting periods, and u_OL = CHO/CR when a meal is announced. CHO corresponds to the carbohydrates ingested, and CR is the insulin-to-CHO ratio (35).

Next, from the three values calculated (IOB_iu, IOB_il and IOB_im), the IOB constraint $\overline{IOB}$ has to be selected. From previous works as (15), it was shown that the sign of the estimated plant-model mismatch provides information about the direction of the physiological variations occurring in the patient. For instance, a negative $\widehat{d}$ suggests variations inducing hypoglycemia, and a positive $\widehat{d}$ suggests variations towards hyperglycemia. To prevent measurement noise from affecting the information provided by d and to smooth the effect of food intake ton d, a moving horizon filter is used as $d^{*}(k)={\Sigma }_{k-N_p}^{k}d(k)/N_p$, where Np is the postprandial window usually between 1-2 hours (15). Here, the idea is to select the IOB constraint of the safety layer by using the information of the plant-model mismatch, but in addition, by including the rate of change of d* (or its first time derivative, μ_d), and its acceleration (or second time derivative, α_d), as well as information of the BG value, and its first μ and second time derivatives α. All these values are obtained from the signal generator block that receives information from the augmented state estimator and produces the derivatives by numerically computation using backward differentiation. The information obtained from these six signals helps to anticipate the variations and glucose behavior, and thus, to take action accordingly. The following rules have been posed to commute between aggressive and conservative control actions:

2.2.2.1 Case 1

The BG value is in the zone of imminent hyperglycemia, and it is detected that its value and rate of change are increasing (140 ≤ BG ≤ 180, μ ≥ 0, and α ≥ 0):

● If in addition the sign of the estimated mismatch has an increasing tendency (μ_d ≥ 0), then a risk of variations inducing hyperglycemia are considered too. Therefore, an aggressive response is established by allowing the delivery of the complete control action u computed by the iZMPC, i.e. the IOB constraint is released (u_f = u).

● If the variations are not increasing, then a less aggressive action is posed, but still allowing high doses of insulin. In this case, the IOB constraint is set as $\overline{IOB}=IO{B}_{iu}$.

2.2.2.2 Case 2

The BG value is not in the zone of imminent hyperglycemia, thus, the control action changes according to the information provided by the estimated mismatch:

● If either d ≤ –ϵ, or μ_d < 0 and d ≤ ϵ, then a risk of variations inducing hypoglycemia is detected, and thus, the lower band is established as the IOB constraint to obtain conservative control actions $\overline{(IOB}=IO{B}_{il})$.

● If d >–ϵ and μ_d ≥0 and α_d ≥0, then a risk of variations towards hyperglycemia is detected, thus, aggressive actions are allowed by changing the IOB constraint to the upper band $\overline{IOB}=IO{B}_{iu}$.

● In any other cases, in which it is detected that the variations are near to a change of direction, then softer control actions than rule 2.2 are allowed by using the middle band $\overline{IOB}=IO{B}_{im}$

A better visualization of the IOB constraint selection, using the cases described, can be seen in Figure 3. It is to clarify that the value of ϵ is set to consider the moments in which the physiological variations reach the higher andlower values. For the results here presented, a population value of ϵ = 0.005 was set from previous analysis on the behavior of d among the virtual subjects. Nevertheless, this parameter could be personalized by running the estimator in an open-loop treatment period to evaluate the magnitude of d. This implies that for the application of this control strategy, the collection of data prior to the activation of the control is a requirement. This practice is common in clinical studies, where data of each patient is collected prior to the clinical study to identify the prediction model, tune the estimator and controller matrices, and other parameters specific of the strategy (7).

FIGURE 3

Figure 3 Illustrative scenario for the IOB constraint selection according to signals of glycemia, the estimated mismatch, and their first and second time derivatives.

After selecting the IOB constraint, the IOB of the subject is estimated using model (3) and a shorter sampling time Δt than the one of the insulin pump T. At each time Δt, the IOB is compared to the selected constraint, constructing the function $\sigma (\Delta t)=\overline{IOB}-IOB(\Delta t)$. According to the sign of σ, the switching signal ω is determined:

$\begin{array}{ll}\omega (\Delta t)=\{\begin{array}{cc}1& if\sigma (\Delta t)\ge 0\\ 0& otherwise\end{array}.& (8)\end{array}$

Given the high frequency commutation between the values of ω, a chattering effect is obtained. To avoid it, the signal ω is smoothed by computed the gain γ(k) as the average of ω:

$$\gamma (k)=\frac{1}{n}{\displaystyle \sum _{t=1}^n\omega (t),}$$

With n = T/Δt the number of samples of w obtained during the pump updating period T. Finally, this gain is the one used to modulate the control action of the iZMPC.

3 Results

To evaluate the control strategy, 10 virtual adult subjects were simulated by identifying the model parameters from the commercially available version of the UVA-Padova T1DM Simulator (36). A 36-hour scenario was selected, the sampling time T was set as 5 minutes, and the daily pattern of carbohydrate intake was 7:00h (55g), 10:00h (20g), 13:00h (90g), and 19:00h (70g). Intra-day variability was introduced in the simulation scenario by modifying three of the model parameters. This was done by considering circadian variability with 30% amplitude of the form:

$$p_i(t)=p_i^{*}\left(1+0.3\sin \left(\frac{2\pi }{1440}\right)t+2\pi RND\right),$$

where $p_i^{*}$, with i = 0,1,3, is the nominal value of the parameters identified for each virtual subject and associated to the hepatic autoregulation, the insulin sensitivity, and the endogenous glucose production, and RND is a randomly uniformly generated number between 0 and 1 (38, 39). It should be noted that the parameter p₃ is 180° out of phase with respect to p₀ and p₁ such that the three parameters induce hypoglycemia or hyperglycemia at the same time and thus prevent the effect of the variations of each parameter from compensating each other.

The iZMPC controller with interval safety layer (iZMPC-ISL) is compared with the iZMPC without the layer, and the iZMPC with a safety layer (iZMPC-SL) that uses the constraint $\overline{IOB}=IO{B}_{im}$ all the time. This with the purpose of putting in evidence the benefits of adapting the IOB constraint to avoid hypoglycemia and react against hyperglycemia episodes.

In Figure 4, an illustrative example of virtual subject #10 is depicted. The Figure shows the glycemia, IOB, insulin, the estimated plant-model mismatch (d), and the gain γ that modulates the control action of the iZMPC. As shown, when using the iZMPC (blue lines) the system presents hypoglycemia episodes because of the delivery of insulin doses higher than the required for the situation of the subject (see that at 7:00h, the variations in the subject induce hypoglycemia). Next, the iZMPC-SL (red lines) is tested, but it can be seen the difficulties of setting one single constraint. If the constraint is too high to counteract variations inducing hypoglycemia, then the insulin doses are higher than the required, and if the constraint is lower than the required to compensate for variations inducing hyperglycemia (see the range between 15:00h and 25:00h) the controller is constrained and thus the insulin doses are not enough to low BG levels. To solve these issues, it can be seen that the iZMPC-ISL (black lines) uses the lower bound IOB_il as IOB constraint when detecting variations towards hypoglycemia, changes to the middle or upper band when detecting variations towards hyperglycemia and release the constraint when the BG levels and its tendency allow to do it (see that the IOB constraint is exceeded at 19:00h). It is also noteworthy that the safety layer seems to almost override the controller during the night, as the control action generated by the iZMPC is highly attenuated. This shows that the safety layer can be decisive at the expense of the iZMPC. However, this situation is desired in diabetes treatment since the control action is precisely canceled when a risk of hypoglycemia is predicted, and therefore, insulin administration must be reduced or suspended. In addition, by using the interval model, the IOB boundary is established at the limit that is actually needed to compensate for current patient BG levels, trend, and detected physiological variations.

FIGURE 4

Figure 4 Comparison of the system evolution with the iZMPC and iZMPC-ISL.

Next, by considering the same simulation scenario, the sensor noise was added to obtain the population outcomes. To that end, the CGM signal based on Dexcom G5 mobile devices was used. Its model and parameters can be seen in detail in (44). For each virtual subject, 10 simulations were performed, thus, 100 different cases were generated. The reported metrics of the population with each evaluated controller can be seen in Table 3. These consist of the mean BG (mg/dl), standard deviation (SD) of BG (mg/dl), coefficient of variation (CV) of BG (%), time percentage of BG in each range (%), and number of events in range. The outcome indexes are reported as mean ± SD for normally distributed data and as median (interquartile range) otherwise (45).

TABLE 3

Table 3 Performance comparison of the iZMPC, iZMPC-SL, iZMPC-OLSL, and iZMPC-ISL.

Figure 5 shows the comparison of the iZMPC and iZMPC-ISL. It is evident how the iZMPC-ISL significantly reduces hypoglycemia events without dangerously increasing BG levels. This is mainly achieved by reducing the insulin doses when required. The iZMPC-ISL vs. iZMPC improves the system performance in terms of time percentage in normoglycemia (92.7 ± 5.8 vs. 83.8 ± 6.7), reduces the time in hypoglycemia (0 ± 2.9 vs. 13.7 ± 6.5), the time in hyperglycemia increases (6.3 pm 5.5 vs. 2.4 ± 3.7) without obtaining events of severe hyperglycemia (BG > 250mg/dl). The complete outcomes can be seen in Table 3.

FIGURE 5

Figure 5 Comparison of the system evolution with the iZMPC and iZMPC-ISL.

In addition, Figure 6 shows the comparison of the iZMPC-ISL with two other strategies. First, the population results obtained with the iZMPC-SL are depicted. It can be seen how this strategy is harmful to compensate for hyperglycemia episodes. In fact, the time in hyperglycemia corresponds to 32.5% (see Table 3). As explained before, this behavior occurs because the IOB is constrained with the same boundary all the time, and thus, when the variations towards hyperglycemia occur(which have a time-random component), the control action is restricted, and the insulin doses delivered are not enough to counteract the variations.

FIGURE 6

Figure 6 Comparison of the system evolution under the iZMPC-SL, iZMPC-OLSL, and iZMPC-ISL.

Secondly, the strategy is compared with the one developed in (10), which uses an open-loop SL (here denoted as iZMPC-OLSL). This strategy was a first attempt of using information of the plant-model mismatch to change the IOB constraint. However, two limitations were detected: (i) the IOB constraint is computed as β times the IOB evolution considering the open-loop treatment, where the β factor is empirically established, and (ii) the rules for adapting the constraint dependedonly on the sign of d, and thus, there is not prediction regarding the changes in the variations (as in the iZMPC-ISL where the first and second time derivatives are used), leading to cases in which aggressive insulin doses are allowed precisely in the moments in which the variations start to decrease, which can cause hypoglycemia. As seen in Figure 6, the iZMPC-OLSL manages to avoid high BG levels and the hypoglycemia episodes in the first meals, nevertheless, for large meals late hypoglycemia occurs. Meanwhile, with the iZMPC-ISL, these events are avoided by using the first and second time derivatives of glycemia and the estimated mismatch to adapt the IOB constraint.

4 Conclusions

A control strategy for the artificial pancreas was introduced using an impulsive offset-free MPC coupled with an interval safety layer. The approach is intended to steer glycemia to the target zone while reducing the risk of hypoglycemia events. To that end, information of the parameter uncertainty is considered by means of an interval IOB model and the estimation of the plant-model mismatch with an augmented state estimator. The interval model establishes possible IOB boundaries to constraint the control action. While the current estimated mismatch, combined with the estimation of glycemia and their first and second time derivatives help to select the most suitable constraint for the current situation of the subject. In this regard, the controller changes between aggressive and conservative control actions to counteract hyperglycemia and hypoglycemia induced variations. This formulation achieved satisfactory results in an adult cohort and was compared with 3 different control approaches that vary the use of the safety layer. Future research includes the use of personalized information of the subject to establish some parameters, and/or the formulation of a function that explicitly relates the IOB constraint selection with the magnitude of the estimated variations.

Data Availability Statement

The original contributions presented in the study are included in the article/supplementary material. Further inquiries can be directed to the corresponding author.

Author Contributions

MV-T developed the code of the control strategy and performed the in-silico tests. MG-J contributed with the theoretical framework of the study about interval models. FL-V and PR conceived the study and analysis of the results, FL-V contributed with the theoretical framework of the safety layer, and PR with the framework of impulsive MPC. All authors contributed to draft the manuscript and approved the submitted version.

Funding

This work was supported by Minciencias (Colombia) with Grant 110180763081.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Haidar A. The Artificial Pancreas: How Closed-Loop Control is Revolutionizing Diabetes. IEEE Control Syst (2016) 36(5):28–47. doi: 10.1109/MCS.2016.2584318

CrossRef Full Text | Google Scholar

2. Hood T, Hovorka R. Coming of Age: The Artificial Pancreas for Type 1 Diabetes. Diabetologia (2016) 59(9):1795–805. doi: 10.1007/s00125-016-4022-4

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Steil GM. Algorithms for a Closed-Loop Artificial Pancreas: The Case for Proportional-Integral-Derivative Control. J Diabetes Sci Technol (2013), 1621–31. doi: 10.1177/193229681300700623

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Fathi AE, Smaoui MR, Gingras V, Boulet B, Haidar A. The Artificial Pancreas and Meal Control: An Overview of Postprandial Glucose Regulation in Type 1 Diabetes. IEEE Control Syst (2018) 7(6):67–85. doi: 10.1109/MCS.2017.2766323

CrossRef Full Text | Google Scholar

5. Dassau E, Pinsker J, Kudva Y, Brown S, Gondhalekar R, Dalla Man C, et al. 12-Week 24/7 Ambulatory Artificial Pancreas With Weekly Adaptation of Insulin Delivery Settings: Effect on Hemoglobin a 1c and Hypoglycemia. Diabetes Care (2017) 40(12):dc171188. doi: 10.2337/dc17-1188

CrossRef Full Text | Google Scholar

6. Fushimi E, Rosales N, De Battista H, Garelli F. Artificial Pancreas Clinical Trials: Moving Towards Closed-Loop Control Using Insulin-on-Board Constraints. Biomed Signal Process Control (2018) 45:1–9. doi: 10.1016/j.bspc.2018.05.009

CrossRef Full Text | Google Scholar

7. Garcia-Tirado J, Brown SA, Laichuthai N, Colmegna P, Koravi CL, Ozaslan B, et al. Anticipation of Historical Exercise Patterns by a Novel Artificial Pancreas System Reduces Hypoglycemia During and After Moderate-Intensity Physical Activity in People With Type 1 Diabetes. Diabetes Technol Ther (2021) 23(4):277–85. doi: 10.1089/dia.2020.0516

PubMed Abstract | CrossRef Full Text | Google Scholar

8. León-Vargas F, Garelli F, De Battista H, Vehí J. Postprandial Blood Glucose Control Using a Hybrid Adaptive PD Controller With Insulin-on-Board Limitation. Biomed Signal Process Control (2013) 8(6):724–32. doi: 10.1016/j.bspc.2013.06.008

CrossRef Full Text | Google Scholar

9. Shi D, Dassau E, Doyle FJ III. Adaptive Zone Model Predictive Control of Artificial Pancreas Based on Glucose- and Velocity-Dependent Control Penalties. IEEE Trans Biomed Eng (2019) 66(4):1045–54. doi: 10.1109/TBME.2018.2866392

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Villa-Tamayo MF, Leon-Vargas F, Garcia-Jaramillo M, Rivadeneira PS. Glycemic Control Strategy Based on an Impulsive MPC With Safety Layer Coupling for IOB Limitation. IEEE Control Syst Lett (2021) 5(5):1669–74. doi: 10.1109/LCSYS.2020.3042651

CrossRef Full Text | Google Scholar

11. Boiroux D, Bátora V, Hagdrup M, Wendt SL, Poulsen NK, Madsen H, et al. Adaptive Model Predictive Control for a Dual-Hormone Artificial Pancreas. J Process Control (2018) 68:105–17. doi: 10.1016/j.jprocont.2018.05.003

CrossRef Full Text | Google Scholar

12. Hajizadeh I, Rashid M, Cinar A. Plasma-Insulin-Cognizant Adaptive Model Predictive Control for Artificial Pancreas Systems. J Process Control (2019) 77:97–113. doi: 10.1016/j.jprocont.2019.03.009

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Turksoy K, Cinar A. Adaptive Control of Artificial Pancreas Systems - a Review. J Healthc Eng (2014) 5(1):1–22. doi: 10.1260/2040-2295.5.1.1

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Grosman B, Dassau E, Zisser H, Jovanovic L, Doyle FJ III. Multi-Zone-MPC: Clinical Inspired Control Algorithm for the Artificial Pancreas. IFAC Proc Vol (IFAC-PapersOnline) (2011) 44(1):7120–5. doi: 10.3182/20110828-6-IT-1002.02215

CrossRef Full Text | Google Scholar

15. Villa-Tamayo MF, Rivadeneira PS. Adaptive Impulsive Offset-Free MPC to Handle Parameter Variations for Type 1 Diabetes Treatment. Ind Eng Chem Res (2020) 59(13):5865–76. doi: 10.1021/acs.iecr.9b05979

CrossRef Full Text | Google Scholar

16. Toffanin C, Visentin R, Messori M, Palma FD, Magni L, Cobelli C. Toward a Run-to-Run Adaptive Artificial Pancreas: In Silico Results. IEEE Trans Biomed Eng (2018) 65(3):479–88. doi: 10.1109/TBME.2017.2652062

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Hajizadeh I, Samadi S, Sevil M, Rashid M, Cinar A. Performance Assessment and Modification of an Adaptive Model Predictive Control for Automated Insulin Delivery by a Multivariable Artificial Pancreas. Ind Eng Chem Res (2019) 58(26):11506–20. doi: 10.1021/acs.iecr.8b06202

CrossRef Full Text | Google Scholar

18. Liu W, Zhang G, Yu L, Xu B, Jin H. Improved Generalized Predictive Control Algorithm for Blood Glucose Control of Type 1 Diabetes. Artif Organs (2019) 43(4):386–98. doi: 10.1111/aor.13350

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Doyle IIIF, Huyett L, Lee J, Zisser H, Dassau E. Closed-Loop Artificial Pancreas Systems: Engineering the Algorithms. Diabetes Care (2014) 37(5):1191–7. doi: 10.2337/dc13-2108

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Ellingsen C, Dassau E, Zisser H, Grosman B, Percival M, Jovanovic L, et al. Safety Constraints in an Artificial Pancreatic βCell: An Implementation of Model Predictive Control With Insulin on Board. J Diabetes Sci Technol (2009) 3(3):536–44. doi: 10.1177/193229680900300319

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Forlenza GP. Relevance of Bolus Calculators in Current Hybrid Closed Loop Systems. Diabetes Technol Ther (2017) 19(7):400–1. doi: 10.1089/dia.2017.0216

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Zisser H, Robinson L, Bevier W, Dassau E, Ellingsen C, Doyle F, et al. Bolus Calculator: A Review of Four ßmart" Insulin Pumps. Diabetes Technol Ther (2008) 10(6):441–4. doi: 10.1089/dia.2007.0284

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Walsh J, Roberts R, Heinemann L. Confusion Regarding Duration of Insulin Action: A Potential Source for Major Insulin Dose Errors by Bolus Calculators. J Diabetes Sci Technol (2014) 8(1):170–8. doi: 10.1177/1932296813514319

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Breton M, Farret A, Bruttomesso D, Anderson S, Magni L, Patek S, et al. Fully Integrated Artificial Pancreas in Type 1 Diabetes: Modular Closed-Loop Glucose Control Maintains Near Normoglycemia. Diabetes (2012) 61(9):2230–7. doi: 10.2337/db11-1445

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Gondhalekar R, Dassau E, Doyle FJ III. Velocity-Weighting & Velocity-Penalty MPC of an Artificial Pancreas: Improved Safety & Performance. Automatica (2018) 91:105–17. doi: 10.1016/j.automatica.2018.01.025

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Lee JB, Dassau E, Gondhalekar R, Seborg DE, Pinsker JE, Doyle FJ. Enhanced Model Predictive Control (Empc) Strategy for Automated Glucose Control. Ind Eng Chem Res (2016) 55(46):11857–68. doi: 10.1021/acs.iecr.6b02718

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Messori M, Incremona GP, Cobelli C, Magni L. Individualized Model Predictive Control for the Artificial Pancreas. IEEE Control Syst Mag (2018) 38(1):86–104. doi: 10.1109/MCS.2017.2766314

CrossRef Full Text | Google Scholar

28. Magdelaine N, Rivadeneira PS, Moog C, Mohammadridha T, Chaillous L, Krempf M, et al. Hypoglycaemia-Free Artificial Pancreas Project. IET Syst Biol (2019) 14(1):16–23. doi: 10.1049/iet-syb.2018.5069

CrossRef Full Text | Google Scholar

29. Hu R, Li C. An Improved Pid Algorithm Based on Insulin-on-Board Estimate for Blood Glucose Control With Type 1 Diabetes. Comput Math Methods Med (2015) 2015:1–8. doi: 10.1155/2015/281589

CrossRef Full Text | Google Scholar

30. Colmegna P, Garelli F, De Battista H, Sánchez-Peña R. Automatic Regulatory Control in Type 1 Diabetes Without Carbohydrate Counting. Control Eng Pract (2018) 74:22–32. doi: 10.1016/j.conengprac.2018.02.003

CrossRef Full Text | Google Scholar

31. Donner T, Sarkar S, et al. Insulin – Pharmacology, Therapeutic Regimens, and Principles of Intensive Insulin Therapy. Feingold KR, Anawalt B, Boyce A, editors. South Dartmouth, MA: Endotext (2019).

Google Scholar

32. León-Vargas F, García-Jaramillo M, Molano A, De Battista H, Garelli F. Glucose Control for T1d Patients Based on Interval Models. Lect Notes Elec Eng (2021) 685:336–44. doi: 10.1007/978-3-030-53021-14

CrossRef Full Text | Google Scholar

33. Calm R, García-Jaramillo M, Bondia J, Sainz M, Vehí J. Comparison of Interval and Monte Carlo Simulation for the Prediction of Postprandial Glucose Under Uncertainty in Type 1 Diabetes Mellitus. Comput Methods Programs Biomed (2011) 104(3):325–32. doi: 10.1016/j.cmpb.2010.08.008

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Villa-Tamayo MF, Caicedo MA, Rivadeneira PS. Offset-Free MPC Strategy for Nonzero Regulation of Linear Impulsive Systems. ISA Trans (2020) 101:91–101. doi: 10.1016/j.isatra.2020.01.005

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Hoyos JD, Villa-Tamayo MF, Builes-Montano CE, Ramirez-Rincon A, Godoy JL, Garcia-Tirado J, et al. Identifiability of Control-Oriented Glucose-Insulin Linear Models: Review and Analysis. IEEE Access (2021) 9:69173–88. doi: 10.1109/ACCESS.2021.3076405

CrossRef Full Text | Google Scholar

36. Man CD, Micheletto F, Lv D, Breton M, B K, Cobelli C. The Uva/Padova Type 1 Diabetes Simulator: New Features. J Diabetes Sci Technol (2014) 8(1):26–34. doi: 10.1177/1932296813514502

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Abuin P, Rivadeneira PS, Ferramosca A, González A. Artificial Pancreas Under Stable Pulsatile MPC: Improving the Closed-Loop Performance. J Process Control (2020) 92:246–60. doi: 10.1016/j.jprocont.2020.06.009

CrossRef Full Text | Google Scholar

38. Herrero P, Bondia J, Adewuyi O, Pesl P, El-Sharkawy M, Reddy M, et al. Enhancing Automatic Closed-Loop Glucose Control in Type 1 Diabetes With an Adaptive Meal Bolus Calculator – in Silico Evaluation Under Intra-Day Variability. Comput Methods Progr Biomed (2017) 146:125–31. doi: 10.1016/j.cmpb.2017.05.010

CrossRef Full Text | Google Scholar

39. Resalat N, El Youssef J, Reddy R, Castle J, Jacobs PG. Adaptive Tuning of Basal and Bolus Insulin to Reduce Postprandial Hypoglycemia in a Hybrid Artificial Pancreas. J Process Control (2019) 80:247–54. doi: 10.1016/j.jprocont.2019.05.018

CrossRef Full Text | Google Scholar

40. García-Jaramillo M, Calm R, Bondia J, Vehí J. Prediction of Postprandial Blood Glucose Under Uncertainty and Intra-Patient Variability in Type 1 Diabetes: A Comparative Study of Three Interval Models. Comput Methods Programs Biomed (2012) 108(1):224–33. doi: 10.1016/j.cmpb.2012.04.003

PubMed Abstract | CrossRef Full Text | Google Scholar

41. González AH, Rivadeneira PS, Ferramosca A, Magdelaine N, Moog CH. Stable Impulsive Zone MPC for Type 1 Diabetic Patients Based on a Long-Term Model. Optim Control Appl Methods (2020) 41(6):2115–36. doi: 10.1002/oca.2647

CrossRef Full Text | Google Scholar

42. Rivadeneira PS, González AH, Ferramosca A. Control Strategies for Non-Zero Set-Point Regulation of Linear Impulsive Systems. IEEE Trans Automat Control (2018) 63(9):2994–3001. doi: 10.1109/TAC.2017.2776598

CrossRef Full Text | Google Scholar

43. Goez-Mora JE, Villa-Tamayo MF, Vallejo M, Rivadeneira PS. Performance Analysis of Different Embedded Systems and Open-Source Optimization Packages Towards an Impulsive Mpc Artifi Cial Pancreas. Front Endocrinol (2021) 12:662348. doi: 10.3389/fendo.2021.662348

CrossRef Full Text | Google Scholar

44. Vettoretti M, Facchinetti A, Sparacino G, Cobelli C. Type-1 Diabetes Patient Decision Simulator for in Silico Testing Safety and Effectiveness of Insulin Treatments. IEEE Trans Biomed Eng (2018) 65(6):1281–90. doi: 10.1109/TBME.2017.2746340

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Maahs D, Buckingham B, Castle J, Cinar A, Damiano E, Dassau E, et al. Outcome Measures for Artificial Pancreas Clinical Trials: A Consensus Report: Table 1. Diabetes Care (2016) 39(7):1175–9. doi: 10.2337/dc15-2716

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Rivadeneira PS, Moog CH. Observability Criteria for Impulsive Control Systems With Applications to Biomedical Engineering Processes. Automatica (2015) 44:125–31. doi: 10.1016/j.automatica.2015.02.042

CrossRef Full Text | Google Scholar

Appendix A1

In this Appendix, the stability, observability, and controllability of the discretized model (1) are analyzed. Furthermore, an analysis on the stability of the interval model is presented.

A1.1 Stability

The discretized system (1) maintains the stability of the original impulsive system (determined by the eigenvalues of A_c), as the discretization is a mapping with the exponential function $(A=e^{A_{c}T})$. Since the eigenvalues of A_c are real and negative, then the eigenvalues of the discrete system go to unit circle of the complex plane. Thus, maintaining the property of stability. It is noteworthy that this property does not depend on the value of T, as for all real and negative eigenvalues of A_c the exponential function maps them between 0 and 1.

A1.2 Observability

To implement the Kalman Filter, model (1) has to be observable. For impulsive systems, the observability criterion has been established in (46), where it is proved that a linear impulsive system is observable on some finite time interval if Rank[C^T A^TC^T …(A^n-1 )^T C^T] = n.

The symbolic observability matrix for model (2) is:

$$\begin{array}{l}Obs\\ =\left[\begin{array}{lllll}1\hfill & -p_o\hfill & p_0^2\hfill & -p_0^3\hfill & p_0^4\hfill \\ 0\hfill & -p_1\hfill & p_o{p}_1+\frac{p_1}{p_4}\hfill & -\frac{(p_0{p}_1+p_1/p_4)}{p_4}-p_0^2{p}_1\hfill & p_0^2(p_0{p}_1+\frac{p_1}{p_4})+\frac{(p_0{p}_1+p_1/p_4)}{p_4^2}\hfill \\ 0\hfill & 0\hfill & -\frac{p_1}{p_4}\hfill & \frac{(p_0{p}_1+p_1/p_4)}{p_4}+\frac{p_1}{p_4^2}\hfill & -2\frac{(p_0{p}_1+p_1/p_4)}{p_4^2}-\frac{p_1}{p_4^3}-\frac{(p_0^2{p}_1)}{p_4}\hfill \\ 0\hfill & p_2\hfill & -p_o{p}_2-\frac{p_2}{p_5}\hfill & \frac{(p_0{p}_2+p_2/p_5)}{p_5}+p_0^2{p}_2\hfill & -p_0^2(p_0{p}_2+\frac{p_2}{p_5})-\frac{(p_0{p}_2+p_2/p_5)}{p_5^2}\hfill \\ 0\hfill & 0\hfill & \frac{p_2}{p_5}\hfill & -\frac{(p_0{p}_2+p_2/p_5)}{p_5}-\frac{p_2}{p_5^2}\hfill & (2(p_0{p}_2+\frac{p_2}{p_5}))/p_5^2+\frac{p_2}{p_5^3}+\frac{(p_0^2(p_2)}{p_5}\hfill \end{array}\right],\end{array}$$

with determinant given by $det=p_1^2{p}_2^2\frac{p_4^4+p_5^4-4p_4{p}_5^3-4p_4^3{p}_5+6p_4^2{p}_5^2}{p_4^5{p}_5^5}$. The determinant satisfies det ≠ 0 except for p₁ = p₂ =0, which never occurs due to the physiological meaning of the parameters.

A1.3 Controllability

To implement the MPC, model (1) has to be controllable. The controllability criterion for impulsive control systems has been discussed in (1), where it is established that a system is impulsively controllable if and only if rank[B AB A²B …A^n-1B] = n. For the diabetes application, it is to clarify that, the state variables x₄. x₅ only depend on the meal disturbance r, and thus, they are not controllable (they are not affected by insulin). Therefore, the controllability of the model is evaluated without considering these states. The symbolic controllability matrix of the model is:

$$Co=\left[\begin{array}{ccc}0& 0& -p_2/p_4^2\\ 0& 1/p_4^2& -2/p_4^3\\ 1/p_4& -1/p_4^2& 1/p_4^3\end{array}\right],$$

which has full rank for all values of the parameters except for p₄ = 0 or p₂ = 0. However, as p₄ represents the time-to-maximum of effective insulin concentration, and p₂ is the carbohydrate bio-availability, their values are never zero.

A1.4 Stability of the Interval Model

The model in (5) can be rewritten in the form x (k + 1) = Ax(k) + Bu(k) when considering a system of 6 state variables ([x_3low x_3upper x_23low x_23upper x_2low x_2upper]). The resulting matrix A of the interval systemis:

$$A=\left[\begin{array}{cccccc}1-\frac{1}{p_1}& 0& 0& 0& 0& 0\\ \frac{1}{p_2}& 1-\frac{1}{p_1}& 0& 0& 0& 0\\ 0& 1& 0& -1& 0& 0\\ 0& 0& 0& 1-\frac{1}{p_2}& 0& 0\\ 0& 0& 0& \frac{1}{p_1}& 1-\frac{1}{p_2}& 0\\ -1& 0& 0& 0& 1& 0\end{array}\right]$$

with eigenvalues [0 0 (p1–1)/p1 (p1–1)/p1 (p2–1)/p2 (p2–1)/p2]'. These eigenvalues are always within the unit circle of the complex plane. Therefore, the interval model is asymptotically stable.