 # Filtering and Prediction of the Blood Glucose Concentration using an Android Smart Phone and a Continuous Glucose Monitor

In this paper we numerically assess the performance of Java linear algebra libraries for the implementation of nonlinear filters in an Android smart phone (Samsung A5 2017). We implemented a linear Kalman filter (KF), an extended Kalman filter (EKF), and an unscented Kalman filter (UKF). These filters are used for state and parameter estimation, as well as fault detection and meal detection in an artificial pancreas. We present the state estimation technologies used for glucose estimation based on a continuous glucose monitor (CGM). We compared three linear algebra libraries: The Efficient Java Matrix Library (EJML), JAMA and Apache Common Math. Overall, EJML provides the best performance for linear algebra operations. We demonstrate the implementation and performance of filtering (KF, EKF and UKF) using real CGM data.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## I Introduction

Kalman filtering is used to estimate the states of a linear state-space model where discrete-time measurements are available. Variations of the Kalman filter also exist for nonlinear systems. The extended Kalman filter (EKF) uses a first order approximation to update the mean and covariance of the state vector. In the unscented Kalman filter (UKF)

[1, 2]

, we select a number of representative points around the mean value (also called sigma points) to numerically evaluate the mean and covariance of the state vector. The particle filter (PF) and ensemble Kalman filter (EnKF) propagate a large number of points to estimate the probability distribution of the states

[3, 4]

. Kolmogorov forward equations (also called Fokker-Planck equations) provide an analytical expression of the state probability density function

.

When the measurement sampling time is fast enough compared to the nonlinearities, the EKF and the UKF provide a suitable approximation of the true system. The PF and the EnKF are more accurate than the EKF and the UKF and can capture non-Gaussian distributions, but suffer from the curse of dimensionality. The Kolmogorov forward equations are exact but need to solve a system of partial differential equations and cannot be numerically solved in a reasonable time when the number of states become too large. The artificial pancreas (AP) is an example of a closed-loop control system. Filtering is important to ensure good performance of the closed-loop algorithm

[6, 7, 8, 9, 10, 11, 12, 13, 14]. The AP gets measurements from a continuous glucose monitor (CGM), computes the optimal insulin dose, and sends this dose to a continuous subcutaneous insulin infusion (CSII) pump. In the AP, filtering and prediction are used for model identification , detection of meals , detection of faulty CGM measurements , as well as in the model predictive control (MPC)  algorithm. The filtering is also important for CGM enabled insulin pen systems [18, 19]. For the AP, EKF and UKF can be used as filtering algorithms and perform better than the PF  due to the fast CGM sampling time, typically every five minutes. More generally speaking, filtering, detection and numerical optimization in the AP involve linear algebra operations for small matrices. In this paper, we evaluate the performance of linear KF and nonlinear filters (EKF and UKF) implemented in an Android platform. We compare three Java linear algebra libraries and demonstrate the implementation of KF, EKF and UKF.

The rest of the paper is structured as follows. Section II presents the linear Kalman filter (KF), and Section III explains the formulation for the extended Kalman filter (EKF), and unscented Kalman filter (UKF). Section IV describes the Java implementation of the filter in an Android smart phone, and the hardware configuration. The results are presented in Section V. Section VI summarizes the main contributions of this paper.

## Ii Linear Model

We use is the Medtronic virtual patient (MVP) type 1 diabetes patient model described in . The MVP is a nonlinear model and we use it for the nonlinear KF. For the linear KF, we identify a second-order transfer function model for insulin to sc glucose dynamics and for the CHO to sc glucose dynamics. The model in the Laplace domain is defined as

 Z(s)=Gu(s)U(s)+Gd(s)D(s), (1a) where Gu(s)=Ku(τus+1)2, (1b) Gd(s)=Kd(τds+1)2. (1c)

The output, , is the sc glucose concentration as the deviation variable, which is the deviation from the steady-state sc glucose concentration (100 mg/dL). is the sc insulin input rate (IU/min), and is the CHO ingestion rate (g/min). and are also deviation variables. The transfer functions, and , are the Laplace transforms of the insulin and CHO impulse responses, respectively. The gains, and , correspond to the steady state change in BG for a unit step in the inputs, and the time constants, and , determine the time to reach the steady state. Mahmoudi et al  present a method for identification of the parameters, , , , and . The identified model is then converted to a linear time-invariant discretized state-space model , which is

 xk+1 =Axk+Buk+E(dk+wk), (2a) zk =Cxk, (2b) yk =zk+vk, (2c)

with

 wk ∼Niid(0,Q), (3) vk ∼Niid(0,R), (4)

where is the subcutaneous insulin infusion rate, is the meal ingestion, is the CGM sensor measurement, and is the sensor noise . For the linear model, We assume that the uncertainty enters the process through the meal ingestion and therefore, we add the process noise, , to the CHO input and use as a tuning parameter . For the linear KF implementation, we compute the stationary state covariance, and the stationary filter gain as follows.

The discrete algebraic Riccati equation (DARE),

 P=APAT+EQET−(APCT)(CPCT+R)−1(APCT)T, (5)

for the KF of the discrete-time stochastic linear model (2) may be used to compute the stationary covariance matrix of the one-step prediction, . The corresponding matrices , , and are computed by

 Re,∞ =CPCT+R, (6a) K∞ =PCTR−1e,∞, (6b) P∞|∞ =P−K∞Re,∞KT∞. (6c) Re,∞ is the innovation covariance. Re,k =CPk|k−1C′+R. (6d)

### Ii-a Kalman Filter

The KF is initialized with .

#### Ii-A1 One-step prediction

The one-step predictions of the states and the covariance are

 ^xk|k−1 =A^xk−1|k−1+Buk−1+Edk−1, (7a)

and the one-step predictions of the outputs and the measurements are

 ^zk|k−1 =C^xk|k−1, (8a) ^yk|k−1 =^zk|k−1. (8b)

#### Ii-A2 Measurement update

When the new measurements, , become available, we compute the innovation,

 ek=yk−^yk|k−1, (9)

as well as the filtered states and the corresponding covariance matrix,

 ^xk|k =^xk|k−1+K∞ek. (10)

## Iii Nonlinear Model

In this section, we describe the EKF and the UKF for stochastic differential equations (SDEs) in the form

 dx(t) =f(x(t),u(t),d(t),p)dt+g(x(t),u(t),d(t),p)dω(t). (11)

The model, , is the MVP model , and in our application the function, , is an additive constant diffusion coefficient and therefore the model reduces to

 dx(t) =f(x(t),u(t),d(t),p)dt+σdω(t). (12)

As in Section II, are the states, are the manipulated inputs, are disturbance variables, and are parameters. is a standard Wiener process, i.e. the incremental covariance of is

. At the initial time, the states are normally distributed:

. The first term in (12) is the drift term which represents the deterministic part of the model equations. The second term is the diffusion term which represents the uncertainty in the process, i.e. the process noise. The process outputs,

 z(t)=h(x(t),p), (13)

are measured at discrete points in time, :

 y(tk)=z(tk)+v(tk). (14)

is the measurement noise.

### Iii-a Extended Kalman Filter

The EKF is initialized using the distribution of the initial states, i.e. and .

#### Iii-A1 One-step prediction

The one-step prediction of the states and their covariance are the solutions to

 ^xk−1(tk−1) =^xk−1|k−1, (15a) Pk−1(tk−1) =Pk−1|k−1, (15b) ddt^xk−1(t) =f(^xk−1(t),u(t),d(t),p), (15c) ddtPk−1(t) =Ak(t)Pk−1(t)+Pk−1(t)Ak(t)′+Qk(t), (15d)

for where

 Ak(t) =∂f∂x(^xk−1(t),u(t),d(t),p), (16a) Qk(t) =σσ′. (16b)

The one-step predictions are

 ^xk|k−1 =^xk−1(tk), (17a) Pk|k−1 =Pk−1(tk), (17b)

and the one-step prediction of the measurements is

 ^yk|k−1 =^zk|k−1=h(^xk|k−1,p). (18)

The covariance of the innovation and the Kalman gain are

 Re,k =CkPk|k−1C′k+R, (19a) Kfx,k =Pk|k−1C′kR−1e,k, (19b)

where is the Jacobian of the sensor model:

 Ck =∂h∂x(^xk|k−1,p). (20)

#### Iii-A2 Measurement update

When the measurements, , become available, we compute the innovation,

 ek =yk−^yk|k−1, (21)

and the filtered estimates of the states and its covariance:

 ^xk|k =^xk|k−1+Kfx,kek, (22a) Pk|k =(I−Kfx,kCk)Pk|k−1(I−Kfx,kCk)′ +Kfx,kRK′fx,k. (22b)

We use the Joseph-stabilized form of the covariance update in order to ensure positive definiteness of the updated covariance matrix as well as numerical stability.

### Iii-B Unscented Kalman Filter

As for the EKF, the UKF is initialized with and .

#### Iii-B1 One- step prediction

First, we introduce the weights as 

 λ=α2(n+k)−n,c=α2(n+k), (23a) W(0)m=λ/(n+λ), (23b) W(0)c=λ/(n+λ)+(1−α2+β), (23c) W(i)m=1/{2(n+λ)},i=1,…,2n, (23d) W(i)c=1/{2(n+λ)},i=1,…,2n, (23e) Wm=[W(0)m…W(2n)m]T, (23f) Wc=[W(0)c…W(2n)c], (23g) where n is the number of states in the state-space model. We use α=0.01, k=0, and β=2 in this work. Next, the sigma points are generated from the previous filtered states and covariance. ^x(0)k−1|k−1=^xk−1|k−1 (23h) ^x(i)k−1|k−1=^xk−1|k−1+√c(√Pk−1|k−1)i,i=1,…,n (23i) ^x(i+n)k−1|k−1=^xk−1|k−1−√c(√Pk−1|k−1)i,i=1,…,n. (23j) We use a Cholesky factorization to compute √Pk−1|k−1, and (√Pk−1|k−1)i denotes the i’th column of √Pk−1|k−1 . For each of the sigma points we solve the ODE, ddt^x(i)k−1(t)=f(^x(i)k−1(t),u(t),d(t),p),i=0,…,2n, t∈[tk−1 tk]. (23k) In parallel, the covariance matrix should be also propagated in time according to dPk−1(t)dt= 2n∑i=0W(i)c(^x(i)k−1(t)−^xk−1(t))(F(i)(t)−F(t))T +2n∑i=0W(i)c(F(i)(t)−F(t))((^x(i)k−1(t)−^xk−1(t))T +σσT,t∈[tk−1 tk], (23l) where, ^xk−1(t)=2n∑i=0W(i)m ^x(i)k−1(t), (23m) F(i)(t)=f(^x(i)k−1(t),u(t),d(t),p), (23n) F(t)=2n∑i=0W(i)m F(i)(t). (23o) The predicted mean and covariance of the states are computed as ^x(i)k|k−1=^x(i)k−1(tk),i=0,…,2n, (23p) ^xk|k−1=2n∑i=0W(i)m ^x(i)k|k−1, (23q) Pk|k−1=Pk−1(tk). (23r)

An intermediate set of sigma points, , can be generated using the predicted state mean and covariance,

 ^~x(0)k|k−1=^xk|k−1, (24a) ^~x(i)k|k−1=^xk|k−1+√c(√Pk|k−1)i,i=1,…,n (24b) ^~x(i+n)k|k−1=^xk|k−1−√c(√Pk|k−1)i,i=1,…,n. (24c)

The predicted outputs corresponding to the sigma points are computed according to

 ^y(i)k|k−1=^z(i)k|k−1=h(^~x(i)k|k−1,p),i=0,…,2n, (24d)

and the mean of the predicted output is

 ^yk|k−1=2n∑i=0W(i)m^y(i)k|k−1. (24e)

Generation of the intermediate set of sigma points can be omitted for the sake of lowering the computational effort . In the UKF implementation, we omitted this step and instead we used the sigma points, , in (23q) to compute in (24d).

The innovation and its covariance are defined as

 Re,k=2n∑i=0W(i)c(^y(i)k|k−1−^yk|k−1)(^y(i)k|k−1−^yk|k−1)T +R. (24f)

The cross-covariance between the predicted states and outputs is

 Rxy,k=2n∑i=0W(i)c (^x(i)k|k−1−^xk|k−1)(^y(i)k|k−1−^yk|k−1)T. (24g)

The filter gain is calculated according to

 Kfx,k=Rxy,k R−1e,k, (24h)

#### Iii-B2 Measurement update

When the measurement is available, we compute the innovation,

 ek=yk−^yk|k−1, (24i)

and the filtered state mean and covariance are

 ^xk|k =^xk|k−1+Kfx,kek, (24j) Pk|k=(I−Kfx,kCk)Pk|k−1(I−Kfx,kCk)′+Kfx,kRK′fx,k. (24k)

As in the EKF, we use the Joseph-stabilized form of the covariance update.

## Iv Implementation

### Iv-a Java implementation

We test three linear algebra libraries: JAMA 1.0.3 (http://math.nist.gov/javanumerics/jama/), EJML 0.36 (http://ejml.org), and Apache Commons Math 3.6.1 (http://commons.apache.org/proper/commons-math/). We evaluate the performance of these libraries by comparison of the computation times for the main operations required in the KF, EKF and UKF: Matrix-matrix multiplications, Cholesky decomposition, matrix scaling, and matrix-matrix addition. The prediction and measurement update in the filters require matrix-matrix multiplication (all filters), Cholesky decomposition (EKF and UKF), matrix scaling (KF and UKF), and matrix-matrix addition (all filters).

Fig. 1 shows a few lines of the Java code from the implementation of the filters. Fig. 1: Java implementation for reading the inputs and filtering the CGM data. The inputs of the filter are the meal and bolus (meal bolus and the extended bolus) information and the CGM data from a specific time in thje past until the current time.

The basal insulin rate is directly added during the one-step prediction. The filtered data is then sent to the plot function. Fig. 2 shows the code snippet for the Cholesky decomposition for the UKF in EJML. Fig. 2: Cholesky decomposition for the UKF in EJML Java library.

Fig. 3 shows the code snippet for computation of the innovation covariance, the cross covariance, and the filter gain for the UKF in EJML. Fig. 3: Computation of innovation covariance, the cross-covariance in (24g) and the filter gain in (24h), for the UKF in EJML Java library.

### Iv-B Hardware and data acquisition

The filters are implemented in a Samsung A5 2017 running Android 8.1.0. We test implementation in Java of the filters in connection to a real CGM. To monitor and filter CGM data, we use AndroidAPS 1.58 (http://github.com/MilosKozak/AndroidAPS). AndroidAPS is an open source Java-based artificial pancreas application developed for Android phones. We added the filter as a plugin to AndroidAPS. The CGM is a FreeStyle Libre from Abbott Laboratories connected to a Blucon Bluetooth transmitter from Ambrosia Systems Inc. Fig. 4 illustrates the hardware configuration. Fig. 5 shows the filter implemented as a general plugin in AndroidAPS. Fig. 4: Hardware used in the study consisting of FreeStyle Libre sensor, BluCon transmitter, and an Android smart phone. The BluCon transmitter is attached to the FreeStyle Libre CGM sensor. Fig. 5: The Filter implemented as a general Plugin in AndroidAPS.

## V Results and discussion

### V-a Comparison of Java linear algebra libraries

Fig. 6 shows the computation time for four matrix algebraic operations in JAMA, EJML, and Commons math libraries in Java. The comparison indicates that EJML is faster than JAMA and Commons math libraries for Cholesky decomposition, and matrix-matrix addition. EJML is relatively comparable with the other two libraries for the matrix-matrix multiplication and matrix scaling. Therefore we choose to implement the filters in EJML.

### V-B Performance of linear and nonlinear Kalman filters

Fig. 7 compares the processing time of the three filters in EJML. For the linear KF, we use a stationary filter, i.e. the state covariance matrix, the Kalman gain and the covariance of the innovation can be computed off-line. We use a forward Euler method to compute the one-step prediction of the state estimate for the UKF, and to compute the one-step prediction of the state covariance for the EKF. The step length is 1 minute, and the sampling time is 5 minutes, i.e. the one-step prediction requires five Euler steps.

Fig. 7 indicates that the linear KF has the smallest computation time for the prediction. This is expected as the one-step prediction in (7) does not need an Euler method implementation. In addition, the filtering step for the linear KF has the smallest computation time among the three filters, because it does not include a covariance update as the stationary state covariance is computed off-line outside the filtering step. The filtering computation time is similar for the EKF and UKF, while the one-step prediction for the UKF is more time-consuming than that for the EKF. For the UKF, it is required to compute the Cholesky factorization of , the one-step predictions for the sigma points and to compute the state covariance matrix, whereas the EKF requires to compute the one-step prediction of states, and the state covariance of size . Fig. 7: Comparison of computation time for prediction step and filtering step, for Kalman filter, extended Kalman filter, and unscented Kalman filter. The values are the average over 10000 runs. (a) CGM filtered by the UKF

Fig. 8 shows the filtered and unfiltered CGM data. The filter is the UKF implemented in the Android smart phone using the EJML Java Matrix library in AndroidAPS platform. The filtered data is shown in the AndroidAPS graph when the Filter plugin in AndroidAPS is enabled (Fig. 5). The unfiltered data are shown in the graph when the Filter plugin in AndroidAPS is disabled. The sensor is connected to a non-diabetes subject, while the model in the UKF is a diabetes patient model. Obviously, the model does not fit the CGM data and that is the reason for the relatively large bias between the unfiltered CGM data and the filtered data. This is of less importance in this study, because the aim of this implementation is to measure the computation time for the filters and solve software/hardware enginering issues with the Java implementation and the sensor connection.

## Vi Conclusions

We presented the linear and nonlinear KF and implemented them in an Android smart phone using the Java programming language. We tested these filters in the AndroidAPS platform and use a FreeStyle Libre sensor as the source of CGM data. We compare the computation times of elementary matrix operations for three numerical linear algebra libraries implemented in Java. These operations are required for linear and nonlinear filtering. The results indicate that the EJML Java library is a suitable package for this purpose. The state-estimator in this study can be used for the implementation in Java of linear and nonlinear parameter and state estimation, model-based controller, model-based sensor fault detection, as well as for meal detection and estimation. These filters and in particular the UKF can be used as part of an artificial pancreas, in model-based bolus advisors for connected insulin pens, and in meal-detection algorithms.

## References

•  S. J. Julier, J. K. Uhlmann, and H. F. Durrant-Whyte, “A new method for the nonlinear transformation of means and covariances in filters and estimators.” IEEE Transactions on Automatic Control, vol. 45, no. 3, pp. 477–482, 2000.
•  S. Särkkä, “On unscented Kalman filtering for state estimation of continuous-time nonlinear systems,” IEEE Transactions on Automatic Control, vol. 52, no. 9, pp. 1631–1641, 2007.
•  G. Burgers, P. Jan van Leeuwen, and G. Evensen, “Analysis scheme in the ensemble kalman filter,” Monthly weather review, vol. 126, no. 6, pp. 1719–1724, 1998.
•  A. Tulsyan, R. B. Gopaluni, and S. R. Khare, “Particle filtering without tears: A primer for beginners,” Computers & Chemical Engineering, vol. 95, pp. 130–145, 2016.
•  S. Challa and Y. Bar-Shalom, “Nonlinear filter design using Fokker-Planck-Kolmogorov probability density evolutions,” IEEE Transactions on Aerospace and Electronic Systems, vol. 36, no. 1, pp. 309–315, 2000.
•  C. Cobelli, E. Renard, and B. Kovatchev, “Artificial pancreas: past, present, future,” Diabetes, vol. 60, pp. 2672 – 2682, 2011.
•  S. Trevitt, S. Simpson, and A. Wood, “Artificial pancreas device systems for the closed-loop control of type 1 diabetes: What systems are in development?” Journal of Diabetes Science and Technology, vol. 10, no. 3, pp. 714–723, 2016.
•  D. Boiroux, D. A. Finan, N. K. Poulsen, H. Madsen, and J. B. Jørgensen, “Meal estimation in nonlinear model predictive control for type 1 diabetes,” IFAC Proceedings Volumes, vol. 43, no. 14, 2010.
•  D. Boiroux, D. A. Finan, J. B. Jørgensen, N. K. Poulsen, and H. Madsen, “Nonlinear model predictive control for an artificial -cell,” in Recent Advances in Optimization and its Applications in Engineering.   Springer, 2010, pp. 299–308.
•  S. Schmidt, D. Boiroux, A. K. Duun-Henriksen, L. Frøssing, O. Skyggebjerg, J. B. Jørgensen, N. K. Poulsen, H. Madsen, S. Madsbad, and K. Nørgaard, “Model-based closed-loop glucose control in type 1 diabetes: The DiaCon experience.” Journal of Diabetes Science and Technology, vol. 7, no. 5, pp. 1255–1264, 2013.
•  D. Boiroux, A. K. Duun-Henriksen, S. Schmidt, K. Nørgaard, S. Madsbad, N. K. Poulsen, H. Madsen, and J. B. Jørgensen, “Overnight glucose control in people with type 1 diabetes,” Biomedical Signal Processing and Control, vol. 39, pp. 503–512, 2018.
•  D. Boiroux, A. K. Duun-Henriksen, S. Schmidt, K. Nørgaard, N. K. Poulsen, H. Madsen, and J. B. Jørgensen, “Adaptive control in an artificial pancreas for people with type 1 diabetes,” Control Engineering Practice, vol. 58, pp. 332–342, 2017.
•  D. Boiroux and J. B. Jørgensen, “A nonlinear model predictive control strategy for glucose control in people with type 1 diabetes,” in 10th IFAC Symposium on Biological and Medical Systems (BMS 2018), 2018, pp. 227–232.
•  ——, “Nonlinear model predictive control and artificial pancreas technologies,” in 57th IEEE Conference on Decision and Control (CDC 2018), 2018, pp. 284–290.
•  D. Boiroux, M. Hagdrup, Z. Mahmoudi, N. K. Poulsen, H. Madsen, and J. B. Jørgensen, “Model identification using continuous glucose monitoring data for type 1 diabetes,” IFAC-PapersOnLine, vol. 49, no. 7, pp. 759–764, 2016.
•  Z. Mahmoudi, F. Cameron, N. K. Poulsen, H. Madsen, B. W. Bequette, and J. B. Jørgensen, “Sensor-based detection and estimation of meal carbohydrates for people with diabetes,” Biomedical Signal Processing and Control, vol. 48, pp. 12–25, 2019.
•  Z. Mahmoudi, K. Nørgaard, N. K. Poulsen, H. Madsen, and J. B. Jørgensen, “Fault and meal detection by redundant continuous glucose monitors and the unscented Kalman filter,” Biomedical Signal Processing and Control, vol. 38, pp. 86–99, 2017.
•  D. Boiroux, T. B. Aradóttir, M. Hagdrup, N. K. Poulsen, H. Madsen, and J. B. Jørgensen, “A bolus calculator based on continuous-discrete unscented Kalman filtering for type 1 diabetics,” IFAC-PapersOnLine, vol. 48, no. 20, pp. 159–164, 2015.
•  D. Boiroux, T. B. Aradóttir, K. Nørgaard, N. K. Poulsen, H. Madsen, and J. B. Jørgensen, “An adaptive nonlinear basal-bolus calculator for patients with type 1 diabetes,” Journal of Diabetes Science and Technology, vol. 11, no. 1, pp. 29–36, 2017.
•  Z. Mahmoudi, S. L. Wendt, D. Boiroux, M. Hagdrup, K. Nørgaard, N. K. Poulsen, H. Madsen, and J. B. Jørgensen, “Comparison of three nonlinear filters for fault detection in continuous glucose monitors,” in IEEE 38th Annual International Conference of the Engineering in Medicine and Biology Society (EMBC 2016), 2016, pp. 3507–3510.
•  S. Kanderian, S. Weinzimer, G. Voskanyan, and G. Steil, “Identification of intraday metabolic profiles during closed-loop glucose control in individuals with type 1 diabetes,” Journal of Diabetes Science and Technology, pp. 1047–1057, 2009.
•  L. Biagi, C. Ramkissoon, A. Facchinetti, Y. Leal, and J. Vehi, “Modeling the error of the medtronic paradigm veo enlite glucose sensor,” Sensors, vol. 17, no. 6, p. 1361, 2017.
•  Z. Mahmoudi, D. Boiroux, and J. B. Jørgensen, “Meal detection for type 1 diabetes using moving horizon estimation,” in 2018 IEEE Conference on Control Technology and Applications (CCTA), 2018, pp. 1674–1679.
•  S. J. Julier, “The scaled unscented transformation,” in Proceedings of the American Control Conference, Anchorage, Alaska, May 2002.
•  D. Simon, Optimal state estimation: Kalman, H, and nonlinear approaches.   John Wiley & Sons, 2006.