Model Predictive Control with Invariant Sets in Artificial Pancreas for Type 1 Diabetes Mellitus

(1)

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 tekniska högskola

Linköpings universitet Linköpings universitet

(2)

(3)

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

(4)

(5)

Avdelning, Institution Division, Department

Avdelningen för reglerteknik Department of Electrical Engineering SE-581 83 Linköping Datum Date 2013-06-20 Språk Language Svenska/Swedish Engelska/English Rapporttyp Report category Licentiatavhandling Examensarbete C-uppsats D-uppsats Övrig rapport

URL för elektronisk version

http://www.ep.liu.se

ISBN — ISRN

LiTH-ISY-EX--13/4699--SE Serietitel och serienummer Title of series, numbering

ISSN —

Titel Title

Modellbaserad prediktionsreglering med invarianta mängder av artificiella pankreas för di-abetes typ 1

Model Predictive Control with Invariant Sets in Artificial Pancreas for Type 1 Diabetes Mel-litus Författare Author Andreas Svensson Sammanfattning Abstract

This thesis deals with Model Predictive Control (mpc) for artificial pancreas for Type 1 Dia-betes 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 dif-ferent ways to incorporate them into the mpc strategy. The work is an extension to the zone mpccontroller 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 vari-ability 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 ap-propriate 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 de-sign, with mpc based on invariant sets.

Nyckelord

(6)

(7)

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 in-varianta 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 ut-vecklats 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 si-mulator (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är-deringarna. 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.

(8)

(9)

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 knowl-edge 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 arti-ficial 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 re-sponses 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 pan-creas control.

Altogether, this thesis pinpoints a possible future direction of artificial pancreas control design, with mpc based on invariant sets.

(10)

(11)

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

(12)

(13)

Notation

General

Symbol Explanation

R Set of real numbers

∀ _{For all}

Abbreviations

Symbol Explanation

arx Auto-Regressive model with eXogenous input cgm Continuous Glucose Monitoring

cvga _{Control-Variability Grid Analysis}

csii _{Continuous Subcutaneous Insulin Infusion} mpc _{Model Predictive Control}

tdi _{Total Daily Insulin} t1dm _{Type 1 Diabetes Mellitus}

zmpc zone-mpc

(16)

(17)

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, meth-ods, 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}

(18)

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 well-functioning 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 im-plemented in the zone-mpc (zmpc) controller for artificial pancreas developed by the ucsb/ Sansum group, and evaluated in silico5 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 systemati-cally, 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 mod-els 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}

(19)

1.2 Outline of the thesis 3

1.2 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

sum-mary 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

com-putations, 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.

(20)

(21)

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 be-haves 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 chap-ter, 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 sys-tems, Rugh [1996] is recommended.

(22)

2.1 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 previous1inputs (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 at2time 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 interpre-tations 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}

(23)

2.1 Dynamical systems 7

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 one4input 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) + a1y(t − 1) + . . . + any(t − n) = b0u(t) + b1u(t − 1) + . . . + bmu(t − m) (2.3) (where the constants a1, . . . , anand b0, . . . , bmmight 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) + a1z −₁ Y (z) + . . . + anz −_n Y (z) = b0U (z) + b1· z −₁ U (z) + . . . + bmy −_m (z) (2.4a)

(where Y (z) is the z-transform of y(t), and similar for U (z)) the equation

Y (z) = b0+ b1z −₁ + . . . + bmz−m 1 + a1z−1+ . . . + anz−n U (z). (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}

(24)

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 in-jected 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

(z − 0.98)(z − 0.965)2U (z) (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 −_2.823 _0.913 1 0 0 0 1 0                 x1(t) x2(t) x3(t)         +         −_0.01 0 0         u(t) (2.6a) y(t) = 0 0 1         x1(t) x2(t) x3(t)         (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}

(25)

2.2 Model Predictive Control 9

2.1.3 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:

ˆ

x(t + 1) =A ˆx(t) + Bu(t) + K[y(t) − C ˆ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].

mpcis 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 mpca 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.

(26)

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 func-tion, 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 + n}_c−₁₎₎ _(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 cur-rent 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.

(27)

2.2 Model Predictive Control 11

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) = A2x(t) + ABu(t)) + Bu(t + 1) (2.9c)

Hence the following optimization problem can be formulated min u(t),u(t+1) X i=t+1,t+2 xT(i)Qx(i) + X i=t,t+1 uT(i)Ru(i) = (2.10a) min u(t),u(t+1)u

T_{(t)Ru(t) + u}T_{(t + 1)Ru(t + 1) + (Ax(t) + Bu(t))}T_{Q(Ax(t) + Bu(t))+}

(A2x(t) + ABu(t) + Bu(t + 1))TQ(A2x(t) + ABu(t) + Bu(t + 1)) = (2.10b)

= min u(t),u(t+1) u(t) u(t + 1) BT _(AB)T 0 BT ! Q _ABB 0_B ! + R 0 0 R !! u(t) u(t + 1) ! + 2 0 BT 0 0 ! Q A A2 ! x(t) !T u(t) u(t + 1) ! + xT(t) I AT Q I A ! x(t) (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

min u(t),u(t+1) 1 2U T_(BT mpcQBmpc+ R)U + (BTmpcQAmpcx(t))TU (2.10d) where U = _{u(t + 1)}u(t) ! , Ampc= _AA2 ! , Bmpc= _ABB 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].

(28)

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].

(29)

2.3 Invariant sets 13

2.3.1 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 exam-ple:

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 equa-tions on state space form, where the two states can be chosen as x1position and x2velocity, 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 there-after 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}

(30)

in the same way. Since there is a damping (e.g., friction) in the system, the max-imum distance to the downward position, as well as the maxmax-imum 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 ini-tial 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 corre-sponding sign convention as for the distances.

Figure 2.6 reveals that the pendulum’s state, once entering a circle9of 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}

(31)

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 con-trolled invariant10sets are defined for systems on the form x(t + 1) = f (x(t), u(t)), where x ∈ X ⊂ Rnand 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 implemen-tation 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

(32)

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 com-putations 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 con-trolled 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 ⊂ Rmand 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 con-trolled 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 exten-sion 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) ∈ U_{t modT} ⊂ Rmt modT and f is periodic such that f (t, x(t), u(t)) =

(33)

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 pe-riodic 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

x(t + 1) = 1.67 −0.71 1 0 ! x(t) + 1 0 ! u(t) + 1 0 0 1 ! d(t) (2.11a) y(t) = 17.3 −₁₉₉ _. _(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

(34)

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-turbance bound for the previous sample12_{, ˜}X_i+1 = {x ∈ R2|x + d ∈ X_i+1∀d ∈ D_i}. The set ˜X_i are thereafter used as a constraint for the corresponding predicted states.

Altogether, the optimization problem solved in each step looks like this: arg min

u(t) u

2_(t) _(2.12)

subject to umin< u(t) < umax

x(t + 1) ∈ ˜X_{t+1 modT}

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} stan-dard 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 pre-diction 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}

(35)

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 con-trolling 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 trajec-tory 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 insulin-glucose 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 distur-bances 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}

(36)

Figure 2.10: Result for a simulation of the system in Example 2.10, trolled by an mpc controller using the robust periodic invariant set as con-straint 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]

(37)

Figure 2.11: Result for a simulation of the system Example 2.10 using a standard mpc controller without knowledge of the disturbances.

(38)

(39)

3

Diabetes and the artificial pancreas

This chapter serves two purposes. First, to give an introduction to Type 1 Dia-betes 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 t1dmis referred to the literature, e.g., Skyler [2012]. For a popular science in-troduction 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 hyper-glycemia1_{. This can be dangerous and cause severe and life-threatening}

compli-cations to the patient if maintained for a long period of time.

A low blood glucose value is known as hypoglycemia2. Hypoglycemia may cause severe and life-threatening complications even if maintained only for a short pe-riod 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.}

(40)

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 to-day 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 asymmet-ric, 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 ar-tificial 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 nei-ther 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 read-ings, but the intravenous blood glucose level. The cgm measurements do not always coincide with the intravenous blood glucose level, for physiological rea-sons (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 mea-surements 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,

(41)

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 pan-creas has been developed at University of Virginia and University of Padova (Ko-vatchev 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 arxmodel) 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) = −5F · 1800 tdi (1 − p1)(1 − p2) 2 1 (z − p1)(z − p2)2 U (z) (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.

(42)

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 small-est 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 pan-creas. An overview of the basic structure of an artificial pancreas is given in Fig-ure 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 tech-nique 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 per-spective) 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 ‘counter-act’ the meal in advance and give a dose of insulin before or at the time the meal is eaten, a so called meal-bolus4. 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}

(43)

3.6 Artificial pancreas 27

Figure 3.2:An overview of the ucsb/Sansum group artificial pancreas sys-tem.

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 fore-hand. Depending on context, only the former may be called announced meal, but they are equivalent from a control algorithm perspective, since the informa-tion 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 pa-tient is having. This is the case considered in this thesis, although extensions to the announced case could be possible.

(44)

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 hav-ing information of the meal, as opposed not to have the information. It is, how-ever, also possible to discuss what the best is from a patient perspective. A com-parison 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, mea-surement errors and an aggressively tuned controller. It can especially be prob-lematic 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 Fig-ure 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.

(45)

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.

(46)

(47)

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 pur-pose of better intuitive understanding of its meaning is presented in Appendix C.)

The computations of the invariant set converges for the particular numerical val-ues 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 implemen-tation in Appendix B). The number of iterations, which also affects the computa-tion 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

(48)

(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 transfor-mations (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 influ-ence 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 in-put 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

(49)

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 measure-ment 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 statis-tical 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 Defini-tion 2.8. A robustificaDefini-tion according to DefiniDefini-tion 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,

(50)

Figure 4.3: The propagated measurement errors, with 90% of the simula-tions in the two-dimensional box.

Example 2.2), the bounds [−30, 30] will apply to those states as well. The robus-tification 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}

(51)

4.3 Implementation of set constraints in MPC 35

4.2.4 Plant-model mismatch

With a description for the plant-model mismatch, an invariant set approach ro-bust against plant-model mismatches could be a very interesting concept to in-vestigate. 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 op-timization, 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 infeasible2problem, 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 con-straints, which means that deviations from the set are allowed, but penalized in the optimization objective function. As opposed to this, a normal constraint im-plementation 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 + n}_c−₁₎₎ _(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 + n}_c−_{1)) + G(}_t+1_{, . . . ,}_t+n

c+1) (4.2)

subject to Ax(t + 1) ≤ b + t+1, . . . , Ax(t + nc) ≤ b + t+nc

t+1, . . . , t+nc≥0

(52)

Figure 4.4: An asymmetricly located rectangle described by Ax ≤ b, com-pared 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 + n}_c−_{1)) + G(θ}_t+1_{, . . . , θ}_t+n

c+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.

Model Predictive Control with Invariant Sets in Artificial Pancreas for Type 1 Diabetes Mellitus

Institutionen för systemteknik

Department of Electrical Engineering

Examensarbete

Model Predictive Control with Invariant Sets in

Artificial Pancreas for Type 1 Diabetes Mellitus

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

Sammanfattning

Abstract

Acknowledgments

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

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

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

3.6.1

Meals

3.6.2

Safety systems and hypoglycemia treatments

3.6.3