Predicting the output error of the suboptimal state estimator to improve the performance of the MPC-based artificial pancreas

The error of single step-ahead output prediction is the information traditionally used to correct the state estimate while exploiting the new measurement of the system output. However, its dynamics and statistical properties can be further studied and exploited in other ways. It is known that in the case of suboptimal state estimation, this output prediction error forms a correlated sequence, hence it can be effectively predicted in real time. Such a suboptimal scenario is typical in applications where the process noise model is not known or it is uncertain. Therefore, the paper deals with the problems of analytical and empirical modeling, identification, and prediction of the output error of the suboptimal state estimator for the sake of improving the output prediction accuracy and ultimately the performance of the model predictive control. The improvements are validated on an empirical model of type 1 diabetes within an in-silico experiment focused on glycemia prediction and implementation of the MPC-based artificial pancreas.


Introduction
In many applications and problems, like the model predictive control of glycemia in diabetic subjects, there is a need to estimate the system state based solely on noisy measurements of the system output.Traditional recursive state estimators/observers, such as the Kalman filter [1], are algorithms based on the prediction and correction of the state estimate by the new output measurement.It means that the state estimate is corrected (innovated) according to the single step-ahead output prediction error produced by the system model.
It is known from the theory of Kalman filtering, that for an optimal state estimator, the sequence of the output error of the state observer (OESO), which is also called the innovation sequence, has the properties of Gaussian white noise.However, in the case of suboptimal state estimation, this innovation sequence is correlated [2,3] and thus can be effectively predicted in real time.
However, why should one even consider using the suboptimal state estimation instead of the optimal one?To answer this question, the suboptimal state estimation is basically not an option of choice, but rather an inevitable consequence that arises from the unknown or poorly estimated parameters of the process noise model.The assumption of suboptimality is because the design of Kalman filter and its optimality relies on the exact knowledge of the process noise model, parameters of which are highly uncertain in many applications, or even the entries of the process noise covariance matrix are often used simply as the tuning parameters.It can be claimed that if there is a mismatch between the noise model and the actual statistical properties of the system, the OESO forms a correlated sequence.
The main practical motivation for studying the dynamics of the OESO is to try to predict it in real time and then correct the predictions of the output variable accordingly, yet the ultimate aim is to involve this prediction in the model predictive control.Therefore, the outlined strategy can be seen as a relatively feasible and cheap way to improve prediction and control performance by effectively compensating for the effect of suboptimal state estimation.
Concerning the target application domain outlined in the title, performing highly accurate predictions yields better forecasting of severe hyper-and hypo-glycemia states, as these are the major risks linked to diabetes and its treatment [4,5].Another application of the proposed strategy is possible within an implementation of the model predictive control-based artificial pancreas [6,7] to control glycemia in subjects with type 1 diabetes by automating the insulin dosing.
The rationale for choosing this application domain to demonstrate the effectiveness of the proposed strategy is supported by the fact that in most available studies, e.g. in [8][9][10][11][12][13][14][15][16][17], the crucial entries of the covariance matrix of the process noise are considered tuning parameters in the Kalman filter design.Therefore, it can be concluded that the state estimator works suboptimally in such scenario.
It should be mentioned for completeness, that there also exist sophisticated and dedicated methods [18][19][20][21][22] for estimating the noise models, which can potentially eliminate the problem of suboptimality of the state estimation.Unfortunately, these covariance matrix estimation methods have limited practical applicability, as they often provide biased estimates and typically require very large datasets to be supplied, what can be infeasible in many applications.
In this paper, we provide important insights into the dynamics of the OESO from the analytical point of view, but also from the practical data-driven perspective where reduced models are considered.The primary motivation for using the reduced structures, such as the autoregressive and the moving average model, rather than the full analytical model, is their feasible prediction and identifiability from the experimentally obtained OESO sequence.
The paper has been divided into the following sections: Sect. 2 introduces the basic preliminaries, formulation of the stochastic state-space model, and equations describing the conventional state observer.In Sect.3, the full analytical model of the OESO is derived to provide some theoretical background.Section 4 considers the reduced autoregressive model and comprises the formulation of the corresponding identification problem and the predictive equation.In Sect.5, the moving average model is studied in a similar way.Section 6 covers the idea of using the identified models to enhance the output prediction accuracy and their inclusion in the model predictive control.The setup of the experiment aimed at prediction and model predictive control of glycemia in subjects with type 1 diabetes is outlined in Sect.7, while its results are discussed in Sect.8.
It should be marginally mentioned that the stochastic nature of glycemia dynamics does not necessarily have to be modeled in the state space by the process noise and the measurement noise, since basic input-output transfer function models like the autoregressive-exogenous, autoregressive moving average with exogenous inputs, and Box-Jenkins model can also be used as reported in [23][24][25][26].The stochastic part of these input-output models is usually estimated in one step together with the deterministic submodels as a result of the system identification procedure from the experimental data, so there is no need to worry about the optimality of the state estimation.
An application of the unconstrained model predictive control with the state-space model along with the state estimation based on the Kalman filter was presented in study [27].A similar control strategy was reported in [28,29], but due to the input-output problem formulation the state estimator was not required.As another example of similar artificial pancreas scheme, in [30] the Kalman filter was used to estimate the state for the linear model predictive control with disturbance rejection.
In our recent work [31], a novel optimal state estimator was proposed as an alternative to the Kalman filter.However, this algorithm was not based on the traditional recursive correction of the state estimate by the OESO, but it used the generalized least squares formulation of the problem.Also in this case, the optimality of the state estimate depended on the exact knowledge of the process noise model.
From the perspective of unique contributions of this paper, it is important to remark that in any of the above referenced studies or in the latest comprehensive survey papers [32][33][34], the strategy of predicting the OESO and compensating for the suboptimality of the state estimation was not proposed or discussed.It can also be concluded that most authors relied on suboptimal state estimation with ad hoc tuning of the process noise model, hence leaving a significant headroom for improving the control performance by targeting this problem.
Unlike the conventional MPC-based artificial pancreas that typically uses the suboptimal state estimator, which normally results in a correlated OESO sequence, we propose a new strategy to effectively compensate for this effect.The proposed modification involves the prediction of the OESO to directly correct the prediction of the system free response.It is worth noting that this modification is easy to embed in already existing MPC schemes [32][33][34] while inducing little more computational cost yet requiring no additional hardware modifications.In other words, the proposed OESO models and the prediction/correction strategy can eventually be retrofitted to any advanced MPC-based artificial pancreas design that utilizes the state estimator at the expense of relatively straightforward structural modification.Note that the other features of the artificial pancreas, such as the safety algorithms and constraints, do not directly interact with the proposed modification.

Model structure and preliminaries
The general discrete-time stochastic state-space empirical model of glycemia dynamics in subjects with type 1 diabetes is postulated as [31] x where k ∈ N is the current sample, the output y [mmol/l] stands for the deviation of glycemia from its steady-state value, the input u [U/min] denotes the deviation of the insulin administration rate from the basal rate, and d [g/min] represents the carbohydrate intake rate input.The state vector of this nth order model is denoted x [n×1], w [n×1] is the process noise vector and the zero-mean uncorrelated random process v ∼ N (0, R) represents the measurement noise of the glucose monitoring device [35].
The parameters of model ( 1) include the state-transition matrix A [n×n], the input matrix B [n×2], and the output vector C [1×n].The state-transition matrix A consists of the submatrices A u , A d and the zero matrices 0 of the conforming dimensions as where matrices A u [n u ×n u ] and A d [n d ×n d ] are in the canonical form and comprise the model coefficients a u d such that The input matrix B is simply where 0 are the zero vectors of conforming dimensions and The output vector C gets where C u [1×n u ] and C d [1×n d ] will comprise the model coefficients c u d as The state vector x also holds the canonical form where n u and n d are the orders of the corresponding submodels, so the overall model order is n = n u +n d .The state variables x u , x d represent the partial effects of insulin administration and carbohydrate intake, respectively.Consider that the process noise w in model ( 1) represents the effect of input uncertainties, so one can write The first random input γ (k) ∼ N 0, σ 2 γ reflects various unmeasurable disturbances, including physiological changes in insulin absorption and action.The second random input δ (k) ∼ N 0, σ 2 δ represents the uncertainty of the meal announcing induced by the patient.
Since all stochastic terms in (1) were defined as zero-mean uncorrelated stationary random processes, the covariance matrix Q [n×n] of the process noise (9) and the variance R of the measurement noise are equal to

State observer
The state vector x of model ( 1) is usually estimated using the recursive state observer where x [n × 1] is the estimated state and K [n × 1] is the gain vector, which is the subject of the observer design.The design is usually based on the optimal Kalman approach [1,36] in the case of stochastic system assumption, or using the pole-placement method if the system is deterministic.
The state estimate residual e [n×1] will be defined as Finally, the single step-ahead model output prediction error , i.e. the OESO, gets

Dynamics of the state observer output error
The model of the state estimate residual e can be derived by substituting x(k+1) from ( 12) and x (k+1) according to (1a) into ( 13) while substituting the output y (k+1) in the terms of 1b as Eq. ( 15) can be transformed by applying the forward timeshift operator z as e (k+1) = ze (k) obtaining The above equation can be reshaped to separate vector e (k) as where I is the unit matrix of the conforming dimensions.
According to ( 14) and ( 17), (k) holds Eq. ( 18) implies that the dynamics of the OESO is represented by a stochastic system with multiple independent noise inputs.One may realize that term C (z I − A+ K C) −1 results in a row vector of rational functions p i (z) s(z) with the common denominator s(z) as the characteristic polynomial of this system.This consideration yields the transfer function model where If Q is the covariance matrix of the process noise according to (10), then w (k) can be written as where 0 [n×1] is the zero vector, η [n + 1×1] is the vector of uncorrelated noise inputs with the unit variance, i.e., cov(η, η) = I , can be replaced by where R is the variance of the measurement noise according to (11).Finally, model ( 19) can be generalized as the sum of n+1 autoregressive-moving-average (ARMA) models by substituting w (k) in the terms of ( 21) and v (k) from ( 23) as Since the process noise (9) has only two nonzero components and diagonal covariance matrix (10), general model ( 24) can be reduced to However, model ( 24) or ( 25) cannot be used to predict the OESO, since the random input vector η as well as the partial outputs r i (z) s(z) are unmeasurable in practice.Another important paradox concerning this analytical model is that since the covariance matrix of the process noise is considered unknown in the case of suboptimal state estimation, analytical model (24) simply cannot be determined.On the contrary, if the covariance matrix of the process noise is known, what implies that the state estimator works optimally, then model (24) will be known, but it will be unnecessary since the OESO is uncorrelated and hence it cannot be predicted.
Concerning the identification of full model ( 24) directly from experimental data, i.e. based on the OESO sequence, this would be hardly possible primarily due to its structure and the large number of parameters to be estimated.
The aforementioned issues with predictability and identifiability are the main motivation for further considering two reduced model structures, particularly the autoregressive and the moving average model.

Optimality test
To test whether the state estimator works optimally or not, the sample autocorrelation function of the OESO sequence has to be analyzed.In the case of optimal state estimation, this autocorrelation function should show a character similar to that of Dirac delta function.
Supposing a finite-length experiment with N samples, the autocorrelation function R (nT s ) is estimated as [38,39] where n ∈ Z satisfies n < N .

Autoregressive model
In this section, the dynamics of the OESO will be approximated by the single-input single-output autoregressive model defined as where η ∼ N 0, σ 2 η is a random process with the properties of zero-mean white noise.
The polynomial q(z) of this n q th order model gets The equivalent difference equation of model ( 27) holds The parameter vector q gets q = q 1 q 2 • • • q n q T .(30)

Identification strategy
According to difference equation ( 29), the corresponding linear regression system considering N available samples takes . . .
or using the shorthand notation where q is the parameter vector ( 30), H AR N ×n q is the regression matrix and [N ×1], η [N ×1] are vectors.The parameter vector q can be estimated as q in a straightforward way using the least squares method with the optimal parameter estimate determined analytically as [40]

Predictive form
For model (27), the explicit prediction formula can be derived based on difference equation (29).The future values of the white noise input are obviously unknown, so assuming that its statistically unbiased prediction is zero, i.e.E η (k+i) = 0, the predictive form considering the prediction horizon n p gets where vectors p n q ×1 and ˆ f n p ×1 are defined as and the matrices M f n p ×n p , M p n p ×n q comprise the elements of vector q (30) as Since the noise input in ( 27) is unmeasurable, by reshaping equation (29), η (k) can be estimated as

Moving average model
In this section, the dynamics of the OESO will be approximated by the moving average model where η ∼ N 0, σ 2 η is the zero-mean white noise input.The polynomial g(z) in the n g th order model (40) gets The difference equation for model (40) can be written as The parameter vector g n g ×1 g i gets

Identification strategy
It is well known that estimating the moving average processes is more difficult than estimating the autoregressive processes [41].Since the input noise η is unmeasurable in practice, the straightforward approach based on the least squares minimization of the model single step-ahead prediction error cannot be directly applied in this case.Therefore, to estimate the coefficient vector (43) using the available OESO sequence, the two-step method of Durbin [41,42] will be adopted.The first step of this method consists of fitting an autoregressive model to the OESO sequence via the ordinary least squares method in the terms of Sect. 4. In the second step, the identified autoregressive model is used to estimate the input noise η by filtering the OESO sequence by the inverse of the estimated autoregressive model according to (39).
The second step uses this estimated input noise sequence η to create the regression system and to estimate the parameters of the moving average process in the least squares sense.
The corresponding regression system then gets . . .
or using the shorthand notation, where g is the parameter vector ( 43), H M A N ×n g is the regression matrix and [N ×1], η [N ×1] are vectors.
The optimal parameter vector g can be estimated as

Predictive form
Having the model parameters estimated, the OESO (36) can be predicted.Assuming that the statistically unbiased prediction of the input zero-mean white noise is zero, the predictive form of the moving average model ( 40) can be derived according to the difference equation (42) as where matrix M η p n p ×n g is formed by the elements of vector g (43) such that and vector ηp n g ×1 comprises the estimated past values of the noise input In practice, the input noise η cannot be measured, so it has to be estimated based on the inverse filtering of according to the difference equation (42) as

Prediction and model predictive control with the OESO compensation
In this section, the algorithm of model predictive control will be adopted from [31], while an important modification that concerns the prediction of the OESO will be proposed.Prediction of the state vector x and the output y can be expressed considering (1) as where k ∈ N is the current sample and i ∈ N gets i = 1 • • • n p , while assuming n p ∈ N is the prediction horizon.Notice that in (52), the OESO ˆ , which was predicted by the identified autoregressive or the moving average model, was taken into account by correcting the output prediction ŷ.This is the most important modification of the traditional prediction and predictive control algorithms as it allows us to effectively compensate for the suboptimality of the state estimation.
The predictive control minimizes the quadratic cost function of the model-based predictions of chosen system variables over the prediction horizon n p .The corresponding quadratic form gets [43,44] where u f [n c ×1] is the vector of future changes of the manipulated variable, while assuming n c is the control horizon.Matrix A [n c ×n c ], vector b [n c ×1] and scalar c are defined as where y r n p ×1 is the reference vector, ŷfree n p ×1 is the system free response vector, and scalar λ u denotes the weight of the manipulated variable increments penalty.
The free response prediction ŷfree in (54) gets Matrices B u n p ×n p and B d n p ×n p get the lower triangular form where matrices were defined by ( 3), ( 5), ( 7), respectively.
The forced response matrix H f n p ×n c in (54) gets The optimization problem ( 53) can be solved by quadratic programming if linear inequalities constraints are considered.For the sake of simplicity, the elements of the reference vector y r are equal to an appropriately chosen constant G t representing the target glycemia.Pursuing the receding horizon strategy, only the first element of the optimal solution u f is actually applied, so one can write Note that the manipulated variable has to be constrained, so the minimal insulin infusion rate u min = 0 U/min, while u max will be adopted from [27].The corresponding linear inequalities system with respect to the decision vector u f can be formed by involving matrix Ψ (57) as in [45] Compared to [31], the crucial modification of the control algorithm is made here by adding the prediction of the OESO ˆ (k+i) to equation (52).
Concerning the safety features of automated insulin therapy, various additional strategies could be considered to enhance the current configuration.To avoid the adverse and dangerous insulin stacking phenomenon, the insulin on board [46][47][48] representing the amount of insulin still active from the previously administered doses can be involved.The dynamics of insulin on board can be represented by simple linear models like in [49], [50] or [51], while this signal can be used to form a special quadratic penalty that is added to the cost function ( 53) of the MPC.Alternatively, hard constraints for the insulin on board signal can also be assumed by extending the linear inequalities system (59) of the quadratic program accordingly.Another type of safety feature is to hard constrain the controlled variable by modifying the inequalities system (59) to prevent the risk of hypoglycemia and hyperglycemia as proposed in [26].As this strategy can potentially induce control infeasibility, soft formulation of the controlled variable constraints should be preferred, as proposed in [52,53].However, further details are beyond the scope of this paper.For information on safety features, see, e.g.[46,48,54].Also note that these safety features do not directly interact with the proposed strategy of predicting the output error of the suboptimal state estimator, which is the main contribution of the paper.

Experimental setup
To validate the proposed strategy and assess its practical effectiveness in application to the problem of prediction and predictive control of glycemia in subjects with type 1 diabetes, a simulation-based experiment was designed and evaluated.
The glycemia response for this experiment was obtained by in-silico approach, simulating the complex physiologybased nonlinear model that was described in [55,56]  i in (7), is not within the scope of this paper, so we suppose that model (1) was identified with parameters (60).However, for more details on this topic, we refer an interested reader to our recent works [25,57].
C d = 0.0320 0.0061 0.0003 0 .( The prediction horizon and the control horizon were assumed as n p = 15, n c = 10, while the sample time was chosen as T s = 10 min.The variance of the measurement noise 11 and the variances of the process noise 10 were empirically adjusted as Note that since the variances of the process noise were just empirically tuned to obtain acceptable performance of the Kalman filter while the actual values cannot be determined because they are not based on any particular physiological mechanisms or characteristics of a diabetic subject, the state estimator will perform only suboptimally in this case, and hence the OESO sequence will be correlated.The observer gain vector K was calculated according to the Kalman filter design [1,36]  Concerning the initial tuning of the proposed empirical models of the OESO, the order of autoregressive model (27) was set as n q = 4, whereas the order of moving average model (40) was chosen as n g = 12.
The experiments were designed to mimic the insulin treatment of a subject with type 1 diabetes during the two-day period.The first investigated problem concerns the prediction of glycemia during standard insulin bolus therapy that was carried out according to the bolus calculator rule (see [58] and the references therein).The second deals with automated insulin dosing managed by the model predictive control algorithm of the artificial pancreas.Both algorithms were modified in the terms of Sect.6.

Discussion
In this section, the results of the outlined simulation experiment will be comprehensively analyzed and discussed.The sequence of the OESO obtained in the terms of equation ( 14) during regular insulin treatment while simultaneously performing the state estimation according to (12) by considering the observer gain (62) is plotted in Fig. 1.This figure suggests that although the state estimate asymptotically converges to the actual state, the character of the OESO sequence is far from ideal uncorrelated noise.
To prove this, the autocorrelation function R of the OESO was estimated according to (26) by processing the sequence from Fig. 1 and is plotted in Fig. 2. Analyzing this autocorrelation function, one can conclude that the OESO sequence is correlated, what confirms that the state estimator works suboptimally due to the empirically adjusted variances of the process noise (61).
Performing the estimation of both reduced models of the OESO by pursuing the strategies presented in Sections 4 and 5, the following coefficients were estimated.
Parameter vector q of the autoregressive model ( 27):  Fig. 3 Estimate of the autocorrelation function R η η(nT s ) of the estimated noise input for the autoregressive model ical models, i.e. by q(z) and 1 g(z) , respectively, the input noise sequence was estimated according to (39) for the autoregressive model and according to (50) for the moving average model.The autocorrelation functions R η η(nT s ) of these sequences for both model structures are depicted in Figs. 3 and 4, which show their Dirac delta function-like nature, proving the estimated models valid.Now follows the prediction of the OESO using the predictive form (34) of the autoregressive model and the predictive form (47) of the moving average model respectively.Relatively accurate predictions for randomly chosen starting points can be observed in Fig. 5.Both model structures showed almost identical performances and could predict the future evolution of the correlated OESO with a satisfying accuracy considering the highly stochastic nature of this signal.
The next comparison concerns the practical impact of correcting the glycemia prediction by predicted OESO as the original contribution of the paper.In Fig. 6, one can see the uncorrected prediction of glycemia ( Ĝ) as the conventional strategy, as well as the predictions that involved the corrections by OESO predicted using the autoregressive ( ĜAR ) and the moving average ( ĜMA ) model.By a basic visual assessment, one can observe an improvement of the prediction accuracy, the improvements concerned primarily the peaks of the response.Keep in mind that such differences between the uncorrected and the corrected prediction can be critical in situations such as decision making with regard to the application of insulin therapy.
In addition to the graphical assessment, the prediction performance will be quantified by the quadratic metric which will provide a better assessment of the prediction performance.
The last part of the experiment is focused on the model predictive control of glycemia in the context of the artificial pancreas implementation, where a positive effect of the proposed predictors on control performance is anticipated.To demonstrate this, Fig. 7 shows the closed loop glycemia response, where involving the predictions of the OESO visibly improved the control performance in terms of tighter control with respect to the reference value, and reduced maximal and minimal observed glycemia, what is especially significant to reduce the risk of hyperglycemia and hypo- Fig. 7 Predictive control of glycemia without the OESO compensation G(t) compared to using the autoregressive model G AR (t) and the moving average model G MA (t) glycemia.It can be concluded that both predictors performed almost identically, but way better than in the original case without compensating for the OESO.Keep in mind that since typical values of the OESO are relatively low compared to the magnitude of the controlled variable (see Fig. 1), the proposed strategy can naturally yield a limited effect.It can also be claimed that the strength of the desired effect is directly related to the performance level of the state observer and thus to the degree of mismatch between the process noise model and the actual statistical properties of the system.Therefore, for systems with an empirically tuned covariance matrix of the process noise, the strategy proposed in this paper is highly recommended.
The control performance will be quantified by the maximal G max and the minimal G min observed glycemia, and by the quadratic metric where G t is the target glycemia.The summary of prediction and control performance metrics obtained during the experiment is documented in Table 1, which also confirm the observations from Figs. 6 and 7.
To investigate the effect of the choice of the tunable parameters of the OESO models, particularly the order of the autoregressive model ( 27) n q and the order of the moving average model ( 40) n g on the resulting performance of the proposed strategy, the experimentation was repeated under various configurations, yielding the results summarized in Tables 2 and 3.
It can be concluded that the moving average model of the OESO performs slightly better in both prediction and predictive control than the autoregressive model, while the choice of the corresponding model orders n q , n g also affected the overall performance, yet not consistently for all the metrics considered.However, all studied models performed better than the original uncompensated configuration.

Conclusions
This study stressed that the OESO sequence in the case of suboptimal state estimation is correlated, while it can be effectively predicted by the autoregressive or the moving average models.These two reduced models could provide good predictability and identifiability from the experimental data.The predicted OESO was then involved to correct the output variable prediction and hence ultimately improve the performance of the model predictive control.It can be concluded that the presented strategy allowed to effectively compensate for the suboptimality of the state estimation caused by the inaccuracy of the process noise model in a relatively inexpensive and feasible way.
We also obtained theoretical results demonstrating that the actual dynamics of the OESO is analytically described as the sum of ARMA processes, while this full structure was approximated by the autoregressive and the moving average models in practice.
A promising application of our results would be possible within an implementation of the artificial pancreas in subjects with type 1 diabetes, where maximizing the performance of the model predictive control is of highest priority.The presented simulation-based experiment demonstrated that the dynamics of the OESO can be effectively predicted by both proposed reduced models, while there were documented positive effects on the accuracy of the glycemia prediction and on the performance of the predictive control.Keep in mind that although qualitative improvements do not appear to be significant at first sight, any feasible improvement matters a lot when managing glycemia and can have significant longterm consequences for patient health.
Compared with the conservative structure of the MPCbased artificial pancreas [6,[32][33][34] that utilizes the Kalman filter state estimation with typically empirically tuned covariance matrix of the process noise, which normally results in its suboptimal performance, the strategy proposed in this paper additionally involves an easily identifiable prediction model of the OESO to effectively compensate for the adverse effect of suboptimality of the state estimation by correcting the system free response prediction.Moreover, this structural modification is not computationally demanding and it can be easily embedded to the existing modular MPC schemes [6,[32][33][34] without involving any hardware adjustments or directly interacting with other features of the artificial pancreas such as constraints.
On the contrary, compared with the alternative solution based on methodological estimation of the covariance matrix of the process noise according to the methods presented in [18][19][20][21][22] to implicitly ensure the optimality of the state estimation, it can be claimed that these methods typically require very large experimental datasets (tens of thousands of samples) to provide reliable and unbiased estimates, whereas the statistical models of the OESO proposed in this paper can be identified from relatively limited experimental data (hundreds of samples).
In a nutshell, the most significant contributions presented in this work include the derivation of the analytical stochastic ARMA model of the OESO and its approximations which can be used to improve the performance of the suboptimal state observer in applications when the process noise model is not exactly known or is uncertain.It can be concluded that the actual effect of the proposed strategy depends primarily on the degree of plant-model mismatch of the noise model.
and the references therein.The basal state of this model was determined with respect to the basal glycemia G b = 6 mmol/l and the corresponding basal insulin administration rate v b = 0.01 U/min.The orders of empirical model (1) were chosen as n u = n d = 4, implying the overall order n = 8.Note that the theory related to estimation of the model parameters a u d i in (3) and c u d while considering the model parameters (60) and the noise model parameters (61) as K = −0.731−0.685 −0.634 −0

Fig. 1
Fig. 1 Evolution of the OESO (k) acquired during the experiment

Fig. 4 5
Fig.4 Estimate of the autocorrelation function R η η(nT s ) of the estimated noise input for the moving average model

Table 1
Comparison of the prediction and control performance metrics

Table 2
Comparison of the prediction and control performance metrics for the autoregressive model

Table 3
Comparison of the prediction and control performance metrics for the moving average model