Control in (Subscript[H, 2])/
          (Subscript[H, ∞]) space via Computer Algebra

Béla Palánc(z^1), Levente Kovác(s^2)

\(\(\^1\) Department\ of\ Photogrammetry\ and\ Geoinformatics\)

\(\^2\)Department of Control Engineering and Information Technology

Budapest University of Technology and Economics, 
      H - 1111 Budapest, Hungary, Muegyetem rkp . 3. 

e - mail : 
        palancz@epito . bme . hu

Abstract

This article presents the (Subscript[H, 2])/(Subscript[H, ∞]) control (disturbance rejection LQ method) of the Bergman minimal model [10] for Type I diabetic patients under intensive care using computer algebra. To design the optimal controller, the disturbance rejection LQ method based on the minimax differential game is applied. The critical, minimax value of the scaling parameter (Subscript[γ, crit]) is determined by using the Modified Riccati Control Algebraic (MCARE) equation employing reduced Groebner basis solution on rational field. The numerical results are in good agreement with those of the Control Toolbox of MATLAB. It turned out, that in order to get positive definite solution stabilizing the closed loop, γ should be greater than (Subscript[γ, crit]). The obtained results are compared with the classical LQ technique on the original non-linear system, using a standard meal disturbance situation. It is also demonstrated that for γ >> (Subscript[γ, crit]), the gain matrix approaches the traditional LQ optimal control design solution. The symbolic and numerical computations were carried out with Mathematica 5.2, and with the CSPS Application 2, as well as with MATLAB 6.5.

Keywords: disturbance rejection LQ control, LQ control, glucose-insulin control, diabetes mellitus, symbolic computations, Mathematica.

Introduction

Diabetes mellitus is one of the most serious diseases which need to be artificially regulated. The newest statistics of the World Health Organization (WHO) predate an increase of adult diabetes population from 4% (in 2000, meaning 171 million people) to 5.4% (366 million worldwide) by the year 2030, [1]. This warns that diabetes could be the “disease of the future”.
In many biomedical systems, external controller provides the necessary input, because the human body could not ensure it. The outer control might be partially or fully automatized. The self-regulation has several strict requirements, but once it has been designed it permits not only to facilitate the patient's life suffering from the disease, but also to optimize (if necessary) the amount of the used dosage.
The blood-glucose control is one of the most difficult control problems to be solved in biomedical engineering. One of the main reasons is that patients are extremely diverse in their dynamics and in addition their characteristics are time-varying. Due to the inexistence of an outer control loop, replacing the partially or totally deficient blood-glucose-control system of the human body, patients are regulating their glucose level manually. Based on the measured glucose levels (obtained from extracted blood samples), they decide on their own what is the necessary insulin dosage to be injected. Although this process is supervised by doctors (diabetologists), mishandled situations often appear. Hyper- (deviation over the basal glucose level) and hypoglycemia (deviation under the basal glucose level) are both dangerous cases, but on short term the latter is more dangerous, leading for example to coma.
Starting from the late Sixties lot of researchers investigate the problem of the glucose-insulin interaction and control. The closed-loop glucose regulation as it was several times formulated [2], [3], [4], requires three components and the current paper focuses on the last component of them:
·	Glucose sensor (already realized even for 10 min. frequent readings - MiniMed [5], Glucowatch [6]);
·	Insulin pump, for insulin injection (MiniMed [7], Disetronic [8]);
·	Control algorithm, which based on the glucose measurements, is able to determine the necessary insulin dosage.
To design an appropriate control, an adequate model is necessary. In the last 50 years several models appeared. The mostly used and also the simplest one proved to be the minimal model of Bergman, [9], [10], but its shortcoming is its big sensitivity to variance in the parameters. Henceforward, the plasma insulin concentration must be known as a function of time. Therefore, extensions of this minimal model have been proposed [11], [12], [13], [14] trying to capture the changes in patient dynamics of the glucose-insulin interaction, particularly with respect to insulin sensitivity. Other more general, but more complicated models appeared in the literature, [15], [16].
Regarding the control strategies applied, the palette is very wide [17], starting form classical control strategies like PID control [18], optimal control [19], to the modern control techniques like adaptive control [11], neuro-fuzzy algorithms [20], [21], model predictive control [2], [22], [23], but also post-modern control strategies, like (Subscript[H, ∞]) control [3], [4], [24], (Subscript[H, 2]) / (Subscript[H, ∞]) control [25], [26], μ-synthesis, [27], Linear Parameter Varying (LPV) technique [28].
The investigations of [2] discourage the use of a low complexity control such as PID, if high level of performance is desired. However, probably the best way to approach the problem is to consider the system model and the applied control technique together, [17], [29].
This article presents the (Subscript[H, 2])/(Subscript[H, ∞]) control (disturbance rejection LQ method) of the Bergman minimal model [10] in symbolic way using Mathematica 5.2 together with its Control System Professional Suite (CSPS). The paper is structured firstly of a brief description of the model, which then is reduced to a two-state model eliminating the unmeasurable slow state variable. Secondly, the LQ and disturbance rejection LQ (minimax control, LQR) methods are presented and the symbolic computations are performed to determine a general solution of the considered model for these control methods. Finally, the obtained results are compared on a standard meal disturbance situation for the original non-linear system and checked under MATLAB 6.5 as well.

Model equations

Several different models of diabetic systems exist in the literature including, for example, the very detailed 21st-order metabolic model of Sorensen [16]. However, to have a system that on one hand, can be readily handled from the point of view of control design, but on the other hand represents the biological process properly, we consider the three- state minimal patient model of Bergman [10],

deq1NL = G'[t] == -p1 G[t] - X[t] (G[t] + GB) + 
            h[t];

deq2NL = X'[t] == -p2 X[t] + p3 Y[t];

deq3NL = Y'[t] == -p4 (Y[t] + YB) + 
            i[t]/VL;

where the three state variables: 
	G(t) - Plasma glucose deviation,  [mg/dL]
	X(t) - Remote compartment insulin utilization, [1/min]
	Y(t) - Plasma insulin deviation, [mU/dL]
The control variable:
	i(t) - Exogenous insulin infusion rate, [mU/min]
The disturbance is,
	h(t) - Exogenous glucose infusion rate, [mg/dL min]
The physical parameters:
	GB - Basal glucose level, [mg/dL]
	YB - Basal insulin level, [mU/dL]
	VL - Insulin distribution volume, [dL]
The model parameters: 
	p1 - [1/min]
	p2 - [1/min]
	p3 - [dL/mU (min^2)]
	p4  - [1/min]

The values of the model and physical parameters are from Furler et al [30]:

NumericalValues = {p1 -> 28/1000, p2 -> 25/1000, 
          p3 -> 13/100000, p4 -> 5/54, GB -> 110, YB -> 1.5, 
          VL -> 120};

Model reduction

Assuming that the unmeasurable variable X(t), remote compartment insulin utilization, is a slow variable, namely X'[t] ≈ 0, then

X[t] = (p3/p2) Y[t];

it can be eliminated by substituting it into the first equation. Now our reduced model is

deq1NML = G'[
            t] == -p1 G[t] -  (p3/p2) Y[t] (G[t] + GB) + 
            h[t];

deq2NML = Y'[t] == -p4 (Y[t] + YB) + 
            i[t]/VL;

To design optimal control, the first step is the linearization of the nonlinear model, which has the following general form,

\(x\& .  = f \((x, u)\)\)
y = g (x, u)

In our case

f = {deq1NML[[2]], 
          deq2NML[[2]]};
x = {G[t], Y[t]};
y = x;
u = {h[t], i[t]};
g = y;

Linearized model

In order to linearize our model, the Control System Applicaton of Mathematica should be loaded

<< ControlSystems`

The linearization will  be carried out at the equilibrium point, namely at (X0, Y0, h0, i0). The linearized model as control object in state space is

ControlObjectSS = 
        Linearize[f, 
          g, 
{{G[t], G0}, {Y[t], Y0}}, {{h[t], h0}, {i[t], 
              i0}}];

ControlObjectSS // EquationForm

The steady state values are

h0 = 0; i0 = p4 YB VL; G0 = 0; Y0 = 0;

Then

ControlObjectSS // EquationForm

or in traditional form

ControlObjectSS // TraditionalForm

The system matrices can be easily separated

A = %[[1]]; MatrixForm[A]

B = %%[[2]]; MatrixForm[B]

Let us check the controllability of the linearized system. The controllability matrix is

MC = ControllabilityMatrix[ControlObjectSS]; MatrixForm[MC]

If the rank of this matrix is equal with the rank of matrix A, then the system can be controlled

MatrixRank[MC] == MatrixRank[A]

or we can use built in function directly

Controllable[ControlObjectSS]

Classical LQ control

It is well-known , that the dynamic of a linear time invariant system can be described in the following way

\(TraditionalForm\`x\& . \)(t) = A x(t) + B u(t)

y(t) = C x(t) + D u(t)

where A,B,C and D are constant matrices.
	Using  classical LQ control, the requirement is to minimize the following quadratic cost functional

\(J \((u \((t)\))\) = \ \(1\/2\) \(∫\_0\%∞\( 
                  y\^T\) \((t)\)\ Q\ y \((t)\)\) + \(u\^T\) \((t)\)\ R\ u \((t)\) dt\)

LQ attempts to find  an optimal control  (SuperStar[u])(t), \(t\  ∈ \ \([0, \ ∞]\)\ \ based\ on\ the\ CARE\ \((Control\ Algebraic\ Riccati\ Equation)\)\ \ such\ that\)

\(J \((\(u\^*\) \((t)\))\) ≤ 
      J \((u \((t)\))\)\)

for all u(t) on t ∈ [0, ∞] under properly choosen R and Q matrices.

The first component of input vector u(t), the exogenous glucose (h(t)) stands for disturbance, therefore its effect on the output should be minimized. As a result, the disturbance should be overweighted. Consequently (Subscript[R, 11]) should be considerable greater than (Subscript[R, 22]). We choose the following matrices, R is

R =  DiagonalMatrix[{(10^3),  (10^(-3))}]; 
    Q =  (10^(-3)) IdentityMatrix[2];  MatrixForm[R]

and considering that in the objective function the portion of the injected insulin is important and the values of the state variables play not an important role, the Q matrix is choosen as

MatrixForm[Q]

Then the optimal gain matrix is

KLQ = LQRegulatorGains[ControlObjectSS /. NumericalValues, Q, R]; 
    MatrixForm[KLQ] // N

The first row in K can be neglected

KLQ = Chop[%, 0.001]; MatrixForm[KLQ]

In order to check this result, it was compared with the KLQ  computed by MATLAB:

(({{0, 0}, {-0.4671, 2.6211}}))

Disturbance rejection LQ control (minimax control)

The disturbance rejection method, LQR (LQ rejection) represents a generalization of the classical LQ method and is based on the minimax criteria, where the system dynamics is generally described as before. However, now the disturbance d(t) is separeted from active control input \(TraditionalForm\`u\&_\)(t) and could be considered as unmeasurable, namely

\(\(x\& . \) \((t)\)\)=A x (t)+B (OverBar[u])(t) +L d (t)

y (t) = 
          C x (t) + D u (t)

Therefore, in this case the quadratic cost functional will be modified with the disturbance explicitly, [31],

J (u (t))= (1/2) (∫_0^∞) 
                    (y^T) (t)  y (t)+((OverBar[u])^T)(t)  \(u\&_\)(t)-\(γ\^2\ \(d\^T\) \((t)\)\ d \((t)\)\ dt\)

Now, the disturbace - while it appears with negative sign - attempts to maximize the cost, while we want to find a control \(TraditionalForm\`\(\(u\&_\)\(\ \)\((t)\)\(\ \)\(that\)\(\ \)\(minimizes\)\(\ \)\(the\)\(\ \)\(maximum\)\(\ \)\(cost\)\(\ \)\)\)achievable by the disturbance (by worst case disturbace). This is a case of so called  "worst case" design and leads to the formulation of a min-max differential game, [32],

(Underscript[(max  ), d (t)]) J((OverBar[u])(t),d (t))-> (min_((OverBar[u])(t))) J (u (t), d (t))

where \(u\&_\)(t) and  d(t) are satisfying the state equation .  It can be demonstrated that the unique solution of the differential game, {\(\(u\&_\^*\) \((t)\), \ \ \(d\^*\) \((t)\)\) }  exists and  satisfies the saddle point condition,

J((SuperStar[(OverBar[u])])(t),d (t))≤ J((OverBar[u])(t),d (t))≤J((OverBar[u])(t),\(\(d\^*\) \((t)\)\))

where \(\(u\&_\^*\) \((t)\)\) is the optimal control and  \(d\^*\)(t) is the worst-case disturbance. These functions can be computed as

(SuperStar[\(u\&_\)])(t)=-(B^T) P x (t)

\(\(d\^*\) \((t)\)\  = \ \(1\/γ\^2\) 
            L\^T\ x \((t)\)\)

where P is the positive definite symmetric solution of the Modified Control Algebraic Riccati Equation (MCARE)

P A  +  (A^T) P + (C^T) C - 
        P (B (B^T) - (1/(γ^2)) L (L^T)) P = 
      0

Solution of MCARE

Now, one of the system matrices, B should be modified and a new one, L will be introduced.

For control input (insulin)

B=(({{0, 0}, {0, (1/VL)}}));

For disturbance (glucose)

L=(({{1, 0}, {0, 0}}));

We are looking for the symmetric solution matrix of the modified Riccati equations in the following form,

P = Table[If[i > j, (Subscript[p, j, i]), (Subscript[p, i, j])], {i, 1, 2}, {j, 1, 2}]; 
    MatrixForm[P]

Keeping in mind, that in our case the C is an identity matrix, the left hand side of the Riccati equation is

RI = 
        P . A + Transpose[A] . P + IdentityMatrix[2] - 
          P . (B . Transpose[B] - (1/(γ^2)) L . Transpose[L]) . 
            P;

MatrixForm[RI]

It is easy to see that RI is a symmetric matrix

Apply[And, 
      Table[RI[[i, j]] == RI[[j, i]], {i, 1, 2}, {j, 1, 2}] // 
        Flatten]

Therefore,to solve the modified Ricatti equation, R I= (({{0, 0}, {0, 0}})),\(only\ three\ equations\ should\ be\ solved\ with\ the\ γ\ as\ parameter\) , namely

eqs = Table[RI[[i, j]] == 0, {i, 1, 2}, {j, i, 2}] // 
          Flatten // FullSimplify

The critical solution of this system belongs to the critical value of the parameter γ. Crossing this critical value with γ the solution, which was real becomes imaginary and vice versa. On one hand finding this critical solution numerically is an ill-conditioned problem, on the other hand to solve this polynomial system in fully symbolic form is hard to manage. Therefore, we used reduced Groebner basis solution on rational field. Substituting the values of the system parameters in rational form, we get

Eqs = eqs /. NumericalValues

To get the solution for  (Subscript[p, 1, 1]), reduced Groebner basis can be employed,

grb21 = 
      GroebnerBasis[
          Eqs, {(Subscript[p, 1, 1]), (Subscript[p, 1, 2]), (Subscript[p, 2, 2])}, {(Subscript[p, 1, 2]), 
            (Subscript[p, 2, 2])}] // Expand

This is a monomial with parameter γ. We can create a function provinding the first root of this monomial, where the independent variable is γ.

p11crit[γ0_] := 
      (Subscript[p, 1, 1]) /. Solve[(grb21[[1]] /. γ -> γ0) == 0, (Subscript[p, 1, 1])][[1]]

The typical graph of the imaginary part of the solution can be seen on Fig.1/a

Plot[Im[p11crit[γ0]], {γ0, -20, 20}, 
        PlotStyle -> {RGBColor[1, 0, 0], Thickness[0.007]}, 
        PlotRange -> All];

Fig .1 /a   The imaginary part of (Subscript[p, 1, 1]) depending on γ

The critical value may computed using NMinimize

γcrit = γ0 /. NMinimize[{Im[
                p11crit[γ0]], γ0 > 
                0}, γ0][[2]]

Plot[Im[p11crit[γ0]], {γ0, -20, 20}, 
        PlotStyle -> {RGBColor[1, 0, 0], Thickness[0.007]}, 
        PlotRange -> All, 
        Epilog -> {AbsolutePointSize[5], RGBColor[0, 0, 1], 
            Point[{γcrit, Im[p11crit[γcrit]]}]}];

Fig .1 /b   Approximate solution for (Subscript[γ, crit])

However, this is just a good approximation for our problem, because

p11crit[γcrit]

see Fig.1/b.

Tuning this approximation a little by trial-error method, we get,

γcrit = 17.117425694;

Indeed, now the imaginary part of the solution is zero,

Im[p11crit[γcrit]]

Employing MATLAB minmax control design we have got the same result, namely, the optimal, minimal γ, (Subscript[(γ ), min])=17.11743 is equal to the critical value (Subscript[γ, crit]) computed above, see Fig.1/c. However at this value γ = (Subscript[(γ ), min]) = (Subscript[γ, crit]), there is no positive definite solution (we are considering only the control input effect feedback!), namely using

Nγ = γcrit

Plot[Im[p11crit[γ0]], {γ0, -20, 20}, 
        PlotStyle -> {RGBColor[1, 0, 0], Thickness[0.007]}, 
        PlotRange -> All, 
        Epilog -> {AbsolutePointSize[5], RGBColor[0, 0, 1], 
            Point[{Nγ, Im[p11crit[Nγ]]}]}];

Fig .1 /c   Proper  solution for (Subscript[γ, crit]) but failing  to ensure positive definite solution

For this Nγ value, the modified Ricatti can be solved numerically

NsolR = 
        Solve[eqs /. NumericalValues /. γ -> 
              Nγ, {(Subscript[p, 1, 1]), (Subscript[p, 1, 2]), 
            (Subscript[p, 2, 2])}];

Considering the eigenvalues, none of the solutions is positive definite,

Map[Eigenvalues[#] &, P /. NsolR] // TableForm // N

However, we choose the second solution like MATLAB does

solPD = P /.  NsolR[[2]] // N

MatrixForm[solPD] // N

MATLAB computed practically the same solution

solPDM=(({{-169.5, 1250.3}, {1250.3, -8950.6}}));

indeed, this is not positive definite solution,

Eigenvalues[solPDM] // N

although this solution satisfies the modified Ricatti equation

Map[
        Norm[RI /. 
                Join[NumericalValues, {(Subscript[p, 1, 1]) -> #[[1, 1]], 
                    (Subscript[p, 1, 2]) -> #[[1, 2]], 
                    (Subscript[p, 2, 2]) -> #[[2, 2]], γ -> 
                      Nγ}]]/Norm[#] &, {solPD}] // First

Considering this solution of the Riccati equation, the corresponding gain matrix can be computed as

KLQR = Transpose[B] . solPD /. NumericalValues // N;

MatrixForm[KLQR]

The MATLAB result is practically the same

(({{0, 0}, {10.4092, -74.5880}}))

Let consider the gain of the closed loop

(A - B . KLQR) /. NumericalValues // N

Eigenvalues[%]

It means that in this case, if we want to control the system only by the control input, this control is not stabilizing the system!

Let us increasing the value of γ up to the value provinding positive definite solution. Increasing somewhat the value of  γ  by  0.25 % , we get

Nγ = γcrit 1.0025

Fig.1/d demonstrates, how sensitive the solution. Visually, the two value of γ represented by the same point on Fig.1/d.

For this γ, the solution is real

Im[p11crit[Nγ]]

and there is a positive definite solution arising

NsolR = 
        Solve[eqs /. NumericalValues /. γ -> 
              Nγ, {(Subscript[p, 1, 1]), (Subscript[p, 1, 2]), 
            (Subscript[p, 2, 2])}];

Map[Eigenvalues[#] &, P /. NsolR] // TableForm // N

Plot[Im[p11crit[γ0]], {γ0, -20, 20}, 
        PlotStyle -> {RGBColor[1, 0, 0], Thickness[0.007]}, 
        PlotRange -> All, 
        Epilog -> {AbsolutePointSize[5], RGBColor[0, 0, 1], 
            Point[{Nγ, Im[p11crit[Nγ]]}]}];

Fig .1 /d   Proper  solution for (Subscript[γ, crit]) providing positive definite solution

Considering the last solution

solPD = P /.  Last[NsolR] // N

MatrixForm[solPD] // N

It satisfies the modified Ricatti equation

Map[
        Norm[RI /. 
                Join[NumericalValues, {(Subscript[p, 1, 1]) -> #[[1, 1]], 
                    (Subscript[p, 1, 2]) -> #[[1, 2]], 
                    (Subscript[p, 2, 2]) -> #[[2, 2]], γ -> 
                      Nγ}]]/Norm[#] &, {solPD}] // First

The control matrix is

KLQR = Transpose[B] . solPD /. NumericalValues // N;

MatrixForm[KLQR]

The closed loop gain matrix

(A - B . KLQR) /. NumericalValues // N

which is now stabilizing the system

Eigenvalues[%]

Comparing classical LQ and disturbance rejection LQ control design

The gain matrix  KLQR,  provided by the disturbance rejection LQ design depends on the actual value of the scaling parameter γ. In case γ -> ∞, we get the gain matrix designed by LQ method [25], [26], namely

(lim_(γ ->  ∞)) KLQR  =  
      KLQ

To demonstrate this we have considered γ = 100 (Subscript[γ, crit]), and we repeated the computations for solving Riccati equation,

Nγ100 = SetPrecision[100 Nγ, 20]

Now the modified Ricatti can be solved numerically

NsolR = 
        Solve[eqs /. NumericalValues /. γ -> 
              Nγ100, {(Subscript[p, 1, 1]), (Subscript[p, 1, 2]), 
            (Subscript[p, 2, 2])}];

The solutions providing positive definite P matrix should be selected

Map[Eigenvalues[#] &, P /. NsolR] // TableForm // N

Therefore, the last solution should be selected

solPD = P /. Last[ NsolR] // N

MatrixForm[solPD] // N

This solution satisfies the modified Ricatti equation

Map[
        Norm[RI /. 
                Join[NumericalValues, {(Subscript[p, 1, 1]) -> #[[1, 1]], 
                    (Subscript[p, 1, 2]) -> #[[1, 2]], 
                    (Subscript[p, 2, 2]) -> #[[2, 2]], γ -> 
                      Nγ100}]]/Norm[#] &, {solPD}] // First

Employing MATLAB minmax control procedure in case of  γ = 100 (Subscript[γ, min]) = 100 17.1174 = 1711.74, the solution of the modified Ricatti equation is:

(({{13.9606, -56.0634}, {-56.0634, 314.6227}}))

which is very close to the result above.

Having the solution of the Riccati equation, the corresponding gain matrix can be computed as

KLQRinf = Transpose[B] . solPD /. NumericalValues;

MatrixForm[KLQRinf] // N

which is indeed very close to KLQ, see above.

KLQ // MatrixForm

Nonlinear system performance with LQ control

The performance of the control is tested by using a standard meal disturbance with about six our duration, modelled by Lehmann and Deutsch, [34]. The following data pairs are considered,

xu = {0, 25, 50, 75, 100, 125, 150, 175, 200, 225, 250, 275, 300, 
          325, 350, 375, 400, 450, 1000};

yu = {0, 0.185, 0.495, 0.765, 0.975, 1.15, 1.07, 0.82, 0.575, 0.335, 
          0.225, 0.145, 0.098, 0.06, 0.035, 0.02, 0.01, 0.005, 
          0};

H = 
      Interpolation[Transpose[{xu, yu}], 
        InterpolationOrder -> 2]

\(\(\(
\)\(p2u = 
        Plot[H[t], {t, 0, 600}, 
          AxesLabel -> {"t[min]", \*""H[\!\(mg\/\(dl\\ min\)\)]""}, PlotStyle -> RGBColor[0, 0, 1]];\)\)\)

Fig. 2.  Exogenous glucose infusion, H(t)

As a result, the disturbance function is

h[t_] = H[t];

The steady state insulin injection is

i0 = i0 /. NumericalValues

Considering negative feedback, the insulin injection rate as control variable is

i[t_] := i0 - KLQ[[2]] . {G[t], Y[t]} /. NumericalValues

Now, we consider the original 3-states nonlinear model, instead of the reduced 2-states linear model which has been employed for control design. The cloosed loop equations are

Clear[X]

deq1NL = G'[
            t] == (-p1 G[t] - X[t] (G[t] + GB) + h[t]) /. 
        NumericalValues

deq2NL = X'[t] == (-p2 X[t] + p3 Y[t]) /. 
        NumericalValues

deq3NL = Y'[t] == (-p4 (Y[t] + YB) + i[t]/VL) /. 
          NumericalValues // N

Now we solve this differential equation system consisting of the three states variables numerically employing implicit backward differentiation formulas,

sol = 
      NDSolve[{deq1NL, deq2NL, deq3NL, X[0] == 0, Y [0] == 0, 
          G[0] == 0}, {X, Y, G}, {t, 0, 1000}, Method -> BDF,  
        MaxSteps -> 50000]

The controlled dynamics of the blood glucose can be seen on Fig. 3

pGLQ = 
        Plot[Evaluate[(G[t] /. sol)], {t, 0, 1000}, 
          PlotStyle -> {RGBColor[0, 0, 1]}, Frame -> True, 
          PlotRange -> All, 
          FrameLabel -> {"t[min]", "G[mg/dl]"}];

Fig. 3 The controlled dynamics of blood glucose concentration, G(t) in case of LQ control

The performance of the corresponding insulin infusion rate, i(t) as control variable is,

piLQ = 
        Plot[Evaluate[i[t]] /. sol, {t, 0, 1000}, 
          PlotStyle -> {RGBColor[0, 0, 1]}, Frame -> True, 
          PlotRange -> All, 
          FrameLabel -> {"t[min]", "i[mU/min]"}];

Fig.4 The corresponding insulin infusion rate, i(t) in case of LQ control

Nonlinear system performance with disturbance rejection LQ control

In our case, it has no physical meaning improving the glucose control by exhausting glucose from the system (some kind of negativ injection, or glucose substraction form the body), see the structure of matrix B in MCARE. In real situation only insuline injection is possible. Therefore we can use only the second row of the control matrix.

Now, the insulin injection rate as control variable is

i[t_] := 
      i0 - KLQR[[2]] . {G[t], Y[t]} /. NumericalValues

Clear[X]

The modell equations are

deq1NL = G'[
            t] == (-p1 G[t] - X[t] (G[t] + GB) + h[t]) /. 
        NumericalValues

deq2NL = X'[t] == (-p2 X[t] + p3 Y[t]) /. 
        NumericalValues

deq3NL = Y'[t] == (-p4 (Y[t] + YB) + i[t]/VL) /. 
          NumericalValues // N

Now we solve this differential equation system consisting of the three states variables numerically employing implicit backward differentiation formulas,

sol = 
      NDSolve[{deq1NL, deq2NL, deq3NL, X[0] == 0, Y [0] == 0, 
          G[0] == 0}, {X, Y, G}, {t, 0, 1000}, Method -> BDF,  
        MaxSteps -> 50000]

The controlled dynamics of the blood glucose can be seen on Fig. 5

pGLQR = 
        Plot[Evaluate[(G[t] /. sol)], {t, 0, 1000}, 
          PlotStyle -> {RGBColor[0, 1, 0]}, Frame -> True, 
          PlotRange -> All, 
          FrameLabel -> {"t[min]", "G[mg/dl]"}];

Fig. 5 The controlled dynamics of blood glucose concentration, G(t) in case of disturbance rejection LQ control

The performance of the corresponding insulin infusion rate, i(t) as control variable is,

piLQR = 
        Plot[Evaluate[i[t]] /. sol, {t, 0, 1000}, 
          PlotStyle -> {RGBColor[0, 1, 0]}, Frame -> True, 
          PlotRange -> All, 
          FrameLabel -> {"t[min]", "i[mU/min]"}];

Fig. 6 The corresponding insulin infusion rate, i(t) in case of disturbance rejection LQ control

Comparing LQ and disturbance rejection LQ control performance

Comparing the results of the two controls, it can be seen that even in this considered case, the disturbance rejection LQ control is more efficient than LQ,

Show[{pGLQ, pGLQR}];

Fig. 7 The comparison of  blood glucose concentration, G(t) in case of LQ and disturbance rejection LQ control

although the latest need higher rate of insulin infusion

Show[{piLQ, piLQR}];

Fig.8  The corresponding insulin infusion rate, i(t) in case of LQ and disturbance rejection LQ control

Conclusions

It turned out that the critical value of γ (where the value for imaginary parts of the solution of MCARE disappear) togehter with the physically realizable interpretation of a control system, will not automatically ensure positive definite solution. It has been demonstrated that "blind" numerical solution of disturbance rejection LQ problem can be misleading, and not always provides positive definite solution. This problem can be detected by using symbolic-numeric solution and increasing the value of γ up to the very value, which already ensures positive eigenvalues of the Ricatti matrix. The other advantage of the symbolic-numeric solution, whenever it is possible, is its robustness concerning round-off errors. However, this problem can be avoided by using another numerical technique for minimax control based on robust control on frequency domain proposed by Helton [33], and its application for glucose-insulin control see [24]. It was also illustrated, that disturbance rejection LQ control could provide better control quality than LQ does, see Fig.7.  Disturbance rejection LQ control is interacting faster as well as employing higher infusion rate (in the considered case) than the classical control, see Fig.8.

References

[1] 	S. Wild, G. Roglic, A. Green, R. Sicree and H. King, “Global Prevalence of Diabetes - Estimates for the year 2000 and projections for 2030”, Diabetes Care, vol. 27 (5), pp. 1047-1053, 2004.
[2] 	N. Hernjak, F. J. Doyle III, “Glucose control Design Using Nonlinearity Assessment Techniques”, AIChE Journal, vol. 51, no. 2, pp. 544-554, Febr. 2005.
[3] 	R. S. Parker, F. J. Doyle III, J. H. Ward, N. A. Peppas, “Robust (Subscript[H, ∞]) Glucose Control in Diabetes Using a Physiological Model”, AIChE Journal, vol. 46 (12), pp. 2537-2549, 2000.
[4] 	E. Ruiz-Vel!azquez, R. Femat, D. U. Campos-Delgado, „Blood glucose control for type I diabetes mellitus: A robust tracking (Subscript[H, ∞]) problem”, Elsevier Control Engineering Practice, vol. 12, pp. 1179-1195, 2004.
[5] 	Mini Med CGMS http://www.minimed.com/patientfam/pf_products_cgms_ov_completetic.shtml.
[6] 	Glucowatch, http://www.glucowatch.com.
[7] 	Mini Med Paradigm Insulin Pump http://www.minimed.com/patientfam/pf_ipt_paradigm_insulin_pump.shtml.
[8] 	Disetronic, http://www.disetronic-usa.com/insulin-pumps.htm.
[9] 	B. N. Bergman, Y. Z. Ider, C. R. Bowden, C. Cobelli, “Quantitive estimation of insulin sensitivity”, American Journal of Physiology, vol. 236, pp. 667-677, 1979.
[10] 	R. N. Bergman, L. S. Philips, C. Cobelli, “Physiologic evaluation of factors controlling glucose tolerance in man”, Journal of Clinical Investigation, vol. 68; pp. 1456-1467, 1981.
[11]	 J. Lin, J. G. Chase, G. M. Shaw, C. V. Doran, C. E. Hann, M. B. Robertson, P. M. Browne, T. Lotz, G. C. Wake, B. Broughton, “Adaptive Bolus-Based Set-Point Regulation of Hyperglycemia in Critical 	Care”, Proc. of 26th Ann. Int. Conf. of IEEE Eng. in Biomedicine Soc., San Francisco, USA, pp. 3463-3466, 2004.
[12] 	M. Fernandez, D. Acosta, M. Villasana, D. Streja, “Enhancing Parameter Precision and the Minimal Modeling Approach in Type I Diabetes”, Proc. of 26th Ann. Int. Conf. of IEEE Eng. in Biomedicine Soc., 	San Francisco, USA, pp. 797-800, 2004.
[13] 	H. C. Morris, B. O’Reilly, D. Streja, “A New Biphasic Minimal Model”, Proc. of 26th Ann. Int. Conf. of IEEE Eng. in Biomedicine Soc., San Francisco, USA, pp. 782-785, 2004.
[14] 	A. de Gaetano, O. Arino, “Some considerations on the mathematical modeling of the Intra-Venous Glucose Tolerance Test”, Journal of Mathematical Biology, vol 40, pp. 136-168, 2000.
[15] 	R. Hovorka, F. Shojaee-Moradie, P. V. Carroll, L. J. Chassin, I. J. Gowrie, N. C. Jackson, R. S. Tudor, A. M. Umpleby, R. H. Jones, “Partitioning glucose distribution/transport, disposal, and endogenous 	production during IVGTT”, American Journal of Physiology - Endocrinology and Metabolism, vol. 282, pp. E992-E1007.
[16] 	J. T. Sorensen, “A Physiologic Model of Glucose Metabolism in Man and Its Use to Design and Assess Improved Insulin Therapies for Diabetes”, PhD Thesis, Dept. of Chemical Engineering, MIT, 		Cambridge, MA, USA, 1985.
[17] 	R. Parker, F. J. Doyle III, N. A. Peppas, “The Intravenous Route to Blood Glucose Control”, IEEE Engineering in Medicine and Biology, vol. 20 (1), pp. 65-73, 2001.
[18] 	F. Chee, T. L. Fernando, A. V. Savkin, V. van Heeden, “Expert PID Control System for Blood Glucose Control in Critically Ill Patients”, IEEE Transactions on Information Technology in Biomedicine, vol. 7 	(4), pp. 419-425, 2003.
[19] 	M. S. Ibbini, M. A. Masadeh, M.M. Bani Amer, “A semiclosed-loop optimal control system for blood glucose level in diabetics”, Journal of Medical Engineering and Technology, vol. 28 (5), pp. 189-195, 	2004.
[20] 	D. Dazzia, F. Taddei, A. Gavarini, E. Uggeri, R. Negro, A. Pezzarossa, “The control of blood glucose in the critical diabetic patient - A neuro-fuzzy method”, Elsevier Journal of Diabetes and its 			Complications, vol. 15, pp. 80-87, 2001.
[21] 	M. Ibbini, “A PI-fuzzy logic controller for the regulation of blood glucose level in diabetic patients”, Journal of Medical Engineering and Technology, vol. 30 (2), pp. 83-92, 2006.
[22] 	R. Hovorka, V. Canonico, L. J. Chassin, U. Haueter, M. Massi-Benedetti, M. Orsini Federici, T. R. Pieber, H. C. Schaller, L. Schaupp, T. Vering, M. E. Wilinska, “Nonlinear model predictive control of glucose 	concentration in subjects with type 1 diabetes”, Physiological measurement, Vol. 25, pp. 905-920, 2004.
[23] 	S. M. Lynch, B. W. Bequette, “Model predictive control of blood glucose in type I diabetics using subcutaneous glucose measurements”, American Control Conference, Anchorage (USA), vol. 5, pp. 4039- 	4043, 2002.
[24] 	L. Kovács, B. Paláncz, B. Benyó, L. Török, Z. Benyó, “Robust Blood-Glucose Control using Mathematica”, EMBC2006 - 28th Annual International Conference of the IEEE Engineering in Medicine and 	Biology Society, New York, USA, pp. 451-454, 2006.
[25] 	L. Kovács, B. Paláncz, Z. Benyó, “Classical and Modern Control Strategies in Glucose-Insulin Stabilization”, 16th IFAC World Congress, Prague, Czech Republic, electronic publication #04165, 2005.
[26] 	L. Kovács, B. Paláncz, Zs. Almássy, Z. Benyó, “Optimal Glucose-Insulin Control in H2 Space”, EMBC2004 - 26th Annual International Conference of the IEEE Engineering in Medicine and Biology Society, 	San Francisco, USA, pp. 762-765, 2004.
[27] 	L. Kovács, B. Kulcsár, Z. Benyó, “On The Use Of Robust Servo Control In Diabetes Under Intensive Care”, SACI2006 - 3rd Romanian-Hungarian Joint Symposium on Applied Computational Intelligence, 	Timisoara, Romania, pp. 236 - 247, 2006.
[28] 	L. Kovács, B. Kulcsár, J. Bokor, Z. Benyó, “LPV Fault Detection of Glucos-Insulin System”, MED06 - 14th Mediterranean Conference on Control and Automation, Ancona, Italy, electronic publication 	TLA2-4, 2006.
[29] 	A. Makroglou, J. Li, Y. Kuang, “Mathematical models and software tools for the glucose - insulin regulatory system and diabetes: an overview”, Elsevier Applied Numerical Mathematics, vol. 56 (3-4), pp. 	559-573, 2006.
[30] 	S. M. Furler, E. W. Kraegen, R. H. Smallwood, D. J. Chisolm, “Blood glucose control by intermittent loop closure in the basal mode: Computer simulation studies with a diabetic model”, Diabetes Care, vol. 8, 	pp. 553-561, 1985.
[31]	K. Zhou, "Robust and Optimal Control", Prentice Hall, New Jersey, pp. 413-497, 1996.
[32]	J. Bokor, "LQ and disturbance rejection", Modern Control Theory II, PhD lecture notes, internal document, Budapest, 2003.
[33] 	J. W. Helton, O. Merino, "Classical Control Using H∞ Methods - Theory, Optimization, and Design", SIAM, pp.17-63, 1998.
[34] 	E. D. Lehmann, T. A. Deutsch, “A physiological model of glucose-insulin interaction in type I diabetes mellitus”, Journal of Biomedical Engineering, vol. 14, pp. 235-242, 1992.