INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN BIOMEDICAL ENGINEERING Int. J. Numer. Meth. Biomed. Engng. (2012) Published online in Wiley Online Library (wileyonlinelibrary.com). DOI: 10.1002/cnm.2466
A genetic algorithm tuned optimal controller for glucose regulation in type 1 diabetic subjects Subhojit Ghosh 1, * ,† and Srihari Gude 2 1 Department 2 Department
of Electrical Engineering, National Institute of Technology, Rourkela, India 769008 of Electrical and Electronics Engineering, Birla Institute of Technology, Deoghar, India 815142
SUMMARY An optimal state feedback controller is designed with the objective of minimizing the elevated glucose levels caused by meal intake in Type 1 diabetic subjects, by the minimal infusion of insulin. The states for the controller based on linear quadratic regulator theory are estimated from noisy data using Kalman filter. The controller designed for a physiological relevant mathematical model is coupled with another model for simulating meal dynamics, which converts meal intake into glucose appearance rate in the plasma. The tuning parameters (weighting matrices) of the controller and the design parameters (noise covariance matrices) of the Kalman filter are optimized using genetic algorithm. The controller based on the combined framework of evolutionary computing and state estimated linear quadratic regulator is found to maintain normoglycemia for meal intakes of varying carbohydrate content. The proposed approach addresses noisy output measurement, modeling error and delay in sensor measurement. Copyright © 2012 John Wiley & Sons, Ltd. Received 27 May 2011; Revised 18 November 2011; Accepted 2 January 2012 KEY WORDS:
diabetes; glucose regulation; insulin infusion; linear quadratic regulation; Kalman filter; genetic algorithm
1. INTRODUCTION Type 1 diabetes mellitus is a metabolic disorder characterized by the retarded or complete absence of insulin secretion by the pancreas. Improper insulin secretion results in glucose levels above the normal range of 80–120 mg/dL, which leads to complications like hypertension, coma, and even death. With diabetes reaching epidemic proportions, recent times have witnessed an increased attention in the development of model-based control algorithms for the maintenance of normoglycemia in diabetic subjects by the external infusion of insulin [1–4]. The algorithms based on the mathematical model of glucose–insulin dynamics, provide an idea of a suitable insulin infusion profile, which when used in conjunction with a proper drug delivery system can maintain glucose levels within desired limits. Along with the type of mathematical model and the control strategy employed, the control algorithms reported in literature differ on their objectives and route of insulin infusion, i.e., intravenous [1] and subcutaneous [2]. The notable control techniques developed in this regard include the conventional proportional-integral-derivative (PID) control [5–7], H1 based robust control [8–10], model predictive control (MPC) [11–14] and linear quadratic regulator (LQR) [15, 16]. Although a majority of these algorithms are for clinical study and validated on mathematical models, many techniques have also been validated in real clinical settings [4, 7, 12, 17–19]. MPC and LQR have proved to be important tools among the state feedback dependent optimal control techniques. These techniques employ weighting matrices for minimizing glucose level and insulin infusion.
*Correspondence to: Subhojit Ghosh, Electrical Engineering, NIT Rourkela, India 769008. † E-mail: [email protected] Copyright © 2012 John Wiley & Sons, Ltd.
S. GHOSH AND S. GUDE
For an observer based state feedback control of a plant corrupted by state and measurement noise, the control action and the appropriateness of the estimated states is heavily dependent on the output and control weighting matrices along with the process and noise covariance matrices. The selection of these parameters is not a trivial problem, and is hence carried out by trial and error. This involves maintaining a trade-off between minimizing the control effort and improving the transient response. An increased weightage towards minimizing glucose level leads to high insulin infusion rate and vice-versa for attaching greater weightage for control effort. Similarly, the covariance matrices also affects the overall control objective by affecting the error dynamics of the observed state vector. In this regard, a genetic algorithm (GA) tuned Kalman filter based LQR theory has been proposed for blood glucose regulation in diabetic subjects in case of a disturbance in the form of a meal. GA [20], which has proved to be a powerful evolutionary technique in recent times, has been utilized for the selection of the weighting matrices of the LQR controller and the noise covariance matrices of the Kalman filter. Applications of genetic algorithms in glucose metabolism have been reported for parameter estimation [21, 22] and quantification of insulin sensitivity [23] from glucose insulin data. For a given patient model, the algorithm searches for the best combination of the tuning parameters that results in the least deviation of glucose from an acceptable basal value. Simulation studies conducted on a physiologically relevant model depicts the effectiveness of the proposed controller in maintaining normal glucose regulation in type 1 diabetes mellitus subjects, using both instantaneous as well delayed glucose measurements. The paper is organized as follows: the next section focuses on a mathematical model of glucose regulation, the Kalman filter based linear quadratic regulator for glucose control is discussed in Section 3, Section 4 gives a brief overview of GA along with its application for tuning the controller and Kalman filter parameters, simulation results for analyzing the controller performance are given in Section 5, and finally Section 6 provides conclusions. 2. MATHEMATICAL MODEL FOR GLUCOSE REGULATION For controller design, the mathematical model proposed by Sturis et al. [24] is considered in the present work. The nonlinear model takes following form: dIp dt dIi dt dG dt dI1 dt dI2 dt dI3 dt
Ip Ii / Vi Vp Ii Ip / D E. Vp Vi D E.
Ip C Uin .t / tp Ii ti
D Din .t / f2 .G/ f3 .G/f4 .Ii / C f5 .I3 / 3 .Ip I1 / td 3 D .I1 I2 / td 3 D .I2 I3 / td
(1)
D
For a normal subject, the following feedback loops are included in the model: pancreatic insulin secretion by stimulation of glucose, insulin stimulates glucose uptake and inhibits hepatic glucose production and glucose enhances its own uptake. However, because the present work involves designing a controller for type 1 diabetic subject, the pancreatic insulin secretion is not considered. The model has three main physiological variables; the concentration of glucose in the plasma G(mg/dL), the concentration of insulin in the plasma Ip (mU/L)and the insulin concentration in the intracellular space Ii (mU/L). In addition to the physiological variables, the model also contains three auxiliary insulin variables I1 , I2 , and I3 . The distribution volume (L) in the plasma insulin, intracellular insulin, and plasma glucose space is represented by Vp , Vi , and Vg , respectively, E(L/min) is the diffusion transfer rate, td (min) is the delay in hepatic glucose production, Ui n .t / Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
GA TUNED OPTIMAL CONTROLLER FOR GLUCOSE REGULATION IN TYPE 1 DIABETES
is the insulin infusion to the plasma compartment, tp (min) and ti (min) are the insulin degradation time constants in the plasma and intercellular space. Functions f2 .G/ and f3 .G/f4 .Ii / represent the insulin independent and insulin dependent glucose utilization, respectively, whereas, function f5 .I3 / represents hepatic glucose production. The controller designed for the above patient model aims to compensate the increase in glucose level caused by the disturbance Di n .t / (mg/dL min1 / by the application of insulin infusion (control input) . In comparison with the intravenous insulin infusion, subcutaneous delivery is more realistic in clinical setting. Hence, it is necessary to augment the above model with the dynamics representing the passage of insulin from the subcutaneous tissue to the plasma insulin compartment. For this purpose, the following model [25] of subcutaneous insulin kinetics has been considered dS1 D k21 S1 C U.t / dt dS2 D k21 S1 .kd C ka /S2 dt
(1a)
where U.t / is the rate of insulin infusion, S1 is the subcutaneous mass where injection takes place, S2 is the subcutaneous insulin mass proximal to the plasma, kd and ka are the degradation constants in the subcutaneous tissue and plasma, respectively, and k21 is the rate constant describing insulin transport. The term ka S2 =Vp corresponds to the rate of insulin infusion to the plasma space from subcutaneous injection. To arrive at a controller for glucose regulation, the integrated model is linearized about the steady state point (basal value) and a new set of variables, i.e., their deviation from the basal values are considered. 3. LINEAR QUADRATIC REGULATOR OPTIMAL FEEDBACK In this section, we briefly explain the observer (Kalman filter) based LQR method [26, 27]. In the present work, the LQR controller is formulated by considering subcutaneous insulin infusion as the control input and glucose appearance rate in the plasma from meal intake as disturbance. The linearized regulator model for controlling glucose levels using state feedback subjected to meal disturbance, can be written in the state space form as dx D Ax C Bu C din dt y D Cx
(2)
where A, B, and C are system matrices and the state vector x is given by x D ŒIp , Ii , G, I1 , I2 , I3 , S1 , S2 , with representing the deviation of the states of (1), from their basal values. The control input is given by u and di n represents the disturbance (di n = [0 0 Di n 0 0 0 0 0]). The LQR method involves the determination of the optimal control input u D kx, that minimizes the following performance index: Z1 JLQR D
.x T Qx C uT Ru/dt
(3)
0
where Q and R are positive definite matrices. These matrices denote the importance attached to the state and control output in the minimization of J . The controller gain K, which minimizes (3), is given by K D R1 B T P ,
(4)
where P is the positive symmetric solution of the following algebraic Riccatti equation PA C AT P C Q PBR1 B T P D 0. Copyright © 2012 John Wiley & Sons, Ltd.
(5)
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
S. GHOSH AND S. GUDE
The selection of Q and R significantly affects the solution of (5) and hence the controller performance. In the conventional approach, Q and R are selected from a number of trial and error iterations. The application of the control input u involves the availability of all the states at the output. However, the online measurement of plasma and intracellular insulin requires a highly invasive and expensive protocol. Also, the auxiliary insulin variables incorporated in the model (1) for representing a delay are nonphysiological in nature. Continuous monitoring of plasma glucose levels, through the subcutaneous route has been made possible by the developments in the field of implantable glucose sensor technology. Hence, of all the states in (2), for real-time control only G is considered to be available in the output. Thus, for applying state feedback control for glucose regulation, an observer is required for the estimation of nonaccessible and nonphysiological state variables. Generalizing the model (2), to include a wider class of systems by including modeling error and measurement noise leads to following representation: dx dt
D Ax C Bu C din C F v y D Cx C ´
(6)
where v is the process noise, which may occur because of modeling error such as noninclusion of nonlinearities and high frequency dynamics, and ´ is the measurement noise. The noises are independent of each other, white with normal probability distribution. In addition, all the state variables are assumed to be affected by the process noise and a single value is considered for the state equation. Kalman filter estimates the state variables.x/of Q the above stochastic plant from the input and noisy output measurement by minimizing the covariance of the conditional error between the actual and the estimated state, given as R D EŒ.x x/ Q T .x x/. Q
(7)
For a predefined noise covariance (Z for measurement noise and V for process noise), the estimated state equation of the Kalman filter is given as d xQ D AxQ C Bu C L.y C x/, Q dt
(8)
where the Kalman gain vector L is given as L D RC T Z 1 .
(9)
The solution of the matrix R is given by the following algebraic Riccatti equation AR C RAT RC T Z 1 CR C F VF T D 0
(10)
Similar to the weighting matrices Q and R for the LQR controller, the covariance matrices Z and V can be considered as tuning parameters for the Kalman filter. Their selection affects the dynamics of the error between actual and estimated state. Since in most cases the power spectral density of the noises are not known in advance, Zand V are also chosen by trial and error. Hence, the two-stage design of an optimal controller for glucose regulation from estimated states involves the design of an optimal regulator that minimizes (3) and a Kalman filter that minimizes (7). The linearity of the system and considering the noises as inputs, allows separating the controller design problem into two isolated processes of estimation and control. The regulator and Kalman filter can be designed separately and combined to form an optimal compensator. The conventional approach of designing a Kalman filter based LQR controller, involves a large number of trial and error iterations for choosing the best combination of the tuning parameters Q, R, Z, and F . The large computational time retards the application of this approach for large-scale systems. In the next section, a technique based on GA is proposed to tune the parameters of the LQR controller and Kalman filter for improving the integrated control action. Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
GA TUNED OPTIMAL CONTROLLER FOR GLUCOSE REGULATION IN TYPE 1 DIABETES
4. GENETIC ALGORITHM BASED CONTROLLER TUNING 4.1. Genetic algorithm Genetic algorithm is an effective search method based on natural evolution with the Darwinian survival of the fittest approach. Unlike conventional optimization algorithms, GA works on coded parameters and not on actual parameters. It does not involve derivation of the objective function and can hence work on discontinuous functions. This advantage allows GA to solve complex multivariable problems, which is difficult or impossible to be solved by traditional methods. GA involves three major operations of reproduction, crossover, and mutation. For an initially randomly selected possible solutions (population) coded as binary bits, during reproduction individual solutions known as chromosomes are selected based on their fitness, and sent to the mating pool. During crossover operation, two strings selected at random from the mating pool undergo crossover with a certain probability at a randomly selected crossover point to generate two new strings. Depending on whether a randomly generated number is larger than a predefined mutation probability or not each bit in the string obtained after crossover is altered (changing 0 to 1 and 1 to 0). The three operations are carried out for a certain number of iterations (generations), or till a chromosome with a certain fitness is obtained. 4.2. Genetic algorithm based linear quadratic regulator controller design S for the present problem, the aim is to control glucose level only, the diagonal elements of matrix Q in (3) for all the variables other than plasma glucose has been constrained to zero. Also, the same covariance has been considered for all the process noises. This results in a total of four tuning
Figure 1. Flowchart for the GA tuned LQR controller. Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
S. GHOSH AND S. GUDE
parameters (one each for Q, R, Z, and V /. GA selects the tuning parameters, so that the controller with the corresponding gains (K and L/ is able to maintain normal glucose levels, in case of a meal disturbance. The detailed procedure for obtaining the parameters is shown in Figure 1. The tuning parameters are encoded into chromosomes of 40 bits (10 bits each for each parameter). The chromosome bits corresponding to Q and R are selected based on the minimization of the following fitness function: Ztf " JRegulat or D 0
G.t / Gbasal
2
U.t / C Ubasal
2 # dt ,
(11)
where, Gbasal and Ubasal represent the basal glucose level and insulin infusion rate, respectively. For a given Q and R, the fitness function is calculated by simulating the closed loop patient model for a given meal disturbance, until the predisturbance glucose level is attained. As mentioned earlier, the selection of Q and R has a direct bearing on the response and magnitude of the control input of the closed loop system. Selecting an appropriate state-weighting matrix Q imposes certain limits on the maximum overshoot of plasma glucose. Similarly, the maximum overshoot limit in the insulin
Figure 2. GA tuned LQR controller response for 50 g meal disturbance. (a) Plasma glucose concentration and (b) insulin infusion rate. Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
GA TUNED OPTIMAL CONTROLLER FOR GLUCOSE REGULATION IN TYPE 1 DIABETES
infusion rate is specified through R. For the Kalman filter, the chromosome bits are selected based on the following objective function: JKalman D
n n X X
Re2 .i, j /.
(12)
i D1 j D1
For each chromosome, the feedback gain K, is determined from the solution of (5), whereas the Kalman gain L, is obtained from the solution of (10). The gain factors are used to calculate JRegulat or and JKalman from the controller response. 5. RESULTS AND DISCUSSION The performance of the GA tuned LQR controller in compensating a meal disturbance from the feedback of noisy glucose measurement is considered in this section. Intravenous route is considered for glucose measurement. The controller performance is evaluated based on the excursion of the plasma glucose level from the basal value and the amount of insulin infused to compensate the disturbance. For simulating the meal dynamics of certain carbohydrate content, the bicompartmental
Figure 3. GA tuned LQR controller response for meal disturbances of 20 g (breakfast), 50 g (lunch) and 40 g (dinner) in a day. (a) Plasma glucose concentration and (b) insulin infusion rate. Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
S. GHOSH AND S. GUDE
model proposed by Lehmann and Deutsch in [28] has been considered. In the model, the conversion from meal intake (in terms of carbohydrate content) to glucose absorption into the plasma compartment is expressed as a trapezoid wave followed by a first-order linear system (Appendix). For tuning the GA, a large number of selection and crossover strategies were tested through pilot runs to determine the most suitable ones. Single point crossover is used in conjunction with elitist strategy based Roulette wheel selection. The elitist strategy, which is able to preserve superior strings, is incorporated in the present work by replacing the 10 worst strings of a particular generation with the 10 best strings of the previous generation. The number of chromosomes in the initial population and the maximum number of generations is set at 50 and 200 in each run. The crossover and mutation probability was fixed at 0.8 and 0.2, respectively. The algorithm was found to converge to the same solution for different intialization of the tuning parameters. The performance of the controller (in terms of plasma glucose level and insulin infusion rate) with the GA derived tuning parameters (Q(3,3)=7.294, R=0.015, V =6.158, Z=0.023), following a meal disturbance of 50 gm carbohydrate content is shown in Figure 2. Because there is no endogenous secretion of insulin, a continuos insulin infusion should always be present even in the absence of a disturbance, to maintain the basal glucose level. The closed loop performance of the diabetic patient for the entire day subjected to breakfast at 8 AM of 20 g, lunch at 1 PM of 50 g and dinner at 8 PM
Figure 4. GA tuned LQR controller response for 50 g meal disturbance using fitness function with glucose deviation only. (a) Plasma glucose concentration and (b) insulin infusion rate. Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
GA TUNED OPTIMAL CONTROLLER FOR GLUCOSE REGULATION IN TYPE 1 DIABETES
of 40 g is shown in Figure 3. Similar plasma glucose and insulin infusion rate profile during the initial stage for different meal intakes is because of the similar rate of glucose appearance into the blood circulation during that period, as obtained from the meal model [28]. However, prolonged and high glucose infusion after that results in larger settling time for meals of high carbohydrate content. Similar to single meal disturbance, for multiple meals also, the insulin infusion rate contains high frequency noise components of the order of ˙0.2 mU/L min approximately. For all the disturbances considered above, the controller is able to maintain the glucose levels within the normal range of 80–120 mg/dL. It should be noted that although the controller is designed based on the linearized model, the controller action results reported here is that on the nonlinear model (1). For many cases, depending on the condition of the subject and the available protocol for insulin infusion, the primary goal might be to keep the glucose level within the normal range without bothering on the rate of insulin infusion. In this regard, it may be advantageous to have periods where insulin infusion deviates a lot from the basal value. Such cases can be handled within the framework of the present work, by not considering the minimization of insulin infusion in the fitness function (11), that is, only the first component of (11) is taken into account. On the basis of this formulation, the response of the controller following a meal disturbance of 50 gm carbohydrate content is shown in Figure 4.
Figure 5. GA tuned LQR controller response for 50 g meal disturbance with a glucose sensor delay of 2 min. (a) Plasma glucose concentration and (b) insulin infusion rate. Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
S. GHOSH AND S. GUDE
An important aspect in controller design for blood glucose regulation is the time delay in glucose sensor measurement. In this paper, until now it has been assumed that the measurement of plasma glucose level, although noisy, is available instantaneously. Along with the degradation of the controller performance, large sensor delay results in an unstable closed loop system. To evaluate the degradation in controller performance, the proposed controller was subjected to measurement delays of 2 and 5 min. The variation in plasma glucose and insulin infusion rate following a meal of 50 g for the above two measurement delays are shown in Figures 5 and 6. For both cases, the glucose levels are well within the acceptable range. However, for large time delays (> 10 min), sustained oscillations are observed in glucose levels with high peak value. For continuous glucose infusion, the standard LQR that provides only proportional control is found to produce a steady state error. In addition to minimizing plasma glucose level and insulin infusion rate, an additional objective of minimizing the integral of the plasma glucose state is required for continuous disturbance. For this, an integral action is added to the standard LQR. The model (2) is augmented to include a new state variable, that is, the integral of the deviation in glucose level from basal value, given as
Figure 6. GA tuned LQR controller response for 50 g meal disturbance with a glucose sensor delay of 5 min. (a) Plasma glucose concentration and (b) insulin infusion rate. Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
GA TUNED OPTIMAL CONTROLLER FOR GLUCOSE REGULATION IN TYPE 1 DIABETES
Zt G. /d
(13)
0
With the additional tuning parameter for the augmented state, the response of the proposed controller with the GA derived tuning parameters (Q(3,3)=0.151, Q(7,7)=0.002, R=23.083, V =0.483, Z=3.92), for continuous glucose infusion of 1 mg/dL min is shown in Figure 7. Unlike the standard LQR controller, where the variation in weighting factors do not effect the time required to reach the steady state, varying the tuning parameter corresponding to the integral of glucose level has a direct effect on the settling time, with higher value resulting in faster convergence to the basal glucose level with oscillations in the steady state. Varying the integral weighting factor is similar to the tuning of integral term in a conventional PID controller. Unlike the regulation for meal disturbance, for constant glucose infusion an insulin infusion rate in excess of the basal value of 3.2 mU/L min is required to maintain normoglycemia at steady state. Oscillations in plasma glucose levels resulting in oscillatory insulin infusion rates should be avoided because the insulin delivery system may not
Figure 7. GA tuned LQR controller response for a continuous glucose infusion of 1 mg/dL min. (a) Plasma glucose concentration and (b) insulin infusion rate. Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
S. GHOSH AND S. GUDE
be able to alter the rate of infusion at a very fast rate. The continuous glucose infusion, although not representing a realistic disturbance to the glucose regulation system, is similar to the situation arising because of changes in the physiological unit processes involved. Changes in insulin sensitivity and retarded insulin independent glucose utilization may result in an offset from the steady state value of the original model. The resulting offset requires integral action to maintain the premodified basal glucose and insulin level. 6. CONCLUSIONS Evolutionary computation has been utilized to design an optimal feedback controller for glucose regulation in type 1 diabetic subjects using noisy output data. The inaccessible and nonphysiological state variables are estimated from the noisy glucose output measurement using a Kalman filter. The selection of the design parameters for the Kalman filter and the feedback controller based on LQR is in general a trial and error process, which involves selection of the best combination from different simulated responses. In the present work, these parameters are optimized using GA. The two-stage algorithm searches for the parameters that minimize the deviation in plasma glucose levels from the basal value and the covariance of the state estimation error. The performance of the controller is evaluated for different meal disturbances with both instantaneous and delayed sensor measurement. The proposed controller is able to maintain glucose levels in type 1 diabetic subjects within the normal range in case of meal disturbances of varying carbohydrate content. The simplicity of the approach and the absence of any online calculation make it easy to store in a digital chip, which can then be embedded with a suitable drug delivery system and an in vivo sensor for automated insulin delivery. Future work in this direction is aimed at testing the efficacy of the proposed strategy using more comprehensive models and studying the robust handling capacity in case of parameter variations. APPENDIX: MEAL DYNAMICS MODEL In the model [28], the absorption of glucose from a meal consists of two stages: one is the emptying of the stomach into the duodenum and the other is the absorption of glucose from the small intestine into the blood circulation. The rate of gastric emptying of carbohydrate is considered as a nonlinear function with saturation characteristics to maintain a constant carbohydrate release from the stomach during intestinal absorption. The rise and fall from the maximum saturated level emptying rate Vm,CHO = 360 mg/min is considered as a linear function spanning over 30 min. These time periods are denoted by Tri se and Tf al l . The period of constant carbohydrate release depends on the carbohydrate content of the meal consumed. Meals with carbohydrate content less than 10.8 g do not reach the saturation level. Larger meals are represented by trapezoid wave. The glucose absorption is modeled by a first-order equation, which takes the form GP gut D kab Ggut C Gempt , Din D kab Ggut
(A1)
where, Ggut is the glucose concentration in the intestine, kab is the rate of glucose absorption into the blood circulation. REFERENCES 1. Parker RS, Doyle III FJ, Peppas NA. The Intravenous Route to Blood Glucose Control. IEEE Engineering in Medicine and Biology 2001; 20(1):65–73. 2. Bellazzi R, Nucci G, Cobelli C. The Subcutaneous Route to Insulin-Dependent Diabetes Therapy. IEEE Engineering in Medicine and Biology 2001; 20(1):54–64. 3. Wayne Bequette B. A critical assessment of algorithms and challenges in the development of a closed-loop artificial pancreas. Diabetes Technology and Therapeutics 2005; 7(1):28–47. 4. Hovorka R. Continuous glucose monitoring and closed-loop systems. Diabetic Medicine 2005; 23(1):1–12. 5. Marchetti G, Barolo M, Jovanovic L, Zisser H, Seborg DE. An improved PID switching control strategy for type 1 diabetes. IEEE Transactions on Biomedical Engineering 2008; 55(3):857–865. Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
GA TUNED OPTIMAL CONTROLLER FOR GLUCOSE REGULATION IN TYPE 1 DIABETES
6. Lam ZH, Hwang KS, Lee JY, Chase JG, Wake GC. Active insulin infusion using optimal and derivative-weighted control. Medical Engineering and Physics 2002; 24(10):663–672. 7. Steil GM, Rebrin K, Darwin C, Hariri F, Saad MF. Feasibility of automating insulin delivery for the treatment of type 1 diabetes. Diabetes 2006; 55:3344–3350. 8. Parker RS, Doyle III FJ, Ward JH, Peppas NA. Robust H1 control in diabetes using a physiological model. AIChE Journal 2000; 46(12):2537–2549. 9. Quiroz G, Femat R. Theoretical blood glucose control in hyper- and hypoglycemic and exercise scenarios by means of an H1 algorithm. Journal of Theoretical Biology 2010; 263:154–160. 10. Ruiz Velazquez E, Femat R, Campos Delgado DU. Blood glucose control for type I diabetes mellitus: a robust H1 tracking problem. Control Engineering Practice 2004; 12:1179–1195. 11. Parker RS, Doyle III FJ, Peppas NA. A model-based algorithm for blood glucose control in type I diabetic patients. IEEE Transactions on Biomedical Engineering 1999; 46(2):148–157. 12. Hovorka R, Canonico V, Chassin LJ, Haueter U, Massi-Benedetti M, Orsini Federici M, Pieber TR, Schaller HC, Schaupp L, Vering T, Wilinska ME. Nonlinear model predictive control of glucose concentration in subjects with type 1 diabetes. Physiological Measurement 2004; 25(4):905–920. 13. Kovatchev BP, Breton MD, Dalla Man C, Cobelli C. In Silico Preclinical Trials: A Proof of Concept in Closed-Loop Control of Type 1 Diabetes. Journal of Diabetes Science and Technology 2009; 3(1):44–55. 14. Abu-Rmileh A, Garcia-Gabin W. Feedforward–feedback multiple predictive controllers for glucose regulation in type 1 diabetes. Computer Methods and Programs in Biomedicine 2010; 99:113–123. 15. Patek SD, Breton MD, Chen Y, Solomon C, Kovatchev B. Linear quadratic Gaussian based closed loop control of type 1 diabetes. Journal of Diabetes Science and Technology 2007; 1(6):834–841. 16. Sanchez Chavez IY, Morales Menendez R, Martinez Chapa SO. Glucose optimal control system in diabetes treatment. Applied Mathematics and Computation 2009; 209:19–30. 17. Chase JG, Shaw GM, Wong XW, Lotz T, Lin J, Hann CE. Model-Based Glycaemic Control in Critical Care-A Review of the State of the Possible. Biomedical Signal Processing and Control 2006; 1:3–21. 18. Palerm CS, Jovanovic L, Zisser H, Doyle III FJ, Bevier WC. Prandial Insulin Dosing Using Run-to-Run ControlApplication of clinical data and medical expertise to define a suitable performance metric. Diabetes Care 2007; 30(5):1131–1136. 19. El-Khatib FH, Jiang J, Damiano ER. A feasibility study of bihormonal closed-loop blood glucose control using dual subcutaneous infusion of insulin and glucagon in ambulatory diabetic swine. Journal of Diabetes Science and Technology 2009; 3(4):789–803. 20. Goldberg DE. Genetic Algorithms in Search, Optimization, and Machine Learning. Addison-Wesley: Reading, MA, 1989. 21. Morbiducci U, Tura A, Grigioni M. Genetic algorithms for parameter estimation in mathematical modeling of glucose metabolism. Computers in Biology and Medicine 2005; 35:862–874. 22. Morbiducci U, Benedetto GD, Kuatzky-Willer A, Pacini G, Tura A. Improved usability of the minimal model of insulin sensitivity based on an automated approach and genetic algorithms for parameter estimation. Clinical Science 2007; 112:257–263. 23. Ghosh S, Maka S. Genetic Algorithm Based NARX model Identification for Evaluation of Insulin Sensitivity. Applied Soft Computing 2011; 11(1):221–226. 24. Sturis J, Polonsky KS, Mosekilde E, Van Cauter E. Computer-model for mechanisms underlying ultradian oscillations of insulin and glucose. American Journal of Physiology 1991; 260:E801–E809. 25. Shimoda S, Nishida K, Sakakida M, Nishida K, Sakakida M, Konno Y, Ichinose K, Uehara M, Nowak T, Shichiri M. Closed-loop subcutaneous insulin infusion algorithm with a short acting insulin analog for long term clinical application of a wearable artificial endocrine pancreas. Frontiers in Medical Biology and Engineering 1997; 8:197–211. 26. Naidu DS. Optimal Control Systems. CRC Press: Boca Raton, Florida, 2003. 27. Crassidis JL, Junkins JL. Optimal Estimation of Dynamic Systems,. Chapman & Hall/CRC: Boca Raton, Florida, 2004. 28. Lehmann ED, Deutsch T. A physiological model of glucose-insulin interaction in type 1 diabetes mellitus. Journal of Biomedical Engineering 1993; 14(3):235–242.
Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
A genetic algorithm tuned optimal controller for glucose regulation in type 1 diabetic subjects Subhojit Ghosh 1, * ,† and Srihari Gude 2 1 Department 2 Department
of Electrical Engineering, National Institute of Technology, Rourkela, India 769008 of Electrical and Electronics Engineering, Birla Institute of Technology, Deoghar, India 815142
SUMMARY An optimal state feedback controller is designed with the objective of minimizing the elevated glucose levels caused by meal intake in Type 1 diabetic subjects, by the minimal infusion of insulin. The states for the controller based on linear quadratic regulator theory are estimated from noisy data using Kalman filter. The controller designed for a physiological relevant mathematical model is coupled with another model for simulating meal dynamics, which converts meal intake into glucose appearance rate in the plasma. The tuning parameters (weighting matrices) of the controller and the design parameters (noise covariance matrices) of the Kalman filter are optimized using genetic algorithm. The controller based on the combined framework of evolutionary computing and state estimated linear quadratic regulator is found to maintain normoglycemia for meal intakes of varying carbohydrate content. The proposed approach addresses noisy output measurement, modeling error and delay in sensor measurement. Copyright © 2012 John Wiley & Sons, Ltd. Received 27 May 2011; Revised 18 November 2011; Accepted 2 January 2012 KEY WORDS:
diabetes; glucose regulation; insulin infusion; linear quadratic regulation; Kalman filter; genetic algorithm
1. INTRODUCTION Type 1 diabetes mellitus is a metabolic disorder characterized by the retarded or complete absence of insulin secretion by the pancreas. Improper insulin secretion results in glucose levels above the normal range of 80–120 mg/dL, which leads to complications like hypertension, coma, and even death. With diabetes reaching epidemic proportions, recent times have witnessed an increased attention in the development of model-based control algorithms for the maintenance of normoglycemia in diabetic subjects by the external infusion of insulin [1–4]. The algorithms based on the mathematical model of glucose–insulin dynamics, provide an idea of a suitable insulin infusion profile, which when used in conjunction with a proper drug delivery system can maintain glucose levels within desired limits. Along with the type of mathematical model and the control strategy employed, the control algorithms reported in literature differ on their objectives and route of insulin infusion, i.e., intravenous [1] and subcutaneous [2]. The notable control techniques developed in this regard include the conventional proportional-integral-derivative (PID) control [5–7], H1 based robust control [8–10], model predictive control (MPC) [11–14] and linear quadratic regulator (LQR) [15, 16]. Although a majority of these algorithms are for clinical study and validated on mathematical models, many techniques have also been validated in real clinical settings [4, 7, 12, 17–19]. MPC and LQR have proved to be important tools among the state feedback dependent optimal control techniques. These techniques employ weighting matrices for minimizing glucose level and insulin infusion.
*Correspondence to: Subhojit Ghosh, Electrical Engineering, NIT Rourkela, India 769008. † E-mail: [email protected] Copyright © 2012 John Wiley & Sons, Ltd.
S. GHOSH AND S. GUDE
For an observer based state feedback control of a plant corrupted by state and measurement noise, the control action and the appropriateness of the estimated states is heavily dependent on the output and control weighting matrices along with the process and noise covariance matrices. The selection of these parameters is not a trivial problem, and is hence carried out by trial and error. This involves maintaining a trade-off between minimizing the control effort and improving the transient response. An increased weightage towards minimizing glucose level leads to high insulin infusion rate and vice-versa for attaching greater weightage for control effort. Similarly, the covariance matrices also affects the overall control objective by affecting the error dynamics of the observed state vector. In this regard, a genetic algorithm (GA) tuned Kalman filter based LQR theory has been proposed for blood glucose regulation in diabetic subjects in case of a disturbance in the form of a meal. GA [20], which has proved to be a powerful evolutionary technique in recent times, has been utilized for the selection of the weighting matrices of the LQR controller and the noise covariance matrices of the Kalman filter. Applications of genetic algorithms in glucose metabolism have been reported for parameter estimation [21, 22] and quantification of insulin sensitivity [23] from glucose insulin data. For a given patient model, the algorithm searches for the best combination of the tuning parameters that results in the least deviation of glucose from an acceptable basal value. Simulation studies conducted on a physiologically relevant model depicts the effectiveness of the proposed controller in maintaining normal glucose regulation in type 1 diabetes mellitus subjects, using both instantaneous as well delayed glucose measurements. The paper is organized as follows: the next section focuses on a mathematical model of glucose regulation, the Kalman filter based linear quadratic regulator for glucose control is discussed in Section 3, Section 4 gives a brief overview of GA along with its application for tuning the controller and Kalman filter parameters, simulation results for analyzing the controller performance are given in Section 5, and finally Section 6 provides conclusions. 2. MATHEMATICAL MODEL FOR GLUCOSE REGULATION For controller design, the mathematical model proposed by Sturis et al. [24] is considered in the present work. The nonlinear model takes following form: dIp dt dIi dt dG dt dI1 dt dI2 dt dI3 dt
Ip Ii / Vi Vp Ii Ip / D E. Vp Vi D E.
Ip C Uin .t / tp Ii ti
D Din .t / f2 .G/ f3 .G/f4 .Ii / C f5 .I3 / 3 .Ip I1 / td 3 D .I1 I2 / td 3 D .I2 I3 / td
(1)
D
For a normal subject, the following feedback loops are included in the model: pancreatic insulin secretion by stimulation of glucose, insulin stimulates glucose uptake and inhibits hepatic glucose production and glucose enhances its own uptake. However, because the present work involves designing a controller for type 1 diabetic subject, the pancreatic insulin secretion is not considered. The model has three main physiological variables; the concentration of glucose in the plasma G(mg/dL), the concentration of insulin in the plasma Ip (mU/L)and the insulin concentration in the intracellular space Ii (mU/L). In addition to the physiological variables, the model also contains three auxiliary insulin variables I1 , I2 , and I3 . The distribution volume (L) in the plasma insulin, intracellular insulin, and plasma glucose space is represented by Vp , Vi , and Vg , respectively, E(L/min) is the diffusion transfer rate, td (min) is the delay in hepatic glucose production, Ui n .t / Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
GA TUNED OPTIMAL CONTROLLER FOR GLUCOSE REGULATION IN TYPE 1 DIABETES
is the insulin infusion to the plasma compartment, tp (min) and ti (min) are the insulin degradation time constants in the plasma and intercellular space. Functions f2 .G/ and f3 .G/f4 .Ii / represent the insulin independent and insulin dependent glucose utilization, respectively, whereas, function f5 .I3 / represents hepatic glucose production. The controller designed for the above patient model aims to compensate the increase in glucose level caused by the disturbance Di n .t / (mg/dL min1 / by the application of insulin infusion (control input) . In comparison with the intravenous insulin infusion, subcutaneous delivery is more realistic in clinical setting. Hence, it is necessary to augment the above model with the dynamics representing the passage of insulin from the subcutaneous tissue to the plasma insulin compartment. For this purpose, the following model [25] of subcutaneous insulin kinetics has been considered dS1 D k21 S1 C U.t / dt dS2 D k21 S1 .kd C ka /S2 dt
(1a)
where U.t / is the rate of insulin infusion, S1 is the subcutaneous mass where injection takes place, S2 is the subcutaneous insulin mass proximal to the plasma, kd and ka are the degradation constants in the subcutaneous tissue and plasma, respectively, and k21 is the rate constant describing insulin transport. The term ka S2 =Vp corresponds to the rate of insulin infusion to the plasma space from subcutaneous injection. To arrive at a controller for glucose regulation, the integrated model is linearized about the steady state point (basal value) and a new set of variables, i.e., their deviation from the basal values are considered. 3. LINEAR QUADRATIC REGULATOR OPTIMAL FEEDBACK In this section, we briefly explain the observer (Kalman filter) based LQR method [26, 27]. In the present work, the LQR controller is formulated by considering subcutaneous insulin infusion as the control input and glucose appearance rate in the plasma from meal intake as disturbance. The linearized regulator model for controlling glucose levels using state feedback subjected to meal disturbance, can be written in the state space form as dx D Ax C Bu C din dt y D Cx
(2)
where A, B, and C are system matrices and the state vector x is given by x D ŒIp , Ii , G, I1 , I2 , I3 , S1 , S2 , with representing the deviation of the states of (1), from their basal values. The control input is given by u and di n represents the disturbance (di n = [0 0 Di n 0 0 0 0 0]). The LQR method involves the determination of the optimal control input u D kx, that minimizes the following performance index: Z1 JLQR D
.x T Qx C uT Ru/dt
(3)
0
where Q and R are positive definite matrices. These matrices denote the importance attached to the state and control output in the minimization of J . The controller gain K, which minimizes (3), is given by K D R1 B T P ,
(4)
where P is the positive symmetric solution of the following algebraic Riccatti equation PA C AT P C Q PBR1 B T P D 0. Copyright © 2012 John Wiley & Sons, Ltd.
(5)
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
S. GHOSH AND S. GUDE
The selection of Q and R significantly affects the solution of (5) and hence the controller performance. In the conventional approach, Q and R are selected from a number of trial and error iterations. The application of the control input u involves the availability of all the states at the output. However, the online measurement of plasma and intracellular insulin requires a highly invasive and expensive protocol. Also, the auxiliary insulin variables incorporated in the model (1) for representing a delay are nonphysiological in nature. Continuous monitoring of plasma glucose levels, through the subcutaneous route has been made possible by the developments in the field of implantable glucose sensor technology. Hence, of all the states in (2), for real-time control only G is considered to be available in the output. Thus, for applying state feedback control for glucose regulation, an observer is required for the estimation of nonaccessible and nonphysiological state variables. Generalizing the model (2), to include a wider class of systems by including modeling error and measurement noise leads to following representation: dx dt
D Ax C Bu C din C F v y D Cx C ´
(6)
where v is the process noise, which may occur because of modeling error such as noninclusion of nonlinearities and high frequency dynamics, and ´ is the measurement noise. The noises are independent of each other, white with normal probability distribution. In addition, all the state variables are assumed to be affected by the process noise and a single value is considered for the state equation. Kalman filter estimates the state variables.x/of Q the above stochastic plant from the input and noisy output measurement by minimizing the covariance of the conditional error between the actual and the estimated state, given as R D EŒ.x x/ Q T .x x/. Q
(7)
For a predefined noise covariance (Z for measurement noise and V for process noise), the estimated state equation of the Kalman filter is given as d xQ D AxQ C Bu C L.y C x/, Q dt
(8)
where the Kalman gain vector L is given as L D RC T Z 1 .
(9)
The solution of the matrix R is given by the following algebraic Riccatti equation AR C RAT RC T Z 1 CR C F VF T D 0
(10)
Similar to the weighting matrices Q and R for the LQR controller, the covariance matrices Z and V can be considered as tuning parameters for the Kalman filter. Their selection affects the dynamics of the error between actual and estimated state. Since in most cases the power spectral density of the noises are not known in advance, Zand V are also chosen by trial and error. Hence, the two-stage design of an optimal controller for glucose regulation from estimated states involves the design of an optimal regulator that minimizes (3) and a Kalman filter that minimizes (7). The linearity of the system and considering the noises as inputs, allows separating the controller design problem into two isolated processes of estimation and control. The regulator and Kalman filter can be designed separately and combined to form an optimal compensator. The conventional approach of designing a Kalman filter based LQR controller, involves a large number of trial and error iterations for choosing the best combination of the tuning parameters Q, R, Z, and F . The large computational time retards the application of this approach for large-scale systems. In the next section, a technique based on GA is proposed to tune the parameters of the LQR controller and Kalman filter for improving the integrated control action. Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
GA TUNED OPTIMAL CONTROLLER FOR GLUCOSE REGULATION IN TYPE 1 DIABETES
4. GENETIC ALGORITHM BASED CONTROLLER TUNING 4.1. Genetic algorithm Genetic algorithm is an effective search method based on natural evolution with the Darwinian survival of the fittest approach. Unlike conventional optimization algorithms, GA works on coded parameters and not on actual parameters. It does not involve derivation of the objective function and can hence work on discontinuous functions. This advantage allows GA to solve complex multivariable problems, which is difficult or impossible to be solved by traditional methods. GA involves three major operations of reproduction, crossover, and mutation. For an initially randomly selected possible solutions (population) coded as binary bits, during reproduction individual solutions known as chromosomes are selected based on their fitness, and sent to the mating pool. During crossover operation, two strings selected at random from the mating pool undergo crossover with a certain probability at a randomly selected crossover point to generate two new strings. Depending on whether a randomly generated number is larger than a predefined mutation probability or not each bit in the string obtained after crossover is altered (changing 0 to 1 and 1 to 0). The three operations are carried out for a certain number of iterations (generations), or till a chromosome with a certain fitness is obtained. 4.2. Genetic algorithm based linear quadratic regulator controller design S for the present problem, the aim is to control glucose level only, the diagonal elements of matrix Q in (3) for all the variables other than plasma glucose has been constrained to zero. Also, the same covariance has been considered for all the process noises. This results in a total of four tuning
Figure 1. Flowchart for the GA tuned LQR controller. Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
S. GHOSH AND S. GUDE
parameters (one each for Q, R, Z, and V /. GA selects the tuning parameters, so that the controller with the corresponding gains (K and L/ is able to maintain normal glucose levels, in case of a meal disturbance. The detailed procedure for obtaining the parameters is shown in Figure 1. The tuning parameters are encoded into chromosomes of 40 bits (10 bits each for each parameter). The chromosome bits corresponding to Q and R are selected based on the minimization of the following fitness function: Ztf " JRegulat or D 0
G.t / Gbasal
2
U.t / C Ubasal
2 # dt ,
(11)
where, Gbasal and Ubasal represent the basal glucose level and insulin infusion rate, respectively. For a given Q and R, the fitness function is calculated by simulating the closed loop patient model for a given meal disturbance, until the predisturbance glucose level is attained. As mentioned earlier, the selection of Q and R has a direct bearing on the response and magnitude of the control input of the closed loop system. Selecting an appropriate state-weighting matrix Q imposes certain limits on the maximum overshoot of plasma glucose. Similarly, the maximum overshoot limit in the insulin
Figure 2. GA tuned LQR controller response for 50 g meal disturbance. (a) Plasma glucose concentration and (b) insulin infusion rate. Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
GA TUNED OPTIMAL CONTROLLER FOR GLUCOSE REGULATION IN TYPE 1 DIABETES
infusion rate is specified through R. For the Kalman filter, the chromosome bits are selected based on the following objective function: JKalman D
n n X X
Re2 .i, j /.
(12)
i D1 j D1
For each chromosome, the feedback gain K, is determined from the solution of (5), whereas the Kalman gain L, is obtained from the solution of (10). The gain factors are used to calculate JRegulat or and JKalman from the controller response. 5. RESULTS AND DISCUSSION The performance of the GA tuned LQR controller in compensating a meal disturbance from the feedback of noisy glucose measurement is considered in this section. Intravenous route is considered for glucose measurement. The controller performance is evaluated based on the excursion of the plasma glucose level from the basal value and the amount of insulin infused to compensate the disturbance. For simulating the meal dynamics of certain carbohydrate content, the bicompartmental
Figure 3. GA tuned LQR controller response for meal disturbances of 20 g (breakfast), 50 g (lunch) and 40 g (dinner) in a day. (a) Plasma glucose concentration and (b) insulin infusion rate. Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
S. GHOSH AND S. GUDE
model proposed by Lehmann and Deutsch in [28] has been considered. In the model, the conversion from meal intake (in terms of carbohydrate content) to glucose absorption into the plasma compartment is expressed as a trapezoid wave followed by a first-order linear system (Appendix). For tuning the GA, a large number of selection and crossover strategies were tested through pilot runs to determine the most suitable ones. Single point crossover is used in conjunction with elitist strategy based Roulette wheel selection. The elitist strategy, which is able to preserve superior strings, is incorporated in the present work by replacing the 10 worst strings of a particular generation with the 10 best strings of the previous generation. The number of chromosomes in the initial population and the maximum number of generations is set at 50 and 200 in each run. The crossover and mutation probability was fixed at 0.8 and 0.2, respectively. The algorithm was found to converge to the same solution for different intialization of the tuning parameters. The performance of the controller (in terms of plasma glucose level and insulin infusion rate) with the GA derived tuning parameters (Q(3,3)=7.294, R=0.015, V =6.158, Z=0.023), following a meal disturbance of 50 gm carbohydrate content is shown in Figure 2. Because there is no endogenous secretion of insulin, a continuos insulin infusion should always be present even in the absence of a disturbance, to maintain the basal glucose level. The closed loop performance of the diabetic patient for the entire day subjected to breakfast at 8 AM of 20 g, lunch at 1 PM of 50 g and dinner at 8 PM
Figure 4. GA tuned LQR controller response for 50 g meal disturbance using fitness function with glucose deviation only. (a) Plasma glucose concentration and (b) insulin infusion rate. Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
GA TUNED OPTIMAL CONTROLLER FOR GLUCOSE REGULATION IN TYPE 1 DIABETES
of 40 g is shown in Figure 3. Similar plasma glucose and insulin infusion rate profile during the initial stage for different meal intakes is because of the similar rate of glucose appearance into the blood circulation during that period, as obtained from the meal model [28]. However, prolonged and high glucose infusion after that results in larger settling time for meals of high carbohydrate content. Similar to single meal disturbance, for multiple meals also, the insulin infusion rate contains high frequency noise components of the order of ˙0.2 mU/L min approximately. For all the disturbances considered above, the controller is able to maintain the glucose levels within the normal range of 80–120 mg/dL. It should be noted that although the controller is designed based on the linearized model, the controller action results reported here is that on the nonlinear model (1). For many cases, depending on the condition of the subject and the available protocol for insulin infusion, the primary goal might be to keep the glucose level within the normal range without bothering on the rate of insulin infusion. In this regard, it may be advantageous to have periods where insulin infusion deviates a lot from the basal value. Such cases can be handled within the framework of the present work, by not considering the minimization of insulin infusion in the fitness function (11), that is, only the first component of (11) is taken into account. On the basis of this formulation, the response of the controller following a meal disturbance of 50 gm carbohydrate content is shown in Figure 4.
Figure 5. GA tuned LQR controller response for 50 g meal disturbance with a glucose sensor delay of 2 min. (a) Plasma glucose concentration and (b) insulin infusion rate. Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
S. GHOSH AND S. GUDE
An important aspect in controller design for blood glucose regulation is the time delay in glucose sensor measurement. In this paper, until now it has been assumed that the measurement of plasma glucose level, although noisy, is available instantaneously. Along with the degradation of the controller performance, large sensor delay results in an unstable closed loop system. To evaluate the degradation in controller performance, the proposed controller was subjected to measurement delays of 2 and 5 min. The variation in plasma glucose and insulin infusion rate following a meal of 50 g for the above two measurement delays are shown in Figures 5 and 6. For both cases, the glucose levels are well within the acceptable range. However, for large time delays (> 10 min), sustained oscillations are observed in glucose levels with high peak value. For continuous glucose infusion, the standard LQR that provides only proportional control is found to produce a steady state error. In addition to minimizing plasma glucose level and insulin infusion rate, an additional objective of minimizing the integral of the plasma glucose state is required for continuous disturbance. For this, an integral action is added to the standard LQR. The model (2) is augmented to include a new state variable, that is, the integral of the deviation in glucose level from basal value, given as
Figure 6. GA tuned LQR controller response for 50 g meal disturbance with a glucose sensor delay of 5 min. (a) Plasma glucose concentration and (b) insulin infusion rate. Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
GA TUNED OPTIMAL CONTROLLER FOR GLUCOSE REGULATION IN TYPE 1 DIABETES
Zt G. /d
(13)
0
With the additional tuning parameter for the augmented state, the response of the proposed controller with the GA derived tuning parameters (Q(3,3)=0.151, Q(7,7)=0.002, R=23.083, V =0.483, Z=3.92), for continuous glucose infusion of 1 mg/dL min is shown in Figure 7. Unlike the standard LQR controller, where the variation in weighting factors do not effect the time required to reach the steady state, varying the tuning parameter corresponding to the integral of glucose level has a direct effect on the settling time, with higher value resulting in faster convergence to the basal glucose level with oscillations in the steady state. Varying the integral weighting factor is similar to the tuning of integral term in a conventional PID controller. Unlike the regulation for meal disturbance, for constant glucose infusion an insulin infusion rate in excess of the basal value of 3.2 mU/L min is required to maintain normoglycemia at steady state. Oscillations in plasma glucose levels resulting in oscillatory insulin infusion rates should be avoided because the insulin delivery system may not
Figure 7. GA tuned LQR controller response for a continuous glucose infusion of 1 mg/dL min. (a) Plasma glucose concentration and (b) insulin infusion rate. Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
S. GHOSH AND S. GUDE
be able to alter the rate of infusion at a very fast rate. The continuous glucose infusion, although not representing a realistic disturbance to the glucose regulation system, is similar to the situation arising because of changes in the physiological unit processes involved. Changes in insulin sensitivity and retarded insulin independent glucose utilization may result in an offset from the steady state value of the original model. The resulting offset requires integral action to maintain the premodified basal glucose and insulin level. 6. CONCLUSIONS Evolutionary computation has been utilized to design an optimal feedback controller for glucose regulation in type 1 diabetic subjects using noisy output data. The inaccessible and nonphysiological state variables are estimated from the noisy glucose output measurement using a Kalman filter. The selection of the design parameters for the Kalman filter and the feedback controller based on LQR is in general a trial and error process, which involves selection of the best combination from different simulated responses. In the present work, these parameters are optimized using GA. The two-stage algorithm searches for the parameters that minimize the deviation in plasma glucose levels from the basal value and the covariance of the state estimation error. The performance of the controller is evaluated for different meal disturbances with both instantaneous and delayed sensor measurement. The proposed controller is able to maintain glucose levels in type 1 diabetic subjects within the normal range in case of meal disturbances of varying carbohydrate content. The simplicity of the approach and the absence of any online calculation make it easy to store in a digital chip, which can then be embedded with a suitable drug delivery system and an in vivo sensor for automated insulin delivery. Future work in this direction is aimed at testing the efficacy of the proposed strategy using more comprehensive models and studying the robust handling capacity in case of parameter variations. APPENDIX: MEAL DYNAMICS MODEL In the model [28], the absorption of glucose from a meal consists of two stages: one is the emptying of the stomach into the duodenum and the other is the absorption of glucose from the small intestine into the blood circulation. The rate of gastric emptying of carbohydrate is considered as a nonlinear function with saturation characteristics to maintain a constant carbohydrate release from the stomach during intestinal absorption. The rise and fall from the maximum saturated level emptying rate Vm,CHO = 360 mg/min is considered as a linear function spanning over 30 min. These time periods are denoted by Tri se and Tf al l . The period of constant carbohydrate release depends on the carbohydrate content of the meal consumed. Meals with carbohydrate content less than 10.8 g do not reach the saturation level. Larger meals are represented by trapezoid wave. The glucose absorption is modeled by a first-order equation, which takes the form GP gut D kab Ggut C Gempt , Din D kab Ggut
(A1)
where, Ggut is the glucose concentration in the intestine, kab is the rate of glucose absorption into the blood circulation. REFERENCES 1. Parker RS, Doyle III FJ, Peppas NA. The Intravenous Route to Blood Glucose Control. IEEE Engineering in Medicine and Biology 2001; 20(1):65–73. 2. Bellazzi R, Nucci G, Cobelli C. The Subcutaneous Route to Insulin-Dependent Diabetes Therapy. IEEE Engineering in Medicine and Biology 2001; 20(1):54–64. 3. Wayne Bequette B. A critical assessment of algorithms and challenges in the development of a closed-loop artificial pancreas. Diabetes Technology and Therapeutics 2005; 7(1):28–47. 4. Hovorka R. Continuous glucose monitoring and closed-loop systems. Diabetic Medicine 2005; 23(1):1–12. 5. Marchetti G, Barolo M, Jovanovic L, Zisser H, Seborg DE. An improved PID switching control strategy for type 1 diabetes. IEEE Transactions on Biomedical Engineering 2008; 55(3):857–865. Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
GA TUNED OPTIMAL CONTROLLER FOR GLUCOSE REGULATION IN TYPE 1 DIABETES
6. Lam ZH, Hwang KS, Lee JY, Chase JG, Wake GC. Active insulin infusion using optimal and derivative-weighted control. Medical Engineering and Physics 2002; 24(10):663–672. 7. Steil GM, Rebrin K, Darwin C, Hariri F, Saad MF. Feasibility of automating insulin delivery for the treatment of type 1 diabetes. Diabetes 2006; 55:3344–3350. 8. Parker RS, Doyle III FJ, Ward JH, Peppas NA. Robust H1 control in diabetes using a physiological model. AIChE Journal 2000; 46(12):2537–2549. 9. Quiroz G, Femat R. Theoretical blood glucose control in hyper- and hypoglycemic and exercise scenarios by means of an H1 algorithm. Journal of Theoretical Biology 2010; 263:154–160. 10. Ruiz Velazquez E, Femat R, Campos Delgado DU. Blood glucose control for type I diabetes mellitus: a robust H1 tracking problem. Control Engineering Practice 2004; 12:1179–1195. 11. Parker RS, Doyle III FJ, Peppas NA. A model-based algorithm for blood glucose control in type I diabetic patients. IEEE Transactions on Biomedical Engineering 1999; 46(2):148–157. 12. Hovorka R, Canonico V, Chassin LJ, Haueter U, Massi-Benedetti M, Orsini Federici M, Pieber TR, Schaller HC, Schaupp L, Vering T, Wilinska ME. Nonlinear model predictive control of glucose concentration in subjects with type 1 diabetes. Physiological Measurement 2004; 25(4):905–920. 13. Kovatchev BP, Breton MD, Dalla Man C, Cobelli C. In Silico Preclinical Trials: A Proof of Concept in Closed-Loop Control of Type 1 Diabetes. Journal of Diabetes Science and Technology 2009; 3(1):44–55. 14. Abu-Rmileh A, Garcia-Gabin W. Feedforward–feedback multiple predictive controllers for glucose regulation in type 1 diabetes. Computer Methods and Programs in Biomedicine 2010; 99:113–123. 15. Patek SD, Breton MD, Chen Y, Solomon C, Kovatchev B. Linear quadratic Gaussian based closed loop control of type 1 diabetes. Journal of Diabetes Science and Technology 2007; 1(6):834–841. 16. Sanchez Chavez IY, Morales Menendez R, Martinez Chapa SO. Glucose optimal control system in diabetes treatment. Applied Mathematics and Computation 2009; 209:19–30. 17. Chase JG, Shaw GM, Wong XW, Lotz T, Lin J, Hann CE. Model-Based Glycaemic Control in Critical Care-A Review of the State of the Possible. Biomedical Signal Processing and Control 2006; 1:3–21. 18. Palerm CS, Jovanovic L, Zisser H, Doyle III FJ, Bevier WC. Prandial Insulin Dosing Using Run-to-Run ControlApplication of clinical data and medical expertise to define a suitable performance metric. Diabetes Care 2007; 30(5):1131–1136. 19. El-Khatib FH, Jiang J, Damiano ER. A feasibility study of bihormonal closed-loop blood glucose control using dual subcutaneous infusion of insulin and glucagon in ambulatory diabetic swine. Journal of Diabetes Science and Technology 2009; 3(4):789–803. 20. Goldberg DE. Genetic Algorithms in Search, Optimization, and Machine Learning. Addison-Wesley: Reading, MA, 1989. 21. Morbiducci U, Tura A, Grigioni M. Genetic algorithms for parameter estimation in mathematical modeling of glucose metabolism. Computers in Biology and Medicine 2005; 35:862–874. 22. Morbiducci U, Benedetto GD, Kuatzky-Willer A, Pacini G, Tura A. Improved usability of the minimal model of insulin sensitivity based on an automated approach and genetic algorithms for parameter estimation. Clinical Science 2007; 112:257–263. 23. Ghosh S, Maka S. Genetic Algorithm Based NARX model Identification for Evaluation of Insulin Sensitivity. Applied Soft Computing 2011; 11(1):221–226. 24. Sturis J, Polonsky KS, Mosekilde E, Van Cauter E. Computer-model for mechanisms underlying ultradian oscillations of insulin and glucose. American Journal of Physiology 1991; 260:E801–E809. 25. Shimoda S, Nishida K, Sakakida M, Nishida K, Sakakida M, Konno Y, Ichinose K, Uehara M, Nowak T, Shichiri M. Closed-loop subcutaneous insulin infusion algorithm with a short acting insulin analog for long term clinical application of a wearable artificial endocrine pancreas. Frontiers in Medical Biology and Engineering 1997; 8:197–211. 26. Naidu DS. Optimal Control Systems. CRC Press: Boca Raton, Florida, 2003. 27. Crassidis JL, Junkins JL. Optimal Estimation of Dynamic Systems,. Chapman & Hall/CRC: Boca Raton, Florida, 2004. 28. Lehmann ED, Deutsch T. A physiological model of glucose-insulin interaction in type 1 diabetes mellitus. Journal of Biomedical Engineering 1993; 14(3):235–242.
Copyright © 2012 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Biomed. Engng. (2012) DOI: 10.1002/cnm
