H1111 Budapest, Műegyetem rkp 3 HUNGARY, Phone: (361) 4632223, Fax: (361) 4633084, EMail: palancz@epito.bme.hu
Abstract – In this case study two symbolicnumeric design and modeling methods are compared for blood glucose control of diabetic patients under intensive care. The analysis is based on a modified two compartment model proposed by [1]. In the case of linear approach, the applied feedback control law decoupling even the nonlinear model leads to a fully symbolic solution of the closed loop equations. In case of nonlinear approach, the nonlinear optimal control law based on the Pontryagin's maximum principle has been developed assuming that glucose injection is the only control variable. The effectivity of the applied symbolic  numeric procedures being mostly builtin the new version of Control System Professional Suite (CSPS) of Mathematica 5, have been demonstrated for controller design in case of a glucose control for treatment of diabetes mellitus. The results are in good agreement with the earlier presented symbolicnumeric analysis by [6], [7].
Keywords: symbolic computation, diabetes mellitus, blood glucose control and modeling with Mathematica, nonlinear control, Pontryagin's maximum principle.
I. 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 bloodglucosecontrol 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. Towards the last 50 years, a variety of models for the interaction between glucose and insulin have been suggested in the literature, [8], [9], [10], [12] and strategies have been designed and applied to the problem [2], [3], [4], [5], [6] and [11].
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 online (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 insulinglucose interaction in human body the following twocompartment model was employed [6], [7]:
deq1_{ }=_{ }p_{1}_{ }X[t]_{ }+_{ }p_{2 }h[t]_{ }==_{ }X’[t]; (1)
deq2_{ }=_{ }(p_{3}X[t])Y[t]+_{ }i[t]+_{ }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.
III. Linear Approach of the GlucoseInsulin Control

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:
=_{ }f(x,u) (2)
y_{ }=_{ }g(x,u) .
In our case the casting is the following. There are two state variables:
{x_{1,}, x_{2}} = {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_{1}, y_{2}} =_{ }{X[t], Y[t]} . (4)
The model has two input variables:
{u_{1}, u_{2}} = {h[t], i[t]} . (5)
Therefore the function f and g can be defined as it follows:
{f_{1}, f_{2}} =_{ }{deq1[[1]], deq2[[1]]} =
= {h[t]_{ }p_{2}+ p_{1}_{ }X[t],_{ }i[t]+_{ }p_{4}_{ }+_{ }(p_{3}X[t])Y[t]}
{g_{1}, g_{2}}_{ }=_{ }{y_{1}, y_{2}} . _{ } (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:
<< ControlSystem` .

Linearization
The linearization should be carried out at a steady state, namely at (X0, Y0, h0, i0), using Taylor approximation:
ControlObjectSS =_{ }Linearize[{f_{1}, f_{2}}, {g_{1}, g_{2}},
{{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:
ControlObjectSS // EquationForm
, (8)
or its representation on frequency domain:
ControlObjectTF_{ }=_{ }TransferFunction[s,
ControlObjectSS]_{ }//_{ }Simplify
TransferFunction[s,{{, 0},
{, }}] . (9)

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:
ControllabilityMatrix[ControlObjectSS]
M_{C }={{p_{2},_{ }0,_{ }p_{1 }p_{2},_{ }0},_{ }{0,1,_{ }Y0_{ }p_{2},_{ }X0_{ }+_{ }p_{3}}} . (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
Therefore, our system is controllable, consequently the system can be stabilized. We get the same result by using the following builtin function:
Controllable[ControlObjectSS] (12)
True

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_{4} + p_{3 }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_{1} and p_{3} < 0).

Control Design
While the linearized system is stable, the control design can be carried out. Let us consider the λ_{1}, λ_{2}, as the eigenvalues of the matrix AKB (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:
K_{ }=_{ }Inverse[B].(ADiagonalMatrix[{λ1,λ2}])
//_{ }Simplify;_{ }MatrixForm[K]
. (16)
We can check this result by computing the eigenvalues:
Eigenvalues[AK.B] // Simplify
{λ1, λ2} . (17)
It is a special case of the control, when λ_{1}=p_{2} and λ_{2}=p_{3}. This is a feasible control, because both of the model parameters are negative, JUHÁSZ (8).

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:
CLControlObjectSS =
StateFeedbackConnect[ControlObjectSS, K] . (18)
The equation form of this closed loop model is:
CLControlObjectSS // EquationForm
. (19)
It is interesting to see, that the closed loop system is decoupled, while each output is depended only on the corresponding input.

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:
OutputResponse[CLControlObjectSS,_{ }{0,_{ }0},
t, InitialConitions_{ }→_{ }{ΔX,_{ }ΔY}]
//_{ }Simplify (20)
{_{ } } .
It can be seen, that greater absolute value of λ_{1}, λ_{2}, 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 bottleneck.
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.

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_{1}X[t]_{ }_{ }(λ1_{ }+_{ }p_{1})X[t]_{ }=_{ }X’[t] , (22)
p_{4}_{}+(p_{3}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λl}X0}} . (24)
Substituting this result into the second equation:
deq2m_{ }=_{ }deq2_{ }/. sol1[[1]]_{ }//_{ }Simplify
, (25)
the solution can be achieved in symbolic form:
sol2_{ }=_{ }DSolve[{deq2m,Y[0]=Y0},Y[t],t]
//FullSimplify
{{Y[t] →_{ }_{ } }} . (26)
Then:
Xc[t_]_{ }=_{ }X[t]_{ }/._{ }sol1[[1]] . (27)
X0
Yc[t_]_{ }=_{ }(Y[t]_{ }/._{ }sol2[[1]]) // Simplify . (28)
_{}
The control variables:
h[t_] = (h[t]/.X[t]_{ }→_{ }Xc[t]) // Simplify
. (29)
i[t_]=(i[t]/.{X[t]_{ }→_{ }Xc[t],_{ }Y[t]_{ }→_{ }Yc[t]) //Simplify
_{} (30)
To finish the analysis, the following numerical values [6] were taken into consideration, demonstrating the designed nonlinear closed loop control:
numericalValues_{ }=_{ }{p_{1} →_{ }0.021151,_{ }p_{2 }→_{ }0.092551,
_{ }p_{3 }→_{ }0.014188,_{ }p_{4 }→_{ }0.077947} . (31)
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. NONLinear Approach of the GlucoseInsulin Control

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:
deq1_{ }=_{ }X’[t]_{ }==_{ }p_{1}_{ }X[t] . (32)
X’[t]_{ }==_{ }p_{1 }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_{1 }< 0. Then substituting this result into the second equation we get a single linear, but nonautonom equation as system equation.
deq2_{ }= Y′[t]==_{ }(p_{3}_{ }_{ }X[t])Y[t]_{ }+ i[t]_{ }+ p_{4} (34)
Y′[t]==_{ }i[t]_{ }+ p_{4} + (X0 + p_{3}) Y[t] .

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:
(35)
.
where α_{1} and α_{2} are weights.
Then the objective function to be minimized is:
(36)
.
where is the termination time of the control.

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):
, (37)
under the constrain of the following differential equations system:
, (38)
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:
, (40)
with the boundary conditions:
, for i = 1,2,…,n. (41)
Moreover, the u(t) must be selected, so that H is a maximum (or minimum as required), which means:
, (42)
and in case of maximum (or minimum):
(or ) , for . (43)

Application of PMP
The Hamiltonian of the problem involving our system equation as a constrain in form of Lagrangian multiplier is:
(44)
+
+(X0 + p_{3})_{ }Y[t])λ_{1}[t] .
Pontryagin's maximum principle states that we should select i(t) control function, which will in our case minimize the Hamiltonian:
eq_{ }= D[H,_{ }i[t]]== 0 (45)
2 i[t]β_{1 }+_{ }λ_{1}[t]_{ }==_{ }0 .
This equation ensures that the stationary point is minimum, if β_{1 }>_{ }0. Then the control function can expressed as:
i[t]_{ }=_{ }i[t] /._{ }Solve[eq,_{ }i[t]][[1]] (46)
 .

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 = λ_{1}′[t]_{ }== λ_{1}[t]_{ }D[deq2[[2]],
_{ }Y[t]]D[Ω,_{ }Y[t]] /. , (47)
λ_{1}′[t]_{ }==(X0 + p_{3})_{ }λ_{1}[t]) .
The boundary conditions of the system are:
bc23 = {λ_{1}[θ]_{ }==_{ }0,_{ }Y[t]==_{ }Y0} . (48)
This is a nonautonom, but linear system, with split boundary values, because deq2 after substituting i(t) is:
deq2 (49)
Y′[t]_{ }== p_{4 }+_{ }(X0 + p_{3})_{ }Y[t]_{ } .
The numerical values were considered the same as in case of the linear approach of the glucoseinsulin control, (31):
numericalValues =_{ }{p_{1} →_{ }0.021151,
p_{2 }→_{ }0.092551, _{ }p_{3 }→_{ }0.014188,
p_{4 }→_{ }0.077947, X0 →_{ }0.1, Y0 →_{ }2.,
θ →_{ }300., α_{ 1} →_{ }0.09, β_{1} →_{ }10}; (50)
The system can be solved numerically:
sol_{ }=_{ }NDSolve[Join[{deq2,_{ }deq3},_{ }bc23]/
/._{ }numericalValues,_{ }{Y[t], λ_{1}[t]},
{t,_{ }0,_{ }300},_{ }Method_{ }→”ImplicitRungeKutta”] (51)
{{Y[t] → InterpolatingFunction[{{0.,300.}},_{ }
<>][t], λ_{1}[t] → InterpolatingFunction [{{0.,_{ }300.}},_{ }<>][t]}} .
However more solutions may exist which satisfy this split boundary problem:
Ys[t_]=_{ }Y[t]/.sol[[1]] . (52)
InterpolatingFunction[{{0.,_{ }300.}},_{ }<>][t]
λ_{1s}[t_]=_{ }λ_{1}[t]/.sol[[1]] . (53)
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 glucoseinsulin 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 online 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 glucoseinsulin 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 flipflop 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.