Road profile estimation for suspension system based on the minimum model error criterion combined with a Kalman filter

2.1. Quarter vehicle suspension model

The quarter vehicle model, which is the most commonly used model for suspension systems, is presented in this section. Due to the nonlinearity of the practical quarter suspension, linearization of the quarter suspension model is used to study the quarter suspension. The linear quarter vehicle model contains the basic vertical information related to the vehicle [3, 19, 27]. The structure of a typical linear quarter vehicle model and road excitation model can be founded in [28, 29].

Based on these models, a linear quarter suspension dynamic equations of the model can be expressed as:

where $c$ represents the suspension damping. $k_s$ and $k_t$ represent the suspension spring stiffness and tire spring stiffness, respectively, and $m_b$, $m_w$ represents the suspension sprung mass and unsprung mass, respectively. The variables,$x_b$, ${\dot{x}}_b$ and ${\ddot{x}}_b$ are the displacement, velocity, and acceleration of the sprung mass. $x_w$, ${\dot{x}}_w$ and ${\ddot{x}}_w$ correspond to the displacement, velocity and acceleration of the unsprung mass, respectively. Finally, $x_r$ is the road unevenness, and the vibration of the suspension system is rooted in the road excitation.

The system state vector and output vector are chosen as:

where the six states variables are the displacement of the sprung mass, the velocity of sprung mass, the displacement of the unsprung mass, the velocity of the unsprung mass, road unevenness and the velocity of the road unevenness. In addition, the output variables correspond to the accelerations.

To apply the road unevenness profile in the estimation observer, it should satisfy the following [30]:

where $a_1$ and $a_2$ are constant values.

The state space equation of the state space variables can thus be expressed as:

where:

where $c$ represents the damping of the suspension, and $\mathbf{w}$ and $\mathbf{v}$ represent the process noise and measure noise, assumed to be independent and Gaussian. Then, $Q\left(k\right)=E\left(w{w}^T\right)$ and $R\left(k\right)=E\left(v{v}^T\right)$ are the associated process noise variance and measurement noise covariance, respectively.

A detailed derivation of the equations can also be found in the literature [31-34]. The discrete-time formulation of the state-space representation may be expressed as:

Substituting Eq. (2) into Eq. (5), and considering the model error, the system equation can be obtained as:

where the $k$ represents the discrete-time instant $kT\text{,}$ and $T$ is the time step. ${\mathbf{E}}_{\mathbf{s}}={\left\{{\epsilon }_1,{\epsilon }_2,{\epsilon }_3,{\epsilon }_4,{\epsilon }_5,{\epsilon }_6\right\}}^T$ is the error vector of the suspension system.

2.2. Augmented suspension system model

The uncertainty effects within suspension system result from the suspension spring characteristic, asymmetric velocity-force characteristics of suspension damping etc. These nonlinear behaviors are accounted for in the state-dependent damping nonlinearity $c\left(k\right)$ and varying spring mass $m_b\left(k\right)$. Since the estimation performance of the state observers for the vertical vehicle dynamics is particularly sensitive to deviations in the vehicle body mass [35, 36], variations in the sprung mass are the main focus in this paper. The equation of motion for the quarter vehicle suspension sprung mass is $m_b\left(k\right)=m_b\left(k\right)+∆m\left(k\right)$. Here $∆m\left(k\right)$ is calculated empirically.

Furthermore, to comprehensively investigate the effect of the KF & MME criterion on the road profile estimation for the suspension system, the linear quarter vehicle model is utilized. The Euler difference method is adopted to obtain the discrete version of a continuous system. Hence, the augmented model of the quarter-vehicle can then be expressed as:

where:

where $f_{1i}$ represents the $i$th component of the minimum and maximum sector bounds and $T$ is time step. $\mathbf{G}$ represents the error propagation matrix, and $\mathbf{d}$ represents the model error matrix (sprung mass error matrix). Referencing the system input, the system disturbance can be approximated using the Gaussian white noise process. Model uncertainties are assumed to be a Gaussian process as well. Each linear sub-model $i$ is computed using bounds and the suspension model parameters listed in Table 1.

The state estimation and observation processes of the discrete system of the 6-state variables and 2-output variables are illustrated in Fig. 1.

Fig. 1State estimation and observation processes of the suspension system

3.1. KF algorithm for road profile estimation

Based on Section 2, the problem of optimal estimation of ${\mathbf{x}}_k$ can be solved by minimization of the loss function [20]:

where ${\widehat{\mathbf{x}}}_{k/k-1}$ represents the prediction of ${\mathbf{x}}_k$, and ${\widehat{\mathbf{x}}}_{k/k}$ represents the prior estimation of ${\mathbf{x}}_k$. Details of the loss function can be found in [20].

A recursive estimation form of ${\mathbf{x}}_k$ may be expressed as:

where ${\widehat{\mathbf{y}}}_{k/k-1}$ denotes the prediction of ${\mathbf{y}}_k$, and ${\mathbf{K}}_k$ is the KF gain. The difference between ${\widehat{\mathbf{y}}}_k$ and ${\mathbf{y}}_k$ is called the filter innovation at the $k$th step. Assuming that the prior estimation ${\widehat{\mathbf{x}}}_{k-1/k-1}$ and the current observation ${\mathbf{y}}_k$ are Gaussian random variables, the optimal solution to the problem is given by the following procedure [20, 21]:

Initialization:

• The initial state and covariance are expressed as:

Time update:

• The prediction of the state and covariance are computed as:

Measurement update:

• The filter gain is determined as:

• The state estimation is deduced as:

• The estimated covariance is calculated as:

where ${\left(.\right)}^T$ denotes the transpose of a vector.

Based on the analysis above, the recursive form of the road profile estimation at the step $k$ may be given as:

where $y\left(k\right)-\widehat{y}(k-1)$ is the correction parameter in the estimation. Combining Eq. (10) and Eq. (14), the flow chart of KF is given in Fig. 2.

Fig. 2The flow chart of the KF algorithm

3.2. MME theory representation

The MME approach solves the system error model under the covariance constraint of the observation Eq. (16). The objective function of the MME compensator is expressed as [23, 26]:

where $\mathbf{W}$ represents the weight matrix of ${\mathbf{d}}_t$, and ${\widehat{\mathbf{y}}}_k$ is estimated using the last time step, and can be expressed as:

where $\partial h/\partial \mathbf{x}$ represents the Jacobi matrix of the observation equation; $L_{\mathbf{G}}\left(\mathbf{h}\right)$ and $T{L}_f\left(\mathbf{h}\right)$ are the 1st Lie derivations of $h\left(\mathbf{x}\right)$ about $\mathbf{G}$ and $f\left(\mathbf{x}\right)$, and are calculated as:

Using the minimum value principle, based on MME, the optimal sprung mass error can be given as:

The key of the MME criterion estimation is to properly choose the weight matrix $\mathbf{W}$. It can be expressed as:

where $m$ is the time scale for approaching $\mathbf{R}$.

Remark: It is difficult to choose $\mathbf{W}$ in practice. Due to the unpredictable white Gauss measurement noise, an ideal $\mathbf{W}$ is strongly time-dependent, unpredictable and even negative definite [26]. Such a result was adopted here, i.e. $\mathbf{W}$ was chosen as a constant.

Fig. 3The flow chart of the estimation process based on MME criterion and KF algorithm

3.3. KF algorithm based MME criterion

The chart of the proposed algorithm in this paper is depicted in Fig. 3. The Fig. 3 can be described as follows. First, the initial data, i.e. system covariance matrix, measurement data etc., are obtained from suspension system. Second, the error of suspension system is calculated using MME criterion after acquiring state estimation at step $k$. Then, combining the MME&KF algorithms, the state at step $k$ is updated to new state at $k+1$. Finally, the estimated information for the road profile at the new step is used as the initial state for the following step. Further details can be found in [23].

As shown in Fig. 3, the integrated estimation MME&KF approach tries to provide a more accurate suspension state estimation.

According to Eq. (6) and Eq. (7), the KF estimation algorithm for the discrete system may be summarized as:

where ${\mathbf{x}}_{k-1/k-1}$_, and ${\mathbf{P}}_{k-1/k-1}$are considered the optimal estimated vector and error covariance matrix of former time step; ${\mathbf{K}}_k^g$ is the KF gain matrix; $\mathbf{I}$ is unit matrix.

The estimation process from $k-1$ to $k$ detailed in Eq. (21) requires three types of information, i.e., the current input and observation information; the estimation results of the last step; and the system and measurement statistical information. Finally, the estimation returns the KF gain, system state and error covariance matrix for the current step.

4.1. Simulation results and analysis

In this work, we applied the MME&KF algorithm to estimate the state of the road profile for a suspension system. The ISO Level-B, ISO Level-C and ISO Level-D were taken as examples to illustrate the method. Also, the ISO level-B, ISO Level-C and ISO Level-D were calculated and used as the road excitation [29]. Note that it was assumed the tire did not lose contact with the ground [36-40].

The estimation accuracy of the MME&KF method proposed in this work was compared to the performance of KF under ISO level-B, ISO Level-C and ISO Level-D road excitations. Details are shown as follows.

4.1.1. Case 1: ISO level-B excitation simulation

The results for the simulation configuration based on an unchanged sprung mass are also illustrated under ISO level-B excitation in Fig. 4. Second, an additional body load $∆m$ ($∆m=$180 kg) was studied [36] under ISO level-B excitation in Fig. 5.

Figs. 5(a) and (b) illustrate the corresponding results of the KF and MME&KF results. Fig. 6(a) and (b) show the estimation error of the road profile for a varying sprung mass under the KF and MME&KF algorithms.

From Figs. 6(a) and (b), we can see that the estimation error of KF clearly fluctuates as the sprung mass changed under Level-B road excitation, and the MME&KF approach can improve the accuracy of the road profile for a varying sprung mass (model error).

Fig. 4Results of road level-B excitation (vveh= 40 km/h): road profile estimation of KF and MME&KF with mb unchanged (mb=mb)

Fig. 5Road profile comparison results of mb changed under road level-B excitation (vveh= 40 km/h)

a) KF and MME&KF achieved with altered $m_b$ ($m_b=m_b+∆m_b$)

b) KF and MME&KF achieved with altered $m_b$ ($m_b=m_b-∆m_b$)

4.1.2. Case 2: ISO level-C excitation simulation

The estimation accuracy of the MME&KF method proposed in this work was compared to the performance of the KF under ISO level-C excitation. The results for the simulation configuration based on an unchanged sprung mass are also illustrated under ISO level-C excitation in Fig. 7. Second, an additional body load $∆m$ ($∆m=$180 kg) was studied [36] under ISO level-C excitation in Fig. 8.

The corresponding results of the KF and MME&KF methods are illustrated in Fig. 7 and show that the higher estimation accuracy can be obtained with MME&KF. Figs. 8 and 9 show the estimation error of the road profile for a varying sprung mass under the KF and MME&KF algorithms. The simulation results showed that a smaller road excitation error is obtained using the proposed method.

Fig. 6Estimation error of KF and MME&KF road profile on road Level-B at vveh= 40 km/h

a)$m_b$ with altered ($m_b=m_b+∆m_b$)

b)$m_b$ with altered ($m_b=m_b-∆m_b$)

Fig. 7Results of road level-C excitation (vveh= 40 km/h): road profile estimation of KF and MME&KF with mb unchanged (mb=mb)

4.1.3. Case 3: ISO level-D excitation simulation

The estimation accuracy of the MME&KF method proposed in this work was compared to the performance of the KF under ISO level-D excitation. The results for the simulation configuration based on an unchanged sprung mass are also illustrated under ISO level-D excitation in Fig. 10. Second, an additional body load $∆m$ ($∆m=$180 kg) was studied [36] under ISO level-D excitation in Fig. 11.

Fig. 8Road profile estimation results of mb changing under road level-C excitation (vveh= 40 km/h)

a) KF and MME&KF achieved with altered $m_b$ ($m_b=m_b+∆m_b$)

b) KF and MME&KF achieved with altered $m_b$ ($m_b=m_b-∆m_b$)

Fig. 9Estimation error of KF and MME&KF road profile on road Level-C at vveh= 40 km/h

a)$m_b$ with altered ($m_b=m_b+∆m_b$)

b)$m_b$ with altered ($m_b=m_b-∆m_b$)

The corresponding results of the KF and MME&KF are illustrated in Fig. 10 and show that higher estimation accuracy can be obtained for MME&KF. Figs. 11 and 12 show the estimation error of the road profile for a varying sprung mass under the KF and MME&KF algorithms. The simulation results show that a smaller road excitation error is obtained using the proposed method.

Fig. 10Results of road level-D excitation (vveh= 40 km/h): road profile estimation of KF and MME& KF with mb unchanged (mb=mb)

Fig. 11Road profile estimation results of mb changing under road level-D excitation (vveh= 40 km/h)

a) KF and MME&KF achieved with altered $m_b$ ($m_b=m_b+∆m_b$)

b) KF and MME&KF achieved with altered $m_b$ ($m_b=m_b-∆m_b$)

In addition, the error values of the estimation standard deviation (STD) were calculated under ISO Level-B, ISO Level-C and ISO Level-D road excitation, and the simulation results are summarized in Table 2. The error of STD is higher with the KF estimation than MME&KF under the varying sprung mass.

Fig. 12Estimation error of KF and MME&KF road profile on road Level-D at vveh= 40 km/h

a)$m_b$ with altered ($m_b=m_b+∆m_b$)

b)$m_b$ with altered ($m_b=m_b-∆m_b$)

Table 2Calculation STD estimation values of different KF variations on road Level-B/C/D Profile at Vvel= 40 km/h

Filter state estimation	Road excitation	Filter mode	$m$ (kg)	Simulation results error STD (%)	Remarks
$x_r$	ISO level-B excitation	KF	$m_b$	7.3	$m_b$ represents the sprung mass of vehicle; $∆m$ represents the changed sprung mass of vehicle
		MME&KF	$m_b$	3.3
		KF	$m_b+∆m$	7.6
		MME&KF	$m_b+∆m$	3.7
		KF	$m_b-∆m$	10.3
		MME&KF	$m_b-∆m$	5.5
	ISO level-C excitation	KF	$m_b$	8.1
		MME&KF	$m_b$	4.6
		KF	$m_b+∆m$	9.5
		MME&KF	$m_b+∆m$	5.2
		KF	$m_b-∆m$	12.1
		MME&KF	$m_b-∆m$	8.3
	ISO level-D excitation	KF	$m_b$	10.5
		MME&KF	$m_b$	6.5
		KF	$m_b+∆m$	12.8
		MME&KF	$m_b+∆m$	7.9
		KF	$m_b-∆m$	16.4
		MME&KF	$m_b-∆m$	11.5

4.2. Experimental results and analysis

Because it is not easy to accurately obtain the time domain road profile, road profile is produced using the excitation equipment [36, 37].

The performance estimation of the KF&MME algorithm was conducted using the available test rig for a quarter suspension system, as pictured in Fig. 13. In the road profile estimation test process, road excitation force on the wheel was measured by a hydraulic ram, and sprung and unsprung mass accelerometers were installed to acquire data from the road excitation. During the experiments, the road excitation reference signal for the estimated quantity was computed off-line using a road excitation model [28, 29] for the test rig.

The road profile estimation accuracy of the proposed MME&KF algorithm was compared to the performance of the KF.

Fig. 13Quarter vehicle suspension test rig for road profile estimation

4.2.1. Case 1: ISO level-B excitation measured

Based the same simulation environment, the experimental and simulation results were compared and summarized in Figs. 14 and 15 under Level-B road excitation.

Fig. 14Road profile results of measured KF and MME&KF estimation on road Level-B at vveh= 40 km/h with mb unchanged (mb=mb)

Figs. 15(a) and (b) illustrate the test and estimation results of KF and MME&KF methods. Figs. 16(a) and (b) show the estimation error of the test road profile for a varying sprung mass under the KF and MME&KF algorithms. The worst estimation accuracy occurred when using the KF algorithm under the varying sprung mass condition, but higher estimation accuracy can be obtained using the MME&KF algorithm under both changing and without changing conditions. This demonstrates that the proposed MME&KF method is most effective under both the varying and without changing sprung mass conditions.

Fig. 15Road profile results of mb changed under road level-B excitation (vveh= 40 km/h)

a) KF and MME&KF achieved with altered $m_b$ ($m_b=m_b+∆m_b$)

b) KF and MME&KF achieved with altered $m_b$ ($m_b=m_b-∆m_b$)

Fig. 16Estimation error of KF and MME&KF road profile on road Level-B at vveh= 40 km/h

a)$m_b$ with altered ($m_b=m_b+∆m_b$)

b)$m_b$ with altered ($m_b=m_b-∆m_b$)

4.2.2. Case 2: ISO level-C excitation measured

The same experiment was performed under the Level-C road excitation. Results of the two estimation methods were compared and summarized in Figs. 17-19 under the Level-C road excitation. Fig. 17 shows the test and estimation results of the KF and MME&KF approaches. Figs. 18 and 19 illustrate that higher estimation accuracy can also be acquired using MME&KF algorithm under a changing and unchanging sprung mass under Level-C excitation conditions, and the proposed MME&KF algorithm is more effective than using only the KF algorithm.

Fig. 17Road profile results of measured KF and MME&KF estimation on road Level-C at vveh= 40 km/h with mb unchanged (mb=mb)

Fig. 18Road profile results of mb changed under road level-C excitation (vveh= 40 km/h)

a) KF and MME&KF achieved with altered $m_b$ ($m_b=m_b+∆m_b$)

b) KF and MME&KF achieved with altered $m_b$ ($m_b=m_b-∆m_b$)

Fig. 19Estimation error of KF and MME&KF road profile on road Level-C at vveh= 40 km/h

a)$m_b$ with altered ($m_b=m_b+∆m_b$)

b)$m_b$ with altered ($m_b=m_b-∆m_b$)

4.2.3. Case 3: ISO level-D excitation measured

The same experiment was performed under the Level-D road excitation. Results of the two estimation methods were compared and summarized in Figs. 20-22. Fig. 20 shows the test and estimation results of the KF and MME&KF approaches. Figs. 21 and 22 illustrate that higher estimation accuracy can be obtained using the MME&KF algorithm under a varying sprung mass or unchanging sprung mass under Level-D excitation conditions, and the proposed MME&KF algorithm is more effective than using only the KF algorithm.

Fig. 20Road profile results of measured KF and MME&KF estimation on road Level-D at vveh= 40 km/h with mb unchanged (mb=mb)

Fig. 21Road profile results of mb changed under road level-D excitation (vveh= 40 km/h)

a) KF and MME&KF achieved with altered $m_b$ ($m_b=m_b+∆m_b$)

b) KF and MME&KF achieved with altered $m_b$ ($m_b=m_b-∆m_b$)

Fig. 22Estimation error of KF and MME&KF road profile on road Level-D at vveh= 40 km/h

a)$m_b$ with altered ($m_b=m_b+∆m_b$)

b)$m_b$ with altered ($m_b=m_b-∆m_b$)

Moreover, the test error values of the STD are listed in Table 3. The error values of the estimation STD were calculated, and the test results are presented in Table 3. Comparing them to the simulation results, the same conclusion can be obtained under the ISO Level-B, ISO Level-C and ISO Level-D road excitations.

Table 3Calculation STD estimation values of different KF variations on road Level-B/C/D Profile at Vvel= 40 km/h

Filter state estimation	Road excitation	Filter mode	$m$ (kg)	Estimation results error STD (%)	Remarks
$x_r$	ISO level-B excitation	KF	$m_b$	7.8	$m_b$ represents the sprung mass of vehicle; $∆m$ represents the changed sprung mass of vehicle.
		MME&KF	$m_b$	4.5
		KF	$m_b+∆m$	8.5
		MME&KF	$m_b+∆m$	4.8
		KF	$m_b-∆m$	12.5
		MME&KF	$m_b-∆m$	6.8
	ISO level-C excitation	KF	$m_b$	9.2
		MME&KF	$m_b$	5.8
		KF	$m_b+∆m$	11.5
		MME&KF	$m_b+∆m$	7.8
		KF	$m_b-∆m$	15.4
		MME&KF	$m_b-∆m$	9.6
	ISO level-D excitation	KF	$m_b$	13.5
		MME&KF	$m_b$	8.5
		KF	$m_b+∆m$	16.8
		MME&KF	$m_b+∆m$	10.9
		KF	$m_b-∆m$	18.4
		MME&KF	$m_b-∆m$	14.5

Based on the analysis above, it can be seen that the estimation data from the MME&KF are closer to the measurement data. It may be concluded that more accurate road profile estimations can be achieved by employing the MME&KF algorithm. It should be noted that because the nonlinear tire characteristics are not considered during the process of road profile estimation; this may lead to a road profile estimation error under varying conditions of the tire.