Kalman and particle filtering methods for full vehicle and tyre identification

ABSTRACT This paper considers identification of all significant vehicle handling dynamics of a test vehicle, including identification of a combined-slip tyre model, using only those sensors currently available on most vehicle controller area network buses. Using an appropriately simple but efficient model structure, all of the independent parameters are found from test vehicle data, with the resulting model accuracy demonstrated on independent validation data. The paper extends previous work on augmented Kalman Filter state estimators to concentrate wholly on parameter identification. It also serves as a review of three alternative filtering methods; identifying forms of the unscented Kalman filter, extended Kalman filter and particle filter are proposed and compared for effectiveness, complexity and computational efficiency. All three filters are suited to applications of system identification and the Kalman Filters can also operate in real-time in on-line model predictive controllers or estimators.


Introduction
Advanced automotive control requires accurate dynamic models to predict vehicle response. In recent years, dynamic modelling has also become increasingly important in vehicle development, due to obvious time and cost saving advantages over traditional prototype testing [1]. Regardless of the application, a compromise between accuracy and complexity needs to be sought, as the ability to run in real time is a major goal in both control and simulation. Hence, the employed models need to be simple and efficient, rather than complex and heavy. The validity of these necessarily simplified models depends on many fixed or estimated parameters. Some of these, such as the centre of gravity (CG) height or the moment of inertia are difficult and expensive to measure. Others, as pointed out in [2] depend on external factors: tyre adherence boundaries for example are heavily influenced by asphalt conditions and precise lab measurements often don't correlate well with the real world. Furthermore, even if these values are physically accurately set, the simplified model can be made to perform better if they are adaptively tuned or completely identified.
Most models in literature are either derived from first principles or generated with multibody techniques. Structured identification, also referred to as grey-box parametrisation can then be employed to identify one or more parameters of this perfectly known mathematical or multibody structure. Models based on Neural Networks are also found [3], but these are black-box and have the disadvantage of not giving any engineering insight into the system.
Common identification tools include the Least Squares, Maximum Likelihood, and Recursive Prediction Error method [4]. And the well-known Kalman Filter [5] which prevails in a large portion of literature, principally in its nonlinear form, the extended Kalman filter (EKF). Several authors have successfully employed EKF to identify a limited number of model parameters, which are concatenated to the state vector and estimated simultaneously with the true states [6,7]. This approach has later been extended to wholly concentrate on parameter estimation [8] and recent findings suggest that Kalman filter methods can be applied to identify all the parameters of any well-conditioned model structure [9]. The Unscented Kalman filter (UKF) has emerged in the last two decades as the main alternative to EKF. First developed for nonlinear state estimation by researchers at Oxford University [10], it is based on the idea that approximating a nonlinear statistical distribution is easier and more accurate than linearising a nonlinear function, as done by the EKF. It avoids the need to calculate Jacobians and is computationally less expensive and easier to implement. It is also easily extended to system identification and dual estimation problems [11,12], in a similar fashion to the EKF. An example of adaptive parameter estimation in vehicle dynamics has been developed in [13] and a recent study [14] includes tyre parameter estimation, though this is carried out using simulation only. Other recent UKF studies consider state estimation only. The particle filter (PF) belongs to the group of recursive Monte Carlo methods and is particularly suited to harsh nonlinearities and non-Gaussian applications [15]. It approximates the posterior probability density function (pdf) of the state vector in a similar way to the UKF, but uses a much larger set of samples, which are randomly selected from an initial uniform distribution. In certain applications the PF has proven to be a promising alternative to gradient-based Kalman Filters [16], since the sampling representation of the pdf is in general a better approximation of non-Gaussian distributions caused by nonlinear model functions. However, examples in literature usually rely on high computational power and low sampling rates, especially for real-time solutions [17].
The focus of this paper is on grey-box automotive parameter identification and three of the most common and promising identification methods are used for the purpose: the EKF, the UKF and the PF. Several past publications have considered the identification of individual parameters of a vehicle model, particularly for tyre and friction coefficient estimation [6,18]. Longitudinal, lateral and vertical dynamics have been considered separately, to estimate a different parameter from each model [19]. And in [20] a method is presented for simultaneous identification of a larger set of parameters, though this method relies on a-priori knowledge of Pacejka tyre coefficients. Here we extend the findings in [21], identifying all the independent parameters of a whole vehicle handling model simultaneously, including longitudinal, yaw and roll freedoms, and with independent combined-slip load dependent tyres. We consider data collected from a test vehicle carrying out medium to high magnitude manoeuvres including wheel-spin and terminal understeer, in order to build a model which is valid over the whole range of the tyres.
In the next sections each filter is first introduced and described in its identifying form. A simplified case is then presented -identification of the front and rear tyre stiffness of a linear single-track model, based on test vehicle data. This example serves as an immediate and straightforward comparison of the different techniques, illustrating how they are tuned for performance and efficiency. The major study of identification and validation of the more complex four degree of freedom full vehicle handling model then completes the paper.

Identifying EKF
The standard state estimating EKF operates on nonlinear system and sensor models f and h, which relate the state vector x, measured sensor set y, known inputs u and model ω represents state propagation and modelling error, υ is the sensor error, and an optimal filter can be derived using estimates or expectations of the error covariance matrices: The EKF also requires model Jacobians to be evaluated at each time step, defined For further details of the full EKF see [6]. The identifying EKF (IEKF) takes advantage of the fact that f and h are general nonlinear functions of x and θ, defining an extended state vector z with extended state derivatives set zero for the parameter states:ż It is computed using a sequence of equations which develop a time-varying estimate of state error covariance, P k and Kalman gain K k ; at each time step of the recorded time histories, compute The final equation in set (6) combines Euler integration of the system using time step T with state and parameter adaptation driven by the output error (known as the innovation sequence).

Unscented Kalman filter
The UKF identifies its own error statistics at each iteration, and hence avoids the need to use Jacobians. According to Julier and Uhlmann [10] and Wu et al. [22], a sample of (2n + 1) so called sigma points χ are selected around the nth order state vector, at each instant k in time: where { √ (n + κ)P k } i is the ith column of the matrix square root of (n + κ)P k (obtained here using Cholesky decomposition) and P k is the current estimate of state error covariance. These sigma points are propagated by the model, here as for IEKF, using Euler approximation: and intermediate estimates for the propagated state and covariance matrix are computed by weighted averages: With W 0 = κ/(n + κ) and for all other i, Similarly, average outputs are obtained according to the output model: The UKF then propagates output error covariance according to the transformed sigma points: and uses this together with a cross correlation estimate to find the Kalman gain by K k+1 = P xy P −1 yy . State and covariance estimates are then updated using the innovation sequence, in a similar way to IEKF: The above UKF derivation is common to all implementations of state and parameter estimation, so the simple substitution of x in the above, with z from Equation (5) provides the identifying filter.

Particle filter
The PF is a recursive Bayesian estimator based on Monte-Carlo simulations. According to Gustafsson et al. [15] and Gordon et al. [17], the posterior density of the state vector is approximated through a set of n p so called particles, initially generated from a uniform distribution on a finite interval [a,b]: So each particle χ i0 is an nth order state vector at time k = 0 with each of its elements selected on a uniform distribution with range [a j , b j ] selected appropriate to the minimum and maximum expected value of that state. Each particle is propagated by the process and output models at every instant k, according to: where Euler integration is used in Equation (15), as for IEKF/UKF. A weight W i is then assigned to each of the particles, based on the error between the propagated output (16) and the true output of the system. These errors are assumed to have a known distribution and zero mean, hence: where n y is the number of outputs and ϕ j is the assumed pdf of the jth output innovation sequence, with zero mean and variance σ j . Here we assume a normal distribution of innovations and hence use The weights are normalised so that their sum equals unity: State and output estimates are then computed by a weighted sum of the propagated particles: The filter then resamples n p particles from the original set using W i as the probability of reselection. In this way, where a new full set of n p particles is generated, these include repeated instances of successful particles and exclude low-weighted particles. The new set of particles for the next iteration is then finalised by randomising the state combinations from the pool of these resampled particles. Equations (15)- (21) are then recursively repeated through the data points.
The PF is typically used for state estimation, for which all n p particles are necessarily propagated, weighted and resampled at least once at each time step k. For identification this has the disadvantage that convergence occurs over a very short section of data (though n p is typically high -around 10 5 -so the process is still computationally intensive). To enable parameter optimisation over a more representative, longer section of identification data (or for continuous identification in real-time) it is sensible to propagate just one particle per time step (i = k). Since the PF operates without a state covariance estimate there is then a further advantage that the true states, x need not be included in the identification set. In converting the PF for identification we therefore select χ ≡ θ so Equations (15)-(16) becomex Identification using batch data then progresses in an obvious way, using multiple passes through the available data to exercise all n p particles. Weighting and resampling operations occur after each n p time steps as in the original method.

Implementation in a simplified case -identification of a linear model
Test data were obtained from a '08 MY Jaguar XF equipped with an OXTS 3200 inertial navigation (IN) device, driven on a dry proving ground. A range of high and low magnitude, separated and combined slip manoeuvres were conducted, exciting both dynamic and steady-state vehicle responses. The controller area network (CAN) data was inter-sampled at a constant 100 Hz to match the IN (setting T = 0.01) and all data was then digitally filtered at 10 Hz to remove higher frequency noise; this is necessary to prevent instability in the filters, as they are propagated by Euler integration. Although here we perform a batch data process, the identifying filters are designed to work equally well in real-time, and realtime parameter tuning could easily be achieved using a combination of dynamic low pass filters and/or higher sampling interval T.
To illustrate execution, convergence and computational effort of the filters, consider first the identification of a linear bicycle handling model with states 60 s of low magnitude random steer data from a constant speed test is applied to each filter using the recorded steer as input, and yaw rate as output Identification is achieved iteratively by 'rinsing' the data repeatedly through each filter, starting with an initial, nominal parameter set. Both Kalman filters require covariance estimates and here Q is fixed throughout as Q = ρI. ρ is the only tuning parameter in both Kalman filters; it weighs the expectation of accuracy of the model f, and particularly the assumptionθ = 0, so it controls the variation rate of the parameters (see [9] for more detail). For any given estimateθ of the identified parameters, R can be obtained numerically from the covariance of υ using Equation (2); here, R is recomputed after each iteration. Finally, P 0 = Q, and for the UKF κ = 1. Convergence and accuracy of the PF is affected by the choice of n p . Figure 1 shows the effect on accuracy and computational effort of varying these core tuning parameters ρ and n p (upper plots). The lower plots then show the convergence of the identified parameters and their corresponding accuracy for a single optimisation using suitable choices, ρ = 10 −8 and n p = 10 6 . Accuracy here is quantified by percentage explanation of any model time historyŷ estimating measurement y: Interestingly, both Kalman Filter methods yield identical results (and we will see indistinguishable difference in the more complex full vehicle case later). Although computed using different techniques, both state covariance and parameter estimates match within a few seconds of filtered data. As the method of computation varies however, so does the processing time; 100 cycles of data takes the UKF 23 s whereas the EKF takes just 9 s (on the mid-range PC used here). The Kalman filters also converge smoothly and definitively at all settings of ρ. E does not rise for ρ < 10 −8 and even at extremely low ρ, these filters yield the best possible accuracy from this data, E = 98.0% in under 3 min. The PF is effective at optimising C α , achieving E = 97.9% in all cases where n p is sufficiently high. However, convergence of the parameters is less definitive, Figure 1(d), and the computational cost becomes high as performance reaches acceptable levels. The final result for n p = 10 5 takes 10 min to converge, and the illustrated n p = 10 6 case takes 90 minutes. For such a small unknown parameter set these are disappointing results; although the PF is effective, it is not competitive with either Kalman filter at delivering well converged, timely results.

Full vehicle model
For full vehicle and tyre identification a suitably complex model structure is the well-known three degree of freedom handling model, simulating yaw, roll, and sideslip using a load dependent, combined-slip Pacejka tyre model. A fourth, longitudinal degree of freedom is required to exercise both longitudinal and lateral tyre slip regimes. This model structure was previously employed in [23]; the principal equations of motion are lateral : roll kinematics :φ = p.
Standard SAE axes are used, fixed relative to the vehicle wheelbase, and the wheels are labelled (1)(2)(3)(4) in ascending order as (front-left, front-right, rear-left, rear-right). Equal halftrack c is assumed, with axles separated distances a and b from the front and rear axles, respectively, and vertical suspension geometry is based on fixed roll centres h R(f,r) and CG height h G such that The forces controlling the vehicle body motion (F x i , F y i ) allow for large steer angles based on lagged tyre forces, where each of the 8 elements are lagged to simulate relaxation within the tyreḞ * The tyre forces (F t x i , F t y i ) are modelled according to a slightly simplified Pacejka magic formula using normalised slip and isotropic similarity scaling [24,25]. The normalised slip vector is where S is the longitudinal slip ratio, and α is the slip angle at each tyre contact patch, based on wheel-oriented velocities and The resulting tyre force vector is then Vertical tyre loads F z are calculated from static weight distribution, modified to accommodate lateral load transfer using separate front/rear distributions according to: and longitudinal load transfer, according to a rigid, zero pitch approximation,

Full vehicle identification
For full vehicle and tyre identification, input data is comprised of steer angle and the four wheel speeds; the outputs are CG longitudinal and lateral accelerations and body roll rate, modelled as All the above measurements with the exception of body roll rate were taken from the vehicle CAN; although roll was measured using the IN, this sensor may already be available in other CAN sets, or would be cheap to add. IN measured vehicle speed, roll angle, yaw rate and lateral velocity are also used for validation of the identified model. One 100 second test including dynamic and steady-state excitations of steering, brake and acceleration inputs was used as the identification set. This test employed random steer inputs carried out on a circuit of the proving ground track, including some extreme high slip inputs in combination and separated from medium to high acceleration and braking applications. A variety of other tests using combined and separated slip were recorded for validation. The only difficult task is the decision of which parameters to adapt in the identification, and which to fix. The filters can identify any number of parameters, but if they are not independent in their influence on the recorded outputs, the parameters will diverge. After various trials the following identification set was established as the minimum non-divergent set: This allows full identification of separate front and rear individual tyre-suspension models but with roll stiffness and damping applied at a constant ratio, known from manufacturer supplied data. The remaining fixed parameters constrain the weight balance, roll inertia and essential geometry (L = wheelbase), and are set Friction is identified via the tyre D parameters, so we set μ = 1. To illustrate the effect of identifying too many parameters, we will also consider the case of adding roll centre heights h Rf and h Rr to the identified set; these are referred to below as the h+ case.
Attempts to identify this full vehicle case using the PF were unsuccessful. Despite lengthy optimisations with high n p the larger set of parameters proved impossible to identify. Results with both Kalman filters are encouraging however; Figure 2 of the identification data. Note how all 14 parameters converge completely in Figure 2(a), as do the values for performance and trace(P k ), in plot (c). The standard IEKF and UKF results are almost identical; it is just possible to make out tiny differences in the parameters as they converge in Figure 2(a). Final parameter values are given in Table 1 and the tyre/suspension model is illustrated in terms of pure slip, in Figure 2(d); note that the absence of suspension derivatives in the simplified model means that the tyre model here is actually mapping the combined tyre/suspension characteristic. All values are in the expected range; the front tyre has a lower stiffness and quicker saturation than the rear, which is consistent with expected front steer compliance seen in an earlier study on the same vehicle [26]. Although it achieves slightly higher output explanations - Figure 2(c), the h+ case results in a less well-conditioned model, with unphysical negative h Rf and high K φ and h G . Note that the result in Table 1 is a snapshot at the 200th iteration; Figure 2(b) highlights the divergent parameters.
As these results optimise parameters for the given simplified model structure, we are not seeking parameters which exactly match manufacturer's data for the test vehicle. However it is a test of model conditioning that the identified parameters appear to be feasible, so manufacturer's vehicle data is also given in Table 1 where direct comparison is meaningful. The only significant difference is in yaw inertia; interestingly the identified model employs a lower front/right difference in tyre forces coupled with lower than expected yaw inertia. A further validation check of steady-state metrics also shows a good match; the model roll rate is 5.2 deg/g compared with manufacturers data of 4.5 deg/g and linear range understeer gradient is 4.9 deg/g which compares favourably with separate steady-state tests carried out on this vehicle found in the range 5-7 deg/g. Figure 3 shows a section of data from the identification set, illustrating the final performance explanations achieved on each of the three output variables. The range of amplitudes in the inputs is also clear here; seeking to characterise the full range of the tyre, over the 100 s longitudinal accelerations range from −7 to +4 m/s 2 , occasionally inducing wheelspin, and lateral accelerations of ±8 m/s 2 cover the full range through saturation. Note also the inclusion of lower amplitudes and the combination of dynamic and steady-state inputs.  Evidence of success in the identification is then given in Figure 4 which shows the response of four additional model outputs not included in the identification filter, over the same section of data. We would expect the excellent longitudinal velocity match, since wheel speeds are inputs to the model, but well matched roll angle and in particular lateral velocity histories show that the parameters actually capture the tyre/suspension and roll modes very well.
The standard, physically appropriate model is also validated using independent test data in Figures 56-7. Note that in earlier testing [21], poorly calibrated instrumentation meant   that lateral velocity drift was present. These validation results were therefore generated after recalibration of the INS, 6 months after the original identification tests, on the same track and in similar weather conditions. These results thus also demonstrate repeatability and robustness of the identified model, over time. Figure 5 shows the result of a step steer manoeuvre which is held while heavy braking is applied approx one second later. Again in addition to accurate results in the identified acceleration outputs, the model is also capable of matching roll angle and lateral velocity very well. Finally, Figures 6 and 7 show further positive results, validating the combined-slip tyre/suspension model in single-slip tests. Figure 6 shows accurate results for a straight line braking and acceleration test under zero steer. The modelled longitudinal acceleration is indistinguishable from the measured data throughout, despite the high k x excursion caused by wheel-spin around 10 s; lateral accelerations track correctly at zero. Finally, Figure 7 sees high accuracy in the single slip lateral case of a constant speed, random steer manoeuvre. Both of these tests explore the tyre in low to medium magnitude and also in over-slip saturation conditions. Both IEKF and UKF identified the same successful model here, so the decision over which filter is 'better' depends only on computational efficiency and complexity of code. UKF has the advantage that Jacobians are not required, so it is simpler to code, though Matlab's Symbolic toolbox means this overhead in complexity is slight. For the simple model in Section 5 the IEKF was clearly faster, but here processing times are only slightly different and the order is reversed; each IEKF iteration takes 266 s whereas each UKF iteration takes 212 s.
The longer processing times are not simply due to the more complex model. As the number of parameters increases, the UKF has more sigma points to process at each time step. Conversely the IEKF has become slower due to the computation of Jacobians F and H; these become highly complex here due to the Pacejka combined slip tyre model. The 'best' Kalman filter in applications on other model structures and for on-line use will therefore depend on the relative complexity of the Jacobians. Provided the number of parameters to be identified is not excessive, the UKF may be quicker to run and is certainly simpler to implement.
In the specific context of on-line application, where the inputs would certainly be less aggressive than considered here and full parameter identification is no longer the objective, a subset of parameters could easily be adapted. For example, even without modification of the sampling time or model conplexity, the UKF will run faster than real-time on the midrange desktop PC used here; if the number of parameters is reduced from 14 to 8, 100 s of data is processed in 84 s.

Conclusions
Three alternative filtering methods have been adapted for system identification and applied to the identification of vehicle handling models. Although effective at identifiying nearoptimal parameters for a simplified model, the PF is very slow to converge and could not be used to carry out full vehicle and tyre identification.
Both UKF and EKF were found to be effective at identifying both simple and complex vehicle models. Although they use different methods for parameter error covariance estimation, both techniques have identical convergence characteristics and yield near-identical models. In applications such as this, where model Jacobians are complex and a relatively small number of model parameters is identified, the UKF is simpler to implement and slightly faster to run. Both Kalman filters can identify a simple linear handling model in under 5 min, while the featured four-dof full vehicle and tyre model takes around 10-15 h to converge fully, depending on the filter used.
A full set of the independent parameters of the nonlinear handling model were identified, including combined-slip tyre/suspension characteristics identified over their full range up to and beyond saturation. This was done using only commonly available CG acceleration and roll rate measurements. The selection of identified parameters and their resulting conditioning has also been considered. The choice of which parameters to identify and which to fix can be critical, but the convergence behaviour of the identified parameters in the running filter make it easy to determine unphysical and/or under-determined combinations.
The resulting full vehicle model has been tested on identification data from a test vehicle and also on independent validation test vehicle data. It was found to be very accurate in matching the measurements used for identification and also in explaining variables that are not readily measured, such as lateral velocity and roll angle. High quality results were seen over a range of validation tests, considering combined and separated longitudinal and lateral slip, at both high and low magnitudes.

Disclosure statement
No potential conflict of interest was reported by the authors.

Funding
This work is supported by Jaguar Land Rover and the UK-EPSRC grant EP/K014102/1 as part of the jointly funded Programme for Simulation Innovation (PSi).