J Murillo
email: murillo@cifasis-conicet.gov.ar

S Guillaume

F Spetale

E Tapia

P Bulacio

Set characterization-selection towards classification based on interaction index

Keywords: Choquet, Subset characterization, HLMS, Generalized Shapley index

In many real world datasets both the individual and coordinated action of features may be relevant for class identification. In this paper, a computational strategy for relevant feature selection based on the characterization of redundant or complementary features is proposed. The characterization is achieved using fuzzy measures and an interaction index computed from fuzzy measure coefficients. Fuzzy measure identification requires raw data to be turned into confidence degrees. This key step is carried out considering the distributions of feature values across all the classes. Fuzzy measure coefficients are then estimated with an improved version of the Heuristic Least Mean Squares algorithm that includes an efficient management of untouched coefficients. Then, a generalization of the Shapley index for an arbitrary number of features is used. Simulations experiments on synthetic datasets are performed to study the behavior of this generalized interaction index. For extreme datasets, containing either redundant or complementary features as well as noise, the index value is defined by mathematical formula. This result is used to motivate feature selection guidelines that take into account feature interactions. Experimental results on benchmark datasets show that the proposal allows for the design of compact, interpretable and competitive classification models.

Introduction

Feature selection is a critical task in classification problems, especially when dealing with noisy data of complex structure. Common feature selection techniques rely on individual feature evaluations [START_REF] Guyon | An introduction to variable and feature selection[END_REF] assuming independence between features. However, this condition is rarely met in practice [START_REF] Grabisch | Aggregation functions: Means[END_REF][START_REF] Kuncheva | Nearest neighbor classifier: Simultaneous editing and feature selection[END_REF]. For example, in micro-array datasets, the expression of certain genes, known to be biologically important for class differentiation, may appear relevant only when evaluated within groups [START_REF] Subramanian | Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles[END_REF]. This suggests that in complex data domains, even if the support of an individual feature may be small, its contribution to the support of feature subsets may be significant.

Diverse techniques have been proposed for evaluating the support of feature subsets in classification problems, including genetic algorithms [START_REF] Da Silva | Improving the ranking quality of medical image retrieval using a genetic feature selection method[END_REF], statistical methods based on mutual information [START_REF] Peng | Feature selection based on mutual information: Criteria of max-dependency, max-relevance, and min-redundancy[END_REF] and fuzzy measures [START_REF] Denneberg | Non-Additive Measure and Integral[END_REF][START_REF] Grabisch | Fundamentals of Uncertainty Calculi, with Applications to Fuzzy Inference[END_REF][START_REF] Wang | Fuzzy measure theory[END_REF]. A main drawback of genetic approaches is that they cannot explain the selection criteria; this drawback is partially solved by mutual information approaches able to explain the selection of feature pairs. On the other hand, fuzzy measures naturally characterize the relevance of feature subsets of arbitrary cardinality for general decision making problems. Hence, they are good candidates for tackling the feature selection and classification problem.

In practice, the estimation of fuzzy measures for n features entails the identification of 2 n coefficients. For small values of n, such identification problem may be solved by a domain expert. However, when n gets moderate or large, machine learning procedures are required. Both genetic algorithms [START_REF] Wang | Genetic algorithms for determining fuzzy measures from data[END_REF][START_REF] Wang | A genetic algorithm for determining nonadditive set functions in information fusion[END_REF][START_REF] Wang | A new type of nonlinear and the computational algorithm[END_REF][START_REF] Wang | Monotone set functions defined by choquet integral[END_REF] and gradient descent approaches [START_REF] Mikenina | Improved feature selection and classification by the 2-additive fuzzy measure[END_REF][START_REF] Grabisch | A new algorithm for identifying fuzzy measures and its application to pattern recognition[END_REF][START_REF] Marichal | Determination of weights of interacting criteria from a reference set[END_REF][START_REF] Wang | Detecting constructions of nonlinear integral systems from input-output data: an application of neural networks[END_REF][START_REF] Jia | Using neural networks to determine sugeno measures by statistics[END_REF] have been considered in literature. The Heuristic Least Mean Squares (HLM S) algorithm [START_REF] Grabisch | A new algorithm for identifying fuzzy measures and its application to pattern recognition[END_REF][START_REF] Grabisch | An empirical study of statistical properties of the choquet and sugeno integrals[END_REF][START_REF] Murillo | Revised hlms: A useful algorithm for fuzzy measure identification[END_REF], based on the Choquet Integral, is a representative of the gradient methods. They exhibit two important properties, the traceability of fuzzy measure estimation due to the lack of random steps and the efficient use of training data [START_REF] Grabisch | The representation of importance and interaction of features by fuzzy measures[END_REF]. They also need much less space. While multiple instances of fuzzy measure coefficients are required by genetic algorithms, HLM S uses only one. Furthermore, the HLM S heuristic based on the error function minimization of individual samples, usually leads to less extreme solutions [START_REF] Grabisch | A review of methods for capacity identification in choquet integral based multi-attribute utility theory: Applications of the kappalab r package[END_REF]. We note, however, that the original HLM S formulation has some convergence problems and thus, its revised HLM Sr version [START_REF] Murillo | Revised hlms: A useful algorithm for fuzzy measure identification[END_REF] is recommended.

Fuzzy measures weight all the possible combinations subsets. Hence, to absolutely characterize a feature subset we need to consider its support 2 across all subsets. This can be accomplished by computing the generalized interaction index proposed in [START_REF] Grabisch | k-order additive discrete fuzzy measures and their representation[END_REF], hereafter called GI. Its computation requires reliable fuzzy measures and thus, special attention must be paid to raw data preprocessing steps and to the fuzzy measures identification process itself.

In this paper, a robust feature selection strategy, based on the GI characterization of feature subsets regarding their relevance to classification problems, is proposed. With this intention, we first analyze GI results when standard raw data conversion procedure is used for the HLM Sr estimation of fuzzy measures. As a result of this study, a raw data conversion procedure based on normalized confidence degrees and an improved version of the HLM Sr algorithm are proposed. Aiming to understand some unexpected GI results previously reported in literature [START_REF] Mikenina | Improved feature selection and classification by the 2-additive fuzzy measure[END_REF][START_REF] Grabisch | k-order additive discrete fuzzy measures and their representation[END_REF], we study the GI behavior on homogeneous datasets, i.e., datasets where, in each class, all features are either redundant or complementary possibly including noisy features. As a result of this study, some guidelines for the efficient GI characterization of feature subsets in homogeneous datasets are introduced.

The outline of the paper is as follows. In Section 2, fuzzy measures, the discrete Choquet integral and GI are briefly reviewed. Then, fuzzy measure identification process with HLM Sr is analyzed, and some modifications to the algorithm are introduced in Section 3. Section 4 deals with the conversion from raw data into confidence degrees. In Section 5, the GI behavior is studied yielding two theorems (proofs are in the appendix). The usefulness of feature subset selection guided by GI characterizations is evaluated in Section 6 with UCI1 benchmark classification problems. Section 7 summarizes the main conclusions and perspectives.

Preliminaries

This section introduces basic concepts related to fuzzy measures, to the discrete Choquet integral and to the GI index. The reader may refer to [START_REF] Grabisch | Fundamentals of Uncertainty Calculi, with Applications to Fuzzy Inference[END_REF][START_REF] Wang | Fuzzy measure theory[END_REF][START_REF] Torra | Modeling Decisions[END_REF] for further details. Let us consider a finite set X = {x 1 , . . . , x i , . . . , x n } and let P(X) denote its power set.

Fuzzy measures and the discrete Choquet integral

A fuzzy measure (FM) is a set function µ : P(X) → [0, 1] fulfilling the following two axioms:

1. Normalization: µ(∅) = 0, µ(X) = 1 2. Monotonicity: A ⊆ B ⊆ X ⇒ µ(A) ≤ µ(B)
While the former allows for fuzzy measure comparisons, the latter ensures that adding any element to a given subset does not make it less informative.

Fuzzy measures µ : P(X) → [0, 1] are used in the definition of discrete Choquet integrals. For a given f : X → + , its discrete Choquet integral C with respect to a fuzzy measure µ : P(X) → [0, 1] is defined as follows:

C µ f (x 1 ), .., f (x n ) n i=1 f (x (i) ) -f (x (i-1) ) µ(A (i) ) (1)
where

x (•) is the rearrangement induced by f (x i ), i = 1, . . . , n, sorted in as- cending order, i.e., f (x (1) ) < • • • < f (x (n) ) and A (i) = {x (i) , x (i+1) , . . . , x (n) }; by convention f (x (0) ) = 0.

The GI characterization of features subsets

In the field of cooperative game theory, the Shapley index can be used to characterize the importance of individual features [25]:

φ i = K⊆X\i (n -|K| -1)!|K|! n! µ(K ∪ {i}) -µ(K) (2)
where | • | indicates the cardinality, and 0!=1 as usual. The Shapley value of µ is the vector φ = [φ 1 • • • φ n ] which has the property to be linear with respect to µ, and to satisfy:

n i=1 φ i = µ(X) = 1 (3)
This index has been generalized, first to characterize the importance of feature pairs [START_REF] Murofushi | Techniques for reading fuzzy measures (iii): interaction index[END_REF] and finally subsets of arbitrary cardinality [START_REF] Grabisch | k-order additive discrete fuzzy measures and their representation[END_REF]:

GI(A) = K⊆Y \A ξ(K) B⊆A (-1) |A|-|B| µ(K ∪ B) (4)
where

ξ(K) = (n -|K| -|A|)!|K|! (n -|A| + 1)!
Note that GI(A) reduces to the Interaction index when only two elements belong to A and it further reduces to Shapley index when A is a singleton. While the Shapley index is known to range in [0, 1] and the Interaction index in [-1, 1], the problem of characterizing the GI range for feature subsets of arbitrary cardinality remains open.

Improving HLM Sr for reliable fuzzy measure identification

To accomplish feature characterization and selection based on GI, a reliable fuzzy measure estimation is needed. A promising alternative is the improved version of HLM S [START_REF] Grabisch | A new algorithm for identifying fuzzy measures and its application to pattern recognition[END_REF] proposed by the authors in [START_REF] Murillo | Revised hlms: A useful algorithm for fuzzy measure identification[END_REF] called HLM Sr. Briefly, HLM Sr is an iterative supervised gradient-based algorithm for the identification of fuzzy measures. The algorithm starts from a set of m training samples involving n features and a target function T . At the beginning, fuzzy measures coefficients u are initialized to the so-called equilibrium state [START_REF] Grabisch | A new algorithm for identifying fuzzy measures and its application to pattern recognition[END_REF]. This initialization strategy reduces the Choquet integral to a simple arithmetic mean. At each iteration step HLM Sr updates the FM coefficient values, for each sample x j with j = 1, ..., m, according to the difference between the target T j and the current Choquet integral (Line 12). A given sample always uses the same coefficients to compute the integral, one for each subset size between 1 and n -1. The coefficient associated to a subset of size l is called u l .

The revised version used in this paper counts with some improvements with respect to the initial version [START_REF] Grabisch | A new algorithm for identifying fuzzy measures and its application to pattern recognition[END_REF]. Firstly, the update formula (Line 12) has been modified to get a true gradient and allows for coefficients to converge to the expected value with synthetic data sets. This modification is discussed in detail in [START_REF] Murillo | Revised hlms: A useful algorithm for fuzzy measure identification[END_REF]. Secondly, the management of untouched coefficients is changed. There might be coefficients which do not participate in any integral computation. As their initial value may limit the evolution of neighboring coefficients, they are identified (Line 6), and not taken into consideration for monotonicity checking (Line 13). The value of an untouched coefficient u l is set at the end of the algorithm (Line 18). The values fulfilling the monotonicity condition verify Eq.( 5):

max{u (l-1) } ≤ u l ≤ min{u (l+1) } (5)
These coefficients are set to the minimal authorized value in order not to influence GI values. Finally, HLM Sr proved sensitive to the default lexicographic order of equal valued features used in the computation of the core Choquet integral. In this case, only the coefficient of the first feature is updated. Since there is no reason to use a lexicographic order of features values, equal feature values are randomly ordered. As a result, features providing the same information are now considered equally relevant. for l ∈ (1 : n -1) do 12:

u l = u l -α × e j emax × (x j (n-l+1) -x j (n-l) ) {Coefficient update} 13:
Neighbors Monotonicity (u l ) {Monotonicity check } 14:

end for 15:

end for

16: E ← 1 m m j=1 (C u (x j ) -T j ) 2 {Error calculation} 17: until Stop Criterion(E) 18: Untouched Monotonicity(D, u) {Final monotonicity correction}
The inputs of HLM Sr are commensurable confidence degrees. The first step of the feature characterization process has to transform raw feature values into confidence degrees.

Raw data to confidence degree conversion

Raw data cannot be used for fuzzy measure identification using HLM Sr. Indeed, what is needed is commensurable information, i.e., the values must be in the same scale and have a common meaning, for instance the higher the more likely to belong to the class (a detailed discussion can be found in [START_REF] Grabisch | On the extension of pseudo-boolean functions for the aggregation of interacting criteria[END_REF]). Many methods for raw data conversion have been proposed in literature, including Parzen Windows, possibilistic histograms and Gaussian densities [START_REF] Grabisch | The representation of importance and interaction of features by fuzzy measures[END_REF]. Among these proposals, Gaussian densities appear as a good option due to their conceptual simplicity and low computational complexity. They have been used, e.g. for image features aggregation [START_REF] Jullien | An attribute fusion system based on the choquet integral to evaluate the quality of composite parts, Instrumentation and Measurement[END_REF].

Using Gaussian densities, the relevance of an individual feature value for the identification of a given class is taken from the probability density distribution of feature values within the class. The transformed feature value is then interpreted as a partial evidence that the corresponding sample belongs to the given class: the higher the transformed feature value, the higher the evidence provided by the feature that the sample belongs to the given class. Overall sample support to a given class can be then obtained by the aggregation of partial evidences, e.g., by means of the discrete Choquet integral. Though useful and intuitive, the Gaussian transformation exhibits some problems when dealing with feature values falling in distribution tails. For example, let us consider the synthetic classification problem shown in Figure 1. This problem involves two relevant features (V1 and V2) and one noisy feature (V3). For class 1, features V1 and V2 are complementary, since both are needed to classify the samples. On the other hand, they are redundant for class 2. Although these observations are acceptable reflected by GI results shown in Table 4, HLM Sr classification results are poor: a detailed inspection of misclassified samples shows that the problem is at the small relevance assigned to feature values falling in the Gaussian tails. Actually, what is happening is that some transformed feature values are being underestimated, i.e., although the absolute relevance of a feature value for a given class may be low, it should be consider relevant if being absolutely irrelevant for remaining classes. This observation suggests that raw data conversion processes must also take into account the distribution of features across all classes. Hence, the following process to convert raw data into confidence degrees is proposed:

xj ik = xj ik { p s=1 xj is } (6)
where xj is is the frequency of the x j raw value in the Gaussian distribution of feature i for class s. xj ik represents the normalized confidence degree over all the classes. This normalization remarkably improved classification results (100% of accuracy for the Atom dataset). However, anomalous GI results were still observed for completely redundant features, i.e., with features having exactly the same values. Although, these type of features rarely occur in practice, they are frequent after raw data conversion into confidence degrees.

Let us consider the data shown in Figure 2. The problem involves two completely relevant redundant features (V 1 and V 2) as well as a noisy one, V 3, not shown). For both classes, any of the V 1 or V 2 features is able to classify the samples. This leads to binary confidence degrees, either {1, 1} or {0, 0} for the pair V 1-V 2. In this case only the pair, and none of the singletons, would be updated, yielding unexpected GI values.

This problem is solved by introducing a small amount of random noise, ε = N (0, 0.01):

xj ik = xj ik { p s=1 xj is } + ε (7)
Thanks to the noise introduction, the GI results, shown in Table 4, are the expected ones and the classification is correctly achieved. This study highlights the importance of this preprocessing step. Now that fuzzy measures seem to be correctly identified, the next step is to analyze the GI behavior. 9

Author-produced version of the article published in Fuzzy Sets and Systems, 2015, N°270, p. 74-89 The original publication is available at http://www.sciencedirect.com http://dx.doi.org/10.1016/j.fss.2014.09.015

Characterization of feature subsets by fuzzy measures and the GI index

To gain insight into the GI behavior, simulations were performed using a synthetic Master dataset with seven normalized and commensurable features, i.e., normalized confidence degrees, A1 to A7, and one binary target T (see Table 3). It can be observed that features A1 and A2 are redundant for class 1 (T=1) since either A1 or A2 can be used for class identification. Features A3 and A4 are complementary for class 1. Finally, features A5 to A7 can be considered noise for both classes. From the Master data, three datasets, Inf , Red and Comp. They correspond to typical situations like noise, redundancy or complementariness. The Inf dataset, A1, A5 to A7 and T , contains only one informative feature, A1. The Red one, A1, A2, A5, A6 and T , includes two redundant features, A1 and A2. In the Comp dataset, A3, A4, A5, A6 and T , A3 and A4 are complementary for class 1. These three datasets are used to analyze the effectiveness of HLM Sr for fuzzy measure estimation and the posterior computation of GI values. In these experiments, HLM Sr was set to work with a maximum of 1000 iterations and a learning rate α = 0.01. Results on the Inf dataset are shown in Table 4. Fuzzy measure coefficients turn to be binary valued, being one for the informative feature A1, and for all subsets containing it, in agreement with the monotonicity constraint. As expected, they are zero for all subsets including only noisy features. Although being informative, fuzzy measure coefficient may be difficult to analyze due to the monotonicity constraint. The GI values make the analysis easy: in this case, the only non null value is the one of the informative feature. 5. As expected, Fuzzy measure coefficients are zero for all combinations of noisy features. For subsets containing features A1 or A2 they seem to indicate the relevance of these features but do not directly designate them. Furthermore, no insight about their interaction can be inferred. GI values help. For singletons, the significant GI values, known to range in [0, 1], are those of features A1 and A2. GI(A2) > GI(A1) means that the evidence brought by A2 is greater than the one by A1, this is in agreement with the data. Only one pair is given a significant GI value, known to range in [-1, 1], A1-A2. The negative sign characterizes this subset as redundant. For subsets of cardinality three, the GI shows some type of interaction in A1-A2-A5. This is in agreement with the residual correlation between A5-A1 and A5-A2 regarding class 1. 6. Fuzzy measure coefficients are one only for all subsets containing the complementary pair A3-A4 and zero for the remaining subsets including singletons A3 and A4 since neither A3 nor A4 can be individually used to identify class 1. GI values are the same for both singletons, indicating their equal relevance, and it is also one for the pair. This positive interaction characterizes the complementariness of the two features. Above mentioned results confirm that the GI is positive valued for complementary feature pairs and negatively valued for complementary ones. In addition, they suggest that GI characterizes as relevant all subsets in it power set. These observations hold for higher cardinality. For instance, if subset A1-A2-A3 is complementary, GI also identify subsets A1-A2, A1-A3, A2-A3 as complementary, and the three singletons as relevant. Moreover, the sign of GI for both complementary and redundant sets of features of cardinality three was positive, suggesting an alternating behavior for redundant feature subsets.

FM set A1 A5 A6 A7 A1-A5 A1-A6 A1-A7 A5-A6 A5-A7 A6-A7 A1-A5-A6 A1-A5-A7 A1-A6-A7 A5-A6-A7 FM value 1 0 0 0 1 1 1 0 0 0 1 1 1 0 GI 1 0 0 0 0 0 0 0 0 0 0 0 0 0
FM set A1 A2 A5 A6 A1-A2 A1-A5 A1-A6 A2-A5 A2-A6 A5-A6 A1-A2-A5 A1-A2-A6 A1-A5-A6 A2-A5-
These trends can be formalized in extreme situations, where a subset of features are fully redundant or complementary while the remaining only bring noise. The two following theorems are proved (the proof is given in the Appendix ):

Theorem 1 -Complementary features

Let N be a dataset with n features among which, for a given class, c of them are fully complementary and n-c are noise. Let C be the subset defined by the c fully complementary features. The GI for power set members of C with cardinality a ∈ {1, . . . , c} is:

GI(a) = 1 c -a + 1
This yields 1/c for each of the c singletons and 1 for the whole feature subset C.

Theorem 2 -Redundant features

Let N be a dataset with n features among which, for a given class, r of them are fully redundant and n -r are noise. Let R be the subset defined by the r fully redundant features. The GI for power set members of R with cardinality a ∈ {1, . . . , r} is:

GI(a) = (-1) a+1 r -a + 1
This yields 1/r for each of the r singletons, 1 for the whole feature subset R when r is odd and -1 when r is even.

The above results are used to propose some feature selection guidelines. Pairs: In this case, the union of complementary feature pairs, i.e., those with a significant positive GI, must be enriched with features coming from redundant feature pairs, i.e., those with a significant negative GI. This can be accomplished by checking that at least one of the features in each redundant feature pair is in the final set of selected features. A reasonable choice for GI thresholds is to set them to the mid value of the GI of fully complementary or redundant feature pairs already known to be ±1, i.e., GI ≥ 0.5 for the complementary case and GI ≤ -0.5 for the redundant one. Alternatively, they can be set according to expert knowledge, as done in [START_REF] Mikenina | Improved feature selection and classification by the 2-additive fuzzy measure[END_REF].

Higher cardinality sets: The pair selection reasoning can be extended to sets of higher cardinality. As follows from Theorems 1 and 2, the GI of fully complementary or redundant feature set of arbitrary cardinality is also ±1 and thus, setting GI threshold to ±0.5 remains valid.

To complete the classification process, once the Choquet integrals are computed for all the classes, using the selected features, the sample is assigned the class label, k, for which the support of the FM is maximum, as shown in Eq.( 8):

k = arg max k∈1,...,p C µ k (8)

Experimental results

The fuzzy measure and interaction index approach for feature selection and classification was evaluated on three UCI benchmark datasets (see Datasets were first pre-processed with the raw data conversion method described in Section 4. Fuzzy measure coefficients required for the computation of the GI were estimated with the modified HLM Sr algorithm described in Section 3. For this purpose, the HLM Sr algorithm was set to work with a maximum of 2000 iterations and a learning rate α = 0.05.

Classification performance was evaluated by means of the classification accuracy in 5-Fold cross validation experiments taking care that the original proportion of samples per class was preserved. Finally, the robustness CI classification under GI feature selection was also evaluated. For this purpose, features were ranked using the sum of their positions in the five GI ranks. This global ranking was then used to measure CI classification accuracy using feature subsets of increasing cardinality. GI feature selection results were compared with those obtained with the SV M -RF E [START_REF] Guyon | Gene Selection for Cancer Classification using Support Vector Machines[END_REF] feature selection technique based on a core SVM classifier set to work with a radial kernel and default constant complexity C=1. Features selected by the GI and SV M -RF E techniques were used to evaluate the performance of SVM and CI classifiers. In these evaluations, the number of GI selected features was used as a cut-off for SV M -RF E.

Iris dataset

As shown in Table 8, fuzzy measure coefficients of V 1, V 2 and V 1-V 2 are null in all classes suggesting that they are irrelevant for the classification task. For class 1, the GI = -1 for coalition V 3-V 4 points out that these features are completely redundant, i.e., only one of them is needed. This is confirmed by the data projection in V 3-V 4, shown in Figure 3, and by the high relevance, GI = 0.5, of individual features V 3 and V 4.

Table 10: Breast Cancer -Features rankings induced by the SV M -RF E and the GI methods. Features relevant to at least one class in the GI method are marked with * . The number of misclassified samples by SV M and CI classifiers is shown at the bottom of each of partition, from P 1 to P 5 .

GI feature selection with a global rank

V 6, V 1, V 3, V 7, V 9, V 2, V 8, V 4, V 5

Wine dataset

As reported in [START_REF] Aeberhard | Comparison of classifiers in high dimensional settings[END_REF], classes in the wine dataset are separable. However, only the regularized discriminant analysis method has been able to accom-Table 11: Wine -Features rankings induced by the SV M -RF E and the GI methods. Features relevant to at least one class in the GI method are marked with * . The number of misclassified samples by SV M and CI classifiers is shown at the bottom of each of partition, from P 1 to P 5 .

The four most significant features according SV M -RF E, (V 7, V 1, V 12 and V 13) belong to set of relevant features selected by GI. As with Breast cancer, GI selected features are also good for SVM classification but the converse does not hold. Stability of CI classification under GI feature selection with a global rank V 13, V 7, V 10, V 11, V 12, V 1, V 2, V 4, V 3, V 5, V 6, V 9, V 8 is shown Fig. 4. The classification error stabilizes after considering 6 out of 13 features. This suggests that the CI classifier is robust to noise, since adding irrelevant features does not affect the classification error.

Conclusions

A computational strategy for the selection of relevant features based on feature interactions has been presented. The strategy is based on the careful characterization of feature interactions by means of fuzzy measures and the posterior computation of the generalized Shapley index, called GI. A detailed analysis of the raw data conversion process and the HLM Sr algorithm required for the identification of fuzzy measures revealed that both tasks needed to be improved. Regarding the raw data conversion process, a normalization step that takes into account the distribution of features across all classes was introduced. The computational complexity of the HLM Sr algorithm was remarkably improved by disregarding monotonicity checking of untouched coefficients during algorithm's evolution. A more reliable estimation of fuzzy measure coefficients has been gained by identifying and solving the HLM Sr sensitivity to equal valued features.

A detailed analysis of the GI index on synthetic datasets containing noisy together with either redundant or complementary features uncovered a surprising GI behavior. The GI turns to be always positive in datasets with complementary features and of alternating sign in those with redundant ones. This result was used to propose some feature selection guidelines. Experimental results on benchmark datasets involving a dozen of features showed that competitive and interpretable fuzzy integral classifiers can be designed in this way. Furthermore, selected features seem to be also useful for the design of accurate SVM classifiers.

Further work is still needed to allow interpretable feature selection in datasets involving hundreds of features. In this regard, studies are first required to investigate the GI behavior involving both complementary and redundant features. Moreover, the reduction of both space and time required by the fuzzy measure identification to make the process tractable with huge datasets is still an open challenge.

Appendix

Mathematical formalization of features characterization for complementary and redundant features using GI

In extreme situations, when the only relevant features are either redundant or complementary, the GI expressions are easy to compute. Theorem 1 -Complementary features. Let N be a dataset with n features among which, for a given class, c of them are fully complementary and n-c are noise. Let C be the subset defined by the c fully complementary features. The GI for power set members of C with cardinality a ∈ {1, . . . , c} is:

GI(a) = 1 c -a + 1 19
This yields 1/c for each of the c singletons and 1 for the whole feature subset C.

Proof : The set N of n features includes a set C ⊂ N of c fully complementary features, while the remaining, N \ C, only bring noise. In this extreme situation, the only coefficients different from zero are the ones of C and the ones needed to keep monotonicity, i.e., the coefficients in C ∪ L with L ∈ P(N \ C). The value of all these coefficients is 1, as shown in Table 6 for c = 2. For given A ⊆ N , its GI is given by the Eq.( 4):

I(A) = K⊆N \A (n -k -a)!k! (n -a + 1)! B⊆A (-1) a-b µ(K ∪ B) (9)
Two cases have to be considered. The first case deals with subsets A entirely contained in C. The second case deals with subsets A having non-null fuzzy measure values only due to the monotonicity condition, i.e., A = C ∪ L with L ∈ P(N \ C).

• A ⊆ C, a ≤ c
In this case, B should be equal to A, otherwise, elements of C would be out of K ∪ B, and the resulting coefficients are 0. As a result, a -b = 0 and (-1) a-b is always positive. The non null coefficients are those of the combinations for which C ⊆ K ∪ A, meaning all the ones included in K = {C \ A} ∪ S where S ⊆ N \ C. As a result, k = c -a + s, µ(K ∪ B) = 1 and the Eq.( 9) becomes:

I(A) = S⊆N \C (n -c -s)!(c -a + s)! (n -a + 1)! (10)
This comes to:

I(A) = 1 (n -a + 1)! n-c s=0 n -c s (n -c -s)!(c -a + s)! (11)
As

n -c s = (n -c)! s!(n -c -s)! : I(A) = (n -c)! (n -a + 1)! n-c s=0 (c -a + s)! s! (12)
Remarking that (c -a + s)! s! = c -a + s s (c -a)!, it becomes:

I(A) = (n -c)!(c -a)! (n -a + 1)! n-c s=0 c -a + s s (13)
It can be proven, by recurrence, that 2 :

n s=0 m + s s = m + n + 1 n , so replacing n-c s=0 c -a + s s by n -a + 1 n -c
finally yields:

I(A) = 1 c -a + 1 (14) • C ⊂ A, a > c
All the sets B to consider must include C, they are of size b = c + p. There are a-c p of size c + p in A which include C. The GI formula can be written as:

I(A) = 1 (n -a + 1)! n-a s=0 (n -a -s)!s! a-c p=0 a -c p (-1) a-c-p (15)
The second sum is the binomial formula:

n p=0 n p (-1) n-p = 1 + (-1) n = 0.
2 Recurrence relation: Hence, the only non-zero GI index values for a set N of n features where c of them are fully complementary and n -c are noise are those belonging to subsets A in P(C). Furthermore, GI values depend just on the cardinality a of such subsets: GI(A) = 1 c -a + 1 .

n+1 s=0 m + s s = m + n + 1 n + m + n + 1 n + 1 = m + n + 2 n + 1
Theorem 2 -Redundant features. Let N be a dataset with n features among which, for a given class, r of them are fully redundant and n -r are noise. Let R be the subset defined by the r fully redundant features. The GI for power set members of R with cardinality a ∈ {1, . . . , r} is:

GI(a) = (-1) a+1 r -a + 1
This yields 1/r for each of the r singletons, 1 for the whole feature subset R when r is odd and -1 when r is even.

Proof : The set N of n features includes a set R ⊂ N of r fully redundant features, while the remaining , N \ R, only bring noise. In this extreme situation, all the coefficients belonging to P(N \ R) are assigned zero while all the subsets including at least one element from R are assigned a value of 1.Table 5 illustrates this trend with 2 partially redundant features. For given A ⊆ N , its GI formula in Eq.( 4) can be written as:

I(A) = K⊆N \A (n -k -a)!k! (n -a + 1)! B⊆A (-1) a-b µ(K ∪ B) (16)
As previously, µ(K ∪ B) is either 0 or 1. The number of elements to sum can be computed as the total number of combinations according to the cardinalities minus the number of null coefficients. µ(K ∪ B) is null if and only if both B and K are included in P(N \ R).

The total number of combinations is given by:

T otC =
Let's define A as an union: A = {Q ⊂ R} ∪ {L ⊂ N \ R}, with a = q + l, and q ≤ r and l ≤ n -r. The sum of null coefficients over all K is computed from the all combinations of N \ R in N \ A. It is given by:

N ulC = n-r-l k=0 n -r -l k (n -k -a)!k! (n -a + 1)! ( 18
)
This sum is non null only for B ⊆ A ⊆ P(N \ R).

Two cases have to be considered:

• A ⊆ R, l = 0
Substituting l by zero and developing the binomial coefficient, Eq.( 18) becomes:

N ulC = (n -r)! (n -a + 1)! n-r k=0 (n -k -a)! (n -r -k)! (19)
This is exactly Eq.( 12) with k = n -r -s.

So, when l = 0, N ulC = 1 r -a + 1 .

The number of subsets of size p ≤ a in A is a p , but the only subset for which this difference is non null is B = ∅. All the others partial sums are the same, with alternate signs. The final GI absolute value is N ulC, and the sign depends on the cardinality, a:

I(A) = (-1) a+1 1 r -a + 1 (20)
• l = 0

Two cases have to be considered. The first, trivial, one is q = 0, A ⊆ N \ R. In this case, all the subsets B give the same sum with alternate signs and the index is zero.

The remaining case includes both kinds of features: redundant and noisy ones. The partial sums are the same for all B ⊆ Q, T otC, and it is also the same, T otC -N ulC, for all B ⊆ L. Thanks to the alternate signs and the sum of the binomial coefficients, the overall index is zero.

For a set of r fully redundant features, while the others only bring noise, the GI values are null for all the combinations except for the power set of the r features. In this case, its value depends on the cardinality of the subset, a and is: (-1) a+1 1 r -a + 1 .

Algorithm 1 for j ∈ examples do 10 :

 110 Revised HLMS 1: Input: Training dataset D with samples (x j , T j ), j = 1, ..., m 2: Output: Fuzzy measure coefficients u = [u p ] , p ∈ P(X) 3: for p ∈ P(X) do {Initialization} e j = C u (x j ) -T j {Individual error calculation} 11:

Figure 2 :

 2 Figure 2: Synthetic Balls dataset with the V1-V2 projection. Class 1 is plotted with triangles (top) and class 2 with circles (bottom).

  FM setA3 A4 A5 A6 A3-A4 A3-A5 A3-A6 A4-A5 A4-A6 A5-A6 A3-A4-A5 A3-A4-A6 A3-A5-A6 A4-

Feature

  selection guidelines using GI Singletons: Relevant individual features are those with GI > 1 n [21].

  is shown Fig.4. The classification error stabilizes after considering 7 out of 9 features.

Figure 4 :

 4 Figure 4: Classification error (%) of CI classifiers under GI feature selection with an increasing number of features in order of relevance for Breast Cancer dataset(left) and Wine dataset (right).

  k -a)!k! (n -a + 1)!

7

  Author-produced version of the article published in Fuzzy Sets and Systems, 2015, N°270, p. 74-89 The original publication is available at http://www.sciencedirect.com http://dx.doi.org/10.1016/j.fss.2014.09.015

		4				Coalition Class 1 Class 2
						V1	0.52	0.22
		2				V2	0.33	0.69
	V2	0			x x x	V3	0.14	0.07
			x	x	x	V1-V2	0.61	-0.41
		-2				V1-V3	0.15	0.06
		-4				V2-V3	-0.24	-0.11
		-4	-2	0	2	4
				V1
	Figure 1: The synthetic Atom
	dataset in the V1-V2 projection.
	Class 1 is plotted with triangles
	(nucleon) and class 2 with circles
	(ring). Misclassified samples are
	marked with crosses.

Table 1 :

 1 GI results for the Atom dataset. For class 1, the GI in-

	dicates that both V1 and V2 are
	relevant and quite complementary.
	For class 2, the GI indicates that
	feature V2 is more relevant than
	V1 and that they are quite redun-
	dant.

Table 2 :

 2 GI results for the Balls dataset. For class 1, the GI indicates that both features V1 and V2 relevant and that they are completely redundant.

Table 3 :

 3 The synthetic Master dataset after relative raw data conversion across 16 samples. Columns A1 to A7 correspond to features values and column T to the target value.

	Sample	A1	A2	A3	A4	A5	A6	A7	T
	1	0.60	0.80	0.80	0.90	0.55	0.33	0.53	1
	2	0.60	0.55	0.80	0.80	0.97	0.01	0.29	1
	3	0.70	0.90	0.75	0.80	0.72	0.27	0.17	1
	4	0.60	0.80	0.70	0.60	0.63	0.30	0.56	1
	5	0.75	0.90	0.80	0.80	0.04	0.68	0.05	1
	6	0.65	0.60	0.90	0.90	0.00	0.47	0.03	1
	7	0.80	0.75	0.70	0.80	0.96	0.65	0.07	1
	8	0.70	0.90	0.80	0.70	0.96	0.81	0.16	1
	9	0.00	0.10	0.10	0.50	0.41	0.17	0.28	0
	10	0.10	0.20	0.40	0.10	0.78	0.35	0.46	0
	11	0.20	0.30	0.50	0.15	0.17	0.64	0.57	0
	12	0.10	0.20	0.10	0.60	0.45	0.73	0.18	0
	13	0.20	0.10	0.50	0.20	0.48	0.68	0.47	0
	14	0.30	0.35	0.60	0.10	0.08	0.49	0.60	0
	15	0.10	0.20	0.70	0.20	0.65	0.92	0.04	0
	16	0.20	0.25	0.10	0.80	0.61	0.25	0.50	0

Table 4 :

 4 Inf

dataset: Fuzzy measure (FM) coefficients and GI values Results on the Red dataset are shown in Table

Table 5 :

 5 Red dataset: Fuzzy measure (FM) coefficients and GI values Results on the Comp dataset are shown in Table

	A6

Table 6 :

 6 Comp dataset: Fuzzy measure (FM) coefficients and GI values

Table 13

 13 Author-produced version of the article published in Fuzzy Sets and Systems, 2015, N°270, p. 74-89 The original publication is available at http://www.sciencedirect.com http://dx.doi.org/10.1016/j.fss.2014.09.015 7). The Iris dataset is a well known dataset, easy to analyze, and allows for a fair comparison with results reported in[START_REF] Grabisch | The representation of importance and interaction of features by fuzzy measures[END_REF].

		#Features	#Samples	#Classes
	Iris	4	150	3
	Breast Cancer	9	683	2
	Wine	13	178	3

Table 7 :

 7 Dataset characteristics 6.1. Experimental protocol

21

  Author-produced version of the article published in Fuzzy Sets and Systems, 2015, N°270, p. 74-89 The original publication is available at http://www.sciencedirect.com http://dx.doi.org/10.1016/j.fss.2014.09.015

22

  Author-produced version of the article published in Fuzzy Sets and Systems, 2015, N°270, p. 74-89 The original publication is available at http://www.sciencedirect.com http://dx.doi.org/10.1016/j.fss.2014.09.015

http://archive.ics.uci.edu/ml/

Author-produced version of the article published in Fuzzy Sets and Systems, 2015, N°270, p. 74-89 The original publication is available at http://www.sciencedirect.com http://dx.doi.org/10.1016/j.fss.2014.09.015

Acknowledgments

The authors are grateful to Tewfik Sari (Irstea, UMR ITAP) for the recurrence relation used in the theorems' proof.

Author-produced version of the article published in Fuzzy Sets and Systems, 2015, N°270, p. 74-89 The original publication is available at http://www.sciencedirect.com http://dx.doi.org/10.1016/j.fss.2014.09.015 For class 3, the GI = -0.56 for coalition V 3-V 4 points out that these features are redundant for the class. For classes 2 and 3, the GI of feature V 4 is higher than that of V 3, pointing out that V 4 is more relevant than V 3 for these classes.

Table 8: Iris -Fuzzy measure coefficients of feature subsets together with their GI values A comparison of these GI results with those reported in [START_REF] Grabisch | The representation of importance and interaction of features by fuzzy measures[END_REF] shows that feature singletons are assigned roughly similar relevance values, except for V 1

Author-produced version of the article published in Fuzzy Sets and Systems, 2015, N°270, p. 74-89 The original publication is available at http://www.sciencedirect.com http://dx.doi.org/10.1016/j.fss.2014.09.015 and V 2 which relevance values are close to zero in our case. When considering feature pairs, the highest absolute GI value in [START_REF] Grabisch | The representation of importance and interaction of features by fuzzy measures[END_REF] is met for coalition V 1-V 3 in class 2, a result quite difficult to explain. Similarly, coalition V 3 -V 4 in class 1 is assigned GI = -0.05 suggesting the lack of interaction between both features in this class, a result quite difficult to explain when observing the data.

Table 9: Iris -Features rankings induced by the SV M -RF E and the GI methods. Features relevant to at least one class in the GI method are marked with * . The number of misclassified samples by SV M and CI classifiers is shown at the bottom of each of partition, from P 1 to P 5 .

GI and SV M -RF E feature selection results are shown in Table 9. In both methods, features V 3 and V 4 are more relevant than V 1 and V 2 for the classification. Although classification errors and selected features are roughly similar for both approaches, in agreement with the 96.7 % reported in [START_REF] Grabisch | The representation of importance and interaction of features by fuzzy measures[END_REF] using all features, an interpretable feature characterization is additionally provided by the GI approach.

Breast Cancer dataset

This dataset has been widely used in scientific literature. To our knowledge, the best classification accuracy results have been obtained with IncNet (97.1%) [START_REF] Duch | A new methodology of extraction, optimization and application of crisp and fuzzy logical rules[END_REF]. With CI classifiers, the best reported classification accuracy is 91.5% [START_REF] Yue | Parameter estimation for choquet fuzzy integral based on takagi-sugeno fuzzy model[END_REF]. As shown in Table 10, GI feature selection and CI classification yield a classification accuracy of 97.5%. Thanks to the untouched coefficient management in HLM Sr around 425 coefficients are identified on average over the different classes (512 for HLM S).

The four most significant features according SV M -RF E (V 1, V 3, V 6 and V 7) belong to set of relevant features selected by GI (Table 10). It is worth noting that GI selected features are also relevant for SV M as indicated by the reduction of misclassified samples.


Xiangjun Qian

Florent Altché

Philipp Bender

Christoph Stiller

Arnaud De

La Fortelle

Optimal Trajectory Planning for Autonomous Driving Integrating Logical Constraints: An MIQP Perspective

This paper considers the problem of optimal trajectory generation for autonomous driving under both continuous and logical constraints. Classical approaches based on continuous optimization formulate the trajectory generation problem as a nonlinear program, in which vehicle dynamics and obstacle avoidance requirements are enforced as nonlinear equality and inequality constraints. In general, gradientbased optimization methods are then used to find the optimal trajectory. However, these methods are ill-suited for logical constraints such as those raised by traffic rules, presence of obstacles and, more generally, to the existence of multiple maneuver variants. We propose a new formulation of the trajectory planning problem as a Mixed-Integer Quadratic Program. This formulation can be solved efficiently using widely available solvers, and the resulting trajectory is guaranteed to be globally optimal. We apply our framework to several scenarios that are still widely considered as challenging for autonomous driving, such as obstacle avoidance with multiple maneuver choices, overtaking with oncoming traffic or optimal lane-change decision making. Simulation results demonstrate the effectiveness of our approach and its real-time applicability.

I. INTRODUCTION

Autonomous driving has been gaining impetus in the last few years, thanks to its foreseen potential of increasing traffic efficiency and reducing the number of road accidents. Recent studies suggest that up to 50% of vehicles may be automated by 2030 [START_REF] De La Fortelle | Network of Automated Vehicles: the AutoNet2030 vision[END_REF].

A challenging research task for autonomous driving is to generate optimal trajectories for certain criteria such as comfort or energy efficiency, while satisfying various constraints arising both from traffic rules and operational limits of the vehicles, some of them being highly specific to the context of driving. A large number of approaches for the motion planning or trajectory planning problem can be found in the literature [START_REF] Paden | A survey of motion planning and control[END_REF], among which one important category is the sampling-based methods. These methods [START_REF] Kuwata | Real-time motion planning with applications to autonomous urban driving[END_REF]- [START_REF] Gu | Focused trajectory planning for autonomous on-road driving[END_REF] typically generate a large set of candidate trajectories by deterministically or stochastically sampling a state space, before executing the best solution among those candidates. The optimality of the resulting trajectory relies on generating a high number of samples, which may not allow real-time applications.

{xiangjun.qian, florent.altche, arnaud.de la fortelle}@mines-paristech.fr {philipp.bender,stiller}@kit.edu 3 Florent Altché is also with Ecole des Ponts ParisTech.

Recently, new approaches based on model predictive control (MPC) [START_REF] Ziegler | Trajectory planning for bertha; a local, continuous method[END_REF]- [START_REF] Nilsson | Receding horizon maneuver generation for automated highway driving[END_REF] have attracted increased attention due to their ability to systematically handle system constraints and quickly find the optimal trajectory. MPC-based approaches rely on iteratively solving a constrained, finite horizon optimal control problem generally using nonlinear optimization techniques. These approaches are intrinsically more efficient than sampling-based ones in finding optimal trajectories, as they exploit useful information like gradient fields to orient the search of optimal trajectories. Real-time implementations of these methods on actual vehicles can be found, e.g., in [START_REF] Ziegler | Trajectory planning for bertha; a local, continuous method[END_REF], [START_REF] Gao | Model Predictive Control for Autonomous and Semiautonomous Vehicles[END_REF], [START_REF] Weiskircher | Predictive adas: A predictive trajectory guidance scheme for advanced driver assistance in public traffic[END_REF].

Even though MPC-based methods have yielded good results for problems with continuous constraints, we argue that there are constraints that are intrinsically logical, in that they are naturally formulated as logical propositions. First, onroad driving is by nature highly constrained by traffic rules and expected driving behaviors, which can be represented naturally as logical propositions, for instance:

(R.1) A vehicle must slow down when approaching a speed bump.

(R.2) A vehicle can only overtake from the left side. (R.3) A vehicle needs to be on the exit lane before exiting on the highway.

Moreover, driving generally involves discrete decisions among multiple maneuver variants, which makes trajectory planning combinatorial in nature [START_REF] Bender | The combinatorial aspect of motion planning: Maneuver variants in structured environments[END_REF], [START_REF] Tas | Integrating Combinatorial Reasoning and Continuous Methods for Optimal Motion Planning of Autonomous Vehicles[END_REF]. Indeed, it has been shown that each maneuver variant can be mapped to a unique homotopy class of trajectories [START_REF] Bender | The combinatorial aspect of motion planning: Maneuver variants in structured environments[END_REF], [START_REF] Tas | Integrating Combinatorial Reasoning and Continuous Methods for Optimal Motion Planning of Autonomous Vehicles[END_REF]. It is generally straightforward to describe these homotopy classes using logical propositions. For example, an obstacle in the middle of the road brings two distinct maneuver choices, as follows:

(R.4) A vehicle should avoid an obstacle, either in a clockwise or in a counterclockwise manner.

By nature, sampling-based approaches can accommodate logical constraints. However, as stated earlier, these methods can yield sub-optimal trajectories even with a large number of samples. Most current MPC-based approaches are illsuited to take these logical rules into account, as most of them require all constraints to be continuous and differentiable. Moreover, most continuous solvers use local, gradientbased optimization methods which can be trapped in a local optimum inside a certain homotopy class. Some heuristics have been proposed to approximate logical constraints by nonlinear functions [START_REF] Gao | Model Predictive Control for Autonomous and Semiautonomous Vehicles[END_REF], [START_REF] Weiskircher | Predictive adas: A predictive trajectory guidance scheme for advanced driver assistance in public traffic[END_REF], [START_REF] Bender | The combinatorial aspect of motion planning: Maneuver variants in structured environments[END_REF] and to initialize the algorithm in a homotopy class that is likely to be the best one [START_REF] Gao | Model Predictive Control for Autonomous and Semiautonomous Vehicles[END_REF], but they provide no guarantee regarding the global optimality of the generated trajectory.

A possible way of handling logical constraints is to formulate the problem as a Mixed-Integer Programming (MIP) problem. General MIP problems are hard to solve, but efficient algorithms exist for special instances of MIP, notably Mixed-Integer Linear Programming (MILP) and Mixed-Integer Quadratic Programming (MIQP). Both have successfully been used for Unmanned Aerial Vehicle (UAV) trajectory planning [START_REF] Schouwenaars | Safe receding horizon path planning for autonomous vehicles[END_REF], [START_REF] Bellingham | Receding horizon control of autonomous aerial vehicles[END_REF] and multi-vehicle collision avoidance problem [START_REF] Borrelli | Mpcbased approach to active steering for autonomous vehicle systems[END_REF]. MILP methods have also been used to design automatic lane change controllers that are proposed in [START_REF] Nilsson | Strategic decision making for automated driving on two-lane, one way roads using model predictive control[END_REF], [START_REF] Du | Autonomous lane-change controller via mixed logical dynamical[END_REF], but they require a simplified vehicle model. To the best of our knowledge, an efficient MIP formulation using realistic vehicle dynamics and capable of handling multiple on-road driving scenarios is not yet available.

The major contribution of this paper is to propose a novel approach for the trajectory generation problem using Mixed Integer Quadratic Programming. This approach allows to seamlessly treat continuous and logical constraints. Furthermore, we apply our framework to several situations which are still widely recognized as challenging for autonomous driving, for instance obstacle avoidance, overtaking in presence of oncoming traffic or optimal lane-change decision making. We demonstrate that our method can intuitively handle these problems and generate globally optimal trajectories.

II. OPTIMAL TRAJECTORY GENERATION AS MIQP

This section proposes an MIQP formulation for the problem of generating optimal trajectories for autonomous driving. We consider that the ego vehicle is driving on a road with lane markings. The ego vehicle uses a Model Predictive Control (MPC) scheme to perform trajectory planning, in which the planned trajectory is updated in successive stages. However, in this paper, we only discuss in detail a single MPC stage. We assume that the road curvature is sufficiently small to consider the road as straight within the horizon of one MPC stage, so that we can model the road in a Cartesian frame (x, y) with x the longitudinal direction of the road and y the lateral direction.

A. Vehicle dynamics and constraints

The ego vehicle is modeled as a third-order point-mass system, with state x and control vector u defined as

x(t) = [x(t), v x (t), a x (t), y(t), v y (t), a y (t)] T , u(t) = [j x (t), j y (t)] T
where x and y denote longitudinal and lateral position in the inertial frame, v the speed and a the acceleration, with subscript x or y respectively indicating their longitudinal and lateral components. The vehicle dynamics are ruled by the following differential equation:

ẋ(t) = A 0 0 A x(t) + B 0 0 B u(t), A =   0 1 0 0 0 1 0 0 0   , B =   0 0 1   , (1)
where 0 denotes zero matrices of proper dimensions.

For practical purposes, we assume that the control u is a piecewise constant function with a time step of τ . The discretized vehicle dynamics can then be approximated as

x(k + 1) = A d 0 0 A d x(k) + B d 0 0 B d u(k), A d =   1 τ 1 2 τ 2 0 1 τ 0 0 1   , B d =   1 6 τ 3 1 2 τ 2 τ   . (2)
where x(k) and u(k) are respectively the state and the control of the vehicle at the beginning of the time interval [kτ, (k + 1)τ ).

Remark that some previous work [START_REF] Nilsson | Receding horizon maneuver generation for automated highway driving[END_REF], [START_REF] Murgovski | Predictive cruise control with autonomous overtaking[END_REF] has used a first or second order point-mass description. Here, we use a third-order model to ensure that the second derivative of the generated trajectory is continuous even with a piecewise constant input. As a result, the yaw rate of the planned trajectory is continuous, and the vehicle can thus smoothly track the trajectory.

Due to the dynamic limitations of the vehicle, bound constraints are enforced on the state and control signals as

x ∈ [x, x], u ∈ [u, u], (3)
with the bounds defined as

x = [0, 0, a x , y, v y , a y ] T , x = [free, v x , a x , y, v y , a y ] T , (4a) u = [j x , j y ] T , u = [j x , j y ] T . (4b)
The above formulation does not consider the nonholonomic constraints of the vehicle, and the longitudinal and lateral dynamics are fully decoupled. The exact coupling between these dynamics involves nonlinear relations [START_REF] Gao | Model Predictive Control for Autonomous and Semiautonomous Vehicles[END_REF]; therefore, we approximate it by two additional constraints which ensure that the generated trajectory is dynamically feasible.

The vehicle heading θ can be reconstructed from the state x as θ = arctan(v y /v x ). We model the coupling of longitudinal and lateral dynamics by enforcing condition θ ∈ [θ, θ], with

v y ∈ [v x tan(θ), v x tan(θ)]. (5)

B. Logical constraints for driving

Winston et al [START_REF] Winston | Introduction to mathematical programming[END_REF] have shown that propositional logic can be further reformulated as a set of linear inequalities with integer variables. Here, we briefly present this approach in the context of autonomous driving.

We define a literal as an atomic statement corresponding to a linear mathematical condition on one of the state variables, for instance: the longitudinal position of the vehicle is larger than 30 m. Literals can be combined using connectors, namely ∧ (and), ∨ (or), ¬ (negation); implications (⇒) and equivalences (⇔) can be formed using the first three connectors.

To illustrate, we consider a speed bump covering the range of longitudinal positions x ∈ [30 m, 50 m], with a speed limit of 10 m/s in this range. Let us define three literals

P 1 = [x(k) ≥ 30], P 2 = [x(k) ≤ 50] and P 3 = [v x (k) ≤ 10
]; rule (R.1) can then be expressed in the form: ∀k ≥ 0, (P 1 ∧P 2 ) ⇒ P 3 , i.e. if 30 ≤ x(k) ≤ 50, then the vehicle's longitudinal speed should be lower than or equal to 10 m/s.

The so-called Big-M method (see [START_REF] Winston | Introduction to mathematical programming[END_REF]) allows to force a binary variable δ(k) to be equal to 1 when a given literal is true, and equal to 0 when the literal is false, using only linear inequalities by introducing a large constant M . In the previous example, we can let δ i (k) = 1 ⇔ P i for i ∈ {1, 2, 3} so that (R.1) can be expressed equivalently as, ∀k ≥ 0,

δ 1 (k) = 1 ⇔ x(k) ≥ 30, δ 2 (k) = 1 ⇔ x(k) ≤ 50, δ 3 (k) = 1 ⇔ v x (k) ≤ 10, -δ 1 (k) + δ 3 (k) ≤ 0, -δ 2 (k) + δ 3 (k) ≤ 0, δ 1 (k) + δ 2 (k) -δ 3 (k) ≤ 1. (6)
Note that these conditions could also be approximated using nonlinear constraints in a continuous solver. However, current state-of-the-art solvers for general nonlinear optimization problems can only find local optima, whereas our mixed-integer formulation ensures global optimality.

C. MIQP formulation

We introduce a new variable δ such that δ(k) = {0, 1} m is a collection of all binary variables resulting from the reformulation of relevant literals as mixed integer linear inequalities. Let x r be the reference trajectory for the vehicle, which can be time-dependent, state-dependent, or dependent on propositions. We also introduce δ r as the reference trajectory for the binary so that we can also express preferences on some binary states, for instance to specify a preferred lane in a multi-lane road. Let K = T /τ be the number of time steps in the prediction horizon. The optimization problem can now be formulated, in its generic form, as

min u,δ K k=0 ||x(k)-x r (k)|| 2 Q + ||δ(k) -δ r (k)|| 2 S + ||u(k)|| 2 R , (7)
subject to

known x(0), (8a)
vehicle dynamics (2), (8b) dynamic constraints (3) and ( 5), (8c)

C   x x r δ   ≤ D. ( 8d
)
where Q, S, R are positive weighting matrices of proper dimensions. Constraint (8d) is the set of all linear inequalities written in matrix form with two matrices C and D of proper dimensions, incorporating all mixed-integer constraints used to enforce driving rules.

x = [0, 0, -4 m/s 2 , 0, -2 m/s, -1 m/s 2 ] T , x = [free, 20 m/s, 3 m/s 2 , 5 m, 2 m/s, 1 m/s 2 ] T , u = [-3 m/s 3 , -2 m/s 3 ] T , u = [3 m/s 3 , 2 m/s 3 ] T , θ = -0.4 rad, θ = 0.4 rad, ω = -0.26 rad/s, ω = 0.26 rad/s, q 1 = 1, q 2 = 2, q 3 = 1, q 4 = 2, q 5 = 4, r 1 = 4, r 2 = 4.

TABLE I: Parameters used for case study

The cost function ( 7) is quadratic if x r and δ r are independent from other variables or are linearly dependent on other variables. The constraints (8) are linear. Therefore, the above optimization problem is an instance of a mixed-integer quadratic program (MIQP). Exact resolution algorithms are known to solve such problems, leading to a globally optimal trajectory which can then be fed to low-level controller charged of tracking this trajectory.

Remark 1: The proposed MIQP formulation is able to cover a vast majority of on-road driving scenarios. For example, Rules (R.1)-(R.4) can be naturally incorporated in our framework. If we must consider non-quadratic cost function, nonlinear vehicle dynamics or constraints, we can formulate a Mixed Integer Nonlinear Program, which, however, is much harder to solve.

III. CASE STUDY

The previous section proposes an MIQP framework that integrates the continuous dynamic of the vehicle and propositional logic. In this section, we explain how this framework can be used to effectively handle various on-road driving scenarios. Note that, although the major purpose of this section is to demonstrate the universality of the framework, the detailed strategies for different driving scenarios are by themselves contributions to some challenging problems for on-road autonomous driving.

Throughout this section, we use the following cost function

J = K k=0 q 1 (v x (k) -v r ) 2 + q 2 a 2 x (k) + q 3 (y(k) -y r (k)) 2 + q 4 v 2 y (k) + q 5 a 2 y (k) + r 1 j 2 x (k) + r 2 j 2 y (k), (9)
such that the vehicle tracks a constant desired speed profile and a potentially time-varying desired lateral deviation, while trying to minimize the control effort. We do not assume any desired binary state, and therefore the term

||δ(k) -δ r (k)|| 2 S
in the generic formulation is ignored. In all cases, the ego vehicle is assumed to start at position (0, 2.5) in the rightmost lane. The initial and desired speeds are equal, with v x (0) = v r = 15 m/s. The width of each lane is 5 m. Other parameters that are common in all cases are recapitulated in Table I. Scenario-specific parameters will be presented respectively in each case study.

We use the commercial solver Gurobi [START_REF]Gurobi optimizer reference manual[END_REF] to compute solutions to our MIQP formulation. Simulation codes are written in Python, and experiments are performed on a laptop with Intel Core i5-5300U CPU clocked at 2.30GHz with 8GB RAM.

A. Speed bump

The first case study considers the speed bump scenario (Fig. 1a) that is used as an example in II-B.

The speed bump conditions are given in [START_REF] Ziegler | Trajectory planning for bertha; a local, continuous method[END_REF]. In the simulation, we use a prediction horizon T = 5 s and τ = 0.25 s. Fig. 1b illustrates the longitudinal speed profile of the planned trajectory with respect to the traveled distance. We observe that vehicles effectively reduces its speed to less than 10 m/s within the interval of [30, 50].

B. Obstacle avoidance

We now consider an obstacle avoidance scenario during on-road driving. The irregular shapes of obstacles are approximated using minimal bounding rectangles. A more complex polygonal modeling is also possible, at the cost of increased computational complexity. For an obstacle ν with bounding rectangle

[x ν (k) -L ν , x ν (k) + L ν ] × [y ν (k) -W ν , y ν (k) + W ν ]
, the set of constraints for collision avoidance is then given as ∀k ≥ 0,

δ ν 1 (k) = 1 ⇔ x(k) ≤ x ν (k) -L ν , (10a) δ ν 2 (k) = 1 ⇔ x(k) ≥ x ν (k) + L ν , (10b) δ ν 3 (k) = 1 ⇔ x(k) ≤ y ν (k) -W ν , (10c) δ ν 4 (k) = 1 ⇔ x(k) ≥ y ν (k) + W ν , (10d) δ ν 1 (k) + δ ν 2 (k) + δ ν 3 (k) + δ ν 4 (k) = 1. ( 10e
)
Note that the formulation allows both moving and still obstacles. The conditions [START_REF] Bender | The combinatorial aspect of motion planning: Maneuver variants in structured environments[END_REF] state that the vehicle must be separated from the obstacle, either by a longitudinal distance L ν or laterally by W ν . For illustration, we consider two identical obstacles centered at (80, 1.5) and (160, 3.5). Parameters are T = 15 s, τ = 1 s, L ν = 10 m and W ν = 2 m. The trajectory labeled "opt" in Fig. 2a is the global optimum found by the MIQP planner; for comparison purposes, we also plot all locally optimal trajectories "nopt-1,2,3" with their respective costs. We observe that our method can effectively find the globally optimal trajectory. Fig. 2b and 2c respectively illustrate the longitudinal and lateral speed and acceleration profiles of the ego vehicle.

C. Overtaking in a two-lane road

This case study considers an overtaking scenario on a twolane road with oncoming traffic. Fig. 3a shows the initial configurations of nearby vehicles, all driving with a constant speed of 10 m/s. This scenario is considered as difficult for both human drivers and autonomous vehicles [START_REF] Nilsson | Receding horizon maneuver generation for automated highway driving[END_REF], [START_REF] Murgovski | Predictive cruise control with autonomous overtaking[END_REF]. Reference [START_REF] Bender | The combinatorial aspect of motion planning: Maneuver variants in structured environments[END_REF] shows the existence of multiple homotopy classes in this scenario and proposes to exhaustively search for the globally optimal solution. Here, we show that we can find the globally optimal solution without explicitly enumerating all homotopy classes.

It is possible to model surrounding vehicles as rectangles as in the previous case study, thus requiring four integer variables for each vehicle and each time step k. However, by introducing the so called ramp barrier [START_REF] Nilsson | Receding horizon maneuver generation for automated highway driving[END_REF], [START_REF] Murgovski | Predictive cruise control with autonomous overtaking[END_REF], the problem can be further simplified by approximating the rectangular obstacle region by a triangle only using two linear constraints as shown in Fig. 3a.

Let ν be a surrounding vehicle; if ν is in the same lane as the ego vehicle, the constraints are given as

δ ν (k) = 0 ⇔ - x(k) -x ν (k) L ν + y(k) -y ν (k) W ν ≥ 1, (11a)
δ ν (k) = 1 ⇔ x(k) -x ν (k) L ν + y(k) -y ν (k) W ν ≥ 1, (11b)
where L ν and W ν are minimal longitudinal and lateral separations during lane change. Similarly, the constraints for an oncoming vehicle ν can be modeled as Using this formulation, only one integer variable δ ν is required per vehicle and per time step.

δ ν (k) = 0 ⇔ x(k) -x ν (k) L ν + y(k) -y ν (k) W ν ≤ -1, (12a)
δ ν (k) = 1 ⇔ - x(k) -x ν (k) L ν + y(k) -y ν (k) W ν ≤ -1. (12b)
x(m) v x (m/s) -2 0 2 4 v y (m/s) v x v y (c) 0 50 100 150 200 250 -3 0 3 x(m) a x (m/s 2 ) -1 0 1 a y (m/s 2 ) a x a y (d)
In the simulation, we adopt a prediction horizon of T = 15 s and we let τ = 1 s. To allow the ego vehicle to temporarily cross the lane border, the upper limit y is relaxed. Fig. 3b illustrates the trajectory of overtaking. We observe that the ego vehicle decides to accelerate slightly so that it can use the space between the first oncoming vehicle and the second oncoming vehicle to perform the overtaking.

For comparison purposes, we reduce the penalty on speed deviation q 1 to 0.5 in the cost function [START_REF] Gao | Model Predictive Control for Autonomous and Semiautonomous Vehicles[END_REF]. The resulting trajectory is shown in Fig. 4: in this case, the ego vehicle chooses not to overtake, as the cost of this maneuver is higher than that of following the slower car, due to the small penalty on the speed deviation. This demonstrates the flexibility of the MIQP formulation: different driving styles can be configured simply by modulating the weighting terms.

D. Lane change

This last case study considers the problem of decision making and trajectory generation for a lane change maneuver: the ego vehicle must decide the objective lane as well as the optimal trajectory to reach this lane, without colliding with surrounding vehicles. The complexity of this problem lies in the multiple discrete choices raised from multiple lanes and multiple vehicles on each lane. References [START_REF] Nilsson | Strategic decision making for automated driving on two-lane, one way roads using model predictive control[END_REF], [START_REF] Du | Autonomous lane-change controller via mixed logical dynamical[END_REF] have considered this problem using MILP formulations; however, their modeling cannot ensure that trajectories are dynamically feasible due to important simplifications of the vehicle dynamics.

We consider a road with N lanes, labeled by γ ∈ {1, ..., N }. We introduce a binary variable δ γ (k) that equals 1 if the ego vehicle is on lane γ at time step k. Let V be the set of surrounding vehicles and V γ be the set of vehicles inside lane γ. We introduce the following logical constraint: ∀k > 0,

δ γ (k) = 1 ⇔ y r (k) = y γ ∧ y(k) ∈ [y γ , y γ ] , (13)
such that, if the ego vehicle is in lane γ, then the vehicle should be within the boundary of lane γ and the reference centerline should be set to the centerline of the lane. Moreover, we add the following collision avoidance constraints: ∀k ≥ 0,

δ γ (k) = 1 ⇔ ∀ν ∈ V γ , δ ν (k) = 1 ⇔ x(k) ≤ x ν (k) -L ν , δ ν (k) = 0 ⇔ x(k) ≥ x ν (k) + L ν , (14)
such that the ego vehicle must avoid collisions with all the vehicles in lane γ.

The ego vehicle is only allowed to be in one lane at any given time, thus we add the following constraint: ∀k ≥ 0,

N γ=1 δ γ (k) = 1. (15)
Fig. 5a shows a highway with three lanes. The ego vehicle starts in the rightmost lane with a speed of 15 m/s. Surrounding vehicles are distributed over three lanes. The vehicle on the leftmost lane drives at a speed of 15 m/s while other surrounding vehicles drive at a speed of 10 m/s. Constraints ( 13), ( 14) and ( 15) are enforced along with other constraints on the formulated MIQP problem. The horizon is set to T = 15 s and the time step τ = 1 s. The constraint on lane boundary is temporarily deactivated. We observe in Fig. 5b that the ego vehicle chooses the left-most lane as the objective lane and plans a dynamically feasible and collisionfree trajectory to reach the lane within the prediction horizon.

Table II presents the computation times for the four case studies, demonstrating the real-time capability of the proposed formulation; future work will focus on further reducing this computation time.

IV. DISCUSSION AND CONCLUSION

We have addressed the problem of optimal trajectory generation integrating both continuous and logical constraints using an MIQP based approach. This approach is wellsuited to the on-road autonomous driving environment as logical constraints can arise from different aspects of driving: traffic rules, on-road obstacles and the existence of multiple maneuver variants. We have showcased the universality and the efficiency of the method by applying it to various challenging driving scenarios. This paper uses a third-order linear vehicle model for trajectory generation. This model is suitable if the longitudinal motion dominates the lateral one, for example when driving on highways or on urban arterial roads. However, for application to low-speed driving, a nonlinear vehicle model might be more desirable. Future work will investigate the applicability of combining the feedback linearization with the MIQP formulation. Another important assumption in this paper is that the vehicle is driving on road segments with small curvature. For large road curvatures, the current model can be imprecise. Future work will investigate this issue.

In the future, we will also integrate the proposed method into a receding horizon framework so that the ego vehicle re-plans regularly to incorporate new information. The algo-rithm will be implemented in the AnnieWay [START_REF] Kammel | Team annieway's autonomous system for the 2007 darpa urban challenge[END_REF] autonomous vehicle to study the performance of our approach in a realworld setting.

2

  Philipp Bender and Christoph Stiller are with Department of Measurement and Control, Karlsruhe Institute of Technology.

Fig. 1 :

 1 Fig. 1: Speed bump scenario: (a) illustration of the scenario, (b) longitudinal speed profile w.r.t. the longitudinal offset.

Fig. 2 :

 2 Fig.2: Obstacle avoidance scenario: (a) the globally optimal trajectory and locally optimal trajectories (respective costs: opt 3.16, nopt-1 36.5, nopt-2 37.2, nopt-3 84.67), (b) speed profiles, (c) acceleration profiles. Obstacles are marked in red.

Fig. 3 :

 3 Fig. 3: Overtaking scenario: (a) illustration of the scenario and the ramp barrier methods, (b) the overtaking trajectory, (c) speed profiles, (d) acceleration profiles. The position of the ego vehicle is shown as blue dots, and rectangles with different colors mark the positions of surrounding vehicles. Time dimension is color-coded, with lighter colors corresponding to instants further away in time.

Fig. 4 :

 4 Fig.4: Effect of reducing the speed deviation weight to q1 = 0.5 in the overtaking scenario. The position of the ego vehicle is shown as blue dots, and rectangles with different colors mark the positions of surrounding vehicles. Time dimension is color-coded, with lighter colors corresponding to instants further away in time.

Fig. 5 :

 5 Fig. 5: Lane change scenario: (a) illustration of the scenario, (b) trajectory of the ego vehicle and the surrounding vehicles, (c) speed profiles, (d) acceleration profiles. The position of the ego vehicle is shown as blue dots, and rectangles with different colors mark the positions of surrounding vehicles. Time dimension is color-coded, with lighter colors corresponding to instants further away in time.

TABLE II :

 II Time statistics

Xiangjun Qian, Florent Altché and Arnaud de La Fortelle are with MINES ParisTech, PSL -Research University, Centre for Robotics.