Institutionen För Systemteknik Department Of Electrical Engineering

Transcript

Institutionen för systemteknik Department of Electrical Engineering Examensarbete Model Predictive Control with Invariant Sets in Artificial Pancreas for Type 1 Diabetes Mellitus Examensarbete utfört i Elektroteknik vid Tekniska högskolan vid Linköpings universitet av Andreas Svensson LiTH-ISY-EX--13/4699--SE Linköping 2013 Department of Electrical Engineering Linköpings universitet SE-581 83 Linköping, Sweden Linköpings tekniska högskola Linköpings universitet 581 83 Linköping Model Predictive Control with Invariant Sets in Artificial Pancreas for Type 1 Diabetes Mellitus Examensarbete utfört i Elektroteknik vid Tekniska högskolan vid Linköpings universitet av Andreas Svensson LiTH-ISY-EX--13/4699--SE Handledare: Ravi Gondhalekar Department of Chemical Engineering, University of Californa Santa Barbara Isak Nielsen isy, Linköpings universitet Examinator: Daniel Axehill isy, Linköpings universitet Linköping, 20 juni 2013 Avdelning, Institution Division, Department Datum Date Avdelningen för reglerteknik Department of Electrical Engineering SE-581 83 Linköping 2013-06-20 Språk Language Rapporttyp Report category ISBN  Svenska/Swedish  Licentiatavhandling ISRN  Engelska/English   Examensarbete   C-uppsats  D-uppsats — LiTH-ISY-EX--13/4699--SE Serietitel och serienummer Title of series, numbering  Övrig rapport  ISSN —  URL för elektronisk version http://www.ep.liu.se Titel Title Modellbaserad prediktionsreglering med invarianta mängder av artificiella pankreas för diabetes typ 1 Model Predictive Control with Invariant Sets in Artificial Pancreas for Type 1 Diabetes Mellitus Författare Author Andreas Svensson Sammanfattning Abstract This thesis deals with Model Predictive Control (mpc) for artificial pancreas for Type 1 Diabetes Mellitus (t1dm) patients. A control strategy exploiting invariant sets in mpc for blood glucose level control is developed, to the authors knowledge for the first time. The work includes various types of invariant sets relevant for the artificial pancreas problem, and different ways to incorporate them into the mpc strategy. The work is an extension to the zone mpc controller for artificial pancreas developed at University of California Santa Barbara and Sansum Diabetes Research Institute (Grosman et al. [2010]; van Heusden et al. [2012]). The evaluation of the proposed control strategy is done in silico in the fda (U.S. Food and Drug Administration) approved metabolic simulator (Kovatchev et al. [2009]). The trials show some promising results in terms of more rapid meal responses and decreased variability between the subjects than the zone mpc. An attempt to robust control employing invariant sets proved to be less promising in the evaluations. The results indicate that the direct application of known robust control techniques is not appropriate, and that more appropriate robust control techniques must be searched for, or developed, more specific to the artificial pancreas control. Altogether, this thesis pinpoints a possible future direction of artificial pancreas control design, with mpc based on invariant sets. Nyckelord Keywords MPC, Invariant sets, Artificial Pancreas, Type 1 Diabetes Mellitus Sammanfattning Det här examensarbetet handlar om modellbaserad prediktionsreglering (mpc) av artificiella pankreas för typ 1 diabetes mellitus (t1dm) patienter. I arbetet utvecklas, till författarens kännedom för första gången, en mpc-strategi med invarianta mängder för reglering av blodsockernivå. I arbetet ingår olika typer av invarianta mängder för artificiella pankreas, och olika sätt att införliva dem i mpc. Arbetet är en utvidgning av den zon-mpc för artificiella pankreas som utvecklats vid University of Californa Santa Barbara och Sansum Diabetes Research Institute (Grosman et al. [2010], van Heusden et al. [2012]). Utvärderingen av den föreslagna reglerstrategin görs in silico i en metabolisk simulator (Kovatchev et al. [2009]) godkänd av U.S. Food and Drug Administration. Försöken visar lovande resultat i form av snabb respons på måltider och minskad variation mellan individer jämfört med nämnda zon-mpc. Ett försök till robust reglering med invarianta mängder visade sig dock vara mindre lovande i utvärderingarna. Resultaten antyder att en rättfram tillämpning av kända metoder för robust reglering inte är ändamålsenlig för den här applikationen, och att metoder mer relevanta för artificiella pankreas behöver hittas eller utvecklas. Sammantaget lyfter det här examensarbetet fram en möjlig framtida inriktning för reglering av artificiella pankreas med mpc baserad på invarianta mängder. iii Abstract This thesis deals with Model Predictive Control (mpc) for artificial pancreas for Type 1 Diabetes Mellitus (t1dm) patients. A control strategy exploiting invariant sets in mpc for blood glucose level control is developed, to the authors knowledge for the first time. The work includes various types of invariant sets relevant for the artificial pancreas problem, and different ways to incorporate them into the mpc strategy. The work is an extension to the zone mpc controller for artificial pancreas developed at University of California Santa Barbara and Sansum Diabetes Research Institute (Grosman et al. [2010]; van Heusden et al. [2012]). The evaluation of the proposed control strategy is done in silico in the fda (U.S. Food and Drug Administration) approved metabolic simulator (Kovatchev et al. [2009]). The trials show some promising results in terms of more rapid meal responses and decreased variability between the subjects than the zone mpc. An attempt to robust control employing invariant sets proved to be less promising in the evaluations. The results indicate that the direct application of known robust control techniques is not appropriate, and that more appropriate robust control techniques must be searched for, or developed, more specific to the artificial pancreas control. Altogether, this thesis pinpoints a possible future direction of artificial pancreas control design, with mpc based on invariant sets. v Acknowledgments First of all, I would like to thank Dr. Ravi Gondhalekar for being a truly excellent supervisor to me. I really appreciate all discussions we have had, your patience, advice and all challenges you have encouraged me to tackle! I would also like to thank Prof. Francis J. Doyle III and Dr. Eyal Dassau for taking me up in the group and letting me work on the artificial pancreas project, and all my fellow Doyle group members for the atmosphere in the lab that made me feel at home. I am also grateful to the people at Sansum for giving me a glimpse of the very interesting world of real (i.e., not simulated) patients and trials. Dr. Daniel Axehill deserves a big thank for initially letting me in touch with Ravi and the Doyle group, and for your advice and proofreading of this thesis together with Isak Nielsen. I would also like to thank Sanna. Although your contributions to this thesis are rather limited, your contributions to my life during this period have been the more important. Finally, I am also very grateful to my family for their support during all my years of studies. Santa Barbara, May 2013 Andreas Svensson vii Contents Notation 1 Introduction 1.1 Problem formulation 1.1.1 Background . 1.1.2 Objective . . . 1.1.3 Methods . . . 1.1.4 Limitations . . 1.2 Outline of the thesis . xi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1 1 2 2 2 3 2 Control theory 2.1 Dynamical systems . . . . . . . . . . . . . . . 2.1.1 Discrete time state space form . . . . . 2.1.2 Discrete time transfer function . . . . 2.1.3 State estimation . . . . . . . . . . . . . 2.2 Model Predictive Control . . . . . . . . . . . . 2.2.1 The standard MPC formulation . . . . 2.2.2 Zone MPC . . . . . . . . . . . . . . . . 2.3 Invariant sets . . . . . . . . . . . . . . . . . . . 2.3.1 Positively invariant sets . . . . . . . . . 2.3.2 Controlled invariant sets . . . . . . . . 2.3.3 Robust controlled invariant sets . . . . 2.3.4 Periodic controlled invariant sets . . . 2.3.5 Maximum invariant sets . . . . . . . . 2.4 Mini example of invariant set in MPC control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 6 6 7 9 9 10 12 12 13 15 16 16 17 17 3 Diabetes and the artificial pancreas 3.1 Diabetes . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Insulin treatment . . . . . . . . . . . . . . . . . . . 3.3 Continuous Glucose Monitoring . . . . . . . . . . . 3.4 Metabolic system simulator . . . . . . . . . . . . . 3.5 System identification for an insulin-glucose model 3.6 Artificial pancreas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 23 24 24 25 25 26 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x CONTENTS 3.6.1 3.6.2 3.6.3 3.6.4 Meals . . . . . . . . . . . . . . . . . . . . . . . Safety systems and hypoglycemia treatments Insulin on board constraints . . . . . . . . . . Night time zone . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 28 28 28 4 MPC with invariant set constraints for an artificial pancreas 4.1 Invariant sets for the insulin-glucose model . . . . . . . . 4.2 Robust invariant set for the insulin-glucose model . . . . 4.2.1 Measurement errors in the state space . . . . . . . 4.2.2 Robustification against measurement errors . . . . 4.2.3 Meals . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.4 Plant-model mismatch . . . . . . . . . . . . . . . . 4.3 Implementation of set constraints in MPC . . . . . . . . . 4.4 In silico trials . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4.1 Test protocol . . . . . . . . . . . . . . . . . . . . . . 4.4.2 Different controller setups in the trials . . . . . . . 4.4.3 Subjects in the trial . . . . . . . . . . . . . . . . . . 4.4.4 In silico results . . . . . . . . . . . . . . . . . . . . 4.4.5 Evaluation of alternative controllers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 31 32 32 33 34 35 35 37 37 37 38 38 42 5 Concluding comments 5.1 Comments and conclusions . 5.1.1 In silico results . . . . 5.1.2 Invariant set approach 5.2 Further work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 47 47 47 48 A Invariant set mathematics A.1 Controlled invariant set . . . . . . . . . . . . . . . . . . . . . . . . . 51 51 B Implementation B.1 Software . . . . . . . . . . . . . . . . . . . . B.2 Matlab code . . . . . . . . . . . . . . . . . . B.2.1 Code for invariant set computations B.2.2 Code for Example 2.10 . . . . . . . . 53 53 53 54 56 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C Experimental invariant set illustration 63 Bibliography 65 Index 69 Notation General Symbol R ∀ Explanation Set of real numbers For all Abbreviations Symbol arx cgm cvga csii mpc tdi t1dm zmpc Explanation Auto-Regressive model with eXogenous input Continuous Glucose Monitoring Control-Variability Grid Analysis Continuous Subcutaneous Insulin Infusion Model Predictive Control Total Daily Insulin Type 1 Diabetes Mellitus zone-mpc xi 1 Introduction Artificial Pancreas for Type 1 Diabetes Mellitus (t1dm) is an evolving field, with numerous projects going on around the world. This thesis deals with invariant sets for robust model predictive control (mpc) for an artificial pancreas. The work behind this thesis is done in the group1 of Prof. Francis J. Doyle III at the University of California Santa Barbara (ucsb), in collaboration with Sansum Diabetes Research Institute2 , where the development of artificial pancreas has been a topic for several years (Grosman et al. [2010]; van Heusden et al. [2012]). 1.1 Problem formulation In this section, an overview of this thesis in terms of background, objective, methods, and limitations is given. 1.1.1 Background Type 1 Diabetes Mellitus (t1dm) is an autoimmune chronic metabolic disease, characterized by the loss of endogenous secretion of insulin from the pancreatic beta-cells. The insulin is crucial in the control of the blood glucose level in the human body. Without treatment, the consequences of the heighten blood glucose level are life-threatening. About 348 million people worldwide had some kind of Diabetes in 2008 (Danaei et al. [2011]), of which approximately 10% (Eiselein et al. [2004]) can be assumed being t1dm. 1 www.thedoylegroup.org 2 www.sansum.org 1 2 1 Introduction The traditional treatment of t1dm involves, since the discovery of the insulin in the 1920s, manually injected insulin (Skyler [2012]), to compensate for the absence of insulin-production in the beta-cells. In the last decades, the use of insulin pumps has evolved as an alternative to manual injections (Skyler [2012]). A fully automatic system that delivers insulin via a pump based on continuous measurement of the blood glucose level is called an artificial pancreas. A wellfunctioning artificial pancreas could improve the quality of life, and lengthen the life-expectancy, of people with t1dm. Besides the practical assistance with the delivery itself, an artificial pancreas could possibly decrease a possible mental stress over blood glucose control for a t1dm patient3 . 1.1.2 Objective The objective for this thesis work is to investigate possible ways to utilize the state space information in the predictions in the mpc for artificial pancreas done in the ucsb/Sansum group4 . This is an attempt to improve the control and handle some of the uncertainties in the artificial pancreas control problem in a mathematically rigorous way. 1.1.3 Methods The methods used is mpc, extended with various invariant set methods from the literature (Blanchini [1999]; Blanchini and Ukovich [1993]; Dórea and Hennet [1999]; Gondhalekar et al. [2013], see Section 2.3 and 2.4). These methods are implemented in the zone-mpc (zmpc) controller for artificial pancreas developed by the ucsb/ Sansum group, and evaluated in silico 5 in the UVa/Padova simulator (Kovatchev et al. [2009], see Section 3.4). 1.1.4 Limitations There are several limitations for this works of which the most important might be: • The controller evaluation is only done in silico • The tuning of the various control parameters is not addressed systematically, but the focus is rather on qualitative aspects of the control behavior than on parameter optimization • Only the zmpc controller is modified, and no tuning or changing of the state estimator, safety systems or other components in the ucsb/Sansum group system is addressed • All work is done on the current insulin-glucose model (van Heusden et al. [2012], see Section 3.5), and no attempts on improving or finding new models are done. 3 See Greenberg [2013] for an anecdotal example. 4 In the current version, only the predicted outputs are used in the control design. 5 i.e., in simulations 1.2 Outline of the thesis 1.2 3 Outline of the thesis The outline of the thesis is as follows • This is Chapter 1 • Chapter 2 gives a brief introduction to control theory and automatic control • A brief introduction to t1dm is given in Chapter 3, together with a summary of related recent work on artificial pancreas • Chapter 4 presents and discusses design choices in the controller design, and presents the in silico trial results • Chapter 5 gives some comments on the results and suggests directions for future work • Some mathematical details on the invariant sets are given in Appendix A • Appendix B provides some details on the implementations of various computations, including Matlab code examples • Appendix C presents an experimental way for illustration of invariant sets for a special class of systems, relevant for the current application. 2 Control theory Figure 2.1: The notion of a ‘system’ The core idea of automatic control is to engineer signals such that a system behaves in a desired way. Some central parts of control theory relevant for this thesis will be presented in this chapter; first the idea and mathematical modeling of a system is briefly presented, followed by a short introduction to the control technique Model Predictive Control (mpc). In the second last section of this chapter, an introduction to invariant sets is found, and an example on how to combine the mpc and the invariant sets methods for improved control is given in the last section. A few references to recommended literature on automatic control are Glad and Ljung [2006] (in Swedish), Åström and Murray [2010] and Glad and Ljung [2000]. For more details on the mathematical modeling and analysis of dynamical systems, Rugh [1996] is recommended. 5 6 2.1 2 Control theory Dynamical systems The term ‘system’ is used to describe situations as shown in Figure 2.1. There are signals which, e.g., can be electronic, medical or chemical, which are named input signals. There are also signals, of possibly different sort than the input signals, which are named output signals. Causal relations between the input and output signals are assumed, where the output signals are depending on the input signals. The terminology can be used even if there is only a limited knowledge of these relations available. If an input signal is used to control the system, it can also be named control signal. If an input signal, on the other hand, is beyond control, it can be named disturbance signal. When a system is subject to control, it is also sometimes called a plant. The notion of dynamical systems is to stress a (possible) dependence between the current output and the previous1 inputs (in contrast to a static system, where the current output depends only on the current, and not the previous, inputs). 2.1.1 Discrete time state space form If the relation between the output at time t and the input at2 time t, t − 1, . . . can be described by a mathematical model with the structure of (2.1), the system is called a linear discrete time system. x(t + 1) =A(t)x(t) + B(t)u(t) (2.1a) y(t) =C(t)x(t) + D(t)u(t) (2.1b) Here, the input is denoted u(t) and the output y(t). The input, the output, the variable x(t) are scalars or vectors. The (possibly time-varying) parameters A(t), B(t), C(t) and D(t) are matrices of appropriate dimension (possibly scalar). The structure of (2.1) is called state space form, which is related to the state variable x(t). The values of x(t) might, but do not have to, have physical interpretations for the system described by the model (such as speed, acceleration, etc.). A crucial point with the state space form is that the state variables x(t) contain all information about the current state of the system, i.e., as long as the current state is known, the input and output histories are of no further use in telling the current situation in the system3 . Consider, e.g., a car. Given only the current position, it is hard to make any claims about the position of the car in the near future. However, given the position and the velocity, the possibilities of making accurate claims about the position of 1 An even more ambitious notion would be causal dynamical systems to stress the dependence on current and previous, but not future, inputs. This will however be implicitly assumed in the sequel. 2 For the sake of a simple and intuitive notation, a sample-period of 1 is assumed. 3 This has big similarities with the markov property of stochastic processes, see, e.g., Yates and Goodman [1999]. 2.1 7 Dynamical systems the car in the near future increases drastically. The velocity and position can be thought of as state variables for this case. The linear state space form in (2.1) can be generalized to cover also systems with nonlinear dynamics. To obtain a more general form of the equations, the (linear) matrix multiplications and additions on the right hand sides in (2.1) are simply replaced by the more general function symbols f and g. x(t + 1) =f (x(t), u(t)) (2.2a) y(t) =g(x(t), u(t)) (2.2b) Of course, the matrix multiplications and additions in (2.1) can be considered as a special case of (2.2). 2.1.2 Discrete time transfer function If there is only one4 input and one output to the system, the input–output relation expressed on the form (2.1) can also be expressed on the following form (see, e.g., Rugh [1996]) eliminating the state variables: y(t) + a1 y(t − 1) + . . . + an y(t − n) = b0 u(t) + b1 u(t − 1) + . . . + bm u(t − m) (2.3) (where the constants a1 , . . . , an and b0 , . . . , bm might be time-varying). From (2.3) it is clear that the output y(t) is calculated using information of the current and old input signals u(t), u(t − 1), . . . . This is to compare with the state space form (2.1), where the current input signal u(t) and the current states x(t) are sufficient for calculating the output y(t). In the case of time-invariant constants a1 , . . . , an , b0 , . . . , bm , the z-transform can be used (see, e.g., Rugh [1996]) to rewrite (2.3) into an equivalent, but slightly more compact, form. Via the step Y (z) + a1 z −1 Y (z) + . . . + an z −n Y (z) = b0 U (z) + b1 · z −1 U (z) + . . . + bm y −m (z) (2.4a) (where Y (z) is the z-transform of y(t), and similar for U (z)) the equation Y (z) = b0 + b1 z −1 + . . . + bm z −m U (z). 1 + a1 z −1 + . . . + an z −n (2.4b) can be obtained. The form (2.4b) is usually referred to as the transfer function from the input U (z) to the output Y (z). The values of z for which the denominator of (2.4b) equals 0 is denoted as the poles of the system. For the same system described on the state space form, all 4 There exists a corresponding theory also for systems with multiple inputs and outputs, see for instance Rugh [1996]. 8 2 Control theory poles are found also as eigenvalues to the matrix5 A. The values of z for which the numerator becomes zero are denoted the zeros of the system. The theory of state space forms (2.1) and transfer functions (2.4b) is extensive and is only very briefly summarized in this section. It is however enough for the scope of this thesis to note that there exists a detailed theory on the relationship and conversion between the forms, and refer to, e.g., Rugh [1996] for more details on the topic. 2.1 Example: Transfer function The relationship between subcutaneous injected insulin and the blood glucose level in the human body can be viewed as a system with one input signal (the injected insulin) and one output signal (the blood glucose level). This relationship has been identified and approximated as (van Heusden et al. [2012]): Y (z) = −0.01 U (z) (z − 0.98)(z − 0.965)2 (2.5) where U (z) is the z-transform of the injected insulin and Y (z) is the z-transform of the blood glucose level6 . In the previous text the sample-period has been set to 1 for the sake of simplicity. For this model, the sample-period is 5 minutes. A more thorough walkthrough of the model can be found in Section 3.5. 2.2 Example: State space description The input–output relation described by the transfer function (2.5) in Example 2.1 can be described on a state space form:   x1 (t + T )   x2 (t + T )  x3 (t + T )     2.91    =  1   0 y(t) =  0 −2.823 0.913 0 0 1 0     x1 (t)  0 1  x2 (t)    x3 (t)       x1 (t)   −0.01    x2 (t)  +  0  u(t)      x3 (t) 0 (2.6a) (2.6b) where T is the sample-period 5 minutes. This particular state space description (there exists infinitely many different state space descriptions equivalent to (2.6)) is sometimes referred to as controllable canonical form (Glad and Ljung [2000]). 5 It is possible to construct a state space description of a given transfer function with not only the poles as the eigenvalues of A, but also other eigenvalues in addition, a so called non-minimal realization. Once again, the reader are referred to the standard literature, e.g., Rugh [1996], for more information on these topics. 6 This model is a linearization around a specific operating point, and U (z) and Y (z) are not absolute values, but deviations from this point. 2.2 Model Predictive Control 2.1.3 9 State estimation For some applications it is not possible to measure all states of the system, but only, e.g., some of them is included in the output. There are, however, situations when knowledge of the states is needed (e.g., in mpc, see Section 2.2). There are however, in absence of access to the current states, numerous ways to estimate the current states, given only the output, input and the state space model. One such technique is the Luenberger observer. The structure of the Luenberger ˆ observer is presented in (2.7), where the estimated states are denoted x: ˆ + 1) =Ax(t) ˆ + Bu(t) + K[y(t) − C x(t)] ˆ x(t (2.7) where the constant vector K is obtained by solving a discrete algebraic Riccati equation, see Glad and Ljung [2000]. 2.2 Model Predictive Control This section will give an overview of the main idea of Model Predictive Control (mpc). mpc is a technique for designing automatic control algorithms, based on future predictions from a model of the controlled system. This section will only give a brief overview of mpc. A good introduction (in Swedish) can be found in Enqvist et al. [2010], and a more thorough treatment can be found in, e.g, Maciejowski [2002]. mpc is only one out of several techniques for designing feedback controllers. The overall structure of feedback control is illustrated in Figure 2.2. In particular is mpc a state feedback controller, as illustrated in Figure 2.3, where the controller in Figure 2.2 is split into a state estimator (see Subsection 2.1.3) and a controller, e.g., mpc. 10 2 Control theory Figure 2.2: The basic structure of a feedback system: Information from the output signals are used for shaping the input signals. Figure 2.3: The structure for a state feedback controller as, e.g., mpc. The controller in Figure 2.2 is here split into two parts; the state estimator and the controller. 2.2.1 The standard MPC formulation The main idea is to formulate the control problem as an optimization problem over nc future samples. The optimization problem contains an objective function, a prediction model and possibly constraints on, e.g., the input, the states or the output. A rather general formulation of the optimization problem (with, for simplicity, only input constraints) is given in (2.8): min u(t),u(t+1),...,u(t+nc −1) F (x(t), u(t), . . . , u(t + nc − 1)) (2.8) subject to umin ≤ u(i) ≤ umax , i = t, t + 1, . . . , t + nc − 1 where the objective function F contains the system model (2.2) for predictions as a function of the optimization variables u(t), u(t + 1), . . . , u(t + nc ) and the current state x(t). The optimization (2.8) is performed, and u(t) is in each sample executed as control signal. An example for nc = 2 and a linear system model is given in Example 2.3. 2.3 Example: Model Predictive Control Assume the system subject to control is described by a state space model of the form (2.1), and the current state is known to the controller. The ultimate objective is to drive the system states to the origin, using a good trade-off between the ‘energy’ in the input signal (i.e., the square of the inputs) and the sum of the squares of the state deviations from the origin. 2.2 11 Model Predictive Control The prediction horizon is chosen to 2. Since x(t + 1) = Ax(t) + Bu(t) (2.9a) it can be noted that x(t + 2) = Ax(t + 1) + Bu(t + 1) = A(Ax(t) + Bu(t)) + Bu(t + 1) (2.9b) = A2 x(t) + ABu(t)) + Bu(t + 1) (2.9c) Hence the following optimization problem can be formulated X X min x T (i)Qx(i) + u T (i)Ru(i) = u(t),u(t+1) i=t+1,t+2 T (2.10a) i=t,t+1 min u (t)Ru(t) + u T (t + 1)Ru(t + 1) + (Ax(t) + Bu(t))T Q(Ax(t) + Bu(t))+ u(t),u(t+1) (A2 x(t) + ABu(t) + Bu(t + 1))T Q(A2 x(t) + ABu(t) + Bu(t + 1)) = (2.10b) = min  u(t),u(t+1) u(t) 0 0  x T (t) 2 ! ! BT (AB)T B 0 u(t + 1) + Q AB B 0 BT ! ! !T ! A u(t) BT + x(t) Q u(t + 1) A2 0 !  I T x(t) Q I A A  0 R R 0 !! u(t) u(t + 1) (2.10c) where Q and R are the relative weights between the cost of deviations from 0 in the input and deviations from the origin in the states. The choice of Q and R is a tuning of the ‘aggressiveness’ of the controller. Noting that the last term is independent of the optimization variable (and hence not affecting the minimizing argument for the problem) and defining a compact matrix notation, the problem can be formulated as 1 T T U (Bmpc QBmpc + R)U + (BTmpc QAmpc x(t))T U u(t),u(t+1) 2 min (2.10d) where U= u(t) u(t + 1) ! , Ampc = A A2 ! , Bmpc = B AB 0 B ! . The problem in (2.10d) is a convex quadratic programming problem and can be subject to various constraints, such as bounds on the input or the states. These problems are a well studied class of problems, see, e.g., Boyd and Vandenberghe [2004]. ! + 12 2 Control theory The optimizing argument of (2.10d) will give the optimal solution for the stated control problem, fulfilling possible constraints. The solution u(t) is then used as the input to the system (u(t + 1), . . . are simply ignored). The entire optimization problem is thereafter solved all over again in the next sample. 2.2.2 Zone MPC If the mpc control objective is not defined as a specific value for the states, but as an interval for the output, the controller is called a zone mpc (zmpc). Instead of penalizing deviations from a set-point, e.g., x(t) = 0 as in Example 2.3, only deviations outside a zone [ymin , ymax ] are penalized. This optimization problem can be formulated in a similar way, using so called slack variables, see, e.g., Boyd and Vandenberghe [2004]. Another possible extension is to allow the penalization of deviations from the zone or set-point to be asymmetric, e.g., to relate a higher cost to deviations in a positive than a negative direction. The costs associated to deviations from a set-point and a zone respectively are exemplified in Figure 2.4. Figure 2.4: Quadratic output costs for a set-point controller (left), and a zone controller with asymmetric costs (right). The horizontal axis illustrates deviations from the set-point and the zone respectively. 2.3 Invariant sets This section will give a brief overview of invariant sets for discrete time systems. For all cases t is assumed to be an integer, and since the focus will be on positively invariant sets of different kinds, t will in most cases be restricted to non-negative integers. The definitions 2.4, 2.7, 2.8 and 2.9 are intended to be equivalent to, but somewhat less technically stated than the definitions in Blanchini and Miani [2007] and Blanchini and Ukovich [1993] respectively. A more detailed introduction to invariant sets is given in the book Blanchini and Miani [2007] (covering also the continuous time case) and the articles Blanchini [1999] and Dórea and Hennet [1999]. 2.3 Invariant sets 2.3.1 13 Positively invariant sets For an autonomous system (i.e., has no inputs) described by x(t + 1) = f (x(t)), where x ∈ X ⊂ Rn , a positively invariant set is defined as follows: 2.4 Definition (Positively invariant set). A set S ⊆ X is positively invariant (w.r.t. the given system) if x(0) ∈ S implies that x(t) ∈ S for all7 t = 1, 2, . . . This definition of a positively invariant set is illustrated with the following example: 2.5 Example: Positively invariant set for a pendulum Consider a pendulum, with some damping and with small deviations from its downward position. The dynamics of a pendulum can be described with equations on state space form, where the two states can be chosen as x1 position and x2 velocity, as sketched in Figure 2.5. Figure 2.5: A pendulum. The dynamics in the pendulum can be caught in the two states x1 , the position, and x2 , the velocity, as illustrated. Starting at an initial position to the right of the downward position and with an initial velocity to the right, the pendulum will eventually come to a point where the velocity is zero and its distance to the downward position is at its maximum. Thereafter, the velocity will start to increase and be directed to the left, and the distance to the downward position will decrease. Eventually the pendulum will pass its downward position, with a maximum speed toward the left, and thereafter reach the left point with maximum distance, and continue back and forth 7 For the definition of an invariant (not only positively invariant) set, the implication has to hold for any integer t, not only t > 0 14 2 Control theory in the same way. Since there is a damping (e.g., friction) in the system, the maximum distance to the downward position, as well as the maximum speed, will decrease for each iteration. This is plotted in Figure 2.6, using x1 and x2 as the axes8 . Figure 2.6: The time evolution of the pendulum states for a certain set of initial conditions. The evolution is indicated by the arrows. The horizontal axis is the position of the pendulum (i.e., its distance to the downward position, taken with a negative sign if the pendulum is to the left of the downward position), and the vertical axis is the velocity of the weight, with the corresponding sign convention as for the distances. Figure 2.6 reveals that the pendulum’s state, once entering a circle9 of a certain radius, it will never leave the circle in the future. Hence, a circle around the origin is an invariant set for this system. An example of such a circle is sketched in Figure 2.7. 8 Such a plot is usually referred to as a phase plane. 9 In fact, an ellipsoid would be more correct. In these plots, however, the axes are scaled in a way such that the ellipsoid becomes a circle. 2.3 Invariant sets 15 Figure 2.7: A positively invariant set (the circle defined via the black dashed line) for a pendulum. The grey lines (all depicting time evolution of the pendulum states for different initial conditions) all enter the set, but never leave it. 2.3.2 Controlled invariant sets For automatic control purposes, systems without inputs are by obvious reasons of limited interest. Therefore, systems with control inputs are considered, and controlled invariant10 sets are defined for systems on the form x(t + 1) = f (x(t), u(t)), where x ∈ X ⊂ Rn and u ∈ U ⊂ Rm : 2.6 Definition (Controlled invariant set). A set S ⊆ X is a controlled invariant set if there for all x(0) ∈ S exists u(t) ∈ U such that x(t) ∈ S, ∀t = 1, 2, . . . In words, this definition says that there exists an (admissible) input signal u(t) such that if applied to the system, once the states x(t) are within the set, they will remain within the set for all future. The algorithm for computations of the (maximum) controlled invariant set is an inversed reachability analysis, and is presented in Appendix A, and an implementation of the algorithm in Matlab is given in Section B.2 in Appendix B. However, in computations of an invariant set, numerical problems may occur and give incorrect results, e.g., indicating points belonging to the set, although they 10 The term ‘positively’ is here implicit. 16 2 Control theory do not. One remedy to avoid such numerical problems is to introduce contractive sets. For some cases, this also improves the convergence properties for the computations of the sets (Blanchini [1991]). A contractive controlled invariant set for the system x(t + 1) = f (x(t), u(t)), where x ∈ X ⊂ Rn and u ∈ U ⊂ Rm , is defined as follows: 2.7 Definition (Controlled contractive set). A set S ⊆ X is a contractive controlled invariant set (w.r.t. the given system and a given λ, 0 < λ < 1) if there exists a u(t) ∈ U such that x(t) ∈ S implies that x(t + 1) ∈ λS, ∀t = 0, 1, . . . . The core idea is still the same as for invariant sets, but instead of requiring a control signal such that the states remain within the set, the existence of a control signal such that the states belong to a slightly shrunken version of the set in the target step. Another way to interpret it is as a margin is added in each time step. 2.3.3 Robust controlled invariant sets So far, full knowledge of the system and all input signals has implicitly been assumed for the definitions to be useful. The next extension of the idea is to introduce robust controlled invariant sets. Consider a system with two inputs, of which one, u(t), is controlled as before, and the other d(t) is an unknown disturbance, known to belong to a bounded set D: x(t + 1) = f (x(t), u(t), d(t)), where x ∈ X ⊂ Rn , u ∈ U ⊂ Rm and d ∈ D ⊂ Rq , (or, for the linear case x(t + 1) = Ax(t) + Bu(t) + Ed(t)). The definition of the robust controlled invariant set for this system reads: 2.8 Definition (Robust controlled invariant sets). A set S ⊆ X is a robust controlled invariant set (w.r.t. the given system) if there for every x(t) ∈ S exists u(t) ∈ U such that x(t + 1) ∈ S, for all d(t) ∈ D and all t = 0, 1, . . . . One way to interpret the robust set is as a game, where you do not know what your opponent (the disturbance) is doing. You want to define the region where you can ‘feel safe’ and can counteract your opponents turns even if he is doing the (from your point of view) worst possible turns. The extension of robust controlled invariant sets to robust contractive controlled invariant sets should be straight forward. 2.3.4 Periodic controlled invariant sets So far, only different kinds of invariant sets for time invariant systems have been considered. The last theoretical extension of the invariant set idea is an extension to cover periodic systems. A periodic controlled invariant set for a system with periodicity T described as x(t + 1) = f (t, x(t), u(t)), where x(t) ∈ Xt modT ⊂ Rnt modT , u(t) ∈ Ut modT ⊂ Rmt modT and f is periodic such that f (t, x(t), u(t)) = f (t + T , x(t), u(t)), is defined as follows: 2.4 Mini example of invariant set in MPC control 17 2.9 Definition (Periodic controlled invariant sets). A finite sequence of sets Sk ⊆ Xk , k ∈ 1, 2, . . . , T is called a periodic controlled invariant (w.r.t. the given system) if for all x(t) ∈ St modT there exists an u(t) ∈ Ut modT such that x(t + 1) ∈ St+1 mod T , ∀t = 0, 1, . . . A more thorough treatment of periodic sets can be found in Blanchini and Ukovich [1993]. Gondhalekar and Jones [2011] and Gondhalekar et al. [2013] both deals with periodic robust invariant sets, and the former also offers an extension to periodic time-varying state space dimensions, used in an example for asynchronous control. 2.3.5 Maximum invariant sets In many cases can several different sets, all fulfilling the definition of an invariant set, be found. In, e.g., the Pendulum example 2.5 are all circles around the origin positively invariant sets. Throughout this thesis, however, will only the largest possible set (i.e., the maximum set, the union of all possible sets) be considered. 2.4 Mini example of invariant set in MPC control To illustrate how to use invariant sets in mpc, an example is given in this section. The entire Matlab code for the example is provided in Appendix B, Section B.2. 2.10 Example: Robust periodic invariant set in mpc control Consider the dynamical system ! ! ! 1 0 1 1.67 −0.71 d(t) u(t) + x(t) + x(t + 1) = 0 1 0 1 0   y(t) = 17.3 −199 . (2.11a) (2.11b) The input u(t) is bounded by [−0.5, 0.5], and the bounds on the output y(t) are periodic with a periodicity of 24 samples. The disturbances d(t) are unknown, but known to be bounded with periodic bounds. The states are known to the controller, and there are no mismatch between the model and the real system. The periodic output bounds are described in Figure 2.8. The disturbance d(t) is a two dimensional signal, only known to be bounded by Di , which is the unit square with a side of 0.02, except at sample i = 7, 12, 20 and i = 16, 22, where it is bounded by a square in the positive quadrant with one corner in the origin and a side of 0.25 and 0.1 respectively (see the green sets in Figure 2.9). The robust periodic contractive set Xi , i = 1, . . . , 24 with11 λ = 0.99 is computed for this problem. 30 iterations are required, and the set is shown in Figure 2.9. The control of the system (2.11) is formulated as an mpc problem. The mpc problem is formulated with prediction horizon nc = 1 and a quadratic input cost 11 The contraction factor, see Definition 2.7. 18 2 Control theory Figure 2.8: Output bounds for system (2.11) in Example 2.10 R for input signals deviating from 0. The input constraints are used as constraints to the optimization. The invariant set Xi , i = 1, . . . , 24 is incorporated into the mpc problem as follows, in order to maintain the output boundaries: Xi is, for each i, shrunk by the dis˜ i+1 = {x ∈ R2 |x + d ∈ Xi+1 ∀d ∈ Di }. turbance bound for the previous sample12 , X ˜ i are thereafter used as a constraint for the corresponding predicted The set X states. Altogether, the optimization problem solved in each step looks like this: arg min u 2 (t) u(t) (2.12) subject to umin < u(t) < umax ˜ t+1 modT x(t + 1) ∈ X This can be compared to the standard mpc formulation (2.8), with nc = 1 and F (x(t), u(t)) = u 2 (t), and the constraint x(t + 1) ∈ X is the extension from standard mpc to mpc with invariant set. Using this mpc control strategy, with (robust periodic contractive) invariant set, the results shown in Figure 2.10 are obtained for a simulation. Note that no output bounds are violated, and the mpc problem solved in each step uses a prediction horizon of only 1. Also note that there are possibilities for various tuning of the controller; the invariant set constraints guarantee the output bounds, but do not imply any optimality of the control. The objective function (2.12) can still be subject to tuning with, e.g., respect to input and output values. As a comparison, the results of an mpc controller not using the invariant set but only the output boundaries as control objective are shown in Figure 2.11. This is the standard mpc as described in Subsection 2.2.1, but with a zone as control objective, as described in Subsection 2.2.2. This controller has no knowledge of 12 Also named the pontryagin difference, X ˜ i+1 = Xi+1 Di 2.4 Mini example of invariant set in MPC control 19 Figure 2.9: Robust periodic contractive set for the system in Example 2.10. Some characteristics in the set shapes is interesting to note: the sizes of the set at sample 7, 12 and 20 are related to the relatively big disturbance bounds at these samples; to guarantee existence of an admissible control signal controlling the system such that it belongs to the set in the consecutive sample, even with the worst case of the unknown disturbances, the states have to be in this rather small set for these samples. The small size of the set in sample 10 is, in a similar way, related to the tightening of the output bounds in the consecutive sample. the disturbances, and uses the prediction13 horizon 3. Clearly the output trajectory violates the output boundaries several times. The tuning of the controller (i.e., the prediction horizon and the input/output costs) is by no means optimal, but gives still an idea of the advantages of robust mpc with invariant set for this example. The system (2.11) used in this example has some similarities with the insulinglucose model briefly presented in Example 2.1; their input-output behavior have similarities, if the sampling time of (2.11) is interpreted as 60 minutes. Also, some of the input, output and disturbances have similarities with the artificial pancreas control problem (e.g., typical meal times and the samples when disturbances occur), although they are not identical. However, as will be apparent in the sequel, the assumptions of perfect knowledge of the system (i.e., the model is correct) and the current states (i.e., no measurement errors and an ‘ideal’ state estimation) are far away from the reality in the artificial pancreas problem. The 13 Since the system in fact is a non-minimum phase system, see, e.g., Glad and Ljung [2000], a shorter prediction horizon will make the system unstable. 20 2 Control theory Figure 2.10: Result for a simulation of the system in Example 2.10, controlled by an mpc controller using the robust periodic invariant set as constraint and a prediction horizon of 1. (Note that the height of the disturbance bars has no direct interpretation, since the disturbances are two dimensional, but gives only an idea of the relative sizes of the disturbances.) entire Matlab code for this example is found in Section B.2. A detailed treatment of robust periodic invariant sets for mpc and an application example to building climate control can be found in Gondhalekar et al. [2013] 2.4 Mini example of invariant set in MPC control 21 Figure 2.11: Result for a simulation of the system Example 2.10 using a standard mpc controller without knowledge of the disturbances. 3 Diabetes and the artificial pancreas This chapter serves two purposes. First, to give an introduction to Type 1 Diabetes Mellitus (t1dm), and second, to give an overview over some relevant work on the artificial pancreas problem. For a more detailed treatment of diabetes and t1dm is referred to the literature, e.g., Skyler [2012]. For a popular science introduction to diabetes in Swedish, e.g., Hellerström [2002] and Ajanki [1999] is recommended. 3.1 Diabetes The blood glucose level in the human body is affected by various organs and processes in the human body. The blood glucose is an energy source for the cells, and the uptake of blood glucose by the metabolic system is partly controlled by the insulin. Insulin is produced in the beta-cells in the pancreas, and the release of insulin from the beta-cells is controlled by an intrinsic system in the body, together with, e.g., the liver. Type 1 Diabetes Mellitus (t1dm) is an autoimmune chronic metabolic disease characterized by the loss of endogenous secretion of insulin from the pancreatic beta-cells. A lack of insulin causes high blood glucose levels, known as hyperglycemia 1 . This can be dangerous and cause severe and life-threatening complications to the patient if maintained for a long period of time. A low blood glucose value is known as hypoglycemia 2 . Hypoglycemia may cause severe and life-threatening complications even if maintained only for a short period of time. 1 Hyperglycemia is typically defined as a blood glucose level over approximately 180 mg/dL 2 Hypoglycemia is typically defined as a blood glucose level below approximately 70 mg/dL. 23 24 3 Diabetes and the artificial pancreas The unit used for blood glucose level in this thesis is milligram per deciliter, mg/dL. The unit mmol/l is also a common unit for blood glucose values, and 1 mmol/l is approximately 18 mg/dL. 3.2 Insulin treatment Insulin is a peptide hormone, and for a patient suffering from t1dm, there is a need for exogenous insulin. The infusion of insulin has to be done carefully, since an over-delivery of insulin can cause hypoglycemia. Since the discovery of insulin in the 1920s, the traditional way to deliver insulin has been through manual injections based on manual blood glucose readings through finger sticks. Continuous subcutaneous insulin infusion (csii) pumps have been developed and have been available for patient use since the 1980s (Skyler [2012]). csii is, however, only a tool for facilitated delivery of exogenous insulin, and does not necessarily imply improved blood glucose level control. A typical use of csii today is in open-loop control with, e.g., pre-programmed insulin delivery profiles. From a control perspective, the effect of the insulin makes the problem asymmetric, since insulin basically drives the blood glucose value down. However, a meal or a snack typically raises the blood glucose value, which in some sense can be used as a control input, see Subsection 3.6.2. One remedy to this asymmetric problem is to use even another peptide hormone, glucagon, with the opposite effect to the insulin, i.e., driving the blood glucose value up. Such bi-hormonal artificial pancreas systems have been designed, see for example Russell et al. [2012]. This solution is however not used by the ucsb/Sansum group, and will hence neither be used in this thesis work. 3.3 Continuous Glucose Monitoring Recently, subcutaneous continuous glucose monitoring (cgm) techniques have been developed, and have during the last decade become available for patient use (Skyler [2012]). Such devices can give almost continuously (e.g., every 5 minutes) measurements of the subcutaneous glucose. The relevant quantity to control for a t1dm patient is however not the cgm readings, but the intravenous blood glucose level. The cgm measurements do not always coincide with the intravenous blood glucose level, for physiological reasons (the blood glucose and the subcutaneous glucose values can be different (Steil et al. [2003])) as well as technical reasons (errors and uncertainties in the measurement method and equipment). From a control perspective the cgm measurements are hence affected by a noise, which is not neglectable (Zisser et al. [2009]; Huyett et al. [2013]). Currently there does not exist any, known to the author, control-relevant description of this noise. The information itself, obtained from the cgm, does of course not improve the control of the blood glucose level. However, in a system with both cgm and csii, 3.4 Metabolic system simulator 25 the blood glucose level control can possibly be improved, as will be discussed in Section 3.6. 3.4 Metabolic system simulator A simulation environment of the metabolic system for evaluation of artificial pancreas has been developed at University of Virginia and University of Padova (Kovatchev et al. [2009]). The simulator is approved by the U.S. Food and Drug Administration (fda) for verification of control algorithms for artificial pancreas before doing clinical trials on t1dm subjects. This simulation environment is in the sequel shortly referred to as the UVa/Padova simulator. The simulator contains a set of different virtual patients with varying parameters, such as weight, age, gender etc, see Subsection 4.4.3. There is also a noise model simulating cgm measurement errors for closed loop control simulations. 3.5 System identification for an insulin-glucose model To obtain a model of the insulin-glucose relation suitable for (linear) mpc, system identification has been done in the UVa/Padova simulator (see Subsection 3.4). The obtained model is a (linear, time-invariant) third order transfer function (an arx model) describing the dynamics of the insulin-glucose relation, as deviations from an input of basal3 rate and an output of 110 mg/dL (van Heusden et al. [2012]): Y (z) = −5 F · 1800 1 (1 − p1 )(1 − p2 )2 U (z) tdi (z − p1 )(z − p2 )2 (3.1) where • p1 = 0.98 • p2 = 0.965 • F is a safety constant tuned by a physician, with a typical value of 1.5 • tdi is the total daily insulin, also prescribed by a physician, with a typical value between 25 and 70 • U is the input signal, i.e., injected insulin, in Insulin Units per hour (U/h). • Y is the output signal, i.e., the blood glucose level, in mg/dL. The sampling time of the model is 5 minutes. Altogether, a typical value of the numerator could be −0.01, which is used for computations in this thesis when nothing else is mentioned. 3 A value prescribed by a physician, related to the tdi value. 26 3 Diabetes and the artificial pancreas The transfer function is a linearization around the patients basal rate (tuned by a physician, which typically is of the size 1 U/h) and a blood glucose value of 110 mg/dL. There is no possibility to give ‘negative’ insulin. But since the model is described as deviations from a working point, a suspension of the pump (which is the smallest possible input) corresponds to minus the basal rate as the value of u(t), i.e., typically −1 U/h. There is also an upper limit of the pump of the size 100 U/h. Altogether, the input is approximately bounded as u(t) ∈ [−1, 100]. There is, based on experience from in silico trials, a non-neglectable plant-model mismatch between this model and the virtual subjects in the simulator. 3.6 Artificial pancreas Figure 3.1: The structure of an artificial pancreas A lot of work has been carried out by different research groups at different places to design closed loop systems for automatic delivery of insulin, i.e., artificial pancreas. An overview of the basic structure of an artificial pancreas is given in Figure 3.1. An artificial pancreas has been developed in the group at ucsb/Sansum (van Heusden et al. [2012]; Harvey et al. [2010]). Especially the zmpc technique has shown some promising results (Dassau et al. [2013]; Grosman et al. [2010]). In the zmpc control algorithm developed in the ucsb/Sansum group, the model (3.1) presented in Subsection 3.5 has been used. An overview of the ucsb/Sansum group system is given in Figure 3.2. 3.6.1 Meals One of the most important and characteristic disturbances (from a control perspective) are the meals. Different approaches can be used to deal with meals: Announced meals With announced meals, the controller is assumed to have knowledge of the time and the size of the meals the patient is having. The controller can hence ‘counteract’ the meal in advance and give a dose of insulin before or at the time the meal is eaten, a so called meal-bolus 4 . The size and the distribution over time of the 4 The idea with the term is to emphasize the similarities with the insulin secretion from normalfunctioning beta cells, which can be divided into bolus and basal parts, see Skyler [2012]. 3.6 Artificial pancreas 27 Figure 3.2: An overview of the ucsb/Sansum group artificial pancreas system. bolus is subject to current research. Announced meal can be designed either as real time inputs from the patient, or forcing the patient to follow a meal schedule known by the controller on forehand. Depending on context, only the former may be called announced meal, but they are equivalent from a control algorithm perspective, since the information on the meal is available to the controller before its effects are seen in the system output. From a control point of view, this can be considered as a feed-forward control design. Unannounced meals With unannounced meals, the controller has no information of the meals a patient is having. This is the case considered in this thesis, although extensions to the announced case could be possible. 28 3 Diabetes and the artificial pancreas Meal detection When having unannounced meals, ideas from the announced meal strategy can still be used as inspiration for the design of the controller. The idea is to detect a meal, basically by a rising cgm value, and give a bolus to treat the meal, as if it was announced. Work on meal detection can be found in Dassau et al. [2008]. By intuition, it is in general possible to achieve a better control performance having information of the meal, as opposed not to have the information. It is, however, also possible to discuss what the best is from a patient perspective. A comparison between announced and unannounced meal for zmpc control is found in Grosman et al. [2010]. 3.6.2 Safety systems and hypoglycemia treatments Since the control problem is asymmetric not only in the input (see Section 3.2) but also since hypoglycemias are more critical than hyperglycemias, safety systems have been developed to react on low cgm readings and, e.g., to give a warning and request the patient to eat some carbohydrates (Dassau et al. [2010]). 3.6.3 Insulin on board constraints To prevent too much insulin to be delivered during a certain time window, the concept of insulin on board constraints has been developed. The problem with over-delivery of insulin can be relevant with, e.g., a plant-model mismatch, measurement errors and an aggressively tuned controller. It can especially be problematic in combination with announced meals where big boluses are given in addition to the action suggested by the control algorithm. The constraints are in principle an integral constraint on the input signal. A detailed description and evaluation of the insulin on board constraints can be found in Ellingsen et al. [2009]. 3.6.4 Night time zone A severe hypoglycemic event during the night could be even more dangerous than during the day, since the possibilities to action from the patient or another person from the household (to, e.g., give carbohydrates to the patient to eat) are more limited during night than day-time. Hence, the concept of a night-time zone is introduced, with a different control objective and tighter constraints on the control signal. A typical value for the night time control objective zone is 110 − 170 mg/dL (instead of the day time control objective of 80 − 140 mg/dL), as shown in Figure 3.3. The night time control objective has been tuned to reduce the risk of hypoglycemia, based on experiences from the clinical trials in the ucsb/Sanum group. The control signals can be limited to typically never exceed, e.g., 1.6 times the basal rate (i.e., u(t) ≤ −0.6 · umin , since u(t) is described as deviations from the basal rate and the basal rate is −umin ) during night time zone. 3.6 Artificial pancreas 29 Figure 3.3: Example of an overnight zone (00:00 - 05:00) with an heightened control objective zone (110-170 mg/dL), for decreased risk of hypoglycemic events during night time. 4 MPC with invariant set constraints for an artificial pancreas In this chapter, the contribution to artificial pancreas control by this thesis, is presented. Various design choices and results are presented and discussed. 4.1 Invariant sets for the insulin-glucose model The controlled invariant set as defined in Definition 2.6 for the insulin-glucose model (3.1) was computed with the Matlab code presented in Section B.2. The set is shown in Figure 4.1a. (An alternative idea to plot the invariant set for the purpose of better intuitive understanding of its meaning is presented in Appendix C.) The computations of the invariant set converges for the particular numerical values used (tdi, input and output limits, etc) after 57 iterations (which corresponds to approximately 10 minutes on an Intel Core i5 2.5 GHz CPU with the implementation in Appendix B). The number of iterations, which also affects the computation time, is highly sensitive to different numerical values. The contractive controlled invariant set as defined in Definition 2.7 for λ = 0.99 and the insulin-glucose model (3.1) computed with Matlab is shown in Figure 4.1b. The computations converges after 63 iterations (approximately 20 minutes on the same CPU). A periodic controlled invariant set for a 24 hour period (i.e., 288 samples with 5 minutes sampling time) for periodic output constraints, can also be computed in a reasonable time. As noted in Subsection 2.1.1, a state space model is not a unique representation of an input-output relation, but there exists infinitely many state space descriptions 31 32 4 MPC with invariant set constraints for an artificial pancreas (a) Controlled invariant set (b) Contractive controlled invariant set Figure 4.1: Controlled invariant and contractive controlled invariant set for the insulin-glucose model (3.1). The sets are very similar, but the contractive controlled invariant set is a subset of the nominal controlled invariant set, and hence slightly smaller. for one input-output relation. However, all (minimal) state space descriptions of the same input-output relation are related through (invertible) linear transformations (Rugh [1996]). Hence the corresponding invariant sets for the different state space descriptions are related through (invertible) linear transformations as well. The structure of the invariant set (number of vertices, etc) can therefore be expected to be independent of the choice of state space description, although the numerical properties of the computations probably can be affected. The influence of different state space descriptions have not been analyzed thoroughly, and the state space description presented in Example 2.1 and used in the controller structure developed by the ucsb/Sansum group was used throughout this thesis. 4.2 Robust invariant set for the insulin-glucose model In addition to the nominal and periodic controlled invariant sets presented in Section 4.1, different robust invariant sets have been considered. 4.2.1 Measurement errors in the state space Since the invariant set is a feature in the state space, the measurement errors (naturally entering in the ‘output space’) have to be converted to a ‘state space feature’ to be able to include in the invariant set framework. The state estimator (estimating the current states from the current and the previous output and input signals, explained in Subsection 2.1.3) is relevant for this conversion. At the time when writing this thesis, the artificial pancreas of the UCSB/Sansum group used a Luenberger observer as state estimator. Simulated measurement errors 4.2 Robust invariant set for the insulin-glucose model 33 were propagated through the state estimator into the state space using a generic method with Monte Carlo simulations. Since no extensive model for the measurement errors was present (Section 3.3), the errors were assumed to be white noise with a realistic variance. Figure 4.2: The propagated measurement errors, from three different angles. According to the simulations, the errors converge to a two-dimensional plane (in the three-dimensional state space), as shown in Figure 4.2. The existence and properties of this plane could not be supported by any theory found by the author in the literature. The existence of the plane was however utilized for fitting a two-dimensional rectangle in this plane to the simulated errors. The 90% of the simulated errors that minimized the circumference of the rectangle were chosen to define the bounds for the rectangle, as shown in Figure 4.3. This rectangle was considered as the measurement error bound in the robustification method described in Subsection 4.2.2. Because of the unknown properties of the errors and the disturbances, the statistical approach with Monte Carlo simulations was considered to be sufficient for this purpose, although an analytical method would be more rigorous. 4.2.2 Robustification against measurement errors The measurement errors propagate, as described in Subsection 4.2.1, through the state observer to the state space. Since the state estimator is a linear system, the propagated measurement errors can be described as an additive error in the state space. The measurement errors do however not affect the system dynamics behavior, as the disturbances do in the disturbed system description in connection to Definition 2.8. A robustification according to Definition 2.8 is therefore not relevant for this problem. Instead, the robustification method chosen for this problem is to interpret the disturbance bounds as a desired robustness margin against the bounds of the states, and compute the invariant set for these tighter bounds. The only bounds on the states are imposed by the constraints ymax ≥ y = Cx and ymin ≤ y = Cx, which for the case without night-time zones reduces to −30 ≤ x3 ≤ 30. However, since x3 (t + T ) = x2 (t), and x2 (t + T ) = x1 (t) for this state space description (see 34 4 MPC with invariant set constraints for an artificial pancreas Figure 4.3: The propagated measurement errors, with 90% of the simulations in the two-dimensional box. Example 2.2), the bounds [−30, 30] will apply to those states as well. The robustification is hence to subtract the disturbance bound (Figure 4.3) from the box [−30, 30] × [−30, 30] × [−30, 30], and compute the invariant set with these new, tighter bounds on the states1 . 4.2.3 Meals The design of a controller that is robust against meals is an interesting challenge. Without any prior knowledge of the meals, and with the present cgm errors, it appears from experience of in silico trials almost impossible to maintain the blood glucose zone when the patient is having a big meal. The lack of a model describing the relationship between a meal and the blood glucose value relevant for control is a bottleneck for further investigation of this theme. However, with a relevant model, robustification according to Definition 2.8 against a smaller meal could probably be interesting to investigate. 1 A comparable, although not equivalent, approach would be to simply shrink the control objective zone by the 90% bound of the measurement errors. 4.3 Implementation of set constraints in MPC 4.2.4 35 Plant-model mismatch With a description for the plant-model mismatch, an invariant set approach robust against plant-model mismatches could be a very interesting concept to investigate. Due to the absence of such a description, this theme is however not looked into further. 4.3 Implementation of set constraints in MPC The invariant set is used as constraints for the predicted states in the mpc optimization, basically as illustrated in Example 2.10, and in particular in (2.12). This is to be compared with the zmpc used, where only the predicted output (and not the predicted states) were directly considered in the optimization. However, to enforce constraints on the predicted states (as done in (2.12)) is in practice not a good idea, due to the risk of formulating an infeasible2 problem, if, e.g., the system has been subject to a major disturbance, and cannot be brought back to the invariant set in one sample. This is highly relevant for the artificial pancreas problem, since the experience suggests that it is hard not to exceed the output 140 mg/dL after a big meal. Hence, the implementation of the constraints must be able to handle deviations from the invariant set. Therefore, the constraints are implemented as soft constraints, which means that deviations from the set are allowed, but penalized in the optimization objective function. As opposed to this, a normal constraint implementation will be named hard constraint for clarity. The soft constraint will cause no issues with infeasible solutions due to the invariant set, since there are no (hard) constraints related to the invariant set to violate. There are, however, various ways to parametrize the penalization of the soft constraints: Assume the invariant set is described as Ax ≤ b. The implementation of this as a hard constraint into the mpc optimization problem (2.8) would then be min u(t),u(t+1),...,u(t+nc −1) F (x(t), u(t), . . . , u(t + nc − 1)) (4.1) subject to Ax(t + 1) ≤ b, . . . , Ax(t + nc ) ≤ b However, implementing the set Ax ≤ b as a soft constraint can either be made as min u(t),u(t+1),...,u(t+nc −1) F (x(t), u(t), . . . , u(t + nc − 1)) + G(t+1 , . . . , t+nc +1 ) subject to Ax(t + 1) ≤ b + t+1 , . . . , Ax(t + nc ) ≤ b + t+nc t+1 , . . . , t+nc ≥ 0 2 This is of course relevant for all constraints, not only constraints related to invariant sets. (4.2) 36 4 MPC with invariant set constraints for an artificial pancreas Figure 4.4: An asymmetricly located rectangle described by Ax ≤ b, compared to Ax ≤ b +  for an  > 0, and Ax ≤ θb with θ (= 2) > 1. Note the different scalings of the rectangle, and compare to the discussed - (4.2) and θ-implementation (4.3) of invariant set constraints. or as min u(t),u(t+1),...,u(t+nc −1) F (x(t), u(t), . . . , u(t + nc − 1)) + G(θt+1 , . . . , θt+nc +1 ) (4.3) subject to Ax(t + 1) ≤ bθt+1 , . . . , Ax(t + nc ) ≤ bθt+nc θt+1 , . . . , θt+nc ≥ 1 where G( · ) is a strict convex function, e.g., a quadratic function. If  > 0 or θ > 1 (and not equal to 0 or 1 respectively), it indicates deviations from the invariant set in the predicted states. θ and  can either be scalar or a matrix or vector with one scalar variable for each inequality in A. The difference between (4.2) and (4.3) may be subtle, but produces significantly different results for this application and is worth a further discussion. The -implementation of the invariant set, (4.2), gives the deviation from the invariant set in absolute terms. It is, however, dependent on the scaling of the rows of the matrix A. If the rows of A are normalized to the length of 1 (which they are with the implementation used for computations of A),  will be in the units of mg/dL, because of the structure of the state space model. The θ-implementation of the invariant set, (4.3), gives the deviation from the invariant set in relative terms, which makes θ unitless and not dependent on the scaling of A. Studying the invariant set in Figure 4.1a, it is clear that some constraints are much closer to the origin than others, which will highly affect the penalization using the θ-implementation and give different ‘weight’ to different inequalities in A. The qualitative difference between the - and θ-implementation is illustrated in Figure 4.4. 4.4 37 In silico trials 4.4 In silico trials The mpc algorithm with invariant set for artificial pancreas has been evaluated in in silico trials. The trials have been done in version 3 of the UVa/Padova simulator (Kovatchev et al. [2009]), see Section 3.4). The trials have been performed as a benchmark against a zmpc controller form the ucsb/Sansum group. This section contains information on the trials, as well as the results. 4.4.1 Test protocol The trial started at 14.00, and the closed loop control of the subjects started 2 hours thereafter. Meals were given according to Table 4.1 and the trials ended after 52 hours. Meal Time Size Dinner 18.30 50 g Breakfast 7.00 40 g Lunch 13.00 50 g Snack 16.00 10 g Dinner 18.30 50 g Breakfast 7.00 40 g Lunch 13.00 50 g Table 4.1: Meal schedule in the simulations. The sizes of the meals is the carbohydrates, e.g., the lunch contains 50 g carbohydrates. To ensure comparability between the results for different controllers, cgm errors were generated using the same random seed. 4.4.2 Different controller setups in the trials Two different controllers were used in the trials. In addition to them, several other controllers have been tested and evaluated in shorter trials, see Section 4.4.5. 1. zmpc: The benchmark controller, tuned with experience from the ucsb/Sansum group: Asymmetric input costs, control objective 80-140 mg/dL, night time zone (see Subsection 3.6.4) between 00:00 and 05:00 with control objective 110-170 mg/dL and an upper bound of control action limited to 1.6 basal rate in the night time zone. 2. MPC with invariant set: The proposed controller; mpc with nominal invariant set as soft constraints with scalar θ-implementation (see Section 4.3). The set used for all controllers was a set computed for a tdi value of 33 and a basal rate of 1. The input costs were asymmetric and tuned iteratively. The controller also had a night time zone with an upper bound of the control action limited to 1.6 basal rate. For all controllers the safety systems (Section 3.6.2) and insulin on board constraints (Section 3.6.3) were inactive, as an attempt to evaluate the efficacy of the control algorithm itself without including the complexities of any safety layers. No information on the meal was present to the controllers, i.e., unannounced meals. 38 4 # Age Weight (kg) tdi (U) Basal (U/h) 1 32 80 42 1.2 MPC with invariant set constraints for an artificial pancreas 2 22 80 43 1.3 3 42 71 52 1.5 4 24 67 35 1.0 5 47 67 40 0.9 6 23 73 72 1.9 7 47 46 43 1.2 8 56 98 52 1.1 9 24 68 34 0.9 10 31 81 47 1.2 Table 4.2: Some parameters for the virtual subjects in the trial 4.4.3 Subjects in the trial The simulator contains a set of 10 virtual subjects. Each subject is associated with individual different physiological parameters. A few of these parameters are presented in Table 4.2, to get an idea of the variety between the subjects in the trial. 4.4.4 In silico results The results from the in silico trial are presented in this section. The results are presented in a few different ways • Tables with ratios of time spent in different blood glucose intervals • Population plots showing the blood glucose mean, standard deviation and min/max values over all subjects, for each time instance • Control Variable Grid Analysis (cvga) plots, showing the 90% confidence bounds on the highest and lowest blood glucose value for each subject during the trials • Entire blood glucose trajectory and control action for 3 different subjects. Tables The in silico results are presented for Controller 1 (the zmpc controller) in Table 4.3a, and for Controller 2 (the mpc controller with invariant sets) in Table 4.3b. Each row adds up to 100, except for rounding effects, and shows the distribution of all blood glucose values for all subjects during the closed loop control in the trials. Population plots The population plots in Figure 4.5a and 4.5b shows the mean, standard deviation and min/max envelope for the blood glucose values in the trials. 4.4 39 In silico trials All day 80-140 mg/dL All day 70-180 mg/dL 22:30-07:00 80-140 mg/dL 22:30-07:00 70-180 mg/dL 00:00-05:00 110-220 mg/dL below 0.37 0.04 0.63 0.00 4.82 within 41.62 74.81 58.40 98.12 95.18 above 58.00 25.15 40.98 1.88 0.00 (a) Controller 1, the benchmark zmpc controller All day 80-140 mg/dL All day 70-180 mg/dL 22:30-07:00 80-140 mg/dL 22:30-07:00 70-180 mg/dL 00:00-05:00 110-220 mg/dL below 0.25 0.00 0.77 0.00 2.74 within 42.27 74.66 69.38 99.81 97.26 above 57.48 25.34 29.84 0.19 0.00 (b) Controller 2, the proposed invariant set controller. Table 4.3: Time in range-statistics for the in silico trials. All values are percentages. CVGA plots The Control Variable Grid Analysis plot (cvga , Magni et al. [2008]) is a plot developed for being a compact way to illustrate trial results. In short, the cvga plot is showing the 95% confidence bound for the highest and lowest blood glucose value for each subject as y- and x-coordinate respectively, and the space is divided into different zones (represented by different colors) to facilitate evaluation. For instance is a subject in the A-zone intended to be preferred over the B-zone. The cvga plots for the trials are found in Figure 4.6a and 4.6b. Note however the nonlinear scaling of the y-axis. General comments In general, the two different controllers performs similarly, but not identically. It is apparent from the population and cvga plots that the proposed Controller 2 results in a smaller variation between the subjects, i.e., the min-max envelope in the population plot is smaller and the points in the cvga are more clustered. It is also interesting to note the more steady and less ‘oscillatory’ control, especially during the night, for Controller 2 (the proposed controller) compared to Controller 1 (the zmpc controller). Subjects controlled by Controller 2 starts their new day, in general, with a less oscillating blood glucose curve, compared to Controller 1. In this comparison, it is however important to remember the different tuning with an overnight zone with a higher objective zone (110-170 mg/dL) for the blood glucose values in Controller 1 than for Controller 2. 40 4 MPC with invariant set constraints for an artificial pancreas (a) Controller 1, the benchmark zmpc controller (b) Controller 2, the proposed invariant set controller Figure 4.5: Population plots. The solid blue line is the mean, the red dashdotted line is 1 standard deviation, and the black dash-dotted line is the min-max envelope, of all subjects. 4.4 41 In silico trials (a) Controller 1, the zmpc controller (b) Controller 2, the invariant set controller Figure 4.6: cvga plots. Each subject is represented by a small black dot, and the mean of all points is represented by a bigger blue dot. 42 4 MPC with invariant set constraints for an artificial pancreas Some interesting subjects Out of the 10 virtual subjects in the simulator, the results for subject 1, 7 and 9 have been chosen to be presented in detail, since they appear to represent characteristics typical for the trial. The blood glucose trajectory for the two different controllers as well as the cgm reading and the controller action for subject 1 are shown in Figure 4.7. Note the more steady overnight control, the somewhat higher top value of the peaks, slightly more damped oscillations after meals and the more aggressive control action for Controller 2 (mpc with invariant set) than for Controller 1 (zmpc). Subject 7, Figure 4.8, is the subject in the trials where the difference between the controllers is largest. The most striking difference is probably the less ‘oscillatory’ behavior in the blood glucose curves for Controller 2, and the higher lows. This behavior seems to be achieved thanks to the more aggressive control action in Controller 2 (than 1) which gives more insulin during the rising slopes and suspends the pump earlier on downward slopes. The results of subject 9, Figure 4.9, are similar to subject 1. Here, however, all peaks are of the same height, but the differences in the lows are larger between the controllers. The blood glucose trajectory for Controller 2 is hardly never less than 110 mg/dL, although Controller 1 goes below 100 mg/dL several times. 4.4.5 Evaluation of alternative controllers Other controllers have also been tested in silico, but are only briefly presented here, since their performance has not been found as promising as Controller 2 in the in silico trials. • mpc with an invariant set for a tdi value close to the actual tdi value3 of the subject. Only very small differences in performance were found, compared to Controller 2, not giving a qualitatively different behavior. However, the tdi value is inverse proportional to the model gain (see (3.1)). Since the effect of the input is scaled by the model gain, and the lower bound on the input is the basal rate, the relevant value for the invariant set shape is in fact not the tdi value nor the basal rate, but the fraction Basal/tdi (assuming that the exact value of the upper bound on the input is of less importance). Looking at Table 4.2, it becomes clear that the range of values for the fraction Basal/tdi is smaller than the range of tdi and basal rates respectively, i.e., there is a correlation between basal rates and tdi values. Hence it is justified to use only one invariant set for all controllers, as done in the in silico trials. • mpc controller using an invariant set robust against cgm errors according to Subsection 4.2.3. However, a robustification against realistic measurement errors yields a controller with a very small invariant set (cf. shrinking the objective zone with the upper bound on measurement errors). With 3 In the in silico trials, a set for a tdi value of 33 was used for the control of all subjects. 4.4 43 In silico trials Blood glucose level 300 Glucose [mg/dL] 250 200 150 100 50 00:00 70−180 mg/dL 80−140 mg/dL 00:00 Meal BG (Invariant Set) CGM (Invariant Set) BG (zMPC) CGM (zMPC) Invariant Set control action. TDI (18:00−18:00) = 31.95 [U] , Basal = 1.20 [U/h] 0.5 Insulin [U] 0.4 0.3 0.2 0.1 0 00:00 00:00 Time of day zMPC control action. TDI (18:00−18:00) = 32.75 [U] , Basal = 1.20 [U/h] 0.5 Insulin [U] 0.4 0.3 0.2 0.1 0 00:00 00:00 Time of day Figure 4.7: Blood glucose trajectory for subject 1. The blood glucose curve for controller 2 (the mpc controller with invariant set) is the solid blue line, and the controller action is the upper graph. The blood glucose for Controller 1 (the zmpc controller) is the dashed black line, and the lower graph with control action. 44 4 MPC with invariant set constraints for an artificial pancreas Blood glucose level 300 Glucose [mg/dL] 250 200 150 100 50 18:00 70−180 mg/dL 00:00 80−140 mg/dL 06:00 12:00 18:00 Meal BG (Invariant Set) 00:00 CGM (Invariant Set) 06:00 12:00 BG (zMPC) 18:00 CGM (zMPC) Invariant Set control action. TDI (18:00−18:00) = 35.60 [U] , Basal = 1.25 [U/h] 0.5 Insulin [U] 0.4 0.3 0.2 0.1 0 18:00 00:00 06:00 12:00 18:00 Time of day 00:00 06:00 12:00 18:00 12:00 18:00 zMPC control action. TDI (18:00−18:00) = 36.00 [U] , Basal = 1.25 [U/h] 0.5 Insulin [U] 0.4 0.3 0.2 0.1 0 18:00 00:00 06:00 12:00 18:00 Time of day 00:00 06:00 Figure 4.8: Blood glucose trajectory for subject 7 4.4 45 In silico trials Blood glucose level 300 Glucose [mg/dL] 250 200 150 100 50 18:00 70−180 mg/dL 00:00 80−140 mg/dL 06:00 12:00 18:00 Meal BG (Invariant Set) 00:00 CGM (Invariant Set) 06:00 12:00 BG (zMPC) 18:00 CGM (zMPC) Invariant Set control action. TDI (18:00−18:00) = 26.30 [U] , Basal = 0.95 [U/h] 0.5 Insulin [U] 0.4 0.3 0.2 0.1 0 18:00 00:00 06:00 12:00 18:00 Time of day 00:00 06:00 12:00 18:00 12:00 18:00 zMPC control action. TDI (18:00−18:00) = 26.80 [U] , Basal = 0.95 [U/h] 0.5 Insulin [U] 0.4 0.3 0.2 0.1 0 18:00 00:00 06:00 12:00 18:00 Time of day 00:00 06:00 Figure 4.9: Blood glucose trajectory for subject 9 46 4 MPC with invariant set constraints for an artificial pancreas a very small target set, the controller becomes similar to a set-point controller, and the advantages of the zmpc idea (Grosman et al. [2010]) becomes less present, and the controller performs poorer than any of the controller in the in silico trials. • mpc controller with invariant set robust against meals or plant-model mismatch. In the absence of a meal model or an uncertainty bound on the model (3.1), no meaningful way to robustify against meals or plantmodel mismatch was found, as described in Subsection 4.2.3 and 4.2.4 respectively. • mpc controller with invariant set using -implementation (see Section 4.3) for penalization of deviations from the invariant set. The result for such a controller are more similar to Controller 1, the zmpc controller, than Controller 2. This could however be an interesting controller if Controller 2 is found to be, e.g., too aggressive. • mpc controller using invariant set with one θ/ variable per constraint (and not one per prediction step, see Section 4.3). This increases the size of the optimization problem drastically, and was not realistic to solve in the framework the implementation was done, using the built-in optimization solvers in Matlab. • mpc controller using contractive invariant set. Due to the very big similarities between the controlled invariant and the controlled contractive (λ = 0.99) invariant sets (see Figure 4.1), is the perfomance of a mpc controller using a contractive set almost identical to a controller using the invariant set controller. • mpc controller with periodic invariant set using overnight zones with a higher control objective (e.g., 110-170 mg/dL) during night time. Such a controller is possible to implement, but since Controller 2 did not show any problems at all with hypoglycemia during overnight, this controller was decided not to be a part of the in silico trials. 5 Concluding comments This chapter gives some comments on the presented work and results. There are also some suggestions on directions for further work. 5.1 Comments and conclusions First, the in silico results will be commented in this section, and the approach with invariant set in general will be discussed thereafter. 5.1.1 In silico results The in silico results looks overall promising and interesting, as seen in Figure 4.5. However, the most interesting evaluation would be a clinical trial, and too comprehensive conclusions should not be drawn from the in silico results. 5.1.2 Invariant set approach Since the same model has been used for prediction in the invariant set mpc controller as in the zmpc controller, and somewhat ‘better’ results have been obtained, it is clear that there exists information in the predicted states that is not fully utilized with the zmpc approach. The concept of invariance is obviously one way to exploit this information, however not proven to be the optimal one. It is not really clear whether the improved results actually are due to the concept of invariance, or if it is of a more coincidental occurrence. There could possibly exists another underlying structure or concept that is more powerful to utilize for the purpose of improved mpc for artificial pancreas control, than the invariant set. 47 48 5 Concluding comments The invariant set approach is relatively ‘model-intensive’, and the results would probably improve with even better models. Even though the arx model (3.1) used for the insulin-glucose relation indeed is relevant to the problem, it does not catch the entire characteristics of the dynamics. Also the shortage of satisfying models for, e.g., meals, cgm errors and plant-model mismatch is apparent. Although the results for the invariant set control are promising, the tuning and the understanding of the control algorithm is more advanced than for the zmpc solution. E.g., the computation of a new invariant set might be necessary in a retuning of an mpc with invariant set, which requires a higher workload than only changing some parameters during a re-tuning of the zmpc. 5.2 Further work In this section, some interesting ways to extend and further investigate the presented work are suggested. • The tuning of Controller 2 in the in silico trials is by no means done systematically, and a systematic tuning of the parameters (i.e., control objective zones, input/output costs, - or θ- implementation, prediction and control horizon, different kinds of invariant sets (nominal, contractive, robust, periodic), etc.) could probably result in an even better control. • The absence of a model describing meals has been a bottleneck for all attempts to develop a controller robust to meals. If, e.g., the insulin-glucose model (3.1) was extended with a second input (describing meals), a robustification against a possible meal at, e.g., the time 13.00 and 18.00 could be an interesting concept to investigate. Such an extension of the model would, however, in the best case not affect the dimensions of the model (3.1), since the computations of invariant sets with the presented methods would probably be a new bottleneck for a 4- (or even higher) dimensional model. • Given knowledge on the mismatch between the model (3.1) and the true system, a controller using invariant sets robust against such plant-model mismatches could be an interesting subject. A parametrization of the mismatch could possibly be a way to better personalization of the control to each subject, as intended by the safety factor F in (3.1) (van Heusden et al. [2012]). • The attempt to develop a controller robust to cgm errors has so far not been successful, see Section 4.4.5. However, with better models of the error dynamics and with more sophisticated theory utilizing not only the error bounds, but also the error dynamics, such an attempt could be interesting. 5.2 Further work 49 • To determine whether the obtained results are due to the concept of invariance or not, as discussed in Section 5.1, it could be interesting to develop an mpc strategy (without invariant set) penalizing not only the absolute values of the cgm readings, but also the rate of change, and see if similar results are possible to obtain with such an approach. A Invariant set mathematics In this chapter, a description of the mathematical algorithm for computations of invariant sets using convex polyhedral sets, is provided. For computations of the robust and periodic case is referred to Blanchini [1999] respectively Blanchini and Ukovich [1993]. A.1 Controlled invariant set A controlled invariant set (Definition 2.6) for the system x(t + 1) = Fx(t) + Gu(t), y = H x(t), x ∈ Rn and u ∈ U (where U is a convex polyhedral set) can be computed via a reversed reachability analysis algorithm. The algorithm is described in detail in the following: 1. Define1 A1 and b1 such that the set X1 = {x ∈ Rn |A1 x ≤ b1 } describes the on given state constraints2 . E.g., output constraints " # the"form y ∈#[ymin , ymax ] H ymax can be transformed to state constraints as x≤ . −H −ymin ¯ i as the set of input signals u ∈ U and states 2. Define the lifted state space X x ∈ X such that" the states in the#following " # " time# step will be in Xi−1 , i.e., x bi ¯ i = {x ∈ Rn+m | Ai−1 F Ai−1 G X ≤ }. 0 Au u bu 1 Although the convention of using the index 0 for this appears to be the most common in the literature, index 1 is here used to conform with the Matlab code. 2 Which here is implicitly assumed to be a convex polyhedral set. 51 52 A Invariant set mathematics 3. Since the invariant set concept only cares about existence of control signals, ¯ is proand not the control law itself, the ‘input dimensions’ of the set X ¯ n ¯ i , n) jected onto the subspace R where the states lives, i.e., Xi = Projection(X ¯ and X , i.e., X = 4. Xi is finally between X i i−1 i " obtained # "by intersection # ¯ ¯ Ai bi ¯ . n ¯ ¯ {x ∈ R | x≤ }, where Ai and bi is describing the set X i Ai−1 bi−1 5. If Xi is the empty set, the invariant set does not exist, and the procedure can be aborted. 6. If Xi = Xi−1 , a fixed-point is reached and Xi is the invariant set. 7. Otherwise, return to Step 2. B Implementation This appendix contains some details on the implementation for various tasks in this thesis. B.1 Software The software used for the computations of the invariant sets are listed in Table B.1. Name Matlab mpt3 Gurobi Organization/ Authors The MathWorks, Inc. Kvasnica et al. [2004] Gurobi Optimization, Inc. Version R2011b 3.0 5.1.0 URL www.mathworks.com/matlab www.tbxmanager.com www.gurobi.com Table B.1: Software used for computations Version 3 of the UVa/Padova simulator (Magni et al. [2008]) was used for the in silico trials (see Section 4.4) . The controller code used for the trials was based on the code described in Grosman et al. [2010]. B.2 Matlab code The Matlab code for some computational tasks are provided in this section. First the code for computing the invariant set for the insulin-glucose model (3.1) on the state space form (2.6) is found. Second the code for Example 2.10 is presented. 53 54 B B.2.1 Implementation Code for invariant set computations The Matlab code in this subsection is divided into two parts. The first part defines the insulin-glucose model (3.1) on state space form (2.6), and calls the second part, the function InvariantSetDiscreteLTISystem.m. This code requires the toolbox mpt3, see Section B.1. The comments about step numbers in the second part refer to the mathematical algorithm presented in Section A.1. main.m K = -0.01; a1 = -0.965*2-0.98; a2 = 2*0.98*0.965+0.965^2; a3 = -0.98*0.965^2; sys.F sys.G sys.H sys.T = = = = [-a1 -a2 -a3; 1 0 0; 0 1 0]; [K; 0; 0]; [0 0 1]; 5*60; % Input constraints sys.umin = -10; sys.umax = 100; % State constraints sys.xBoundaries = Polyhedron('A',[eye(3);-eye(3)],'b',30*[1 1 1 1 1 1]); % Output constraints sys.ymin = -Inf; sys.ymax = Inf; max_it = 100; lambda = 1; time = cputime; [InvSet, N] = InvariantSetDiscreteLTI(sys,lambda,max_it); display(['Total time elapsed ', num2str(cputime-time)]) InvariantSetDiscreteLTISystem.m function [ InvSet, N ] = InvariantSetDiscreteLTI( sys, lambda, max_it ) %INVARIANTSETDISCRETELTI Calculates the invariant (or contractive) set for %the system sys % Calculate the invariant or contractive set for the discrete LTI system % sys. % % sys - Description of discrete LTI system, with dynamics % x(t+1) = F*x(t) + G*u(t), y(t) = H*x(t). % The struct sys must be on the following form: % sys.F % sys.G % sys.H % sys.umin % sys.umax B.2 % % % % % % % % % % % % 55 Matlab code sys.xBoundaries sys.ymin sys.ymax sys.T MPT Polyhedron object Sampling time lambda - contraction paramterer when calculating contractive sets. lambda = 1 gives an invariant set. max_it - number of iterations before the calculation is stopped (if no convergence before) Andreas Svensson, 2013 F = sys.F; G = sys.G; % Define dimensions for the problem dim = size(F,1); %Preallocate structs X = cell([1,max_it]); X_ = cell([1,max_it]); X__ = cell([1,max_it]); A = cell([1,max_it]); b = cell([1,max_it]); A{1} = [sys.xBoundaries.A; sys.H; -sys.H]; b{1} = [sys.xBoundaries.b; sys.ymax; -sys.ymin]; X{1} = Polyhedron('A',A{1},'b',b{1}); %STEP 1 X{1}.plot; title('Iteration number 1') UcA = [eye(size(G,2)); -eye(size(G,2))]; Ucb = [sys.umax; -sys.umin]; i = 2; exit = 0; while (i > max_it) + exit == 0 %Exit conditions: Either fix point, or too many iterations tic [¬,n] = size(A{i-1}); X_{i} = Polyhedron('A', [A{i-1}*F A{i-1}*G; zeros(2,n), UcA],... 'b',[lambda*b{i-1};lambda*Ucb]); %STEP 2, augmentation of the state space and the input space X__{i} = X_{i}.projection(1:dim); %STEP 3, Projection X{i} = X{i-1}.intersect(X__{i}); %STEP 4, Intersection A{i} = X{i}.A; b{i} = X{i}.b; X{i}.plot; title(['Iteration number ', num2str(i)]) pause(0.05) exit = X{i-1}==X{i}; %STEP 6 Check if fix point is reached T = toc; 56 B Implementation display(['Iteration number ', num2str(i), ': ', num2str(T)]) i = i+1; end N = i-1; InvSet = X{1}; for i = 2:N InvSet = InvSet.intersect(X{i}); end if exit title([num2str(i-1), ', Invariant set found! :-)']) display('Invariant set found! :-)') else title([num2str(i-1), ', Terminated early. :-(']) display('Terminated early. :-(') end end B.2.2 Code for Example 2.10 This Matlab computes the robust periodic contractive set described in Example 2.10, and performs a simulated mpc control using this. The code automatically produces Figure 2.10 and a simplified version of Figure 2.9. The code is, however, relatively general in its notation, and should be straight forward to modify for related problems. This code requires mpt3 and Gurobi, see Section B.1, clear time = cputime; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% Define system struct 'sys' %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % % % % % System dynamics on the form x(t+1) = sys.F*x(t) + sys.G*u(t) + sys.E*d(t), y(t) = sys.H*x(t) where y(t) is the system output, u(t) is the control input and d(t) are the disturbances. sys.F sys.G sys.E sys.H = = = = [1.67 -0.71; 1 0]; [1; 0]; -1*eye(2); [17.3 -199]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% Define problem struct 'prob' %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % The struct prob defines the entire problem. Each feature (except prob.P) B.2 57 Matlab code % is an array with prob.P cells, one for each corresponding time step %%% Define periodicity of period P prob.P = 24; %%% Define state bounds on the form Ax*x < bx, for each time step % Upper and lower bound on y; lby < y < uby uby = [10 30 lby = [10 30 10 30 10 30 10 30 10 30 10 30 10 30 10 30 10 30 10 30 10 30 30 30 30 30 5 30 0 30 30 30 10 30 30 30 5 0 30 30 30 30 30 ... 30]; 30 ... 30]; % Distribute uby and lby to corresponding time step, using y = H*x % together with the ultimate bounds on x [-1e5, 1e5]*[-1e5, 1e5] prob.Ax = cell(prob.P,1); [prob.Ax{1:prob.P}] = deal([sys.H; -sys.H; eye(2); -eye(2)]); prob.bx = cell(prob.P,1); for i = 1:prob.P prob.bx{i} = [uby(i); lby(i); 1e5*[1; 1; 1; 1]]; end %%% Define input bounds for each step on the form Au*u < bu prob.Au = cell(prob.P,1); [prob.Au{1:prob.P}] = deal([1; -1]); prob.bu = cell(prob.P,1); [prob.bu{1:prob.P}] = deal([0.5; 0.5]); %%% Define disturbance bounds for each step on the form Ad*d < bd prob.Ad = cell(prob.P,1); [prob.Ad{1:prob.P}] = deal([eye(2); -eye(2)]); prob.bd = cell(prob.P,1); [prob.bd{1:prob.P}] = deal(0.01*[1; 1; 1; 1]); prob.bd{7} = deal(0.25*[1; 1; 0; 0]); prob.bd{12} = deal(0.25*[1; 1; 0; 0]); prob.bd{16} = deal(0.1*[1; 1; 0; 0]); prob.bd{20} = deal(0.25*[1; 1; 0; 0]); prob.bd{22} = deal(0.1*[1; 1; 0; 0]); %%% Define contraction factor prob.lambda = 0.99; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% Define simulation parameters %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% Number of simulated periods sim.n = 3; %%%%%%%%%%%%%%%%%%%% %%% End of input %%% %%%%%%%%%%%%%%%%%%%% % Add some state space model information in struct sys 58 B Implementation sys.n = size(sys.F,2); sys.un = size(sys.G,2); sys.dn = size(sys.E,2); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% Calculate invariant set %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Initialize for h = prob.P:-1:1 S{h,1} = Polyhedron('A',prob.Ax{h},'b',prob.bx{h}); end set_found = false; % Initialization is_empty = false; % Initialization max_it = 300; % Maximum number of iterations before ending prematurely k = -prob.P-1; %Initialize step counter while -k < max_it + prob.P tic && ¬set_found && ¬is_empty lcurrent = floor( ( - k - 1) / prob.P ) + 1; % Index of periodic round hcurrent = k + lcurrent * prob.P + 1; % Index within periodic round ktarget = k + 1; % Target step ltarget = floor( ( - ktarget - 1 ) / prob.P ) + 1; % Target index htarget = ktarget + ltarget * prob.P + 1; % Target index Atarget = S{htarget,ltarget}.A; % Get the set matrices for the target btarget = S{htarget,ltarget}.b; % step Ap = S{hcurrent,lcurrent-1}.A; % Get the set matrices for the step in bp = S{hcurrent,lcurrent-1}.b; % the previous periodic round [m,n] = size(Atarget); [mp,np] = size(Ap); % Calculate the erosion for robustification against disturbances erosion = zeros(m,1); for j = 1:m prob.A = sparse(prob.Ad{hcurrent}); prob.lb = -Inf*ones(sys.dn,1); prob.rhs = prob.bd{hcurrent}; prob.sense = '<'; prob.obj = -Atarget(j,:)*sys.E; params.outputflag = 0; y = gurobi(prob,params); erosion(j) = -y.objval; end % Define the augmented set S_ S_ = Polyhedron('A',... [Atarget*sys.F Atarget*sys.G; ... zeros(2*sys.un,sys.n), prob.Au{hcurrent}],... 'b', [prob.lambda*btarget-erosion;prob.lambda*prob.bu{hcurrent}]); B.2 59 Matlab code % Project S_ on the dimensions belonging to the state space Sp = S_.projection(1:sys.n); if Sp.isEmptySet is_empty = true; end % Intersect current step with corresponding step from previous round S{hcurrent,lcurrent} = Polyhedron('A',[Ap; Sp.A],'b',[bp; Sp.b]); % Define set for control constraints C{htarget,ltarget} = Polyhedron('A',S{htarget,ltarget}.A,... 'b',S{htarget,ltarget}.b-erosion); figure(1) S{hcurrent,lcurrent}.plot title([num2str(hcurrent),', ',num2str(lcurrent),... ' (', num2str(k+prob.P),')']) pause(0.05) if S{hcurrent,lcurrent}.isEmptySet is_empty = true; elseif lcurrent > 2 && S{hcurrent,lcurrent-1} == S{hcurrent,lcurrent} set_found = true; end T = toc; display([num2str(hcurrent),', ',num2str(lcurrent),... ' (', num2str(k+prob.P),') ', num2str(T), ' s.']) k = k-1; end k = k+1; if set_found display('Invariant set found!') elseif is_empty display('Empty set; no periodically invariant set exists!') clf else display('Terminated early') end display(['Total time elapsed ', num2str(cputime-time)]) if ¬set_found return end % End script %%%%%%%%%%%%%%%%%%%%%%%%%% %%% Plot invariant set %%% %%%%%%%%%%%%%%%%%%%%%%%%%% % Size of subplot plot p = ceil(sqrt(prob.P)); q = p - (p*(p-1) ≥ prob.P); 60 B Implementation figure(2), clf for n = k:1:k+prob.P-1 lp = floor((-n-1)/prob.P)+1; hp = n + lp*prob.P + 1; subplot(p,q,prob.P+k-n) dim = S{hp,lp}.Dim; if dim == 2 % Plot set (MPT plot function is not supported in subplots) hull = convhull(S{hp,lp}.minVRep.V(:,1),S{hp,lp}.minVRep.V(:,2)); fill(S{hp,lp}.minVRep.V(hull,1),S{hp,lp}.minVRep.V(hull,2),'r'); elseif dim == 3 % Plot set (MPT plot function is not supported in subplots) M = S{hp,lp}.minVRep.V; hull = convhull(M(:,1),M(:,2),M(:,3)); trisurf(hull,M(:,1),M(:,2),M(:,3),'Facecolor','r') end title([num2str(hp),', ',num2str(lp), ' (', num2str(n+prob.P),')']) end InvSet = cell(1,prob.P); for j = 1:prob.P subplot(p,q,j) if isempty(S{j,end}) InvSet{j} = S{j,end-1}; else InvSet{j} = S{j,end}; end dim = InvSet{j}.Dim; if dim == 2 % Plot set (MPT plot function is not supported in subplots) hull = ... convhull(InvSet{j}.minVRep.V(:,1),InvSet{j}.minVRep.V(:,2)); fill(InvSet{j}.minVRep.V(hull,1),InvSet{j}.minVRep.V(hull,2),'r'); elseif dim == 3 % Plot set (MPT plot function is not supported in subplots) M = InvSet{j}.minVRep.V; hull = convhull(M(:,1),M(:,2),M(:,3)); trisurf(hull,M(:,1),M(:,2),M(:,3),'Facecolor','r') end title(num2str(j),'FontWeight','bold') end % Extract the set aimed for control constraints SetConstr = cell(1,prob.P); for j = 1:prob.P if isempty(C{j,end}) SetConstr{j} = C{j,end-1}; else SetConstr{j} = C{j,end}; end end %%%%%%%%%%%%%%%%%%% %%% Simulations %%% %%%%%%%%%%%%%%%%%%% B.2 Matlab code 61 N = sim.n*prob.P; % Set random seed rng(1,'v5uniform') % Preallocations x = zeros(sys.n,N+1); y = zeros(1,N); u = zeros(1,N); d = zeros(2,N); % Initializatiopns x(:,1) = [0,0]; t = 1:N; for i = t j = mod(i,prob.P); jplus1 = j+1; j = j + (j==0)*10; % Define index % Generate disturbances according to description % (assuming given on a form with Ad = [eye(2); -eye(2)]) b = prob.bd{j}; d(:,i) = (b(1:2)+b(3:4)).*(rand(2,1)-0.5) + (b(1:2)-b(3:4))/2; % Calculate output y(1,i) = sys.H*x(:,i); % MPC problem with prediction horizon 1 and no output pr.obj = 0; pr.Q = sparse(1); pr.A = sparse([SetConstr{jplus1}.A*sys.G; prob.Au{j}]); pr.sense = '<'; pr.rhs = [SetConstr{jplus1}.b-SetConstr{jplus1}.A*(sys.F*x(:,i));... prob.bu{j}]; pr.lb = -Inf; pr.ub = Inf; params.outputflag = 0; s = gurobi(pr,params); % Simulate system u(i) = s.x; x(:,i+1) = sys.F*x(:,i) + sys.G*s.x + sys.E*d(:,i); end % Plotting figure(3) subplot(2,1,1) limh = plot(t,repmat(uby,1,sim.n),'Color',0.8*[1 1 1],'LineWidth',2); hold on; plot(t,-repmat(lby,1,sim.n),'Color',0.8*[1 1 1],'LineWidth',2); yh = plot(t,y,'k'); disth = stem(t,150*(d(1,:)+d(2,:))-100,'k-','fill','marker','.'); baseline_handle = get(disth,'BaseLine'); set(baseline_handle,'BaseValue',-100) hold off; axis([t(1) t(end) -100 100]) legend([yh disth limh],'Output','Disturbances','Output boundaries') title('Output trajectory and disturbances'); xlabel('Time') subplot(2,1,2) stem(t,u,'k','fill'); axis([t(1) t(end) -0.4 0.4]) title('Control input'); xlabel('Time') C Experimental invariant set illustration It is, at least for the author, not easy to get an intuitive feeling for the meaning of the shape of an invariant set of dimension N when it is displayed as a N -dimensional figure, as in, e.g., the 2-dimensional set in Figure 2.10 or the 3dimensional ser in Figure 4.1a. The plotting of an invariant set of dimension N is also limited, by obvious reason, to N ≤ 3 However, the states of the chosen state space description of the model, presented in Example 2.1, are possible to interpret as future outputs, i.e., x3 (t) = y(t), x2 (t) = y(t + T ) and x1 (t) = y(t + 2T ), according to the model structure. This idea of interpreting the states as consecutive outputs is utilized in Figure C.1. The plotting shows, for one given x3 value at a time (with tick marker x3), the projection of the x1 - x2 plane onto the x2 axis (with tick marker x2), and finally, for every x2 value, the projection on the x1 axis of the intersection of the given x2 and x3 values (with tick marker1 x1). 1 Since this is plotted for a lot of different x values at the same x tick on the axis, the association 2 1 between the x1 and x2 values are lost. This plot hence only gives an idea of the, in some sense, extreme bounds on the invariant set. Since the set is convex, however, the extreme values of x2 are associated with the extreme values of x1 , which guarantees that the plot is not misleading in that sense. Figure C.1: Experimental illustration of the invariant set 63 64 C Experimental invariant set illustration Three consecutive output values at the times t, t + T and t + 2T from the model are possible to interpret as the state of the model at time t. Although the practical use of this might be limited, it can still provide an intuitive idea of how steep the output curves can go2 , according to the invariant set. This way of plotting an invariant set does, indeed, not give full information of the invariant set, since only a sample of it is plotted. 2 However, for the states at time t + T to still be within the set, if they were at time t, the control action taken at time t had to be appropriate. Bibliography T Ajanki. Historien om diabetes och insulinets upptäckt. Svenska diabetesförb., Sundbyberg, 1999. Cited on page 23. F Blanchini. Ultimate boundedness control for uncertain discrete-time systems via set-induced Lyapunov functions. In Decision and Control, 1991., Proceedings of the 30th IEEE Conference on, pages 1755–1760. IEEE, 1991. Cited on page 16. F Blanchini. Set invariance in control. Automatica, 35(11):1747–1767, 1999. Cited on pages 2, 12, and 51. F Blanchini and S Miani. Set-theoretic methods in control. Birkhäuser Boston, 2007. Cited on page 12. F Blanchini and W Ukovich. Linear programming approach to the control of discrete-time periodic systems with uncertain inputs. Journal of optimization theory and applications, 78(3):523–539, 1993. Cited on pages 2, 12, 17, and 51. S Boyd and L Vandenberghe. Convex optimization. Cambridge university press, 2004. Cited on pages 11 and 12. G Danaei, MM Finucane, Y Lu, G M Singh, M J Cowan, C J Paciorek, J K Lin, F Farzadfar, YH Khang, GA Stevens, et al. National, regional, and global trends in fasting plasma glucose and diabetes prevalence since 1980: systematic analysis of health examination surveys and epidemiological studies with 370 country-years and 2· 7 million participants. The Lancet, 378(9785):31–40, 2011. Cited on page 1. E Dassau, B W Bequette, B A Buckingham, and F J Doyle III. Detection of a meal using continuous glucose monitoring implications for an artificial β-cell. Diabetes Care, 31(2):295–300, 2008. Cited on page 28. E Dassau, F Cameron, H Lee, B W Bequette, H Zisser, L Jovanovič, H P Chase, D M Wilson, B A Buckingham, and F J Doyle III. Real-time hypoglycemia prediction suite using continuous glucose monitoring. Diabetes Care, 33(6): 1249–1254, 2010. Cited on page 28. 65 66 Bibliography E Dassau, H Zisser, R A Harvey, M W Percival, B Grosman, W Bevier, E Atlas, S Miller, R Nimri, L Jovanovič, et al. Clinical evaluation of a personalized artificial pancreas. Diabetes care, 36:801–809, 2013. Cited on page 26. C E T Dórea and J C Hennet. (A, B)-invariant polyhedral sets of linear discretetime systems. Journal of optimization theory and applications, 103(3):521–542, 1999. Cited on pages 2 and 12. L Eiselein, H J Schwartz, and J C Rutledge. The challenge of type 1 diabetes mellitus. ILAR Journal, 45(3):231–236, 2004. Cited on page 1. C Ellingsen, E Dassau, H Zisser, B Grosman, M W Percival, L Jovanovič, and F J Doyle III. Safety constraints in an artificial pancreatic β cell: An implementation of model predictive control with insulin on board. Journal of diabetes science and technology (Online), 3(3):536, 2009. Cited on page 28. M Enqvist, T Glad, S Gunnarsson, P Lindskog, L Ljung, J Löfberg, T McKelvey, A Stenman, and J E Strömberg. Industriell reglerteknik kurskompendium. Division of Automatic Control, Department of Electrical Engineering, Linköping University, Sweden, 2010. Cited on page 9. T Glad and L Ljung. Control theory. CRC Press, 2000. Cited on pages 5, 8, 9, and 19. T Glad and L Ljung. Reglerteknik : grundläggande teori. Studentlitteratur, Lund, 4., [omarb.] uppl. edition, 2006. Cited on page 5. R Gondhalekar and C N Jones. MPC of constrained discrete-time linear periodic systems-a framework for asynchronous control: Strong feasibility, stability and optimality via periodic invariance. Automatica, 47(2):326–333, 2011. Cited on page 17. R Gondhalekar, F Oldewurtel, and C N Jones. Least-restrictive robust periodic model predictive control applied to room temperature regulation. Automatica, 2013. doi: 10.1016/j.automatica.2013.05.009. Cited on pages 2, 17, and 20. R Greenberg. Turning diabetes over to the bionic pancreas. Huffington Post, 04 2013. URL http://www.huffingtonpost.com/riva-greenberg/ diabetes-clinical-trial_b_3110140.html. Cited on page 2. B Grosman, E Dassau, H Zisser, L Jovanovič, and F J Doyle III. Zone model predictive control: A strategy to minimize hyper- and hypoglycemic events. Journal of diabetes science and technology (Online), 4(4):961–975, 2010. Cited on pages iii, v, 1, 26, 28, 46, and 53. R A Harvey, Y Wang, B Grosman, M W Percival, W Bevier, DA Finan, H Zisser, DE Seborg, L Jovanovic, F J Doyle III, et al. Quest for the artificial pancreas: combining technology with treatment. Engineering in Medicine and Biology Magazine, IEEE, 29(2):53–62, 2010. Cited on page 26. C Hellerström. Diabetes : forskningen, framstegen, framtiden. Vetenskapsrådet, Stockholm, 2002. Cited on page 23. Bibliography 67 L M Huyett, R A Harvey, W Bevier, H Zisser, L Jovanovic, E Dassau, and F J Doyle III. Closed-loop control and quality of continuous glucose monitor data based on prospective closed-loop clinical trial, March 2013. Presented at the 6th International Conference on Advanced Technologies and Treatments for Diabetes. Cited on page 24. B P Kovatchev, M Breton, C Dalla Man, and C Cobelli. In silico preclinical trials: a proof of concept in closed-loop control of type 1 diabetes. Journal of diabetes science and technology, 3:44–55, January 2009. Cited on pages iii, v, 2, 25, and 37. M Kvasnica, P Grieder, M Baotić, and M Morari. Multi-parametric toolbox (MPT). Hybrid Systems: Computation and Control, pages 121–124, 2004. Cited on page 53. J M Maciejowski. Predictive control with constraints. Pearson education, 2002. Cited on page 9. L Magni, D M Raimondo, C Dalla Man, M Breton, S Patek, G De Nicolao, C Cobelli, and B P Kovatchev. Evaluating the efficacy of closed-loop glucose regulation via control-variability grid analysis. Journal of diabetes science and technology (Online), 2(4):630, 2008. Cited on pages 39 and 53. K J Åström and R M Murray. Feedback Systems: An Introduction for Scientists and Engineers. Princeton University Press, 2010. URL http://www.cds. caltech.edu/~murray/amwiki/index.php/Main_Page. Cited on page 5. W J Rugh. Linear system theory. Prentice-Hall, Inc., 1996. Cited on pages 5, 7, 8, and 32. S J Russell, F H El-Khatib, DM Nathan, K L Magyar, J Jiang, and E R Damiano. Blood glucose control in type 1 diabetes with a bihormonal bionic endocrine pancreas. Diabetes care, 35(11):2148–2155, 2012. Cited on page 24. J Skyler, editor. Atlas of diabetes. Springer Verlag, 5th edition, 2012. Cited on pages 2, 23, 24, and 26. G M Steil, K Rebrin, J Mastrototaro, B Bernaba, and M F Saad. Determination of plasma glucose during rapid glucose excursions with a subcutaneous glucose sensor. Diabetes technology & therapeutics, 5(1):27–31, 2003. Cited on page 24. K van Heusden, E Dassau, H C Zisser, D E Seborg, and F J Doyle III. Controlrelevant models for glucose control using a priori patient characteristics. Biomedical Engineering, IEEE Transactions on, 59(7):1839–1849, 2012. Cited on pages iii, v, 1, 2, 8, 25, 26, and 48. R D Yates and D J Goodman. Probability and stochastic processes. John Wiley & Sons, 2nd edition, 1999. Cited on page 6. 68 Bibliography H C Zisser, T S Bailey, S Schwartz, R E Ratner, and J Wise. Accuracy of the seven® continuous glucose monitoring system: Comparison with frequently sampled venous glucose measurements. Journal of diabetes science and technology (Online), 3(5):1146, 2009. Cited on page 24. Index artificial pancreas, 2, 26 arx, 25 abbreviation, xi asymmetric costs, 12 input signals, 6 insulin on board, 28, 37 insulin-glucose model, 25 invariant sets, 12–17 computation, 51–52, 54–56 constraints, 35–36 contractive, 16 controlled, 15 glucose-insulin model, 31 glucose-insulin model, robust, 32– 34, 46 maximum, 17 periodic controlled, 17 positively, 13 robust controlled, 16 basal rate, 25 bolus, 26 cgm, 24 abbreviation, xi errors, see measurement errors computations, 53–61 control signals, 6 csii, 24 abbreviation, xi cvga, 39 abbreviation, xi limitations, 2 Luenberger observer, 9 disturbance signals, 6 -implementation, 35–36, 46 Matlab, 53 Matlab code, 53–61 fda, see U.S. Food and Drug Adminis- meals tration announced, 26 feedback control, 9 detection, 28 robustification, 34, 48 glucagon, 24 unannounced, 27 Gurobi, 53 measurement errors, 24, 32–34, 48 model mismatch, 48 hard constraint, 35 mpc, 9–12 hyperglycemia, 23 abbreviation, xi hypoglycemia, 23, 28 example, 17–21 zone, see zmpc implementation, 53 mpt3, 53 in silico, 2, 37–46 69 70 night time zone, 28, 37 objective, 2 output signals, 6 pendulum, 13 phase plan, 14 plant, 6 plant-model mismatch, 26, 35 poles, 7 results, 37–46 safety system, 37 safety systems, 28 simulations, 25, 37–46 simulator, 25, 53 soft constraint, 35 state estimation, 9 state space form, 6 t1dm, 1, 23–25 abbreviation, xi tdi, 25 abbreviation, xi θ-implementation, 36–37 transfer function, 7 U.S. Food and Drug Administration, 25 ucsb/Sansum group, 1 UVa/Padova simulator, see simulator z-transform, 7 zeros, 8 zmpc, 12 abbreviation, xi Index Upphovsrätt Detta dokument hålls tillgängligt på Internet — eller dess framtida ersättare — under 25 år från publiceringsdatum under förutsättning att inga extraordinära omständigheter uppstår. Tillgång till dokumentet innebär tillstånd för var och en att läsa, ladda ner, skriva ut enstaka kopior för enskilt bruk och att använda det oförändrat för ickekommersiell forskning och för undervisning. Överföring av upphovsrätten vid en senare tidpunkt kan inte upphäva detta tillstånd. All annan användning av dokumentet kräver upphovsmannens medgivande. För att garantera äktheten, säkerheten och tillgängligheten finns det lösningar av teknisk och administrativ art. Upphovsmannens ideella rätt innefattar rätt att bli nämnd som upphovsman i den omfattning som god sed kräver vid användning av dokumentet på ovan beskrivna sätt samt skydd mot att dokumentet ändras eller presenteras i sådan form eller i sådant sammanhang som är kränkande för upphovsmannens litterära eller konstnärliga anseende eller egenart. För ytterligare information om Linköping University Electronic Press se förlagets hemsida http://www.ep.liu.se/ Copyright The publishers will keep this document online on the Internet — or its possible replacement — for a period of 25 years from the date of publication barring exceptional circumstances. The online availability of the document implies a permanent permission for anyone to read, to download, to print out single copies for his/her own use and to use it unchanged for any non-commercial research and educational purpose. Subsequent transfers of copyright cannot revoke this permission. All other uses of the document are conditional on the consent of the copyright owner. The publisher has taken technical and administrative measures to assure authenticity, security and accessibility. According to intellectual property law the author has the right to be mentioned when his/her work is accessed as described above and to be protected against infringement. For additional information about the Linköping University Electronic Press and its procedures for publication and for assurance of document integrity, please refer to its www home page: http://www.ep.liu.se/ © Andreas Svensson