- Email: [email protected]
Identification of a nonlinear model for a glucoregulatory benchmark problem
Biomedical Signal Processing and Control 13 (2014) 168–173
Contents lists available at ScienceDirect
Biomedical Signal Processing and Control journal homepage: www.elsevier.com/locate/bspc
Identification of a nonlinear model for a glucoregulatory benchmark problem Laurent Vanbeylen a,∗ , Anne Van Mulders a , Amjad Abu-Rmileh b a b
Department ELEC, Vrije Universiteit Brussel, Belgium Department of Brain and Cognitive Sciences, Ben-Gurion University of the Negev, Israel
a r t i c l e
i n f o
Article history: Received 15 January 2014 Received in revised form 1 April 2014 Accepted 24 April 2014 Keywords: Insulin-glucose modelling Artificial pancreas System identification Black-box modelling Best linear approximation Nonlinear models Nonlinear fractional representation
a b s t r a c t Recently, a novel identification method for a nonlinear dynamic model, called nonlinear Linear Fractional Representation (NL-LFR) model, has been developed. The model, composed of a static nonlinearity (SNL) surrounded by linear dynamics, can account for both nonlinear feed-forward and nonlinear feedback effects. Using two classical frequency response measurements, the SNL is automatically recovered in a user-friendly and efficient (non-iterative) way. In this contribution, the method is illustrated on a glucoregulatory benchmark dataset (insulin–glucose relationship of the human body). The research on insulin–glucose models is essential to develop methodologies to control the blood glucose level in diabetes patients. The obtained results outperform earlier results on the same benchmark data, while providing an excellent accuracy-complexity tradeoff. © 2014 Elsevier Ltd. All rights reserved.
1. Introduction To be able to describe measured data with an increased accuracy, e.g. over a wide amplitude range, it is important to explicitly incorporate nonlinearities into the model. Until now, a range of nonlinear models have been proposed [1–9]. However, most of them either require a (typically) quite large number of parameters (i.e., they are not parsimonious enough) or are not well suited to represent ubiquitous nonlinear feedback effects, such as amplitude-dependent resonance frequencies and dampings (i.e., they are not flexible enough). The nonlinear Linear Fractional Representation (NL-LFR) model belongs to the popular class of nonlinear block-oriented models. It offers a good flexibility-parsimony tradeoff, and is therefore applicable to a wide range of scientific domains. Block-oriented nonlinear models [1,7] can be defined as interconnections of linear dynamic elements (a.k.a. Linear Time Invariant parts or LTI blocks) and static nonlinear elements (a.k.a. Static NonLinearities or SNL blocks). The most studied block-oriented model structures are cascade structures of the form LTI-SNL, SNL-LTI and LTI-SNLLTI (known as Wiener, Hammerstein and Wiener–Hammerstein model, respectively).
∗ Corresponding author. Tel.: +32 26292979. E-mail addresses: [email protected] (L. Vanbeylen), [email protected] (A. Van Mulders), [email protected] (A. Abu-Rmileh). http://dx.doi.org/10.1016/j.bspc.2014.04.007 1746-8094/© 2014 Elsevier Ltd. All rights reserved.
The NL-LFR model, to be discussed in Section 3, is more general than these cascades and surrounds the SNL by arbitrary dynamics. We refer to [10] for a complete description of these matters and of the benefits of the model. In this paper, an insulin–glucose benchmark problem (see Section 2) is presented as an application of the NL-LFR model. The results and comparison to earlier results with different model structures are shown in Sections 4–6. Section 7 concludes the paper. 2. Insulin–glucose modelling problem 2.1. Aim To enhance the lives of (type 1) diabetes (mellitus) patients, research is conducted on the regulation of the blood glucose level [11]. E.g., glucose control schemes for an artificial pancreas have been developed for the insulin–glucose metabolism [12], based on a Wiener (LTI-SNL) model and nonlinear sliding mode control. More sophisticated simulation insulin–glucose models built from first principles exist [13,14], but are avoided for control due to a large computational effort involved with the implementation and model parameters are very difficult to estimate from the available patient data. In this work, as in [15], simplified nonlinear models, for use in control in a later stage, represent the system as a replacement for the sophisticated models. To estimate a simplified nonlinear model, a data-driven, measurement point of view is taken here: the insulin is treated as a measured input and the glucose concentration
L. Vanbeylen et al. / Biomedical Signal Processing and Control 13 (2014) 168–173
Fig. 1. (Concatenated) input and output signals for estimation (dataset A): randomphase multisines at 12 operating points.
as a measured output of the glycoregulatory system. In this paper, estimation results for the NL-LFR model are presented, and are compared with the Wiener and other nonlinear models. The data from which the simplified nonlinear models are estimated, are priorly generated through the Dalla-Man model [14]. The same data as in [15] are used: partly with a multisine1 insulin input (dataset A) and partly a multisine input with band-limited pulses superimposed, intended to simulate a more realistic situation (dataset B). The data are available at 5 multisine amplitudes and 12 DC-levels. Please note that the estimation of a (simplified) nonlinear model is based on simulated data. The simple reason is that performing real measurements on patients can be very expensive, and that high-quality results on simulated data may be required to convince the medical experts that performing a real measurement campaign makes sense. 2.2. Benchmark datasets The benchmark datasets used consist of steady-state2 measurements of the glucose response of the Dalla-Man model to an insulin input at different basal (DC) glucose levels, while keeping the meal input zero [15]. Note that, for both datasets, only the data at a single input std3 -value is taken as estimation data and another (slightly lower) std-value is taken as validation data, only intended to assess the model quality. The input (insulin4 ) and output (glucose) signals for datasets A and B are depicted in Figs. 1 and 2, respectively. The properties of both datasets are listed in Table 1. To compare the results, the following measure for the fit quality is used on the validation data:
Fit% =
y − yˆ 2 1− ¯ 2 y − y
169
Fig. 2. (Concatenated) input and output signals for estimation (dataset B): bandlimited pulses superimposed on a random-phase multisine at 12 operating points.
measured output (constant signal), respectively. Distinct values for the fit quality for the different operating points will be considered.
3. Nonlinear LFR model The focus of this paper is to identify an NL-LFR model on the insulin–glucose benchmark data. The model’s input and output are represented by u and y, as shown in Fig. 3. The SNL has an input z and output v. It has no memory effects, and can therefore be considered as a simple nonlinear mapping fSNL : R → R operating on each input sample z(k) (with discrete time index k):
v(k) = fSNL (z(k))
(1)
The linear part (LTI) connects the SNL and the model’s input and output: the MIMO LTI converts the two inputs, u and v, to two outputs, y and z. Signals v and z are internal and not accessible for measurement. The MIMO-LTI part can, in general, be described as: x(k + 1) = Ax(k) + Bu(k)
(2)
y(k) = Cx(k) + Du(k)
× 100
with y, yˆ , and y¯ the steady-state measured output signal, the steadystate simulated model output signal and the mean value of the
1 The use of multisine signals (among other types of variations), to represent low amplitude slow fluctuations around a basal insulin level is known in the field of artificial pancreas, mainly for the identification of control-relevant models [16]. The use of such signals taking into account patients safety limits, would be possible. 2 In the context of this paper, “steady-state” means “after die-out of the transient phenomena”, i.e. the periodic regime response is considered. 3 std = standard deviation 4 It can be remarked that, at some places, negative insulin values are visible. This is not realistic. However, the idea behind using negative values is to implement fluctuations below the basal (viz., average) insulin levels.
Fig. 3. The NL-LFR model. The signals u, y, z and v represent the model input and output and the SNL’s input and output, respectively (the identification method has no access to measurements of z and v).
170
L. Vanbeylen et al. / Biomedical Signal Processing and Control 13 (2014) 168–173
Table 1 Properties of both datasets.
Kind of excitation Sampling period Ts Number of data samples (per period) Number of periods Bandwidth of the input [1/min] Input std range [pmol/min] Input std for estimation [pmol/min] Input std for validation [pmol/min] Input o/pa range [pmol/min] Total number of experiments a
Dataset A
Dataset B
Odd random-phase multisine 20 min 2000 2 0.01 26–176 105 60 115.0–539.0 5 std-values × 12 DC-values
Band-limited pulses superimposed on an odd random-phase multisine 1 min 40,000 2 0.33 115–316 212 177 154.0–539.0 5 std-values × 12 DC-values
o/p = operating point (DC-offset) = mean value = basal level.
with states x(k) ∈ Rn (order n), input vector u(k) =
vector y(k) =
D=
Dyu Dzu
Dyv Dzv
y(k) z(k)
. The matrices B =
Bu
u(k) v(k)
Bv , C =
, output Cy Cz
4.1. Brief summary and
are splitted in subblocks. This results into:
x(k + 1) = Ax(k) + Bu u(k) + Bv v(k) y(k) = Cy x(k) + Dyu u(k) + Dyv v(k)
(3)
z(k) = Cz x(k) + Dzu u(k) Note that, as in [10], the coefficient Dzv , is set to zero to avoid algebraic loops when performing the simulation (computation of the response for a given signal u, SNL and state-space matrices). Discussion. Through the interconnection of the SNL (1) and the MIMO-LTI (3), the NL-LFR model features both nonlinearities and dynamics. Moreover, the dynamics surround the SNL in a flexible way, such that, in principle, any finite-dimensional system with a localised SNL can be captured by an NL-LFR model. Besides, please note that an alternative, equivalent representation for the model, shown in Fig. 4 exists, where the MIMO-LTI is split in four individual contributions. The transfer function matrix u y of the MIMO-LTI from inputs u = to outputs y = , conz v taining each of the four contributions, is connected as follows to the state-space representation:
Gyu (Z) =
GD (Z)
GH (Z)
GW (Z)
GFB (Z)
= D + C(ZI − A)−1 B.
4. Dataset A: identification methodology and results
(4)
where Z is the z-transform variable. The subscripts of the four contributions are associated with their positions w.r.t. the SNL in the block-scheme of Fig. 4: GD denotes the direct linear contribution from u to y, this signal path is not affected by the SNL; GFB denotes the feed-back which forms a nonlinear feedback loop together with the SNL; this loop is preceded and followed by the dynamics GW and GH , respectively. This corresponds to the positions of the linear dynamics of a Wiener and a Hammerstein system w.r.t. the SNL (viz., before and after, respectively).
Fig. 4. The NL-LFR model (LTI decomposed in four contributions). Note that the output y is here shown in the bottom left corner.
The identification of the NL-LFR model presented in the previous section can be carried out by following the lines of [10]. In a nutshell, the approach is based on data at different input conditions (e.g., amplitudes or operating points), exciting each time the nonlinearity differently. From each part of the data the Best Linear Approximation (BLA) is extracted. By solving a Riccati equation, the MIMO-LTI part (viz., the common linear dynamics in the NL-LFR) can be obtained. In a last step, the internal signals z(k) and v(k) are reconstructed, from which the SNL is estimated. The MIMO-LTI and the SNL together constitute the initial nonlinear model. 4.2. Detailed description The NL-LFR identification method starts with the identification of the BLA, the linear dynamics optimal in mean-square sense [17] for two different input conditions. Here, the data at input operating points 2 and 5 are selected. First, the two BLAs are estimated nonparametrically (see Fig. 5 for their Frequency Response Function (FRF) estimates). Next, each FRF is fitted as a parametric model, more precisely a 3rd order linear state-space model. This fit is obtained by means of a frequency-domain subspace technique. It is shown in Fig. 6. Given the two state-space models, a Nonsymmetric Algebraic Riccati Equation (NARE) is solved. This can be done very fast (see [18,19]) and in general leads to several possible solutions of the MIMO dynamics (here, there were 4 solutions). The input and output of the SNL are reconstructed for each of these solutions in the frequency domain (the data from all operating points are used). Next, a scatter plot between these reconstructed signals is generated (see Fig. 7), and the linearly parameterized SNL is fitted in one least-squares step (here, a polynomial of degree 4). The best fit
Fig. 5. Nonparametric FRFs for operating points 2 (grey dots) and 5 (black dots).
L. Vanbeylen et al. / Biomedical Signal Processing and Control 13 (2014) 168–173
171
Fig. 6. Nonparametric FRFs (black dots) and their parametric fit (black line) with the difference between both (grey dots). Top plot: operating point 2; bottom plot: operating point 5.
Fig. 8. 4th degree polynomial fit on the SNL for which the residual of the fit was the lowest. Top: nonparametric SNL (dots) and polynomial fit (black line); bottom: residuals of the fit.
Fig. 7. The nonparametric SNLs corresponding to the 4 possible solutions of the NARE. Light grey to black for operating points 1–12.
(best out of the 4 solutions shown in Fig. 7) of the SNL is given in Fig. 8. This solution (MIMO dynamics and SNL) is retained and used as an initial value in a further nonlinear-least-squares optimisation of the model parameters (matrix entries in (3) and polynomial coefficients). The FRFs of GD , GW , GH and GFB (corresponding to the final model) are given in Fig. 9.
Fig. 9. Four dynamic parts of the estimated MIMO-LTI block.
5. Dataset A: validation and comparison This section describes the validation result and compares it with other identification approaches (different model types). Fig. 10 shows the output signal of the validation data and the errors made by the two linear models used in the identification method. It can be seen that the linear models perform very well at the operating point they were fitted on, but fail completely outside this region.5 Moreover, the DC-errors (mean error values) were removed in these plots. Hence, the actual results are worse. The result of the NL-LFR model is shown in Fig. 11. Note already the difference in the order of magnitude of the error signals (see vertical axis): a relative error of the order of a few percent can be expected. The nonlinear model performs well, and this on all
5
Note also that a linear model fitted on the entire dataset would not produce good results: it would be poor everywhere.
Fig. 10. Dataset A. Top: (concatenated) output signal of the validation data; middle: error of the BLA constructed at operating point 2; bottom: error of the BLA constructed at operating point 5. Note that DC-errors are removed in this plot.
172
L. Vanbeylen et al. / Biomedical Signal Processing and Control 13 (2014) 168–173
Fig. 11. Dataset A. Top: (concatenated) output signal of the validation data; bottom: error of the NL-LFR model. No DC-correction was applied.
operating points. DC-levels are estimated automatically with the model: the error levels have no DC-correction. The fit-percentages for each of the 12 operating points of the models (linear and nonlinear) are reported in Table 2. In contrast to the linear models, the NL-LFR model can describe the entire operating region with a relatively high fit, in the range from 89.8% to 97.3%. It can be observed that the fit is not increasing for each operating point during the optimisation. This is because the optimisation balances the errors over the different set points due to the absence of weighting in the least-squares objective function optimised. The results of other nonlinear models are given in [15]. A schematic overview of these results, including the more recent result of this paper, is given in Fig. 12. The NL-LFR model presented in this paper clearly has a very good flexibility-parsimony tradeoff and outperforms the other results. 6. Dataset B: results and comparison In [15], the identification was redone for a (slightly) more realistic input signal (see Fig. 2), in which band-limited pulses are applied, superimposed by a multisine signal. This should mimic insulin shots administered to the patient. Without going into the details of the identification or results, this section provides a schematic overview of the results (see Fig. 13). The new dataset contained so many data points, that the nonlinear state-space (NLSS) models could no longer be obtained. The method presented in this paper did not suffer from this. It took only a few minutes to obtain an initial estimate with average fit 87.5%, Table 2 Dataset A. Fit% for each operating point. First two columns: results of the two BLAs; last two columns: results of the NL-LFR model before and after optimisation. The fits higher than 80% are indicated in boldface. o/p 1 2 3 4 5 6 7 8 9 10 11 12
BLA o/p 2
BLA o/p 5
NL-LFR init
NL-LFR
74.0 94.0 64.2 22.1 −10.3 −70.5 −141.0 −219.3 −303.4 −393.0 −485.1 −579.1
34.2 46.6 63.1 83.3 94.3 71.2 37.7 0.2 −40.1 −83.1 −127.4 −172.7
95.6 95.3 92.9 93.0 91.5 98.2 97.7 96.4 94.0 91.1 87.7 83.8
96.4 97.3 95.7 91.9 91.7 96.4 95.0 96.5 96.7 92.6 91.5 89.8
Fig. 12. Dataset A. Fit% as function of the number of independent model parameters (model complexity). The average fit (over all operating points) is indicated by a dot and the range from minimal to maximal fit quality (depending on the operating point) is given by a box. NN, (P)NLSS and SA denote neural network, (polynomial) nonlinear state-space, and state-affine model, respectively. The BLA depicted here corresponds to the operating point 7, located close to the center.
Fig. 13. Dataset B. Fit% as function of the number of independent model parameters. The same notation as in Fig. 12 is used.
minimal fit 78.6% and maximal fit 91.2%. Starting from this fit, only four iterations were performed, but the result is already quite good compared to the other results. The MIMO dynamic part of the NLLFR model has order n = 4 and the SNL has degree d = 3. 7. Conclusion In this paper, the aim was to obtain a simplified representation for a benchmark dataset representing the human glycoregulatory system (insulin–glucose behaviour), priorly simulated by means of the so-called Dalla-Man model. High-quality fits, outperforming earlier results, have been obtained with the nonlinear LFR model, consisting of a MIMO-LTI part and a polynomial SNL interconnected with each other. An excellent accuracy-complexity tradeoff is achieved with fit percentages of 90–97% and only a couple of tens of parameters. The model’s high quality can be attributed to its flexibility: the SNL is surrounded by arbitrary linear dynamics, taking both nonlinear feedforward and nonlinear feedback effects into account (which is not the case for Wiener-type models). The NLLFR model combines efficiently the power of the linear framework and the simplicity of block-oriented nonlinear structures. Acknowledgments Laurent Vanbeylen and Anne Van Mulders’s work was supported by the Fund for Scientific Research (FWO-Vlaanderen), the Flemish Government (Methusalem Grant METH-1), the Belgian Program on Inter-university Poles of Attraction (IAP VII/19 – Dysco), and
L. Vanbeylen et al. / Biomedical Signal Processing and Control 13 (2014) 168–173
by the ERC advanced grant SNLSID, under contract 320378. Amjad Abu-Rmileh acknowledges the support of the University of Girona through the BR-UdG research grant. References [1] S. Billings, S. Fakhouri, Identification of systems containing linear dynamic and static nonlinear elements, Automatica 18 (1) (1982) 15–26. [2] S. Boyd, L. Chua, Fading memory and the problem of approximating nonlinear operators with Volterra series, IEEE Trans. Circuits Syst. 32 (11) (1985) 1150–1161. [3] I. Hunter, M. Korenberg, The identification of nonlinear biological systems: Wiener and Hammerstein cascade models, Biol. Cybern. 55 (2) (1986) 135– 144. [4] S. Chen, S. Billings, Representation of non-linear systems: the NARMAX model, Int. J. Control 49 (3) (1989) 1012–1032. [5] G. Giannakis, E. Serpedin, A bibliography on nonlinear system identification, Signal Process. 81 (3) (2001) 533–580. [6] I. Goethals, K. Pelckmans, J.A. Suykens, B. De Moor, Identification of MIMO Hammerstein models using least squares support vector machines, Automatica 41 (7) (2005) 1263–1272. [7] F. Giri, E. Bai, Block-oriented Nonlinear System Identification, Springer, Berlin, Germany, 2010. [8] J. Paduart, L. Lauwers, J. Swevers, K. Smolders, J. Schoukens, R. Pintelon, Identification of nonlinear systems using polynomial nonlinear state space models, Automatica 46 (4) (2010) 647–656.
173
[9] D.T. Westwick, J. Schoukens, Initial estimates of the linear subsystems of Wiener–Hammerstein models, Automatica 48 (11) (2012) 2931–2936. [10] L. Vanbeylen, Nonlinear LFR block-oriented model: potential benefits and improved, user-friendly identification method, IEEE Trans. Instrum. Meas. 62 (12) (2013) 3374–3383. [11] B. Bequette, Challenges and recent progress in the development of a closed-loop artificial pancreas, Annu. Rev. Control 36 (2) (2012) 255–266. [12] A. Abu-Rmileh, W. Garcia-Gabin, Wiener sliding-mode control for artificial pancreas: a new nonlinear approach to glucose regulation, Comput. Methods Program Biomed. 107 (2) (2012) 327–340. [13] R. Hovorka, V. Canonico, L.J. Chassin, U. Haueter, M. Massi-Benedetti, M.O. Federici, T.R. Pieber, H.C. Schaller, L. Schaupp, T. Vering, et al., Nonlinear model predictive control of glucose concentration in subjects with type 1 diabetes, Physiol. Meas. 25 (4) (2004) 905–920. [14] C. Dalla Man, R.A. Rizza, C. Cobelli, Meal simulation model of the glucose–insulin system, IEEE Trans. Biomed. Eng. 54 (10) (2007) 1740–1749. [15] A. Marconato, M. Schoukens, K. Tiels, W. Widanage, A. Abu-Rmileh, J. Schoukens, Comparison of several data-driven nonlinear system identification methods on a glucoregulatory system example, IET Control Theory Applicat. (2013), under review. [16] H. Lee, B.W. Bequette, A closed-loop artificial pancreas based on model predictive control: human-friendly identification and automatic meal disturbance rejection, Biomed. Signal Process. Control 4 (4) (2009) 347–354. [17] J. Schoukens, R. Pintelon, T. Dobrowiecki, Y. Rolain, Identification of linear systems with nonlinear distortions, Automatica 41 (3) (2005) 491–504. [18] J. Potter, Matrix quadratic solutions, SIAM J. Appl. Math. 14 (3) (1966) 496–501. [19] K. Mårtensson, On the matrix Riccati equation, Inform. Sci. 3 (1) (1971) 17–49.
Contents lists available at ScienceDirect
Biomedical Signal Processing and Control journal homepage: www.elsevier.com/locate/bspc
Identification of a nonlinear model for a glucoregulatory benchmark problem Laurent Vanbeylen a,∗ , Anne Van Mulders a , Amjad Abu-Rmileh b a b
Department ELEC, Vrije Universiteit Brussel, Belgium Department of Brain and Cognitive Sciences, Ben-Gurion University of the Negev, Israel
a r t i c l e
i n f o
Article history: Received 15 January 2014 Received in revised form 1 April 2014 Accepted 24 April 2014 Keywords: Insulin-glucose modelling Artificial pancreas System identification Black-box modelling Best linear approximation Nonlinear models Nonlinear fractional representation
a b s t r a c t Recently, a novel identification method for a nonlinear dynamic model, called nonlinear Linear Fractional Representation (NL-LFR) model, has been developed. The model, composed of a static nonlinearity (SNL) surrounded by linear dynamics, can account for both nonlinear feed-forward and nonlinear feedback effects. Using two classical frequency response measurements, the SNL is automatically recovered in a user-friendly and efficient (non-iterative) way. In this contribution, the method is illustrated on a glucoregulatory benchmark dataset (insulin–glucose relationship of the human body). The research on insulin–glucose models is essential to develop methodologies to control the blood glucose level in diabetes patients. The obtained results outperform earlier results on the same benchmark data, while providing an excellent accuracy-complexity tradeoff. © 2014 Elsevier Ltd. All rights reserved.
1. Introduction To be able to describe measured data with an increased accuracy, e.g. over a wide amplitude range, it is important to explicitly incorporate nonlinearities into the model. Until now, a range of nonlinear models have been proposed [1–9]. However, most of them either require a (typically) quite large number of parameters (i.e., they are not parsimonious enough) or are not well suited to represent ubiquitous nonlinear feedback effects, such as amplitude-dependent resonance frequencies and dampings (i.e., they are not flexible enough). The nonlinear Linear Fractional Representation (NL-LFR) model belongs to the popular class of nonlinear block-oriented models. It offers a good flexibility-parsimony tradeoff, and is therefore applicable to a wide range of scientific domains. Block-oriented nonlinear models [1,7] can be defined as interconnections of linear dynamic elements (a.k.a. Linear Time Invariant parts or LTI blocks) and static nonlinear elements (a.k.a. Static NonLinearities or SNL blocks). The most studied block-oriented model structures are cascade structures of the form LTI-SNL, SNL-LTI and LTI-SNLLTI (known as Wiener, Hammerstein and Wiener–Hammerstein model, respectively).
∗ Corresponding author. Tel.: +32 26292979. E-mail addresses: [email protected] (L. Vanbeylen), [email protected] (A. Van Mulders), [email protected] (A. Abu-Rmileh). http://dx.doi.org/10.1016/j.bspc.2014.04.007 1746-8094/© 2014 Elsevier Ltd. All rights reserved.
The NL-LFR model, to be discussed in Section 3, is more general than these cascades and surrounds the SNL by arbitrary dynamics. We refer to [10] for a complete description of these matters and of the benefits of the model. In this paper, an insulin–glucose benchmark problem (see Section 2) is presented as an application of the NL-LFR model. The results and comparison to earlier results with different model structures are shown in Sections 4–6. Section 7 concludes the paper. 2. Insulin–glucose modelling problem 2.1. Aim To enhance the lives of (type 1) diabetes (mellitus) patients, research is conducted on the regulation of the blood glucose level [11]. E.g., glucose control schemes for an artificial pancreas have been developed for the insulin–glucose metabolism [12], based on a Wiener (LTI-SNL) model and nonlinear sliding mode control. More sophisticated simulation insulin–glucose models built from first principles exist [13,14], but are avoided for control due to a large computational effort involved with the implementation and model parameters are very difficult to estimate from the available patient data. In this work, as in [15], simplified nonlinear models, for use in control in a later stage, represent the system as a replacement for the sophisticated models. To estimate a simplified nonlinear model, a data-driven, measurement point of view is taken here: the insulin is treated as a measured input and the glucose concentration
L. Vanbeylen et al. / Biomedical Signal Processing and Control 13 (2014) 168–173
Fig. 1. (Concatenated) input and output signals for estimation (dataset A): randomphase multisines at 12 operating points.
as a measured output of the glycoregulatory system. In this paper, estimation results for the NL-LFR model are presented, and are compared with the Wiener and other nonlinear models. The data from which the simplified nonlinear models are estimated, are priorly generated through the Dalla-Man model [14]. The same data as in [15] are used: partly with a multisine1 insulin input (dataset A) and partly a multisine input with band-limited pulses superimposed, intended to simulate a more realistic situation (dataset B). The data are available at 5 multisine amplitudes and 12 DC-levels. Please note that the estimation of a (simplified) nonlinear model is based on simulated data. The simple reason is that performing real measurements on patients can be very expensive, and that high-quality results on simulated data may be required to convince the medical experts that performing a real measurement campaign makes sense. 2.2. Benchmark datasets The benchmark datasets used consist of steady-state2 measurements of the glucose response of the Dalla-Man model to an insulin input at different basal (DC) glucose levels, while keeping the meal input zero [15]. Note that, for both datasets, only the data at a single input std3 -value is taken as estimation data and another (slightly lower) std-value is taken as validation data, only intended to assess the model quality. The input (insulin4 ) and output (glucose) signals for datasets A and B are depicted in Figs. 1 and 2, respectively. The properties of both datasets are listed in Table 1. To compare the results, the following measure for the fit quality is used on the validation data:
Fit% =
y − yˆ 2 1− ¯ 2 y − y
169
Fig. 2. (Concatenated) input and output signals for estimation (dataset B): bandlimited pulses superimposed on a random-phase multisine at 12 operating points.
measured output (constant signal), respectively. Distinct values for the fit quality for the different operating points will be considered.
3. Nonlinear LFR model The focus of this paper is to identify an NL-LFR model on the insulin–glucose benchmark data. The model’s input and output are represented by u and y, as shown in Fig. 3. The SNL has an input z and output v. It has no memory effects, and can therefore be considered as a simple nonlinear mapping fSNL : R → R operating on each input sample z(k) (with discrete time index k):
v(k) = fSNL (z(k))
(1)
The linear part (LTI) connects the SNL and the model’s input and output: the MIMO LTI converts the two inputs, u and v, to two outputs, y and z. Signals v and z are internal and not accessible for measurement. The MIMO-LTI part can, in general, be described as: x(k + 1) = Ax(k) + Bu(k)
(2)
y(k) = Cx(k) + Du(k)
× 100
with y, yˆ , and y¯ the steady-state measured output signal, the steadystate simulated model output signal and the mean value of the
1 The use of multisine signals (among other types of variations), to represent low amplitude slow fluctuations around a basal insulin level is known in the field of artificial pancreas, mainly for the identification of control-relevant models [16]. The use of such signals taking into account patients safety limits, would be possible. 2 In the context of this paper, “steady-state” means “after die-out of the transient phenomena”, i.e. the periodic regime response is considered. 3 std = standard deviation 4 It can be remarked that, at some places, negative insulin values are visible. This is not realistic. However, the idea behind using negative values is to implement fluctuations below the basal (viz., average) insulin levels.
Fig. 3. The NL-LFR model. The signals u, y, z and v represent the model input and output and the SNL’s input and output, respectively (the identification method has no access to measurements of z and v).
170
L. Vanbeylen et al. / Biomedical Signal Processing and Control 13 (2014) 168–173
Table 1 Properties of both datasets.
Kind of excitation Sampling period Ts Number of data samples (per period) Number of periods Bandwidth of the input [1/min] Input std range [pmol/min] Input std for estimation [pmol/min] Input std for validation [pmol/min] Input o/pa range [pmol/min] Total number of experiments a
Dataset A
Dataset B
Odd random-phase multisine 20 min 2000 2 0.01 26–176 105 60 115.0–539.0 5 std-values × 12 DC-values
Band-limited pulses superimposed on an odd random-phase multisine 1 min 40,000 2 0.33 115–316 212 177 154.0–539.0 5 std-values × 12 DC-values
o/p = operating point (DC-offset) = mean value = basal level.
with states x(k) ∈ Rn (order n), input vector u(k) =
vector y(k) =
D=
Dyu Dzu
Dyv Dzv
y(k) z(k)
. The matrices B =
Bu
u(k) v(k)
Bv , C =
, output Cy Cz
4.1. Brief summary and
are splitted in subblocks. This results into:
x(k + 1) = Ax(k) + Bu u(k) + Bv v(k) y(k) = Cy x(k) + Dyu u(k) + Dyv v(k)
(3)
z(k) = Cz x(k) + Dzu u(k) Note that, as in [10], the coefficient Dzv , is set to zero to avoid algebraic loops when performing the simulation (computation of the response for a given signal u, SNL and state-space matrices). Discussion. Through the interconnection of the SNL (1) and the MIMO-LTI (3), the NL-LFR model features both nonlinearities and dynamics. Moreover, the dynamics surround the SNL in a flexible way, such that, in principle, any finite-dimensional system with a localised SNL can be captured by an NL-LFR model. Besides, please note that an alternative, equivalent representation for the model, shown in Fig. 4 exists, where the MIMO-LTI is split in four individual contributions. The transfer function matrix u y of the MIMO-LTI from inputs u = to outputs y = , conz v taining each of the four contributions, is connected as follows to the state-space representation:
Gyu (Z) =
GD (Z)
GH (Z)
GW (Z)
GFB (Z)
= D + C(ZI − A)−1 B.
4. Dataset A: identification methodology and results
(4)
where Z is the z-transform variable. The subscripts of the four contributions are associated with their positions w.r.t. the SNL in the block-scheme of Fig. 4: GD denotes the direct linear contribution from u to y, this signal path is not affected by the SNL; GFB denotes the feed-back which forms a nonlinear feedback loop together with the SNL; this loop is preceded and followed by the dynamics GW and GH , respectively. This corresponds to the positions of the linear dynamics of a Wiener and a Hammerstein system w.r.t. the SNL (viz., before and after, respectively).
Fig. 4. The NL-LFR model (LTI decomposed in four contributions). Note that the output y is here shown in the bottom left corner.
The identification of the NL-LFR model presented in the previous section can be carried out by following the lines of [10]. In a nutshell, the approach is based on data at different input conditions (e.g., amplitudes or operating points), exciting each time the nonlinearity differently. From each part of the data the Best Linear Approximation (BLA) is extracted. By solving a Riccati equation, the MIMO-LTI part (viz., the common linear dynamics in the NL-LFR) can be obtained. In a last step, the internal signals z(k) and v(k) are reconstructed, from which the SNL is estimated. The MIMO-LTI and the SNL together constitute the initial nonlinear model. 4.2. Detailed description The NL-LFR identification method starts with the identification of the BLA, the linear dynamics optimal in mean-square sense [17] for two different input conditions. Here, the data at input operating points 2 and 5 are selected. First, the two BLAs are estimated nonparametrically (see Fig. 5 for their Frequency Response Function (FRF) estimates). Next, each FRF is fitted as a parametric model, more precisely a 3rd order linear state-space model. This fit is obtained by means of a frequency-domain subspace technique. It is shown in Fig. 6. Given the two state-space models, a Nonsymmetric Algebraic Riccati Equation (NARE) is solved. This can be done very fast (see [18,19]) and in general leads to several possible solutions of the MIMO dynamics (here, there were 4 solutions). The input and output of the SNL are reconstructed for each of these solutions in the frequency domain (the data from all operating points are used). Next, a scatter plot between these reconstructed signals is generated (see Fig. 7), and the linearly parameterized SNL is fitted in one least-squares step (here, a polynomial of degree 4). The best fit
Fig. 5. Nonparametric FRFs for operating points 2 (grey dots) and 5 (black dots).
L. Vanbeylen et al. / Biomedical Signal Processing and Control 13 (2014) 168–173
171
Fig. 6. Nonparametric FRFs (black dots) and their parametric fit (black line) with the difference between both (grey dots). Top plot: operating point 2; bottom plot: operating point 5.
Fig. 8. 4th degree polynomial fit on the SNL for which the residual of the fit was the lowest. Top: nonparametric SNL (dots) and polynomial fit (black line); bottom: residuals of the fit.
Fig. 7. The nonparametric SNLs corresponding to the 4 possible solutions of the NARE. Light grey to black for operating points 1–12.
(best out of the 4 solutions shown in Fig. 7) of the SNL is given in Fig. 8. This solution (MIMO dynamics and SNL) is retained and used as an initial value in a further nonlinear-least-squares optimisation of the model parameters (matrix entries in (3) and polynomial coefficients). The FRFs of GD , GW , GH and GFB (corresponding to the final model) are given in Fig. 9.
Fig. 9. Four dynamic parts of the estimated MIMO-LTI block.
5. Dataset A: validation and comparison This section describes the validation result and compares it with other identification approaches (different model types). Fig. 10 shows the output signal of the validation data and the errors made by the two linear models used in the identification method. It can be seen that the linear models perform very well at the operating point they were fitted on, but fail completely outside this region.5 Moreover, the DC-errors (mean error values) were removed in these plots. Hence, the actual results are worse. The result of the NL-LFR model is shown in Fig. 11. Note already the difference in the order of magnitude of the error signals (see vertical axis): a relative error of the order of a few percent can be expected. The nonlinear model performs well, and this on all
5
Note also that a linear model fitted on the entire dataset would not produce good results: it would be poor everywhere.
Fig. 10. Dataset A. Top: (concatenated) output signal of the validation data; middle: error of the BLA constructed at operating point 2; bottom: error of the BLA constructed at operating point 5. Note that DC-errors are removed in this plot.
172
L. Vanbeylen et al. / Biomedical Signal Processing and Control 13 (2014) 168–173
Fig. 11. Dataset A. Top: (concatenated) output signal of the validation data; bottom: error of the NL-LFR model. No DC-correction was applied.
operating points. DC-levels are estimated automatically with the model: the error levels have no DC-correction. The fit-percentages for each of the 12 operating points of the models (linear and nonlinear) are reported in Table 2. In contrast to the linear models, the NL-LFR model can describe the entire operating region with a relatively high fit, in the range from 89.8% to 97.3%. It can be observed that the fit is not increasing for each operating point during the optimisation. This is because the optimisation balances the errors over the different set points due to the absence of weighting in the least-squares objective function optimised. The results of other nonlinear models are given in [15]. A schematic overview of these results, including the more recent result of this paper, is given in Fig. 12. The NL-LFR model presented in this paper clearly has a very good flexibility-parsimony tradeoff and outperforms the other results. 6. Dataset B: results and comparison In [15], the identification was redone for a (slightly) more realistic input signal (see Fig. 2), in which band-limited pulses are applied, superimposed by a multisine signal. This should mimic insulin shots administered to the patient. Without going into the details of the identification or results, this section provides a schematic overview of the results (see Fig. 13). The new dataset contained so many data points, that the nonlinear state-space (NLSS) models could no longer be obtained. The method presented in this paper did not suffer from this. It took only a few minutes to obtain an initial estimate with average fit 87.5%, Table 2 Dataset A. Fit% for each operating point. First two columns: results of the two BLAs; last two columns: results of the NL-LFR model before and after optimisation. The fits higher than 80% are indicated in boldface. o/p 1 2 3 4 5 6 7 8 9 10 11 12
BLA o/p 2
BLA o/p 5
NL-LFR init
NL-LFR
74.0 94.0 64.2 22.1 −10.3 −70.5 −141.0 −219.3 −303.4 −393.0 −485.1 −579.1
34.2 46.6 63.1 83.3 94.3 71.2 37.7 0.2 −40.1 −83.1 −127.4 −172.7
95.6 95.3 92.9 93.0 91.5 98.2 97.7 96.4 94.0 91.1 87.7 83.8
96.4 97.3 95.7 91.9 91.7 96.4 95.0 96.5 96.7 92.6 91.5 89.8
Fig. 12. Dataset A. Fit% as function of the number of independent model parameters (model complexity). The average fit (over all operating points) is indicated by a dot and the range from minimal to maximal fit quality (depending on the operating point) is given by a box. NN, (P)NLSS and SA denote neural network, (polynomial) nonlinear state-space, and state-affine model, respectively. The BLA depicted here corresponds to the operating point 7, located close to the center.
Fig. 13. Dataset B. Fit% as function of the number of independent model parameters. The same notation as in Fig. 12 is used.
minimal fit 78.6% and maximal fit 91.2%. Starting from this fit, only four iterations were performed, but the result is already quite good compared to the other results. The MIMO dynamic part of the NLLFR model has order n = 4 and the SNL has degree d = 3. 7. Conclusion In this paper, the aim was to obtain a simplified representation for a benchmark dataset representing the human glycoregulatory system (insulin–glucose behaviour), priorly simulated by means of the so-called Dalla-Man model. High-quality fits, outperforming earlier results, have been obtained with the nonlinear LFR model, consisting of a MIMO-LTI part and a polynomial SNL interconnected with each other. An excellent accuracy-complexity tradeoff is achieved with fit percentages of 90–97% and only a couple of tens of parameters. The model’s high quality can be attributed to its flexibility: the SNL is surrounded by arbitrary linear dynamics, taking both nonlinear feedforward and nonlinear feedback effects into account (which is not the case for Wiener-type models). The NLLFR model combines efficiently the power of the linear framework and the simplicity of block-oriented nonlinear structures. Acknowledgments Laurent Vanbeylen and Anne Van Mulders’s work was supported by the Fund for Scientific Research (FWO-Vlaanderen), the Flemish Government (Methusalem Grant METH-1), the Belgian Program on Inter-university Poles of Attraction (IAP VII/19 – Dysco), and
L. Vanbeylen et al. / Biomedical Signal Processing and Control 13 (2014) 168–173
by the ERC advanced grant SNLSID, under contract 320378. Amjad Abu-Rmileh acknowledges the support of the University of Girona through the BR-UdG research grant. References [1] S. Billings, S. Fakhouri, Identification of systems containing linear dynamic and static nonlinear elements, Automatica 18 (1) (1982) 15–26. [2] S. Boyd, L. Chua, Fading memory and the problem of approximating nonlinear operators with Volterra series, IEEE Trans. Circuits Syst. 32 (11) (1985) 1150–1161. [3] I. Hunter, M. Korenberg, The identification of nonlinear biological systems: Wiener and Hammerstein cascade models, Biol. Cybern. 55 (2) (1986) 135– 144. [4] S. Chen, S. Billings, Representation of non-linear systems: the NARMAX model, Int. J. Control 49 (3) (1989) 1012–1032. [5] G. Giannakis, E. Serpedin, A bibliography on nonlinear system identification, Signal Process. 81 (3) (2001) 533–580. [6] I. Goethals, K. Pelckmans, J.A. Suykens, B. De Moor, Identification of MIMO Hammerstein models using least squares support vector machines, Automatica 41 (7) (2005) 1263–1272. [7] F. Giri, E. Bai, Block-oriented Nonlinear System Identification, Springer, Berlin, Germany, 2010. [8] J. Paduart, L. Lauwers, J. Swevers, K. Smolders, J. Schoukens, R. Pintelon, Identification of nonlinear systems using polynomial nonlinear state space models, Automatica 46 (4) (2010) 647–656.
173
[9] D.T. Westwick, J. Schoukens, Initial estimates of the linear subsystems of Wiener–Hammerstein models, Automatica 48 (11) (2012) 2931–2936. [10] L. Vanbeylen, Nonlinear LFR block-oriented model: potential benefits and improved, user-friendly identification method, IEEE Trans. Instrum. Meas. 62 (12) (2013) 3374–3383. [11] B. Bequette, Challenges and recent progress in the development of a closed-loop artificial pancreas, Annu. Rev. Control 36 (2) (2012) 255–266. [12] A. Abu-Rmileh, W. Garcia-Gabin, Wiener sliding-mode control for artificial pancreas: a new nonlinear approach to glucose regulation, Comput. Methods Program Biomed. 107 (2) (2012) 327–340. [13] R. Hovorka, V. Canonico, L.J. Chassin, U. Haueter, M. Massi-Benedetti, M.O. Federici, T.R. Pieber, H.C. Schaller, L. Schaupp, T. Vering, et al., Nonlinear model predictive control of glucose concentration in subjects with type 1 diabetes, Physiol. Meas. 25 (4) (2004) 905–920. [14] C. Dalla Man, R.A. Rizza, C. Cobelli, Meal simulation model of the glucose–insulin system, IEEE Trans. Biomed. Eng. 54 (10) (2007) 1740–1749. [15] A. Marconato, M. Schoukens, K. Tiels, W. Widanage, A. Abu-Rmileh, J. Schoukens, Comparison of several data-driven nonlinear system identification methods on a glucoregulatory system example, IET Control Theory Applicat. (2013), under review. [16] H. Lee, B.W. Bequette, A closed-loop artificial pancreas based on model predictive control: human-friendly identification and automatic meal disturbance rejection, Biomed. Signal Process. Control 4 (4) (2009) 347–354. [17] J. Schoukens, R. Pintelon, T. Dobrowiecki, Y. Rolain, Identification of linear systems with nonlinear distortions, Automatica 41 (3) (2005) 491–504. [18] J. Potter, Matrix quadratic solutions, SIAM J. Appl. Math. 14 (3) (1966) 496–501. [19] K. Mårtensson, On the matrix Riccati equation, Inform. Sci. 3 (1) (1971) 17–49.