Buletinul Stiintific al Universitatii "Politehnica" din Timisoara, romania

I. INTRODUCTION

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 [4], [13], the authors orientated on [6] and [7], considering as the best appropriate model.

Symbolic computation was used to design and compare linear and nonlinear multivariable model control based on the space representation of a verified nonlinear model. In this way it was able to design a control law which can be evaluated rapidly, therefore it can be applied in on-line (real time) applications. Computations were carried out and the article was written in Mathematica Version 5 and presented as a live worksheet.
II. Model equations
To simulate the insulin-glucose interaction in human body the following two-compartment model was employed [6], [7]:

deq1 = p₁ X[t] + p₂ h[t] == X’[t]; (1)

deq2 = (p₃-X[t])Y[t]+ i[t]+ p₄ == 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.

1. 
   Nonlinear model
   

In order to handle the problem with Control System Professional Suite (CSPS) Application of Mathematica, the following general form of a nonlinear model, where x(t) is the state variable, u(t) and y(t) are the input and output variables, should be considered:

In our case the casting is the following. There are two state variables:

{x_1,, x₂} = {X[t], Y[t]} . (3)

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.

{y₁, y₂} = {X[t], Y[t]} . (4)

The model has two input variables:

{u₁, u₂} = {h[t], i[t]} . (5)

Therefore the function f and g can be defined as it follows:

{f₁, f₂} = {deq1[[1]], deq2[[1]]} =

= {h[t] p₂+ p₁ X[t], i[t]+ p₄ + (p₃-X[t])Y[t]}

{g₁, g₂} = {y₁, y₂} . (6)

To design the control for this system, the nonlinear system should be linearized, while numerical solution methods are available for nonlinear systems. Now, the CSPS Application can be loaded:

2. 
   Linearization
   

The linearization should be carried out at a steady state, namely at (X0, Y0, h0, i0), using Taylor approximation:

{{X[t], X0}, {Y[t], Y0}}, {{h[t], h0},

{i[t], i0}}] . (7)

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

or its representation on frequency domain:

TransferFunction[s,{{, 0},

{, }}] . (9)

3. 
   Controllability Test
   

In order to check whether if the system can be stabilized at all, first we test the controllability of the linearized system. The controllability matrix is:

M_C ={{p₂, 0, p₁ p₂, 0}, {0,1, -Y0 p₂, -X0 + p₃}} . (10)

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

Rank[M_C] = Rank[ControlObjectSS[[1]]] (11)

True

True

4. 
   Steady state values
   

As it was mentioned before, the linearization should be carried out at the steady state, namely at X[0]=0, h[0]=0, i[0]=0 and Y[0]=Y0.

We need to compute Y0. The first model equation is identity ( 0 ≡ 0 ), but the second one is:

eq = deq2[[1]] /. t → 0

p₄ + p₃ Y[0] = 0 . (13)

Therefore, it can be easily solved for Y[0]:

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

{{Y[0]}} . (14)

Now, the linearized model is:

ControlObjectSS // EquationForm

. (15)

Introducing A and B for the matrices as usual, the equilibrium is stable, if the poles of the system are on the left half plane ( p₁ and p₃ < 0).

5. 
   Control Design
   

While the linearized system is stable, the control design can be carried out. Let us consider the λ₁, λ₂, as the eigenvalues of the matrix A-KB (these are exactly the poles of the closed loop system too), where K is the gain matrix, what we want to design. Then K can be computed as:

// Simplify; MatrixForm[K]

We can check this result by computing the eigenvalues:

{λ1, λ2} . (17)

It is a special case of the control, when λ₁=p₂ and λ₂=p₃. This is a feasible control, because both of the model parameters are negative, JUHÁSZ (8).

6. 
   Closed Loop System
   

Now, we can implement the designed control in our model, resulting a closed loop. The state space form of our closed loop model is:

The equation form of this closed loop model is:

It is interesting to see, that the closed loop system is decoupled, while each output is depended only on the corresponding input.

7. 
   Simulation of the Linear Closed Loop Model
   

As disturbance, let us perturbate the initial conditions ∆X = X - X0 > 0 and ∆Y = Ys - Y0 > 0. Then we can compute the system responses even in symbolic form:

// Simplify (20)

{ } .

It can be seen, that greater absolute value of λ₁, λ₂, make the system reach the steady state faster. So 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 CSPS. Now, we continue our study with further symbolic computations.

8. 
   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 consideration the steady values. So, the control vector is:

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

{, } . (21)

Then our model equations are:

p₁X[t] - (-λ1 + p₁)X[t] = X’[t] , (22)

p₄+(p₃-X[t])Y[t]==

== Y’[t] . (23)

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

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

{{X[t] → e^tλlX0}} . (24)

Substituting this result into the second equation:

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

, (25)

the solution can be achieved in symbolic form:

{{Y[t] → }} . (26)

Then:

Xc[t_] = X[t] /. sol1[[1]] . (27)

X0

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

The control variables:

To finish the analysis, the following numerical values [6] were taken into consideration, demonstrating the designed nonlinear closed loop control:

We consider a hypoglykaemic episode with initial values X0=0.1 and Y0=2, so the eigenvalues were: ; .

The simulation results can be seen in section IV (Fig.1 – Fig.4), together with the simulation results of the nonlinear control using Pontryagin's maximum principle, which is discussed in the following.

III. NON-Linear Approach of the Glucose-Insulin Control

1. 
   Model Reduction
   

Assuming that we do not use glucose injection (h(t)  0) - which is a very reasonable assumption, because the reference value of the glucose concentration is zero - and the only control variable is the insulin inlet, the first equation of our model can be reduced as:

X’[t] == p₁ X[t]

This can be solved with initial condition X0:

X[t_] = X[t] /. DSolve[{deq1, X[0] == X0},

X[t], t][[1]] , (33)

X0 .

This state variable X[t] is exponentially decreasing towards the steady state, because in our model the parameter p₁ < 0. Then substituting this result into the second equation we get a single linear, but nonautonom equation as system equation.

Y′[t]== i[t] + p₄ + (-X0 + p₃) Y[t] .

2. 
   Objective Function
   

The kernel of our objective function represents the deviation of the actual state from the steady state and the consumed insulin of the process:

where α₁ and α₂ are weights.

Then the objective function to be minimized is:

(36)

.

where  is the termination time of the control.

3. 
   The Pontryagin's Maximum Principle (PMP)
   

Let us suppose we have a y functional to be maximized or minimized according to the control function u(t):

under the constrain of the following differential equations system:

with the initial condition X(0) = X0.

To find the optimal control u(t), the following Hamiltonian can be defined employing Lagrangian multipliers:

. (39)

The set of necessary condition for the existence of an optimum, the adjoint functions (Lagrangian multipliers) can then be written from the Euler -Lagrange equations for the modified objective function as:

with the boundary conditions:

Moreover, the u(t) must be selected, so that H is a maximum (or minimum as required), which means:

and in case of maximum (or minimum):

4. 
   Application of PMP
   

The Hamiltonian of the problem involving our system equation as a constrain in form of Lagrangian multiplier is:

+(-X0 + p₃) Y[t])λ₁[t] .

2 i[t]β₁ + λ₁[t] == 0 .

This equation ensures that the stationary point is minimum, if β₁ > 0. Then the control function can expressed as:

i[t] = i[t] /. Solve[eq, i[t]][[1]] (46)

- .

5. 
   Solution of Split Boundary Problem
   

The steady state of the system can be easily computed. From the first equation we get X_st = 0, and from the second one, assuming that at steady state i = 0, .

The adjoint function (Lagrangian multiplier) considering that the steady state for Y is , becomes:

deq3 = λ₁′[t] ==- λ₁[t] D[deq2[[2]],

Y[t]]-D[Ω, Y[t]] /. , (47)

λ₁′[t] ==-(-X0 + p₃) λ₁[t]) .

The boundary conditions of the system are:

bc23 = {λ₁[θ] == 0, Y[t]== Y0} . (48)

This is a nonautonom, but linear system, with split boundary values, because deq2 after substituting i(t) is:

Y′[t] == p₄ + (-X0 + p₃) Y[t] - .

The numerical values were considered the same as in case of the linear approach of the glucose-insulin control, (31):

numericalValues = {p₁ → -0.021151,

p₂ → 0.092551, p₃ → -0.014188,

p₄ → 0.077947, X0 → 0.1, Y0 → 2.,

θ → 300., α ₁ → 0.09, β₁ → 10}; (50)

The system can be solved numerically:

{{Y[t] → InterpolatingFunction[{{0.,300.}},

However more solutions may exist which satisfy this split boundary problem:

InterpolatingFunction[{{0., 300.}}, <>][t]

InterpolationFunction[{{0., 300.}}, <>][t]

The boundary conditions are satisfied fairly well:

{Ys[0], λ_1s[300]} (54)

{2.01953, 0.} .

The steady state {0, -} is nearly reached:

{X[300]/. numericalValues, Ys[300]} (55)

{0.000175498, 5.48098} ,

because

-/. numericalValues (56)

5.49387 .

IV. SIMULATION RESULTS
Using the same numerical values ((31), (50)) in both control cases, the simulations point out the differences between the two control strategies: linear and nonlinear approach of the glucose-insulin control

In case of nonlinear approach using Pontryagin's maximum principle, the glucose concentration is stabilized faster (Fig.3), while the insulin concentration reaches the steady state value slower (Fig.4). Due to the fact, that the priority of the control is the stabilization of the glucose concentration, we can tell that the nonlinear approach using Pontryagin’s maximum principle gives better results.

This conlusion is also demonstrated by the evolution of the input signals. In case of nonlinear approach, the variation of the exogenous insulin input is smoother than in case of linear approach (Fig.1).

However, it must be taken into consideration that the glucose injection was considered zero using Pontryagin’s maximum principle, but this is not a big constraint, while phisiologically speaking we inject only insulin. This is also demonstrated by the very small variation of the glucose input in case of linear approach (Fig.2).

V. Conclusions
The main advantage of the selected model, in comparison with the other mentioned models is that it is an on-line adaptive model, which operates on strong control theoretical foundations as well as describes the physiological system properly. In the literature there are some articles dealing with more sophisticated models of glucose-insulin interaction, but they have not been applied for solving control problems [12], [14]. Until now, as we know, the existing adaptive models worked in a flip-flop way, selecting the control command between two values, but it goes without saying they have neighter physiological nor control theoretical background.

Using the presented control model, now the control command can be selected from an interval, changing its value continuously, in every moment. Moreover the model is not complicated (it uses only two differential equations), so it gives a big advantage in case of practical implementation.

Results show, that it is possible to carry out the design algorithm in a fully symbolical way using the CSPS of Mathematica.

The nonlinear optimal control employing only one control variable proved to be more effective than the control strategy using both glucose and insulin injection as control variables, based on the linearized model [6], [7]. However, the existence of multiply solution may jeopardize its practical application. We plan to investigate the control of this process employing more complex dynamical models in the future.

However, the model is not implemented yet in a real (practical) application, but after the necessary further verifications could provide a useful help to control of blood glucose level in diabetics under intensive care, and to the optimization process of diabetic administration.

ACKNOWLEDGMENTS
This research has been supported by Hungarian National Research Fund, Grants No. OTKA T029830, T042990 and by Hungarian Ministry of Education Grant No. FKFP 200/2001.

REFERENCES

1. 
   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.
   
2. 
   B. Candas, J. Radziuk, “An adaptive controller of glycemia based on a physiological model of insulin-dependent glucose removal”, Annual International Conference of the IEEE Engineering in Medicine and Biology Society, vol. 13, no. 5, pp. 2285 - 2286, 1991.
   
3. 
   M.E. Fischer, K.L. Teo, “Optimal Insulin Infusion Resulting from a Mathematical Model of Blood Glucose Dynamics”, IEEE Transactions on Biomedical Engineering, vol. 36, no. 3, pp. 479 - 486, 1989.
   

4. 
   Cs. Juhász, B. Asztalos, “AdASDiM: An Adaptive Control Approach to Diabetic Management”, Innovation et Technologie en Biologie et Medicine, vol. 17, no. 1, 1996.
   
5. 
   A. Sano, “Adaptive and optimal schemes for control of blood glucose levels”, Biomedical Measurements, Informatics and Control, vol. 1, No. 1, pp. 16 – 22, 1986.
   
6. 
   Z. Benyó, B. Paláncz, Cs. Juhász, P. Várady, “Design of Glucose Control via Symbolic Computation”, Proceedings of the 20^th Annual International Conference of the IEEE Engineering in Medicine and Biology Society, vol. 20, no. 6, pp. 3116 - 3119, 1998.
   
7. 
   B. Benyó, Z. Benyó, B. Paláncz, L. Kovács, and L. Szilágyi, “A Fully Symbolic Design and Modelling of Nonlinear Glucose Control with Control System Professional Suite (CSPS) of Mathematica”, Proceedings of the World Congress on Medical Physics and Biomedical Engineering, Sydney, Australia, Mathematica notebook in e-publication form Mathsource 5043, Wolfram Research, no. 34, 2003.
   
8. 
   J. Sturis, K.S. Polonsky, E. Mosekilde, E. van Cauter, “Computer model for mechanisms underlying ultradian oscillations of insulin and glucose”, Am. J. Physiol., vol. 260, pp. 439 - 445, 1991.
   
9. 
   E.D. Lehmann, T. Deutsch, “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.
   
10. 
    I.M. Tolic, E. Mosekilde, J. Sturis, “Modeling the Insulin-Glucose Feedback System: The Signification of Pulsatile Insulin Selection”, Journal of Theoretical Biology, vol. 207-3, pp. 361 - 375, 2000.
    
11. 
    Ogunye, A.B. “Process Control and Symbolic Computation: An Overview with Maple V”, MapleTech, vol. 3, no. 1, pp. 94 - 103, 1996.
    
12. 
    D.L. Benett, S.A. Gourley, “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.
    
13. 
    Cs. Juhász, Medical Application of Adaptive Control, Supporting Insulin-Therapy in case of Diabetes Mellitus, PhD thesis, Budapest, 1993.
    
14. 
    A Mukhopadhyay, et al, “Modelling the Intra-Venous Glucose Tolerance Test: a global study for a single-distributed-delay model” American Institute of Mathematical Sciences, Discrete and Continuous Dynamical Systems–Series B, vol. 4, no. 2, pp. 407 - 417, 2004.