*Article* **Estimation of Vehicle Dynamic Parameters Based on the Two-Stage Estimation Method**

**Wenfei Li 1,2,3,\*, Huiyun Li 1,2,3, Kun Xu 1,2,3, Zhejun Huang 1,2,3, Ke Li 1,2,3 and Haiping Du <sup>4</sup>**


**Abstract:** Vehicle dynamic parameters are of vital importance to establish feasible vehicle models which are used to provide active controls and automated driving control. However, most vehicle dynamics parameters are difficult to obtain directly. In this paper, a new method, which requires only conventional sensors, is proposed to estimate vehicle dynamic parameters. The influence of vehicle dynamic parameters on vehicle dynamics often involves coupling. To solve the problem of coupling, a two-stage estimation method, consisting of multiple-models and the Unscented Kalman Filter, is proposed in this paper. During the first stage, the longitudinal vehicle dynamics model is used. Through vehicle acceleration/deceleration, this model can be used to estimate the distance between the vehicle centroid and vehicle front, the height of vehicle centroid and tire longitudinal stiffness. The estimated parameter can be used in the second stage. During the second stage, a single-track with roll dynamics vehicle model is adopted. By making vehicle continuous steering, this vehicle model can be used to estimate tire cornering stiffness, the vehicle moment of inertia around the yaw axis and the moment of inertia around the longitudinal axis. The simulation results show that the proposed method is effective and vehicle dynamic parameters can be well estimated.

**Keywords:** vehicle dynamic parameters; Unscented Kalman Filter; multiple-model

### **1. Introduction**

Nowadays, modern road vehicles are using an increasing number of active systems to improve vehicle safety, passenger comfort, vehicle performance and energy efficiency. Advanced Driver Assistance Systems (ADAS), as well as Automated Driving (AD) technologies, are being increasingly implemented in vehicles, aiming for improved driving safety and passenger comfort [1,2]. In addition, the autonomous driving test rig is also an important method to test autonomous driving control algorithms (as shown in Figure 1, it is an autonomous driving test rig proposed by our research group) [3–5]. The implementation of these fields greatly depends on accurate vehicle dynamic parameters. Vehicle dynamic parameters are also important for vehicle modeling. Thus, vehicle dynamic parameters are important for vehicle design and testing. The vehicle dynamic parameters (VDPs), such as the vehicle mass, moment of inertia and position of the vehicle centroid, affect the closed-loop behavior of active safety systems and play an important role [6]. It is necessary to determine the VDPs to obtain real vehicle responses. Some of the VDPs can be easily measured such as the mass, the track width or the wheelbase. However, other parameters are unknown and difficult to be measured directly, such as the distance from vehicle centroid to the front axis. The moment of inertia around each axis can be measured by

**Citation:** Li, W.; Li, H.; Xu, K.; Huang, Z.; Li, K.; Du, H. Estimation of Vehicle Dynamic Parameters Based on the Two-Stage Estimation Method. *Sensors* **2021**, *21*, 3711. https:// doi.org/10.3390/s21113711

Academic Editor: Felipe Jiménez

Received: 20 April 2021 Accepted: 22 May 2021 Published: 26 May 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

special equipment which is extremely costly. In contrast, the estimation method is a less intrusive and expensive way to obtain VDPs. The VDPs can be estimated by combining the estimation algorithm with some cheap sensors such as Inertial Measurement Unit (IMU), Global Positioning System (GPS), wheel speed sensors and steering angle sensor [7].

**Figure 1.** Autonomous driving test rig.

To obtain VDPs, many different methods have been proposed. In [8], a novel modelbased parameter identification approach using optimized excitation trajectory is proposed to identify the VDPs. However, this method needs test rigs, which is a huge cost. In addition, a variety of algorithms for VDPs estimation have been presented in works of literature [9–39]. The influence of VDPs on vehicle dynamics often involves coupling. Most of the papers only study the estimation part of some parameters and the other parameters are treated as being easily measured or obtained. In actual applications, this strategy is not feasible. In real applications, all VDPs need to be obtained through simple sensors and estimation strategies. Since the VDPs are always coupled with the vehicle states, the state-parameter joint and dual estimation methods [9,10] have become increasingly prevalent and have been studied by many researchers. Some researchers use the Dual Kalman Filter (DKF) to identify the VDPs and the vehicle states simultaneously. Besides, VDPs estimation is usually classified based on the parameters of interest and the vehicle dynamics model used. In [11], common onboard sensors which are able to measure the lateral acceleration and yaw rate and a non-linear vehicle model are used. Augmented Extended Kalman Filtering is used to estimate motion states and tire cornering stiffness based on a non-linear vehicle model and sensor. Sideslip and roll angles of electric are estimated using lateral tire force sensors through RLS and the Kalman Filter based on the Single-track model in [12]. Sprung mass, yaw moment of inertia and longitudinal position of the center of gravity are identified through a dual unscented Kalman Filter in [13]. In [14],

a four-wheel nonlinear vehicle model with roll dynamics and a correlation between the inertial parameters is used for a dual Unscented Kalman Filter to simultaneously identify the inertial parameters and the vehicle state. A local observability analysis on the nonlinear vehicle model is used to activate and deactivate different modes of the proposed algorithm. A Dual Extended Kalman Filter (DEKF) is used to estimate both vehicle states and vehicle parameters such as the vehicle mass, moment of inertia about the vertical axis and distance between the center of gravity and the front axle [15]. An extended Kalman Filter-based estimator adopting a dynamic vehicle model for determining the vehicle's longitudinal and lateral velocity as well as the yaw rate is proposed in [16]. In [17], a novel approach based on combined *H*∞ and extended Kalman Filter (*H*∞-EKF) is used to estimate the center of gravity position of electric vehicles. To implement this estimation algorithm, a simplified vehicle dynamics model is applied to the filter formulation. The *H*∞ estimator is employed to filter states by means of minimizing the influence of unexpected noise, whose statistics are unknown. Simultaneously, the other EKF estimator uses the states derived by the former filter to identify the position of the vehicle centroid. A methodology based on multiple-models and a switching method for real-time estimation of the position of vehicle centroid is proposed in [13]. The method uses the well-known simple linear vehicle models for lateral and roll dynamics and assumes the availability of lateral acceleration, the yaw rate, velocity, and steering angle measurements. As mentioned in previous research, the existing estimation methods are either expensive or only portions of the VDPs can be estimated. However, vehicle dynamics modeling needs to completely determine the completed VDPs, while the cost of VDPs acquisition should be as small as possible. Thus, a method that can obtain completed VDPs at low cost urgently needs to be proposed.

In order to obtain completed VDPs at a low cost, we propose a two-stage estimation method consisting of multiple-models and Unscented Kalman Filter to estimate VDPs. In the first stage, the vehicle is set to accelerate/decelerate and the longitudinal vehicle model is used. During this stage, the height of the vehicle centroid, tire longitudinal stiffness and the longitudinal position of the vehicle centroid are estimated by the Unscented Kalman Filter. After these parameters are estimated, these estimated parameters can be used in the second stage. In the second stage, a Single-track with roll dynamics vehicle model is adopted and the vehicle is set to continuous steering. Through vehicle steering, this model can be used to estimate tire cornering stiffness, the vehicle moment of inertia around the yaw axis and the moment of inertia around the longitudinal axis. After the two-stage estimation, all VDPs are estimated. The rest of the paper is organized as follows: vehicle dynamics model are shown in Section 2. The method used to estimate VDPs is provided in Section 3. Section 4 shows and discusses the simulation results. Finally, Section 5 delivers the conclusions and points towards future work.

#### **2. Vehicle Model**

The vehicle model used in this paper is a multiple-model approach which is based on a longitudinal vehicle dynamics model (as shown in Figure 2a) and a single-track with roll dynamics vehicle model (as shown in Figure 2b,c), which comprises: the motion in the longitudinal direction *x*, the longitudinal velocity; the motion in the lateral direction *y* or lateral velocity; the yaw around the vertical axis *z*, described by the yaw rate and roll with regard to the longitudinal axis *x*; and the roll rate [13]. Figure 2 illustrates the vehicle model adopted in this paper. The whole motion of the vehicle is a direct result of the forces (the aerodynamic forces and rolling resistance are neglected in this paper) that are generated between the road and tires. As shown in Figure 2b, the four-wheel vehicle dynamics model can be simplified as a single-track model. Other states that depend directly on these states can be derived, such as longitudinal and lateral accelerations. The tire states, such as the wheel slip angle, slip ratio and rotational velocities are also important. Tire-road friction force can be obtained based on tire states and tire stiffness. The vehicle states are also largely dependent on VDPs. VDPs include vehicle mass, moments of inertia around each axis and the position of the vehicle centroid.

**Figure 2.** Vehicle model: (**a**) Longitudinal vehicle dynamics model; (**b**) Single-track vehicle model; (**c**) Vehicle roll dynamics.

The vehicle dynamic model can be described by differential equations. The vehicle model implemented here can be obtained from [17,18]. When the vehicle was accelerating or decelerating along the longitudinal direction, the longitudinal vehicle dynamics model was adopted. As shown in Figure 2a, the longitudinal vehicle dynamics model was built with the longitudinal motion, as well as the front and rear wheel rotations

$$
\dot{m}\dot{v}\_x = F\_{xf} + F\_{xr} \tag{1}
$$

$$J\dot{\omega}\_f = T\_f - rF\_{xf} \tag{2a}$$

$$J\dot{\omega}\_r = T\_r - rF\_{\text{x}r} \tag{2b}$$

where *m* is vehicle total mass, *vx* represents the longitudinal vehicle velocity *Fxf* and *Fxr* represent the longitudinal forces of the front and rear tires. *J* is the wheel's moment of inertia. *r* is the equivalent radius of the front and rear tires. *Ti*, *ω<sup>i</sup>* (*i* = *f*,*r*) represent the wheel torque and angular speed. The load distribution can be expressed by the vertical forces that act on each of the four wheels. These can be calculated as follows:

$$F\_{zf} = mg\frac{l\_r}{L} - ma\_x \frac{h}{L} \tag{3a}$$

$$F\_{zr} = mg\frac{l\_f}{L} + ma\_x \frac{h}{L} \tag{3b}$$

where *Fzf* and *Fzr* are vertical force of the front and rear wheels. *ax* is the longitudinal accelerations, *g* is the gravitational constant, *lf* is the distance between the vehicle centroid and vehicle front axis, *lr* is the distance between the vehicle centroid and vehicle rear axis and *h* denotes the height of the vehicle centroid. *L* is the distance between the front axis and rear axis.

When the vehicle was being steered, the single-track with roll dynamics vehicle model was adopted. As shown in Figures 1c and 2b, the differential equations for the calculation of longitudinal and lateral acceleration are as follows:

.

$$
\dot{\boldsymbol{w}}\_{\boldsymbol{x}} = \boldsymbol{a}\_{\boldsymbol{x}} + \boldsymbol{v}\_{\boldsymbol{y}} \dot{\boldsymbol{\psi}} \tag{4}
$$

$$
\dot{\upsilon}\_y = a\_y + \upsilon\_\chi \dot{\psi} \tag{5}
$$

$$a\_{\chi} = \frac{1}{m} \left( F\_{\chi\emptyset} \cos \delta + F\_{y\emptyset} \sin \delta + F\_{xr} \right) \tag{6}$$

$$a\_y = \frac{1}{m} \left( F\_{xf} \sin \delta + F\_{yf} \cos \delta + F\_{yr} \right) \tag{7}$$

Yaw and roll motion can be obtained from:

$$
\ddot{\psi} = \frac{\Gamma}{I\_z} \tag{8}
$$

$$I\_x \ddot{\Phi} = mh \left( a\_y + \mathfrak{g}\phi \right) - k\_\phi \phi - c\_\phi \dot{\phi} \tag{9}$$

where . *<sup>ψ</sup>* is the yaw rate, . *φ* is the roll rate, *Iz* is the moment of inertia around the yaw axis, *Ix* is the moment of inertia around the longitudinal axis, *k<sup>φ</sup>* is the roll stiffness, *c<sup>φ</sup>* is the roll damping and *ay* is lateral acceleration. **Γ** can be calculated as follows:

$$
\Gamma = l\_f (F\_{xf} \cos \delta + F\_{yf} \sin \delta) - l\_l F\_{yr} \tag{10}
$$

where *δ* is the wheel steer angle while *Fyr* represents the lateral forces of the rear tires. There are many different approaches for achieving tire force, such as the so-called 'Magic Formula' by Pacejka [19], the tire model by Fiala [20] or the 'TMeasy' tyre model [21]. When the acceleration/deceleration strength of the vehicle is small and the steering angle is small, the tire force can be calculated as follows:

$$F\_{\rm xl} = C\_{\sigma} s\_l \tag{11a}$$

$$F\_{y\bar{l}} = \mathcal{C}\_a \mathfrak{a}\_{\bar{l}} \tag{11b}$$

where *Fxi*, *Fyi* (*i* = *f*,*r*) represent the longitudinal and lateral tire forces, *C<sup>α</sup>* denotes the tire cornering stiffness, *C<sup>σ</sup>* denotes the tire longitudinal stiffness, *si* is the slip ratio and *α<sup>i</sup>* is the slip angle. *α<sup>i</sup>* can be presented as follows:

$$\mathbf{a}\_f = \frac{\mathbf{v}\_y - l\_f \dot{\psi}}{\mathbf{v}\_x} - \boldsymbol{\mathcal{S}} \tag{12a}$$

$$\alpha\_r = \frac{v\_y - l\_r \dot{\psi}}{v\_x} \tag{12b}$$

the slip ratio *si* (*i* = *f*,*r*) can be presented as follows:

$$s\_l = \frac{\omega\_l r}{v\_\chi} - 1\tag{13}$$

When the acceleration/deceleration strength of the vehicle was small, the tire-road friction coefficient was proportional to the slip ratio rate [21]. Then the longitudinal tire force can also be presented as follows:

$$F\_{\rm xl} = F\_{\rm zl} \mathcal{C}\_{\rm K} \mathbf{s}\_{\rm l} \tag{14}$$

where *CK* is the slip ratio rate. It is a constant value related to the road surface. When the road surface was different, *CK* changed as well. From Equations (11a) and (14), it can be seen that tire longitudinal stiffness can be calculated based on the slip ratio rate and vertical force of the wheel. This means that the tire longitudinal stiffness can be obtained when the slip ratio rate is estimated.

#### **3. Estimation Method**

To adapt to non-linear problems in vehicle dynamics estimation, EKF is widely used for estimating different vehicle states. However, the accuracy of EKF-based estimation cannot be guaranteed due to linearization errors with Jacobian matrices when approximating non-linear systems [20–26]. More recently, additional attention has been paid to UKF estimation, which uses a set of sigma points to conduct non-linear transformation so that it can deal with strong non-linear estimation problems for vehicle dynamics systems [31]. The UKF, developed by Julier et al. [32] and refined by Wan and van der Merwe et al. [33] provides a new estimation approach. Unlike the EKF, the UKF approximates the probability density function of system states by implementing the Unscented Transformation (UT) instead of the system dynamics model. The UT captures the mean and covariance of the Gaussian random vector (GRV) to at least second-order accuracy through the use of a set of sample points. UKF is an effective method to estimate the states or the parameters of a discrete dynamic system. In this paper, we use UKF to estimate VDPs through a two-stage method. The frame diagram of the two-stage estimation method is shown in Figure 3. In this paper, we assume that the velocity of the vehicle can be measured by GPS, and the vehicle mass is known. The driving or braking torques of vehicle (*Ti*) can be obtained. Longitudinal acceleration *ax*, lateral acceleration *ay* and the yaw rate . *ψ* can be measured by IMU. Rolling stiffness *kφ* and roll damping *cφ* are given by the manufacturer. The relevant parameters are listed in Table 1.


**Table 1.** Nomenclature.

**Figure 3.** Two-stage estimation method.

As shown in Figure 3, the longitudinal vehicle dynamics model was adopted during the first stage. For the first stage, the longitudinal vehicle dynamics model could be described by Equations (1)–(3), (11a) and (14). To make sure the longitudinal vehicle dynamic model is able to reflect the real state of the vehicle, the absolute value of the frontwheel steering angle needed to be smaller than **0.62 deg** and the absolute value of yaw rate needed to be smaller than **1 deg/s** (as shown in Table 1). The inputs of the longitudinal vehicle dynamics model were wheel torque *Ti*(*i* = *f*,*r*). The states of the longitudinal vehicle dynamics model included the angular speed *ωi*(*i* = *f*,*r*) and longitudinal vehicle velocity *vx*. The measurable outputs were the angular speed *ωi*(*i* = *f*,*r*) and longitudinal vehicle velocity *vx*. As shown in Figure 3, the distance between the vehicle centroid and vehicle front *lf*, the height of vehicle centroid *h* and the tire longitudinal stiffness *CK* are estimated parameters. To enable the VDPs to be estimated (persistent excitation requirement), specific command signals needed to be given to activate the corresponding parameters (As shown in Table 1).

When *lf* , *h* and *CK* were estimated during the first stage, these estimated VDPs could be used in the second stage (as shown in Figure 3). During the second stage, a singletrack with roll dynamics vehicle model was used and is described by Equations (4)–(12). As shown in Figure 3, the inputs of the single-track with the roll dynamics vehicle model the wheel steer angle *<sup>δ</sup>*, lateral accelerations *ay*, yaw rate . *<sup>ψ</sup>* and the roll rate . *φ*. The states of this model include longitudinal vehicle velocity *vx*, lateral velocity *vy*, longitudinal accelerations *ax*, lateral acceleration *ay*, yaw rate . *<sup>ψ</sup>* and the roll rate . *φ*. The estimated parameters are the moment of inertia around the yaw axis *Iz*, the moment of inertia around the longitudinal axis *Ix* and the tire cornering stiffness *Cα*. When the conditions of the second stage in Table 2 are met, VDPs (*Iz*, *Ix*, *Cα*) can be estimated.

**Table 2.** Two-stage estimation condition requirements.


As shown in Figure 3, UKF was used in both the first stage and second stage. UT is one of the most important parts of UKF. First, we introduce UT here. UT is shown in Table 2 [29].

When a system function is given as *y* = *f*(*x*), *x* is the state and the dimension of *x* is *L* (as shown in Table 2). Given an L-dimensional GRV *x* with mean *x*ˆ and covariance *Px*, the statistics of *y* = *f*(*x*) were approximated by the selection of **2***L*+**1** discrete sample points {*χi*}**2***<sup>L</sup> <sup>i</sup>*=**<sup>0</sup>** <sup>=</sup> % *x and* ˆ *x*ˆ ± *σj*, *j* = **1**,..., *L* & where *σ<sup>j</sup>* is the *i th* column of the matrix )(*L* + *λ*)*Px*. *λ* is a scaling parameter and depends on *α*, *κ* and *L*. The constant *α* determines the spread of sigma points about the mean *x*ˆ. The constant *κ* is generally set to **3** − *L*. The constant *β* was used to incorporate prior knowledge of the distribution. In this paper, *α* = **0.01**, *β* = **2**.

As shown in Table 3, *ω* represents VDPs; *x* represents the states of dynamics; *d* represents the measured vector; *u* represents the input vector of the dynamic system. *R<sup>e</sup> <sup>k</sup>* is the measurement noise covariance. *R<sup>r</sup> <sup>k</sup>* is the processing noise covariance. The corresponding

parameters are shown in Table 4 and Section 4. At different stages, these variables represented different parameters. The local observability was demonstrated by investigating the rank of the observability matrix [20]. If the observability matrix had the full column rank, it was said to be locally observable. Using the continuous state-space representation, the discretized state-space representation was written by the Euler's forward discretization in (15).


**Table 4.** UKF for VDPs estimation.


$$d(\mathbf{x}\_k) = \mathbf{x}\_{k-1} + T\_s G(\mathbf{x}\_{k-1}, \omega\_{\mathbf{k}-1}, \mathbf{u}\_{k-1}) \tag{15}$$

where *Ts* is the sampling time. The observability matrix is the Jacobian of measurement vector *d*, with respect to the parameter vector *ω*. During the first stage, the Longitudinal vehicle dynamics model was used. According to Equations (1)–(3), (11a) and (14), the dynamic functions could be rewritten as:

$$m\dot{v}\_x = \left(m\lg\frac{l\_r}{L} - ma\_x\frac{h}{L}\right)\mathbf{C}\_K\mathbf{s}\_f + \left(m\lg\frac{l\_f}{L} + ma\_x\frac{h}{L}\right)\mathbf{C}\_K\mathbf{s}\_r\tag{16}$$

$$J\dot{\omega}\_f = T\_f - r \left(mg\frac{l\_r}{L} - ma\_x\frac{h}{L}\right) \mathbb{C}\_{K^\text{S}f} \tag{17}$$

$$J\dot{\omega}\_{r} = T\_{r} - r \left( m g \frac{l\_{f}}{L} + m a\_{x} \frac{h}{L} \right) \mathbb{C}\_{K^{\mathfrak{S}\_{T}}} \tag{18}$$

and the measurement vector *d* was written as:

$$\mathbf{d} = \mathbf{G}\_{\mathbf{1}}(\mathbf{x}, \omega\_{\mathbf{1}}, \mathbf{u}) = \begin{bmatrix} \mathbf{g}\_{\mathbf{1}1} \\ \mathbf{g}\_{\mathbf{1}2} \\ \mathbf{g}\_{\mathbf{1}3} \end{bmatrix} = \begin{bmatrix} \left( \mathbf{g}\_{\mathbf{1}}^{\mathsf{L}} - a\_{\mathbf{x}} \frac{\mathsf{h}}{\mathsf{L}} \right) \mathsf{C}\_{\mathsf{K}} \mathsf{s}\_{\mathsf{f}} + \left( \mathbf{g}\_{\mathbf{1}}^{\mathsf{L}} + a\_{\mathbf{x}} \frac{\mathsf{h}}{\mathsf{L}} \right) \mathsf{C}\_{\mathsf{K}} \mathsf{s}\_{\mathsf{r}} \\\ \frac{T\_{\mathsf{f}} - r \left( mg \frac{\mathsf{h}}{\mathsf{L}} - ma\_{\mathbf{x}} \frac{\mathsf{h}}{\mathsf{L}} \right) \mathsf{C}\_{\mathsf{K}} \mathsf{s}\_{\mathsf{f}}}{I} \\\ \frac{T\_{\mathsf{r}} - r \left( mg \frac{\mathsf{h}}{\mathsf{L}} + ma\_{\mathbf{x}} \frac{\mathsf{h}}{\mathsf{L}} \right) \mathsf{C}\_{\mathsf{K}} \mathsf{s}\_{\mathsf{r}}}{I} \end{bmatrix} \tag{19}$$

where *ω***<sup>1</sup>** represents a constant vector of the vehicle inertial parameters. During the first stage, *ω***<sup>1</sup>** = ' *lf h Ck* (*T* . The observability matrix was defined as the Jacobian of *G***<sup>1</sup>** with respect to the parameter vector *ω***1**. The Jacobian matrix, *C***<sup>1</sup>** = ∇*ω***1***G***1**, was represented as:

$$\mathbf{C}\_{1} = \nabla \,\omega\_{1} \,\mathbf{G}\_{1} = \begin{bmatrix} \frac{\partial \mathbf{g}\_{11}}{\partial \omega\_{11}} & \frac{\partial \mathbf{g}\_{11}}{\partial \omega\_{12}} & \frac{\partial \mathbf{g}\_{11}}{\partial \omega\_{13}} \\ \frac{\partial \mathbf{g}\_{12}}{\partial \omega\_{11}} & \frac{\partial \mathbf{g}\_{12}}{\partial \omega\_{12}} & \frac{\partial \mathbf{g}\_{12}}{\partial \omega\_{13}} \\ \frac{\partial \mathbf{g}\_{13}}{\partial \omega\_{11}} & \frac{\partial \mathbf{g}\_{13}}{\partial \omega\_{12}} & \frac{\partial \mathbf{g}\_{13}}{\partial \omega\_{13}} \end{bmatrix} \tag{20}$$

where

*∂g*<sup>11</sup> *∂ω*<sup>11</sup> *<sup>=</sup>*<sup>−</sup> *<sup>g</sup> <sup>L</sup>CKsf+ <sup>g</sup> <sup>L</sup>CKsr, ∂g*<sup>11</sup> *∂ω*<sup>12</sup> *<sup>=</sup>*−*ax <sup>L</sup> CKsf+ax <sup>L</sup> CKsr, ∂g*<sup>11</sup> *∂ω*<sup>13</sup> *=gsf*−*g lf <sup>L</sup>sf*−*ax h <sup>L</sup>sf+g lf <sup>L</sup>sr+ax h Lsr ∂g*<sup>12</sup> *<sup>∂</sup>ω*<sup>11</sup> *<sup>=</sup>rmg JL CKsf, ∂g*<sup>12</sup> *<sup>∂</sup>ω*<sup>12</sup> *<sup>=</sup>rmaxCKsf JL , ∂g*<sup>12</sup> *∂ω*<sup>13</sup> *<sup>=</sup>*−*rmg*<sup>1</sup> *<sup>J</sup> sf+rmg lf JL sf+maxr <sup>h</sup> JL sf ∂g*<sup>13</sup> *<sup>∂</sup>ω*<sup>11</sup> *<sup>=</sup>*−*rmgCKsr JL , ∂g*<sup>13</sup> *<sup>∂</sup>ω*<sup>12</sup> *<sup>=</sup>*−*rmaxCKsr JL , ∂g*<sup>13</sup> *∂ω*<sup>13</sup> *<sup>=</sup>*−*r***(***mg lr JL +max <sup>h</sup> JL* **)***sr*

Then *C*<sup>1</sup> could be written as:

$$\begin{array}{l} \mathbf{C}\_{1} = \nabla\_{\omega\_{1}} \mathbf{G}\_{1} = \\ \begin{bmatrix} -\frac{\mathbf{g}}{L} \mathbf{C}\_{K} \mathbf{s}\_{f} + \frac{\mathbf{g}}{L} \mathbf{C}\_{K} \mathbf{s}\_{r} & -\frac{a\_{\mathbf{x}}}{L} \mathbf{C}\_{K} \mathbf{s}\_{f} + \frac{a\_{\mathbf{x}}}{L} \mathbf{C}\_{K} \mathbf{s}\_{r} & \mathbf{g} \mathbf{s}\_{f} - \mathbf{g} \frac{\mathbf{l}\_{f}}{L} \mathbf{s}\_{f} - a\_{\mathbf{x}} \frac{\mathbf{l}\_{f}}{L} \mathbf{s}\_{f} + \mathbf{g} \frac{\mathbf{l}\_{f}}{L} \mathbf{s}\_{r} + a\_{\mathbf{x}} \frac{\mathbf{l}\_{f}}{L} \mathbf{s}\_{r} \\ \hline \frac{r \mathbf{m}\_{\mathbf{x}} \mathbf{C}\_{K} \mathbf{s}\_{f}}{L} & -r \mathbf{m}\_{\mathbf{x}} \mathbf{c}\_{\mathbf{x}} \mathbf{c}\_{f} & -r (\mathbf{m}\_{f} \mathbf{c}\_{f} + r \mathbf{m}\_{\mathbf{x}} \mathbf{c}\_{f}) \mathbf{s}\_{f} \\ \hline \end{bmatrix} \\ \end{array} \tag{21}$$

As shown in the above equation, the observability matrix was able to meet the requirement of the full column rank as long as the acceleration *ax* was properly selected.

During the second stage, the measurement vector consisted of the longitudinal vehicle velocity, yaw rate and roll rate. According to Equations (4)–(14), the dynamic functions could be rewritten as:

$$\dot{w}\_{\mathcal{X}} = \frac{1}{m} \left( \mathbb{C}\_{\mathcal{C}} s\_f \cos \delta + \mathbb{C}\_{\mathcal{A}} a\_f \sin \delta + \mathbb{C}\_{\mathcal{C}} s\_r \right) + v\_{\mathcal{Y}} \dot{\psi} \tag{22}$$

$$\ddot{\psi} = \frac{l\_f \left( \mathbb{C}\_v s\_f \cos \delta + \mathbb{C}\_u a\_f \sin \delta \right) - l\_r \mathbb{C}\_u a\_r}{l\_z} \tag{23}$$

$$\ddot{\phi} = \frac{h\left(\mathbb{C}\_{\upsilon}s\_f\sin\delta + \mathbb{C}\_{a}a\_f\cos\delta + \mathbb{C}\_{a}a\_r\right) + mgh\phi - k\_{\phi}\phi - c\_{\phi}\dot{\phi}}{I\_x} \tag{24}$$

then the measurement vector *d* was written as:

$$\begin{aligned} d &= G\_2(\mathbf{x}, \omega\_2) = \begin{bmatrix} \mathcal{S}\_{21} \\ \mathcal{S}\_{22} \\ \mathcal{S}\_{23} \end{bmatrix} \\ &= \begin{bmatrix} \frac{1}{m} (\mathcal{C}\_{\sigma} \mathbf{s}\_f \cos \delta + \mathcal{C}\_{\mathbf{a}} \mathbf{a}\_f \sin \delta + \mathcal{C}\_{\sigma} \mathbf{s}\_f) + \upsilon\_y \dot{\Psi} \\\ \frac{l\_f (\mathcal{C}\_{\sigma} \mathbf{s}\_f \cos \delta + \mathcal{C}\_{\mathbf{a}} \mathbf{a}\_f \sin \delta) - l\_r \mathcal{C}\_{\mathbf{a}} \mathbf{a}\_r}{I\_z} \\\ \frac{h (\mathcal{C}\_{\sigma} \mathbf{s}\_f \sin \delta + \mathcal{C}\_{\mathbf{a}} \mathbf{a}\_f \cos \delta + \mathcal{C}\_{\mathbf{a}} \mathbf{a}\_r) + \mathfrak{w} \mathbf{g} \mathbf{l} \phi - k\_\theta \boldsymbol{\phi} - \mathcal{c}\_\theta \dot{\phi}}{\dot{\Phi}} \end{bmatrix} \end{aligned} \tag{25}$$

where *ω*<sup>2</sup> represents a constant vector of the vehicle inertial parameters. During the second stage, *ω*2*=***[** *Iz Ix C<sup>α</sup>* **]** *T* . The observability matrix was defined as the Jacobian of *G*<sup>2</sup> with respect to the parameter vector *ω*2. The Jacobian matrix, *C*2*=*∇*ω*2*G*2, was represented as:

$$\mathbf{C}\_{2} = \nabla\_{\omega\nu\_{2}} \mathbf{G}\_{2} = \begin{bmatrix} \frac{\partial \mathbf{g}\_{21}}{\partial \omega \omega\_{21}} & \frac{\partial \mathbf{g}\_{21}}{\partial \omega \omega\_{22}} & \frac{\partial \mathbf{g}\_{21}}{\partial \omega \omega\_{23}} \\ \frac{\partial \mathbf{g}\_{22}}{\partial \omega \omega\_{21}} & \frac{\partial \mathbf{g}\_{22}}{\partial \omega \omega\_{22}} & \frac{\partial \mathbf{g}\_{22}}{\partial \omega \omega\_{23}} \\ \frac{\partial \mathbf{g}\_{23}}{\partial \omega \omega\_{21}} & \frac{\partial \mathbf{g}\_{23}}{\partial \omega \omega\_{22}} & \frac{\partial \mathbf{g}\_{23}}{\partial \omega \omega\_{23}} \end{bmatrix} \tag{2.6}$$

where

$$\frac{\partial g\_{21}}{\partial I\_z} = 0, \frac{\partial g\_{21}}{\partial I\_x} = 0, \frac{\partial g\_{21}}{\partial C\_a} = \frac{a\_f \sin \delta}{m}$$

$$\frac{\partial g\_{22}}{\partial I\_z} = -\frac{I\_f (C\_a, s\_f \cos \delta + C\_a a\_f \sin \delta) - I\_r C\_a a\_r}{I\_z^2}, \frac{\partial g\_{22}}{\partial I\_x} = 0, \frac{\partial g\_{22}}{\partial C\_a} = \frac{I\_f a\_f \sin \delta - I\_r a\_r}{I\_z}.$$

$$\frac{\partial g\_{21}}{\partial I\_z} = 0, \frac{\partial g\_{22}}{\partial I\_x} = -\frac{h (C\_a, s\_f \sin \delta + C\_a a\_f \cos \delta + C\_a a\_r) \sin \delta \exp(-k\_\phi \phi - c\_\phi \dot{\phi})}{I\_z^2}, \frac{\partial g\_{22}}{\partial C\_a}$$

$$= \frac{h a\_f \cos \delta + h a\_r}{I\_x}$$

Then *C*<sup>2</sup> could be written as:

$$\begin{aligned} \mathbf{C\_{2}} &= \nabla\_{\omega\_{2}} \mathbf{G\_{2}} \\ \mathbf{I} &= \begin{bmatrix} 0 & 0 & \frac{\mathbf{a}\_{r} \sin \delta}{m} \\ -\frac{\mathbf{l}\_{\parallel} (\mathbf{C\_{r}} \boldsymbol{\varphi} \cos \delta + \mathbf{C\_{a}} \boldsymbol{\varphi} \sin \delta) - \mathbf{l}\_{r} \mathbf{C\_{a}} \mathbf{a\_{r}}}{\mathbf{I}\_{z}^{2}} & 0 & \frac{\frac{\mathbf{a}\_{r} \sin \delta}{m} - \mathbf{l}\_{r} \mathbf{a\_{r}}}{\mathbf{I}\_{z}} \\ 0 & -\frac{\mathbf{l} (\mathbf{C\_{r}} \boldsymbol{\varphi} \sin \delta + \mathbf{C\_{a}} \boldsymbol{\varphi} \cos \delta + \mathbf{C\_{a}} \mathbf{a\_{r}}) + \mathbf{m} \mathbf{l} \mathbf{g} \boldsymbol{\phi} - \mathbf{k}\_{\theta} \boldsymbol{\phi} - \boldsymbol{\varepsilon}\_{\theta} \dot{\boldsymbol{\phi}}}{\mathbf{I}\_{z}^{2}} & \frac{\mathbf{h} \mathbf{a\_{r}} \cos \delta + \mathbf{h} \mathbf{a\_{r}}}{\mathbf{I}\_{z}} \end{bmatrix} \end{aligned} (27)$$

As shown in the above equation, the observability matrix was able to meet the requirement of full column rank as long as the steering angle *δ* was properly selected. Based on the above analysis, the VDPs can be estimated when the acceleration *ax* and steering angle *δ* are designed according to Table 2.

In order to compare the estimation performance of the method proposed in this paper, we used the commonly used extended Kalman algorithm for comparison. The extended Kalman algorithm is shown in Table 5.


The meanings of relevant parameters in the Table 5 are same as the meanings of relevant parameters in Table 4.

#### **4. Simulation Results**

The parameters of the vehicle model are shown in Table 6.

**Table 6.** Model parameters and definitions.


As shown in Figure 3, the whole parameter estimation process was divided into two parts. The second stage of the estimation could only start after the first stage of the estimation was completed. Due to space limitations, the control of the vehicle is not discussed here. The vehicle control can refer to reference [36–39]. First, we estimated *Iz*, *Ix*, *Cα*. As shown in Table 1, the vehicle needed to be continuously accelerated and braked for VDPs estimation. For this paper, the vehicle speed command signal was set as shown in Figure 4. It is a sinusoidal signal with a period of 12.5 s and an amplitude of 10. It includes acceleration/deceleration and can meet the requirements of the first stage (as shown in Table 1).

**Figure 4.** Vehicle longitudinal speed command.

VDPs can be estimated by making the car follow the command signal (as shown in Figure 4) to run for 4 cycles and 50 s. The simulation results are shown in Figure 5.

**Figure 5.** First stage VDPs estimation: (**a**) slip ratio rate estimation; (**b**) the height of vehicle centroid estimation; (**c**) estimation of the distance between the vehicle centroid and vehicle front axis.

As shown in Figure 5, the slip ratio rate *Cσ*,the height of vehicle centroid *h* and the distance between the vehicle centroid and vehicle front axis *lf* were estimated and the estimated values approximated the real values in a short time (about 10 s) through the proposed method in this paper. Compared with the UKF used in this paper, the estimation error of EKF was larger (as shown in Figure 5). When *Cσ*, *h* and *lf* were estimated, they were used in the second stage. To meet the requirement of the second stage (as shown in Table 2), the vehicle steering angle command signal was set as Figure 6.

**Figure 6.** Vehicle steering angle command signal.

As shown in Figure 6, the vehicle steering angle command signal is a sawtooth wave with a period of 10 s and an amplitude of 5. Additionally, the vehicle operated at a speed of 10 m/s. The simulation results are shown in Figure 7.

**Figure 7.** *Cont*.

**Figure 7.** Second stage VDPs estimation. (**a**) the tire cornering stiffness estimation (**b**) the moment of inertia around the longitudinal axis estimation (**c**) the moment of inertia around the yaw axis estimation.

As shown in Figure 7, tire cornering stiffness *Cα*, the moment of inertia around the longitudinal axis *Ix* and the moment of inertia around the yaw axis *Iy* can be well estimated and the estimated error is small through the method proposed by us. However, the estimation error by EKF becomes larger compared with the first stage (as shown in Figure 7). The simulation results show that the proposed method is very capable of estimating the VDPs, and thus proves the effectiveness of the proposed method. This is mainly caused by two reasons. First, the estimated parameters with larger errors in the first stage are used in the second stage. Second, the parameter estimation in the second stage is a non-linear estimation (This can be seen in Equation (27)). The two simulation results prove that the method proposed in this paper can more accurately estimate the VDPs.

#### **5. Discussion and Conclusions**

In this paper, a new method is proposed to estimate VDPs. Different from other studies that only estimated portions of VDPs, the proposed two-stage estimation method which combines multiple-models and the Unscented Kalman Filter is able to estimate more VDPs. Because the states of a vehicle are affected by the tire stiffness, the tire stiffness is difficult to measure. The proposed estimation method is able to estimate VDPs and tire stiffness. The proposed two-stage estimation method also solves the problem that VDPs have a coupling effect on vehicle motion, which makes the VDPs difficult to estimate. For comparison, EKF is used. The simulation results prove that the proposed method not only can estimate VDPs but also that the estimation errors are small.

The proposed two-stage estimation method in this paper can obtain all the VDPs needed for vehicle dynamics modeling at one time. It is useful for vehicle modeling, control and autonomous driving control algorithm tests on a test rig. More and more artificial intelligence technologies are being widely used in autonomous driving. However, most intelligent control algorithms are trained using vehicle kinematics models. An intelligent control algorithm trained with the kinematics model cannot accurately reflect the state of a real vehicle on the road. In order to ensure the effectiveness of the intelligent control algorithm, the vehicle dynamics model needs to be used in the algorithm training process. However, VDPs provided by most vehicle and devices manufacturers are not complete. The method proposed in this paper can estimate most of the VDPs required for vehicle dynamics modeling. Then it can be used to develop intelligent control algorithms for autonomous vehicles.

Our current research work verifies the effectiveness of the method proposed in this paper from the simulation. It verifies the program in advance for the next step of real vehicle test verification. In addition, it is assumed that some vehicle states can be measured directly in this paper. However, they are difficult to obtain in real scenarios. In the future work, we will use a Dual Unscented Kalman Filter to estimate the unmeasurable states and VDPs simultaneously.

**Author Contributions:** Conceptualization, W.L.; methodology, W.L.; validation, Z.H. and K.L.; formal analysis, K.L.; investigation, W.L.; writing—original draft preparation, W.L.; writing—review and editing, Z.H.; supervision, H.L. and H.D.; project administration, H.L.; funding acquisition, H.L. and K.X. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China, Grant No. 62073311 and the Key Program of Natural Science Foundation of Shenzhen, Grant Nos. JCYJ20200109 115403807, JCYJ20200109115414354, Science and Technology Development Fun, Macao S.A.R. (FDCT), No.0015/2019/AKP.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** This work was supported by grant from National Natural Science Foundation of China (Grant No. 62073311), the Key Program of Natural Science Foundation of Shenzhen (Grant Nos. JCYJ20200109115403807, JCYJ20200109115414354) and Science and Technology Development Fun, Macao S.A.R. (FDCT)(No.0015/2019/AKP). The authors would like to thank for support from National Natural Science Foundation of China, the Key Program of Natural Science Foundation of Shenzhen and Science and Technology Development Fun, Macao S.A.R. (FDCT).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

