Optimal  control in (Subscript[H, 2]) space

B.  Palánc\(TraditionalForm\`\(z\^*\)\) and  L. Kovác\(TraditionalForm\`\(s\^ ** \)\)

\(TraditionalForm\`\(\^*\)\) Department of Photogrammetry and Geoinformatics
\(TraditionalForm\`\(\( . \^\(\(**\)\(\ \)\)\)\(Department\)\)\)of Control Engineering and Information Technology
 Budapest University of Technology  and Economics, 
 H-1111 Budapest, Hungary, Müegyetem rkp. 3.
 
 e-mail:  palancz@epito.bme.hu

Abstract

In this case study, an optimal  control in  (Subscript[H, 2]) space is presented for glucose - insulin system of diabetic patients under intensive care. The analysis is based on a modified two-compartment model. To design optimal controller, the disturbance rejection,  LQR method  based on minimax  differential game is applied. The critical, minimax value of the scaling parameter  (Subscript[γ, crit]),  is determined by symbolic solution of the modified Ricatti equation. The numeric evaluation of the symbolic computation for γ > (Subscript[γ, crit]) leads two  different solutions, but  the norms of the vectors {(Subscript[λ, 1]), (Subscript[λ, 2])} formed by the eigenvalues of the pair of the gain matrices are the same. The numerical results are in good agreement with that of  the μ-Toolbox of MATLAB, however this latest gives only one of the two solutions. One of the gain matrices with increasing γ,  approaches the gain matrix computed with the traditional LQ optimal control design. The symbolic and numerical computations were carried out with Mathematica 5. and with the CSPS Application 2. as well as with  MATLAB 6.5

Introduction

From engineering point of view, treatment of diabetes mellitus can be represented by outer control loop to replace the partially or totally failing blood-glucose control system of human body. To maintain the glucose level in a diabetic patient under intensive care is currently an actively researched topic in the field of Biomedical Engineering. To design an appropriate control,  an adequate model  is necessary. Towards the last 50 years,  a variety of models  for the interaction  between glucose and insulin have been suggested in the  literature  [1], [2], [3], [4] and [5] as well as control strategies have been designed and applied to the problem [6], [7], [8], [9], [10] and [11]. Most of the models were realized for "artificial pancreas" function, in such conditions, where the patient blood glucose level being monitored and insulin injection are performed continuously during surgery.
	The aim of the authors was to select a model which is not so complicated, however it is able to properly describe the performance of the physiological system to be controlled. Therefore comparing the mentioned models, the authors choose a modified two- compartment model [2], considering it as the best appropriate model.
	Improving the control strategy an optimal glucose-insulin control in Hardy (Subscript[H, 2])space (LQ, LQR control) has been designed and  simulation the effect of food (sugar) intake carried out.

Model equations

To simulate the insulin-glucose interaction in human body, the following two-compartment model was employed [2]

deq1 = X'[t] == (Subscript[p, 1]) X[t] + (Subscript[p, 2]) h[t] ;

deq2 = Y'[t] == ((Subscript[p, 3]) - X[t]) Y[t] + i[t] + 
            (Subscript[p, 4]);

The terms h(t) and i(t) stand for  exogenous glucose and insulin as the inputs of the system. The state variables X(t) and Y(t) are the concentration of glucose in the plasma and that of the insulin remote from plasma, respectively. In our case X(t) and Y(t) represent both the state and the output of the system, because the dynamic of the measurement and actuator devices are considerably faster than that of the system itself.  The constants (Subscript[p, i]),  i = 1, 2,...,4 are the model parameters.
	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 = {deq1[[2]], deq2[[2]]}; x = {X[t], Y[t]}; y = x; 
    u = {h[t], i[t]}; g = y;

Linearized model

In order to linearize our model, Control Systems Professional for Mathematica should be loaded

<< ControlSystems`

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

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

Or as is usual, in equation form

ControlObjectSS // EquationForm

Equilibrium state

Considering, that at t = 0 the system has no input,

i(0 ) = 0,  h(0) = 0

then from the first state equation,  we get

X(0) = 0

therefore the equilibrium values are

data = {i[0] -> 0, h[0] -> 0, Y[0] -> Y0, 
          X[0] -> 0};

The equilibrium value, Y0 can computed from the second state equation,

eq = 0 == (f[[2]] /. t -> 0) /. data

solving this for Y0,

sol = Solve[eq, Y0]

Inserting the equilibrium values  into our linearized model,

X0 = 0; Y0 = Y0 /. sol[[1]];

we get

ControlObjectSS // EquationForm

or in traditional form

ControlObjectSS // TraditionalForm

The system matrices A and B are

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

and

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

Controllability Test

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]

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\ \ \(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 chosen R and Q matrices. 
Hereafter the system considered is characterized by the following parameter values [12],

numericalValues = {(Subscript[p, 1]) -> -0.021151, (Subscript[p, 2]) -> 0.092551, 
          (Subscript[p, 3]) -> -0.014188, (Subscript[p, 4]) -> 0.077947};

The first component of u(t), the exogenous glucose, h(t) stands for disturbance, therefore it should be eliminated from the LQ control. 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 the Q matrix

MatrixForm[Q]

Then the  optimal gain matrix is

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

The first row in K can be neglected

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

The eigenvalues of the closed loop are

Eigenvalues[(A - B . KLQ) /. numericalValues]

They are negative, therefore the control stabilize the system.

LQR control

The disturbance rejection method, LQR (LQ rejection) represents a generalization of the classical LQ method and based on the minimax criteria in Hardy (Subscript[H, 2]) space, where the system dynamics is generally described as before. However, now the disturbance d(t) is separated from active control input \(TraditionalForm\`u\&_\)(t) and  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,

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 disturbance -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 disturbance). This is a case of so called  "worst case" design and leads to the formulation of a min-max differential game [13],

(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.  According to [13], 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) [13]

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

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

For disturbance

L=(({{(Subscript[p, 2]), 0}, {0, 0}}));

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

P=(({{P11, P12}, {P12, P22}}));

Keeping in mind,  that in our case the C is 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

RI[[1, 2]] == RI[[2, 1]]

Therefore,to solve the modified Ricatti equation, R I= (({{0, 0}, {0, 0}})),only three equations to be solved, namely

eq1 = RI[[1, 1]] == 0

eq2 = RI[[1, 2]] == 0

eq3 = RI[[2, 2]] == 0

However, to solve these three nonlinear equations in symbolic form is not an easy task. We use a trick, now solving only two equations for  unknowns P11 and P22, with  the third one, P12 as a parameter:

solR = Solve[{eq1, eq2}, {P11, P22}]

Because the elements of matrix P should be real numbers, the discriminant (it is the same in every solution)

dRI = 
        4 (γ^4) (Subscript[p, 1]^2) (Subscript[p, 3]^2) - 
          4 (Subscript[p, 2]^2) (Subscript[p, 3]) ((γ^2) (Subscript[p, 3]) - 
                (P12^2) (γ^2) (Subscript[p, 3]) + 
                2 P12 (γ^2) (Subscript[p, 4]));

should be positive, dRI ≥ 0. Now we consider the worst case, namely dRI = 0

solP12 = Solve[dRI == 0, P12]

Using the same technique, the critical value of γ  can be computed from dP12 = 0, where

dP12 = -(γ^2) (Subscript[p, 1]^2) (Subscript[p, 2]^2) (Subscript[p, 3]^2) + 
          (Subscript[p, 2]^4) (Subscript[p, 3]^2) + (Subscript[p, 2]^4) (Subscript[p, 4]^2);

Solving this equation for γ

solγcrt = Solve[dP12 == 0, γ]

The value of γ should be positive, therefore only the first solution can be considered (model parameter (Subscript[p, 2]) is positive),

therefore

γcrit = 
      Simplify[(( (Sqrt[(Subscript[p, 2]^2) (Subscript[p, 3]^2) + (Subscript[p, 2]^2) (Subscript[p, 4]^2)]))/((Subscript[p, 1]) (Subscript[p, 3]))), Assumptions -> (Subscript[p, 2]) > 0]

In case of γ = (Subscript[γ, crit]), P = (Subscript[P , crit]). The elements of the critical solution matrix are

Clear[P11, P12, P22, P11crit, P12crit, P22crit]

P12crit = ((P12 /. 
                solP12[[1]]) /. γ -> γcrit) // 
        Simplify

and then the critical solution for P11

P11crit = (P11 /. solR[[1]]) /. {γ -> γcrit, 
            P12 -> P12crit} // Simplify

but for  P22,  the critical solution is blowing up:

P22crit = (P22 /. solR[[1]]) /. {γ -> γcrit, 
            P12 -> P12crit} // Simplify

which means, that the solution of the modified Riccati equation at γ = (Subscript[γ, crit]) is a singular solution.

The optimal value of γ is just greater than (Subscript[γ, crit]), because substituting Eq. (11) into Eq. (7) in order to achieve Eq. (8), the smallest  admissible γ should be chosen.

The modified Ricatti equation for matrix P with numerical values of the model parameters is,

eqs[γ_] = 
      Map[(# /. numericalValues) &, {eq1, eq2, eq3}]

Employing the numerical values of the model parameter, the critical value of γ

Nγcrit = 
      SetPrecision[γcrit /. numericalValues, 20]

Let us consider  γ =( 1 + ϵ)  (Subscript[γ, crit]), where ϵ is a small number, for example,  ϵ = (10^(-10))

ϵ = (10^(-10));

then the corresponding admissible γ value is,

Nγ = 
      SetPrecision[(1 + ϵ) Nγcrit, 20]

Now the modified Ricatti can be solved numerically

NsolR = Solve[eqs[Nγ], {P11, P12, P22}]

The solutions  providing  positive definite P matrix should be selected

solPD =  
      Select[{{P11, P12}, {P12, P22}} /. 
          NsolR, (Eigenvalues[#][[1]] > 
                0 ∧ Eigenvalues[#][[2]] > 0) &]

indeed these solutions are positive definite,

Map[Eigenvalues[#] &, solPD]

In order to use matrix norm, the following standard package should be loaded

<< LinearAlgebra`MatrixManipulation`

These solutions satisfy the modified Ricatti equation

Map[
      MatrixNorm[
            RI /. Join[
                numericalValues, {P11 -> #[[1, 1]], 
                  P12 -> #[[1, 2]], 
                  P22 -> #[[2, 2]], γ -> Nγ}]]/
          MatrixNorm[#] &, solPD]

The L matrix is

NL = L /. numericalValues

Having the solutions of the Riccati equation, the corresponding gain matrices can be computed as

K = 
      Map[(Transpose[B] . # + (1/N)(γ^2) Transpose[NL] . #) &, 
        solPD]; Map[MatrixForm[#] &, K]

These two solutions are fairly close to each other, so we consider the average value of them

KLQR = (K[[1]] + K[[2]])/2; MatrixForm[KLQR]

The eigenvalues of the closed loop are

Eigenvalues[(A /. numericalValues) - B . KLQR]

We should mentioned that the μ-Toolbox of MATLAB provides only numerical value for (Subscript[γ, crit]), and only one solution for the gain matrix, which is very closed to the average gain KLQR, namely the MATLAB solution for the gain matrix is [14]

(({{0.2284, -0.00083}, {-5.4159, 0.9861}}))

and the eigenvalues are

(({{-1.00031, -0.02115}}))

These values show that the average of symbolic solutions of MCARE and the numeric solution of the μ-Toolbox of MATLAB are in very good agreement.

Comparing LQ and LQR control design

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

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

Considering γ = 100 (Subscript[γ, crit]), we repeat the computations for solving Riccati equation,

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

The solutions of the MCARE are

NsolR = Solve[eqs[Nγ], {P11, P12, P22}]

For the selected positive definite solutions, we get

solPD =  
      Select[{{P11, P12}, {P12, P22}} /. 
          NsolR, (Eigenvalues[#][[1]] > 
                0 ∧ Eigenvalues[#][[2]] > 0) &]

Checking these solutions, whether they satisfy the Riccati equation,

Map[
      MatrixNorm[
            RI /. Join[
                numericalValues, {P11 -> #[[1, 1]], 
                  P12 -> #[[1, 2]], 
                  P22 -> #[[2, 2]], γ -> Nγ}]]/
          MatrixNorm[#] &, solPD]

Then the gain matrices are

K = 
      Map[(Transpose[B] . # + (1/N)(γ^2) Transpose[NL] . #) &, 
        solPD]; Map[MatrixForm[#] &, K]

The eigenvalues of these gain matrices are

Map[Eigenvalues[#] &, K]

However, the eigenvalues are different, the norm of the eigenvalues vector {(Subscript[λ, 1]), (Subscript[λ, 2])}practically are the same

Map[VectorNorm[#] &, %]

Now, we consider the first solution

KLQR∞ = K[[1]]; MatrixForm[KLQR∞]

The eigenvalues of the closed loop

Eigenvalues[(A /. numericalValues) - 
        B . KLQR∞]

These results are very close to the KLQ  solution, as we have seen

KLQ // MatrixForm

and

Eigenvalues[(A - B . KLQ) /. numericalValues]

The solution of the μ-Toolbox of MATLAB is  [14]

(({{(10^(-5)), -8.2 (10^(-8))}, {-5.3037, 0.98591}}))

and the eigenvalues are

(({{-1.0001, -0.02115}}))

which are in very  good agreement with the results of the LQ and LQR ( γ = 100 (Subscript[γ, crit]) ) designs.

Nonlinear system performance without control

First, the dynamical performance of the nonlinear system is simulated in case of food intake as a jump type increase of glucose input h(t) , considering as disturbance

δ[t_] := 
      0.01 (UnitStep[t - 10] - UnitStep[t - 40])

Plot[δ[t], {t, 0, 300}, Frame -> True, 
        PlotLabel -> "Food (sugar) intake simulation", 
        FrameLabel -> {"time(sec)", "Intake rate μU/kg \
min"}, PlotStyle -> RGBColor[1, 0.5, 0], 
        Background -> RGBColor[0, 0, 0]];

Therefore the glucose inlet is

h[t_] := δ[t]

The insulin inlet is zero

i[t_] := 0;

The differential equations of the model in numerical form are

deq1NL = deq1 /. numericalValues

deq2NL = deq2 /. numericalValues

The numerical solution of this system is

solNL = 
      NDSolve[{deq1NL, deq2NL, X[0] == 0. , 
          Y[0] == Y0 /. numericalValues}, {X[t], Y[t]}, {t, 0, 
          300}]

The solution for the glucose concentration

XNLNC = X[t] /. solNL[[1]];

which can be plotted

pXNC = 
        Plot[(XNLNC /. t -> z), {z, 0, 300},  Frame -> True, 
          PlotLabel -> "Glucose concentration without control", 
          FrameLabel -> {"time(sec)", "Concentration g/100 ml \
"}, PlotStyle -> RGBColor[1, 0.5, 0], 
          Background -> RGBColor[0, 0, 0]];

The solution for the insulin concentration

YNLNC = Y[t] /. solNL[[1]];

pYNC = 
        Plot[(YNLNC /. t -> z), {z, 0, 300},  Frame -> True, 
          PlotLabel -> "Insulin concentration without control", 
          FrameLabel -> {"time(sec)", "Concentration \
μU/ml"}, PlotRange -> {3, 6}, 
          PlotStyle -> RGBColor[1, 0.5, 0], 
          Background -> RGBColor[0, 0, 0]];

Nonlinear system performance with LQ control

Now, we employ insulin injection as a control action based on LQ control. The insulin inlet computed as a feedback of the deviance of the actual state from the equilibrium

i[t_] := -KLQ[[2]] . {X[t], 
            Y[t] - Y0 /. numericalValues}

This insulin inlet term exists only in the second differential equation, therefore we reset only this equation

deq2NL = deq2 /. numericalValues

The numerical solution of this system

solNL = 
      NDSolve[{deq1NL, deq2NL, X[0] == 0. , 
          Y[0] == Y0 /. numericalValues}, {X[t], Y[t]}, {t, 0, 
          300}]

then the concentrations are

XNLLQ = X[t] /. solNL[[1]];

the glucose concentration

pXLQ = 
        Plot[(XNLLQ /. t -> z), {z, 0, 300},  Frame -> True, 
          PlotRange -> All, 
          PlotLabel -> "Glucose concentration with LQ control", 
          FrameLabel -> {"time(sec)", "Concentration g/100 ml"}, 
          PlotStyle -> RGBColor[0, 0, 1], 
          Background -> RGBColor[0, 0, 0]];

similarly, the insulin concentration

YNLLQ = Y[t] /. solNL[[1]];

pYLQ = 
        Plot[(YNLLQ /. t -> z), {z, 0, 300},  Frame -> True, 
          PlotLabel -> "Insulin concentration with LQ control", 
          FrameLabel -> {"time(sec)", "Concentration μU/kg min\
"}, PlotRange -> All, PlotStyle -> RGBColor[0, 0, 1], 
          Background -> RGBColor[0, 0, 0]];

Nonlinear system performance with LQR control

Now, let us compute the system performance for LQR control. In this case, we have insulin  as well as glucose injection action

the glucose inlet is

h[t_] := δ[t] - 
        KLQR[[1]] . {X[t], Y[t] - Y0 /. numericalValues}

the insulin inlet is

i[t_] := -KLQR[[2]] . {X[t], 
            Y[t] - Y0 /. numericalValues}

The model equations are

deq1NL = deq1 /. numericalValues

deq2NL = deq2 /. numericalValues

The numerical solution of this system is

solNL = 
      NDSolve[{deq1NL, deq2NL, X[0] == 0. , 
          Y[0] == Y0 /. numericalValues}, {X[t], Y[t]}, {t, 0, 
          300}]

The solution for the glucose concentration

XNLLQR = X[t] /. solNL[[1]];

pXLQR = 
        Plot[(XNLLQR /. t -> z), {z, 0, 300},  Frame -> True, 
          PlotRange -> All, 
          PlotLabel -> "Glucose concentration with LQR control", 
          FrameLabel -> {"time(sec)", "Concentration g/100 ml"}, 
          PlotStyle -> RGBColor[1, 0, 0], 
          Background -> RGBColor[0, 0, 0]];

and for the insulin concentration

YNLLQR = Y[t] /. solNL[[1]];

pYLQR = 
        Plot[(YNLLQR /. t -> z), {z, 0, 300},  Frame -> True, 
          PlotRange -> {5.4925, 5.494}, 
          PlotLabel -> "Insuline concentration with LQR control", 
          FrameLabel -> {"time(sec)", "Concentration μU/kg min\
"}, PlotRange -> All, PlotStyle -> RGBColor[1, 0, 0], 
          Background -> RGBColor[0, 0, 0]];

Comparing dynamical performance controlled by LQ and LQR controller

Plotting the results of the two control strategies, the higher performance of the LQR control is clear.

Show[{pXLQ, pXLQR}, 
        PlotLabel -> "Glucose concentration with LQ and LQR \
control"];

Show[{pYLQ, pYLQR}, 
        PlotLabel -> "Insulin concentration with LQ and LQR \
control"];

Fig. 8. and   Fig.9.   show that the control designed by LQR  for  γ = \(γ\_crit\ is\ better\ \ than\ that\ of\ LQ\ control\ \ or\ \ which\ is\ the\ same\) , the control designed by  LQR for γ  =  
        100 (Subscript[γ, crit]).

Conclusions

Figures 8. and 9. show that the LQR controller design  provides a better control performance than the LQ controller design.

Employing symbolic computation, it was possible to determine the critical value of the scaling parameter γ as function of the model parameter:

(Subscript[γ, crit]) = (((Subscript[p, 2]) (Sqrt[(Subscript[p, 3]^2) + (Subscript[p, 4]^2)]))/((Subscript[p, 1]) (Subscript[p, 3])))

\(It\ turned\ out\ that\ for\ γ\  > \ γ\_crit\ there\ are\ more\ than\ one\ positive\ definite\ solutions\ of\ the\ modified\ Ricatti\ equation\), see Table 1. and  they are in good agreement with the result of  μ-Toolbox of MATLAB, however this latest gives only one of solutions, (K^a) .  As an illustration, the  (Subscript[K, 11]) (γ) 
        is computed for different γ.

({{(γ/(Subscript[γ, crit])), 1.1, 1.2, 1.3, 1.4, 1.5, 2., 3., 5.}, {((Subscript[K, 11])^a), 0.2284, 0.1021, 0.0824, 0.0685, 0.0580, 0.0305, 0.0130, 0.0046}, {((Subscript[K, 11])^b), 0.3237, 0.3548, 0.3745, 0.3884, 0.3988, 0.4264, 0.4439, 0.4524}})

Plotting these values together with the singular solution (γ ≈(Subscript[γ, crit])), we get

data = {{1. , 0.225278, 0.231658}, {1.1, 0.2284, 0.3237}, {1.2, 
            0.1021, 0.3548}, {1.3, 0.0824, 0.3745}, {1.4, 0.0685, 
            0.3884}, {1.5, 0.058, 0.3988}, {2. , 0.0305, 0.4264}, {3. , 
            0.013, 0.4439}, {5. , 0.0046, 0.4524}, {10. , 0.0011, 
            0.4559}, {100. , 0. , 0.457}};

\(\(p1 = 
        ListPlot[Map[Drop[#, \(-1\)] &, data], Frame -> True, 
          PlotJoined -> True, 
          PlotLabel -> \*""Gaim Matrix Element (Subscript[K, 11])(γ)"", 
          FrameLabel -> {\*""(γ/(Subscript[γ, crit]))"", \*""(Subscript[K, 11])""}, PlotStyle -> RGBColor[1, 0.5, 0], 
          Background -> RGBColor[0, 0, 0], 
          DisplayFunction -> Identity];\)\)

\(\(p2 = 
        ListPlot[Map[Drop[#, {2, 2}] &, data], Frame -> True, 
          PlotJoined -> True, 
          PlotLabel -> \*""Gaim Matrix Element (Subscript[K, 11])(γ)"", 
          FrameLabel -> {\*""(γ/(Subscript[γ, crit]))"", \*""(Subscript[K, 11])""}, PlotStyle -> RGBColor[1, 0.5, 0], 
          Background -> RGBColor[0, 0, 0], 
          DisplayFunction -> Identity];\)\)

Show[{p1, p2, 
          Graphics[{PointSize[0.04], RGBColor[0, 0, 1], 
              Point[Drop[data[[1]], -1]]}]}, 
        DisplayFunction -> $DisplayFunction];

It was also revealed that, the norms of the vectors, which elements are the eigenvalues of the pairs of gain matrices are practically the same for the different values of the scaling parameter.

References

1. Sturis K.S., Polonsky E., Mosekilde van Cauter, "Computer model for mechanisms underlying ultradian oscillations of insulin and glucose", American Journal of Physiology. Vol. 260. pp. 439-445, 1991.

2. Juhasz Cs. Medical Application of Adaptive Control Supporting Insulin Therapy in case Of Diabetes mellitus, Ph.D. dissertation, electrical Eng. Program.,Budapest Univ. of technology and Economics, Budapest, hungary, 1993.

3.  Lehmann, E.D., Deutsch T., "Compartmental models for glycaemic prediction and decision support in clinical diabetes care: promise and reality", Elsevier, Computer Methods and Programs in Biomedicine, Vol. 56., pp. 193-204, 1998. 

4. Tolic I.M., Mosekilde E., Sturis J. "Modeling the Insulin Glucose Feedback System: Signification of Pulsatile Insulin Selection", Journal of Theoretical Biology, Vol. 3. No. 1., pp.361-375, 2000.

5. Benett D.L., Gourley S.A. "Asymptotic properties of a delay differential equation model for the interaction of glucose with plasma and interstitial insulin", Elsevier, Applied Mathematics and Computation, in press 2003.

6. Sano A. "Adaptive and optimal schemes for control of blood glucose levels", Biomedical Measurements, Informatics and Control, Vol.1.No.1. pp. 16-22, 1986.

7. Fischer M.E., Teo K.L. "Optimal Insulin Infusion Resulting from a Mathematical Model of Blood Glucose Dynamics", IEEE Trans. Biomed.Eng. Vol. 36. No.3. pp.479-486,1989.

8.  Candas B., Radziuk J. "An adaptive controller of glycemia removal", Annu. Int. Conf. of IEEE Engng. in Medicine and Biology Soc. Vol.13. No.5. pp.2285-2286,1991.

9.  Benyó Z., Paláncz B.,Juhász Cs., Várady P. "Design of Glucose Control via Symbolic Computation", Proc. of the 20th Annu. Int. Conf. of IEEE Engng. in Medicine and Biology Soc. Vol.20. No.6. pp.3116-3119,1998.

10. Benyó B., Benyó Z., Paláncz B., Kovács L.  and Szilágyi L. "A Fully symbolic Design and Modelling of Nonlinear Glucose Control with Control system professional Suite (CSPS) of Mathematica", Proc. of the world Congress on Medical Physics and Biomedical Engng.  Available: http//library.wolfram.com/infocenter/Mathsource/5043/

11. Kovács L., Paláncz B. "Linear and Nonlinear Approach of the Glucose-Insulin Control using Mathematica", Periodica Politechnica Trans. On Automatic Control and Computer Science,Vol. 49 (63). No. 2. pp.65-71, 2004.

12. Juhász Cs., Asztalos B. "AdASDiM: An Adaptive Control Approach to Diabetic Management", Innovation et technologie en Bilologie et Medicine, Vol.'7.No.1, 1996.

13. Zhou K. Robust and Optimal Control, Prentice Hall, New Jersey, pp. 376-412, 1996.

14. Kovács L., Paláncz B., Almássy Zs. and Benyó Z. "Optimal Glucose-Insulin Control in (Subscript[H, 2])space", 26th Annual Int. Conf. of IEEE Engng. in Biomedicine Soc., San Francisco, (accepted for Session  1.5.5,  ID 554), 2004.