Fully Symbolic Design of Glucose Control  with CSP2 of Mathematica

Benyó (SuperStar[B])., Benyó, Z.,Paláncz, B. Kovacs,L. and Szil ágyi, L.  
Szechenyi Istvan Universit(SuperStar[y]), Hungary 
Budapest Technical University,Hungary
e-mail:   palancz@epito.bme.hu
Budapest, 2003

This notebook is an electronic version of the lecture delivered at the   World Congress on Medical and Biomedical Engineering  24-29  August 2003,  in Sydney, Australia

Abstract

In this case study a fully symbolic design and modelling method are  presented for blood glucose control  of diabetic patiens under intensive  care. The analysis is based on a modified two compartment model proposed by  BERGMAN et al. (1). The applied feedback control law decoupling even the  nonlinear model leads to a fully symbolic solution of the closed loop  equations. The effectivity of the applied symbolic procedures being mostly  built-in the new version of Control System Professional Suite (CSPS) Application of  Mathematica - have been demonstrated for controller design in case of a  glucose control for treament of diabetes mellitus. The results are in good  agreement with the earlier presented symbolic-numeric analysis by BENYO et  al. (7).
Keywords: symbolic computation, diabetes mellitus, blood glucose control  and modeling with Mathematica

Introduction

From the 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 the human body. To maintain the glucose level  in a diabetic under intensive care is currently an actively researched topic  in the field of Biomedical Engineering. Many different models and strategies  have been designed and applied to the problem SANO (6), FISCHER (3), CANDAS  (2), JUHASZ (4) and BENYO et al. (7). The authors orientated to BENYO (7),  considering as the best appropiate model.
In this paper symbolic computation  was used to design multivariable modal control based on the space  representation of a verified nonlinear model. The computations were carried  out with Mathematica Version 4.2, and the article was written in Mathematica and presented as a live worksheet

Model equations

To simulate the insulin-glucouse interaction in human body the following 
two-compartment model was employed:

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

and

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

The terms h(t) and i(t) as exogenous insulin and glucose take the impacts on 
glucose level into consideration, X(t) and Y(t) stand for the concentration 
of glucose in the plasma and that of the insulin remote from plasma, 
respectively.

Nonlinear model

In order to handle the problem with CSPS Application the following general form of a nonlinear model should be  considered:

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

where x(t) is the state variable, u(t) and y(t) are the input and output variables, respectively. In our case the casting is the following. There are two state variables

{(Subscript[x, 1]), (Subscript[x, 2])} = {X[t], Y[t]};

The two output variables supposed to be the same as the state variables 
provided the dynamical performances of the measurement and actuator devices 
are considerably faster than that of the system itself.

{(Subscript[y, 1]), (Subscript[y, 2])} = {X[t], Y[t]};

The model has two input variables:

{(Subscript[u, 1]), (Subscript[u, 2])} = {h[t], i[t]};

The function f and g can be defined as:

{(Subscript[f, 1]), (Subscript[f, 2])} = {deq1[[1]], deq2[[1]]}

{(Subscript[g, 1]), (Subscript[g, 2])} = {(Subscript[y, 1]), (Subscript[y, 2])};

To design the control for the system , it should be linearized. Now, we  load the CSP2 Application:

<< ControlSystems`

Linearization

The linearization should be carried out at a steady state, namely at (X0,  (Subscript[Y0, , ]) h0,  i0):

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

From the state representation, one can get the equation form of the 
linearized system:

ControlObjectSS // EquationForm

or its representation on the frequency domain:

ControlObjectTF = 
      TransferFunction[s, ControlObjectSS] // Simplify

Controllability test

First, we test the controllability of the linearized system.The controllability matrix is

ControllabilityMatrix[ControlObjectSS]

In case of controllable system the rank of this matrix must be equal with 
that of the system matrix

Rank[%] == Rank[ControlObjectSS[[1]]]

Therefore the system controllable. This can be checked by the following 
function, too:

Controllable[ControlObjectSS]

Steady state values

The linearization should be carried out at steady state, X[0] = 0, h[0] = 0,  i[0] = 0 and Y[0] = Y0.

X0 = X[0] = 0; h0 = h[0] = 0; i0 = i[0] = 0; 
    Y0 = Y[0];

We need to computeY0.The first model equation is identity,the second one is

eq = deq2[[1]] /. t -> 0

It can be easily solved for Y[0]

sol = Solve[eq == 0, Y[0]]

Y[0] = Y[0] /. sol[[1]]

Now, the linearized model is

ControlObjectSS // EquationForm

Introducing A and B for the matrices as usual, the equillibrum is stable, because the model  parameters p1 and p3 are negative.

A = ControlObjectSS[[1]]; MatrixForm[A]

B = ControlObjectSS[[2]]; MatrixForm[B]

Control Design

Let us consider the λ1, λ2 as the eigenvalues of the matrix A-KB, where K is the gain matrix. Then K can be computed as:

K = 
      Inverse[B] . (A - DiagonalMatrix[{λ1, λ2}]) // 
        Simplify; MatrixForm[K]

We can check this result by computing the eigenvalues:

Eigenvalues[A - K . B] // Simplify

This is a feasible control, because both of the model parameters ((Subscript[p, 1])and (Subscript[p, 3])) are negative.

Closed loop system

The state space form of our closed loop model is:

CLControlObjectSS = 
        StateFeedbackConnect[ControlObjectSS, K];

The equation form of this closed loop model is

CLControlObjectSS // EquationForm

It is interesting to see, that the closed loop system is decoupled.

Simulation of the linear closed loop model

As disturbance, let the initial conditions Xs = X - X0 > 0 and Ys - Y0  > 0, 
then we can compute the system responses even in symbolic form:

OutputResponse[CLControlObjectSS, {0, 0}, t, 
        InitialConditions -> {Xs, Ys}] // Simplify

It can be seen, that greater absolute value of ë1,  ë2, make the system reach the steady state faster. With other  words, the quality of the control will be improved with the increase of the  absolute value of the ë's, however in real cases, the dynamical  performance ability of the actuators can be the bottle-neck.
So, we were  able to get symbolic results for control design of our system and could draw  certain conclusions concerning the control performance, which demonstrate one  of the unique features of the Control System Professional Application! However this is not the end, but just the end of the beginning. Now, we  continue our study with further symbolic computations!

Nonlinear closed loop model

The variables of the linearized model represent the deviations from the 
steady state instead of the total values. Therefore to get the nonlinear 
closed loop model, one should take into considaration the steady values. So, 
the control vector is:

{h[t], i[t]} = {h0, i0} - K . {X[t] - X0, Y[t] - Y0}

Then our model equations are

deq1

deq2

The first equation is independent from the second one, therefore it can be 
solved as a single equation:

Clear[X, X0]

sol1 = DSolve[{deq1, X[0] == X0}, X[t], t]

Substituting this result into the second equation:

deq2m = deq2 /. sol1[[1]] // Simplify

the solution can be achieved in symbolic form:

Clear[Y, Y0]

sol2 = 
      DSolve[{deq2m, Y[0] == Y0}, Y[t], t] // FullSimplify

Then the time functions of the state variables of the closed loop system 
are

Xc[t_] = X[t] /. sol1[[1]]

Yc[t_] = (Y[t] /. sol2[[1]]) // Simplify

The control variables:

h[t_] = (h[t] /. X[t] -> Xc[t]) // Simplify

i[t_] = (i[t] /. {X[t] -> Xc[t], Y[t] -> Yc[t]}) // 
        Simplify

Now, let us finish our analysis with a numerical example for the nonlinear 
closed loop control simulation.

Numerical simulation of the nonlinear closed loop model

Let us consider the following numerical values JUHASZ et.al. (1996):

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

We consider a hypoglykaemic episode with initial values X0 = 0.1 and Y0 = 2.

X0 = 0.1;  Y0 = 2. ;

Let the eigenvalues be λ1=0.7(Subscript[p, 1])and λ2=1.05(Subscript[p, 3]).

The dynamic performance of the  state variables with control:

g1 = 
        Plot[Xc[
                t] /. {λ1 -> 0.7 (Subscript[p, 1]), λ2 -> 
                  1.05 (Subscript[p, 3])} /. numericalValues, {t, 0, 300}, 
          PlotStyle -> RGBColor[1, 0, 0],  
          AxesLabel -> {"t", "X[t]"}];

Fig.1 The first state variable in time with control.

g2 = 
        Plot[Yc[
                t] /. {λ1 ->  0.7 (Subscript[p, 1]), λ2 ->  
                  1.05 (Subscript[p, 3])} /. numericalValues, {t, 0, 300}, 
          PlotStyle -> RGBColor[1, 0, 0],  
          AxesLabel -> {"t", "Y[t]"}, 
          PlotRange -> All];

Fig.2 The second state variable in time with control.

Comparing controlled and not controlled dynamics

In order to qualify the effect of the control on the dynamical  performance, we carry out the simulation without control. This can be done by  only numerically. The state equations of the system without control (h(t) ≡ 0 and i(t) ≡ 0) are

deq1 = (Subscript[p, 1]) X[t] == X'[t]

and

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

The numerical solution of this system { XnL(t), YnL(t) } is

sol = 
        NDSolve[{deq1 /. numericalValues, deq2 /. numericalValues, 
            X[0] == X0, Y[0] == Y0}, {X[t], Y[t]}, {t, 0, 
            300}];

{XnL[t_], YnL[t_]} = {X[t], Y[t]} /. sol[[1]];

g3 = 
        Plot[XnL[t], {t, 0, 300}, PlotStyle -> RGBColor[0, 0, 1],  
          AxesLabel -> {"t", "X[t]"}, 
          DisplayFunction -> Identity];

g4 = 
        Plot[YnL[t], {t, 0, 300}, PlotStyle -> RGBColor[0, 0, 1],  
          AxesLabel -> {"t", "Y[t]"}, 
          DisplayFunction -> Identity];

Show[{g1, g3}, 
        DisplayFunction -> $DisplayFunction];

Fig.3 The first state variable in time with and without control.

Show[{g2, g4}];

Fig.4 The second state variable in time with and without control.

The figures clearly show, that the control considerably improves the 
dynamical performance of the system, although in case of the first variable 
the response is somewhat better without control.

Conclusions

Although, the control design was based on the linearized model, the control 
could improve considerably the system response. The symbolic results shows, 
that it is possible to carry out the design algorithm fully symbolical way!

According to the first results, the system is expected - after the necessary 
further verifications - to provide a useful help to control of blood glucose 
level in diabetics under intensive care, and to the optimisation process of 
diabetic administration.

References

[1] Bergman, B.N.-Ider, Y.Z.-Bowden, C.R.-Cobelli C. (1979). Quantitive 
estimation of insulin sensitivity., American Journal of Physiology, Vol. 236, 
pp. 667-677.

[2] Candas, B.-Radziuk, J. (1991). An adaptive controller of glycemia based 
on a physiological model of insulin-dependent glucose removal. Ann.Int.Conf. 
of the IEEE Eng. In Medicine and Biology Soc., Vol. 13, No. 5, pp. 2285-2286.

[3] Fischer, M.E.-Teo, K.L. (1989). Optimal Insulin Infusion Resulting from a 
Mathematical Model of Blood Glucose Dynamics. IEEE Trans. On BME, Vol. 36, 
No. 3, pp. 479-486.

[4] Juhász, Cs. - Asztalos, B. (1996). AdASDiM: An  

Adaptive Control Approach to Diabetic Management. Innovation et Technologie 
en Biologie et Medicine, Vol. 17, No. 1.

[5] Ogunye, A.B. (1996). Process Control and Symbolic Computation: An 
Overview with Maple V. MapleTech, Vol. 3, No. 1, pp. 94-103.

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

[7] Benyó, Z., Paláncz, B., Juhász, Cs., 
Várady, P. (1998). Design of Glucose Control via Symbolic 
Computation., Proceedings of the 20th Annual International Conference of the 
IEEE Engineering in Medicine and Biology Society, Vol. 20, No. 6, pp. 
3116-3119.