Model Order Reduction of Glucose-Insulin Homeostasis Using Empirical Gramians and Balanced Truncation*

Proceedings of the 20th World Congress The International Federation of Congress Automatic Control Proceedings of the 20th World The International Federation of Congress Automatic Control Proceedings of the 20th9-14, World Toulouse, France, July 2017 The International Federation of Automatic Control Available online at www.sciencedirect.com Toulouse, France,Federation July 9-14, 2017 The International of Automatic Control Toulouse, France, July 9-14, 2017 Toulouse, France, July 9-14, 2017

ScienceDirect

IFAC PapersOnLine 50-1 (2017) 14735–14740 Model Order Reduction of Glucose-Insulin Model Order Reduction of Glucose-Insulin Model Order Reduction of Model Order Reduction of Glucose-Insulin Glucose-Insulin Homeostasis Using Empirical Gramians Homeostasis Using Empirical Gramians Homeostasis Using Empirical Gramians  Homeostasis Using Empirical Gramians and Balanced Truncation and Balanced Truncation and Balanced Balanced Truncation Truncation  and ∗ ∗∗

Christian Christian Tolks Tolks ∗ Christoph Christoph Ament Ament ∗∗ Christian Tolks ∗∗ Christoph Ament ∗∗ Christian Tolks Christoph Ament ∗∗ ∗ ∗ University of Augsburg, 86159 Augsburg, Germany (e-mail: University of Augsburg, 86159 Augsburg, Germany (e-mail: ∗ of Augsburg, 86159 Augsburg, Germany (e-mail: [email protected]). ∗ University [email protected]). University of Augsburg, 86159 Augsburg, Germany (e-mail: ∗∗ [email protected]). of ∗∗ University of Augsburg, Augsburg, 86159 86159 Augsburg, Augsburg, Germany Germany (e-mail: (e-mail: [email protected]). ∗∗ University of Augsburg, 86159 Augsburg, Germany (e-mail: [email protected]) ∗∗ University [email protected]) University of Augsburg, 86159 Augsburg, Germany [email protected]) (e-mail: [email protected]) Abstract: In this contribution, Abstract: In this contribution, a a novel novel method method for for model model order order reduction reduction is is presented. presented. Since Since Abstract: In this contribution, a novel method for model order reduction is presented. Since biological systems are often highly nonlinear or having a high number of parameters, the biological systems are often highlya nonlinear or having a highorder number of parameters, the concept concept Abstract: In this contribution, novel method for model reduction is presented. Since biological systems are often highly nonlinear or having aand highobservability number of parameters, the allows concept of Gramians is for analysis. aa of empirical empirical Gramians is considered considered for controllability controllability and observability analysis. This This allows biological systems are often highly nonlinear or having aand high numberby oftruncation. parameters, themotivate concept of empirical Gramians isof considered forwhich controllability observability analysis. This allows a balanced representation the system, offers model reduction We balanced representation of the system, which offers model reduction by truncation. We motivate of empirical Gramians isofconsidered forwhich controllability andreduction observability analysis. This allows a balanced the system, offers model bythat, truncation. We Gramian motivate to the states of parameters. For loadability to transfer transferrepresentation the idea idea of of balancing balancing states to to of balancing balancing parameters. For that, loadability Gramian balanced representation of the system, which offers model reduction by truncation. We motivate to transfer thefirst, idea of balancing states to of balancing parameters. For that, loadability Gramian is characterizing controllability Gramian for Second, is introduced introduced first, characterizing controllability Gramian for parameters. parameters. Second, parameters parameters to transfer the idea of balancing states to of balancing parameters. For that, loadability Gramian is introduced first, characterizing controllability Gramian for parameters. Second, parameters were treated as constant states, which allows observability-based identification. It is then possible were treated asfirst, constant states, which allows observability-based identification. It is then possible is introduced characterizing controllability Gramian for parameters. Second, parameters were treated as constant states, which allows observability-based identification. It is thenbalanced possible to obtain a balanced realization of the parameter-space of the system, which enables to obtain a balanced realization of the parameter-space of the system, which enables balanced were treated as constant states, which allows observability-based identification. It is then possible to obtain a This balanced realization of the parameter-space of the system, which enables balanced truncation. method is to of homeostasis. truncation. This method is applied applied to systems systems of glucose-insulin glucose-insulin homeostasis. to obtain a This balanced realization of the parameter-space of the system, which enables balanced truncation. method is applied to systems of glucose-insulin homeostasis. truncation. method isFederation applied to systems ofControl) glucose-insulin © 2017, IFACThis (International of Automatic Hosting byhomeostasis. Elsevier Ltd. All rights reserved. Keywords: Keywords: Glucose-Insulin Glucose-Insulin Homeostasis, Homeostasis, Empirical Empirical Gramians, Gramians, Balanced Balanced truncation, truncation, Keywords: Glucose-Insulin Homeostasis, Empirical Gramians, Balanced truncation, Controllability, Observability, Identifiability, Hankel singular values, Model order Controllability, Observability, Identifiability, Hankel singular values, Model order reduction reduction Keywords: Glucose-Insulin Homeostasis, Empirical Balanced truncation, Controllability, Observability, Identifiability, HankelGramians, singular values, Model order reduction Controllability, Observability, Identifiability, Hankel singular values, Model order reduction 1. ciple 1. INTRODUCTION INTRODUCTION ciple component component analysis analysis provides provides aa way way of of finding finding this this 1. INTRODUCTION ciple component analysis provides a way of findingthose this best low-dimensional approximation. In other words, best low-dimensional approximation. In other words, those 1. INTRODUCTION ciple component analysis provides aInway ofinput-output findingthose this best low-dimensional approximation. other words, states with a small influence on the system’s In general, biomedical processes can be described by orstateslow-dimensional with a small influence on the system’s input-output In general, biomedical processes can be described by or- best approximation. In other words, those states with a small influence on the system’s input-output can be In general, biomedical processes can be described by or- behavior dinary differential equations, characterizing their behavior can be neglected. neglected. dinary differential equations, characterizing their system system with a subset small influence on the system’s input-output In general, biomedical processes cantechniques be described bylead or- states behavior can be neglected. An optimal can dinary differential equations, characterizing their system behavior. Advances in experimental have An optimal subset can be be found found if if the the system system is is balbalbehavior. Advancesequations, in experimental techniques have lead behavior can be neglected. dinary differential characterizing their system An optimal subset can be found if the system is balanced. That means, controllability and observability in behavior. Advances in experimental techniques have lead to an increasing number of mathematical models with a anced. That means, controllability and observability in to an increasing number of mathematical models with a An optimal subset can be found ifand the system is balbehavior. Advances in experimental techniques have lead anced. That means, controllability observability in this representation are weighted equally and the states are to an increasing number of mathematical models with a raising level of accuracy, that usually have a high count of this representation are weighted equally and the states are raising level of accuracy, that usually have a high count of anced. That means, controllability and observability in to an increasing number ofnot mathematical models with ofa this representation are weighted equally and the states are decoupled. We suggest to transfer this concept of state raising level of accuracy, that usually have a high count states and parameters. But all of them are measurable decoupled. We suggest to transfer thisand concept of state states and parameters. But notusually all of them are measurable this representation are weighted equally the states are raising level of accuracy, that have a high count of decoupled. We suggest to transfer this concept of state order reduction to of parameters. For states and parameters. notFor all of them are or known priori. instance, it is reduction to the the reduction reduction of model model parameters. For or otherwise otherwise known a a But priori. For instance, it measurable is impossiimpossi- order decoupled. Wewesuggest toparameters transfer this concept of apply state states and parameters. But not all of them are measurable order reduction to the reduction of model parameters. For this purpose, consider as states and or otherwise known a priori. For instance, it is impossible to measure all relevant states, e.g. concentrations of this purpose, we consider parameters as states and apply ble to measure all relevant states, e.g. concentrations of order reduction to the reduction of model parameters. For or otherwise known a priori. For instance, it iscontroller impossithis purpose, we consider parameters as states and apply the same algorithm for parameter reduction. ble to measure all relevant states, e.g. concentrations of substances in a living subject for monitoring, the same algorithm for parameter reduction. substances in a living subject for monitoring, controller this purpose, we consider parameters as states and apply ble to measure relevant states, e.g. concentrations of the same algorithm for parameter reduction. and To obtain information about substances in a all living subject forAnyhow, monitoring, controller design state estimation. it always obtain information about controllability controllability and observobservdesign or or online online state estimation. Anyhow, it is is not not always To the same algorithm for parameter reduction. substances in a living subject for monitoring, controller To obtain information about controllability and observability of nonlinear systems, such as biomedical ones, design or online state estimation. Anyhow, it is not always necessary to all from ability of nonlinear systems, such as biomedical ones, necessary to identify identify all parameters parameters fromit experimental experimental To obtain information about controllability and observdesign or online state estimation. Anyhow, is not always ability of nonlinear systems, such as biomedical ones, analytical approaches are often not applicable. Here, we necessary to identify all parameters from experimental data. analytical approaches are often not applicable. Here, we data. ability of nonlinear systems, such as biomedical ones, necessary to identify allfield parameters from experimental analytical approaches are often not applicable. Here, we consider the use of empirical Gramians, which allows indata. Reduced models in the of systems biology, here esconsider the use of empirical Gramians, which allows inReduced models in the field of systems biology, here esanalytical approaches are often not applicable. Here, we data. consider the use of empirical Gramians, which allows insights in local system behavior in aa region around nominal Reduced models in the field of systems biology, here especially in the area of metabolic disorders, can play a sights in local system behavior in region around nominal pecially in the area of field metabolic disorders, canhere playes-a consider the use of empirical Gramians, which allows inReduced models in the of systems biology, in local system behavior in a region around nominal operation points. pecially in Models the areawith of metabolic disorders, play a sights great reduced of states and points. great role. role. Models with reduced number number of can states and sights in local system behavior inhas a region around nominal pecially in lower the area of metabolic disorders, can play a operation operation points. But first, an appropriate great role. Models with reduced of states and number parameters computational effort on simulating, parBut first, an appropriate model model has to to be be defined. defined. parameters lower computational effort on simulating, paroperation points. great role. Models with reduced number of states and But first, an appropriate model has to be defined. parameters lower computational effort on simulating, particularly if these models run on mobile devices. A reticularly if these models run on mobile devices. A reBut first, an appropriate model has to be defined. parameters effort on simulating, 2. ticularly if lower these computational models run on mobile devices. Aparreduced space also simplifies design of estimators 2. MODEL MODEL OF OF THE THE GLUCOSE-INSULIN GLUCOSE-INSULIN duced state state space also simplifies design of state state estimators ticularly if these models run on mobile devices. A re2. MODEL OF THE GLUCOSE-INSULIN HOMEOSTASIS duced state space design of state like filter for online and prediction HOMEOSTASIS 2. MODEL OF THE GLUCOSE-INSULIN like Kalman Kalman filteralso for simplifies online monitoring monitoring and estimators prediction duced state space also simplifies design or of state estimators HOMEOSTASIS like Kalman filter for online monitoring and prediction purposes (Eberle and Ament (2011)) the building of HOMEOSTASIS purposes (Eberle and Ament (2011)) or the building of like Kalman filter pumps for Ament online monitoring andloop prediction purposes (2011)) or the building of With automatic insulin acting in control With aa growing growing number number of of affected affected people people living living with with DiaDiaautomatic(Eberle insulin and pumps acting in a a closed closed loop control purposes (Eberle and Ament (2011)) or the building of With a growing number of affected people living with Diabetes (International Diabetes Federation (2015)), underautomatic insulin pumps acting in a closed loop control (Russell et al. (2016)). betes (International Diabetes Federation (2015)), underWith a growing number of affected people living with Dia(Russell et insulin al. (2016)). automatic pumps acting in a closed loop control betes (International Diabetes Federation (2015)), understanding mechanisms of this metabolic disorder through (Russell et al. (2016)). Model order reduction can be a method of choice to adstanding mechanisms of this metabolic disorder through betes (International Diabetes Federation (2015)), underModel order reduction can be a method of choice to ad(Russell et al. (2016)). The mechanisms of this homeostasis metabolic disorder through modeling the is Model order reduction can be method of choice to ad- standing dress problems. goal to a modeling the glucose-insulin glucose-insulin homeostasis is aa indispensindispensstanding mechanisms of this metabolic disorder through dress these these problems. The goalaa is is to derive derive a low-order low-order Model order reduction can be method of choice to admodeling the glucose-insulin homeostasis is a indispensable step, as it provides the possibility of simulation, dress these problems. The goal is to in low-order approximation of system such that step,the as glucose-insulin it provides thehomeostasis possibility isof asimulation, modeling indispensapproximation of the the original original system inderive such aa way, way, that able dress these problems. The goal is to derive low-order able step, as it model provides the controller possibility design of simulation, system analysis, driven or approximation the original system in such a way, of behavior the of reduced model does not analysis, driven design or state state able step, as it model provides the controller possibility of of simulation, the input-output input-output behavior of the the reduced model doesthat not system approximation of behavior the original system in such a way, that system analysis, model driven controller design or state estimation procedures, without the necessity taking a the input-output of the reduced model does not differ significantly from that of the original system. This estimation procedures, without the necessity of taking a system analysis, model driven controller design or state differ significantly from that of the original system. This the input-output behavior of subspace thethe reduced model does not estimation procedures, without the necessityto ofblack-box taking a high amount of measurements. In contrast differ significantly This from that of original system. can be done by utilizing a of the state space high amount of measurements. In contrast to black-box estimation procedures, without the necessity of taking a can be done by utilizing a subspace of the state space differ significantly from that of the original system. This high amount of measurements. In contrast to black-box modeling approaches like ARMAX models or neural nets, can be which done by subspace of the state space system, only contains important dynamics. Prinapproaches like ARMAX models ortoneural nets, high amount of measurements. In contrast black-box system, which onlyutilizing containsaathe the important dynamics. Prin- modeling can be done by utilizing subspace of the state space approaches ARMAX or neural nets, a physiological model approach is preferred in to system, which only contains the important dynamics. Prin- modeling physiological modellike approach is models preferred in order order to approaches like ARMAX models or neural nets, system, which contains theAutomotive importantEngineering dynamics.GmbH, Prin- amodeling  agather physiological model approach is preferred in order to information about fluxes and their transportation work was only supported by IAV  This gather information about fluxes and their transportation This work was supported by IAV Automotive Engineering GmbH, a physiological model approach is preferred in order to  gather information about fluxes and their transportation within the body. In Bergman et al. (1981), one of the Berlin, Germany This work was supported by IAV Automotive Engineering GmbH, within the body. Inabout Bergman etand al. their (1981), one of the Berlin, Germany  This work gather information fluxes transportation was supported by IAV Automotive Engineering GmbH, within the body. In Bergman et al. (1981), one of the Berlin, Germany within the body. In Bergman et al. (1981), one of the Berlin, Germany

Copyright © 2017 IFAC 15300 Copyright 2017 IFAC 15300 2405-8963 © 2017, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Copyright ©under 2017 responsibility IFAC 15300 Peer review of International Federation of Automatic Control. Copyright © 2017 IFAC 15300 10.1016/j.ifacol.2017.08.2576

Proceedings of the 20th IFAC World Congress 14736 Christian Tolks et al. / IFAC PapersOnLine 50-1 (2017) 14735–14740 Toulouse, France, July 9-14, 2017

first models of glucose-insulin homeostasis was described and introduced as minimal model. It has been extended and improved since then, see Bergman (2005) or Cobelli et al. (2009), and is widely excepted in clinical practice. The model depict the main physiological effects with three states. Two describe concentrations of glucose and insulin in plasma, denoted by G(t) and I(t), respectively. Remote insulin, as the third state, is denoted by X(t). These states were extended in Eberle and Ament (2012) by two dynamics describing first- and second-phase insulin response, denoted by v1 (t) and v2 (t), respectively. Model of the gastroenterological system is achieved by a second-order delay structure, denoted by qG1 (t) and qG2 (t), respectively. To allow applying subcutaneous insulin bolus, input uis (t) and a first-order delay system qI (t) were introduced. To gather information about tissue glucose concentrations, a new interstitial compartment, modeled as a first-order delay system with output Y (t), was added. The system is modeled as followed with nine states and 16 parameters in total, summarized in Table 1. 1 ˙ G(t) = − G(t) − kx X(t) [G(t) + Gb ] + . . . TG 1 1 uGV (t) qG2 (t) + VG BW ˙ = − 1 I(t) + v1 (t) + v2 (t) + qI (t) + uIV (t) I(t) TI 1 1 ˙ X(t) + I(t) X(t) =− TX TX 1 1 Y (t) + G(t) Y˙ (t) = − TY TY 1 kG1 v˙ 1  (t) = − v1  (t) − 2 G(t), T1 T1 k G1 G(t) v 1 = v1  + T1 1 kG2 G(t) v˙ 2 (t) = − v2 (t) + T2 T2 1 1 1 1 uGO (t) qG1 (t) + q˙G1 (t) = − TO1 TO1 VG BW 1 1 qG2 (t) + qG1 (t) q˙G2 (t) = − TO2 TO2 1 1 uIS (t) q˙I (t) = − qI (t) + TS TS

3. METHODS 3.1 Review of Gramians First introduced by Moore (1981), Gramian matrices are used for controllability and observability analysis of linear state space systems. In the nonlinear case, these matrices can be extended to empirical Gramians, which were first introduced by Lall et al. (1999) and improved by Hahn and Edgar (2002). A survey can be found in Baur et al. (2014). In Himpe and Ohlberger (2013) a software framework for empirical Gramians was presented. Since controllability and observability depends on the system’s trajectory, the empirical Gramians are computed by calculating the average over various simulations or experimental data following a special design of experiment. For controllability, input signals are deflected around nominal excitation and influence on the system’s states is considered. For observability, initial states are perturbed and the corresponding output trajectory is measured. Extending this concept to balancing of parameters provides a second way to further reduce state space. Analog to the controllability Gramian for states, the pseudo loadability Gramian for parameters is introduced (section 3.4). For the identifiability procedure, parameters were converted to constant states and observability Gramian can be applied. In stable linear systems, controllability and observability matrices WC and WO can be found by using linear inputto-state map (5) and state-to-output map (7), respectively, which represent the energy transfer from inputs to states and from states to outputs, respectively.

(1)

x(t) = eAt x0  ∞ EC = Tr x(t) xT (t) dt 0  ∞ T WC = eAt BBT eA t dt

x(t) = [G(t), I(t), X(t), Y (t), v1 (t), v2 (t), . . . T

(4) (5)

0

y(t) = Cx = CeAt x0  ∞ EO = yT (t) y(t) dt 0  ∞ T = xT0 eA t CT CeAt dt x0 0  

Eq. (1) builds a nonlinear state space system that can be summarized into state vector

qG1 (t), qG2 (t), qI (t)] .

(3)

(2)

(6)

(7)

WO

With A, B, and C denoting state, input and output matrices of linear state space system, respectively. Preserved energies are denoted by EC and EO . Given a nonlinear system equation of the form

The model has three inputs. Signals uGV (t) and uGO (t) are defined such that an intravenous glucose tolerance test (IVGTT) and an oral glucose tolerance test (OGTT) can be modeled. Input uIS (t) takes into account supplying of insulin subcutaneously. Note that all states are described as differences to their basal values Gb and Ib , to simplify model equations. Here, only uGO (t) is used to simulate ingestion of glucose. Three outputs can be measured: concentrations of plasma and interstitial glucose, G(t) and Y (t) respectively, as well as insulin concentration in plasma I(t).

˙ x(t) = f (x(t), u(t), p) , y(t) = g (x(t)) ,

(8)

with states x(t) ∈ Rn , inputs u(t) ∈ Rm , outputs y(t) ∈ Ro and parameters p(t) ∈ Rp empirical Gramian can be considered. Therefore one has first to define a set of input perturbations given by {Eu , Ru , Qu } and a set of state perturbations, given by {Ex , Rx , Qx }. This allows rotation and scaling of inputs and initial states, details can be found in Himpe and Ohlberger (2014):

15301

Proceedings of the 20th IFAC World Congress Toulouse, France, July 9-14, 2017 Christian Tolks et al. / IFAC PapersOnLine 50-1 (2017) 14735–14740

14737

Table 1. Compartments, maximum states and parameters of glucose-insulin model (values taken from Eberle and Ament (2012)), described by a set of nine first-order differential equations (only compartment 1 has nonlinear state relations). Maximal value

Compartment

Plasma glucose G(t) in mg dl−1

100

Parameters Name TG kX Gb VG BW TI Ib

Unit min ml/µU/min mg/dl dl/kg kg min µU/ml

Nominal values 1 32.4 507 × 10−6 85 1.88 (Ref. Dalla Man et al. (2007)) 70 3.33 9

Plasma insulin I(t) 50 in µU/ml Interstitial insulin X(t) TX min 47.78 10 in µU/ml Interstitial glucose Y (t) TY min 16 50 in mg dl−1 First-phase glucose T1 min 4 20 response v1 (t) in mg dl−1 min−1 kG1 dl/mg/µU/ml 1.08 2 Second-phase glucose T2 min 12 10 response v2 (t) in mg dl−1 min−1 kG2 dl/mg/µU/ml/min 0.107 Oral glucose flow qG1 (t) TO1 min 15 10 (1) in mg dl−1 min−1 Oral glucose flow qG2 (t) TO2 min 40 10 (2) in mg dl−1 min−1 Subcutaneous insulin TS min 20 30 flow qI (t) in mg dl−1 min−1 1 : Nominal values match with a non-pathological metabolism 2: k G1 is multiplied by factor 20 in case of oral excitation to model incretin effect, Eberle and Ament (2012)

  Eu = ei ∈ Rj ; ei  = 1; i = 1, . . . , m ,   Ru = Si ∈ Rj×j ; ST S = I; i = 1, . . . , s , Qu = {ci ∈ R; ci > 0; i = 1 . . . , q} , Ex = {fi ∈ Rn ; fi  = 1; i = 1, . . . , n} ,   Rx = Ti ∈ Rn×n ; TT T = I; i = 1, . . . , t , Qx = {di ∈ R; di > 0; i = 1, . . . , r} .

3.3 Empirical observability Gramian

(9)

3.2 Empirical controllability Gramian Empirical controllability Gramian WC can be obtained using the sets Eu , Ru and Qu , as well as input signal u(t) ¯ (t) and monitoring the around nominal input trajectory u ¯ (¯ corresponding state trajectory x u):  ∞ q s m 1  1 WC = Ψhij (t)dt, 2 qs c 0 (10) h=1 i=1 j=1 h  hij   hij T hij n×n ¯ x (t) − x ¯ ∈R Ψ (t) = x (t) − x .

Empirical observability Gramian characterizes how well an output can be driven by the states and is calculated by perturbing the initial states of the system and observing the output trajectory y(t) for each numerical experiment ¯ (¯ and comparing with nominal output y x):  ∞ r  t  1 T Ψuv (t)dt TTv , (13) WO = 2 v r t d 0 u u=1 v=1   uva ¯ (t)]T yuvb (t) − y ¯ (t) . (t) − y Ψuv ab (t) = [y The initial state vector x0 is perturbed as follows: ¯ + du Tv fa , xuva =x 0

describing a combination of δ(t)-impulses around nominal ¯ . Furthermore, the scalar constant du describes the value x scaling of amplitudes within the maximum range given by diagonal matrix T. For a system with discrete sampling time ∆T and time steps k, WO can be obtained by:

q  s  m 

k end 

1 1 Ψhij (k)∆T, 2 qs c h=1 i=1 j=1 h k=1  hij  T hij ¯ xhij (k) − x ¯ ∈ Rn×n . Ψ (k) = x (k) − x WC =

(12)

t r  

k

end  1 T Ψuv (k)∆T TTv , v 2 r t d u u=1 v=1 k=1  T  uvb uva ¯ ¯ (k) . y (k) − y (k) = [y (k) − y (k)] Ψuv ab

Inputs for a specific set of input vector u(t) is generated by configuration ¯ (t) + ch Si ej u(t). (11) uhij (t) = u For a system with sampling time ∆T and discrete time steps k, WC can be obtained by:

(14)

WO =

(15)

3.4 Loadability Gramian

Here it is introduced the pseudo Gramian WL , transferring the principle of controllability of states to that of parameters as illustrated in Fig. 1. Note the difference to the sensitivity Gramian (for details see Himpe and

15302

Proceedings of the 20th IFAC World Congress 14738 Christian Tolks et al. / IFAC PapersOnLine 50-1 (2017) 14735–14740 Toulouse, France, July 9-14, 2017

Note that Ψuv (t) and Ti now have dimension Rnˆ ׈n . Using  sampling time ∆T with discrete time steps k, WO is

Dynamic system Inputs  u1     um 

States  x1      xn 

 y1      yo 

Sensitivity Parameter  p1       pp   

Parameter not controllable

Observability and identifiability

Fig. 1. Decomposition of Gramians: Concept of controlla06.03.20 Seite 0 bility and observability of states and loadability and 15 identifiability of parameters. Ohlberger (2013) or Sun and Hahn (2006)) that handle parameters as additional inputs to a system and hence, give information about how parameter changes influence the states of a system. It is pondered a balanced realization of parameters, since they can also be scaled freely just as states. One can, e.g. by altering the units, improve observability of a system. But in return, controllability will deteriorate. Since parameters are not controllable, energy in the parameters is considered:   EL = Tr ppT (16) WL = p pT .

(17)

It is obtained the loadability Gramian matrix WL ∈ Rp×p similar to controllability Gramian WC for balancing purposes. 3.5 Empirical identifiability Gramian Empirical identifiability is about observing the system’s parameters and is an excepted method in biological (Streif et al. (2006)) or chemical (Sun and Hahn (2006)) processes. Identifiability can be done by augmenting the nominal state vector with parameters that are now assumed to be constant states over time. One gets the following ˆ (t): augmented system x       ˙ x x(t) f (x(t), u(t), p) ˙x ˆ (t) = ˆ 0 = 0 , (18) = , x p0 p˙ 0 ˆ (t) = g (ˆ y x(t)) , (19) ˆ 0 containing the ˆ (t) ∈ R , n ˆ = n + p and x with x initial values of prior state vector x0 and the parameter values. Now, observability of the extended system can be obtained. Applying Eq. (13) on the extended state  ˆ (t), augmented observability Gramian WO vector x can be calculated as followed:  ∞ r  t  1  WO = T Ψuv (t)dt TTv , (20) 2 v r t d 0 u u=1 v=1   uva ¯ (t)]T yuvb (t) − y ¯ (t) . (t) − y Ψuv ab (t) = [y n ˆ

k

(21)

 holds information about observability of Gramian WO states in the upper left sub matrix WO ∈ Rn×n , and information about observability of parameters in the lower right sub matrix WP ∈ Rp×p :   WO WOP  WO = (22) WP O WP

Controllability

Initial „loading“ of parameter

r  t 

end  1 T Ψuv (k)∆T TTv , 2 v r t d u u=1 v=1 k=1  T  uvb uv uva ¯ (k)] y (k) − y ¯ (k) . Ψab (k) = [y (k) − y

 = WO

Outputs

Gramian identifiability matrix WI can then be extracted via the Schur complement: (23) WI = WP − WP O WP−1 WOP . 3.6 Balancing transformation After obtaining empirical Gramians WC , WO , WL , and WI , balancing can be conducted directly from the state space representation, see e.g. Laub et al. (1987). Note that the empirical Gramians only give a scope over a limited region around the nominal system trajectory and not a global approximation. Although balancing breaks the physical meaning of original states and parameters, but gives linear transformation matrices allowing to switch between former and balanced realization. For balancing system (8), it must be controllable, observ˜C = able and stable. It is principle axis balanced, if W ˜ O = Σ. With Σ being a diagonal matrix containing the W Hankel singular values σi in descending order: (24) Σ = diag (σ1 , . . . , σn ) . This means, one searches a transformation matrix V, such that states, which are difficult to reach are equal to those, which are difficult to observe, (see Antoulas (2005), pp. 209–211), and the following holds: ˜ C = VWC VT , (25) W ˜ O = V−T WO V−1 , (26) W with V donating the transposed inverse of V. Finding V is done by computing the Cholesky factorization of WC followed by decomposition of the eigenvalues: (27) WC = UUT , −T

KΣ2 KT = UT WO U. (28) A principle axis transformation V can then be calculated as follows: (29) V = Σ1/2 KT U−1 , V−1 = UKΣ−1/2 . Using V, general model (8) can be transformed balanced representation by   ˜˙ (t) = Vf V−1 x ˜ (t), u(t), p , x   ˜ (t) , ˜ (t) = g V−1 x y ˜ = Vx denoting the balanced state vector. with x

(30) into (31) (32)

3.7 Model reduction Model reduction by truncating simply means cutting off those states from the balanced realization, which do not

15303

200

Table 2. RMSE of nominal output trajectories compared to state (upper panel) and parameter reduced models (bottom panel).

50 0 60

120

180

240

Simulation time (min)

Fig. 2. Simulation (∆T = 1 min) of OGTT with ingestion of 2600 mg/min glucose for 30 min. Concentrations of glucose in plasma G(t) and tissue Y (t), and plasma insulin concentration I(t), were observed.

9 0.0 10−5 10−4

8 0.0 0.0 0.0

7 11.1 9.7 3.8

6 6.1 4.0 18.5

5 6.2 7.0 22.8

4 4.2 5.7 21.2

Output G(t) Y (t) I(t)

14 10−12 10−12 10−11

13 10−12 10−12 10−11

12 5.5 5.3 5.9

11 6.0 5.8 9.8

10 6.4 6.2 8.9

9 3.5 3.4 8.8

Input

4. RESULTS

First, system (1) was taken to obtain controllability and observability Gramians. As nominal input trajectory, an OGTT was fulfilled as shown in Fig. 2 (OGTT is typical for daily system excitation). Inputs and initial states were scaled to a maximum of 10% of their maximal values shown in Table 1. For balancing, rotation matrix S from set Ru and matrix T from set Rx were restricted to {±I}. That means, only one state/ parameter is perturbed at once, independently from all others. The model was reduced from 9 states to the order of 4 with  = 0.019 as can be seen in Fig. 3 (upper panel). In Fig. 4, output signals are shown for the balanced model as well as for the state reduced models using same configuration 3 2 1 0

0

1

2

3

4

5

6

7

8

Number of balanced states

1000

9

10

500 0

0

2

4

6

8

10

12

14

16

Number of balanced parameter

 Fig. 3. Hankel singular values σi (σ = diag(Σ2 )) of balanced model (1) of state vector (upper panel) and of parameter vector (lower panel).

Plasma glucose (mg/dl) Interstitial glucose (mg/dl) Plasma insulin ( U/ml)

150 100 50 200 150 100 Balanced model Reduced 8th order model Reduced 7th order model Reduced 4th order model

50 200 100 0 0

60

120

180

240

Simulation time (min)

k=r+1

with r denoting the remaining order of the system, such   that P = I(r×r) O(r×n−r) ∈ Rr×n .

Ingested glucose (2600 mg/min)

200

have a strong influence on the system’s input-output behavior. These states are characterized by low Hankel singular values in Σ. Since they are ordered in decreasing singular values, a Galerkin projection matrix P can be used for truncation, see Lall et al. (1999). The balanced ˇ = PVx and reduced model, represented by state vector x is then given by   ˇ˙ (t) = PVf V−1 PT x ˇ (t), u(t), p , x (33)  −1 T  ˇ (t) = g V P x(t) . y (34) To decide which states to truncate threshold  is used: n  Σk,k , (35) ≤2

Hankel singular values

Output G(t) Y (t) I(t)

100

0

Hankel singular values

States

Ingested glucose Plasma glucose Interstitial glucose Plasma insulin

150

14739

Parameter

Glucose concentration (mg/dl) Insulin concentration (µU/ml)

Proceedings of the 20th IFAC World Congress Toulouse, France, July 9-14, 2017 Christian Tolks et al. / IFAC PapersOnLine 50-1 (2017) 14735–14740

Fig. 4. Output trajectories during simulation of an OGTT with (state) balanced and reduced models. Root mean squared error is shown in upper field of Table 2. as displayed in Fig. 2. It can be seen that balanced and 8th order model have an insignificant small error, due to the fact that one Hankel singular value is very small. Error in plasma insulin output (bottom panel in Fig. 4) is relative high. This could be explained by the fact, that the insulin subsystem was not triggered during observability and controllability investigation and therefore, this part of the model is not represented adequately. Next, all parameters in model (1) were transformed into constant states, independently from any possible a priori knowledge. Observability analysis was utilized scaling the parameters to a maximum of 10% of their nominal values, also listed in Table 1. It was shown that parameter TS is not identifiable, because its corresponding state qI (t) remained zero and insulin input was not triggered in this case. Also, basal insulin concentration Ib is not identifiable, since it is only an offset to the output signal I(t). Hence, 14 of 16 parameters are identifiable and were taken into account for balancing and model reduction. Hankel singular values are depicted in Fig. 3 (lower panel). Since two parameters are not observable, reduction begins with order 14 and is stable until 7 parameters are left. RMSE calculation between nominal output trajectories and the outputs of the reduced models are listed in Table 2 (bottom panel). It can be stated that the error is comparably low down to 9 parameters. As a last test, another input configuration is considered for nominal and state reduced model of order 4. The system is now stimulated with a combined glucose ingestion and

15304

Proceedings of the 20th IFAC World Congress 14740 Christian Tolks et al. / IFAC PapersOnLine 50-1 (2017) 14735–14740 Toulouse, France, July 9-14, 2017

Inputs

insulin bolus administration, as it is usually done by patients with Diabetes. Note that the parameter set used here, only correspond to non-diabetic case. In Fig. 5, it can be seen that the reduced model is able to track this new input setup, too. Nevertheless, further investigations have to be done, e.g. with a parameter set describing a diabetic case. Insulin bolus (180 U/ml/min) Ingested glucose (2600 mg/min)

Plasma glucose (mg/dl)

140 120 100 80

Interstitial glucose (mg/dl)

60 140 120 100 80

Insulin ( U/ml)

60 150 Nominal model 9th order

100

Reduced model 4th order

50 0 0

60

120

180

240

Simulation time (min)

Fig. 5. Output trajectories using a combined input of insulin bolus at t = 0 min and glucose ingestion at t = 10 min for nominal and reduced model of 4th order, respectively.

5. CONCLUSION The method of balancing both, state-subspace and parametersubspace enables an easy way for model order reduction by balanced truncation. We pointed out that state reduction to half of the original number of states is possible and parameter reduction is applicable to models of glucose-insulin homeostasis. The reduced model can then be a starting point for observer design of states and parameters, system identification or adaptation. However, to get best controllability and observability, an optimal combination of suitable input excitations (OGTT, IVGTT, or insulin bolus) and signal measurements has to be found. Also, other dynamical phenomena have to be investigated to validate methods. REFERENCES Antoulas, A.C. (2005). Approximation of large-scale dynamical systems. Advances in design and control. Society for Industrial and Applied Mathematics. Baur, U., Benner, P., and Feng, L. (2014). Model order reduction for linear and nonlinear systems: A systemtheoretic perspective. Archives of Computational Methods in Engineering, 21(4), 331–358. Bergman, R.N. (2005). Minimal model: Perspective from 2005. Hormone Research, 64(3), 8–15. Bergman, R.N., Phillips, L.S., and Cobelli, C. (1981). Physiologic evaluation of factors controlling glucose tolerance in man: Measurement of insulin sensitivity and

beta-cell glucose sensitivity from the response to intravenous glucose. Journal of Clinical Investigation, 68, 1456–1467. Cobelli, C., Dalla Man, C., Sparacino, G., Magni, L., de Nicolao, G., and Kovatchev, B.P. (2009). Diabetes: Models, signals, and control. IEEE Reviews in Biomedical Engineering, 2, 54–96. Dalla Man, C., Rizza, R.A., and Cobelli, C. (2007). Meal simulation model of the glucose-insulin system. IEEE Transactions on Biomedical Engineering, 54(10), 1740– 1749. Eberle, C. and Ament, C. (2011). The unscented kalman filter estimates the plasma insulin from glucose measurement. Biosystems : Journal of Biological and Information Processing Sciences, 103(1), 67–72. Eberle, C. and Ament, C. (2012). Real-time state estimation and long-term model adaptation: A two-sided approach toward personalized diagnosis of glucose and insulin levels. Journal of diabetes science and technology, 6(5), 1148–1158. Hahn, J. and Edgar, T.F. (2002). An improved method for nonlinear model reduction using balancing of empirical gramians. Computers & Chemical Engineering, 26(10), 1379–1397. Himpe, C. and Ohlberger, M. (2013). A unified software framework for empirical gramians. Journal of Mathematics, 2013(9), 1–6. Himpe, C. and Ohlberger, M. (2014). Cross-gramian-based combined state and parameter reduction for large-scale control systems. Mathematical Problems in Engineering, 2014(2), 1–13. International Diabetes Federation (2015). IDF Diabetes Atlas. Lall, S., Marsden, J.E., and Glavaski, S. (1999). Empirical model reduction of controlled nonlinear systems. In Proceedings of the 14th IFAC World Congress, 473–478. Pergamon, Oxford and England. Laub, A., Heath, M., Paige, C., and Ward, R. (1987). Computation of system balancing transformations and other applications of simultaneous diagonalization algorithms. IEEE Transactions on Automatic Control, 32(2), 115– 122. Moore, B.C. (1981). Principal component analysis in linear systems: Controllability, observability, and model reduction. IEEE Transactions on Automatic Control, 26(1), 17–32. Russell, S.J., Hillard, M.A., Balliro, C., Magyar, K.L., Selagamsetty, R., Sinha, M., Grennan, K., Mondesir, D., Ehklaspour, L., Zheng, H., Damiano, E.R., and ElKhatib, F.H. (2016). Day and night glycaemic control with a bionic pancreas versus conventional insulin pump therapy in preadolescent children with type 1 diabetes: A randomised crossover trial. The Lancet Diabetes & Endocrinology, 4(3), 233–243. Streif, S., Findeisen, R., and Bullinger, E. (2006). Relating cross gramians and sensitivity analysis in systems biology. In Proceedings of the 17th International Symposium on Mathematical Theory of Networks and Systems (MTNS 2006), 437–442. Sun, C. and Hahn, J. (2006). Parameter reduction for stable dynamical systems based on hankel singular values and sensitivity analysis. Chemical Engineering Science, 61(16), 5393–5403.

15305