Keywords: friction potential estimation, Kalman filters, parametric model learning, state-space representation, noise covariance matrix, MCMC estimation

The tire friction potential is the quantity characterizing the amount of friction remaining before the tire begins to skid on the road. Knowing this quantity during a travel turns out to be particularly advantageous, especially for the development of autonomous vehicles and also to significantly improve the driver assistance systems performances. Unfortunately, under standard driving conditions, the longitudinal tire force provided by a tire is low compared to the normal load applied on it. Therefore, the grip potential is particularly difficult to estimate under standard driving conditions. Thus, the goal of this PhD thesis is to carry out a method enabling to estimate the grip potential under standard driving conditions. This method should work using only the sensors fitted on production vehicle. The introduced method is divided into two main steps. Firstly, friction data points are combined with an adaptive Monte-Carlo Markov Chain method in order to predict the grip potential value. In this first main step, the friction data points are measurements coming from simulations or tests carried out on standard tires with a tire testing machine. However, in the practical case where the only accessible measurements are the ones provided by sensors fitted on production vehicle, these friction data points measurements are not available. Therefore, in a second main step, an extended Kalman filter is used in order to obtain friction data points estimates by using exclusively signals provided by the sensors fitted on production vehicle. This step is validated on simulated data representative of a straight line maneuver.

Résumé

Le potentiel d'adhérence d'un pneu est la grandeur caractérisant la quantité d'adhérence restante à un pneu avant que celui-ci ne se mette à patiner sur la route. La connaissance de cette grandeur lors d'un trajet pourrait se révéler particulièrement avantageuse, notamment pour le développement des véhicules autonomes et également pour améliorer sensiblement les performances des systèmes d'aide à la conduite. Malheureusement, lors d'une manoeuvre standard, les efforts longitudinaux générés par le pneu restent faibles par rapport à la charge subie par celui-ci. En conséquence, le potentiel d'adhérence est particulièrement difficile à estimer. Ainsi, l'objectif de ces travaux est de mettre en oeuvre une méthode permettant d'estimer le potentiel d'adhérence lors d'une manoeuvre standard en utilisant uniquement les capteurs équipant les véhicules de série. La méthode mise en oeuvre se scinde en deux grandes étapes. Tout d'abord, une méthode de type Monte-Carlo adaptative utilisant des points de friction est mise en oeuvre afin de déterminer le potentiel d'adhérence. Dans cette première grande partie, les points de friction utilisés sont des mesures provenant de simulations ou d'essais effectués avec des machines de test sur des pneus standards. Cependant, dans le cas pratique où les seules mesures accessibles sont celles fournies par les capteurs présents sur les véhicules de série, ces mesures de points de friction ne sont pas accessibles. Ainsi, dans une seconde grande étape, un filtre de Kalman étendu est utilisé afin d'obtenir des estimations des points de friction en ne se servant que des signaux fournis par les capteurs présents sur les véhicules de série. Cette étape est validée sur des données de simulation représentatives d'une manoeuvre en ligne droite.

Mots clés: estimation du potentiel d'adhérence, filtres de Kalman, estimation paramétrique, représentation d'état, matrice de covariance de bruit, estimation MCMC. vi humeur même si celle-ci aborde des sujets qui ne font pas forcément partie de ses domaines de prédilection. Je le remercie d'avoir pris le temps de faire un suivi de thèse très régulier au milieu d'un emploi du temps plus que bien rempli. A ce titre je lui témoigne mon admiration pour arriver à jongler entre tous ses projets et être toujours partant pour ''faire de la science" entre deux présentations PowerPoint. Je lui exprime également ma gratitude pour toutes les connaissances qu'il m'a transmises avec passion et les discussions toujours intéressantes. Enfin je tiens à lui présenter mes excuses si quelque fois je l'ai parfois fait rester tard au bureau ou si mes rédactions lui ont demandé des relectures nocturnes.

Je remercie également chaleureusement Jérémy Vayssettes, sans qui cette thèse n'aurait tout simplement pas vu le jour. Bien qu'il n'ait pas pu co-encadrer ma thèse, il a toujours su se montrer disponible quand j'ai rencontré des difficultés et il a toujours pris le temps de répondre à mes questions de façon détaillé. Je le remercie également pour sa générosité et ses mots d'encouragement à chacune de nos rencontres. Savoir que je pouvais compter sur lui en cas d'obstacle majeur a toujours été une source de réconfort.

Je tiens également à remercier l'équipe ETN qui m'a accueilli dans de bonnes conditions pour mener à bien mes travaux de thèse. Au cours de mon doctorat, j'ai eu la chance de passer quelque mois au laboratoire le temps de prodiguer des enseignements aux élèves de l'ENSI Poitiers. Je profite donc de ces pages pour adresser mes remerciements à l'ensemble du personnel du LIAS pour leur accueil, leur disponibilité et leur bonne humeur.

Finalement, bien qu'un travail de thèse soit un travail professionnel, il me semble impossible de le mener à bien sans un équilibre personnel épanouissant. C'est pourquoi, je ne peux terminer ces remerciements sans remercier mes proches à commencer par mes parents et mes frères pour leur confiance inébranlable. Ils ont toujours été à mes côtés dans les bons moments comme dans les tempêtes. Comment ne pas citer également celle qui m'accompagne depuis quelques années déjà. Margaux, ces moments passés avec toi sont les plus précieux qui soient. Je te remercie pour ton soutien indéfectible et la joie de vivre que tu me transmets au quotidien.

CHAPTER 1 Introduction and state of the art 1.1 Introduction

Because traffic accidents are responsible for a huge number of deaths every year all over the world, the passenger's safety is one of the priority for vehicle manufacturers [START_REF] Boufadene | Vehicle online parameter estimation using a nonlinear adaptive observer[END_REF][START_REF] Boudali | Unified dynamic and geometrical vehicle guidance strategy to cope with the discontinuous reference trajectory[END_REF]. In order to ensure this safety, passenger vehicles are nowadays equipped with control chassis systems such as the anti-lock braking system (ABS) or the electronic stability programs (ESP) [START_REF] Burton | Effectiveness of ABS and vehicle stability control systems[END_REF][START_REF] Murgovski | Predictive cruise control with autonomous overtaking[END_REF]. These tools act on vehicle inputs to correct the vehicle trajectory or to prevent the wheels from locking up in dangerous situations. These dangerous situations often appear when the car is exposed to severe external conditions such as slippery roads or sudden avoidance manoeuvres at high speed. With the rise of autonomous vehicles and all the requirements necessary to create a vehicle that should drive itself without the intervention of a human driver, the driver assistance systems, commonly called Advanced Driver Assistance Systems (ADAS) [START_REF] Acosta | Road friction virtual sensing: A review of estimation techniques with emphasis on low excitation approaches[END_REF], are particularly studied nowadays [START_REF] Piao | Advanced driver assistance systems from autonomous to cooperative approach[END_REF][START_REF] Nilsson | Strategic decision making for automated driving on two-lane, one way roads using model predictive control[END_REF]. In the case of autonomous vehicles, the ADAS should not only have corrective actions but also preventive ones in order to ensure the traffic fluidity. One of the path considered in the literature to develop these future ADAS goes through the acquisition of essential information influencing the vehicle behavior [START_REF] Martinez | A safe longitudinal control for adaptive cruise control and stop-and-go scenarios[END_REF]. Among the information which could improve the ADAS, the maximum friction coefficient [START_REF] Lee | Real-time slip-based estimation of maximum tire-road friction coefficient[END_REF] (also called the grip potential, friction potential or maximum tire-road friction coefficient [START_REF] Acosta | Virtual tyre force sensors: An overview of tyre model-based and tyre model-less state estimation techniques[END_REF]) is paramount. The maximum friction coefficient represents the maxi-mum effort a tire can transmit to the ground [START_REF] Rajamani | Tire-road friction-coefficient estimation[END_REF] before it begins to slip on the road. This grip potential mainly depends on three different components, (i) the tire (model, dimensions, pressure, wear, to name a few), (ii) the ground (e.g., type of covering, state of the road) and (iii) the presence of a third element at the interface between the tire and the road such as snow, water or dust [START_REF]The tyre grip[END_REF]. Furthermore, knowledge of the grip potential value turns out to be advantageous in various situations such as during emergency manoeuvres because it allows the ADAS to brake with the optimal brake pressure and thus, reduce the stopping distance or anticipate the right trajectory. In addition to these cases, the tire-road friction coefficient can also be used to detect the low grip area and thus reveals problems such as worn roads, poor rainwater drainage or need for snow removal [START_REF] Singh | Estimation of tire-road friction coefficient and its application in chassis control systems[END_REF].

As illustrated with this list of (non exhaustive) applications involving grip potential values, it is essential for car and tire manufacturers to develop reliable and accurate grip potential estimators. While several solutions are now available in the literature as described in Section 1.2, there is still room for improvement as explained in Section 1.2.2 where the main limits of these solutions are listed. Hence, the main objective of this thesis is to develop new friction potential estimators which bypass the difficulties encountered with the methods introduced in the literature so far.

State of the art

Several studies have already been carried out to estimate the maximum friction coefficient. Most of the developed methods are gathered in different review articles such as [START_REF] Acosta | Road friction virtual sensing: A review of estimation techniques with emphasis on low excitation approaches[END_REF][START_REF] Acosta | Virtual tyre force sensors: An overview of tyre model-based and tyre model-less state estimation techniques[END_REF][START_REF] Singh | Literature review and fundamental approaches for vehicle and tire state estimation[END_REF][START_REF] Rajamani | Tire-road friction-coefficient estimation[END_REF][START_REF] Khaleghian | A technical survey on tire-road friction estimation[END_REF][START_REF] Müller | Estimation of the maximum tire-road friction coefficient[END_REF][START_REF] Guo | A review of estimation for vehicle tire-road interactions toward automated driving[END_REF]. In these articles, different method classifications are proposed. However, the type of classification and the names employed vary from one article to another, which may confuse the reader. For this reason, a classification of the different method is suggested hereafter.

In this chapter, the state of the art focuses specifically on the different friction potential estimation methods. However, as can be seen in the following sections, solving the grip potential estimation problem often requires to solve intermediate problems such as tire forces estimation or modeling a suitable vehicle model. For this reason, more detailed states of the art are presented in the different chapters solving intermediate problems.

Friction potential estimation method classification

The suggested classification is highly inspired by the ones indicated in [START_REF] Acosta | Road friction virtual sensing: A review of estimation techniques with emphasis on low excitation approaches[END_REF] and [START_REF] Müller | Estimation of the maximum tire-road friction coefficient[END_REF]. This one is illustrated in Figure 1.1. The grip potential estimation methods can be separated into two main categories. The first class of solutions is the cause-based method, also called experimentbased method in [START_REF] Khaleghian | A technical survey on tire-road friction estimation[END_REF]. These techniques infer a grip potential value by measuring quantities affecting the grip potential such as the surface roughness or the water film [START_REF] Acosta | Road friction virtual sensing: A review of estimation techniques with emphasis on low excitation approaches[END_REF]. Once these quantities are measured, they are passed into a friction model in order to infer a grip potential value. This class of solutions has many benefits. The first one concerns the high level of accuracy reachable. For example, in [START_REF] Eichhorn | Prediction and monitoring of tyre/road friction[END_REF], a method based on an optical sensor detecting the surface roughness succeeds in estimating the grip potential with an uncertainty of 0.1 in 92% of the treated cases. Another significant advantage is the excitation level necessary. This one varies with the type of sensor used. In the case of optical sensor, since the quantity influencing the grip potential is the surface roughness, the required friction excitation level is close to zero.

Therefore this class of methods can be used under standard driving conditions [START_REF] Müller | Estimation of the maximum tire-road friction coefficient[END_REF]. However, these methods have also drawbacks. Most of them require additional sensors to detect the quantities affecting the grip potential. These sensors can be costly and thus inconsistent in an industrial context. Besides, as mentioned previously, once the sensor has measured the quantity of interest, this one should be passed into a model. Elaborating a physical model linking the measured quantity with the grip potential may be troublesome. One solution is to resort to a machine learning model (e.g., a neural networks [START_REF] Graber | Application of neural networks to external parameter estimation for nonlinear vehicle models[END_REF]). However, this solution requires a large database and can lead to a poor accuracy of the friction potential estimate in situations not supported by the training database. For these reasons, this class of solutions is difficult to apply in an industrial context where the use of costly sensor is prohibitive.

The second main category gathers the effect-based methods, also called model-based methods in [START_REF] Khaleghian | A technical survey on tire-road friction estimation[END_REF]. This class of solutions determines a friction potential by examining the impact of a grip potential variation on some quantities. For example, in icy conditions where the grip potential reduces drastically, the tire can no longer transmit all the engine torque to the ground and thus the longitudinal tire force F x decreases. Therefore one method to estimate the grip potential is to monitor the amount of the longitudinal tire force. Among the effect-based techniques, the authors distinguish three sub-classes (see Figure 1.1).

• The tire tread approach estimates the friction potential by measuring the deformation of the tire tread in the contact patch. This type of solution has the advantage to predict the friction potential before reaching it. However, it requires an additional sensor to measure the tire tread deformation.

• The vibration-based methods aim at inferring a grip potential estimate by examining the noise frequency content at the tire-road interface (high frequency method) or by establishing correlation between the grip potential and the resonance frequency of a vehicle subsystem (low frequency method).

• Tire slip-based methods estimate the grip potential by using the link between the normalized tire forces and the slips. In the longitudinal case, this link is characterized by the friction curve de-picted in Figure 1.2. Under standard driving conditions, the available friction information is represented by the points located at the bottom of the friction curve. Here, these points are called the friction points. Therefore, the main challenge of slip-based methods is to infer a grip potential estimate by using the friction points. Among the tire slip-based methods, we can distinguish between the black-box methods [START_REF] Rogers | A first course in machine learning[END_REF] and the grey-box methods which include a physical model representing the tire behavior. For example, in [START_REF] Berntorp | Bayesian tire-friction learning by Gaussian-process state-space models[END_REF], the author employs Gaussian processes to learn the friction curve from production car sensors. Thus, this black-box technique enables to determine not only the grip potential but also the entire friction curve as a function of the slip ratio. However, this method requires, like many other methods, high-friction data points to give reliable results.

Finally, as shown in Figure 1.1, an additional distinction between the methods is made depending on the vehicle dynamics. Some methods are exclusively designed for the longitudinal dynamics, i.e., when the vehicle drives in a straight line. On the contrary, lateral dynamics-based methods are effective when these dynamics occur, e.g., when the vehicle turns a corner or when it is in a roundabout. The coupled dynamicsbased methods work in both longitudinal and lateral situations. As illustrated in Figure 1.1, the method considered in this manuscript belongs to the grey-box tire slip-based methods. For this reason, the different methods presented subsequently mainly belong to the this class of methods as well. The reasons why we select the tire-slip based category are exposed in Section 1.3.

Article review

The different review articles [START_REF] Acosta | Road friction virtual sensing: A review of estimation techniques with emphasis on low excitation approaches[END_REF][START_REF] Acosta | Virtual tyre force sensors: An overview of tyre model-based and tyre model-less state estimation techniques[END_REF][START_REF] Singh | Literature review and fundamental approaches for vehicle and tire state estimation[END_REF][START_REF] Rajamani | Tire-road friction-coefficient estimation[END_REF][START_REF] Khaleghian | A technical survey on tire-road friction estimation[END_REF][START_REF] Müller | Estimation of the maximum tire-road friction coefficient[END_REF][START_REF] Guo | A review of estimation for vehicle tire-road interactions toward automated driving[END_REF] provide an important basis of all the works already carried out to solve the friction potential estimation problem. This section summarizes the main results of these different articles.

According to the different review articles, one main constraint preventing the friction potential estimation is the excitation level. Under standard driving conditions, the engine torque does not reach high values. Consequently, the excitation level might be insufficient to apply most of the grip potential estimation methods [START_REF] Acosta | Road friction virtual sensing: A review of estimation techniques with emphasis on low excitation approaches[END_REF][START_REF] Singh | Literature review and fundamental approaches for vehicle and tire state estimation[END_REF][START_REF] Rajamani | Tire-road friction-coefficient estimation[END_REF][START_REF] Müller | Estimation of the maximum tire-road friction coefficient[END_REF].

One way to overcome this constraint is to resort to cause-based methods. As indicated in [START_REF] Müller | Estimation of the maximum tire-road friction coefficient[END_REF] and reminded in Section 1.2.1, one of the main advantages of this type of methods is the low excitation level necessary. However, most of the cause-based methods include data-driven models (usually neural networks) and thus, require training sets. Besides, these methods provide less accurate results when situations not handled by the training sets occur [START_REF] Khaleghian | A technical survey on tire-road friction estimation[END_REF][START_REF] Müller | Estimation of the maximum tire-road friction coefficient[END_REF]. Finally, the efficiency of these methods needs to be proven on production vehicle where few sensors are available [START_REF] Singh | Literature review and fundamental approaches for vehicle and tire state estimation[END_REF]. For these reasons, a large number of elaborated methods belongs to the effect-based methods. As underlined in [START_REF] Acosta | Road friction virtual sensing: A review of estimation techniques with emphasis on low excitation approaches[END_REF][START_REF] Khaleghian | A technical survey on tire-road friction estimation[END_REF][START_REF] Müller | Estimation of the maximum tire-road friction coefficient[END_REF][START_REF] Guo | A review of estimation for vehicle tire-road interactions toward automated driving[END_REF], among this class of solutions, tire slip-based methods are predominant.

According to the different authors, except for the slip-slope method (see Section 1.2.3), most of the slip-based methods require a signifi-cant excitation level, i.e., friction points close to the grip potential, in order to achieve friction potential inference. The slip-slope method is a slip-based method introduced by F. Gustafsson in [START_REF] Gustafsson | Slip-based tire-road friction estimation[END_REF] which assumes the existence of a link between the slope of the friction curve and the maximum. This method is described in more details in Section 1.2.3. However, it is important to underline that the idea of a link between the slope and the friction potential is not unanimously accepted [START_REF] Acosta | Road friction virtual sensing: A review of estimation techniques with emphasis on low excitation approaches[END_REF]. Besides, this hypothetical link holds while other quantities such as the tire inflation pressure and tire temperature are fixed. Therefore, in practice, the need to keep constant some quantities that are difficult to control causes problems of test repeatability [START_REF] Khaleghian | A technical survey on tire-road friction estimation[END_REF].

On another note, one method called self-aligning torque (SAT) method might provide an accurate grip potential estimation by using friction points around µ = 0.3 [START_REF] Acosta | Road friction virtual sensing: A review of estimation techniques with emphasis on low excitation approaches[END_REF][START_REF] Singh | Literature review and fundamental approaches for vehicle and tire state estimation[END_REF]. This method consists in predicting the grip potential by looking at the curve of the self-aligning torque as a function of the slip angle. An example of SAT curve is illustrated in begins to decrease. Physically, a decrease of the SAT means that the tire is close to its friction limit. Besides, as depicted in Figure 1.3, the SAT begins to decrease before the grip potential is reached. Therefore, the main advantage of the SAT is that it can bring information on the friction grip potential with a reduced excitation level. However, this method holds for lateral dynamics only. Besides, it requires to model the vehicle suspensions and steering kinematics which can be troublesome. These constraints lead to the necessity to introduce another class of estimation methods, especially in the longitudinal case.

According to Acosta et al. [START_REF] Acosta | Road friction virtual sensing: A review of estimation techniques with emphasis on low excitation approaches[END_REF], one promising alternative to the slipbased methods is the vibration-based method and especially the high vibration-based method. As reminded in Section 1.2.1, this class of solutions examines the noise frequency content at the tire-road interface and then, deduces a friction potential. In order to apply this type of methods, it is necessary to find how tire-road interactions physically influence the grip potential. For example, it is accepted that the micro and macro roughnesses impact the grip potential, however, the literature lacks of theoretical models providing an explicit link between these quantities. Therefore, this class of methods is mostly applied with datadriven techniques such as neural networks. Consequently, similarly to the cause-based approach, these techniques require a high amount of data and may deliver poorly results in situations outside of the training set.

Finally, according to the different review articles, the slip-based techniques seem to be the only techniques able to provide a friction potential estimate without additional sensors. For this reason, as explained further in Section 1.3, the approach suggested in this manuscript belongs to this class of solutions. In particular, the approach designed in this manuscript is mainly inspired of F. Gustafsson's [START_REF] Gustafsson | Slip-based tire-road friction estimation[END_REF] and L. Ray's works [START_REF] Ray | Nonlinear tire force estimation and road friction identification: Simulation and experiments[END_REF]. For this reason, in the following, these two methods are described in further details.

Slip-slope method

As mentioned in review articles, the slip-slope method is frequently highlighted when the problem of grip potential estimation raises. The slip-slope method has been introduced by F. Gustafsson in [START_REF] Gustafsson | Slip-based tire-road friction estimation[END_REF]. This method is based on the assumption that the slip-slope K x , defined by the initial slope of the friction curve, depends on the grip potential. As a consequence, under this assumption, if the slip-slope is estimated, which is possible with measurements available under standard condi-tions, the grip potential can be determined. This concept is illustrated in Figure 1. [START_REF] Acosta | Virtual tyre force sensors: An overview of tyre model-based and tyre model-less state estimation techniques[END_REF] where the red lines are the slip-slope of the different friction curves. In his paper, F. Gustafsson estimates a front and a rear slip-slope. However, in order to simplify the approach idea, the following explanations take into account of a global slip-slope K x only. A scheme summarizing the slip-slope method is illustrated in Figure 1.5. The assumed where δ is a wear-related offset. This relation can be straightforwardly rewritten as follows

s = 1 K x µ + δ. (1.2)
In order to obtain accurate estimates of δ and K x (or more precisely its inverse 1 Kx ), F. Gustafsson implements a Kalman filter such as the one indicated in Figure 1.5. In addition, in order to improve the time response of the state observer, F. Gustafsson adds a CUSUM detection algorithm [START_REF] Basseville | Detection of abrupt changes: theory and application[END_REF]. This algorithm should detect strong variations and correct the Kalman filter coefficients accordingly. Once K x has been estimated, its value is used with a classifier in order to determine a type of road and thus, a grip potential estimate.

One of the main advantages of the slip-slope method is that it does not require any additional sensors. Indeed, as represented in Figure 1.5, the method only needs measurements provided by the vehicle CANBUS. These ones are used to compute friction points {s, µ}. However, the slipslope method has been particularly discussed and the literature contains contradictory conclusions [START_REF] Acosta | Road friction virtual sensing: A review of estimation techniques with emphasis on low excitation approaches[END_REF]. For some authors, the method is applicable with small friction measurements (lower than 0.3) [START_REF] Singh | Literature review and fundamental approaches for vehicle and tire state estimation[END_REF]. For others, it requires higher friction (around 0.4 at least) and additional investigations to conclude on the validity of the method [START_REF] Müller | Estimation of the maximum tire-road friction coefficient[END_REF][START_REF] Acosta | Road friction virtual sensing: A review of estimation techniques with emphasis on low excitation approaches[END_REF]. In addition, as pointed out by F. Gustafsson in [START_REF] Gustafsson | Slip-based tire-road friction estimation[END_REF], the slip-slope depends on several factors such as the tire inflation pressure, the tire temperature or the tire wear, to name a few. Therefore, the method is effective only if all these factors are known. Thus, this method necessitates a calibration to be effective. This point is however not the main issue of this method. The main problem of the slip-slope method is the fundamental assumption of a link between the slip-slope and the grip potential [START_REF] Acosta | Road friction virtual sensing: A review of estimation techniques with emphasis on low excitation approaches[END_REF]. As mentioned in [START_REF] Gustafsson | Slip-based tire-road friction estimation[END_REF], this link has not been proven theoretically yet. Tests to establish the existence of this link are performed in [START_REF] Andrieux | New results on the relation between tyre-road longitudinal stiffness and maximum available grip for motor car[END_REF] but the authors conclude that a link between the slope and the tire-road friction coefficient cannot be confirmed. Hence the slip-slope method needs more examinations to verify if it can provide a reliable estimate of the grip potential in all situations.

Tire forces observer and Bayesian selection

Another method allowing to infer the grip potential is introduced by L. Ray in [START_REF] Ray | Nonlinear tire force estimation and road friction identification: Simulation and experiments[END_REF] and is summarized by Rajamani et al. in [START_REF] Rajamani | Tire-road friction-coefficient estimation[END_REF]. Figure 1.6 illustrates the concept of the method. As shown in Figure 1.6, Figure 1.6: L. Ray's approach concept adapted from [START_REF] Ray | Nonlinear tire force estimation and road friction identification: Simulation and experiments[END_REF].

the fundamental idea of L. Ray's method is to infer a grip potential value by statistically comparing a friction estimate provided by an Extended Kalman filter (EKF) with friction curves obtained with a tire model. In order to simplify the explanations, we restrict the basic idea to the longitudinal case. However, note that L. Ray's work was applied with both longitudinal and lateral dynamics.

First of all, L. Ray implements an EKF with a vehicle model representing both longitudinal and lateral dynamics as well as roll dynamics. The observer used by L. Ray provides at each time k state and tire forces estimates.

Remark 1.1. One important contributions in L. Ray's work consists in considering the longitudinal and lateral tire forces as random walks and thus, simplifying their integration in the state vector.

Once the state and tire forces are estimated, these ones are used to compute, at each time k, an estimated friction point {ŝ EKF , μEKF } represented in red in Figure 1.7. Besides, by using the tire model described in [START_REF] Szostak | Analtical modeling of driver response in crash avoidance maneuvering. volume 2: An interactive tire model for driver/vehicle simulation[END_REF], L. Ray computes off-line different friction curves. These curves are plotted in blue in Figure 1.7. In order to estimate a grip potential value, L. Ray compares the friction curves and the EKF estimates in a Bayesian way (see the bloc "Bayesian hypothesis selection" in Figure 1.6). More precisely, the basic idea is to translate the distance between the EKF estimate (the red point) and the different friction curves (the green points) into a probability that the EKF estimate belongs to a given friction curve. Hence, once we know at which friction curve the EKF estimate belongs to, we can infer a grip potential value straightforwardly by examining the maximum of this friction curve. For the sake of clarity, this idea is explained in a mathematical form. Let p(µ max T j |μ EKF k ) be the probability that the EKF friction point estimated at time k {ŝ EKF k , μEKF k } belongs to the friction curve T j which has the maximum µ max T j . According to the Bayes formula [START_REF] Theodoridis | Machine learning: a Bayesian and optimization perspective[END_REF], this probability is given by

p(µ max T j |μ EKF k ) = p(μ EKF k |µ max T j )p(µ max T j |μ EKF k-1 ) J i=1 p(μ EKF k |µ max T i )p(µ max T i |μ EKF k-1 ) . (1.3)
Therefore, it is necessary to compute a probability p(μ EKF |µ max T ). In order to compute this probability, L. Ray uses the Gaussian formulation

p(μ EKF |µ max T ) = 1 (2π) n/2 S 1/2 exp - 1 2 (μ EKF -µ T ) ⊤ S -1 (μ EKF -µ T ) , (1.4)
where µ T is the ordinate of the friction point belonging to the friction curve T , i.e., the ordinate of the green point of the friction curve T . Once p(µ max T j |μ EKF k ) is computed for all the friction curves, the grip potential is estimated by taking the weighted average value

μmax k = J j=1 p(µ max T j |μ EKF k )µ max T j . (1.5)
L. Ray's method has several advantages. First of all, it accurately reconstructs the states and especially the tire forces in various situations. Besides, it does not require additional sensors other than the ones already fitted on production vehicle. Finally, it provides accurate estimates of the grip potential as long as the excitation is high enough. However, for a low excitation level, the algorithm encounters difficulties in estimating a reliable grip potential value due to the abacus friction curves which are too close between them for low slip ratio (see Figure 1.7). In addition, the tests made to validate the method have been conducted on dry asphalt only. Thus, tests on others surfaces are needed in order to validate the method performance.

Current methodology

According to the different information provided by the state of the art, two constraints should be faced if we want to solve the problem of estimating the grip potential. Firstly, because the grip potential is an information which is intended to be used in an industrial context, it is necessary to develop a cost effective method which does not require any additional sensor. Therefore, the estimation should work by using only the measurements provided by the sensors fitted on production vehicles. Secondly, the grip potential should be estimated under standard driving conditions. This area corresponds to situations where the friction data points are lower than 0.3. These two main constraints and the review of the different methods exposed in the literature should guide us all along the elaboration of the estimation approach.

First of all, it is necessary to choose the dynamics taken into account. According to the state of the art [START_REF] Acosta | Road friction virtual sensing: A review of estimation techniques with emphasis on low excitation approaches[END_REF][START_REF] Acosta | Virtual tyre force sensors: An overview of tyre model-based and tyre model-less state estimation techniques[END_REF][START_REF] Singh | Literature review and fundamental approaches for vehicle and tire state estimation[END_REF], in the lateral case, the SAT method seems to provide accurate estimate of the slip potential with a low excitation level. However, in the purely longitudinal case, a curve such as the SAT curve illustrated in Figure 1.3 does not exist. Therefore, the required excitation level is much higher and does not correspond to standard driving conditions. An exception is the slip-slope method but, as mentioned previously, this method needs further test in order to conclude on its validity. For all these reasons, the present study will focus only on the longitudinal dynamics since important accuracy gain are expected.

Due to both constraints, i.e., the constraint on the available sensor measurements and the friction points limit, solving the global problem of the grip potential estimation can be split into the following two subproblems:

• how to estimate a friction potential from the knowledge of friction data points?

• How to infer friction data points by using only sensors fitted on production vehicle?

These two problems are represented in the master plan (see Figure 1.8) elaborated in order to estimate the grip potential.

The first problem handled in this manuscript concerns the estimation of a friction potential value from friction data points. In this part, we make the assumption that we have access to noisy friction point measurements. This assumption is not strong because even if the friction data points are reconstructed and not directly measured, they are still inferred from measurements coming from real sensors. Thus, they are still corrupted by noise. Besides, because standard driving conditions must be accounted for, the friction measurements are almost always lower than 0.3. These constraints on the measurements guide our choice on the use of grey-box methods. Indeed, although black-box methods are particularly efficient for interpolation, they might be less suitable for extrapolation, and particularly to estimate points far from the available measurements. The method suggested in this manuscript is a friction curve fitting method based on a combination of the Maximum Likelihood (ML) method and an adaptive Monte-Carlo Markov Chain (MCMC) method [START_REF] Luengo | A survey of Monte Carlo methods for parameter estimation[END_REF]. The MCMC method is suggested because it satisfies well our constraints on the measurements. More precisely, this class of solutions has the advantage to work well with few noisy points by giving the opportunity to the user to take into account of the noise properties [START_REF] Rogers | A first course in machine learning[END_REF]. This method requires the use of a tire model which should represent the friction curve accurately. One tire model commonly used by tire manufacturers is the Pacejka's magic formula tire model [START_REF] Pacejka | Tire and vehicle dynamics[END_REF]. This one is described in more detail in Chapter 2. Hence, the idea suggested to solve the first sub-problem consists in estimating the unknown parameters of the Pacejka's model by initializing an MCMC method with ML estimates. The part of the master plan dealing with the first problem is the one framed in red in Figure 1.9 with the bloc "Friction potential estimation" detailed in Figure 1.10.

Remark 1.2.

At this stage, it can be interesting to compare our approach with the slip-slope method. First of all, similarly to the slip-slope method, our method resort to friction points located at the bottom of the friction curve. However, unlike F. Gustafsson's approach, in our method, the friction potential is not inferred from an estimate of the slip-slope K x but from the estimates of a tire model parameters describing the entire friction curve. In other words, rather than placing our confidence in the slip-slope, we transfer our confidence in a tire model commonly used by tire manufacturers [START_REF] Pacejka | Tire and vehicle dynamics[END_REF]. Once the first sub-problem has been tackled, it remains to solve the second one, i.e., how to infer friction data points by using only sensors fitted on production vehicle. This problem corresponds to the part framed in red in Figure 1.11. One approach providing interesting results is L. Ray's method [START_REF] Ray | Nonlinear tire force estimation and road friction identification: Simulation and experiments[END_REF]. In her paper, L. Ray uses a double-track model with an EKF in order to determine state and tire forces estimates. Then, these ones are used to compute slip ratio estimates and normal load estimates.

The approach suggested hereafter is highly inspired of L. Ray's method. In our case, the idea is to also use an EKF but to modify the considered as well as roll dynamics. Although this vehicle model is accurate, it involves many quantities and thus, necessitates a large state vector which make difficult the setting of the noise covariance matrices. Besides, in our case, we focus our study on the longitudinal behavior. For these reasons, the considered model is a single track model. Furthermore, the friction data points are computed from estimates of the slip ratio s, the longitudinal tire force F x and the normal load F z . Hence, these quantities should be estimated with accuracy. In particular, during a travel, the load applied on the front and the rear part of the vehicle evolves and thus, modifies the normal load F z . This phenomena called load transfer is not taken into account by L. Ray in her study. However, it should be taken into account if we want to accurately estimate the normal load. Therefore, the vehicle model used here is extended with a suspension model.

Finally, another difference with L. Ray's work is the setting of the EKF noise covariance matrices. When we use a Kalman filter, it is necessary to tune noise covariance matrices. The common way employed to set these matrices is to proceed to an error-trial phase based on our confidence on the measurements and the model accuracy. However, in practice, if the state vector is high dimensional, finding the different matrices coefficients can be quite troublesome. In order to simplify the setting, different methods were established. Most of them are described in the survey [START_REF] Duník | Noise covariance matrices in state-space models: A survey and comparison of estimation methods-part i[END_REF]. Although some of these methods simplify the tuning of the noise covariance matrices, they still require the tuning of other parameters. Therefore, a method which do not require any setting is introduced in this manuscript. This method is based on subspace identification.

Structure of the thesis and road map

The thesis manuscript in organized as follows (see Figure 1.12). Chapter 2 is dedicated to the first sub-problem, namely, how to estimate the grip potential by having access to friction data points. In this chapter, an MCMC procedure is combined with a Maximum likelihood approach in order to get a grip potential estimate. This approach is applied to both simulated and real friction data points. Chapter 3 describes the theoretical aspect of the method used herein in order to estimate friction data points from classical measurements provided by sensors fitted on production vehicle. In particular, the employed method is based on an EKF. Consequently, it appears the question of the settings of the noise covariance matrices. Chapter 4 presents a method designed in order to set the noise covariance matrices. This method based on subspace identification does not requires any tuning from the user. Chapter 5 is devoted to the application of the methods exposed in Chapter 3 and 4 with simulated data coming, on the one hand, from a simplified vehicle model and, on the other hand, from a realistic vehicle simulator. Finally Chapter 6 concludes this manuscript.

As reported in Figure 1.12, given that the problematic of grip potential estimation is divided into two sub-problems, the reader can read Chapter 3-5 without having to read Chapter 2.

CHAPTER 2

Grip potential estimation

Introduction

One essential part of the work done during this thesis is the development and the application of a procedure allowing to estimate the friction potential from friction points. As mentioned in Chapter 1, the friction points cannot be measured directly with sensor fitted on production vehicle. For this reason, they have to be estimated and it is the main subject of Chapter 3. However, in Chapter 2, the friction points are assumed to be available. Hence, the objective of this current chapter is to determine a grip potential estimate from friction point measurements (see Figure 2.1).

In this chapter, the grip potential estimation procedure should operate in realistic situations, i.e., with friction points corresponding to standard driving conditions. As a consequence, the used friction data points are corrupted by noise and are almost always lower than 0.3 (see Figure 2.2). These constraints on the measurements can make difficult the use of black-box methods which are particularly efficient for interpolation, but might be less appropriate to estimate points far from the measurement area as it is the case here. Hence, a parametric model-based approach is considered hereafter to estimate the tire friction potential. The method suggested in this chapter is a tire curve fitting method based on a Maximum Likelihood (ML) method and an adaptive Monte-Carlo Markov Chain (MCMC) method [START_REF] Luengo | A survey of Monte Carlo methods for parameter estimation[END_REF], respectively. The MCMC method is suggested because it has the advantage to work well with few points and enables the user to take into consideration the noise properties effectively [START_REF] Rogers | A first course in machine learning[END_REF]. This method requires the use of a tire model which should represent the tire behavior in various situations accurately. One tire model fulfilling these conditions is the Pacejka's magic formula tire model [START_REF] Pacejka | Tire and vehicle dynamics[END_REF]. This model is commonly used by tire manufacturers because it is a parsimonious model which represents accurately the tire behavior. Therefore, the approach suggested in this chapter consists in estimating the unknown parameters of the Pacejka's model by initializing an MCMC method with ML estimates. Figure 2.3 illustrates the grip potential estimation procedure by using friction point measurements.

The rest of the chapter is organized as follows. Section 2 is dedicated to the description of the problem and the notations used in this chapter. Section 3 describes the method used herein to estimate the tireroad friction coefficient and underlines the different settings necessary to maximize the method performance. Section 4 presents the results obtained when the method is tested with simulated data and real data, respectively. Finally, Section 5 concludes this chapter.

Problem formulation and notations

In this chapter, we aim at determining the maximum of the tire-road friction from noisy and short data sets acquired with standard sensor signals. This maximum value, called in the literature the tire friction potential or the grip potential [START_REF] Acosta | Virtual tyre force sensors: An overview of tyre model-based and tyre model-less state estimation techniques[END_REF], is the maximum force a tire can generate before sliding. More precisely, when only longitudinal dynamics is considered 1 , the tire friction µ is defined as follows [START_REF] Canudas-De-Wit | Dynamic friction models for road/tire longitudinal interaction[END_REF] and thus, the grip potential is

µ = F x F z , (2.1)
µ max = max F x F z , (2.2)
where F x and F z are the longitudinal tire force and the tire normal load, respectively.

As shown, e.g., in [START_REF] Devinney | Factors affecting tire traction[END_REF], the friction µ, and by extension the grip potential µ max , is impacted by several quantities such as the tire inflation pressure, the tire temperature, the tire load or the tire tread depth to name a few. Among all these quantities, the slip ratio s defined by (see Table 2.1 for the definition of ω, R rol and v x , respectively)

s = ωR rol -v x max(ωR rol , v x ) , (2.3)
plays a central role [START_REF] Rajamani | Vehicle Dynamics and Control[END_REF]. The basic idea of the solution developed in this chapter is to determine the friction dependency in the slip ratio s with a (parametric) model

µ = f (s), (2.4)
and to compute its maximum straightforwardly. Figure 2.4 illustrates standard but significant characteristics of this link between the tire friction µ and the slip ratio s. Such curves clearly show that the mapping µ = f (s) is a nonlinear function having a sharp initial slop, a maximum, then a gentle decrease with larger values of s. Besides, the different friction curves indicate that µ and thus µ max depend of the road characteristics as well. For obvious practical reasons, it is essential that the model learning strategy introduced hereafter gives access to µ max accurately during "standard" driving conditions, i.e., when only low values of µ (lower than 0.3) and s are observed. Said differently, the solutions developed in this chapter to determine the grip potential µ max accurately must estimate this unknown mapping µ = f (s) from

• data sets acquired for low values of s, i.e., far from the ones corresponding to the tire friction potential,

• real measurements acquired with standard sensors, i.e., from noisy data sets.

In order to reach this goal, a data driven curve fitting, also named model learning solution, is suggested hereafter.

In addition to the aforementioned noisy data, the second ingredient of our model learning strategy is the choice of the mapping f (•). Different model structures have been suggested in the literature [START_REF] Pacejka | Tire and vehicle dynamics[END_REF][START_REF] Dugoff | Tire performance characteristics affecting vehicle response to steering and braking control inputs[END_REF][START_REF] Canudas-De-Wit | A new model for control of systems with friction[END_REF][START_REF] Van Zanten | Measurement and simulation of transient tire forces[END_REF][START_REF] Szostak | Analtical modeling of driver response in crash avoidance maneuvering. volume 2: An interactive tire model for driver/vehicle simulation[END_REF]. In this contribution, a specific attention is paid to a parametric modelling approach because it can take into account the a priori knowledge about the tire behavior. More specifically, among the different tire models available in the literature [START_REF] Pacejka | Tire and vehicle dynamics[END_REF][START_REF] Dugoff | Tire performance characteristics affecting vehicle response to steering and braking control inputs[END_REF][START_REF] Canudas-De-Wit | A new model for control of systems with friction[END_REF][START_REF] Van Zanten | Measurement and simulation of transient tire forces[END_REF][START_REF] Szostak | Analtical modeling of driver response in crash avoidance maneuvering. volume 2: An interactive tire model for driver/vehicle simulation[END_REF], the very popular Pacejka's model [START_REF] Pacejka | Tire and vehicle dynamics[END_REF] is selected because (i) it can represent the tire dynamics under various practical conditions accurately (ii) it is a parsimonious model, i.e. it has a limited number of parameters, which facilitates its use for parametric estimation. Such a model satisfies the following equation (see Figure 2.5 for a geometrical interpretation of the so-called "magic formula" parameters [START_REF] Pacejka | Tire and vehicle dynamics[END_REF])

µ = D sin (C arctan [Bx -E(Bx -arctan(Bx))]) + s v , (2.5a)
x = s + s h , (2.5b)
or, in a compact way,

µ = f (s, θ) with θ = B C D E s h s v ⊤ , (2.6)
where the corresponding "magic formula" parameters composing the unknown parameter vector θ ∈ R 6×1 are defined in Table 2.2. • lower and upper bounds for each component of θ, our main goal is to estimate the parameter vector θ accurately. As illustrated in Section 2.4, such a goal can be reached by resorting to the data driven model learning solutions introduced in the next section.

Pacejka's model learning: combining maximum likelihood and Monte Carlo Markov Chain

When parameter estimation of nonlinear static function from noisy data comes into play, the dominant solution in the literature consists in resorting to the maximum likelihood (ML) framework [START_REF] Seber | Nonlinear Regression[END_REF]. The ML solutions indeed benefit from interesting statistical properties like their asymptotic unbiasedness, efficiency and consistency under mild conditions [START_REF] Theodoridis | Machine learning: a Bayesian and optimization perspective[END_REF]Chapter 3]. Such asymptotic properties can indeed help the user deliver reliable parameter estimates with quantification of the estimation accuracy [START_REF] Van Den Bos | Parameter estimation for scientists and engineers[END_REF]. As shown, e.g., in [START_REF] Garatti | Assessing the model quality in system identification: the asymptotic theory revisited and application to iterative control[END_REF], using asymptotic results to quantify model uncertainties from finite data sets can yield unreliable results when small size sets of data are used. This is the main reason why the Bayesian framework is nowadays often suggested to tackle parameter estimation problems when short data sets are handled [START_REF] Luengo | A survey of Monte Carlo methods for parameter estimation[END_REF]. The techniques involved in this framework, and more specifically the Monte Carlo (MC) random sampling methods [START_REF] Robert | Introducing Monte Carlo methods with R[END_REF], have the attractive property to give access to samples of a desired distribution (instead of a parameter vector as the ML solutions do) from user defined density function proposals by resorting to easy-to-implement algorithms mainly based on accept or reject conditions [START_REF] Spall | Estimation via Markov Chain Monte Carlo[END_REF]. Interesting when short data sets are considered, as adviced in [6, Section 1], "it is important that we do not treat them as black boxes in order to obtain the best results out of this class of algorithms". This practical observation is probably the main reason why we suggest hereafter combining ML and MC algorithms to benefit from their respective advantages, i.e., we tackle the problem of data driven model learning by using first a ML approach for get reliable initial guesses to tune, in a second step, a Monte Carlo Markov Chain (MCMC) algorithm. These two steps are more precisely introduced in the next sections.

Maximum likelihood: a short review

Because the tire friction is measured with (standard) sensors, the measured tire friction samples are assumed to be noisy. As a direct consequence, the model we consider from now on has the following form

µ i = f (s i , θ) + e i , i ∈ {1, • • • , N }, (2.7)
where e i , i ∈ {1, • • • , N }, stands for a sample of the realization e N 1 = (e i ) i∈{1,••• ,N } of a stochastic sequence (e t ) t∈Z characterized by a specific probability density function (pdf) p e (e) [START_REF] Leon-Garcia | Probability, Statistics, and Random Processes for Electrical Engineering[END_REF]. Such an additive noise in the model description implies that the friction µ should also been considered as a random variable. Said differently, the acquired samples µ i , i ∈ {1, • • • , N }, are now seen as components of the realization µ N 1 = (µ i ) i∈{1,••• ,N } of the stochastic sequence (µ t ) t∈Z . This observation is at the basis of the ML estimation solutions. More specifically, given the pdf of the random sequence (µ t ) t∈Z parameterized by θ, i.e., given

p µ (µ; s, θ), (2.8)
the maximum likelihood techniques aim at determining the unknown parameter vector θ which "moves" this pdf such that it is maximum for our observed data or realization (µ i ) i∈{1,••• ,N } . This probability density function p µ (µ; s, θ) evaluated at our observed data

µ N 1 = (µ i ) i∈{1,••• ,N } acquired for the slip ratios s N 1 = (s i ) i∈{1,••• ,N } , i.e., p µ (µ = µ N 1 ; s = s N 1 , θ)
, is called the likelihood function and is denoted herein as follows

ℓ(µ N 1 ; s N 1 , θ).
(2.9)

Thus, the maximum likelihood solutions estimate the parameter vector

θ ML = arg max θ ℓ(µ N 1 ; s N 1 , θ). (2.10)
The solution of this maximization problem is strongly linked to the analytic equation of the likelihood function [START_REF] Van Den Bos | Parameter estimation for scientists and engineers[END_REF]. Because of the model structure given in Eq. (2.7), the likelihood function ℓ(µ N 1 ; s N 1 , θ) is directly related to the noise probability density function p e (e). More specifically, if it is assumed that the noise at each data point is independently distributed and zero mean, we have

ℓ(µ N 1 ; s N 1 , θ) = i∈{1,••• ,N } P(f (s i , θ), σ 2 i ), (2.11)
where

σ 2 i , i ∈ {1, • • • , N }, stands for the variance of component e i , i ∈ {1, • • • , N }, while P(• 1 , • 2
) is a generic notation for the noise probability density function of mean • 1 and variance • 2 , respectively. Different popular continuous density functions can be suggested for P(• 1 , • 2 ) [START_REF] Leon-Garcia | Probability, Statistics, and Random Processes for Electrical Engineering[END_REF]. In many practical cases, considering noise components which are • conditionally independent,

• normally distributed with zero mean and with the same standard deviation σ, is a reliable assumption. This is the main reason why we consider herein that

ℓ(µ N 1 ; s N 1 , θ) = i∈{1,••• ,N } N (f (s i , θ), σ 2 ), (2.12)
where N (• 1 , • 2 ) stands for the normal probability density function of mean • 1 and variance • 2 , respectively. Because the natural logarithm is a monotonic increasing function, the maximum of ℓ(µ N 1 ; s N 1 , θ) occurs at the same value of θ as the maximum of the log-likelihood defined as follows

L(µ N 1 ; s N 1 , θ) = log(ℓ(µ N 1 ; s N 1 , θ)) = -N log(σ) - N 2 log(2π) - 1 2σ 2 V (θ), (2.13
) with

V (θ) = N i=1 (µ i -f (s i , θ)) 2 .
(

Such a direct relation between the (log-)likelihood and the least squares cost function V (θ) proves that maximizing ℓ(µ N 1 ; s N 1 , θ) with respect to θ boils down to minimizing V (θ) with respect to θ. The minimization of the cost function V (θ) requires the use of nonlinear optimization algorithms because of the nonlinear dependency of f with respect to θ. Several algorithms are available in the literature for minimizing V (θ) under such practical conditions [START_REF] Nocedal | Numerical optimization[END_REF][START_REF] Rao | Engineering Optimization: Theory and Practice[END_REF]. In this chapter, the interiorreflective Newton method described in [START_REF] Coleman | An interior trust region approach for nonlinear minimization subject to bounds[END_REF] (see also [START_REF] Nocedal | Numerical optimization[END_REF]Chapter 19]) is selected to reach a local optimum of V (θ) in order to take into the aforementioned lower and upper bounds on θ explicitly during this numerical optimization.

Because the ML estimator is generated from noisy data, the estimated parameter vector θ ML can be considered as a realization of a random vector θ ML . As any random vector, θ ML can be characterized by standard (centered) moments [START_REF] Leon-Garcia | Probability, Statistics, and Random Processes for Electrical Engineering[END_REF] from which uncertainty domain or estimation quality certificates can be generated [START_REF] Kay | Fundamentals of Statistical Signal Processing: Estimation Theory[END_REF]. As proved, e.g., in [START_REF] Van Den Bos | Parameter estimation for scientists and engineers[END_REF]Chapter 5], the asymptotic covariance matrix of the ML parameter estimator is

cov(θ ML ) = σ 2 (F ⊤ θ ML F θ ML ) -1 , (2.15)
where F θ ML is the Jacobian matrix defined as

F θ ML =     ∂f (s 1 ,θ) ∂θ 1 • • • ∂f (s 1 ,θ) ∂θn θ . . . . . . . . . ∂f (s N ,θ) ∂θ 1 • • • ∂f (s N ,θ) ∂θn θ     θ=θ ML ∈ R N ×n θ . (2.16)
Furthermore, an unbiased estimate of σ 2 is given by [95, Chapter 5]

σ2 = V (θ ML ) N -n θ .
(2.17)

These statistical results will be used in Section 2.4 for characterizing the estimation quality.

Adaptive Metropolis algorithms

While the minimization of V (θ) is feasible whatever the friction of data samples, the asymptotic statistical results of the ML estimator introduced previously may be not reliable when only low-friction data sets are handled (see Section 2.4.1 for an illustration of this claim). This is the main reason why the ML class of solutions should be combined with another class of methods when, as shown in Section 2.4, the friction data points remain under 0.3. The Bayesian framework, and more specifically the MC random sampling techniques, are used hereafter as a second step of our model learning solution because of their efficiency under such practical constraints [START_REF] Luengo | A survey of Monte Carlo methods for parameter estimation[END_REF]. By using again the fact that the unknown parameter vector θ can be treated as a realization of a random vector θ, the Bayesian inference solutions aim at determining the conditional or a posteriori probability density function p θ,µ (θ | µ). More precisely, these solutions determine samples which form a distribution that asymptotically approaches p θ,µ (θ | µ) [START_REF] Spall | Estimation via Markov Chain Monte Carlo[END_REF]. Indeed, as shown hereafter, once these samples with common pdf are generated accurately, empirical means or covariance matrices can be computed effectively (see also [START_REF] Rogers | A first course in machine learning[END_REF]). In order to reach this goal, the heart of the Bayesian approach is the following pdf based Bayes formula [START_REF] Theodoridis | Machine learning: a Bayesian and optimization perspective[END_REF]

p θ,µ (θ | µ) = p θ,µ (µ | θ)p θ (θ) p µ (µ) , (2.18)
given two jointly distributed random vectors µ and θ, where [START_REF] Kay | Fundamentals of Statistical Signal Processing: Estimation Theory[END_REF] • p θ (θ) stands for the prior pdf of θ, i.e., the belief we have on θ before seeing data,

• p θ,µ (µ | θ) stands for the conditional pdf of µ given θ, i.e., our knowledge provided by µ conditioned on knowing θ,

• p θ,µ (θ | µ) stands for the posterior pdf we are interested in,

• p µ (µ) stands for the marginal distribution of µ, i.e.,

p µ (µ) = • • • p θ,µ (µ | θ)p θ (θ)dθ. (2.19)
As far as our model learning problem is concerned, we assume the access to lower and upper bounds for each parameter. This prior is translated into a uniform flat prior p θ (θ) ∝ 1 within these bounds, i.e.,

θ i ∼ U([θ i lower , θ iupper ]), i ∈ {1, • • • , n θ }. (2.20)
In addition, by assuming again zero mean independent and normally distributed measurement errors with a (known) constant variance σ 2 , the conditional pdf p θ,µ (µ | θ) can be written as follows

p θ,µ (µ | θ) = 1 (2πσ 2 ) N/2 exp - V (θ) 2σ 2 , (2.21)
or, equivalently,

p θ,µ (µ | θ) ∝ exp - V (θ) 2σ 2 . (2.22)
Remark 2.1. As pointed out, e.g., in [START_REF] Kay | Fundamentals of Statistical Signal Processing: Estimation Theory[END_REF], this pdf is identical in form to the likelihood function introduced previously. In the Bayesian case, the pdf is however a conditional pdf while, in the ML approach, the pdf is an unconditional pdf parameterized by θ.

In a nutshell, for our model learning problem,

p θ,µ (θ | µ) ∝ exp -V (θ) 2σ 2 p µ . ( 2

.23)

The challenging step in the Bayesian approach is the determination of p µ . Indeed, in most of practical cases, the multiple integrals involved in the definition of p µ cannot be computed analytically while standard numerical integration methods [START_REF] Walter | Numerical methods and optimization: A consumer guide[END_REF] fail when n θ is larger than a few [START_REF] Andrieu | An introduction to MCMC for machine learning[END_REF]. Fortunately, the class of Metropolis-Hastings (MH) algorithms [START_REF] Robert | Monte Carlo statistical methods[END_REF] allows to sample a pdf on condition that this distribution is known up to a multiplicative constant. Indeed, with the MH algorithm, new samples θ (i) are generated from a user defined proposal distribution, then accepted or rejected according to a simple rule involving the acceptance probability ratio defined as follows

α (i) = p θ,µ (θ (i) | µ) p θ,µ (θ (i-1) | µ) = p θ,µ (µ | θ (i) )p θ (θ (i) ) p θ,µ (µ | θ (i-1) )p θ (θ (i-1) )
.

(2.24)

Because α (i) is a ratio of posterior distributions, it does not involve p µ . As a direct consequence for our model learning problem, a MH algorithm can be used hereafter to generate samples from this posterior distribution. Among the MH algorithms available in the literature (see, e.g., [START_REF] Luengo | A survey of Monte Carlo methods for parameter estimation[END_REF] for a recent list of MH solutions), the starting point of our solution is the random walk Metropolis algorithm [START_REF] Luengo | A survey of Monte Carlo methods for parameter estimation[END_REF]Section 3.1.1]. This specific choice is dictated by the facts that

• for the random walk Metropolis algorithm, the aforementioned proposal distribution can be a Gaussian pdf centered on the current candidate, i.e., the candidate parameter vector can be easily generated with standard random number generators once an initial parameter covariance matrix is selected by the user,

• this algorithm satisfies convergence constraints to guarantee that, asymptotically, the selected samples approach the correct expected values, i.e., share the same pdf and, more importantly, this pdf is • "thinning" [81, Section 9.4], i.e., removing samples in order to keep a chain of uncorrelated samples.

p θ,µ (θ | µ) (
Such steps are summed in Algorithm 1.

Algorithm 1 MCMC based estimation method

Inputs:

• θ lower , θ upper : parameter bounds

• M : chain lengths

• n mc : number of chains

• n M L : number of initial guesses for ML estimation Outputs:

• n mc estimates θ MCMC

• n mc estimates of µ max 1: Generate n M L initial guesses θ (0) by picking values randomly according to a uniform distribution between the parameter bounds 2: Estimate n M L parameter vectors θ ML with the ML method 3: Keep the parameter vector θ ML giving the smallest residual norm V (θ) defined in Eq. (2.14) and estimate the covariance matrix Σ ML defined in Eq. (2.15) 4: Initialize n mc chains at the value θ ML and the initial covariance matrix Σ (0) at the value Σ ML 5: Initialize S (0) to be the lower-triangular Cholesky factor of the initial covariance Σ (0) 6: for i = 1 to M do 7:

Sample a new candidate, i.e., θ * = θ (i-1) + S (i-1) r (i) , where r (i) ∼ N (0, I)

8:
Compute the acceptance probability α (i)

9:

Sample a uniform random variable u ∼ U(0, 1) 10:

If u ≤ α (i) , accept the sample and set θ (i) = θ * , otherwise, reject the sample and set θ (i) = θ (i-1)

Core of MH algorithm Maximum Likelihood 11:

Compute a lower-triangular matrix S (i) with positive diagonal elements satisfying the equation

S (i) S (i) ⊤ = S (i-1) I + η (i) (α (i) -ᾱ * ) r (i) r (i) ⊤ ||r (i) || 2 S (i-1) ⊤ ,
where η (i) ∈ [0, 1] is an adaptation step size sequence decaying to zero (e.g.

η (i) = i -γ with γ ∈ [0.5, 1]
) and ᾱ * = 0.234 the ideal acceptance probability 12: Thin the different chains to keep chains of uncorrelated samples 13: Compute an estimate θ MCMC for each chain using Equation (2.25) 14: Estimate n mc friction curves with the n mc estimates θ MCMC by using Eq. (2.26) 15: Determine the maximum of each friction curves Covariance adaptation from [START_REF] Särkkä | Bayesian Filtering and Smoothing[END_REF] Remark 2.2. As indicated in Algorithm 1, the random walk MH algorithm is performed on n mc Markov chains. By doing so, the method provides n mc estimates that can be gathered into histograms. However, after assessing that the different chains converge to the same distribution, it is common to select one chain and thus to have one mean estimation. One way to select the chain is to keep the chain which has the acceptance rate the closest to the ideal acceptance rate for the random walk MH defined by ᾱ * = 0.234 [START_REF] Särkkä | Bayesian Filtering and Smoothing[END_REF]. In the case treated here, the estimation procedure is applied with low friction measurements. For this reason, the estimates obtained are more likely to vary from one simulation to another. Thus it can be difficult to select one chain among all. Therefore, all the n mc Markov chains are kept in this study.

Once reliable chain samples are generated, the parameters estimation can be done using the expression [START_REF] Rogers | A first course in machine learning[END_REF]

θ MCMC = E p θ,µ (θ|µ) (θ | µ) = • • • θp θ,µ (θ | µ)dθ ≈ 1 M n mc M nmc j=1 θ (j) .
(2. [START_REF] Dugoff | Tire performance characteristics affecting vehicle response to steering and braking control inputs[END_REF] In others words, the parameters estimation is made by averaging the available samples. Finally, the prediction on an unseen point {s M +1 , µ M +1 } can be carried out as follows [START_REF] Rogers | A first course in machine learning[END_REF]

µ M +1 = E p θ,µ (θ|µ) (f (θ, s M +1 ) | µ) ≈ 1 M n mc M nmc j=1 f (θ (j) , s M +1 ).
(2.26) Because the slip ratio corresponding to the tire-road friction µ max is not known a priori, Equation (2.26) cannot be used directly to estimate µ max . However, this equation can be used with many slip ratio to draw the entire friction curve. Once the curve is obtained, µ max is determined straightforwardly by taking the maximum of the friction curve or by examining when the derivative of the friction curve with respect to the slip ratio is equal to zero.

Case studies

In this section, the MCMC model learning technique performance is assessed with both simulated data (Section 2.4.1) and real data (Section 2.4.2). Because in Section 2.4.1, the measurements are simulated, we have access to the real parameter vector θ. Thus the estimated values obtained with the MCMC3 method can be compared to the actual values directly. For this reason, Section 2.4.1 focuses on the ability of the method to estimate the parameter vector θ rather than the friction potential. If θ is estimated accurately, the estimated grip potential determined, e.g., by computing the derivative ∂f ∂s , should indeed be reliable as well. In Section 2.4.2, we use on the contrary real measurements of the friction potential generated with a flat track tire testing machine [START_REF] Cabrera | A versatile flat track tire testing machine[END_REF]. Under such practical conditions, we do not have access to the real parameter vector θ. Hence, the method performance is assessed by comparing the estimated grip potential with a grip potential extracted from the real measurements only.

Results with simulated data

In order to evaluate the effectiveness of the estimation method introduced in this chapter, the MCMC model learning method is first tested with simulated data. The results obtained with the MCMC model learning method are compared with the ones obtained with the ML method carried out on the same simulated data. The simulated data is generated with the Pacejka's tire model given in Eq. (2.5) disturbed with an output zero mean white Gaussian noise characterized by a standard deviation σ equal to 2.53e-2 (see Figure 2

Results with the Maximum Likelihood method

The ML method is first tested with the aforementioned simulated training sets. In the different cases treated here, the ML method is performed 10000 times (n M L = 10000 in Algorithm 1) with 10000 different initial guesses. As indicated in Algorithm 1, the different initial guesses θ (0) are generated randomly according to a uniform distribution between the parameter bounds. Here, the bounds, given in Table 2.3, are taken so that the Pacejka's model can generate realistic friction curves corresponding to any kind of situations such as a vehicle driving on dry asphalt, wet asphalt or even icy road [START_REF] Olazagoitia | Identification of tire model parameters with artificial neural networks[END_REF]. As illustrated in Figure 2.7, 2.8 and 2.9, the means of the 10000 estimates are different from the minimizers of the loss function. In addition, the mean estimates are closer to the actual parameter values than the minimizers of the loss function. However, in the first two cases (measurements limited by µ < 0.3 and µ < 0.6, respectively), the mean estimates and the minimizer of the loss function are far from the actual values. On the contrary, when all the data points are considered, the ML method provides accurate estimates of the parameters. Besides, in this case, the measurement standard deviation supplied by the ML method is equal to the actual value used to simulate the data (see Table 2.5). Furthermore, in the first two cases, the Jacobian matrices F θ ML given in Eq. (2.16) are ill-conditioned. For this reason, the theoretical covariance matrices (see Eq. (2.15)) computed by using the Jacobian matrices F θ ML , contain huge values as reported in Table 2.6. On the opposite, when the last training set is considered, the Jacobian matrix F θ ML is not ill-conditioned. Thus, in this case, the theoretical covariance matrix has realistic values and is almost equal to the covariance matrix computed with the 10000 ML estimates.

All these results indicate that the ML method is suitable when data points describing the entire friction curve are available. However, when only training sets containing low friction measurements are at disposal, the ML method provides poor results. For this reason, another estimation method is needed under standard driving conditions.

Results with the MCMC method

The MCMC based estimation method is now tested with the different training sets described in the paragraph before the Section 2.4.1. As indicated in Algorithm 1, the MCMC model learning method requires the initialization of some quantities. Namely, the chain lengths, the chain numbers and the parameter bounds. Herein, the MCMC method developed in Section 3 is performed with 100 chains to obtain 100 estimates of θ. Each chain has a length of 50 000 samples. Besides the parameter bounds are the ones given in Table 2.3. Algorithm 1). However, in our case, the Jacobian matrices F θ ML obtained with the ML method carried out on the two first training sets are ill-conditioned (see Section 2.4.1). As a consequence, the theoretical ML covariance matrices estimated in Section 2.4.1 are not suitable for initializing the covariance matrix of the MCMC proposal distribution in these two cases. The proposal covariance matrix should indeed be set so that the generated samples browse the parameters space properly [START_REF] Rogers | A first course in machine learning[END_REF]. The most important factor influencing the sample generation is the variance of the parameters. In the proposal covariance matrix, the parameter variances are representing by the diagonal elements. For this reason, the proposal covariance matrix is chosen diagonal. Besides, the magnitude of the diagonal elements is set regarding the ML covariance matrix estimated when all the data points are selected. Finally, the diagonal elements are set to higher values than the ML matrix. The diagonal values are higher so that the generated samples cover a 2.7 correspond to the mean of the 100 estimates. In addition, Gaussian pdf have been fitted to the MCMC estimates. The means and the standard deviations used to fit the pdf correspond to the means and the standard deviations of the 100 estimates. As illustrated in Figures 2.10, 2.11 and 2.12 (see also Table 2.7), the MCMC model learning method outperforms the ML method for all scenarii. For the three different cases, the MCMC estimates are indeed closer to the true values than the ones obtained with the ML method. Besides, we can note that the MCMC estimates obtained with a friction Table 2.6: Parameters covariance matrices obtained with the ML method. The theoretical column corresponds to the covariance matrices computed with Eq. (2.15). The empirical covariance matrices corresponds to the covariance matrices computed with the 10000 estimates.

Remark 2.3. At this point, the previous settings should be sufficient for performing the MCMC based estimation method (see

Theoretical Empirical

Friction limit : µ = 0.30

        59679083
               
22.2905 0.0366 -0.6287 -0.0111 -0.0048 -0.0122 0.0366 0.1500 -0.0891 0.0894 0.0011 -0.0297 -0.6287 -0.0891 0.0823 -0.0642 -0.0005 0.0166 -0.0111 0.0894 -0.0642 0.4112 0.0041 -0.0799 -0.0048 0.0011 -0.0005 0.0041 0.0001 -0.0013 -0.0122 -0.0297 0.0166 -0.0799 -0.0013 0.0278

        Friction limit : µ = 0.60         215526.
               
25.4821 -0.1198 -1.2810 -2.5468 -0.0190 0.3863 -0.1198 0.0119 -0.0040 0.0066 0.0000 -0.0008 -1.2810 -0.0040 0.0754 0.1309 0.0010 -0.0212 -2.5468 0.0066 0.1309 0.2759 0.0018 -0.0359 -0.0190 0.0000 0.0010 0.0018 0.0000 -0.0003 0.3863 -0.0008 -0.0212 -0.0359 -0.0003 0.0069

        Without friction limit         2.3744
0.0952 -0.1696 0.5560 -0.0084 0.1685 0.0952 0.0039 -0.0070 0.0224 -0.0003 0.0069 -0.1696 -0.0070 0.0124 -0.0396 0.0006 -0.0124 0.5560 0.0224 -0.0396 0.1335 -0.0020 0.0393 -0.0084 -0.0003 0.0006 -0.0020 0.0000 -0.0006 0.1685 0.0069 -0.0124 0.0393 -0.0006 0.0123

               
2.1736 0.0774 -0.1844 0.6428 -0.0088 0.1841 0.0774 0.0028 -0.0066 0.0229 -0.0003 0.0066 -0.1844 -0.0066 0.0157 -0.0546 0.0008 -0.0156 0.6428 0.0229 -0.0546 0.1905 -0.0026 0.0545 -0.0088 -0.0003 0.0008 -0.0026 0.0000 -0.0007 0.1841 0.0066 -0.0156 0.0545 -0.0007 0.0156

       
limit of 0.3 are closer to the actual values (except for B) than the ones obtained with a friction limit of 0.6.

To conclude this part, the MCMC model learning technique provides interesting results to estimate the parameters of the Pacejka's tire model and outperforms the ML method. Thus, it is now time to test the ML and the MCMC methods on real data with realistic friction measurement sets.

Results with real data

In this section, the method described in Section 3 is tested with real data coming from a flat track tire testing machine. A flat track tire testing machine is a machine allowing the user to conduct experiments on a tire by applying specific constraints like a desired normal load or by rolling the tire at a desired speed [START_REF] Cabrera | A versatile flat track tire testing machine[END_REF]. Therefore, a flat track tire testing machine offers the possibility to collect useful information on the tire behavior by studying its dynamical responses. Among the signals available thanks to the flat track tire testing machine, we can have access to the slip ratio s as well as the tire forces F x , F y and F z respectively. Hence, by using Eq. (2.1), we obtained the friction measurements given in Figure 2.13. In Figure 2.13, it can be pointed out that the friction measurements depict the entire friction curve which is not available under standard practical driving conditions. Here, all these data points are accessible because the flat track offers the possibility to impose high slip ratios to the tire and thus allows the acquisition of the tire response for a large range of slip ratio. As a consequence, these representative measurements offer the opportunity to test the MCMC based estimation method and the ML method on different part of the friction curve.

Tests with full measurements

In the first scenario, the method is applied using the full data sets for s ∈ [0, 0.2]. Although this situation does not correspond to classical driving conditions, this test enables to check if the method can provide good results under ideal situations. In the same way as we proceeded in the case with simulated data, the MCMC based estimation method is applied with a high number of chains to ensure a satisfying efficiency of the method. Here the number of chains is set to 300 which imply 300 estimates of the grip potential. Each chain has a length of 30000 samples. In addition, as stated in Algorithm 1, the different chains are initialized with the parameter vector obtained with the ML method. Herein, the ML method is carried out 10000 times similarly to what have been done with the simulated data. Furthermore, the parameter bounds used are still the ones indicated in Table 2.3.

Remark 2.4. In the case treated here, the estimated theoretical covariance matrix obtained with the ML method is

ĉov(θ ML ) =         950
        , (2.28)
which contains high diagonal elements. As a consequence, this matrix is not suitable for initializing the covariance matrix of the proposal distribu- ML method have been carried out 10000 times with 10000 different initializations. However, in this case, the resulting grip potential estimates are almost always the same. For this reason, the grip potential obtained with the minimizer of the loss function and the mean of the estimates are approximately equals and thus, are not distinguishable.

For the same reason, no ML pdf have been plotted in Figure 2.15. The friction curve estimated with the ML method is represented in Figure 2.14. The means and the standard deviations of the different ML and MCMC estimates are summarized in Table 2.8. As illustrated in Figure 2.14, the MCMC estimates match well to the friction measurements. On top of that, the mean of these estimates is almost similar to the one obtained with the ML method. This assertion is reinforced by Table 2.8 and the histograms in Figure 2.15 where the fric- tion potential are estimated between 1.15 and 1.3 with a mean equal to 1.222 while the ML value is 1.223. As shown in the second histogram in Figure 2.15, the slip ratio corresponding to the grip potential estimated with both ML and MCMC methods are close to each other. Indeed, the slip ratio estimated with the MCMC method are situated between 0.03 and 0.055 for a mean value of 0.0425 while the ML value is 0.0397. These results indicate that the MCMC based estimation method works well to estimate the grip potential in the case treated here. More precisely, it delivers mean results similar to the ML method.

Finally note that simulating one chain of 30000 samples has taken around 7.5 seconds. Experiments were conducted with Matlab R2018b on a computer with an Intel core i7 processor running at 2.0 GHz.

Because the ML and the MCMC methods provide an accurate similar estimate of the friction curve when all data points are taken into account, this estimate is now used as a reference. This reference will be used in the following to compare the results obtained with different training sets.

Tests with realistic scenarii

The ML and MCMC estimation methods are now tested with different training sets composed of the data points represented in Figure 2.13. Each training set is restricted by a different friction limit. The different friction limits start with µ = 0.05 and grow until µ = 0.7 by step of 0.05 for a total of 14 different training sets. Under these practical conditions, the smallest training set contains 24 points and the largest 112 points. Once the training sets are created, it remains to proceed to the initialization of the different methods. For both methods, the parameter bounds are still the ones given in Table 2.3. As in the previous cases, the ML method is first performed by picking initials guesses randomly between the parameter bounds. In the different cases treated here, the ML method is performed with 1000 different initial guesses. Furthermore, among all the 1000 ML estimates, the one kept is the parameter vector θ ML giving the smallest residual norm. For the MCMC model learning technique, it remains to define the chain numbers and chain lengths. Here, because the training sets are reduced compared to the one considered in the previous case, the estimates obtained with the MCMC method are more likely to vary. For this reason, the number of chains, and thus the estimates, has been set to 1000 with each chain having a length of 100000 samples. Each chain is initialized with the ML estimate θ ML .

Under such conditions, the MCMC method gives 1000 estimates for each of the 14 training sets. The method performance is assessed by comparing, for the 14 training sets, the relative error between the grip potential of the mean of the MCMC estimates and the grip potential reference. The grip potential reference is the one obtained in the Section 2.4.2 when the ML and the MCMC methods are used with all the available measurements. The same comparison is made with the slip ratio corresponding to the grip potential. Figure 2.16 shows the evolution of the estimated grip potential with the different estimation methods. Figure 2.16 also includes the evolution of the relative error of the estimates with the reference value. Figure 2.17 contains the same type of results about the slip ratio corresponding to the grip potential, i.e. the evolution of the estimates and the relative errors. In addition to these curves, Figures 2.16 and 2.17 contain a third curve named "MCMC estimate with physical prior" where physical prior is added on MCMC estimates. More precisely, in the situation considered here, we know that the corresponding slip ratio value of the grip potential cannot physically be too high. For this reason, after making the 1000 MCMC estimates, we keep the ones indicating a grip potential with a corresponding slip ratio lower than 0.1 and we consider the mean of these estimates. In this way, the physical prior can be combined with the MCMC based estimation method to improve the final grip potential estimation. Figure 2.18 indicates the percentage of MCMC estimates rejected after adding physical prior. Figures 2.19 and 2.20 show the estimated friction curves obtained in two practical cases. More precisely, when the friction measurements available remain under µ < 0.2 and µ < 0.3, respectively. Finally, Figures 2.21 and 2.22 display plots of iterations versus sampled values for each variable in the MCMC chain when added physical prior is considered. The plots are the ones obtained when the used friction measurements remain under µ < 0.2 and µ < 0.3, respectively.

As indicated in Figure 2.16 and 2.17, adding physical prior to the MCMC estimates improves the results. Besides, once the training set friction limit is greater than 0.15, the MCMC estimates with physical prior bring more accurate results than the ML method. Indeed, once this friction limit is reached, the MCMC based estimation method provides an estimate of the grip potential with a relative error lower than 20% while the relative error of the estimate obtained with the ML method is almost always greater than 20%. This result is interesting because as mentioned in Section 1, most of the methods used to estimate the grip potential are effective once the grip consumption level is greater than 0.3 [START_REF] Acosta | Road friction virtual sensing: A review of estimation techniques with emphasis on low excitation approaches[END_REF]. Therefore, estimating the grip potential with a lower friction consumption represents an improvement. Besides, as shown in Figure 2.18, the percentage of MCMC estimates rejected after adding physical prior decreases with the friction limit. This result suggests that the estimates obtained are more reliable once the friction limit increase. Furthermore, as depicted in Figures mix well [START_REF] Andrieu | An introduction to MCMC for machine learning[END_REF]. More precisely, on the one hand, the samples are not stuck in a small part of the parameter space. On the other hand, the samples are not far from each other. Thus, these traceplots indicate that the generated samples explore the parameter space very well. In order to evaluate the method performance more deeply, the histograms of the MCMC estimates of the grip potential and its corresponding slip ratio have been plotted in Figure 2.23 and 2.24, respectively. Figure 2.23 indicates that the histograms begin to be centered on the reference once the friction limit is equal to 0.2 which corroborates the fact that the method provides satisfactory results from this friction limit measurement value. However, the histogram corresponding to a friction limit of 0.2 in Figure 2.24 is not clearly centered on the reference value and it does not contain values that emerge from all the estimates obtained. Thus, with this friction limit, the grip potential estimate can be accurate but not necessary the associated slip ratio. This situation is not necessarily an issue because we try to estimate the grip potential in priority.

Finally simulating one chain of 100000 samples have taken around 7.5 seconds for the smallest training set and 10 seconds for the largest one.

Characterization of real measurements

As indicated in the different cases treated in this part, with real data, the MCMC estimates are closer to the reference values than the estimates obtained with the ML method. Looking at the results, we can wonder why the MCMC method works better than the ML method. One reason that could explain these results is the validity of the assumptions made for both methods.

When we described the ML method in Section 2.3.1, we stated that the stochastic part of the measurements can be modelled with i.i.d. Gaussian stochastic components. Now that we have access to real mea- surements and a reliable estimate of the friction curve (given by ML and MCMC), the residuals can be computed (see Figure 2.25), then analysed in order to validate the aforementioned assumptions. Two tests are conducted on the residuals, a normality test [START_REF] Thode | Testing for normality[END_REF] and an autocorrelation test (ACF) [START_REF] Rogers | A first course in machine learning[END_REF]. These tests are carried out on two different training sets in order to ensure that all the measurements have the same statistical properties. The first training set includes the residuals corresponding to friction measurements lower than 0.6. The second data set contains all the residuals.

Firstly, the normality test is performed to verify if the residuals can be model with a Gaussian distribution. The normality test used here is the Kolmogorov-Smirnov test at a significance level of 5% [START_REF] Thode | Testing for normality[END_REF]. This one indicates that the null hypothesis cannot be rejected for both residual sets. Besides, this test indicates a p-value of 0.8960 when only a part of the residuals are considered and a p-value of 0.5339 when all the residuals are selected. Thus, considering the noise normally distributed is a suitable assumption.

Secondly, the ACF test is carried out to verify if the residuals can be considered as independent. As illustrated in Figure 2.26, it appears that no matter which residual sets is considered, the data points are not i.i.d.. Therefore, modelling the stochastic part of the measurements with i.i.d. Gaussian stochastic components might not be the most appropriate choice. Furthermore, because the i.i.d. assumptions is not valid, it can explained why the ML method provides estimates far from the actual value. The MCMC method seems less sensitive to these assumptions on the measurements. This difference in the sensitivity could explained why the MCMC method provides better estimates than the ML method.

Conclusions

In this chapter a method combining the Maximum Likelihood approach and an Adaptive Metropolis-Hastings MCMC procedure is applied to determine the tire friction potential. In this method, the ML estimates are used to initialize an Adaptive Metropolis-Hastings MCMC algorithm. This method is tested with different training sets coming from both simulated data and real data. The real data is obtained from experiments conducted with a flat track tire testing machine. When data points describing the entire friction curve are selected (µ ∈ [0, 1.2]), the ML and MCMC methods give accurate predictions of the Pacejka's model coefficients. However, under standard driving conditions, only friction measurements lower than 0.3 are available. For this reason, the estimation methods are assessed with realistic training sets containing friction measurements under 0.3 only. Under these conditions, the results reveal that the MCMC solution outperforms the classical Maximum Likelihood approach. Indeed, this approach produces accurate estimates as soon as friction measurements reach a value of approximately 0.2. Thus, this result makes the method applicable during standard driving conditions where this range of measurements is accessible. These results are obtained by assuming that the friction data points are actually measured. However, in practice, the friction data points are not measured but estimated. For this reason, it is necessary to adapt this solution when the friction data is estimated (see Figure 2.27).

CHAPTER 3 Tire force estimation

As mentioned in the past sections, the objective of this thesis is to determine the grip potential defined in Eq. (2.1) under standard driving conditions. In order to reach this goal, it is necessary to use friction points. Although friction points measurements were used in Chapter 3, these ones are not available with sensors fitted on production vehicle. Therefore, they should be estimated and this is the main objective of this chapter (See Figure 3.1). The friction data points are deduced through estimates of the longitudinal tire force F x , the normal load F z and their corresponding slip ratio s. These estimates are obtained by resorting to a state observer [START_REF] Simon | Optimal state estimation: Kalman, H ∞ and nonlinear approaches[END_REF]. To achieve this, it is necessary (i) to choose a vehicle model (ii) to define an observer structure. Hence, the chapter is organized as follows. Section 3.1 describes the problem tackled in this chapter, then introduces the main notations used in the next paragraphs. Section 3.2 is devoted to the description of the used vehicle model. Section 3.3 is dedicated to the observer description and its settings.

Problem formulation and notations

In this chapter, we aim at extracting friction information from available measurements. Because friction data points cannot be measured directly with the available sensors [START_REF] Acosta | Road friction virtual sensing: A review of estimation techniques with emphasis on low excitation approaches[END_REF], the friction data information must be reconstructed from measurements. Among the different solutions available in the literature [START_REF] Acosta | Virtual tyre force sensors: An overview of tyre model-based and tyre model-less state estimation techniques[END_REF], a specific attention is paid to a state observer because it requires only a model of the considered system and measurements coming from sensors fitted on production vehicle. In control theory, a state observer is indeed a process combining the inputs u, the measurements z and a dynamical system representation in order to infer estimates of the internal states x of the system [START_REF] Simon | Optimal state estimation: Kalman, H ∞ and nonlinear approaches[END_REF][START_REF] Kailath | Linear Estimation[END_REF]. This concept is illustrated in Figure 3.2. As indicated in Figure 3.2, the state observer uses a theoretical description of the system to estimate the state. Therefore, in order to apply a state observer, it is necessary to define • a model representation of the considered system, i.e., a mathematical description of the system dynamics we want to estimate,

• a state observer structure.

In the case treated here, the studied system is a vehicle car. Therefore, the considered model should describe the dynamics of a car over the time. In automotive field, such a model is called a vehicle model. Besides, in our case, the objective is to determine friction points through estimates of the longitudinal tire force F x , the normal load F z and the slip ratio s. Thus, the selected vehicle model should include these quantities or quantities related to them. The used vehicle model and the observer structure are described in Section 3.2 and in Section 3.3, respectively.

Notations: In the following, the tire forces are noted with the form F uij where the subscript u ∈ {x, y, z} designates the longitudinal, lateral and normal tire force, respectively. The subscript i ∈ {f, r} stands for front and rear respectively and the subscript j ∈ {l, r} stands for left and right respectively. Thus, for example, with this notation, the quantity F xf l indicates the front left longitudinal tire force. More generally, in this manuscript, others quantities are denoted with the subscripts i ∈ {f, r} and j ∈ {l, r}. In any case, i and j make the distinction between front and rear and left and right, respectively.

Description of the vehicle model

A vehicle model is a set of mathematical equations used to mimic the behavior of a vehicle. Depending on what part of the vehicle is studied, the model can be more or less complex. For this reason, many different types of model exist in the literature [START_REF] Guiggiani | The Science of Vehicle Dynamics[END_REF], [START_REF] Rajamani | Vehicle Dynamics and Control[END_REF], [START_REF] Kiencke | Automotive control systems: for engine, driveline, and vehicle[END_REF]. Among the different vehicle models, we distinguish the kinematic models from the dynamic models. The kinematic models do not include the forces acting on the vehicle while the dynamic models include them [START_REF] Guiggiani | The Science of Vehicle Dynamics[END_REF]. Besides, depending on what sort of vehicle behavior is studied (longitudinal, lateral or coupled), the vehicle model is viewed as a longitudinal model, lateral model or coupled model, respectively.

In the present study, the main objective is to estimate the longitudinal tire force F x , the normal load F z and the slip ratio s. Thus, the vehicle model should include these quantities or quantities related to them in the set of equations. As a result, the model taken into account is a dynamical longitudinal model. Moreover, in this study, we want a model sufficiently accurate to get reliable estimates of the tire forces and at the same time sufficiently simple to avoid a heavy computational cost. For these reasons, the model chosen in this study is a single-track model [START_REF] Guiggiani | The Science of Vehicle Dynamics[END_REF] considering only the longitudinal dynamics extended with a suspension model. This choice of vehicle model should represent a good trade-off between accuracy and complexity. Indeed, on the one hand, the single-track model is used to describe the longitudinal tire force F x and on the other hand, the suspension model is used to represent the normal load F z and especially to capture the variations of F z . Because the single-track is a simplification of the double-track model [START_REF] Guiggiani | The Science of Vehicle Dynamics[END_REF], it is necessary to introduce the double-track model before presenting the single-track model. Consequently, the double-track model is presented in Section 3.2.1, then the single-track model is presented in Section 3.2.2. The suspension model is presented in Section 3.2.3. Besides, we are also interested in the slip ratio. As presented in Section 3.2.4, this quantity depends of the effective tire radius model. Thus, Section 3.2.5 is dedicated to the effective tire radius model. Finally, all the equations of the vehicle model considered are summarized in Section 3.2.6.

Double-track model

The double-track model is a vehicle model commonly used to describe the vehicle dynamics [START_REF] Guiggiani | The Science of Vehicle Dynamics[END_REF][START_REF] Ray | Nonlinear tire force estimation and road friction identification: Simulation and experiments[END_REF][START_REF] Goldfain | Autorally: An open platform for aggressive autonomous driving[END_REF]. In this model, the four wheels of the vehicle are represented separately (see Figure 3.3 and Table 3.1) in order to get closer to a real car. As a consequence, each wheel is related to three different tire forces F x , F y and F z that should be taken into consideration. However, in the common double-track model, the pitch and roll dynamics are neglected [START_REF] Rajamani | Vehicle Dynamics and Control[END_REF]. Therefore, the different normal loads F zif applied on each wheel are considered as constant values.

Model equations

The model equations of the double-track model are obtained by applying the Newton's second law to the system composed of the vehicle body and the four wheels [START_REF] Guiggiani | The Science of Vehicle Dynamics[END_REF][START_REF] Rajamani | Vehicle Dynamics and Control[END_REF]. In this case, we obtained the equations

vx = v y ψ + 1 m [X f + X r -F aero ] , vy = -v x ψ + 1 m [Y f + Y r ] , ψ = 1 I z L f Y f -L r Y r + t f 2 ∆X f + t r 2 ∆X r , (3.1a) (3.1b) (3.1c)
f = (F xf l cos(δ f l ) + F xf r cos(δ f r )) -(F yf l sin(δ f l ) + F yf r sin(δ f r )) , X r = (F xrl cos(δ rl ) + F xrr cos(δ rr )) -(F yrl sin(δ rl ) + F yrr sin(δ rr )) , Y f = (F xf l sin(δ f l ) + F xf r sin(δ f r )) + (F yf l cos(δ f l ) + F yf r cos(δ f r )) , Y r = (F xrl sin(δ rl ) + F xrr sin(δ rr )) + (F yrl cos(δ rl ) + F yrr cos(δ rr )) , ∆X f = (F xf r cos(δ f r ) -F xf l cos(δ f l )) + (F yf l sin(δ f l ) -F yf r sin(δ f r )) , ∆X r = (F xrr cos(δ rr ) -F xrl cos(δ rl ))
+ (F yrl sin(δ rl ) -F yrr sin(δ rr )) .

(3.2a)

(3.2b) (3.2c) (3.2d) (3.2e) (3.2f)
As indicated in Equations (3.1a)-(3.2f), the double-track model is a vehicle model involving many different quantities. However, in order to carry out state estimation, it is necessary to dispose of wheel measurements (see Section 3.3). Therefore, using this model requires four different sensors only for the wheels since all the four vehicle wheels are represented separately. Besides, modelling the four wheels separately results in a high dimensional state vector since it should contain at least the four longitudinal tire forces and the four wheel speeds (see Section 3.3). As indicated in [START_REF] Berntorp | Trajectory tracking for autonomous vehicles on varying road surfaces by friction-adaptive nonlinear model predictive control[END_REF] and [START_REF] Fang | Nonlinear bayesian estimation: From kalman filtering to a broader horizon[END_REF], the double-track model is an accurate model which describes the vehicle dynamics with a high accuracy but this amount of accuracy is not necessarily required in the case of state estimation. For these reasons, the double-track model might not be a suitable choice for the problematic treated here and another less complex model should be investigated.

Single-track model

The single-track model (see Figure 3.4 and Table 3.2) is another vehicle model particularly used to represent the vehicle dynamics [START_REF] Guiggiani | The Science of Vehicle Dynamics[END_REF], [START_REF] Rajamani | Vehicle Dynamics and Control[END_REF]. It is a simplification of the double-track model. In addition to the double-track model assumption (pitch and roll effect neglected), the single-track model requires the following simplifying assumptions

• the left and right front steering angles are equal ⇒ δ f l = δ f r = δ,

• the rear steering angles are equal to zero ⇒ δ rl = δ rr = 0.

The first assumption on the equal left and right steering angles leads to the equality between the left and right tire forces. Consequently, in the single-track model, the left and right front, respectively rear, wheels are gathered to form one front, respectively rear, wheel. As a result, the tire forces are brought back to each axle which means that, instead of estimating four longitudinal tire forces (one for each wheel), we will only estimate two longitudinal tire forces, one for the front axle, F xf and one for the rear axle, F xr . Besides, in the purely longitudinal singletrack model, the axle tire forces are the sum of the left and right tire forces. As a consequence, we obtain the equations

F xf = F xf l + F xf r , F xr = F xrl + F xrr , F yf = F yf l + F yf r , F yr = F yrl + F yrr , F zf = F zf l + F zf r , F zr = F zrl + F zrr . (3.3a) (3.3b) (3.3c) (3.3d) (3.3e) (3.3f)
Hence, by making the assumption that the tire forces are the same in the left and right side, the axle tire forces are enough to estimate the tire forces on each wheel.

Model equations

The equations of the single-track model are obtained by applying Newton's second law to the different parts composing the single-track model, that is, the vehicle body and the front and rear tires. In order to apply Newton's second law, it is necessary to described the forces acting on the vehicle. These forces are represented in Figures 3.4 As indicated in [START_REF] Rajamani | Vehicle Dynamics and Control[END_REF], the equations of the single-track model consid-ering only the longitudinal dynamics are

vx = 1 m [F xf + F xr -F aero -(R xf + R xr )] , ωf = 1 2I wf [T f -R f load F xf ] , ωr = 1 2I wr [T r -R rload F xr ] , (3.4a)
(3.4b) (3.4c)
where

F aero = 1 2 ρ a S a C x (v wind + v x ) 2 , R xf + R xr = f RR (F zf + F zr ) , R f load = R 0 - F zf K zz , R rload = R 0 - F zr K zz . (3.5a) (3.5b) (3.5c) (3.5d)
The signals and parameters involved in the model equations are described in Table 3.4. Because the pitch dynamics are not considered, the

F zf = F zf 0 = mg L r L f + L r , F zr = F zr0 = mg L f L f + L r . (3.6a) (3.6b) Thus F zf + F zr = mg and R xf + R xr = f RR mg.
Besides, because the wind speed is generally unknown, it is neglected in the Equation (3.5a). Under these practical conditions, the model equations become

vx = 1 m F xf + F xr - 1 2 ρ a S a C x v 2 x -f RR mg , ωf = 1 2I wf [T f -R f load F xf ] , ωr = 1 2I wr [T r -R rload F xr ] . (3.7a) (3.7b) (3.7c)

Suspension model

As mentioned in Section 3.2.2, in the single-track model, the pitch and roll effects are neglected. This assumption leads to static expressions of the normal loads F zi . However, the goal of this chapter is to estimate friction data points which depend of the normal load. For this reason, a great importance should be paid to the estimation of F zi in order to obtain an accurate estimate of the friction µ. Thus, it is necessary to consider the normal loads with dynamic expressions and not static ones. One way to represent these dynamics is to extend the single-track model with a suspension model which will be used to characterize the load transfer. The load transfer represents the evolution of the normal load which occurs during a travel. For example, during acceleration and braking phases, a load transfer between the front and the rear parts of the vehicle takes place. This load transfer is due to the pitch effect (see Figure 3.6 for an illustration of the pitch effect). Thus, during these phases, the normal load applied on the tire varies and these variations are described by using a suspension model.

Simplified suspension model

In order to include the pitch dynamics, it is necessary to introduce a suspension model. Many models exist in the literature to represent the vehicle suspensions [START_REF] Poussot-Vassal | Robust LPV multivariable automotive global chassis control[END_REF][START_REF] Doumiati | Road profile estimation using an adaptive youla-kučera parametric observer: Comparison to real profilers[END_REF]. One simplified passive suspension model is the one represented in Figure 3.7. The parameters used in this model are summarized in Table 3.5.

∆F zf = K f L f sin(κ) + B f L f κ cos(κ), ∆F zr = K r L r sin(κ) + B r L r κ cos(κ). (3.8a) (3.8b)
As a result, the normal loads become

F zf = F zf 0 + [K f L f sin(κ) + B f L f κ cos(κ)] , F zr = F rf 0 -[K r L r sin(κ) + B r L r κ cos(κ)] , (3.9a) (3.9b) with F zf 0 = mg L r L f + L r , F zr0 = mg L f L f + L r . (3.10a) (3.10b)
As indicated in Eq. (3.9a) and (3.9b), respectively, it is necessary to describe the evolution of the pitch angle κ and the pitch speed κ in order to compute the normal loads. By applying Newton's second law at the CoG in the model represented in Figure 3.7, we can establish the following equation governing the pitch dynamics

I y κ = (F xf + F xr -R xf -R xr )h G -F zf L f + F zr L r . ( 3
κ + κ cos(κ) I y [B f L f (f RR h G + L f ) + B r L r (L r -f RR h G )] + sin(κ) I y [K f L f (f RR h G + L f ) + K r L r (L r -f RR h G )] = h G I y m vx + 1 2 ρ a S a C x v 2 x .
(3.12)

Thus, with the suspension model represented in Figure 3.7, the pitch dynamics are governed by a second order nonlinear differential equation with nonconstant second member. This equation can be simplified by making assumptions on the pitch angle and the spring and damper parameters. During a regular travel, the pitch angle of a car stays relatively small (κ < 10 • ) [START_REF] Berntorp | Trajectory tracking for autonomous vehicles on varying road surfaces by friction-adaptive nonlinear model predictive control[END_REF]. In this case, sin(κ) ≃ κ and cos(κ) ≃ 1. Besides, by adding the assumption of equilibrium of load transfer between front and rear axles, we have

K f L f = K r L r = k, B f L f = B r L r = b. (3.13a) (3.13b)
Hence, the normal loads become

F zf = F zf 0 + ∆F z , F zr = F zr0 -∆F z , (3.14a) (3.14b)
with ∆F z = kκ + b κ. Furthermore, the equation governing the pitch dynamics becomes

κ + Lb I y κ + Lk I y κ = h G I y m vx + 1 2 .ρ a S a C x v 2 x . (3.15)
Thus, with the additional assumptions of small pitch angle and equilibrium of load transfer, Equation (3.15), which governs the pitch dynamics, becomes an ordinary second order linear differential equation with nonconstant second member. Therefore, it can be put into a canonical form

κ + ω 0 Q κ + ω 2 0 κ = h G I y m vx + 1 2 .ρ a S a C x v 2 x , (3.16)
with

ω 0 = Lk I y , Q = 1 b I y k L . ( 3
κ(t) = κ 1 (t) + κ 2 (t), (3.18)
where

κ 1 = exp -ω 0 t 2Q (A cos(Ωt) + B sin(Ωt)) , Ω = ω 0 1 - 1 4Q 2 , κ 2 = h G 2Lk ρ a S a C x v 2 x , (3.19a) (3.19b) (3.

19c)

with A and B given by the initials conditions. If the initial conditions are κ(0) = κ ini and κ(0) = κini , the constants A and B are

A = κ ini -κ 2 , B = 1 Ω κini + ω 0 2Q A . (3.20a) (3.20b)

Slip ratio expression

In order to determine friction data points, it is necessary to estimate the slip ratio s corresponding to the friction µ of the different points. The definitions of the front and rear slip ratios are given by [START_REF] Rajamani | Vehicle Dynamics and Control[END_REF]

s f = ω f R rol f -v x max(ω f R rol f , v x ) , s r = ω r R rolr -v x max(ω r R rolr , v x ) . (3.21a) (3.21b)
Therefore, computing the slip ratio requires the knowledge of the effective tire radius R rol i (i ∈ {f, r}), the longitudinal vehicle speed v x and the wheel speed ω i (i ∈ {f, r}). While the vehicle speed and the wheel speed are common information available on production cars, this is not the case for the effective tire radius. For this reason, it is necessary to introduce an effective tire radius model.

Effective tire radius model

During a travel, the tire is distorted by the normal load applied on it. As a result, this deformation induces a modification of the effective tire radius. Different models exist to represent the evolution of the effective tire radius [START_REF] Kiencke | Automotive control systems: for engine, driveline, and vehicle[END_REF], [START_REF] Rajamani | Vehicle Dynamics and Control[END_REF], [START_REF] Pacejka | Tire and vehicle dynamics[END_REF]. Among them, two effective tire radius models are particularly employed, (i) the sinc effective tire radius model [START_REF] Kiencke | Automotive control systems: for engine, driveline, and vehicle[END_REF], (ii) the Pacejka effective tire radius model [START_REF] Pacejka | Tire and vehicle dynamics[END_REF]. The sinc model is a model giving an expression of the effective tire radius as a function of the free-tire radius R 0 , the normal load F z applied on the tire and the stiffness K zz . This tire model is based on the relation

v x = R rol i ω i , i ∈ {f, r}. (3.22)
However, this relation corresponds to situations where the slip ratio is equal to zero. Indeed, by inserting Equation (3.22) in Eq. (3.21a) and Eq. (3.21b), the resulting front and rear slip ratios are equal to zero. However, the goal of this chapter is to estimate the friction points defined by a friction µ and a slip ratio s. Therefore, the effective tire radius model used to estimate the slip ratio should operate in situations where the slip ratio is different from zero. For this reason, the sinc effective tire radius model is not chosen here.

The second effective tire radius model well documented in the literature is the semi-empirical model developed by Pacejka in [START_REF] Pacejka | Tire and vehicle dynamics[END_REF]. In addition to geometrical interpretation, this effective tire radius model is also based on practical measurement data of the tire radius deformation during a travel. Therefore, this tire radius model describes the evolution of the effective tire radius in various situations and especially in the cases where the slip ratio is different from zero. Different formulations of the Pacejka's tire radius model are given in [START_REF] Pacejka | Tire and vehicle dynamics[END_REF]. One of them is

R rol = R ω - F nomin K zz F ref f F z F nomin + D ref f arctan B ref f F z F nomin , (3.23)
with

R ω = R 0 q reo + q V 1 R 0 ω V 0 2 . ( 3.24)
In Eq. (3.23) and (3.24), F ref f , B ref f , D ref f , q reo and q V 1 are empirical parameters fitted to measurement data. Pacejka also provides an expression of the quantity Fz F nomin as a function of physical quantities such as the tire forces F x , F y and the tire inflation pressure dp i 3.6. The different parameters q re0 , q V 1 , q V 2 , q F cx1 , q F cy1 , q F z1 , q F z2 and p F z1 are parameters allowing the user to set the sensitivity of the model with the different physical quantities such as the tire forces or the tire inflation pressure. As for F ref f , B ref f and D ref f , these parameters are empirical and require to be fitted to measurement data.

F z F nomin = 1 + q V 2 |ω| R 0 V 0 -q F cx1 F x F nomin 2 -q F cy1 F y F nomin 2 q F z1 ρ R 0 + q F z2 ρ 2 R 2 0 (1 + p F z1 dp i ) , (3.25) with ρ = max (R ω -R load , 0) , (3.26)
R load = R 0 - F z K zz , (3.27)
K zz0 = F nomin R 0 q F z1 + 2q F z2 ρ 0 R 0 , (3.28)
K zz = K zz0 (1 + p F z1 dp i ) . ( 3
At this point, the Pacejka effective tire radius model seems to be a very accurate model because it includes many physical quantities influencing the tire radius. However, this model requires to carry out experiments on the tire to set all the coefficients

F ref f , B ref f , D ref f , q re0 , q V 1 , q V 2 , q F cx1
, q F cy1 , q F z1 , q F z2 and p F z1 . Nevertheless, as mentioned by Pacejka [START_REF] Pacejka | Tire and vehicle dynamics[END_REF], this model can still be used if only the free tire radius R 0 , the vertical stiffness K zz and the four parameters

F ref f , D ref f et B ref f are known.
In this case, the parameters q V 1 , q V 2 , q F cx1 , q F cy1 , q F z2 , p F z1 should be set to zero and q re0 should be equal to one. Consequently, Equation (3.23) becomes

R rol = R 0 - F nomin K zz [F ref f ρ + D ref f arctan (B ref f ρ)] , (3.30)
with ρ = max (R 0 -R load , 0) . (3.31)

Summary of the vehicle model equations and assumptions

The vehicle model considered in this study is a single-track model extended with a simplified suspension model. The assumptions made to establish this model are • the roll dynamics are neglected,

• the left and right front steering angles are equal ⇒ δ f l = δ f r = δ,

• the rear steering angles are equal to zero ⇒ δ rl = δ rr = 0,

• variations of the CoG height neglected,

• equilibrium of the load transfer between the front and rear parts of the vehicle.

The model equations including the tire forces F xi , F zi and the slip ratio are reminded in Equation (3.32a)-(3.32j).

vx = 1 m F xf + F xr - 1 2 ρ a S a C x v 2 x -f RR mg , ωf = 1 2I wf [T f -R f load F xf ] , ωr = 1 2I wr [T r -R rload F xr ] , κ + ω 0 Q κ + ω 2 0 κ = h G I y m vx + 1 2 .ρ a S a C x v 2 x , F zf = mg L r L f + L r + (kκ + b κ) , F zr = mg L f L f + L r -(kκ + b κ) , s f = ω f R rol f -v x max(ω f R rol f , v x ) , s r = ω r R rolr -v x max(ω r R rolr , v x ) , R rol f = R 0 - F nomin K zz [F ref f ρ f + D ref f arctan (B ref f ρ f )] , R rolr = R 0 - F nomin K zz [F ref f ρ r + D ref f arctan (B ref f ρ r )] . (3.32a)
(3.32b)

(3.32c) (3.32d) (3.32e) (3.32f) (3.32g) (3.32h) (3.32i) (3.32j)

Description of the observer

State-space representation

In order to perform state estimation, it is necessary to put the used model in the form of a state-space representation

ẋ(t) = f (x(t), u(t), t) , z(t) = h (x(t), t) . (3.33)
where x ∈ R nx represent the state vector, u ∈ R nu the input vector and z ∈ R nz the vector of measurements.

The measurement vector z contains the quantities that can be measured with the sensors fitted on production vehicles. Thus, the measurement vector contains the longitudinal vehicle speed, the wheel speeds and the pitch speed, i.e.,

z(t) = [v x (t), ω f (t), ω r (t), κ(t)] T . (3.34)
The input vector u contains the inputs of the system. In the case of a vehicle car, the inputs are mainly the torques applied on the different parts of the vehicle. Therefore,

u(t) = [T f (t), T r (t)] T . (3.35)
The quantities included in the state vector x are the quantities that are estimated by the observer. Hence, x should at least include the signals we want to estimate or quantities related to them. In our case, we want to estimate the tire forces F xi , F zi and the slip ratios s i . In addition, in order to be able to put the model in a state-space representation form (see Eq.(3.33)), it is necessary to determine relation characterizing the derivative of the state as a function of the state x, the input u and the time.

In our case, the dynamical parts of F zf and F zr are characterized by κ and κ (see Eq. (3.32e) and (3.32f)). Therefore, instead of including the normal loads in the state, the choice is made here to include only the pitch dynamics κ and κ. This choice has several advantages. Firstly, in this case, it is not necessary to find a relation linking the derivative of the normal loads to the state. It is only necessary to find this type of relation for the pitch dynamics such as the one given in Eq (3.32d). Secondly, avoiding to insert the normal loads in the state has the advantage to keep a reasonable dimension of the state. As indicated in [START_REF] Fang | Nonlinear bayesian estimation: From kalman filtering to a broader horizon[END_REF], having a low dimensional state reduces the computation cost, especially in the case where the used observer is an extended Kalman filter.

Furthermore, as reminded in Eq. (3.32g) and (3.32h), the slip ratios depend on the longitudinal speed v x , the wheel speeds ω f and ω r and the normal loads F zf and F zr . Thus, the choice is made here to include v x , ω f and ω r in the state and then, to use their estimates with the normal load estimate in order to determine slip ratios estimates.

Therefore, it only remains to find a relation characterizing the derivative of the longitudinal tire forces. Different solutions exist in the literature [START_REF] Wielitzka | State and maximum friction coefficient estimation in vehicle dynamics using ukf[END_REF][START_REF] Ray | Nonlinear tire force estimation and road friction identification: Simulation and experiments[END_REF][START_REF] Rajamani | Tire-road friction-coefficient estimation[END_REF]. One solution is to introduce a tire force model expression extracted from a tire model such as the Pacejka tire model or the brush model [START_REF] Pacejka | Tire and vehicle dynamics[END_REF][START_REF] Dugoff | Tire performance characteristics affecting vehicle response to steering and braking control inputs[END_REF][START_REF] Canudas-De-Wit | A new model for control of systems with friction[END_REF][START_REF] Van Zanten | Measurement and simulation of transient tire forces[END_REF][START_REF] Szostak | Analtical modeling of driver response in crash avoidance maneuvering. volume 2: An interactive tire model for driver/vehicle simulation[END_REF]. However, this solution has the disadvantage to require a tire model and thus maybe others signals and parameters. Another solution introduced by Ray in [START_REF] Ray | Nonlinear tire force estimation and road friction identification: Simulation and experiments[END_REF] is to consider the longitudinal tire forces as a random walk, also called Gauss-Markov process [START_REF] Sri-Jayantha | Determination of nonlinear aerodynamic coefficients using the estimation-before-modeling method[END_REF]. In this case, the derivative of the longitudinal tire force is considered constant. Consequently, the longitudinal tire forces and its derivatives can be represented by Ḟxi Fxi = 0 1 0 0

F xi Ḟxi + w F x , i ∈ {f, r} (3.36)
where w F x is a zero mean white noise. This way to represent the longitudinal tire forces has the advantage to not require specific tire model for describing the tire forces. This advantage is valuable because, if the tire force expressions depend on parameters related to the tire, it is necessary to adapt the tire force expression for each tested tire. Furthermore, representing the tire forces with random walks was already used with success [START_REF] Ray | Nonlinear tire force estimation and road friction identification: Simulation and experiments[END_REF]. For these reasons, the longitudinal tire forces are here represented with the random walk describe in Equation (3.36). Thus, the selected state is given by

x(t) = [v x (t), ω f (t), ω r (t), F xf (t), F xr (t), Ḟxf (t), Ḟxr (t), κ(t), κ(t)] T . (3.37)
Once all the necessary relations have been established, the vehicle model can be put into the form of a state-space representation with

x(t) = [v x (t), ω f (t), ω r (t), F xf (t), F xr (t), Ḟxf (t), Ḟxr (t), κ(t), κ(t) T , u(t) = [T f (t), T r (t)] T , z(t) = [v x (t), ω f (t), ω r (t), κ(t)] T , (3.38a)
(3.38b) (3.38c)
and

ẋ(t) =                vx (t) ωf (t) ωr (t) Ḟxf (t) Ḟxr (t) Fxf (t) Fxf (t) κ(t) κ(t)                =                     1 m [F xf (t) + F xr (t) -1 2 ρ a S a C x v 2 x (t) -f RR mg] 1 2I wf [T f (t) -R f load (t)F xf (t)] 1 2Iwr [T r (t) -R rload (t)F xr (t)] Ḟxf (t) Ḟxr (t) 0 0 κ(t) h G Iy m vx (t) + 1 2 .ρ a S a C x v 2 x (t) -ω 0 Q κ(t) -ω 2 0 κ(t)                     , (3.39)
z(t) =     v x (t) ω f (t) ω r (t) κ(t)     =    
1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1

                   v x (t) ω f (t) ω r (t) F xf (t) F xr (t) Ḟxf (t) Ḟxr (t) κ(t) κ(t)                . (3.40)
After describing the state-space representation and the different models used hereafter, we can illustrate with blocs the different steps leading to the friction points estimates. These step are given in Figure 3.9. As indicated in Figure 3.9, vehicle inputs and measurements are combined with a vehicle model and an observer to estimate the state of the model. Among the state estimates, the pitch dynamics are combined with a suspension model to determine normal load estimates. This normal load estimate is on one side combined with the estimated longitudinal tire force to determine friction estimate. On the other side, the normal load is combined with an effective tire radius model and estimated speeds in order to determine slip ratios.

Kalman filters

In addition to a vehicle model, the second essential element necessary to perform state estimation is to define an observer structure. The choice of the observer structure relies on the linear or nonlinear behavior of the vehicle model and the existence of noise and disturbances in the model and measurements used. In our case, the measurements used are noisy measurements. Besides, the used vehicle model is an approximation of the real model describing the system behavior. For this reason, it is necessary to introduce a process noise in order to take into account of the differences between the approximate used model and the real model. Well-known state observers for such stochastic systems are the Kalman filters [START_REF] Kailath | Linear Estimation[END_REF][START_REF] Simon | Optimal state estimation: Kalman, H ∞ and nonlinear approaches[END_REF]. However, as described by Eq. (3.39), the model used in this work model is nonlinear with respect to the state variables. Since the classical Kalman filter is not designed to handle non-linear models [START_REF] Simon | Optimal state estimation: Kalman, H ∞ and nonlinear approaches[END_REF], an alternative designed for nonlinear models will be used : the Extended Kalman Filter (EKF) [START_REF] Simon | Optimal state estimation: Kalman, H ∞ and nonlinear approaches[END_REF]. In order to simplify the description of the EKF, let us first introduce the classic Kalman filter.

Kalman filter

As mentioned earlier, the Kalman filter is a stochastic estimator providing a state estimate of a model. This estimator can be formulated with a continuous-time or a discrete-time formulation [START_REF] Simon | Optimal state estimation: Kalman, H ∞ and nonlinear approaches[END_REF]. Given that the real system measurements are mostly discrete, the Kalman filter described here is a discrete one.

Let us consider a model described by a discrete-time linear time state-space representation

x k = F x k-1 + Gu k-1 + w k-1 , z k = Hx k + v k . (3.41a) (3.41b)
where w and v are white, zero-mean and uncorrelated noise. w k ∈ R nx and v k ∈ R nz are respectively called the process noise and the measurement noise and are characterized by their covariance matrices W k and V k . In addition, the covariance matrix of the state x k at time k is denoted P k .

The goal of the Kalman filter is to provide an estimate of the state x k and its covariance matrix P k by using the model description (Eq. (3.41a) and (3.41b)) and noisy measurements z k . In the Kalman filter description, this objective is reached in two steps. First, a prediction step where the model equations are used to predict the state value at time k. This estimate is here denoted xk|k-1 . Then, a correction step is performed on which the prediction xk|k-1 of x k is modified to take into account of the measurement vector at time k, z k . This estimate is here denoted xk|k . It is important to note that xk|k-1 and xk|k are both estimates of x k . The only difference between the two estimates is that in xk|k , we take into account of the measurement vector z k . An illustration of the Kalman filter is depicted in Figure 3.10.

During the prediction step, the estimated state and its covariance matrix are propagated with the model equations (Eq. (3.41a)). The

xk|k-1 = F xk-1|k-1 + Gu k-1 , P k|k-1 = F P k-1|k-1 F ⊤ + W k-1 . (3.42a) (3.42b)
Then, during the correction step, the measurements z k are taken into account in the state estimation which result according to [START_REF] Simon | Optimal state estimation: Kalman, H ∞ and nonlinear approaches[END_REF] in

K k = P k|k-1 H ⊤ HP k|k-1 H ⊤ + V k -1 , xk|k = xk|k-1 + K k z k -H xk|k-1 , P k|k = F P k|k-1 F ⊤ -F K k HF ⊤ + W k . (3.43a) (3.43b) (3.43c)
Thus, the Kalman filter provides estimate of the state xk and its covariances matrices P k . Here, it is applied with a linear model (See Eq. (3.41a) and (3.41a)). However, during this thesis, the vehicle model considered is nonlinear with respect to the state (See Eq. (3.39)). For this reason, the Kalman filter should be adpated to be employed with nonlinear model. Fortunately, different formulation of the Kalman filter exist to estimate the state of nonlinear system [START_REF] Simon | Optimal state estimation: Kalman, H ∞ and nonlinear approaches[END_REF][START_REF] Särkkä | Bayesian Filtering and Smoothing[END_REF][START_REF] Allik | Nonlinear estimators for censored data: A comparison of the ekf, the ukf and the tobit kalman filter[END_REF][START_REF] Wan | The unscented kalman filter for nonlinear estimation[END_REF]. One of them is the Extended Kalman filter. Remark 3.4. In Section 3.3.2, the considered state-space representations are time-invariant (see Eq. (3.41a) and (3.41b)). We decide to present the Kalman filter in a time-invariant framework because the used vehicle model is represented with a nonlinear time-invariant state-space representation (see Eq. (3.39) and (3.40)). However, the Kalman filter equations also applied with time-invariant state-space representation (see for example [START_REF] Simon | Optimal state estimation: Kalman, H ∞ and nonlinear approaches[END_REF]. In this case, the different matrices F, G and H are replace by time-varying matrices F k , G k and H k .

Extended Kalman filter

Since the Kalman filter equations (see Eq. (3.42a)-(3.43c)) involve linear time state-space representations only, the idea of the EKF is to linearise the basic model around the current state estimate, then to apply the Kalman filter equations to the linearised model. In a mathematical form, let us consider a discrete-time nonlinear model

x k = f (x k-1 , u k-1 ) + w k-1 , z k = h (x k ) + v k , (3.44a) (3.44b)
where v k ∈ R nz and w k ∈ R nx are respectively the measurement noise and process noise. By using a Taylor expansion around the current state estimate xk-1 , Equation (3.44a) becomes

xk = F x k-1 + ũk-1 + w k-1 , (3.45)
with ũk-1 = f ( xk-1 , u k-1 ) -F xk-1 , F = ∂f ∂x xk-1 . (3.46a) (3.46b)
Another Taylor expansion is performed around the state estimate xk|k-1 to linearise the Equation (3.44b),

ẑk = Hx k + zk + v k , (3.47)
with zk = h xk|k-1 -H xk|k-1 , H = ∂h ∂x xk|k-1 . (3.48a) (3.48b)
As a result, Eq. (3.45) and (3.47) form a discrete-time linear state space representation on which the Kalman filter can be applied. In this case, the extended Kalman filter equations become

P k|k-1 = F P k-1|k-1 F ⊤ + W k-1 xk|k-1 = f xk-1|k-1 , u k-1 , K k = P k|k-1 H ⊤ HP k|k-1 H ⊤ + V k -1 , xk|k = xk|k-1 + K k z k -h xk|k-1 , P k|k = (I n -K k H) P k|k-1 .
(3.49a)

(3.49b) (3.49c) (3.49d) (3.49e)
These equation can be used to estimate the state of discrete-time nonlinear model. However, in our case, the system studied is model with continuous-time equations (see Eq. (3.39)). Therefore, the EKF equations must be slightly modified in order to be applied to the continuoustime model.

Discrete-time EKF applied to continuous-time model

Let us consider the nonlinear continuous-time model

ẋ(t) = f c (x(t), u(t), t) + w(t), z(t) = h c (x(t), t) + v(t). (3.50)
In order to apply the discrete-time EKF, ẋ(t) can be discretized using an Euler explicit method with a sampling period dt,

x k = x k-1 + dt f c (x k-1 , u k-1 ) , x k = f (x k-1 , u k-1 ) . (3.51) with f (x k-1 , u k-1 ) = x k-1 + dt f c (x k-1 , u k-1
). In these conditions, the Jacobian F is given by

F = ∂f ∂x xk-1 = ∂ ∂x [x k-1 + dt f c (x k-1 , u k-1 )] xk-1 = I nx + dt F c ( xk-1 , u k-1 ) (3.52)
with

F c = ∂fc ∂x xk-1
.

Consequently, the Kalman filter equations became

xk|k-1 = xk-1 + dt f c ( xk-1 , u k-1 ) , P k|k-1 = F P k-1|k-1 F ⊤ + W k-1 , K k = P k|k-1 H ⊤ HP k|k-1 H ⊤ + V k -1 , xk|k = xk|k-1 + K k z k -h c xk|k-1 , P k|k = (I n -K k H) P k|k-1 , (3.53a)

Conclusion

In this chapter, the methodology for state observer estimation have been described in order to estimate the friction data points. This method requires to introduce a vehicle model and an observer structure. The considered vehicle model is a single-track model extended with a suspension model in order to take into account of the load transfer and thus, the variation of normal load applied on a tire during a travel. Because the chosen vehicle model is an approximation of the real system, process noise is added on the vehicle model which results in a stochastic model. Consequently the observer structure chosen is a Kalman filter since it is design to handle such model. Finally, because the used vehicle model is nonlinear, the selected Kalman filter is an extended Kalman filter.

Until now, the question of the Kalman filter settings has not been addressed. However, the Kalman filter requires the setting of four different quantities. In details, an initial state vector x 0 , its corresponding covariance matrix P 0 and more importantly, values for W and V . While prior information can be used to set the initial state value of the Kalman filter, setting the W and V matrices can be troublesome. For this reason, the following chapter is dedicated to the setting of the noise covariance matrices.

Introduction

Kalman filtering [START_REF] Stengel | Optimal Control and Estimation[END_REF][START_REF] Kailath | Linear Estimation[END_REF][START_REF] Simon | Optimal state estimation: Kalman, H ∞ and nonlinear approaches[END_REF] is probably the first solution any engineer would suggest when dynamical system state estimation problems come into play. This matter of fact can be illustrated by the myriads of Kalman filtering solutions commonly used, e.g., for guidance, navigation or control of any kind of vehicles [START_REF] Bar-Shalom | Estimation with applications to tracking and navigation: theory, algorithms and software[END_REF][START_REF] Zarchan | Fundamentals of Kalman filtering: a practical approach[END_REF]. One standard practical case is the use of Kalman filters for vehicle state estimation in Advanced Driver Assistance Systems [START_REF] Doumiati | Vehicle Dynamics Estimation using Kalman Filtering[END_REF] where signals like the rear and front axle side forces cannot be measured directly but must be reconstructed from measured signals like the yaw rate or the body slip angle at the center of gravity, respectively. The Kalman filter success for state estimation is probably due to (i) its easiness of implementation [START_REF] Verhaegen | Filtering and system identification: a least squares approach[END_REF], (ii) its capability to give access to the minimum variance state estimate under mild practical conditions [START_REF] Kailath | Linear Estimation[END_REF][START_REF] Simon | Optimal state estimation: Kalman, H ∞ and nonlinear approaches[END_REF], (iii) its recursive update-predict structure which makes it still efficient even if the aforementioned practical conditions [START_REF] Simon | Optimal state estimation: Kalman, H ∞ and nonlinear approaches[END_REF] are not perfectly satisfied.

As shown in Section 4.2, the standard Kalman filter equations involve a deterministic stochastic state space model [START_REF] Van Overschee | Subspace identification for linear systems. Theory, implementation, applications[END_REF] of the system to track. The stochastic part is usually characterized via additive output and process noises. In Kalman filtering theory, this stochastic contribution not only accounts for measurement noise effects but also for dynamics not embedded into the deterministic model components. These output and process noises, more specifically their mean values and their covariance matrices, thus play a central role in the Kalman filter efficiency because

• they are (with the deterministic state space matrices) the main inputs of the standard Kalman filter algorithms [START_REF] Stengel | Optimal Control and Estimation[END_REF][START_REF] Simon | Optimal state estimation: Kalman, H ∞ and nonlinear approaches[END_REF],

• they quantify the confidence the user has into the deterministic model description of the system to track,

• their mispecification may lead to divergence issues [START_REF] Simon | Optimal state estimation: Kalman, H ∞ and nonlinear approaches[END_REF].

The dependence of the Kalman filter performance on the availability of reliable output and process covariance matrices may explain the reasons why a significant part of the literature dedicated to Kalman filter tuning is devoted to the development of techniques and algorithms for the estimation of these covariance matrices. As recently described in the survey [START_REF] Duník | Noise covariance matrices in state-space models: A survey and comparison of estimation methods-part i[END_REF] of Dunik et al., most of the noise covariance matrix estimation methods can be classified into two main families:

• the "feedback methods" where the noise covariance matrices and the unknown state signals are estimated in one shot,

• the "feedback free methods" where the state variables are first estimated (e.g., with a non optimal estimator), then used for noise characteristics estimation.

In this chapter, in order to avoid nonlinear optimization issues [START_REF] Bianchi | Process noise covariance estimation via stochastic approximation[END_REF] usually encountered with the "feedback methods", a "feedback free method" using linear least squares only is suggested. More specifically, instead of resorting to a suboptimal filter as required, e.g., by the well known innovation correlation methods [START_REF] Mehra | On the identification of variances and adaptive kalman filtering[END_REF][START_REF] Mehra | Approaches to adaptive filtering[END_REF] developed by R. Mehra in the 70s, the method introduced in this chapter relies on the comparison of an estimated model of the deterministic system dynamics (determined from the available input-output data sets) and the discrete time state space representation involved in the Kalman filter equations. By doing so, contrary to most of the "feedback free methods", no specific tuning parameter is required. By assuming that the selected model identification method is reliable enough to yield accurate and consistent discrete time linear time invariant state space models under the assumptions required by the Kalman filter, the basic idea of the solution developed herein consists in translating the discrepancy between the estimated model and the dynamical state space model used by the Kalman filter into noise covariance matrix estimates. As we said previously, in Kalman filtering theory, even if the system to track is not necessarily truly stochastic, the stochasticity brought by the output and process noises is used to describe the model uncertainties as closely as possible. Accurately determining the modeling error induced by the deterministic state space model of the Kalman filter should be a smart way to quantify the model uncertainties and, by extension, the noise covariance matrices.

As shortly explained in the former paragraph, the noise covariance matrix estimation technique introduced in this chapter includes a data based model learning step. The system identification literature [START_REF] Söderström | System identification[END_REF][START_REF] Johansson | System modeling and identification[END_REF][START_REF] Ljung | System identification. Theory for the user[END_REF] is full of algorithms and model learning techniques dedicated to discrete time state space model estimation. This is all the more true when linear time invariant models are sought. Herein, a class of subspace model identification methods [START_REF] Van Overschee | Subspace identification for linear systems. Theory, implementation, applications[END_REF][START_REF] Katayama | Subspace methods for system identification[END_REF][START_REF] Verhaegen | Filtering and system identification: a least squares approach[END_REF][START_REF] Mercere | Regression techniques for subspace-based black-box state-space system identification: an overview[END_REF] is adapted to estimate a discrete time linear time invariant state space model used afterwards for noise covariance matrix estimation. The main reason why this class of methods is selected is the strong and inherent link between the N4SID or CCA class [START_REF] Katayama | Subspace methods for system identification[END_REF] of subspace identification methods and the Kalman filter [START_REF] Van Overschee | N4SID: subspace algorithms for the identification of combined deterministic stochastic systems[END_REF][START_REF] Van Overschee | Subspace identification for linear systems. Theory, implementation, applications[END_REF]. Indeed, as shown first by P. van Overschee and B. De Moor [START_REF] Van Overschee | N4SID: subspace algorithms for the identification of combined deterministic stochastic systems[END_REF], N4SID-like methods "introduce a bank of non steady state Kalman filters generating a sequence of state estimates [START_REF] Van Overschee | Subspace identification for linear systems. Theory, implementation, applications[END_REF]". As far as model error modeling is concerned, the main drawback of the standard subspace model identification methods is their ability to give access to state space models up to a similarity transformation only. As shown in Section 4.3, such a weak point of subspace model identification can be bypassed by resorting to prior available when Kalman filtering is the final goal of the estimation procedure.

The chapter is organized as follows. Section 4.2 describes the problem tackled in this chapter, then introduces the main notations used in the next paragraphs. Section 4.3 is devoted to the main contribution of this chapter, i.e., (i) the detailed description of the algorithmic steps leading to accurate estimates of the noise covariance matrices necessary for Kalman filtering, (ii) the proof of the asymptotic unbiasedness of these estimates. The validation of this approach is carried out in Sec-tion 4.4 where three different simulation scenarii are considered. Section 4.5 concludes this chapter.

Problem formulation and notations

In this chapter, the goal is to estimate the state vector trajectory of a dynamical system, the behavior of which can be described by the following linear time invariant continuous time state space representation

ẋ(t) = Ax(t) + Bu(t), (4.1a)
y(t) = Cx(t), (4.1b)
where x(t) ∈ R nx stands for the state vector, y(t) ∈ R ny is the output vector and u(t) ∈ R nu is the input vector, respectively. In order to reach this goal, we assume that a Kalman filter can be designed. More specifically, it is assumed that the state at time kT s (k ∈ N and T s being the sampling period) can be determined from the current and previous measurements by using the famous update-predict equations [START_REF] Simon | Optimal state estimation: Kalman, H ∞ and nonlinear approaches[END_REF], e.g., the following set of equations (starting with x+ 0 = E[x 0 ] and

P + 0 = E[(x 0 -x+ 0 )(x 0 -x+ 0 ) ⊤ ]) x- k = A d x+ k-1 + B d u k-1 , (4.2a)
P - k = A d P + k-1 A ⊤ d + W , (4.2b
)

K k = P - k C ⊤ d + S C d P - k C ⊤ d + V + C d S + S ⊤ C ⊤ d -1 , (4.2c
)

x+ k = x- k + K k (y k -C d x- k ), (4.2d)
P + k = P - k -K k C d P - k + S ⊤ , (4.2e)
where the notations -and + stand for "before and after new measurements". These Kalman filter equations are valid and lead to an optimal [START_REF] Stengel | Optimal Control and Estimation[END_REF][START_REF] Kailath | Linear Estimation[END_REF][START_REF] Simon | Optimal state estimation: Kalman, H ∞ and nonlinear approaches[END_REF] state estimator when the tracked system is described by the following discrete time linear time invariant state space representation

x k+1 = A d x k + B d u k + w k , (4.3a
)

y k = C d x k + v k , (4.3b)
where the output and process noises v k ∈ R ny and w k ∈ R nx , respectively, are zero mean random variables with covariance matrices satisfying

E v i w i v ⊤ j w ⊤ j = V S S ⊤ W δ ij , (4.4)
where δ ij is the Kronecker delta function. In addition to initial values for x + 0 and X + 0 , the Kalman equations require the knowledge of the state space matrices A d , B d and C d as well as the covariance matrices S, V and W . While the matrices A d , B d and C d can be generated from the matrices A, B and C given in Eq. (4.1) as follows [START_REF] Stengel | Optimal Control and Estimation[END_REF] A d = e ATs , (4.5a)

B d = e ATs
Ts 0 e -Aτ Bdτ, (4.5b)

C d = C, (4.5c)
by assuming that the input signals are constant in the interval of integration, the matrices S, V and W cannot be directly extracted from prior knowledge. The problem of estimating these matrices is thus now tackled by resorting to a specific model learning technique inspired by the famous subspace model identification methods [START_REF] Van Overschee | Subspace identification for linear systems. Theory, implementation, applications[END_REF][START_REF] Katayama | Subspace methods for system identification[END_REF][START_REF] Verhaegen | Filtering and system identification: a least squares approach[END_REF]. Before describing our new solution for noise covariance matrix estimation into details, let us introduce the following notations. For any vector r k ∈ R nr and parameters M , i and ℓ ∈ N + * , we define

r i,M =      r i r i+1 . . . r i+M -1      ∈ R M nr , (4.6)
R i,M = r i r i+1 • • • r i+M -1 ∈ R nr×M , (4.7)
and the block Hankel matrix as follows

R i,ℓ,M =      r i r i+1 • • • r i+M -1 r i+1 r i+2 • • • r i+M . . . . . . . . . . . . r i+ℓ-1 r i+ℓ • • • r i+M +ℓ-2      ∈ R ℓnr×M . ( 4

.8)

With matrices A, B, C and D of appropriate dimensions, for ℓ ≥ n x , the extended controllability matrix can be defined as follows

Ω ℓ (A, B) = A ℓ-1 B • • • AB B . (4.9)
We also define the extended observability matrix

Γ ℓ (A, C) =      C CA . . . CA ℓ-1      , (4.10)
and the block lower triangular Toeplitz matrix

∆ ℓ (A, B, C, D) =      D 0 • • • 0 CB D • • • 0 . . . . . . . . . . . . CA ℓ-2 B • • • CB D      . (4.11)

Noise covariance matrix estimation with a modified subspace model identification approach

As we said previously, in the Kalman filter equations, the stochastic part of the state space model (4.3) is introduced mainly to describe the uncertainties in the model dynamics, i.e., the confidence we have in the values of the matrices A d , B d and C d to characterize the system dynamics. Said differently, the matrices S, V and W , which characterize the randomness of the zero mean noises (v k ) k∈Z and (w k ) k∈Z , must be chosen to quantify the confidence we have in the deterministic part of the state space representation (4.3). In order to reach this goal, it is necessary to determine if the matrices A d , B d and C d chosen a priori by the user describe the real system dynamics accurately. By assuming the access to a sufficiently rich data set, i.e., by assuming that the samples {u k , y k } N k=1 used by the Kalman filter embed the system dynamics of interest, the solution considered in this chapter consists in (i) using system identification for extracting an accurate discrete time model from the available data sets, (ii) comparing the estimated state space with the deterministic part of the state space form (4.3) we know a priori in order to get reliable model misfit measurements, (iii) transforming this discrepancy measurements into covariance matrix estimates.

Among all of the state space model identification techniques available in the literature, a specific attention is paid herein to the subspace model identification method class because of its capability to give access to state variable estimates with linear least squares only.

Remark 4.1. In model learning, when subspace model identification comes into play, the user must choose an upper bound of the system order a priori. Because, in the Kalman filtering framework, the state space matrices A d , B d and C d are chosen a priori, the system order is in the list of priors. Thus, the technique developed in this chapter does not require any user action all along the estimation procedure.

The starting point of this model learning approach is the innovation state space model [START_REF] Kailath | Linear Estimation[END_REF] (which is equivalent to the Kalman filter given in Eq. (4.2) as proved, e.g., by [START_REF] Kailath | Linear Estimation[END_REF] T. Kailath et al.)

xk+1 = A d xk + B d u k + Ke k , (4.12a)
y k = C d xk + e k , (4.12b)
where K is the steady state Kalman gain while e k = y k -C d xk , called the innovation term or innovations [START_REF] Kailath | Linear Estimation[END_REF], stands for the part of the measurements which contains new information [START_REF] Kailath | Linear Estimation[END_REF]. Then, by using the notations introduced in Section 4.2, for N , f and ℓ ∈ N + * , straightforward recursions lead to the following data equation [START_REF] Mercere | Regression techniques for subspace-based black-box state-space system identification: an overview[END_REF][START_REF] Mercère | Innovation-based subspace identification in open-and closed-loop[END_REF]

Y f,ℓ,N = Γ f (A d , C d ) Xf,N + ∆ u f U f,ℓ,N + N f,ℓ,N , (4.13)
where

∆ u f = ∆ f (A d , B d , C d , 0), (4.14a
)

N f,ℓ,N = ∆ f (A d , K, C d , I ny ) ∆ e f E f,ℓ,N . (4.14b)
By explicitly using the definition of the innovation term into Eq. (4.12), with Ãd = A d -KC d , the predictor state space form [START_REF] Qin | An overview of subspace identification[END_REF]

xk+1 = Ãd xk + B d u k + Ky k , (4.15a
)

y k = C d xk + e k , (4.15b
xk = Ω p ( Ãd , K)y k-p,p + Ω p ( Ãd , B d )u k-p,p , (4.17)
can be viewed as the the optimal linear estimate of xk (in the mean square error sense [START_REF] Jansson | A linear regression approach to state space subspace system identification[END_REF]) given y k-p,p and u k-p,p . Of course, this approximation becomes better with increasing p. By using this state approximation, the quantity Xf,N involved in Eq. (4.13) becomes

Xf,N = xf xf+1 • • • xf+N-1 ≃ xf xf+1 • • • xf+N-1 . (4.18)
Furthermore, by using Eq. (4.17), we have

xf xf+1 • • • xf+N-1 = Ω p ( Ãd , B d ) Ω p ( Ãd , K) u f -p,p u f -p+1,p • • • u f +N -1-p,p y f -p,p y f -p+1,p • • • y f +N -1-p,p . (4.19) Thus, by taking f = p = ℓ, we get Xf,N ≃ Ω p ( Ãd , B d ) Ω p ( Ãd , K) u 0,ℓ u 1,ℓ • • • u N -1,ℓ y 0,ℓ y 1,ℓ • • • y N -1,ℓ ≃ Ω p ( Ãd , B d ) Ω p ( Ãd , K) U 0,ℓ,N Y 0,ℓ,N . (4.20)
Substituting Eq. (4.20) into Eq. (4.13) leads to

Y f,ℓ,N =Γ f (A d , C d ) Ω p ( Ãd , B d ) Ω p ( Ãd , K) U 0,ℓ,N Y 0,ℓ,N + ∆ u f U f,ℓ,N + N f,ℓ,N . (4.21)
In this equation, the unknown quantities are

Γ f (A d , C d )Ω p ( Ãd , B d ), Γ f (A d , C d )Ω p ( Ãd , K), ∆ u
f and, of course, the term representing the noise N f,ℓ,N . One of the main interests of this data equation is its linearity with respect to the unknown matrices

Γ f (A d , C d )Ω p ( Ãd , B d ), Γ f (A d , C d )Ω p ( Ãd , K) and ∆ u
f . In this chapter, a specific attention is paid to the first two block matrices ,K) is equivalent to solve a linear least-squares problem. This one can be solve by resorting to the following QR factorization [START_REF] Meyer | Matrix analysis and applied linear algebra[END_REF]

Γ f (A d , C d )Ω p ( Ãd , B d ) and Γ f (A d , C d )Ω p ( Ãd , K), respectively. As represented in Eq. (4.21), estimate Γ f (A d , C d )Ω p ( Ãd , B d ) and Γ f (A d , C d )Ω p ( Ãd
    U f,ℓ,N U 0,ℓ,N Y 0,ℓ,N Y f,ℓ,N     =   L 11 0 0 L 21 L 22 0 L 31 L 32 L 33     Q 1 Q 2 Q 3   . (4.22)
Indeed, it can be proved [START_REF] Peternell | Statistical analyss of novel subspace identification methods[END_REF] that, if

• the innovation sequence (e k ) k∈Z is a zero mean white noise sequence with covariance matrix R e ,

• the pair (A d , C d ) is observable [START_REF] Kailath | Linear Estimation[END_REF] and the pair A d , B d , KR

1/2 e is reachable [START_REF] Kailath | Linear Estimation[END_REF],

• the input signals u are quasi stationary and exciting of sufficient order [START_REF] Ljung | System identification. Theory for the user[END_REF],

we have

lim N →∞ L 32 L -1 22 = Γ f (A d , C d ) Ω p ( Ãd , B d ) Ω p ( Ãd , K) . (4.23)
By recalling that,

Ω p ( Ãd , B d ) Ω p ( Ãd , K) U 0,ℓ,N Y 0,ℓ,N = Xf,N , (4.24)
Eq. (4.23) leads to

lim N →∞ L 32 L -1 22 U 0,ℓ,N Y 0,ℓ,N = Γ f (A d , C d ) Xf,N . (4.25)
This last equation shows that the SVD [59]

L 32 L -1 22 U 0,ℓ,N Y 0,ℓ,N = U ΣV ⊤ (4.26)
can be used to approximate the column space of Γ f (A d , C d ) and the row space of the state sequence Xf,N as follows

Γf (A d , C d ) = U Σ 1/2 , (4.27a) Xf,N = Σ 1/2 V ⊤ . (4.27b)
The main issue with this SVD is the difficulty to guarantee that the row space of the state sequence Xf,N is generated in the system state basis, the decomposition of U ΣV ⊤ into two blocks being valid up to a similarity transformation [START_REF] Meyer | Matrix analysis and applied linear algebra[END_REF] T . Therefore, it is necessary to determine an estimate T of the similarity transformation in order to move the estimated state sequence Xf,N into the system state basis. This is the place where our prior on A

Γ f (A d , C d )T = Γf (A d , C d ). (4.28)
Once T is estimated with a Moore Penrose pseudo inverse [START_REF] Meyer | Matrix analysis and applied linear algebra[END_REF], the state sequence Xf,N can be moved into the "good" state basis, i.e., the system state basis, as follows

Xf,N = T Σ 1/2 V ⊤ . (4.29)
Instead of using this state sequence estimate for determining Âd , Bd and Ĉd as usually performed with the N4SID-like algorithms [START_REF] Van Overschee | Subspace identification for linear systems. Theory, implementation, applications[END_REF][START_REF] Katayama | Subspace methods for system identification[END_REF][START_REF] Verhaegen | Filtering and system identification: a least squares approach[END_REF], the prior on A d and C d is used herein to quantify the discrepancy between this prior information and the information brought by the data sets via the estimated state sequence Xf,N . More specifically, in order to determine accurate estimates of the covariance matrices S, V and W , respectively, knowing A d and C d and having access to an estimate of the state sequence Xf,N in the system state basis, the sought discrepancy measurement can be generated by computing the residuals [99]

Ŵf,1,N-1 Vf,1,N-1 = Xf+1,N Y f,1,N -1 - A d B d C d 0 Xf,N-1 U f,1,N -1 . (4.30)
These residuals can indeed be used to estimate S, V and W as follows [99]

V Ŝ Ŝ⊤ Ŵ = lim N →∞ 1 N Ŵf,1,N-1 Vf,1,N-1 Ŵ⊤ f,1,N -1 V⊤ f,1,N -1 . (4.31)
On top of that, it can be proved that these estimates are asymptotically unbiased estimates of V , W and S, respectively.

Theorem 4.1. By assuming that

• the system to be identified is linear time invariant,

• the deterministic input (u k ) k∈{1,...,N -1} is uncorrelated with the process and measurement noises (w k ) k∈Z and (v k ) k∈Z , respectively,

• the input (u k ) k∈{1,...,N -1} is persistently exciting [97, Definition 5] of order 2f ,

• the number of measurements N goes to infinity,

• the process and measurement noises are not identically zero, we have

lim N →∞ 1 N Ŵf,1,N-1 Vf,1,N-1 Ŵ⊤ f,1,N -1 V⊤ f,1,N -1 = V S S ⊤ W , (4.32)
i.e., the estimated covariance matrices V , Ŵ and Ŝ are asymptotically unbiased.

Proof. By referring to Section 4.2.1 of the book [START_REF] Van Overschee | Subspace identification for linear systems. Theory, implementation, applications[END_REF] entitled "Subspace Identification for Linear Systems-Theory, Implementation, Applications", the estimated state sequence Xf,N as well as the estimated observability matrix Γf (A d , C d ) generated from the SVD given in Eq. (4.26) are both asymptotically unbiased under the aforementioned assumptions. Thus, the similarity transformation relating Γf (A d , C d ) and Γ f (A d , C d ) can be determined accurately when number of measurements goes to infinity. Once T is available, the generation of the state sequence in the good state basis is direct thanks to Eq. (4.29). The generation of the residuals Ŵf,1,N-1 and Vf,1,N-1 follows directly and, thanks to the strong law of large numbers [START_REF] Papoulis | Probability, random variables, and stochastic processes[END_REF], the time average covariance matrix generated from the residuals (as in Eq. (4.31)) converges to the ensemble covariance matrix with probability one as long as N tends to infinity (see also Section 10.6 of the book [START_REF] Katayama | Subspace methods for system identification[END_REF] entitled "Subspace Methods for System Identification").

As clearly shown with the former equations,

• this covariance matrix estimation procedure involves linear least squares and robust linear algebra tools only,

• neither specific user action nor supplementary prior is required,

• the modeling error due to the initial choice of A d , B d and C d is first determined by comparing, in an indirect way, the state trajectory of an estimated model (by trusting the available data sets only) and the state sequence of the innovation form, then used for estimating S, V and W by generating residuals quantifying the modeling uncertainties.

It is now time to test the efficiency of this procedure by considering different numerical simulations.

Numerical validation 4.4.1 A toy example

In order to test the capabilities of our noise covariance matrix estimation solution under ideal conditions, let us first consider the following toy example used, e.g., in the book [START_REF] Van Overschee | Subspace identification for linear systems. Theory, implementation, applications[END_REF] written by P. van Overschee and B. de Moor in 1996 on subspace model identification. This discrete time linear time invariant system written in an innovation form (see Eq.(4.12)) is characterized by the following matrices The data length used in this simulation is N = 1000. The input signals are independent and zero mean Gaussian noise signals of unit variances while the innovation sequences are zero mean Gaussian noises with covariance matrix R e . Thus, the theoretical noise covariance matrices are V , Ŵ and Ŝ components, respectively. For comparison, the same parameters have been estimated (i) by computing the sample covariance matrices using the 1000 samples of each noise realization, (ii) with the ICM [START_REF] Mehra | On the identification of variances and adaptive kalman filtering[END_REF] for which the tuning parameters are chosen equal to 4 and K (see (4.33d)), (iii) with the DCM [START_REF] Odelson | A new autocovariance least-squares method for estimating noise covariances[END_REF] for which the tuning parameter is chosen equal to 5, (iv) with the CMM [START_REF] Myers | Adaptive sequential estimation with unknown noise statistics[END_REF] for which the initial matrices for V , W and P are chosen to be equal to identity matrices of correct dimensions, respectively. Theo. value 0.0183 -0.0045 0.0131 -0.0172 -0.0311 0.0063 -0.0234 0.0314 Sample cov. avg. 0.0183 -0.0045 0.0131 -0.0172 -0.0311 0.0063 -0.0234 0.0314 std. 0.0008 0.0002 0.0006 0.0008 0.0014 0.0003 0.0010 0.0014 New meth. avg. 0.0181 -0.0039 0.0137 -0.0169 -0.0302 0.0061 -0.0239 0.0309 std. 0.0011 0.0007 0.0010 0.0010 0.0018 0.0011 0.0015 0.0016

A d = 0.
V = R e = 0.0176 -0.0267 -0.0267 0.0497 , (4.34a)
Even if the estimated values yielded by our new method are slightly biased, the figures available in Table 4.1, 4.2 and 4.3, respectively, clearly show that (i) this new technique significantly outperforms methods available in the literature [START_REF] Mehra | On the identification of variances and adaptive kalman filtering[END_REF][START_REF] Myers | Adaptive sequential estimation with unknown noise statistics[END_REF][START_REF] Odelson | A new autocovariance least-squares method for estimating noise covariances[END_REF], (ii) the variances of the estimates obtained with this method are comparable to the variances computed from noise samples directly, i.e., of the same order of magnitude as what we could expect if we had access to the noise realizations. These values demonstrate that, under ideal practical conditions, our noise covariance matrix estimation technique is efficient.

A mass spring damper system

Let us now tackle the problem of state estimation. More specifically, let us test the efficiency of our new technique for tracking the state signals (position and velocity) of a mass spring damper system excited by a pseudo random binary sequence [START_REF] Ljung | System identification. Theory for the user[END_REF]. The dynamics of such a system can be described by a second order constant coefficient differential equation which can be compactly written with a continuous time linear time invariant state space representation like Eq. (4.1) with

A = 0 1 -k m -b m , B = 0 g m , C = 1 0 , (4.35)
where m is the mass of the moving box, k is the spring constant, b is the damping constant and g is the position sensor gain. Herein, we select m = 100 kg, g = 2, k = 35 kN/m and b = 25 kN.s/m. The input force is a pseudo random binary sequence tuned as suggested by J. Chen and C. Yu [START_REF] Chen | Optimal input design using generalized binary sequence[END_REF] (see Figure 4.1 for a realization of this input signal). The noise free output signal generated by simulating this continuous time linear time invariant system is perturbed, in a second step, by output measurement disturbances built from a zero mean white Gaussian noise with a variance selected so that the signal to noise ratio defined as follows (y nf standing for the noise free output signal)

SNR = 20 log σ 2 y nf σ 2 v , ( 4

Passenger car lateral dynamics

The performance of the noise covariance matrix estimation procedure is finally tested with simulated data sets used in the context of pas- senger vehicle guidance control. More specifically, we tackle the problem of estimating the state variables involved in the description of the lateral dynamics of passenger vehicles. As we said in Section 4.1 and the aforementioned references [START_REF] Bar-Shalom | Estimation with applications to tracking and navigation: theory, algorithms and software[END_REF][START_REF] Zarchan | Fundamentals of Kalman filtering: a practical approach[END_REF][START_REF] Doumiati | Vehicle Dynamics Estimation using Kalman Filtering[END_REF], this important task used, e.g., for controlling autonomous cars, can be performed, in the first place, by resorting to a Kalman filter. This third simulation example is introduced in this chapter to test the efficiency of our estimation solution when unmodeled dynamics start playing a central role in the Kalman filter tuning step.

In the literature [START_REF] Rajamani | Vehicle Dynamics and Control[END_REF][START_REF] Guiggiani | The Science of Vehicle Dynamics[END_REF], several models have been suggested to describe the lateral dynamics of passenger cars. Herein, in order to make simulations easier without losing dynamics modeling capabilities, a single track model is used, i.e., a model for which (i) the left and right steering angles are assumed to be equal, (ii) the rear steering angle is assumed to be null, (iii) the pitch and roll effects are neglected. As shown in Figure 4.6, such a vehicle model has only one front and one rear wheel, thus is often called the bicycle model [START_REF] Rajamani | Vehicle Dynamics and Control[END_REF][START_REF] Guiggiani | The Science of Vehicle Dynamics[END_REF]. Table 4.4 gathers the different signals involved in this bicycle model. Because we focus herein on the lateral dynamics only, the longitudinal vehicle speed can be assumed to be constant, i.e., vx (t) = 0, ∀ t ≥ 0. Thanks to standard laws of mechanics and trigonometric rules, the following set of equations can be used to describe the bicycle dynamical behavior (see Table 4.5 for the definition of the parameters involved in these equations)

ma y (t) = f r (t) + f f (t) cos(δ(t)), (4.37a)
a y (t) = vy (t) + v x ψ(t), (4.37b) The next step for the construction of a model of the lateral vehicle dynamics consists in modeling the lateral tire forces f f and f r which act on the vehicle. Many tire models have been introduced in the literature [START_REF] Khaleghian | A technical survey on tire-road friction estimation[END_REF]. Most of them try to describe the link between the lateral tire forces f i and the side slip angles α i , i ∈ {f, r}, via curve fitting techniques, i.e., by suggesting physical models or black box equations mimicking experimental data sets as those available in Figure 4.7. Under normal usage conditions, i.e., for a y < 0.5 -0.6g, the tire is used in its linear part as shown in Figure 4.7. Thus, a standard assumption consists in modeling the efforts by the following linear function of the side slip angle [START_REF] Doumiati | Vehicle Dynamics Estimation using Kalman Filtering[END_REF]

I z ψ(t) = ℓ f f f (t) -ℓ r f r (t), (4.37c) tan(β(t)) = v y (t) v x , (4.37d) α f (t) = δ(t) -arctan v y (t) + ℓ f ψ(t) v x , (4.37e) α r (t) = arctan ℓ r ψ(t) -v y (t) v x . ( 4
f f = c y f α f , f r = c yr α r , (4.38)
where c y i stands for the cornering stiffness, i ∈ {f, r}. These lateral force models are easy to use, describe real measurements efficiently when small angle values are considered but are, by construction, quasistatic only. In reality, tire forces are not generated instantaneously when the side slip angle changes due to the smoothness of the rubber. Thus, tire models involving such transient behaviors should be introduced. As suggested, e.g., by M. Doumiati et al. [START_REF] Doumiati | Vehicle Dynamics Estimation using Kalman Filtering[END_REF], a first order model can be used, i.e.,

τ f ḟf (t) + f f (t) = c y f α f (t), (4.39)
τ r ḟr (t) + f r (t) = c yr α r (t), (4.40)
where τ i stands for a relaxation time constant, i ∈ {f, r}. This time constant can be approximated as follows

τ i = c y i k i v x = σ α i v x , i ∈ {f, r}, (4.41)
where k i , i ∈ {f, r}, stands for an equivalent tire lateral stiffness while σ α i , i ∈ {f, r}, is a relaxation length which approximates the distance needed to build up tire forces. By selecting

x(t) =     β(t) ψ(t) f y f (t) f yr (t)     , u(t) = δ swa (t), y(t) =   β(t) ψ(t) a y (t)   , (4.42)
the former set of equations can be written in a compact way, leading to the following nonlinear state space representation

ẋ(t) = f (x(t), u(t)), (4.43a
)

y(t) = g(x(t), u(t)), (4.43b)
where f (•) and g(•) stand for the nonlinear equations linking the signals ẋ(t), x(t), y(t) and u(t), respectively. By assuming that [START_REF] Jazar | Vehicle dynamics: theory and application[END_REF] the data sets are generated so that the experiments are carried out in the linear range (up to 5 m/s 2 of lateral acceleration), the former nonlinear state space form can be linearized, leading to the following continuous time linear time invariant state space representation

ẋ(t) = Ax(t) + Bu(t), (4.44a
)

y(t) = Cx(t), (4.44b)
with

x(t) =     β(t) ψ(t) f y f (t) f yr (t)     , A =       0 -1 1 mvx 1 mvx 0 0 ℓ f Izz -ℓr Izz -c y f vx σα f -c y f ℓ 1 σα f -vx σα f 0 -c yr vx σα r c yr ℓ 2 σα r 0 -vx σα r       , (4.45a) v(t) = δ swa (t), B ⊤ = 0 0 vxcy f nσα f 0 , (4.45b)
y(t) =   β(t) ψ(t) a y (t)   , C =   1 0 0 0 0 1 0 0 0 0 1 m 1 m   , (4.45c)
where, again, the signals δ swa (t), β(t), ψ(t), f y f (t), f yr (t), a y (t) and the parameters m, v x , ℓ f , ℓ r , I zz , c y f , c yr , σ α f , σ αr , n are explicitly defined in Table 4.4 and Table 4.5, respectively. Both linear and nonlinear state space representations are going to be used in the sequels. While the linear model given in Eq. (4.44) is used, after discretization, to form the matrices A d , B d and C d used by the Kalman filter, the nonlinear model given in Eq. (4.43) is used to generate the simulation data sets.

By considering the car parameters available in Table 4.5 and by selecting T s = 0.01 s, let us first generate a noise free data set by exciting the nonlinear model given in Eq. (4.43) with the input signal given in Figure 4.8 assuming zero initial conditions. Remark 4.2. In order to reproduce practical conditions encountered when real cars are driven, the input signal is chosen to be a real signal, i.e, a chirp clear that this maximum value of 5 Hz is slightly larger than the frequency bandwidths of each SISO dynamical system (see the Bode plots given in Figure 4.9). This practical limit will play a crucial role as far as model learning accuracy is concerned.

Second, in order (again) to mimic real practical conditions encountered on sensors used for acquiring the output signals y(t) on real passenger vehicles, the noise free data is disturbed by adding, on each output signal, a zero mean white Gaussian noise with a variance selected so that the signal to noise ratio equals 15 dB for each output. One sample of each noisy output signal is given in Figure 4.10.

Given a noisy data set {u k , y k } N k=0 like the one in Figure 4.8 and 4.10, respectively, the noise covariance matrix estimation method described in Section 4.3 is tested again. Because the state sequence estimate Xf,N is not optimal [START_REF] Van Overschee | Subspace identification for linear systems. Theory, implementation, applications[END_REF] especially when short, little informative and noisy data sets are handled, we observe in Figure 4.11 that the reconstructed state variables yielded by our N4SID-like algorithm mimic the main variations of the real state signals only. Such a discrepancy can be attributed to both the noise acting on the system and the modeling error induced by the linear time invariant model given in Eq. (4.44). In order to account for this modeling error in the Kalman filter equations, the effect of these unmodeled dynamics is translated into noise covariance matrix estimates as suggested by Eq. (4.31), leading to the Kalman filter state trajectories given in Figure 4.12. The good performance of the Kalman filter estimates illustrate the efficiency of our approach to transform model uncertainties into reliable noise covariance matrix estimates.

Conclusions

In this chapter, a new method is introduced to solve a common problem hindering the use of a Kalman filter, namely the accurate estimation of the noise covariance matrices appearing in the Kalman filter equa-

CHAPTER 5

Tire force estimation: practical application

Introduction

The objective of this chapter is to apply the estimation method introduced in Chapter 3 in order to infer friction points. Since the estimation method is based on an extended Kalman filter (EKF), it is necessary to find a way to set the noise covariance matrices. In this chapter, this setting is made by using the subspace identification method introduced in Chapter 4. Therefore, this chapter is the application of the methods introduced in Chapter 3 and 4 to solve the problem of estimating friction points (see Figure 5.1). Furthermore, the results obtained in this chapter are used to establish if the estimated friction points can be combined with the MCMC method introduced in Chapter 2 in order to estimate a grip potential value.

In this chapter, the EKF performance is assessed with simulated data. Because the data is simulated, we have access to the real states. Consequently, we can evaluate the performance of the observer by comparing the estimated values to the true values. Two types of simulated data are used, namely:

• data coming from simulations made with an extended single-track model,

• data coming from simulations made with a realistic real-time vehicle simulator called VI-CarRealTime (VI-CRT) [START_REF]Vi-carrealtime 19[END_REF]. On the one hand, the results on the data coming from an extended single-track model allows to check if modeling the longitudinal tire forces with random walks is a suitable choice. On the other hand, the purpose to the test on VI-CRT data is to verify if the EKF can estimate signals coming from a complex vehicle model accurately.

The rest of the chapter is organized as follows. Section 5.2 describes a way to generate data with the extended single-track model. Section 5.3 is dedicated to the EKF estimation of the states generated with the extended single-track model. Section 5.4 describes how data is generated with VI-CRT, then shows the reconstructed signals obtained with the EKF. Section 5.5 concludes this chapter.

Generation of data with an extended singletrack model

This part is dedicated to the generation of simulated data allowing to test the EKF performance. Here, the data is generated by using a single-track model extended with a suspension model such as the one used by the observer (see Chapter 3). The main difference between the model used to generate simulated data and the one used by the EKF lies in the longitudinal tire force expressions. In the EKF, F xf and F xr are modeled with random walks. However, these representations have no physical meaning, thus, are not suitable to generate consistent data. For this reason, F xf and F xr are simulated by using the extended singletrack model only. The equations introduced in Chapter 3 are reminded in Eq. (5.1a)-(5.1l), respectively.

vx = 1 m F xf + F xr - 1 2 ρ a S a C x v 2 x -f RR mg , (5.1a) ωf = 1 2I wf [T f -R f load F xf ] , (5.1b
)

ωr = 1 2I wr [T r -R rload F xr ] , (5.1c) κ + ω 0 Q κ + ω 2 0 κ = h G I y m vx + 1 2 .ρ a S a C x v 2 x , (5.1d)
F zf = mg L r L f + L r + (kκ + b κ) , (5.1e)
F zr = mg L f L f + L r -(kκ + b κ) , (5.1f
)

s f = ω f R rol f -v x max(ω f R rol f , v x ) , (5.1g
)

s r = ω r R rolr -v x max(ω r R rolr , v x ) , (5.1h)
R rol f = R 0 - F nomin K zz [F ref f ρ f + D ref f arctan (B ref f ρ f )] , (5.1i
)

R rolr = R 0 - F nomin K zz [F ref f ρ r + D ref f arctan (B ref f ρ r )] , (5.1j) R f load = R 0 - F zf K zz , (5.1k
)

R rload = R 0 - F zr K zz . (5.1l)
In order to generate simulated data with the extended single-track model (Eq. (5.1a)-(5.1l)), the following assumptions are introduced.

• The longitudinal vehicle speed and acceleration are inputs of the vehicle data simulator.

• The front slip ratio is constant.

• All the engine power is transmitted to the front wheels, thus, T r = 0.

The first assumption is used to generate the data by setting a vehicle speed instead of a front axle torque. This choice is motivated by our knowledge of the order of magnitude of the vehicle speed. The second assumption of constant front slip ratio is used to simplify the data generation. This situation is not representative of the reality because, in real life situation, the longitudinal tire force is generated in the contact patch by the slip ratio [START_REF] Guiggiani | The Science of Vehicle Dynamics[END_REF]. However, even if this situation is not realistic, the generated data can still be used to test the EKF performance as well as the random walk representation. Finally, the last assumption on the engine torque power repartition is a standard assumption used to depict what happens with a traction vehicle [START_REF] Guiggiani | The Science of Vehicle Dynamics[END_REF]. Besides, this assumption leads to the simplifying equation v x = R rolr ω r connecting the vehicle speed v x and the rear wheel speed ω r . Indeed, if T r = 0, the rear longitudinal tire force F xr should be close to zero. However, because the rear tire force is driven by the rear slip ratio s r , F xr ≈ 0 leads to s r ≈ 0. Therefore, setting s r = 0 in Eq. (5.1h) results in the expression v x = R rolr ω r .

After making these simplifying assumptions, the inputs and outputs of the vehicle data simulator are the ones reported in Table 5.1. These signals are also represented in Figure 5.2 In Table 5.1, the derivatives Table 5.1: Inputs and outputs of the vehicle data simulator.

Inputs T r , v x vx , s f , s r Outputs T f , F xf , F xr , F zf , F zr , ω f ωf , ω r ωr , κ κ κ , R rol f , R rolr , R f load , R rload
of some quantities such as v x and ω r are represented in column, e.g., ωr is indicated below ω r . This representation is used to highlight that all the quantities and their derivatives are considered here as one single input or output quantity. More precisely, once one quantity of a column is known, its derivative or its integral is determined by using numerical derivation or integration. Since the data signals are generated without noise, resorting to numerical derivation does not raise any numerical issue. For example, in our case, v x is an input of the simulator. Therefore, we first define a signal v x , then, deduce a signal vx by proceeding to numerical derivation of the the signal v x . Another example can be given with ω r . First, the output signal ω r is inferred from the relation ω r = vx R rolr , then, ωr is deduced from ω r thanks to numerical derivation. Consequently, we are in a situation where we have 12 equations (Eq. (5.1a)-(5.1l)) and 12 outputs (see Table 5.1). In this case, because the number of equations is equal to the number of output signals, data can be generated. The output signals are simulated by using Eq. (5.1a)-(5.1l) according to the plan indicated in Figure 5.3. Besides, the used parameters are the ones indicated in Table 5.2.

In this section, v x is set to correspond to a realistic situation. More precisely, the represented situation consists of different acceleration and deceleration phases between 30 km/h and 130 km/h. Once v x is defined, vx is deduced by numerical derivation of v x . As mentioned previ- ously, because the vehicle data simulator represents a traction vehicle, T r = 0 and s r = 0, respectively. Finally, the front slip ratio is set such as, when the acceleration vx is positive, s f = 0.01 and when vx is negative, s f = -0.01. Therefore when the tires generate a positive traction force F xf , s f is positive and when F xf < 0, s f < 0. As previously explained, the signal representing the front slip ratio s f is not representative of a true slip ratio signal. In the case exposed here, because the signal v x varies, s f should vary accordingly. However, even if this signal is not correct from a physical point of view, it does not prevent to test the EKF performance. The used v x , vx and s f are represented in Finally, in order to have simulated measurements close to the ones obtained with real sensors fitted on production vehicle, a noise is added on the quantities included in the measurement vector of the EKF. In our case, the measurement quantities are the vehicle speed v x , the wheel speeds ω f and ω r and the pitch speed κ. In agreement with real vehicle sensor accuracy [START_REF] Fleming | Overview of automotive sensors[END_REF], the noise used to corrupt the simulated data is an additive zero-mean white Gaussian noise with a signal-to-noise ratio

Extended Kalman filter with extended singletrack model simulated data

As indicated in Section 5.2, the purpose of this step is to test the EKF performance and especially to settle if representing the longitudinal tire forces with a random walk is a suitable choice. After generating the data with the vehicle model, the noisy measurements (see Figure 5.13) are used with the EKF to reconstruct the model state. In order to perform the EKF estimation, it is necessary to set the values of the noise covariance matrices W and V , respectively. Here, these matrices are set with the subspace approach designed in Chapter 4. In this case, the subspace approach provides the following W sub and V sub V sub =

W sub =       5.42e-4 1.13e-3 1.78e-3 -7.70e-3 ... 1.13e-3 1.02e-2 3.52e-3 2.22 ... 1.78e-3 3.52e-3 5.97e-3 -1.23e-1 ... -7.70e-3 2.22 -1.23e-1 6.84e3 ... -1.17e-2 -3.29e-1 3.44e-1 7.08e2 ... 1.85 -1.59e1 -2.51e1 -5.26e5 ... -5.46e-1 -1.90e1 -3.36e1 -1.37e5 ... -2.10e-7 -1.36e-6 -1.06e-6 -2.49e-3 ... 1.11e-8 2.70e-6 -1.02e-6 -1.24e-3 ... ... -1.17e-2 1.85 -5.46e-1 -2.10e-7 1.11e-8 ... -3.29e-1 -1.59e1 -1.90e1 -1.36e-6 2.70e-6 ... 3.44e-1 -2.51e1 -3.36e1 -1.06e-6 -1.02e-6 ... 7.08e2 -5.26e5 -1.37e5 -2.49e-3 -1.24e-3 ... 1.48e3 -1.47e5 -1.50e5 -2.63e-3 -3.25e-3 ... -1.47e5 6.19e7 1.64e7 4.15e-1 5.48e-1 ... -1.50e5 1.64e7 1.67e7 3.27e-1 2.66e-1 ... -2.63e-3 4.15e-1 3.27e-1 1.43e-7 -4.91e-8 ... -3.25e-3 5.48e-1 2.66e-1 -4.91e-8 3.91e-8       . ( (5.2)
By setting the noise covariance matrices W and V equal to W sub and V sub respectively, the EKF provides the estimated state signals indicated in Figures 5.14-5.20. In addition, in order to provide estimates of W and V , the subspace approach proceed to a state estimation with the N4SID method. These estimated states are also indicated in Figures 5.14-5.20. Then, these two types of estimates (the ones given by the EKF and the others given by N4SID) are used to compute approximations of the normal loads, the frictions, the slip ratios, the loaded tire radii and the effective tire radii. These estimates are represented in Figures 5. 21-5.30.

In order to evaluate the EKF performance, the errors between the estimated and the true values are plotted for each quantity inferred. Twotype of error representations are illustrated: (i) the temporal evolution of the residuals, (ii) error histograms. Besides, a criterion quantifying the quality of the estimates is computed. This goodness of fit criterion (GF C) quantifies the error between the estimated signals and the true signals according to the equation

GF C = 1 - ∥x true -x∥ 2 2 ∥x true -mean (x true )∥ 2 2 . (5.4)
Hence, according to Eq. (5.4), the GF C defined here varies between -∞ and 1. The closer to 1 the GF C is, the best the estimates are. The evaluated GF C are available in the figure titles. According to Figures 5.14-5.20, except for F xr , the EKF is able to reconstruct the dynamics of the state signals. This result is confirmed through the GF C values with a GF C greater than 66% for F xf and greater than 88% for the other state variables. In our case, F xr is badly reconstructed. However, the data generated in Section 5.2 mimics the behavior of a traction vehicle. Consequently, the dynamics of F xr are very low which can explain the reasons why the EKF poorly estimates F xr . Furthermore, as depicted in Figures 5.14, 5.15, 5.16 and 5.20, the EKF provides accurate estimates of the noisy measured quantities v x , ω f , ω r and , κ. Therefore, because the state values are accurately estimated, the normal loads and the different radii are well estimated as well with a GF C > 88% (see Figures 5.21, 5.22 and 5.27-5.30). Because F xf is not so well estimated, the front friction µ f is accordingly not as well estimated as the other variables with a GF C ≈ 64% (see Figure 5.23). Naturally, µ r is poorly estimated because F xr is badly estimated. Finally, the slip ratios s f and s r are also wrongly estimated. These results can be explained by the consistency of the used data. Indeed, as mentioned in Section 5.2, s f is a constant piecewise signal and s r = 0. However, because the longitudinal tire forces are physically generated by the slip ratios, once the tire forces vary, the slip ratios should vary accordingly. This inconsistency can be one reason of the poorly slip ratio estimates.

Regarding the estimates provided by the N4SID algorithm, their accuracies vary from one estimate to another. First of all, the different speeds v x , ω f and ω r are as well estimated as the EKF estimates with a GF C > 98% (see Figures 5. 14-5.16). Although the pitch speed κ is also quite well estimated with a GF C > 86%, we can observe a reduction of the accuracy compared to the EKF estimate (see Figure 5.20). Finally, as shown in Figure 5.17, the longitudinal tire force F xf estimated with N4SID is less accurate than the one provided by the EKF. More precisely, the N4SID estimate is less noisy than the EKF estimates, however, it is not centered on the true value contrary to the EKF estimate.

Furthermore, some comments might be made on the histograms (Figures 5. 14-5.30). Firstly, the sampling period T s is equal to 1 ms and the simulation lasts 100 s. Consequently, the different signals are formed of 100000 samples and thus, the histograms contain 100000 samples. Therefore, although the histograms take into account of only one realization, they still include a large number of samples and thus the strong law of large number can be assumed to hold. Besides, as reminded in Chapter 4, as long as the sample number tends to infinity, the time average covariance matrices generated from the residuals converges to the ensemble covariance matrices V and W with probability one. Hence, a high number of samples increases the confidence on the estimated noise covariance matrices V sub and W sub . Secondly, in the case treated here, the noise added on the measurements is zero mean and Gaussian. Therefore, the residuals should be Gaussian too. However, the histograms of the N4SID residuals are not always Gaussian. This result suggests that the N4SID method is not enough to completely reconstruct the states and thus it is necessary to introduce an EKF in order to improve the estimates.

Although the longitudinal tire forces are quite well estimated with a GF C ≈ 66%, we can wonder if the EKF might provide more accurate estimates with a different setting of the noise covariance matrices W and V . Indeed, because the longitudinal tire forces are modelled with random walks, the subspace identification procedure can encounter some difficulties to evaluate the coefficients related to F xf , F xr , Ḟxf and Ḟxr . More precisely, the subspace identification procedure needs to know a linearised state-space representation of the vehicle model used in order to move the transformation matrix T into the right basis (see Chapter 4). However, because the longitudinal tire force model has no physical meaning, the corresponding part of the linearised model has no physical meaning as well. Consequently, the approach may fail to move T into the right basis. This assumption is confirmed by the N4SID estimates of the longitudinal tire forces given in Figures 5.17 In order to suggest another setting W mix and V mix for W and V , the idea is to start from the estimates W sub and V sub , then to modify the coefficients of W sub corresponding to the state poorly estimated with our N4SID-like algorithm. Since the EKF successes to reduce the measurement noise as mentioned previously, V mix = V sub . In the same vein for the process noise covariance matrix W , the plan is to keep the diagonal elements of W sub except the ones corresponding to F xf , F xr , Ḟxf and Ḟxr and to set the missing coefficients with an error-trial phase. Under these practical conditions, the number of coefficients to set is only equal to four. This plan leads to the following process noise covariance matrix

W mix = diag             5.42e-04 1.02e-02 5.97e-03 1.00e+02 1.00e-02 1.00e+04 1.00e+00 1.43e-7 3.91e-8             .
(5.5)

By using the EKF with W mix and V mix , the observer provides the longitudinal tire forces and frictions indicated in Figures 5. [START_REF] Garatti | Assessing the model quality in system identification: the asymptotic theory revisited and application to iterative control[END_REF] In order to compare the results obtained with both settings of the EKF, Table 5.3 gathering the GF C obtained with both methods is shown. In addition, Table 5.4 and Table 5.5 comparing respectively the means and medians of the absolute errors estimated are also provided.

As shown in Figure 5.31 and Tables 5.3-5.5, setting the EKF with W mix and V mix increases the accuracy of the estimated F xf significantly with a reduction of 83% of the absolute error median and a GF C of 88%. Consequently, the precision of the estimated front friction increases accordingly with the same values for the GF C and the decreasing of the absolute error median (see Figure 5.33 and Table 5.5). Furthermore, although the rear longitudinal tire force signal has low dynamics, the EKF set with W mix and V mix is able to track them quite well as depicted in Figure 5.32. More precisely, the GF C increases to 61% and the median absolute error decreases by 97.8% when compared to the estimates obtained with the past EKF settings (see Tables 5.3 and 5.5).

Finally, the results obtained in this part indicate that it is possible to estimate tire forces and thus friction by using the random walk model. By using the subspace approach described in Chapter 4, we can set most of the coefficients of the EKF noise covariance matrices. Consequently, it only remains to adjust the coefficients of the states which are not well estimated by the N4SID-like algorithm, i.e., only the four coefficients corresponding to the longitudinal tire forces and their derivatives, respectively.

Extended Kalman filter with VI-CRT simulated data

Data generation with VI-CRT

The next step for the validation of the vehicle state observer is to evaluate its performance using more realistic simulation data sets. These simulations are carried out with the software VI-CarRealTime (VI-CRT) developed by VI-Grade company [START_REF]Vi-carrealtime 19[END_REF]. VI-CRT is a virtual vehicle modeling and a real-time vehicle simulation software. With this software, it is possible to make simulations of vehicles in various real life situations. For example, this software allows the user to select a type of car and tire, or to create different types of roads. The different vehicle models used by VI-CRT are much more complete than the single-track and doubletrack models introduced beforehand. Besides, the vehicle and tire parameters are obtained using vehicle and tire test benches and elaborated fitting methods. Hence, VI-CRT gives access to very realistic data describing the vehicle dynamics. In addition, it should be noted that VI-CRT provides data signals for the four vehicle wheels. Figure 5.35 shows the VI-CRT main interface on which is selected the choice of the vehicle subsystem such as the tires and the suspensions. VI-CRT is used here as a black-box simulator.

In order to generate data with VI-CRT, it is necessary to select a vehicle model, a tire model and to define a maneuver. The chosen vehicle model is the CompactCar model provided by VI-CRT. This model is selected because it represents the behavior of compact class vehicles which are vehicles commonly used by drivers. Besides, it must be pointed out that this vehicle model is a traction model. Therefore, the driving torque is applied exclusively to the front axle of the vehicle. As a consequence, the rear tires generate only small residual tire forces. The selected tire model is a Pacejka's tire model designed for the Compact-Car model which is provided by VI-CRT.

Once the vehicle and the tire model have been selected, it remains to define the maneuver settings such as • the vehicle direction (exclusively straight line or presence of turns),

• the vehicle speed or acceleration,

• the experience duration.

In this manuscript, we focus on the longitudinal vehicle behavior. For this reason, the simulated maneuver consists of a succession of acceleration and deceleration phases in a straight line. The vehicle acceleration is set such as the vehicle speed increases from 80 km/h to 120 km/h during 6 s and then, decreases to 80 km/h during around 16 s. This succession of acceleration and deceleration phases is repeated seven times for a total duration time of 200 s. The acceleration and deceleration phases are repeated in order to verify if the estimates provided by the observer gain in precision as the different phases are executed. The different signals genetared by VI-CRT are given in Figures 5. 37-5.45. Especially, Figure 5.38 shows the evolution of the vehicle speed during the different phases. Besides, in Figure 5.42, it can be pointed out that the front friction signal µ f reaches values around 0.5 while common values are around 0.3. This amount of friction is due to the maneuver in which the speed increases significantly (from 80 km/h to 120 km/h) during a short amount of time (6 s).

Finally, such as in Section 5.2, in order to be closer to real measurements provided by sensors fitted on production vehicles, a noise is added on the quantities included in the measurement vector of the EKF. Once

Results of EKF estimates

After generating data with VI-CRT, the signals corresponding to the observer inputs (the front and rear torques T f and T r ) and the measurements (see Figures 5.46) are combined with the EKF in order to proceed to state estimation. Once again, the EKF is tuned by applying the subspace identification procedure described in Chapter 4. This one provides the following results W sub and V sub V sub = 7.80e-02 -5.50e-02 2.09e-02 -4.34e-04 -5.50e-02 6.89e+00 -1.40e+00 4.42e-02 2.09e-02 -1.40e+00 1.19e+00 -1.04e-02 -4.34e-04 4.42e-02 -1.04e-02 3.26e-04

,

(5.6) nals and the GF C > 88%. Furthermore, the pitch angle κ is also well estimated by the EKF with a GF C > 89% (see Figure 5.52). However, as indicated in Figure 5.53, the pitch speed κ is surprisingly poorly estimated by the EKF. Especially, the estimated signal is less accurate than the measurement signal. This result indicates an incorrect setting of the coefficients corresponding to κ in W sub and V sub . This statement is confirmed by the wrongly pitch speed κ estimated with N4SID. Nonetheless, this lack of precision on κ does not have a significant impact on the EKF estimated normal loads. Indeed, as shown in Figures 5.54 and 5.55, F zf and F zr are accurately estimated with a GF C > 84% and errors around 50 N for an order of magnitude around 9000 N and 6500 N , respectively. Moreover, with the noise covariance matrices set with W sub and V sub , the front longitudinal tire force F xf is estimated by the EKF with accuracy (see Figure 5.50) with a GF C > 92%. This result constitutes a valuable progress compared to the previous results obtained by using data generated with the extended single-track model. This result may seems surprising because the longitudinal tire forces estimated with the N4SID-like algorithm are quite far from the actuals signals. However, we can found an explanation by looking to the process noise covari-ance matrix coefficients. According to Table 5.6, while the coefficients of W sub corresponding to F xf and F xr are closed in the cases where the EKF is tested on the extended single-track model generated data and the VI-CRT data, the coefficients corresponding to Ḟxf and Ḟxr differ considerably. In the first treated case, the order of magnitudes of these coefficients is in 10 7 while in the second case it is in 10 3 . Therefore, this gap in the coefficient values could explain why setting W to the estimate W sub provides accurate results of the longitudinal tire forces. Furthermore, in the VI-CRT data case, the coefficients of W sub related to Ḟxf and Ḟxr are close to the ones proposed to improve the EKF estimates in Section 5.3 (see Eq. (5.5)). This observation reinforces the idea that the difference in the estimates accuracy comes from the coefficients related to Ḟxf and Ḟxr . As a consequence, the front friction µ f is accordingly well estimated with a GF C > 92% (see Figure 5.56). However, as shown in Figures 5.50 and 5.56, we can observe peak values on the absolute errors curves. These ones appear each time the vehicle dynamic abruptly changes. These results indicate a small delay of the EKF response. Nonetheless, the combination of the subspace approach and the EKF provides accurate friction estimates.

Another important quantity necessary to estimate here is the slip ratio s f . Unfortunately, as represented in Figure 5.58, this one is not well estimated with a GF C ≈ 49%. More precisely, the reconstructed signal is able to reproduce the main variations of the true signal but it is corrupted by an important noise. This lack of precision is due to the way s f is inferred. Here, s f is computed by using Eq. (5.1g). Therefore, it is necessary to divide ω f R rol f -v x (a small quantity) by ω f R rol f or v x which are noisy estimates. Even if these ones are really well estimated (GF C > 96%), the fact that it is necessary to divide a small quantity by a noisy signal result in a poor estimate of s f . For this reason, if we want to improve the precision of the estimated s f , we should either change the way we estimate it or use slip ratio measurements. As reported in Figure 5.58, the slip ratio estimated with the N4SID-like algorithm is quite accurate. Therefore, using a black-box method such as N4SID could be a solution to increase the estimated slip ratio.

On another side, because the simulated vehicle is a traction vehicle, the rear dynamics are really weak (see Figures 5.51,5.57 and 5.59). For this reason, the EKF is unable to reconstruct the signals corresponding to the rear dynamics with precision.

Finally, the sampling period of the VI-CRT data is T s = 0.01 s and the simulation last 216 s. Therefore, the different signals are constituted of 21600 samples. Once again, the residuals computed with the N4SID estimates do not describe Gaussian curves (see .

Proposition of another setting of the noise covariance matrices

As in the Section 5.3, another setting of the covariances matrices V mix and W mix is proposed in order to see if we can improve the results on quantities of interest, i.e., s and µ. This new setting consists in keeping the matrix V given by the subspace approach (V mix = V sub ) and modifying some coefficients of the matrix W . More precisely, although the front longitudinal tire force is estimated with a high level of accuracy (GF C > 92%), the four coefficients corresponding to the longitudinal tire forces are slightly modified in order to check if a different setting could improve even more the estimates of the longitudinal tire forces. In addition, because the pitch speed is not really well estimated, the corresponding coefficient of W is also modified. Under these conditions, W mix is given by

W mix = diag              
              (5.8)
Hence, only 5 coefficients of the process noise matrix are modified. The remaining diagonal coefficients of W mix and the measurement noise matrix V mix are tuned with the results of the subspace identification approach.

By using the new setting W mix and V mix , the EKF provides the estimated pitch speed indicated in Figure 5.64, the longitudinal tire forces represented in Figures 5. [START_REF] Nocedal | Numerical optimization[END_REF] Besides, as in Section 5.3, Tables reporting the GF C obtained with both methods (see Table 5.7) and comparing the means (see Table 5.8) and medians (see Table 5.9) of the estimated absolute errors are also provided.

As shown in Figure 5.64 and Tables 5.7-5.9, setting the EKF with the proposed W mix and V mix increases significantly the accuracy of κ with a GF C rising from 12% to 67% and a reduction of more than 85% of the median absolute error. Besides, this new setting has also a positive influence on the estimated F xr . As represented in Figure 5.66, the EKF estimate of F xr is less noisy and follows slightly better the signal dynamics. However, these benefits on the pitch speed and the rear dynamics estimates do not improve the results on the front longitudinal tire force F xf and more importantly on the front friction µ f and slip ratio s f (see Figure 5.65, 5.67, 5.69 and Tables 5.7-5.9). Hence, although the new setting of the noise covariance matrices W and V increases slightly some EKF estimates, it is not sufficient to enhance the accuracy of the estimated friction and slip ratio which compose the friction data points.

Setting of the noise covariance matrices including the noise statistics

Finally, because we are in a simulation case, we have not only access to noisy measurements but also to the true signals and thus, to the measurement noise statistics. Therefore, in this case, we can propose another setting W mix2 and V mix2 of the noise covariance matrices based on these statistics. This new setting consists in modifying both matrices W and V . More precisely, V is set by using the standard deviations of the noise added on the measurements. 1 Note that in real life situations, this information is not available. Therefore, the new setting

V mix2 =
2.76e-01 0.00e+00 0.00e+00 0.00e+00 0.00e+00 9.34e-01 0.00e+00 0.00e+00 0.00e+00 0.00e+00 9.20e-01 0.00e+00 0.00e+00 0.00e+00 0.00e+00 2.22e-04

(5.9)

Because the measurement noises are not involved in the process noise covariance matrix, W mix2 is tuned to the previous proposed setting. Thus, W mix2 = W mix . This new setting W mix2 and V mix2 leads to the estimated pitch speed indicated in Figure 5.71, the longitudinal tire forces represented in Figures 5.72 Tables reporting the GF C obtained with the EKF estimates set with W sub , V sub and W mix2 , V mix2 , respectively (see Table 5.10) and com- paring the means (see Table 5.10) and medians (see Table 5.12) of the estimated absolute errors are also provided.

As shown in Figure 5.71 and Tables 5.10-5.12, setting the EKF with W mix2 and V mix2 increases even more the accuracy of κ than W mix and V mix . In this case, the GF C increasing from 12% to 90% and we observe a reduction of more than 94% of the median absolute error.

However, compare to the results obtained with the EKF set with W mix and V mix , this new setting does not improve the results on the quantities of interest, i.e., the friction µ f and slip ratio s f (see 5.74, 5.76 and Tables 5.10-5.12). Thus, setting the noise covariance matrices W and V to W mix and V mix seems to represent the best compromise. On the one hand, this setting is mainly based on the subspace identification procedure. Therefore, the user needs only to change 5 coefficients. On the other hand, except for the pitch speed, this setting provides equivalent estimates than the EKF set with the knowledge of the noise statistics.

Conclusion

In this Chapter, the EKF introduced in Chapter 3 have been applied to simulated data generated with an extended single-track model and the real-time software VI-CRT. By setting the noise covariance matrices thanks to the subspace approach introduced in Chapter 4, the EKF is able to reconstruct accurately tire force signals and thus to estimate the friction with a high level of accuracy. Therefore, although the subspace approach has been design for a linear Kalman filter, performing this approach with an EKF provides attractive results. Besides, this result pointed out that modeling the longitudinal tire forces with random walks constitutes a suitable choice. Hence, combining the subspace approach and the EKF forms an accurate estimation method which have the main advantage to require almost no adjustments from the user. However, the levels of accuracy on the estimated slip ratios reveal to be insufficient. As a consequence, the estimated friction points do not coincide with the true friction curve (see Figure 5.78). If the slip ratio were perfectly estimated, the friction points obtained would be the ones indicated in Figure 5.79. Therefore, before proceeding to MCMC estimation with the estimated friction points, it is necessary to either use slip ratio measurements or to find a way to improve the estimated slip ratio.

CHAPTER 6

Conclusion and perspectives

Conclusion

The objective of this thesis has been to estimate the tire friction potential under standard driving conditions by using the measurements provided by the sensors fitted on production vehicles. After establishing, in Chapter 1, a state of the art about the different methods existing in the literature to solve this problem, it has been decided to design an effect-based method belonging to the tire slip-based class of solutions. In addition, because it exists an important margin of progress to improve the results in the purely longitudinal case, the study conducted has been focused on these dynamics. In order to simplify the resolution of the global issue, this one has been divided into two sub-problems. Firstly, the estimation of a grip potential value from friction point measurements restricted to standard driving conditions. Secondly, the estimation of friction data points with measurements fitted on production vehicle.

The first sub-problem has been treated in Chapter 2. In this part, a solution based on a combination of a classical Maximum Likelihood (ML) approach and an adaptive MCMC approach has been suggested. The results have been validated with both simulations and real data. In particular, real measurements coming from a flat track tire testing machine have been used. These results have shown that, when only measurements corresponding to standard conditions are considered, the suggested MCMC technique delivers more accurate results than the tra-ditional ML method. More precisely, the MCMC method has been able to infer an accurate value of the grip potential once friction measurements of approximately 0.2 were available.

The second sub-problem has been tackled is Chapters 3-5. Chapter 3 has been mainly focused on the theoretical aspects of the method developed to estimate friction data points. In particular, this method has been based on the use of an Extended Kalman Filter (EKF). Therefore, in this part, it has been necessary to introduce a vehicle model and to put this one under a state-space representation form. The vehicle model considered was a single-track model extended with a suspension model. This one has been chosen because it represents a good trade-off between complexity and accuracy. Besides, based on L. Rays's work presented in Chapter 1, the longitudinal tire forces have been modeled with random walks. Furthermore, because the estimation method involves the use of an EKF, it has naturally raised the question of the setting of the noise covariance matrices V and W . This question has been important in our case because the state vector used here contains nine quantities. Therefore, it could be complicated to set the nine coefficients of the process noise matrix by hand during a trial-error phase. For this reason, a method which has the advantage to require no tuning from the user has been introduced in Chapter 4. This method has been based on subspace identification. More precisely, with this method, a discrete time model estimated with a subspace identification approach has been compared with the deterministic part of a Kalman filter model. Then, the comparison between these two models has led to residuals which have been used to estimate the noise covariance matrices. Although this method has been designed for Kalman filter, i.e., for linear Kalman filter, this approach was found to work well with an EKF. Finally, Chapter 5 has been dedicated to the application of the methods developed in Chapters 3 and 4 in order to estimate friction data points. These methods have been assessed with simulated data coming, on the one hand, from an extended single-track model and, on the other hand, from the realistic vehicle simulator VI-CRT. The results have shown that combining the subspace identification procedure with the EKF leads to accurate estimate of the friction signals. In particular, the use of the subspace identification method developed in Chapter 4 has simplified drastically the setting of the noise covariance matrices. Sometimes, it has been necessary to adapt some coefficients of the noise covariance matrices corresponding to the tire forces. However, in any cases, the subspace approach has allowed to reduce significantly the number of coefficients needed to be set. Nevertheless, although the EKF has provided accurate friction signal, the estimated slip ratio was found to be corrupted by noise and thus, unusable. Hence, others methods such as the ones suggested in Section 6.2 should be investigated in order to improve these results.

Perspectives

In view of the results obtained during this thesis, several perspectives might be considered. Since this work have been realized in an industrial context, these perspectives are classified between short, medium and long term targets from an industrial point of view.

Short term perspectives

The first short term perspective concerns the validation steps required in order to enhance the efficiency of the proposed approach. The MCMC based approach introduced in Chapter 2 needs to be tested on different type of roads such as wet road or icy road. In particular, it would be interesting to assess this method with friction data points similar to the ones provided by the flat track (see Chapter 2), but corresponding to a wet road where µ max ≈ 0.6. This test is important because it would allow to know if the MCMC-based method is able to capture a reduction of the grip potential under standard driving conditions.

Another important validation step concerns the estimation of friction data points. In Chapter 5, the tests were carried out with simulation data only. Therefore, the next stage is to test the observer with real measurements coming from the CANbus of a vehicle. This test is essential because, in real life situations, these measurements are the only ones available. Furthermore, this step raises questions about result validation because in the CANbus, the true friction point values are not known. Thus, in this case, one important challenge is to find a way to conclude on the reliability of the results provided by the EKF.

One of the main concern of manufacturers is the real-time estimation of the friction potential. Therefore, an important perspective it to find a way to reduce the computational time of the MCMC based method. As a reminder, depending on the size of the training set, simulating 1000 chains of 100000 samples have taken between 2h and 2h45 (see Chapter 2). With the MCMC procedure, most of the time is spent during the repetition of the sample generation for each chain. Consequently, an idea to reduce the computational time would be to generate samples depicting faster the target distribution. One solution provided in the literature is to resort to Hamiltonian Monte-Carlo (HMC) [START_REF] Betancourt | A conceptual introduction to hamiltonian monte carlo[END_REF][START_REF] Gelman | Bayesian data analysis[END_REF][START_REF] Särkkä | Bayesian Filtering and Smoothing[END_REF][START_REF] Rogers | A first course in machine learning[END_REF]. In theory, this method should generate samples closer to the target distribution and thus, decrease the number of samples necessary and thus, the time computation. If the HMC technique does not provide interesting results, another idea is to resort to a recursive approach where the result given by the MCMC method is adapted at each new measurement.

Finally, the slip ratio estimation may be the subject of another major short term perspective. As pointed out in Chapter 5, the estimated slip ratio provided by the EKF is noisy and quite far from the true value. However, since the friction data points are composed of the slip ratio and the friction, a wrong estimation of the slip ratio signal is a major problem. One solution to this issue might be to resort to slip ratio measurements. However, in some cases, even slip ratio measurements are strongly corrupted by noise. An example of this type is illustrated in Figure 6.1 where a curve representing the longitudinal tire force F x as a function of the slip ratio s is drawn. The curve represented in Figure 6.1 is plotted with longitudinal tire force measurements and slip ratio measurements coming from sensors fitted on a tractor. As can be seen, for a given F x , numerous different slip ratio values are given. In these conditions, it is difficult to extract an accurate curve representing the real evolution of the the longitudinal tire force as a function of the slip ratio. As a consequence, using slip ratio measurements would not necessary help to infer accurate slip ratio values and another way should be examined. Based on the results obtained in Chapter 5, one solution could be to resort to black-box methods. As a reminder, in Chapter 5, contrary to the EKF, the N4SID algorithm provides a slip ratio estimated signal with little noise (see Figure 5.56).

Mid term perspectives

Once short term perspectives have been investigated, mid term perspectives which require further developments can be mentioned. One of them concerns the setting of the noise covariance matrices V and W . In Chapter 4, the proposed method to set the noise covariance matrices provides constant estimates V and Ŵ . However, during a travel, the considered vehicle model is more or less appropriate depending on the situations. Therefore, it should be interesting to consider time varying noise covariance matrices V and W . One solution would be to estimate V and Ŵ in a recursive way. More precisely, the idea of the method would be similar to the one presented in Chapter 4. However, instead of using the N4SID algorithm, we would use a recursive subspace method such as the ones presented in [START_REF] Mercère | A new recursive method for subspace identification of noisy systems: EIVPM[END_REF][START_REF] Mercère | Recursive subspace identification based on instrumental variable unconstrained quadratic optimization[END_REF][START_REF] Mercère | Convergence analysis of instrumental variable recursive subspace identification algorithms[END_REF]. Hence, estimates of V and W could be inferred at each time step.

A second mid term perspective involving some new developments consists in identifying some vehicle and tire parameters instead of fixing them. More precisely, in Chapter 5, the used vehicle and tire parameters correspond to the ones used by the vehicle simulator VI-CRT. Therefore, these parameters are the ones of a complex vehicle model designed by VI-CRT. However, the extended single-track model used to model the vehicle behavior is a simplified model compared to the VI-CRT model. Consequently, although some parameters are common and should not vary between both models, e.g., the vehicle mass m, other parameters might change. For example, since the suspension description is quite different between the extended single-track model and the VI-CRT model, the corresponding stiffness and damping coefficients should be different too. Hence, one idea to improve the results is first to identify some parameters with an identification procedure and then, use these estimates with the EKF in order to determine friction data points.

Long term perspectives

Finally, long term perspectives may be cited. These ones would require further investigations and important modifications of the approach introduced in this manuscript. A first long term perspective would be to resort to another tire model taking into account other effects affecting the friction curve such as the tire internal temperature T i or the tire inflation pressure, to name a few. As pointed out by F. Gustafsson in [START_REF] Gustafsson | Slip-based tire-road friction estimation[END_REF], the tire internal temperature influences the bottom of the friction curve and thus, the friction data points situated in this region. This temperature effect is illustrated in Figure 6.2 where friction curves corresponding to two different type of roads and two different internal temperatures have been plotted. As depicted in Figure 6.2, for one type of road, the tire internal temperature has a significant impact on the bottom of the friction curve and a low impact on the friction potential value. Thus, including the tire internal temperature could improve the reliability of the grip potential estimated value. However, many challenges need to be faced in order to achieve this perspective. First of all, considering the temperature effect requires the use of a tire model describing the friction curve and taking into account the tire internal temperature. However, this type of tire model is not common and most of the tire models designed in the literature are purely mechanical models [START_REF] Svendenius | Tire modeling and friction estimation[END_REF]. Nevertheless, some thermomechanical tire models exist such as the model TameTire [START_REF] Durand-Gasselin | Assessing the thermo-mechanical tametire model in offline vehicle simulation and driving simulator tests[END_REF][START_REF] Pearson | Tametire: Introduction to the model[END_REF] developed by Michelin company. However, using the TameTire model is not as easy as using the Pacejka's model. For instance, one of the main advantages of the Pacejka's tire model is its parsimony, i.e., it can accurately describe the friction curve by using a few number of parameters. On the contrary, the TameTire model involved more than 170 parameters. Among them, approximately 20 parameters have a more or less strong influence on the friction curve. Therefore, one essential point would be to classify the influence of the different parameters on the friction curve depending on the situations.

Another long term perspective would be to adapt the MCMC-based estimation method for different dynamics. In particular, in the purely lateral dynamics case, a curve similar to the used friction curve exist.

Instead of representing the ratio Fx

Fz as a function of the slip ratio s, this similar curve represents the ratio Fy Fz as a function of the slip angle s a . Consequently, the MCMC-based method can be adapted to the purely lateral case. One interesting point of this perspective is the possibility to combine the MCMC-based approach with the SAT method (see Chapter 1) which is one of the most promising method to estimate the grip potential in a purely lateral case. Nevertheless, adapting the method to this case requires some modifications such as the vehicle model considered in order to represent the steering kinematics.

In the same spirit, another perspective would be to extend the MCMC-based approach to coupled dynamics. This step would make the method more general and thus more applicable. However, as for the adaptation to the lateral dynamics, extending the MCMC method to the coupled dynamics requires to model the steering dynamics. In addition, in this case, the definition of the grip potential changes and begins

µ max = max   F 2 x + F 2 y F z   . (6.1)
Furthermore, in the longitudinal case, we work with a friction curve representing the normalized longitudinal tire forces as a function of the slip ratio. As mentioned previously, in the lateral case, we have a similar curve representing the normalized lateral tire force as a function of the slip angle. Therefore, in the coupled dynamics case, two dimensions should be considered and thus, we should not work with friction curve but with friction surface.

Résumé étendu

Introduction et état de l'art Introduction

De nos jours, les accidents de la circulation causent toujours de nombreuses victimes chaque année. Pour cette raison, l'amélioration de la sécurité des passagers est une des priorités pour les constructeurs automobiles et les manufacturiers. C'est dans cette optique que de nombreux systèmes d'aide à la conduite tel que le système antiblocage des roues (ABS) et le correcteur électronique de trajectoire (ESP) ont été développés. Ces systèmes permettent de corriger la trajectoire du véhicule en cas de situations dangereuses où encore permettent d'éviter le blocage des roues. Actuellement, avec le développement des véhicules autonomes, les systèmes d'aide à la conduite qualifiés d'avancés (ADAS) sont de plus en plus étudiés. Une des méthodes permettant d'améliorer les ADAS est de leur fournir des informations détaillées sur certaines grandeurs influençant le comportement du pneu. Une de ces grandeurs est le potentiel d'adhérence du pneu. Cette quantité représente l'effort maximum qu'un pneu peut transmettre au sol avant de se mettre à patiner sur la route. Le potentiel d'adhérence peut se révéler particulièrement intéressant pour développer les ADAS car il permettrait à ceux-ci de freiner avec la quantité optimale de pression de freinage. De plus, il pourrait permettre également de détecter les portions de routes ayant une faible adhérence et mettrait donc en évidence les routes usées ou encore celles ayant une mauvaise évacuation de l'eau.

Cette liste non exhaustive des applications possibles grâce à la connaissance du potentiel d'adhérence justifie la nécessité pour les manu-facturiers de développer des méthodes d'estimation permettant de calculer le potentiel d'adhérence avec précision. Bien qu'il existe de nombreuses méthodes élaborées dans la littérature pour estimer le potentiel d'adhérence, celles-ci ne résolvent pas toutes les difficultés liées à l'estimation du potentiel d'adhrence. C'est pourquoi l'objectif principale de cette thèse est de développer une méthode d'estimation du potentiel d'adhérence qui permet de s'affranchir des difficultés rencontrées par les méthodes décrites dans la littérature.

Etat de l'art

La littérature contient de nombreuses méthodes permettant d'estimer le potentiel d'adhérence comme le montre les nombreux articles de synthèses [START_REF] Acosta | Road friction virtual sensing: A review of estimation techniques with emphasis on low excitation approaches[END_REF][START_REF] Acosta | Virtual tyre force sensors: An overview of tyre model-based and tyre model-less state estimation techniques[END_REF][START_REF] Singh | Literature review and fundamental approaches for vehicle and tire state estimation[END_REF][START_REF] Rajamani | Tire-road friction-coefficient estimation[END_REF][START_REF] Khaleghian | A technical survey on tire-road friction estimation[END_REF][START_REF] Müller | Estimation of the maximum tire-road friction coefficient[END_REF][START_REF] Guo | A review of estimation for vehicle tire-road interactions toward automated driving[END_REF]. Ces articles proposent différentes classifications des méthodes qui varient en fonction de l'article de synthèse sélectionné. Ainsi, la classification illustrée par la Figure 6.3 est proposée dans ce manuscrit.

Approche mise en oeuvre

En examinant la littérature, deux contraintes font principalement obstacles aux différentes méthodes d'estimation du potentiel d'adhérence. Tout d'abord, étant donné que cette grandeur est destinée à être utilisée dans un contexte industriel, c'est-à-dire sur des véhicules de série, il est primordial de développer une méthode économique qui ne requiert pas d'autre capteurs autres que ceux déjà présent sur les véhicules de série. De plus, le potentiel d'adhérence doit être déterminé dans des conditions standards de conduite. Comme illustré sur le Figure 6.4, ces conditions de conduite correspondent à des points de friction ne dépassant pas 0.3. En examinant la littérature certaines méthodes comme celles se basant sut le couple d'auto-alignement (méthode SAT) semblent fournir des résultats prometteurs en dépits des deux principales contraintes précédemment citées. Cependant, un équivalent de la méthode SAT n'est pas possible dans le cas longitudinal. Par ailleurs, il semble qu'il existe une importante marge d'amélioration des résultats pour ce type de dynamique. Ainsi, la méthode considérée ici se concentre sur les dynamiques longitudinales. En outre, afin de résoudre le problème d'estimation du potentiel d'adhérence en respectant les deux contraintes majeures, la solution employée ici propose de diviser le problème en deux parties :

• comment estimer le potentiel d'adhérence à partir de points de friction ?

• Comment déterminer des points de friction en utilisant uniquement les mesures fournies par les capteurs équipant les véhicules de série ?

Cette simplification du problème principale est représentée sur la Figure 6.5.

Estimation du potentiel d'adhérence

Cette section concerne la première étape du plan générale mis en place pour estimer une valeur du potentiel d'adhérence µ max dans des conditions standards de conduite, c'est-à-dire avec des points de friction situés en-dessous de la limite µ = 0.3. Ainsi, dans un premier temps, notre étude va se focaliser sur les méthodes d'estimation du potentiel d'adhérence en supposant que l'on a accès à des mesures des points de friction (voir Figure 6.6).

La méthode mise en oeuvre pour estimer le potentiel d'adhérence est illustrée sur la Figure 6.7. L'idée de la méthode est de combiner • des mesures des points de friction,

• un modèle pneu décrivant la courbe de friction ,

• une méthode d'estimation paramétrique.

Le modèle pneu choisi est le modèle semi-empirique Pacejka. Ce modèle parcimonieux à l'avantage de représenter fidèlement la courbe de friction en utilisant simplement six paramètres.

La méthode d'estimation mise en oeuvre est une méthode de Monte-Carlo par chaines de Markov (MCMC). Cette méthode est basée sur l'algorithme de Metropolis-Hasting en marche aléatoire. De plus, les différentes chaines de l'algorithme sont initialisées avec une estimation fournie par la méthode du maximum de vraisemblance (ML). Enfin, afin d'améliorer encore plus la méthode, on y ajoute un a priori physique. Cet a priori consiste à ne garder que les courbes qui ont un glissement correspondant au maximum de la courbe réaliste. Finalement, la méthode mise en oeuvre est appliquée avec des mesures venant d'une machine de test de pneu. Cette machine à l'avantage de nous fournir des mesures de la courbe de friction complète. Ainsi, ces mesures permettent de valider que la courbe estimée colle bien à toutes les mesures de la courbe. En utilisant les trois ingrédients cités précédemment avec des mesures des points de friction correspondant à une conduite standard (µ < 0.3), la combinaison de la méthode MCMC et du modèle Pacejka fournit les courbes représentées sur la Figure 6.8. Sur la Figure 6.8, on peut observer que la méthode MCMC avec ajout de l'a priori physique est très proche de la courbe de référence. En particulier, cette courbe est bien plus proche que celle obtenue avec la méthode ML. Maintenant que la méthode a prouvé qu'elle pouvait fournir des estimations précises du potentiel d'adhérence, il reste à résoudre le problème de l'estimation des points de friction à partir des signaux mesurables avec les capteurs équipant les véhicules de série.

Estimation des efforts pneus

Comme indiqué sur la Figure 6.9, cette section s'intéresse plus particulièrement à l'estimation des points de friction à partir des signaux directement disponibles sur les véhciules de série. Afin de ne pas inclure trop de grandeurs dans l'état du système, les différentes grandeurs permettant d'estimer les points de friction sont estimées suivant le schéma illustré sur la Figure 6.12. Tout d'abord, on es-Figure 6.12: Schéma détaillant l'estimation des points de fricton. time les états avec l'EKF. Ici, l'état est composé de la vitesse du véhicule, de la vitesse des roues, des dynamiques de tangage et des efforts longitudinaux. Une fois les états estimés, on utilise différents modèles (voir Figure 6.12) pour remonter jusqu'à l'adhérence et au glissement afin de calculer les points de friction. Etant donnée que cette méthode s'appuie sur l'utilisation d'un EKF, il vient le problème du réglage des matrices de covariances des bruits. Afin de simplifier cette étape, une méthode ne nécessitant aucun réglage de la part de l'utilisateur a été mise au point.

Réglage des matrices de covariances des bruits du filtre de Kalman

Le problème du réglage des matrices de covariances V et W est fréquemment rencontré lors de l'utilisation d'un filtre de Kalman. En effet, régler ces matrices avec des valeurs incorrectes peut réduire considérablement les performances du filtre de Kalman. Malheureusement, en pratique, les statistiques des bruits d'état et de sortie permettant de calculer les matrices de covariances ne sont pas accessibles. Ainsi, dans la majeure partie des cas, le filtre de Kalman est réglé lors d'une phase d'essai erreur. Durant cette phase, chacun des coefficients des matrices V et W est ajusté un par un. Ainsi pour des systèmes ayant un nombre élevé d'état, cette phase peut rapidement devenir fastidieuse pour l'utilisateur. Afin de palier à ce problème, certaines méthodes ont été développées dans la littérature afin d'obtenir des estimations des matrices de covariances V et W (voir l'article de synthèse [START_REF] Duník | Noise covariance matrices in state-space models: A survey and comparison of estimation methods-part i[END_REF]). Toutefois, la plupart de ces méthodes nécessitent le réglage par l'utilisateur de certains paramètres afin de pouvoir fonctionner correctement. Or, il n'est pas toujours évident de régler ces paramètres. Ainsi, une méthode permettant d'estimer des valeurs des matrices de covariances sans réglage est mise au point dans ce manuscrit. Cette méthode se base sur la théorie de l'identification à partir des méthodes des sous-espaces. Plus précisément, avec cette méthode, des estimations des matrices de covariances sont déduites en comparant une représentation d'état estimée et la représentation d'état discrète utilisée par le filtre de Kalman. L'estimation de la représentation d'état est obtenue en utilisant une méthode des sous-espaces basée sur l'algorithme N4SID. D'un point de vue plus formel, considérons la représentation d'état

x k+1 = A d x k + B d u k + w k , (6.2a
) Après avoir décrit le fonctionnement de la méthode, celle-ci est testée sur différents exemples de simulations. Les résultats obtenus avec ceuxci montrent que la méthode développée permet d'obtenir des estimations précises des matrices de covariances V et W . De plus, les estimations obtenues sont bien plus précises que celles fournies par d'autres méthodes plus standards rencontrées dans la littérature. Finalement, bien que cette méthode soit faite pour fonctionner avec des systèmes linéaires à temps invariant ayant des matrices de covariances V et W constantes, la méthode est aussi appliquée dans un cas ou le système considéré est non linéaire. Les expériences réalisées sur ce cas-ci ont montré que la méthode fournissait des résultats prometteurs. Plus précisément, les estimations obtenues permettent à l'EKF de reconstruire correctement les états.

y k = C d x k + v k . ( 6

Application des méthodes d'estimation des efforts pneus

Afin de tester les performances de l'EKF, celui est testé sur des données de simulation venant d'un modèle bicycle étendu et sur des données générées avec le logiciel de simulation dynamique VI-CarRealTime (VI-CRT). Dans le cas des données VI-CRT, les résultats fournis par N4SID et l'EKF réglé avec la méthode des sous-espaces sont illustrés sur les Figures 6.13-6.14. Comme on peut le voir sur la Figure 6.13, l'EKF permet de reconstruire avec précision le signal de l'adhérence contrairement à l'algorithme N4SID qui ne parvient pas à reproduire avec précision µ. Cependant, comme illustré sur la Figure 6.14, l'estimation du glissement fournit par EKF parvient à reproduire les dynamiques du vrai signal de glissement mais l'estimation est très bruitée. A l'inverse, l'algorithme N4SID fourni une estimation qui est très peu bruité.

Afin de vérifier si on ne peut pas avoir de meilleurs résultats avec l'EKF, un second réglage des matrices de covariances est proposé. Ce réglage consiste à modifier les coefficients de la matrice de covariances de bruit d'état qui correspondent aux efforts pneus et aux dynamiques

Conclusion et perspectives Conclusion

Dans ce manuscrit, une méthode est proposée pour estimer le potentiel d'adhérence d'un pneu dans des conditions normales de conduite et en utilisant uniquement les capteurs présents sur les véhicules de série. La méthode mise au point se scinde en deux parties. Tout d'abord le problème d'estimation du potentiel d'adhérence à partir de mesures des points de friction est étudié. En particulier, la méthode utilisée est une méthode de type MCMC. Cette méthode délivre des résultats promet- teurs en respectant la contrainte des conditions normales de conduite. Dans un second temps, un EKF appliqué à un modèle bicycle étendu avec un modèle de suspensions est utilisé afin d'obtenir des estimations des points de friction. Afin de régler l'EKF, une méthode se basant sur les méthodes des sous-espaces a été développée. Contrairement à la plupart des méthodes rencontrées dans la littérature, la méthode développée ne requiert aucun réglage de la part de l'utilisateur. Dans ce manuscrit, l'EKF est testé sur différente données de simulations. Tout d'abord sur des données générées avec un modèle bicycle étendu avec un modèle de suspensions. Puis, avec des données générées avec le logiciel de simulation dynamique VI-CRT. Les résultats obtenus avec l'EKF montre que celui-ci est capable de reconstruire avec précision le signal d'adhérence mais qu'il rencontre des difficultés pour estimer le signal de glissement.

Perspectives

Les travaux réalisés dans ce manuscrit peuvent faire l'objet de nombreuses perspectives. A court terme, il serait intéressant de valider l'estimation des points de friction avec des mesures venant du bus CAN d'un véhicule. Ce test est important car les mesures provenant du bus CAN sont les seules mesures accessibles sur les véhicules de série. De plus il serait aussi intéressant d'appliquer la méthode MCMC avec des mesures de points de friction correspondant à une route mouillée où le potentiel d'adhérence est plus faible. Ce test est important car il permettrait de déterminer si la méthode MCMC est capable d'observer une réduction du potentiel d'adhérence dans des conditions normales de conduite.

Par ailleurs, il apparait essentiel d'améliorer l'estimation du glissement. Une des solutions pourrait être d'utiliser des mesures du glissement. Cependant comme le montre la Figure 6.17, les mesures de glissement relevées sont parfois fortement bruitées. Sur la Figure 6.17

Figure 1 . 1 :

 11 Figure 1.1: Grip potential estimation method classification. The path marked in red indicates the type of estimation method elaborated in this manuscript.

Figure 1 . 2 :

 12 Figure 1.2: Example of tire friction curve. The blue points are the data points available under standard driving conditions.

Figure 1 . 3 :

 13 Figure 1.3: Example of SAT curve and tire friction curve.

Figure 1 . 4 :

 14 Figure 1.4: Illustration of the slip-slope concept. The red lines are the slip-slope of the different friction curves.

Figure 1 . 5 :

 15 Figure 1.5: Details of the slip-slope method.

Figure 1 . 7 :

 17 Figure 1.7: Illustration of L. Ray's method in a longitudinal case.

Figure 1 . 8 :

 18 Figure 1.8: Grip potential estimation master plan.

Figure 1 . 9 :

 19 Figure 1.9: Part of the master plan handled in the first problem.

Figure 1 . 10 :

 110 Figure 1.10: Grip potential estimation plan.

Figure 1 . 11 :

 111 Figure 1.11: Friction points estimation plan.

Figure 1 . 12 :

 112 Figure 1.12: Manuscript organization.

Figure 2 . 1 :

 21 Figure 2.1: Part of the master plan handled in this chapter. The addressed part corresponds to the part inside the dashed red rectangle. In this plan the shadow part corresponding to the estimation of friction point is replaced by friction points measurements.

Figure 2 . 2 :

 22 Figure 2.2: Example of tire friction curve. The blue points are the data points accessible under standard driving conditions.

Figure 2 . 3 :

 23 Figure 2.3: Grip potential estimation plan.

Figure 2 . 5 :

 25 Figure 2.5: Example of tire friction curve calculated with Pacejka's formula and geometrical interpretation of the model parameters.

Figure 2 . 7 :

 27 Figure 2.7: Histograms of the Pacejka's parameters estimated with the ML method (for µ ≤ 0.3). The red line represents the true value used to generate the friction measurements. The green line indicates the mean of the ML estimates. The magenta dash-dot line indicates the minimizer of the loss function. The green dash-dot line represents a Gaussian pdf fitted to the 10000 ML estimates.

Figure 2 . 8 :

 28 Figure 2.8: Histograms of the Pacejka's parameters estimated with the ML method (for µ ≤ 0.6). The red line represents the true value used to generate the friction measurements. The green line indicates the mean of the ML estimates. The magenta dash-dot line indicates the minimizer of the loss function. The green dash-dot line represents a Gaussian pdf fitted to the 10000 ML estimates.

Figure 2 . 9 :

 29 Figure 2.9: Histograms of the Pacejka's parameters estimated with the ML method (all data points selected). The red line represents the true value used to generate the friction measurements. The green line indicates the mean of the ML estimates. The magenta dash-dot line indicates the minimizer of the loss function. The green dash-dot line represents a Gaussian pdf fitted to the 10000 ML estimates.

Figure 2 . 10 :

 210 Figure 2.10: Histograms of the Pacejka's parameters estimated with the MCMC method (for µ ≤ 0.3). The red line represents the true value used to generate the friction measurements. The black dashed line indicates the mean of the MCMC estimates. The magenta dash-dot line indicates the minimizer of the loss function. The green line indicates the mean of the ML estimates. The blue line represents a Gaussian pdf fitted to the MCMC estimates. The green line represents a Gaussian pdf fitted to the 10000 ML estimates.

Figure 2 . 11 :

 211 Figure 2.11: Histograms of the Pacejka's parameters estimated with the MCMC method (for µ ≤ 0.6). The red line represents the true value used to generate the friction measurements. The black dashed line indicates the mean of the MCMC estimates. The magenta dash-dot line indicates the minimizer of the loss function. The green line indicates the mean of the ML estimates. The blue line represents a Gaussian pdf fitted to the MCMC estimates. The green line represents a Gaussian pdf fitted to the 10000 ML estimates.

15 . 5

 155 -69.4 199 -1.71 69.3 15.5 0.25 -1.13 3.26 -0.03 1.13 -69.4 -1.13 5.01 -14.5 0.13 -5.06 199 3.26 -14.5 41.8 -0.36 14.5 -1.71 -0.03 0.13 -0.36 0.003 -0.13 69.3 1.13 -5.06 14.5 -0.13 5.06

Figure 2 . 12 :Figure 2 . 13 :

 212213 Figure 2.12: Histograms of the Pacejka's parameters estimated with the MCMC method (all data points selected). The red line represents the true value used to generate the friction measurements. The black dashed indicates the mean of the MCMC estimates. The magenta dashdot line indicates the minimizer of the loss function. The green line indicates the mean of the ML estimates. The blue line represents a Gaussian pdf fitted to the MCMC estimates. The green line represents a Gaussian pdf fitted to the 10000 ML estimates.

Figure 2 . 14 :

 214 Figure 2.14: Results of ML and MCMC curve fitting on all the data available. The shadowed area represents the 300 estimates and the red curve the mean of these estimates. The dashed blue line indicates the ML estimate.

Figure 2 . 15 :

 215 Figure 2.15: Histograms of grip potential estimates (upper plot) and their corresponding slip ratios (lower plot). The red line represents the values obtained with the ML method. The dashed black vertical line represents the mean of the MCMC estimates. The blue line represents a Gaussian pdf fitted to the MCMC estimates.

Figure 2 . 16 :

 216 Figure 2.16: The upper plot represents the evolution of the estimated friction potential obtained with the ML and MCMC methods, respectively. The lower plot represents the evolution of the relative error of the estimated friction potential with the reference value.

Figure 2 . 17 :Figure 2 . 18 :Figure 2 . 19 :

 217218219 Figure2.17: The upper plot represents the evolution of the estimated slip ratio corresponding to the friction potential obtained with the ML and MCMC methods, respectively. The lower plot represents the relative errors with the reference of the estimated slip ratio corresponding to the friction potential.

Figure 2 . 20 :

 220 Figure 2.20: Results of ML and MCMC curve fitting (µ < 0.3). The red line represents the reference friction curve. The dashed-dot green line indicates the ML estimate. The dashed-dot magenta line indicates the MCMC estimate. The dashed blue line indicates the MCMC with physical prior estimate.

Figure 2 . 21 :

 221 Figure 2.21: Traceplot of θ M CM C (µ < 0.20) in the case where physical prior is added. The blue curve represents the value of the samples kept. The horizontal red line indicates the reference value.

Figure 2 . 22 :

 222 Figure 2.22: Traceplot of θ M CM C (µ < 0.30) in the case where physical prior is added. The blue curve represents the value of the samples kept. The horizontal red line indicates the reference value.

Figure 2 . 23 :

 223 Figure 2.23: Evolution of the histograms representing the friction potential estimated with the MCMC based estimation method with physical prior. The red line represents the reference value. The dashed black line represents the mean of the MCMC estimates.

Figure 2 . 24 :

 224 Figure 2.24: Evolution of the histograms representing the slip ratio corresponding to the friction potential estimated with the MCMC based estimation method with physical prior. The red line represents the reference value. The dashed black line represents the mean of the MCMC estimates.

Figure 2 . 25 :

 225 Figure 2.25: Measurement residuals. The upper plot represents the evolution of the residuals as a function of the slip ratio. The lower plot is an histogram of the residuals.

Figure 2 . 26 :

 226 Figure 2.26: ACF of the residuals. The left figure represents the ACF of the residuals corresponding to the friction measurements lower than 0.6. The right figure represents the ACF of all the residuals. In both figures, the dashed black lines represent the standard error of the ACF with a confidence interval of 95%.

Figure 2 . 27 :

 227 Figure 2.27: Master plan handled in this thesis. Contrary to what have been done in this chapter, the friction point measurements are not measured but estimated.

Figure 3 . 1 :

 31 Figure 3.1: Part of the master plan handled in this chapter. The addressed part corresponds to the part inside the dashed red rectangle.

Figure 3 . 2 :

 32 Figure 3.2: Concept of state observer.

Figure 3 . 3 :

 33 Figure 3.3: Double-track model with longitudinal and lateral tire forces.

Figure 3 . 4 :

 34 Figure 3.4: Single-track model with longitudinal and lateral tire forces.

  and 3.5 and described in Table3.3.

Figure 3 . 5 :

 35 Figure 3.5: Longitudinal forces acting on a vehicle during a travel.

Figure 3 . 6 :

 36 Figure 3.6: Illustration of pitch effect during a braking phase.

Figure 3 . 7 :

 37 Figure 3.7: Simplified passive suspension model.

. 11 )

 11 By introducing Eq. (3.7a), (3.9a) and (3.9b) to Eq. (3.11), the equation governing the pitch dynamics becomes

. 29 )

 29 The signals and physical parameters involved in Eq. (3.23)-(3.29) are summarized in Table

Remark 3 . 3 .Figure 3 . 8 :

 3338 Figure 3.8: Example of tire model.

Figure 3 . 9 :

 39 Figure 3.9: Friction data points estimation plan.

Figure 3 . 10 :

 310 Figure 3.10: Illustration of the discrete-time Kalman filter.

  d and C d plays a central role. Knowing A d and C d , the expected column space Γ f (A d , C d ) can be determined a priori. Then, as soon as f ≥ n x , the similarity transformation T between Γ f (A d , C d ) known a priori and the estimated column space Γf (A d , C d ) can be computed by noticing that

W

  = KR e K ⊤ = 0.0202 -0.0045 0.0149 -0.0198 -0.0045 0.0012 -0.0031 0.0041 0.0149 -0.0031 0.0111 -0.0148 -0.0198 0.0041 -0.0148 0.0199 , (4.34b) S = R e K ⊤ = 0.0183 -0.0045 0.0131 -0.0172 -0.0311 0.0063 -0.0234 0.0314 (4.34c)

Figure 4 . 1 :

 41 Figure 4.1: Pseudo Random Binary Sequence exciting the mass spring damper system.

Figure 4 . 2 :

 42 Figure 4.2: Noisy realization of the mass spring damper system output.

Figure 4 . 3 :

 43 Figure 4.3: Frequency responses of the sytem and the estimated model.

Figure 4 . 4 :

 44 Figure 4.4: Comparison of the real and estimated state variable trajectories (position (top) and velocity (bottom) of the mass, respectively).N4SID model learning.

Figure 4 . 5 :

 45 Figure 4.5: Comparison of the real and estimated state variable trajectories (position (top) and velocity (bottom) of the mass, respectively).Kalman filtering.

Figure 4 . 6 :

 46 Figure 4.6: Single-track model for lateral dynamics description.

Figure 4 . 7 :

 47 Figure 4.7: Lateral tire forces f i vs. side slip angles α i for different loads.

Figure 4 . 8 :Figure 4 . 9 :

 4849 Figure 4.8: Steering wheel angle generated by a pilot. This signal is a chirp generated by a pilot driving a real passenger vehicle.

Figure 4 . 10 :

 410 Figure 4.10: Noisy output signal samples.

Figure 4 . 11 :

 411 Figure 4.11: Comparison of the real and estimated state variable trajectories. N4SID model learning.

Figure 5 . 1 :

 51 Figure 5.1: Part of the master plan handled in this chapter. The addressed part corresponds to the part inside the dashed red rectangle.

Figure 5 . 2 :

 52 Figure 5.2: Macroscopic plan of the vehicle data simulator.

Figure 5 . 3 :

 53 Figure 5.3: Descriptive plan used to generate simulated data with the extended single-track model. The quantities framed in red are the user defined quantities.

Figure 5 . 4 .

 54 Then, the resulting data signals generated are the ones depicted in Figures 5.5-5.12.

Figure 5 . 4 :Figure 5 . 5 :Figure 5 . 6 :

 545556 Figure 5.4: Simulation of the known quantities.

Figure 5 . 7 :

 57 Figure 5.7: Pitch angle and pitch velocity simulated data.

5. 2 Figure 5 . 8 :

 258 Figure 5.8: Front and rear axle longitudinal tire forces simulated data.

Figure 5 . 9 :Figure 5 . 10 :Figure 5 . 11 :

 59510511 Figure 5.9: Front and rear axle normal loads simulated data.

(

  SNR) equal to 25 dB. The resulting noisy signals are represented in Figure 5.13.

Figure 5 . 12 :

 512 Figure 5.12: Front and rear effective tire radii simulated data.

Figure 5 . 13 :

 513 Figure 5.13: Noisy and noise-free signals of the measured quantities.

Figure 5 . 14 :

 514 Figure 5.14: Vehicle speed estimates. The upper plot shows the true signal (in black), its estimates with EKF (in green) and N4SID (in blue) and the measurements (in yellow). The bottom left figures illustrate the evolution of the residuals. The bottom right plots display the error histograms.

Figure 5 . 15 :

 515 Figure 5.15: Front wheel speed estimates. The upper plot shows the true signal (in black), its estimates with EKF (in green) and N4SID (in blue) and the measurements (in yellow). The bottom left figures illustrate the evolution of the residuals. The bottom right plots display the error histograms.

Figure 5 . 16 :

 516 Figure 5.16: Rear wheel speed estimates. The upper plot shows the true signal (in black), its estimates with EKF (in green) and N4SID (in blue) and the measurements (in yellow). The bottom left figures illustrate the evolution of the residuals. The bottom right plots display the error histograms.

Figure 5 . 17 :

 517 Figure 5.17: Front longitudinal tire force estimates. The upper plot shows the true signal (in black) and its estimates with EKF (in green) and N4SID (in blue). The bottom left figures illustrate the evolution of the residuals. The bottom right plots display the error histograms.

Figure 5 . 18 :

 518 Figure 5.18: Rear longitudinal tire force estimates. The upper plot shows the true signal (in black) and its estimates with EKF (in green) and N4SID (in blue). The bottom left figures illustrate the evolution of the residuals. The bottom right plots display the error histograms.

Figure 5 . 19 :

 519 Figure 5.19: Pitch angle estimates. The upper plot shows the true signal (in black) and its estimates with EKF (in green) and N4SID (in blue). The bottom left figures illustrate the evolution of the residuals. The bottom right plots display the error histograms.

Figure 5 . 20 :

 520 Figure 5.20: Pitch speed estimates. The upper plot shows the true signal (in black), its estimates with EKF (in green) and N4SID (in blue) and the measurements (in yellow). The bottom left figures illustrate the evolution of the residuals. The bottom right plots display the error histograms.

Figure 5 . 21 :

 521 Figure 5.21: Front normal load estimates. The upper plot shows the true signal (in black) and its estimates with EKF (in green) and N4SID (in blue). The bottom left figures illustrate the evolution of the residuals.The bottom right plots display the error histograms.

Figure 5 . 22 :

 522 Figure 5.22: Rear normal load estimates. The upper plot shows the true signal (in black) and its estimates with EKF (in green) and N4SID (in blue). The bottom left figures illustrate the evolution of the residuals.The bottom right plots display the error histograms.

Figure 5 . 23 :

 523 Figure 5.23: Front friction estimates. The upper plot shows the true signal (in black) and its estimates with EKF (in green) and N4SID (in blue). The bottom left figures illustrate the evolution of the residuals.The bottom right plots display the error histograms.

Figure 5 . 24 :

 524 Figure 5.24: Rear friction estimates. The upper plot shows the true signal (in black) and its estimates with EKF (in green) and N4SID (in blue). The bottom left figures illustrate the evolution of the residuals.The bottom right plots display the error histograms.

Figure 5 . 25 :

 525 Figure 5.25: Front slip ratio estimates. The upper plot shows the true signal (in black) and its estimates with EKF (in green) and N4SID (in blue). The bottom left figures illustrate the evolution of the residuals.The bottom right plots display the error histograms.

Figure 5 . 26 :

 526 Figure 5.26: Rear slip ratio estimates. The upper plot shows the true signal (in black) and its estimates with EKF (in green) and N4SID (in blue). The bottom left figures illustrate the evolution of the residuals.The bottom right plots display the error histograms.

Figure 5 . 27 :

 527 Figure 5.27: Front loaded radius estimates. The upper plot shows the true signal (in black) and its estimates with EKF (in green) and N4SID (in blue). The bottom left figures illustrate the evolution of the residuals.The bottom right plots display the error histograms.

Figure 5 . 28 :

 528 Figure 5.28: Rear loaded radius estimates. The upper plot shows the true signal (in black) and its estimates with EKF (in green) and N4SID (in blue). The bottom left figures illustrate the evolution of the residuals.The bottom right plots display the error histograms.

Figure 5 . 29 :

 529 Figure 5.29: Front effective tire radius estimates. The upper plot shows the true signal (in black) and its estimates with EKF (in green) and N4SID (in blue). The bottom left figures illustrate the evolution of the residuals. The bottom right plots display the error histograms.

Figure 5 . 30 :

 530 Figure 5.30: Rear effective tire radius estimates. The upper plot shows the true signal (in black) and its estimates with EKF (in green) and N4SID (in blue). The bottom left figures illustrate the evolution of the residuals. The bottom right plots display the error histograms.

  -5.18.

  -5.32 and Figures 5.33-5.34, respectively.

Figure 5 .

 5 Figure 5.31: F xf estimates. The upper plot shows the true signal (in black) and its estimates with W sub , V sub (in green) and W mix , V mix (in red). The bottom left figures illustrate the evolution of the residuals.The bottom right plots display the error histograms.

Figure 5 . 32 :

 532 Figure 5.32: F xr estimates. The upper plot shows the true signal (in black) and its estimates with W sub , V sub (in green) and W mix , V mix (in red). The bottom left figures illustrate the evolution of the residuals.The bottom right plots display the error histograms.

Figure 5 .

 5 Figure 5.33: µ f estimates. The upper plot shows the true signal (in black) and its estimates with W sub , V sub (in green) and W mix , V mix (in red). The bottom left figures illustrate the evolution of the residuals.The bottom right plots display the error histograms.

Figure 5 .

 5 Figure 5.34: µ r estimates. The upper plot shows the true signal (in black) and its estimates with W sub , V sub (in green) and W mix , V mix (in red). The bottom left figures illustrate the evolution of the residuals.The bottom right plots display the error histograms.

Figure 5 .

 5 [START_REF] Guiggiani | The Science of Vehicle Dynamics[END_REF] illustrates an example of simulation made with the CompactCar vehicle model of VI-CRT. However, although VI-CRT gives access to a large amount of data, it does not give access to the vehicle model equations. Therefore

Figure 5 .

 5 Figure 5.35: VI-CRT main interface.

Figure 5 . 36 :

 536 Figure 5.36: Example of VI-CRT simulation.

Figure 5 . 37 :

 537 Figure 5.37: Front and rear axle torques VI-CRT simulated data.

Figure 5 . 38 :Figure 5 . 39 : 20 Figure 5 . 40 :

 53853920540 Figure 5.38: Front and rear wheel speeds VI-CRT simulated data.

Figure 5 . 41 :

 541 Figure 5.41: Front and rear axle normal loads VI-CRT simulated data.

Figure 5 . 42 :

 542 Figure 5.42: Front and rear friction VI-CRT simulated data.

Figure 5 . 43 :

 543 Figure 5.43: Front and rear loaded tire radii VI-CRT simulated data.

Figure 5 . 44 :Figure 5 . 45 :

 544545 Figure 5.44: Front and rear effective tire radii VI-CRT simulated data.

Figure 5 . 7 )

 57 Figure 5.46: VI-CRT simulated data corresponding to the measurements used by the observer.

Figure 5 . 47 :

 547 Figure 5.47: Vehicle speed estimates. The upper plot shows the true signal (in black), its estimates with EKF (in green) and N4SID (in blue) and the measurements (in yellow). The bottom left figures illustrate the evolution of the residuals. The bottom right plots display the error histograms.

Figure 5 . 48 :

 548 Figure 5.48: Front wheel speed estimates. The upper plot shows the true signal (in black), its estimates with EKF (in green) and N4SID (in blue) and the measurements (in yellow). The bottom left figures illustrate the evolution of the residuals. The bottom right plots display the error histograms.

Figure 5 . 49 :

 549 Figure 5.49: Rear wheel speed estimates. The upper plot shows the true signal (in black), its estimates with EKF (in green) and N4SID (in blue) and the measurements (in yellow). The bottom left figures illustrate the evolution of the residuals. The bottom right plots display the error histograms.

Figure 5 . 50 :

 550 Figure 5.50: Front longitudinal tire force estimates. The upper plot shows the true signal (in black) and its estimates with EKF (in green) and N4SID (in blue). The bottom left figures illustrate the evolution of the residuals. The bottom right plots display the error histograms.

Figure 5 . 51 :

 551 Figure 5.51: Rear longitudinal tire force estimates. The upper plot shows the true signal (in black) and its estimates with EKF (in green) and N4SID (in blue). The bottom left figures illustrate the evolution of the residuals. The bottom right plots display the error histograms.

Figure 5 . 52 :

 552 Figure 5.52: Pitch angle estimates. The upper plot shows the true signal (in black) and its estimates with EKF (in green) and N4SID (in blue). The bottom left figures illustrate the evolution of the residuals. The bottom right plots display the error histograms.

Figure 5 . 53 :

 553 Figure 5.53: Pitch speed estimates. The upper plot shows the true signal (in black), its estimates with EKF (in green) and N4SID (in blue) and the measurements (in yellow). The bottom left figures illustrate the evolution of the residuals. The bottom right plots display the error histograms.

Figure 5 . 54 :

 554 Figure 5.54: Front normal load estimates. The upper plot shows the true signal (in black) and its estimates with EKF (in green) and N4SID (in blue). The bottom left figures illustrate the evolution of the residuals.The bottom right plots display the error histograms.

Figure 5 . 55 :

 555 Figure 5.55: Rear normal load estimates. The upper plot shows the true signal (in black) and its estimates with EKF (in green) and N4SID (in blue). The bottom left figures illustrate the evolution of the residuals.The bottom right plots display the error histograms.

Figure 5 . 56 :

 556 Figure 5.56: Front friction estimates. The upper plot shows the true signal (in black) and its estimates with EKF (in green) and N4SID (in blue). The bottom left figures illustrate the evolution of the residuals.The bottom right plots display the error histograms.

Figure 5 . 57 :

 557 Figure 5.57: Rear friction estimates. The upper plot shows the true signal (in black) and its estimates with EKF (in green) and N4SID (in blue). The bottom left figures illustrate the evolution of the residuals.The bottom right plots display the error histograms.

Figure 5 . 58 :

 558 Figure 5.58: Front slip ratio estimates. The upper plot shows the true signal (in black) and its estimates with EKF (in green) and N4SID (in blue). The bottom left figures illustrate the evolution of the residuals.The bottom right plots display the error histograms.

Figure 5 . 59 :

 559 Figure 5.59: Rear slip ratio estimates. The upper plot shows the true signal (in black) and its estimates with EKF (in green) and N4SID (in blue). The bottom left figures illustrate the evolution of the residuals.The bottom right plots display the error histograms.

Figure 5 . 60 :

 560 Figure 5.60: Front loaded radius estimates. The upper plot shows the true signal (in black) and its estimates with EKF (in green) and N4SID (in blue). The bottom left figures illustrate the evolution of the residuals.The bottom right plots display the error histograms.

Figure 5 . 61 :

 561 Figure 5.61: Rear loaded radius estimates. The upper plot shows the true signal (in black) and its estimates with EKF (in green) and N4SID (in blue). The bottom left figures illustrate the evolution of the residuals.The bottom right plots display the error histograms.

Figure 5 . 62 :

 562 Figure 5.62: Front effective tire radius estimates. The upper plot shows the true signal (in black) and its estimates with EKF (in green) and N4SID (in blue). The bottom left figures illustrate the evolution of the residuals. The bottom right plots display the error histograms.

Figure 5 . 63 :

 563 Figure 5.63: Rear effective tire radius estimates. The upper plot shows the true signal (in black) and its estimates with EKF (in green) and N4SID (in blue). The bottom left figures illustrate the evolution of the residuals. The bottom right plots display the error histograms.

  -5.66, frictions and slip ratios indicated in Figures 5.67

  -5.68 and Figures 5.69

  -5.70, respectively.

Figure 5 .

 5 Figure 5.64: κ estimates. The upper plot shows the true signal (in black), the measurements (in yellow) and its estimates with W sub , V sub (in green) and W mix , V mix (in red). The bottom left figures illustrate the evolution of the residuals. The bottom right plots display the error histograms.

Figure 5 .

 5 Figure 5.65: F xf estimates. The upper plot shows the true signal (in black) and its estimates with W sub , V sub (in green) and W mix , V mix (in red). The bottom left figures illustrate the evolution of the residuals.The bottom right plots display the error histograms.

Figure 5 . 66 :

 566 Figure 5.66: F xr estimates. The upper plot shows the true signal (in black) and its estimates with W sub , V sub (in green) and W mix , V mix (in red). The bottom left figures illustrate the evolution of the residuals.The bottom right plots display the error histograms.

Figure 5 .

 5 Figure 5.67: µ f estimates. The upper plot shows the true signal (in black) and its estimates with W sub , V sub (in green) and W mix , V mix (in red). The bottom left figures illustrate the evolution of the residuals.The bottom right plots display the error histograms.

Figure 5 .

 5 Figure 5.68: µ r estimates. The upper plot shows the true signal (in black) and its estimates with W sub , V sub (in green) and W mix , V mix (in red). The bottom left figures illustrate the evolution of the residuals.The bottom right plots display the error histograms.

Figure 5 .

 5 Figure 5.69: s f estimates. The upper plot shows the true signal (in black) and its estimates with W sub , V sub (in green) and W mix , V mix (in red). The bottom left figures illustrate the evolution of the residuals.The bottom right plots display the error histograms.

Figure 5 .

 5 Figure 5.70: s r estimates. The upper plot shows the true signal (in black) and its estimates with W sub , V sub (in green) and W mix , V mix (in red). The bottom left figures illustrate the evolution of the residuals.The bottom right plots display the error histograms.

  -5.73, frictions and slip ratios indicated in Figures 5.74-5.75 and Figures 5.76-5.77, respectively.

Figure 5 .

 5 Figure 5.71: κ estimates. The upper plot shows the true signal (in black), the measurements (in yellow) and its estimates with W sub , V sub (in green) and W mix2 , V mix2 (in red). The bottom left figures illustrate the evolution of the residuals. The bottom right plots display the error histograms.

Figure 5 . 72 :

 572 Figure 5.72: F xf estimates. The upper plot shows the true signal (in black) and its estimates with W sub , V sub (in green) and W mix2 , V mix2 (in red). The bottom left figures illustrate the evolution of the residuals.The bottom right plots display the error histograms.

Figure 5 . 73 :

 573 Figure 5.73: F xr estimates. The upper plot shows the true signal (in black) and its estimates with W sub , V sub (in green) and W mix2 , V mix2 (in red). The bottom left figures illustrate the evolution of the residuals.The bottom right plots display the error histograms.

Figure 5 .

 5 Figure 5.74: µ f estimates. The upper plot shows the true signal (in black) and its estimates with W sub , V sub (in green) and W mix2 , V mix2 (in red). The bottom left figures illustrate the evolution of the residuals.The bottom right plots display the error histograms.

Figure 5 .

 5 Figure 5.75: µ r estimates. The upper plot shows the true signal (in black) and its estimates with W sub , V sub (in green) and W mix2 , V mix2 (in red). The bottom left figures illustrate the evolution of the residuals.The bottom right plots display the error histograms.

Figure 5 .

 5 Figure 5.76: s f estimates. The upper plot shows the true signal (in black) and its estimates with W sub , V sub (in green) and W mix2 , V mix2 (in red). The bottom left figures illustrate the evolution of the residuals.The bottom right plots display the error histograms.

Figure 5 .

 5 Figure 5.77: s r estimates. The upper plot shows the true signal (in black) and its estimates with W sub , V sub (in green) and W mix2 , V mix2 (in red). The bottom left figures illustrate the evolution of the residuals.The bottom right plots display the error histograms.

Figure 5 . 78 :

 578 Figure 5.78: Friction curves. The black curve represents the true friction curve computed with VI-CRT signals. The red points represent the friction points computed with the estimated friction and slip ratio.

Figure 5 . 79 :

 579 Figure 5.79: Friction curves. The black curve represents the true friction curve computed with VI-CRT signals. The blue points represent the friction points computed with the estimated friction and true slip ratio.

Figure 6 . 1 :

 61 Figure 6.1: Longitudinal tire force and slip ratio measurements coming from sensors fitted on a tractor.

Figure 6 . 2 :

 62 Figure 6.2: Visualization of the temperature effect and road effect on the friction curve.

Figure 6 . 3 :

 63 Figure 6.3: Classification des méthodes d'estimation du potentiel d'adhérence. Le chemin représenté en rouge indique la type de méthode développée dans ce manuscrit.

Figure 6 . 4 :

 64 Figure 6.4: Exemple de courbe de friction. Les points bleus représentent le niveau d'ahérence consommé dans des conditions standards de conduite.

Figure 6 . 5 :

 65 Figure 6.5: Plan principal pour estimer le potentiel d'adhérence.

Figure 6 . 6 :

 66 Figure 6.6: Partie du plan principal étudiée dans le Chapitre 2. La partie considérée est celle encadrée dans le rectangle en pointillé rouge. Dans ce plan, la partie grisée est remplacée par des mesures des points de friction.

Figure 6 . 7 :Figure 6 . 8 :

 6768 Figure 6.7: Plan pour estimer le potentiel d'adhérence à partir de mesures des points de friction.

Figure 6 . 9 :

 69 Figure 6.9: Partie du plan principal étudiée dans le Chapitre 3. La partie considérée est celle encadrée dans le rectangle en pointillé rouge.

Figure 6 . 11 :

 611 Figure 6.11: Modèle de suspensions passif simplifié.

  .2b)Cette représentation d'état peut être réécrite de façon équivalentew k = x k+1 -A d x k -B d u k ,(6.3a)v k = y k -C d x k . (6.3b) Or, les matrices A d , B d et C d ainsi que u k et y k sont connues. Ainsi, l'idée de la méthode est d'utiliser la méthode N4SID pour reconstruire un état xk et d'utiliser cette estimation pour déduire des estimations vk et ŵk . Ces estimations sont par la suite utilisées afin de déduire des estimations V et Ŵ des matrices de covariances.

Figure 6 . 13 :

 613 Figure 6.13: Estimation de la friction avant. Le tracé effectué en haut montre le vrai signal (en noir), ses estimations avec EKF (en verte) et N4SID (en bleue). Les figures situées en bas à gauche illustrent les évolutions temporelles des résidus. Les figures situées en bas à droite montrent les histogrammes des résidus.

Figure 6 . 14 :Figure 6 . 15 :

 614615 Figure 6.14: Estimation du glissement avant. Le tracé effectué en haut montre le vrai signal (en noir), ses estimations avec EKF (en verte) et N4SID (en bleue). Les figures situées en bas à gauche illustrent les évolutions temporelles des résidus. Les figures situées en bas à droite montrent les histogrammes des résidus.

Figure 6 . 16 :

 616 Figure 6.16: Estimation du glissement avant. Le tracé effectué en haut montre le vrai signal (en noir), ses estimations avec EKF réglé avec W sub , V sub (en verte) et avec W mix , V mix (en rouge). Les figures situées en bas à gauche illustrent les évolutions temporelles des résidus. Les figures situées en bas à droite montrent les histogrammes des résidus.

Figure 6 . 17 :

 617 Figure 6.17: Mesures des efforts longitudinaux et du glissement provenant d'un capteur équippant un tracteur.

Table 2 .

 2 1: Signals and parameters used for defining the slip ratio.

	Symbol	Definition	Units
	ω	wheel speed	rad/s
	R rol	effective tire radius	m
	v x	longitudinal vehicle speed m/s

1 Figure 2.4: Example of curves representing the friction µ as a function of the slip ratio s.

Table 2 .

 2 2: Magic formula parameters.

	Symbol	Definition
	B	stiffness factor
	C	shape factor
	D	peak value
	E	curvature factor
	s v	vertical shift
	s

h horizontal shift

In a nutshell, by assuming that the user has access to

Table 2 .

 2 3: Intervals of Pacejka's parameters.

	Parameter Lower and upper bound
	B	[5, 30]
	C	[0.5, 2]
	D	[0.2, 2]
	E	[-2, 0]
	s h	[-0.05, 0.05]
	s v	[-0.3, 0.3]
	Under these conditions, the ML method provides 10000 estimates of the
	parameter vector θ. These estimates are represented in Figures 2.7, 2.8
	and 2.9, respectively. In addition, these Figures also include a Gaussian
	pdf fitted to the 10000 estimates. More precisely, the means and the
	standard deviations of the pdf are determined by computing the means
	and the standard deviations of the 10000 estimates obtained with each
	training set. Among the 10000 estimates, different estimated values are
	also extracted. Firstly, the mean of the 10000 estimates. Secondly, the
	estimated parameter vector θ ML giving the smallest residual norm V (θ)
	defined in Eq. (2.14). In the following, this second estimate is referred
	as "Minimizer of the loss function". These two estimated values are rep-
	resented in Figure 2.7, 2.8 and 2.9 and are summarized in Table 2.4.
	Furthermore, the ML method provides estimates of the measurement
	standard deviation (see Eq. (2.17)). These estimates are gathered and

Table 2 .

 2

	5: Measurement standard deviation estimated with the ML
	method. The relative error lines give the relative errors between the
	different estimates and the true value.
		σ
	True value	2.53e -2
	Friction limit : µ = 0.30
	ML	3.32e -2
	(Relative error)	30.9%
	Friction limit : µ = 0.60
	ML	2.96e -2
	(Relative error)	16.8%
	Without friction limit
	ML	2.53e -2
	(Relative error)	0%
	wider space. As a consequence, the initial covariance matrix of the proposal
	distribution is set as follows
	Σ (0) = diag (7, 0.43, 0.3, 0.3, 0.005, 0.01) .	(2.27)
	Under such practical conditions, the MCMC based estimation method
	provides the results given in Figures 2.10, 2.11 and 2.12 and sum up in

Table 2 .

 2 7. The MCMC estimated values gathered in Table

  65 -10463.68 800.89 16267.38 -48.61 952.37

	-10463.68	527.43	-47.37 -789.87	2.16	-42.26
	800.89	-47.37	6.69	60.58	-0.09	1.80
	16267.38	-789.87	60.58	1235.65 -3.68	72.16
	-48.61	2.16	-0.09	-3.68	0.01	-0.26
	952.37	-42.26	1.80	72.16	-0.26	5.04

Table 2 .

 2 7: MCMC estimates with simulated data. The relative error lines give to the relative errors between the different estimates and the true parameter values.

		B	C	D	E	s h	s v
	True value	15.4	1.60	8.71e -1	-1.09	0	0
	Friction limit : µ = 0.30
	MCMC	17.5	1.24	1.03	-1.03	2.41e -3 -1.67e -2
	(Relative error)	13.7% 22.1%	18.2%	5.97%	-	-
	Friction limit : µ = 0.60
	MCMC	16.3	1.21	1.21	-9.80e -1 4.20e -4	6.25e -3
	(Relative error)	6.03% 24.0%	38.7%	10.2%	-	-
	Friction limit : all points
	MCMC	16.6	1.62	8.09e -1 -9.31e -1 -2.14e -3 5.29e -2
	(Relative error)	7.70% 1.39%	7.19%	14.68%	-	-

Table 2 .

 2 8: ML and MCMC estimates with real data.

	Estimation method µ max	µ max	arg(µ max )	arg(µ max )
		mean standard deviation	mean	standard deviation
	ML	1.223	2.10e -6	3.970e -2	3.67e -15
	MCMC	1.222	2.65e -2	4.250e -2	4.23e -3

Table 3 .

 3 1: Signals and parameters used to define the double-track model.

	Symbol	Units	Definition
	v x , v y	m.s -1 Longitudinal/Lateral vehicle speed
	v g	m.s -1	Vehicle speed at the CoG
	ψ	rad.s -1	Yaw rate
	F xij	N	Longitudinal tire forces
	F yij	N	Lateral tire forces
	δ ij	rad	Steering angle
	β	rad	Sideslip angle
	t f , t r	m	Front/Rear track
	L f , L r	m	Distance between the vehicle CoG
			and the front/rear axle
	with
		X

Table 3 . 2 :

 32 Signals and parameters used to define the single-track model.

	Symbol	Units	Definition
	v

x , v y m.s -1 Longitudinal/lateral vehicle speed ψ rad.s -1 Yaw rate F xf , F xr N Front/Rear axle longitudinal tire forces F yf , F yr N Front/Rear axle lateral tire forces α f , α r rad Front/Rear slip angle δ rad Steering angle L f , L r m Distance between the vehicle CoG and the front/rear axle

Table 3 . 3 :

 33 Forces acting in a vehicle during a maneuver.

	Symbol Units	Definition
	F zf , F zr	N	Front/Rear normal load
	P	N	Vehicle weight
	R xf , R xr	N	Front/Rear rolling resistance
	F xf , F xr	N	Front/Rear longitudinal tire forces
	F aero	N	Aerodynamic drag force
	Remark 3.1. Although the Figure 3.4 represents also the lateral forces
	acting on a vehicle, in our case, we only focus on the longitudinal dynamics.
	For this reason, the equations modelling the lateral dynamics of a vehicle
	are omitted.

Table 3 .

 3

	Symbol	Units	Definition
	m	kg	Vehicle mass
	g	m.s -2	Acceleration of the gravity
	ρ a	kg.m -3	Air density
	S a	m 2	Frontal area
	C x	-	Aerodynamic drag coefficient
	f RR	-	Rolling resistance coefficient
	ω f , ω r	rad.s -1	Front/Rear wheel rotational speed
	T f , T r	N.m	Front/Rear driving and braking torque
	I wf , I wr	kg.m 2	Front/Rear wheel moment of inertia
	R 0	m	Free tire radius
	R f load , R rload	m	Front/Rear loaded tire radius
	K zz	N.m -1	Tire radial stiffness

4: Signals and parameters involved in the torque balance applied to the vehicle wheels.

Table 3 .

 3 5: Simplified suspension model parameters.

	Symbol	Units	Definition
	B f , B r N.s.m -1 Front/Rear suspension damping coefficient
	K f , K r	N.m -1	Front/Rear suspension stiffness
	h G	m	Height of the CoG
	κ	rad	Pitch angle of the sprung mass
	By assuming no variation of the center of gravity (CoG) height, the
	suspension model considered in Figure 3.7 gives us the following load
	transfer

  If the velocity v x is constant, the second member of the differential equation(3.16) becomes constant and an analytical solution can be computed. In this case, one solution of the Equation (3.16) is

	.17a)
	(3.17b)
	Remark 3.2.

Table 3 .

 3 6: Signals and parameters involved in the Pacejka's tire radius model.

	Symbol	Units	Definition
	R 0	m	Free tire radius ( = unload tire radius)
	F nomin	N	Nominal normal load
	K zz	N.m -1	Normal tire stiffness
	K zz0	N.m -1 Normal tire stiffness at the nominal normal load
	F z V 0	N m.s -1	Normal load Reference velocity (=	√	gR 0 )
	ω	rad.s -1	Wheel speed
	F x	N	Longitudinal tire force
	F y	N	Lateral tire force
	dp i	-	Normalized change in inflation pressure

  = |λ max (A d -KC d )| < 1,we have ∥ Ãd ∥ p < ρ p 0 which decreases towards 0 when p increases. Therefore, by assuming that p is taken sufficiently large, the quantity Ãd

	can be generated from Eq. (4.12) straightforwardly. Then, with standard
	recursions again, we have, for any user defined p ∈ N * ,
	xk = Ãd	p xk-p + Ω p ( Ãd , K)y k-p,p + Ω p ( Ãd , B d )u k-p,p .	(4.16)
	By assuming that the discrete time linear time invariant state space
	form (4.12) is strict minimum phase 1 , i.e.,
	ρ 0
			)

p xk-p can be neglected. Hence the state sequence approximation xk defined as follows

Table 4 .

 4

1, 4.2 and 4.3 gather the estimation results obtained by running 10 3 different realizations of the noise signals. More specifically, each table reports the average and the standard deviation of the main

Table 4 .

 4 1: Estimates of the elements of V .

		v11	v12	v22
	Theo. value	0.0176 -0.0267 0.0497
	Sample cov. avg. 0.0176 -0.0267 0.0497
		std. 0.0008 0.0012 0.0022
	ICM	avg. 0.0588 -0.0750 0.0888
		std. 0.0412 0.0149 0.0057
	DCM	avg. 0.031	-0.041	0.028
		std. 0.0041 0.0013	0.006
	CMM	avg. 0.023	-0.073	0.031
		std. 0.0087 0.0066 0.0092
	New meth. avg. 0.0198 -0.0272 0.0516
		std. 0.0011 0.0013 0.0024

Table 4 .

 4

			ŵ11	ŵ12	ŵ13	ŵ14	ŵ22	ŵ23	ŵ24	ŵ33	ŵ34	ŵ44
	Theo. value		0.0202	-0.0045	0.0149	-0.0198	0.0012	-0.0031	0.0041	0.0111	-0.0148	0.0199
	Sample cov. avg.	0.0202	-0.0045	0.0149	-0.0198	0.0012	-0.0031	0.0041	0.0111	-0.0148	0.0199
		std. 0.8886e-03 0.2054e-03 0.6563e-03 0.8750e-03 0.0509e-03 0.1487e-03 0.1967e-03 0.4877e-03 0.6520e-03 0.8725e-03
	ICM	avg.	0.0526	-0.0150	0.0355	-0.0454	0.0041	-0.0103	0.0133	0.0238	-0.0303	0.0386
		std.	0.0146	0.0079	0.0066	0.0067	0.0042	0.0036	0.0038	0.0029	0.0030	0.0030
	DCM	avg.	0.0113	-0.0058	0.0186	-0.0285	0.003	-0.0103	0.0033	0.0138	-0.0188	0.0223
		std.	0.004	0.0033	0.0068	0.0067	0.0039	0.0033	0.0034	0.0028	0.0031	0.0030
	CMM	avg.	0.0170	-0.0041	0.0124	-0.0234	0.0041	-0.0043	0.0041	0.0138	-0.0245	0.0224
		std.	0.0097	0.0082	0.0064	0.0062	0.0052	0.0028	0.0036	0.0031	0.0024	0.0030
	New meth. avg.	0.0196	-0.0041	0.0145	-0.0190	0.0015	-0.0026	0.0038	0.0117	-0.0147	0.0195
		std.	0.0017	0.0006	0.0011	0.0011	0.0004	0.0005	0.0007	0.0011	0.0011	0.0013

2: Estimates of the elements of W .

Table 4 .

 4 3: Estimates of the elements of S for the new method only (ICM, DCM and CMM do not give access to estimates for S indeed).

	ŝ11	ŝ12	ŝ13	ŝ14	ŝ21	ŝ22	ŝ23	ŝ24

Table 4 .

 4 4: Bicycle model variables.

	Symbol	Definition	Units
	f y f	front axle side force	N
	f yr	rear axle side force	N
	β	body slip angle at the CoG	rad
	ψ	yaw rate	rad/s
	ν	vehicle speed	km/h
	δ swa	steering wheel angle	rad
	α f	front side slip angle	rad
	α r	rear side slip angle	rad
	a y	lateral acceleration	m/s 2

Table 4 .

 4 5: Bicycle model parameters.

	Symbol	Definition	Units	Rough nominal Values
	v x	vehicle speed	m/s	25
	m	vehicle mass	Kg	1500
	c y f	front axle cornering stiffness	N/rad	200000
	c yr	rear axle cornering stiffness	N/rad	250000
	ℓ f	front axle distance to the CoG	m	1.5
	ℓ r	rear axle distance to the CoG	m	0.9
	I zz	inertia about the z axis	Kg/m 2	3000
	n	steering ratio	-	15
	σ α f	front axle relaxation length	m	0.7
	σ αr	rear axle relaxation length	m	0.2

  Traditionally, the noise covariance matrices are considered as tunable parameters to be adjusted by the user in order to make the Kalman filter work. However, accurately determining the noise covariance matrices can be complicated in practice, especially when multivariable systems are considered. Thus, this chapter introduces a new solution to estimate the covariance matrices of linear time invariant systems with stationary random disturbances. Obtaining such a solution is possible by using the well known subspace model identification theory. More specifically, our new solution consists in (i) using subspace model identification for extracting an accurate discrete time model from the available data sets, (ii) comparing the estimated state space model with the Kalman filter deterministic state space form known a priori in order to get reliable model misfit measurements, (iii) transforming this discrepancy measurements into covariance matrix estimates. By estimating the noise covariances matrices this way, no tuning of any parameters is re-

	tions.
	0.01
	0.008
	0.006
	0.004
	0.002
	0
	-0.002
	-0.004
	-0.006
	-0.008
	-0.01
	2	4	6	8	10	12	14	16	18	20
	0.25
	0.2
	0.15
	0.1
	0.05
	0
	-0.05
	-0.1
	-0.15
	-0.2
	-0.25
	2	4	6	8	10	12	14	16	18	20
	4000
	3000
	2000
	1000
	0
	-1000
	-2000
	-3000
	-4000
	-5000
	2	4	6	8	10	12	14	16	18	20
	6000
	4000
	2000
	0
	-2000
	-4000
	-6000
	2	4	6	8	10	12	14	16	18	20
	Figure 4.12: Comparison of the real and estimated state variable trajec-
		tories. Kalman filtering.

quired by the user. After describing the theory, the method has been tested with several simulation examples. These simulation results have shown that (i) this new solution can significantly outperforms standard solutions like the innovation correlation method, (ii) can yield reliable covariance matrix estimates when realistic Kalman filtering problems come into play. These results confirm that our new method can be efficiently combined with the Kalman filter equations to give access to accurate estimated state trajectories.

Table 5 .

 5 2: Parameters values used to generate simulated data.

	Symbol	Value	Units	Definition
	B f	2.50e4 N.s.m -1	Front suspension damping coefficient
	C x	0.4238	-	Aerodynamic drag coefficient
	f RR	0.0104	-	Rolling resistance coefficient
	g	9.806	m.s -2	Acceleration of the gravity
	h G	0.563	m	Height of the CoG
	I wf	2.6498	kg.m 2	Front wheel moment of inertia
	I wr	2.6498	kg.m 2	Rear wheel moment of inertia
	I y	1542.46	kg.m 2	Vehicle moment of inertia toward y-axis
	K f	1597.7	N.m -1	Front suspension stiffness
	K zz	2e5	N.m -1	Tire radial stiffness
	L	2.58	m	Wheelbase
	L f	1.02	m	Distance front axle to CoG
	L r	1.56	m	Distance rear axle to CoG
	m	1548.38	kg	Vehicle mass
	R 0	0.312	m	Free tire radius
	S a	2	m 2	Frontal area
	ρ a	1.22	kg.m -3	Air density
	F nomin	4000	N	Nominal normal load
	B ref f	6.1	-	Pacejka's effective tire radius model coefficient
	D ref f	0.45	-	Pacejka's effective tire radius model coefficient
	F ref f	0.01	-	Pacejka's effective tire radius model coefficient

Table 5 .

 5 3: GFC of the estimates obtained by setting the EKF with W sub , V sub and W mix and V mix , respectively.

	Symbol	Fit of the subspace EKF Fit of the mixed EKF
	v x (m.s -1 )	0.9919	0.9919
	ω f (rad.s -1 )	0.9854	0.9888
	ω r (rad.s -1 )	0.9855	0.9912
	F xf (N )	0.6606	0.8778
	F xr (N )	-8.0549	0.6051
	κ (rad)	0.8841	0.8959
	κ (rad.s -1 )	0.9721	0.9778
	F zf (N )	0.8846	0.8966
	F zr (N )	0.8846	0.8966
	µ f (-)	0.6437	0.8779
	µ r (-)	-7.4177	0.6132
	s f (-)	-0.0971	0.0538
	s r (-)	-Inf	-Inf
	R f load (m)	0.8846	0.8966
	R rload (m)	0.8846	0.8966
	R rol f (m)	0.8799	0.8961
	R rolr (m)	0.8927	0.8940

Table 5 .

 5 4: Means of the absolute errors obtained by setting the EKF with W sub , V sub and W mix and V mix , respectively.

	Symbol	Absolute error Absolute error Decreasing between
		means with	means with	absolute
		W sub and V sub W mix and V mix	errors means
	v x (m.s -1 )	6.14e -02	6.14e -02	0.1%
	ω f (rad.s -1 )	3.71e -01	2.82e -01	24.1%
	ω r (rad.s -1 )	3.65e -01	2.23e -01	38.8%
	F xf (N )	3.08e + 02	6.26e + 01	79.7%
	F xr (N )	2.95e + 02	9.13e + 00	96.9%
	κ (rad)	8.89e -04	5.60e -04	37.1%
	κ (rad.s -1 )	2.94e -04	2.31e -04	21.7%
	F zf (N )	2.27e + 01	1.41e + 01	37.7%
	F zr (N )	2.27e + 01	1.41e + 01	37.7%
	µ f (-)	3.35e -02	6.73e -03	79.9%
	µ r (-)	4.93e -02	1.56e -03	96.8%
	s f (-)	7.43e -03	6.42e -03	13.7%
	s r (-)	7.36e -03	5.57e -03	24.4%
	R f load (m)	5.68e -05	3.53e -05	37.7%
	R rload (m)	5.68e -05	3.53e -05	37.7%
	R rol f (m)	3.68e -06	2.28e -06	38.2%
	R rolr (m)	7.72e -06	4.90e -06	36.6%

Table 5 .

 5 5: Medians of the absolute errors obtained by setting the EKF with W sub , V sub and W mix and V mix , respectively.

	Symbol	Absolute error Absolute error Decreasing between
		medians with	medians with	absolute
		W sub and V sub W mix and V mix	errors medians
	v x (m.s -1 )	5.19e -02	5.20e -02	-0.2%
	ω f (rad.s -1 )	3.15e -01	2.37e -01	24.8%
	ω r (rad.s -1 )	3.08e -01	1.89e -01	38.8%
	F xf (N )	2.61e + 02	4.33e + 01	83.4%
	F xr (N )	2.50e + 02	5.54e + 00	97.8%
	κ (rad)	7.51e -04	4.07e -04	45.8%
	κ (rad.s -1 )	2.48e -04	1.94e -04	21.9%
	F zf (N )	1.92e + 01	1.02e + 01	46.5%
	F zr (N )	1.92e + 01	1.02e + 01	46.5%
	µ f (-)	2.83e -02	4.70e -03	83.4%
	µ r (-)	4.17e -02	9.17e -04	97.8%
	s f (-)	5.20e -03	4.49e -03	13.7%
	s r (-)	5.17e -03	3.93e -03	24.0%
	R f load (m)	4.79e -05	2.56e -05	46.5%
	R rload (m)	4.79e -05	2.56e -05	46.5%
	R rol f (m)	3.11e -06	1.66e -06	46.4%
	R rolr (m)	6.49e -06	3.46e -06	46.7%

Table 5 .

 5 6: Coefficients of W sub related to the longitudinal tire forces.

	Symbol Extended single-track generated data VI-CRT data
	F xf	6.84e3	1.02e3
	F xr	1.48e3	9.85e2
	Ḟxf	6.19e7	9.65e3
	Ḟxr	1.67e7	5.49e3

Table 5 .

 5 7: GFC of the estimates obtained by setting the EKF with W sub , V sub and W mix and V mix , respectively.

	Symbol	Fit of the subspace EKF Fit of the mixed EKF
	v x (m.s -1 )	0.9897	0.9898
	ω f (rad.s -1 )	0.9670	0.9762
	ω r (rad.s -1 )	0.9780	0.9811
	F xf (N )	0.9213	0.9290
	F xr (N )	-2.5881	-1.4744
	κ (rad)	0.8907	0.9074
	κ (rad.s -1 )	0.1236	0.6698
	F zf (N )	0.8410	0.8489
	F zr (N )	0.8674	0.8763
	µ f (-)	0.9240	0.9296
	µ r (-)	-2.7969	-1.5400
	s f (-)	0.4901	0.6103
	s r (-)	-27.4398	-23.4418
	R f load (m)	0.8410	0.8489
	R rload (m)	0.8674	0.8763
	R rol f (m)	0.8475	0.8549
	R rolr (m)	0.8567	0.8673

Table 5 .

 5 8: Absolute error means obtained by setting the EKF with W sub , V sub and W mix and V mix , respectively.

	Symbol	Absolute error Absolute error Decreasing between
		means with	means with	absolute
		W sub and V sub W mix and V mix	errors means
	v x (m.s -1 )	1.80e -03	1.94e -03	-8.3%
	ω f (rad.s -1 )	7.46e -03	4.62e -04	93.8%
	ω r (rad.s -1 )	7.97e -05	1.15e -03	-1346.9%
	F xf (N )	1.66e + 01	1.54e + 01	7.3%
	F xr (N )	6.39e -01	1.74e + 00	-171.9%
	κ (rad)	5.41e -05	7.09e -05	-31.1%
	κ (rad.s -1 )	5.88e -05	2.40e -07	99.6%
	F zf (N )	4.31e + 01	4.24e + 01	1.7%
	F zr (N )	-2.27e + 01	-2.19e + 01	3.3%
	µ f (-)	1.42e -03	1.35e -03	5.2%
	µ r (-)	4.37e -04	6.51e -04	-49.0%
	s f (-)	-1.75e -04	-2.90e -04	-65.2%
	s r (-)	-6.27e -05	-9.43e -05	-50.3%
	R f load (m)	-1.08e -04	-1.06e -04	1.7%
	R rload (m)	5.67e -05	5.48e -05	3.3%
	R rol f (m)	-7.20e -06	-7.06e -06	1.8%
	R rolr (m)	7.47e -06	7.22e -06	3.3%

Table 5 .

 5 9: Absolute error medians obtained by setting the EKF with W sub , V sub and W mix and V mix , respectively.

	Symbol	Absolute error Absolute error Decreasing between
		medians with	medians with	absolute
		W sub and V sub W mix and V mix	errors medians
	v x (m.s -1 )	1.52e -03	2.02e -03	-32.5%
	ω f (rad.s -1 )	1.00e -02	4.75e -03	52.6%
	ω r (rad.s -1 )	1.46e -03	3.23e -04	77.9%
	F xf (N )	1.96e + 01	1.70e + 01	13.4%
	F xr (N )	6.26e + 00	2.44e + 01	-289.0%
	κ (rad)	7.93e -05	1.00e -04	-26.4%
	κ (rad.s -1 )	5.93e -05	8.67e -06	85.4%
	F zf (N )	4.24e + 01	4.22e + 01	0.5%
	F zr (N )	-2.16e + 01	-2.11e + 01	2.7%
	µ f (-)	2.05e -03	1.86e -03	8.9%
	µ r (-)	1.05e -03	4.13e -03	-293.1%
	s f (-)	-1.60e -04	-2.46e -04	-53.6%
	s r (-)	-8.92e -05	-1.01e -04	-13.7%
	R f load (m)	-1.06e -04	-1.06e -04	0.5%
	R rload (m)	5.41e -05	5.27e -05	2.7%
	R rol f (m)	-7.00e -06	-6.99e -06	0.1%
	R rolr (m)	7.06e -06	6.70e -06	5.1%
	called V mix2 is given by

Table 5 .

 5 10: GFC of the estimates obtained by setting the EKF with W sub , V sub and W mix2 and V mix2 , respectively.

	Symbol	Fit of the subspace EKF Fit of the mixed EKF
	v x (m.s -1 )	0.9897	0.9920
	ω f (rad.s -1 )	0.9670	0.9751
	ω r (rad.s -1 )	0.9780	0.9813
	F xf (N )	0.9213	0.9267
	F xr (N )	-2.5881	-1.4808
	κ (rad)	0.8907	0.9063
	κ (rad.s -1 )	0.1236	0.9061
	F zf (N )	0.8410	0.8491
	F zr (N )	0.8674	0.8766
	µ f (-)	0.9240	0.9273
	µ r (-)	-2.7969	-1.5464
	s f (-)	0.4901	0.6044
	s r (-)	-27.4398	-22.1690
	R f load (m)	0.8410	0.8491
	R rload (m)	0.8674	0.8766
	R rol f (m)	0.8475	0.8551
	R rolr (m)	0.8567	0.8675

Table 5 .

 5 11: Absolute error means obtained by setting the EKF with W sub , V sub and W mix2 and V mix2 , respectively.

	Symbol	Absolute error	Absolute error	Decreasing between
		means with	means with	absolute
		W sub and V sub W mix2 and V mix2	errors means
	v x (m.s -1 )	4.02e -02	3.09e -02	23.0%
	ω f (rad.s -1 )	3.66e -01	3.24e -01	11.6%
	ω r (rad.s -1 )	2.87e -01	2.43e -01	15.3%
	F xf (N )	7.58e + 01	6.85e + 01	9.6%
	F xr (N )	6.02e + 01	4.24e + 01	29.6%
	κ (rad)	5.56e -04	4.78e -04	14.0%
	κ (rad.s -1 )	2.46e -03	1.65e -04	93.3%
	F zf (N )	4.68e + 01	4.50e + 01	3.8%
	F zr (N )	3.68e + 01	3.25e + 01	11.8%
	µ f (-)	8.67e -03	8.06e -03	7.1%
	µ r (-)	9.76e -03	6.79e -03	30.5%
	s f (-)	4.54e -03	3.95e -03	13.0%
	s r (-)	3.75e -03	3.07e -03	17.9%
	R f load (m)	1.17e -04	1.12e -04	3.8%
	R rload (m)	9.21e -05	8.13e -05	11.8%
	R rol f (m)	7.81e -06	7.52e -06	3.7%
	R rolr (m)	1.21e -05	1.06e -05	12.1%

Table 5 .

 5 [START_REF] Betancourt | A conceptual introduction to hamiltonian monte carlo[END_REF]: Absolute error medians obtained by setting the EKF with W sub , V sub and W mix2 and V mix2 , respectively.

	Symbol	Absolute error	Absolute error	Decreasing between
		medians with	medians with	absolute
		W sub and V sub W mix2 and V mix2	errors medians
	v x (m.s -1 )	3.38e -02	2.56e -02	24.3%
	ω f (rad.s -1 )	2.51e -01	2.66e -01	-5.9%
	ω r (rad.s -1 )	2.44e -01	2.04e -01	16.6%
	F xf (N )	4.76e + 01	4.23e + 01	11.1%
	F xr (N )	5.05e + 01	3.12e + 01	38.3%
	κ (rad)	4.84e -04	4.31e -04	11.0%
	κ (rad.s -1 )	1.82e -03	9.52e -05	94.8%
	F zf (N )	4.31e + 01	4.26e + 01	1.2%
	F zr (N )	3.08e + 01	2.54e + 01	17.4%
	µ f (-)	5.49e -03	4.98e -03	9.2%
	µ r (-)	8.35e -03	5.23e -03	37.4%
	s f (-)	3.13e -03	3.13e -03	-0.1%
	s r (-)	3.09e -03	2.54e -03	17.7%
	R f load (m)	1.08e -04	1.07e -04	1.2%
	R rload (m)	7.69e -05	6.36e -05	17.4%
	R rol f (m)	7.12e -06	7.07e -06	0.8%
	R rolr (m)	1.01e -05	8.18e -06	19.3%

1.2 State of the art

Notice that, most of the time, a car drives straightforward[START_REF] Singh | Literature review and fundamental approaches for vehicle and tire state estimation[END_REF].

As shown in[START_REF] Andrieu | An introduction to MCMC for machine learning[END_REF], "in order to obtain best results out of this class of algorithms", the user must "incorporate as much domain specific knowledge as possible into their design".

In the following, the MCMC method refers to the method described in Algorithm 1, i.e., the adaptive MCMC initialized with ML estimate.

This assumption is satisfied[START_REF] Kailath | Linear Estimation[END_REF] when K is a Kalman filter gain.

This setting of the matrix V have been tried in the previous case where the data signals are generated with the extended single-track model. However, in this case, this setting does not improve the accuracy of the results.

Remerciements