- Email: [email protected]
Adaptive Model Predictive Control for Nonlinearity in Biomedical Applications
6th IFAC IFAC Conference Conference on on Nonlinear Nonlinear Model Model Predictive Predictive Control Control 6th Madison, WI, 19-22, 6th IFAC Conference on Nonlinear Model Predictive Control Available online at www.sciencedirect.com Madison, WI, USA, USA, August August 19-22, 2018 2018 6th IFAC Conference on Nonlinear Model Predictive Control Madison, WI, USA, August 19-22, 2018 Madison, WI, USA, August 19-22, 2018
ScienceDirect
IFAC PapersOnLine 51-20 (2018) 368–373
Adaptive Model Predictive Control for Adaptive Model Predictive Control for Adaptive Model Predictive Control for Adaptive Model Predictive Control for Nonlinearity in Biomedical Applications Nonlinearity in Biomedical Applications Nonlinearity in Biomedical Applications Nonlinearity in† Biomedical Applications † ‡
† Mudassir Rashid † Mert Sevil ‡ Iman Iman Hajizadeh Hajizadeh Mert Sevil ‡ † Mudassir Rashid † † Mert Sevil ‡‡ Iman Hajizadeh Rashid † Rachel Brandt ‡‡ Sediqeh Sediqeh Samadi Nicole Hobbs ‡ † Mudassir † Rachel Brandt Samadi Iman Hajizadeh Mudassir Rashid MertHobbs Sevil ‡ ‡ † Nicole †,‡,1 Rachel Brandt ‡ Sediqeh Samadi †,‡,1 † Nicole Hobbs ‡ Ali Cinar Ali Cinar Rachel Brandt Sediqeh Samadi Nicole Hobbs †,‡,1 Ali Cinar †,‡,1 Ali Cinar † † Department of Chemical and Biological Engineering, † Department of Chemical and Biological Engineering, of Chemical and Biological Illinois Institute of Technology, Chicago, IL, † Department Illinois Institute of Technology, Chicago,Engineering, IL, USA USA Department of Chemical and Biological Engineering, ‡ Institute Illinois ofof Chicago, IL, USA ‡ Department ofTechnology, Biomedical Engineering, Department Biomedical Engineering, Illinois ‡ Institute of Technology, Chicago, IL, USA BiomedicalChicago, Engineering, Illinois Institute IL, ‡ Department Illinois Institute of ofofTechnology, Technology, IL, USA USA BiomedicalChicago, Engineering, IllinoisDepartment Institute ofofTechnology, Chicago, IL, USA Illinois Institute of Technology, Chicago, IL, USA Abstract: Abstract: An An adaptive adaptive model model predictive predictive control control (MPC) (MPC) formulation formulation is is proposed proposed in in this this work work Abstract: An adaptive model predictive control (MPC) formulation is proposed in work for optimal insulin dosing decisions in artificial pancreas (AP) systems. To this end, a recursive for optimal An insulin dosingmodel decisions in artificial pancreas (AP) systems.isTo this end,ina this recursive Abstract: adaptive predictive control (MPC) formulation proposed this work for optimal insulin dosing decisions in approach artificial pancreas systems. the To this end, a dynamics recursive subspace-based system identification approach is used used to to(AP) characterize the transient dynamics subspace-based system identification is characterize transient for optimal insulin dosing decisions in artificial pancreas (AP) systems. To this end, a recursive subspace-based systemspecifically identification used to characterize the transient dynamics of biological the metabolic processes involved Subsequent to of biological systems, systems, specifically the approach metabolicis involved in in diabetes. diabetes. Subsequent to subspace-based systeman identification approach isprocesses used to characterize therecursively transient dynamics of biological systems, specifically the metabolic processes involved in diabetes. Subsequent to system identification, adaptive MPC algorithm is designed using the identified system identification, an adaptive MPC algorithm is designed using the recursively identified of biological systems,compute specifically the metabolic processes involved in the diabetes. Subsequent to system identification, an adaptive MPC algorithm is designed using recursively identified models to effectively the optimal insulin delivery for AP systems. A feature extraction models to effectively compute the optimal insulin delivery for APusing systems. A feature extraction system identification, an adaptive MPC algorithm is designed the recursively identified models to effectively compute the optimal insulin delivery for AP systems. A feature extraction method based based on on glucose glucose measurements measurements is is used used to to detect detect rapid rapid deviations deviations from from the the desired desired method models effectively compute thedisturbances optimal is insulin APmodify systems. A constraints feature extraction method tobased onby glucose measurements used to detectforrapid deviations from the desired set-point caused significant and delivery subsequently the of set-point caused by significant disturbances modify the constraints of the the and subsequently method based onby glucose measurements is used to detect rapid deviations from the desired set-point caused significant disturbances and subsequently modify the constraints of the optimization problem for negotiating between the aggressiveness and robustness of the controller optimization problem for negotiating between the aggressiveness and robustness of the controller set-point caused by significant disturbances andaggressiveness subsequently modify the constraints of the optimization problem for negotiating and robustness of the controller to suggest required amount insulin. The efficacy of of the the proposed adaptive MPC to suggest the the required amount of of between insulin. the proposed adaptive MPC is is The efficacy optimization problem for negotiating between the aggressiveness and robustness of the controller to suggest the required amount of insulin. The efficacy of the proposed adaptive MPC is demonstrated using using simulation simulation case case studies. studies. demonstrated to suggest the required amount of insulin. The efficacy of the proposed adaptive MPC is demonstrated using simulation case studies. demonstrated using simulation case of studies. © 2018, IFAC (International Federation Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Keywords: Keywords: Model Model predictive predictive control, control, Subspace Subspace methods, methods, Recursive Recursive system system identification, identification, Keywords:pancreas Model predictive control, Subspace methods, Recursive system identification, Artificial pancreas Artificial Keywords: Model predictive control, Subspace methods, Recursive system identification, Artificial pancreas Artificial pancreas 1. INTRODUCTION INTRODUCTION for for the the computably computably constrained constrained hardware hardware settings settings prevapreva1. 1. INTRODUCTION for computably constrainedTo hardware settings prevalentthe in biomedical biomedical applications. To address this this limitation, lent in applications. address limitation, 1. INTRODUCTION for the computably constrained hardware settings prevalent in biomedical applications. To address this limitation, advanced step MPC is proposed that solves a full nonlinear Biological processes processes are are complex complex dynamical dynamical systems systems with with advanced step MPC is proposed To that solves this a fulllimitation, nonlinear Biological lent in biomedical applications. address advanced step MPC is proposed thatthough solves athe fullsolution nonlinear at each sampling instance, is Biological are complex dynamical systems with significant processes uncertainties and exogenous exogenous disturbances that program program at each sampling instance, solution is significant uncertainties and disturbances that advanced step MPC is proposed thatthough solves athe full nonlinear Biological processes areand complex dynamical systems with program at each sampling instance, though the solution is computed between sampling instances, with sensitivitiessignificant uncertainties and exogenous disturbances that render their regulation control a challenging problem computed between sampling instances, with sensitivitiesrender their regulation and control a challenging problem program at each sampling instance, though the solution is significant uncertainties and exogenous disturbances that computed between sampling sensitivitiesbased updates updates of the the solutioninstances, performedwith on-line. Another render their regulation and control a challenging problem (Doyle III et al., 2011). Model predictive control (MPC) based of solution performed on-line. Another (Doyle III etregulation al., 2011).and Model predictive controlproblem (MPC) computed between sampling instances, with sensitivitiesrender their control a challenging updates of the solution performedburden on-line.of Another approach to the nonlin(Doyle III et often al., 2011). Model predictive control (MPC) formulations, often used in in biomedical applications, partic- based approach to overcome overcome the computation computation burden nonlinformulations, used biomedical applications, particbased updates of the solution performed on-line.of Another (Doylefor IIIcontrolled et often al., 2011). Model predictive control (MPC) approach to overcome the computation burden of nonlinear programming involves the implementation a formulations, used in biomedical applications, particularly drug delivery, are effective in rejecting ear programming involves the implementation of a realrealularly for controlled drug delivery, are effective in rejecting to overcome the nonlinear computation burden of nonlinformulations, oftenmitigating used indelivery, biomedical applications, partic- approach ear involves the implementation ofwhere a real-a timeprogramming iteration scheme for optimization where a controlled drug are ularly for effective in rejecting disturbances and their effects through receding time iteration scheme for nonlinear optimization disturbances and mitigating their effects through receding ear programming involves the implementation of a real-a ularly for controlled drug delivery, are effective in rejecting time iteration scheme forproblem nonlinear where quadratic programming problem is optimization solved at at each each iteradisturbances and mitigating their effects through receding horizon control, while robust and offset-free control formuquadratic programming is solved iterahorizon control, robust their and offset-free controlreceding formu- time iteration scheme forproblem nonlinear where disturbances andwhile mitigating effects through programming is optimization solved atand each itera-a tion, allowing for active changes ensuring horizon are control, while robust anduncertainty offset-free control formu- quadratic lations are proposed to handle uncertainty in the the models models tion, allowing for multiple multipleproblem active set set changes and ensuring lations proposed to handle in quadratic programming is solved at each iterahorizon control, while robust and offset-free control formuallowing for multiple active set changes ensuring that the algorithm performs no than linear MPC. lations to 2002; handleRawlings uncertainty the models (Muske are andproposed Badgwell, 2002; Rawlings and in Mayne, 2009; tion, that the algorithm performs no worse worse than aa and linear MPC. (Muske and Badgwell, and Mayne, 2009; tion, allowing for multiple active set changes and ensuring lations are proposed to handle uncertainty in the models that the algorithm performs no worse than a linear MPC. In aa different different direction direction involving involving shifting shifting the the computations computations (Muske and Biegler, Badgwell, 2002; Rawlings andFirst-principles Mayne, 2009; In Zavala and and Biegler, 2009; Mayne, 2014). First-principles Zavala 2009; Mayne, 2014). the algorithm performs no worse than linearare MPC. (Muske and and Biegler, Badgwell, 2002; Rawlings andFirst-principles Mayne, 2009; that In a different direction involving shifting theacomputations to an off-line procedure, explicit MPC approaches Zavala 2009; Mayne, 2014). physiological models are often developed to characterize to an off-line procedure, explicit MPC approaches are proprophysiological models 2009; are often developed to characterize In a different direction involving shifting the computations Zavala and Biegler, Mayne, 2014). First-principles an off-line procedure, explicit MPC approaches are proposed that compute the solution of the quadratic programphysiological models oftentime-varying developed toand characterize biological systems dueare to their their time-varying and nonlinear to posed that compute the solution of the quadratic programbiological systems due to nonlinear to an off-line procedure, explicit MPC approaches are prophysiological oftentime-varying developed toand characterize that compute the solution of the of quadratic programming explicitly as the biological systems dueare tovariabilities. their nonlinear dynamics andmodels temporal variabilities. Nevertheless, deriv- posed ming problem problem explicitly as aa function function of the parametrized parametrized dynamics and temporal Nevertheless, derivposed that compute the solution of the quadratic programbiological systems due to their time-varying and nonlinear ming problem explicitly as a function of the parametrized initial state, state, thereby thereby circumventing circumventing the the need need to to solve solve the the dynamics temporal variabilities. Nevertheless, ing models modelsand from first-principles fundamental laws derivoften initial ing from first-principles fundamental laws often ming problem explicitly as a functionthe of need the parametrized dynamics and temporal variabilities. Nevertheless, derivinitial state, thereby circumventing to solve the problems on-line. Using multiparametric programming to ing models from first-principles fundamental laws often yields nonlinear dynamic systems where convex optimizaproblems on-line. Using multiparametric programming to dynamic systems fundamental where convexlaws optimizayields nonlinear initial state, thereby circumventing theaction need to solve the ing models from first-principles often problems on-line. Using multiparametric programming to find explicit control laws, the control can then be yields nonlinear dynamic systems where convex optimization techniques techniques typically typically cannot cannot be be readily readily applied. applied. find explicit control laws, the control action can then be tion on-line. Using multiparametric programming to yieldstechniques nonlinear typically dynamic cannot systemsbewhere convex optimiza- problems find explicit control laws, the control action can then be implemented on-line in the form of expedited calculations readily applied. tion implemented on-line laws, in thethe form of expedited calculations explicit control control action can then be Linear MPC with with convex cannot optimization problems are abunabun- find tion techniques typically be readily applied. implemented on-line in the form of expedited calculations Linear MPC convex optimization problems are involving aa lookup lookup table table of of the the appropriate appropriate function function for for the the implemented on-line in the form of expedited calculations Linear MPC with convex optimization problems aresolvers abun- involving dant because because efficient and reliable interior interior point solvers involving a observed lookup table of the appropriate function for the particular state and reference vectors. Despite dant efficient and reliable point Linear MPC with convex optimization problems are abunparticular observed state and reference vectors. Despite a observed lookup table of the appropriate function for the dant because efficientthat andcan reliable interior point for solvers are widely widely available that can be easily easily compiled for em- involving particular state andefficient reference vectors. all the in solution of matheare available be compiled emdantwidely because efficient andcan reliable interior point for solvers all the improvements improvements in the the efficient solution of Despite matheparticular observed state and reference vectors. Despite are available that be easily compiled embedded hardware and possibly also for field-programmable the improvements in the efficient solution physiologof mathematical programming the bedded hardware andthat possibly also for field-programmable are widely available can be easily compiled em- all matical programming problems, problems, the nonlinear nonlinear physiologall improvements in the efficient solution of mathebedded hardware and possibly also for field-programmable gate arrays. In contrast contrast to linear linear MPC, where for convex matical programming problems, the nonlinear physiologgate arrays. In to MPC, where convex icalthe models are not not tractable tractable enough for implementation bedded hardware and possibly also for field-programmable ical models are enough for implementation matical programming problems, the nonlinear physiologgate arrays. In contrast to linear MPC, where convex quadratic programming programs are solved at each sammodels are not tractable enoughhardware, for implementation on wearable devices with while quadratic programming areMPC, solvedwhere at each sam- ical gate arrays. In nonlinear contrastprograms to linear convex wearable devices with embedded embedded hardware, while the the ical modelsmodels are nothave tractable enough for implementation quadratic programming programs are solved at each sampling instance, (and possible nonconvex) opti- on on wearable devices with embedded hardware, while the nonlinear various time-varying and unique pling instance, nonlinear (and possible nonconvex) optiquadratic programming programs are solved at each samnonlinear models have various time-varying and unique wearable devices with embedded hardware, while the pling instance, nonlinear (and possible opti- on mization problems involving MPC with nonconvex) nonlinear fundafundanonlinear models have various time-varying and unique parameters that must be estimated on-line (Pannocchia mization problems involving MPC with nonlinear pling instance, nonlinear (and possible nonconvex) optiparametersmodels that must be estimated on-line (Pannocchia have various time-varying and unique mizationfirst-principles-based problems involving models MPC with funda- nonlinear mental first-principles-based models are nonlinear not as as tractable tractable parameters that must be estimated on-line (Pannocchia mental are not et al., al., 2007; 2007; Wang Wang and and Boyd, Boyd, 2010; 2010; Bernardini Bernardini and and BemBemmizationfirst-principles-based problems involving models MPC with funda- et parameters that must be estimated on-line (Pannocchia mental are nonlinear not as tractable et al., 2007; Wang and Boyd, 2010;Kouramas Bernardini and Bemporad, 2012; Rivotti et al., 2012; et al., 2013; 1 mental first-principles-based models are not as tractable porad, 2012; Rivotti et al., 2012; Kouramas et al., 2013; 1 Correspondence concerning this article should be addressed to A. et al., et 2007; Wang and Boyd, 2010; Bernardini and BemCorrespondence concerning this article should be addressed to A. porad, 2012; Rivotti et al., 2012; Kouramas et al., 2013; Jerez al., 2014). For many biomedical control applica1 Jerez et al., 2014). For many biomedical control applicaCinar at [email protected] Correspondence concerning this article should be addressed to A. porad, RivottiFor et many al., 2012; Kouramas et al., 2013; Cinar at [email protected] 1 Jerez et2012; al., 2014). biomedical control applicaCorrespondence concerning this article should be addressed to A. Cinar at [email protected] Jerez et al., 2014). For many biomedical control applicaCinar at [email protected]
2405-8963 © 2018, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Copyright 2018 451 Copyright © under 2018 IFAC IFAC 451 Control. Peer review© responsibility of International Federation of Automatic Copyright © 2018 IFAC 451 10.1016/j.ifacol.2018.11.061 Copyright © 2018 IFAC 451
2018 IFAC NMPC Madison, WI, USA, August 19-22, 2018
Iman Hajizadeh et al. / IFAC PapersOnLine 51-20 (2018) 368–373
tions, an identified linear model with constrained inputs is utilized that is well suited for control design, thus circumvent the onerous requirements of first-principles modeling. The general linear models employed for controller synthesis do, however, suffer from a drawback; they do not well describe the nonlinear transient dynamics that are inherent in biological systems. Instead of estimating the important parameters of physiological models on-line, system identification techniques facilitate the recursive updating of the model whenever new data becomes available. Adaptive models also allow for characterizing the time-varying and nonlinear transient behavior of biological systems. Among the computationally viable alternatives to nonlinear models is the subspace-based state-space system identification technique, where structured block Hankel matrices are constructed using measured input/output data and computationally efficient projection operations retrieve certain subspaces related to the system matrices. Recent developments in system identification include a recursive predictor-based subspace identification (PBSID) method that avoids the computationally burdensome procedure of computing or updating singular value decomposition (SVD) factorizations (Houtzager et al., 2012). The adaptive identification allows the realized models to be valid over a diverse range of conditions without requiring in-depth physiological knowledge on the biological systems (Hajizadeh et al., 2017b, 2018a). An example of MPC in biomedical applications is the artificial pancreas (AP) systems that automatically administer the required amount of insulin dose to patients living with the chronic disease of type 1 diabetes to regulate their glucose concentrations (Hovorka et al., 2004; Gondhalekar et al., 2016; Turksoy et al., 2017; Chakrabarty et al., 2018; Messori et al., 2018). As the core of the artificial pancreas is a control algorithm that computes the required dose throughout diverse daily life conditions and across various patients, the use of linear models is of interest to design computational efficient control algorithms. Nevertheless, as physiological characteristics among people differ substantially, and also vary significantly within individuals over time, a general linear model is not sufficient for effective control. To overcome this limitation, personalized and recursive identification techniques can employed to better characterize the time-varying nonlinear dynamics of the biological processes. Motivated by the above considerations, this work proposes an adaptive MPC algorithm for AP systems. To this end, a recursive subspace-based system identification approach is used to characterize the transient dynamics of biological systems, specifically the metabolic processes involved in diabetes. Subsequent to system identification, an adaptive MPC algorithm is designed using the recursively identified models to effectively compute the optimal insulin delivery for AP systems. A feature extraction method based on glucose measurements is used to detect rapid deviations from the desired set-point caused by significant disturbances and subsequently modify the constraints of the optimization problem for negotiating between the aggressiveness and robustness of the controller to suggest the required amount of insulin. The efficacy of the proposed adaptive MPC is demonstrated using the UVa/Padova metabolic simulator (Dalla Man et al., 2014). 452
369
2. PRELIMINARIES In this section, a brief overview of the glucose–insulin dynamic model and the formulation of the nonlinear observer for parameter estimation is provided, followed by a review of the PBSID algorithm for the identification of linear, time-varying state-space models. 2.1 Adaptive Personalized State and Parameter Observer The Hovorka’s model, a glucose–insulin dynamic model, is used for designing an estimator of the insulin present in the bloodstream, called plasma insulin concentration (PIC) (Hovorka et al., 2004; Hajizadeh et al., 2017a). The model has different compartments to characterize the blood glucose dynamics, the subcutaneous insulin infusion, and the glucose transport from plasma to interstitial tissues. Incorporating uncertainty in the dynamics of the discrete-time system (referred to as process noise) and noise in the measurements, the model can be written in the form xk+1 = f (xk , uk , θk ) + wk , wk ∼ N (0, Q(x) ) (1)
(2) yk = g (xk , uk , θk ) + vk , vk ∼ N (0, R(v) ) n m p where x ∈ R , u ∈ R and y ∈ R denote the vectors of state, input, and output variables, respectively, with the output as the continuous glucose monitoring (CGM) measurements and the input variable is the infused insulin. The vectors w ∈ Rn and v ∈ Rp denote the vectors of process and measurement noises, respectively, and θ denotes the vector of uncertain time-varying model parameters. In the simultaneous state and parameter estimation approach, the unscented Kalman filter (UKF) algorithm is used to recursively compute both the state and parameter ˆ estimates x ˆ, θ at each sampling instance (Kol˚ as et al., 2009). To achieve the simultaneous estimation, the original model is transformed by treating the parameters to be T estimated as additional states with x k := xTk θkT the augmented state vector including the uncertain model parameters. The new system for simultaneous state and parameter estimation is f (xk , uk , θk ) (3) f x k , uk := θk g x k , uk := g (xk , uk , θk ) (4) for the where w and v denote the vectors process and ,(x) := (x) (θ) observation noise, Q and R ,(v) := diag Q , Q
R(v) denote the covariance matrix and variance of the process and measurement noises, respectively, with Q(θ) the covariance of the normally distributed random stochastic noise in the parameters. The augmentedstate estimation (θ) (x) error covariance is Pk := diag Pk , Pk . Bounds on the states and parameters can be similarly augmented as x L ≤ T T T T and x U := xTU θU x k ≤ x U , where x L := xTL θL denote the lower and upper bounds on the augmented state vector. Employing the state and parameter estimation approach with Hovorka’s model allows for simultaneously estimating the PIC, which is a state variable in the model, along with other time-varying parameters. A nonlinear observer for the simultaneous estimation of the state variables and the time-varying parameters can be designed as
2018 IFAC NMPC 370 Madison, WI, USA, August 19-22, 2018
Iman Hajizadeh et al. / IFAC PapersOnLine 51-20 (2018) 368–373
ˆk , uk + Kk (yk − yˆk ) x ˆ k+1 = f x ˆk , uk yˆk = g x
(5)
where x ˆ denotes the estimated augmented state vector and K is the observer gain chosen such that the observer error dynamics are asymptotically stable. The proposed PIC estimator can be individualized by appropriately initializing the time-varying parameters using partial least squares regression (Hajizadeh et al., 2017a). The proposed PIC estimator is able to capture the inter- and intrasubject variabilities to provide accurate information on the amount of insulin present in the body, which is used to determine the constraints on the MPC formulation.
and
u u u θk−p θk−p+1 · · · θk−1 .. u u θk−p . θk−2 0 ΓL = . .. . .. ··· ··· u 0 0 · · · θk−f
y y y θk−p θk−p+1 · · · θk−1 .. y y θk−p . θk−2 0 ΓK = . .. . .. ··· ··· y 0 0 · · · θk−f
where f is the future window length. 2.2 Recursive Subspace-Based System Identification A stable, reliable, and computationally tractable dynamic model is essential for the design of model-based predictive control algorithms in AP systems. To this end, the predictor-based subspace identification approach is implemented to track the time-varying linear system and is coupled with a constrained optimization solver to guarantee the stability of the model (Houtzager et al., 2012; Hajizadeh et al., 2017b, 2018a). Consider a vector autoregressive with exogenous variables (VARX) model yˆk|k−1 =
p
u θk−i uk−i +
i=0
p
y θk−i yk−i
(6)
i=1
where yˆk|k−1 is the predicted output for the kth sampling instance using the past inputs uk , . . . , uk−p and outputs yk−1 , . . . , yk−p . The VARX model parameter p is the length of the past window of data considered when predicting the future outputs. The coefficient matrices θu and θy are readily estimated though recursive least square (RLS) techniques at each sampling time. Furthermore, the stacked vector yk−p,p is defined with respect to the past T T T T window of length p as yk−p,p = yk−p yk−p+1 . . . yk−1 . Stacked vector uk−p,p is also similarly defined. Recognizing that the predicted state x ˆk is given by p
x ˆk = A x ˆk−p + Luk−p,p + Kyk−p,p
(7)
x ˆk = Luk−p,p + Kyk−p,p
(8)
Γx ˆk = Γ Luk−p,p + Γ Kyk−p,p
(9)
where Land K denote the extended controllability matri- ces L = Ap−1 B . . . AB B and K = Ap−1 K . . . AK K and assuming that the state transition matrix is nilpotent with degree p, that is the contribution of the initial state x ˆk−p is negligible for sufficiently large p (Ap = 0), the predicted state can be expressed as Pre-multiplying the predicted state by the observability matrix Γ gives with
Γ =
C CA .. . CAp−1
The product of the matrices Γ L and Γ K can be constructed from the VARX model coefficient matrices as 453
Therefore, after estimating the VARX coefficient matrices, the estimated coefficient matrices θu and θy can be used to determine all quantities on the right-hand side of eq. 9, and a SVD factorization can be used to readily obtain a low-rank approximation of the state sequence. For recursive identification, a selection matrix S of appropriate dimensions can be determined such that the basis of the state estimation is consistent at each sampling time as (10) x ˆk = Sk Wk Γ Luk−p,p + Γ Kyk−p,p where Wk is a pre-defined weight matrix and the selection matrix S can be recursively updated through the projection approximation subspace tracking (PAST) method (Chiuso, 2007). The estimated state sequence is then employed along with the inputs and measured outputs to estimate the system matrices by solution of RLS problems. Specifically, after computing an estimate of the state sequence x ˆk , two RLS problems that ensure stability of the estimated system are used to determine the state-space matrices, thus yielding the identified model x ˆk+1 = Ak x ˆ k + B k u k + Kk e k (11) yk = C k x ˆ k + D k uk + e k where A, B, C, and D are the system matrices, K is the Kalman gain matrix, and ek = yk − yˆk for feedback correcting the state variables. The data-driven model determined from the subspacebased system identification technique is integrated with a physiological insulin compartment model derived from Hovorka’s model (Hajizadeh et al., 2018b). Considering a delay of order d and assuming no direct feedthrough term in the state-space model (that is Dk = 0), we can write ˆk + Bk Ik−d + Kk ek x ˆk+1 = Ak x (12) yk = C k x ˆ k + ek The PIC variable in Eq. 12, denoted I, is a filtered insulin to handle time-delays and discrete or impulse variations, which can be obtained from the compartment model with injected basal and bolus insulin as the input variable. To this end, the compartment model in Eq. 1 is integrated with the recursively identified model Eq. 12 yielding ˜ k uk ˜ˆk + B x ˜ˆk+1 = A˜k x (13) ˆ˜k yˆk = C˜k x with the new state vector of T ˆ˜k = x (14) x ˆTk+3+d Ik+2 and the output
yˆk = yˆk+3+d .
(15)
2018 IFAC NMPC Madison, WI, USA, August 19-22, 2018
Iman Hajizadeh et al. / IFAC PapersOnLine 51-20 (2018) 368–373
3. ADAPTIVE MPC ALGORITHM In this section, the glycemic and plasma insulin risk indexes, PIC bounds, the feature extraction for manipulating constraints and the adaptive MPC formulation are presented. Glycemic Risk Index A glycemic risk index (GRI) is used to determine the weighting matrix, denoted Q (ˆ yk+d+3 ), for penalizing the deviations of the outputs from their nominal set-point (Rashid et al., 2018). The glycemic risk index asymmetrically increases the set-point tracking weight when the off-line predicted CGM yˆk+d+3 diverges from the target range. Plasma Insulin Risk Index A plasma insulin penalty index (PIRI), denoted γk+d+3 , is defined to manipulate the weighting matrix for penalizing the amount of input actuation (aggressiveness of insulin dosing) depending on the estimated PIC, thus suppressing the infusion rate if sufficient insulin is present in the bloodstream (Rashid et al., 2018). To this end, the time-varying positive definite weighting matrix Rk := R (γk+d+3 ) is developed from the PIRI as R (γk+d+3 ) := γk+d+3 · R, with R as a nominal weight and γk+d+3 defined as 2 PICk+d+3 γk+d+3 := (16) PICbasal,k+d+3 where
Idb,k+d+3 (17) VI · ke,k and Idb is the patient specific (possibly time-varying) basal insulin rate that is known in practice, and VI and ke are parameters of Hovorka’s model. PICbasal,k+d+3 :=
Plasma Insulin Concentration Bounds In the proposed MPC, the estimated future PIC is dynamically bounded depending on the value of the CGM measurements. For instance, if the CGM measurement values are elevated, the bounds on the PIC are increased to ensure sufficient insulin is administered to regulate the glucose concentration. Furthermore, the PIC bounds also constrain the search space in the optimization problem. The PIC bounds are determined based on the CGM measurements as XPIC := yk ), where XPIC defines the lower and upper (1 + Pmeal )X (ˆ bounds and a desired target for the normalized PIC through the estimated CGM yˆk and Pmeal is a parameter which modifies the bounds when there is a rapid increase in the CGM measurements. 3.1 Feature Extraction for Manipulating Constraints One of the main disturbances that an AP system faces is meal consumption. The blood glucose rises drastically after a meal and it is necessary to counteract this significant disturbance with the required amount of insulin infusion in a correct time. In order to detect meal consumption and give the insulin dose in a timely manner, several meal detection algorithms are developed (Samadi et al., 2017). In this work, we build upon these to handle the effects of meals and other disturbances that result in a rapid and persistent increase in the blood glucose. The new approach modifies the constraints if there is any rapid and significant 454
371
increase in the glucose measurements. Meal detection can be based on the qualitative description of segments of glucose time series data (Samadi et al., 2017). In this approach, features from the glucose measurement data are extracted using the first and second derivatives of a secondorder polynomial fitted to the last few measurement data. If the first-order derivative is positive and greater than a threshold while taking into consideration the sign of second derivative, the constraints are modified to make the controller more aggressive to suggest a sufficient insulin dose. At every sampling time, the approach estimates the coefficients of the following polynomial function Yk = pn1 × tk 2 + pn2 × tk 1 + pn3 .
(18)
where pn1 , pn2 and pn3 are the parameters of a secondorder polynomial fitted to the last few glucose measurement and tk is the sampling time. At each sampling instance, by updating Eq. 18, the Pmeal parameter is calculated by the first and second derivatives as follows: Yk/T dm,1 if Y ≥ T d m,11 , Yk ≥ 0 k Pmeal = (19) Yk/T dm,2 if Y ≥ T d m,11 , Yk < 0 k
where the Yk and Yk are the first and second order derivatives at the last sampling time, the T dm,1 , T dm,2 and T dm,11 are patient-specific threshold parameters. 3.2 Adaptive MPC Formulation In this subsection, we propose a novel adaptive MPC algorithm cognizant of the PIC for computing the optimal insulin infusion rate. The proposed MPC formulation employs the glycemic and PIC risk indexes that manipulate the penalty weighting matrices in the cost function. To this end, the MPC computes the optimal insulin infusion over a finite horizon using the identified time-varying subspacebased models by solving at each ith sampling instance the following quadratic programming problem nP P := arg min JnP ,i yˆi+d+3 , γk+d+3 , {uk }nk=0 {u∗k }k=0 u∈U ˆ˜k+1 = Ak x ˆ˜k + Bk uk−d−3 , ∀k ∈ ZnP −1 x 0 ˆ˜k , ∀k ∈ ZnP s.t. yˆk = Ck x 0 x ˆ˜0 = x ˆ˜i x ∈ X (20) with the objective function nP nP JnP ,i (·) := eTk Qi ek + uTk Ri uk k=0
k=0
ˆ˜ and yˆ denote the predicted states and outwhere x puts, respectively, for the prediction/control horizon nP , u ∈ Rm denotes the vector of constrained input variables, taking values in a nonempty convex set U ⊆ Rm with U := {u ∈ Rm : umin ≤ u ≤ umax }, umin ∈ Rm and umax ∈ Rm denote the lower and upper bounds on the manipulated input, respectively, and ek := yˆk − ysp . The index Zn0 P represents all integers in a set as Zn0 P := n ˜ {0, . . . , n˜nP }. The nonemptyconvex set Xn˜ ⊆ R with X := x ∈ R : xmin ≤ x ≤ xmax , xmin ∈ R and xmax ∈ Rn˜ denote the lower and upper bounds on the state variables, respectively, with one of the states as the estimated PIC that is constrained through the PIC bounds. Furthermore, ˆ˜i provides an initialization of the state vector, Q ≥ 0, x
2018 IFAC NMPC 372 Madison, WI, USA, August 19-22, 2018
Iman Hajizadeh et al. / IFAC PapersOnLine 51-20 (2018) 368–373
Table 1. Closed-loop Simulation Results Subject Adult Adult Adult Adult Adult Adult Adult Adult Adult Adult
1 2 3 4 5 6 7 8 9 10
Average
< 55
[55, 70)
> 250
Mean
Statistics SD Min
Max
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
Percent of time in range [80, 140] [70, 180) [180, 250) 43.5 41.0 59.6 38.3 43.4 51.2 44.0 39.5 43.7 48.4
89.7 65.9 80.5 61.9 61.6 73.2 76.6 71.4 83.2 82.2
10.3 34.1 19.5 38.1 35.2 26.8 23.3 28.6 16.8 17.8
0.0 0.0 0.0 0.0 3.3 0.0 0.1 0.0 0.0 0.0
138.9 150.4 137.2 153.9 163.4 144.7 150.1 150.6 141.3 141.0
35.1 46.1 37.7 45.9 47.4 44.3 41.9 42.0 36.6 37.4
83.0 83.0 76.0 79.0 74.0 80.0 98.0 90.0 87.0 85.0
214.0 236.0 224.0 233.0 287.0 236.0 252.0 236.0 218.0 233.0
0.0
0.0
45.3
74.6
25.0
0.3
147.1
41.4
83.5
236.9
Qk := Q (ˆ yk+d+3 ) is a positive semi-definite symmetric matrix used to penalize the deviations of the outputs from their nominal set-point, and R > 0, Rk := R (γk+d+3 ) is a strictly positive definite symmetric matrix to penalize the manipulated input variables. 4. RESULTS The efficacy of the proposed PIC cognizant MPC with the CGM feature extraction method is demonstrated using the UVa/Padova metabolic simulator. Overall 10 subjects are simulated for 12 days with varying times and quantities of meals consumed on each day. The times and amounts of meals are selected from specified ranges to evaluate the proposed MPC under real challenging realistic scenarios. The meal information is not utilized in the MPC algorithm as the controller is designed to regulate the (blood glucose concentration) BGC in the presence of significant disturbances such as unannounced meals. In order to handle the meal effects efficiently, the CGM feature extraction method is designed to make the MPC more aggressive. The controller set-point is specified to be 110 mg/dL. The model is designed with a delay of order d = 1 with regards to the effect of PIC on CGM. The quantitative evaluation of the closed-loop results based on the proposed PIC cognizant MPC algorithm using the CGM feature extraction method are presented in Table 1, which gives the percentage of samples in defined glycemic ranges and selected statistics for the glucose measurements. It is readily observed that no hypoglycemia occurs as the BGC is never below 70 mg/dL. The average percentage of time spent in the target range (BGC ∈ [70, 180] mg/dL) and the higher tier of BGC ∈ [180, 250] mg/dL are 74.6% and 25.0%, respectively. The BGC is thus tightly controlled to be within the safe range for a significant duration of time (Fig. 1). Although the BGC is tightly controlled to be within the safe range in both cases (Fig. 2), the results with the MPC that uses the feature extraction method show a superior performance especially during day time when the AP system is challenged by meals. The average percentage of time spent in the target range (BGC ∈ [70, 180] mg/dL) is improved by 12% and the mean and standard deviation of BGC are reduced by 10 and 5 mg/dL, respectively. Furthermore, the average maximum value of BGC is reduced by 17 mg/dL without any significant change in the average minimum value of BGC (a decrease of 0.5 mg/dL) and the injected insulin (manipulated variable) (an increase of 0.5 Unit of insulin per day). 5. DISCUSSION In this work, a physiological model drawn from the insulin compartment of Hovorka’s model is incorporated with a 455
Fig. 1. Closed-loop results for 10 adult subjects using the feature extraction method
Fig. 2. Comparison of MPC results without (gray) and with (white) the feature extraction method. The bottom and top of the boxes are the first and third quartiles and the line inside the box is the median. The whiskers ends represents minimum and maximum values and + indicates mean values. data-driven model developed using a recursive subspace identification technique. The subspace identification technique with recursive modeling renders the identified models valid over a diverse range of daily activities without requiring onerous and obscure information such as the amount of carbohydrate consumption. In addition to the PIC constraints, in order to better compensate for the effects of meal consumption, the rate of change of the glucose values (as determined by the derivatives of the CGM measurements) are computed in real time to modify the aggressiveness of the controller and give the required quantity of insulin in an appropriate amount of time. The proposed approach involving adaptation of the model, the penalty weights of the optimization cost function, and the constraints is necessary as the nonlinear and transient operating behavior of biological systems is difficult to characterize precisely. The adaptive models can well characterize the nonlinear and time-varying dynamics across various operating conditions by recursively updating the identified
2018 IFAC NMPC Madison, WI, USA, August 19-22, 2018
Iman Hajizadeh et al. / IFAC PapersOnLine 51-20 (2018) 368–373
models, thereby tailoring the models to the particular conditions of any given situation. This adaptive nature of the proposed technique results in good control performance, as demonstrated through the tight glucose control using the metabolic simulator. 6. CONCLUSION An adaptive MPC is proposed in this work for optimal insulin dosing decisions in AP systems using a recursive subspace-based system identification approach to characterize the transient dynamics of processes involved in diabetes. Subsequent to system identification, an adaptive MPC algorithm is designed using the recursively identified models to effectively compute the optimal insulin delivery for AP systems. A feature extraction method based on glucose measurements is used to detect rapid deviations from the desired set-point caused by significant disturbances and subsequently modify the constraints of the optimization problem for negotiating between the aggressiveness and robustness of the controller to suggest the required amount of insulin. REFERENCES Bernardini, D. and Bemporad, A. (2012). Stabilizing model predictive control of stochastic constrained linear systems. IEEE Trans Automat Contr, 57(6), 1468–1480. Chakrabarty, A., Zavitsanou, S., Doyle, F.J., and Dassau, E. (2018). Event-triggered model predictive control for embedded artificial pancreas systems. IEEE Trans. Biomed. Eng., 65(3), 575–586. Chiuso, A. (2007). The role of vector autoregressive modeling in predictor-based subspace identification. Automatica, 43(6), 1034–1048. Dalla Man, C., Micheletto, F., Lv, D., Breton, M., Kovatchev, B., and Cobelli, C. (2014). The UVa/Padova type 1 diabetes simulator: new features. J. Diabetes Sci. Technol., 8(1), 26–34. Doyle III, F.J., Bequette, B.W., Middleton, R., Ogunnaike, B., Paden, B., Parker, R.S., and Vidyasagar, M. (2011). Control in biological systems. In T. Samad, & A. M. Annaswamy (Eds.), The impact of control technology, part 1. IEEE Control Systems Society. Gondhalekar, R., Dassau, E., and Doyle, F.J. (2016). Periodic zone-MPC with asymmetric costs for outpatientready safety of an artificial pancreas to treat type 1 diabetes. Automatica, 71, 237 – 246. Hajizadeh, I., Rashid, M., and Cinar, A. (2018a). Ensuring stability and fidelity of recursively identified controlrelevant models. In The 18th IFAC Symposium on System Identification (SYSID), 927–932. Hajizadeh, I., Rashid, M., and Cinar, A. (2018b). Integrating compartment models with recursive system identification. In American Control Conference (ACC), 3583–3588. Hajizadeh, I., Rashid, M., Turksoy, K., Samadi, S., Feng, J., Frantz, N., Sevil, M., Cengiz, E., and Cinar, A. (2017a). Plasma insulin estimation in people with type 1 diabetes mellitus. Ind. Eng. Chem. Res., 56(35), 9846– 9857. Hajizadeh, I., Rashid, M., Turksoy, K., Samadi, S., Feng, J., Sevili, M., Frantz, N., Lazaro, C., Maloney, Z., Littlejohn, E., and Cinar, A. (2017b). Multivariable 456
373
recursive subspace identification with application to artificial pancreas systems. IFAC-PapersOnLine, 909 – 914. Houtzager, I., van Wingerden, J.W., and Verhaegen, M. (2012). Recursive predictor-based subspace identification with application to the real-time closed-loop tracking of flutter. IEEE Trans. Control Syst. Technol., 20(4), 934–949. 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. Physiol. Meas., 25(4), 905. Jerez, J.L., Goulart, P.J., Richter, S., Constantinides, G.A., Kerrigan, E.C., and Morari, M. (2014). Embedded online optimization for model predictive control at megahertz rates. IEEE Trans Automat Contr, 59(12), 3238–3251. Kol˚ as, S., Foss, B., and Schei, T. (2009). Constrained nonlinear state estimation based on the UKF approach. Comput. Chem. Eng., 33(8), 1386–1401. Kouramas, K., Panos, C., Fasca, N., and Pistikopoulos, E. (2013). An algorithm for robust explicit/multiparametric model predictive control. Automatica, 49(2), 381 – 389. Mayne, D.Q. (2014). Model predictive control: Recent developments and future promise. Automatica, 50(12), 2967 – 2986. Messori, M., Incremona, G.P., Cobelli, C., and Magni, L. (2018). Individualized model predictive control for the artificial pancreas: In silico evaluation of closed-loop glucose control. IEEE Control Syst. Mag., 38(1), 86–104. Muske, K.R. and Badgwell, T.A. (2002). Disturbance modeling for offset-free linear model predictive control. J. Process Control, 12(5), 617 – 632. Pannocchia, G., Rawlings, J.B., and Wright, S.J. (2007). Fast, large-scale model predictive control by partial enumeration. Automatica, 43(5), 852 – 860. Rashid, M., Hajizadeh, I., and Cinar, A. (2018). Plasma insulin cognizant predictive control for artificial pancreas. In American Control Conference (ACC), 3589–3594. Rawlings, J.B. and Mayne, D.Q. (2009). Model predictive control: Theory and design. Nob Hill Pub. Rivotti, P., Lambert, R.S., and Pistikopoulos, E.N. (2012). Combined model approximation techniques and multiparametric programming for explicit nonlinear model predictive control. Comput. Chem. Eng., 42, 277 – 287. Samadi, S., Turksoy, K., Hajizadeh, I., Feng, J., Sevil, M., and Cinar, A. (2017). Meal detection and carbohydrate estimation using continuous glucose sensor data. IEEE J. Biomed. Health Inform., 21(3), 619–627. Turksoy, K., Hajizadeh, I., Samadi, S., Feng, J., Sevil, M., Park, M., Quinn, L., Littlejohn, E., and Cinar, A. (2017). Real-time insulin bolusing for unannounced meals with artificial pancreas. Control Eng. Pract., 59, 159–164. Wang, Y. and Boyd, S. (2010). Fast model predictive control using online optimization. IEEE Trans Control Syst Technol, 18(2), 267–278. Zavala, V.M. and Biegler, L.T. (2009). The advanced-step NMPC controller: Optimality, stability and robustness. Automatica, 45(1), 86 – 93.
ScienceDirect
IFAC PapersOnLine 51-20 (2018) 368–373
Adaptive Model Predictive Control for Adaptive Model Predictive Control for Adaptive Model Predictive Control for Adaptive Model Predictive Control for Nonlinearity in Biomedical Applications Nonlinearity in Biomedical Applications Nonlinearity in Biomedical Applications Nonlinearity in† Biomedical Applications † ‡
† Mudassir Rashid † Mert Sevil ‡ Iman Iman Hajizadeh Hajizadeh Mert Sevil ‡ † Mudassir Rashid † † Mert Sevil ‡‡ Iman Hajizadeh Rashid † Rachel Brandt ‡‡ Sediqeh Sediqeh Samadi Nicole Hobbs ‡ † Mudassir † Rachel Brandt Samadi Iman Hajizadeh Mudassir Rashid MertHobbs Sevil ‡ ‡ † Nicole †,‡,1 Rachel Brandt ‡ Sediqeh Samadi †,‡,1 † Nicole Hobbs ‡ Ali Cinar Ali Cinar Rachel Brandt Sediqeh Samadi Nicole Hobbs †,‡,1 Ali Cinar †,‡,1 Ali Cinar † † Department of Chemical and Biological Engineering, † Department of Chemical and Biological Engineering, of Chemical and Biological Illinois Institute of Technology, Chicago, IL, † Department Illinois Institute of Technology, Chicago,Engineering, IL, USA USA Department of Chemical and Biological Engineering, ‡ Institute Illinois ofof Chicago, IL, USA ‡ Department ofTechnology, Biomedical Engineering, Department Biomedical Engineering, Illinois ‡ Institute of Technology, Chicago, IL, USA BiomedicalChicago, Engineering, Illinois Institute IL, ‡ Department Illinois Institute of ofofTechnology, Technology, IL, USA USA BiomedicalChicago, Engineering, IllinoisDepartment Institute ofofTechnology, Chicago, IL, USA Illinois Institute of Technology, Chicago, IL, USA Abstract: Abstract: An An adaptive adaptive model model predictive predictive control control (MPC) (MPC) formulation formulation is is proposed proposed in in this this work work Abstract: An adaptive model predictive control (MPC) formulation is proposed in work for optimal insulin dosing decisions in artificial pancreas (AP) systems. To this end, a recursive for optimal An insulin dosingmodel decisions in artificial pancreas (AP) systems.isTo this end,ina this recursive Abstract: adaptive predictive control (MPC) formulation proposed this work for optimal insulin dosing decisions in approach artificial pancreas systems. the To this end, a dynamics recursive subspace-based system identification approach is used used to to(AP) characterize the transient dynamics subspace-based system identification is characterize transient for optimal insulin dosing decisions in artificial pancreas (AP) systems. To this end, a recursive subspace-based systemspecifically identification used to characterize the transient dynamics of biological the metabolic processes involved Subsequent to of biological systems, systems, specifically the approach metabolicis involved in in diabetes. diabetes. Subsequent to subspace-based systeman identification approach isprocesses used to characterize therecursively transient dynamics of biological systems, specifically the metabolic processes involved in diabetes. Subsequent to system identification, adaptive MPC algorithm is designed using the identified system identification, an adaptive MPC algorithm is designed using the recursively identified of biological systems,compute specifically the metabolic processes involved in the diabetes. Subsequent to system identification, an adaptive MPC algorithm is designed using recursively identified models to effectively the optimal insulin delivery for AP systems. A feature extraction models to effectively compute the optimal insulin delivery for APusing systems. A feature extraction system identification, an adaptive MPC algorithm is designed the recursively identified models to effectively compute the optimal insulin delivery for AP systems. A feature extraction method based based on on glucose glucose measurements measurements is is used used to to detect detect rapid rapid deviations deviations from from the the desired desired method models effectively compute thedisturbances optimal is insulin APmodify systems. A constraints feature extraction method tobased onby glucose measurements used to detectforrapid deviations from the desired set-point caused significant and delivery subsequently the of set-point caused by significant disturbances modify the constraints of the the and subsequently method based onby glucose measurements is used to detect rapid deviations from the desired set-point caused significant disturbances and subsequently modify the constraints of the optimization problem for negotiating between the aggressiveness and robustness of the controller optimization problem for negotiating between the aggressiveness and robustness of the controller set-point caused by significant disturbances andaggressiveness subsequently modify the constraints of the optimization problem for negotiating and robustness of the controller to suggest required amount insulin. The efficacy of of the the proposed adaptive MPC to suggest the the required amount of of between insulin. the proposed adaptive MPC is is The efficacy optimization problem for negotiating between the aggressiveness and robustness of the controller to suggest the required amount of insulin. The efficacy of the proposed adaptive MPC is demonstrated using using simulation simulation case case studies. studies. demonstrated to suggest the required amount of insulin. The efficacy of the proposed adaptive MPC is demonstrated using simulation case studies. demonstrated using simulation case of studies. © 2018, IFAC (International Federation Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Keywords: Keywords: Model Model predictive predictive control, control, Subspace Subspace methods, methods, Recursive Recursive system system identification, identification, Keywords:pancreas Model predictive control, Subspace methods, Recursive system identification, Artificial pancreas Artificial Keywords: Model predictive control, Subspace methods, Recursive system identification, Artificial pancreas Artificial pancreas 1. INTRODUCTION INTRODUCTION for for the the computably computably constrained constrained hardware hardware settings settings prevapreva1. 1. INTRODUCTION for computably constrainedTo hardware settings prevalentthe in biomedical biomedical applications. To address this this limitation, lent in applications. address limitation, 1. INTRODUCTION for the computably constrained hardware settings prevalent in biomedical applications. To address this limitation, advanced step MPC is proposed that solves a full nonlinear Biological processes processes are are complex complex dynamical dynamical systems systems with with advanced step MPC is proposed To that solves this a fulllimitation, nonlinear Biological lent in biomedical applications. address advanced step MPC is proposed thatthough solves athe fullsolution nonlinear at each sampling instance, is Biological are complex dynamical systems with significant processes uncertainties and exogenous exogenous disturbances that program program at each sampling instance, solution is significant uncertainties and disturbances that advanced step MPC is proposed thatthough solves athe full nonlinear Biological processes areand complex dynamical systems with program at each sampling instance, though the solution is computed between sampling instances, with sensitivitiessignificant uncertainties and exogenous disturbances that render their regulation control a challenging problem computed between sampling instances, with sensitivitiesrender their regulation and control a challenging problem program at each sampling instance, though the solution is significant uncertainties and exogenous disturbances that computed between sampling sensitivitiesbased updates updates of the the solutioninstances, performedwith on-line. Another render their regulation and control a challenging problem (Doyle III et al., 2011). Model predictive control (MPC) based of solution performed on-line. Another (Doyle III etregulation al., 2011).and Model predictive controlproblem (MPC) computed between sampling instances, with sensitivitiesrender their control a challenging updates of the solution performedburden on-line.of Another approach to the nonlin(Doyle III et often al., 2011). Model predictive control (MPC) formulations, often used in in biomedical applications, partic- based approach to overcome overcome the computation computation burden nonlinformulations, used biomedical applications, particbased updates of the solution performed on-line.of Another (Doylefor IIIcontrolled et often al., 2011). Model predictive control (MPC) approach to overcome the computation burden of nonlinear programming involves the implementation a formulations, used in biomedical applications, particularly drug delivery, are effective in rejecting ear programming involves the implementation of a realrealularly for controlled drug delivery, are effective in rejecting to overcome the nonlinear computation burden of nonlinformulations, oftenmitigating used indelivery, biomedical applications, partic- approach ear involves the implementation ofwhere a real-a timeprogramming iteration scheme for optimization where a controlled drug are ularly for effective in rejecting disturbances and their effects through receding time iteration scheme for nonlinear optimization disturbances and mitigating their effects through receding ear programming involves the implementation of a real-a ularly for controlled drug delivery, are effective in rejecting time iteration scheme forproblem nonlinear where quadratic programming problem is optimization solved at at each each iteradisturbances and mitigating their effects through receding horizon control, while robust and offset-free control formuquadratic programming is solved iterahorizon control, robust their and offset-free controlreceding formu- time iteration scheme forproblem nonlinear where disturbances andwhile mitigating effects through programming is optimization solved atand each itera-a tion, allowing for active changes ensuring horizon are control, while robust anduncertainty offset-free control formu- quadratic lations are proposed to handle uncertainty in the the models models tion, allowing for multiple multipleproblem active set set changes and ensuring lations proposed to handle in quadratic programming is solved at each iterahorizon control, while robust and offset-free control formuallowing for multiple active set changes ensuring that the algorithm performs no than linear MPC. lations to 2002; handleRawlings uncertainty the models (Muske are andproposed Badgwell, 2002; Rawlings and in Mayne, 2009; tion, that the algorithm performs no worse worse than aa and linear MPC. (Muske and Badgwell, and Mayne, 2009; tion, allowing for multiple active set changes and ensuring lations are proposed to handle uncertainty in the models that the algorithm performs no worse than a linear MPC. In aa different different direction direction involving involving shifting shifting the the computations computations (Muske and Biegler, Badgwell, 2002; Rawlings andFirst-principles Mayne, 2009; In Zavala and and Biegler, 2009; Mayne, 2014). First-principles Zavala 2009; Mayne, 2014). the algorithm performs no worse than linearare MPC. (Muske and and Biegler, Badgwell, 2002; Rawlings andFirst-principles Mayne, 2009; that In a different direction involving shifting theacomputations to an off-line procedure, explicit MPC approaches Zavala 2009; Mayne, 2014). physiological models are often developed to characterize to an off-line procedure, explicit MPC approaches are proprophysiological models 2009; are often developed to characterize In a different direction involving shifting the computations Zavala and Biegler, Mayne, 2014). First-principles an off-line procedure, explicit MPC approaches are proposed that compute the solution of the quadratic programphysiological models oftentime-varying developed toand characterize biological systems dueare to their their time-varying and nonlinear to posed that compute the solution of the quadratic programbiological systems due to nonlinear to an off-line procedure, explicit MPC approaches are prophysiological oftentime-varying developed toand characterize that compute the solution of the of quadratic programming explicitly as the biological systems dueare tovariabilities. their nonlinear dynamics andmodels temporal variabilities. Nevertheless, deriv- posed ming problem problem explicitly as aa function function of the parametrized parametrized dynamics and temporal Nevertheless, derivposed that compute the solution of the quadratic programbiological systems due to their time-varying and nonlinear ming problem explicitly as a function of the parametrized initial state, state, thereby thereby circumventing circumventing the the need need to to solve solve the the dynamics temporal variabilities. Nevertheless, ing models modelsand from first-principles fundamental laws derivoften initial ing from first-principles fundamental laws often ming problem explicitly as a functionthe of need the parametrized dynamics and temporal variabilities. Nevertheless, derivinitial state, thereby circumventing to solve the problems on-line. Using multiparametric programming to ing models from first-principles fundamental laws often yields nonlinear dynamic systems where convex optimizaproblems on-line. Using multiparametric programming to dynamic systems fundamental where convexlaws optimizayields nonlinear initial state, thereby circumventing theaction need to solve the ing models from first-principles often problems on-line. Using multiparametric programming to find explicit control laws, the control can then be yields nonlinear dynamic systems where convex optimization techniques techniques typically typically cannot cannot be be readily readily applied. applied. find explicit control laws, the control action can then be tion on-line. Using multiparametric programming to yieldstechniques nonlinear typically dynamic cannot systemsbewhere convex optimiza- problems find explicit control laws, the control action can then be implemented on-line in the form of expedited calculations readily applied. tion implemented on-line laws, in thethe form of expedited calculations explicit control control action can then be Linear MPC with with convex cannot optimization problems are abunabun- find tion techniques typically be readily applied. implemented on-line in the form of expedited calculations Linear MPC convex optimization problems are involving aa lookup lookup table table of of the the appropriate appropriate function function for for the the implemented on-line in the form of expedited calculations Linear MPC with convex optimization problems aresolvers abun- involving dant because because efficient and reliable interior interior point solvers involving a observed lookup table of the appropriate function for the particular state and reference vectors. Despite dant efficient and reliable point Linear MPC with convex optimization problems are abunparticular observed state and reference vectors. Despite a observed lookup table of the appropriate function for the dant because efficientthat andcan reliable interior point for solvers are widely widely available that can be easily easily compiled for em- involving particular state andefficient reference vectors. all the in solution of matheare available be compiled emdantwidely because efficient andcan reliable interior point for solvers all the improvements improvements in the the efficient solution of Despite matheparticular observed state and reference vectors. Despite are available that be easily compiled embedded hardware and possibly also for field-programmable the improvements in the efficient solution physiologof mathematical programming the bedded hardware andthat possibly also for field-programmable are widely available can be easily compiled em- all matical programming problems, problems, the nonlinear nonlinear physiologall improvements in the efficient solution of mathebedded hardware and possibly also for field-programmable gate arrays. In contrast contrast to linear linear MPC, where for convex matical programming problems, the nonlinear physiologgate arrays. In to MPC, where convex icalthe models are not not tractable tractable enough for implementation bedded hardware and possibly also for field-programmable ical models are enough for implementation matical programming problems, the nonlinear physiologgate arrays. In contrast to linear MPC, where convex quadratic programming programs are solved at each sammodels are not tractable enoughhardware, for implementation on wearable devices with while quadratic programming areMPC, solvedwhere at each sam- ical gate arrays. In nonlinear contrastprograms to linear convex wearable devices with embedded embedded hardware, while the the ical modelsmodels are nothave tractable enough for implementation quadratic programming programs are solved at each sampling instance, (and possible nonconvex) opti- on on wearable devices with embedded hardware, while the nonlinear various time-varying and unique pling instance, nonlinear (and possible nonconvex) optiquadratic programming programs are solved at each samnonlinear models have various time-varying and unique wearable devices with embedded hardware, while the pling instance, nonlinear (and possible opti- on mization problems involving MPC with nonconvex) nonlinear fundafundanonlinear models have various time-varying and unique parameters that must be estimated on-line (Pannocchia mization problems involving MPC with nonlinear pling instance, nonlinear (and possible nonconvex) optiparametersmodels that must be estimated on-line (Pannocchia have various time-varying and unique mizationfirst-principles-based problems involving models MPC with funda- nonlinear mental first-principles-based models are nonlinear not as as tractable tractable parameters that must be estimated on-line (Pannocchia mental are not et al., al., 2007; 2007; Wang Wang and and Boyd, Boyd, 2010; 2010; Bernardini Bernardini and and BemBemmizationfirst-principles-based problems involving models MPC with funda- et parameters that must be estimated on-line (Pannocchia mental are nonlinear not as tractable et al., 2007; Wang and Boyd, 2010;Kouramas Bernardini and Bemporad, 2012; Rivotti et al., 2012; et al., 2013; 1 mental first-principles-based models are not as tractable porad, 2012; Rivotti et al., 2012; Kouramas et al., 2013; 1 Correspondence concerning this article should be addressed to A. et al., et 2007; Wang and Boyd, 2010; Bernardini and BemCorrespondence concerning this article should be addressed to A. porad, 2012; Rivotti et al., 2012; Kouramas et al., 2013; Jerez al., 2014). For many biomedical control applica1 Jerez et al., 2014). For many biomedical control applicaCinar at [email protected] Correspondence concerning this article should be addressed to A. porad, RivottiFor et many al., 2012; Kouramas et al., 2013; Cinar at [email protected] 1 Jerez et2012; al., 2014). biomedical control applicaCorrespondence concerning this article should be addressed to A. Cinar at [email protected] Jerez et al., 2014). For many biomedical control applicaCinar at [email protected]
2405-8963 © 2018, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Copyright 2018 451 Copyright © under 2018 IFAC IFAC 451 Control. Peer review© responsibility of International Federation of Automatic Copyright © 2018 IFAC 451 10.1016/j.ifacol.2018.11.061 Copyright © 2018 IFAC 451
2018 IFAC NMPC Madison, WI, USA, August 19-22, 2018
Iman Hajizadeh et al. / IFAC PapersOnLine 51-20 (2018) 368–373
tions, an identified linear model with constrained inputs is utilized that is well suited for control design, thus circumvent the onerous requirements of first-principles modeling. The general linear models employed for controller synthesis do, however, suffer from a drawback; they do not well describe the nonlinear transient dynamics that are inherent in biological systems. Instead of estimating the important parameters of physiological models on-line, system identification techniques facilitate the recursive updating of the model whenever new data becomes available. Adaptive models also allow for characterizing the time-varying and nonlinear transient behavior of biological systems. Among the computationally viable alternatives to nonlinear models is the subspace-based state-space system identification technique, where structured block Hankel matrices are constructed using measured input/output data and computationally efficient projection operations retrieve certain subspaces related to the system matrices. Recent developments in system identification include a recursive predictor-based subspace identification (PBSID) method that avoids the computationally burdensome procedure of computing or updating singular value decomposition (SVD) factorizations (Houtzager et al., 2012). The adaptive identification allows the realized models to be valid over a diverse range of conditions without requiring in-depth physiological knowledge on the biological systems (Hajizadeh et al., 2017b, 2018a). An example of MPC in biomedical applications is the artificial pancreas (AP) systems that automatically administer the required amount of insulin dose to patients living with the chronic disease of type 1 diabetes to regulate their glucose concentrations (Hovorka et al., 2004; Gondhalekar et al., 2016; Turksoy et al., 2017; Chakrabarty et al., 2018; Messori et al., 2018). As the core of the artificial pancreas is a control algorithm that computes the required dose throughout diverse daily life conditions and across various patients, the use of linear models is of interest to design computational efficient control algorithms. Nevertheless, as physiological characteristics among people differ substantially, and also vary significantly within individuals over time, a general linear model is not sufficient for effective control. To overcome this limitation, personalized and recursive identification techniques can employed to better characterize the time-varying nonlinear dynamics of the biological processes. Motivated by the above considerations, this work proposes an adaptive MPC algorithm for AP systems. To this end, a recursive subspace-based system identification approach is used to characterize the transient dynamics of biological systems, specifically the metabolic processes involved in diabetes. Subsequent to system identification, an adaptive MPC algorithm is designed using the recursively identified models to effectively compute the optimal insulin delivery for AP systems. A feature extraction method based on glucose measurements is used to detect rapid deviations from the desired set-point caused by significant disturbances and subsequently modify the constraints of the optimization problem for negotiating between the aggressiveness and robustness of the controller to suggest the required amount of insulin. The efficacy of the proposed adaptive MPC is demonstrated using the UVa/Padova metabolic simulator (Dalla Man et al., 2014). 452
369
2. PRELIMINARIES In this section, a brief overview of the glucose–insulin dynamic model and the formulation of the nonlinear observer for parameter estimation is provided, followed by a review of the PBSID algorithm for the identification of linear, time-varying state-space models. 2.1 Adaptive Personalized State and Parameter Observer The Hovorka’s model, a glucose–insulin dynamic model, is used for designing an estimator of the insulin present in the bloodstream, called plasma insulin concentration (PIC) (Hovorka et al., 2004; Hajizadeh et al., 2017a). The model has different compartments to characterize the blood glucose dynamics, the subcutaneous insulin infusion, and the glucose transport from plasma to interstitial tissues. Incorporating uncertainty in the dynamics of the discrete-time system (referred to as process noise) and noise in the measurements, the model can be written in the form xk+1 = f (xk , uk , θk ) + wk , wk ∼ N (0, Q(x) ) (1)
(2) yk = g (xk , uk , θk ) + vk , vk ∼ N (0, R(v) ) n m p where x ∈ R , u ∈ R and y ∈ R denote the vectors of state, input, and output variables, respectively, with the output as the continuous glucose monitoring (CGM) measurements and the input variable is the infused insulin. The vectors w ∈ Rn and v ∈ Rp denote the vectors of process and measurement noises, respectively, and θ denotes the vector of uncertain time-varying model parameters. In the simultaneous state and parameter estimation approach, the unscented Kalman filter (UKF) algorithm is used to recursively compute both the state and parameter ˆ estimates x ˆ, θ at each sampling instance (Kol˚ as et al., 2009). To achieve the simultaneous estimation, the original model is transformed by treating the parameters to be T estimated as additional states with x k := xTk θkT the augmented state vector including the uncertain model parameters. The new system for simultaneous state and parameter estimation is f (xk , uk , θk ) (3) f x k , uk := θk g x k , uk := g (xk , uk , θk ) (4) for the where w and v denote the vectors process and ,(x) := (x) (θ) observation noise, Q and R ,(v) := diag Q , Q
R(v) denote the covariance matrix and variance of the process and measurement noises, respectively, with Q(θ) the covariance of the normally distributed random stochastic noise in the parameters. The augmentedstate estimation (θ) (x) error covariance is Pk := diag Pk , Pk . Bounds on the states and parameters can be similarly augmented as x L ≤ T T T T and x U := xTU θU x k ≤ x U , where x L := xTL θL denote the lower and upper bounds on the augmented state vector. Employing the state and parameter estimation approach with Hovorka’s model allows for simultaneously estimating the PIC, which is a state variable in the model, along with other time-varying parameters. A nonlinear observer for the simultaneous estimation of the state variables and the time-varying parameters can be designed as
2018 IFAC NMPC 370 Madison, WI, USA, August 19-22, 2018
Iman Hajizadeh et al. / IFAC PapersOnLine 51-20 (2018) 368–373
ˆk , uk + Kk (yk − yˆk ) x ˆ k+1 = f x ˆk , uk yˆk = g x
(5)
where x ˆ denotes the estimated augmented state vector and K is the observer gain chosen such that the observer error dynamics are asymptotically stable. The proposed PIC estimator can be individualized by appropriately initializing the time-varying parameters using partial least squares regression (Hajizadeh et al., 2017a). The proposed PIC estimator is able to capture the inter- and intrasubject variabilities to provide accurate information on the amount of insulin present in the body, which is used to determine the constraints on the MPC formulation.
and
u u u θk−p θk−p+1 · · · θk−1 .. u u θk−p . θk−2 0 ΓL = . .. . .. ··· ··· u 0 0 · · · θk−f
y y y θk−p θk−p+1 · · · θk−1 .. y y θk−p . θk−2 0 ΓK = . .. . .. ··· ··· y 0 0 · · · θk−f
where f is the future window length. 2.2 Recursive Subspace-Based System Identification A stable, reliable, and computationally tractable dynamic model is essential for the design of model-based predictive control algorithms in AP systems. To this end, the predictor-based subspace identification approach is implemented to track the time-varying linear system and is coupled with a constrained optimization solver to guarantee the stability of the model (Houtzager et al., 2012; Hajizadeh et al., 2017b, 2018a). Consider a vector autoregressive with exogenous variables (VARX) model yˆk|k−1 =
p
u θk−i uk−i +
i=0
p
y θk−i yk−i
(6)
i=1
where yˆk|k−1 is the predicted output for the kth sampling instance using the past inputs uk , . . . , uk−p and outputs yk−1 , . . . , yk−p . The VARX model parameter p is the length of the past window of data considered when predicting the future outputs. The coefficient matrices θu and θy are readily estimated though recursive least square (RLS) techniques at each sampling time. Furthermore, the stacked vector yk−p,p is defined with respect to the past T T T T window of length p as yk−p,p = yk−p yk−p+1 . . . yk−1 . Stacked vector uk−p,p is also similarly defined. Recognizing that the predicted state x ˆk is given by p
x ˆk = A x ˆk−p + Luk−p,p + Kyk−p,p
(7)
x ˆk = Luk−p,p + Kyk−p,p
(8)
Γx ˆk = Γ Luk−p,p + Γ Kyk−p,p
(9)
where Land K denote the extended controllability matri- ces L = Ap−1 B . . . AB B and K = Ap−1 K . . . AK K and assuming that the state transition matrix is nilpotent with degree p, that is the contribution of the initial state x ˆk−p is negligible for sufficiently large p (Ap = 0), the predicted state can be expressed as Pre-multiplying the predicted state by the observability matrix Γ gives with
Γ =
C CA .. . CAp−1
The product of the matrices Γ L and Γ K can be constructed from the VARX model coefficient matrices as 453
Therefore, after estimating the VARX coefficient matrices, the estimated coefficient matrices θu and θy can be used to determine all quantities on the right-hand side of eq. 9, and a SVD factorization can be used to readily obtain a low-rank approximation of the state sequence. For recursive identification, a selection matrix S of appropriate dimensions can be determined such that the basis of the state estimation is consistent at each sampling time as (10) x ˆk = Sk Wk Γ Luk−p,p + Γ Kyk−p,p where Wk is a pre-defined weight matrix and the selection matrix S can be recursively updated through the projection approximation subspace tracking (PAST) method (Chiuso, 2007). The estimated state sequence is then employed along with the inputs and measured outputs to estimate the system matrices by solution of RLS problems. Specifically, after computing an estimate of the state sequence x ˆk , two RLS problems that ensure stability of the estimated system are used to determine the state-space matrices, thus yielding the identified model x ˆk+1 = Ak x ˆ k + B k u k + Kk e k (11) yk = C k x ˆ k + D k uk + e k where A, B, C, and D are the system matrices, K is the Kalman gain matrix, and ek = yk − yˆk for feedback correcting the state variables. The data-driven model determined from the subspacebased system identification technique is integrated with a physiological insulin compartment model derived from Hovorka’s model (Hajizadeh et al., 2018b). Considering a delay of order d and assuming no direct feedthrough term in the state-space model (that is Dk = 0), we can write ˆk + Bk Ik−d + Kk ek x ˆk+1 = Ak x (12) yk = C k x ˆ k + ek The PIC variable in Eq. 12, denoted I, is a filtered insulin to handle time-delays and discrete or impulse variations, which can be obtained from the compartment model with injected basal and bolus insulin as the input variable. To this end, the compartment model in Eq. 1 is integrated with the recursively identified model Eq. 12 yielding ˜ k uk ˜ˆk + B x ˜ˆk+1 = A˜k x (13) ˆ˜k yˆk = C˜k x with the new state vector of T ˆ˜k = x (14) x ˆTk+3+d Ik+2 and the output
yˆk = yˆk+3+d .
(15)
2018 IFAC NMPC Madison, WI, USA, August 19-22, 2018
Iman Hajizadeh et al. / IFAC PapersOnLine 51-20 (2018) 368–373
3. ADAPTIVE MPC ALGORITHM In this section, the glycemic and plasma insulin risk indexes, PIC bounds, the feature extraction for manipulating constraints and the adaptive MPC formulation are presented. Glycemic Risk Index A glycemic risk index (GRI) is used to determine the weighting matrix, denoted Q (ˆ yk+d+3 ), for penalizing the deviations of the outputs from their nominal set-point (Rashid et al., 2018). The glycemic risk index asymmetrically increases the set-point tracking weight when the off-line predicted CGM yˆk+d+3 diverges from the target range. Plasma Insulin Risk Index A plasma insulin penalty index (PIRI), denoted γk+d+3 , is defined to manipulate the weighting matrix for penalizing the amount of input actuation (aggressiveness of insulin dosing) depending on the estimated PIC, thus suppressing the infusion rate if sufficient insulin is present in the bloodstream (Rashid et al., 2018). To this end, the time-varying positive definite weighting matrix Rk := R (γk+d+3 ) is developed from the PIRI as R (γk+d+3 ) := γk+d+3 · R, with R as a nominal weight and γk+d+3 defined as 2 PICk+d+3 γk+d+3 := (16) PICbasal,k+d+3 where
Idb,k+d+3 (17) VI · ke,k and Idb is the patient specific (possibly time-varying) basal insulin rate that is known in practice, and VI and ke are parameters of Hovorka’s model. PICbasal,k+d+3 :=
Plasma Insulin Concentration Bounds In the proposed MPC, the estimated future PIC is dynamically bounded depending on the value of the CGM measurements. For instance, if the CGM measurement values are elevated, the bounds on the PIC are increased to ensure sufficient insulin is administered to regulate the glucose concentration. Furthermore, the PIC bounds also constrain the search space in the optimization problem. The PIC bounds are determined based on the CGM measurements as XPIC := yk ), where XPIC defines the lower and upper (1 + Pmeal )X (ˆ bounds and a desired target for the normalized PIC through the estimated CGM yˆk and Pmeal is a parameter which modifies the bounds when there is a rapid increase in the CGM measurements. 3.1 Feature Extraction for Manipulating Constraints One of the main disturbances that an AP system faces is meal consumption. The blood glucose rises drastically after a meal and it is necessary to counteract this significant disturbance with the required amount of insulin infusion in a correct time. In order to detect meal consumption and give the insulin dose in a timely manner, several meal detection algorithms are developed (Samadi et al., 2017). In this work, we build upon these to handle the effects of meals and other disturbances that result in a rapid and persistent increase in the blood glucose. The new approach modifies the constraints if there is any rapid and significant 454
371
increase in the glucose measurements. Meal detection can be based on the qualitative description of segments of glucose time series data (Samadi et al., 2017). In this approach, features from the glucose measurement data are extracted using the first and second derivatives of a secondorder polynomial fitted to the last few measurement data. If the first-order derivative is positive and greater than a threshold while taking into consideration the sign of second derivative, the constraints are modified to make the controller more aggressive to suggest a sufficient insulin dose. At every sampling time, the approach estimates the coefficients of the following polynomial function Yk = pn1 × tk 2 + pn2 × tk 1 + pn3 .
(18)
where pn1 , pn2 and pn3 are the parameters of a secondorder polynomial fitted to the last few glucose measurement and tk is the sampling time. At each sampling instance, by updating Eq. 18, the Pmeal parameter is calculated by the first and second derivatives as follows: Yk/T dm,1 if Y ≥ T d m,11 , Yk ≥ 0 k Pmeal = (19) Yk/T dm,2 if Y ≥ T d m,11 , Yk < 0 k
where the Yk and Yk are the first and second order derivatives at the last sampling time, the T dm,1 , T dm,2 and T dm,11 are patient-specific threshold parameters. 3.2 Adaptive MPC Formulation In this subsection, we propose a novel adaptive MPC algorithm cognizant of the PIC for computing the optimal insulin infusion rate. The proposed MPC formulation employs the glycemic and PIC risk indexes that manipulate the penalty weighting matrices in the cost function. To this end, the MPC computes the optimal insulin infusion over a finite horizon using the identified time-varying subspacebased models by solving at each ith sampling instance the following quadratic programming problem nP P := arg min JnP ,i yˆi+d+3 , γk+d+3 , {uk }nk=0 {u∗k }k=0 u∈U ˆ˜k+1 = Ak x ˆ˜k + Bk uk−d−3 , ∀k ∈ ZnP −1 x 0 ˆ˜k , ∀k ∈ ZnP s.t. yˆk = Ck x 0 x ˆ˜0 = x ˆ˜i x ∈ X (20) with the objective function nP nP JnP ,i (·) := eTk Qi ek + uTk Ri uk k=0
k=0
ˆ˜ and yˆ denote the predicted states and outwhere x puts, respectively, for the prediction/control horizon nP , u ∈ Rm denotes the vector of constrained input variables, taking values in a nonempty convex set U ⊆ Rm with U := {u ∈ Rm : umin ≤ u ≤ umax }, umin ∈ Rm and umax ∈ Rm denote the lower and upper bounds on the manipulated input, respectively, and ek := yˆk − ysp . The index Zn0 P represents all integers in a set as Zn0 P := n ˜ {0, . . . , n˜nP }. The nonemptyconvex set Xn˜ ⊆ R with X := x ∈ R : xmin ≤ x ≤ xmax , xmin ∈ R and xmax ∈ Rn˜ denote the lower and upper bounds on the state variables, respectively, with one of the states as the estimated PIC that is constrained through the PIC bounds. Furthermore, ˆ˜i provides an initialization of the state vector, Q ≥ 0, x
2018 IFAC NMPC 372 Madison, WI, USA, August 19-22, 2018
Iman Hajizadeh et al. / IFAC PapersOnLine 51-20 (2018) 368–373
Table 1. Closed-loop Simulation Results Subject Adult Adult Adult Adult Adult Adult Adult Adult Adult Adult
1 2 3 4 5 6 7 8 9 10
Average
< 55
[55, 70)
> 250
Mean
Statistics SD Min
Max
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
Percent of time in range [80, 140] [70, 180) [180, 250) 43.5 41.0 59.6 38.3 43.4 51.2 44.0 39.5 43.7 48.4
89.7 65.9 80.5 61.9 61.6 73.2 76.6 71.4 83.2 82.2
10.3 34.1 19.5 38.1 35.2 26.8 23.3 28.6 16.8 17.8
0.0 0.0 0.0 0.0 3.3 0.0 0.1 0.0 0.0 0.0
138.9 150.4 137.2 153.9 163.4 144.7 150.1 150.6 141.3 141.0
35.1 46.1 37.7 45.9 47.4 44.3 41.9 42.0 36.6 37.4
83.0 83.0 76.0 79.0 74.0 80.0 98.0 90.0 87.0 85.0
214.0 236.0 224.0 233.0 287.0 236.0 252.0 236.0 218.0 233.0
0.0
0.0
45.3
74.6
25.0
0.3
147.1
41.4
83.5
236.9
Qk := Q (ˆ yk+d+3 ) is a positive semi-definite symmetric matrix used to penalize the deviations of the outputs from their nominal set-point, and R > 0, Rk := R (γk+d+3 ) is a strictly positive definite symmetric matrix to penalize the manipulated input variables. 4. RESULTS The efficacy of the proposed PIC cognizant MPC with the CGM feature extraction method is demonstrated using the UVa/Padova metabolic simulator. Overall 10 subjects are simulated for 12 days with varying times and quantities of meals consumed on each day. The times and amounts of meals are selected from specified ranges to evaluate the proposed MPC under real challenging realistic scenarios. The meal information is not utilized in the MPC algorithm as the controller is designed to regulate the (blood glucose concentration) BGC in the presence of significant disturbances such as unannounced meals. In order to handle the meal effects efficiently, the CGM feature extraction method is designed to make the MPC more aggressive. The controller set-point is specified to be 110 mg/dL. The model is designed with a delay of order d = 1 with regards to the effect of PIC on CGM. The quantitative evaluation of the closed-loop results based on the proposed PIC cognizant MPC algorithm using the CGM feature extraction method are presented in Table 1, which gives the percentage of samples in defined glycemic ranges and selected statistics for the glucose measurements. It is readily observed that no hypoglycemia occurs as the BGC is never below 70 mg/dL. The average percentage of time spent in the target range (BGC ∈ [70, 180] mg/dL) and the higher tier of BGC ∈ [180, 250] mg/dL are 74.6% and 25.0%, respectively. The BGC is thus tightly controlled to be within the safe range for a significant duration of time (Fig. 1). Although the BGC is tightly controlled to be within the safe range in both cases (Fig. 2), the results with the MPC that uses the feature extraction method show a superior performance especially during day time when the AP system is challenged by meals. The average percentage of time spent in the target range (BGC ∈ [70, 180] mg/dL) is improved by 12% and the mean and standard deviation of BGC are reduced by 10 and 5 mg/dL, respectively. Furthermore, the average maximum value of BGC is reduced by 17 mg/dL without any significant change in the average minimum value of BGC (a decrease of 0.5 mg/dL) and the injected insulin (manipulated variable) (an increase of 0.5 Unit of insulin per day). 5. DISCUSSION In this work, a physiological model drawn from the insulin compartment of Hovorka’s model is incorporated with a 455
Fig. 1. Closed-loop results for 10 adult subjects using the feature extraction method
Fig. 2. Comparison of MPC results without (gray) and with (white) the feature extraction method. The bottom and top of the boxes are the first and third quartiles and the line inside the box is the median. The whiskers ends represents minimum and maximum values and + indicates mean values. data-driven model developed using a recursive subspace identification technique. The subspace identification technique with recursive modeling renders the identified models valid over a diverse range of daily activities without requiring onerous and obscure information such as the amount of carbohydrate consumption. In addition to the PIC constraints, in order to better compensate for the effects of meal consumption, the rate of change of the glucose values (as determined by the derivatives of the CGM measurements) are computed in real time to modify the aggressiveness of the controller and give the required quantity of insulin in an appropriate amount of time. The proposed approach involving adaptation of the model, the penalty weights of the optimization cost function, and the constraints is necessary as the nonlinear and transient operating behavior of biological systems is difficult to characterize precisely. The adaptive models can well characterize the nonlinear and time-varying dynamics across various operating conditions by recursively updating the identified
2018 IFAC NMPC Madison, WI, USA, August 19-22, 2018
Iman Hajizadeh et al. / IFAC PapersOnLine 51-20 (2018) 368–373
models, thereby tailoring the models to the particular conditions of any given situation. This adaptive nature of the proposed technique results in good control performance, as demonstrated through the tight glucose control using the metabolic simulator. 6. CONCLUSION An adaptive MPC is proposed in this work for optimal insulin dosing decisions in AP systems using a recursive subspace-based system identification approach to characterize the transient dynamics of processes involved in diabetes. Subsequent to system identification, an adaptive MPC algorithm is designed using the recursively identified models to effectively compute the optimal insulin delivery for AP systems. A feature extraction method based on glucose measurements is used to detect rapid deviations from the desired set-point caused by significant disturbances and subsequently modify the constraints of the optimization problem for negotiating between the aggressiveness and robustness of the controller to suggest the required amount of insulin. REFERENCES Bernardini, D. and Bemporad, A. (2012). Stabilizing model predictive control of stochastic constrained linear systems. IEEE Trans Automat Contr, 57(6), 1468–1480. Chakrabarty, A., Zavitsanou, S., Doyle, F.J., and Dassau, E. (2018). Event-triggered model predictive control for embedded artificial pancreas systems. IEEE Trans. Biomed. Eng., 65(3), 575–586. Chiuso, A. (2007). The role of vector autoregressive modeling in predictor-based subspace identification. Automatica, 43(6), 1034–1048. Dalla Man, C., Micheletto, F., Lv, D., Breton, M., Kovatchev, B., and Cobelli, C. (2014). The UVa/Padova type 1 diabetes simulator: new features. J. Diabetes Sci. Technol., 8(1), 26–34. Doyle III, F.J., Bequette, B.W., Middleton, R., Ogunnaike, B., Paden, B., Parker, R.S., and Vidyasagar, M. (2011). Control in biological systems. In T. Samad, & A. M. Annaswamy (Eds.), The impact of control technology, part 1. IEEE Control Systems Society. Gondhalekar, R., Dassau, E., and Doyle, F.J. (2016). Periodic zone-MPC with asymmetric costs for outpatientready safety of an artificial pancreas to treat type 1 diabetes. Automatica, 71, 237 – 246. Hajizadeh, I., Rashid, M., and Cinar, A. (2018a). Ensuring stability and fidelity of recursively identified controlrelevant models. In The 18th IFAC Symposium on System Identification (SYSID), 927–932. Hajizadeh, I., Rashid, M., and Cinar, A. (2018b). Integrating compartment models with recursive system identification. In American Control Conference (ACC), 3583–3588. Hajizadeh, I., Rashid, M., Turksoy, K., Samadi, S., Feng, J., Frantz, N., Sevil, M., Cengiz, E., and Cinar, A. (2017a). Plasma insulin estimation in people with type 1 diabetes mellitus. Ind. Eng. Chem. Res., 56(35), 9846– 9857. Hajizadeh, I., Rashid, M., Turksoy, K., Samadi, S., Feng, J., Sevili, M., Frantz, N., Lazaro, C., Maloney, Z., Littlejohn, E., and Cinar, A. (2017b). Multivariable 456
373
recursive subspace identification with application to artificial pancreas systems. IFAC-PapersOnLine, 909 – 914. Houtzager, I., van Wingerden, J.W., and Verhaegen, M. (2012). Recursive predictor-based subspace identification with application to the real-time closed-loop tracking of flutter. IEEE Trans. Control Syst. Technol., 20(4), 934–949. 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. Physiol. Meas., 25(4), 905. Jerez, J.L., Goulart, P.J., Richter, S., Constantinides, G.A., Kerrigan, E.C., and Morari, M. (2014). Embedded online optimization for model predictive control at megahertz rates. IEEE Trans Automat Contr, 59(12), 3238–3251. Kol˚ as, S., Foss, B., and Schei, T. (2009). Constrained nonlinear state estimation based on the UKF approach. Comput. Chem. Eng., 33(8), 1386–1401. Kouramas, K., Panos, C., Fasca, N., and Pistikopoulos, E. (2013). An algorithm for robust explicit/multiparametric model predictive control. Automatica, 49(2), 381 – 389. Mayne, D.Q. (2014). Model predictive control: Recent developments and future promise. Automatica, 50(12), 2967 – 2986. Messori, M., Incremona, G.P., Cobelli, C., and Magni, L. (2018). Individualized model predictive control for the artificial pancreas: In silico evaluation of closed-loop glucose control. IEEE Control Syst. Mag., 38(1), 86–104. Muske, K.R. and Badgwell, T.A. (2002). Disturbance modeling for offset-free linear model predictive control. J. Process Control, 12(5), 617 – 632. Pannocchia, G., Rawlings, J.B., and Wright, S.J. (2007). Fast, large-scale model predictive control by partial enumeration. Automatica, 43(5), 852 – 860. Rashid, M., Hajizadeh, I., and Cinar, A. (2018). Plasma insulin cognizant predictive control for artificial pancreas. In American Control Conference (ACC), 3589–3594. Rawlings, J.B. and Mayne, D.Q. (2009). Model predictive control: Theory and design. Nob Hill Pub. Rivotti, P., Lambert, R.S., and Pistikopoulos, E.N. (2012). Combined model approximation techniques and multiparametric programming for explicit nonlinear model predictive control. Comput. Chem. Eng., 42, 277 – 287. Samadi, S., Turksoy, K., Hajizadeh, I., Feng, J., Sevil, M., and Cinar, A. (2017). Meal detection and carbohydrate estimation using continuous glucose sensor data. IEEE J. Biomed. Health Inform., 21(3), 619–627. Turksoy, K., Hajizadeh, I., Samadi, S., Feng, J., Sevil, M., Park, M., Quinn, L., Littlejohn, E., and Cinar, A. (2017). Real-time insulin bolusing for unannounced meals with artificial pancreas. Control Eng. Pract., 59, 159–164. Wang, Y. and Boyd, S. (2010). Fast model predictive control using online optimization. IEEE Trans Control Syst Technol, 18(2), 267–278. Zavala, V.M. and Biegler, L.T. (2009). The advanced-step NMPC controller: Optimality, stability and robustness. Automatica, 45(1), 86 – 93.