**1. Introduction**

The automotive industry has made significant improvements in vehicle safety and driving performance in the last decades thanks to active control systems. These control systems rely on the measurement and estimation of several parameters and signals such as the wheel slip, sideslip angle. Tire normal forces can also be used to improve the vehicle safety and performance. Indeed, the longitudinal and lateral tire forces are coupled with the wheel loads [1]. Cho [2] showed that the estimate of tire forces, including tire loads, can be used to implement Global Chassis Control (GCC) systems and to further improve vehicle stability. Another practical application of tire normal force estimation is roll-over avoidance and understeering or oversteering prevention [3].

Tire loads are permanently changing when the vehicle is moving. Loads are transferred between wheels during accelerating, braking, and cornering. The position of the center of gravity, the road grade and irregularities on the road also impact the distribution of wheel normal forces making the estimation task a complex problem [4]. Due to the lack of low-cost sensors to measure the vehicle vertical force, a common approach is to consider the normal forces as constant parameters or to use an algebraic expression based on the vehicle longitudinal and lateral accelerations [5]. These open-loop estimation schemes, albeit simple, are not able to give a precise representation of the normal forces.

Doumiati et al. [6,7] provided a cascaded observer to estimate the tire normal forces. The first step of the algorithm estimates the lateral load transfer using a linear Kalman filter from the suspension deflection and accelerometer measurements and an estimate of the vehicle mass. A second observer is then used to infer the normal tire force from the lateral load transfer and a formulation of the normal force algebraic expression with coupling

**Citation:** Filipozzi, L.; Assadian, F.; Kuang, M.; Johri, R.; Velazquez Alcantar, J. Estimation of Tire Normal Forces including Suspension Dynamics. *Energies* **2021**, *14*, 2378. https://doi.org/10.3390/en14092378

Academic Editors: Wiseman Yair and Francesco Bellotti

Received: 3 March 2021 Accepted: 20 April 2021 Published: 22 April 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/).

between longitudinal and lateral acceleration. This yields an Extended Kalman Filter. Jiang [8] extended the application of this estimation framework by adding the vehicle pitch dynamics to take in account the road angle and the road irregularities.

Ozkan [9] used a Controller Output Observer (COO) to estimate lateral and normal tire forces. The estimation model is a half-car model which includes the vehicle lateral, heave, roll, and yaw motion. The model does not include pitch dynamics. Moreover, it assumes a constant vehicle velocity and a perfectly flat road. The COO has been successfully used in other automotive applications, e.g., to estimate the vehicle states [10] or to estimate longitudinal tire forces [11].

The goal of the paper is to provide a framework capable of delivering a reliable estimate of the wheel normal forces using low-calibration estimation methods. Contrary to previous work, the estimation of the wheel force will not be derived from algebraic expressions which estimate the normal forces from the vehicle acceleration. Instead, we intend to integrate the suspension dynamics in the normal force estimator to truly capture the tire force generation using suspension deflection sensors and accelerometers. This manuscript also illustrates the application of a recently developed estimation methodology called Youla Controller Output Observer (YCOO). This new approach is compared to the established method for signal estimation: Kalman filtering. Since the estimation problem is not formulated as a state estimation but an input estimation problem, the Unbiased Minimum Variance Filter (UMVF) is used in place of the standard state-estimation Kalman Filter. The comparison between the two estimation methods is based on several criteria: estimation performance, robustness to uncertainties, and ease of design. The performance is evaluated based on simulation results. In these simulations, we cover the major ways normal force can be generated (load transfer, disturbance from the road profile, and inclined road). The robustness is analyzed by conducting a sensitivity study of the suspension parameters, by introducing discrepancies between the estimation model and the simulation model, and by introducing noise in the measurements. Finally, by considering the different requirements to implement the YCOO and UMVF and the different approaches used by each observer (the UMVF uses a stochastic approach in the time domain while the YCOO relies on a deterministic approach in the frequency domain), we aim to highlight that the YCOO estimation framework is easily implementable with a low-calibration burden to guarantee robustness.

The following section analyzes a vehicle model and provides a model that can be used for closed-loop estimation. Sections 3 and 4 introduce the YCOO and UMVF estimation frameworks. As the normal force estimation is dependent on the vehicle mass estimation problem, Section 5 explains the computation of the vehicle mass estimate, which is used by the two frameworks. Finally, both the YCOO and UMVF are tested in simulation.

#### **2. System Modeling**

The generation of vertical forces of a rigid vehicle is linked with the compliance of the suspensions [12]. To develop the estimation frameworks and evaluate the effectiveness, a vehicle model with realistic heave, roll, and pitch dynamics is indispensable. Fortunately, the literature is filled with such vehicle models. Shim [13] described a 14 degrees of freedom vehicle model and validated it against the commercial vehicle models Carsim and ADAMS/Car. Figure 1 shows a schematic of the vehicle model. The degrees of freedom are the longitudinal, lateral, heave, roll, pitch, and yaw motion of the chassis, the vertical dynamics of each unsprung masses, and the wheels spin. The model assumes rigid bodies for the sprung and unsprung masses and neglects the compliance between the chassis and the unsprung mass in the vehicle longitudinal and lateral directions.

The model described in Figure 1 gives an accurate representation of the vehicle. However, it is highly complex which limits its usage to implement vehicle observers. Before beginning the observer development, a linear system analysis of the vehicle is performed. Its purpose is to determine the cross-coupling effects between the input–output

pairs to facilitate the observer design. The inputs of the model are the tire loads and the outputs are the suspension deflections.

**Figure 1.** Fourteen degrees of freedom vehicle model. The *x*, *y*, and *z* axes indicate, respectively, the longitudinal, lateral, and vertical directions of the vehicle. Parameters *ks*, *bs*, and *kt* denote the suspension stiffness, damping, and the tire stiffness.

The coupling between input-output pairs of the plant *G* is analyzed using the Relative Gain Array (RGA) <sup>Λ</sup> = (*G*(0)−1)*<sup>T</sup>* × *<sup>G</sup>*(0) where *<sup>G</sup>*(0) is the plant gain [14]. Since RGA is a linear analysis tool, the model mapping the tire loads *fFLz*, *fFRz*, and *fRLz*, *fRRz* to the suspension deflections *qFLs*, *qFRs*, *qRLs*, and *qRRs* must first be linearized. The operating point is chosen to be a steady-state cornering such that the vehicle velocity is *vx* = 90 km h−<sup>1</sup> and the vehicle lateral acceleration is *ay* = 0.4 g (i.e., the vehicle lateral velocity is *vy* = 0.2 m s<sup>−</sup>1; the heave velocity is *vz* = 0ms<sup>−</sup>1; and the roll, pitch, and roll angular velocities are *wx* = *wy* = 0 rad s−<sup>1</sup> and *wz* = 0.15 rad s<sup>−</sup>1). The matrices Λ<sup>1</sup> and Λ<sup>2</sup> correspond, respectively, to the RGA of a vehicle without an anti-roll bar and with an anti-roll bar on the rear axle. The inputs are *u<sup>T</sup>* = *fFLz fFRz fRLz fRRz* .


The matrix Λ<sup>1</sup> is almost equivalent to an identity matrix which indicates that there is almost no coupling between the four corners of the vehicle without anti-roll bar. Concerning the vehicle equipped with anti-roll bar, the two axles decoupled as the bottom left and upper right 2 × 2 submatrices of Λ<sup>2</sup> are almost zero. However, the anti-roll bar introduces a coupling between the left and right rear normal loads *fRLz* and *fRRz*, as can be seen in the bottom right corner of Λ2.

We assume that the vehicle is not equipped with an anti-roll bar. Thus, the wheel load of each corner of the vehicle can be estimated individually. A simpler model (Figure 2) is introduced to model the suspension of each wheel of the vehicle. The estimation model is a quarter-car model [15] with two inputs: a force applied on the sprung mass to represent the load transfer and another force which represents the tire load. In practice, the anti-roll bar of the vehicle should be considered and the coupling between the left and right normal forces should not be ignored. This coupling would appear as a load transfer from one side to the other when the two wheel loads are not equal, e.g., during cornering. The anti-roll bar would prevent the use of a quarter-car model and require a half-car model.

**Figure 2.** Quarter-car model of the suspension with wheel normal force and load transfer as inputs.

The bond graph of the model in Figure 2 yields the equation of motion of the quartercar model

$$
\dot{p}\_s = k\_s q\_s + b\_s \dot{q}\_s - m\_s g - \Delta f\_z \tag{2}
$$

$$
\dot{p}\_u = -k\_s q\_s - b\_s \dot{q}\_s - m\_u g + f\_z \tag{3}
$$

$$
\dot{q}\_s = \frac{p\_u}{m\_u} - \frac{p\_s}{m\_s} \tag{4}
$$

where *ps* and *pu* are the sprung mass and unsprung mass momentum, *qs* is the suspension deflection, *fz* is the wheel normal load, Δ*fz* represents the load transfer applied to the wheel, *ks* is the suspension stiffness, *bs* is the damper coefficient, and *ms* and *mu* are the sprung and unsprung masses of the corner of the vehicle. The states of the model are *x<sup>T</sup>* = *ps pu qs* . It is assumed that the vehicle is equipped with suspension deflection sensors and accelerometers on the sprung mass, hence the measurements are *y<sup>T</sup>* = " *p*˙*<sup>s</sup> ms qs* # . The model inputs are *u<sup>T</sup>* = *fz* Δ*fz* .

The equation of motions of the quarter-car model (Equations (2)–(4)) are not linear but affine. The constant term due to gravity can be eliminated by translating the states, inputs, and outputs of the system as follows

$$\begin{aligned} \text{i.i}^T &= \begin{bmatrix} p\_\ $ & p\_u & q\_\$  - \frac{m\_\xi q}{k\_s} \end{bmatrix}, \quad \text{ii}^T = \begin{bmatrix} f\_z - (m\_\xi + m\_u)g & \Delta f\_z \end{bmatrix}, \text{and } \mathfrak{j}^T = \begin{bmatrix} \frac{p\_\flat}{m\_\varkappa} - \mathcal{g} & q\_\flat - \frac{m\_\varkappa q}{k\_s} \end{bmatrix} \tag{5}$$

#### **3. Controller Output Observer**

The Youla Controller Output Observer, based on the COO framework [9], is a modelbased estimation technique that uses a controller to minimize the error between the measurement and the virtual measured signals of an estimation model. Contrary to the Luenberger observer, the YCOO does not assume that all system inputs are known; it is therefore well suited for an input estimation problem. Instead of designing a static gain controller via pole placement similarly to the Luenberger observer, the YCOO uses a dynamic controller designed with Youla parameterization [16]. This technique allows including information about the sensor dynamics and its noise content in the frequency domain to ensure good robustness and performance. A block diagram of the estimation concept is given in Figure 3. Measurements *y* are fed to the YCOO to provide an estimate *u*ˆ of the signal *u*. The YCOO is decomposed into two components: an estimation model *G*ˆ *<sup>p</sup>* that maps the estimated signal *u*ˆ to the virtual measurement *y*ˆ and a controller *Gc* that is responsible to follow those measurements [17].

The transfer function from the the true signal *u* and the estimated signal *u*ˆ is given by (*I* + *Lu*)−1*GcGp* where *Lu* = *GcG*ˆ *<sup>p</sup>* is the return ratio. If there is no discrepancy between the estimated plant and the actual one (*Gp* = *G*ˆ *<sup>p</sup>*), then this transfer function corresponds to the closed-loop transfer function *Tu*. If the plant has multiplicative uncertainty such that *Gp* = *G*ˆ *<sup>p</sup>*(*I* + Δ), the relation becomes *u*ˆ = (*Tu* + *Y*Δ)*u* where *Y* is the Youla transfer function. This shows that the YCOO relies on an accurate model of the system. Indeed, to guarantee good tracking of the measured quantities, *Tu* ≈ *I* is needed at low frequencies. This condition on *Tu* also constrains the gain of *Y* to the inverse of the plant gain at low frequencies since *Tu* = *YGp*. If the model is not correct, however, the term *Y*Δ introduces a steady-state error that cannot be compensated by the observer unless the plant gain is very high at low-frequency. Moreover, for the nominal system *Gp* = *G*ˆ *<sup>p</sup>*, the transfer function mapping the sensor noise to the estimation error is also given by *Y*. Thanks to its loop shaping approach, the YCOO directly addresses the trade-off between noise rejection, bandwidth, and robustness to high-frequency multiplicative uncertainties. Indeed, a higher bandwidth would increase the gain of *Y* at higher frequency, making the estimation more sensitive to noise and less robust to multiplicative uncertainties.

**Figure 3.** Block diagram of the YCOO estimation concept.

The quarter-car model described in the last section is used as the estimation model *G*ˆ *<sup>p</sup>*. The Youla parameterization technique is applied to design the controller *Gc* from the estimation model. The plant model can be written as a transfer function *G*ˆ *<sup>p</sup>* = *<sup>P</sup> <sup>δ</sup>* mapping the signals *u*ˆ to *y*ˆ.

$$\hat{\mathbf{G}}\_{\mathcal{P}} = \begin{bmatrix} \frac{p\_s}{f\_z} & \frac{p\_s}{\Delta f\_z} \\ \frac{q\_s}{f\_z} & \frac{q\_s}{\Delta f\_z} \end{bmatrix} = \frac{1}{m\_s m\_u s^2 + (m\_s + m\_u)b\_s s + (m\_s + m\_u)k\_s} \begin{bmatrix} k\_s + b\_s s & -(m\_u s^2 + b\_s s + k\_s) \\ m\_s & m\_u \end{bmatrix} \tag{6}$$

The first step in deriving a controller using the Youla parameterization technique is to find the Smith–McMillan form *MP* of the plant *Gp* such that *MP* = *ULGpUR*. The Smith–McMillan form [18] is useful in multi-variable control as it gives a realization of the plant in a basis where the plant is decoupled (i.e., its transfer function matrix is a diagonal matrix). The poles and transmission zeros of the Smith–McMillan form correspond to the poles and zero of the original system and the unimodular matrices *UL* and *UR* describe the transformation from the original basis to the basis used by the Smith–McMillan form. The Smith–McMillan form of the plant and its unimodular matrices *UL* and *UR* are

$$M\_P = \begin{bmatrix} \frac{1}{m\_s m\_u s^2 + (m\_s + m\_u)b\_s s + (m\_s + m\_u)k\_s} & 0\\ 0 & \frac{1}{m\_s m\_u} \end{bmatrix}, \\ \mathcal{U}\_L = \begin{bmatrix} 0 & 1\\ m\_u & m\_u s^2 + b\_s s + k\_s \end{bmatrix}, \\ \mathcal{U}\_R = \begin{bmatrix} 0 & \frac{1}{m\_s m\_u} \\ \frac{1}{m\_u} & -\frac{1}{m\_u^2} \end{bmatrix} \tag{7}$$

The controller is designed such that the decoupled system is a second-order Butterworth filter of unit gain with additional poles to make the controller proper, see Equation (8). The damping ratio *ζ* is set to <sup>√</sup><sup>1</sup> <sup>2</sup> as it offer good trade-off between fast transient and small oscillations. A large enough bandwidth is necessary for the wheel load estimate to be used by the control system. Moreover, the frequency response of a suspension mapping the road disturbance to tire force is shaped as a band-pass filter [12] whose high cutoff frequency is the wheelhop frequency (typically located at 10 Hz). Thus, to capture the tire force response, the closed-loop bandwidth should be faster than the wheelhop frequency. The controller is designed such that *ζ* = 1/ <sup>√</sup><sup>2</sup> and the bandwidth of the closed-loop system is 30 Hz. Singular values of the closed-loop transfer function and of the controller are given in Figure 4. At frequencies below the bandwidth, *Tu* is 0 dB and *Su* has low gain, ensuring a good tracking. At higher frequencies, the gain of *Tu* decreases to reject sensor noise and make the estimate robust against high-frequencies model mismatch.

$$M\_T = \frac{\omega\_0^2}{s^2 + 2\zeta\omega\_0 s + \omega\_0^2} \frac{1}{(\tau s + 1)^2} \tag{8}$$

Let *MY* such that *MT* = *MYMP*, the closed-loop transfer function and the controller transfer function matrix is obtained from the following equations which are derived from [17]

$$T\_{\rm II} = \mathcal{U}R\_{\rm R}\mathcal{U}\mathcal{U}\_{\rm R}^{-1} \tag{9}$$

$$S\_u = I - T\_u \tag{10}$$

$$Y = \mathcal{U}\_R M\_Y \mathcal{U}\_L \tag{11}$$

$$G\_{\mathfrak{C}} = S\_{\mathfrak{u}}^{-1} Y \tag{12}$$

This yields the controller

**Figure 4.** Singular values of the closed-loop transfer function *Tu*, *Su*, and *Y* and of the return ratio *Lu*, *Gc*, and *Gp*.

The transfer function from the measured signal *y* to the estimated input *u*ˆ is given

$$Y = \frac{\omega\_0^2}{(\tau s + 1)^2 (s^2 + 2\zeta\omega\_0 s + \omega\_0^2)} \begin{bmatrix} m\_{\text{ll}} & m\_{\text{ll}}s^2 + b\_{\text{s}}s + k\_{\text{s}} \\ -m\_{\text{s}} & b\_{\text{s}}s + k\_{\text{s}} \end{bmatrix} \tag{14}$$

Hence, the estimate of the normal force given by the YCOO is

$$\hat{f}\_z(s) = \left[ m\_{\text{il}}s \times \left( \frac{p\_s(s)}{m\_s} + sq\_s(s) \right) + b\_s sq\_s(s) + k\_s q\_s(s) \right] \frac{\omega\_0^2}{(s^2 + 2\zeta\omega\_0 s + \omega\_0^2)(\tau s + 1)^2} \tag{15}$$

The same equation can be obtained by combining (3) and (4) and by adding a filter with unit gain. Moschuk et al. [19] patented a concept to estimate wheel normal force using only suspension deflection sensors. The invention uses derivative filters to compute the suspension deflection velocity and the unsprung mass velocity (assuming the sprung mass vertical acceleration is null). It then uses damper and spring force maps to compute the tire normal force. Writing the derivative filter as *Fd*, the estimation in Reference [19] is

$$\hat{f}\_z(\mathbf{s}) = m\_\mathcal{U} F\_d \left( F\_d(q\_\mathcal{s}) \right) + b\_\mathcal{s} \left( F\_d(q\_\mathcal{s}) \right) + k\_\mathcal{s} (q\_\mathcal{s}) \tag{16}$$

This is similar to the estimate of the wheel normal force given by the YCOO in (15).

#### **4. Unbiased Minimum Variance Filtering**

The Unbiased Minimum Variance Filter is a variation of the Kalman filter for systems with unknown inputs. It gives an unbiased (zero-mean error) estimate of the model states and unknown inputs [20] with the assumption that the system is strongly observable. Consider the discrete Linear Time-Invariant (LTI) system

$$\begin{aligned} \mathbf{x}\_{k+1} &= A\mathbf{x}\_k + Bu\_k + He\_k\\ \mathbf{y}\_k &= \mathbf{C}\mathbf{x}\_k + Du\_k + Ge\_k \end{aligned} \tag{17}$$

where *xk* are the model states, *uk* the known inputs, and *ek* the unknown inputs. The states and inputs of the estimation model are *x<sup>T</sup> <sup>k</sup>* = " *ps pu qs* <sup>−</sup> *msg ks* # , *u<sup>T</sup> <sup>k</sup>* <sup>=</sup> <sup>∅</sup>, and *<sup>e</sup><sup>T</sup> <sup>k</sup>* = *fz* − (*ms* + *mu*)*g* Δ*fz* . The LTI system (17) is strongly observable if the matrix Ψ has full column rank [21] (*Gd* is the *G* matrix considering only feedthrough unknown inputs, i.e., with all zero-columns removed).

$$\Psi = \begin{bmatrix} \mathbf{C} & \mathbf{G} \\ \mathbf{CA} & \mathbf{CH} & \ddots \\ \vdots & \vdots & \ddots & \mathbf{G} \\ \mathbf{CA^{n-1}} & \mathbf{CA^{n-2}H} & \dots & \mathbf{CH} & \mathbf{G}\_d \end{bmatrix} \tag{18}$$

Since the first column of Ψ corresponds to the observability matrix, observability is a necessary condition for strong observability. Unfortunately, the system representing the quarter-car model with unknown force inputs (Equations (2)–(4)) is not observable. Indeed, similarity transformation shows that the state associated to the direction (*ms p*˜*<sup>s</sup>* + *mu p*˜*u*) does not produce any observable output. The system is reduced to eliminate the unobservable states. Moreover, it is necessary to use additional measurements to make the system strongly observable. The measured signals used by the UMVF is *y*˜*<sup>T</sup>* = " *p*˙*<sup>s</sup> ms* − *g qs* − *msg*/*ks q*˙*<sup>s</sup>* # . Note that suspension deflection sensors such as linear variable transformers which defines an electrical signal based on the position of an objected it is connected to can only measure deflection [22]. The measurement of the suspension relative velocity *q*˙*<sup>s</sup>* needed by the UMVF requires differentiating the signal *qs*, which requires additional signal processing.

Similar to the Kalman filter, the estimated states is computed in two steps. First, the estimated signals are computed based on the plant model.

$$\pounds\_{k+1\mid k} = A\pounds\_{k\mid k} + Bu\_k \tag{19}$$

$$P\_{k+1|k} = AP\_{k|k}A^T + Q\_k \tag{20}$$

Second, the gain *Lk*<sup>+</sup><sup>1</sup> is computed to guarantee an unbiased estimate of the model states.

$$\mathcal{R}\_{k+1} = \mathbb{C}P\_{k+1|k}\mathbb{C}^T + \mathbb{R}\_{k+1} \tag{21}$$

$$
\Phi\_{k+1} = \begin{bmatrix} -G & \mathcal{C}H \end{bmatrix} \tag{22}
$$

$$
\Omega\_{k+1} = \begin{bmatrix} 0\_{n \times p} & H \end{bmatrix} - P\_{k+1|k} \mathbb{C}^T \mathbb{R}\_{k+1}^{-1} \Phi\_{k+1} \tag{23}
$$

$$L\_{k+1} = P\_{k+1|k} \mathbb{C}^T \tilde{R}\_{k+1}^{-1} - \Omega\_{k+1} (\Phi\_{k+1}^T \tilde{R}\_{k+1}^{-1} \Phi\_{k+1})^{-1} \Phi\_{k+1}^T \tilde{R}\_{k+1}^{-1} \tag{24}$$

$$\pounds\_{k+1|k+1} = \pounds\_{k+1|k} + L\_{k+1}(y\_{k+1} - \mathbb{C}\pounds\_{k+1|k} - Du\_{k+1})\tag{25}$$

$$P\_{k+1|k+1} = L\_{k+1} \mathcal{R}\_{k+1} L\_{k+1}^T - P\_{k+1|k} \mathcal{C}^T L\_{k+1}^T - L\_{k+1} \mathcal{C} P\_{k+1|k}^T + P\_{k+1|k} \tag{26}$$

with *n* the number of states and *p* the number of unknown input *ek*. Given an unbiased estimate *x*ˆ*k*|, Palanthandalam-Madapusi [23] showed that an unbiased estimate of unknown inputs can be obtained from the following equations

$$\pounds\_k = H^\dagger L\_{k+1} (y\_{k+1} - \mathbb{C} \pounds\_{k+1|k} - D u\_{k+1}) \tag{27}$$

$$\pounds\_k = G^\dagger (y\_k - \text{Cat}\_{k|k} - Du\_k) \tag{28}$$

where † denotes the Moore–Penrose pseudo-inverse. Computing the unknown inputs using Equations (27) and (28) guarantees that *E*[*e*ˆ*k*] = *G*†*GE*[*ek*] and *E*[*e*ˆ*k*] = *H*†*HE*[*ek*], respectively. If *H* or *G* have full column rank, then this translates to *E*[*e*ˆ*k*] = *E*[*ek*]. In the general case where both *H* and *G* are not full column rank, it is necessary to combine Equations (27) and (28) to compute an unbiased estimate *e*ˆ*k*. In this article, we propose to compute the unknown input by solving a linear system. Let *V<sup>T</sup> <sup>H</sup>* and *<sup>V</sup><sup>T</sup> <sup>G</sup>* be the matrices of left eigenvectors associated to non-zero eigenvalues of *H*†*H* and *G*†*G*. The unknown input is the solution of:

$$
\begin{bmatrix} V\_H^T \\ V\_G^T \end{bmatrix} \hat{e}\_k = \begin{bmatrix} V\_H^T H^\dagger L\_{k+1} (y\_{k+1} - \mathbb{C} \pounds\_{k+1|k} - Du\_{k+1}) \\ V\_G^T G^\dagger (y\_k - \mathbb{C} \pounds\_{k|k} - Du\_k) \end{bmatrix} \tag{29}
$$

Without loss of generality, we can assume that rank *H<sup>T</sup> GT* = *p*, and it is possible to obtain at least *p* left eigenvectors of *H*†*H* and *G*†*G* associated to non-zero eigenvalues. Therefore, the matrix \$ *V<sup>T</sup> H V<sup>T</sup> G* % has full column rank. Taking the mean of (29) yields

$$
\begin{bmatrix} V\_H^T \\ V\_G^T \end{bmatrix} E[\hat{\mathbf{e}}\_k] = \begin{bmatrix} V\_H^T \\ V\_G^T \end{bmatrix} E[\mathbf{e}\_k] \tag{30}
$$

Since the matrix has full column rank, this guarantee that *E*[*e*ˆ*k*] = *E*[*ek*], i.e., *e*ˆ*<sup>k</sup>* is an unbiased estimate of *ek*.

#### **5. Vehicle Mass Estimation**

Both observers require knowledge of the vehicle mass and the location of the center of gravity. Indeed, the estimate of signals *u*˜ from the observer, as given in (5), directly depends on the vehicle mass. The static wheel load must be added to this estimate to compute the wheel load. This section presents a simple algorithm to obtain the static load of each wheel. More elaborate algorithms could be applied [24].

Algebraic expressions can be used to describe the wheel load distribution in quasisteady-state.

$$f\_{\rm ijz} = f\_{\rm ijz}^0 \pm \Delta f\_{\rm j}^x a\_{\rm x} \pm \Delta f\_{\rm i}^y a\_{\rm y} \; (i, j) \in \{F, L\} \times \{L, R\} \tag{31}$$

where the static load and the load transfer terms are

$$f\_{ijz}^{0} = \frac{m g (L - l\_j)(\mathcal{W} - w\_i)}{L \mathcal{W}} \tag{32}$$

$$
\Delta f\_j^x = \frac{mh(\mathcal{W} - w\_j)}{LW} \tag{33}
$$

$$
\Delta f\_i^y = \frac{mh(L - l\_i)}{L\mathcal{W}}\tag{34}
$$

Variables *wL* and *wR* denote the distance from the center of gravity of the vehicle to the left and ride sides; *W* = *wL* + *wR* is the track width; *lF* and *lR* denote the distance from the center of mass to the front and rear axles; *L* = *lF* + *lR* is the wheelbase; *m* is the total vehicle mass; *g* is the acceleration of gravity; *h* is the height of the center of mass; and *ax* and *ay* are the longitudinal and lateral acceleration due to vehicle acceleration, which also include the gravity component on the vehicle longitudinal and lateral axes.

The vehicle mass estimation algorithm is run when the vehicle is not moving. The wheel load *fijz* corresponds to the force produced by the suspension deflection ignoring the weight of the unsprung mass. Hence, at steady-state, the wheel load in Equation (31) is replaced by *fijz* = *ksqijs*. This yields

$$k\_{\bar{s}}q\_{\bar{i}\bar{j}\bar{z}} = f^{0}\_{\bar{i}\bar{j}\bar{z}} \pm \Delta f^{\bar{x}}\_{\bar{j}}a\_{\bar{x}} \pm \Delta f^{\bar{y}}\_{\bar{i}}a\_{\bar{y}} \tag{35}$$

where *ax* and *ay* are the longitudinal and lateral accelerations measured by the sensors and correspond to the acceleration of gravity in those directions if the vehicle is parked on a slope.

This corresponds to a system of four equations with four unknowns, *m*, *lF* , *wL*, and *h* (replacing *lR* and *wR* by *l* − *lF* and *w* − *wL* with *l* the vehicle wheelbase and *w* the axle track width). Thus, the position of the center of mass and the vehicle mass can be obtained when the vehicle is not moving.

#### **6. Simulation Results**

In this section, the YCOO and UMVF observers developed in Section 3 and 4 are tested in simulation. The full-car model shown in Figure 1 is assumed to represent the actual vehicle dynamics and the measured ground-truth signals are extracted from the 14 degrees of freedom vehicle model. The driving scenarios presented aim to cover all possible ways to redistribute tire loads, i.e., longitudinal or lateral load transfer, road irregularities, and sloped road. The results are also compared to the algebraic expression for normal force.

Figure 5a shows the estimates during a braking step of 3000 N m at 1 s from an initial velocity of 90 km h<sup>−</sup>1. The two observers provide better estimates than the algebraic expression which suffers from a steady-state error. Moreover, the two observers are intentionally not initialized; both observers converge in approximately 0.1 s. Figure 5b shows the estimate during a double lane change maneuver with a constant velocity of 90 km/h and with maximum lateral acceleration of 0.6 g. Both estimators provide a good estimate of the tire vertical force, whereas the estimation from algebraic expression does not capture the transient response. Figure 5a,b validates the two estimators for situations where the load transfer is due to longitudinal or lateral acceleration.

**Figure 5.** Vertical tire force estimation on maneuver with longitudinal and lateral only acceleration.

lateral acceleration.

Figure 6 shows the estimate during a bounce sine sweep test. The vehicle velocity is maintained at 20%. The road profile corresponds to sinusoidal bumps of decreasing wavelength with decreasing amplitude. The minimum wavelength is 1.6 m. Thus, the road excites the suspension over the frequency range 0 Hz to 3.5 Hz. The YCOO and the UMVF are able to estimate the wheel loads. Both observers reproduce the frequency response of the suspension: the wheel load amplitude increases when the road excitation get closer to the the suspension frequency (1 Hz obtained when *t* ≈ 13 s) and remains constant at frequencies between the suspension and wheelhop frequencies. Since the longitudinal and lateral accelerations during this maneuver are almost zero, the algebraic expression is not able to provide an accurate estimation of the wheel loads. The effect of road slope on the estimation scheme is investigated in Figure 7. The vehicle is driven from a flat road to a slope of 20% gradient. The YCOO and the UMVF capture the load transfer due to the road gradient and provide a good estimate during transient.

**Figure 6.** Vertical tire force estimation during a bounce sine sweep test. Vehicle speed is constant at 20 km h<sup>−</sup>1. The minimum wavelength is 1.6 m at *t* = 20%. The bottom figure shows the road profile.

**Figure 7.** Vertical tire force estimation when driving on a 20 % slope.

It is not practical to assume that the estimation model is a perfect representation of the actual suspension. Figure 8 evaluates the sensitivity of the wheel load estimation against the suspension stiffness. Uncertainties over this parameter result in an offset between the real and estimated wheel load. This is due to the wrong calibration of the mass estimation strategy. The load transfer estimate also suffers from uncertainties in the suspension stiffness. Indeed, without any uncertainty, both observers yield a correct load transfer of 700 N, but with a 50% stiffer suspension the load transfer estimate is only 450 N. The robustness against the damping coefficient *bs* is investigated in Figure 9a. The YCOO and the UMVF provide the same estimate, thus only the estimate given by the UMVF is shown in Figure 9a. The estimation is not robust against the damping coefficient in the transient but it does not affect the steady-state estimation. Similarly, nonlinearities in the damper map affect the transient of the wheel load estimate when the suspension operates in the region approximated by the linear damper map. The linear and nonlinear damper maps are given in Figure 9b.

**Figure 8.** Robustness against suspension stiffness during a braking maneuver. Solid lines show the ground-truth signals and dashed lines show the estimated ones.

(**a**) Robustness against damping coefficient (**top**) and damping map (**bottom**).

(**b**) Suspension linear and nonlinear damper map.

**Figure 9.** Robustness against uncertainties in the damping map during a braking maneuver. Solid lines show the ground-truth signals and dashed lines show the estimated ones.

Finally, Figure 10 shows the estimated signals obtained with the YCOO and the UMVF when Gaussian white noise of time correlation 10 ms and of power spectral density 10−<sup>4</sup> and 10−<sup>9</sup> is, respectively, added to the sprung mass vertical acceleration and to the suspension deflection measurements. The YCOO offers better noise rejection than the UMVF.

**Figure 10.** Estimation with noisy measurements.
