Practical implementation of extended Kalman filtering in chemical systems with sparse measurements

1 Introduction

It is commonly accepted that the state observation in chemical systems is often based on the available measurements of some parameters (depending on the utilized technology) coupled with computation of remaining (not measurable) variables by means of appropriate software. The latter is referred to as the Software Sensor (SS) in [6]. The SS implies that a proper mathematical model of chemical kinetics is to be utilized in parallel with this or that nonlinear filter implemented for estimating the state of considered chemical system. The dynamic behaviour of chemical reactions is simulated effectively with use of conservation laws and presented conveniently in the form of an Ordinary Differential Equation (ODE). However, the corresponding measurement equation is assumed to be discrete due to the modern practice. The process and measurement models are supplied with stochastic noise terms because of possible random disturbances and uncertainties in the process and measurement as well as due to existing plant-model mismatch. Eventually, it is accepted in chemical research and engineering that chemical phenomena are modelled adequately by continuous-time stochastic dynamic state-space systems written in the form of Itô-type Stochastic Differential Equation (SDE)

where x(t) ∈ ℝ^n₁ is the n₁-dimensional vector of system state at time t, u(t) ∈ ℝ^n₂ is a measurable input at time t (further explanation of function u(t) is given in [24]), F : ℝ^n₁ × ℝ^n₂ → ℝ^n₁ is a sufficiently smooth drift function representing reaction kinetics, G(t) is a time-variant diffusion matrix of size n₁ × q and {w(t), t > 0} is a zero-mean Brownian process with square covariance matrix Q(t) > 0 of size q. The initial state x₀ of SDE (1.1) is supposed to be a random variable. More precisely, x₀ ~ 𝒩(x̅₀, Π₀) with Π₀ > 0, where the notation 𝒩(x̅₀, Π₀) stands for the normal distribution with mean x̅₀ and covariance Π₀.

In addition, some measurement information arrives discretely and in equidistant time intervals of size δ = t_k – t_k-1. This interval δ is called the sampling period (or waiting time) in filtering theory. The relation of measurement y_k to the system state x_k is usually fixed by the formula

Here, k is a discrete time index (i.e., x_k means x(t_k)), y_k ∈ ℝ^m is the information available at time point t_k, h : ℝ^n₁ → ℝ^m is a differentiable function, and the measurement noise {υ_k, k ⩾ 0} is a zero-mean Gaussian white-noise process with covariance matrix R_k > 0. Also, all realizations of the noises w(t), υ_k and the initial state x₀ are taken from mutually independent Gaussian distributions.

To treat nonlinear continuous-discrete stochastic models of form (1.1), (1.2), a special technology called the Extended Kalman Filter (EKF) has been designed [8, 9, 15, 30, 38]. Historically, the first EKF variant was grounded in the Euler-Maruyama method [18] and implemented as follows. The mentioned method is applied to SDE (1.1) on a time interval [t_k-1, t_k]. This results in the discretization

where δ = t_k − t_k−1 and the random variable W̃(t_k−1) ~ 𝒩(0, δQ(t_k-1)). It follows from the discrete-time stochastic model (1.3) that taking the expectation yields

with E{x(t_k)} = x̂(t_k) and E{x(t_k-1)} = x̂(t_k-1). The state vector in the stochastic process (1.3) is independent of the driving noise. Therefore the associated covariance is determined by the formula

The EKF implies that equations (1.4) and (1.5) are solved approximately in each sampling interval [t_k-1, t_k ] with use of the first-order Taylor expansion of the nonlinear drift function F(x(t), u(t)) around the filtering estimate x_k-1\k-1 at the time t_k-1:

where the Jacobian J_k-1 = ∂F(x̂_k-1|k-1, u(t_k-1))/∂x̂(t) is evaluated at the chemical system state (x̂_k-1|k-1, u(t_k-1)), and the notation HOT stands for higher-order terms of this Taylor expansion. Now substituting formula (1.6) into the moment equations (1.4) and (1.5) yields the time-update step of the EKF method:

where P_k|k-1 = E{(x_k - x̂_k|k-1)(x_k - x̂_k|k-1)^T} is the predicted covariance at time t_k, while P_k-1|k-1 = E{(x_k-1 - x̂_k-1|k-1)(x_k-1 - x̂_k-1|k-1)^T} is the filtering covariance at time t_k-1 and I_n1 stands for the identity matrix of size n₁. After arrival of a new measurement y_k, the measurement-update step of this EKF is:

where the matrix H_k = dh(x̂_k|k-1)/dx_k is the Jacobian of the right-hand side function in the measurement model (1.2) evaluated at the predicted state expectation x̂_k|k-1 from the time-update step (1.7), (1.8), and e_k ~ 𝒩 (0, R_e,k) are innovations of the KF. The described filter is the simplest but successful state estimator that has been utilized by practitioners for decades [8, 9, 15, 30, 38]. Theoretical justifications of this method in both deterministic and stochastic settings can be found in [4,34].

We stress that the notion of EKF does not imply the unique technique, as presented above, but refers to a class of various methods with different properties. For instance, Frogerais et al. [7] consider and examine five EKF implementations on two nonlinear test problems. They elaborate at least two approaches for constructing EKF schemes, which can lead to a great variety of practical filtering algorithms based on ordinary or stochastic differential equation numerical methods. Thus, some criticism published on performance of the EKF in chemical systems for offline models and industrial applications [6,12,33,35,36,39,40] does not mean that all other implementations of this method will also fail, as shown below.

It should be taken into account that many reasons, namely, poor initialization, sparse measurements, large discretization errors, strong nonlinearities, linearization errors, poor modelling (i.e., models that do not enforce physical constraints), poor tuning of noise covariance, and so on, may influence the EKF performance in practical state estimation tasks. For example, Haseltine and Rawlings [12] report that their EKF fails for a batch reactor with poor initialization. Thus, we aim at addressing the issue of accurate state estimation in stochastic chemical systems (1.1), (1.2) with poor initialization and sparse measurements in detail.

Certainly, the useful model correction technology elaborated by Schneider and Georgakis [37] improves the EKF performance for the batch reactor considered in [12]. However, it can resolve the problem of poor initialization of the traditional EKF, but does not help for chemical systems with sparse measurements. This is illustrated by numerical tests below. The reason is that the above-discussed EKF implementation poorly approximates chemical kinetics for long sampling times, and this may fail the EKF. On the other hand, adaptive ODE solvers with automatic global error control used instead of the fixed-stepsize stochastic Euler-Maruyama method resolve the mentioned difficulty and lead to a more advanced state estimation technology, which is referred to as the Accurate Continuous-Discrete Extended Kalman Filter (ACD-EKF) [22, 24,25, 26, 28]. This filter is outlined in the next section.

The mentioned ACD-EKF is designed for treating chemical systems with short and long waiting times (including variable sampling periods) in the same manner, i.e., without manual tuning. The latter is important because, for instance, Soroush [39] writes: ‘In the chemical/petrochemical and biochemical industries, there are many processes wherein the choice of sampling rate is limited by the availability of the output measurements. For example, composition analyzers such as gas chromatographs have a cycle time say 5-10 min compared to a desired control interval of say 0.1-1 min. If the control interval is increased to match the availability of measurements then control performance deteriorates significantly’. The cited inconsistency can be attacked by the accurate extended Kalman filtering technique constructed in [22, 24, 25, 26, 28]. In the present paper, we study its capacity for estimating states in chemical SDE models with poor initialization and infrequent measurements and compare it with the aforementioned traditional EKF method, which is commonly utilized for estimation purposes in chemistry research and engineering.

2 Software sensors for stochastic systems with sparse measurements

There exist two principal variants of the EKF implementation, namely, the continuous-discrete and discrete- discrete EKFs [7,22]. The continuous-discrete variant is the most suited one for treating chemical engineering tasks due to the continuous-time matter of equation (1.1). This technology is based on replacement of the predicted values of the state mean and covariance matrix determined in the time-update step (1.7), (1.8) of the EKF with the values satisfying the Moment Differential Equations (MDEs)

where J(x̂ (t), u (t)) denotes the Jacobian of the drift function F(x̂ (t), u (t)) in SDE (1.1) (i.e., J (x̂(t), u (t)) = ∂F (x̂(t), u(t))/∂x̂(t)) evaluated at the state mean trajectory, G(t) is the diffusion matrix in the process noise term, Q(t) is the covariance matrix of the zero-mean Gaussian white-noise process w(t), x̂(t) is the statemean of the random system state vector x(t) at time t (i.e., x(t) is a solution to SDE (1.1)), and u(t) is the measurable input (i.e., a known function of time).

More formally, we use filtering values of the state mean and covariance matrix calculated in the previous sampling time as the initial values of MDEs (2.1), (2.2)for solving them in the next sampling interval [t_k-1, t_k], i.e., x̂(t_k-1) = x̂_k-1\k-1, P(t_k-1) = P_k-1|k-1. Then, their predicted values are taken to be: x̂_k|k-1 = x̂(t_k) and P_k|k-1 = P(t_k). After arrival of a new measurement y_k, the measurement-update step of this filter is conducted in line with the above formulas (1.9)-(1.10).

The main difficulty associated with the continuous-discrete filtering grounded in formulas (1.9)-(2.2) is that the MDEs arisen in chemical engineering are nonlinear and, hence, should be treated numerically. Strictly speaking, any computed predicted state mean x̂_k|k-1 and covariance matrix P_k|k-1 do not further satisfy their MDEs, because these values are always determined with some error, which is referred to as the discretization error. The magnitude of the discretization error depends on the quality of implemented MDE solver. Unfortunately, large discretization errors of inaccurate numerical solutions may fully ruin the proper performance of the EKF method (1.9)-(2.2) and result in wrong solutions to the mathematical model (1.1), (1.2), especially in the case of infrequent measurements [7, 22]. Such a situation often occurs when MDEs (2.1), (2.2) are discretized by a numerical scheme for one step of the size δ. The classical EKF algorithm presented by formulas (1.7)-(1.10) is the case. Its inefficiency is clearly shown for chemical systems with sparse measurement information, below. This is quite obvious because the theory of numerical methods for ODEs [5,10,11,14] says that small discretization errors are ensured only for sufficiently small sizes of the sampling period δ, and this is not the case for chemical systems with long waiting times. On the other hand, Soroush [39] explains that such an inaccurate numerical solution seriously limits the applied potential of the traditional EKF method because it does not succeed in chemical models (1.1), (1.2) with infrequent measurements, and the short sampling periods may be technically (or by any other reason) impossible or too expensive in practice.

In this paper, we focus on the issue of large discretization errors committed in continuous-discrete stochastic systems with strong nonlinearities and sparse measurements. The above-mentioned obstacle can be overcome by means of the variable-stepsize ODE solvers with automatic error control, i.e., for example, by commonly used MATLAB codes [13, Section 12.2]. However, it is shown in [16, 23, 29,32] that usually EKF techniques grounded in general purpose ODE solvers, including MATLAB software, lose the EKF implementations based on specially created methods for solving MDEs (2.1), (2.2). This is because all MATLAB codes and general purpose ODE integrators are not intended for retaining the positive semi-definiteness of computed covariance matrix and treat the MDEs as conventional ODE systems. We recall that the matrix P(t) in equation (2.2) has the physical meaning of being the variance of the state prediction error, i.e., x(t) - x̂(t), and, theoretically, has to be positive semi-definite. Thus, they compromise the mathematical theory underlying the EKF and, hence, may reduce its practical performance.

The second shortcoming of all MATLAB solvers outlined in [13, Section 12.2] as well as of other methods developed and tested in [19, 32] (all these are utilized for implementation of the above-cited continuous- discrete EKF) is that they exploit only local error control and, hence, cannot ensure a preassigned accuracy of numerical integration of MDEs (2.1), (2.2). This follows from the fact that the local error regulated in such filters is an artificial notion, which hardly relates to the actual accuracy of numerical solution. So the committed numerical integration error (also termed as the global/discretization error) is unpredictable and can be of any magnitude, there. Thus, the global error control seems to be the necessary option of a robust filter. It is also worthwhile to mention that, in contrast to the EKF with only local error control, the methods with global error control regulate the true error of numerical integration of MDEs (2.1), (2.2), automatically. The user needs only to set the required accuracy of numerical integration (a single number) and the code does all remaining work for calculating a numerical solution to MDEs (2.1), (2.2) with the error corresponding to the set accuracy level. The latter property is the reason to name any such state estimator as the Accurate Continuous-Discrete Extended Kalman Filter (ACD-EKF) in [22, 23, 24, 25, 26, 28].

Among the ACD-EKF versions developed in the cited papers, the most promising one is based on the hybrid triple NIRK6(4)M2 and presented with all technical particulars in [24,26]. This is because of the sixth- order convergence of the output numerical solution to the nonlinear state expectation equation (2.1) exploited there instead of the fourth-order convergence enjoyed in the other hybrid triple NIRK4(2)M2 built in [23, 25, 26, 28]. The mentioned numerical scheme applied to MDE (2.1) results in the discretization

with the fixed coefficients b₁ = b₃ := 5/18, b₂ := 4/9, c12:=3−3/6,c22:=3+3/6,a112=a222:=1/2+23/9,a122=a212:=1/2−23/9,d112=−d222:=3+3/36,d122=−d212:=−3+3/36,c13:=5−15/10,c23:=1/2,c33:=5+15/10,a113=a323:=125+3915/250,a123=a313:=125−3915/250,a213=a223:=1/2,d113=−d323:=7+215/200,d123=−d313:=(−7+215)/200,d133=−d343:=1815+153/1000,d143=−d333:=1815−153/1000,d213=−d223:=1/32,d233=−d243:=33/32 , obtained at each node t_l, l = 0,1, 2,..., L - 1, of some mesh

introduced in the sampling interval [t_k-1, t_k]. Notice that the presented discrete-time equations are nonlinear and should be solved for the unknown vector x̂_l+1 approximating the state mean at time t_l+1. The stage value x^lji of method (2.3) implies an approximation to the state mean x^(tlji) evaluated at the time instant tlji=tl+cjiτl, and the measurable input ulji=u(tlji). We recall that the vector u(t) is a known function of time. The function F(‧) represents the right-hand side of MDE (2.1) or, this is the same, the drift function in SDE (1.1).

The high convergence rate is not the only advantage of the scheme under consideration. It is also A-stable, stiffly accurate, symmetric and conjugate to a symplectic method up to order 6 at least. So, discretization (2.3) has sufficiently high stage and classical orders, which are 3 and 6, respectively, and some stage values calculated in a step of this scheme are rather accurate and allow for a sixth-order dense output (see a further discussion of the presented method in [21, 27]). Most importantly, the nonlinear problem (2.3) admits a sufficiently cheap iteration for finding the unknown vector x̂_l+1, which is presented with all necessary particulars in [24]. Also, to generate automatically a good mesh (2.4) in each sampling interval [t_k-1, t_k], our ACD-EKF exploits the local and global error control mechanisms designed for method (2.3) in [20, Algorithm 3.2]. These controls allow the committed numerical integration error to be made negligible in automatic mode.

Next, the covariance equation (2.2) is discretized on the above mesh (2.4) generated by the aforementioned variable-stepsize solver grounded in formulas (2.3) by means of the modified implicit mid-point rule [32]. The cited discretization reads

where t_l+1/2 = t_l + τ_l/2 and K_l+1/2 and M_l+1/2 are defined as follows:

with the Jacobian matrix J_l+1/2 = ∂F(x̂_l+1/2, u(t_l+1/2))/∂x̂(t) evaluated at the point (x_l+1/2, u(t_l+1/2)). In formulas (2.6), the vector x̂_l+1/2 is taken to be equal to the stage value x^l23 from the discretization method (2.3) (see an additional explanation in [24, Section 2.3]). The error of approximation (2.5), (2.6) is not regulated in our filter and considered to be negligible because of the linearity of MDE (2.2).

Mazzoni’s method (2.5), (2.6) is A-stable and convergent of order 2 [32]. However, the main advantage of using this method in practical state estimations is its symmetry and positive semi-definiteness. However, the latter properties are preserved in exact arithmetic, only. So, round-off operations existing in reality might compromise the symmetry and positive semi-definiteness of the computed covariance P_l+1 and fail the ACD- EKF. On the other hand, it is known that square-root implementations of Kalman-like filtering methods based on stable orthogonal rotations allow the committed round-off errors to be reduced dramatically (see, for instance, [1, 2, 3, 9, 17]). That is why we utilize the square-root ACD-EKF variant published in [24, Appendix A] for showing the superiority of the filtering technique (1.9)-(2.2) obtained in the continuous-discrete approach in comparison to the classical EKF method (1.7)-(1.10) built within the discrete-discrete filtering, in particular for chemical stochastic SDE systems with sufficiently long sampling periods. In other words, we address the inconsistency stated for chemical continuous-time stochastic models with sparse measurements by Soroush [39].

We emphasize that the ACD-EKF presented in [24, Appendix A] does not require any manual tuning, maybe, except for altering the value of accuracy parameter ε_g (its default value ε_g = 10^–4 is accepted in the numerical experiments below). This method is expected to work successfully for short and sufficiently long sampling intervals δ (see more explanation in [22, 24, 25, 26]). In the next section, our argument is demonstrated numerically on a chemical system from [12], for which the classical EKF can fail, in order to exhibit differences in performance of the traditional state estimation technology and the contemporary ACD-EKF.

3 Numerical simulation and discussion

Following Haseltine and Rawlings [12], we consider the gas-phase reversible reaction with three species denoted as A, B, and C:

in which the fixed coefficients k₁ = 0.5, k₂ = 0.05, k₃ = 0.2 and k₄ = 0.01. The stoichiometric matrix v of reaction (3.1) and the reaction rates r are chosen to be

where c_A, c_B, and c_C denote the concentrations of the species A, B, and C in moles per liter, respectively. Therefore, the state of this chemical system is defined by the vector

The measurement equation is given in the following form:

Here, x_k means an approximation to the state vector x(t_k), R is the ideal gas constant, T is the reactor temperature in Kelvin, and the measurement noise is υ_k ~ 𝒩 (0, R_k) with R_k = 0.25² in each measurement information y_k. Following [12], we set RT = 32.84. Below, we consider the task of state estimation for a Batch Reactor (BR) discussed in the cited paper.

The well-mixed, constant-volume, isothermal BR is modeled by the SDE

with the initial values

where the matrix υ and the vectors r, x(t) are defined in (3.2) and (3.3). The matrix G is assumed to be constant and diagonal, i.e., G= diag{1,1,1} in (3.5). The continuous-time process noise w(t) is the Brownian motion with the fixed covariance matrix Q= diag{10^–6/δ,10^–6/δ,10^–6/δ} because the covariance matrix of the discrete-time zero-mean Gaussian process is Q̃= diag{0.001², 0.001², 0.001²} in [12]. This process noise is added to simulate possible random disturbances and uncertainties in the BR and also due to existing plant- model mismatch in reality. The sampling period δ = t_k+1 - t_k = 0.25 is set in the first experiment. The initial information for the classical EKF and for the novel ACD-EKF method is supposed to be the same and poor, i.e., x̄₀ = [0,0,4]^T and Π₀ = diag{0.5², 0.5², 0.5²} (compare x̄₀ with the exact initial state given by formula(3.6)). We stress that the numerical simulation results presented for the traditional EKF in [12] exhibit convergence to a wrong steady-state and even negative concentrations. Now we repeat the same experiment, but with our own codes. The state estimation methods and all computations are implemented in MATLAB.

To simulate a measurement history, we integrate SDE (3.5) by the stochastic Euler-Maruyama solver [18] with the small fixed step size equal to 0.001 and apply formula (3.4) with the additive noise υ_k to calculate a measurement information y_k at every sampling instant t_k. Then, we solve the reverse problem, i.e., we treat SDE (3.5) by means of the traditional EKF and our ACD-EKF methods with the same above-set initial data and simulated measurements. The results of this numerical experiment (i.e., the estimated concentrations c_A, c_B, and c_C) are shown in Fig. 1, where plots (a) and (b) display one typical behaviour (Regime 1, i.e., Convergence Regime) of the state estimators, when the convergence to the no noise solution is clearly observed, and plots (c) and (d) exhibit the divergence and failure (Regime 2, i.e., Divergence Regime) of both EKF implementations, i.e., when the wrong steady-state and even negative concentrations are calculated after convergence. We stress that the no noise solution has been computed by the ODE solver NIRK6(4) [20] (with the global error control) applied to the underlying ODE (3.5), (3.6), i.e., without the noise term Gdw(t). The latter solution is shown additionally in all figures of this paper and marked as ‘No noise’. The steady-state values of the no noise, EKF, and ACD-EKF estimated concentrations in the convergence and divergence regimes are gathered in Tables 1 and 2, respectively.

Fig. 1

Actual (no noise), EKF and ACD-EKF estimated concentrations (with initial values [0, 0, 4]^T and the initial covariance Π₀) for the batch reactor: Regimes 1 (convergence, top plots) and 2 (divergence, bottom plots).

Having compared the derived results to those published in [12] we see a principal difference. Our simulation shows that about 25% of EKF runs and 17% of ACD-EKF runs are successful. The exact numbers of fails of these state estimators out of 1000 runs are presented for this numerical experiment in the first column of Table 6, below. Note that Schneider and Georgakis [37] also confirm the existence of these two performance regimes for the batch reactor (3.5), (3.6) and the inexact initial state vector x̄₀ = [0,0, 4]^T and the initial covariance Π₀.

In Regime 1, the EKF and ACD-EKF estimated concentrations converge to the actual steady-state in contrast to the results published in [12], where the convergence only to the wrong value is reported. Our Tables 1 and 2 exhibit good accuracies of the estimated steady-states in Regime 1. Moreover, the traditional EKF algorithm, given by formulas (1.7)-(1.10), and the advanced ACD-EKF method from [24, Appendix A] do not produce non-physical data (i.e., negative concentrations) in this performance regime. Certainly, negative values are possible during the convergence (transient) period (see Figs. 1a and 1b) because any filtering method is just looking for the true solution in this stage. However, both filters converge to the true solution of the BR (3.5), (3.6) in Regime 1 and, then, all the concentrations become positive, as predicted by theory.

It is worthwhile to recall that altering the true initial values (3.6) to the untrue ones given by the vector x̄₀ = [0,0,4]^T and utilized by the filter obliges the latter to solve another, say, untrue mathematical model. Besides, the EKF does not know that the model is wrong and follows its wrong state trajectory, which has no relation to the original system (3.5), (3.6). On the other hand, at every sampling time, the filter receives some information on the true state of the batch reactor (i.e., true measurements). Then, this information is used to amend the calculated wrong state in the measurement-update step of the filtering procedure. In other words, after each measurement-update, the filter moves from one state expectation trajectory of the process model to the other, which is expected to be closer to the true solution of the considered problem. A certain amount of true measurements allow the state expectation trajectory of the BR (3.5), (3.6) to be identified correctly. The traditional EKF and ACD-EKF methods expose good accuracies of their convergence regimes after reasonable time, i.e., when t_k ⩾ 10 in this numerical experiment. We stress that the filter convergence time (the transient period) is not prefixed by the model and depends on many factors including the frequency of measurements, i.e., the sampling rate δ. Thus, the wrong states computed in the transient period of the filter (i.e., for t_k < 10) are not the concentrations of the BR under examination. Their values may be positive or negative because they do not relate to the original mathematical model (3.5), (3.6). That is why they must be merely ignored and excluded from consideration. The states estimated by the converged filter, when t_k ⩾ 10, are only of the interest. We stress that the states estimated in the transient period of filtering are of no sense, but the convergence time itself is an important characteristic of any practical SS. A further discussion on convergence issues of the Kalman filter, also known as the transient behaviour, for linear and nonlinear stochastic models can be found, for instance, in [15,31].

Eventually, the EKF and ACD-EKF estimated steady-states in Regime 1 are positive and correspond well to its no noise value (see Tables 1 and 2) in contrast to the data published in [12]. On the other hand, Regime 2 is in line with the outcome reported in the cited paper. Thus, if we take now into account the numbers of fails presented in the first column of Table 6 we will conclude that the EKF and the ACD-EKF work poorly for the BR model under consideration. The reason is clearly identified in [37]. This is because the continuous-discrete stochastic system (3.5), (3.4) with the above discussed initial state x̂₀ and the initial covariance Π₀ is ill-posed.

In other words, the matrix Π₀, which has the physical meaning of being the variance of the initial state error, i.e., x₀ -x̄₀, is too far from the exact covariance corresponding to the initial state given by formula (3.6). Thus, altering this matrix to the modified initial covariance Π0mod = diag{0.5², 0.05²,4²} recommended in [37] fully resolves all the problems and, then, the EKF and ACD-EKF state estimators work well and without a single fail out of 1000 independent runs for the inexact initial state vector x₀ = [0, 0, 4]^T (see Fig. 2 and Tables 3 and 6). In addition, both filters demonstrate their swift convergence in this experiment.

Fig. 2

Actual (no noise), EKF and ACD-EKF estimated concentrations (with initial values [0, 0, 4]^T and the modified initial covariance Π0mod) for the batch reactor.

Next, we show that the initialization technique elaborated in [37] is not capable of addressing the problem of infrequent measurements discovered in [39] in the frame of the traditional EKF, and the advanced ACD-EKF helps a lot. For that, we assume that the BR (3.5) is to be estimated by the same EKF and ACD-EKF methods with the initial data x̄₀ = [0, 0, 4]^T and Π0mod = diag{0.5², 0.05², 4²}, but for the longer sampling period δ= 2.0. The behaviour of estimated concentrations and their steady-states are displayed in Fig. 3 and Table 4, respectively.

Fig. 3

Actual (no noise), EKF and ACD-EKF estimated concentrations (with initial values [0, 0, 4]^T and the modified initial covariance Π0mod) for the batch reactor.

The latter numerical experiment exhibits that the ACD-EKF clearly outperforms the traditional EKF. Figure 3 and Table 4 show that despite the slightly worsened convergence to the true solution (because the ACD- EKF needs a certain amount of measurements to converge, and these have been reduced in our experiment), the modern method computes the true state estimates to the BR (3.5) and its steady-state with good accuracy. In contrast, the classical state estimator is not able to capture correctly the BR dynamic behaviour because of poor approximation of the chemical kinetics and, hence, it produces wrong trajectories of the concentrations and their untrue steady-states in the majority of runs. We have done 1000 independent runs of our codes and revealed no fail for the ACD-EKF, whereas the classical EKF has been successful in less than 1% of runs (see the precise numbers of fails of these state estimators in the third column of Table 6). Here, we do not show the convergence regime of the EKF for the BR under examination because the method working successfully with the probability of less than 1% is of no practical significance. In addition to the accuracy of state estimation, the ACD-EKF is flexible and insensitive to the sampling rate δ. We stress that exactly the same ACD-EKF code has been run for δ = 0.25 and for δ= 2.0.So this filter accommodates changes in the length of the sampling interval, automatically [22,24, 25, 26]. That is why it can also be efficient in chemical systems with irregular sampling.

Next, we show that even the exact initial data do not ensure the accurate state estimation by the traditional EKF method when the waiting time of the chemical system is sufficiently long. To simulate such a situation, we merely solve SDE (3.5) assuming that its exact initial value (3.6) is known, i.e., the initial values for the EKF and ACD-EKF techniques are now taken to be x̄₀ = [0.5, 0.05,0]^T and Π0exact = diag{0.01², 0.01², 0.01²}. Note that the latter correspond well to the theoretical stability requirement of EKF [34]. At the same time, we suppose that the measurement information arrives discretely and in the interval δ = 4.0. The outcome of our last experiment is displayed in Fig. 4 and Table 5.

Fig. 4

Actual (no noise), EKF and ACD-EKF estimated concentrations (with the initial values [0.5, 0.05, 0]^T and the initial covariance Π0exact) for the batch reactor.

Again, we observe a great difference in the behaviour of the traditional EKF and modern ACD-EKF methods. We see that the classical state estimator has not converged yet in the interval [0,30] (see Fig. 4a). So, its ‘steady-state’ shown in Table 5 is not actually a steady-state, but only the value calculated at the last sampling point in the simulation interval. At the same time, the ACD-EKF works well and reconstructs accurately the concentration trajectories (see Fig. 4b) and the chemical system steady-state (see Table 5). Moreover, the described picture retains for all 1000 independent runs of the codes (see the last column in Table 6). All this confirms the superiority of the modern state estimation technique presented in [24, Appendix A] towards the traditional EKF method (1.7)-(1.10). The latter experiment also explains that the method of model tuning [37] may be useless when the state estimator is not able to reconstruct correctly the true chemical system dynamics due to sparse measurements. In contrast, the ACD-EKF is flexible and insensitive to the size δ of the utilized sampling intervals. The presented computation results have been obtained by the same code, i.e., without any manual tuning for a specific δ. This is a clear evidence that our novel filter can be useful for addressing the contradiction between the measurement time available in reality and the optimum sampling rate of the filter stated in [39]. In other words, the ACD-EKF works for much longer waiting periods without significant deterioration of the accuracy and robustness of state estimation in chemical systems. The latter is important when frequent measurements are technically (or by any other reason) impossible or too expensive in practice.

In conclusion, Table 6 summarizing performances of two state estimators applied to SDE (3.5) with various initial data demonstrates that when the given mathematical model is ill-posed both EKF implementations work poorly. However, the initialization technique recommended in [37] fully resolves this difficulty and all the methods expose the excellent performance with no fail. Besides, when the sampling rate δ is reduced only the ACD-EKF retains accurate. The latter confirms the promising applied potential of the ACD-EKF and creates the necessary background for its further practical testing, including in industrial environment.

4 Conclusion

This paper investigates the issue of the EKF failure reported for the batch reactor in [12]. For resolving this problem, two different approaches have been examined above. The first one is grounded in the model correction technology recommended in [37]. It improves the EKF performance by addressing the issue of poor initialization of the state estimator. In other words, Schneider and Georgakis suggest reformulation of the chemical system such that the conventional EKF method works well for the batch reactor studied in [12]. The sole shortcoming of this approach is that the technique discussed in [37] may fail at least for some chemical systems with infrequent measurements. The reason is that the ‘textbook’ EKF given by formulas (1.7)-(1.10) poorly approximates swiftly changing chemical kinetics for long sampling intervals, and this may also fail the EKF. The second approach is based on adaptive ODE solvers with automatic global error control used instead of the stochastic Euler-Maruyama method. It reconstructs accurately chemical model dynamics and resolves the problem of sparse (and even irregular) measurements. Thus, this approach underlies advanced Accurate Continuous-Discrete Extended Kalman Filters (ACD-EKF) designed in [22, 23, 24, 25, 28].

Actually, our research and the technique considered in [37] deal with different aspects of the problem of EKF failure in chemical research and engineering. We do not change the ‘model’ (i.e., the chemical system with initial values and noises), but improve the ‘method’ (the state estimation technique) so that it works for long waiting times. In contrast, Schneider and Georgakis [37] improve the ‘model’ such that even the classical EKF treats the examples of the EKF failure published in [12], successfully. Obviously, these two approaches are equally important and only their combination will lead to most effective state estimation techniques for difficult offline chemical models and industrial applications in future.