# Handling Parameter Variations during the Treatment of Type 1 Diabetes Mellitus: In Silico Results

Authors: María F. Villa Tamayo; Michelle A. Caicedo Alvarez; Pablo S. Rivadeneira
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

Type 1 diabetic patients need a strict treatment to regulate blood glucose concentration in a target range. Despite the development of different control strategies, the model parameter variations, given by physiological changes, can generate an inaccurate treatment and in consequence hyperglycemia and hypoglycemia episodes. Therefore, it is necessary to use control techniques that compensate such effects and maintain the control goals. Here, the effect of parametric variations is examined by the sensitivity analysis from which the most influential parameters in glycemia dynamics are detected. Based on that, an offset-free MPC strategy for impulsive systems is given for the first time in literature and simulated for type 1 diabetes treatment. This scheme along with the impulsive zone MPC with artificial variables reestablishes the normoglycemia behavior since the parameter variations are adequately rejected. However, only parametric variations up to 50% from their nominal values are well compensated, which suggests that more robust formulations are needed to ensure a greater rejection of physiological variations.


---

## Body

## 1. Introduction

Type 1 diabetes mellitus (T1DM) is an autoimmune disease in which the pancreaticβ-cells are destroyed, causing inability to secrete insulin and regulate blood glucose (BG) in the body. Additionally, it leads to secondary pathophysiological alterations in many systems and can be the cause of nephropathy, blindness, or even nontraumatic amputations of lower extremities [1]. To counteract this dysfunction, one of the most common treatments is functional insulin therapy. It consists of daily insulin injections according to glycemia measurements and carbohydrate intake with the objective of maintaining normoglycemia (70 mg/dl ≤ BG ≤ 180 mg/dl). Another treatment option is based on the development of the project called Artificial Pancreas (AP). AP attempts to emulate the natural behavior of the pancreas by the use of an insulin pump, continuous glucose monitoring, and closed-loop control strategies. This system aims to avoid hypoglycemia (BG < 70 mg/dl) and hyperglycemia (BG > 180 mg/dl) events, which result in complications in the patient with T1DM [2].AP has been extensively studied the last decade, and several control strategies have been proposed to close the loop, ranging from PID up to more complex strategies such as model predictive control. Recently, Medtronic has released to the market the first commercial device using a closed-loop strategy producing astonishing results. It consists of a proportional-integral-derivative with insulin feedback of a model-predicted insulin profile (PID-IFB) with safety constraints, to calculate the basal insulin delivery, and a manual option for bolus insulin [3]. This indicates that it is necessary to maintain the effort to get fully automated strategies. In this regard, one of the most studied control strategies is precisely model predictive control (MPC), which has shown adequate performances in clinical and in silico environments [4–6].The prediction capacity of the MPC strategy provides the selection of the optimal insulin administration that drives glycemia to the objective according to a cost function. Additionally, the zone MPC (ZMPC), which allows BG regulation in the target zone without requiring an extra correction to take it to a set-point, has shown a significant reduction of hypoglycemia and hyperglycemia events [7]. Recently, an incursion of impulsive control strategies applied to biomedical systems has been done. Due to the characteristics of T1DM treatment, the input can be approximated as an impulse when the insulin bolus is injected because of its short duration in relation to the sample time [8]. Many biomedical applications, such as HIV, malaria, and diabetes, are better modeled and controlled as impulsive systems [8]. In [9], an initial application of impulsive ZMPC (iZMPC) to T1DM patient dynamics is introduced. Furthermore, an iZMPC with artificial variables (iZMPC-AV) is developed in [10], where the use of these variables guarantees the convergence to an objective set and provides an enlarged domain of attraction.An essential component in MPC is the prediction model that describes the interaction between glycemia, insulin, and carbohydrate absorption. Several models have been developed: maximal models which seek a detailed representation of a diabetic patient and the control-oriented models. The most used model for control is the Bergman minimal model [11], which has been modified to include different features such as meal absorption [12]. Based on this results, advances as the improvement of the UVa/Padova simulator [13] and the implementation with control strategies as MPC [14] have been done. Recently, many control-oriented models have been introduced; in [15, 16], different long-term models are presented, with a validation of 2 and 12 days, respectively. These have realistic equilibriums and are simple due to their affine linear forms. Finally, in [17, 18], a modification of the last two models is developed, considering different time constants in the compartments of absorption of insulin and carbohydrates and the endogenous glucose. The model in [18] is used in this paper.So far, control techniques assume that the model adequately represents the process dynamics, which is not completely true due to the inherent mismatch implicated by the variability on human physiology caused by physical activity, hormonal changes, stress, the dawn phenomenon, and other factors. As a consequence of this assumption, when the resulting control is applied to the patient, hypoglycemia or hyperglycemia events may occur. With the purpose of understanding the plant-model mismatch, here, a sensitivity analysis is carried out as explained in [19] to determine the parameters that have more influence in glycemia. Subsequently, tests are conducted with different variation scenarios of the parameters.To overcome the inconvenience of the mismatch evidenced in the tests, the offset-free MPC strategy developed in [20, 21] is adopted. Here, an extension of this strategy is performed for impulsive systems for the first time in literature. To assess the performance of the impulsive MPC strategies, a large comparison using six different formulations is developed. In conclusion, the iZMPC-AV has the best performance. Based on that, this strategy is adapted to the offset-free scheme. The adaptation consists of an augmented system with a disturbance model and a state estimator (the Kalman estimator is chosen). The information provided by the estimator is used by the MPC control to achieve the target zone and, in clinical terms, to reduce the possibility of hypoglycemia and hyperglycemia episodes, even when there are parametric variations.The outline of the paper is as follows: In the Methodology, a description of the selected model for glucose-insulin-carbohydrate dynamics is presented. Then, a subsection is devoted to introduce the iZMPC-AV formulation. In addition, a sensitivity analysis of the closed-loop system is addressed to determine the parameters that have more influence in glycemia. The section concludes with the description of the offset-free control strategy. In the Results and Discussion, a performance comparison between the MPC strategies for discrete systems and the MPC strategies for impulsive systems, applied to T1DM treatment problems, is carried out. Afterwards, the issue shown by the mismatch is illustrated by the sensitivity analysis, and the results and performance evaluation of the offset-free control strategy are presented. Perspectives and conclusions are discussed at the end.


## 2. Methodology

The development of the offset-free MPC strategy applied to T1DM treatment consists of four fundamental subsections: the model to describe the interaction between glucose, insulin, and carbohydrates; the appropriate MPC formulation to maintain glycemia in the target zone; the sensitivity analysis of the influence of the parameters in the glycemia dynamics; and an extended MPC formulation considering the offset-free problem.


### 2.1. Glucose-Insulin-Carbohydrate Model

Different models have been developed to describe the glucose-insulin-carbohydrate interaction. Among them, the ones presented in [15, 16] are advantageous for having a long-term validity, and their equilibriums adequately describe the physiology of a patient with type 1 diabetes. Two recent models are those introduced in [17, 18]. Both models present a great similarity, differentiating themselves from the previously mentioned, mainly by including a term associated with glucose self-regulation mechanism, which allows the glucose level to be driven toward a basal state, stabilizing the system, and by using different time constants for each compartment of absorption of insulin. In this work, the selected model is the one in [18], since it also allows the use of different time constants for each compartment of absorption of carbohydrates. The model is formed by five state variables, and its state space representation is (1)x˙t=Axt+Buut+Brrt+E,A=-p0-p10p200-1p41p50000-1p600000-1p81p90000-1p11,Bu=001p700,Br=00001p10,E=p30000.The five state variables are as follows: x1, the glycemia-blood glucose concentration (mg/dl); x2 and x3, the delivery rates of insulin in the blood and interstitial space compartments, respectively (U/min); and x4 and x5, the delivery rates of carbohydrates in the stomach and duodenum compartments, respectively (g/min). The inputs are u, the exogenous insulin (U), and r, the intake of carbohydrates due to meals (g). The output of the system is the glycemia, i.e., y=10000x.The initial conditions and model parameters were identified for the 33 virtual patients from the UVa/Padova simulator [13] and also validated with clinical data of 42 real patients. The identification algorithm and details can be seen in [18]. One important feature is that all parameters have physiological meaning; their description and units are exhibited in Table 1.Table 1
Parameters description of the model.
Parameter Description Units p 0 Fractional rate of glucose self-regulation. 1/min p 1 Insulin action effectiveness. mg/U/min p 2 Carbohydrate bioavailability. - p 3 Net balance between endogenous glucose production and insulin independent glucose consumption. mg/dl/min p 4 ,  p 6 Diffusion time constants in the insulin compartments. min p 5 , p 7 Insulin diffusion effectiveness in each compartment. - p 8 ,  p 10 Diffusion time constants in the two digestion compartments. min p 9 , p 11 Meal absorption effectiveness in each compartment. -


### 2.2. Impulsive Zone Model Predictive Control with Artificial Variables Strategy

A recent approach for the control of the system is by considering the input of insulin as an impulse. Insulin doses are administered as small spaced pulses, rather than a continuous input or a discrete one. Then, the consideration of it as an impulsive system is appropriate to emulate the real T1DM treatment [9, 22]. Affine linear impulsive systems are described as follows:(2)x˙t=Axt+Brrt+E,x0=x0,t≠τk,xτk+=xτk+Buuτk,k∈N,where x∈X⊆Rn is the constrained state vector, u∈U⊆Rm denotes the constrained impulsive control input, and r⊆Rmr denotes the disturbances due to meals. Both X and U are compact sets and contain the origin. The impulsive-time instants are denoted by τk,k∈N, and τk+ refers to the time instant after τk. As explained in [23], from the relationship between the state after the free response xτk and the state after the jump xτk+ given by the application of the impulsive input, it is possible to expand the characterization of (2) by means of two underlying discrete-time subsystems at time instants τk and τk+, respectively: (3)x•k+1=A•x•k+Bu•u•k+Br•r•k+E•,x•0=xτ0,(4)x∘k+1=A∘x∘k+Bu∘u∘k+Br∘r∘k+E∘,x∘0=xτ0+,where A•=eAT,A∘=eAT,Bu•=eATBu,Bu∘=Bu, Br•=Br∘=∫0TeAsdsBr, E•=eATE, E∘=E,u∘k=u•k, and T=τk+1-τk is a constant period of impulses. The state xτk is denoted as x•(k), and xτk+ is x∘(k). The purpose is to control the whole system by the manipulation of (3).Afterward, the control objective is to drive the linear impulsive system to a nonempty target setXs•Tar⊆Xs⊆X, with its associated input target set UsTar⊆Us⊆U, where Xs and Us are the equilibrium sets of the system and maintain the system inside Xs•Tar indefinitely. Note that the equilibrium control is not applied continuously as it happens in discrete-time systems since now the control system is impulsive, and it cannot be applied in this form. The equilibrium control is also sent impulsively at each time τk, which leads to a generalized equilibrium. This concept was studied at [10] including its dynamical characterization.All MPC algorithms have three essential elements for their formulation: the prediction model needed to compute the predicted output, a cost function in a prediction horizonN, and the constraints associated with the nature of the system [24–26]. The iZMPC-AV formulation was developed and substantiated in [10]. Here, that formulation is applied to the T1DM treatment problem since there is a required zone to maintain glycemia, and the iZMPC-AV is a stable formulation. This strategy is chosen due to its performance in comparison to other MPC formulations for discrete and impulsive systems, as it will be shown later. The performance comparison of these formulations is established in the Results and Discussion.DenoteXs•⊆Xs as the equilibrium set to the subsystem (3) and its associated set of equilibrium inputs Us. The cost function of the iZMPC-AV is then(5)VNx•,r•,Xs•Tar,UsTar;u,ua,xa=Vdynx•,r•;u,xa,ua+VfXs•Tar,UsTar;xa,uawhere Vdyn and Vf are(6)Vdynx•,r•;u,xa,ua=∑i=0N-1x•k+i∣k-xaQ2+∑i=0N-1u•k+i∣k-uaR2(7)VfXs•Tar,UsTar;xa,ua=PdistXs•Tarxa+distUsTarua,with Q, R, P being semidefinite positive matrices. The term Vdyn is meant to steer the system to a certain artificial equilibrium given by ua,xa∈UsxXs•. The term Vf is a terminal cost function; it represents the deviation between the artificial equilibrium point and any other point (ut,xt) in the objective equilibrium set (UsTar,Xs•Tar).The optimization problem to be solved by the MPC at timek is (8)minu,ua,xaVNx•,r•,Xs•Tar,UsTar;u,ua,xas.t.x•0=x•k∣k,x•k+i∣k=A•x•k+i-1∣k+Bu•u•k+i-1∣k+Br•r•k+i-1∣k+E•,yk+i∣k=Cx•k+i∣k,u•k+i∣k∈U,x•k+i∣k∈X,x•k+N∣k=xa,xa=A•xa+Bu•ua+E•,xa∈Xs•,ua∈Uswhere the constraint x•(k+N∣k)=xa drives the terminal state at the end of the horizon to reach the artificial equilibrium xa, and the constraint xa=A•xa+Bu•ua+E• forces the pair (ua,xa) to fulfill the equilibrium condition. Notice that the information about the disturbance due to the meals, Br•r•∙, is entered into the formulation, but the control strategy will assume that the intake of meals is unannounced.


### 2.3. Sensitivity Analysis

To expose the dependence ofx to the parameters and their variations, a sensitivity analysis is performed. This denotes which of the parameters are the most influential in the change of blood glucose, for example. Assuming the identified parameters as the nominal value pn of the parameters p, the sensitivity function is defined [19]:(9)Sxi,pj=∂xi∂pj,(10)S˙xi,pjt=ASxi,pjt+∂fx,u,pn∂p,Sxi,pj0=0,where Sx,p represents the sensitivity of the state x to the parameters p=p0p2…p11 of model (1). Then, the analysis of variations of the most influential parameters is done to retrieve information about the effect on the control performance when there is a plant-model mismatch.


### 2.4. The Offset-Free iZMPC-AV Strategy

To compensate the plant-model mismatch given by variation of parameters, the offset-free MPC strategy introduced in [20, 21] is adapted here for the first time to the impulsive formulation stated in Section 2.2. The idea of the strategy is to have some knowledge about the mismatch and provide this information to the iZMPC-AV. For that purpose, the affine impulsive system (2) is augmented with an additional model of the disturbance to capture the plant-model mismatch in the steady state. It is assumed that the perturbation is constant, i.e., d˙t=0, and it affects the state in the impulse instants. Then, the augmented system has the form(11)x˙td˙t=A000xtdt+Br0rt+E0,x0d0=x00,t≠τk,xτk+=xτk+Buuτk+Bddτk,k∈N,yτk=Cxτk+Cddτk,in which d∈Rnd is the constant disturbance, and its associated matrices are Bd∈Rnxxnd and Cd∈Rnyxnd. Therefore, the underlying discrete system described in (3) and the output equation are rewritten as(12)x~•k+1=A~•x~•k+B~u•u•k+B~r•r•k+E~•+wk,x~•0=x~τ0,y~k=C~x~•k+vk,where x~•(k)=x•kd•k′ denotes the augmented state, y~ is the augmented output, and the augmented matrices of the new system (12) are given by A~•=A•Bd•0I, B~u•=B~u•0, B~r•=B~r•0, E~•=E•0, C~=CCd, and Bd•=eATBd. Additionally, the term w(k) denotes a zero-mean and Qe-covariance Gaussian process noise entering into the states, and v(k) is a zero-mean and Re-covariance Gaussian process noise affecting the measure of the system. Afterward, the condition in which the state and the disturbance can be estimated is that the pair C~,A~• is observable. It is satisfied if and only if(13)rankC~C~A~•⋮C~A~•nx+nd-1=nx+nd,and, then, the state and disturbance estimator is given by(14)x~^•k+1=A~•x~^•k+B~u•u•k+B~r•r•k+E~•+Key~k-C~x~^•k,in which the gain matrix Ke=K1K2′ is chosen so that the estimator is stable. The selected method to estimate the augmented state is the Kalman estimator [27]. This estimator consists of a prediction or a priori stage (15) denoted by x~^•-, which is an estimate of x~• before the measurement y~ is taken into account, and a correction stage or a posteriori stage (16) denoted by x~^•+, which is the estimate of x~• after the measurement is processed. Pe denotes the covariance of the estimation error, Ke is the Kalman filter gain, and Re and Qe are the covariance matrices associated with measurement noise and model noise, respectively.(15)x~^•-k+1=A~•x~^•+k+B~u•uk+B~r•r•k+E~•,Pe-k+1=A~•Pe+kA~•T+Qe,(16)Kek+1=Pe-k+1C~′C~Pe-k+1C~′+Re-1,x~^•+k+1=x~^•-k+1+Kek+1y~k+1-C~x~^•-k+1,Pe+k+1=I-Kek+1C~Pe-k+1I-Kek+1C~′+Kek+1ReKe′k+1Assuming that the estimator is stable and the number of disturbances is equal to the number of outputs, in steady state, it satisfies(17)A•-IBu•C0x•^∞u•∞=-E•-Bd•d^∞y∞-Cdd^∞,where the subscript ∞ denotes the steady state values. Given the current augmented state x~^•(k), the optimization problem coupled with the observation problem to solve in each time step k is (18)minu,ua,xaVNx~•,r•,Xs•Tar,UsTar;u,ua,xas.t.x~•0=x~^•k∣k,x~•k+i∣k=A~•x~•k+i-1∣k+B~u•u•k+i-1∣k+B~r•r•k+i-1∣k+E~•,y~k+i∣k=C~x~•k+i∣k,u•k+i∣k∈U,x•k+i∣k∈X,x•k+N∣k=xa,xa=A•xa+Bu•ua+Bd•d+E•,xa∈Xs•,ua∈Us,ya=Cxa+Cddk+i∣k,where VNx~•,r•,Xs•Tar,UsTar;u,ua,xa is the augmented version of the cost function in (5)-(7) and the difference between the plant and the model is included in the calculation of the artificial equilibrium xa. Let ya,∞ be the artificial reference in steady state. To assure that the variables will have zero offset, the steady state values must satisfy the observer condition (17) and the artificial values must satisfy the equilibrium constraint of the controller. (19)A•-IBu•C0xa,∞ua,∞=-E•-Bd•d^∞ya,∞-Cdd^∞In (18) the terms xa=A•xa+Bu•ua+Bd•d+E• and ya=Cxa+Cddk+i∣k refer to this equilibrium constraint of the controller. Note that, as in the nominal MPC formulation (see Section 2.2), the intake of meals is available information for the control strategy, but it will not be used since a secondary objective of the control is to be able to counteract the unknown disturbances, and then the meals will be considered as unannounced.


## 2.1. Glucose-Insulin-Carbohydrate Model

Different models have been developed to describe the glucose-insulin-carbohydrate interaction. Among them, the ones presented in [15, 16] are advantageous for having a long-term validity, and their equilibriums adequately describe the physiology of a patient with type 1 diabetes. Two recent models are those introduced in [17, 18]. Both models present a great similarity, differentiating themselves from the previously mentioned, mainly by including a term associated with glucose self-regulation mechanism, which allows the glucose level to be driven toward a basal state, stabilizing the system, and by using different time constants for each compartment of absorption of insulin. In this work, the selected model is the one in [18], since it also allows the use of different time constants for each compartment of absorption of carbohydrates. The model is formed by five state variables, and its state space representation is (1)x˙t=Axt+Buut+Brrt+E,A=-p0-p10p200-1p41p50000-1p600000-1p81p90000-1p11,Bu=001p700,Br=00001p10,E=p30000.The five state variables are as follows: x1, the glycemia-blood glucose concentration (mg/dl); x2 and x3, the delivery rates of insulin in the blood and interstitial space compartments, respectively (U/min); and x4 and x5, the delivery rates of carbohydrates in the stomach and duodenum compartments, respectively (g/min). The inputs are u, the exogenous insulin (U), and r, the intake of carbohydrates due to meals (g). The output of the system is the glycemia, i.e., y=10000x.The initial conditions and model parameters were identified for the 33 virtual patients from the UVa/Padova simulator [13] and also validated with clinical data of 42 real patients. The identification algorithm and details can be seen in [18]. One important feature is that all parameters have physiological meaning; their description and units are exhibited in Table 1.Table 1
Parameters description of the model.
Parameter Description Units p 0 Fractional rate of glucose self-regulation. 1/min p 1 Insulin action effectiveness. mg/U/min p 2 Carbohydrate bioavailability. - p 3 Net balance between endogenous glucose production and insulin independent glucose consumption. mg/dl/min p 4 ,  p 6 Diffusion time constants in the insulin compartments. min p 5 , p 7 Insulin diffusion effectiveness in each compartment. - p 8 ,  p 10 Diffusion time constants in the two digestion compartments. min p 9 , p 11 Meal absorption effectiveness in each compartment. -


## 2.2. Impulsive Zone Model Predictive Control with Artificial Variables Strategy

A recent approach for the control of the system is by considering the input of insulin as an impulse. Insulin doses are administered as small spaced pulses, rather than a continuous input or a discrete one. Then, the consideration of it as an impulsive system is appropriate to emulate the real T1DM treatment [9, 22]. Affine linear impulsive systems are described as follows:(2)x˙t=Axt+Brrt+E,x0=x0,t≠τk,xτk+=xτk+Buuτk,k∈N,where x∈X⊆Rn is the constrained state vector, u∈U⊆Rm denotes the constrained impulsive control input, and r⊆Rmr denotes the disturbances due to meals. Both X and U are compact sets and contain the origin. The impulsive-time instants are denoted by τk,k∈N, and τk+ refers to the time instant after τk. As explained in [23], from the relationship between the state after the free response xτk and the state after the jump xτk+ given by the application of the impulsive input, it is possible to expand the characterization of (2) by means of two underlying discrete-time subsystems at time instants τk and τk+, respectively: (3)x•k+1=A•x•k+Bu•u•k+Br•r•k+E•,x•0=xτ0,(4)x∘k+1=A∘x∘k+Bu∘u∘k+Br∘r∘k+E∘,x∘0=xτ0+,where A•=eAT,A∘=eAT,Bu•=eATBu,Bu∘=Bu, Br•=Br∘=∫0TeAsdsBr, E•=eATE, E∘=E,u∘k=u•k, and T=τk+1-τk is a constant period of impulses. The state xτk is denoted as x•(k), and xτk+ is x∘(k). The purpose is to control the whole system by the manipulation of (3).Afterward, the control objective is to drive the linear impulsive system to a nonempty target setXs•Tar⊆Xs⊆X, with its associated input target set UsTar⊆Us⊆U, where Xs and Us are the equilibrium sets of the system and maintain the system inside Xs•Tar indefinitely. Note that the equilibrium control is not applied continuously as it happens in discrete-time systems since now the control system is impulsive, and it cannot be applied in this form. The equilibrium control is also sent impulsively at each time τk, which leads to a generalized equilibrium. This concept was studied at [10] including its dynamical characterization.All MPC algorithms have three essential elements for their formulation: the prediction model needed to compute the predicted output, a cost function in a prediction horizonN, and the constraints associated with the nature of the system [24–26]. The iZMPC-AV formulation was developed and substantiated in [10]. Here, that formulation is applied to the T1DM treatment problem since there is a required zone to maintain glycemia, and the iZMPC-AV is a stable formulation. This strategy is chosen due to its performance in comparison to other MPC formulations for discrete and impulsive systems, as it will be shown later. The performance comparison of these formulations is established in the Results and Discussion.DenoteXs•⊆Xs as the equilibrium set to the subsystem (3) and its associated set of equilibrium inputs Us. The cost function of the iZMPC-AV is then(5)VNx•,r•,Xs•Tar,UsTar;u,ua,xa=Vdynx•,r•;u,xa,ua+VfXs•Tar,UsTar;xa,uawhere Vdyn and Vf are(6)Vdynx•,r•;u,xa,ua=∑i=0N-1x•k+i∣k-xaQ2+∑i=0N-1u•k+i∣k-uaR2(7)VfXs•Tar,UsTar;xa,ua=PdistXs•Tarxa+distUsTarua,with Q, R, P being semidefinite positive matrices. The term Vdyn is meant to steer the system to a certain artificial equilibrium given by ua,xa∈UsxXs•. The term Vf is a terminal cost function; it represents the deviation between the artificial equilibrium point and any other point (ut,xt) in the objective equilibrium set (UsTar,Xs•Tar).The optimization problem to be solved by the MPC at timek is (8)minu,ua,xaVNx•,r•,Xs•Tar,UsTar;u,ua,xas.t.x•0=x•k∣k,x•k+i∣k=A•x•k+i-1∣k+Bu•u•k+i-1∣k+Br•r•k+i-1∣k+E•,yk+i∣k=Cx•k+i∣k,u•k+i∣k∈U,x•k+i∣k∈X,x•k+N∣k=xa,xa=A•xa+Bu•ua+E•,xa∈Xs•,ua∈Uswhere the constraint x•(k+N∣k)=xa drives the terminal state at the end of the horizon to reach the artificial equilibrium xa, and the constraint xa=A•xa+Bu•ua+E• forces the pair (ua,xa) to fulfill the equilibrium condition. Notice that the information about the disturbance due to the meals, Br•r•∙, is entered into the formulation, but the control strategy will assume that the intake of meals is unannounced.


## 2.3. Sensitivity Analysis

To expose the dependence ofx to the parameters and their variations, a sensitivity analysis is performed. This denotes which of the parameters are the most influential in the change of blood glucose, for example. Assuming the identified parameters as the nominal value pn of the parameters p, the sensitivity function is defined [19]:(9)Sxi,pj=∂xi∂pj,(10)S˙xi,pjt=ASxi,pjt+∂fx,u,pn∂p,Sxi,pj0=0,where Sx,p represents the sensitivity of the state x to the parameters p=p0p2…p11 of model (1). Then, the analysis of variations of the most influential parameters is done to retrieve information about the effect on the control performance when there is a plant-model mismatch.


## 2.4. The Offset-Free iZMPC-AV Strategy

To compensate the plant-model mismatch given by variation of parameters, the offset-free MPC strategy introduced in [20, 21] is adapted here for the first time to the impulsive formulation stated in Section 2.2. The idea of the strategy is to have some knowledge about the mismatch and provide this information to the iZMPC-AV. For that purpose, the affine impulsive system (2) is augmented with an additional model of the disturbance to capture the plant-model mismatch in the steady state. It is assumed that the perturbation is constant, i.e., d˙t=0, and it affects the state in the impulse instants. Then, the augmented system has the form(11)x˙td˙t=A000xtdt+Br0rt+E0,x0d0=x00,t≠τk,xτk+=xτk+Buuτk+Bddτk,k∈N,yτk=Cxτk+Cddτk,in which d∈Rnd is the constant disturbance, and its associated matrices are Bd∈Rnxxnd and Cd∈Rnyxnd. Therefore, the underlying discrete system described in (3) and the output equation are rewritten as(12)x~•k+1=A~•x~•k+B~u•u•k+B~r•r•k+E~•+wk,x~•0=x~τ0,y~k=C~x~•k+vk,where x~•(k)=x•kd•k′ denotes the augmented state, y~ is the augmented output, and the augmented matrices of the new system (12) are given by A~•=A•Bd•0I, B~u•=B~u•0, B~r•=B~r•0, E~•=E•0, C~=CCd, and Bd•=eATBd. Additionally, the term w(k) denotes a zero-mean and Qe-covariance Gaussian process noise entering into the states, and v(k) is a zero-mean and Re-covariance Gaussian process noise affecting the measure of the system. Afterward, the condition in which the state and the disturbance can be estimated is that the pair C~,A~• is observable. It is satisfied if and only if(13)rankC~C~A~•⋮C~A~•nx+nd-1=nx+nd,and, then, the state and disturbance estimator is given by(14)x~^•k+1=A~•x~^•k+B~u•u•k+B~r•r•k+E~•+Key~k-C~x~^•k,in which the gain matrix Ke=K1K2′ is chosen so that the estimator is stable. The selected method to estimate the augmented state is the Kalman estimator [27]. This estimator consists of a prediction or a priori stage (15) denoted by x~^•-, which is an estimate of x~• before the measurement y~ is taken into account, and a correction stage or a posteriori stage (16) denoted by x~^•+, which is the estimate of x~• after the measurement is processed. Pe denotes the covariance of the estimation error, Ke is the Kalman filter gain, and Re and Qe are the covariance matrices associated with measurement noise and model noise, respectively.(15)x~^•-k+1=A~•x~^•+k+B~u•uk+B~r•r•k+E~•,Pe-k+1=A~•Pe+kA~•T+Qe,(16)Kek+1=Pe-k+1C~′C~Pe-k+1C~′+Re-1,x~^•+k+1=x~^•-k+1+Kek+1y~k+1-C~x~^•-k+1,Pe+k+1=I-Kek+1C~Pe-k+1I-Kek+1C~′+Kek+1ReKe′k+1Assuming that the estimator is stable and the number of disturbances is equal to the number of outputs, in steady state, it satisfies(17)A•-IBu•C0x•^∞u•∞=-E•-Bd•d^∞y∞-Cdd^∞,where the subscript ∞ denotes the steady state values. Given the current augmented state x~^•(k), the optimization problem coupled with the observation problem to solve in each time step k is (18)minu,ua,xaVNx~•,r•,Xs•Tar,UsTar;u,ua,xas.t.x~•0=x~^•k∣k,x~•k+i∣k=A~•x~•k+i-1∣k+B~u•u•k+i-1∣k+B~r•r•k+i-1∣k+E~•,y~k+i∣k=C~x~•k+i∣k,u•k+i∣k∈U,x•k+i∣k∈X,x•k+N∣k=xa,xa=A•xa+Bu•ua+Bd•d+E•,xa∈Xs•,ua∈Us,ya=Cxa+Cddk+i∣k,where VNx~•,r•,Xs•Tar,UsTar;u,ua,xa is the augmented version of the cost function in (5)-(7) and the difference between the plant and the model is included in the calculation of the artificial equilibrium xa. Let ya,∞ be the artificial reference in steady state. To assure that the variables will have zero offset, the steady state values must satisfy the observer condition (17) and the artificial values must satisfy the equilibrium constraint of the controller. (19)A•-IBu•C0xa,∞ua,∞=-E•-Bd•d^∞ya,∞-Cdd^∞In (18) the terms xa=A•xa+Bu•ua+Bd•d+E• and ya=Cxa+Cddk+i∣k refer to this equilibrium constraint of the controller. Note that, as in the nominal MPC formulation (see Section 2.2), the intake of meals is available information for the control strategy, but it will not be used since a secondary objective of the control is to be able to counteract the unknown disturbances, and then the meals will be considered as unannounced.


## 3. Results and Discussion

### 3.1. Glycemia Control in T1DM Patients: Comparison of Six MPC Strategies

In this subsection, different MPC strategies for discrete and impulsive systems, applied to the T1DM treatment problem, are compared. These strategies are the standard discrete MPC (dMPC) [24], the discrete ZMPC (dZMPC) with slack variable [26], the discrete ZMPC with artificial variables (dZMPC-AV) [28, 29], and the adaptation of the same formulations for impulsive systems: iMPC [30], iZMPC [22], and iZMPC-AV [10]; the description of each formulation can be seen in the appendix.The three mentioned techniques for discrete and impulsive systems are applied using model (1) and assuming that the plant and the model are equal, i.e., the obtained results are from nominal conditions. As stated before, the formulation of each of the control strategies does not consider the announcement of meals; instead these are considered as unknown disturbances. Therefore, the term Brrt is not considered in the prediction model. The simulated scenario is based on 33 virtual patients from the UVa/Padova simulator, which has 10 children, 10 adolescents, and 10 adults with their average patients, respectively. The initial condition for x1 is set in 120 mg/dl. The objective is to maintain the system in the target zone, even when food intake is entered into the system. The time of the simulation is 24 hours, and meals considered are 50 g of breakfast at 7:00h, lunch of 80g at 12:00h, and dinner of 60g at 21:00h. The prediction horizon N and the matrices Q, R, P are adjusted for each patient using the Nelder-Mead method for autotuning to maximize the time inside the target zone [31]. This tuning procedure is appropriate for in silico environments but other strategies have to be considered for use in real-life conditions such as the one presented in [32]. To all patients, the constraint sets are 0,0,0≤X≤1000,10,10 and 0≤U≤30. In addition, the reference for all strategies is set in 100 mg/dl. The target zone of the glycemia is set as 80≤x1≤140 for all controllers, and the zone called normoglycemia is defined as 70≤x1≤180. As detailed in [33], a value above that range is considered as level 1 hyperglycemia and a value below it as level 1 hypoglycemia, and the critical events x1< 54 mg/dl and x1> 250 mg/dl are denominated as level 2 hypoglycemia and level 2 hyperglycemia, respectively.For all virtual patients evaluated, the performance statistics are computed for every strategy as recommended in [34], taking into account the fact that the results presented here are obtained from simulation and not from clinical trials. Table 2 includes the percentage of time in which glycemia is found both in and out of the normoglycemia zone and the amount of hypoglycemia and hyperglycemia events. Furthermore, the mean and standard deviation (SD) of BG with each controller are calculated. For each outcome the median and the interquartile range (IQR) of the population are shown, except for the mean and SD of BG concentration which are normally distributed so that the average and SD of the 11 patients are shown. The glycemia evolution for the 33 patients using the six different control strategies can be seen in the appendix, which shows that all controllers maintain the patients in the safe zone most of the time without risk of level 2 hypoglycemia. From the results, it is inferred that the iZMPC-AV presents a better performance than the other strategies. It has the highest percentage of precise control, achieving an 83.9% (IQR 11.8%) of time inside the target zone with a median of 1 (IQR 2) event that exceeds 180 mg/dl corresponding to 1.3% (7.9%) of the time. According to this, the iZMPC-AV is selected as the controller to be used in the offset-free control strategy, as was anticipated in early sections.Table 2
Performance comparison of the MPC strategies for discrete and impulsive systems. Values are medians (interquartile range) unless stated otherwise.
Strategy iMPC dMPC iZMPC dZMPC iZMPC-AV dZMPC-AV Mean (SD) BG (mg/dl) 127.8 (10.7) 121.6 (12.2) 129.2 (4.9) 125.5 (8.8) 112.5 (8.7) 112.3 (8.4)  SD (SD) (mg/dl) 36.2 (12.7) 30.9 (11.9) 27.9 (12.1) 31.7 (10.9) 27.5 (13.0) 29.0 (12.8)  Time percentage of BG (%)  < 7 0 mg/dl 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) ∈ 70,140 mg/dl 73.2 (11.4) 76.9 (11.9) 74.1 (7.8) 75.2 (11.3) 83.9 (11.8) 83.8 (12.8) ∈ 70,180 mg/dl 91.5 (14.6) 94.9 (14.8) 94.9 (9.5) 93.9 (11.4) 98.7 (7.9) 97.8 (9.3) > 180 mg/dl 8.5 (14.6) 5.1 (14.8) 5.0 (9.5) 6.1 (11.4) 1.3 (7.9) 2.2 (9.3) > 250 mg/dl 0 (2.0) 0 (1.1) 0 (0) 0 (1.6) 0 (0) 0 (0) > 300 mg/dl 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0)  Number of events BG (-)  < 70 mg/dl 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) > 180 mg/dl 2(2) 2(3) 1(2) 1(2) 1(2) 1(2) > 250 mg/dl 0(1) 0(1) 0 (0) 0(1) 0 (0) 0 (0) > 300 mg/dl 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0)


### 3.2. Sensitivity Analysis

The sensitivity function defined in (9) and given by the differential equation (10) is solved simultaneously with x to obtain its time evolution. The proposed simulation scenario consists of two days, only one of which includes the intake of carbohydrates, with the purpose of demonstrating the sensitivity to parameters during fasting and when there is meal intake. The considered meals are 50g at 7:00h, 80g at 12:00h, and 60g at 21:00h. In Figure 1, the magnitude of the sensitivity of glycemia to pj, denoted as Sx1,pj, is plotted for the average adult as illustration. The parameters are categorized according to the behavior of the sensitivity magnitude, noticing similar responses to p0,p1,p2, and p3 associated with glycemia; p4,p5,p6, and p7 with insulin; and p8,p9,p10, and p11 with intake of carbohydrates.Figure 1
Glycemia sensitivity to each parameter of the model. Top: parameters associated with glycemia. Lower left: parameters associated with insulin. Lower right: parameters associated with intake of carbohydrates.According toSx1,pj, the parameters p0,p1,p2,p3,p5, and p8 are selected as the most influential in blood glucose concentration. Therefore, the simulation of the plant-model mismatch is done when variations in these parameters are considered. The closed-loop system with the iZMPC-AV control strategy is simulated for variations of 10%, 30%, and 50% in the selected parameters and illustrated using the average patient. The behavior of glycemia in the second day confirms the offset provoked by these variations. In Figure 2, the effect of each of these is depicted and the following statements are established:(i)
By increasing the parameterp0 in the plant, less amount of insulin is required. As the prediction model has a lower value of p0, the controller overestimates the necessary dose of insulin and thus the glycemia decreases. If the value of p0 is reduced, the required insulin is higher, but the controller underestimates the amount of insulin, which causes hyperglycemia.(ii)
Variations above the nominal value ofp1 steer to lower levels of BG, obtaining hypoglycemia events. This is because the patient is more sensitive to insulin and therefore requires a smaller amount in each dose, whereas, as the prediction model has a lower value of the parameter, the calculated insulin dose is higher than the actually needed. When insulin action effectiveness parameter is lower than the nominal value, the glycemia increases and hyperglycemia events may occur.(iii)
The variations in parameterp2, corresponding to carbohydrate bioavailability, directly reduce or increase BG levels at the moment of the carbohydrate intake. When the real value of p2 increases, the patient gets more sensitive to meal disturbances and episodes of hyperglycemia in the postprandial period are of greater magnitude. When the real value of p2 is lower than the identified parameter, the effect of meals in the BG levels is reduced.(iv)
Variations inp3 affect glycemia at the time of meals and at steady state. An increase in this parameter leads to an increment of glycemia and it tends to hyperglycemia events. A decrease in the parameter results in a reduction of BG levels and then in hypoglycemia events.(v)
Changes in the parameterp5 are more critical at the time of meals and when its real value is less than the nominal value, leading to hypoglycemia episodes. When the effectiveness of insulin diffusion is greater, the glycemia rises.(vi)
Variations inp8 have an effect on glycemia only at the time of meals. When it is augmented, the glycemia is more sensitive to meal disturbances due to a longer time of diffusion of carbohydrates, leading the patient to have hyperglycemia episodes. By decreasing the value of the parameter, glycemia is reduced in the postprandial period. There is no offset presented because of its dependence on carbohydrate intake.Figure 2
Variation of parametersp0, p1, p2, p3, p5, and p8 in 10%, 30%, and 50% of their nominal values.


### 3.3. Offset-Free Strategy for Glycemia Control with Plant-Model Mismatch

In this subsection, the iZMPC-AV strategy is affected by parameter variations. These variations normally appear in association with exercise, the dawn phenomenon, production of endogenous glucose, among others [35]. As was illustrated in the above subsection, depending on the magnitude of the parameter variation, the iZMPC-AV strategy cannot compensate the disturbance to maintain the adequate control. For that reason, the offset-free iZMPC-AV (iZMPC-AV-OF) strategy is designed. Its performance is evaluated and compared with the iZMPC-AV without the offset-free strategy (iZMPC-AV-NOF) under the next two different scenarios:Scenario 1.
Variation of the parameters to induce hyperglycemia: the parametersp2,p3,p5, and p8 are increased with respect to their nominal values and parameters p0 and p1 are decreased by the same percentages.Scenario 2.
Variation of the parameters to induce hypoglycemia: the parametersp0 and p1 are increased with respect to their nominal values and parameters p2,p3,p5, and p8 are decreased by the same percentages.Asx1 is the only measured state that corresponds to glycemia, the state variables x2, x3, x4, and x5 must be estimated, as well as the disturbance d to correct the steady state mismatch. The simulation consists of 36 hours, starting at 0:00h. The meals considered are a breakfast of 50g at 7:00h, a lunch of 80g at 12:00h, and a dinner of 60g at 21:00h. The control objective is to drive glycemia to the zone between 70 mg/dl and 180 mg/dl and remain there as long as possible. The control strategies do not consider meal announcements.To illustrate the resulting effect when the offset-free strategy is used, Figure3 depicts a comparison of the evolution of the glycemia and the insulin doses for the average adult patient under Scenario 2 considering variations of 0% (nominal situation), 30%, and 50%. The dashed lines are produced when the plant is controlled with the iZMPC-AV-NOF strategy, and the solid lines when the iZMPC-AV-OF strategy is applied. The maximum variation in parameters before the occurrence of hypoglycemia events (or prolonged hyperglycemia events) is established at 50% of the nominal values of the parameters. Meals are only given on the first day to expose the correction at steady state of the second day.Figure 3
Comparison of glycemia evolution and insulin doses for the average adult patient under Scenario2, considering variations of 0%, 30%, and 50%. Dashed lines (blue, purple, and green) represent iZMPC-AV-NOF while solid lines (yellow, red, and black) represent iZMPC-AV-OF.It is noticed that with the iZMPC-AV-NOF the glycemia levels are adequately regulated when the plant and the model are equal, but they decrease to critical BG levels when there are parameter variations of 30% and 50%, maintaining glycemia below 50 mg/dl at meal times. At steady state (from 30:00h to 36:00h), the offset is well appreciated since the glycemia is out of the target zone. This behavior shows the need to apply a strategy that compensates the plant-model mismatch. On the other hand, the iZMPC-AV-OF manages to maintain the glycemia in the desired zone despite the variations. For variations of 30%, the glycemia remains above 80 mg/dl, and for variations of 50% the glycemia remains between 70 and 180 mg/dl and the effect of meal disturbances is reduced. When the iZMPC-AV-OF is applied to the nominal model, the behavior is very similar to that obtained with the iZMPC-AV-NOF; there are negligible differences at meal times, but at steady state the behavior is the same.Regarding the control action for both strategies, when there is a plant-model mismatch and the iZMPC-AV-NOF is applied, it is observed that the offset at steady state is given by a miscalculation of the doses of insulin that are administered. The basal insulin is greater than that needed for the plant with the mismatch. In contrast, when applying the iZMPC-AV-OF strategy, the amount of basal insulin administered is corrected, and the dose is reduced. Finally, for the nominal model, it is seen that the basal has the same magnitude for the iZMPC-AV-NOF and the iZMPC-AV-OF.The performances of iZMPC-AV-NOF and iZMPC-AV-OF, both in Scenarios1 and 2, are compared using the Control Variability Grid Analysis (CVGA) [36]. Figure 4 shows the CVGA for the 10 adult patients and the average adult patient. Black color corresponds to the extreme glucose excursions for the 11 patients when the iZMPC-AV-NOF is used, and white color corresponds to the extreme excursions when the iZMPC-AV-OF is used. For Scenario 1, variations of 30% and 50% are represented by the triangle and diamond marker, respectively. For Scenario 2, 30%, and 50% variations are represented by the square and the circle marker, respectively.Figure 4
CVGA comparison of the iZMPC-AV-NOF (black) and iZMPC-AV-OF (white) strategies for 11 adult patients for each variation scenario.Scenario1: the triangle mark is used for the 30% variation, and the diamond for 50% variation.Scenario2: the square mark is used for the 30% variation, and the circle for 50% variation. The average patient in all cases is represented for the filled form.It is evident that the iZMPC-AV-OF achieves better results than the iZMPC-AV-NOF for each case. For Scenario1, it can be seen that with the iZMPC-AV-NOF most patients have excursions in level 2 hyperglycemia, obtaining excursions greater than 400 mg/dl for variations of 50%. Meanwhile, with the iZMPC-AV-OF strategy, the population is moved to lower levels of BG; for variations of 30%, the iZMPC-AV-OF maintains all patients with benign control deviations (zone B), and for 50% variations, 55% of the cases show benign deviations into hyperglycemia (zone upper B). In Scenario 2, more critical results are obtained with the iZMPC-AV-NOF strategy, since there are only 4 cases for variations of 30% that do not present level 2 hypoglycemia. When applying the iZMPC-AV-OF strategy, the BG levels are augmented, moving the population to a safe control (zone A+B), and for variations of 50% there is a safe control in 64% of the cases and an overcorrection of hyperglycemia (zone lower C) for 36% of the cases.Figure5 shows the median of the 11 patients when applying the iZMPC-AV-OF strategy for the same simulation scenarios described above. In addition, the region between the 25th and 75th percentiles of the glycemia is displayed. From this range, it can be observed that in all cases the BG is regulated within the target zone at steady state. The minimum glycemia curves show that there are no level 2 hypoglycemia events. For variations to induce hypoglycemia (Scenario 2), the BG is maintained in the normoglycemia zone for all patients, and then the risk of postprandial hyperglycemia and induced hypoglycemia is eliminated.Figure 5
Glycemia evolution of adult patients when variations in the parameters of 30% and 50% in Scenario1 (to induce hyperglycemia) and Scenario 2 (to induce hypoglycemia) are presented and the iZMPC-AV-OF strategy is used. Continuous lines represent the 50th percentile of the glycemia (median), and the shaded regions are the 25th to 75th percentiles of each case.The glucose control metrics of the 11 patients in each variation scenario to induce hyperglycemia and hypoglycemia are shown in Table3. From the results with the iZMPC-AV-NOF it is possible to reaffirm the need for an extra strategy to compensate the plant-model mismatch. With each variation, the percentage of time in normoglycemia decreases significantly, obtaining 30.0% (IQR 3.3%) of the time in normoglycemia with the variations of 50% in Scenario 1 and 38.2% (45.3%) in Scenario 2, and driving BG to level 2 hyperglycemia and level 2 hypoglycemia during a high percentage of time, 43.4% (10.2%) and 37.5% (35.4%), respectively. On the other hand, a greater percentage of time can be noticed in the target when the iZMPC-AV-OF is applied. In Scenario 1, for variations of the parameters to induce hyperglycemia, it is possible to maintain the glycemia in the zone 70-140 mg/dl for a percentage of time greater than 70%. Only 3.6% (4.2%) of the time is obtained in level 2 hyperglycemia for variations up to 50%, reaching a maximum of 306.04 mg/dl in a patient. In Scenario 2, when variations to induce hypoglycemia are made, a significant increase in the percentage of time in the same zone is noticed, maintaining glycemia in the target more than 90% of the time. In this scenario there are no episodes of level 2 hyperglycemia or hypoglycemia. For greater variations of the parameters in Scenario 1, the percentage of time in the target zone is less than 70%, so that it is considered that variations up to 50% are adequately compensated with the iZMPC-AV-OF strategy.Table 3
Performance comparison of the iZMPC-AV-NOF and iZMPC-AV-OF strategies. Values are medians (interquartile range) unless stated otherwise.
Scenario1: variations to induce hyperglycemia Scenario2: variations to induce hypoglycemia Nominal 10% 30% 50% 10% 30% 50%  Strategy NOF OF NOF OF NOF OF NOF OF NOF OF NOF OF NOF OF  Mean (SD) BG (mg/dl) 114.8 (14.9) 114.0 (14.6) 138.2 (13.1) 119.5 (11.3) 159.5 (25.7) 118.3 (9.9) 253.1 (34.2) 126.3 (8.2) 88.6 (4.6) 110.1 (5.4) 69.5 (8.7) 100.1 (9.2) 51.9 (20.8) 95.4 (4.8)  SD (SD) (mg/dl) 21.2 (10.8) 24.3 (6.8) 33.9 (9.1) 27.6 (6.3) 58.9 (24.0) 35.8 (6.7) 99.9 (29.5) 46.94 (9.5) 17.9 (4.6) 17.0 (3.7) 24.3 (9.0) 13.9 (2.9) 33.4 (15.5) 12.1 (3.7)  Time percentage of BG (%)  < 54 mg/dl 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 22.6 (30.0) 0 (0) 37.5 (35.4) 0 (0) < 6 0 mg/dl 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (3.4) 0 (0) 26.0 (29.9) 0 (0) 47.2 (41.6) 0 (0) < 7 0 mg/dl 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 7.9 (20.4) 0 (0) 38.1 (31.9) 0.3 (0) 61.8 (45.3) 0 (7.3) 70 - 140 mg/dl 83.9 (9.6) 84.3 (6.8) 68.4 (7.9) 78.3 (9.2) 56.6 (14.5) 75.0 (7.4) 7.2 (24.6) 73.0 (8.4) 91.9 (19.6) 91.5 (11.5) 61.9 (31.9) 97.4 (3.3) 38.2 (45.3) 96.8 (7.8) 70 - 180 mg/dl 96.0 (4.0) 98.9 (3.4) 87.9 (13.0) 94.4 (5.3) 68.9 (20.2) 93.6 (6.1) 30.0 (3.3) 82.2 (6.6) 92.1 (20.4) 99.46 (0.9) 61.9 (31.9) 99.7 (0) 38.2 (45.3) 100 (7.3) > 180 mg/dl 4.0 (4.0) 1.1 (3.4) 12.1 (13.0) 5.6 (5.3) 31.1 (20.2) 6.3 (6.8) 70.0 (3.3) 17.8 (6.6) 0 (0) 0 (0.9) 0 (0) 0 (0) 0 (0) 0 (0) > 250 mg/dl 0 (0) 0 (0) 0 (1.3) 0 (0) 7.8 (13.0) 0 (0) 43.4 (10.2) 3.6 (4.2) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) > 300 mg/dl 0 (0) 0 (0) 0 (0) 0 (0) 2.2 (5.7) 0 (0) 32.8 (17.3) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0)  Number of events BG (-)  < 54 mg/dl 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 2(3) 0 (0) 3(2) 0 (0) < 60 mg/dl 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0(1) 0 (0) 3(1) 0 (0) 3(1) 0 (0) < 7 0 mg/dl 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 1(2) 0 (0) 4(1) 0 (0) 3(1) 0(2) > 180 mg/dl 0(1) 1(1) 2(1) 2(1) 3 (0) 1(1) 3 (0) 3 (0) 0 (0) 0(1) 0 (0) 0 (0) 0 (0) 0 (0) > 250 mg/dl 0 (0) 0 (0) 0(1) 0 (0) 1(1) 0 (0) 3 (0) 1(1) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) > 300 mg/dl 0 (0) 0 (0) 0 (0) 0 (0) 1(1) 0 (0) 3(1) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0)To complement the results obtained with the offset-free strategy, a comparison between the dZMPC without offset-free strategy (dZMPC-NOF) and the offset-free dZMPC (dZMPC-OF) is done in order to illustrate the effect in a discrete scheme. For that goal, in Figure6 the same simulation scenarios used in Figure 3 are taken into account. Then it shows a comparison of the evolution of the glycemia and the insulin doses for the average adult patient.Figure 6
Comparison of glycemia evolution and insulin doses for the average adult patient under Scenario2, considering variations of 0%, 30%, and 50%. Dashed lines (blue, purple, and green) are used for dZMPC-NOF while solid lines (yellow, red, and black) are used for dZMPC-OF.As observed with the iZMPC-AV-NOF, the dZMPC-NOF adequately regulates the glycemia when there is no plant-model mismatch, but BG levels fall below 54 mg/dl when there are parameter variations of 30% and 50%. Contrarily, with the dZMPC-OF, the glycemia remains in the desired zone despite the variations; it completely eliminates the offset at steady state and reduces the effect of meal intake. However, it is observed that for 30% variations, BG levels fall almost to the limit of normoglycemia due to the late postprandial effect, which was not evidenced with the iZMPC-AV-OF in Figure3 for the same percentage of variation. Furthermore, it can be seen that the offset in steady state with the dZMPC-NOF is given by a miscalculation of insulin doses, administering a greater amount of basal insulin than that needed for the patient; when applying the offset-free strategy, the amount of basal insulin administered is reduced. Finally, from Figures 3 and 6, the difference of the control action of a discrete system and an impulsive system can be recognized. In the impulsive case, the magnitude of each infusion of insulin is greater than that in the discrete case, because insulin is not delivered continuously but in every period of time.The performance metrics obtained with the dZMPC-NOF and dZMPC-OF are reported in Table4. The improvement obtained by applying the dZMPC-OF in relation to the dZMPC-NOF is evident. For the variations of 50%, in Scenario 1, 8.5% (IQR 32.9%) of the time in normoglycemia results when the dZMPC-NOF is used compared to 71.1% (13.0%) of the time when the dZMPC-OF is used; there is also a big difference in Scenario 2, since there is only 29.7% (15.6%) of time in normoglycemia with the dZMPC-NOF compared to 95.8% (10.3%) with the dZMPC-OF. It is possible to see how the mean of BG is reduced in Scenario 1 and increased in Scenario 2, reaching the desired zone, when applying the dZMPC-OF for each variation.Table 4
Performance comparison of the dZMPC-NOF and dZMPC-OF strategies. Values are medians (interquartile range) unless stated otherwise.
Scenario1: variations to induce hyperglycemia Scenario2: variations to induce hypoglycemia Nominal 10% 30% 50% 10% 30% 50%  Strategy NOF OF NOF OF NOF OF NOF OF NOF OF NOF OF NOF OF  Mean (SD) BG (mg/dl) 111.9 (4.9) 115.6 (6.1) 128.2 (11.8) 127.6 (9.1) 180.4 (31.9) 128.9 (14.2) 291.0 (56.4) 143.5 (10.1) 89.8 (10.4) 109.1 (10.5) 55.4 (27.5) 97.0 (9.1) 21.0 (52.4) 92.5 (7.4) SD (SD) (mg/dl) 21.8 (5.5) 25.7 (4.9) 29.6 (5.7) 31.1 (6.1) 56.6 (12.4) 34.6 (10.1) 109.0 (24.8) 53.1 (7.5) 26.2 (10.6) 18.4 (3.6) 36.3 (32.2) 12.3 (2.3) 54.3 (57.7) 12.5 (2.3)  Time percentage of BG (%)  < 54 mg/dl 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (17.7) 0 (0) 29.8 (33.6) 0 (0) 57.1 (22.1) 0 (0) < 6 0 mg/dl 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 1.2 (23.4) 0 (0) 44.3 (22.4) 0 (0) 67.2 (15.9) 0 (0) < 7 0 mg/dl 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (2.4) 0 (0) 3.5 (7.0) 18.1 (29.2) 0 (0) 56.2 (16.8) 2.8 (3.8) 70.3 (15.6) 4.2 (10.3) 70 - 140 mg/dl 80.7 (8.0) 80.7 (7.2) 70.7 (12.2) 74.4 (8.2) 4.3 (48.4) 68.8 (31.1) 2.1 (27.7) 55.3 (18.2) 76.4 (29.6) 91.5 (11.0) 43.8 (16.8) 96.9 (4.1) 29.7 (15.6) 95.8 (10.3) 70 - 180 mg/dl 96.7 (3.3) 96.9 (4.6) 90.4 (6.3) 91.8 (3.8) 54.4 (24.0) 90.2 (14.5) 8.5 (32.9) 71.1 (13.0) 81.9 (29.2) 100 (0) 43.8 (16.8) 97.3 (3.8) 29.7 (15.6) 95.8 (10.3) > 180 mg/dl 3.3 (3.3) 3.1 (4.6) 8.7 (8.7) 8.2 (3.8) 45.6 (24.0) 9.8 (13.1) 91.5 (32.9) 25.3 (8.8) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) > 250 mg/dl 0 (0) 0 (0) 0 (0) 0 (0) 13.6 (16.1) 0 (0) 69.2 (21.9) 5.6 (5.1) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) > 300 mg/dl 0 (0) 0 (0) 0 (0) 0 (0) 4.8 (5.8) 0 (0) 46.7 (31.2) 0 (0.4) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0)  Number of events BG (-)  < 54 mg/dl 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0(2) 0 (0) 2(2) 0 (0) 3(2) 0 (0) < 60 mg/dl 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 1(2) 0 (0) 3 (0) 0 (0) 4(1) 0 (0) < 7 0 mg/dl 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0(1) 0 (0) 1(1) 2(3) 0 (0) 3(1) 1(1) 4(1) 1(2) > 180 mg/dl 0(1) 1(1) 2(2) 2(1) 3 (0) 2(2) 3 (0) 3 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) > 250 mg/dl 0 (0) 0 (0) 0 (0) 0 (0) 2(1) 0 (0) 3 (0) 1 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) > 300 mg/dl 0 (0) 0 (0) 0 (0) 0 (0) 1(1) 0 (0) 3 (0) 0(1) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0)Lastly, in Figure7 the median of the 11 patients and the region between the 25th and 75th percentiles of the glycemia is displayed when applying the iZMPC-AV-OF strategy considering a robustness test with different meals, parameter variations of 30% in Scenario 2, and sensor noise (i.e., the CGM signal is used instead of the glycemia). The meals consist of 50g at 6:00h, 20g at 9:00h, 90g at 13:00h, 20g at 17:00h, and 70g at 22:00h. The accuracy and dynamics of continuous glucose sensors have been studied in [37]. The noise is then modeled as an ARMA process with a non-Gaussian noise following the Johnson distribution. The percentage of time in the target zone corresponds to 81.56%; it has 15 events of level 1 hypoglycemia and 2 events of level 2 hypoglycemia corresponding to 6.82% and 0.25% of the time, respectively; there is one event of hyperglycemia with a percentage of time of 2.84%. It is noticed that, during the second day, the glycemia reaches its steady state without offset (it achieves the target zone).Figure 7
Glycemia evolution of adult patients when the iZMPC-AV-OF strategy is applied. There are variations of 30% in Scenario2, sensor noise is considered, and five meals are provided. Continuous lines represent the median, and the shaded regions are the 25th to 75th percentiles of each case.


## 3.1. Glycemia Control in T1DM Patients: Comparison of Six MPC Strategies

In this subsection, different MPC strategies for discrete and impulsive systems, applied to the T1DM treatment problem, are compared. These strategies are the standard discrete MPC (dMPC) [24], the discrete ZMPC (dZMPC) with slack variable [26], the discrete ZMPC with artificial variables (dZMPC-AV) [28, 29], and the adaptation of the same formulations for impulsive systems: iMPC [30], iZMPC [22], and iZMPC-AV [10]; the description of each formulation can be seen in the appendix.The three mentioned techniques for discrete and impulsive systems are applied using model (1) and assuming that the plant and the model are equal, i.e., the obtained results are from nominal conditions. As stated before, the formulation of each of the control strategies does not consider the announcement of meals; instead these are considered as unknown disturbances. Therefore, the term Brrt is not considered in the prediction model. The simulated scenario is based on 33 virtual patients from the UVa/Padova simulator, which has 10 children, 10 adolescents, and 10 adults with their average patients, respectively. The initial condition for x1 is set in 120 mg/dl. The objective is to maintain the system in the target zone, even when food intake is entered into the system. The time of the simulation is 24 hours, and meals considered are 50 g of breakfast at 7:00h, lunch of 80g at 12:00h, and dinner of 60g at 21:00h. The prediction horizon N and the matrices Q, R, P are adjusted for each patient using the Nelder-Mead method for autotuning to maximize the time inside the target zone [31]. This tuning procedure is appropriate for in silico environments but other strategies have to be considered for use in real-life conditions such as the one presented in [32]. To all patients, the constraint sets are 0,0,0≤X≤1000,10,10 and 0≤U≤30. In addition, the reference for all strategies is set in 100 mg/dl. The target zone of the glycemia is set as 80≤x1≤140 for all controllers, and the zone called normoglycemia is defined as 70≤x1≤180. As detailed in [33], a value above that range is considered as level 1 hyperglycemia and a value below it as level 1 hypoglycemia, and the critical events x1< 54 mg/dl and x1> 250 mg/dl are denominated as level 2 hypoglycemia and level 2 hyperglycemia, respectively.For all virtual patients evaluated, the performance statistics are computed for every strategy as recommended in [34], taking into account the fact that the results presented here are obtained from simulation and not from clinical trials. Table 2 includes the percentage of time in which glycemia is found both in and out of the normoglycemia zone and the amount of hypoglycemia and hyperglycemia events. Furthermore, the mean and standard deviation (SD) of BG with each controller are calculated. For each outcome the median and the interquartile range (IQR) of the population are shown, except for the mean and SD of BG concentration which are normally distributed so that the average and SD of the 11 patients are shown. The glycemia evolution for the 33 patients using the six different control strategies can be seen in the appendix, which shows that all controllers maintain the patients in the safe zone most of the time without risk of level 2 hypoglycemia. From the results, it is inferred that the iZMPC-AV presents a better performance than the other strategies. It has the highest percentage of precise control, achieving an 83.9% (IQR 11.8%) of time inside the target zone with a median of 1 (IQR 2) event that exceeds 180 mg/dl corresponding to 1.3% (7.9%) of the time. According to this, the iZMPC-AV is selected as the controller to be used in the offset-free control strategy, as was anticipated in early sections.Table 2
Performance comparison of the MPC strategies for discrete and impulsive systems. Values are medians (interquartile range) unless stated otherwise.
Strategy iMPC dMPC iZMPC dZMPC iZMPC-AV dZMPC-AV Mean (SD) BG (mg/dl) 127.8 (10.7) 121.6 (12.2) 129.2 (4.9) 125.5 (8.8) 112.5 (8.7) 112.3 (8.4)  SD (SD) (mg/dl) 36.2 (12.7) 30.9 (11.9) 27.9 (12.1) 31.7 (10.9) 27.5 (13.0) 29.0 (12.8)  Time percentage of BG (%)  < 7 0 mg/dl 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) ∈ 70,140 mg/dl 73.2 (11.4) 76.9 (11.9) 74.1 (7.8) 75.2 (11.3) 83.9 (11.8) 83.8 (12.8) ∈ 70,180 mg/dl 91.5 (14.6) 94.9 (14.8) 94.9 (9.5) 93.9 (11.4) 98.7 (7.9) 97.8 (9.3) > 180 mg/dl 8.5 (14.6) 5.1 (14.8) 5.0 (9.5) 6.1 (11.4) 1.3 (7.9) 2.2 (9.3) > 250 mg/dl 0 (2.0) 0 (1.1) 0 (0) 0 (1.6) 0 (0) 0 (0) > 300 mg/dl 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0)  Number of events BG (-)  < 70 mg/dl 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) > 180 mg/dl 2(2) 2(3) 1(2) 1(2) 1(2) 1(2) > 250 mg/dl 0(1) 0(1) 0 (0) 0(1) 0 (0) 0 (0) > 300 mg/dl 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0)


## 3.2. Sensitivity Analysis

The sensitivity function defined in (9) and given by the differential equation (10) is solved simultaneously with x to obtain its time evolution. The proposed simulation scenario consists of two days, only one of which includes the intake of carbohydrates, with the purpose of demonstrating the sensitivity to parameters during fasting and when there is meal intake. The considered meals are 50g at 7:00h, 80g at 12:00h, and 60g at 21:00h. In Figure 1, the magnitude of the sensitivity of glycemia to pj, denoted as Sx1,pj, is plotted for the average adult as illustration. The parameters are categorized according to the behavior of the sensitivity magnitude, noticing similar responses to p0,p1,p2, and p3 associated with glycemia; p4,p5,p6, and p7 with insulin; and p8,p9,p10, and p11 with intake of carbohydrates.Figure 1
Glycemia sensitivity to each parameter of the model. Top: parameters associated with glycemia. Lower left: parameters associated with insulin. Lower right: parameters associated with intake of carbohydrates.According toSx1,pj, the parameters p0,p1,p2,p3,p5, and p8 are selected as the most influential in blood glucose concentration. Therefore, the simulation of the plant-model mismatch is done when variations in these parameters are considered. The closed-loop system with the iZMPC-AV control strategy is simulated for variations of 10%, 30%, and 50% in the selected parameters and illustrated using the average patient. The behavior of glycemia in the second day confirms the offset provoked by these variations. In Figure 2, the effect of each of these is depicted and the following statements are established:(i)
By increasing the parameterp0 in the plant, less amount of insulin is required. As the prediction model has a lower value of p0, the controller overestimates the necessary dose of insulin and thus the glycemia decreases. If the value of p0 is reduced, the required insulin is higher, but the controller underestimates the amount of insulin, which causes hyperglycemia.(ii)
Variations above the nominal value ofp1 steer to lower levels of BG, obtaining hypoglycemia events. This is because the patient is more sensitive to insulin and therefore requires a smaller amount in each dose, whereas, as the prediction model has a lower value of the parameter, the calculated insulin dose is higher than the actually needed. When insulin action effectiveness parameter is lower than the nominal value, the glycemia increases and hyperglycemia events may occur.(iii)
The variations in parameterp2, corresponding to carbohydrate bioavailability, directly reduce or increase BG levels at the moment of the carbohydrate intake. When the real value of p2 increases, the patient gets more sensitive to meal disturbances and episodes of hyperglycemia in the postprandial period are of greater magnitude. When the real value of p2 is lower than the identified parameter, the effect of meals in the BG levels is reduced.(iv)
Variations inp3 affect glycemia at the time of meals and at steady state. An increase in this parameter leads to an increment of glycemia and it tends to hyperglycemia events. A decrease in the parameter results in a reduction of BG levels and then in hypoglycemia events.(v)
Changes in the parameterp5 are more critical at the time of meals and when its real value is less than the nominal value, leading to hypoglycemia episodes. When the effectiveness of insulin diffusion is greater, the glycemia rises.(vi)
Variations inp8 have an effect on glycemia only at the time of meals. When it is augmented, the glycemia is more sensitive to meal disturbances due to a longer time of diffusion of carbohydrates, leading the patient to have hyperglycemia episodes. By decreasing the value of the parameter, glycemia is reduced in the postprandial period. There is no offset presented because of its dependence on carbohydrate intake.Figure 2
Variation of parametersp0, p1, p2, p3, p5, and p8 in 10%, 30%, and 50% of their nominal values.


## 3.3. Offset-Free Strategy for Glycemia Control with Plant-Model Mismatch

In this subsection, the iZMPC-AV strategy is affected by parameter variations. These variations normally appear in association with exercise, the dawn phenomenon, production of endogenous glucose, among others [35]. As was illustrated in the above subsection, depending on the magnitude of the parameter variation, the iZMPC-AV strategy cannot compensate the disturbance to maintain the adequate control. For that reason, the offset-free iZMPC-AV (iZMPC-AV-OF) strategy is designed. Its performance is evaluated and compared with the iZMPC-AV without the offset-free strategy (iZMPC-AV-NOF) under the next two different scenarios:Scenario 1.
Variation of the parameters to induce hyperglycemia: the parametersp2,p3,p5, and p8 are increased with respect to their nominal values and parameters p0 and p1 are decreased by the same percentages.Scenario 2.
Variation of the parameters to induce hypoglycemia: the parametersp0 and p1 are increased with respect to their nominal values and parameters p2,p3,p5, and p8 are decreased by the same percentages.Asx1 is the only measured state that corresponds to glycemia, the state variables x2, x3, x4, and x5 must be estimated, as well as the disturbance d to correct the steady state mismatch. The simulation consists of 36 hours, starting at 0:00h. The meals considered are a breakfast of 50g at 7:00h, a lunch of 80g at 12:00h, and a dinner of 60g at 21:00h. The control objective is to drive glycemia to the zone between 70 mg/dl and 180 mg/dl and remain there as long as possible. The control strategies do not consider meal announcements.To illustrate the resulting effect when the offset-free strategy is used, Figure3 depicts a comparison of the evolution of the glycemia and the insulin doses for the average adult patient under Scenario 2 considering variations of 0% (nominal situation), 30%, and 50%. The dashed lines are produced when the plant is controlled with the iZMPC-AV-NOF strategy, and the solid lines when the iZMPC-AV-OF strategy is applied. The maximum variation in parameters before the occurrence of hypoglycemia events (or prolonged hyperglycemia events) is established at 50% of the nominal values of the parameters. Meals are only given on the first day to expose the correction at steady state of the second day.Figure 3
Comparison of glycemia evolution and insulin doses for the average adult patient under Scenario2, considering variations of 0%, 30%, and 50%. Dashed lines (blue, purple, and green) represent iZMPC-AV-NOF while solid lines (yellow, red, and black) represent iZMPC-AV-OF.It is noticed that with the iZMPC-AV-NOF the glycemia levels are adequately regulated when the plant and the model are equal, but they decrease to critical BG levels when there are parameter variations of 30% and 50%, maintaining glycemia below 50 mg/dl at meal times. At steady state (from 30:00h to 36:00h), the offset is well appreciated since the glycemia is out of the target zone. This behavior shows the need to apply a strategy that compensates the plant-model mismatch. On the other hand, the iZMPC-AV-OF manages to maintain the glycemia in the desired zone despite the variations. For variations of 30%, the glycemia remains above 80 mg/dl, and for variations of 50% the glycemia remains between 70 and 180 mg/dl and the effect of meal disturbances is reduced. When the iZMPC-AV-OF is applied to the nominal model, the behavior is very similar to that obtained with the iZMPC-AV-NOF; there are negligible differences at meal times, but at steady state the behavior is the same.Regarding the control action for both strategies, when there is a plant-model mismatch and the iZMPC-AV-NOF is applied, it is observed that the offset at steady state is given by a miscalculation of the doses of insulin that are administered. The basal insulin is greater than that needed for the plant with the mismatch. In contrast, when applying the iZMPC-AV-OF strategy, the amount of basal insulin administered is corrected, and the dose is reduced. Finally, for the nominal model, it is seen that the basal has the same magnitude for the iZMPC-AV-NOF and the iZMPC-AV-OF.The performances of iZMPC-AV-NOF and iZMPC-AV-OF, both in Scenarios1 and 2, are compared using the Control Variability Grid Analysis (CVGA) [36]. Figure 4 shows the CVGA for the 10 adult patients and the average adult patient. Black color corresponds to the extreme glucose excursions for the 11 patients when the iZMPC-AV-NOF is used, and white color corresponds to the extreme excursions when the iZMPC-AV-OF is used. For Scenario 1, variations of 30% and 50% are represented by the triangle and diamond marker, respectively. For Scenario 2, 30%, and 50% variations are represented by the square and the circle marker, respectively.Figure 4
CVGA comparison of the iZMPC-AV-NOF (black) and iZMPC-AV-OF (white) strategies for 11 adult patients for each variation scenario.Scenario1: the triangle mark is used for the 30% variation, and the diamond for 50% variation.Scenario2: the square mark is used for the 30% variation, and the circle for 50% variation. The average patient in all cases is represented for the filled form.It is evident that the iZMPC-AV-OF achieves better results than the iZMPC-AV-NOF for each case. For Scenario1, it can be seen that with the iZMPC-AV-NOF most patients have excursions in level 2 hyperglycemia, obtaining excursions greater than 400 mg/dl for variations of 50%. Meanwhile, with the iZMPC-AV-OF strategy, the population is moved to lower levels of BG; for variations of 30%, the iZMPC-AV-OF maintains all patients with benign control deviations (zone B), and for 50% variations, 55% of the cases show benign deviations into hyperglycemia (zone upper B). In Scenario 2, more critical results are obtained with the iZMPC-AV-NOF strategy, since there are only 4 cases for variations of 30% that do not present level 2 hypoglycemia. When applying the iZMPC-AV-OF strategy, the BG levels are augmented, moving the population to a safe control (zone A+B), and for variations of 50% there is a safe control in 64% of the cases and an overcorrection of hyperglycemia (zone lower C) for 36% of the cases.Figure5 shows the median of the 11 patients when applying the iZMPC-AV-OF strategy for the same simulation scenarios described above. In addition, the region between the 25th and 75th percentiles of the glycemia is displayed. From this range, it can be observed that in all cases the BG is regulated within the target zone at steady state. The minimum glycemia curves show that there are no level 2 hypoglycemia events. For variations to induce hypoglycemia (Scenario 2), the BG is maintained in the normoglycemia zone for all patients, and then the risk of postprandial hyperglycemia and induced hypoglycemia is eliminated.Figure 5
Glycemia evolution of adult patients when variations in the parameters of 30% and 50% in Scenario1 (to induce hyperglycemia) and Scenario 2 (to induce hypoglycemia) are presented and the iZMPC-AV-OF strategy is used. Continuous lines represent the 50th percentile of the glycemia (median), and the shaded regions are the 25th to 75th percentiles of each case.The glucose control metrics of the 11 patients in each variation scenario to induce hyperglycemia and hypoglycemia are shown in Table3. From the results with the iZMPC-AV-NOF it is possible to reaffirm the need for an extra strategy to compensate the plant-model mismatch. With each variation, the percentage of time in normoglycemia decreases significantly, obtaining 30.0% (IQR 3.3%) of the time in normoglycemia with the variations of 50% in Scenario 1 and 38.2% (45.3%) in Scenario 2, and driving BG to level 2 hyperglycemia and level 2 hypoglycemia during a high percentage of time, 43.4% (10.2%) and 37.5% (35.4%), respectively. On the other hand, a greater percentage of time can be noticed in the target when the iZMPC-AV-OF is applied. In Scenario 1, for variations of the parameters to induce hyperglycemia, it is possible to maintain the glycemia in the zone 70-140 mg/dl for a percentage of time greater than 70%. Only 3.6% (4.2%) of the time is obtained in level 2 hyperglycemia for variations up to 50%, reaching a maximum of 306.04 mg/dl in a patient. In Scenario 2, when variations to induce hypoglycemia are made, a significant increase in the percentage of time in the same zone is noticed, maintaining glycemia in the target more than 90% of the time. In this scenario there are no episodes of level 2 hyperglycemia or hypoglycemia. For greater variations of the parameters in Scenario 1, the percentage of time in the target zone is less than 70%, so that it is considered that variations up to 50% are adequately compensated with the iZMPC-AV-OF strategy.Table 3
Performance comparison of the iZMPC-AV-NOF and iZMPC-AV-OF strategies. Values are medians (interquartile range) unless stated otherwise.
Scenario1: variations to induce hyperglycemia Scenario2: variations to induce hypoglycemia Nominal 10% 30% 50% 10% 30% 50%  Strategy NOF OF NOF OF NOF OF NOF OF NOF OF NOF OF NOF OF  Mean (SD) BG (mg/dl) 114.8 (14.9) 114.0 (14.6) 138.2 (13.1) 119.5 (11.3) 159.5 (25.7) 118.3 (9.9) 253.1 (34.2) 126.3 (8.2) 88.6 (4.6) 110.1 (5.4) 69.5 (8.7) 100.1 (9.2) 51.9 (20.8) 95.4 (4.8)  SD (SD) (mg/dl) 21.2 (10.8) 24.3 (6.8) 33.9 (9.1) 27.6 (6.3) 58.9 (24.0) 35.8 (6.7) 99.9 (29.5) 46.94 (9.5) 17.9 (4.6) 17.0 (3.7) 24.3 (9.0) 13.9 (2.9) 33.4 (15.5) 12.1 (3.7)  Time percentage of BG (%)  < 54 mg/dl 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 22.6 (30.0) 0 (0) 37.5 (35.4) 0 (0) < 6 0 mg/dl 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (3.4) 0 (0) 26.0 (29.9) 0 (0) 47.2 (41.6) 0 (0) < 7 0 mg/dl 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 7.9 (20.4) 0 (0) 38.1 (31.9) 0.3 (0) 61.8 (45.3) 0 (7.3) 70 - 140 mg/dl 83.9 (9.6) 84.3 (6.8) 68.4 (7.9) 78.3 (9.2) 56.6 (14.5) 75.0 (7.4) 7.2 (24.6) 73.0 (8.4) 91.9 (19.6) 91.5 (11.5) 61.9 (31.9) 97.4 (3.3) 38.2 (45.3) 96.8 (7.8) 70 - 180 mg/dl 96.0 (4.0) 98.9 (3.4) 87.9 (13.0) 94.4 (5.3) 68.9 (20.2) 93.6 (6.1) 30.0 (3.3) 82.2 (6.6) 92.1 (20.4) 99.46 (0.9) 61.9 (31.9) 99.7 (0) 38.2 (45.3) 100 (7.3) > 180 mg/dl 4.0 (4.0) 1.1 (3.4) 12.1 (13.0) 5.6 (5.3) 31.1 (20.2) 6.3 (6.8) 70.0 (3.3) 17.8 (6.6) 0 (0) 0 (0.9) 0 (0) 0 (0) 0 (0) 0 (0) > 250 mg/dl 0 (0) 0 (0) 0 (1.3) 0 (0) 7.8 (13.0) 0 (0) 43.4 (10.2) 3.6 (4.2) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) > 300 mg/dl 0 (0) 0 (0) 0 (0) 0 (0) 2.2 (5.7) 0 (0) 32.8 (17.3) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0)  Number of events BG (-)  < 54 mg/dl 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 2(3) 0 (0) 3(2) 0 (0) < 60 mg/dl 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0(1) 0 (0) 3(1) 0 (0) 3(1) 0 (0) < 7 0 mg/dl 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 1(2) 0 (0) 4(1) 0 (0) 3(1) 0(2) > 180 mg/dl 0(1) 1(1) 2(1) 2(1) 3 (0) 1(1) 3 (0) 3 (0) 0 (0) 0(1) 0 (0) 0 (0) 0 (0) 0 (0) > 250 mg/dl 0 (0) 0 (0) 0(1) 0 (0) 1(1) 0 (0) 3 (0) 1(1) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) > 300 mg/dl 0 (0) 0 (0) 0 (0) 0 (0) 1(1) 0 (0) 3(1) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0)To complement the results obtained with the offset-free strategy, a comparison between the dZMPC without offset-free strategy (dZMPC-NOF) and the offset-free dZMPC (dZMPC-OF) is done in order to illustrate the effect in a discrete scheme. For that goal, in Figure6 the same simulation scenarios used in Figure 3 are taken into account. Then it shows a comparison of the evolution of the glycemia and the insulin doses for the average adult patient.Figure 6
Comparison of glycemia evolution and insulin doses for the average adult patient under Scenario2, considering variations of 0%, 30%, and 50%. Dashed lines (blue, purple, and green) are used for dZMPC-NOF while solid lines (yellow, red, and black) are used for dZMPC-OF.As observed with the iZMPC-AV-NOF, the dZMPC-NOF adequately regulates the glycemia when there is no plant-model mismatch, but BG levels fall below 54 mg/dl when there are parameter variations of 30% and 50%. Contrarily, with the dZMPC-OF, the glycemia remains in the desired zone despite the variations; it completely eliminates the offset at steady state and reduces the effect of meal intake. However, it is observed that for 30% variations, BG levels fall almost to the limit of normoglycemia due to the late postprandial effect, which was not evidenced with the iZMPC-AV-OF in Figure3 for the same percentage of variation. Furthermore, it can be seen that the offset in steady state with the dZMPC-NOF is given by a miscalculation of insulin doses, administering a greater amount of basal insulin than that needed for the patient; when applying the offset-free strategy, the amount of basal insulin administered is reduced. Finally, from Figures 3 and 6, the difference of the control action of a discrete system and an impulsive system can be recognized. In the impulsive case, the magnitude of each infusion of insulin is greater than that in the discrete case, because insulin is not delivered continuously but in every period of time.The performance metrics obtained with the dZMPC-NOF and dZMPC-OF are reported in Table4. The improvement obtained by applying the dZMPC-OF in relation to the dZMPC-NOF is evident. For the variations of 50%, in Scenario 1, 8.5% (IQR 32.9%) of the time in normoglycemia results when the dZMPC-NOF is used compared to 71.1% (13.0%) of the time when the dZMPC-OF is used; there is also a big difference in Scenario 2, since there is only 29.7% (15.6%) of time in normoglycemia with the dZMPC-NOF compared to 95.8% (10.3%) with the dZMPC-OF. It is possible to see how the mean of BG is reduced in Scenario 1 and increased in Scenario 2, reaching the desired zone, when applying the dZMPC-OF for each variation.Table 4
Performance comparison of the dZMPC-NOF and dZMPC-OF strategies. Values are medians (interquartile range) unless stated otherwise.
Scenario1: variations to induce hyperglycemia Scenario2: variations to induce hypoglycemia Nominal 10% 30% 50% 10% 30% 50%  Strategy NOF OF NOF OF NOF OF NOF OF NOF OF NOF OF NOF OF  Mean (SD) BG (mg/dl) 111.9 (4.9) 115.6 (6.1) 128.2 (11.8) 127.6 (9.1) 180.4 (31.9) 128.9 (14.2) 291.0 (56.4) 143.5 (10.1) 89.8 (10.4) 109.1 (10.5) 55.4 (27.5) 97.0 (9.1) 21.0 (52.4) 92.5 (7.4) SD (SD) (mg/dl) 21.8 (5.5) 25.7 (4.9) 29.6 (5.7) 31.1 (6.1) 56.6 (12.4) 34.6 (10.1) 109.0 (24.8) 53.1 (7.5) 26.2 (10.6) 18.4 (3.6) 36.3 (32.2) 12.3 (2.3) 54.3 (57.7) 12.5 (2.3)  Time percentage of BG (%)  < 54 mg/dl 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (17.7) 0 (0) 29.8 (33.6) 0 (0) 57.1 (22.1) 0 (0) < 6 0 mg/dl 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 1.2 (23.4) 0 (0) 44.3 (22.4) 0 (0) 67.2 (15.9) 0 (0) < 7 0 mg/dl 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (2.4) 0 (0) 3.5 (7.0) 18.1 (29.2) 0 (0) 56.2 (16.8) 2.8 (3.8) 70.3 (15.6) 4.2 (10.3) 70 - 140 mg/dl 80.7 (8.0) 80.7 (7.2) 70.7 (12.2) 74.4 (8.2) 4.3 (48.4) 68.8 (31.1) 2.1 (27.7) 55.3 (18.2) 76.4 (29.6) 91.5 (11.0) 43.8 (16.8) 96.9 (4.1) 29.7 (15.6) 95.8 (10.3) 70 - 180 mg/dl 96.7 (3.3) 96.9 (4.6) 90.4 (6.3) 91.8 (3.8) 54.4 (24.0) 90.2 (14.5) 8.5 (32.9) 71.1 (13.0) 81.9 (29.2) 100 (0) 43.8 (16.8) 97.3 (3.8) 29.7 (15.6) 95.8 (10.3) > 180 mg/dl 3.3 (3.3) 3.1 (4.6) 8.7 (8.7) 8.2 (3.8) 45.6 (24.0) 9.8 (13.1) 91.5 (32.9) 25.3 (8.8) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) > 250 mg/dl 0 (0) 0 (0) 0 (0) 0 (0) 13.6 (16.1) 0 (0) 69.2 (21.9) 5.6 (5.1) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) > 300 mg/dl 0 (0) 0 (0) 0 (0) 0 (0) 4.8 (5.8) 0 (0) 46.7 (31.2) 0 (0.4) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0)  Number of events BG (-)  < 54 mg/dl 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0(2) 0 (0) 2(2) 0 (0) 3(2) 0 (0) < 60 mg/dl 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 1(2) 0 (0) 3 (0) 0 (0) 4(1) 0 (0) < 7 0 mg/dl 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0(1) 0 (0) 1(1) 2(3) 0 (0) 3(1) 1(1) 4(1) 1(2) > 180 mg/dl 0(1) 1(1) 2(2) 2(1) 3 (0) 2(2) 3 (0) 3 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) > 250 mg/dl 0 (0) 0 (0) 0 (0) 0 (0) 2(1) 0 (0) 3 (0) 1 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) > 300 mg/dl 0 (0) 0 (0) 0 (0) 0 (0) 1(1) 0 (0) 3 (0) 0(1) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0) 0 (0)Lastly, in Figure7 the median of the 11 patients and the region between the 25th and 75th percentiles of the glycemia is displayed when applying the iZMPC-AV-OF strategy considering a robustness test with different meals, parameter variations of 30% in Scenario 2, and sensor noise (i.e., the CGM signal is used instead of the glycemia). The meals consist of 50g at 6:00h, 20g at 9:00h, 90g at 13:00h, 20g at 17:00h, and 70g at 22:00h. The accuracy and dynamics of continuous glucose sensors have been studied in [37]. The noise is then modeled as an ARMA process with a non-Gaussian noise following the Johnson distribution. The percentage of time in the target zone corresponds to 81.56%; it has 15 events of level 1 hypoglycemia and 2 events of level 2 hypoglycemia corresponding to 6.82% and 0.25% of the time, respectively; there is one event of hyperglycemia with a percentage of time of 2.84%. It is noticed that, during the second day, the glycemia reaches its steady state without offset (it achieves the target zone).Figure 7
Glycemia evolution of adult patients when the iZMPC-AV-OF strategy is applied. There are variations of 30% in Scenario2, sensor noise is considered, and five meals are provided. Continuous lines represent the median, and the shaded regions are the 25th to 75th percentiles of each case.


## 4. Conclusions

A performance comparison of the different MPC strategies for discrete and impulsive systems, applied to the T1DM treatment problem, was done in a simulation environment by means of the UVa/Padova simulator. All these strategies were tested against meal disturbances, without announcement to the control strategy. The use of ZMPC strategies seems to be adequate for the T1DM treatment, since it applies enough insulin to maintain the blood glucose in the target zone, eliminating the need to administer greater amount of insulin in each dose to reach a set-point and fulfil constraints. The iZMPC-AV presents the least percentage of time out of the normoglycemia zone. One of its main benefits is that it has an enlarged domain of attraction due to the artificial intermediary variables and the application of the basal as small boluses instead of a continuous dose because of the impulsive scheme.The main goal of this work is to propose a control strategy based on the iZMPC-AV for the glycemic regulation problem in type 1 diabetic patients, including a plant-model mismatch caused by variations in the parameters. The iZMPC-AV-OF strategy is then used to treat the problem. The main achievement of this control application is to maintain the blood glucose levels within the target zone of normoglycemia, decreasing hypoglycemia and hyperglycemia events caused by parameter variations and carbohydrate intake.By simulating two scenarios for 11 adult virtual patients, satisfactory performances were achieved. It was shown that the dZMPC-OF and iZMPC-AV-OF strategies maintain the glycemia in the normal range (70-180 mg/dl) for variations up to 50%. Hypoglycemia episodes were successfully avoided and hyperglycemia events were minimized. In Scenario1, when the variations are made to induce hyperglycemia, the percentage time in the target exceeds 70% by applying the iZMPC-AV-OF strategy, while in Scenario 2, when variations to induce hypoglycemia are made, the percentage of time increases over 90%. This is because of two main reasons: (i) variations in the parameters, in Scenario 2, associated with carbohydrates favor the reduction of hyperglycemia peaks in postprandial times, and (ii) the controller manages to reduce the episodes of hypoglycemia, maintaining glycemia in the target zone for approximately 90% of the time.This work presents a first step to counteract the plant-model mismatch inherent in a biomedical problem such as the treatment of type 1 diabetes mellitus. The results are obtained in a limited simulation environment built with the UVa/Padova simulator, which tend to have more contra-regulation on low glycemic values than that observed in real subjects. Although the results are very good, it is necessary to remember that the simulated scenarios do not exactly describe real-life conditions. For further works more realistic scenarios will be taken into account, including temporary disturbances in the parameters given by exercise, the dawn phenomenon, insulin sensitivity, and the time constant of insulin absorption. In addition, an analysis of the noise given by the glucose sensor must be considered in more detail, which, having a nonzero-mean, can alter the elimination of the offset.



---

*Source: 2640405-2019-05-28.xml*