Fuzzy functional observer for the control of the glucose-insulin system

Abstract

This paper presents an algorithm for the design of a Generalized Dynamic Functional Observer (GDFO) for Takagi-Sugeno systems. A set of conditions are provided for the existence of the observer, which can be applied to the proportional functional observer that is presented as a particular case of the proposed generalized observer. Lyapunov function is used to find stability condition of the observer. Stability conditions are given in terms of linear matrix inequalities (LMIs). A glucose-insulin system example illustrates the performance and effectiveness of the proposed approach.

Keywords

Generalized dynamic functional observer Takagi-Sugeno systems glucose-insulin system feedback control

1. Introduction

Fuzzy theory represents a solution to control complex systems in their whole operation range. Fuzzy theory has been gradually constituted as a powerful technique of control [1 –4]. Fuzzy sets and fuzzy logic are a useful mathematical tools for modelling and controlling systems in industry, humanity, and nature (see [5 –8]).

The Takagi-Sugeno (T-S) system is a class of fuzzy systems, it has been widely studied over the past few years because it can effectively represent the system dynamics of nonlinear systems. The T-S systems are based on fuzzy rules IF-THEN, which represent linear local relationships of input and output of a nonlinear system. The main characteristic of a T-S model is to express the local dynamics of each fuzzy rule by a linear model, valid locally in a corresponding region. The global diffuse model is achieved through a combination of linear local models. Over decades, stability analysis, optimization, observer construction, and controller design for Takagi-Sugeno fuzzy systems have been proposed as in [9 –12].

Observers are mathematical algorithms that have become a useful tool because of their ability to estimate the state vector from measured inputs and outputs. Sometimes, state estimation is important for a wide range of applications in control, such as the design of a controller for output feedback and the detection of system failures. Often, in control applications, only a linear combination of state variables is required.

The idea of an observer to directly estimate a linear function of the state vector was initially proposed by [13], this observer was named as Functional observer. From this research, some works have been developed, as [14] where the authors present the necessary and sufficient conditions for the existence and stability of functional observers. Other works that address the design of functional observers reported in the literature are [15] among others. This type of observer has two main objectives, to estimate an specific set or state variables, or to estimate a control law that establishes the system from the input data and the measured outputs. Another important characteristic of the functional observer that has been widely studied is the condition of observability. In [16] they introduced the conditions of functional observability, and later, [17] simplified those conditions.

On the other hand, a new form of observer called generalized dynamic observer (GDO) which structure is formed by a dynamical part and a static part, it has been introduced for descriptor systems [18], linear time invariant systems [19], discrete-time systems [20] and continuous-time linear parameter variant (LPV) systems [21]. This observer structure is based on [22, 23], where the principal idea is to add dynamics to increase its degrees of freedom, with the purpose of achieving steady state accuracy and improve robustness in estimation error against disturbances and parametric uncertainties.

In this work, a Generalized Dynamic Functional Observer (GDFO) for Takagi-Sugeno type systems is proposed. The main contribution of this work is the structure of the GDFO, this observer presents an alternative structure that can be configured in a proportional functional observer (PFO) which is only a particular case. The main advantage of employing a functional observer lies in its ability to directly estimate a given linear function of the state vector without having to estimate all individual states, whereas a state observer design scheme can not (see Fig. 1).

Fig. 1

Comparison of the classic regulation approach and the functional approach.

The case of study presented in this paper, is the glucose-insulin system, this system requires a control in the supply of insulin, an equivocal supply of this component in patients can take them to hyperglycemic or hypoglycemic levels. The portable device that supplies insulin is of limited dimensions [24], so information storage is also restricted, in this sense, the control design based on a functional observer presents an alternative of regulation of insulin supply, since it can generate directly an appropriate control law from the input data and the measured outputs, that regulates glucose levels in face of disturbances associated with food intake. On the other hand, with the classical approach of control based observer, it is first necessary to obtain the estimation of all the system state variables, and then with these estimates generate the control law, this task requires greater computational demand and therefore storage.

Notation: The symbol (*) denotes the transpose element on the symmetric position. I_n is an identity matrix of dimension n × n, ones_n,m means a matrix of dimension n × m with all its elements one, A^⊥ denotes a maximal row rank matrix such that A^⊥A = 0, by convention A^⊥ = 0 when A is of full row rank. The symbol A⁺ denotes any generalized inverse of the matrix A, satisfying AA⁺A = A. The symbol $\dot{x} (t)$ represents the first derivative of x (t) with respect to time, i.e. $\frac{dx (t)}{dt}$ and $\hat{x} (t)$ represents the estimate of x (t).

2. Preliminaries

The proposed observer is developed for continuous-time Takagi-Sugeno systems, which can be expressed in the following form $\begin{matrix} \dot{x} (t) = \sum_{i = 1}^{m} ω_{i} (ρ) (A_{i} x (t)) + Bu (t) + Δ \\ y (t) = Cx (t) \\ z (t) = Lx (t) \end{matrix}$ (1) where $x (t) \in ℝ^{n}$ is the state vector of the system, $u (t) \in ℝ^{s}$ is the input vector, $y (t) \in ℝ^{p}$ is the measured output vector and $z (t) \in ℝ^{q}$ is the linear function that is required to estimate. A_i, B, Δ, C and L are known constant matrices of appropriate dimensions, m is the number of local models, ω_i are convex weighting functions depending on the premise variable vector ρ (t). The convex weighting functions satisfy the convex set sum property $\sum_{i = 1}^{m} ω_{i} (ρ) = 1, 0 \leq ω_{i} (ρ) \leq 1, \forall i = 1, \dots, m$ (2)

The following definition introduces the concept of functional observability.

Definition 1. [16] Functional observability is a test of whether linear functions of states of a system can be inferred by knowledge of external outputs and inputs of the system.

Lemma 1. There exists a linear functional observer with arbitrary poles such that $\hat{z} (t)$ converges asymptotically to z (t) if the following two conditions are satisfied:

The matrix L must satisfy $rank [\begin{matrix} {LA}_{i} \\ {CA}_{i} \\ C \\ L \end{matrix}] = rank [\begin{matrix} {CA}_{i} \\ C \\ L \end{matrix}], \forall i = 1, \dots, m$

The matrix L must satisfy, for all $s \in ℂ$

$rank [\begin{matrix} sL - {LA}_{i} \\ {CA}_{i} \\ C \end{matrix}] = rank [\begin{matrix} {CA}_{i} \\ C \\ L \end{matrix}]$

Let us recall the following definition and lemma:

Definition 2. [25] The pairs (A_i, B) ∀i = 1, …, m are controllable if $rank [\begin{matrix} B & A_{i} B & \dots & A_{i}^{n - 1} B \end{matrix}] = n$ where n is the dimension of A_i, then the system (1) is controllable.

Lemma 2. [26] Let matrices $B$ and $Q_{i}$ be given. Then the following statements are equivalent.

There exists a matrix $X_{i}$ satisfying $B X_{i} + {(B X_{i})}^{T} + Q_{i} < 0$

The following condition holds: $B^{⊥} Q_{i} B^{⊥ T} < 0 or B B^{T} > 0$

Suppose that the above statements hold and further assume that

B^{T} B > 0

. Then all matrices

X_{i}

in statement i are given by

X_{i} = - σ B^{T} + \sqrt{σ} L ϑ_{i}^{1 / 2}

where

L

is any matrix such that

| | L | | < 1

y σ > 0 is any scalar such that

ϑ_{i} = σ B B^{T} - Q_{i} > 0

Assumption 1. System (1) is controllable (see Definition 2).

Assumption 2. The conditions of Lemma 1 are verified, this implies that it is possible to estimate the function z (t) with the inputs and outputs measurements.

3. Problem statement

Let us consider the following generalized dynamic functional observer (GDFO) for system (1): $\begin{matrix} \dot{ζ} (t) & = \sum_{i = 1}^{m} ω_{i} (ρ) (N_{i} ζ (t) + H_{i} ν (t) + F_{i} y (t)) + \\ Ju (t) + T Δ \end{matrix}$ (3) $\dot{ν} (t) = \sum_{i = 1}^{m} ω_{i} (ρ) (S_{i} ζ (t) + G_{i} ν (t) + M_{i} y (t))$ (4) $\hat{z} (t) = P ζ (t) + Qy (t)$ (5) where $ζ (t) \in ℝ^{q_{0}}$ represents the state vector of the observer, $ν (t) \in ℝ^{q_{1}}$ is an auxiliary vector and $\hat{z} (t) \in ℝ^{q}$ is the estimate of z (t). N_i, H_i, F_i, J, S_i, G_i, M_i, P and Q are constant matrices of appropriate dimensions to be determined such that $\hat{z} (t)$ asymptotically converges to z (t).

Lemma 3. There exists an observer of the form (3)-(5) for system (1) if matrices $\begin{matrix} [\begin{matrix} N_{i} & H_{i} \\ S_{i} & G_{i} \end{matrix}], \forall i = 1, \dots, m \end{matrix}$ are Hurwitz and if there exists a matrix T such that the following conditions are satified:

N_iT + F_iC - TA_i = 0

J = TB

S_iT + M_iC = 0

PT + QC = L

Proof. Let $T \in ℝ^{q_{0} \times n}$ be a parameter matrix and define the transformed error vector ɛ (t) = ζ (t) - Tx (t), such that its derivative is given by $\dot{ɛ} (t) = \sum_{i = 1}^{m} ω_{i} (ρ) (N_{i} ɛ (t) + H_{i} ν (t) +$ (6) $(N_{i} T + F_{i} C - {TA}_{i}) x (t)) + (J - TB) u (t)$

By using the definition of ɛ (t) in (4) and (5), they can be rewritten as $\begin{matrix} \dot{ν} (t) & = \sum_{i = 1}^{m} ω_{i} (ρ) (S_{i} ɛ (t) + (S_{i} T + M_{i} C) x (t) + \\ L_{i} ν (t)) \end{matrix}$ (7) $\hat{z} (t) = P ɛ (t) + (PT + QC) x (t)$ (8) Considering that conditions a)–c) are satisfied, then the dynamics error, formed by (6) and (7), can be written as: $\underset{\dot{φ} (t)}{\underset{︸}{[\begin{matrix} \dot{ɛ} (t) \\ \dot{ν} (t) \end{matrix}]}} = \sum_{i = 1}^{m} ω_{i} (ρ) \underset{A_{i}}{\underset{︸}{[\begin{matrix} N_{i} & H_{i} \\ S_{i} & G_{i} \end{matrix}]}} \underset{φ}{\underset{︸}{[\begin{matrix} ɛ (t) \\ ν (t) \end{matrix}]}}$ (9) and using condition d) we have the following expression from (8) $\hat{z} (t) - z (t) = e_{z} (t) = P ɛ (t)$ (10) in this case if matrices $A_{i}$ are Hurwitz ∀i = 1, …, m we have that the $lim_{t \to \infty} φ (t) = 0$ , then $lim_{t \to \infty} e_{z} (t) = 0$ .□

4. Generalized dynamic functional observer design

4.1. Parameterization of the observer matrices

Consider a matrix $E \in ℝ^{q_{0} \times n}$ of full row rank with arbitrary elements, such that $Σ = [\begin{matrix} E \\ C \end{matrix}]$ be a full column rank matrix.

Conditions c) and d) can be written as: $[\begin{matrix} S_{i} & M_{i} \\ P & Q \end{matrix}] [\begin{matrix} T \\ C \end{matrix}] = [\begin{matrix} 0 \\ L \end{matrix}]$ (11) The necessary and sufficient condition for (11) to have a solution is $rank [\begin{matrix} T \\ C \end{matrix}] = rank [\begin{matrix} T \\ C \\ 0 \\ L \end{matrix}] = n$ . Once this condition is satisfied, there exist two matrices

$T \in ℝ^{q_{0} \times n}$ and $K \in ℝ^{q_{0} \times p}$ such that $T + KC = E$ (12) which can be written as $[\begin{matrix} T & K \end{matrix}] Ω = E$ (13) where $Ω = [\begin{matrix} I_{n} \\ C \end{matrix}]$ , since the $rank (Ω) = rank [\begin{matrix} Ω \\ E \end{matrix}]$ the general solution of (13) is given by $[\begin{matrix} T & K \end{matrix}] = E Ω^{+} - Z (I_{n + p} - Ω Ω^{+})$ (14) which can be decomposed in $T = T_{1} - {ZT}_{2}$ (15) $K = K_{1} - {ZK}_{2}$ (16) where $T_{1} = E Ω^{+} [\begin{matrix} I_{n} \\ 0 \end{matrix}],$

$T_{2} = (I_{n + p} - Ω Ω^{+}) [\begin{matrix} I_{n} \\ 0 \end{matrix}]$ , $K_{1} = E Ω^{+} [\begin{matrix} 0 \\ I_{p} \end{matrix}],$ $K_{2} = (I_{n + p} - Ω Ω^{+}) [\begin{matrix} 0 \\ I_{p} \end{matrix}]$ and Z is a matrix of appropriate dimension with arbitrary elements.

Considering matrix T from equation (12), condition a) of Lemma 3 can be written as: $N_{i} E + {\tilde{K}}_{i} C = {TA}_{i}$ (17) where ${\tilde{K}}_{i} = F_{i} - N_{i} K$ . Equation (17) can be written as $[\begin{matrix} N_{i} & {\tilde{K}}_{i} \end{matrix}] Σ = {TA}_{i}$ (18) The general solution of (18), considering the parameterization of matrix T is given by $N_{i} = N_{1, i} - {ZN}_{2, i} - Y_{1, i} N_{3}$ (19) ${\tilde{K}}_{i} = {\tilde{K}}_{1, i} - Z {\tilde{K}}_{2, i} - Y_{1, i} {\tilde{K}}_{3}$ (20) where $N_{1, i} = T_{1} A_{i} Σ^{+} [\begin{matrix} I_{q_{0}} \\ 0 \end{matrix}],$

$N_{2, i} = T_{2} A_{i} Σ^{+} [\begin{matrix} I_{q_{0}} \\ 0 \end{matrix}],$

$N_{3} = (I_{q_{0} + p} - Σ Σ^{+}) [\begin{matrix} I_{q_{0}} \\ 0 \end{matrix}]$ , ${\tilde{K}}_{1, i} = T_{1} A_{i} Σ^{+} [\begin{matrix} 0 \\ I_{p} \end{matrix}],$ ${\tilde{K}}_{2, i} = T_{2} A_{i} Σ^{+} [\begin{matrix} 0 \\ I_{p} \end{matrix}]$ and ${\tilde{K}}_{3} = (I_{q_{0} + p} - Σ Σ^{+}) [\begin{matrix} 0 \\ I_{p} \end{matrix}]$ and Y_1,i are matrices of appropriate dimensions with arbitrary elements.

From equation (17) we can determine that

$F_{i} = {\tilde{K}}_{i} + N_{i} K$ , which can be written as $\begin{matrix} F_{i} & = {\tilde{K}}_{1, i} + N_{1, i} K - Z ({\tilde{K}}_{2} + N_{2, i} K) - \\ Y_{1, i} ({\tilde{K}}_{3} + N_{3} K) \end{matrix}$ (21) $F_{i} = F_{1, i} - {ZF}_{2, i} - Y_{1, i} F_{3}$ (22) where $F_{1, i} = T_{1} A_{i} Σ^{+} [\begin{matrix} K \\ I_{p} \end{matrix}]$ , $F_{2} = (T_{2} A_{i} Σ^{+}) [\begin{matrix} K \\ I_{p} \end{matrix}]$ and $F_{3} = (I_{q_{0} + p} - Σ Σ^{+}) [\begin{matrix} K \\ I_{p} \end{matrix}]$ .

On the other hand, by using (12) we can obtain the following expression $[\begin{matrix} T \\ C \end{matrix}] = [\begin{matrix} I_{q_{0}} & - K \\ 0 & I_{p} \end{matrix}] Σ$ (23) replacing (23) in (11) we get $[\begin{matrix} S_{i} & M_{i} \\ P & Q \end{matrix}] [\begin{matrix} I_{q_{0}} & - K \\ 0 & I_{p} \end{matrix}] Σ = [\begin{matrix} 0 \\ L \end{matrix}]$ (24) which leads to the following solution

$\begin{matrix} [\begin{matrix} S_{i} & M_{i} \\ P & Q \end{matrix}] & = ([\begin{matrix} 0 \\ L \end{matrix}] Σ^{+} - [\begin{matrix} Y_{2, i} \\ Y_{3} \end{matrix}] (I_{q_{0} + p} - Σ Σ^{+})) \times \\ [\begin{matrix} I_{q_{0}} & K \\ 0 & I_{p} \end{matrix}] \end{matrix}$ (25) where we have used the inverse as ${[\begin{matrix} I_{q_{0}} & - K \\ 0 & I_{p} \end{matrix}]}^{- 1} = [\begin{matrix} I_{q_{0}} & K \\ 0 & I_{p} \end{matrix}]$ , with matrices Y_2,i and Y₃ as matrices of appropriate dimensions with arbitrary elements.

Then, matrices S_i, M_i, P and Q can be written as $S_{i} = - Y_{2, i} N_{3}$ (26) $M_{i} = - Y_{2, i} F_{3}$ (27) $P = P_{1} - Y_{3} N_{3}$ (28) $Q = Q_{1} - Y_{3} F_{3}$ (29) where $P_{1} = L Σ^{+} [\begin{matrix} I_{q_{0}} \\ 0 \end{matrix}]$ and $Q_{1} = L Σ^{+} [\begin{matrix} K \\ I_{p} \end{matrix}]$ .

By using (19) and (26) the observer error dynamics (9) can be written as $\dot{φ} (t) = \sum_{i = 1}^{m} ω_{i} (ρ) ((A_{i} - Y_{i} A_{2}) φ (t))$ (30) where $φ (t) = [\begin{matrix} ɛ (t) \\ ν (t) \end{matrix}]$ , $A_{i} = [\begin{matrix} N_{1, i} - {ZN}_{2} & 0 \\ 0 & 0 \end{matrix}]$ , $A_{2} = [\begin{matrix} N_{3} & 0 \\ 0 & - I_{q_{0}} \end{matrix}]$ and $Y_{i} = [\begin{matrix} Y_{1, i} & H_{i} \\ Y_{2, i} & G_{i} \end{matrix}]$ .

The problem has been now reduced to find matrices $Y_{i}$ and Z such that (30) is stable. This can be reached by using the linear matrix inequality (LMI) approach.

4.2. Observer stability analysis

In this section a method to design the observer is presented. The following theorem gives the LMI conditions that allow the obtention of all the observer matrices.

Theorem 1. There exists parameter matrices $Y_{i}$ and Z such that system (30) is asymptotically stable if there exists a matrix $X = [\begin{matrix} X_{1} & X_{2} \\ X_{2}^{T} & X_{3} \end{matrix}] > 0$ such that the following LMI is verified

$\begin{matrix} N_{3}^{T ⊥} & [N_{1, i}^{T} X_{1} - N_{2, i}^{T} W^{T} + \\ X_{1} N_{1, i} - {WN}_{2, i}] N_{3}^{T ⊥ T} < 0 \end{matrix}$ (31) where $Z = X_{1}^{- 1} W$ , by using Lemma 2 matrix $Y_{i}$ is parameterized as $Y_{i} = X^{- 1} {(- σ B^{T} + \sqrt{σ} L ϑ_{i}^{1 / 2})}^{T}$ (32) where $L$ is a matrix of arbitrary elements such that $| | L | | < 1$ and σ > 0 is any scalar that satisfy: $ϑ_{i} = σ B B^{T} - Q_{i} > 0$ (33) with $Q_{i} =$

$[\begin{matrix} {(N_{1, i} - {ZN}_{2, i})}^{T} X_{1} + X_{1} (N_{1, i} - {ZN}_{2, i}) & (*) \\ X_{2}^{T} (N_{1, i} - {ZN}_{2, i}) & 0 \end{matrix}]$ and $B = [\begin{matrix} - N_{3}^{T} & 0 \\ 0 & I_{q_{0}} \end{matrix}]$ .

Proof. Consider the following Lyapunov function $V (φ (t)) = φ (t)^{T} X φ (t) > 0$ (34) with $X = X^{T} = [\begin{matrix} X_{1} & X_{2} \\ X_{2}^{T} & X_{3} \end{matrix}] > 0$ . Its derivative along the trajectory of (30) is given by

$\begin{matrix} \dot{V} (φ (t)) = & \sum_{i = 1}^{m} ω_{i} (ρ) φ (t)^{T} [{(A_{i} - Y_{i} A_{2})}^{T} X + \\ (A_{i} - Y_{i} A_{2})] φ (t) < 0 \end{matrix}$ (35) The asymptotic stability of system (30) is guaranteed if and only if $\dot{V} (φ (t)) < 0$ . This leads to the following LMI $\sum_{i = 1}^{m} ω_{i} (ρ) [{(A_{i} - Y_{i} A_{2})}^{T} X + X (A_{i} - Y_{i} A_{2})] < 0$ (36) which can be written as $\sum_{i = 1}^{m} ω_{i} (ρ) (B X_{i} + {(B X_{i})}^{T} + Q_{i}) < 0$ (37) where $B = - A_{2}^{T}$ , $Q_{i} = A_{i}^{T} X + X A_{i}$ and $X_{i} = Y_{i}^{T} X$ .□

According to Lemma 2 there exist matrices $X_{i}$ satisfying (37) if and only if the following condition holds: $B^{⊥} Q_{i} B^{⊥ T} < 0$ (38) with $B = [\begin{matrix} - N_{3}^{T ⊥} & 0 \end{matrix}]$ and $Q_{i}$ ∀i ∈ 1, …, m defined in (46). Considering matrices $B$ and $Q_{i}$ we obtain inequality (31). If (38) BO or its equivalence (31) is verified, then matrices $Y_{i}$ are obtained as shown in (32).

The following algorithm summarize the observer design to obtain the corresponding matrices.

Algorithm 1:

Choose a matrix $E \in ℝ^{q_{0} \times n}$ such that matrix Σ is of full column rank.

Compute the matrices N_1,i, N_2,i, N₃, T₁, T₂, K₁,K₂, P₁ and Q₁ as shown in section 4.1.

Solve the LMI of the equation (31) to find matrices X and Z.

Choose a matrix $L$ , such that $| | L | | < 1$ and a scalar σ > 1 such that ϑ_i > 0, then determine the parameter matrix $Y_{i}$ as in the equation (32).

Compute all the matrices of the GDFO (3)-(5), using (19) to compute N_i, (32) to compute H_i and G_i, (26)-(29) to compute S_i, M_i, P y Q taking matrix Y₃ = 0. F_i is given in (22) and matrix J is given by the condition b).

5. Particular case

5.1. Proportional functional observer

In order to obtain a Proportional Functional Observer (PFO) from the structure of the GDFO, the following considerations must be taken into account: H_i = 0, S_i = 0, G_i = 0 and M_i = 0, to obtain the following observer: $\begin{matrix} \dot{ζ} (t) & = \sum_{i = 1}^{m} ω_{i} (ρ) (N_{i} ζ (t) + F_{i} y (t)) + Ju (t) + T Δ \\ \hat{z} (t) & = P ζ (t) + Qy (t) \end{matrix}$ (39) and the observer error dynamics (30) becomes: $\dot{ɛ} (t) = ({\tilde{A}}_{i} - {\tilde{Y}}_{i} {\tilde{A}}_{2}) ɛ (t)$ (40) where ${\tilde{A}}_{i} = [\begin{matrix} N_{1, i} - {ZN}_{2, i} \end{matrix}]$ , ${\tilde{A}}_{2} = - N_{3}$ and

${\tilde{Y}}_{i} = Y_{1, i}$ . Consequently, matrices $Q_{i}$ , $B$ and $X_{i}$ of Theorem 1. become $\begin{matrix} Q_{i} = & {(N_{1, i} - {ZN}_{2, i})}^{T} X + X (N_{1, i} - {ZN}_{2, i}) \\ B = - N_{3} and X_{i} = Y_{i}^{T} X \end{matrix}$ With these matrices the observer design can be obtained directly by applying Algorithm 1 .

6. Application to a study case

6.1. Model of the glucose-insulin system

Diabetes mellitus is a group of metabolic disorders, whose main characteristic is the presence of high concentrations of glucose in the blood persistently or chronically, due either to a defect in the production of insulin, to a resistance to the action of it to use the glucose, to an increase in the production of glucose, or to a combination of these causes, these factors determine the type of diabetes either type 1 or type 2.

The model of the glucose-insulin that is shown in [27] is a modification of the minimal model of Bergman in [28] that represents the system of glucose-insulin of a patient with diabetes mellitus type 1. This model is described by the following equations: $\begin{matrix} \frac{d G_{Δ} (t)}{dt} & = - p_{1} G_{Δ} (t) - \\ X (t) [G_{Δ} (t) + G_{b}] + D (t) \end{matrix}$ (41) $\frac{dX (t)}{dt} = - p_{2} X (t) + p_{3} I_{Δ} (t)$ (42) $\frac{{dI}_{Δ} (t)}{dt} = - N [I_{Δ} (t) + I_{b}] + \frac{U (t)}{V_{I}}$ (43) where U (t) is the infusion of exogenous insulin, G_Δ and I_Δ represent the differences in plasma glucose concentration and plasma insulin concentration from their basal values G_b and I_b, i.e. G (t) = G_b + G_Δ (t) and I (t) = I_b + I_Δ (t). D (t) is a disturbance, which represents the effect of glucose resulting from a meal, which is given in the following form [28] $D (t) = {Be}^{- drate \cdot t}$ (44) In the literature there is another way to represent the disturbance D (t), which can be calculated using the Rayleigh probability density function as follow: $D (t) = \frac{kt}{b^{2}} e^{- t^{2} / 2 b^{2}}$ (45) with different values for k depending on the carbohydrates amount, and for the constant b [29].

The parameters definition and values are given in Table 1. It should be noted that for diabetic patients, the basal value of the insulin concentration in plasma I_b is not a natural value, but it must be interpreted as the objective value for the insulin infusion program.

Table 1

Parameters definition

Parameter	Value	Definition
G _b	81[mg/dL]	Subject’s baseline glycemia
I _b	2.5[mU/dL]	Subject’s baseline insulinemia
V _I	120[dL]	Insulin distribution volume
p ₁	$0.05 [\min^{- 1}]$	Glucose mass action rate constant
p ₂	0.1 $[\min^{- 1}]$	Rate constant expressing the spontaneous decrease of tissue glucose uptake ability
p ₃	0.00065 $[\frac{dL}{mU} \min^{2}]$	Insulin-dependent rate of increase in tissue glucose uptake ability, per unit of insulin concentration excess over baseline insulin
N	5/54 $[\min^{- 1}]$	Fractional disappearance rate of insulin
drate	0.05 $[\frac{mg}{dL} \min]$	Decomposition rate of meal disturbance
B	15g, 19g, 18g	Amount of glucose ingested during the meal.
k	15 $[\frac{g}{dL}]$ , 18 $[\frac{g}{dL}]$ , 20 $[\frac{g}{dL}]$	Constant depending on the size of the meal.
b	80	Constant for the Rayleigh probability density function.

6.2. Takagi-Sugeno formulation

The nonlinear sector approach is one of the most used for the construction of Takagi-Sugeno models and for the design of fuzzy control, since an exact representation of a nonlinear system can be obtained in a compact set of state space [30].

By following the methodology to represent a nonlinear model in a Takagi-Sugeno system using the nonlinear sector approach, the system (41)-(43) can be expressed as: $\begin{matrix} \dot{x} (t) & = A (x (t)) x (t) + Bu (t) + Dw (t) + Δ \\ y (t) & = Cx (t) \end{matrix}$ (46) where $A = [\begin{matrix} - p_{1} & - (G_{Δ} (t) + G_{b}) & 0 \\ 0 & - p_{2} & p_{3} \\ 0 & 0 & - N \end{matrix}], B = [\begin{matrix} 0 \\ 0 \\ \frac{1}{V_{I}} \end{matrix}],$

$Δ = [\begin{matrix} 0 \\ 0 \\ - {NI}_{b} \end{matrix}], x (t) = [\begin{matrix} G_{Δ} (t) \\ X (t) \\ I_{Δ} (t) \end{matrix}],$ $D = [\begin{matrix} 1 \\ 0 \\ 0 \end{matrix}],$

$C = [\begin{matrix} 1 & 0 & 0 \end{matrix}], u (t) = U (t)$ and w (t) = D (t). In this case, the nonlinearity of the system is given by

$ρ (t) = - (G_{Δ} (t) + G_{b})$ , there is a single nonlinearity, so then m = 2^p local models are obtained, where p is the number of nonlinearities, is this case two local models are defined.

It is known that ρ (t) varies in a bounded region $ρ (t) = - (G_{Δ} (t) + G_{b})$ , where $\underline{ρ}$ and $\bar{ρ}$ are the lower and upper limit of variation of ρ (t).

In such way, that the membership functions are given by $ω_{1} (ρ (t)) = η_{0}^{1} = \frac{\bar{ρ} - ρ (t)}{\bar{ρ} - \underline{ρ}},$ (47) $ω_{2} (ρ (t)) = η_{1}^{1} = 1 - η_{0}^{1}$ (48) so that, the following Takagi-Sugeno representation is obtained

$\begin{matrix} \dot{x} (t) & = \sum_{i = 1}^{m} ω_{i} (ρ) (A_{i} x (t)) + Bu (t) + Dw (t) + Δ \\ y (t) & = Cx (t) \end{matrix}$ (49) where $A_{1} = [\begin{matrix} - p_{1} & \underline{ρ} & 0 \\ 0 & - p_{2} & p_{3} \\ 0 & 0 & - N \end{matrix}],$ $A_{2} = [\begin{matrix} - p_{1} & \bar{ρ} & 0 \\ 0 & - p_{2} & p_{3} \\ 0 & 0 & - N \end{matrix}],$ matrices B, D, Δ and C are defined in (46).

7. Regulation of the glucose-insulin system

This section presents the regulation of the glucose-insulin system by a functional observer to demonstrate the results obtained in this paper. The objective of the test is stabilizing the system using a functional observer, to compare the response of the GDFO and the PFO in presence of a disturbance due to food intake.

Considering the parameters of Table 1 the following matrices of the Takagi-Sugeno system (49) are obtained.

$A_{1} = [\begin{matrix} - 0.05 & - 220.6742 & 0 \\ 0 & - 0.1 & 0.0006 \\ 0 & 0 & - 0.0926 \end{matrix}]$ ,

$B = [\begin{matrix} 0 \\ 0 \\ 0.0083 \end{matrix}]$ ,

$A_{2} = [\begin{matrix} - 0.05 & - 81 & 0 \\ 0 & - 0.1 & 0.0006 \\ 0 & 0 & - 0.0926 \end{matrix}]$ , $D = [\begin{matrix} 1 \\ 0 \\ 0 \end{matrix}]$ , $Δ = [\begin{matrix} 0 \\ 0 \\ - 0.2315 \end{matrix}]$

and $C = [\begin{matrix} 1 & 0 & 0 \end{matrix}]$ .

The aim is to estimate a linear function that allows to stabilize the system, by using a functional observer. In this case, the matrix L of the function to be estimated, is determined by state feedback u (t) = - Kx (t) by assigning the poles of the system in a circular LMI region with center q = -0.11 and radio r = 0.13, to give us $L_{m} atrix L = - K = [\begin{matrix} - 0.1692 & 593.7892 & 8.7418 \end{matrix}] .$

By following Algorithm 1 the GDFO can be designed as follows:

1. Chose matrix $E = [\begin{matrix} 3 & 1 & 2 \\ 1 & 2 & 3 \\ 2 & 1 & 4 \end{matrix}]$ such that rank (Σ) =3 .

2. Obtain matrices

$N_{1, 1} = [\begin{matrix} - 59.1923 & - 241.1751 & 210.4313 \\ - 19.7136 & - 80.4889 & 70.1544 \\ - 39.4136 & - 160.7799 & 140.1993 \end{matrix}]$ , $T_{1} = [\begin{matrix} 1.5 & 1 & 2 \\ 0.5 & 2 & 3 \\ 1 & 1 & 4 \end{matrix}]$ ,

$N_{1, 2} = [\begin{matrix} - 21.7300 & - 88.5511 & 77.2322 \\ - 7.2262 & - 29.6142 & 25.7546 \\ - 14.4387 & - 59.0306 & 51.3999 \end{matrix}]$ , $T_{2} = [\begin{matrix} 0.5 & 0 & 0 \\ 0 & 0 & 0 \\ - 0.5 & 0 & 0 \end{matrix}],$

$N_{2, 1} = [\begin{matrix} - 19.7382 & - 80.3743 & 70.1499 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ 19.7382 & 80.3743 & - 70.1499 \end{matrix}]$ $K_{1} = [\begin{matrix} 1.5 \\ 0.5 \\ 1 \end{matrix}]$ ,

$N_{2, 2} = [\begin{matrix} - 7.2508 & - 29.4997 & 25.7502 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ 7.2508 & 29.4997 & - 25.7502 \end{matrix}]$ $K_{2} = [\begin{matrix} - 0.5 \\ 0 \\ 0 \\ 0.5 \end{matrix}]$ ,

$N_{3} = [\begin{matrix} 0.1656 & - 0.0662 & - 0.0331 \\ - 0.0662 & 0.0265 & 0.0132 \\ - 0.0331 & 0.0132 & 0.0066 \\ - 0.3642 & 0.1457 & 0.0728 \end{matrix}],$ $Q_{1} = [\begin{matrix} 7.5113 \end{matrix}]$

and $P_{1} = [\begin{matrix} - 104.2037 & - 431.6087 & 373.6230 \end{matrix}]$ .

3. By using the YALMIP toolbox, to solve the LMI (31) finding

$X = [\begin{matrix} 0.5629 & - 0.1955 & - 0.2561 & 0.1348 & - 0.0268 & - 0.0379 \\ - 0.1955 & 0.7449 & - 0.0125 & - 0.0268 & 0.1643 & - 0.0065 \\ - 0.2561 & - 0.0125 & 0.6154 & - 0.0379 & - 0.0065 & 0.1444 \\ 0.1348 & - 0.0268 & - 0.0379 & 0.1000 & 0 & 0 \\ - 0.0268 & 0.1643 & - 0.0065 & 0 & 0.1000 & 0 \\ - 0.0379 & - 0.0065 & 0.1444 & 0 & 0 & 0.1000 \end{matrix}],$

$Z = [\begin{matrix} 3.2081 & 2.1809 & 2.1809 & 1.6512 \\ 1.9593 & 1.2983 & 1.2983 & 0.9606 \\ 2.9068 & 1.7770 & 1.7770 & 1.2714 \end{matrix}]$ .

4. Consider matrix $L = {ones}_{7, 6} \times 0.1$ such that $| | L | | < 1$ , and choose σ = 1000 such that ϑ_i > 0 then from (32) we obtain

$Y_{1} = [\begin{matrix} 0.004 & - 0.067 & - 0.057 & - 0.157 & 0.372 & - 0.001 & 0.015 \\ - 0.041 & - 0.039 & - 0.040 & - 0.040 & 0.009 & 0.323 & - 0.036 \\ - 0.035 & - 0.050 & - 0.048 & - 0.070 & 0.020 & - 0.042 & 0.346 \\ 0.071 & 0.161 & 0.148 & 0.276 & - 1.391 & 0.172 & 0.202 \\ 0.167 & 0.143 & 0.147 & 0.114 & 0.187 & - 1.434 & 0.186 \\ 0.149 & 0.145 & 0.145 & 0.139 & 0.213 & 0.182 & - 1.397 \end{matrix}] \times 10^{4}$

$Y_{2} = [\begin{matrix} 0.003 & - 0.067 & - 0.057 & - 0.158 & 0.372 & - 0.001 & 0.015 \\ - 0.042 & - 0.040 & - 0.041 & - 0.038 & 0.008 & 0.322 & - 0.037 \\ - 0.034 & - 0.050 & - 0.047 & - 0.069 & 0.021 & - 0.042 & 0.347 \\ 0.072 & 0.161 & 0.149 & 0.276 & - 1.391 & 0.172 & 0.202 \\ 0.168 & 0.145 & 0.148 & 0.115 & 0.188 & - 1.432 & 0.187 \\ 0.148 & 0.143 & 0.144 & 0.138 & 0.211 & 0.181 & - 1.398 \end{matrix}] \times 10^{4}$

5. Finally, we compute all the matrices of the observer, taking for simplicity Y₃ = 0.

$N_{1} = [\begin{matrix} - 669.5788 & 140.4060 & 229.4387 \\ - 104.5084 & 41.5880 & 20.9941 \\ - 253.9907 & 69.4045 & 74.8495 \end{matrix}], F_{1} = [\begin{matrix} 441.1821 \\ 71.5546 \\ 169.5688 \end{matrix}],$

$N_{2} = [\begin{matrix} - 653.2811 & 214.5128 & 165.7099 \\ - 107.1697 & 42.7226 & 21.4738 \\ - 247.3392 & 87.1156 & 58.2404 \end{matrix}], F_{2} = [\begin{matrix} 440.9948 \\ 73.3868 \\ 167.7864 \end{matrix}],$

$H_{1} = [\begin{matrix} 3720.7685 & - 6.5388 & 150.5313 \\ 86.3528 & 3230.4601 & - 361.4838 \\ 200.7622 & - 422.5666 & 3465.0814 \end{matrix}], J = [\begin{matrix} 0.0167 \\ 0.0250 \\ 0.0333 \end{matrix}],$

$H_{2} = [\begin{matrix} 3714.9867 & - 12.3206 & 144.7495 \\ 77.3700 & 3221.4773 & - 370.4667 \\ 207.8029 & - 415.5259 & 3472.1221 \end{matrix}],$

$S_{1} = [\begin{matrix} 1041.9288 & - 416.7715 & - 208.3858 \\ 281.5928 & - 112.6371 & - 56.3186 \\ 404.6216 & - 161.8486 & - 80.9243 \end{matrix}], M_{1} = [\begin{matrix} - 713.5534 \\ - 192.8457 \\ - 277.1006 \end{matrix}],$

$S_{2} = [\begin{matrix} 1043.6117 & - 417.4447 & - 208.7223 \\ 285.8274 & - 114.3310 & - 57.1655 \\ 400.7721 & - 160.3088 & - 80.1544 \end{matrix}], M_{2} = [\begin{matrix} - 714.7058 \\ - 195.7458 \\ - 274.4643 \end{matrix}],$

$G_{1} = [\begin{matrix} - 1.3912 & 0.1718 & 0.2018 \\ 0.1867 & - 1.4337 & 0.1858 \\ 0.2126 & 0.1818 & - 1.3971 \end{matrix}] \times 10^{4}, P = {[\begin{matrix} - 104.2037 \\ - 431.6087 \\ 373.6230 \end{matrix}]}^{T},$

$G_{2} = [\begin{matrix} - 1.3912 & 0.1718 & 0.2018 \\ 0.1867 & - 1.4337 & 0.1858 \\ 0.2126 & 0.1818 & - 1.3971 \end{matrix}] \times 10^{4}, Q = 7.5113$ and $T = [\begin{matrix} 0.7216 & 1 & 2 \\ 0.0006 & 2 & 3 \\ 0.1823 & 1 & 4 \end{matrix}] .$

In order to provide a comparison between the performance of the GDFO, we have designed a PFO. By considering matrices $E = I_{3}, L = {ones}_{4, 3} \times 0.1$ and σ = 10 the following PFO matrices are obtained:

$N_{1} = [\begin{matrix} - 0.1472 & - 110.3367 & 0 \\ 0.0001 & - 0.4601 & 0.0006 \\ 0 & - 0.0009 & - 0.0926 \end{matrix}], F_{1} = [\begin{matrix} 0.2287 \\ 0.0006 \\ 0 \end{matrix}],$

$N_{2} = [\begin{matrix} - 0.1472 & - 40.4999 & 0 \\ 0.0001 & - 0.2322 & 0.0006 \\ 0 & - 0.0003 & - 0.0926 \end{matrix}], F_{2} = [\begin{matrix} 0.1147 \\ 0.0003 \\ 0 \end{matrix}],$

$J = [\begin{matrix} 0 \\ 0 \\ 0.0083 \end{matrix}], T = [\begin{matrix} 0.5 & 0 & 0 \\ 0.0016 & 1 & 0 \\ 0 & 0 & 1 \end{matrix}], Q = 1.0959$

and $P = [\begin{matrix} 0.0846 & - 593.7982 & - 8.7418 \end{matrix}] .$

8. Simulation results

In the following section, two simulations are performed to regulate the glucose level in the presence of disturbances due to food intake. In the simulation 1, the disturbance considered is described by equation (44), which represents the rate at which glucose enters the blood from intestinal absorption following a meal. Three disturbances are considered with the amplitude of 15, 19 and 18g.

According to the literature a better way to represent the exogenous glucose input D (t) is by the Rayleigh probability density function, as is described in equation (45), where the value of k depends upon the amount of food consumed. Three disturbances are considered with k = 15 for 10g, k = 18 for 12g and k = 20 for 18g.

Simulation 1

In order to evaluate the performance of the proposed approach, a comparison of the glucose-insulin system response is made considering a constant supply of insulin (open loop) and considering the regulation made by the functional observer, as shown in Fig. 2

Fig. 2

Glucose-insulin system in open and closed-loop.

Figure 3 show disturbance w (t) (Eq. (44)) that is associated with food intake during the course of the day.

Fig. 3

Disturbance due to food intake w (t) from (44).

Figure 4 shows the exogenous insulin supply. The black line represents an open loop system where a constant supply of insulin is applied, with the value of u (t) =27.7778 mU/dL. The blue line is the exogenous insulin estimated by the GDFO and the red one is the control law estimated by the PFO.

Fig. 4

Exogenous insulin signal u (t).

In Fig. 4 it can be seen how the control law compensates the effects of the disturbance in the system. Since the structure of the GDFO is more robust, the control law is greater, therefore it is expected that the GDFO regulate better the glucose level.

Figure 5 shows the regulation of the blood glucose level.

Fig. 5

Blood glucose concentration G (t).

In the first figure, it can be seen that the GDFO and the PFO have almost the same behaviour, this is due to, both observer are designed using the same approach of Theorem 1. In the second figure, which is a zoom of the first disturbance, we can see that the glucose levels increase due to a disturbance by food intake. The insulin supply, shown in Fig. 4 has the task of bring these levels to baseline glucose values. Both observers achieve the regulation objective, however the GDFO reaches a lower glucose value. The reason is that the the GDFO has additional degrees of freedom in its structure, which allows it to obtain a control law (insulin supply) that minimizes the effect of the disturbance better.

Simulation 2

Figure 6 shows the the perturbation w (t) which is taken from equation (45).

Fig. 6

Disturbance due to food intake w (t) from (45).

Figure 7 shows the control laws of exogenous insulin, which are computed by the functional observers. The behaviour of the control laws respond to the supply of insulin to compensate the presence of external glucose due to feeding.

Fig. 7

Control law u (t).

The black line represents the constant supply of insulin (open loop), the blue line is the insulin provided by the GDFO and the red one the insulin estimated by the PFO. In this case, it is important to remark that, because the shape of the disturbance has changed, the amplitude of the control laws are less than in the previous simulation, even if the disturbances are almost the same.

The results presented in Fig. 8 are similar to those of the previous simulation. As the amplitude of the control laws are less, there is a higher level of glucose in the bloodstream. In the zoom figure, it can be seen that, for the smaller disturbance the glucose level in open loop (black line) is close to the hyperglycemia level (greater than $250 [\frac{mg}{dL}]$ ), whereas the control laws in closed loop allow to regulate better the levels of glucose preventing reaching at hyperglycemic levels.

Fig. 8

Comparison of the glucose levels G (t) in open and closed loop.

9. Conclusion

In this paper, an alternative structure of a functional observer for Takagi-Sugeno systems is proposed. The necessary conditions for the existence of this observer are provided, and its stability is proved through the use of linear matrix inequalities. The presented approach unifies also the design for full ans reduced order observers for T-S systems. It was also shown that this observer is a generalization of the existing PFO; this observer is only a particular case. A glucose-insulin system was considered to show the performances of the observers designed. In simulation results, it was proved that GDFO is able to estimate a control law that regulates the glucose-insulin system. An algorithm is provided with the main steps involved in the design of the functional observer.

References

Miao ,

You and

Zhen , A new method for fuzzy system representation and implementation, in Fifth International Conference on Fuzzy Systems and Knowledge Discover, Shandong, China, 2008, pp. 1–5.

Sun ,

Mou ,

Qiu ,

Wang and

Gao , Adaptive fuzzy control for nontriangular structural atochastic switched nonlinear systems with full state constraints, IEEE Transactions on Fuzzy Systems 27(8) (2018), 1587–1601.

Li ,

Zhou and

Xu , Fuzzy control system design via fuzzy lyapunov functions, IEEE Transactiona on systems man, and cybernetics-Part B: Cybernetics 38(6) (2008), 1657–1661.

Qiu ,

Sun ,

Wang and

Gao , Observer-based fuzzy adaptive event-triggered control for pure-feedback nonlinear systems with prescribed performance, IEEE Transactions on Fuzzy Systems (2019), 1.

Pamuar ,

Vasin ,

Atanasković and

Miliić , Planning the city logistics terminal location by applying the green p-median model ans type-2 neurofuzzy network, Computational Intelligence and Neuroscience ID 6972818 (2016), 1–15.

Pamuar ,

Ljubojević ,

Kostadinović and

Dorović , Cost and risk aggregation in multi-objective route planning for hazardous materials transportation - a neuro-fuzzy and artificial bee colony approach, Expert Systems with Applications 65 (2016), 1–15.

Lukovac ,

Pamuar ,

Popović and

Dorović , Portfolio model for analyzing human resources: An approach based on neuro-fuzzy modeling and the simulated annealing algorithm, Expert Systems with Applications 90 (2017), 318–331.

Su ,

Xia ,

Liu and

Wu , Event-triggered fuzzy control of nonlinear systems with its application to inverted pendulum systems, Automatica 94 (2018), 236–248.

Tanaka ,

Ikeda and

H.O.

Wang , Fuzzy regulators and fuzzy observers: Relaxed stability conditions and lmi-based designs, IEEE Transactions on Fuzzy Systems 6(2) (1998), 250–265.

10.

Tanaka ,

Taniguchi and

H.O.

Wang , Model-based fuzzy control of tora system: Fuzzy regulator and fuzzy observer design via lmis that represent decay rate, disturbance rejection, robustness, optimality, in Fuzzy Systems Proceedings, 1998. IEEE World Congress on Computational Intelligence, The 1998 IEEE International Conference on, IEEE, vol. 1, 1998, pp. 313–318.

11.

Tanaka ,

Yoshida ,

Ohtake and

H.O.

Wang , A sumof-squares approach to modeling and control of nonlinear dynamical systems with polynomial fuzzy systems, IEEE Transactions on Fuzzy Systems 17(4) (2009), 911–922.

12.

Chen ,

Shi ,

C.-C.

Lim , et al., Robust constrained control for mimo nonlinear systems based on disturbance observer, IEEE Trans Automat Contr 60(12) (2015), 3281–3286.

13.

Luenberger , Observers for multivariable systems, IEEE Transactions on Automatic Control 11(2) (1966), 190–197.

14.

Darouach , Existence and design of functional observers for linear systems, IEEE Transactions on Automatic Control 45(5) (2000), 940–943.

15.

Moreno , Quasi-unknown input observers for linear systems, in Control Applications, 2001 (CCA’01) Proceedings of the 2001 IEEE International Conference on, IEEE, 2001, pp. 732–737.

16.

T.L.

Fernando ,

H.M.

Trinh and

Jennings , Functional observability and the design of minimum order linear functional observers, IEEE Transactions on Automatic Control 55(5) (2010), 1268–1273.

17.

L.S.

Jennings ,

T.L.

Fernando and

H.M.

Trinh , Existence conditions for functional observability from an eigenspace perspective, IEEE Transactions on Automatic Control 56(12) (2011), 2957–2961.

18.

G.-L.

Osorio-Gordillo ,

Darouach and

C.-M.

Astorga-Zaragoza , h1 dynamical observers design for linear descriptor systems. application to state and unknown input estimation, European Journal of Control 26 (2015), 35–43.

19.

Gao ,

Darouach ,

Voos and

Alma , New unified h1 dynamic observer design for linear systems with unknown inputs, Automatica 65 (2016), 43–52.

20.

Gao ,

Darouach and

Alma , h1 dynamic observer design for discrete-time systems, in Proceedings of the 20th IFAC World Congress, IFAC, Toulouse, France, IEEE, 2017, pp. 2811–2816.

21.

A.-J.

Pérez-Estrada ,

G.-L.

Osorio-Gordillo ,

Darouach and

V.-H.

Olivares-Peregrino , Generalized dynamic observer design for quasi-lpv systems, At-Automatisierungstechnik 66(3) (2018), 225–233.

22.

H.-J.

Marquez , A frequency domain approach to state estimation, Journal of the Franklin Institute 340 (2003), 147–157.

23.

J.-K.

Park ,

D.-R.

Shin and

T.-M.

Chung , Dynamic observer for linear time-invariant systems, Automatica 38(6) (2002), 1083–1087.

24.

Cobelli ,

Renard and

Kovatchev , Artificial pancreas: Past, present, future, Diabetes 60(11) (2011), 2672–2682.

25.

Yang ,

Michel ,

Wang and

Chen , Hybrid controllability of linear switched systems, in IEEE Conference Decicion on Control, 1998, pp. 4350–4355.

26.

Skelton ,

Iwasaki and

Grigoriadis , A unified algebraic approach to linear control design, 1998.

27.

De Gaetano ,

Di Martino ,

Germani and

Manes Mathematical models and state observation of the glucoseinsulin homeostasis, in IFIP Conference on System Modeling and Optimization, Springer, 2003, pp. 281–294.

28.

M.E.

Fisher , A semiclosed-loop algorithm for the control of blood glucose levels in diabetics, IEEE Transactions on Biomedical Engineering 38(1) (1991), 57–61.

29.

C.L.

Chen and

H.W.

Tsai , Modeling the physiological glucose-insulin system on normal and diabetic subjects, Computer Methods and Programs in Biomedicine 97 (2010), 130–140.

30.

Lendek ,

T.M.

Guerra ,

Babuska and

De Schutter , Stability analysis and nonlinear observer design using Takagi-Sugeno fuzzy models, Springer, 2011.