Alexandre Barachant

Stéphane Bonnet

Marco Congedo

Christian Jutten

Keywords:

come

Multi-class Brain Computer Interface Classification by Riemannian Geometry

Alexandre Barachant, Stéphane Bonnet, Marco Congedo, and Christian Jutten

Abstract-This paper presents a new classification framework for Brain Computer Interface (BCI) based on motor imagery. This framework involves the concept of Riemannian geometry in the manifold of covariance matrices. The main idea is to use spatial covariance matrices as EEG signal descriptors and to rely on Riemannian geometry to directly classify these matrices using the topology of the manifold of Symmetric and Positive Definite (SPD) matrices. This framework allows to extract the spatial information contained in EEG signals without using spatial filtering. Two methods are proposed and compared with a reference method (multi-class Common Spatial Pattern (CSP) and Linear Discriminant Analysis (LDA)) on the multi-class dataset IIa from the BCI competition IV. The first method, named Minimum Distance to Riemanian Mean (MDRM), is an implementation of the Minimum Distance to Mean (MDM) classification algorithm using Riemannian distance and Riemannian mean. This simple method shows comparable results with the reference method. The second method, named Tangent Space LDA (TSLDA), maps the covariance matrices onto the Riemannian tangent space where matrices can be vectorized and treated as Euclidean objects. Then, a variable selection procedure is applied in order to decrease dimensionality and a classification by LDA is performed. This latter method outperforms the reference method increasing the mean classification accuracy from 65.1% to 70.2%.

I. INTRODUCTION

A Brain Computer Interface (BCI) aims at translating brain signals into commands. On one hand, it will be useful for people with severe motor impairments in order to restore communication and movement. On the other hand, it could be a new interface for healthy people, e.g., in gaming applications.

In the literature several electro-physiological sources of BCI control have been investigated. This paper focuses on Motor Imagery (MI), although the method we propose herein may be applied to other types of BCIs. In MI-based EEG-BCI, the standard operation mode consists of a cue-based (or synchronous) calibration stage followed by an asynchronous operation mode [START_REF] Gouy-Pailler | Nonstationary brain source separation for multiclass motor imagery[END_REF]. The first stage is often devoted to build subject-dependent frequency and spatial filters discriminating between EEG datasets corresponding to two different classes of MI [START_REF] Blankertz | Optimizing spatial filters for robust EEG Single-Trial analysis[END_REF]. Such spatial filters perform linear combination of the EEG signals in order to create new signals with maximal variance in one condition and minimal variance in the other condition [START_REF] Blankertz | Optimizing spatial filters for robust EEG Single-Trial analysis[END_REF]. Once these spatial filters have been designed, the (log-)variance of the spatially-filtered signals are used A. Barachant as features by a supervised classification algorithm. Linear Discriminant Analysis (LDA) is often used to perform this processing [START_REF] Blankertz | Single-trial analysis and classification of ERP components -a tutorial[END_REF].

In this work, we propose an entirely different approach by exploiting directly the covariance structure of the data as the feature of interest. Since the covariance matrices contain the spatial information embedded in EEG signal, this approach aims to merge the spatial filtering and the classification procedure into one unique step. This idea has already been proposed in [START_REF] Farquhar | A linear feature space for simultaneous learning of spatiospectral filters in BCI[END_REF] where the spatial filters and the classifier are trained together in an optimization process. However, the particular structure of covariance matrices, which belong to the Riemannian manifold of the symmetric and positive definite (SPD) matrices, has to be treated carefully. In this context, the Riemannian geometry provides a rich framework to manipulate these matrices. In addition, the properties of the manifold of SPD matrices allow to have explicit formulas for the main operations in the Riemannian manifold, which lead to an easier implementation of algorithms.

A similar approach has also been recently developed in the context of EEG for sleep states detection [START_REF] Li | EEG signal classification based on a riemannian distance measure[END_REF] that focused on the exploitation of the frequential information through the estimation of the power spectral covariance matrices. Working with covariance matrices in their manifold has proven useful in other fields such as radar image processing [START_REF] Barbaresco | Innovative tools for radar signal processing based on cartan's geometry of SPD matrices & information geometry[END_REF], diffusion tensor imaging [START_REF] Fletcher | Principal geodesic analysis on symmetric spaces: Statistics of diffusion tensors[END_REF] and image processing [START_REF] Tuzel | Pedestrian detection via classification on Riemannian manifolds[END_REF]. This methodology has also been used for pedestrian detection in image processing [START_REF] Tuzel | Pedestrian detection via classification on Riemannian manifolds[END_REF].

The first idea of our paper is to manipulate EEG spatial covariance matrices in their native space and to make use of the Riemannian distance between them in this space. Spatial filtering is no more necessary, since the spatial information is embedded in the spatial covariance matrices. Therefore all distance-based classification algorithms can be extended by applying this new Riemannian distance.

The second idea we propose is to use the concept of tangent space to apply more sophisticated classification algorithms. The operation named tangent space mapping sends the covariance matrices, belonging to a manifold, into an Euclidean space where they can be treated as vectors. This mapping operation allows the use of state-of-the-art classifiers within the Riemannian framework.

The proposed methods are developed based on an earlier work in [START_REF] Barachant | Riemannian geometry applied to BCI classification[END_REF]. Here we address the multi-class case and the high dimensionality of the tangent space This article is organized as follows. First we discuss about the use of the covariance matrices in motor imagery-based BCI. Then, in Section III we introduce the basic concepts of Riemannian geometry. Section IV presents the two proposed classification methods and Section V is dedicated to the related results on a BCI dataset. Finally, we conclude on the proposed approach in Section VI.

II. ON THE USE OF COVARIANCE MATRICES IN BCI

By definition, motor imagery consists in the imagination of a limb movement like hand or foot. The topographical representation and band power change of such mental tasks are well-known since different body parts are represented in different area of the motor cortex (roughly, right hand under C 3 electrode, left hand under C 4 , foot under C z , etc.). The associated frequential phenomena are known as Event-Related Desynchronisation (ERD) and Event-Related Synchronisation (ERS) [START_REF] Pfurtscheller | Event-related EEG/MEG synchronization and desynchronization: basic principles[END_REF].

The EEG recording of the brain activity is usually cut into trials, i.e., short-time windows. During the calibration operation mode, the trials are supervised, in the sense that the MI mental task or class is known. In the functional operation mode, single-trial classification must be performed to recognize the different mental tasks.

Spatial covariance matrices are used in this context. For instance, the popular Common Spatial Pattern (CSP) algorithm is entirely based on their estimation, from which spatial filters are derived to enhance class separability [START_REF] Blankertz | Optimizing spatial filters for robust EEG Single-Trial analysis[END_REF].

In the following, let x t ∈ R n denote the EEG signal vector at a specific time point t, with n denoting the number of recording channels. The spatial covariance matrix is formally defined by

Σ = E (x t -E {x t }) (x t -E {x t })
T , where superscript T denotes vector or matrix transposition and E {.} the expected value.

In BCI we consider short-time segments of EEG signal, or trials, in the form of a matrix X i = [x t+Ti . . . x t+Ti+Ts-1 ] ∈ R n×Ts which corresponds to the i-th trial of imagined movement started at time t = T i . Here T s denotes the number of sampled time points in each trial. Here and hereafter we suppose that each channel measurement has been previously time-centered by a high-pass filtering operation. For the ith trial, the spatial covariance matrix is estimated using the Sample Covariance Matrix (SCM) P i ∈ R n×n such as :

P i = 1 T s -1 X i X T i . ( 1
)
The SCM is known to be an unbiased estimator of the covariance matrix Σ provided that the number of observations T s is much larger than the number of variables n [START_REF] Smith | Covariance, subspace, and intrinsic Cramèr-Rao bounds[END_REF].

III. RIEMANNIAN GEOMETRY

In this section we briefly define the concepts and tools of Riemannian geometry that are needed in the proposed approach. The main point is that spatial SCM belongs to a particular manifold and classification can be directly achieved in this Riemannian space, as it will be shown in Section IV.

A. Notations

Denote by S(n) = {S ∈ M (n), S T = S} the space of all n × n symmetric matrices in the space of square real matrices M (n) and P (n) = {P ∈ S(n), u T Pu > 0, ∀ u ∈ R n } the set of all n × n symmetric positive-definite (SPD) matrices. A SPD matrix is always diagonalizable with strictly real positive eigenvalues. Finally, denote by Gl(n) the set of all n × n invertible matrices in M (n).

The Frobenius norm of a matrix is given by : A 2 F = Tr AA T = |A ij | 2 where Tr(.) denotes the trace operator.

For a vector a, the L 2 norm is denoted by a 2 .

For SPD matrices in P (n), the exponential matrix of P is obtained using the eigenvalue decomposition of P :

P = U diag (σ 1 , . . . , σ n ) U T ,
where σ 1 > σ 2 > . . . σ n > 0 are the eigenvalues and U the matrix of eigenvectors of P. It reads :

exp(P) = U diag (exp(σ 1 ), . . . , exp(σ n )) U T .
The inverse operation is the logarithm of a SPD matrix :

log(P) = U diag (log(σ 1 ), . . . , log(σ n )) U T .
We also have the following properties :

• ∀ P ∈ P (n), det(P) > 0 • ∀ P ∈ P (n), P -1 ∈ P (n) • ∀ (P 1 , P 2 ) ∈ P (n) 2 , P 1 P 2 ∈ P (n) • ∀ P ∈ P (n), log(P) ∈ S(n) • ∀ S ∈ S(n), exp(S) ∈ P (n)
Finally, notation P 1/2 defines a symmetric matrix A that fulfils the relation AA = P.

B. Riemannian natural metric

The space of SPD matrices P (n) is a differentiable Riemannian manifold M [START_REF] Förstner | A metric for covariance matrices[END_REF]. The derivatives at a matrix P on the manifold lies in a vector space T P , which is the tangent space at that point. The tangent space is lying in the space S(n). The manifold and the tangent space are m = n(n+1)/2 dimensional.

Each tangent space has an inner product , P that varies smoothly from point to point over the manifold. The natural metric on the manifold of SPD matrices is defined by the local inner product :

S 1 , S 2 P = Tr S 1 P -1 S 2 P -1 .
(

) 2
The inner product induces a norm for the tangent vectors on the tangent space, such that, S 2 P = S, S P = Tr SP -1 SP -1 . We note that, at Identity matrix, such norm simplifies into the Frobenius norm, i.e. S, S I = S 2 F . Working in the manifold M of the SPD matrices has several advantages. Due to the properties of SPD matrices, there exists explicit formulae for the major concepts of Riemannian geometry. Therefore they can be applied easily in the context of signal processing.

C. Riemannian Geodesic distance

Let Γ(t) : [0, 1] → P (n) be any (differentiable) path from Γ(0) = P 1 to Γ(1) = P 2 . The length of Γ(t) is given by

L(Γ(t)) = 1 0 Γ(t) Γ(t) dt, (3)
with the norm defined previously. The minimum length curve connecting two points on the manifold is called the geodesic, and the Riemannian distance between the two points is given by the length of this curve. The natural metric (2) induces the geodesic distance [START_REF] Moakher | A differential geometric approach to the geometric mean of symmetric Positive-Definite matrices[END_REF] δ R (P 1 , P 2 ) = log P -1

1 P 2 F = n i=1 log 2 λ i 1/2 , (4)
where λ i , i = 1 . . . n are the real eigenvalues of P -1 1 P 2 . The main properties of the Riemannian geodesic distance are as follows :

• δ R (P 2 , P 1 ) = δ R (P 1 , P 2 ) • δ R (P -1 1 , P -1 2 ) = δ R (P 1 , P 2 ) • δ R (W T P 1 W, W T P 2 W) = δ R (P 1 , P 2 ) ∀ W ∈ Gl(n).
The third property is very important, since it implies that the space of SPD matrices is invariant by projection. It allows us to manipulate such space with tools like PCA without incidence on the distance.

D. Exponential map

For each point P ∈ P (n), we can thus define a tangent space composed by the set of tangent vectors at P. Each tangent vector S i can be seen as the derivative at t = 0 of the geodesic Γ i (t) between P and the exponential mapping P i = Exp P (S i ), defined as

Exp P (S i ) = P i = P 1 2 exp P -1 2 S i P -1 2 P 1 2 . ( 5
)
The inverse mapping is given by the logarithmic mapping defined as

Log P (P i ) = S i = P 1 2 log P -1 2 P i P -1 2 P 1 2 . ( 6
)
This procedure is described geometrically in Fig. 1.

Figure 1. Tangent space at point P, S i a tangent vector at P and Γ i (t) the geodesic between P and P i .

Equivalent definitions of the Riemannian distance are :

δ R (P, P i ) = Log P (P i ) P = S i P = upper P -1 2 Log P (P i ) P -1 2 2 = s i 2 . ( 7
)
where the upper(.) operator consists in keeping the upper triangular part of a symmetric matrix and vectorizing it by applying unity weight for diagonal elements and √ 2 weight for out-of-diagonal elements [START_REF] Tuzel | Pedestrian detection via classification on Riemannian manifolds[END_REF]. Here s i is the m-dimensional vector upper P -1 2 Log P (P i ) P -1 2 of the normalized tangent space.

Under some conditions over P and the P i [8], (7) leads to an approximation in terms of distance between the tangent space and the Riemannian manifold, such as ∀i, j δ R (P i , P j ) ≈ s i -s j 2 .

(

) 8
In order to verify these conditions the P i must be locally distributed into the manifold, i.e., located into a small part of the manifold and P must be the mean of the P i .

E. Mean of SPD Matrices

1) Euclidean mean: Using the Euclidean distance on M(n), δ E (P 1 , P 2 ) = P 1 -P 2 F , it is possible to define the Euclidean mean of I ≥ 1 SPD matrices by

A (P 1 , . . . , P I ) = arg min P∈P (n) I i=1 δ 2 E (P, P i ) = 1 I I i=1 P i . (9)
Such mean is also referred to as the arithmetic mean in the literature [START_REF] Moakher | A differential geometric approach to the geometric mean of symmetric Positive-Definite matrices[END_REF] and it has been used ubiquitously so far.

2) Riemannian mean: Using Riemannian geodesic distance (4) the Riemannian mean of I ≥ 1 SPD matrices is given by G (P 1 , . . . , P I ) = arg min

P∈P (n) I i=1 δ 2 R (P, P i ) , ( 10
)
This mean is also referred to as the geometric mean. Indeed, if we consider

1 × 1 SPD matrices {x i > 0} 1≤i≤I , this defini- tion gives G (x 1 , . . . , x I ) = arg min x>0 I i=1 log 2 (x i /x) = I √ x 1 . . . x I .
For a manifold of non-positive sectional curvature like P (n), such local minimum exists and is unique [START_REF] Karcher | Riemannian center of mass and mollifier smoothing[END_REF]. However, there is no closed-form expression to compute the mean and optimisation algorithms must be employed. An efficient iterative algorithm to compute the Riemannian mean of I SPD matrices is given in [START_REF] Fletcher | Principal geodesic analysis on symmetric spaces: Statistics of diffusion tensors[END_REF].

IV. MI-BASED BCI CLASSIFICATION USING RIEMANNIAN FRAMEWORK

During the calibration operation mode, a set of annotated trials is obtained for each class of motor imagery. Suppose we have a set of trials obtained for the k-th condition X i , i ∈ I (k) with I (k) the set of indices of the trials corresponding to the k-th condition. The SCM of each trial can be computed using [START_REF] Gouy-Pailler | Nonstationary brain source separation for multiclass motor imagery[END_REF]. Since SCMs are SPD matrices, they do belong to the manifold M. Let P denote the SCM of the unlabel trial X. It is possible, using the results in Section III, to derive several efficient classification algorithms to estimate the unknown label of trial X. G . The condition with minimum-distance is affected to the unknown trial X. For each new trial X, this procedure amounts to estimating K intra-class covariance matrices and then compute K eigenvalue-decompositions of (4).

Algorithm 1 Minimum Distance to Riemannian Mean

Input: a set of trials X i of K different known classes. Input: X an EEG trial of unknown class. Input : I (k) the set of indices of the trials corresponding to the k-th condition. Output: k the estimated class of test trial X.

1: Compute SCMs of X i to obtain P i , (1). 2: Compute SCM of X to obtain P, (1).

3: for k = 1 to K do 4: Many popular and efficient classification algorithms [START_REF] Lotte | A review of classification algorithms for EEG-based brain-computer interfaces[END_REF] (LDA, SVM, Neural Network) are based on projections into hyperplanes. As such, they cannot be implemented directly in Riemannian manifold. The implementation of more sophisticated classification algorithms can be readily obtained with the use of the tangent space located at the geometric mean of the whole set of trials : P G = G (P i , i = 1 . . . I). Each SCM P i is then mapped into this tangent space, to yields the set of m = n(n + 1)/2 dimensional vectors :

P (k) G = G P i , i ∈ I (k) , ( 10
s i = upper P -1 2 G Log P G (P i ) P -1 2 G . ( 11
)
Algorithm 2 Tangent Space Mapping

Input: a set of I SPD matrices P i ∈ P (n) Output: a set of I vectors s i . 1: P G = G (P i , i = 1 . . . I) {Compute Riemannian mean of the whole set, (10). } 2: for i = 1 to I do 3:

s i = upper P -1 2 G Log P G (P i )P -1 2 G
4: end for 5: return s i 1) Variable selection: Since the tangent space is a m = n(n + 1)/2 dimensional space, the number of dimensions may now exceed the number of trials of each class. Regularised classification algorithms are usually needed in this case to tackle this problem [START_REF] Blankertz | Single-trial analysis and classification of ERP components -a tutorial[END_REF], [START_REF] Friedman | Regularization paths for generalized linear models via coordinate descent[END_REF].

An alternative solution that we have explored is to consider a variable selection procedure to decrease the space dimensionality. We use a one-way ANOVA to select the most discriminant variables. To account for any correlation between variables, the vectors s i are first orthogonalized using a Singular Value Decomposition (SVD) :

S = UΛV T ,
where S = [s 1 . . . s I ] ∈ R m×I , U ∈ R m×m and V ∈ R I×I two orthogonal matrices and Λ ∈ R m×I a diagonal matrix composed by the singular values of S. The tangent space S is projected using the orthogonal basis U :

S o = U T S,
to obtain the orthogonalized tangent space S o where all variables are uncorrelated. This operation is known as Principal Component Analysis (PCA).

Then, for each variable of S o , a one-way ANOVA is applied. All variables are ranked according to their p-values and a weighted False Discovery Rate (FDR) [START_REF] Benjamini | Multiple hypotheses testing with weights[END_REF] is applied to select automatically the minimal number of variables. The singular values are used as weight for the FDR procedure in order to give priority to variables that best represent the structure of the data.

2) Tangent Space Linear Discriminant Analysis (TSLDA): After the variable selection, we can apply any kind of classical classification algorithms since the tangent space is Euclidean. In this article, the classify function of Matlab (The MathWorks, inc, Natick MA) will be used. This function performs the linear discriminant analysis of Fisher between each pair of class. Then, the K(K-1) 2 classifiers are combined to obtain the classification results.

V. RESULTS

A. Description of data

We analyse the dataset IIa from the BCI competition IV [START_REF] Leeb | BCI competition 2008 graz data set b[END_REF]. It contains EEG data from 9 subjects who perform four kinds of motor imagery (right hand, left hand, foot and tongue imagined movements). EEG signals are recorded using 22 electrodes. For each subject, a training set and a test set are available. In this paper both sets are concatenated and performances are assessed by means of a cross-validation procedure. A total of 576 trials per subject are available (144 trials per class). The same pre-processing step is applied on the whole dataset. The EEG signals are bandpass filtered by a 5-th order Butterworth filter in the 8 -30 Hz frequency band. The time interval is restricted to the time segment comprised between 0.5 and 2.5s after the cue instructing the user to perform the mental task.

B. Results

1) Minimum Distance to Riemannian Mean:

Let us observe the results of the MDRM (Algorithm 1) algorithm. For illustration purpose we consider the case of the subject S1, who achieves fair performances. In order to show interpretable figures, we reduce the dataset to a two-class case, considering the data from the right hand and left hand MI realisations. Fig. 2 shows the distance of each trial from the two classrelated mean covariance matrices P G . On the abscissa, we have the distance from the the right hand mean, i.e., δ R (P

G , P i ). On the ordinate, we have the distance from the left hand mean, i.e., δ R (P (2) G , P i ). The dashed line represent the decision border whereupon δ R (P (1)

G , P i ) = δ R (P (2) G , P i ).
In this instance, only a small percentage of trials are misclassified and the use of Riemannian distance to mean is effective. To understand the benefit of the Riemannian framework, the same process is applied using both Euclidean distance and Euclidean mean, as expressed in [START_REF] Barachant | Riemannian geometry applied to BCI classification[END_REF]. Fig. 3 illustrates the obtained results. The trials are distributed around the dashed line regardless of their class membership. This fact implies that the relevant information about class membership is not accessible using the Euclidean distance and the associated mean. 2) Tangent space: We need to check the approximation, in terms of distance between the Riemannian manifold and the tangent space, as stated by [START_REF] Tuzel | Pedestrian detection via classification on Riemannian manifolds[END_REF]. It can be empirically verified looking at the normalized pairwise errors given by

ij = |δ R (P i , P j ) -s i -s j 2 | δ R (P i , P j ) . ( 12
)
The average and standard deviation of ij across the 9 subjects are 2% and 0.6%, respectively. We conclude that approximation ( 8) is verified in these data.

Then we apply Algorithm 2 to map the data in the Tangent space. It is interesting to observe the distribution of each class across trials. Fig. 4 illustrates these distributions for the most discriminant variable, selected using an ANOVA over the m variables of the tangent space. Clearly, all classes are not separated equally well. We can distinguish two groups. On one side, the trials corresponding to both hands and on the other side the trials corresponding to the foot and tongue. This phenomenon is due to physiological reasons, as it was discussed in [START_REF] Pfurtscheller | Mu rhythm (de)synchronization and EEG single-trial classification of different motor imagery tasks[END_REF]. Also, for this variable, the trials corresponding to the left hand are almost fully separated from those corresponding to the tongue movement. Obviously, increasing the number of variables will bring supplementary information useful to treat the other cases. The CSP algorithm is widely used in BCI. If we observe the distributions of features corresponding to the best CSP spatial filter (Fig. 5) and selected using the methods described in [START_REF] Grosse-Wentrup | Multiclass common spatial patterns and information theoretic feature extraction[END_REF], we found the same pattern of distributions as observed in Fig. 4. This is not surprising, according to our work in [START_REF] Barachant | Common spatial pattern revisited by riemannian geometry[END_REF]; we have shown that the CSP feature space can be viewed as a filtered and approximated representation of the Riemannian manifold. In a similar way, the tangent space is also an approximated representation of the Riemannian manifold as explained in III-D. In this example the feature corresponding to the best CSP filter is able to better discriminate the different classes as compared to the best variable in tangent space. However, we need to keep in mind that CSP is a supervised technique and the spatial filters are build using the knowledge of class label in order to maximize the discrimination of the different classes. On the other hand, the tangent space mapping is a fully unsupervised operation. In this sense, the variables obtained with the tangent space mapping are not subject to over-fitting. In this article, we use a supervised variable selection in order to reveal the variables of interest. However, this step is not mandatory and depends on the used classification algorithms.

The number of electrodes is n = 22, thus the tangent space is of dimension m = 253, which is not much lower than the number of trials (576). Albeit possible it is not efficient to use standard classification algorithms like LDA. This is due to the difficulty to compute an unbiased estimation of the features covariance matrix. For this reason we perform a selection of discriminant variables according to the procedure described in Section IV-B1.

After the orthogonalization of the tangent space, the ANOVA sorts the variables in terms of discriminability between the classes. Fig. 6 shows the p-value evolution along the sorted variables. One can observe a quick increase of the pvalue, suggesting that only a small subset of variables is useful. The minimal number of variables is set using the weighted FDR algorithm with q = 0.05 (in the FDR procedure q is the expected proportion of false rejections with respect to all rejections). For this user, the number of selected variables is 10. 3) Classification: Fig. 7 shows the evolution of the classification error rate against the number of selected variables. The classifier is a multi-class implementation of LDA, thus it is very sensitive to over-fitting: the error on the training set decreases to zero while the test set error rate increases after 10 variables. The minimum classification error rate is given for 10 variables. The weighted FDR procedure allows to select the optimal number of variables. Retaining the best 10 variables the error rate amounts to 19.5% for this user. A graphical representation of the confusion matrix is given in Fig. 8. As expected, most ambiguities are found between the two pairs of hands and foot vs. tongue.

C. Classification results

We have evaluated the performance of the proposed methods by means of classification accuracy. For all methods, we applied a 30-fold cross-validation procedure [START_REF] Graimann | Brain-Computer Interfaces: Revolutionizing Human-Computer Interaction[END_REF]. Performance is measured in terms of classification accuracy, thus the higher the score the better the performance. Since there are four classes, the chance level is 25%.

The two proposed Riemannian methods are compared with the classical approach consisting in a spatial filtering by multiclass CSP [START_REF] Gouy-Pailler | Nonstationary brain source separation for multiclass motor imagery[END_REF], [START_REF] Grosse-Wentrup | Multiclass common spatial patterns and information theoretic feature extraction[END_REF] followed by a LDA classification [START_REF] Blankertz | Optimizing spatial filters for robust EEG Single-Trial analysis[END_REF] on the feature space composed by the log variance of the spatially filtered signals. The multi-class CSP used in this paper is an implementation of [START_REF] Gouy-Pailler | Nonstationary brain source separation for multiclass motor imagery[END_REF] that performed approximate joint diagonalization of all the class-related mean covariance matrices. This algorithm is more efficient than the standard CSP used in a one-versus-rest manner for the multi-class case. The number of CSP spatial filter is set to 8 as proposed in the reference papers [START_REF] Gouy-Pailler | Nonstationary brain source separation for multiclass motor imagery[END_REF], [START_REF] Grosse-Wentrup | Multiclass common spatial patterns and information theoretic feature extraction[END_REF]. For a more fair comparison, we also presented the results for the CSP method (denoted CSP*+LDA) whereby the optimal number of filters was selected for each subject according to a weighted FDR criterion that is similar to our method.

No parameters need to be set for both Riemannian methods. For the TSLDA, the procedure is fully automatic since the weighted FDR choose the optimal number of variables.

The results are shown in Table I. First of all, the subject specific selection of the number of spatial filters did not improve significantly the results of the reference method. However, this selection makes the procedure fully automatic. The MDRM method offers a slightly worse result as compared to the reference methods. This is not surprising since this method does not perform any kind of de-noising. More interesting results are given by the TSLDA method. For all subjects except subject S1, this method outperforms the reference methods. It is also interesting to see that the improvement is more important for the difficult cases. The average number of selected variables is small (∼ 13.5) in comparison to the total number of variables.

A more sophisticated analysis can be done using the confusion matrices given in Tables II-III-IV-V. It appears that the reference methods have more difficulties to treat the left hand and foot classes. We observe the same behaviour with the MDRM method. On the other hand, the TSLDA approach seems not affected by this problem, displaying around 70% correct classification rate for all classes. In summary, the improvements brought upon by TSLDA are mainly due to a better handling of critical cases, resulting in a 7% improvement classification for both left hand and foot classes.

Finally, Table VI gives the results in terms of kappa values, as it was done for the BCI competition IV. The described methods are applied in the conditions of the Competition in order to compare performances. The MDRM method achieved a mean performance of 0.52 which ranks this method to the second place of the competition. The TSLDA method achieved a mean performance of 0.567 which is close to the score of the winner. Moreover, the scores are obtained without exploiting the frequential information, contrary to the winner. Higher performances can be reached using an optimisation of frequency filters.

VI. CONCLUSION

This article presents a new framework to classify single trials in BCI applications based on motor imagery. It relies on the covariances matrices as descriptors of EEG signal. Two methods based on Riemannian geometry have been proposed to classify directly the covariance matrices. The first method, named MDRM, is an implementation of the Minimum Distance to Mean (MDM) algorithm using Riemannian Distance and Riemannian Mean instead of the classical Euclidean equivalent. This algorithm is simple, effective and offers results close to those obtained with the reference method (consisting in a multi-class CSP followed by a LDA). The second method, named TSLDA, is based on a powerful Riemannian concept: the tangent space. The covariance matrices are mapped onto a higher dimensional space where they can be vectorized and treated as Euclidean objects. In this tangent space a variable selection procedure is applied in order to decrease dimensionality and a classification by LDA is performed. Significant better results have been achieved with the TSLDA method compared to the reference one. This improvement is mainly due to a better handling of difficult (noisy) cases where there is a strong ambiguity between the class memberships.

Both methods do not need spatial filtering anymore. In the case of the MDRM classification method, the only necessary signal processing steps are a temporal filtering, an estimation of the mean covariance matrices and a comparison of Riemannian distances. No parameters need to be set and the obtained results appear satisfactory. In addition, the tangent space mapping is an unsupervised operation that extracts spatial information with a quality comparable to the stateof-the-art CSP. This result may lead to the development of efficient unsupervised training algorithms in BCI applications.

Yet, there are several issues that need to be investigated in details. On one hand, spatial filtering by CSP can be viewed as a way to reduce the dimensionality of the problem and to give a criterion based on eigenvalues in order to select the best components. On the other hand, the tangent space mapping allows gathering more spatial information by making a higher dimensional space. This high dimensionality of the tangent space could bring difficulties for some algorithms, causing overfitting or bias in statistical estimations. For this reason, the unconditional replacement of spatial filtering by the tangent space mapping is not recommended in all cases. However, several solutions are available to address this problem. In this article we have chosen to use a variable selection procedure but a regularization strategy may be efficient as well. Future works should try to reduce the computational cost of these methods. Since a computation of a Riemannian distance involves a matrix inversion and a matrix diagonalization, the computational time dramatically increases with the number of electrodes. This effect could be reduced using an optimized implementation of eigenvalue decomposition or by using other metrics approximating the Riemannian metric [START_REF] Arsigny | Geometric means in a novel vector space structure on symmetric Positive-Definite matrices[END_REF].

Finally, it is difficult to interpret the results physiologically, unlike the observation of the spatial patterns in CSP, which gives information about the brain areas involved in the MI. Even if selected variables in the tangent space correspond to geodesics in the Riemann manifold, these geodesics can not be trivially used to visualize physiological principles, i.e., electrode positions or weights. We will investigate this point in a future work.

Beyond these few limitations the proposed framework appears very promising. This work opens the door to a new family of BCI algorithms wherein investigations have just started.

A. Classification in the Riemannian manifold 1 )

 1 Minimum Distance to Riemannian Mean (MDRM): The intra-class covariance matrices for each condition P (k) G , where k = [1 : K] denotes the class indices, can be computed using the results in Section III-E. A very simple supervised classification algorithm consists in computing the Riemannian distance between the unknown SCM and every intra-class covariance matrix P (k)

  ). 5: end for 6: k = arg min k δ R (P, P (k) G ), (4). 7: return k B. Classification in the Riemannian tangent space

Figure 2 .

 2 Figure 2. Riemannian distance to Riemannian mean for 2 class-related mean covariance matrices -case of right hand vs. left hand

Figure 3 .

 3 Figure 3. Euclidean distance to Euclidean mean for 2 class-related mean covariance matrices -case of right hand vs. left hand

Figure 4 .

 4 Figure 4. Class distribution across trials for the most discriminant variable in Tangent space

Figure 5 .

 5 Figure 5. Class distribution across trials for the best CSP spatial filter

Figure 6 .

 6 Figure 6. Evolution of the p-value (straight line) and number of selected variables (dashed line)

Figure 7 .

 7 Figure 7. Classification error vs. number of selected variables

Figure 8 .

 8 Figure 8. Confusion matrix in 30-fold cross validation for the subject S1

  and S. Bonnet are with the CEA-LETI, MINATEC Campus, 38054 Grenoble, France.

M. Congedo and C. Jutten are with the GIPSA-lab, CNRS, Grenoble University, 38402 Saint Martin d'Hères, France. The third and fourth author have been partially supported by ANR (French Natioanal research Agency) project OpenViBE2 and RoBIK.

Table VI RESULTS

 VI OF THE BCI COMPETITION IN KAPPA VALUES

		mean kappa	S1	S2	S3	S4	S5	S6	S7	S8	S9
	1 st	0.57	0.68	0.42 0.75	0.48	0.40	0.27	0.77	0.75	0.61
	TSLDA	0.567	0.74	0.38 0.72	0.50	0.26	0.34	0.69	0.71	0.76
	MDRM	0.52	0.75	0.37 0.66	0.53	0.29	0.27	0.56	0.58	0.68
	2 nd	0.52	0.69	0.34 0.71	0.44	0.16	0.21	0.66	0.73	0.69
	3 rd	0.31	0.38	0.18 0.48	0.33	0.07	0.14	0.29	0.49	0.44


Moustapha Doumiati

Olivier Sename
email: olivier.sename@gipsa-lab.grenoble-inp.fr

John Jairo

Martinez Molina

Luc Dugard

Peter Gaspar
email: gaspar@sztaki.hu

Zoltan Szabo

Jozsef Bokor

J Martinez

Vehicle yaw control via coordinated use of steering/braking systems

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Under critical driving circumstances, it is difficult for a driver to stabilize the vehicle, and dangerous accidents could happen. Safety of ground vehicles requires the improvement of yaw stability by active control. The basic idea is to assist the vehicle handling to be close to a linear vehicle handling characteristic that is familiar to the driver, and to restrain the vehicle lateral dynamics to be within the stable handling region in aggressive maneuvers.

By applying additional steering angle to front/rear wheels, Active Steering control (AS) can enhance steerability and lateral vehicle dynamics by directly regulating tire slip angles and thus the lateral tire forces [START_REF] Guven | Robust two degree-of-freedom vehicle steering controller design[END_REF], [START_REF] Rajamani | Vehicle dynamics and control[END_REF]). Unfortunately, AS becomes less effective when the vehicle reaches the handling limits due to the tire saturations. In order to maintain vehicle stability under critical driving conditions, an alternative approach utilizes differential braking forces between the left and the right sides of the vehicle to produce the required corrective yaw moment [START_REF] Rajamani | Vehicle dynamics and control[END_REF], [START_REF] Park | H ∞ direct yaw-moment control with brakes for robust performance and stability of vehicles[END_REF], [START_REF] Boada | Fuzzy-logic applied to yaw moment control for vehicle stability[END_REF]). This braking based technique is referred to as Direct Yaw moment Control (DY C). Note that DY C is effective in both vehicle linear/nonlinear region, however, it is not desirable in normal conditions because of the direct influence of the control action on the longitudinal vehicle dynamics (i.e it causes the vehicle to slow down significantly). Consequently, these two control techniques are optimized individually in specific handling regions, and the maximum benefit can be gained through the coordinated use of both methods of corrective yaw motion generation in the control strategy. In this field, several solutions have been proposed in recent years [START_REF] He | Coordination of active steering, driveline, and braking for integrated vehicle dynamics control[END_REF] and [START_REF] Yang | Coordinated control of AFS and DYC for vehicle handling and stability based on optimal guaranteed cost theory[END_REF]), and this topic is still an object of intense research activities from both industrial and academic sides. It is worth noting that some studies also involve suspension actuators to improve handling performances, like in (Poussot-Vassal et al. [2010]), but this is out of the scope of this paper. This work deals with the design of a new vehicle chassis control scheme that integrates and coordinates braking and front steering in order to enhance vehicle handling and yaw stability. The proposed control design is worked out on the basis of a 2-DOF (Degree-Of-Freedom) linear planar vehicle model, and it is formulated as a Linear Parameter Varying (LP V ) controller structure. The exogenous scheduling control parameter is function of the sideslip angle dynamics, and it activates the braking control only in critical maneuvers. The controller is synthesized within the Linear Matrix Inequalities (LM I) framework, while guaranting H ∞ performances. This new control methodology constitutes the major contribution of this paper. It aims at simplifying engineer design, reducing development time in making actuators cooperate, guaranteeing robustness properties with respect to model uncertainties, and internal stability. Besides, the proposed strategy does not involve any online optimization process, that means it could be suitable for real-time applications.

The response of the vehicle with the proposed control scheme has been evaluated with MATLAB computer simulations using a complex full vehicle nonlinear model subject to various driving situations. The obtained results confirm the effectiveness of the proposed control.

The rest of this paper is structured as follows. Section II first introduces the global control scheme, and then develops and synthesizes the vehicle dynamical stability controller. Performance analysis is done in Section III through time simulations. Conclusions and discussions are given in Section IV.

Paper notations:

Throughout the paper, the following notations will be

CONTROL SYSTEM DESIGN

The control system is shown in the block diagram of figure 1. The main goal of the proposed control system is to make the actual yaw rate, ψ, to follow the desired yaw rate, ψd . The reference model is adopted to keep the vehicle within the linear region that is familiar to the driver. Another purpose is to limit the vehicle sideslip angle, β, to be within an acceptable region to prevent vehicle spin.

As seen in figure 1, the proposed controller responds to yaw rate error, and its outputs are the active steer angle and the wheel brake torque. In the following, each block of the control scheme will be described in details.

Reference model

To obtain a desired vehicle response, it is necessary that the yaw rate follows its target value ψd . In this work, ψd is calculated using a 2-DOF classical linear bicycle model depicted in figure 2. The use of this model is explained in details in [START_REF] Dugoff | An analysis of tire traction properties and their influence on vehicle dynamic performance[END_REF]). Roll, pitch, and longitudinal dynamics are neglected to simplify the vehicle dynamics. The equations governing the lateral and yaw motions in this model can be expressed as:

• Equation of lateral motion:

mv β -ψ = C f δ -β -l f ψ v + Cr -β + lr ψ v (1)
• Equation of yaw motion:

Iz ψ = C f -l f ψ v -β -δ l f + Cr β -lr ψ v lr (2)
Consequently, ψd is function of the driver steering wheel angle, δ d , and the vehicle speed, v.

Since the lateral acceleration of the vehicle cannot exceed the maximum friction coefficient µ, the desired yaw rate must be limited by the following value [START_REF] Rajamani | Vehicle dynamics and control[END_REF]):

ψd,max ≤ |µg/v| (3)

Yaw controller

The function of the yaw control system is to restore the yaw rate of the vehicle as close as possible to the nominal motion expected by the driver. Consequently, the controller is designed so that the vehicle follows the reference yaw rate through driving the tracking error between the actual and desired yaw rate to zero. The controller structure is hierarchical and designed in 2 layers:

a) The upper controller provides the active steer angle, δ * , and the corrective yaw moment, M * z , needed to track the target yaw rate, and thus ensures the vehicle handling. Besides, it is supposed to reject disturbances that may affect the lateral motion of the vehicle. Note that when the vehicle is in its linear region, the controller ensures steerability, and only steering is used to follow the desired response. However, when the vehicle reaches the handling limits, steering and braking act together to maintain the vehicle stability. b) The lower controller converts the stabilizing yaw moment generated by the upper controller into effective braking torque, and it decides which wheel must be braked to counteract the undesired yaw motion.

a) Upper controller: LP V /H ∞ controller design

The 2-DOF linear bicycle model is used as the control model. Although the bicycle model is relatively simple, it includes the important features of the lateral vehicle dynamics. Taking into account the controller structure, this model is extended to include: 1) the direct yaw moment input M * z , and 2) a lateral disturbance force F dy . In the following, the extended linear bicycle model given in system (4) is used for synthesis. To synthesize the upper control, the H ∞ control performance is used. For more information about the H ∞ theory, reader can refer to [START_REF] Skogestad | Multivariable feedback control, analysis and design[END_REF]. In the following, the generalized synthesis plant together with the performance weighting functions, called here g , is presented and illustrated in figure 3. Dynamics of the actuators are neglected during the controller design process. g is given thereafter:

ψ β =    - l 2 f C f + l 2 r C r I z v l r C r -l f C f I z 1 + l r C r -l f C f mv 2 - C f + C r mv    ψ β +    l f C f I z C f mv    δ * + 1 I z 0 M * z + 0 1 mv F dy (
Σ g : ẋ z y = A B 1 B 2 C 1 D 11 D 12 C 2 0 0 x w u (5)
where x includes the state variables of the system and of the weighting functions, w = [ ψd , F dy ] T is the exogenous input vector, u = [δ * , M * z ] T represents the control input signals, y = [e ψ ] T is the measurement vector, and z = [z 1 , z 2 , z 3 ] T contains the weighted controlled outputs which have to be as small as possible.

Weighting functions:

In order to formulate the standard structure for the H ∞ controller defined in figure 3, the weighting functions W 1 , W 2 , and W 3 are defined to characterize the performance objectives and the actuator limitations (actuator descriptions are given in Subsection 2.4):

• W 1 weights the yaw rate error signal:

W 1 = 7 s/10πf 1 + 1 s/2πf 1 + 1 . ( 6
)
where f 1 = 8 Hz is the cut-off frequency of the low pass filter. W 1 is shaped in order to reduce the yaw rate error. • W 2 weights the braking control signal according to a scheduling parameter ρ:

W 2 = ρ × 10 -3 s/70f 2 + 1 s/700f 2 + 1 , (7)
where f 2 = 10 Hz is the braking actuator cut-off frequency. W 2 is linearly parameterized by the considered varying parameter ρ(.), where ρ ∈ ρ ≤ ρ ≤ ρ (with ρ = 0.1 and ρ = 10). Then, when ρ = ρ, the braking input is penalized, on the contrary, when ρ = ρ, the braking control signal is relaxed. • W 3 weights the steering control signal:

W 3 = G 0 δ (s/2πf 3 + 1)(s/2πf 4 + 1) (s/α2πf 4 + 1) 2 G 0 δ = (∆ f /α2πf 4 + 1) 2 (∆ f /2πf 3 + 1)(∆ f /2πf 4 + 1) ∆ f = 2π(f 4 + f 3 )/2, ( 8
)
where f 3 = 1 Hz is the lower limit of the actuator intervention, and f 4 = 10 Hz is the steering actuator bandwidth. This filter is designed in order to allow the steering control to act only in [f 3 , f 4 ] frequency range, where the driver cannot act. Thus, it is ensured that the steering action is comfortable for the driver. This filter design is inspired from [START_REF] Guven | Robust two degree-of-freedom vehicle steering controller design[END_REF].

These weighting functions are recalled thereafter in the sensitivity function plots as upper bounds limits, 1/W i (see figure 5).

Problem resolution: LM

I based LP V /H ∞ :
The H ∞ problem consists in finding a stabilizing controller, named S(ρ) (see figure 3), scheduled by ρ, of the form:

S(ρ) : ẋc u = A c (ρ) B c (ρ) C c (ρ) 0 x c y (9)
that minimizes the H ∞ norm of the closed-loop LP V system formed by the interconnection of equations ( 5) and ( 9), where u = [δ , M * z ] T and y = [e ψ ] T . As the parameter dependency, ρ, enters in a linear way in the system definition, the LP V polytopic approach solution consists in finding a common Lyapunov function at each vertex ρ, ρ , of the polytope defined by system (5). Thus, a LM I problem has to be solved, minimizing the attenuation level γ [START_REF] Scherer | Multiobjective output-feedback control via LMI optimization[END_REF]). Using Yalmip/Sedumi solver [START_REF] Lofberg | YALMIP: a toolbox for modeling and optimization in MATLAB[END_REF], [START_REF] Sturm | Using seDuMi 1.02, a Matlab toolbox for optimization over symmetric cones[END_REF]), one obtains the optimal attenuation level γ opt = 0.89. Then the applied controller is a convex combination of the controllers synthesized at the vertices ρ, ρ (for more details about the polytopic approach, refer to [START_REF] Apkarian | A convex characterization of gain scheduled H ∞ controllers[END_REF]).

Figure 4 shows the steering and braking controller output according to ρ. As the braking weight has been described as parameter dependent, it is shown that when ρ = 10 braking signal is attenuated, and conversely, when ρ = 0.1 braking gain is larger. As a consequence, when ρ is low (resp. high), braking is activated (resp. disabled). Intermediate values will give intermediate behaviors. Remember that, for any ρ ∈ ρ ≤ ρ ≤ ρ , the closed loop stability is guaranted thanks to the LP V design and the polytopic approach.

According to the sensitivity functions Bode diagrams illustrated in figure 5, it is to make the following deductions:

• e ψ ψd
emphasizes that the yaw rate tracking performance satisfies the required specifications (see figure 5 The desired yaw moment M * z determined by the upper controller can be generated by applying brake torque to appropriate wheels. In this study, to avoid overlapping with front steering actuators, only rear wheels are concerned in the control law:

T * brj = 2RM * z t r (10)
In the following, we develop the control law of the lower controller (assuming counterclockwise positive).

There are two situations that accompany yaw instability:

• for the understeer condition ( ψ ≺ ψd ), the control moment, M * z , is generated by applying braking torque on the inner rear wheel.

• for the oversteer condition ( ψ ≻ ψd ), the control yaw moment, M * z , is generated by applying braking torque on the outer rear wheel.

Monitor: coordination strategy of steering and braking

As the brake-based DYC technique is not desirable in normal driving situations because of its direct influence on the longitudinal dynamics, the aim of the monitor is to minimize the usage of the braking. Consequently, the braking actuators must only be used when the vehicle goes toward instability. Since vehicle stability is directly related to the sideslip motion of the vehicle, judging the vehicle stability region is derived from the phase-plane (β -β). A stability bound defined in [START_REF] He | Coordination of active steering, driveline, and braking for integrated vehicle dynamics control[END_REF] is used here, that is formulated as:

χ < 1, (11)
where χ = 2.49 β + 9.55β is the "Stability Index". Moreover, the control task is also supposed to provide a seamless introduction of the direct yaw moment control, when it is required. Hence, the scheduling parameter ρ(χ) can be defined as:

/W 1 , 1/W 2 , 1/W 3 . ρ :=        ρ if χ ≤ χ (steering control) χ -χ χ -χ ρ + χ -χ χ -χ ρ if χ < χ < χ (steering+braking) ρ if χ ≥ χ (steering+full braking) (12)
where χ = 0.8 (χ is user defined) and χ = 1. To calculate the actual stability index χ, the block called "Sideslip dynamics observer" (see figure 3) evaluates β and β in real-time:

• β can be reconstructed using available sensors, according to the following relationship:

β = a y v x -ψ, (13)
where a y is the lateral acceleration, and v x is approximated by the mean of the rear wheel velocities. • β is not available using standard sensors, and thus, it must be estimated. "β-estimation" is widely discussed in the literature, and many papers deal with that topic (see [START_REF] You | New adaptive approaches to real-time estimation of vehicle sideslip angle[END_REF]). This estimation problem was also studied by the authors in a previous publication [START_REF] Doumiati | Observers for vehicle tyre/road forces estimations: experimental validation[END_REF]).

Actuator models

The corrective steer angle and rear braking torques control signals can be generated via actuation systems. In this particular research, let us consider the following actuators:

• Steer-by-wire Active Steering (AS) system providing an additional steering angle. This actuator is modeled as: δ+ = 2πκ(δ *δ + ) (14) where, κ = 10Hz is the actuator cut-off frequency, δ * and δ + are the steering controller and actuator outputs respectively. This actuator is bounded between

[-5 • , +5 • ].
• Brake-by-wire Electro Mechanical Braking (EMB) actuators providing a continuously variable braking torque. The EMB model is given by:

Ṫ + brj = 2πf 2 (T * brj -T + brj ) (15)
where, f 2 = 10 Hz is the actuator cut-off frequency, T * brj and T + bij are the local braking controller and actuator outputs respectively. This actuator control is limited between [0, 1200] Nm.

Full vehicle model

The block labelled "Vehicle simulation model" represents a full nonlinear vehicle model validated on a real "Renault Mégane Coupé" car. This model will be used in simulation for validation purpose only. Note that the main interest in using the full vehicle model is that it allows to take into account nonlinear load transfer, suspension dynamics, slipping and sideslip angles that are essential factors which play a major role on the global chassis dynamics, especially in critical driving situations. For more information about the vehicle model, see Poussot-Vassal et al. [2010].

SIMULATION RESULTS AND ANALYSIS

To clarify the effects of the proposed control scheme, both the vehicle dynamics with and without controllers are checked and compared. Different simulations were carried out on the nonlinear vehicle model platform briefly introduced in Subsection 2.5. Due to space limitations, only a double-lane-change test on a dry road (µ = 0.9) is reported in this paper. The vehicle is driven at 100 km/h. Figure 6 shows the yaw rate response versus the steering input. We can deduce that the uncontrolled vehicle becomes rather unstable as the amplitude of the steering input becomes larger. On the other hand, the controlled output of the yaw rate has nearly converged to the output of the desired linear model. Comparisons between the sideslip angles and the trajectories of the uncontrolled and controlled vehicles are illustrated in figure 7. Note that the sideslip angle is not specifically controlled in this study, it is simply the effect of the yaw rate control. According to these results, it is obvious that the handling performances are much improved by the proposed controller. Figure 8 illustrates how the stability index and the dependency parameter ρ evolve according to the driving situations. As stated before, when the stability index, χ, is below 0.8, only steering control is involved to enhance the handling performances. Therefore, ρ is equal to 10 and the corrective yaw moment is penalized. On the other hand, when χ exceeds 0.8, the braking system starts to work with the active steering to keep the vehicle stable. When χ becomes superior to 1, ρ takes the value 0.1, and braking is fully activated.

Figure 9 shows the generated corrective steering angle and the brake torques to enhance the lateral vehicle control. It is worth noting, that despite this test agressivity, actuators are far from saturation that may lead to instability.

CONCLUSION

The focus of this work is on presenting the vehicle yaw stabilizing problem in the framework of a control scheme.

The intent is to enhance vehicle steerability and stability during dynamic vehicle handling maneuvers. In this context, a new LP V /H ∞ controller that coordinates steering and braking actuators, is designed in this report. This LP V controller is designed in an original way and ensures that: 1)Steerability is enhanced in normal driving condition, and 2) Braking is used only when the vehicle tends to instability. Simulations of critical driving situations that compare the responses of a controlled vehicle with respect to a passive vehicle show the effectiveness of the proposed control design. Future work may consist to implement suspension actuators in the control scheme to enhance stability and prevent rollover situations.

  Fig. 1. Global control scheme.

Fig. 3 .

 3 Fig. 3. Generalized plant model for synthesis.

Fig. 4 .

 4 Fig. 4. Bode diagrams of the controller outputs δ * and M * z . approach, the closed loop stability is guaranteed for any ρ ∈ ρ, ρ . • The steering control is activated especially in the specified frequency range [1Hz, 10Hz] where the driver cannot act (see figure 5(c)). b) Lower controller: braking control scheme

  Closed loop transfer functions between M * and exogenous inputs.

  loop transfer functions between δ * and exogenous inputs.

Fig. 5 .

 5 Fig. 5. Closed loop transfer functions: LP V (red dashed (ρ = 0.1), or green dashed (ρ = 10)) synthesis results; Inverse of weighting functions (black dashed) of 1/W 1 , 1/W 2 , 1/W 3 .

Fig. 6 .

 6 Fig. 6. Response of the yaw rates versus steering wheel angle: Controlled (red dash-dot), Uncontrolled (blue dot), Reference model (black solid)

Fig. 7 .

 7 Fig. 7. Sideslip angles and trajectories of the controlled and uncontrolled vehicles: Controlled (red dash-dot), Uncontrolled (blue dash)

Fig. 8 .

 8 Fig. 8. M * z and ρ variations according to χ for the double lane-change maneuver.

Fig. 9 .

 9 Fig. 9. Control signals generated by the controller for the double lane-change maneuver.

Table 1 .

 1 Notations and vehicle parameters. adopted: index i = {f, r} and j = {l, r} are used to identify vehicle front, rear and left, right positions respectively. Table1summarizes the vehicle parameters, notations and values used in the paper.

	Symbol Value	Unit	Signification
	m	1535	kg	vehicle mass
	mr	648	kg	vehicle rear mass
	Iz	2149	kg.m 2	vehicle yaw inertia
	C f	40000	N/rd	cornering stiffness of front tires
	Cr	40000	N/rd	cornering stiffness of rear tires
	l f	1.4	m	distance COG -front axle
	lr	1	m	distance COG -rear axle
	tr	1.4	m	rear axle length
	h	0.5	m	height of the center of gravity
	µ	[2/5; 1]	-	tire/road contact friction interval
	v	[50; 130]	km/h	vehicle velocity interval
	g	9.81	m/s 2	gravitational acceleration

ACKNOWLEDGMENTS

This work is supported by the national french project INOVE/ ANR 2010 BLAN 0308.