Lisa Serir
email: lisa.serir@femto-st.fr

Emmanuel Ramasso
email: emmanuel.ramasso@femto-st.fr

Noureddine Zerhouni
email: noureddine.zerhouni@femto-st.fr

E2GK: Evidential Evolving Gustafsson-Kessel Algorithm For Data Streams Partitioning Using Belief Functions

A new online clustering method, called E2GK (Evidential Evolving Gustafson-Kessel) is introduced in the theoretical framework of belief functions. The algorithm enables an online partitioning of data streams based on two existing and efficient algorithms: Evidantial c-Means (ECM) and Evolving Gustafson-Kessel (EGK). E2GK uses the concept of credal partition of ECM and adapts EGK, offering a better interpretation of the data structure. Experiments with synthetic data sets show good performances of the proposed algorithm compared to the original online procedure.

Introduction

Given a set of N data points, clustering refers to a wide variety of algorithms that aim at discovering c groups (clusters) ω 1 , ..., ω c whose members are similar in some way. The purpose is to summarize the data or to verify an existing structure of the data. In most cases, a cluster is defined as a subset of data for which the similarity between data within this subset is larger than the similarity with the data in other subsets. In many cases, the Euclidean distance between data is used as a dissimilarity measure.

A wide variety of clustering methods has been developed. The most commonly used methods are divided into two main categories: hierarchical and nonhierarchical methods. Among the latter, the K-means algorithm [START_REF] Bishop | Neural networks for pattern recognition[END_REF] is the most commonly used. The idea of K-means algorithm is to randomly create K clusters and to assign each data point to the closest one in an iterative way, reallocating points until a convergence criterion is satisfied.

Using hard partitioning methods, data are grouped in an exclusive way, i.e., data can't belong to two (or more) different clusters. In fuzzy partitioning, each data can belong to more than one cluster with different membership degrees. The most popular fuzzy partitioning method is Bezdek's Fuzzy C-means (FCM) algorithm [START_REF] Bezdek | Pattern Recognition with fuzzy objective function algorithms[END_REF]. One can also mention the Gustafson-Kessel fuzzy clustering algorithm [START_REF] Gustafson | Fuzzy clustering with a fuzzy covariance matrix[END_REF] that is capable of detecting hyper-ellipsoidal clusters of different sizes and orientations by adjusting the covariance matrix of data.

Another concept of partition, introduced in [START_REF] Denoeux | Evclus: Evidential clustering of proximity data[END_REF], is the credal partition based on belief functions theory. A credal partition extends the existing concepts of hard, fuzzy (probabilistic) and possibilistic partition by allocating, for each data, a mass of belief, not only to single clusters, but also to any subset of Ω = {ω 1 , ..., ω c }. This particular representation allows coding all the situations, from certainty to total ignorance of membership to clusters. In the Evidential c-Means (ECM) algorithm [START_REF] Masson | ECM: An evidential version of the fuzzy c-means algorithm[END_REF], the credal partition is in particular exploited for outliers detection.

Online clustering is an important problem that frequently arises in many fields, such as pattern recognition and machine learning [START_REF] Duda | Pattern Classification[END_REF]. Numerous techniques have been developed for clustering data in a static environment [START_REF] Bishop | Neural networks for pattern recognition[END_REF]. However, in many real-life applications, non-stationary data (i.e., with time-varying parameters) are commonly encountered. The task of online clustering is to group incoming data into clusters in a temporal sequence. Also called incremental clustering in machine learning [START_REF] Janichen | Acquisition of concept description by conceptual clustering[END_REF], online clustering, is generally unsupervised and has to manage recursive training in order to incorporate new information gradually and to take into account model evolutions over time.

In this paper, we propose the Evidential Evolving Gustafson Kessel algorithm (E2GK) which permits to adapt a credal partition matrix as data gradually arrive. This clustering algorithm is introduced in the theoretical framework of belief functions, and more precisely of Smets' Transferable Belief Model (TBM, [START_REF] Smets | The transferable belief model[END_REF]). E2GK is composed of two main steps, both performed online:

1. Determination of clusters' prototypes (also called centers), either by moving existing prototypes or by creating new ones. To do so, we use some results from the Evolving Gustafson-Kessel algorithm (EGK) proposed in [START_REF] Georgieva | Gustafson-kessel algorithm for evolving data stream clustering[END_REF]. 2. Allocation of the belief masses to the different subsets of classes. This step is based on some results of the Evidential c-means algorithm (ECM) [START_REF] Masson | ECM: An evidential version of the fuzzy c-means algorithm[END_REF].

E2GK benefits from two efficient algorithms: EGK and ECM, by dealing with -in an online manner -doubt between clusters and outliers. Doubt is generally encountered in data transition and can be useful to limit the number of clusters in the final partition. Moreover, outliers are well managed using the conflict degree explicitly emphasized in the TBM framework.

In Section 2, we present GK and ECM algorithms as well as some tools of the theory of belief functions giving the necessary background for Section 3 in which we introduce E2GK. Some results are finally presented in Section 4.

Background

Let the data be in the form of a collection {x 1 , . . . , x k , . . . , x N } of feature vectors x k ∈ q , and c the number of clusters, each of them characterized by a prototype (or a center) v i ∈ q .

Gustafson-Kessel Algorithm

Clustering algorithms based on an optimization process aim at minimizing a suitable fuction J that represents the fitting error of the clusters regarding the data:

J(V, U ) = c i=1 N k=1 (u ik ) β d 2 ik , (1)
where u ik is the membership degree of point k to the i-th prototype (cluster center),

-U = [u ij ] is the resulting partition matrix with dimension c × N , -V = [v i ]
is the c × q matrix of prototypes, -d ik is the distance between the k-th data point x k and the i-th prototype, -Paramater β > 1 is a weighting exponent that controls the fuzziness of the partition (it determines how much clusters may overlap).

The distance d ik used in the GK algorithm is a squared inner-product distance norm (Mahalanobis) that depends on a positive definite symmetric matrix A i defined by:

d 2 ik = x k -v i 2 Ai = (x k -v i )A i (x k -v i ) T . (2)
This adaptive distance norm is unique for each cluster as the norm inducing matrix A i , i = 1...c, is calculated by estimates of the data covariance

A i = [ρ i det(F i )] 1/q F -1 i , (3)
where ρ i is the cluster volume of the i-th cluster and F i is the fuzzy covariance matrix calculated as follows:

F i = N k=1 (u ik ) β (x k -v i ) T (x k -v i ) N k=1 (u ik ) β . ( 4
)
The objective function is minimized using an iterative algorithm, which alternatively optimizes the cluster centers and the membership degrees:

v i = N k=1 (u ik ) β x k N k=1 (u ik ) β , i = 1 • • • c, k = 1 • • • N , (5)
and

u ik = 1 c j=1 (d ik /d jk ) 2/β-1 , i = 1 • • • c, k = 1 • • • N . (6)
The GK algorithm has the great advantage to adapt the clusters according to their real shape.

Belief Functions and Credal partition

Dempster-Shafer theory of evidence, also called belief functions theory, is a theoretical framework for reasoning with partial and unreliable information. It was first introduced by A. P. Dempster (1968), then developed by G. Shafer (1976). Later, Ph. Smets proposed a general framework, the Transferable Belief Model (TBM) [START_REF] Smets | The transferable belief model[END_REF], for uncertainty representation and combination of various pieces of information without additional priors. Considering a variable ω taking values in a finite set called the frame of discernment Ω, the belief of an agent in subsets of Ω can be represented by a basic belief assignment (BBA), also called belief mass assignment:

m : 2 Ω → [0, 1] A → m(A) , (7)
with A⊆Ω m(A) = 1. A belief mass can not only be assigned to a singleton (|A| = 1), but also to a subset (|A| > 1) of variables without any assumption concerning additivity. This property permits the explicit modeling of doubt and conflict, and constitutes a fundamental difference with probability theory. The subsets A of Ω such that m(A) > 0, are called the focal elements of m. Each focal element A is a set of possible values of ω. The quantity m(A) represents a fraction of a unit mass of belief allocated to A. Complete ignorance corresponds to m(Ω) = 1, whereas perfect knowledge of the value of ω is represented by the allocation of the whole mass of belief to a unique singleton of Ω, and m is then said to be certain. In the case of all focal elements being singletons, m boils down to a probability function and is said to be bayesian.

A positive value of m(∅) is considered if one accepts the open-world assumption stating that the set Ω might not be complete, and thus ω might take its values outside Ω. This value represents the degree of conflict and is then interpreted as a mass of belief given to the hypothesis that ω might not lie in Ω. This interpretation is useful in clustering for outliers detection [START_REF] Masson | ECM: An evidential version of the fuzzy c-means algorithm[END_REF].

Belief functions theory is largely used in clustering and classification problems [START_REF] Denoeux | A k-nearest neighbor classification rule based on dempster-shafer theory[END_REF][START_REF] Kim | Evidential reasoning approach to multisource-data classification in remote sensing[END_REF]. Recently (2003) was proposed the use of belief functions for cluster analysis. Similar to the concept of fuzzy partition but more general, the concept of Credal Partition was introduced. It particularly permits a better interpretation of the data structure. A credal partition is constructed by assigning a BBA to each possible subset of clusters. Partial knowledge regarding the membership of a datum i to a class j is represented by a BBA m ij on the set Ω = {ω 1 , . . . , ω c }. This particular representation makes it possible to code all situations, from certainty to total ignorance.

A ∅ ω1 ω2 {ω1, ω2} ω3 {ω1, ω3} {ω2, ω3} {ω1, ω2, ω3} m1(A) 0 1 0 0 0 0 0 0 m2(A) 0 0 0 0 0 0 0 1 m3(A) 0 0 0 0 0.2 0.5 0 0.3 m4(A) 1 0 0 0 0 0 0 0 2.3 ECM: Evidential C-Means algorithm
Our approach for developing E2GK (Evidential Evolving GK algorithm) is based on the concept of credal partition as described in ECM [START_REF] Masson | ECM: An evidential version of the fuzzy c-means algorithm[END_REF] where the objective function was defined as:

J ECM (M, V ) = N k=1 {i/Ai =∅,Ai⊆Ω} |A i | α m β ki d 2 ki + N k=1 δ 2 m k (∅) β , (8)
subject to {i/Ai =∅,Ai⊆Ω}

m ki + m k (∅) = 1 ∀k = 1, . . . , N , (9)
where:

α is used to penalize the subsets of Ω with high cardinality, -β > 1 is a weighting exponent that controls the fuzziness of the partition, -d ki denotes the Euclidean distance between datum k and prototype v i , -δ controls the amount of data considered as outliers.

The N × 2 c partition matrix M is derived by determining, for each datum k, the BBAs m ki = m k (A i ) , A i ⊆ Ω such that m ki is low (resp. high) when the distance d ki between datum k and focal element A i is high (resp. low). The matrix M is computed by the minimization of criterion [START_REF] Duda | Pattern Classification[END_REF] and was shown to be [START_REF] Masson | ECM: An evidential version of the fuzzy c-means algorithm[END_REF], ∀k = 1 . . . N , ∀i/A i ⊆ Ω, A i = ∅:

m ki = |A i | -α/(β-1) d -2/(β-1) ki A l =∅ |A l | -α/(β-1) d -2/(β-1) kl + δ -2/(β-1) , (10)
and

m k (∅) = 1 -Ai =∅ m ki .
The distance between a datum and any non empty subset A i ⊆ Ω is then defined by computing the center of each subset A i . The latter is the barycenter v i of the clusters' centers (obtained by minimizing criterion (8)) composing A i .

Deriving E2GK

GK algorithm [START_REF] Gustafson | Fuzzy clustering with a fuzzy covariance matrix[END_REF] has the great advantage to adapt the clusters according to their real shape. The resulting clusters are hyper-ellipsoids with arbitrary orientation and are well suited for a variety of practical problems. However, GK is not able to deal with streams of data (relies on an iterative optimization scheme). Moreover, it assumes that the number of clusters is known in advance.

In [START_REF] Georgieva | Gustafson-kessel algorithm for evolving data stream clustering[END_REF], an online version of GK clustering algorithm (EGK) was developed to enable online partitioning of data streams based on a similar principle to the one used in the initial GK algorithm [START_REF] Gustafson | Fuzzy clustering with a fuzzy covariance matrix[END_REF]. In particular, online updating of the fuzzy partition matrix relies on the same formula [START_REF] Denoeux | A k-nearest neighbor classification rule based on dempster-shafer theory[END_REF]. Rules were then proposed to decide whether a new cluster has to be created or existing prototypes should evolve.

E2GK: Evidential Evolving Gustafsson-Kessel algorithm

The adaptation of the EGK algorithm to belief functions is introduced in this section. The E2GK algorithm is presented in Tab. 2. It relies on some parts developed in [START_REF] Georgieva | Gustafson-kessel algorithm for evolving data stream clustering[END_REF] and the proposed adaptations are emphasized in bold characters.

Step 1 -Initialization: At least one cluster's center should be provided. Otherwise, the first point is chosen as the first prototype. If more than one prototype is assumed in the initial data, GK or ECM algorithm can be applied to identify an initial partition matrix. The result of the initialization phase is a set of c prototypes v i and a covariance matrix1 F i .

Step 2 -Decision making: The boundary of each cluster is defined by the cluster radius r i , defined as the medium distance between the cluster center v i and the points belonging to this cluster with membership degrees larger or equal to a given threshold u h :

r i = median ∀xj ∈ i-th cluster and Pji>u h v i -x j Ai . ( 11
)
where P ij is the confidence degree that point j belongs to ω i ∈ Ω and can be obtained by three main processes: either by using the belief mass m j (ω i ), or the pignistic transformation [START_REF] Smets | The transferable belief model[END_REF] that converts a BBA into a probability distribution, or by using the plausibility transform [START_REF] Cobb | On the plausibility transformation method for translating belief function models to probability models[END_REF]. We propose here to choose the pignistic transformation. The median value is used (instead of the maximum rule in EGK) to reduce the sensitivity to extreme values. Moreover, the minimum membership degree u h -initially introduced in [START_REF] Georgieva | Gustafson-kessel algorithm for evolving data stream clustering[END_REF] and requiring to decide whether a data point belongs or not to a cluster -can be difficult to assess. It may depend on the density of the data as well as on the level of cluster overlapping. We rather set u h automatically to 1/c in order to reduce the number of parameters while ensuring a natural choice for its value.

Step 3 -Computing the partition matrix: Starting from the resulting set of clusters at a given iteration, we build the partition matrix M (10) using the Mahalanobis distance (2)(3). We assumed that each cluster volume ρ i = 1 as in standard GK algorithm.

Step 4 -Adapting the structure: Given a new data point x k , two cases are considered: Case 1: x k belongs to an existing cluster, thus a clusters' update has to be performed. Data point x k is assigned to the closest cluster p if d pk ≤ r p . Then, the p-th cluster is updated:

v p,new = v p,old + θ • (x k -v p,old ) , (12)
and

F p,new = F p,old + θ • (x k -v p,old ) T (x k -v p,old ) -F p,old , (13)
where θ is a learning rate, v p,new and v p,old denote respectively the new and old values of the center, and F p,new and F p,old denote respectively the new and old values of the covariance matrix. Case 2: x k is not within the boundary of any existing cluster (i.e. d pk > r p ), thus a new cluster may be defined and a clusters' update has to be performed. The number of clusters is thus incremented: c = c + 1. Then, the incoming data x k is accepted as the center v new of the new cluster and its covariance matrix F new is initialized with the covariance matrix of the closest cluster F p,old .

In the initial EGK algorithm [START_REF] Georgieva | Gustafson-kessel algorithm for evolving data stream clustering[END_REF], a parameter P i was introduced to assess the number of points belonging to the i-th cluster. The authors suggested a threshold parameter P tol to guarantee the validity of the covariance matrices and to improve the robustness. This (context-determined) parameter corresponds to the desired minimal amount of points falling within the boundary of each cluster. The new created cluster is then rejected if it contains less than P tol data points.

After creating a new cluster, the data structure evolves. However, the new cluster may contain data points previously assigned to another cluster. Thus, the number of data points in previous clusters could change. We propose an additional step to verify, after the creation of a new cluster, that all clusters have at least the required minimum amount of data points (P tol or more). If not, the cluster with the lowest number of points is deleted. Therefore, compared to the initial EGK algorithm, in which the number of clusters only increases, E2GK is more flexible because the structure can change either by increasing or decreasing the number of clusters.

The overall algorithm is presented in Tab. 2 where the proposed adaptation appears in bold.

Application of E2GK

To illustrate the ability of the proposed algorithm, let consider the following synthetic data randomly generated from five different bivariate gaussian distributions with parameters as given in Tab. 3.

Initial clusters (Fig. 1) of N = 15 data points each, of type G 1 and G 2 , were identified by batch GK procedure with u h = 0.5, P tol = 20 and θ = 0.1. To test the updating procedure, we gradually (one point at a time) added the following data points (in this given order): 1) 15 data points of type G 1 , 2) 15 data points Each new incoming data point leads to a new credal partition. Figure 2 shows the final resulting partition. The center of gravity of each cluster is marked by a big star (the notation ω ij stands for {ω i , ω j }). A data point falling in a subset ω ij means that this point could either belong to ω 1 or ω 2 . The points represented in circles are those with the highest mass given to the empty set and considered as outliers. It can be seen that a meaningful partition is recovered and that outliers are correctly detected.

The online adaptation of the clusters is illustrated in Figure 3. One can see how E2GK assigns each new data point to the desired cluster or subset. The figure depicts the evolution of the partition regarding the order of arrival of the data (like mentionned before). The first 30 points are used to initialize clusters ω 1 and ω 2 . Then, from t = 31 to 45 points are assigned by E2GK to cluster ω 2 . The next 15 points are assigned to ω 1 then to ω 4 , ω 3 (30 points) and to ω 4 . The next points correspond to noise and are mainly assigned to subsets, for example point 160 to ω 134 . Figure 4 also depicts the structure evolution, that is the number of clusters at each instant. The scenario given at the begining of this section is recovered: at t = 76 data from group G 3 arrive but still, not enough data are available to create clustrs while a cluster is created at t = 93 and t = 110 for group G 4 and G 3 respectively. "Noise" and atypical points arriving from t = 181 to t = 211 do not affect the structure. This figure does not illustrate clusters' removing because this operation is made within the algorithm. Figure 5 describes the dataset partitioning after decision making by applying the pignistic transformation [START_REF] Smets | The transferable belief model[END_REF] on the final credal partition matrix. Datatips provide the center coordinates, which are close to the real parameters (Tab 3). In comparison, we also provide in Figure 6 the centers obtained by EGK algorithm with parameters P tol = 20, u h =1/c and θ = 0.1 (the same as in E2GK).

Conclusion

To our knowledge, only one incremental approach to clustering using belief functions has been proposed [START_REF] Ben-Hariz | IK-BKM: An incremental clustering approach based on intra-cluster distance[END_REF]. However, in this approach the number of clusters is known in advance so this is not adapted for online applications. Moreover, data are described by a given number of attributes, each labeled by a mass of belief provided by an expert. This prior information is generally not available in pattern recognition problems. E2GK algorithm, described in this paper, is an evolving clustering algorithm using belief functions theory, which relies on the credal partition concept. This type of partition allows a finer representation of datasets by emphasizing doubt between clusters as well as outliers. Doubt is important for data streams analysis from real systems because it offers a suitable representation of gradual changes in the stream. E2GK relies on some parts of EGK algorithm [START_REF] Georgieva | Gustafson-kessel algorithm for evolving data stream clustering[END_REF], initially based on a fuzzy partition, to which we bring some modifications:

using the median operator to calculate cluster radius (vs. max. for EGK), using the credal partitioning (vs. fuzzy for EGK), changing the partitionning structure by adding or removing clusters (vs.

adding only in EGK).

Simulation results show that E2GK discovers relatively well the changes in the data structure. A thorough analysis of parameters' sensitivity (P tol and θ) is now required to properly and automatically set them.

Example 1 .

 1 Considering N = 4 data and c = 3 classes, Tab. 1 gives an example of a credal partition. BBAs for each datum in Tab. 1 illustrate various situations: datum 1 certainly belongs to class 1, whereas the class of datum 2 is completely unknown. Partial knowledge is represented for datum 3. As m 4 (∅) = 1, datum 4 is considered as an outlier, i.e., its class does not lie in Ω.

  5 2.5] 2 of type G 2 , 3) 15 data points of type G 3 , 4) 30 data points of type G 4 , 5) 15 data points of type G 3 , 6) 90 data points of type "noise", 7) 6 data points at the following positions: [10.1 3.2], [10.1 -3.2], [-4.1 -3.1], [-2.3 8.3], [8.6 -3.1] and [6.2 9.2]. E2GK parameters were set to: P tol = 20, θ = 0.1, δ = 10, α = 1 and β = 2.

1 Fig. 1 .

 11 Fig. 1. Initialization of E2GK algorithm using some data from two clusters. Centers are represented by stars.

Fig. 2 .

 2 Fig.2. Credal partition with δ = 10, α = 1, β = 2, θ = 0.1, P tol = 20. Big stars represent centers. We also displayed the centers corresponding to subsets, e.g. ω123, and atypical data (dots) are well detected.

Fig. 3 .

 3 Fig.3. Structure adaptation: a datum arrives at each instant (x-axis) and is assigned to one of all possible subsets (y-axis). The set of possible subsets also evolves with the number of clusters.

Fig. 4 .

 4 Fig. 4. Structure evolution: the number of clusters at each instant varies as data arrive.

Fig. 5 .

 5 Fig.5. Decision on clusters for each point based on the pignistic probabilities obtained from the credal partition (Fig.2) using E2GK algorithm. Also are displayed the coordinates of the centers found by E2GK.

5 Fig. 6 .

 56 Fig.6. Decision on clusters for each point based on the maximum degree of membership from the fuzzy partition using GK algorithm. Also are displayed the coordinates of the centers found by EGK. The parameter u h was set to 1/c and the other parameters are the same as in E2GK (θ = 0.1 and P tol = 20).

Table 1 .

 1 Example of a credal partition.

Table 2 .

 2 E2GK algorithmInitialization1. Take the first point as a center or apply the off-line GK or ECM algorithm to get the initial number of clusters c and the corresponding centers V and covariances Fi, i = 1 • • • c 2. Calculate vj, the barycenter of the clusters' centers composing Aj ⊆ Ω 3. Calculate the credal partition M , using (10) Updating Repeat for each new data point x k 4. Find the closest cluster p 5. Decision-making: Calculate the radius rp of the closest cluster using (11) with the median value If d pk ≤ rp 6. Update the center vp (12) 7. Update the covariance matrix Fp (13) else 8. Create a new cluster: vc+1 := x k ; Fc+1 := Fp end 9. Recalculate the credal partition M using (10) 10. Check the new structure: remove the cluster with the minimum number of data points if less than P tol

Table 3 .

 3 Parameters of the synthetic data

	type	µ	σ
	G1	[0 5] 0.3
	G2	[0 0] 0.3
	G3	[6 6] 0.6
	G4	[6 0] 0.6
	noise [2.

To obtain a covariance matrix from ECM, one can also use the Mahalanobis distance as proposed in[START_REF] Antoine | Cecm -adding pairwise constraints to evidential clustering[END_REF].


Emmanuel Ramasso

Sylvie Jullien

Parameter identification in Choquet Integral by the Kullback-Leibler divergence on continuous densities with application to classification fusion

Keywords: Information fusion, Fuzzy measures, Relative Entropy, Health assessment, Classification

Classifier fusion is a means to increase accuracy and decision-making of classification systems by designing a set of basis classifiers and then combining their outputs. The combination is made up by non linear functional dependent on fuzzy measures called Choquet integral. It constitues a vast family of aggregation operators including minimum, maximum or weighted sum. The main issue before applying the Choquet integral is to identify the 2 M -2 parameters for M classifiers. We follow a previous work by Kojadinovic and one of the authors where the identification is performed using an informationtheoritic approach. The underlying probability densities are made smooth by fitting continuous parametric and then the Kullback-Leibler divergence is used to identify fuzzy measures. The proposed framework is applied on widely used datasets.

Introduction

In most of pattern recognition tasks, a first step consists in extracting relevant features bringing information on classes of interest. Features are then transformed into membership degrees according to the different classes by entities called classifiers.

A classifier is a system that takes as inputs a Qdimensional vector X t = [x 1 . . . x Q ] (also called features or attributes) and generates a degree of confidence in the statement "X belongs to class ω j " for all classes in Ω = {ω 1 . . . ω K }.

Multiple Classifier Systems (MCS) [START_REF] Kuncheva | Decision templates for multiple classifier fusion[END_REF] are designed when complementary (and sometimes redundant) information sources (here classifiers) are used in order to improve classification accuracy and decision-making. MCS can be also viewed as information fusion systems where inputs are classifiers. MCS can take several forms and among them the parallel one which takes as input M × K partial degrees of confidence and generates an output, called global confidence degree, made up of K degrees of confidence (one for each class). We denote by φ m,j (X) the degree of confidence delivered by a classifier m ∈ {1 . . . M } for class ω j ∈ Ω given the observation X.

Usual combinations of classifier outputs include product, naive Bayes and decision templates among others [START_REF] Kuncheva | Decision templates for multiple classifier fusion[END_REF] but most of them can be used provided each output represents an independent source of information. However, the independance assumption is not always satisfied. To face this problem, an approach that considers interactions among classifier outputs such as fuzzy integrals and in particular the Choquet Integral [START_REF] Grabisch | The application of fuzzy integrals in multicriteria decision making[END_REF][START_REF] Grabisch | A new algorithm for identifying fuzzy measures and its application to pattern recognition[END_REF] can be used. The explicit interaction coefficients (in the 2-additive form) provide very interesting information on complementarity and redundancy of the fused data which can also be used for subset selection [START_REF] Grabisch | Fuzzy integral for classification and feature extraction[END_REF].

A fuzzy integral is a type of non-linear functional dependent on fuzzy measures which constitutes a vast family of aggregation operators including many widely used operators (minimum, maximum, weighted sum, ordered weighted sum and so on) [START_REF] Grabisch | Fuzzy aggregation of numerical preferences[END_REF]. In order to be combined by the Choquet Integral, the commensurability [START_REF] Kojadinovic | Unsupervised aggregation of commensurate correlated attributes by means of the choquet integral and entropy functionals[END_REF] of classifier outputs must be satisfied. That means the classifier outputs must be defined on the same measurement scale.

The combination of all partial confidence degrees provided by classifiers is thus made up by a Choquet Integral which is described in the next section.

Choquet capacities and Choquet Integral

Let the M classifiers (sources) be denoted by Θ = {θ 1 , θ 2 . . . θ M }. A fuzzy measure µ k for a given class ω k weighs the importance of a subset of sources S ⊆ Θ and is defined by [START_REF] Grabisch | The application of fuzzy integrals in multicriteria decision making[END_REF][START_REF] Grabisch | A new algorithm for identifying fuzzy measures and its application to pattern recognition[END_REF]:

µ k : 2 Θ → [0, 1] S → µ k (S) (1)
satisfying the following constraints:

• µ k (∅) = 0 and µ k (Θ) = 1 • S ⊆ T ⇒ µ k (S) ≤ µ k (T ) (monotonicity)
The fuzzy measure is said:

• additive when µ k (S ∪ T ) = µ k (S) + µ k (T ), ∀S, T ⊆ Θ/S ∩ T = ∅ (probability measure), • super-additive when µ k (S∪T ) ≥ µ k (S)+µ k (T ), ∀S, T ⊆ Θ/S ∩ T = ∅, • sub-additive when µ k (S ∪ T ) ≤ µ k (S) + µ k (T ), ∀S, T ⊆ Θ/S ∩ T = ∅.
In classification problems, the fuzzy measure is used in order to take into account interactions between sources. One fuzzy measure is tuned for each class and each discrete Choquet Integral aggregates the information provided by the sources as follows [START_REF] Grabisch | The application of fuzzy integrals in multicriteria decision making[END_REF][START_REF] Grabisch | A new algorithm for identifying fuzzy measures and its application to pattern recognition[END_REF]:

C(φ 1 , . . . , φ M ) = M i=1 (φ (i) -φ (i-1) ) • µ k (S (i) ) (2)
where µ k (S (i) ) is the importance of subset of sources S (i) = {θ (i) , ..., θ (M ) } and the value φ (i) is provided by source θ (i) . The notation (.) indicates a permutation of indices according to the values provided by the sources such as φ(1) ≤ φ(2) ≤ ... ≤ φ(M ) ≤ 1 (and by convention φ (0) = 0). The Choquet integral thus coincides to the weighted arithmetic mean when the fuzzy measure is additive. One approximation of Eq. 2 called 2-additive Choquet Integral is often used and consists in considering a 2-order additive capacity which takes into account both the weights of each source and the interaction between pairs. The weight ν i of a source θ i (for the detection of class ω k ) and the coefficient I ij of interaction between both sources θ i and θ j can be obtained from the fuzzy measure µ k by [START_REF] Grabisch | The application of fuzzy integrals in multicriteria decision making[END_REF][START_REF] Grabisch | A new algorithm for identifying fuzzy measures and its application to pattern recognition[END_REF]:

ν i = T ⊆Θ\i (M -|T | -1)!|T |! M ! × µ k (T ∪ {i}) -µ k (T ) (3a)
I ij = T ⊆Θ\{i,j} (M -|T | -2)!|T |! (M -1)! × µ k (T ∪ {i, j})- µ k (T ∪ {i}) -µ k (T ∪ {j}) + µ k (T ) (3b)
These parameters are interesting for interpreting the fuzzy measure and also to highlight which sources are important and how they interact. When interactions between two sources are positive, the sources are said complementary while they are said redundant when interactions are negative. The problem of Choquet Integral parameters identification was treated by several authors [START_REF] Grabisch | Fuzzy integral for classification and feature extraction[END_REF][START_REF] Grabisch | A review of methods for capacity identification in choquet integral based multi-attribute utility theory, applications of the R package[END_REF]. In the context of classification as considered here (where the classes are known), the method proposed by Grabisch [START_REF] Grabisch | A new algorithm for identifying fuzzy measures and its application to pattern recognition[END_REF] and called Heuristic Least Mean Square (HLMS) is often used. However it requires the global scores (the real output of the Choquet Integral) to perform the optimization of the fuzzy measure. Recently, two information theoritic methods based on entropy [START_REF] Kojadinovic | Estimation of the weights of interacting criteria from the set of profiles by means of information-theoretic functionals[END_REF][START_REF] Kojadinovic | Unsupervised aggregation of commensurate correlated attributes by means of the choquet integral and entropy functionals[END_REF] and on relative entropy [START_REF] Jullien | Identification of choquet integral's parameters based on relative entropy and applied to classification of tomographic images[END_REF] were proposed. The former is purely unsupervised and requires only degrees of confidence of classifiers while the latter requires the ground truth, i.e. the real class of each pattern. The relative entropy-based approach is supervised but requires less prior information than HLMS.

Identification of Choquet Integral parameters based on discrete relative entropy

A probabilistic view

Each fuzzy value µ k (S) expresses the relative importance of a subset S for distinguishing class ω k from the others [START_REF] Kojadinovic | Estimation of the weights of interacting criteria from the set of profiles by means of information-theoretic functionals[END_REF][START_REF] Kojadinovic | Unsupervised aggregation of commensurate correlated attributes by means of the choquet integral and entropy functionals[END_REF]. In order to identify them, the authors in [START_REF] Jullien | Identification of choquet integral's parameters based on relative entropy and applied to classification of tomographic images[END_REF] proposed to use the relative entropy, also called Kullback-Leibler divergence (KL) [START_REF] Kullback | On information and sufficiency[END_REF], which is a measure of divergence between two densities. It could be interpreted as the expected discrimination information between two hypotheses and thus appears very natural for the identification of fuzzy measures.

To compute KL, one needs first to compute:

• the distribution (say P S k ) of confidence degrees in class ω k conditional to class ω k ,

• and the distribution (say P S k ) of confidence degrees in class ω k conditional to the other classes (ω k = Ω\ω k ), both given a subset of sources S. These distributions characterize the input data (confidence degrees) and the greater is the difference (calculated by KL) between them, the higher is the discrimination power.

Identifying fuzzy measure using a probabilistic approach was introduced in [START_REF] Kojadinovic | Estimation of the weights of interacting criteria from the set of profiles by means of information-theoretic functionals[END_REF][START_REF] Kojadinovic | Unsupervised aggregation of commensurate correlated attributes by means of the choquet integral and entropy functionals[END_REF] where the author proposed an unsupervised entropy-based method. When the class is known for each input pattern, the KL-based approach proposed in [START_REF] Jullien | Identification of choquet integral's parameters based on relative entropy and applied to classification of tomographic images[END_REF] should be used. It fully exploits the available information provided by the training dataset and, as expected, increases the discrimination power.

Relative entropy

We assume all confidence degrees to be commensurable values in [0, 1] which is generally true in classification.

Let

P Θ k (Y ) with Y = (φ 1,k (X), φ 2,k (X) . . . φ M,k (X)) ∈ [0, 1] |Θ| (resp. P Θ k (Y )) be the probability that classifiers 1, 2 . . . M jointly provide the values φ 1,k (X),φ 2,k (X),. . . and φ M,k (X) given the ground truth is class ω k (resp.
given ω k ) and observation X. In [START_REF] Jullien | Identification of choquet integral's parameters based on relative entropy and applied to classification of tomographic images[END_REF], the distributions were assumed discrete and the relative entropy (KL) of both distributions was thus given by:

D(P Θ k ||P Θ k ) = Y P Θ k (Y ) log P Θ k (Y ) P Θ k (Y ) (4)
For the sake of simplicity, we will denote by R k (S) the KL value given by D(P S k ||P S k ) for a given subset of sources S ⊆ Θ:

R k (S) ≡ D(P S k ||P S k ) (5)
Note that when the distributions P Θ k and P Θ k are computed, the distributions P S k and P S k for S ⊂ Θ are obtained by marginalizing out the components θ ∈ Θ, θ / ∈ S. To compute Eq. 4, the support of the distribution P S k must be included in the support of the distribution P S k otherwise the relative entropy diverges towards infinity. In order to respect this constraint, the skew divergence was used in [START_REF] Jullien | Identification of choquet integral's parameters based on relative entropy and applied to classification of tomographic images[END_REF].

From relative entropy to Choquet capacities

The relative entropy has to satisfy the conditions presented Section 2 in order to be interpreted as a Choquet capacity. For that, the relative entropy R k (S) for a subset S is normalized as in Kojadinovic's method [START_REF] Kojadinovic | Estimation of the weights of interacting criteria from the set of profiles by means of information-theoretic functionals[END_REF][START_REF] Kojadinovic | Unsupervised aggregation of commensurate correlated attributes by means of the choquet integral and entropy functionals[END_REF] by the entropy of the whole set of sources R k (Θ):

µ k (S) = R k (S) R k (Θ) (6)
Moreover, the relative entropy is zero when the set S is empty but also when both distributions are identical. Therefore, a source that provides the same degrees of support for a sought-after class ω k and for the other classes ω k is assigned a low importance value since it can not distinguish class ω k from the others. This is exactly the mean of discrimination power.

The relative entropy has also to satisfy the monotonicity constraint (Section 2), i.e. given two sources θ i and θ j , the relative entropy has to satisfy the following equations:

µ({θ i , θ j }) ≥ µ({θ i }) (7a) µ({θ i , θ j }) ≥ µ({θ j }) (7b)
In order to check these constraints, one can rewrite the relative entropy as [START_REF] Cover | Elements of information theory[END_REF][START_REF] Kojadinovic | Estimation of the weights of interacting criteria from the set of profiles by means of information-theoretic functionals[END_REF][START_REF] Kojadinovic | Unsupervised aggregation of commensurate correlated attributes by means of the choquet integral and entropy functionals[END_REF]:

R k ({θ i , θ j }) = R k ({θ i }) + R k ({θ j | θ i }) (8)
that is always positive and therefore has a monotonic behavior [START_REF] Cover | Elements of information theory[END_REF][START_REF] Kojadinovic | Estimation of the weights of interacting criteria from the set of profiles by means of information-theoretic functionals[END_REF][START_REF] Kojadinovic | Unsupervised aggregation of commensurate correlated attributes by means of the choquet integral and entropy functionals[END_REF][START_REF] Jullien | Identification of choquet integral's parameters based on relative entropy and applied to classification of tomographic images[END_REF]:

R k ({θ i , θ j }) ≥ R k ({θ i }) (9a) R k ({θ i , θ j }) ≥ R k ({θ j }) (9b)
This reasoning can be extended easily to larger subsets of sources. Therefore, the normalized relative entropy satisfies all the constraints in order to be interpreted as a Choquet capacity.

Modeling positive and negative interactions

When sources θ i and θ j , that provide distributions P S k and P S k , are independent, the relative entropy has an additive behavior [START_REF] Cover | Elements of information theory[END_REF][START_REF] Kojadinovic | Estimation of the weights of interacting criteria from the set of profiles by means of information-theoretic functionals[END_REF][START_REF] Kojadinovic | Unsupervised aggregation of commensurate correlated attributes by means of the choquet integral and entropy functionals[END_REF]:

R k ({θ i , θ j }) = R k ({θ i }) + R k ({θ j }) (10)
When sources θ i and θ j are interacting one each other, the relative entropy can be expressed by:

R k ({θ i , θ j }) = R k ({θ i }) + R k ({θ j }) + R k ({θ j |θ i }) -R k ({θ j }) (11)
where the last term can be negative or positive according to sources θ i and θ j implying that the identified Choquet capacities can be super-additive or sub-additive. Therefore the proposed method is able to model and identify both positive and negative interactions whereas Kojadinovic's approach can only identify negative ones.

Identification of Choquet Integral parameters based on continuous relative entropy

On using a continuous approach

The core of the KL-based method is the evaluation of the multidimensional probability distributions (P S k and P S k ). In [START_REF] Kojadinovic | Estimation of the weights of interacting criteria from the set of profiles by means of information-theoretic functionals[END_REF][START_REF] Kojadinovic | Unsupervised aggregation of commensurate correlated attributes by means of the choquet integral and entropy functionals[END_REF][START_REF] Jullien | Identification of choquet integral's parameters based on relative entropy and applied to classification of tomographic images[END_REF], the distributions were computed using discretization of confidence degrees (and histograms). We rather propose to remain in the continuous space (the space of the degrees of confidence) and to use parametric continuous densities for confidence degrees modelling. These densities allow to:

• Ensuring an infinite support for the distributions and therefore avoiding using artificial methods to solve the problem of minimum support. • Avoiding the necessity of finding the optimal number of bins for the histograms. This could be a serious problem for high dimensional data such as in image processing or in complex systems diagnosis. • Obtaining a more precise paving of the input space and therefore generating smooth distributions and improving the computation of the relative entropy by summing over more data points sampled from the continuous densities. • Simplifying the computation of marginalizations (according to the family of densities).

Fig. 1 depicts the problem of finding the number of bins for discrete histograms. We drew 100 points from a mixture of five Gaussians with equal probability and with means 0, 10, 25, 35, 50 and unit variance. We then computed histograms (depicted in the first three figures) with 10, 50 and 100 bins and we also run an EM in order to identify automatically the parameters of a continuous density made up of five components. The results are very different with a preference given to the last figure.

Modelling inputs by continuous densities

Modelling

We assume that the joint probability density function related to P Θ k (and similarly for P Θ k ) has a continuous and parametric form. For example, we consider mixtures of Gaussians which are very general and have interesting properties:

P Θ k (Y ) = L k a=1 c k,a • f k,a (Y ) (12)
where

f k,a (Y ) = N (Y |α k,a , Σ k,a ) (13)
with

Y = (φ 1,k (X), φ 2,k (X) . . . φ M,k (X)) ∈ [0, 1] |Θ| the joint
N (Y |α k,a , Σ k,a ) = exp -1 2 (Y -α k,a ) T Σ -1 k,a (Y -α k,a ) (2π) M/2 |Σ k,a | 1/2 (14)
The dimension of the parameters is the same as Θ which is M . The same expression holds for P Θ k :

P Θ k (Y ) = L b b=1 c k,b • g k,b (Y ) (15)
with different parameters indexed by subscripts (k, b).

Learning parameters of densities

The parameters of the densities can be estimated automatically by standard methods such as the Expectation-Maximization algorithm (EM [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF]) where L, the number of components, can also be estimated.

When the parameters of the distributions P Θ k and P Θ k have been specified, it is easy to compute P S k and P S k for subsets S ⊂ Θ by marginalization. In case the joint density related to P Θ k is represented by a mixture of multivariate Gaussians, the marginal is also a mixture of multivariate Gaussians where some components (marginalized out) have been eliminated. In particular, the |S| components of the mean vector of the marginal are the means of the variables in S and its covariance matrix is composed of the pairwise covariances of the same variables.

Continuous relative entropy

For two unimodal multivariate normal densities f k and g k (with L a = L b = 1), the KL has an exact closed form [START_REF] Hershey | Approximating the kullbackleibler divergence between gaussian mixture models[END_REF]:

R ex k (S; f k ; g k ) = 1 2 log |Σ k | |Σ k | + . . . Tr Σ -1 k Σ k -M + . . . µ k -µ k t Σ -1 k µ k -µ k (16)
When densities are multimodal, the continuous relative entropy is obtained by integrating on the support of P S k , Supp(P S k ) = {Y :

P S k (Y ) > 0}: R k (S) = Y ∈Supp(P S k ) P S k (Y ) log P S k (Y ) P S k (Y ) dY
(17) To evaluate this expression, several methods can be used [START_REF] Hershey | Approximating the kullbackleibler divergence between gaussian mixture models[END_REF]. In this paper, we have used Monte Carlo sampling (MC) and variational approximation (VA).

The MC method consists in drawing samples from the mixture associated to P S k . For that, a component is chosen randomly using the distribution c k,. . A continuous sample is then drawn from the associated Gaussian component and the density is evaluated. Given {Y i , i = 1 . . . N } the set of i.i.d. sampled points, we can approximate the integral (17) by its MC estimate:

R MC k (S) = 1 N i log P S k (Y i ) P S k (Y i ) → D(P S k ||P S k ) (18)
The precision of the evaluation of KL depends obviously on the number of simulations.

In the VA method, the integral is approximated by the following expression [START_REF] Hershey | Approximating the kullbackleibler divergence between gaussian mixture models[END_REF]:

R VA k (S) = La a=1 c k,a •log        La a =1 c k,a • e -R ex k (S;f k,a ;f k,a ) L b b=1 c k,b • e -R ex k (S;f k,a ;g k,b )        (19) where R ex k (S; f k,a ; g k,b
) is the exact value of KL between component a of f k and component b of g k given by Eq. 16.

Final algorithm

The overall algorithm for computing the fuzzy measure is as follows: R k (S) ← Apply Eq. 17 with P S k and P S k 9:

µ k (S) ← R k (S) R k (Θ) (Eq. 6)
10: end for From µ k , one can compute the weight of each source (i.e. classifier) and their interactions using Eq. 3b. These values can help an end-user or any people interested in knowing which classifiers contribute to the final results as well as how they interact.

Experiments

A toy example is first presented. Then, the proposed method is evaluated on two datasets from UCI [START_REF] Frank | UCI machine learning repository[END_REF]: vehicle and image segmentation. Classifiers used were the following: Evidential Neural Network (EvNN) [START_REF] Denoeux | An evidence-theoric neural network classifier[END_REF] (with 4 prototypes for each class), Evidential Nearest Neighborhood (EvKNN) [START_REF] Denoeux | A k-nearest neighbor classification rule based on Dempster-Shafer theory[END_REF] (with K = 5) and Support Vector Machines (SVM) [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF] (with a Gaussian Kernel of size 2.2). Classifiers were learnt using 1-vs-1 strategy for each class, and the final scores were obtained by using a weighted vote. SVM scores were transformed into probabilities using a sigmoid transfer function. Note that classifier parameters were not "optimized" for each dataset, since the goal is here to assess the fusion process. The KL was assessed using the MC method using 1e6 samples. Figure 3: A density P Θ k .

A toy example

As an example, let consider two classifiers θ 1 and θ 2 with confidence degrees in ω k , given the ground truth is ω k , being distributed according to Fig. 2, and according to Fig. 3 for confidence degrees in ω k given the ground truth is another class ω k . From these densities, we are looking for characterizing the importance of the coalition {θ 1 , θ 2 } in distinguishing ω k from ω k .

In Fig. 2, the confidence degrees of classifier 1 are globally close to unity given class ω k . That means classifier 1 often provides high scores for ω k given the ground truth is ω k . Classifier 2 however seems to provide some results close to 0.5 meaning classifier 2 is frequently not certain about the predicted class. Given the ground truth is ω k (Fig. 3), classifier outputs are globally close to 0 for ω k . That means classifiers generally provide low values for ω k when the ground truth is ω k as expected.

In order to quantify the importance µ k ({θ 1 , θ 2 }) of the coalition {θ 1 , θ 2 } given ω k , we compute the divergence between both distributions. The higher is the divergence, the higher is the importance of {θ 1 , θ 2 } for distinguishing ω k from the other classes. In this example, densities were obtained using two mixtures with the following parameters: With these parameters, Eq. 17 leads to R k ({θ 1 , θ 2 }) ≈ 9.44 (with N = 1.10 6 ).

• given ω k , α k,1 = [0.1, 0.1], Σ k,1 = [.

Vehicle dataset

The UCI's vehicle dataset is a four-classes problem composed of 946 examples almost uniformly distributed between classes. The goal is to classify data into one of the following types of vehicle: OPEL (ω 1 ), SAAB (ω 2 ), BUS (ω 3 ) and VAN (ω 4 ). The half of the dataset was used for classifier training, and the other half for testing. Also are displayed the Area Under the Curve (AUC) which reflects the efficiency in detecting the class. These curves can be compared with ROC curves of Figure 7 computed from the results of the fusion process proposed in this paper. Table 1 also gives the obtained fuzzy measures, while interaction and classifier weights computed from them (as detailed previously) are provided in Tables 2 and3.

ROC curves clearly show the complementarity of individual classifiers. For example, class ω 1 is better recognized using EvKNN (Fig. 5) with almost 98% of good classification. Class ω 2 is better recognized using SVM (Fig. 6) with an accuracy close to 83%. Class ω 3 is well recognized (about 93%) by EvNN (Fig. 4) and EvKNN (Fig. 5). Lastly, class ω 4 is better detected by SVM (Fig. 6). 1.

As shown in Figure 7, the proposed fusion process draw benefits of all these classifiers, providing AUCs close to 99%, 84%, 97% and 88% for class ω 1 , ω 2 , ω 3 and ω 4 respectively (improvement close to 10%). Interaction indexes can explain this result. Indeed, class ω 1 , that is well detected by all classifiers, is represented by a fuzzy measure with negative interactions because of redundancy. The highest redundancy is detected for class ω 2 between EvNN and EvKNN (I 12 = -0.20) while the highest complementarity is detected for class ω 3 between EvNN and SVM classifiers (I 23 = +0.33). Weights are also the highest for classifiers with the best accuracies, except for class ω 2 .

Image segmentation dataset

The UCI's image segmentation dataset is a sevenclasses problem composed of 210 training examples and 2100 for testing. Initially, features are 19dimensional but we reduced them to 6 dimensions and we kept features [1 3 4 8 9 17]. The goal is thus to classify data into one of the following types of vehicle: Brickface (ω 1 ), Sky (ω 2 ), Foliage (ω 3 ), Cement (ω 4 ), Window (ω 5 ), Path (ω 6 ) and Grass (ω 7 ).

For this dataset, we present the results in the form of confusion matrices ( Table 7 shows the confusion matrix of the fusion process result. This matrix clearly shows that the proposed method is able to draw benefits from individual classifiers. 9. For classes ω 1 , ω 3 and ω 5 interactions between classifiers emphasize complementarity. In particular, classifiers EvNN and EvKNN present the strongest complementarities. This is represented in confusion matrix where these classifiers mix sometimes several classes while classifier 3 confuses between class ω 7 and the others.

In general, an efficient classifier has also a relatively high weight (Tab. 10), and if it is not the case, interaction values provide compensation.

Conclusion

We proposed an information-theoritic approach relying on Kullback-Leibler divergence for fuzzy measures identification in the context of supervised classification. The use of well known parametric and continuous functions for the representation of confidence degrees allows to simplifying the estimation of joint densities and marginalization. We shown its application on widely used datasets where the fuzzy measure brought a lot of useful information concerning classifier importance and interactions. Results also emphasized that the proposed fusion process allows, on the one hand, to improve classification results and, on the other hand, to be robust to classifiers mistakes.

Further investigations concern the study of algorithms used for learning distribution parameters which are of key of importance.

Figure 1 :

 1 Figure 1: Role of the number of bins (see comments in the text).

  observation of degrees of confidence, c k,a the mixing coefficient of component a (among L k , with a c k,a = 1) and N (Y |α k,a , Σ k,a ) a multidimensional normal density with mean α k,a and covariance Σ k,a (positive definite):

5 :

 5 Require: L k the set of confidence degrees in class ω k of M classifiers given the ground truth is class ω k Require: L k the set of confidence degrees in class ω k of M classifiers given the ground truth are classes different from ω k Ensure: the fuzzy measure µ k for class ω k 1: P Θ k ← Estimate the parameters of the densities on L k 2: P Θ k ← Estimate the parameters of the densities on L k 3: µ k (Θ) ← 1, µ k (∅) ← 0 4: R k (Θ) ← Apply Eq. 17 with P Θ k and P Θ k for all S ⊂ Θ do

Figure 2 :

 2 Figure 2: A density P Θ k .

  01, 8; 8, .01], α k,2 = [0.2, 0.4] and Σ k,2 = [.02, .01; .01, .02] (with equal mixing coefficients c k,a ). • given ω k , α k,1 = [0.8, 0.8], Σ k,1 = [8, 1; 1, 8], α k,2 = [0.55, 0.9], Σ k,2 = [.02, 8; 8, .02], α k,3 = [0.95, 0.97], Σ k,3 = [4, 1; 1, 4] (with equal mixing coefficients c k,a ).

Figure 4 :

 4 Figure 4: EvNN classification results. Top: classes 1 and 2, bottom: classes 3 and 4.

Figures 4 -

 4 Figures 4-6 pictorially described ROC curves computed for each class given results of classifiers. Also are displayed the Area Under the Curve (AUC) which reflects the efficiency in detecting the class.These curves can be compared with ROC curves of Figure7computed from the results of the fusion process proposed in this paper. Table1also gives the obtained fuzzy measures, while interaction and classifier weights computed from them (as detailed previously) are provided in Tables2 and 3.ROC curves clearly show the complementarity of individual classifiers. For example, class ω 1 is better recognized using EvKNN (Fig.5) with almost 98%

Figure 5 :

 5 Figure 5: EvKNN classification results. Top: classes 1 and 2, bottom: classes 3 and 4.

Figure 6 :

 6 Figure 6: SVM classification results. Top: classes 1 and 2, bottom: classes 3 and 4.

Figure 7 :

 7 Figure 7: Fusion by Choquet Integral. Top: classes 1 and 2, bottom: classes 3 and 4.

Table 1 :

 1 The fuzzy measure for each class learnt by the proposed algorithm for the "vehicle" dataset.

	set	ω1	ω2	ω3	ω4
	1	0.22 0.39 0.23 0.24
	2	0.66 0.57 0.34 0.37
	12 0.75 0.70 0.44 0.65
	3	0.36 0.29 0.20 0.44
	13 0.50 0.60 0.43 0.62
	23 0.89 0.84 0.83 0.65
	set	ω1	ω2	ω3	ω4
	I12 -0.08 -0.20 -0.10 +0.10
	I13 -0.03 -0.02 +0.03 +0.01
	I23 -0.08 +0.03 +0.33 -0.10

Table 2 :

 2 Interaction values associated to the fuzzy measures of Table1.

	set	ω1	ω2	ω3	ω4
	ν1	0.15 0.26 0.19 0.28
	ν2	0.56 0.47 0.45 0.35
	ν3	0.29 0.28 0.37 0.37

Table 3 :

 3 Weight values associated to the fuzzy measure of Table

Table 4 -

 4 [START_REF] Grabisch | A review of methods for capacity identification in choquet integral based multi-attribute utility theory, applications of the R package[END_REF] where the ground truth is on columns while results of classifiers are

	class	ω1	ω2	ω3	ω4	ω5	ω6	ω7
	ω1	0	0	0	6	18	0	0
	ω2	67	264	0	67	15	0	0
	ω3	0	0	41	8	7	0	0
	ω4	233	36	259 219 256	0	0
	ω5	0	0	0	0	3	0	0
	ω6	0	0	0	0	1	300	1
	ω7	0	0	0	0	0	0	299

Table 4 :

 4 Confusion matrix of EvNN classifier (global acc.: 53.6%, degraded on purpose).

	class	ω1	ω2	ω3	ω4	ω5	ω6	ω7
	ω1	152	9	26	69	61	0	0
	ω2	38	223	0	55	17	0	0
	ω3	22	1	176	12	45	0	0
	ω4	45	55	13	140	6	30	0
	ω5	43	12	85	24	171	0	14
	ω6	0	0	0	0	0	270	0
	ω7	0	0	0	0	0	0	286

Table 5 :

 5 Confusion matrix of EvKNN classifier (global acc.: 67.5%)

	class	ω1	ω2	ω3	ω4	ω5	ω6	ω7
	ω1	98	0	1	3	4	0	0
	ω2	4	78	0	3	0	0	0
	ω3	5	0	59	0	0	0	0
	ω4	18	6	3	49	1	0	0
	ω5	0	0	14	0	94	0	0
	ω6	0	0	0	0	0	111	0
	ω7	175 216 223 245 201 189 300

Table 6

 6

		: Confusion matrix of SVM classifier (global
	acc.: 37.6%)
	on lines. Tables 4-6 are confusion matrices of indi-
	vidual classifiers. We here degraded on purpose the
	results of the first classifier (EvNN) on classes ω 1 ,
	ω 3 , ω 5 and ω 7 (by adding noise on the parameters
	trained by the algorithm). As a result, the confusion
	matrix presents low detection rate for these classes
	(close to 0 for some of them). We then applied the
	fusion process.
	class	ω1	ω2	ω3	ω4	ω5	ω6	ω7
	ω1	97	3	3	1	11	0	0
	ω2	55	234	0	58	15	0	0
	ω3	0	0	156	6	20	0	0
	ω4	148	63	62	231 113	3	0
	ω5	0	0	79	4	141	0	7
	ω6	0	0	0	0	0	297	0
	ω7	0	0	0	0	0	0	293

Table 7 :

 7 Confusion matrix after fusion using MC sampling (global acc.: 69%)

Table 8 :

 8 Table 8 is the confusion matrix obtained by the fusion process based on the variational approximation of the KL divergence which shows that, for this application, the approximation Confusion matrix after fusion using variational approximation (global acc.: 66%) is satisfying.Interaction values obtained in this application are shown in Table

	class	ω1	ω2	ω3	ω4	ω5	ω6	ω7
	ω1	95	3	8	3	25	0	0
	ω2	75	226	0	68	21	0	0
	ω3	5	0	145	6	25	0	0
	ω4	125	71	81	222 100	3	0
	ω5	0	0	66	1	129	0	5
	ω6	0	0	0	0	0	297	2
	ω7	0	0	0	0	0	0	293

Table 9 :

 9 Interaction values for application 2.

	class	ω1	ω2	ω3	ω4	ω5	ω6	ω7
	I12	0.32	-0.29 0.15	-0.22 0.03	-0.17 -0.28
	I13	0.06	0.00	0.06	-0.04 0.03	-0.02 -0.03
	I23	-0.09 0.07	0.02	-0.09 0.05	0.05	-0.02
	class ω1	ω2	ω3	ω4	ω5	ω6	ω7
	ν1	0.32 0.54 0.35 0.69 0.33 0.52 0.48
	ν2	0.54 0.28 0.50 0.24 0.50 0.33 0.48
	ν3	0.14 0.18 0.15 0.07 0.17 0.15 0.04

Table 10 :

 10 Weight values for application 2.


P Mondalek
email: pamela.mondalek@mines-paristech.fr

Prof L M Silva

M Bellet

Dr ] P Mondalek

Bellet

A numerical model for powder densification by SPS technique**

Spark Plasma Sintering belongs to a class of sintering techniques that employs electric current to assist compaction. This technology seems very promising to obtain net-shape components made of intermetallic alloys with a fine microstructure. However the SPS process is difficult to stabilize because of density heterogeneities arising from non homogeneous temperature or stress in the powder. This motivates the development of a three-dimensional finite element simulation in order to understand the distribution of current, temperature and porosity. The model couples three physical problems: electrical, thermal and mechanical. The numerical implementation is based on a monolithic formulation consisting in solving the different conservation equations on a single mesh including the specimen and the tooling. The general set of equations is described; the effect of the powder physical properties and of the geometry of the set up on the distribution of electrical current and temperature is discussed. A macroscopic Abouaf constitutive model is used to simulate powder densification; first results involving porosity evolution are presented.

Powder sintering technique assisted by electrical current, which is also called Spark Plasma Sintering (SPS), is a revolutionary high speed powder consolidation technology. This process presents many benefits such as reduction of sintering times and capacity of producing

materials with excellent properties. In fact, the speed of the process allows powder densification with nanostructured materials. The SPS is used to produce different materials like ceramics, composites, metals and also intermetallic alloys that are difficult to form with conventional techniques due to their different melting temperatures. [1] Application of SPS to high performance intermetallics should allow their cost-effective introduction in aero-engines, in which these low density materials have great potential to improve fuel efficiency and reduce the amount of pollution. [2] The main characteristic of SPS is that heat is generated internally by Joule effect, in contrast to conventional hot pressing, where heat is provided by conduction and/or radiation. More precisely, the process consists in applying a pulsed continuous current and a uniaxial pressure to obtain nanostructured materials due to very high heating rates and short sintering time (about 15 min). The powder material undergoes different mechanisms leading to full compaction, such as plastic deformation, sintering electromigration and surface diffusion induced by surface tension. During compaction, the particle size and current porosity value influence those different contributions. [3][4] In addition, there is a strong influence of temperature which is directly connected to the distribution of current. [5] Therefore, to control the microstructure homogeneity, especially for complex components, it is important to understand current flow, heat transfer, densification mechanisms and how they are coupled. A 3D numerical model has been developed using the CimLib library (implemented at CEMEF) to simulate the coupled electrical thermal mechanical problems. The finite element modeling and the equations are presented in the first section. The monolithic approach, on which the numerical code is based, is then presented.

Finally, results concerning the distribution of temperature and density are discussed.

Finite element modeling

Governing equations: coupled problems

The finite element modeling of the SPS process is based on a set of governing equations. [6- 7-8] The powder is placed inside a graphite mold and is surrounded by two graphite pistons on the top and bottom. An electrical current is then applied. Since all the tools are conductive, the electric current goes through the whole assembly and generates heat. Simultaneously, loads are applied to assist sintering. As a consequence, a compaction of the powder density is observed. Therefore, due to the interaction between the electric current, temperature and densification behaviour, three coupled problems need to be solved: electrical, thermal and mechanical.

First, the electrical modeling in the assembly (powder and tooling) is based on the charge conservation law:

0 = ⋅ ∇ J r (1)
Applying Ohm's law, the current density J r can be written:

U σ σ e e ∇ - = = E J r r (2)
( ) 0 = ∇ - ⋅ ∇ U σ e (3)
As a consequence of electric flow, heat is generated in the assembly by Joule effect. The associated source term in the heat equation is written:

E J v r ⋅ = e q
. The temperature distribution is calculated by solving the heat equation:

e p q dt dT ρc T k = + ∇ - ⋅ ∇ ) ( ( 4
)
where T is the temperature, k is the thermal conductivity, ρ the density, c p the specific heat and q e the heat source. In this first approach, these parameters do not depend on the porosity. Furthermore, the powder is subjected to mechanical loads exerted by the pistons. Stresses are generated inside the powder and tools. Here, as a first approach, tooling components are considered rigid. The mechanical problem in the powder (considered here as a continuum) is governed by the momentum equation, in which gravity and inertia effects are ignored:

0 = ⋅ ∇ σ (5)
where σ σ σ σ is the stress tensor in the powder.

The densification of the powder during electrically activated sintering results from different physical phenomena: [3] plastic deformation of grains, mass transport between grains by diffusion under load, electro-migration, or by surface tension. Therefore, the global strain rate of the compressible continuum is the sum of four contributions:

st em dl pl ε ε ε ε ε + + + = (6)
It is assumed here that plastic deformation of powder particles is the dominant phenomenon and consequently, in a first approach, only the contribution pl ε & is considered. A macroscopic model is used to follow the powder density evolution. More precisely, the powder is considered as a continuum characterized by the rate of its relative density, which is defined as the ratio of the material volume by the apparent volume (ρ r =V m /V a ). The density of the continuum is then expressed as ρ= ρ r ρ 0 , where ρ 0 the density in the dense state. The framework is the same as the one proposed by Shima & Oyane complemented by the thermo dynamical justifications of Abouaf et al.. [9][10] For a given metallic alloy, supposing that the dense material obeys a one-dimensional constitutive equation ( )

,T σ ε F = &
, where σ and ε & denote respectively the von Mises generalized stress and the generalized plastic strain rate, then the powder continuum composed of particles of the same alloy satisfies the same relation, but with the following expressions for invariants σ and ε & , proposed by Green: [11] ( )

2 / 1 2 tr ) ( ) dev( : ) dev( ) ( 2 3       + = σ σ σ r r ρ f ρ c σ (7) 2 / 1 2 ) (tr ) ( 9 2 ) ( 9 1 ) ( 3 2               - + = ε ε : ε r r r ρ c ρ f ρ c ε (8)
Parameters c and f are two decreasing functions of the relative density ρ r . In fact, c and f get their minimum value, 1 and 0 respectively, when densification is completed. During SPS compaction, high temperatures are reached and maintained for a few minutes (typically 800°C to 1200°C, in 10 to 20 min), so a viscoplastic model is considered based on a Norton-Hoff formulation. The viscoplastic potential φ is written:

( ) m m 3 K σ T K m m 1 1 +       + = ϕ ( 9
)
where K(T) is the material consistency and m the strain rate sensitivity index. According to Abouaf, [10] ε & in the porous state should be expressed as in the dense state:

σ ε ∂ ∂ = ϕ & (10)
Consequently the stress tensor can be separated into deviatoric and volumetric contributions and we write respectively the deviatoric stress tensor and the hydrostatic pressure as:

( ) ( ) ( ) v ε v e v e v s vp m r ε ρ c dev ) ( with ) ( ) ( 3 ) ( ) ( 2 1 = = - (11) ( ) ( ) ( ) ( ) ( ) ( ) v ε v σ vp r m ρ f ε T K p & & tr ) ( 3
1 3 tr 3 1 1 - - = - = (12)
In this study, and as a first approach, friction between powder and tooling is neglected, perfect slip conditions are considered.

The mechanical problem includes another equation from which the relative density is deduced. More precisely, during compaction, closing of pores causes the variation of the apparent volume and so the relative density. Since ρ=ρ r ρ 0 , the evolution of ρ r is calculated by solving mass conservation equation:

0 = ⋅ ∇ + ⋅ ∇ + ∂ ∂ v v r r r ρ ρ t ρ (13)
where v denotes the velocity field of the powder continuum. Once calculated, f and c functions are updated using the following expressions: [12] 1 and 1

+ = - - = f c ρ ρ ρ k f c c r r f α (14)
Where k f , ρ c and α c are material parameters. All physical properties are dependent on the temperature and are implemented as functions of the temperature calculated at each time step.

Numerical approach

The electric and thermal problems are solved using a monolithic method, through the CimLib library developed at CEMEF. The assembly is represented by a single mesh (Fig. 1)

and a single equation is solved for each problem on the entire geometry. As a result, the different boundary conditions between the different components of the SPS assembly, like heat conduction, are avoided but on the other hand, specific methods are used in order to characterize the subdomains occupied by the different materials. The concept behind is based on using a function to differentiate the materials and on a mesh adaptation technique providing mesh control at the interfaces.

Let Ω be the computational domain, divided into a certain number of subdomains Ω i where each subdomain refers to the domain occupied by the powder or by the different materials constituting the parts of the assembly: dies, spacers and pistons. These subdomains Ω i are described with level set functions α i (x), each one being the signed distance to the interface Γ i boundary of Ω i . [13] Therefore, the function α i (x) is expressed as follows:

     Γ ∈ Ω ∉ Γ - Ω ∈ Γ = i i i i i i x x x x x x if 0 if ) , d( if ) , d( ) ( α (15)
Once calculated, the level set allows us to define a presence function of the subset Ω i . A "smoothed" Heaviside function H(α i ) is used in this work, where a thickness ε m is fixed in the surrounding of the interface in order to prevent the discontinuous transition in the region:

       > ≤                 + + - < = m i m i m i m i m i i ε ε ε ε ε ε ε ε ε ε ε ε π π π π π π π π ε ε ε ε ε ε ε ε α if 1 α if α sin 1 α 1 2 1 α if 0 ) H(α (16)
In calculations, the value of ε m is chosen accordingly to the mesh size. Furthermore, the objective being to calculate the physical properties on the entire geometry, mixing laws are introduced to express the properties on Ω. These laws are defined as functions of the level set.

Different expressions can be used, either approximating at nodes or at elements (P1 or P0type). As an example, for the density and specific heat, linear interpolations are applied at the neighborhood of the interface between the two subdomains i and j:

j i i i ρ ρ ρ )) H(α 1 ( ) H(α - + = j p i i p i p c c c )) H(α 1 ( ) H(α - + = (17)
Besides, same approach is used for the conductivities; but it can be shown that harmonic average provides better results:

j e i i e i e σ σ σ , , )) H(α 1 ( ) H(α 1 - + = j i i i k k k )) H(α 1 ( ) H(α 1 - + = (18)
Whereas, concerning the consistency K and the strain rate sensivity coefficient m, discontinuous mixing laws show better results, following expressions are used: Ω e being the element on which the mixing is applied, ψ i and ψ j the respective properties of Ω i and Ω j .

In order to enhance the accuracy around the interface, the level set method and mixing laws are associated with an anisotropic adaptive meshing. More precisely, the main point is to keep an isotropic mesh in the area far from the interface, and to generate a refined mesh with anisotropic elements in its neighbourhood. For that reason, the metric used for the mesh is dependent on the gradient of the level set function:

α α α where otherwise 2 α if ∇ ∇ ⊗ ∇ =      +       - > = T rem rem rem ε ε e N e ε A Ι A Ι M (20)
where N is the number of elements generated in the thickness 2e and in the direction of the gradient of the level set, and ε rem is related to the background mesh size. In practice, the mesh is generated using the MTC mesher developed by Coupez. [14][15] Simulations Figure 1 presents the 3D geometry considered in the simulation of the electric-thermalmechanical coupled problems. Due to symmetry reasons, only a quarter of the set up is simulated.

An electrical current density is applied on the upper graphite edge and a constant zero potential is applied on the lower graphite edge. The current density is continuously controlled by a PID algorithm ("Proportional-Integral-Derivative"). The temperature is controlled in the mould 3 mm away from the specimen surface. At each time step, the injected current density J is calculated as a function of the error between calculated and prescribed temperature if K p , K i and K d are the proportional, integral and derivative constants J is written:

Two different material samples are chosen: TiAl and Al 2 O 3 . The mould, pistons and spacers are made of graphite. To minimize heating rates in the SPS setup, two Inconel spacers are in contact with the graphite edges and are water cooled (Fig. 1). As Inconel pieces are not included in our numerical domain, it is necessary to express the thermal boundary condition directly onto the boundaries of the top and bottom graphite spacers. This can be done as follows. Assuming a quasi-stationary state, the heat flux at the Inconel boundary is equal to the heat flux through the Inconel piece, and to the heat flux through the interface between Inconel and graphite. As a consequence, using these two expressions, the heat flux on the graphite boundary can be expressed as a function of the graphite temperature along this boundary and the water temperature:

( )

water graphite inconel inconel cond inconel inconel cond T T e k h e k h - + = φ ( 23
)
here h cond = 880 W.m -2 .K -1 is a constant heat transfer coefficient, k inconel the Inconel thermal conductivity, e inconel its thickness and T water the water temperature (assumed constant 296 K).

Moreover, the process taking place in a vacuum chamber, heat losses by convection and conduction are neglected. All lateral surfaces have heat losses by radiation towards chamber walls, which are held at room temperature (300 K), a radiative heat flux is considered:

( 4 4 wall r B radiation T T ε σ - = φ φ φ φ ( 24
)
where σ B = 5.67x10 -8 W.m -2 .K -4 is the Stefan-Boltzmann's constant, r the emissivity (assumed to be equal to 0.8 in this study), T the local temperature of the vertical boundary surfaces and T wall the temperature of the chamber walls. On the other hand, a constant voltage is applied during the simulations. Furthermore, all contact surfaces between the parts of the SPS assembly are considered ideal: no contact resistance is taken into account. Physical parameters are taken from Molenat et al.. [8] The powder follows the viscoplastic law for porous materials as mentioned previously.

Besides, tools are assumed to be rigid (stresses are not calculated inside graphite) and a frictionless contact is considered between tools and powder. Moreover, a constant velocity of 6x10 -3 mm/s is applied to the upper tooling, whereas the mould and lower tooling are considered fixed. The initial porosity in the TiAl powder is considered to be equal to 30%.

The consistency K(T) and the strain rate sensitivity m(T) are thermal dependent, their values are deduced from literature: measurements carried out on Ti-47Al-2Cr-2Nb. [16] Results

Electrical potential and temperature distribution

The distribution of the temperature during the SPS process is directly related to the distribution of current which is firstly affected by the geometry and the physical properties of the different components and materials involved in the whole set-up. The choice of the two samples, with two different electrical conductivities, shows effect on the distribution of the current and consequently on the distribution of the temperature. In Fig. 2, the norm of the current density is compared. Results show that the highest current density is in the region of the pistons in both cases, due to their smaller section. Besides, in the case of TiAl, most of the current flows through the sample, the distribution of the current density is more important at the lateral surface of the conducting specimen and decreases in the radial direction towards the centre during the process. Whereas in the case of alumina the current is, as expected, relatively absent inside the sample, because of its low electrical conductivity at low temperature. In fact, at the beginning, alumina is heated only by conduction and not directly by Joule effect, its electrical conductivity increases with temperature, allowing the current to go through after 15 min, as it can be seen in figure 2.

Since the homogeneity of the microstructure in the SPS parts is directly dependent on the temperature homogeneity, it is interesting to understand the evolution of the temperature distribution during the process. This was studied in several works. [5][6] Pistons are heated first due to the concentration of electrical current in their reduced section and present the highest temperatures along the process, that way increasing heating of the powder by conduction.

Besides, the gradient of temperature inside the sample is directly related to the electrical conductivity since the current is the only heat source. [6] Actually, the electrical conductivity of the powder affects the distribution of the current in the graphite that surrounds it (Fig. 2), consequently the distribution of the current in graphite affects the distribution of the source term and so the temperature gradient inside the sample. As shown in the graphics of figure 3, the temperature gradient inside the alumina specimen varies with time. The increase of its electrical conductivity with temperature changes the current distribution and so the heating phenomenon.

Moreover, geometry of the setup could be a major cause of heterogeneities. When the load is applied, the symmetry in the vertical direction might be lost. Coupling with the mechanical problem allows controlling this source of homogeneity. This can be seen in figure 4, where a displacement is imposed in the upper tools to compact the powder while the lower tools are kept fixed.

In figure 5, the axial temperature profiles after 25 s and 900 s show that once symmetry is lost, heating is no longer similar on both sides. When the piston enters into the mould, the height of the restricted section is lower. Consequently, Joule effect is less important since the concentration of the current flux covers a smaller region.

The temperature heterogeneities in the sample affect directly the mechanical problem. Since the constitutive law is temperature dependant, the deformation resulting from the pressure applied during the SPS process will be heterogeneous. Densification won't be achieved identically inside the specimen. In order to better display the evolution of porosity, a different geometry is chosen to simulate the mechanical problem (figure 6). Same boundary conditions are used as previously. The fully coupled algorithm (electrical, thermal and mechanical) is used. Then all level set functions are updated as well as remeshing to conserve interfaces. The norm of the current density is presented in Fig. 6, the current is concentrated in the bottom of the powder where the radius is smaller. However, the source term of the heat equation is not affected significantly and so the gap of temperature in the powder doesn't exceed 40 °C.

In that case, pressure gradients result mainly from the geometry. The pressure is the highest in the upper section of the specimen and decreases towards the centre and the bottom. Since the consistency decreases with temperature, the pressure decreases with time as seen in figure 7.

The heterogeneity of pressure induces the heterogeneity during densification; the porosity is lower when pressure is high. The initial relative density is set to 0.7, after 16 min it reaches 0.99 in the upper region where pressure is the highest, as seen in figure 8.

Hence, simulations allow predicting densification during the SPS process in order to optimize the compacted specimen. The position of the specimen in the mould and its orientation is one of the major issues that affect homogeneity.

Conclusions

A 3D finite element analysis has been presented to simulate the electrical thermal and mechanical problems during the SPS process. It is based on CimLib using a monolithic approach. Its application to two samples (TiAl and alumina) has proven that distribution of current strongly depends on geometry and physical properties. It was shown how, in the case of alumina, the increase of its electrical conductivity causes a variation of the temperature distribution inside the powder. Besides by coupling with the mechanical problem, the distribution of the temperature in the axial direction was shown to be no longer symmetrical due to the motion of the upper tools.

In a second step, mechanical results were shown for a more complex geometry in which the pressure values are heterogeneous, causing quite different local densification kinetics. The preceding results show that the coupled electrical-thermal-mechanical finite element simulation may be of great help in understanding the temperature distribution and the densification during the SPS process. Such simulations should allow engineers to better determine the nature and geometry of the involved tooling materials, to better control microstructure formation. Forthcoming developments should aim at enriching the finite

  field and e σ the electric conductivity. The solution of the electrical problem consists in solving the following Poisson type equation for the electrical potential U :

Fig. 1 .Fig. 2 .Fig. 3 .

 123 Fig. 1. Geometry of the SPS set-up with the generated mesh used in the simulation and a zoom in the region of the sample showing refined elements around the sample/matrix interface (radius of the sample 18 mm, of the matrix 35 mm ).

Fig. 4 .

 4 Fig. 4. Vertical motion of the upper tools causing non symmetrical distribution of the temperature in SPS machine. The extreme right picture better shows the effect by filtering temperature values in the central region.

Fig. 5 .Fig. 6 .Fig. 7 .Fig. 8 .

 5678 Fig. 5. Temperature variation along the vertical axis for the TiAl sample for different time steps showing the effect of the symmetry.

& & & & &

& & & &

& & & &

element model through the elimination of the different simplifying assumptions mentioned along the paper.

This study has been conducted in the framework of the cooperative project "IRIS" supported by the French Agence Nationale de la Recherche (ANR).