Genetic algorithm tuning of artificial pancreas MPC with individualized models

(1)

Genetic algorithm tuning of artificial pancreas MPC with individualized models

Olov Sehlin

Engineering Physics and Electrical Engineering, master's level 2019

Luleå University of Technology

Department of Computer Science, Electrical and Space Engineering

(2)

(3)

L

ULEÅ UNIVERSITY OF TECHNOLOGY

DEPARTMENT OFCOMPUTERSCIENCE, ELECTRICAL AND SPACE ENGINEERING

ENGINEERING PHYSICS ANDELECTRICAL ENGINEERING

ELECTRONIC SYSTEM ANDCONTROL ENGINEERING

Genetic algorithm tuning of artificial pancreas MPC with individualized models

Author:

Olov Sehlin

Supervisor:

Simone Del Favero

Examiner:

Khalid Tourkey Atta

June 30, 2019

(4)

(5)

Foreword

When I was five years old I got Type 1 Diabetes. I have never felt that I would like to use an insulin pump since they always been manual and clumsy. When Damiano Varagnolo told me about Simone Del Favero and the research they are doing in Padova it was not much to think about. To do my master thesis in this subject was a dream for me. Combining my life long experience with diabetes and my knowledge about control theory to help millions in their battle against diabetes.

This was the perfect end for my five years at Luleå University of Technology.

Why the thesis was made in Padova and together with Simone had three main reasons. In Padova, they have good knowledge about both diabetes and control theory, much collected data and the UVA/Padova Simulator. It also had personal reasons, for me it was an opportunity to expanding my contact network, exploring the world and of course help myself in my own battle against diabetes.

I would like to thank Damiano Varagnolo and Khalid Atta for being very supportive and patient at all times. Always pushing me to my limit to increase my confidence and knowledge and making me believe in my self.

Simone del Favero for letting me come to Padova and supervise my thesis. Supporting me in our shared interests but also my own interests. Placing me in a lab together with other great minds that made me feel welcome from day one.

Last but not least, my family and especially my wife. This would not have been possible without Emmelie by my side in Padova. The most supportive person I know who left Sweden and her job to just to live with me in Italy and support me.

(6)

(7)

Abstract

Diabetes is a growing chronic disease and a worldwide problem. Without any available cure in sight for the public other methods needs to be applied to increase the life quality of diabetic patients. Artificial Pancreas (AP), a concept of having a closed loop system to control the glucose level on Type 1 Diabetes (T1D) patients has been introduced and is under development.

In this thesis, Model Predictive Control (MPC) has been re implemented from scratch in MAT- LAB/SIMULINK with associated Kalman filter and prediction function. It was implemented in the latest version of the UVA/Padova Simulator which is a tool approved by FDA for simulating diabetes treatment in order to speed up the AP development.

Different MPC cost functions where tested together with integral action on a simplified system using a linear approximation of a population model. It was implemented and tuned with a new simulation tuning method using Genetic Algorithm (GA). It showed that the quadratic cost function without integral action was the best with respect to performance and time efficiency. 3 hours was the best prediction horizon and was used for the individualized tuning using the University of Virginia (UVA)/Padova simulator.

For the individualized MPC, models identified by the University of Padova were used. These simulations showed that an individualized model could be used for improved T1D treatment compared to an average population model even though the results were mixed. Almost all of the patients got improved treatment with the closed treatment and non hypoglycemic event occurred.

The identification of better models is a great challenge for the future development of the AP MPC due to the excitation problems.

(8)

(9)

List of Figures

2.1 Block diagram of general treatment of T1D. . . 4 2.2 The model of the U.S Food and Drug Administration (FDA) accepted simulator.

The white block are the from the s2008 version, the gray were updated in the s2013 version and the black were completely new to the s2013 version. . . 6 2.3 The structure of the latest version (s2017) of the UVA/Padova simulator. . . 7 2.4 Concept plot of CVGA. . . 16 2.5 The argument of the Hyper + Hypo² weighted Mean (HH²_wM) cost function in

Eq.2.32 for the GA tuning. . . 16 3.1 MPC verification with three different weights at time k and prediction horizon

N= 5 hours and a meal of 50g of carbohydrates. Also showing the concept of MPC. 21 3.2 MPC state weights found for all prediction horizons using the both tuning fitness

functions. . . 22 3.3 Fitness values obtained with the weights in Fig. 3.2 and all prediction horizons

and applied on its own and the other fitness function. . . 23 3.4 Mean glucose with the obtained weights for all prediction horizon for and both

fitness functions. . . 24 3.5 Minimum and maximum glucose with the obtained weights for all prediction hori-

zon for and both fitness functions. . . 24 3.6 Weights found for all prediction horizons using the both tuning fitness functions

and integral action. . . 25 3.7 Fitness values for all weights and prediction horizons with integral action and

fitness values obtained from its own fitness function. . . 26 4.1 Downsample strategy. . . 27 4.2 Downsampled example with the blue dot original signal, blue split signal and red

downsampled signal. . . 28 5.1 Three day simulation for patient 23 with glucose and insulin from all three methods. 34 5.2 The Hyper + Hypo²weighted CVGA (HH²_{wCV GA}) fitness value from all three cases

where a score below the target line means 100 % time in target. . . 35 5.3 Parallel coordinate comparison of the HH²_{wCV GA}fitness value at an individual level. 37

(12)

LIST OF FIGURES

5.4 Percentage spent in target for the three methods and all 53 patients. . . 37 5.5 Distributions of all glucose values from the three days. . . 38 5.6 Mean glucose value for all 53 patients and all three methods. . . 38 E.1 One day simulation with all prediction horizons and weights obtained with both

fitness functions. . . . G.1 Three day simulation for patient 1 with glucose and insulin from all three methods.

G.2 Three day simulation for patient 2 with glucose and insulin from all three methods.

(13)

List of Tables

4.1 Changes made in units and signals in the simulator to fit the used individualized models. . . 28 5.1 The tests that were preformed and the number of passing models for each test.

After failing a test they were directly removed from the set. . . 33 5.2 Selected subjects for simulation plots. . . 34 5.3 The number of subjects out of the 53 with 100% time in target and the number

of the subjects that got improved treatment vs open loop and the average model tuning. The "100% in target" and "outside target" is refereed to the open loop treatment results. . . 36 D.1 The states used in the average model with description and units. . . . D.2 Diagonal elements of process noise covariance matrix QKFand measurement noise matrix RKFused in the kalman filter, grouped by units. . . .

(14)

(15)

Acronyms

AP Artificial Pancreas

ARMAX AutoRegressive Moving Average eXogenous CGM Continuous Glucose Monitoring

CVGA Control Variability Grid Analysis FDA U.S Food and Drug Administration

FL Fuzzy Logic

GA Genetic Algorithm

HH²_wM Hyper + Hypo²weighted Mean HH²_{wCV GA} Hyper + Hypo²weighted CVGA IG Interstitial Glucose

in silico computer simulation in reference to biological experiments LQR Linear Quadratic Regulator

MIMO Multiple Input Multiple Output MPC Model Predictive Control

NICE National Institute for Health and Care Excellence PID Proportional Integral Derivative

T1D Type 1 Diabetes T2D Type 2 Diabetes UVA University of Virginia

(16)

(17)

Chapter 1 Introduction

Diabetes is a common and growing chronic disease characterized by elevated levels of glucose in the bloodstream. Diabetes can generally be divided into two different classes: T1D which is characterized by reduced or non existing production of insulin while Type 2 Diabetes (T2D) is characterized by being resistant to insulin in combination with reduced production. T2D is the most common [1]. In 2014 the estimated number of people living with diabetes was 422 million

5

worldwide compared to 108 million in 1980. 8.5 percent of the adult population in 2014 which has nearly doubled from 4.7 percent in 1980. In 2012, 1.5 million deaths were directly connected to diabetes and the same year, diabetes was estimated to be the eighth leading cause of death and the fifth leading cause of death among women [2].

Without any cure for diabetes for the public in sight in the nearest future, other methods are

10

needed to improve the quality of life for those who are living with diabetes. AP, a concept in- vented to have a closed loop system built to mimic the insulin delivery of a non diabetic person [3]. Several method had been proposed earlier, including various version Proportional Integral Derivative (PID), Fuzzy Logic (FL) and MPC control schemes [4].

Regardless of the chosen control algorithm, the tuning is very important for robust control at

15

an individual level. Finding the optimal prediction horizon and weights for the MPC cost function in clinical trials might harm or even kill the patient. Instead, computer simulation in reference to biological experiments (in silico) trials can be made on a personalized level to build a base of parameters that can be used in a regression analysis to find parameters for clinical trials [5]. These simulations can be made in the UVA/Padova simulator which is a FDA approved simulator for

20

testing T1D treatment. The simulator emulates meals and uses the carbohydrate and insulin to blood glucose dynamics at an individualized level to simulate T1D treatment.

Finding the in silico tuning parameters is an optimization problem and the Control Variability Grid Analysis (CVGA) metric is a one well known and used method for tuning and assessing the performance [6] of glucose control. In this thesis, GA will be used to find the optimal tuning of the

25

MPC to see if the current individualized models identified by the University of Padova can be used to improve T1D treatment using the latest version of the UVA/Padova simulator. The MPC will be implemented from scratch using the theory to create own functions instead of using the already existing SIMULINK functions. Different versions of the MPC cost functions such as quadratic,

1

(18)

CHAPTER 1. INTRODUCTION

linear and zone together with integral action will be tested to see if the glucose control could be improved.

In the background chapter, the background of diabetes, artificial pancreas, and the UVA/Padova simulator will be presented together with the theory regarding MPC, kalman filter and the tuning of the controller.

5

In the third chapter, "Testing the MPC on a linear approximation of the average patient" the theory from the background section was implemented and tested on the linear average diabetes model. The MPC, kalman filter and the human system were all based on the same linear model and no noise were used. The GA was used to tune the MPC.

Then the changes made in the UVA/Padova simulator together with the implementation veri-

10

fication of the MPC are presented in the fourth chapter, "UVA/Padova simulator MPC implementation".

Chapter five, MPC on the virtual non linear patients using the UVA/Padova simulator, is the chapter where the UVA/Padova simulator was used to tune and simulate patients on an individual level. The individualized models were compared to individual tuning of the average model and

15

open loop treatment.

Finally, in the last chapter, conclusion, the thesis is summarized and suggestions for future development are discussed.

2

(19)

Chapter 2 Background

This chapter will cover all the background regarding this thesis. From diabetes and the treatment of the disease to the theory of MPC and Kalman filter and also the tuning of the controller.

A new way of in silico tuning using optimization algorithms will also be presented together with the well known CVGA metric and a newly proposed metric based on the authors own experience with diabetes.

5

2.1 Diabetes

Diabetes is a chronic disease where the body is not able to control the level of glucose in the bloodstream by itself. When not treated properly with insulin, diabetes is characterized by elevated levels of glucose in the blood (hyperglycemia, above 180 mg/dl) and can over time cause serious complications on the heart, blood vessels, eyes, kidneys, nerves and limb losses. Insulin is a

10

hormone and works like a key for the glucose into the cells and allows the body to use the energy from meals. If to much insulin is injected the glucose levels get lower (hypoglycemia, below 70 mg/dl) and life threatening situations can occur. The disease is usually divided into two categories:

Type 1 Diabetes (T1D) and Type 2 Diabetes (T2D), where the second is the most common. In T2D, the body is getting resistant to insulin or not producing enough insulin. This is often caused by an

15

unhealthy lifestyle and most common among adults. T1D, on the other hand, is an autoimmune disease and the opposite of T2D. In T1D, the immune system of the diseased person has attacked the beta cells in the islets of Langerhans, placed on the pancreas, causing a decreased or non existing production of insulin [7–9]. T1D is sometimes referred to as juvenile diabetes but it can occur in all ages.

20

2.1.1 Treatment

Diabetes is a very individual disease which means that every patient has to use a treatment method that works in their life and lifestyle. One common thing for every patient when treating T1D is the need to measure the glucose levels in the blood to know how to dose the insulin. This is recommended before and after every meal but also when expecting risks of hypo or hyperglycemia.

25

Different methods that can be used, either the old classical method where the patient pick a finger, 3

(20)

CHAPTER 2. BACKGROUND

Human system

Blood glucose Disturbance

Glucose monitoring Insulin

Insulin decision

Meals Meal

Measured glucose

Fig. 2.1: Block diagram of general treatment of T1D.

arm or another part of your body with a sharp needle. Then put a drop of blood on the test stick. The other common and more modern method is using a Continuous Glucose Monitoring (CGM) system which continuously can monitor your glucose levels. These systems are usually less accurate but show trends and patterns to help the patient understanding the glucose levels better [10].

5

The other common part for all T1Dpatients is the need of manually injecting insulin. This can be via injections or with help of an insulin pump. In the first case with injection, there are some different kinds of insulin that are being used. Fast acting insulin in connections with meals or when needed, called bolus. Medium or long acting insulin once or twice a day to keep a stable level during the whole day, called basal insulin. This means between 4 − 8 injections a normal day. The

10

insulin pump works different, it only uses fast acting insulin and is continuously injecting small doses of insulin that works similarly as the medium or long acting insulin [11, 12]. Fig. 2.1 shows a block diagram a general treatment of diabetes. The disturbance could be sensitivity variability or some sort of activity that the person performs. The insulin decision could either be made by the patient via manual injections, an insulin pump or an AP.

15

2.2 Artificial Pancreas

The traditional treatment of T1D is indeed a challenge in the everyday life for the patient. It requires many important decisions every day and could potentially be very stressful. Stress could possibly make the treatment even harder were other hormones could affect the glucose levels [13].

To reduce the burden of the patient, a new and more effective method is under development but

20

also some products are already approved for commercial distribution [14]. This method is called Artificial Pancreas (AP) which is a closed loop system built to mimic how the body works for a non-diabetic person. The AP is based on three main part:

• Continuous Glucose Monitoring: As mentioned earlier, it continuously measures the glucose values in the blood. This information is used for feedback for the control algorithm.

25

• Control algorithm: This is what makes this method for glucose control unique and more effective. Instead of constant decision taking from the patient, a control algorithm is used

4

(21)

to control the glucose levels. Different methods have been tested to control the glucose levels. Non model required controller such as PID has been tested but due to the complex dynamics of the body with different impulse responses for insulin and carbohydrate and the hard constraints of the system [4, 15]. Instead, model based controller such as MPC, which uses the model dynamics to predict the future, can be used. In any of the cases, the closed

5

loop algorithm can be combined with open loop meal boluses [16]. After calculating the next control insulin action, the control algorithm sends the information to the insulin pump.

• Insulin pump: Instead of having a constant and manually set basal insulin the AP uses the instructions from the control algorithm to control the glucose levels in the bloodstream.

2.3 UVA/Padova Simulator

10

During the last few years, the development of artificial pancreas control algorithms has been ac- celerating thanks to computer simulations. It has provided a more cost and time efficient way of preclinical trials that also avoid animal trials.

The first version of the UVA / Padova T1D simulator (s2008) was approved by FDA in 2008 and it was the first computer based model accepted as a substitute of preclinical trials of closed

15

loop insulin treatments. It consisted of 300 in silico patients (100 adults, 100 adolescents, 100 children) where each virtual patient was presented by a parameter vector which was randomly extracted from an appropriate joint parameter distribution. The simulator emulate meals and included features like glucose, glucose rate of appearance, glucose utilization, renal extraction, and insulin rate of insulin appearance from the subcutaneous tissue and insulin degradation [17]. This

20

version can be summarized by the white and grey blocks in Fig. 2.2

The first update of the simulator (s2013) was released in 2013 and accepted by FDA the same year. The grey blocks in Fig. 2.2 was improved from the s2008 version and the black blocks were added. Some important improvements were added to the s2013 version such as the joint parameter distribution, the definition of relevant clinical parameters and strategy for the virtual

25

patient generation. Improved rules of determining carbohydrate ratio and the correction factor was implemented to better comply with clinical trials. Also, non linearity response for hypoglycemia glucose treatment was added [18].

The latest and most recent version of the simulator is the s2017 version. It incorporates new models of time varying parameters such as intraday variability of insulin sensitivity and dawn

30

phenomenon. It also has new insulin deliver routes and models of the Dexcom^©G5 Mobile CGM device and the Bayer^©Contour Next USB self monitoring blood glucose device. The most important features of the s2017 simulator are shown in Fig. 2.3.

The latest version also has the same structure as the general treatment in Fig. 2.1. The insulin decision block is divided into two parts. First, it calculates a meal bolus based on the current

35

glucose value, the size of the meal and insulin on board. Then the meal bolus can be used as a feed forward in the closed loop algorithm which is the other part. This feed forward is also referred to as open loop treatment.

In the glucose monitoring block, different monitoring methods can be selected. It includes the CGM model and also the Interstitial Glucose (IG) which the CGM is based on and also pure blood

40

5

(22)

28

Journal of Diabetes Science and Technology 8(1)

was modified as follows (see appendix for the complete set of model equations).

Glucagon secretion and kinetics. Glucagon kinetics are

described with a 1-compartment linear model

⁸

:

H t( )= − ⋅n H t( )+SR t_H( ) H( )0 =H_b

(5) with H(t) plasma hormone concentration, SR

_H(t) glucagon

secretion (SR

_H^b

its basal value), and n clearance rate.

Glucagon secretion is described as the sum of 2 components:

SR t_H( )=SR t_H^s ( )+SR t_H^d( )

(6) with

SRHs (t)

SRHs (t) max 2 G h G(t) SRHb ,0

=

− ⋅^ρ _^ −

(

^σ ⋅[ ^t − ]+

)

^_^{if G t}^{( ) G}≥Gb

t f G t

− ⋅ − ⋅ −

+ +

[ ]



 





 



ρ σ

SRHs (t) max G h G(t)

I(t) 1 SRHb ,0 ( ))<







 ^Gb

(7)

with G(t) plasma glucose (G

_b

its basal value), I(t) plasma insulin concentration (I

_b

its basal value), σ and σ

₂

alpha-cell responsivity to glucose level, 1/ρ delay between static gluca- gon secretion and plasma glucose. In this way, static secre- tion is stimulated when G<G

_b

(but modulated by insulin) and inhibited when G≥G

_b.

The second component is proportional to glucose rate of change,

SR_H^d( ) max ( ) ,

t dG t

= ⋅ − dt

 



δ 0

(8) with dG(t)/dt glucose rate of change, and δ alpha-cell respon- sivity to glucose rate of change.

Model parameters σ and δ also reflect that glucagon response declines with age of diabetes (see also the “New Joint Parameter Distribution” section below).

Of note, in real life, glucagon secretion is almost certainly dependent on insulin level in the alpha cells (paracrine effect), not in the circulation. However, it is very difficult to model the intrapancreatic levels, so the use of plasma insulin was the best we could do, even if not perfectly physiologic.

Glucagon action. The model of glucagon action is based

on the EGP model reported in equation (4), but assumes that EGP is stimulated by the above basal glucagon concentration with some delay,

EGP t( )=k_p₁−k_p₂⋅G t_p( )−k_p₃⋅X t^L( )+ ⋅ξ X^H( )t

I t’( )= − ⋅k_i

[

I t’( )−I t( )

]

I’( )⁰ ⁼I_b

X t^L( )= − ⋅k_i X t^L( )−I t’( ) X^L( )0 =I_b

X t k X t k H t H

X

H

H H

H b

H

( ) ( ) max ( ) ,

( )

= − ⋅ + ⋅ _

(

−

)

_

=

0

0 0

(9)

with H(t) plasma glucagon concentration, X

^H(t) delayed

glucagon action on EGP, ξ liver responsivity to glucagon, and 1/k

_H

delay between glucagon concentration and action.

Glucose utilization in hypoglycemia. The model of glucose uti-

lization reported in equations (2) and (3) is unable to describe well the hypoglycemic range likely due to an inadequate description of insulin action, which paradoxically increases when glucose decreases under a given threshold. This phe- nomenon has been described during hyperinsulemic clamps in T1DM.

^9,10

Based on these observation we developed a new

Figure 1. Scheme of the model included in the FDA-accepted

T1DM simulator. White blocks are the unit processes of S2008⁴ (gastro-intestinal tract, glucose kinetics and insulin kinetics); gray blocks are those that have been updated in the S2013 to account for counterregulation (liver, muscle, and adipose tissue); black blocks are new (alpha cell, glucagon kinetics, and delivery).

Fig. 2.2: The model of the FDA accepted simulator. The white block are the from the s2008 version, the gray were updated in the s2013 version and the black were completely new to the s2013 version.

glucose.

The meal block emulating meal which has settings such as carbohydrate counting error to give the insulin decision block another meal size than the human system. It also has feedback from the glucose value for hypoglycemic treatments.

Lastly the human system block. It has the differential equation that describes insulin and

5

carbohydrate responses. The disturbance in the simulator could be some sensitivity variability.

2.4 Model Predictive Control

Model Predictive Control (MPC) is today a well used model based control method. It originates from the late seventies and has been under constant development since then. MPC has many

6

(23)

Visentin et al 275

subjects during early morning, due to both an increased endogenous glucose production (EGP) and an increased insulin requirements.^29,30 Recently, Mallad et al²⁹ quantified an almost 30 mg/dL increase in glucose concentration from 3:00 am to 7:00 am and observed an increase of EGP of about 1.5 mg/kg/min in the same interval. Thus, we modeled the EGP variation as a linear increase of EGP at zero glucose and insulin (represented in the model by parameter kp1 of equation (A5)) from 3:00 ^am to 7:00 ^am. Similarly, the increased insulin requirement is described by a decrease in insulin-dependent glucose utilization (through a parameter k_ir of equation (A10)). Intersubject variability of both time interval and magnitude of k_p1 increase and k_ir decrease are obtained as random modulation of the averages reported in Mallad et al.²⁹

Updated Glucagon Model. With respect to S2013, S2017 incorporates two updates in the model of glucagon secretion, kinetics and action: first, glucagon clearance, originally expressed as absolute rate (min^–1), is now fractional, that is, per unit of distribution volume (mL/kg/min), in order to pro- vide an easy comparison with the literature;³¹ second, no static secretion is assumed when blood glucose is over its basal value (see equation (A25)), in order to avoid nonphysi- ological oscillations in simulated glucagon concentration.

Insulin Delivery Models

Subcutaneous administration of fast-acting insulin analogs is part of conventional insulin therapy in T1D, and is thus a fundamental component of the T1D simulator. However, in order to accommodate alternative insulin delivery routes, S2017 also incorporates models of intradermal and inhaled insulins (Figure 2).

Subcutaneous Fast-acting Insulin. S2017 incorporates an improved model of subcutaneous insulin kinetics, which is able to reproduce the commercially fast-acting insulin analogs.³² This model has been identified on a large dataset of 116 T1D subjects and introduces, with respect to that included in S2008 and S2013, a subject-specific delay in insulin absorption (τ), accounting for both removal of auxil- iary substances for insulin storage in the vial and diffusion processes occurring in the subcutis (see equations A16, A17).

Intradermal Insulin. The model of intradermal insulin route is based on Lv et al.³³ This model (equations A18, A19) has been identified on 10 healthy subjects, and assumes that insulin deliv- ered at the intradermal site is then transported to the bloodstream via two independent routes, that is, a diffusion-like process and a combination of a diffusion-like processes followed by an additional compartment before entering the blood.

Figure 2. Scheme of T1D S2017 highlighting the introduced upgrades. Specifically, the model incorporates time-varying parameters describing intraday S_I variability (k_p3, V_mx) and dawn phenomenon (k_p1, k_ir); insulin delivery routes account for administration of subcutaneous fast-acting insulin, intradermal and inhaled insulins; glucose monitoring devices include models of both CGM and SMBG.Fig. 2.3: The structure of the latest version (s2017) of the UVA/Padova simulator.

advantages and these are some that stand out.

• It uses an optimized control action based on the predicted future and unlike the PID controller which only uses information from the past. This can easily be explained by con- trolling the speed and steering of a car. The driver uses the information seen through the windscreen to plan the driving. If an animal or pedestrian walks out on the road ahead, the

5

driver makes a new prediction/plan to act on while a PID only can see through the driving mirrors and what is happening on the road behind the car. This is also shown in Fig. 3.1 where different weights yielding different control signals and different outcomes.

• A broad use of applications. From process industries to a wide sense of robotics, to chemical processes or as in this case, medical applications.

10

• Can handle Multiple Input Multiple Output (MIMO) system.

• Takes constraints into account when optimizing the control action.

• Reference tracking when the future reference is known and also handling known disturbances.

7

(24)

However, it also has some drawbacks. The main drawback is the computational power that is needed for the optimization problem. However, the development of cheaper, smaller and more powerful computers is constantly decreasing this drawback. It also acquires more complex deriva- tions than other control methods and full state feedback [19].

This section will deal with the concepts and theory behind MPC.

5

2.4.1 Prediction

As stated in the name of MPC, the prediction part is the main idea behind the method. It allows the control algorithm to see the future and make decisions based on the past, present and future.

This can be done in continuous time for non linear systems as well as for discrete time linear systems and all combinations of them [20, 21]. It will be explained further for the discrete time

10

linear case.

Let the discrete time linear system be described by:

x[k + 1] = Ax[k] + Bu[k] + Md[k]

y[k] = Cx[k] (2.1)

where x ∈ Rⁿ is the system states, u ∈ R^m^u is the system inputs, y ∈ R^p is the system outputs and d ∈ R^m^d is a known disturbance at the current time instance but possible even for future time instances as well. It is assumed that the system is fully controllable and observable. Also, it is

15

assumed that the system does not have any direct input-output relation, hence no D matrix in the Eq.(2.1).

The future state predictions can be described by:

x[k + N] = A^Nx[k] +

N

∑

i=1

A^Ni(Bu[k + i1] + Md[k + i1]) (2.2)

It can be proved by induction. To solve the optimization problem the prediction formula in Eq.(2.2) can be written on a matrix form:

20

X =Asx[k] +BsU+MsD (2.3)

where

X =x[k + 1]^T, x[k + 2]^T. . . x[k + N − 1]^T, x[k + N]^TT

∈ R^Nn U =u[k]^T, u[k + 1]^T. . . u[k + N − 2]^T, u[k + N − 1]^TT

∈ R^Nm^u D=d[k]^T, d[k + 1]^T. . . d[k + N − 2]^T, d[k + N − 1]^TT

∈ R^Nm^d

(2.4)

8

(25)

and

As=





 A A²

... A^N−1

A^N







Bs=







B 0 0 . . . 0

AB B 0 . . . 0

A²B AB B . . . 0

... . .. ...

A^N−1B A^N−2B . . . AB B







Ms=







M 0 0 . . . 0

AM M 0 . . . 0

A²M AM M . . . 0

... . .. ...

A^N−1M A^N−2M . . . AM M







(2.5)

whereAs∈ R^Nnxn,Bs∈ R^NnxNm^u andDs∈ R^NnxNm^d. For output prediction, the same state space equation can be used together with Eq.(2.1) and the following equation is obtained:

y[k + N] = CA^Nx[k] +

N

∑

i=1

CA^Ni(Bu[k + i1] + Mu[k + i1]) (2.6) with the same matrix form as before:

Y =Ax[k] +BU+MD (2.7)

where

5

A=





 CA CA²

... CA^N−1

CA^N







B=







CB 0 0 . . . 0

CAB CB 0 . . . 0

CA²B CAB CB . . . 0

... . .. . ..

CA^N−1B CA^N−2B . . . CAB CB







M=







CM 0 0 . . . 0

CAM CM 0 . . . 0

CA²M CAM CM . . . 0

... . .. . ..

CA^N−1M CA^N−2M . . . CAM CM







(2.8)

whereA_{∈ R}^{N p},B_{∈ R}^{N pxNm}^u andD_{∈ R}^{N pxNm}^d. 9

(26)

2.4.2 Cost function

The MPC algorithm is an optimization problem where a cost function based on the receding horizon prediction is minimized. At each time instance the future control sequence in Eq.(2.4) is computed and minimized according to the cost function design. It can take different shapes de- pending on the objectives and use cases. Usually, a quadratic cost function is used and in the easy

5

case with just regulation to zero without any known disturbances the cost function is:

J=

N

∑

k=1

y^T[k]qy[k] + u^T[k]ru[k] (2.9)

or in the matrix form:

J=Y^T[k]QY[k] +U^T[k]RU[k] (2.10)

andY is the same as in Eq. (2.7) without the known disturbance and X andU from Eq.(2.4) andA andB from Eq.(2.8). The weights q and r are the design parameters together with the prediction horizon N. q weights the error in reference tracking. Large q allows big control action

10

in the input to full fill the tracking requirement and vice versa. In the matrix case the weights will be as follows:

Q=







q 0 0 . . . 0 0 q 0 . . . 0

... ... . .. ...

0 0 . . . q 0 0 0 . . . 0 q







R=







r 0 0 . . . 0 0 r 0 . . . 0

... ... . .. ...

0 0 . . . r 0 0 0 . . . 0 r







(2.11)

With Eq.(2.10) and Eq.(2.11) the cost function can be written as:

J= (AX +BU)^TQ(AX+BU) +U^TRU

= x^T[k]A^TQAx[k] + x^T[k]A^TQBU+U^TB^TQBU+U^TB^TQAx[k] +U^TRU

= 2x^T[k]A^TQBU+U^T(B^TQB+ R)U

(2.12)

In the exact same way the cost function for reference tracking, feed forward input and with known disturbances can be derived starting from:

15

J=

N

∑

k=1

(x[k] − x₀[k])^Tq(x[k] − x₀[k]) + (u[k] − u₀[k])^Tr(u[k] − u₀[k]) (2.13)

10

(27)

This time, lets use the state prediction instead of output prediction from Eq.(2.6) but for output prediction. This will yield new weighting matrices and introduce a new weight matrix ˆT:

Qˆ =







C^TqC 0 0 . . . 0

0 C^TqC 0 . . . 0

... ... . .. . ..

0 0 . . . C^TqC 0

0 0 . . . 0 C^TqC







Rˆ=







r 0 0 . . . 0 0 r 0 . . . 0

... ... . .. ...

0 0 . . . r 0 0 0 . . . 0 r







Tˆ =







C^Tq 0 0 . . . 0 0 C^Tq 0 . . . 0

... ... . .. ...

0 0 . . . C^Tq 0 0 0 . . . 0 C^Tq







(2.14)

resulting in the final form cost function for reference tracking:

J= 2((x^T[k]A^T_s +D^TM^T_s) ˆQBs+Y₀T^ˆBs+U₀R)^ˆ U

+U^T(B_s^TQ^ˆB_s+ ˆR)U ^(2.15)

where the reference output vector is:

Y₀=y₀[k + 1]^T, y₀[k + 2]^T. . . y₀[k + N − 1]^T, y₀[k + N]^TT

∈ R^{N p} (2.16) 2.4.3 Quadratic programming

5

The final form of the cost function can be recognized as a quadratic programming optimization problem which is formulated as this in its general form:

minimize

x f(x) := x^TQx+ c^Tx subject to A_eqx= beq

Ax≤ b lb≤ x ≤ ub

x∈ R^Nm

(2.17)

Finding the closed form solution for this problem is easy without the presence of constraints.

∂ f (x)

∂ x = 0 ⇒ x = −c^T/Q⁻¹= K1x[k] + K2 (2.18) which will just be normal state feedback with a additional part from the the known disturbance, reference tracking and feed forward. When dealing with constraints which is one of the advantages

10

11

(28)

with MPC, more advanced optimization algorithms are needed. In MATLAB, quadprog can be used [22].

U= quadprog(QQP, cQP, A, b, Aeq, beq, lb, ub) (2.19) In both cases,U =u[k]^T, u[k + 1]^T. . . u[k + N − 2]^T, u[k + N − 1]^TT

will be calculated but only the first will be used as the control action at each time instance k.

2.4.4 Integral action

5

The standard MPC does not have any integral control action unless the system has an integrator itself with the argument that it is not necessarily optimal [23]. To get the desired offset free reference tracking an integrator has to be augmented into the system. There are some different ways of doing this but some might not always give the expected results. However, two methods will be derived here to show two cases where it will work and not. In both methods, the state space

10

model in Eq.(2.1) will be considered.

Method 1

The first method is using ∆u = u[k] − u[k − 1] as input and with that input the original state space model will be without the known disturbance.

x[k + 1] = Akx[k] + Bku[k] = Akx[k] + Bk(u[k − 1] + ∆u[k]) (2.20) In order to write it on a complete state space form the previous input u[k − 1] must be added as a

15

state as follows:

x_m1[k] =x[k]^T, u[k − 1]^TT

(2.21) The augmented state space model using the first method can now be expressed as follows with the known disturbance matrix added as well:





 x[k + 1]

u[k]





=







A B

0 I







| {z }

A˜





 x[k]

u[k − 1]





+





 B

I







| {z }

B˜

∆u +





 M

0







| {z }

M˜

d

y=

C 0

| {z }

C˜





 x[k]

u[k − 1]







(2.22)

It can be showed from the eigenvalues of the new state space matrix ˜A that system has been augmented with a discrete pole in 1 which was the objective. However, it can be showed in the

20

closed loop system without the presence of an observer that the integrator will disappear, hence it will have direct feedback around the integrator. Why it works with the observer in the loop can be explained by the estimated ˆu[k − 1] not necessarily being equal to the actual u[k − 1] in the presence of unknown disturbances [23].

12

(29)

Method 2

Instead of differentiating the input, the state vector can be differentiated:

x[k + 1] − x[k] = Ak(x[k] − x[k − 1]) + Bk(u[k] − u[k − 1]) = Ak∆x[k] + Bk∆u[k] (2.23) and note that the known disturbance is not added yet. The same is done for the output.

∆y = y[k + 1] − y[k] = CAk∆x[k] + CBk∆u[k] (2.24) In order to write it on a complete state space form the current output y[k] must be added as a state as follows:

5

x_m2[k] =

∆x[k]^T, y[k]^TT

(2.25) The augmented state space model using the second method can now be expressed as follows with the known disturbance matrix added as well:







∆x[k + 1]

y[k + 1]





=







A 0

CA I







| {z }

A˜







∆x[k]

y[k]





+





 B CB







| {z }

B˜

∆u +





 M CM







| {z }

M˜

∆d

y=

0 I

| {z }

C˜







∆x[k]

y[k]







(2.26)

In the same way for this method it can be showed that the augmented system now has a discrete pole in 1 which again was the objective. This method will work both with and without an observer in the closed loop. It should be noted that the known disturbance ∆d now also are differentiated

10

[24].

2.5 Kalman filter

As stated in Eq.(2.2), the MPC algorithm need full state feedback. This is usually not possible, therefore an observer is needed in the closed loop. It was assumed in 2.4.1 that the used model is fully observable. Since most of the models have uncertainties and noise or incomplete mea-

15

surements, a kalman filter is a good choice for state estimation [25, 26]. For the kalman filter, the following state space model is considered:

x[k + 1] = Akx[k] + Bku[k] + Mkd[k] + Nkw[k]

y[k] = Ckx[k] + v[k] (2.27)

Where w[k] ∈ Rⁿ and v[k] ∈ R^p are zero mean white noise of the process and measurements respectively with the following covariances:

E











 w[k]

v[k]





w^T[k], v^T[k]





=





 Q_k S_k S_k^T R_k





 (2.28)

13

(30)

The presence of correlated process and measurement noise Sk 6= 0 makes the kalman filter sub optimal [27]. It can be showed [27, 28] that formulating the prediction step as follows:

ˆ

x[k|k − 1] = A_k−1x[k − 1|k − 1] + Bˆ _k−1u[k − 1] + M_k−1d[k − 1]

+ ¯G_k−1(y[k − 1] −Ckx[k − 1|k − 1])ˆ

P_k|k−1= (A_k−1− ¯G_k−1Ck)P_k−1|k−1(A_k−1− ¯GkCk)^T + ¯Q_k−1− N_k−1S_k−1R¯⁻¹_k−1S^T_k−1N_k−1^T

(2.29)

where

Q¯k= NkQkN_k^T R¯k= Rk

G¯k= NkSkR¯⁻¹_k

(2.30)

will make the process and measurement noise uncorrelated and optimal estimation is achieved again. The kalman filter algorithm then follows by the measurement update.

5

K_k= P_k|k−1C_k^T(CkP_k|k−1C_k^T+ ¯R_k)⁻¹

˜

y[k] = y −Ckx[k|k − 1]ˆ ˆ

x[k|k] = ˆx[k|k − 1] + Kky[k]˜ P[k|k] = (I − KkC_k)P_k|k−1

(2.31)

where ˆx[k|k] can be used for state feedback.

2.6 MPC tuning

The tuning of the MPC can be made in several ways. In the most simple cases where only one parameter, like the state weight q, is to be tuned methods such as line search can be successfully used in a time efficient way. When the complexity increases with more tuning parameters the

10

usage of optimization algorithms can be used. There are several optimization algorithms available in MATLAB, GA, Particle Swarm, fmincon and quadprog to mention some of them. Different optimization algorithm has different strengths and weaknesses and needs to be selected for the specific problem. Using this kind of tuning is only possible in in silico trials since they might be very dangerous and also time consuming if they were not dangerous. The tuning also needs to be

15

assessed and evaluated with a metric designed for the specific case.

In this case, where the MPC will be used for T1D control, the tuning at an individual level is very important since every patient is different. Some patients might be able to use very aggressive control while others might get harmed with aggressive control.

2.6.1 Genetic algorithm

20

Genetic Algorithm (GA) is one of many optimization methods based on evolutionary biology [29].

It is a combination of survival of the fittest, a central concept Darwinian evolutionary theories describing the natural selection and some random elements added to that [30]. The main idea behind the method is to have a selected number of children searching for the global minimum of

14

(31)

the selected fitness function. The first generation is randomized. The next generation is then based on the three main parts of the GA, the selection, crossover and mutation and they are explained as follows:

• Selection: This is the part of the survival of the fittest strategy. In each generation, a selected number of "elite children" are chosen. These are the children with the best fitness value and

5

their chromosomes will be copied to the next generation.

• Crossover: This occurs between randomly selected parents and genetic information is combined to generate a new child for the next generation. How much genetic information that will be used from each of the parents is randomized from a selected probability distribution.

• Mutation: To generate diversity between one generation and another, mutated chromo-

10

somes are used. This gives a new, completely random element to the search of the global minimum of the fitness function.

This will run for a fixed number of generations or until a stall criterion is fulfilled. Exactly how to set parameters such as crossover rate, stall conditions, numbers of generations, etc in the MATLAB implementation can be found at [31].

15

2.6.2 Tuning metrics

Previously the Control Variability Grid Analysis (CVGA) metric has been widely used for tuning the artificial pancreas controllers. The idea behind this method is to minimize the distance between the highest and lowest glucose values. A concept plot of the CVGA is shown in Fig. 2.4 where the distance d should be minimized. The x-axis shows the minimum glucose value and the y-axis

20

shows the maximum glucose value. It should be noted that the x-axis is reverted and that the y-axis is not linear. It is fitted to a 3rd degree polynomial to penalize a maximum of 180 mg/dl equally as a minimum of 90 mg/dl and the same for 300 and 70 and so on. This means that abnormalities on the minimum glucose are penalized more than on the maximum glucose. Even if the CVGA is widely used, it only considers the extreme point and not the time spent in different regions. A

25

good CVGA does not necessarily mean that the patient is having a good mean glucose and could potentially spend much time in hyperglycemia. As an alternative to the CVGA tuning method, a new fitness function was created. This function was created to maximize the time in the target and will be called Hyper + Hypo²weighted Mean (HH²_wM):

HH²_wM= hyper + hypo²+ mean(Glucose)W (2.32) Where, hyper = time samples in hyperglycemia, hypo = time samples in hypoglycemia. This

30

method has similarities with the CVGA since both penalize low values more since they can potentially be life threatening. The argument of this fitness function can be explained by Fig. 2.5.

Where the both red lines having the same time in target but the since the lower spent time in the hypoglycemia zone it will be penalized more. The both blue lines have the same time in target but to distinguish between them the element of the mean value was added as a final decider. The

35

mean value was multiplied by W = 1/100 dl/mg since the target range of glucose is in the range of 10²mg/dl and then the mean value element will just act like one more sample outside the target

15

(32)

>110 90 70 <50

<110 180 300

>400

A Upper B Upper C

Lower B B Upper D

Lower C Lower D

E

A Upper B Upper C

Lower B B Upper D

Lower C Lower D

E

Minimum BG

MaximumBG

A=0 B=1 C=0 D=0 E=0

d

Fig. 2.4: Concept plot of CVGA.

0 2 4 6 8 10 12 14 16 18

50 100 150 200

Samples[units]

Glucose[mg/dl]

HH²_wMcost function design

Target range

Fig. 2.5: The argument of the HH²_wM cost function in Eq.2.32 for the GA tuning.

16

(33)

range and not affect the result in the big picture. It should be noted that this is an empirical fitness function based on the diabetic experience of the author with the main objective to maximize the time in target while at the same time avoid hypoglycemic events.

17

(34)

18

(35)

Chapter 3 Testing the MPC on a linear approximation of the average patient

This chapter will present the average model and the results found using GA as the tuning method.

This stage were implemented as intermediate step between the theory and using the UVA/Padova simulator. It was used to investigate the possibility of different cost functions for the MPC and the usage of integral action. Also to test if the GA tuning method was working and and learning how T1D control could behave before moving on to the simulator and more complex setups.

5

The MPC and the subject system were both based on the average linear model and no noise were used to give an easier first step of T1D control using MPC.

3.1 Average linear model of T1D population

The average linear model was derived from 100 adults where parameters such as age, weight, basal insulin and glucose, carbohydrate ratio among others were averaged to generate an average

10

non linear model [32]. Then, the non linear model was linearized around an average equilibrium point of the population. Finally, it was discretized at a sample time of Ts= 5 minutes and put on the following state space form:

x[k + 1] = Ax[k] + Bu[k] + Md[k]

y[k] = Cx[k] (3.1)

where

• x[k] is the state vector presented in Table. D.1 in appendix D.

15

• u[k] = i[k] - u₀is the difference between the bolus insulin i[k] and the basal insulin u₀[k] and measured in pmol/min/kg, normalized over the patient body weight.

• d[k] is the meal acting as a disturbance and measured in mg/min.

• y[k] = CGM[k] - y0 is the difference between the measured values from the CGM and the basal glucose level y₀and is measured in mg/dl.

20

19

(36)

CHAPTER 3. TESTING THE MPC ON A LINEAR APPROXIMATION OF THE AVERAGE PATIENT Since only the linear average model will be used, it will be referred to as the average model.

The kalman filter covariances used for the average model can be found in Table. D.2.

3.2 MPC cost function

The first cost function considered for the MPC was the following quadratic cost function:

J=

N−1

∑

k=1

q(y[k] − y₀[k])²+ (u[k] − u₀[k])²+ ||x[N]||²_P (3.2) where P is the unique non negative solution if the discrete time Riccati equation:

5

P= A^TPA+ qC^TC+ A^TPB(1 + B^TPB)⁻¹B^TPA (3.3) and r = 1 to simplify the tuning. This structure of the cost has previously been used in diabetes MPC [6]. Using this different weight on the last state in the prediction horizon can be explained by comparing MPC to Linear Quadratic Regulator (LQR) and is used to increase stability for short prediction horizon [32]. To penalize inputs that would leave the system in a state where it can recover from. As mentioned in 2.4.3, this optimization can easily be solved with the quadprog

10

function in MATLAB.

The second cost function that was considered was the following linear cost function:

J=

N

∑

k=1

q|y[k] − y₀[k]| + |u[k] − u₀[k]| (3.4) This could potentially be solved using the linprog function in MATLAB but due to absolute value it gets very complicated. Instead fmincon was used to optimize the cost function. Since fmincon is not designed for this kind of optimization like quadprog is for the quadratic cost

15

function, the optimization took too much time. Therefore, the linear cost function was excluded.

3.3 Implementation verification

To verify that the implementation of the MPC and kalman filter was correct, a prediction test with different weights was made. One yielding the wanted result, one less aggressive in the state deviation and one more aggressive. A long prediction horizon of 5 hours was used for easier visual

20

inspection of the test. The result of the test can be found in Fig. 3.1 showing that the three different weights gave three different routes and that the prediction was good.

The MPC and kalman filter were also compared with the built in SIMULINK functions and they yielded the same results.

3.4 Simulation setup for the average model

25

The scenario used for simulating and tuning the average model: Three meals, breakfast at 07 : 00 with 50g of carbohydrates, lunch at 12 : 00 with 60g and dinner at 18 : 00 with 80g and the simulation was running over 24 hours starting at 04 : 00. Similar meal setups has been used before [6].

20

(37)

CHAPTER 3. TESTING THE MPC ON A LINEAR APPROXIMATION OF THE AVERAGE PATIENT

0

Time [steps]

Glucose[mg/dl]

MPC weight scenarios

Target Glucose Target weight Low weight High weight Past glucose

k + N

k k k + N

k + N k

0

Time [steps]

Insulin[pmol/min/kg]

Predicted insulin Low weight

0

Time [steps]

Predicted insulin Target weight

0

Time [steps]

Predicted insulin High weight Past Prediction horizon Unknown future

k + N k

Fig. 3.1: MPC verification with three different weights at time k and prediction horizon N = 5 hours and a meal of 50g of carbohydrates. Also showing the concept of MPC.

The meal bolus was calculated only from using the average model carbohydrate ratio. The meal was announced at the same time to the MPC and the virtual patient and the bolus was also given at the same time as the meal.

The MPC settings for the GA were set as follow: To find the best prediction horizon for the MPC, horizon lengths from 0.5 hours to 5 hours with steps of 30 minutes were used. The lower

5

bound of the prediction horizon was set shorter than necessary just to make sure all interesting horizons were tested but still longer than just one time sample. The upper boundary of 5 hours was decided by inspecting an impulse response with both insulin and carbohydrates showing that the states were affected for about that long. The boundaries for the MPC state weight q was just set to 10⁻⁸ as the lower bound a no upper bound was used. The lower bound was used to make sure

10

no negative values was tested and low values as 10⁻⁸behaved just as open loop. The lower bound for integral action was set to 10⁻¹⁰since it showed integral action preforms better with lower cost q. No upper bound was used since the GA can some times get stuck on the boundaries.

The GA settings were set as follows: The standard population size was used which is 50 when the number of variables ≤ 5. 50 generation was selected just to give a good result without taking

15

to much time. The standard setting for elite rate = 5% and crossover rate = 80% and no stall condition was used.

21

Genetic algorithm tuning of artificial pancreas MPC with individualized models

Genetic algorithm tuning of artificial pancreas MPC with individualized models

Olov Sehlin

L

Genetic algorithm tuning of artificial pancreas MPC with individualized models

Foreword

Table of contents

List of Figures

List of Tables

Acronyms

Chapter 1

Introduction

Chapter 2

Background

2.1 Diabetes

2.2 Artificial Pancreas

2.3 UVA/Padova Simulator

28

was modified as follows (see appendix for the complete set of model equations).

described with a 1-compartment linear model

:

(5) with H(t) plasma hormone concentration, SR

secretion (SR

its basal value), and n clearance rate.

Glucagon secretion is described as the sum of 2 components:

(6) with

(

)

(7)

with G(t) plasma glucose (G

its basal value), I(t) plasma insulin concentration (I

its basal value), σ and σ

alpha-cell responsivity to glucose level, 1/ρ delay between static gluca- gon secretion and plasma glucose. In this way, static secre- tion is stimulated when G<G

(but modulated by insulin) and inhibited when G≥G

The second component is proportional to glucose rate of change,

(8) with dG(t)/dt glucose rate of change, and δ alpha-cell respon- sivity to glucose rate of change.

Model parameters σ and δ also reflect that glucagon response declines with age of diabetes (see also the “New Joint Parameter Distribution” section below).

on the EGP model reported in equation (4), but assumes that EGP is stimulated by the above basal glucagon concentration with some delay,

[

]

(

)

(9)

with H(t) plasma glucagon concentration, X

glucagon action on EGP, ξ liver responsivity to glucagon, and 1/k

delay between glucagon concentration and action.

Based on these observation we developed a new

2.4 Model Predictive Control

Insulin Delivery Models

∑

∑

∑

∑

2.5 Kalman filter

2.6 MPC tuning

Chapter 3

Testing the MPC on a linear approximation of the average patient

3.1 Average linear model of T1D population

3.2 MPC cost function

∑

∑

3.3 Implementation verification

3.4 Simulation setup for the average model