A mathematical model for the robust blood glucose tracking

Weijiu Liu; Weijiu Liu

doi:10.3934/mbe.2019036

Mathematical Biosciences and Engineering

2019, Volume 16, Issue 2: 759-781. doi: 10.3934/mbe.2019036

Previous Article Next Article

Research article

A mathematical model for the robust blood glucose tracking

Weijiu Liu ^,

Department of Mathematics, University of Central Arkansas, 201 Donaghey Avenue, Conway, AR 72035, USA

Received: 20 August 2018 Accepted: 23 November 2018 Published: 15 January 2019

In this paper, we study the problem of the robust blood glucose tracking. Tracking here means that the error between a state variable of a system under control and its desired time-varying reference converges to zero over time. Robustness here means that a controller designed for a system can tolerate a small variation of the system parameters. Since the parameters in the blood glucose regulation system differ in people, such a robust controller is useful in the insulin pump technology: an insulin pump equipped with such a robust controller could be used in a group of people. Thus, in our study, parameter uncertainties are introduced into a mathematical model of the blood glucose regulation system. Using an actual blood glucose level as feedback and an exogenous glucose input and a desired glucose reference as feedforward, we design a robust feedback and feedforward controller, which drives the blood glucose to track the desired time-varying glucose reference for any small uncertainties. Numerical simulations with published experimental blood glucose data are conducted to further confirm our theoretical results.

Keywords:

Citation: Weijiu Liu. A mathematical model for the robust blood glucose tracking[J]. Mathematical Biosciences and Engineering, 2019, 16(2): 759-781. doi: 10.3934/mbe.2019036

Related Papers:

[1]	Lela Dorel . Glucose level regulation via integral high-order sliding modes. Mathematical Biosciences and Engineering, 2011, 8(2): 549-560. doi: 10.3934/mbe.2011.8.549
[2]	Danilo T. Pérez-Rivera, Verónica L. Torres-Torres, Abraham E. Torres-Colón, Mayteé Cruz-Aponte . Development of a computational model of glucose toxicity in the progression of diabetes mellitus. Mathematical Biosciences and Engineering, 2016, 13(5): 1043-1058. doi: 10.3934/mbe.2016029
[3]	Lingmin Lin, Kailai Liu, Huan Feng, Jing Li, Hengle Chen, Tao Zhang, Boyun Xue, Jiarui Si . Glucose trajectory prediction by deep learning for personal home care of type 2 diabetes mellitus: modelling and applying. Mathematical Biosciences and Engineering, 2022, 19(10): 10096-10107. doi: 10.3934/mbe.2022472
[4]	Xingjia Li, Jinan Gu, Zedong Huang, Chen Ji, Shixi Tang . Hierarchical multiloop MPC scheme for robot manipulators with nonlinear disturbance observer. Mathematical Biosciences and Engineering, 2022, 19(12): 12601-12616. doi: 10.3934/mbe.2022588
[5]	Xiaohan Yang, Yinghao Cui, Zhanhang Yuan, Jie Hang . RISE-based adaptive control of electro-hydraulic servo system with uncertain compensation. Mathematical Biosciences and Engineering, 2023, 20(5): 9288-9304. doi: 10.3934/mbe.2023407
[6]	Kimberly Fessel, Jeffrey B. Gaither, Julie K. Bower, Trudy Gaillard, Kwame Osei, Grzegorz A. Rempała . Mathematical analysis of a model for glucose regulation. Mathematical Biosciences and Engineering, 2016, 13(1): 83-99. doi: 10.3934/mbe.2016.13.83
[7]	Balázs Csutak, Gábor Szederkényi . Robust control and data reconstruction for nonlinear epidemiological models using feedback linearization and state estimation. Mathematical Biosciences and Engineering, 2025, 22(1): 109-137. doi: 10.3934/mbe.2025006
[8]	Hugo Flores-Arguedas, Marcos A. Capistrán . Bayesian analysis of Glucose dynamics during the Oral Glucose Tolerance Test (OGTT). Mathematical Biosciences and Engineering, 2021, 18(4): 4628-4647. doi: 10.3934/mbe.2021235
[9]	Huazhong Yang, Wang Li, Maojin Tian, Yangfeng Ren . A personalized multitasking framework for real-time prediction of blood glucose levels in type 1 diabetes patients. Mathematical Biosciences and Engineering, 2024, 21(2): 2515-2541. doi: 10.3934/mbe.2024111
[10]	Xingjia Li, Jinan Gu, Zedong Huang, Wenbo Wang, Jing Li . Optimal design of model predictive controller based on transient search optimization applied to robotic manipulators. Mathematical Biosciences and Engineering, 2022, 19(9): 9371-9387. doi: 10.3934/mbe.2022436

Abstract

1. Introduction

Blood glucose needs to be maintained within a narrow range, for example, from 3.6 to 6.0 mmol/l (65 to 108 mg/dl) for the fasting people without diabetes. This narrow range differs in people and in time. A high blood glucose concentration outside the range is a symptom of diabetes, one of the prevalent diseases affecting millions of people worldwide. A peptide hormone insulin is required to interact with glucose to achieve this goal.

The interaction between blood glucose and insulin involves multiple feedback control mechanisms and is of extremely high complexity. Here we give a very basic description about it. Glucose comes from food and liver, and utilized by brain and nerve cells (insulin-independent) via the glucose transporter 3 (GLUT3) or by tissue cells such as muscle, kidney, and fat cells (insulin-dependent) via the glucose transporter 4 (GLUT4). Glucose is transported into and out of liver cells by the concentration-driven GLUT2, which is insulin-independent^[17]. In response to a low blood glucose level, the $\alpha$ cells of the pancreas produce the hormone glucagon. The glucagon initiates a series of activations of kinases, and finally leads to the activation of the glycogen phosphorylase, which catalyzes the breakdown of glycogen into glucose^[14]. In addition, the series of activations of kinases also result in the inhibition of glycogen synthase and then stop the conversion of glucose to glycogen. In response to a high blood glucose level, the $\beta$ cells of the pancreas secrete insulin^[15,35,36]. Insulin triggers a series of reactions to activate glycogen synthase, which catalyzes the conversion of glucose into glycogen^[16,21]. Insulin also initiates a series of activations of kinases in tissue cells to lead to the redistribution of the glucose transporter 4 (GLUT4) from intracellular storage sites to the plasma membrane^[13]. Once at the cell surface^[8,32], GLUT4 transports glucose into the muscle or fat cells.

Mathematical models for the blood glucose regulation system are required for integrating a glucose monitoring system into insulin pump technology to form a closed-loop insulin delivery system, the so-called "artificial pancreas" (see Hovorka^[18], Panteleon et al.^[33], Steil et al.^[38]). To make this artificial pancreas as close as possible to the natural pancreas, many mathematical models about the interaction between blood glucose and insulin have been proposed. In many of these models, the compartmental modeling approach has been used. In the 1960's, Ackerman et al.^[1,2] proposed a simplified compartmental model of two linear ordinary differential equations for the glucose-tolerance test. This model lumped the large number of kinetic parameters into a much smaller number which can, at least in part, characterize the human glucose regulatory system. In the 1970's, Bergman et al.^[4,5,6,7] proposed a minimal compartmental model of nonlinear ordinary differential equations. It was assumed that the blood glucose and hormone insulin are contained in two different compartments and interact with each other. The plasma insulin is secreted from the endocrine system under the stimulation of the blood glucose and enters a "remote compartment", where it is active in accelerating glucose utilization. Then the dynamics of the blood glucose, plasma insulin, and the insulin in the remote compartment was modeled by a system of three nonlinear ordinary differential equations. In the 1990's, to determine whether the oscillations of insulin and blood glucose could result from the feedback loops between insulin and glucose, Sturis et al.^[39,40] developed a parsimonious model including the major mechanisms involved in glucose regulation. This model comprises two major negative feedback loops describing the effects of insulin on glucose utilization and glucose production, respectively, and both loops include the stimulatory effect of glucose on insulin secretion. Then this model consists of six nonlinear ordinary differential equations, three of which are used to model a time delay for the effect of insulin on glucose production and a sluggish effect of insulin on glucose utilization. On the basis of Sturis et al.'s model, Li et al.^[23] introduced two explicit time delays and proposed a more robust alternative model for better understanding the glucose-insulin endocrine metabolic regulatory system and the ultradian insulin secretory oscillations for the cases of continuous enteral nutrition and constant glucose infusion. More complicated compartmental models proposed by Hovorka et al.^[18], Liu et al.^[24,25,27], Man et al.^[29,30,31], and Sorensen^[37] have considered labeling more compartments for better understanding the behavior of different parts of the body. To assess the liver ability to extract insulin, Toffolo et al.^[41,42], proposed a minimal model of insulin delivery and kinetics by combining the classical C-peptide minimal model. In the 2010's, Vahidi et al.^[45,46] developed a model for type Ⅱ diabetes, based on a detailed compartmental physiological model proposed by Sorensen^[37] for a healthy human body. There have been a huge set of references on the modeling research and it is hopeless to mention all of them.

In all these models, no parameter uncertainties were introduced. The parameter uncertainties here means that a parameter may vary in a neighborhood of its nominal value. This can happen in practice. For instance, an estimation error may be introduced when a model parameter is estimated. In addition, parameter values may vary slightly from person to person. Therefore, it is of practical interest to introduce parameter uncertainties into a model. As far as the robust tracking is concerned, the robustness against the parameter uncertainties means that a feedback and feedforward controller should be designed in a way such that it can tolerate the parameter variation in a small neighborhood of their nominal values. That is, even though the parameters have a small change, the controller can still achieve its tracking goal. Such a robust problem may be solved by only a robust control design method and other methods such as the sensitivity analysis, bootstrapping to quantify parameter ranges, and stochastic modeling might not be relevant.

As a first attempt in this direction, we introduce parameter uncertainties into the model proposed by Bergman et al.^[4,5]. Although this model is simple and consists of only three differential equations, one equation for the glucose dynamics, one equation for the effect of the remote insulin on glucose, and the other equation for the hormone insulin dynamics, it can simulate well the major biological dynamics about the interaction between glucose and insulin, as evidenced in the numerous literature (see. e.g., ^[4,5,6,7,11]). It has been widely used and further updated later (see, e.g., ^{[10,18,20,23,24,39,40,41]}).

Furthermore, we also include a time-dependent exogenous glucose input $J$ from food through the intestine in this model and assume that the desired glucose reference $r$ to be tracked is time-dependent, as described in . Then we design a robust feedback (using an actual glucose concentration $g$ ) and feedforward (using the exogenous glucose input $J$ and the desired glucose reference $r$ ) controller $u(g, r, J)$ such that the glucose concentration asymptotically tracks the time-varying glucose reference for any small uncertainties. Numerical simulations conducted with published experimental blood glucose data show that our controller is effective in achieving the blood glucose tracking and agree with our theoretical results.

Figure 1. Feedback and feedforward control flow diagram. For a detailed description, see the text.

DownLoad: Full-Size Img PowerPoint

2. Formulation of a model with parameter uncertainties

In the 1970's, Bergman et al.^[4,5] proposed a mathematical model for the blood glucose regulation system. To establish such a model, they made a simplified assumption. The blood glucose is produced from food through the intestine. The blood glucose and hormone insulin are contained in two different compartments and interact with each other. The blood glucose can be converted into the liver glycogen and used by the tissue metabolism with or without the help of insulin. The plasma insulin is secreted from the endocrine system under the stimulation of the blood glucose and enters a "remote compartment", where it is active in accelerating glucose utilization. It is used up by the insulin metabolism. Under this assumption, the dynamics of the blood glucose and insulin was modeled by the following system of nonlinear differential equations ^[5]:

$\begin{eqnarray} \frac{dg}{dt} & = & -m_1 g- h_a g+J , \end{eqnarray}$

(2.1)

$\begin{eqnarray} \frac{dh_a}{dt} & = & -m_2 h_a+m_3h , \end{eqnarray}$

(2.2)

$\begin{eqnarray} \frac{dh}{dt} & = & -m_4 h+u . \end{eqnarray}$

(2.3)

In the above equations, $g$ and $h$ denote concentrations of blood glucose and plasma insulin, respectively, $h_a$ is the effect of the remote insulin on glucose, $J$ is a rate of the exogenous glucose input from the intestine, $u$ is a rate of insulin secreted from the endocrine system or infused externally, and $m_1, m_2, m_3, m_4$ are positive rate constants. The biological meaning of these rate constants are explained in the Table 1. Because the glucose needs to be monitored and regulated, we introduce an output equation:

Table 1. Values of parameters of the model (3.32)-(3.36).

Parameter	Value	Description
$m_1$	0.000296 (/min)	Glucose utilization rate constant
$m_2$	0.045 (/min)	Insulin effect rate constant
$m_3$	0.000012 (/(pmol/l)/min $^2$ )	Insulin effect rate constant
$m_4$	0.268 (/min)	Insulin metabolism constant
$G_b$	0.3218 (mM/min)	Basal glucose input
$\bar{g}$	6 (mM)	Glucose equilibrium
$\bar{h}$	200 (pmol/l)	Insulin equilibrium
$\bar{h}_a$	0.0533 (/min)	Insulin effect equilibrium
$\bar{u}$	53.63 (pmol/l/min)	Insulin input equilibrium

| Show Table

DownLoad: CSV

$\begin{equation} y = g . \end{equation}$

(2.4)

To maintain the blood glucose at a normal level of about 6 mmol/l, a basal exogenous glucose input $G_b$ is required^[24,47]. Thus the blood glucose input $J$ could be split into two parts:

$\begin{equation} J = G_b+G_d, \end{equation}$

(2.5)

where $G_d$ can be treated as a glucose disturbance.

In what follows, $(\bar{g}, \bar{h}_a, \bar{h}, \bar{u})$ denotes the equilibrium of the system (2.1)-(2.3) without the glucose disturbance. That is, they are solutions of the following steady state system:

$\begin{eqnarray} 0 & = & -m_1 \bar{g}- \bar{g}\bar{h}_a +G_b , \end{eqnarray}$

(2.6)

$\begin{eqnarray} 0 & = & -m_2 \bar{h}_a+m_3\bar{h} , \end{eqnarray}$

(2.7)

$\begin{eqnarray} 0 & = & -m_4 \bar{h}+\bar{u} . \end{eqnarray}$

(2.8)

In practice, an estimation error may be introduced when a model parameter is estimated. In addition, the values of the parameters $m_1, m_2, m_3, m_4$ may vary slightly from person to person. Therefore, it is reasonable to introduce parameter uncertainties into the system (2.1)-(2.3):

$\begin{eqnarray} \frac{dg}{dt} & = & -(m_1+w_1) g- gh_a+ w_1 \bar{g} + G_b+G_d, \end{eqnarray}$

(2.9)

$\begin{eqnarray} \frac{dh_a}{dt} & = & - (m_2 +w_2) h_a+(m_3 +w_3) h+ w_2 \bar{h}_a-w_3 \bar{h} , \end{eqnarray}$

(2.10)

$\begin{eqnarray} \frac{dh}{dt} & = & - (m_4 +w_4) h+w_4 \bar{h}+u . \end{eqnarray}$

(2.11)

In the above equations, $m_1, m_2, m_3$ and $m_4$ represent the nominal part of the parameters and $w_1, w_2, w_3, w_4$ represent the uncertain part. Since the equilibrium $(\bar{g}, \bar{h}_a, \bar{h}, \bar{u})$ would not depend on the uncertainties, the terms $w_1 \bar{g}, w_2 \bar{h}_a, w_3 \bar{h}$ and $w_4\bar{h}$ are included such that this is guaranteed.

Because the blood glucose needs to be maintained within a narrow range from 3.6 to 6.0 mmol/l (65 to 108 mg/dl), we need to develop an insulin infusion rate $u$ in response to a glucose surge to achieve this goal. This problem can be mathematically formulated as follows. Given a glucose reference $r(t)$ (a constant 6 mmol/l or fluctuating around 6 mmol/l), we want to design a controller $u = u(g, J, r)$ such that for all sufficiently small uncertainties $w_1, w_2, w_3, w_4$

$\lim\limits_{t\to \infty} (g(t)- r(t)) = 0.$

3. A static feedback and feedforward controller

To solve this control problem, we need to express $G_d$ and $r$ by using solutions of linear differential systems. The glucose produced from food via the intestine is not exactly known and might be random. This glucose disturbance could be of periodic nature and then can be approximated by a Fourier polynomial

$\begin{eqnarray} G_d& = & a_0 +a_1\sin (\omega_d t) +b_1\cos(\omega_d t)+a_2\sin (2\omega_d t)+b_2\cos(2\omega_d t)+\cdots\\ &&+ a_n\sin (n\omega_d t)+b_n\cos(n\omega_d t). \end{eqnarray}$

(3.1)

The number $1$ and functions cosine's and sine's can be generated by the linear differential system

$\frac{d\mathbf{v}_d}{dt} = \mathbf{A}_{d} \mathbf{v}_d,$

where

$\mathbf{A}_{d} = \left[\begin{array}{cccccccc} 0&0&0&0&0& \cdots& 0&0\\ 0& 0& \omega_d& 0 & 0& \cdots& 0&0 \\ 0&-\omega_d &0&0&0& \cdots& 0&0\\ 0&0&0&0&2\omega_d& \cdots& 0&0\\ 0&0&0&-2\omega_d&0& \cdots& 0&0\\ \vdots&\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\ 0&0&0&0&0& \cdots& 0&n\omega_d\\ 0&0&0&0&0& \cdots& -n\omega_d&0 \end{array} \right].$

Then $G_d$ can be expressed as

$G_d = \mathbf{C}_d\mathbf{v}_d,$

where $\mathbf{C}_d = [a_0, a_1, b_1, a_2, b_2, \cdots, a_n, b_n]$ is a constant vector.

Since the glucose reference $r$ is of periodic nature, too, it can be generated in the same way as follows:

$r = \bar{g}+\mathbf{C}_r\mathbf{v}_r,$

where $\bar{g}$ is the glucose equilibrium, $\mathbf{C}_r$ is a constant vector and $\mathbf{v}_r$ is governed by

$\frac{d\mathbf{v}_r}{dt} = \mathbf{A}_{r} \mathbf{v}_r$

with $\mathbf{A}_{r}$ being similar to $\mathbf{A}_{d}$ :

$\mathbf{A}_{r} = \left[\begin{array}{cccccccc} 0&0&0&0&0& \cdots& 0&0\\ 0& 0& \omega_r& 0 & 0& \cdots& 0&0 \\ 0&-\omega_r &0&0&0& \cdots& 0&0\\ 0&0&0&0&2\omega_r& \cdots& 0&0\\ 0&0&0&-2\omega_r&0& \cdots& 0&0\\ \vdots&\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\ 0&0&0&0&0& \cdots& 0&m\omega_r\\ 0&0&0&0&0& \cdots& -m\omega_r&0 \end{array} \right].$

Let

$\mathbf{v} = \left[\begin{array}{c} \mathbf{v}_r\\ \mathbf{v}_d \end{array} \right],\quad \mathbf{A}_e = \left[\begin{array}{cc} \mathbf{A}_{r}&0\\ 0& \mathbf{A}_{d} \end{array} \right].$

Then glucose reference and disturbance can be lumped together as follows:

$\begin{equation} \frac{d\mathbf{v} }{dt} = \mathbf{A}_e \mathbf{v}. \end{equation}$

(3.2)

We introduce the tracking error:

$e = g-r = g-\bar{g}-\mathbf{C}_r\mathbf{v}_r.$

Combining the exosystem (3.2) with the model (2.9)-(2.11), we obtain a composite system

$\begin{eqnarray} \frac{dg}{dt} & = & -(m_1+w_1) g- gh_a+ w_1 \bar{g} + G_b+[\mathbf{0},\mathbf{C}_d]\mathbf{v}, \end{eqnarray}$

(3.3)

$\begin{eqnarray} \frac{dh_a}{dt} & = & - (m_2 +w_2) h_a+(m_3 +w_3) h+ w_2 \bar{h}_a-w_3 \bar{h} , \end{eqnarray}$

(3.4)

$\begin{eqnarray} \frac{dh}{dt} & = & - (m_4 +w_4) h+w_4 \bar{h}+u , \end{eqnarray}$

(3.5)

$\begin{eqnarray} \frac{d\mathbf{v} }{dt}& = &\mathbf{A}_{e} \mathbf{v}, \end{eqnarray}$

(3.6)

$\begin{eqnarray} e & = &g- \bar{g}-[\mathbf{C}_r, \mathbf{0}]\mathbf{v}. \end{eqnarray}$

(3.7)

When the disturbance $G_d = 0$ and $r = \bar{g}$ , we can readily show that $g, h_a, h$ and $u$ converge to their constant equilibrium, respectively. When $G_d$ and $r$ are time-dependent, we can expect that $g, h_a, h$ and $u$ would converge to their time-varying "equilibrium" that would depend on $\mathbf{v}$ . In fact, such time-varying "equilibrium" called a stable center manifold is guaranteed by the center manifold theory (see, e.g., ^[9]). Thus we introduce the following $\mathbf{v}$ -dependent variable transform:

$\begin{equation} g = \hat{g}+\bar{g} +G(\mathbf{v}),\quad h_a = \hat{h}_a+\bar{h}_a+H_a(\mathbf{v}),\quad h = \hat{h}+\bar{h} +H(\mathbf{v}),\quad u = \hat{u}+\bar{u} +U( \mathbf{v}) \end{equation}$

(3.8)

to translate the time-varying "equilibrium" to the origin. Substituting the above transformation into the equations (3.3)-(3.7), we obtain

$\begin{eqnarray} \frac{d\hat{g}}{dt}+\frac{\partial G}{\partial \mathbf{v}} \mathbf{A}_e \mathbf{v} & = & -a_1(\mathbf{w}) (\hat{g}+G(\mathbf{v}))- \bar{g} (\hat{h}_a+H_a(\mathbf{v}))\\ &&- [ \hat{g}\hat{h}+G(\mathbf{v}) \hat{h}_a+\hat{g}H_a(\mathbf{v})] - G(\mathbf{v})H_a(\mathbf{v})\\ && +[\mathbf{0},\mathbf{C}_d]\mathbf{v} , \end{eqnarray}$

(3.9)

$\begin{eqnarray} \frac{d\hat{h}_a}{dt}+\frac{\partial H_a}{\partial \mathbf{v}} \mathbf{A}_e \mathbf{v} & = & -a_2(\mathbf{w}) (\hat{h}_a+H_a(\mathbf{v}))+a_3(\mathbf{w}) (\hat{h}+H(\mathbf{v})) , \end{eqnarray}$

(3.10)

$\begin{eqnarray} \frac{d\hat{h}}{dt}+\frac{\partial H}{\partial \mathbf{v}} \mathbf{A}_e \mathbf{v} & = & -a_4(\mathbf{w}) (\hat{h}+H(\mathbf{v}))+\hat{u} +U(\mathbf{v}) , \end{eqnarray}$

(3.11)

$\begin{eqnarray} e & = &\hat{g}+G(\mathbf{v}) -[\mathbf{C}_r, \mathbf{0}]\mathbf{v} , \end{eqnarray}$

(3.12)

where $\mathbf{w} = [w_1, w_2, w_3, w_4]$ and

$\begin{eqnarray} a_1(\mathbf{w}) & = & m_1+w_1+ \bar{h}_a, \end{eqnarray}$

(3.13)

$\begin{eqnarray} a_2(\mathbf{w}) & = & m_2 +w_2 , \end{eqnarray}$

(3.14)

$\begin{eqnarray} a_3(\mathbf{w}) & = & m_3 +w_3 , \end{eqnarray}$

(3.15)

$\begin{eqnarray} a_4(\mathbf{w}) & = & m_4 +w_4. \end{eqnarray}$

(3.16)

In order to eliminate the terms that do not contain either $\hat{g}$ or $\hat{h}$ , we set

$\begin{eqnarray} \frac{\partial G}{\partial \mathbf{v}} \mathbf{A}_e \mathbf{v}& = & -a_1(\mathbf{w}) G -\bar{g} H_a- GH_a + [\mathbf{0},\mathbf{C}_d]\mathbf{v} , \end{eqnarray}$

(3.17)

$\begin{eqnarray} \frac{\partial H_a}{\partial \mathbf{v}} \mathbf{A}_e \mathbf{v} & = & -a_2(\mathbf{w}) H_a+a_3(\mathbf{w}) H , \end{eqnarray}$

(3.18)

$\begin{eqnarray} \frac{\partial H}{\partial \mathbf{v}} \mathbf{A}_e \mathbf{v} & = & -a_4(\mathbf{w}) H+U , \end{eqnarray}$

(3.19)

$\begin{eqnarray} 0 & = &G-[\mathbf{C}_r, \mathbf{0}]\mathbf{v}. \end{eqnarray}$

(3.20)

Notice that the equations (3.17)-(3.19) are the result of the center manifold theory (see^[9]). In this way, the tracking problem (3.3)-(3.7) is converted to the following stabilization problem:

$\begin{eqnarray} \frac{d\hat{g}}{dt} & = & -a_1(\mathbf{w}) \hat{g} -\bar{g} \hat{h}_a- ( \hat{g}\hat{h}+G(\mathbf{v}) \hat{h}_a+\hat{g}H_a(\mathbf{v})), \end{eqnarray}$

(3.21)

$\begin{eqnarray} \frac{d\hat{h}_a}{dt} & = & -a_2(\mathbf{w}) \hat{h}_a + a_3(\mathbf{w}) \hat{h} , \end{eqnarray}$

(3.22)

$\begin{eqnarray} \frac{d\hat{h}}{dt} & = & -a_4(\mathbf{w}) \hat{h} + \hat{u} , \end{eqnarray}$

(3.23)

$\begin{eqnarray} e & = &\hat{g} . \end{eqnarray}$

(3.24)

The regulator equations (3.17)-(3.20) can be solved explicitly as follows:

$\begin{eqnarray} G(\mathbf{v}) & = & [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v}, \end{eqnarray}$

(3.25)

$\begin{eqnarray} F(\mathbf{v}) & = & [\mathbf{0},\;\mathbf{C}_d]\mathbf{v}-a_1(\mathbf{w}) [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v} -[\mathbf{C}_r, \; \mathbf{0}] \mathbf{A}_e \mathbf{v} , \end{eqnarray}$

(3.26)

$\begin{eqnarray} H_a(\mathbf{v}) & = &\frac{ F(\mathbf{v}) } {\bar{g}+ [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v}} , \end{eqnarray}$

(3.27)

$\begin{eqnarray} H(\mathbf{v}) & = & \frac{ \frac{\partial H_a}{\partial \mathbf{v}} \mathbf{A}_e \mathbf{v}+a_2(\mathbf{w}) H_a }{a_3(\mathbf{w})} \end{eqnarray}$

(3.28)

$\begin{eqnarray} & = & \frac{ ( \bar{g}+ [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v}) \frac{\partial F}{\partial \mathbf{v}} \mathbf{A}_e \mathbf{v} -F(\mathbf{v})[\mathbf{C}_r, \; \mathbf{0}] \mathbf{A}_e \mathbf{v} } {a_3(\mathbf{w})[\bar{g}+ [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v}]^2} \\ &&+ \frac{ a_2(\mathbf{w})F(\mathbf{v})} {a_3(\mathbf{w})( \bar{g}+ [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v})}, \end{eqnarray}$

(3.29)

$\begin{eqnarray} U(\mathbf{v}) & = & \frac{\partial H}{\partial \mathbf{v}} \mathbf{A}_e \mathbf{v}+a_4(\mathbf{w}) H \\ & = & \frac{ ( [\mathbf{C}_r, \; \mathbf{0}] \mathbf{A}_e \mathbf{v}) \frac{\partial F}{\partial \mathbf{v}} \mathbf{A}_e \mathbf{v} +( \bar{g}+ [\mathbf{C}_r, \; \mathbf{0}] \mathbf{v})\frac{\partial F}{\partial \mathbf{v}} \mathbf{A}_e^2 \mathbf{v} -\frac{\partial F}{\partial \mathbf{v}} \mathbf{A}_e \mathbf{v}[\mathbf{C}_r, \; \mathbf{0}] \mathbf{A}_e \mathbf{v}-F(\mathbf{v})[\mathbf{C}_r, \; \mathbf{0}] \mathbf{A}_e^2 \mathbf{v} } {a_3 (\mathbf{w})[\bar{g}+ [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v}]^2} \\ &&- \frac{ 2 [( \bar{g}+ [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v}) \frac{\partial F}{\partial \mathbf{v}} \mathbf{A}_e \mathbf{v} -F(\mathbf{v})[\mathbf{C}_r, \; \mathbf{0}] \mathbf{A}_e \mathbf{v} ][\bar{g}+ [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v}]a_3(\mathbf{w}) [\mathbf{C}_r, \; \mathbf{0}]\mathbf{A}_e\mathbf{v}} {a_3^2(\mathbf{w})[\bar{g}+ [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v}]^4} \\ &&+ \frac{ a_2(\mathbf{w})\frac{\partial F}{\partial \mathbf{v}} \mathbf{A}_e \mathbf{v} a_3 (\mathbf{w})( \bar{g}+ [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v}) - a_2(\mathbf{w})F(\mathbf{v})a_3 (\mathbf{w}) [\mathbf{C}_r, \; \mathbf{0}]\mathbf{A}_e\mathbf{v})} {a_3^2(\mathbf{w})( \bar{g}+ [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v})^2}\\ &&+a_4(\mathbf{w})\frac{ ( \bar{g}+ [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v}) \frac{\partial F}{\partial \mathbf{v}} \mathbf{A}_e \mathbf{v} -F(\mathbf{v})[\mathbf{C}_r, \; \mathbf{0}] \mathbf{A}_e \mathbf{v} } {a_3(\mathbf{w})[\bar{g}+ [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v}]^2} \\ &&+ \frac{ a_2(\mathbf{w})a_4(\mathbf{w})F(\mathbf{v})} {a_3(\mathbf{w})( \bar{g}+ [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v})}. \end{eqnarray}$

(3.30)

It is easy to see that the equilibrium $(0, 0, 0)$ of the system (3.21)-(3.24) at $\mathbf{v} = \mathbf{0}$ can be stabilized exponentially and locally by a static output feedback

$\hat{u} = k_g\hat{g}.$

Therefore, we have constructed a static output feedback

$\begin{eqnarray} u & = & \bar{u}+\hat{u}+U(\mathbf{v})\\ & = & \bar{u}+k_g\hat{g}+U(\mathbf{v})\\ & = & \bar{u}+k_g(g-\bar{g}- [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v}))+U(\mathbf{v}) \end{eqnarray}$

(3.31)

to regulate the blood glucose to its desired reference $\bar{g}+ [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v}$ . Incorporating this feedback controller into the system (3.3)-(3.7), we obtain a closed-loop system:

$\begin{eqnarray} \frac{dg}{dt} & = & -(m_1+w_1) g- gh_a+ w_1 \bar{g} + G_b+[\mathbf{0},\;\mathbf{C}_d]\mathbf{v}, \end{eqnarray}$

(3.32)

$\begin{eqnarray} \frac{dh_a}{dt} & = & - (m_2 +w_2) h_a+(m_3 +w_3) h+ w_2 \bar{h}_a-w_3 \bar{h} , \end{eqnarray}$

(3.33)

$\begin{eqnarray} \frac{dh}{dt} & = & - (m_4 +w_4) h+w_4 \bar{h}+ \bar{u}+k_g(g-\bar{g}- [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v}))+U(\mathbf{v}) , \end{eqnarray}$

(3.34)

$\begin{eqnarray} \frac{d\mathbf{v} }{dt}& = &\mathbf{A}_{e} \mathbf{v}, \end{eqnarray}$

(3.35)

$\begin{eqnarray} e & = &g- \bar{g}- [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v}. \end{eqnarray}$

(3.36)

We summarize this result in the following theorem.

Theorem 3.1. Suppose that $k_g$ is a constant to exponentially stabilize the equilibrium $(0, 0, 0)$ of the system (3.21)-(3.24). Then the closed-loop nonlinear system (3.32)-(3.36) asymptotically tracks the blood glucose reference $\bar{g} +[\mathbf{C}_r, \; \mathbf{0}] \mathbf{v}(t)$ , that is,

$\lim\limits_{t\to \infty} [g(t)-\bar{g}-[\mathbf{C}_r, \; \mathbf{0}] \mathbf{v}(t)] = 0.$

In numerical simulations below, we use the software Matlab and the experimental data from Sturis et al.'s ^[39]. We use the software Engauge Digitizer 4.1 to read the data.

To conduct a numerical simulation, we use the experimental data from Sturis et al.'s ^[39] to determine $\mathbf{C}_d$ , $\mathbf{C}_r$ , and $\mathbf{A}_e$ . These data were obtained from normal men without a history of diabetes and the detailed experimental methods and procedures are referred to^[39]. Since the data contain 24 hours profiles of glucose and plasma insulin, they are appropriate for testing whether the closed-loop nonlinear system (3.32)-(3.36) asymptotically tracks the blood glucose reference. The data were read from ^[39] by using the software Engauge Digitizer 4.1. According to the second panel of Figure 2 of Sturis et al.'s ^[39], the exogenous glucose infusion shown in Figure 2 is given by

Figure 2. Glucose infusion shown in the second panel of Figure 2 of Sturis et al.'s ^[39]. The unit cc per hour is converted into mM per minute.

DownLoad: Full-Size Img PowerPoint

$J = 125[1+0.33\cdot\sin(\pi t/72 )]\; (\mbox{cc/hr}) .$

Because the glucose concentration unit is mM, we need to covert the unit of cc per hour of the infused solution to the unit of mM per minute of the infusion rate inside the body. According to Sturis et al.'s paper^[39], the infused glucose solution is 20% solution. This could mean that one cc solution contains 0.2 grams of glucose. The blood volume of a person with the weight of 70 kg is about 6 liter. Using the conversion 1 mM = 18 mg/dl, we obtain

$J = \frac{125\cdot 0.2\cdot 1000}{60\cdot 18\cdot 60} [1+0.33\cdot\sin(\pi t/72 )] = 0.3858 [1+0.33\cdot\sin(\pi t/72 )]\; (\mbox{mM/min}) .$

In this simulation, we find that $G_b = 0.3218$ mM/min. Thus we obtain

$\mathbf{C}_d = [0.3858 - G_b, 0.1273, 0 ],$

and

$\mathbf{A}_{d} = \left[\begin{array}{ccc } 0&0&0 \\ 0& 0& \pi /72 \\ 0&-\pi /72 &0 \end{array} \right] .$

Fitting the experimental glucose data shown in the second panel of Figure 2 of Sturis et al.'s ^[39] into the Fourier polynomial:

$g = a_0+\sum\limits_{i = 1}^{17}[a_i\sin(i\pi t/750 )+b_i\cos(i\pi t/750)],$

we obtain

$\begin{eqnarray*} a_0 & = & 7.599 , \\ a_1 & = & -0.2394 , \quad b_1 = -1.024,\\ a_2 & = & 0.328 , \quad b_2 = -0.2217,\\ a_3 & = & -0.2314 , \quad b_3 = -0.1684 ,\\ a_4 & = & -0.1663 , \quad b_4 = -0.05297,\\ a_5 & = & -0.1167 , \quad b_5 = 0.1937,\\ a_6 & = & 0.07061 , \quad b_6 = -0.1443,\\ a_7 & = & 0.3512 , \quad b_7 = -0.08337,\\ a_8 & = & 0.02829 , \quad b_8 = -0.0705,\\ a_9 & = & 0.2508 , \quad b_9 = -0.259,\\ a_{10} & = & 0.6841 , \quad b_{10} = -0.8653,\\ a_{11} & = & -0.3087 , \quad b_{11} = 0.5024,\\ a_{12} & = & -0.2117 , \quad b_{12} = 0.2057,\\ a_{13} & = & -0.07391 , \quad b_{13} = 0.02245,\\ a_{14} & = & -0.01738 , \quad b_{14} = -0.1047,\\ a_{15} & = & -0.1031 , \quad b_{15} = -0.05237,\\ a_{16} & = & -0.001116 , \quad b_{16} = 0.004296,\\ a_{17} & = & -0.04128 , \quad b_{17} = 0.003007. \end{eqnarray*}$

The unit of the coefficients $a_0$ , $a_i$ and $b_i$ is mM, where $i = 1, 2, \cdots, 17$ . The glucose equilibrium is set to $\bar{g} = 6$ mM. Thus we obtain

$\mathbf{C}_r = [a_0 - \bar{g}, a_1, b_1 , a_2, b_2 , \cdots a_{17}, b_{17} ],$

and

$\mathbf{A}_{r} = \left[\begin{array}{cccccccc} 0&0&0&0&0& \cdots& 0&0\\ 0& 0& \pi/750& 0 & 0& \cdots& 0&0 \\ 0&-\pi/750 &0&0&0& \cdots& 0&0\\ 0&0&0&0&2\pi/750& \cdots& 0&0\\ 0&0&0&-2\pi/750&0& \cdots& 0&0\\ \vdots&\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\ 0&0&0&0&0& \cdots& 0&17\pi/750\\ 0&0&0&0&0& \cdots& -17\pi/750&0 \end{array} \right] .$

Figure 3 shows that the data can be fitted into this Fourier polynomial very well.

Figure 3. Fourier polynomial fitting of the blood glucose data shown in the second panel of Figure 2 of Sturis et al.'s ^[39]. The data were read by using the software Engauge Digitizer 4.1.

DownLoad: Full-Size Img PowerPoint

The estimation of the parameters $m_1, m_2, m_3$ , and $m_4$ is challenging. Because the expression of the controller $u$ is too big to be expressed in Matlab, we failed to use the Matlab function "sbioparamestim" to estimate them. Thus, starting with the values from Bergman et al.'s paper^[5]: $m_1 = 0.0296$ (/min), $m_2 = 0.0186$ (/min), $m_3 = 0.0000065$ (/( $\mu$ U/ml)/min $^2$ ), $m_4 = 0.238$ (/min), we manually adjust them to achieve as a best result as we can. These estimated values are presented in . Compared with Bergman et al.'s original values, only $m_1$ has a large change. The uncertainties $w_1, w_2, w_3,$ and $w_4$ are taken to be 10% of their nominal values. The proportional constant $k_g$ is taken to be $0.01$ (/min).

The model is solved numerically with Matlab. We first run the model 1500 minutes to achieve a stable state and then set that time as the initial time. shows that the blood glucose $g$ is asymptotically tracking the desired time-varying reference determined by the experimental data from Sturis et al.'s ^[39]. Although the simulated insulin profile does not exactly match with the experimental data, the simulated and experimental pulses occur simultaneously. The dynamics of the controller $u$ and insulin effect $h_a$ follows the dynamics of the glucose infusion.

Figure 4. Glucose tracking under the static output feedback controller (3.31). The blood glucose

$g$ is tracking the desired time-varying reference determined by the experimental data from Sturis et al.'s ^[39]. The data were read by using the software Engauge Digitizer 4.1.

DownLoad: Full-Size Img PowerPoint

Because the static output feedback (3.31) contains the uncertainties $\mathbf{w}$ , it is not robust. If the uncertainties $\mathbf{w}$ are set to 0 in the static output feedback controller (3.31), shows that the glucose $g$ fails to track its reference. In fact, it has been shown that there are no static state feedback controls that can solve the robust tracking problem (see^[19]).

Figure 5. Glucose fails to track its reference when the uncertainties are set to 0 in the static output feedback controller (3.31).

DownLoad: Full-Size Img PowerPoint

4. A robust dynamical feedback and feedforward controller

We now design a robust dynamical feedback and feedforward controller via linearization. The linearization of the system (3.3)-(3.7) at its equilibrium is given by

$\begin{eqnarray} \frac{d\tilde{g}}{dt} & = & -a_1(\mathbf{w}) \tilde{g} - \bar{g} \tilde{h}_a +[\mathbf{0},\;\mathbf{C}_d]\mathbf{v}, \end{eqnarray}$

(4.1)

$\begin{eqnarray} \frac{d\tilde{h}_a}{dt} & = & - a_2(\mathbf{w}) \tilde{h}_a + a_3(\mathbf{w}) \tilde{h} , \end{eqnarray}$

(4.2)

$\begin{eqnarray} \frac{d\tilde{h}}{dt} & = & - a_4(\mathbf{w}) \tilde{h} +\tilde{u} , \end{eqnarray}$

(4.3)

$\begin{eqnarray} \frac{d\mathbf{v} }{dt}& = &\mathbf{A}_{e} \mathbf{v}, \end{eqnarray}$

(4.4)

$\begin{eqnarray} e & = &\tilde{g}- [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v}, \end{eqnarray}$

(4.5)

where

$\begin{equation} \tilde{g} = g-\bar{g},\; \tilde{h}_a = h_a-\bar{h}_a,\; \tilde{h} = h-\bar{h},\; \tilde{u} = u-\bar{u}. \end{equation}$

(4.6)

To construct a robust control law, we introduce the following dynamic compensator:

$\begin{equation} \frac{d\mathbf{z} }{dt} = \mathbf{G}_1 \mathbf{z}+\mathbf{G}_2(\tilde{g} - [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v}), \end{equation}$

(4.7)

where the matrix pair $(\mathbf{G}_1, \mathbf{G}_2)$ is determined below. Since the state variables $\tilde{h}_a$ and $\tilde{h}$ are not available for feedback, we also need to use the Luenburger observer:

$\begin{equation} \frac{d \mathbf{o}}{dt} = \mathbf{A} \mathbf{o}+ \mathbf{B}\tilde{u}+\mathbf{L}(\tilde{g}-\mathbf{C}\mathbf{o}) , \end{equation}$

(4.8)

where the vector $\mathbf{L}$ is to be designed and

$\begin{equation} \mathbf{A} = \left[\begin{array}{ccc} -a_1(\mathbf{0})& -\bar{g} &0\\ 0& -a_2(\mathbf{0}) &a_3(\mathbf{0})\\ 0& 0&-a_4(\mathbf{0}) \end{array} \right],\quad \mathbf{B} = \left[\begin{array}{c} 0\\ 0\\ 1 \end{array} \right],\quad \mathbf{C} = [1, 0, 0]. \end{equation}$

(4.9)

We then introduce the following dynamic output feedback and feedforward controller:

$\begin{equation} \tilde{u} = -\mathbf{K}_z \mathbf{z} -\mathbf{K}_o \mathbf{o}, \\ \end{equation}$

(4.10)

where the vectors $\mathbf{K}_o$ and $\mathbf{K}_z$ are to be designed.

A minimal polynomial of a matrix $\mathbf{A}$ is a polynomial $P(\lambda)$ with the least degree such that $P(\mathbf{A}) = \mathbf{0}$ . If $\omega_r = \omega_d$ , then the minimal polynomial of $\mathbf{A}_{e}$ is given by

$\begin{eqnarray*} p(\lambda) & = &\lambda(\lambda^2+\omega_d^2)(\lambda^2+(2\omega_d)^2)\cdots (\lambda^2+(l\omega_d)^2)\\ & = &\lambda^{2l+1}+c_1\lambda^{2l}+\cdots+c_{2l}\lambda +c_{2l+1}, \end{eqnarray*}$

where $l = \max(m, n)$ . If $\omega_r\neq\omega_d$ , then the minimal polynomial of $\mathbf{A}_{e}$ is given by

$\begin{eqnarray*} p(\lambda) & = &\lambda(\lambda^2+\omega_d^2)(\lambda^2+(2\omega_d)^2)\cdots (\lambda^2+(n\omega_d)^2)(\lambda^2+\omega_r^2)(\lambda^2+(2\omega_r)^2)\cdots (\lambda^2+(m\omega_r)^2)\\ & = &\lambda^{2l+1}+c_1\lambda^{2l}+\cdots+c_{2l}\lambda +c_{2l+1}, \end{eqnarray*}$

where $l = m+n$ . Then $\mathbf{G}_1$ is defined to be a companion matrix of the minimal polynomial and $\mathbf{G}_2$ is set to be a matrix such that the pair $(\mathbf{G}_1, \mathbf{G}_2)$ is controllable as follows:

$\begin{equation} \mathbf{G}_1 = \left[\begin{array}{ccccc} 0&1&0&\cdots&0\\ 0& 0& 1& \cdots & 0\\ \vdots&\vdots &\ddots&\vdots&\vdots\\ 0&0&0&\cdots&1\\ -c_{2l+1}&-c_{2l} &-c_{2l-1}&\cdots&-c_1 \end{array} \right],\quad \mathbf{G}_2 = \left[\begin{array}{c} 0\\ 0 \\ \vdots\\ 0\\ 1 \end{array} \right]. \end{equation}$

(4.11)

In control theory, such a matrix pair $(\mathbf{G}_1, \mathbf{G}_2)$ is called to incorporate a 1-copy internal model of the matrix $\mathbf{A}_{e}$ . In the following proof, we can see the reason why the pair $(\mathbf{G}_1, \mathbf{G}_2)$ is set up this way.

To convert this tracking problem to a stabilization problem, we introduce the coordinate transformation

$\begin{equation} \tilde{g} = \hat{g} +\mathbf{G} \mathbf{v},\quad \tilde{h}_a = \hat{h}_a +\mathbf{H}_a \mathbf{v},\quad \tilde{h} = \hat{h} +\mathbf{H} \mathbf{v},\quad \quad \mathbf{z} = \hat{\mathbf{z}} +\mathbf{Z} \mathbf{v} ,\quad \quad \mathbf{o} = \hat{\mathbf{o}} +\mathbf{O} \mathbf{v}. \end{equation}$

(4.12)

Substituting this transformation (4.12) and the feedback (4.10) into the composite system (4.1)-(4.5) and the equations (4.7) and (4.8), we obtain

$\begin{eqnarray*} \frac{d\hat{g}}{dt} + \mathbf{G} \mathbf{A}_{e}\mathbf{v}& = & - a_1(\mathbf{w}) (\hat{g} +\mathbf{G} \mathbf{v})- \bar{g} (\hat{h}_a +\mathbf{H}_a \mathbf{v}) +[\mathbf{0},\;\mathbf{C}_d]\mathbf{v}, \\ \frac{d\hat{h}_a}{dt}+ \mathbf{H}_a \mathbf{A}_{e}\mathbf{v} & = & - a_2(\mathbf{w}) (\hat{h}_a +\mathbf{H}_a \mathbf{v}) + a_3(\mathbf{w}) (\hat{h} +\mathbf{H} \mathbf{v}),\\ \frac{d\hat{h}}{dt}+ \mathbf{H} \mathbf{A}_{e}\mathbf{v} & = & - a_4(\mathbf{w}) (\hat{h} +\mathbf{H} \mathbf{v}) - \mathbf{K}_z (\hat{\mathbf{z}}+\mathbf{Z}\mathbf{v})- \mathbf{K}_o (\hat{\mathbf{o}}+\mathbf{O} \mathbf{v}) , \\ \frac{d\hat{\mathbf{z}}}{dt}+ \mathbf{Z} \mathbf{A}_{e}\mathbf{v} & = & \mathbf{G}_1 (\hat{\mathbf{z}}+\mathbf{Z} \mathbf{v})+ \mathbf{G}_2(\hat{g}+\mathbf{G} \mathbf{v} -[\mathbf{C}_r, \; \mathbf{0}]\mathbf{v}),\\ \frac{d\hat{\mathbf{o}}}{dt}+ \mathbf{O} \mathbf{A}_{e}\mathbf{v} & = & \mathbf{A} (\hat{\mathbf{o}}+\mathbf{O} \mathbf{v})- \mathbf{B}[\mathbf{K}_z (\hat{\mathbf{z}}+\mathbf{Z}\mathbf{v})+ \mathbf{K}_o (\hat{\mathbf{o}}+\mathbf{O} \mathbf{v}) ]+\mathbf{L}(\hat{g} +\mathbf{G} \mathbf{v}-\mathbf{C}(\hat{\mathbf{o}}+\mathbf{O} \mathbf{v})),\\ e & = &\hat{g} +\mathbf{G} \mathbf{v} - [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v}. \end{eqnarray*}$

In order to eliminate $\mathbf{v}$ from the above system, we set

$\begin{eqnarray} \mathbf{G} \mathbf{A}_{e} & = & - a_1(\mathbf{w}) \mathbf{G} - \bar{g} \mathbf{H}_a +[\mathbf{0},\;\mathbf{C}_d], \end{eqnarray}$

(4.13)

$\begin{eqnarray} \mathbf{H}_a \mathbf{A}_{e} & = & - a_2(\mathbf{w}) \mathbf{H}_a +a_3(\mathbf{w}) \mathbf{H} , \end{eqnarray}$

(4.14)

$\begin{eqnarray} \mathbf{H} \mathbf{A}_{e} & = & - a_4(\mathbf{w}) \mathbf{H} -\mathbf{K}_z \mathbf{Z} -\mathbf{K}_o \mathbf{O}, \end{eqnarray}$

(4.15)

$\begin{eqnarray} \mathbf{Z} \mathbf{A}_{e} & = & \mathbf{G}_1 \mathbf{Z} + \mathbf{G}_2( \mathbf{G} -[\mathbf{C}_r, \; \mathbf{0}] ), \end{eqnarray}$

(4.16)

$\begin{eqnarray} \mathbf{O} \mathbf{A}_{e} & = & \mathbf{A} \mathbf{O} -\mathbf{B}\mathbf{K}_z \mathbf{Z} -\mathbf{B}\mathbf{K}_o \mathbf{O}+\mathbf{L}(\mathbf{G}-\mathbf{C}\mathbf{O}), \end{eqnarray}$

(4.17)

$\begin{eqnarray} 0 & = & \mathbf{G} - [\mathbf{C}_r, \; \mathbf{0}] \end{eqnarray}$

(4.18)

and then

$\begin{eqnarray} \frac{d\hat{g}}{dt} & = & - a_1(\mathbf{w}) \hat{g}- \bar{g} \hat{h}_a , \end{eqnarray}$

(4.19)

$\begin{eqnarray} \frac{d\hat{h}_a}{dt} & = & - a_2(\mathbf{w}) \hat{h}_a + a_3(\mathbf{w}) \hat{h} , \end{eqnarray}$

(4.20)

$\begin{eqnarray} \frac{d\hat{h}}{dt} & = & - a_4(\mathbf{w}) \hat{h} - \mathbf{K}_z \hat{\mathbf{z}} - \mathbf{K}_o \hat{\mathbf{o}} , \end{eqnarray}$

(4.21)

$\begin{eqnarray} \frac{d\hat{\mathbf{z}}}{dt} & = & \mathbf{G}_1 \hat{\mathbf{z}} + \mathbf{G}_2 \hat{g} , \end{eqnarray}$

(4.22)

$\begin{eqnarray} \frac{d\hat{\mathbf{o}}}{dt} & = & \mathbf{A} \hat{\mathbf{o}}- \mathbf{B}[\mathbf{K}_z \hat{\mathbf{z}}+ \mathbf{K}_o \hat{\mathbf{o}} ]+\mathbf{L}(\hat{g} -\mathbf{C}\hat{\mathbf{o}}), \end{eqnarray}$

(4.23)

$\begin{eqnarray} e & = &\hat{g} . \end{eqnarray}$

(4.24)

If the regulator equations (4.13)-(4.18) have a solution, then the tracking problem is converted to the problem of stabilization of the system (4.19)-(4.24).

We first show that the equilibrium $(0, 0, 0, \mathbf{0}, \mathbf{0})$ of the system (4.19)-(4.24) can be exponentially stabilized at $\mathbf{w} = \mathbf{0}$ and then robustly stabilized for sufficiently small $\mathbf{w}$ . For this, we set $\hat{\mathbf{x}} = [\hat{g}, \hat{h}_a, \hat{h}]^T$ and rewrite the system (4.19)-(4.24) at $\mathbf{w} = \mathbf{0}$ as follows:

$\begin{eqnarray} \frac{d\hat{\mathbf{x}}}{dt} & = & \mathbf{A} \hat{\mathbf{x}}- \mathbf{B}[\mathbf{K}_z \hat{\mathbf{z}}+ \mathbf{K}_o \hat{\mathbf{o}} ], \end{eqnarray}$

(4.25)

$\begin{eqnarray} \frac{d\hat{\mathbf{z}}}{dt} & = & \mathbf{G}_1 \hat{\mathbf{z}} + \mathbf{G}_2 \mathbf{C}\hat{\mathbf{x}} , \end{eqnarray}$

(4.26)

$\begin{eqnarray} \frac{d\hat{\mathbf{o}}}{dt} & = & \mathbf{A} \hat{\mathbf{o}}- \mathbf{B}[\mathbf{K}_z \hat{\mathbf{z}}+ \mathbf{K}_o \hat{\mathbf{o}} ]+\mathbf{L}(\mathbf{C}\hat{\mathbf{x}} -\mathbf{C}\hat{\mathbf{o}}). \end{eqnarray}$

(1.27)

Setting $\hat{\mathbf{q}} = \hat{\mathbf{x}}-\hat{\mathbf{o}}$ , we derive from the above system that

$\begin{eqnarray} \frac{d\hat{\mathbf{x}}}{dt} & = & ( \mathbf{A} - \mathbf{B}\mathbf{K}_o)\hat{\mathbf{x}}- \mathbf{B}[\mathbf{K}_z \hat{\mathbf{z}}- \mathbf{K}_o \hat{\mathbf{q}} ], \end{eqnarray}$

(4.28)

$\begin{eqnarray} \frac{d\hat{\mathbf{z}}}{dt} & = & \mathbf{G}_1 \hat{\mathbf{z}} + \mathbf{G}_2 \mathbf{C}\hat{\mathbf{x}} , \end{eqnarray}$

(4.29)

$\begin{eqnarray} \frac{d\hat{\mathbf{q}}}{dt} & = & (\mathbf{A} -\mathbf{L}\mathbf{C})\hat{\mathbf{q}}, \end{eqnarray}$

(4.30)

and then

$\begin{equation} \frac{d}{dt}\left[\begin{array}{c} \hat{\mathbf{x}} \\ \hat{\mathbf{z}}\\ \hat{\mathbf{q}} \end{array} \right] = \left[\begin{array}{ccc} \mathbf{A} - \mathbf{B}\mathbf{K}_o& - \mathbf{B}\mathbf{K}_z & \mathbf{B}\mathbf{K}_o\\ \mathbf{G}_2 \mathbf{C}&\mathbf{G}_1 &\mathbf{0}\\ \mathbf{0}&\mathbf{0} &\mathbf{A} -\mathbf{L}\mathbf{C} \end{array} \right]\left[\begin{array}{c} \hat{\mathbf{x}} \\ \hat{\mathbf{z}}\\ \hat{\mathbf{q}} \end{array} \right]. \end{equation}$

(4.31)

It is easy to see that $(\mathbf{A}, \mathbf{C})$ is observable and the pair

$\begin{equation} \mathbf{\Phi} = \left[\begin{array}{cc} \mathbf{A} & \mathbf{0} \\ \mathbf{G}_2 \mathbf{C}&\mathbf{G}_1 \end{array} \right],\quad \mathbf{\Psi} = \left[\begin{array}{ccc} \mathbf{B} \\ \mathbf{0} \end{array} \right] \end{equation}$

(4.32)

is stabilizable (see, e.g., ^[19]). So there exist constant vectors $\mathbf{K}_o, \mathbf{K}_z, \mathbf{L}$ such that all eigenvalues of the coefficient matrix of the system

$\begin{equation} \mathbf{\Omega} = \left[\begin{array}{ccc} \mathbf{A} - \mathbf{B}\mathbf{K}_o& - \mathbf{B}\mathbf{K}_z & \mathbf{B}\mathbf{K}_o\\ \mathbf{G}_2 \mathbf{C}&\mathbf{G}_1 &\mathbf{0}\\ \mathbf{0}&\mathbf{0} &\mathbf{A}-\mathbf{L}\mathbf{C} \end{array} \right] \end{equation}$

(4.33)

have negative real parts. That is, the equilibrium $(0, 0, 0, \mathbf{0}, \mathbf{0})$ of the system (4.19)-(4.24) can be exponentially stabilized at $\mathbf{w} = \mathbf{0}$ .

So it remains to show that the regulator equations (4.13)-(4.18) have a solution. To this end, we set $\mathbf{X} = [\mathbf{G}, \mathbf{H}_a, \mathbf{H}]^T$ and rewrite the equations (4.13)-(4.17) as one matrix equation:

$\begin{equation} \left[\begin{array}{c } \mathbf{X} \\ \mathbf{Z}\\ \mathbf{Z} \end{array} \right] \mathbf{A}_e = \mathbf{Q}\left[\begin{array}{c } \mathbf{X} \\ \mathbf{Z} \\ \mathbf{O} \end{array} \right]+ \left[\begin{array}{c } \mathbf{C}^T [0,\;\mathbf{C}_d] \\ \mathbf{G}_2[\mathbf{C}_r, \; 0]\\ \mathbf{0} \end{array} \right], \end{equation}$

(4.34)

where

$\begin{equation} \mathbf{Q} = \left[\begin{array}{ccc} \mathbf{A} & - \mathbf{B}\mathbf{K}_z & -\mathbf{B}\mathbf{K}_o\\ \mathbf{G}_2 \mathbf{C}&\mathbf{G}_1 &\mathbf{0}\\ \mathbf{L}\mathbf{C} &- \mathbf{B}\mathbf{K}_z &\mathbf{A} -\mathbf{L}\mathbf{C}-\mathbf{B}\mathbf{K}_o \end{array} \right]. \end{equation}$

(4.35)

Since $\mathbf{Q}$ and $\mathbf{\Omega}$ are equivalent, all eigenvalues of $\mathbf{Q}$ also have negative real parts. Because all eigenvalues of $\mathbf{A}_e$ have zero real parts, the above matrix equation has a unique solution (see^[19]). To show the solution also satisfies the equation (4.18), we let $\mathbf{Z} = [\mathbf{Z}_{1}, \mathbf{Z}_{2}, \cdots, ,\mathbf{Z}_{2l+1}]^T$ . It then follows from the equation (4.16) that

$\begin{eqnarray} \mathbf{Z}_{2} & = &\mathbf{Z}_{1} \mathbf{A}_e, \end{eqnarray}$

(4.36)

$\begin{eqnarray} \mathbf{Z}_{3} & = &\mathbf{Z}_{2} \mathbf{A}_e = \mathbf{Z}_{1} \mathbf{A}_e^2, \end{eqnarray}$

(4.37)

$\begin{eqnarray} &\vdots& \\ \mathbf{Z}_{2l+1} & = &\mathbf{Z}_{2l} \mathbf{A}_e = \mathbf{Z}_{1} \mathbf{A}_e^{2l}, \end{eqnarray}$

(4.38)

$\begin{eqnarray} \mathbf{G} & = &[\mathbf{C}_r, \; 0]+\mathbf{Z}_{2l+1} \mathbf{A}_e+c_{2l+1}\mathbf{Z}_{1}+c_{2l}\mathbf{Z}_{2}+\cdots+c_1\mathbf{Z}_{2l+1} \\ & = &[\mathbf{C}_r, \; 0]+\mathbf{Z}_{1} \mathbf{A}_e^{2l+1} +c_{2l+1}\mathbf{Z}_{1}+c_{2l}\mathbf{Z}_{1}\mathbf{A}_e+\cdots +c_1\mathbf{Z}_{1} \mathbf{A}_e^{2l} \\ & = &[\mathbf{C}_r, \; 0]+\mathbf{Z}_{1} ( c_{2l+1}\mathbf{I}+c_{2l} \mathbf{A}_e+\cdots +c_1 \mathbf{A}_e^{2l}+\mathbf{A}_e^{2l+1}) \\ & = &[\mathbf{C}_r, \; 0]+\mathbf{Z}_{1}\mathbf{0} \\ & = &[\mathbf{C}_r, \; 0] . \end{eqnarray}$

(4.39)

Here we have used the property that $p(\lambda)$ is the minimal polynomial of $\mathbf{A}_e$ , that is, $p(\mathbf{A}_e) = \mathbf{0}$ . The above proof of the equation (4.18) gives the reason why the pair $(\mathbf{G}_1, \mathbf{G}_2)$ is constructed as in (4.11).

Using the translation (4.6), we obtain a dynamic output feedback and feedforward controller for the original nonlinear system (3.3)-(3.7):

$\begin{eqnarray} u& = &\bar{u}+\tilde{u}\\ & = & \bar{u}-\mathbf{K}_z \mathbf{z}-\mathbf{K}_o \mathbf{o}, \end{eqnarray}$

(4.40)

$\begin{eqnarray} \frac{d\mathbf{z} }{dt} & = &\mathbf{G}_1 \mathbf{z}+\mathbf{G}_2( g-\bar{g} - [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v}), \end{eqnarray}$

(4.41)

$\begin{eqnarray} \frac{d \mathbf{o}}{dt} & = & \mathbf{A} \mathbf{o}- \mathbf{B}(\mathbf{K}_z \mathbf{z}+\mathbf{K}_o \mathbf{o})+\mathbf{L}(g-\bar{g}-\mathbf{C}\mathbf{o}). \end{eqnarray}$

(4.42)

For the convenience of reference, we collect all relevant equations together to form a closed-loop system:

$\begin{eqnarray} \frac{dg}{dt} & = & -(m_1+w_1) g- gh_a+ w_1 \bar{g} +G_b+[\mathbf{0},\;\mathbf{C}_d]\mathbf{v}, \end{eqnarray}$

(4.43)

$\begin{eqnarray} \frac{dh_a}{dt} & = & - (m_2 +w_2) h_a+ (m_3 +w_3) h +w_2 \bar{h}_a-w_3 \bar{h} , \end{eqnarray}$

(4.44)

$\begin{eqnarray} \frac{dh}{dt} & = & - (m_4 +w_4) h+w_4 \bar{h}+\bar{u}-\mathbf{K}_z \mathbf{z}-\mathbf{K}_o \mathbf{o}, \end{eqnarray}$

(4.45)

$\begin{eqnarray} \frac{d\mathbf{z} }{dt} & = &\mathbf{G}_1 \mathbf{z}+\mathbf{G}_2( g-\bar{g} - [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v}), \end{eqnarray}$

(4.46)

$\begin{eqnarray} \frac{d \mathbf{o}}{dt} & = & \mathbf{A} \mathbf{o}- \mathbf{B}(\mathbf{K}_z \mathbf{z}+\mathbf{K}_o \mathbf{o})+\mathbf{L}(g-\bar{g}-\mathbf{C}\mathbf{o}), \end{eqnarray}$

(4.47)

$\begin{eqnarray} \frac{d\mathbf{v} }{dt}& = &\mathbf{A}_{e} \mathbf{v}, \end{eqnarray}$

(4.48)

$\begin{eqnarray} e & = &g- \bar{g}- [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v}. \end{eqnarray}$

(4.49)

Theorem 4.1. Suppose that $\mathbf{K}_o, \mathbf{K}_z, \mathbf{L}$ are constant vectors such that the matrix

$\begin{equation*} \mathbf{\Omega} = \left[\begin{array}{ccc} \mathbf{A} - \mathbf{B}\mathbf{K}_o& - \mathbf{B}\mathbf{K}_z & \mathbf{B}\mathbf{K}_o\\ \mathbf{G}_2 \mathbf{C}&\mathbf{G}_1 &\mathbf{0}\\ \mathbf{0}&\mathbf{0} &\mathbf{A} -\mathbf{L}\mathbf{C} \end{array} \right] \end{equation*}$

is exponentially stable. Then the blood glucose of the closed-loop nonlinear system (4.43)-(4.49) asymptotically tracks its reference $\bar{g} + [\mathbf{C}_r, \; \mathbf{0}] \mathbf{v}(t)$ for any sufficiently small uncertainty $\mathbf{w}$ , that is,

$\begin{equation*} \lim\limits_{t\to \infty} [g(t)-\bar{g}- [\mathbf{C}_r, \; \mathbf{0}] \mathbf{v}(t)] = 0. \end{equation*}$

We use the same data as used in the above section to conduct a numerical simulation for this robust approximate controller (4.40). In this simulation, we first used an internal model 1-copy $\mathbf{G}_1$ of the matrix

$\begin{equation*} \mathbf{A}_{e} = \left[\begin{array}{ccccccccc} 0&0&0& \cdots& 0&0&0&0&0\\ 0& 0& \pi/750& \cdots& 0&0&0&0&0 \\ 0&-\pi/750 &0& \cdots& 0&0&0&0&0\\ 0&0&0& \cdots& 0&0&0&0&0\\ 0&0&0& \cdots& 0&0&0&0&0\\ \vdots&\vdots&\vdots& \ddots&\vdots&\vdots&\vdots&\vdots&\vdots\\ 0&0&0& \cdots& 0&17\pi/750&0&0&0\\ 0&0&0& \cdots& -17\pi/750&0&0&0&0\\ 0&0&0& \cdots& 0&0&0&0&0\\ 0&0&0& \cdots& 0&0&0&0&\pi/72\\ 0&0&0& \cdots& 0&0&0&-\pi/72&0 \end{array} \right] . \end{equation*}$

Unfortunately, Matlab failed to solve the closed-loop control system (4.43)-(4.49) with this high dimensional $\mathbf{G}_1$ . Thus we have to reduce the dimension to obtain an internal model 1-copy $\mathbf{G}_1$ of the matrix

$\begin{equation*} \tilde{\mathbf{A}}_{e} = \left[\begin{array}{cccccccccc} 0& 0& 0& 0& 0& 0&0&0&0&0\\ 0& 0& \pi/750& 0& 0& 0&0&0&0&0 \\ 0&-\pi/750 &0& 0& 0& 0&0&0&0&0\\ 0& 0& 0& 0& 2\pi/750&0&0&0&0&0\\ 0& 0& 0& - 2\pi/750&0 & 0&0&0&0&0\\ 0& 0& 0& 0& 0&0&3\pi/750&0&0&0\\ 0& 0& 0& 0& 0&-3\pi/750&0&0&0&0\\ 0&0&0& 0& 0&0&0&0&0&0\\ 0&0&0& 0& 0&0&0&0&0&\pi/72\\ 0&0&0& 0& 0&0&0&0&-\pi/72&0 \end{array} \right] . \end{equation*}$

In this case, the matrices $\mathbf{G}_1$ and $\mathbf{G}_2$ are given by

$\begin{eqnarray*} \mathbf{G}_1 & = & \left[\begin{array}{ccccccccc} 0& 1& 0& 0& 0& 0&0&0&0\\ 0& 0& 1& 0& 0& 0&0&0&0\\ 0& 0& 0& 1& 0& 0&0&0&0\\ 0& 0& 0& 0& 1& 0&0&0&0\\ 0& 0& 0& 0& 0& 1&0&0&0\\ 0& 0& 0& 0& 0& 0&1 &0&0\\ 0& 0& 0& 0& 0& 0&0&1&0\\ 0& 0& 0& 0& 0& 0&0&0&1\\ 0 &-3.70\cdot 10^{-16} & 0& -2.89\cdot 10^{-11} & 0 &-4.83\cdot 10^{-7} & 0 & -0.0021495 & 0 \end{array} \right] ,\\ \mathbf{G}_2 & = & \left[\begin{array}{c} 0\\ 0 \\ 0\\ 0\\ 0\\ 0 \\ 0\\ 0\\ 1 \end{array} \right]. \end{eqnarray*}$

The vectors $\mathbf{K}_o$ and $\mathbf{K}_z$ are determined in the way such that the eigenvalues of the matrix $\mathbf{\Phi}- \mathbf{\Psi}\mathbf{K}$ are: $-0.0500, -0.0909, -0.1318, -0.1727, -0.2136, -0.2545, -0.2955, -0.3364,$ $-0.3773, -0.4182, -0.4591, -0.5000,$ where $\mathbf{K} = [\mathbf{K}_o, \mathbf{K}_z]$ . Using the algorithm for computing the feedback gain matrix from Section 3.4 of the book^[26], we obtain

$\begin{eqnarray*} \mathbf{K}_o & = & [ -55379.84, 379273.25, 2.93],\\ \mathbf{K}_z & = & [ -0.00028, -0.019, -0.56, -9.41, -100.48, -728.94, -3701.22, -13180, -33587.83]. \end{eqnarray*}$

The vector $\mathbf{L}$ is determined in the way such that the eigenvalues of the matrix $\mathbf{A}- \mathbf{L}\mathbf{C}$ are: $-0.01, -0.02, -0.03$ . Using the Matlab function place, we obtain

$\begin{equation*} \mathbf{L } = [ -0.30679, -0.011386,211.93]^T. \end{equation*}$

According to our simulation studies, the eigenvalues should be between $-0.5$ and $-0.01$ to achieve a good performance of the controller. All other parameter values are the same as in the previous simulations.

Comparing Figure 6 with Figure 4, we found it interesting that the dynamics of all state variables under the robust controller (4.40) is identical to the one under the static controller (3.31).

Figure 6. Glucose tracking under the robust controller (4.40). The dynamics of the blood glucose tracking is almost identical to the one under the static controller (3.31).

DownLoad: Full-Size Img PowerPoint

The static controller (3.31) and the robust dynamic controller (4.40) have their own advantages and disadvantages. The static controller (3.31) is easy to implement, but is not flexible, no parameters being used to adjust to improve its performance. Since it is not robust, a small change in the model parameters results in a deviation of blood glucose from its target reference. On the other hand, the robust dynamical controller (4.40) can tolerate certain small changes in the model parameters and its performance can be improved by adjusting certain controller parameters such as $\mathbf{K}_o, \mathbf{K}_z$ and $\mathbf{L}$ . However, it is more difficult to implement because the dimension of the closed-loop control system is greatly increased through the dynamic compensator (4.7). As discussed above, Matlab failed to solve the closed-loop control system numerically if the original $\mathbf{G}_1$ for the matrix $\mathbf{A}_e$ is used. Furthermore, its performance critically depends on the gain vectors $\mathbf{K}_o, \mathbf{K}_z$ and $\mathbf{L}$ . It is a great challenging problem to determine them to achieve the best performance of the controller in a way such that the simulated both glucose and insulin dynamics agrees with the real experimental data. Since the controller performance issue itself is an independent research area^[12,22,28], this problem is beyond the scope of this work.

5. Discussion and conclusion

In this work, we have solved a robust glucose tracking problem. Parameter uncertainties were introduced into the mathematical model proposed by Bergman et al.^[5]. Using an actual blood glucose level as feedback and an exogenous glucose input and a desired glucose reference as feedforward, we have designed a robust feedback and feedforward controller such that the blood glucose tracks the desired time-varying glucose reference for any small uncertainties. Numerical simulations conducted with published experimental blood glucose data show that our controller is effective in achieving the blood glucose tracking and agree with our theoretical results.

Limitations exist in this work. Although our robust feedback and feedforward controllers might have a potential application in the insulin pump technology, more comprehensive tests and improvements with experimental and clinic data are needed to be done before they might be implemented in an insulin pump. The interaction between glucose and insulin is complex. Our current simple model missed considering many other important control mechanisms such as the glucose regulation by the hormone glucagon and the organ liver. Therefore we may need to consider more complex mathematical models, such as the model proposed by Sturis et al.^[40] and the model proposed by Liu et al.^[24], and introduce parameter uncertainties into them. In all these cases, the regulator equations are more complex and it might be impossible to solve them analytically. Then a different approach such as used in^[3] is needed to design a controller.

Acknowledgments

The author heartily appreciates the reviewers' constructive comments, which greatly enhance the quality of the paper to meet the journal standard.

Conflict of interest

The author declares no conflicts of interest in this paper.

References

[1]	E. Ackerman, J.W. Rosevear andW. F. McGuckin, A mathematical model of the glucose-tolerance test, Phys. Med. Biol., 9 (1964), 203–213.
[2]	E. Ackerman, L. C. Gatewood, J. W. Rosevear and G. D. Molnar, Model studies of blood glucose regulation, Bull. Math. Biophys., 27 (1965), 21–37.
[3]	B. Ashley and W. Liu, Asymptotic tracking and disturbance rejection of the blood glucose regulation system, Math. Biosci., 289 (2017), 78–88.
[4]	R. N. Bergman, Y. Z. Ider, C. R. Bowden and C. Cobelli, Quantitative estimation of insulin sensitivity, Am. J. Physiol. Endocrinol. Metab., 236 (1979), E667–E677.
[5]	R. N. Bergman, L. S. Phillips and C. Cobelli, Measurement of insulin sensitivity and β-cell glucose sensitivity from the response to intraveous glucose, J. Clin. Invest., 68 (1981), 1456–1467.
[6]	R. N. Bergman, D. T. Finegood and M. Ader, Assessment of insulin sensitivity in vivo, Endocrine Reviews, 6 (1985), 45–86.
[7]	R. N. Bergman, Toward physiological understanding of glucose tolerance, Minimal-model approach, Diabetes, 38 (1989), 1512–1527.
[8]	A. Bertoldo, R. R. Pencek, K. Azuma, J. C. Price, C. Kelley, C. Cobelli and D. E. Kelley, Interactions between delivery, transport, and phosphorylation of glucose in governing uptake into human skeletal muscle, Diabetes, 55 (2006), 3028–3037.
[9]	J. Carr, Applications of Center Manifold Theory, Applied Mathematical Sciences 35, Springer, New York, 1981.
[10]	C. Cobelli, G. Federspil, G. Pacini, A. Salvan and C. Scandellari, An integrated mathematical model of the dynamics of blood glucose and its hormonal control, Math. Biosci., 58 (1982), 27– 60.
[11]	K. Fessel, J. B. Gaither, J. K. Bower, G. Gaillard and K. Osei, Mathematical analysis of a model for glucose regulation, Mathematical Biociences and Engineering, 13(2016), 83–90.
[12]	L. B. Freidovich and H. K. Khalil, Performance recovery of feedback-linearization-based designs, IEEE Trans. Automat. Contr., 53 (2008), 2324–2334.
[13]	W. T. Garvey, L. Maianu, J. H. Zhu, G. Brechtel-Hook, P. Wallace and A. D. Baron, Evidence for defects in the tra cking and translocation of GLUT4 glucose transporters in skeltal muscle as a cause of human insulin resistance, J. Clin. Invest., 101 (1998), 2377–2386.
[14]	C. J. Goodner, B. C. Walike, D. J. Koerker, J. W. Ensinck, A. C. Brown, E. W. Chideckel, J. Palmer and L. Kalnasy, Insulin, glucagon, and glucose exhibit synchronous, sustained oscillations in fasting monkeys, Science, 195 (1977), 177–179.
[15]	O. I. Hagren and A. Tengholm, Glucose and insulin synergistically activate phosphatidylinositol 3- kinase to trigger oscillations of phosphatidylinositol 3,4,5-trisphosphate in β-cells, J. Biol. Chem., 281 (2006), 39121–39127.
[16]	B. C. Hansen, K. C. Jen, S. B. Pek and R. A.Wolfe, Rapid oscillations in plasma insulin, glucagon, and glucose in obese and normal weight humans. J. Clin. Endocr. Metab., 54 (1982), 785–792.
[17]	A. Klip and M. Vranic, Muscle, liver, and pancreas: Three Musketeers fighting to control glycemia, Am. J. Physiol. Endocrinol. Metab., 291 (2006), E1141–E1143.
[18]	R. Hovorka, Continuous glucose monitoring and closed-loop systems, Diabetic Med., 23 (2006), 1–12.
[19]	J. Huang, Nonlinear output regulation, theory and applications, Society for Industrial and Applied Mathematics, Philadelphia, 2004.
[20]	H. Kang, K. Han and M. Choi, Mathematical model for glucose regulation in the whole-body system, Islets, 4 (2012), 84–93.
[21]	D. A. Lang, D. R. Matthews, J. Peto, and R. C. Turner, Cyclic oscillations of basal plasma glucose and insulin concentrations in human beings, New Engl. J. Med., 301 (1979), 1023–1027.
[22]	J. Lee, R. Mukherjee and H.K. Khalil, Output feedback performance recovery in the presence of uncertainties, Syst. Control Lett., 90 (2016), 31–37.
[23]	J. Li, Y. Kuang and C. C. Mason, Modeling the glucose-insulin regulatory system and ultradian insulin secretory oscillations with two explicit time delays, J. Theor. Biol., 242 (2006), 722–735.
[24]	W. Liu and F. Tang, Modeling a simplified regulatory system of blood glucose at molecular levels, J. Theor. Biol., 252 (2008), 608–620.
[25]	W. Liu, C. Hsin, and F. Tang, A molecular mathematical model of glucose mobilization and uptake, Math. Biosciences, 221 (2009), 121–129.
[26]	W. Liu, Elementary Feedback Stabilization of the Linear Reaction Diffusion Equation and the Wave Equation, Mathematiques et Applications, Vol. 66, Springer, 2010.
[27]	W. Liu, Introduction to Modeling Biological Cellular Control Systems, Modeling, Simulation and Applications, Vol. 6, Springer, 2012.
[28]	K. Ma, H. K. Khalil and Y. Yao, Guidance law implementation with performance recovery using an extended high-gain observer, Aerosp. Sci. Technol., 24 (2013), 177–186.
[29]	C. D. Man, A. Caumo, R. Basu, R. A. Rizza, G. Toffolo and C. Cobelli, Minimal model estimation of glucose absorption and insulin sensitivity from oral test: validation with a tracer method, Am. J. Physiol. Endocrinol. Metab., 287 (2004), E637–E643.
[30]	C. D. Man, M. Campioni, K. S. Polonsky, R. Basu, R. A. Rizza, G. Toffolo and C. Cobelli, Twohour seven-sample oral glucose tolerance test and meal protocol: minimal model assessment of β-cell responsivity and insulin sensitivity in nondiabetic individuals, Diabetes, 54 (2005), 3265– 3273.
[31]	C. D. Man, R. A. Rizza and C. Cobelli, Meal simulation model of the glucose-insulin system, IEEE Trans. Biomed. Eng., 54 (2007), 1740–1749.
[32]	H. Nishimura, F. Pallardo, G. A. Seidner, S. Vannucci, I. A. Simpson and M. J. Birnbaum, Kinetics of GLUT1 and GLUT4 glucos transporters expressed in Xenopus oocytes, J. Biol. Chem., 268 (1993), 8514–8520.
[33]	A. E. Panteleon, M. Loutseiko, G. M. Steil and K. Rebrin, Evaluation of the effect of gain on the meal response of an automated closed-loop insulin delivery system, Diabetes, 55 (2006), 1995– 2000.
[34]	A. R. Sedaghat, A. Sherman and M. J. Quon, A mathematical model of metabolic insulin signaling pathways, Am. J. Physiol. Endocrinol. Metab., 283 (2002), E1084–E1101.
[35]	E. T. Shapiro, H. Tillil, K. S. Polonsky, V. S. Fang, A. H. Rubenstein and E. V. Cauter, Oscillations in insulin secretion during constant glucose infusion in normal man: relationship to changes in plasma glucose, J. Clin. Endocr. Metab., 67 (1988), 307–314.
[36]	C. Simon, G. Brandenberger and M. Follenius, Ultradian oscillations of plasma glucose, insulin, and C-peptide in man during continuous enteral nutrition, J. Clin. Endocr. Metab., 64 (1987), 669–674.
[37]	J. T. Sorensen, A Physiological Model of Glucose Metabolism in Man and its Use to Design and Assess Improved Insulin Therapies for Diabetes, PhD Thesis, Massachusetts Institute of Technology, 1985.
[38]	G. M. Steil, K. Rebrin, C. Darwin, F. Hariri and M. F. Saad, Feasibility of automating insulin delivery for the treatment of type 1 diabetes, Diabetes, 55 (2006), 3344–3350.
[39]	J. Sturis, E. V. Cauter, J. D. Blackman and K. S. Polonsky, Entrainment of pulsatile insulin secretion by oscillatory glucose infusion, J. Clin. Invest., 87 (1991), 439-445.
[40]	J. Sturis, K. S. Polonsky, E. Mosekilde and E. V. Cauter, Computer model for mechanisms underlying ultradian oscillations of insulin and glucose, Am. J. Physiol. Endocrinol. Metab., 260 (1991), E801–E809.
[41]	G. Toffolo and C. Cobelli, The hot IVGTT two-compartment minimal model: an improved version, Am. J. Physiol. Endocrinol. Metab., 284 (2003), E317–E321.
[42]	G. Toffolo, M. Campioni, R. Basu, R. A. Rizza and C. Cobelli, A minimal model of insulin secretion and kinetics to assess hepatic insulin extraction, Am. J. Physiol. Endocrinol. Metab. 290 (2006), E169–E176.
[43]	I. M. Tolic, E. Mosekilde and J. Sturis, Modeling the insulin-glucose feedback system: the significance of pulsatile insulin secretion, J. Theor. Biol., 207 (2000), 361–375.
[44]	R. C. Turner, R. R. Holman, D. Matthews, T. D. Hockaday and J. Peto, Insulin deficiency and insulin resistance interaction in diabetes: estimation of their relative contribution by feedback analysis from basal plasma insulin and glucose concentrations, Metabolism, 28 (1979), 1086– 1096.
[45]	O. Vahidi, K. E. Kwok, R. B. Gopaluni and L. Sun, Developing a physiological model for type II diabetes mellitus, Biochem. Eng. J., 55 (2011), 7-16.
[46]	O. Vahidi1, K. E. Kwok, R. B. Gopaluni and F. K. Knop, A comprehensive compartmental model of blood glucose regulation for healthy and type 2 diabetic subjects, Med. Biol. Eng. Comput., 54 (2016), 1383-1398.
[47]	R. R. Wolfe, J. R. Allsop and J. F. Burke, Glucose metabolism in man: Responses to intravenous glucose infusion, Metabolism, 28 (1979), 210–220.

This article has been cited by:

1.	Weijiu Liu, The exponential tracking and disturbance rejection for the nonlinear wave equation via boundary control without any geometrical conditions on the domain, 2022, 508, 0022247X, 125863, 10.1016/j.jmaa.2021.125863
2.	Weijie Wang, Shaoping Wang, Yixuan Geng, Yajing Qiao, Teresa Wu, An OGI model for personalized estimation of glucose and insulin concentration in plasma, 2021, 18, 1551-0018, 8499, 10.3934/mbe.2021420

Reader Comments

Your name:*

Email:*
© 2019 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)

通讯作者: 陈斌, bchen63@163.com

1.
沈阳化工大学材料科学与工程学院沈阳 110142

Mathematical Biosciences and Engineering

4.4

Metrics

Article views(4912) PDF downloads(968) Cited by(2)

Preview PDF

Download XML

Export Citation

Article outline

Show full outline

Figures and Tables

Figures(6) / Tables(1)

Mathematical Biosciences and Engineering

A mathematical model for the robust blood glucose tracking

Related Papers:

Abstract

1. Introduction

2. Formulation of a model with parameter uncertainties

3. A static feedback and feedforward controller

4. A robust dynamical feedback and feedforward controller

5. Discussion and conclusion

Acknowledgments

Conflict of interest

References

This article has been cited by:

Reader Comments

通讯作者: 陈斌, bchen63@163.com

Metrics

Figures and Tables

Other Articles By Authors

Catalog

Mathematical Biosciences and Engineering

A mathematical model for the robust blood glucose tracking

Related Papers:

Abstract

1. Introduction

2. Formulation of a model with parameter uncertainties

3. A static feedback and feedforward controller

4. A robust dynamical feedback and feedforward controller

5. Discussion and conclusion

Acknowledgments

Conflict of interest

References

This article has been cited by:

Reader Comments

通讯作者: 陈斌, bchen63@163.com

Metrics

Figures and Tables

Other Articles By Authors

Related pages

Tools

Export File

Citation

Format

Content

Catalog