Real-time performance of DAE and ODE based estimators evaluated on a diesel engine

Computation and sampling time requirements for real-time implementation of observers is studied. A common procedure for state estimation and observer design is to have a system model in continuous time that is converted to sampled time with Euler forward method and then the observer is designed and implemented in sampled time in the real time system. When considering state estimation in real time control systems for production there are often limited computational resources. This becomes especially apparent when designing observers for stiff systems since the discretized implementation requires small step lengths to ensure stability. One way to reduce the computational burden, is to reduce the model stiffness by approximating the fast dynamics with instantaneous relations, transforming an ordinary differential equations (ODE) model into a differential algebraic equation (DAE) model. Performance and sampling frequency limitations for extended Kalman filter (EKF)’s based on both the original ODE model and the reduced DAE model are here analyzed and compared for an industrial system. Furthermore, the effect of using backward Euler instead of forward Euler when discretizing the continuous time model is also analyzed. The ideas are evaluated using measurement data from a diesel engine. The engine is equipped with throttle, exhaust gas recirculation (EGR), and variable geometry turbines (VGT) and the stiff model dynamics arise as a consequence of the throttle between two control volumes in the air intake system. The process of simplifying and modifying the stiff ODE model to a DAE model is also discussed. The analysis of the computational effort shows that even though the ODE, for each time-update, is less computationally demanding than the resulting DAE, an EKF based on the DAE model achieves better estimation performance than one based on the ODE with less computational effort. The main gain with the DAE based EKF is that it allows increased step lengths without degrading the estimation performance compared to the ODE based EKF.


Introduction
Model based estimation and learning for systems with both fast and slow dynamics is a wide and active research area, with different application areas such as combustion engines [1,2], and electrochemical and reactive distillation processes [3]. In the automotive industry, there is a constant pursuit of designs for cost effective systems that are able to meet stricter emission legislations and consumer demands on cheap and reliable platforms. In particular this is true for the competitive market of heavy duty diesel engines where concepts like intake throttle, exhaust gas recirculation (EGR) and variable geometry turbines 2 Background and problem motivation Differential-algebraic models appear naturally in different applications; see for example [22] for a process industry example, [23] for a power electronics example, and [3] for a chemical process example. Here, an automotive application is considered and a main reason to consider DAE's is stiffness properties of the model. Numerically efficient and stable observer solutions for embedded systems are desirable and an often limiting fact is the base frequencies in which the control system is scheduled, such as 20 or 100 Hz. This can cause problems when implementing EKF's for models with fast dynamics, e.g., stiff models, due to their need for a small step length to ensure solver stability [13].
If the step length in the numerical method is increased from a stable condition the numerical solution often starts to oscillate when the step length comes close to the stability limit and if it is increased further the solution becomes unstable. This illustrates that the combination of system dynamics and the integration method can result in numerical solutions that are unstable when the step length is too long; this situation will be referred to as unstable estimates.
Stiffness in engine models. In the automotive industry, and especially in the field of engine control a typical example, where stiff models is a problem, is when a throttle is used [10]. A throttle in the air intake system, between the compressor and engine cylinders, connects the intake manifold and intercooler volumes with a variable restriction; see Figure 1. A common model for filling and emptying dynamics of volumes, is to consider mass conservation of ideal gases in fixed control volumes [24] p CV = dp CV dt = R CV T CV V CV (W CV,in (·) − W CV,out (·)), Höckerdal Figure 1 Schematic of the diesel engine model [10] with throttle, EGR and VGT, showing model states (p im , pem, p ic , ωt, and Tem), inputs (uegr, uvgt, u δ , u th , and ne), and flows between the different components (Wc, W th , Wegr, W ei , Weo, and Wt). Rectangles with rounded corners represent control volumes.
where p CV and T CV denote the control volume pressure and temperature, R CV the ideal gas constant for the gas in the control volume, V CV the volume, and W CV,in /W CV,out are the mass-flows in/out of the control volume. The flows depend on pressures and temperatures in the surrounding systems [8], but here these dependencies are omitted for the sake of clarity. An ODE model for the intake manifold and intercooler pressures, presented in Figure 1, iṡ In operating points with fully open, or nearly open throttle, the filling and emptying dynamics become fast, resulting in oscillating estimates if the discretization method and/or step length is inadequate, i.e., too long. For heavy duty diesel engines, these stiff operating conditions are common and need to be addressed.
The fast dynamics around the throttle and p ic results in that the flows W c (·) and W th (·) follow each other closely. Therefore, one way to handle this is to remove the intercooler control volume and consider the flows in and out of that volume as an algebraic constraint. By doing this, the original ODE (1) is transformed into the following DAE: with p ic now becoming an algebraic variable. When a control volume is removed, the dynamics of the total system is affected, especially operating conditions with wide open throttle, where the dynamics becomes faster than that of the original model. Since most operation of diesel engines is with wide open throttle, it is desirable to have a model, and an EKF, for these operating conditions. Figure 2 shows a simulation of the default ODE model, defined by Intercooler, and intake and exhaust manifold pressures Intercooler, and intake and exhaust manifold pressures and together with simulations of the DAE model, with modifications of (3d) and (3e) according to (2). The measurement equation (4) is split in two, where (4a) is used for feedback in the EKF and (4b) is used as an evaluation output in Section 4 only. Two different control volume sizes, V DAE im = V im and V DAE im = V im + V ic , are evaluated. Especially it shows that the smaller volume gives undesired overshoots in the simulations of the pressures, turbine speed, and mass-flow for open throttle operation (see solid ellipses in Figure 2(a)), while the closed throttle operation is not significantly degraded (see dotted ellipses in Figure 2(b)). Therefore, the volume remaining in the DAE is increased to the sum of the original ODE's control volumes, i.e.,

DAE observer
Before presenting the EKF algorithm used for the DAE observer in Subsection 3.2, a review of previous results and algorithms is conducted and presented in Subsection 3.1, while observability of the DAE model is the subject of Subsection 3.3.

Observer design -EKF for DAE's
For continuous time ODE's, such asẋ where w and v are zero mean white noise processes with covariance intensities Q and R, EKF implementations are well known [25]. For general DAE'ṡ where x and z denote differential and algebraic variables/states, respectively this is not the case. EKF's for DAE's are treated in for example [3,22,23], and all consider semi-explicit DAE's of index 1, i.e., ∂g ∂z has full column rank. The difference between the EKF in [22] and [3] is that the latter includes the algebraic subsystem (6b) in the standard EKF algorithm. This is done by linearizing (6) with respect to x, z, and w and differentiating the algebraic part with respect to time, In the method [3], the partial derivatives , are assumed constant, and therefore there are no inner derivatives in the righthand side of (7).
Another issue with EKF's for DAE's is that system noise must be introduced with care for the problem to be well posed [26]. In the approach that follows system noise is introduced only in the differential part of the DAE (6a), which may be more restrictive than necessary but guarantees well posedness.

EKF algorithm
The algorithm used is presented here although any estimation approach could be used. The used approach has a lot in common with the one presented in [3]. Among the differences are that the state prediction is performed with the explicit forward Euler, or implicit backward Euler, and that the covariance matrix is predicted using the corresponding discrete-time linearization (8) instead of using the transition matrix φ = eĀ t|t Ts . Note that the solution using forward Euler is supplemented with an equation solver to ensure consistency of the algebraic states after prediction.
The algorithm used is summarized in Algorithm 1. For the diesel engine f is composed by (3a)-(3d), g by (3e), and h by (4a), where modifications of (3d) and (3e) are made according to (2). In Algorithm 1, G t|t is a result of the differentiation in (7) and implies that system noise is present in the differential variables of the DAE only. x 0|0 where the implication indicates solvingx t+1|t andẑ t+1|t from the system of equations, Ts the sampling time, andḠ 3. Measurement update: 4. Let t = t + 1 and repeat from 2.

Observability of DAE models
To ensure that the observer estimates converge to the true states, observability of the underlying model along the studied trajectory is central. In contrast to standard state space descriptions, for DAE's there exist several concepts of observability [27][28][29][30]. The concepts are C-, R-and I-observability, denoting complete observability, observability within the reachable set and impulse observability, respectively. Generally, descriptor systems, or DAE's, are not C-observable since they contain algebraic constraints that force the solution, and output, onto a specific manifold. For this reason, the observability within the reachable set -R-observability has been defined (see e.g., [27]), which is the concept used here. This concept needs an appropriate projection of the dynamical part of the system, sometimes referred to as the slow sub-system, onto a manifold defined by the algebraic equations, or fast sub-system. For linear time-invariant systemsĒẋ =Āx, that are regular, i.e., that det αĒ − βĀ = 0 for some (α, β) ∈ C 2 ; this can be described in terms ofĒ, A, andH.
Definition 1 (R-observable). The system (9) is called observable within the reachable set if the zero output of the descriptor system with u = 0 implies that all solutions of this system satisfy P rx = 0, where P r denotes the projection onto the right deflating sub-space corresponding to the finite eigenvalues of λĒ −Ā. -" indicate that the step length was too large for stable simulation with the current discretization method. "Bold face" entries denote the largest step length with non-oscillating and stable simulation results for each model and discretization method. b) Oscillating estimates, see Figure 3(b).
The concept of R-observability can, in the linear time-invariant case (9), be characterized algebraically according to Theorem 1 [28].

Evaluation with respect to step length
The evaluation is performed by running the two EKF's based on the default ODE model and the modified DAE model. Forward and backward Euler are used for the discretization of both observers, and the evaluation uses measurement data from an engine in a test cell; see Appendix A for more detailed information. The EKF is tuned to the system and measurement noise with the covariances Q and R. For low dimension systems, there exist methods for off-line computation of these tuning variables. However, normally at least Q is found by manual tuning while R could be found using sensor characteristics provided by the sensor manufacturer. Here, the measurement noise is affected by both sensor characteristics as well as sensor installation, and is therefore treated in the same manner as Q, i.e., through manual tuning. The observers are tuned equally with respect to model and measurement uncertainty with Q = diag(10 −4 , 0.1, 1, 1, 1), R = diag(100, 100, 1, 10 −8 ), corresponding to (3a)-(3d) and (4a), respectively. Stability limits for these settings are ∼15 ms for the ODE observers and ∼125 ms for the DAE observers which shows that a significant increase in step length can be achieved for the DAE based EKF compared to the ODE EKF; see Table 1.
The performance measure used in the evaluations is the root mean square error, whereŷ t is the estimate and y t the measurement. Since several of the model outputs are used for feedback (4a) in the EKF's, an extra model output (4b) not used in the observer, is used for the evaluation, i.e., p ic .
The evaluation highlights three aspects of the real-time performance of DAE and ODE based EKF's. It starts by analyzing the effects of stiff models in estimation, then the forward and backward Euler for the discretization are analyzed, and it finishes up by studying the possibility of increased step lengths in the ODE and DAE based EKF's with maintained estimation accuracy. Figure 3 presents the state estimates of both the ODE and DAE based EKF's, using forward Euler for the prediction, for a 100 s segment of the WHTC where two different step lengths, 3 and 10 ms, have been used. The main difference between Figure 3(a), representing simulations with 3 ms step length which is the largest step size for which the ODE based EKF estimates are stable and without oscillations, and Figure 3(b), representing a step length of 10 ms, is the oscillating estimates of the ODE observer during high mass-flow operation. These oscillating properties, observed as a band of noise in Figure 3(b) of the ODE observer, can be explained by the stiffness of the ODE in operating points with wide open throttle, i.e., high compressor air mass-flow.

Effects of stiff ODE dynamics
The estimated probability density functions (PDF) of the estimation errors for the two observers with two different step lengths are presented in Figure 4. In Figure 4(a), presenting the estimation error PDF's for a step length of 3 ms, the appearance of the two observers are rather similar indicating that for a step length of 3 ms an EKF based on an ODE would be preferable due to its less demanding computational complexity. This benefit of the ODE based observer is due to the fact that the EKF update needs to evaluate the system model only once while the DAE needs to iterate over the algebraic equation. However, for a step length of 10 ms the degradation of the ODE estimation error PDF's for the outputs closely connected to the throttle, i.e., p im , p ic and W c , is apparent in Figure 4(b).

Influence of discretization method
A comparison of selected discretization step lengths for forward and backward Euler as discretization methods is presented in Table 1. It shows that, even though there is a gain in stability using backward Euler, e.g., non-oscillating estimates for the ODE model, the upper limit of possible step lengths for stable estimation is not affected for either the ODE or DAE based EKF's. That is, the ODE and DAE based EKF's lose stability for the same step length, independent of the discretization method, the ODE  Figure 4 (Color online) Estimated PDF's of the estimation errors for p im , p ic , and Wc from observers based on ODE and DAE models with different step lengths. The PDF's of the estimation errors for simulation step lengths of (a) 3 ms are similar for both the DAE and the ODE based EKF's. For step lengths of (b) 10 ms the ODE based EKF estimation error PDF's are wider than those of the DAE based EKF which agrees with the noisy estimates of the ODE based EKF observed in Figure 3 based EKF for 15 ms and the DAE based EKF for 125 ms. A key result here is that, while some gain in stability can be achieved by using a discretization method with better stability properties, the main gain is achieved by transforming the stiff ODE model into a DAE which results in a significant increase in step length without loss of stability or estimation accuracy.
It can be noted that the implementation of backward Euler utilizes a limited number of Newton iterations for finding the solution to the differential equation. The limited number of iterations ensures an upper execution time limit for each time step which is tuned for a 3 ms step length. From Table 1, studying the ODE columns for 3, 10, and 15 ms, it can be seen that the relative execution times between backward and forward Euler are 2.5, 3.2, and 3.9, respectively. The increasing trend comes from the fact that more Newton iterations are needed in backward Euler when the step length is increased. The limiting factor for the accuracy and performance is the stability limit of the integration method. The stability region for Euler backward is much larger so it can take longer and fewer steps which in total gives a lower computation burden even though each individual step is more computationally demanding. Another interesting remark is that, considering the total system with EKF and backward Euler, a benefit not utilized in the implementation is the fact that the EKF could provide the Jacobian for the Newton iterations, which would reduce the computational complexity further for backward Euler. This is out of the scope of this paper so it is left as a future direction. Figure 5 presents the normalized RMSE of the two observers as a function of step length for both forward and backward Euler. The normalization is with respect to the RMSE of the estimates of the ODE observer with 3 ms step length. The trend of the performance as a function of the step length is clear for all signals, i.e., the RMSE for the ODE observer increases more than the RMSE for the DAE observer as the step Step length (ms) N. RMSE  length increases. Even though the DAE observer estimates are not the best for all step lengths, for example W c at 3 ms step length, the RMSE is equal or better than for the ODE for all step lengths from 10 ms and up, which represents a common timing loop in modern truck ECU's. The reason for the better performance of the ODE observer for W c and 3 ms step length originates from the change of model dynamics that was introduced when the intercooler pressure dynamics was removed, which also affected the throttle dependent dynamics, discussed in the last paragraph of Subsection 2.1. The performance of the DAE observer with a step length of 125 ms is comparable to the performance of the ODE observer with a step length of ∼10 ms. So even if the DAE observer is computationally more demanding than the ODE observer in each time step the gain in step length results in a gain in total computational effort.

Step length analysis
It can be noted that the solver used to compute the state prediction in the DAE based EKF is not optimized with respect to efficient evaluation and there exist several improvements that would reduce the computation time, see for example [13,Subsection 5.2]. Table 2 contains the data of Figure 5(a) and shows the absolute values of the RMSE for the estimates. In particular it shows the relation of the RMSE between the estimated variables that is not shown in Figure 5. For example, it shows, for the DAE observer, that no particular variable estimate distinguishes itself as better, or worse, than any other including p ic . The data of Figure 5(b) is consistent with the observations made in Table 2 for Figure 5(a) and is omitted.

Conclusion
The benefit of using DAE's instead of ODE's for EKF design in real-time applications with fixed discretization step lengths has been analyzed. The analysis is performed using an ODE model of a heavy-duty Scania diesel engine with throttle, EGR and VGT, and a reduced DAE model obtained through removal of one of the volumes surrounding the throttle causing stiff model dynamics at wide open throttle operation. In this simplification, it is also shown how the volumes can be modified to describe the original system's dynamics. Criteria for observability for DAE's are given and an analysis established that both the DAE and ODE are observable.
The step length required to get equal performance of the two observers is investigated using experimental data from an engine that has been operated to follow the WHTC. It is shown that, even though the computational complexity for each time step of the DAE based EKF is higher than for the ODE based EKF, the possibility to use more than 10 times larger step lengths for the DAE, compared to the ODE, results in a more computational efficient implementation with maintained estimation performance. The effect of discretization method for solving the differential equations is also studied. By comparing forward and backward Euler, representing an explicit and an implicit differential equation solver, it is shown that, even if there is a significant gain in increased discretization step length for backward Euler, the main gain is obtained by using a DAE based EKF.

Appendix A Engine model and data
The default model, on which the method is applied, is a 5th order non-linear state space model of a six cylinder Scania diesel engine with EGR, VGT, and intake throttle. The model states and inputs are presented in Tables A1 and A2, respectively. It is based on a model originally developed in [34] and extended with intake throttle in [10]. The modifications are that the states for the intake and exhaust manifold oxygen concentrations and actuator dynamics are removed.  The data is collected in an engine test cell at Scania CV AB in Södertälje, Sweden. The data is from a six cylinder Scania diesel engine with EGR, VGT, and throttle and was collected during a WHTC [33]. The sensor signals used are presented in Table A3. All signals are collected with a sampling rate of 100 Hz. It can be noted that the turbo speed sensor is unable to measure rotational speeds below 20000 rpm (∼2094 rad/s) and those measurements are therefore excluded in the turbo speed measurements.