Design of Luenberger observer for insulin-glocose control via Mathematica

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

Many articles dealing with insulin - glucose control have been published in the last decade, and all of them assumed that all the system state variables are available for feedback. However this is not usually the case, because mainly the blood plasma glucose measurable in the practise. In this study the use of the reduced-order estimator (also known as the Luenberger observer) is considered to reconstruct those state variables that cannot be recovered directly from the system outputs. An observer for state variables - remote compartment insulin and plasma insulin - is designed in case of extended minimal Bergman model, employing Polynomial Control System Application of Mathematica. The results of nonlinear closed loop simulations with the classical  LQ control showed that an observer, which about ten times faster than the system itself can provide a very good state recovery performance.

System Equations

Several different models of diabetic systems exist in the literature including, for example, the very detailed 21st-order metabolic model of  Sorensen (1985). 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 (1981),

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 (1985):

numPar = {p1 -> 0.028, p2 -> 0.025, p3 -> 0.00013, GB -> 110, 
          YB -> 1.5, p4 -> 5. /54, VL -> 120};

Steady State

In order to linearize the system, we need its steady-state values. The steady-state values of the three state variables are,

G0 = X0 = Y0 = 0;

The steady-state value of the control variable can be computed from the third differential equation,

Solve[i0/VL == p4 YB, i0] /. numPar

which means

i0 = i0 /. %[[1]]

The steady-state value of the disturbance variable is,

h0 = 0;

Now, we save the state of the actual context path, and then load the Control System Professional Suite (CSPS) application,

Linearization

<< ControlSystems`

The linearization can be carried out in the vicinity of this steady state, using new variables denoting the deviations from the equilibrium,

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

In our model the only measurable output is the G(t) - plasma glucose deviation,

{(Subscript[y, 1])} = {G[t]};

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

{(Subscript[f, 1]), (Subscript[f, 2]), (Subscript[f, 3])} = {deq1NL[[2]], deq2NL[[2]], 
          deq3NL[[2]]};

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

The linearized model, as a control object in CSPS,

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

Its state-space form is

ControlObjectSS // EquationForm

or giving the model matrices,

ss = ControlObjectSS // TraditionalForm

Steady-State Stability

It is useful to check the stability of the open-loop system. The system matrix can be extracted as

A = ss[[1, 1]]; MatrixForm[A]

and its eigenvalues are

Eigenvalues[A]

or

% /. numPar

which means the system is stable.

Controllability and Observability

The controllability and observability can be tested by using the corresponding built - in functions,

Controllable[ControlObjectSS]

Observable[ControlObjectSS]

State Reconstruction

In our model only the G(t) - plasma glucose deviation can be measured, the other two state variables, the X(t) - remote compartment insulin utilization and 	the Y(t) - plasma insulin deviation,  and their corresponding output variables, (Subscript[y, 2]) and (Subscript[y, 3]) should be reconstructed in order to provide feedback signal for the controller.
Here, the use of the reduced-order estimator (also known as the Luenberger observer) is considered, as this can be used to reconstruct those state variables that cannot be recovered directly from the system outputs. The reduced-order state estimator described in Patel and Munro (1982) is based on the Luenberger observable canonical form of the system,

Given an (n^"th")-order linear system in state-space form

(OverDot[x])= A x+B u
y= C x+D u

with q outputs and p inputs, then a reduced-order dynamical system of order n-q, called an estimator, can be constructed with state vector z such that z approximates a linear combination of n-q of the system states given by

z=Lx

where L has dimensions (n-q) × n. An observer of order n-q is then defined for the system above by

(OverDot[z])=Fz+Gu+Hy

where F, G, and H are matrices of constants with dimensions (n-q) × (n-q), (n-q) × p, and (n-q) × q, respectively. The desired estimate (OverHat[x]) of the system state vector x is then recovered from z and the system outputs y as

(OverHat[x])=(Subscript[K, 1])y+(Subscript[K, 2])z

The design of the observer requires the solution of the observer equation

LA-FL=HC

such that the matrix L has maximum rank, in order to guarantee the recovery of the unmeasurable state variables. Thus, not only must L have maximum rank, but the matrix [({{C}, {L}})] must also have maximum rank, since (Subscript[K, 1]) and (Subscript[K, 2]) are determined from

[(Subscript[K, 1]) (Subscript[K, 2])] =(([({{C}, {L}})])^(-1))

The reduced-order estimator algorithm developed by Munro (Patel and Munro (1982)) is based on the Luenberger observable canonical form of the system.
To determine a reduced-order estimator for a given system, you simply have to specify the dynamics of the observer, that is, the n-q eigenvalues of the matrix F. Here, these do not need to be distinct, as is the case with other reduced-order observer algorithms, but must be chosen to be faster than the fastest eigenvalue of the system matrix A. An advantage of this algorithm is that the designer can specify exactly the speed with which the observer tracks the system states. It can also deal with systems containing symbols.

Using the Polynomial Control Systems Application of Mathematica, the function ReducedOrderEstimator implements this approach and returns the state-space system

Fig.1 System with observer

where the observer,with inputs(({y,u})^T)and outputs(OverHat[x]),is connected to the system as shown in Fig .  
          1

Let us load the Application

<< PolynomialControlSystems`

Let us consider the n - q = 3 - 1 = 2 eigenvalues of the F matrix, {(Subscript[p, 1]), (Subscript[p, 2])}, then the reduced order estimator is

observer = 
      Simplify[ReducedOrderEstimator[
          ControlObjectSS, {(Subscript[p, 1]), (Subscript[p, 2])}]]

The linear form of the original system plus observer, see Fig. 1, can be constructed as he following:

The number of input is,

ins = CountInputs[ControlObjectSS]

Let us connects a dummy identity block I in parallel with the original system, to allow the system inputs to be connected to the observer,

ss1 = 
      ParallelConnect[TransferFunction[IdentityMatrix[ins]], ControlObjectSS, 
        Table[{i, i}, {i, ins}], {}]

Then let us connect the observer

ss2 = SeriesConnect[ss1, observer]

or in equation form

ss2 // EquationForm

Testing the observer via nonlinear simulation

Let us extract the observer matrices from the observer object, considering its structure

[({{F1, G1, H1}, {K2, 0, K1}})]

F1 = observer[[1]]; MatrixForm[F1]

<< LinearAlgebra`MatrixManipulation`

G1 = SubMatrix[observer[[2]], {1, 1}, {2, 2}]; 
    MatrixForm[G1]

H1 = SubMatrix[observer[[2]], {1, 3}, {2, 1}];; 
    MatrixForm[H1]

K2 = observer[[3]]; MatrixForm[K2]

K1 = SubMatrix[observer[[4]], {1, 3}, {3, 1}]; 
    MatrixForm[K1]

Then the differential equations for the observer for z(t)

deqz1 = z1'[t] == 
        First[F1[[1, 1]] z1[t] + F1[[1, 2]] z2[t] + 
            G1[[1, 1]] h[t] + G1[[1, 2]] h[t] i[t] + 
            H1 [[1]] G[t]]

deqz2 = z2'[t] == 
        First[F1[[2, 1]] z1[t] + F1[[2, 2]] z2[t] + 
            G1[[2, 1]] h[t] + G1[[2, 2]] h[t] i[t] + 
            H1 [[2]] G[t]]

The initial conditions for  z(t) can be computed from the initial condition of the states to be recovered, see equations for  (OverHat[x]).  

For z1(0) corresponding to X(0)

eqin1 = 
      X0 == 
        K1[[2]] G0 + K2[[2, 1]] z10 + 
          K2[[2, 2]] z20

and for z2(0) corresponding to Y(0)

eqin2 = 
      Y0 == K1[[3]] G0 + K2[[3, 1]] z10 + 
          K2[[3, 2]] z20

Solve[{eqin1, eqin2}, {z10, z20}]

The performance of the test a standard meal disturbance with about six our duration, modelled by Lehman and Deutsch (1992) is used. 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 nonlinear model equations using the numerical parameters, are

deq1 = deq1NL /. numPar

deq2 = deq2NL /. numPar

deq3 = deq3NL /. numPar

Let the observer roughly ten times faster the the system itself, then the eigenvalues for F matrix, {(Subscript[p, 1]), (Subscript[p, 2])}={-1, -1}

deq4 = deqz1 /. numPar /. {(Subscript[p, 1]) -> -1, 
          (Subscript[p, 2]) -> -1}

deq5 = deqz2 /. numPar /. {(Subscript[p, 1]) -> -1, 
          (Subscript[p, 2]) -> -1}

We employ a sinus signal for the insulin input on the top of its steady state

i[t_] := i0 + 5 Sin[ t/50];

Plot[i[t], {t, 0, 1000}, PlotStyle -> RGBColor[0, 0, 1], 
        Frame -> True, 
        FrameLabel -> {"t[min]", "i[mU/min]"}];

Fig.3  Insulin infusion

Now we solve this differential equation system consisting of the three states variables and the six controller states, simultaneously

sol = 
      NDSolve[{deq1, deq2, deq3, deq4, deq5, X[0] == 0, 
          Y [0] == 0, G[0] == 0, z1[0] == 0, 
          z2[0] == 0}, {X, Y, G, z1, z2}, {t, 0, 1000}, 
        Method -> BDF,  MaxSteps -> 50000]

The dynamics of the remote compartment insulin utilization of the nonlinear system

P1 = 
        Plot[Evaluate[(X[t] /. sol)], {t, 0, 1000}, 
          PlotStyle -> {RGBColor[0, 0, 1]}, Frame -> True, 
          PlotRange -> All, 
          FrameLabel -> {"t[min]", "X[1/min]"}, 
          DisplayFunction -> Identity];

The dynamics of the plasma insulin deviation of the nonlinear system

P2 = 
        Plot[Evaluate[(Y[t] /. sol)], {t, 0, 1000}, 
          PlotStyle -> {RGBColor[0, 0, 1]}, Frame -> True, 
          PlotRange -> All, 
          FrameLabel -> {"t[min]", "Y[mU/dL]"}, 
          DisplayFunction -> Identity];

The recovered states are

Xrec[
        t_] := K1[[2]] (G[t] /. sol) + 
            K2[[2, 1]] (z1[t] /. sol) + 
            K2[[2, 2]] (z2[t] /. sol) /. 
          numPar /. {(Subscript[p, 1]) -> -1, (Subscript[p, 2]) -> -1}

Yrec[
        t_] := K1[[3]] (G[t] /. sol) + 
            K2[[3, 1]] (z1[t] /. sol) + 
            K2[[3, 2]] (z2[t] /. sol) /. 
          numPar /. {(Subscript[p, 1]) -> -1, (Subscript[p, 2]) -> -1}

The recovered X(t) state

P3 = 
        Plot[Xrec[t], {t, 0, 1000}, PlotStyle -> {RGBColor[1, 0, 0]}, 
          Frame -> True, PlotRange -> All, 
          FrameLabel -> {"t[min]", "Xrec[1/min]"}, 
          DisplayFunction -> Identity];

The recovered Y(t) state

P4 = 
        Plot[Yrec[t], {t, 0, 1000}, PlotStyle -> {RGBColor[1, 0, 0]}, 
          Frame -> True, PlotRange -> All, 
          FrameLabel -> {"t[min]", "Yrec[mU/dL]"}, 
          DisplayFunction -> Identity];

We can compare the original (blue) and recovered (red) states

Show[{P1, P3}, Frame -> True, PlotRange -> All, 
        FrameLabel -> {"t[min]", "X and Xrec[1/min]"}, 
        DisplayFunction -> $DisplayFunction];

Fig.4 The original (blue) and recovered (red) states of remote compartment insulin utilization

and for plasma insulin

Show[{P2, P4}, Frame -> True, PlotRange -> All, 
        FrameLabel -> {"t[min]", "Y and Yrec[mU/dL]"}, 
        DisplayFunction -> $DisplayFunction];

Fig.5 The original (blue) and recovered (red) states of  plasma insulin deviation

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.

We choose the following matrices, R is

R =  DiagonalMatrix[{(10^4),  (10^(-4))}]; 
    Q =  (10^(-3)) IdentityMatrix[3];  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 /. numPar, Q, R]; 
    MatrixForm[KLQ] // N

The first row in K can be neglected

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

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

Then the  control law is

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

Nonlinear closed loop simulation for ideal case

Now we solve this differential equation system consisting of the three states variables supposing that all of the three states can be feeded back to the controller

sol = 
      NDSolve[{deq1, deq2, deq3, 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

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

Fig. 6 Controlled plasma glucose deviation in  case of ideal state feedback

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

Fig. 7 Controlled remote compartment insulin utilization in  case of ideal state feedback

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

Fig. 8  Controlled plasma insulin deviation in case of ideal state feedback

PC1 = 
        Plot[Evaluate[(i0 - 
                KLQ[[2]] . {G[t] /. sol, X[t] /. sol, 
                    Y[t] /. sol})], {t, 0, 1000}, 
          PlotStyle -> {RGBColor[0, 0, 1]}, Frame -> True, 
          PlotRange -> All, 
          FrameLabel -> {"t[min]", "i[mU/min]"}];

Fig. 9. Exogenous insulin infusion rate in case of ideal state feedback

Simulation with observer

Now in the  control law,  we should use the recovered states for X(t) and Y(t) estimated by the observer, namely

Xrec[t_] := 
      First[K1[[2]] (G[t]) + K2[[2, 1]] (z1[t]) + 
              K2[[2, 2]] (z2[t]) /. numPar /. {(Subscript[p, 1]) -> -1, 
            (Subscript[p, 2]) -> -1}]

Yrec[t_] := 
      First[(K1[[3]] (G[t]) + K2[[3, 1]] (z1[t]) + 
                  K2[[3, 2]] (z2[t]) /. 
                numPar /. {(Subscript[p, 1]) -> -1, (Subscript[p, 2]) -> -1}) // 
          Flatten]

then

i[t_] := i0 - KLQ[[2]] . {G[t], Xrec[t], Yrec[t]}

We solve the state equations with the observer equations, simultaneously

sol = 
      NDSolve[{deq1, deq2, deq3, deq4, deq5, X[0] == 0, 
          Y [0] == 0, G[0] == 0, z1[0] == 0, 
          z2[0] == 0}, {X, Y, G, z1, z2}, {t, 0, 1000}, 
        Method -> BDF,  MaxSteps -> 50000]

The controlled dynamics of the state variables in case of observer

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

P5 = 
        Plot[Evaluate[(X[t] /. sol)], {t, 0, 1000}, 
          PlotStyle -> {RGBColor[1, 0, 0]}, Frame -> True, 
          PlotRange -> All, 
          FrameLabel -> {"t[min]", "X[1/min]"}, 
          DisplayFunction -> Identity];

P6 = 
        Plot[Evaluate[(Y[t] /. sol)], {t, 0, 1000}, 
          PlotStyle -> {RGBColor[1, 0, 0]}, Frame -> True, 
          PlotRange -> All, 
          FrameLabel -> {"t[min]", "Y[mU/dL]"}, 
          DisplayFunction -> Identity];

PC2 = 
        Plot[Evaluate[(i0 - 
                KLQ[[2]] . {G[t] /. sol, Xrec[t] /. sol, 
                    Yrec[t] /. sol})], {t, 0, 1000}, 
          PlotStyle -> {RGBColor[1, 0, 0]}, Frame -> True, 
          PlotRange -> All, 
          FrameLabel -> {"t[min]", "i[mU/min]"}, 
          DisplayFunction -> Identity];

Now we can compare the dynamical performance of  state variables concerning the two simulations

Show[{P1, P4}, 
        DisplayFunction -> $DisplayFunction];

Fig. 10 Controlled plasma glucose deviation in  case of  observer (red) and in case of ideal state feedback (blue)

Show[{P2, P5}];

Fig.11 Controlled remote compartment insulin utilization in case of observer (red) and in case of ideal state feedback (blue)

Show[{P3, P6}];

Fig. 12  Controlled plasma insulin deviation in case of observer (red) and in case of ideal state feedback (blue)

Show[{PC1, PC2}];

Fig.13. Exogenous insulin infusion rate in case of observer (red) and in case of ideal state feedback (blue)

There is only a slight difference in the dynamical performances.

Conclusions

Numerical simulations show, that ten times faster observer - means, that the absolute value of the greatest eigenvalue of the observer (F matrix) is ten times greater than that of the system - can ensure satisfactory state recovery for the unmeasurable state variables, remote compartment insulin utilization and plasma insulin deviation. This observer provide good performance in case of standard meal disturbance (exogenous glucose infusion) for classical LQ as well as for disturbance rejection control. Symbolic form of the observer makes easy to change parameter values without designing the observer again.

References

Patel, R. V. and N. Munro, Multivariable System Theory and Design, Oxford: Pergamon Press, 1982.

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.

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.

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.