*Article* **Multi-sensor Fusion Road Friction Coe**ffi**cient Estimation During Steering with Lyapunov Method**

**Letian Gao 1,2, Lu Xiong 1,2,\*, Xuefeng Lin 1,2, Xin Xia 1,2, Wei Liu 1,2, Yishi Lu 1,2 and Zhuoping Yu 1,2**


Received: 2 July 2019; Accepted: 3 September 2019; Published: 4 September 2019

**Abstract:** The road friction coefficient is a key parameter for autonomous vehicles and vehicle dynamic control. With the development of autonomous vehicles, increasingly, more environmental perception sensors are being installed on vehicles, which means that more information can be used to estimate the road friction coefficient. In this paper, a nonlinear observer aided by vehicle lateral displacement information for estimating the road friction coefficient is proposed. First, the tire brush model is modified to describe the tire characteristics more precisely in high friction conditions using tire test data. Then, on the basis of vehicle dynamics and a kinematic model, a nonlinear observer is designed, and the self-aligning torque of the wheel, lateral acceleration, and vehicle lateral displacement are used to estimate the road friction coefficient during steering. Finally, slalom tests and DLC (Double Line Change) tests in high friction conditions are conducted to verify the proposed estimation algorithm. Test results showed that the proposed method performs well during steering and the estimated road friction coefficient converges to the reference value rapidly.

**Keywords:** road friction coefficient; tire model; nonlinear observer; self-aligning torque; lateral displacement; Lyapunov method

#### **1. Introduction**

Vehicle safety-related state estimation [1–3] and control systems [4–6] have received much attention in recent decades. Vehicle dynamic control systems, such as the ASR (Anti-slip Regulation System), ESC (Electronic Stability Control), and AEB (Autonomous Emergency Brake), are realized by controlling the driving forces or braking forces so that the forces exerted by the road on the tires can be changed to maintain the stability of the wheels and the vehicle. The road friction coefficient is a key parameter for vehicle dynamic control systems [7,8], because it can reflect the dynamic motion limitations to a certain extent [9]. For human-operated vehicles, drivers can estimate the motion limitation of the vehicle and adapt their driving style using experience to prevent the vehicle from driving into critical conditions. However, with the development of intelligent vehicles, progressively more ADAS (Advanced Driving Assistant System) functions are being implemented by automated systems, which means that driving and braking forces and steering angles need to be calculated and controlled by control units, such as ACC (Adaptive Cruise Control) and LKA (Lane Keep Assistant). Therefore, an accurate road friction coefficient provides the automated system with the current motion limitation of the vehicle. Furthermore, for highly automated vehicles, the road friction coefficient is critical for decision-making, trajectory planning, and trajectory tracking.

Road friction coefficient estimation methods can be divided into two general types: cause-based methods and effect-based methods. The state-of-the-art methods of road friction estimation have been reviewed [10]. The principle of cause-based methods is the direct determination of road surface characteristics by special sensors, such as cameras, laser scanners, optical sensors, and so on. Alonso, J. et al. [11] proposed an asphalt status classification system based on real-time acoustic analysis of the tire-road interaction, but only wet and dry asphalt states were covered. Roychowdhury, S. et al. [12] proposed a two-stage method based on images captured by the front camera. A convolutional neural network model was applied to learn the road characteristics, and then the road states were divided into three types according to a rule-based strategy. Although cause-based methods can accurately characterize road states with special sensors, only a rough road friction coefficient is estimated, and the interval may vary within a very large range. The road friction coefficient describes the interaction effects between the road and tires [13], which means it cannot be estimated accurately by cause-based methods without considering the vehicle and tire characteristics. Effect-based methods estimate the road friction coefficient from the dynamic or kinematic motion responses [14] of the vehicle or wheels due to the tire forces caused by the interaction between the tires and the road. Normally, effect-based methods estimate the road friction coefficient by using models, and the estimation results are more accurate than those of cause-based methods. Effect-based methods fall into two categories: methods based on longitudinal dynamics and methods based on lateral dynamics. In order to obtain the road friction coefficient in all conditions, Ahn, C. et al. [15] divided driving conditions according to the slip ratio, the sideslip angle, and lateral acceleration. Then, different estimation methods were applied to estimate road friction coefficients under different conditions. Using longitudinal dynamics, Castillo Aguilar, J.J. et al. [16] applied a fuzzy logic algorithm to estimate the road friction according to the variation curve of the relationship between the road friction coefficient and the longitudinal slip ratio, and the algorithm was utilized in the hydraulic pressure control system of the EHB (Electric Hydraulic Brake) [17]. Enisz, K. et al. [18] designed an augmented vehicle model with the road friction coefficient and slip ratio, and the vehicle speed, wheel rotation speed, slip ratio, and road friction coefficient were simultaneously estimated by the EKF (Extended Kalman Filter). Taking advantage of the fact that the wheel torque of distributed drive electric vehicles is available and can be controlled precisely, Xia, X. et al. [19] proposed a road friction coefficient estimation algorithm under driving conditions using a nonlinearobserver. The observer performed well with strong longitudinal excitation, and its stability was proved. Moreover, many studies have focused on road friction coefficient estimation methods based on lateral dynamic characteristics. Using the relationship between lateral force and the sideslip angle, Wang, R. et al. [20] proposed a road friction coefficient estimation method that is effective when the vehicle has enough lateral excitation. Qi, Z. et al. [21] designed a Kalman Filter to estimate the longitudinal and lateral forces of each tire, as well as the derivatives of the two forces, using a 4 DOF(Degree of Freedom) vehicle model, and then the road friction coefficient was estimated according to the estimated tire lateral force. Compared with lateral force, self-aligning torque enters the nonlinear region earlier, so less lateral excitation is needed to estimate the road friction coefficient. Therefore, the relationship between the sideslip angle and self-aligning torque has been applied in many studies to obtain an accurate road friction coefficient. Luque, P. et al. [22] used a Kalman Filter to estimate longitudinal and lateral tire forces, and the self-aligning torque of the tire was calculated by a pretrained neural network. Then, the road friction coefficient was obtained by the relation curves between self-aligning torque and the sideslip angle in different road states. Matsuda, T. et al. [23] considered the road friction coefficient as a state and designed an EKF using a nonlinear 2 DOF single-track vehicle model, and self-aligning torque was measured to update the states. With varying road friction, Hsu, Y.-H.J. et al. [24] estimated the road friction coefficient from the relationships between (i) self-aligning torque and the tire trail and (ii) the tire trail and the sideslip angle. Ahn, C et al. [25] used a Kalman Filter to estimate the self-aligning torque of tires on the basis of the steering system and designed a nonlinear observer to estimate the road friction coefficient using self-aligning torque and a nonlinear vehicle model, and the stability and robustness of the nonlinear observer were proved. Shao, L. and Jin, C. [26] adopted a novel strategy to estimate the front axle lateral force. Then, combined with an indirect measurement based on total aligning torque estimation, a nonlinearadaptive observer was designed to estimate

the road friction coefficient with guaranteed stability. The self-aligning torque of the front axle is coupled with the sideslip angle; so, to precisely calculate self-aligning torque, we must know the current sideslip angles of each tire. Therefore, an accurate sideslip angle contributes to improvement in the estimation accuracy of the road friction coefficient using self-aligning torque-based estimation methods. With the development of intelligent vehicles, besides conventional onboard sensors, such as steering wheel angle sensors, wheel speed sensors, and inertial measurement units, information from new sensors equipped on intelligent vehicles can also be used to estimate the vehicle state. Yoon, J. and Peng, H. [27] used velocity measurements from two GPS receivers to estimate the sideslip angle. To reduce the cost, they took advantage of the direction measurement using a magnetometer and proposed a sideslip angle estimation method that integrated a magnetometer with a GPS [28]. Wang, Y. et al. [29] proposed a combined vehicle and vision model to increase the robustness of the body-slip-angle estimation to uncertain vehicle parameters, and multi-rate and time-delay issues were explained. Furthermore, camera-aided estimation of the lateral state for the integrated control of automated vehicles was discussed in Reference [30].

Since new sensors equipped on intelligent vehicles facilitate the estimation of vehicle states, they could be useful for improving the accuracy of the results of road friction coefficient estimation. In this paper, we introduced vehicle lateral displacement, which is based on the relationship between road friction and the self-aligning torque of the front axle, to the framework of road friction coefficient estimation. On the one hand, lateral displacement information contributes to improvement in the estimation accuracy of the vehicle's sideslip angle so that tire forces can be estimated more precisely. For intelligent vehicles, vehicle lateral displacement information can be obtained from cameras, GNSS (Global Navigation Satellite System) and maps, or V2X (Vehicle to Everything) systems. We acquired this information from a high-accuracy GNSS and a pre-established lane line map. On the other hand, compared with methods based on the relationship between road friction and the longitudinal or lateral forces of the tires, the self-aligning torque-based method requires fewer excitations, so the road friction coefficient can be estimated before the vehicle drives into critical conditions. We adjusted the tire brush model to fit the tire test data more accurately. Then, by integrating lateral displacement information, self-aligning torque measurement, lateral acceleration measurement, the tire model, and the vehicle model, we developed a nonlinear observer for road friction coefficient estimation. The stability of the observer was proved, and the observer's robustness was analyzed.

The main contributions of this paper are summarized as follows:


The remainder of this paper is organized as follows. In Section 2, the vehicle model and modified tire model are introduced. In Section 3, the nonlinear observer for road friction coefficient estimation is proposed, and its stability and robustness are analyzed. Section 4 presents experiments that were conducted to prove the proposed estimation method, and the experimental results are discussed. Finally, the conclusions are summarized in Section 5.

#### **2. Vehicle and Tire Model**

#### *2.1. Vehicle Model*

A nonlinear 2 DOF vehicle model is introduced to express the vehicle lateral dynamics, as shown in Figure 1. Both longitudinal and lateral load transfer are considered, and the dynamic model is expressed as:

$$\dot{\omega} = \frac{l\_f}{I\_z} (F\_{yfl} \cos \delta\_{fl} + F\_{yfr} \cos \delta\_{fr}) - \frac{l\_r}{I\_z} (F\_{yrl} + F\_{yrr}) \tag{1}$$

$$
\dot{\beta} = \frac{a\_y}{v\_x} - w
\tag{2}
$$

where ω is the yaw rate; *lf* and *lr* are the distance between the COG (Center of Gravity) and the front and rear axles, respectively; *Iz* is the vehicle yaw moment of inertia; *Fyfl*, *Fyfr*, *Fyr f* , and *Fyrr* are the lateral forces of the front left tire, front right tire, rear left tire, and rear right tire, respectively; δ*fl* and δ*f r* are the steering angles of the front left wheel and front right wheel, respectively; β is the sideslip angle; *ay* is the lateral acceleration of the vehicle; *vx* is the longitudinal speed of the vehicle.

**Figure 1.** Vehicle model.

Lateral force *Fy* is calculated by the tire model, and it corresponds to the road friction coefficient μ and sideslip angle α of each tire and vertical load *Fz*. The sideslip angles of tires are:

$$\begin{cases} \alpha\_{fl} = \frac{v\_{y} + l\_{f}\omega}{v\_{x} - \frac{b}{\tau}\omega} - \delta\_{fl} & \alpha\_{rl} = \frac{v\_{y} - l\_{r}\omega}{v\_{x} - \frac{b}{\tau}\omega} \\\ a\_{fr} = \frac{v\_{y} + l\_{f}\omega}{v\_{x} + \frac{b}{\tau}\omega} - \delta\_{fr} & \alpha\_{rr} = \frac{v\_{y} - l\_{r}\omega}{v\_{x} + \frac{b}{\tau}\omega} \end{cases} \tag{3}$$

where *vy* is the lateral speed of the vehicle, and *b* is track base. Given the load transfer, the vertical load of each tire is:

⎧ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨

⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩

$$\begin{cases} F\_{zf1} = mg \frac{l\_r}{2(l\_f + l\_r)} - ma\_X \frac{h}{2(l\_f + l\_r)} - ma\_y \frac{hl\_r}{b(l\_f + l\_r)}\\ F\_{zfr} = mg \frac{l\_r}{2(l\_f + l\_r)} - ma\_x \frac{h}{2(l\_f + l\_r)} + ma\_y \frac{hl\_r}{b(l\_f + l\_r)}\\ F\_{zf1} = mg \frac{l\_r}{2(l\_f + l\_r)} + ma\_x \frac{h}{2(l\_f + l\_r)} + ma\_y \frac{h}{b(l\_f + l\_r)}\\ F\_{zf1} = mg \frac{l\_r}{2(l\_f + l\_r)} + ma\_x \frac{h}{2(l\_f + l\_r)} - ma\_y \frac{hl\_f}{b(l\_f + l\_r)} \end{cases} \tag{4}$$

where *Fzfl*, *Fzfr*, *Fzrf* , and *Fzrr* are the vertical forces of the front left tire, front right tire, rear left tire, and rear right tire, respectively; *m* is the vehicle mass; *ax* is longitudinal acceleration; *h* is the height of the COG.

#### *2.2. Tire Model*

A novel modified tire brush model is applied to describe lateral force *Fy* and the self-aligning torque of the tire *Mz*. Compared with the traditional tire brush model, the proposed modified tire model has a simpler form and fits the tire test data better, so it is more convenient for road friction coefficient estimation during normal steering conditions. In the proposed modified tire model, the relationship between <sup>α</sup> *F*0.15 *z* and *Fy F*0.81 *z* and the relationship between <sup>α</sup> *F*0.45 *z* and *Mz F*1.85 *z* are established, as shown in Equations (5) and (6), respectively. With the new mapping relationships, the tire model can describe the variation in lateral force and self-alignment with different vertical loads more precisely.

$$\frac{F\_y}{F\_z^{0.81}} = -\text{sign}(a) \bullet 3d\_1\mu \frac{c\_1}{3\mu} \left(\frac{|a|}{F\_z^{0.15}}\right) \left\{ 1 - \frac{c\_1}{3\mu} \left(\frac{|a|}{F\_z^{0.15}}\right) + \frac{1}{3} \left[\frac{c\_1}{3\mu} \left(\frac{|a|}{F\_z^{0.15}}\right)\right]^2 \right\} \tag{5}$$

$$\frac{M\_z}{F\_z^{1.85}} = \text{sign}(\alpha) \bullet \frac{d\_2}{2} \mu \frac{c\_2}{3\mu} \left(\frac{|\alpha|}{F\_z^{0.45}}\right) \left[1 - \frac{c\_2}{3\mu} \left(\frac{|\alpha|}{F\_z^{0.45}}\right)\right]^3 \tag{6}$$

where *d*1, *d*2, *c*1, and *c*<sup>2</sup> are parameters, and μ is the road friction coefficient.

Tire tests were conducted on a tire test bench to verify the proposed tire model. The raw test data are shown in Figure 2. If we normalize the lateral tire force with the vertical load using the traditional tire brush model, the curves for different tire loads do not line up very well, as shown in Figure 3, which means that *Fy* and *Mz* are not directly correlative with *Fz*. With the new relationship proposed in the modified tire brush model, the test data are normalized, as shown in Figure 4. For lateral tire force, Figure 4a reveals that the relationships between <sup>α</sup> *F*0.15 *z* and *Fy F*0.81 *z* calculated by the modified tire model are nearly the same for different tire loads. Similarly, the test results in Figure 4b show that the relationships between <sup>α</sup> *F*0.45 *z* and *Mz F*1.85 *z* are the same for different vertical loads. Therefore, the proposed modified tire brush model can calculate the lateral tire force more precisely with tire load variation.

**Figure 2.** Row tire test data: (**a**) lateral force; (**b**) self-aligning torque.

**Figure 3.** Normalization of tire test data with the original tire brush model: (**a**) lateral force; (**b**) self-aligning torque.

**Figure 4.** Normalization of tire test data with the modified tire brush model: (**a**) lateral force; (**b**) self-aligning torque.

The tire test data in Figure 4 prove the function form of the modified tire model, so the next step is to use the test data normalized by the vertical load to fit the tire model function. The fitting results are shown as the black line in Figure 4, and the proposed modified tire model can be expressed as:

$$F\_{\mathcal{Y}}(a, F\_{z}, \mu) = -\text{sign}(a) \bullet 3d\_{1}\mu F\_{z}^{0.81} \frac{c\_{1}}{3\mu} \left(\frac{|a|}{F\_{z}^{0.15}}\right) \left\{1 - \frac{c\_{1}}{3\mu} \left(\frac{|a|}{F\_{z}^{0.15}}\right) + \frac{1}{3} \left[\frac{c\_{1}}{3\mu} \left(\frac{|a|}{F\_{z}^{0.15}}\right)\right]^{2}\right\} \tag{7}$$

$$M\_{\overline{z}}(\alpha, F\_{\overline{z}}, \mu) = \text{sign}(\alpha) \bullet \frac{d\_2}{2} \mu F\_{\overline{z}}^{1.85} \frac{c\_2}{3\mu} \left(\frac{|\alpha|}{F\_{\overline{z}}^{0.45}}\right) \left[1 - \frac{c\_2}{3\mu} \left(\frac{|\alpha|}{F\_{\overline{z}}^{0.45}}\right)\right]^3 \tag{8}$$

#### **3. Nonlinear Observer Design for Road Friction Coe**ffi**cient Estimation**

#### *3.1. NonlinearObserver Design*

Assuming that the road friction coefficient μ is piecewise constant and using the vehicle dynamic model (2) we have:

$$
\dot{\mu} = 0 \tag{9}
$$

$$
\dot{\upsilon}\_y = a\_y(\mu\_\prime \upsilon\_{y\prime} F\_z) - \alpha \upsilon\_x \tag{10}
$$

where *ay* μ, *vy*, *Fz* is lateral acceleration, which is:

$$\begin{array}{lcl}a\_{g}\{\mu,\nu\_{g},F\_{z}\} &= \frac{1}{m} \Big[\boldsymbol{F}\_{gf\boldsymbol{f}}\cos\delta\_{fl} + \boldsymbol{F}\_{gf\boldsymbol{f}}\cos\delta\_{fr} + \boldsymbol{F}\_{gf\boldsymbol{f}} + \boldsymbol{F}\_{g\boldsymbol{rr}}\\ &= \frac{1}{m} \Big[\boldsymbol{F}\_{\mathcal{Y}}\big(\mu,\alpha\_{f\boldsymbol{\alpha}},\boldsymbol{F}\_{z\boldsymbol{f}\boldsymbol{f}}\big)\cos\delta\_{fl} + \boldsymbol{F}\_{\mathcal{Y}}\big(\mu,\alpha\_{fr},\boldsymbol{F}\_{z\boldsymbol{f}\boldsymbol{f}}\big)\cos\delta\_{fr} + \boldsymbol{F}\_{\mathcal{Y}}\big(\mu,\alpha\_{fl},\boldsymbol{F}\_{z\boldsymbol{f}\boldsymbol{f}}\big) + \boldsymbol{F}\_{\mathcal{Y}}\big(\mu,\alpha\_{r\boldsymbol{f}},\boldsymbol{F}\_{z\boldsymbol{f}\boldsymbol{f}}\big)\Big] \end{array} \tag{11}$$

Since the lane line map information is available, we can measure (i) the distance between the COG of the vehicle and the lane line *yl* and (ii) the angle between the lane line and vehicle heading ϕ. *yl* and ϕ can either be obtained by a camera installed on the vehicle or calculated through the location information from a GNSS receiver and lane map, which is a priori knowledge. The lateral displacement can also be obtained through V2X technology; for example, through the UWB (Ultra-Wideband) localization technique, the distance between the vehicle and infrastructure along the road can be calculated. According to the kinematic relationships of vehicle motion, the dynamics of the distance between the COG and the lane line can be expressed as:

$$\dot{y}\_l = v\_\mathbf{x} \sin \varphi + v\_y \cos \varphi \tag{12}$$

From the system defined by Equations (10)–(12), the corresponding nonlinear observer isdesigned as: .

$$\dot{\hat{\mu}} = k\_1 \text{sign}(f\_1)[\mathbf{M}\_k - f\_k(\hat{\mu}, \hat{\alpha}, F\_z)] + k\_2 \text{sign}(f\_2)[a\_y - a\_y(\hat{\mu}, \hat{\alpha}, F\_z)] \tag{13}$$

$$\dot{\mathcal{g}}\_l = \upsilon\_x \sin \varrho + \vartheta\_y \cos \varrho + k\_3 (\mathcal{y}\_l - \mathcal{y}\_l) \tag{14}$$

$$\dot{\mathfrak{H}}\_{\mathcal{Y}} = a\_{\mathcal{Y}} \Big( \hat{\mathfrak{mu}}\_{\mathcal{Y}}, \hat{\mathfrak{y}}\_{\mathcal{Y}} F\_{\mathcal{Z}} \Big) - r v\_{\mathcal{X}} + k\_{4} (y\_{\mathcal{Y}} - \hat{\mathfrak{y}}\_{\mathcal{Y}}) + k\_{5} \Big[ a\_{\mathcal{Y}} - a\_{\mathcal{Y}} \Big( \hat{\mathfrak{mu}}, \hat{\mathfrak{y}}\_{\mathcal{Y}} F\_{\mathcal{Z}} \Big) \Big] + k\_{6} \int\_{0}^{t} \Big[ a\_{\mathcal{Y}} - a\_{\mathcal{Y}} \Big( \hat{\mathfrak{mu}}, \hat{\mathfrak{y}}\_{\mathcal{Y}} F\_{\mathcal{Z}} \Big) \Big] dt \tag{15}$$

where the superscript ^ denotes an estimated value; *k*1, *k*2, *k*3, *k*4, *k*5, and *k*<sup>6</sup> are parameters; *k*1, *k*2, *k*<sup>3</sup> and *k*<sup>6</sup> are positive, *k*<sup>4</sup> = cosϕ. From Equation (13)*f*<sup>1</sup> is defined as:

$$\begin{cases} f\_1 > 0 & \alpha\_f > 0, \hat{\alpha}\_f > 0 \\ f\_1 < 0 & \alpha\_f < 0, \hat{\alpha}\_f < 0 \\ f\_1 = 0 & \text{else} \end{cases} \tag{16}$$

and *f*<sup>2</sup> is defined as:

$$\begin{cases} f\_2 < 0 & \alpha\_f > 0, \left. \alpha\_f > 0, \alpha\_l > 0, \left. \alpha\_l > 0 \right. \\ f\_2 > 0 & \alpha\_f < 0, \left. \alpha\_f < 0, \alpha\_r < 0, \left. \alpha\_r < 0 \right. \\ f\_2 = 0 & \text{else} \end{cases} \right. \tag{17}$$

*Mk* is the self-aligning torque at the kingpin, and it can be calculated as:

$$M\_k = F\_{kl} L\_l(\delta\_{f1}) - F\_{kr} L\_r(\delta\_{fr}) \tag{18}$$

where *Fkl* and *Fkr* are tie rod forces on the left and right sides, respectively. *Ll*(δ*fl*) and *Lr*(δ*lr*) are the distances from the steering rods to the kingpin on the left and right sides, respectively. It has to be mentioned that, in our approach, *Fkl* and *Fkr* are measured by tension and compression force sensors installed at the left and right tie rods. This measurement can be provided by the steer-by-wire system of intelligent vehicles. It can also be estimated by the steering system if the algorithm of the EPS (Electric Power Steering) system is available. *fk* μˆ, *v*ˆ*y*, *Fz* is the self-aligning torque at the kingpin estimated by the wheel steering model and tire model, which can be expressed as:

*Sensors* **2019**, *19*, 3816

$$\begin{aligned} f\_k(\mu, \nu\_y, F\_z) &= M\_{zfl} + M\_{zfr} + l\_{\text{ml}} F\_{yfl} + l\_{\text{mr}} F\_{yfr} \\ &= M\_z(\mu, \alpha\_{fl}, F\_{zfl}) + M\_z(\mu, \alpha\_{fr}, F\_{zfr}) + l\_{\text{ml}} F\_y(\mu, \alpha\_{fl}, F\_{zfl}) + l\_{\text{mr}} F\_y(\mu, \alpha\_{fr}, F\_{zfr}) \end{aligned} \tag{19}$$

where *lml* and *lmr* are the mechanical trails of the left and right front tires, respectively.

#### *3.2. Stability Analysis*

According to the system dynamics in (9), (10), (12) and nonlinear observer in (13), (14), (15), the corresponding error dynamics is:

$$\dot{\tilde{\mu}} = -k\_1 \text{sign}(f\_1) \left[ f\_k(\mu, \upsilon\_{y^\*}, F\_z) - f\_k(\mu, \upsilon\_{y^\*}, F\_z) \right] - k\_2 \text{sign}(f\_2) \left[ a\_y(\mu, \upsilon\_{y^\*}, F\_z) - a\_y(\mu, \upsilon\_{y^\*}, F\_z) \right] \tag{20}$$

$$\dot{\overline{y}}\_l = -k\_3 \overline{y}\_l + \overline{v}\_y \cos \varphi \tag{21}$$

$$\begin{array}{rcl}\dot{\widetilde{\boldsymbol{w}}}\_{\mathcal{Y}} &=& (1-k\mathfrak{z}) \Big[a\_{\mathcal{Y}} \Big(\mu, \boldsymbol{\upsilon}\_{\mathcal{Y}}, \boldsymbol{F}\_{\mathcal{Z}}\big) - a\_{\mathcal{Y}} \Big(\hat{\mu}, \boldsymbol{\loz}\_{\mathcal{Y}}, \boldsymbol{F}\_{\mathcal{Z}}\big)\Big] - k\_{\mathcal{S}} \int\_{0}^{t} \Big[\Big[a\_{\mathcal{Y}} \big(\mu, \boldsymbol{\upsilon}\_{\mathcal{Y}}, \boldsymbol{F}\_{\mathcal{Z}}\big) - r\boldsymbol{\upsilon}\_{\mathcal{Z}}\big] - \Big[\mathbb{A}\_{\mathcal{Y}} \big(\hat{\mu}, \boldsymbol{\loz}\_{\mathcal{Y}}, \boldsymbol{F}\_{\mathcal{Z}}\big) - r\boldsymbol{\upsilon}\_{\mathcal{Z}}\big] \Big] dt - k\_{\mathcal{S}} \widetilde{y}\_{l} \\ &=& -k\_{\mathcal{S}} \widetilde{\boldsymbol{w}}\_{\mathcal{Y}} - (k\_{\mathcal{Y}} - 1) \Big[a\_{\mathcal{Y}} \big(\mu, \boldsymbol{\upsilon}\_{\mathcal{Y}}, \boldsymbol{F}\_{\mathcal{Z}}\big) - a\_{\mathcal{Y}} \big(\hat{\mu}, \boldsymbol{\upsilon}\_{\mathcal{Y}}, \boldsymbol{F}\_{\mathcal{Z}}\big) \Big] - k\_{\mathcal{A}} \widetilde{y}\_{l} \end{array} \tag{22}$$

where the superscript ~denotes the error between the estimated value and the real value.

The Lyapunov function is defined as:

$$V = \frac{1}{2}\widetilde{\mu}^2 + \frac{1}{2}\widetilde{y\_l}^2 + \frac{1}{2}\widetilde{v\_y}^2\tag{23}$$

Then, we have:

$$\begin{array}{l}\dot{V} = \widetilde{\mu} \Big[ -k\_{\rm l} \text{sign}(f\_{\rm l}) \Big] f\_{\rm l} \big( \mu, \upsilon\_{y}, F\_{z} \big) - f\_{\rm l} \big( \widehat{\mu}, \widehat{\upsilon}\_{y}, F\_{z} \big) \Big] - k\_{\rm l} \text{sign}(f\_{\rm l}) \Big[ a\_{\rm y} \big( \mu, \upsilon\_{y}, F\_{z} \big) - a\_{\rm y} \big( \widehat{\mu}, \widehat{\upsilon}\_{y}, F\_{z} \big) \Big] \\ + \widetilde{y}\_{\rm l} \big[ -k\_{\rm l} \widetilde{y}\_{\rm l} + (\cos \varphi) \widetilde{\upsilon}\_{y} \big] + \widetilde{\upsilon}\_{y} \big[ -k\_{\rm l} \widetilde{\upsilon}\_{y} - (k\_{\rm l} - 1) \big[ a\_{\rm y} \big( \mu, \upsilon\_{y}, F\_{z} \big) - a\_{\rm y} \big( \widehat{\mu}, \widehat{\upsilon}\_{y}, F\_{z} \big) \Big] - k\_{\rm l} \widetilde{y}\_{\rm l} \end{array} \tag{24}$$

Using the mean value theorem, we have:

$$f\_k(\mu, \upsilon\_{y\prime}, F\_z) - f\_k(\hat{\mu}, \upsilon\_{y\prime}, F\_z) = \frac{\partial \overline{f}\_k}{\partial \mu} \overline{\mu} + \frac{\partial \overline{f}\_k}{\partial \upsilon\_y} \overline{\upsilon}\_y \tag{25}$$

$$a\_y \left(\mu\_\prime \upsilon\_{y\prime} F\_z\right) - a\_y \left(\hat{\mu}\_\prime \upsilon\_{y\prime} F\_z\right) = \frac{\partial \overline{a}\_y}{\partial \mu} \overline{\mu} + \frac{\partial \overline{a}\_y}{\partial \upsilon\_y} \overline{\upsilon}\_y \tag{26}$$

where

$$\begin{cases} \frac{\partial \overline{f}\_k}{\partial \mu} = \frac{\partial f\_k}{\partial \mu} \left( \overline{\mu}, \overline{v}\_y \right) \\ \frac{\partial \overline{f}\_k}{\partial \overline{v}\_y} = \frac{\partial f\_k}{\partial v\_y} \left( \overline{\mu}, \overline{v}\_y \right) \\ \frac{\partial \overline{v}\_y}{\partial \mu} = \frac{\partial a\_y}{\partial \mu} \left( \overline{\mu}, \overline{v}\_y \right) \\ \frac{\partial \overline{a}\_y}{\partial \overline{v}\_y} = \frac{\partial a\_y}{\partial v\_y} \left( \overline{\mu}, \overline{v}\_y \right) \end{cases} \tag{27}$$

In Equation (27) μ is a median between μ and μˆ; *vy* is a median between *vy* and *v*ˆ*y*. Substituting (25), (26) and (27) into (24) we have:

$$\begin{split} \dot{V} &= -\begin{bmatrix} k\_1 \text{sign}(f\_1) \frac{\partial \overline{f}\_k}{\partial \mu} + k\_2 \text{sign}(f\_2) \frac{\partial \overline{u}\_y}{\partial \mu} \Big| \overline{\mu}^2 - k\_3 \overline{y}\_l^2 - \Big[ (k\_5 - 1) \frac{\partial \overline{u}\_y}{\partial v\_y} + k\_6 \right] \overline{v}\_y^2 \\\ -\begin{bmatrix} k\_1 \text{sign}(f\_1) \frac{\partial \overline{f}\_k}{\partial v\_y} + k\_2 \text{sign}(f\_2) \frac{\partial \overline{u}\_y}{\partial v\_y} + (k\_5 - 1) \frac{\partial \overline{u}\_y}{\partial \mu} \Big] \overline{\mu v\_y} + (-k\_4 + \cos \rho) \overline{y\_l} \overline{v\_y} \\\ -\begin{bmatrix} \ \overleftarrow{\mu} & \ \overleftarrow{y\_l} & \ \overleftarrow{v}\_y \end{bmatrix} \begin{bmatrix} \overline{\mu} \\ \overline{\overline{y}\_l} \\ \overline{\overline{v}\_y} \end{bmatrix} \end{split} \tag{28}$$

*Sensors* **2019**, *19*, 3816

where

$$A = \begin{bmatrix} A\_{11} & 0 & A\_{13} \\ 0 & k\_3 & A\_{23} \\ A\_{31} & A\_{32} & A\_{33} \end{bmatrix} \tag{29}$$

$$A\_{11} = k\_1 \text{sign}(f\_1) \frac{\partial \overline{f}\_k}{\partial \mu} + k\_2 \text{sign}(f\_2) \frac{\partial \overline{a}\_y}{\partial \mu} \tag{30}$$

$$A\_{13} = A\_{31} = \frac{1}{2}k\_1 \text{sign}(f\_1) \frac{\partial \overline{f}\_k}{\partial v\_y} + \frac{1}{2}k\_2 \text{sign}(f\_2) \frac{\partial \overline{a}\_y}{\partial v\_y} + \frac{1}{2}(k\_5 - 1)\frac{\partial \overline{a}\_y}{\partial \mu} \tag{31}$$

$$A\_{23} = A\_{32} = -\frac{1}{2}k\_4 + \frac{1}{2}\cos\varphi\tag{32}$$

$$A\_{33} = (k\_5 - 1)\frac{\partial \overline{a}\_y}{\partial v\_y} + k\_6 \tag{33}$$

If *<sup>A</sup>* is a positive definite matrix, then . *V* < 0 holds, which means that the estimation system is stable, and the estimation error will converge to zero as time *t* → ∞. To ensure that *A* is a positive definite matrix, all sequential principal minors of *A* need to be positive. According to the modified tire model in (7) and (8) and the symbolic function defined in (16) and (17) it can be deduced that:

$$\text{sign}(f\_1)\frac{\partial \overline{f}\_k}{\partial \mu} \ge 0 \tag{34}$$

$$\text{sign}(f\_2)\frac{\partial \overline{a}\_y}{\partial \mu} \ge 0 \tag{35}$$

Therefore, if we choose *k*1, *k*2, and *k*<sup>3</sup> as positive constants, the first-order and second-order sequential principal minors of *A* are positive. If we choose *k*<sup>4</sup> = cosϕ, then the third-order sequential principal minor of *A* is:

$$\begin{array}{ll}|A| = & k\_3 \Big[ k\_1 \text{sign}(f\_1) \frac{\partial \tilde{f}\_k}{\partial \mu} + k\_2 \text{sign}(f\_2) \frac{\partial \mathbb{E}\_y}{\partial \mu} \Big] \Big| (k\_5 - 1) \frac{\partial \mathbb{E}\_y}{\partial v\_y} + k\_6 \Big] \\ & - \frac{k\_3}{4} \Big[ k\_1 \text{sign}(f\_1) \frac{\partial \tilde{f}\_k}{\partial v\_y} + k\_2 \text{sign}(f\_2) \frac{\partial \mathbb{E}\_y}{\partial v\_y} + (k\_5 - 1) \frac{\partial \mathbb{E}\_y}{\partial \mu} \Big]^2 \end{array} \tag{36}$$

Since <sup>∂</sup> *<sup>f</sup> <sup>k</sup>* ∂*vy* , ∂*ay* ∂*vy* , and <sup>∂</sup>*ay* ∂μ are bounded according to the modified tire model, there exists a parameter *k*<sup>5</sup> such that:

$$\left[k\_1 \text{sign}(f\_1) \frac{\partial \overline{f}\_k}{\partial v\_y} + k\_2 \text{sign}(f\_2) \frac{\partial \overline{a}\_y}{\partial v\_y} + (k\_5 - 1) \frac{\partial \overline{a}\_y}{\partial \mu}\right]^2 = \left[A\_{11} + (k\_5 - 1) \frac{\partial \overline{a}\_y}{\partial \mu}\right]^2 = 0\tag{37}$$

If the chosen value of *<sup>k</sup>*<sup>6</sup> is large enough, then <sup>|</sup>*A*<sup>|</sup> <sup>&</sup>gt; 0 holds. Therefore, . *V* < 0 holds, and it can be deduced that:

$$
\begin{bmatrix}
\frac{\overline{\mu}}{\overline{y}\_l} \\
\frac{\overline{\nu}\_y}{\overline{\nu}\_y}
\end{bmatrix} \to \begin{bmatrix}
0 \\
0 \\
0
\end{bmatrix} \text{ as } t \to \infty \tag{38}
$$

and the estimation error can converge to zero exponentially.

#### *3.3. Robustness Analysis*

The uncertainties of the tire and vehicle models or measurements from sensors introduce perturbance to the system. It is necessary to analyze the performance of the estimator with bounded external excitation. According to the error dynamics of the system in (20), (21) and (22) without external inputs, the error dynamics of the system with inputs can be expressed as

$$\dot{\overline{\mu}} = -k\_1 \text{sign}(f\_1) \left( \frac{\partial \overline{f}\_k}{\partial \mu} \overline{\mu} + \frac{\partial \overline{f}\_k}{\partial v\_y} \overline{v\_y} \right) - k\_2 \text{sign}(f\_2) \left( \frac{\partial \overline{a}\_y}{\partial \mu} \overline{\mu} + \frac{\partial \overline{a}\_y}{\partial v\_y} \overline{v\_y} \right) + u\_1 \tag{39}$$

$$\dot{\overline{y}}\_l = -k\_3 \overline{y}\_l + \overline{v}\_y \cos q + \nu\_2 \tag{40}$$

$$\dot{\overline{\overline{\upsilon}}}\_{y} = -k\_{6}\overline{\overline{\upsilon}}\_{y} - (k\_{5} - 1) \Big( \frac{\partial \overline{\overline{a}}\_{y}}{\partial \mu} \overline{\mu} + \frac{\partial \overline{\overline{a}}\_{y}}{\partial \upsilon\_{y}} \overline{\overline{\upsilon}}\_{y} \Big) - k\_{4}\overline{y}\_{l} + \mu\_{3} \tag{41}$$

where *u*1, *u*2, and *u*<sup>3</sup> are external bounded inputs. *u*<sup>1</sup> includes the uncertainties of the tire model and the uncertainties in the estimation results of *vy*. *u*<sup>2</sup> includes the uncertainties in the measurements of lateral distance between the COG and the lane line and uncertainties in the estimation results of *vy*.*u*<sup>3</sup> includes the uncertainties in lateral distance measurements and uncertainties in the estimation results of μ.

The Lyapunov function is chosen and defined in Equation (23) thus, according to (39), (40), and (41) we have:

. *<sup>V</sup>* <sup>=</sup> μ 2 −*k*1sign(*f*1) ∂ *f <sup>k</sup>* ∂μ μ + <sup>∂</sup> *<sup>f</sup> <sup>k</sup>* ∂*vy vy* − *k*2sign(*f*2) ∂*ay* ∂μ <sup>μ</sup> + <sup>∂</sup>*ay* ∂*vy vy* + *u*<sup>1</sup> 3 +*yl* <sup>−</sup>*k*3*yl* <sup>+</sup> (cosϕ)*vy* + *u*<sup>2</sup> <sup>+</sup>*vy* 2 <sup>−</sup>*k*6*vy* − (*k*<sup>5</sup> − 1) ∂*ay* ∂μ <sup>μ</sup> + <sup>∂</sup>*ay* ∂*vy vy* <sup>−</sup> *<sup>k</sup>*4*yl* + *u*<sup>3</sup> 3 <sup>=</sup> <sup>−</sup>*A*11<sup>μ</sup><sup>2</sup> <sup>−</sup> *<sup>k</sup>*3*yl* <sup>2</sup> <sup>−</sup> *<sup>A</sup>*33*vy* 2 − 0 *A*<sup>11</sup> + (*k*<sup>5</sup> − 1) ∂*ay* ∂μ 1 μ*vy* <sup>+</sup> (cos<sup>ϕ</sup> <sup>−</sup> *<sup>k</sup>*4)*vyyl* <sup>+</sup> <sup>μ</sup>*u*<sup>1</sup> <sup>+</sup> *ylu*<sup>2</sup> <sup>+</sup>*vyu*<sup>3</sup> (42)

Since *k*<sup>4</sup> = cosϕ, substituting (37) into (42) we have

. *<sup>V</sup>* <sup>=</sup> <sup>−</sup>*A*11<sup>μ</sup><sup>2</sup> <sup>−</sup> *<sup>k</sup>*3*yl* <sup>2</sup> <sup>−</sup> *<sup>A</sup>*33*vy* <sup>2</sup> <sup>+</sup> <sup>μ</sup>*u*<sup>1</sup> <sup>+</sup> *ylu*<sup>2</sup> <sup>+</sup>*vyu*<sup>3</sup> ≤ −*A*<sup>11</sup> μ 2 − *k*<sup>3</sup> *yl* 2 − *A*<sup>33</sup> *vy* 2 + μ |*u*1<sup>|</sup> <sup>+</sup> *yl* |*u*2<sup>|</sup> <sup>+</sup> *vy* |*u*3<sup>|</sup> = 0 −*A*11(1 − θ1) μ 2 − *A*11θ<sup>1</sup> μ 2 + μ |*u*1<sup>|</sup> 1 + 0 −*k*3(1 − θ2) *yl* 2 − *k*3θ<sup>2</sup> *yl* 2 + *yl* |*u*2<sup>|</sup> 1 + 0 −*A*33(1 − θ3) *vy* 2 − *A*11θ<sup>1</sup> *vy* 2 + *vy* |*u*3<sup>|</sup> 1 = −*A*11(1 − θ1) μ 2 − *k*3(1 − θ2) *yl* 2 − *A*33(1 − θ3) *vy* 2 − *A*11θ<sup>1</sup> μ 2 − μ |*u*1<sup>|</sup> − *k*3θ<sup>2</sup> *yl* 2 − *yl* |*u*2<sup>|</sup> − *A*11θ<sup>1</sup> *vy* 2 − *vy* |*u*3<sup>|</sup> (43)

where 0 < θ<sup>1</sup> < 1, 0 < θ<sup>2</sup> < 1, and 0 < θ<sup>3</sup> < 1. Therefore, if the bounded inputs satisfy:

$$\begin{cases} |\mu\_1| < A\_{11}\theta\_1 \\ |\mu\_2| < k\_3\theta\_2 \\ |\mu\_3| < A\_{33}\theta\_3 \end{cases} \tag{44}$$

then we have: .

$$\dot{V} \le -A\_{11}(1-\theta\_1) \left| \overleftarrow{\mu} \right|^2 - k\_3(1-\theta\_2) \left| \overleftarrow{y} \right|^2 - A\_{33}(1-\theta\_3) \left| \overleftarrow{v}\_y \right|^2 \tag{45}$$

If we define:

$$A\_{11}(1 - \theta\_1) \left| \overleftarrow{\mu} \right|^2 + k\_3(1 - \theta\_2) \left| \overleftarrow{y} \right|^2 + A\_{33}(1 - \theta\_3) \left| \overleftarrow{v}\_y \right|^2 = \mathcal{W}(\mathbf{x}) \tag{46}$$

where *<sup>x</sup>* = <sup>μ</sup> *yl vy T* , then inequalities (47) (48) hold:

$$\frac{1}{4} \|\mathbf{x}\|^2 \le V \le \|\mathbf{x}\|^2\tag{47}$$

$$
\frac{\partial V}{\partial t} + \frac{\partial V}{\partial \mathbf{x}} \dot{\mathbf{x}} \le -\mathcal{W}(\mathbf{x}).\tag{48}
$$

By applying Theorem 4.19 from [31], we can reason that the system expressed in (39), (40), and (41) is input-state stable; thus, if the estimation system is interfered by bounded inputs, then the system will still stay stable.

#### **4. Experimental Validation**

Tests based on an electric vehicle were conducted to verify the proposed road friction coefficient estimation algorithm. The experimental setup and results are discussed in this section.

#### *4.1. Experimental Setup*

#### 4.1.1. Test Vehicle

The test vehicle is an electric vehicle and shown in Figure 5a, and the vehicle parameters are listed in Table 1. Information about the wheel speed and steering wheel angle was obtained through CAN-Bus. The GNSS receiver is Novatel 718D, which provides the absolute position and heading angle of the vehicle. It is necessary to point out that the lane lines were modeled in advance in the navigation coordinates so that the distance between the vehicle and lane line could be calculated in real time. ADIS16495 is the IMU (Inertial Measurement Unit), which measures the acceleration and angular velocities. The angle between the vehicle heading and lane line can be calculated by integrating the yaw rate of the vehicle. The steering tie rods are cut off and two Kistler tension and compression force sensors 9321B are installed on the left and right tie rods, respectively, as shown in Figure 5b. The tension and compression force sensors measure the force at the tie rods so that the self-aligning torque of the wheel can be measured indirectly.

**Figure 5.** Test vehicle implementation: (**a**) test vehicle; (**b**) steering tie rod with a tension and compression force sensor.


**Table 1.** Vehicle parameters.

#### 4.1.2. Test Road

To verify the proposed road friction coefficient estimation algorithm, slalom tests and DLC (Double Line Change) tests were conducted on the road shown in Figure 6. The white lane lines shown in Figure 6 were mapped in advance. According to a large number of emergency braking experiments, the real road friction coefficient is considered to be around 0.8.

**Figure 6.** Test road.

#### *4.2. Experimental Results and Analysis*

#### 4.2.1. Slalom Test

Slalom test results are shown in Figure 7. Figure 7a shows the vehicle speed measured by the GNSS receiver. The blue line and green line show the steering wheel angle and yaw rate during the test, respectively. Figure 7c shows the accelerations recorded by the IMU, and we can see that the maximum lateral acceleration reaches 8 *ms*<sup>−</sup>2, which means that the vehicle has enough lateral excitation. Figure 7d shows the estimated self-aligning torque at the kingpin according to the tension and compression force sensors installed at the steering tie rods. Figure 7e shows the estimated lateral distance between the COG of the vehicle and the lane line. The reference is calculated by the absolute position and the lane line map. The road friction coefficient estimation result is shown in Figure 7f. Since the road friction coefficient is related not only to the road surface but also to the tires, an absolutely precise road friction coefficient could not be obtained. From braking tests, we know that the real friction coefficient is about 0.8; therefore, we set 0.75–0.85 as the reference region. The initial road friction coefficient is set at 1. From Figure 8f, we can see that at around 9 s, the estimated road coefficient converges to the reference region, and the convergence time was about 3 s with continuous lateral excitation. After 9 s, the estimated value remains within the reference region, although there is a slight fluctuation, which means that the nonlinear estimator performs well. If the reference value of the road friction is considered as 0.8, then the estimation accuracy was about 97.2%. From 16 s onward, the vehicle drives straightly, and the road friction coefficient estimation algorithm stops without lateral excitation.

**Figure 7.** *Cont.*

**Figure 7.** Results of the slalom test: (**a**) vehicle speed; (**b**) steering wheel angle and yaw rate; (**c**) longitudinal and lateral acceleration; (**d**) self-aligning torque at the kingpin; (**e**) distance between the vehicle and the left lane line; (**f**) road friction coefficient.

**Figure 8.** *Cont.*

**Figure 8.** Results of the DLC (Double Line Change) test: (**a**) vehicle speed; (**b**) steering wheel angle and yaw rate; (**c**) longitudinal and lateral acceleration; (**d**) self-aligning torque at the kingpin; (**e**) distance between the vehicle and the left lane line; (**f**) road friction coefficient.

#### 4.2.2. DLC Test

Figure 8 shows the experimental results of DLC maneuvering. Figure 8a shows the vehicle speed during the test, and in Figure 8b represents the steering wheel angle and yaw rate of the vehicle by the blue line and green line, respectively. Figure 8c shows the longitudinal and lateral accelerations of the vehicle, and the maximum lateral acceleration is between 8 and 9 *ms*<sup>−</sup>2. The estimated aligning torques of the left and right kingpins are shown in Figure 8d. Figure 8e shows the estimated lateral distance between the vehicle and the lane line, and the estimated value tracks the reference value with little error. The road friction coefficient estimation results are shown in Figure 8f. If the reference value of the road friction is considered as 0.8, then the estimation accuracy was about 97.8%. Since the nonlinear estimator only works during steering, the estimation holds if the vehicle's lateral acceleration is relatively small during DLC, for example, from 7 to 8 s. Compared with the slalom test results, the road friction estimation results dose not fluctuate because the lateral excitation is not continuous. From the experimental results, we can see that the road friction coefficient rapidly converges to the reference value.

#### **5. Conclusions**

In this paper, a nonlinear observer for the road friction coefficient during steering based on the self-aligning torque characteristics of the tires aided by vehicle lateral displacement information was proposed. A modified tire brush model was established according to the tire test data, and the model describes the tire characteristics more precisely than the original model. A nonlinear observer using vehicle lateral displacement information was designed, and the stability and robustness were analyzed. Experiments were conducted to verify the proposed road friction coefficient estimation algorithm. The test results demonstrate that the proposed method performs well during vehicle steering, and the estimated road friction coefficient converges to the reference value very rapidly.

#### **6. Future Work**

We have modified the tire brush model according to the tire test data, and the results show that the modified model describes the tire characteristics properly. However, the tire tests were done with only one type tire in high friction condition, which is not sufficient to verify that the modified tire model is suitable for other tires with different sizes or types. Tire tests with more tires in different friction conditions should be conducted to verify the modified tire brush model.

The experiments were conducted only in high friction condition due to the test condition limitation. The algorithm should also be validated in low friction condition and high to low or low to high friction transition conditions in the future.

**Author Contributions:** methodology, L.G., L.X. and X.X.; software, L.G. and X.L.; validation, X.X., W.L. and Y.L.; resources, L.X. and Z.Y.; data curation, L.G. and X.X.; writing—original draft preparation, L.G.; writing—review and editing, W.L., X.X. and Y.L.; supervision, L.X. and Z.Y.; funding acquisition, L.X. and Z.Y.

**Funding:** This research is supported by the National Key Research and Development Program of China (grant no. 2016YFB0100901) and Shanghai Scientific Research Project (grant no. 16DZ1100700).

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

#### **References**


© 2019 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 (http://creativecommons.org/licenses/by/4.0/).
