Proceedings of the 8th IFAC Symposium on Robust Control Design Proceedings of the 8th IFAC Symposium on Robust Control Design July 8-11, 2015. Bratislava, Republic Proceedings of the the 8th IFAC IFACSlovak Symposium on Robust Robust Control Control Design Design Proceedings of 8th Symposium on July 8-11, 2015. Bratislava, Slovak Republic Available online at www.sciencedirect.com July 8-11, 2015. Bratislava, Slovak July 8-11, 2015. Bratislava, Slovak Republic Republic
ScienceDirect IFAC-PapersOnLine 48-14 (2015) 082–087
Robust Gain-Scheduled Controller Design Robust Gain-Scheduled Controller Design Robust Gain-Scheduled Controller Design Robust Gain-Scheduled Controller Design for T1DM Individualised Model for T1DM Individualised Model for for T1DM T1DM Individualised Individualised Model Model Adrian Ilka, Tom´ aˇ s Ludwig, Ivan Ottinger, Mari´ an T´ arn´ık, Adrian Ilka, Tom´ a ˇ ss Ludwig, Ivan Ottinger, Mari´ a n T´ a rn´ık, Adrian Ilka, Tom´ a ˇ Ludwig, Ivan Ottinger, a n a Mikloviˇ cov´ a, Vojtech Vesel´ y, J´ anMari´ Murgaˇ sT´ AdrianEva Ilka, Tom´ a ˇ s Ludwig, Ivan Ottinger, Mari´ a n T´ arn´ rn´ık, ık, Eva Mikloviˇ c ov´ a , Vojtech Vesel´ y , J´ a n Murgaˇ s Eva Mikloviˇ c ov´ a , Vojtech Vesel´ y , J´ a n Murgaˇ s Eva Mikloviˇ cov´ a, Vojtech Vesel´ y, J´ an Murgaˇ s Institute of Robotics and Cybernetics, Institute of Robotics and Cybernetics, Institute of Cybernetics, Faculty of Electrical Engineering andand Information Technology, Slovak Institute of Robotics Robotics and Cybernetics, Faculty of Electrical Engineering and Information Technology, Slovak Faculty of of of Electrical Engineering and Information Information Technology, Slovak University Technology in Bratislava, Ilkoviˇcova 3,Technology, 812 19 Bratislava, Faculty Electrical Engineering and Slovak University of Technology in Bratislava, Ilkoviˇ c ova 3, 812 19 Bratislava, University Technology in Bratislava, Ilkoviˇ ccova 3, 812 19 Bratislava, e-mail: of {adrian.ilka, tomas.ludwig, ivan.ottinger, marian.tarnik, University of Technology in Bratislava, Ilkoviˇ ova 3, 812 19 Bratislava, e-mail: {adrian.ilka, tomas.ludwig, ivan.ottinger, marian.tarnik, e-mail: {adrian.ilka, marian.tarnik, eva.miklovicova, vojtech.vesely,ivan.ottinger, jan.murgas}@stuba.sk e-mail: {adrian.ilka, tomas.ludwig, tomas.ludwig, ivan.ottinger, marian.tarnik, eva.miklovicova, vojtech.vesely, jan.murgas}@stuba.sk eva.miklovicova, vojtech.vesely, jan.murgas}@stuba.sk eva.miklovicova, vojtech.vesely, jan.murgas}@stuba.sk Abstract: Novel original robust gain-scheduled controller design for individualized type 1 Abstract: Novel original robust gain-scheduled controller design for individualized type 1 Abstract: Novel original robust controller type diabetes mellitus (T1DM) subject model is presented in thisdesign paper.for Forindividualized controller design an11 Abstract: Novel (T1DM) original subject robust gain-scheduled gain-scheduled controller design for individualized typean diabetes mellitus model is presented in this paper. For controller design diabetes mellitus (T1DM) subject model is issubsystem presentedbased in this this paper. For controller controller designwith an LPV model is created from subject insulin-glucose on paper. Bergman’s minimal model diabetes mellitus (T1DM) model presented in For design an LPV model is created from insulin-glucose subsystem based on Bergman’s minimal model with LPV model is created from insulin-glucose subsystem based on Bergman’s minimal model with subcutaneous insulin absorption and absorption of digested carbohydrates. Identification was LPV model is created from insulin-glucose subsystem based on Bergman’s minimal model with subcutaneous insulin absorption and absorption of digested carbohydrates. Identification was subcutaneous insulin and absorption carbohydrates. Identification was performed on pharmacokinetics and pharmacodynamics characteristics of administered insulin subcutaneous insulin absorption absorptionand andpharmacodynamics absorption of of digested digested carbohydrates. Identification was performed on pharmacokinetics characteristics of administered insulin performed on pharmacokinetics and pharmacodynamics characteristics of administered insulin and data collected from continuous glucose monitoring (CGM) system. The controller design performed on pharmacokinetics and glucose pharmacodynamics characteristics ofThe administered insulin and data collected from continuous monitoring (CGM) system. controller design and from monitoring system. controller design approach guarantees the continuous closed-loopglucose stability and cost (CGM) for all scheduled parameter changes. and data data collected collected from continuous glucose monitoring (CGM) system. The The controller design approach guarantees the closed-loop stability and cost for all scheduled parameter changes. approach guarantees the closed-loop stability and cost for all scheduled parameter changes. Simulation results show the benefits of the proposed approach. approach guarantees the closed-loop stability and cost for all scheduled parameter changes. Simulation results the benefits of the proposed approach. Simulation results show show benefits the approach. Simulation show the the benefitsofof ofAutomatic the proposed proposed approach. © 2015, IFACresults (International Federation Control) Hosting by Elsevier Ltd. All rights reserved. Keywords: LPV system, Robust controller, Gain-scheduling, Output feedback, Quadratic Keywords: LPV system, Robust controller, Gain-scheduling, Output feedback, Quadratic Keywords: LPV controller, stability, Type diabetesRobust mellitus model Gain-scheduling, Keywords:Type LPV11 system, system, Robust controller, Gain-scheduling, Output Output feedback, feedback, Quadratic Quadratic stability, diabetes mellitus model stability, stability, Type Type 11 diabetes diabetes mellitus mellitus model model 1. INTRODUCTION controller design for nonlinear systems is nowadays a very 1. INTRODUCTION INTRODUCTION controller design for nonlinear systems is nowadays a very 1. controller design nonlinear systems is determinative andfor of research. 1. INTRODUCTION controller design forimportant nonlinearfield systems is nowadays nowadays a a very very determinative and important field of research. and important field of research. Computer modeling of type 1 diabetes mellitus (T1DM) determinative determinative and important field of research. Computer modeling modeling of of type type 11 diabetes diabetes mellitus mellitus (T1DM) (T1DM) Gain-scheduling is one of the most common used controller Computer is one of the most common used controller has attracted considerable attention in mellitus the past(T1DM) decade. Gain-scheduling Computer modeling of typeattention 1 diabetes Gain-scheduling is the common used controller approaches for of nonlinear a wide has attracted considerable in the past decade. design is one one of the most mostsystems commonand usedhas controller has attracted in the past design approaches for nonlinear systems and has a wide Patients with considerable T1DM sufferattention from high levels of decade. glucose Gain-scheduling has attracted considerable attention in the past decade. design approaches for nonlinear systems and has a wide range of use in industrial applications. Many of the Patients with T1DM suffer from high levels of glucose approaches for nonlinear systems and of has a early wide Patients with T1DM suffer from high levels of glucose range of use in industrial applications. Many the early concentration due to defective insulin secretion. The lack design Patients with T1DM suffer from high levels of glucose range of ofwere use in in industrialwith applications. Many of the theetearly early articles associated flight control (Adams al., concentration due to defective insulin secretion. The lack range use industrial applications. Many of concentration due to defective insulin secretion. The lack articles were associated with flight control (Adams et al., of insulin is preventing glucose uptake and utilisation by concentration due to defective insulin secretion. The lack articles were associated with flight control (Adams et al., 1992) and aerospace (Hyde and Glover, 1993). Then, gradof insulin is preventing glucose uptake and utilisation by articles were associated with flight control (Adams et al., of insulin is preventing glucose uptake and utilisation by 1992) and aerospace (Hyde and Glover, 1993). Then, gradcells. Long-term high glucose concentration results in sevof insulin is preventing glucose uptake and results utilisation by ually, 1992) and aerospace (Hyde and Glover, 1993). Then, gradthis approach has been used almost everywhere in cells. Long-term high glucose concentration in sev1992) and aerospace (Hyde and Glover, 1993). Then, gradcells. Long-term high concentration results in this approach has been used almost everywhere in eral health complications. The most common intensified cells.health Long-term high glucose glucose concentration results in sevsev- ually, ually, this approach has been used almost everywhere in control engineering, which was greatly helped with the eral complications. The most common intensified ually, this approach has been used almost everywhere in eral health complications. most control engineering, which was greatly helped with the insulin therapy nowadays isThe based oncommon manual intensified exogenous eral health complications. The most common intensified control engineering, engineering, which wasLinear greatlyparameter-varying helped with with the the introduction of LPV systems. insulin therapy nowadays is based on manual exogenous control which was greatly helped insulin therapy nowadays is based on manual exogenous introduction of LPV systems. Linear parameter-varying dosing to either keep the level of basal insulin or insulin dosing therapytonowadays is based on manual of Linear parameter-varying systems are time-varying plants whose space matrices insulin either keep keep the level level of basal basal exogenous insulin or or introduction introduction of LPV LPV systems. systems. Linearstate parameter-varying insulin dosing to insulin are time-varying plants whose state space matrices to suppress glycemic excursions after meal. The patient insulin dosing to either either keep the theafter levelaa of of basal insulin or systems systems are time-varying plants whose state space matrices are fixed functions of some vector of varying parameters to suppress glycemic excursions meal. The patient systems are time-varying plants whose state space matrices to suppress glycemic aa meal. The patient are fixed functions of some vector of varying parameters needs to take severalexcursions fingerstickafter blood glucose measureto suppress glycemic excursions after meal. The patient are fixed functions of some vector of varying parameters θ(t). They were introduced first by Jeff S. Shamma in needs to take several fingerstick blood glucose measureare fixed functions of some vector of Jeff varying parameters needs to take several fingerstick blood glucose measureθ(t). They were introduced first by S. Shamma in ments a day and make decisions on insulin doses. A closedneeds atoday take several fingerstickonblood glucose measureθ(t). They They weregain-scheduling. introduced first firstToday by Jeff Jeff S.LPV Shamma in 1988 to model the (Linear ments and make decisions insulin doses. A closedθ(t). were introduced by S. Shamma in ments a day and make decisions on insulin doses. A closed1988 to model gain-scheduling. Today the LPV (Linear loop blood glucose control would dramatically improve mentsblood a day glucose and make decisions on insulin doses. Aimprove closed- 1988 to model gain-scheduling. Today the LPV (Linear Parameter-Varying) paradigm has become a standard forloop control would dramatically 1988 to model gain-scheduling. Today the LPV (Linear loop blood glucose control would improve paradigm has become a standard forthe life of T1DM subjects. Despite the fast development looplife blood glucosesubjects. control Despite would dramatically dramatically improve Parameter-Varying) Parameter-Varying) paradigm become standard formalism in systems and controlshas with lot of aaresearches the of T1DM T1DM the fast fast development development Parameter-Varying) paradigm has become standard and forthe life of Despite the malism in systems and controls with lot of researches and of insulin pumps subjects. and continuous glucose measurement the life of T1DM subjects. Despite the fast development malism in systems and controls with lot of researches and articles devoted to analysis, controller design and system of insulin pumps and continuous glucose measurement malism in systems and controls with lotdesign of researches and of insulin pumps and continuous glucose measurement articles devoted to analysis, controller and system systems, a fully autonomous control of glycemia has not of insulina pumps and continuous measurement articles devoted devoted to analysis, analysis, controller design and system system of these modelscontroller (Shamma,design 2012).and systems, fully autonomous autonomous controlglucose of glycemia glycemia has not not identification to systems, aa fully control of identification of these models (Shamma, 2012). been introduced in a commercially available devicehas yet.not articles systems, fully autonomous control of glycemia has identification of these models (Shamma, 2012). been introduced in a commercially available device yet. identification of these models (Shamma, 2012). been introduced in commercially available device yet. The main motivation of our paper were our previous rebeen robust introduced in aa theory commercially device main motivation of our paper were our previous reThe control is well available established for yet. linear The The motivation our previous results main in gain-scheduling y and were Ilka, our 2013), (Ilka and The robust control theory is well established for linear The main motivation of of(Vesel´ our paper paper were our previous reThe robust control theory is well established for linear sults in gain-scheduling (Vesel´ y and Ilka, 2013), (Ilka and systems but almost all real processes are more or less The robust control theory is well established for linear sults in gain-scheduling (Vesel´ y and Ilka, 2013), (Ilka and Vesel´ y , 2014a), (Vesel´ y and Ilka, 2014) and the results from systems but almost all real processes are more or less sults in gain-scheduling (Vesel´ y and Ilka, 2013), (Ilka and systems all processes are or y , 2014a), (Vesel´ y and Ilka, 2014) and the results from nonlinear. If almost the plant operating region ismore small, one Vesel´ systems but butIf almost all real real processes are is more or less less Vesel´ y 2014a), y and Ilka, 2014) and the results from (T´arn´ et al., arn´ ık et al., 2014) nonlinear. the plant plant operating region small, one T1DM Vesel´ y,, research 2014a), (Vesel´ (Vesel´ yık and Ilka,2013), 2014)(T´ and the results from nonlinear. the operating region is small, one T1DM research (T´ a rn´ ık et al., 2013), (T´ a rn´ ık et al., 2014) can use theIf robust control approaches to design a linear nonlinear. If the plant operating region is small, one T1DM research (T´ a rn´ ık et al., 2013), (T´ a rn´ ık et al., 2014) and (Ottinger et al., 2015). In this paper a novel robust can use the robust control approaches to design a linear T1DM researchet(T´ arn´2015). ık et al., 2013), (T´arn´ ıknovel et al.,robust 2014) can use the robust control approaches to design a linear and (Ottinger al., In this paper a robust controller where the nonlinearities are treated as can usecontroller the robustwhere control to design a linear (Ottinger et In paper novel robust discrete gain-scheduling controller Bergman’s robust theapproaches nonlinearities are treated treated as and and (Ottinger et al., al., 2015). 2015). In this this design paper aafor novel robust robust controller where the nonlinearities are as discrete gain-scheduling controller design for Bergman’s model uncertainties. However, for real nonlinear processes, robust controller where the nonlinearities are treated as discrete gain-scheduling controller design for Bergman’s minimal model of glucose-insulin dynamics coupled with model uncertainties. However, for real nonlinear processes, discrete gain-scheduling controller design for Bergman’s model real nonlinear processes, model of glucose-insulin dynamics coupled with where the operating However, region is for large, the above mentioned model uncertainties. uncertainties. However, for realthe nonlinear processes, minimal minimal model of glucose-insulin dynamics coupled with insulin and carbohydrates absorption subsystems is prowhere the operating region is large, above mentioned model of glucose-insulin dynamics coupled with where the operating region is large, the mentioned insulin and carbohydrates absorption subsystems is controller synthesis is inapplicable. Forabove this reason the minimal where the operating region is large, the above mentioned insulin and and carbohydrates carbohydrates absorption absorption subsystems subsystems is is proproposed. controller synthesis is inapplicable. For this reason the insulin procontroller synthesis is For reason the controller is inapplicable. inapplicable. For this this reasonGrant the posed. posed. work synthesis has been supported by the Slovak Scientific The m×n posed. The work has been supported by the Slovak Scientific Grant Our notations are standard, D ∈ Rm×n denotes the set of Agency VEGA, 1/1241/12 andSlovak Grant Scientific No. 1/2256/12. The has been supported by Grant Our notations are D denotes the set m×n The work work has Grant been No. supported by the the Grant Agency VEGA, Grant No. 1/1241/12 andSlovak Grant Scientific No. 1/2256/12. m×n denotes Our m×n notations are standard, standard, D∈ ∈R Ridentity the setZof of real matrices. Im is an m×m matrix and m The paper is one of the outcomes of the research work for the Our notations are standard, D ∈ R denotes the set of Agency VEGA, Grant No. 1/1241/12 and Grant No. 1/2256/12. real m×n matrices. IIm is an m×m identity matrix and Z Agency VEGA, Grant No.outcomes 1/1241/12 1/2256/12. m The paper is one of the of and the Grant researchNo. work for the real m×n matrices. is an m×m identity matrix and Z denotes a zero matrix. If the size can be determined from m m project entitled ”Research center for severe diseases and related The paper is one of the outcomes of the research work for the real m×n matrices. I is an m×m identity matrix and Z m If the size can be determined from m The paper is one of the outcomes the research the denotes aa zero matrix. project entitled ”Research center forof severe diseaseswork and for related denotes zero matrix. If the size can be determined from the context, we will omit the subscript. P > 0 (P ≥ 0) is complications”, ”ITMS: 26240120038”. ”This project is being coproject entitled ”Research center for severe diseases and related denotes a zero matrix. If the size can be determined from project entitled ”Research center for severe diseases and related the context, we will omit the subscript. P > 0 (P ≥ 0) is complications”, ”ITMS: 26240120038”. ”This project is being cofinanced by the European Union. We support research activities in the context, we will omit the subscript. P > 0 (P ≥ 0) complications”, ”ITMS: 26240120038”. ”This project is being coa real symmetric, positive definite (semidefinite) matrix. the context, we will omit the subscript. P > 0 (Pmatrix. ≥ 0) is is complications”, ”This project being cofinanced by the ”ITMS: European26240120038”. Union. We support research isactivities in aa real symmetric, positive definite (semidefinite) Slovakia”. financed by the European Union. We support research activities in real symmetric, symmetric, positive positive definite definite (semidefinite) (semidefinite) matrix. matrix. financed a real Slovakia”.by the European Union. We support research activities in
Slovakia”. Slovakia”. Copyright 82 Hosting by Elsevier Ltd. All rights reserved. 2405-8963 ©©2015 2015,IFAC IFAC (International Federation of Automatic Control) Copyright © 2015 IFAC 82 Copyright © 82 Peer review underIFAC responsibility of International Federation of Automatic Copyright © 2015 2015 IFAC 82 Control. 10.1016/j.ifacol.2015.09.438
IFAC ROCOND 2015 July 8-11, 2015. Bratislava, Slovak Republic
Adrian Ilka et al. / IFAC-PapersOnLine 48-14 (2015) 082–087
83
1 1 (2c) S2 (t) TI VI where parameter TI [min] is a time constant of the subsystem, kI [1/min] is a decay rate of insulin in plasma and parameter VI [dl/kg] represents a insulin distribution volume per kilogram of body weight. Input v(t) [µU/kg/min] is insulin subcutaneous infusion rate, S1 (t) and S2 (t) [µU/kg] represent the amount of insulin in compartments of the subsystem.
Organisation of the paper is following. Section 2 includes problem formulation and some preliminaries are given. In Section 3 sufficient stability conditions in the form of BMI and/or LMI are given for the design of a robust discrete gain-scheduled controller. In Section 4 the obtained results are illustrated on the T1DM model.
˙ = −kI I(t) + I(t)
2. PROBLEM FORMULATION AND PRELIMINARIES
Our aim was to adjust the parameters of the proposed model so that the output of the model fits the continuous glucose monitoring (CGM) data of a particular T1DM subject. For identification of specific model parameters we used pharmacokinetics (PK) and pharmacodynamics (PD) measurements (as published in (Novo-Nordisk, 2002; Mudaliar et al., 1999)) of the particular insulin prescribed to the patient. The information about ingested carbohydrates was also recorded during data acquisition.
Third subsystem describes the glucose absorption from gastrointestinal tract, i.e. output of the subsystem is the signal Ra(t) [mg/kg/min]. The subsystem is described as follows 1 1 ˙ D(t) = − D(t) + AG d(t) (3a) T T D D 1 1 ˙ Ra(t) + D(t) (3b) Ra(t) =− TD TD where parameter TD [min] is a time constant and AG [dimensionless] is a friction of ingested carbohydrates which are effectively absorbed. Input d(t) [mg/kg/min] is the rate of carbohydrate ingestion at meal time, i.e. signal d(t) is an impulse with a width of one sampling period while the impulse area corresponds to the amount of ingested carbohydrates.
2.1 T1DM model
2.2 Identification of model parameters
Bergman’s minimal model consists of two differential equations in the form ˙ X(t) = −p2 X(t) + p2 SI (I(t) − Ib ) (1a) 1 ˙ G(t) = − SG + X(t) G(t) + SG Gb + Ra(t) (1b) VG where SG [1/min] is the rate constant which gives the rate of change of glucose caused by deviation from the basal glucose concentration Gb [mg/dl], parameter SI [ml/µU/min] is known as the insulin sensitivity index and p2 [1/min] is a rate constant. Parameter VG [dl/kg] represents the glucose distribution volume per kilogram of body weight BW [kg]. G(t) [mg/dl] is the blood glucose concentration and signal X(t) [1/min] represents the insulin in remote compartment. Values Ib [µU/ml] and Gb [mg/dl] are the basal insulin concentration and the basal glucose concentration respectively. In a basal steady state we have X(0) = 0 and G(0) = Gb .
For identification of model parameters we used data collected from a male T1DM subject aged 14, with BW = 64.6 [kg] and using fast-acting insulin NovoRapid (insulin Aspart) from an insulin pump.
In this section we briefly describe the mathematical model of a T1DM subject, which was based on Bergman’s minimal model of insulin-glucose interaction (Bergman et al., 1979). Later in this work the model will be used as a base for controller design and as a patient simulator for verification of the controller.
Insulin absorption subsystem: The first step in model identification was identifying of insulin absorption subsystem based on pharmacokinetics data of the used insulin. PK data from (Mudaliar et al., 1999) were used. An average basal insulin infusion rate vb [µU/kg/min] of the subject during the day is known since data from insulin pump are available. Signal v(t) is a sum of bolus part vB (t) and basal part vb . The aim is to identify the vector of unknown parameters Θ1 = [TI kI VI ] so that the error between simulated insulin concentration I(t) and PK data is minimized. In basal (steady) state for a given vb we get the basal insulin concentration Ib as the output and I(t) response after a bolus administration. We used the nonlinear least-squares optimization to identify the vector Θ1 .
Inputs of the model (1) are plasma insulin concentration I(t) [µU/ml] and glucose rate of appearance Ra(t) [mg/kg/min]. Signal Ra(t) can have in general two sources – the absorption of glucose from gastro-intestinal tract (modeled as a subsystem) and direct intravenous glucose administration.
Insulin sensitivity index and insulin action time: In the next step we identified the parameter related to insulin sensitivity SI and insulin action time p2 . These parameters determine dynamics of remote insulin signal X(t).
Insulin absorption is modeled as a separate subsystem where the output is insulin concentration I(t) (Herrero et al., 2013; Hovorka et al., 2004). The subsystem has the form 1 S1 (t) + v(t) (2a) S˙ 1 (t) = − T I 1 1 S2 (t) + S1 (t) (2b) S˙ 2 (t) = − TI TI
The measuring principle of pharmacodynamics is to maintain glycemia at basal concentration after bolus administration by intravenous glucose infusion. This glucose infusion corresponds to the signal Ra(t) in the equation (1b). If the equation (1b) is written in the form 1 ˙ G(t) = −SG (G(t) − Gb ) + Ra(t) − VG X(t)G(t) (4) VG 83
IFAC ROCOND 2015 84 July 8-11, 2015. Bratislava, Slovak Republic
Adrian Ilka et al. / IFAC-PapersOnLine 48-14 (2015) 082–087
furthermore
Table 1. T1DM identified model parameters TI 44.55
kI
VI
0.1645
138.8
SI 0.00159
p2 0.0106
SG 0.032
p1 (θ) = p1 0 +
TD
Gb
AG
1.467
8.5
0.95
p1 i θi , p3 (θ) = p3 0 +
i=1 p
33.474
a(θ) = a0 +
Table 2. Other fixed parameters of the model VG
p
ai θi , b(θ) = b0 +
i=1
p
p3 i θ i ,
i=1 p
bi θ i
i=1
The coefficient a(θ) is used to cover the nonlinear part of (1b) X(t)G(t) → x1 x2 in the following way (6) a(θ) = G(t) ⇒ a0 + a1 θ1 = x1 ⇒ a(θ)x2 = x1 x2 y−a0 where θ1 (t) = a1 . The coefficients a0 and a1 were calculated so as to maintain the scheduling parameter θ1 in the range −1, 1 min(y) + max(y) min(y) − max(y) a0 = , a1 = (7) 2 2 Note, the coefficients ai , i = 2, 3, 4, 5 are equal to zero. Similarly, the coefficient b(θ) is calculated in the following way b(θ)y = Ra (t) ⇒ b(θ) = b0 + b2 θ2 + b3 θ3 (8) where coefficients b0 and b2 are calculated so as to maintain the scheduling parameter θ2 in the range −1, 1 min(Ra /y) + max(Ra /y) b0 = 2 min(Ra /y) − max(Ra /y) b2 = 2 Ra y − b0 θ2 = b1 Furthermore b3 = 5 % of avarage Ra (t) (uncertainty) and bi = 0, i = 1, 4, 5 as well as θ3 (t) ∈ −1, 1 is unknown but constant parameter describing uncertainty.
˙ and we measure the PK, i.e. G(t) = 0, then G(t) ≈ Gb ∀t. It is obvious that parameter SG has minor influence in ˙ order to achieve G(t) = 0, so we assume SG = 0 during this step of parameter identification. The aim is to identify vector of unknown parameters Θ2 = [SI p2 ] so that the error between G(t) and Gb is minimized. Signal Ra(t) is given by PD data. Finalizing the model: At last, remaining parameters SG and TD are identified based on CGM data. The data containing both basal and bolus insulin dosing together with the amount of ingested carbohydrates were used as inputs to the model. Now we are identifying a vector of unknown parameters Θ3 = [SG TD ] so that the error between measured CGM data and the simulator output is minimized. Again, nonlinear least–square optimization was used. All identified parameters are reported in table 1. For the extended description of the identification process, please refer to our preliminary work (Ottinger et al., 2015)
For parameters p1 and p3 we also considered an uncertainty (±5%) (9) p1 (θ) = p1 0 + p1 4 θ4 , p3 (θ) = p3 0 + p3 5 θ5 where p1 0 = p1 , p1 4 = 5 % of p1 , p3 0 = p3 , p3 5 = 5 % of p3 , p1 i = 0, i = 1, 2, 3, 5, p3 i = 0, i = 1, 2, 3, 4 and θ4 , θ5 ∈ −1, 1 are unknown but constant parameters.
3. LPV-BASED ROBUST GAIN-SCHEDULED CONTROLLER DESIGN In this section a new LPV model is presented on the base of the nonlinear Bergman’s minimal model, which is then used to design a robust discrete LPV-based gain-scheduled controller for T1DM.
For the robust discrete LPV-based gain-scheduling controller design the model (5) is transformed to discrete time-space and W (θ) is neglected, because has no effect on stability.
3.1 LPV model of T1DM The Bergman’s model (1) with the insulin absorbation model (2) can be transfomed to the following LPV model with substitutions x1 (t) = G(t), x2 (t) = X(t), x3 (t) = S1 (t), x4 (t) = S2 (t), x5 (t) = I(t) and u(t) = v(t) x(t) ˙ = A(θ)x(t) + Bu(t) + W (θ) (5) y(t) = Cx(t)
3.2 Robust gain-scheduled controller design The output feedback gain-scheduled control law is considered for discrete-time PID (often denoted as PSD) controller in the form k e(i) u(k) = KP (θ|k)e(k) + KI (θ|k) (10)
where θ(t) ∈ Ω is a vector of scheduled parameters and 0 0 −p1 (θ) + b(θ) −a(θ) 0 0 −p2 0 0 p3 (θ) 1 0 0 0 0 − A(θ) = Ti 1 1 0 0 0 Ti − Ti 0 0 0 Ti1Vi −ki 0 x1 p1 (θ)Gb 0 x2 −p3 (θ)Ib I W (θ) = , B = 1 , x = x3 x 0 0 4 0 0 x5 C = [1 0 0 0 0]
i=0
+ KD (θ|k)(e(k) − e(k − 1)) where e(k) = y(k)−w(k) is control error, w(k) is reference value and gain matrices KP (·), KI (·), KD (·) are controller parameter matrices 1 (indexes P, I, D means proportional, sum (integral) and first difference (derivative), respectively) in the form p KP (θ|k) = KP 0 + i=1 KP i θi (k) p KI (θ|k) = KI 0 + i=1 KI i θi (k) p KD (θ|k) = KD 0 + i=1 KD i θi (k) 1
84
For SISO systems they are scalars.
IFAC ROCOND 2015 July 8-11, 2015. Bratislava, Slovak Republic
Adrian Ilka et al. / IFAC-PapersOnLine 48-14 (2015) 082–087
Note that the number of controller gain matrices is only 2 (for θ1 and θ2 ), the rest 3 (uncertainty) is equal to zero. Because the reference signal w(k) does not influence the closed-loop stability, we assume that it is equal to zero. For w(k) = 0, the control law (10) can be rewritten as u(k) = KP (θ|k)y(k) + KI (θ|k)
k
y(i)
and the value of closed-loop cost function (15) satisfies J ≤ J ∗ then J ∗ is said to be a guaranteed cost and u∗ guaranteed cost control law for system (5). Substituting the control law (12) to the quadratic cost function (15) we can obtain ˜T Q(θ|k) + C T F (θ|k)T R F (θ|k) C x ˜ (16) Jd (θ|k) = x Definition 2. (Apkarian et al., 1995) The linear closedloop system (14) for θ ∈ Ω is quadratically stable if and only if there exist a symmetric positive definite matrix P > 0 and for the first difference of Lyapunov function V (k) = xT P x along the trajectory of closed-loop system (14) holds (17) ∆V (θ|k) = Ac (θ|k)T P Ac (θ|k) + P < 0
(11)
i=0
+ KD (θ|k)(y(k) − y(k − 1)) State space description of PID controllers can be derived in the following way (Vesel´ y and Rosinov´ a, 2013). We can extend the system with two state variables z(k) = k−2 [z1T (k), z2T (k)]T where z1 (k) = i=0 y(i) and z2 (k) = k−1 y(i), then y(k − 1) = z (k) − z 2 1 (k). Substituting to i=0 (11) one obtains u(k) = (KP (θ|k) + KI (θ|k) + KD (θ|k)) y(k) (12) + KI (θ|k)z2 (k) − KD (θ|k)(z2 (k) − z1 (k))
From LQ theory we introduce the well known results. Lemma 1. (Kunceviˇc and Lyˇcak, 1977) Consider the closed-loop system (14). Closed-loop system (14) is affinely quadratically stable with guaranteed cost if and only if the following inequality holds (18) Be = min {∆V (θ|k) + Jd (θ|k)} ≤ 0
Control law (12) can be transformed to matrix form u(t) = F (θ|k)˜ y (k) (13)
u
where y˜ = [y(k), z1 (k), z2 (k)]T is the extended measurement output vector and KP (θ|k) + KI (θ|k) + KD (θ|k) KD (θ|k) F (θ|k)T = KI (θ|k) − KD (θ|k) Substituting the control law (13) to the discrete uncertain LPV system the closed-loop system is obtained in the form ˜(k) (14) x ˜(k + 1) = Ac (θ|k) x
for all θ(k) ∈ Ω. The main result of this section, the robust discrete gainscheduled controller design procedure, relies on the concept of multi-convexity, that is, convexity along each direction θi of the parameter space. The implications of multiconvexity for scalar quadratic functions are given in the next lemma (Gahinet et al., 1996). Lemma 2. Consider a scalar quadratic function of θ ∈ Rp . p p p p f (θ1 , . . . , θp ) = a0 + ai θi + bij θi θj + ci θi2
where x ˜(k) = [x(k), z1 (k), z2 (k)]T , Ac (θ|k) = Ar (θ|k) + Br (θ|k)F (θ|k)Cr (θ|k) and A(θ|k) 0 0 B(θ|k) 0 0 I , Br (θ|k) = 0 , Ar (θ|k) = 0 C 0 I C 0 0 Cr (θ|k) = 0 I 0 0 0 I
i=1
= p
i=1
Using Lemma 1 and 2 the following theorem is obtained Theorem 1. Closed-loop system (14) is quadratically stable with guaranteed cost if a positive defined P > 0 for all θ(k) ∈ Ω exists, matrices Qi , R, i = 1, 2, . . . p and gainscheduled controller matrices F (θ(k)) satisfy M (θ(k)) < 0; θ(k) ∈ Ω (19) (20) Mii ≥ 0; i = 1, 2, . . . p where p p p p Mi θ i + Mij θi θj + Mii θi2 (21) M (θ) = M0 +
To access the performance quality a quadratic cost function (Engwerda and Weeren, 2008) known from LQ theory is used in this paper, where weighting matrices depends on scheduling parameters (Ilka and Vesel´ y, 2014b). Using this approach we can affect performance quality in each operating point separately. The quadratic cost function is then in the form ∞ Jdf (θ) = (˜ x(k)T Q(θ|k) x ˜(k) + u(k)T R u(k) ∞
i=1 j>i
and assume that f (θ1 , . . . , θp ) is multi-convex, that is ∂ 2 f (θ) = 2ci ≥ 0 for i = 1, 2, . . . , p. Then f (θ) is negative ∂θi2 for all θ ∈ Ω if and only if it takes negative values at the corners of θ.
Remark 1. The controller’s filter of the derivative (differential) part can be included in the system model.
k=0
85
i=1
furthermore
(15)
Jd (θ|k)
M0 =
k=0
where Q(θ|k) = Q0 + i=1 Qi θi , Qi = QTi ≥ 0 where Q0 , Qi ∈ Rn×n , R ∈ Rm×m are symmetric positive definite (semidefinite) and definite matrices, respectively. The concept of guaranteed cost control is used in a standard way. Definition 1. Consider the system (5) with control algorithm (10). If there exists a control law u∗ and a positive scalar J ∗ such that the closed-loop system (14) is stable
Mi
=
Mij = Ac0 = Aci = Acij = 85
i=1 j>i
i=1
−P + Q0 C T F0T AcT0 −1 F0 C −R 0 −1 0 −P Ac 0 T T T Qi C Fi Aci Fi C 0 0 0 0 Aci 0 0 AcTij 0 0 AcTii , Mii = 0 0 0 0 0 0 Acii 0 0 Acij 0 0 Ar 0 + Br 0 F0 Cr Ar i + Br i F0 Cr + Br 0 Fi Cr Br i Fj Cr + Br j F iCr , Acii = Br i Fi Cr
IFAC ROCOND 2015 86 July 8-11, 2015. Bratislava, Slovak Republic
Adrian Ilka et al. / IFAC-PapersOnLine 48-14 (2015) 082–087
Kp (θ) = −101.1243 − 250.7895 θ1 + 239.9897 θ2 Ki (θ) = −5.1342 × 10−8 − 1.2379 × 10−7 θ1 −1.0598 × 10−5 θ2 Kd (θ) = −1012.1475 − 5998.7879 θ1 + 0.1235 θ2 For the illustration propose simulation experiment results are shown in Fig. 2. During manual administration of insulin by the T1DM subject, the measured glycemia has been higher than 10 mmol/l during 45% of the monitored time. In the case of automatic dosing controlled by the proposed gain-scheduled algorithm, the time when glycemia reached the level of 10 mmol/l or more was reduced to 9.9% of the simulation time.
Proof 1. Proof is based on Lemma 1 and 2. From (18) we can obtain M (θ(k)) = Ac (θ(k))T P Ac (θ(k)) + P (22) + Q + C T F (θ(k))T RF (θ(k))C < 0 Using Schur complement we obtain T T W11 W21 W31 T M (θ(k)) = W21 W22 W32 < 0 (23) W31 W32 W33 where W11 = −P + Q(θ(k)) W22 = −R−1 W21 = F (θ(k))C W32 = 0 W31 = Ac(θ(k)) W33 = −P −1 After we extend (23) to affine form we obtain (19) and (20) which proofs the Theorem 1.
5. CONCLUSION The paper addresses the problem of the robust discrete gain-scheduled controller design for Bergman’s minimal model of glucose-insulin dynamics coupled with insulin absorption subsystem and carbohydrates absorption subsystem. The obtained design procedure can be implemented easily to the standard LMI or BMI approaches. It can be used in low-cost micro-controllers, which is specifically beneficial for systems where we need to save the operation energy. The presented theory opens new possibilities for further research and study in this area.
Note that Theorem 1 in its presented form is in the form of BMI. One can use a free and open source BMI solver PENLAB or we can linearize the nonlinear part of (19) to use LMI solver (LMILAB or SEDUMI ). lin(−P −1 ) ≤ X −1 (P − X)X −1 − X −1 (24) where in each iteration pores X = P . Using this linearization, the element obtaining the nonlinear part (M0 ) become as follows AcT0 −P + Q0 C T F0T −1 M0 = F0 C −R 0 Ac0 0 X −1 (P − X)X −1 − X −1
REFERENCES Adams, R., Sparks, A.G., and Banda, S. (1992). A gain scheduled multivariable design for a manual flight control system. In First IEEE Conference on Control Applications, 584–589. Apkarian, P., Gahinet, P., and Becker, G. (1995). SelfScheduled H∞ Control of Linear Parameter-Varying Systems: A Design Example. Automatica, 31(9), 1251– 1261. Bergman, R., Ider, Y., Bowden, C., and Cobelli, C. (1979). Quantitative estimation of insulin sensitivity. Am. J. Physiol. Endocrinol. Metab, 236(6), E667–E677. Engwerda, J. and Weeren, A. (2008). A result on output feedback linear quadratic control. Automatica, 44(1), 265–271. Gahinet, P., Apkarian, P., and Chilali, M. (1996). Affine parameter-dependent Lyapunov functions and real parametric uncertainty. IEEE Transactions on Automatic Control, 41(3), 436–442. Herrero, P., Georgiou, P., Oliver, N., Reddy, M., Johnston, D., and Toumazou, C. (2013). A composite model of glucagon-glucose dynamics for in silico testing of bihormonal glucose controllers. Journal of Diabetes Science and Technology, 7(4), 941–951. Hovorka, R., Canonico, V., Chassin, L.J., Haueter, U., Massi-Benedetti, M., Federici, M.O., Pieber, T.R., Schaller, H.C., Schaupp, L., Vering, T., and Wilinska, M.E. (2004). Nonlinear model predictive control of glucose concentration in subjects with type 1 diabetes. Physiological Measurement, 25(4), 905. Hyde, R. and Glover, K. (1993). The application of scheduled h∞ controllers to a VSTOL aircraft. IEEE Transactions on Automatic Control, 38(7), 1021–1039. Ilka, A. and Vesel´ y, V. (2014a). Discrete gain-scheduled controller design: Variable weighting approach. In 15th International Carpathian Control Conference (ICCC), 2014, 186–191.
4. SIMULATION EXPERIMENTS In this section the proposed robust discrete gain-scheduled controller is verified using the individualized T1DM Bergman’s model served as a patient. For controller synthesis the LPV model described in section 3.1 with parameters presented in section 2.2, transformed to discrete time-space with sample time Ts = 5 min was used. The response of Bergman’s model with the discrete LPV model to an insulin bolus is shown in Fig 1. The disturbance has been considered in the form of mixed meal ingestion. The main objective was to keep the blood glucose concentration levels within normal glycemic range (3.8 – 10 mmol/l).
Fig. 1. Bergman’s model with the discrete LPV model The obtained model was extended for robust discrete gainscheduled PID controller design (14). Then using Theorem 1 with weighting matrices Q = qi I, q0 = 1 × 10−5 , q1 = 1 × 10−1 , q2 = 1 × 10−2 , q3 = q4 = q5 = 0 R = rI, r = 1 and ξL ≤ P (θ) ≤ ξU , ξU = 1 × 108 , ξL = 1 × 10−5 , Ts = 5 min with LMILAB one can obtain robust discrete gain-scheduled controller in the form (10) where 86
IFAC ROCOND 2015 July 8-11, 2015. Bratislava, Slovak Republic
Adrian Ilka et al. / IFAC-PapersOnLine 48-14 (2015) 082–087
87
Fig. 2. Simulation results for time period of 4 days T´arn´ık, M., Mikloviˇcov´a, E., Murgaˇs, J., Ottinger, I., and Ludwig, T. (2014). Model reference adaptive control of glucose in type 1 diabetics: A simulation study. In Proceedings of the 19th IFAC World Congress, 2014, Cape Town International Convention Centre, Cape Town, South Africa. T´arn´ık, M., Murgaˇs, J., Mikloviˇcov´a, E., and Farkas, v. (2013). Adaptive control of time-delayed systems with application for control of glucose concentration in type 1 diabetic patients. In 11th IFAC International Workshop on Adaptation and Learning in Control and Signal Processing - July 3-5 2013, Caen, France. Vesel´ y, V. and Ilka, A. (2013). Gain-scheduled PID controller design. Journal of Process Control, 23(8), 1141–1148. Vesel´ y, V. and Ilka, A. (2014). PID robust gain-scheduled controller design. In 2014 European Control Conference (ECC), 2756–2761. Vesel´ y, V. and Rosinov´a, D. (2013). Robust PID-PSD Controller Design: BMI Approach. Asian Journal of Control, 15(2), 469–478.
Ilka, A. and Vesel´ y, V. (2014b). Gain-Scheduled Controller Design: Variable Weighting Approach. Journal of Electrical Engineering, 65(2), 116–120. Kunceviˇc, V. and Lyˇcak, M. (1977). Control system design using Lyapunov function approach (in Russian). Nauka, Moskva. Mudaliar, S., Lindberg, F.A., Joyce, M., Beerdsen, P., Strange, P., Lin, A., and Henry, R.R. (1999). Insulin aspart (b28 asp-insulin): A fast-acting analog of human insulin. Diabetes Care, 22(9), 1501–1506. Novo-Nordisk (2002). NovoRapid (insulin aspart) - Product Monograph. Watermeadow Medical, Two Rivers House, Station Lane, Witney, Oxfordshire OX28 4BH, UK on behalf of Novo Nordisk A/S. Ottinger, I., Ludwig, T., Mikloviˇcov´ a, E., B´ atora, V., Murgaˇs, J., and T´ arn´ık, M. (2015). Individualized t1dm simulator for verification of adaptive controller. Process Control 2015 (Accepted). Shamma, J.S. (2012). Control of Linear Parameter Varying Systems with Applications, chapter An overview of LPV systems, 3–26. Springer.
87