Hachem Kadri
email: hachem.kadri@inria.fr

Philippe Preux
email: philippe.preux@inria.fr

Alain Rakotomamonjy
email: alain.rakotomamonjy@insa-rouen.fr

Functional Regularized Least Squares Classification with Operator-valued Kernels

Although operator-valued kernels have recently received increasing interest in various machine learning and functional data analysis problems such as multi-task learning or functional regression, little attention has been paid to the understanding of their associated feature spaces. In this paper, we explore the potential of adopting an operatorvalued kernel feature space perspective for the analysis of functional data. We then extend the Regularized Least Squares Classification (RLSC) algorithm to cover situations where there are multiple functions per observation. Experiments on a sound recognition problem show that the proposed method outperforms the classical RLSC algorithm.

Introduction

Following the development of multi-task and complex output learning methods [START_REF] Pontil | Open House on Multi-Task and Complex Outputs Learning[END_REF], operator-valued kernels have recently attracted considerable attention in the machine learning community (Micchelli & Pontil, 2005b;[START_REF] Reisert | Learning equivariant functions with matrix valued kernels[END_REF][START_REF] Caponnetto | Universal multi-task kernels[END_REF]. It turns out that these kernels lead to a new class of algorithms well suited for learning multi-output functions. For example, in multi-task learning contexts, standard single-task kernel learning methods such as support vector machines (SVM) and regularization networks (RN) are extended to deal with several dependent tasks at once by learning a vector-valued function using multi-task kernels [START_REF] Evgeniou | Learning multiple tasks with kernel methods[END_REF]. Also, in functional data analysis (FDA) where observed continuous data are measured over a densely sampled grid and then represented by real-valued functions rather than by discrete finite dimensional vectors, function-valued reproducing kernel Hilbert spaces (RKHS) are constructed from nonnegative operator-valued kernels to extend kernel ridge regression from finite dimensions to the infinite dimensional case [START_REF] Lian | Nonlinear functional models for functional responses in reproducing kernel Hilbert spaces[END_REF][START_REF] Kadri | Nonlinear functional regression: a functional RKHS approach[END_REF].

While most recent work has focused on studying operator-valued kernels and their corresponding RKHS from the perspective of extending Aronszajn's pioneering work (1950) to the vector or functionvalued case (Micchelli & Pontil, 2005a;[START_REF] Carmeli | Vector valued reproducing kernel Hilbert spaces and universality[END_REF][START_REF] Kadri | Nonlinear functional regression: a functional RKHS approach[END_REF], in this paper we pay special attention to the feature space point of view [START_REF] Schölkopf | Input space vs. feature space in kernel-based methods[END_REF]. More precisely, we provide some ideas targeted at advancing the understanding of feature spaces associated with operator-valued kernels and we show how these kernels can design more suitable feature maps than those associated with scalar-valued kernels, especially when input data are complex, infinite dimensional objects like curves and distributions.

In many experiments, the observations consist of a sample of random functions or curves. So, we adopt in this paper a functional data analysis point of view in which each curve corresponds to one observation. This is an extension of multivariate data analysis where observations consist of vectors of finite dimension. For an introduction to the field of FDA, the two monographs by [START_REF] Ramsay | Functional Data Analysis[END_REF]2002) provide a rewarding and accessible overview on foundations and applications, as well as a collection of motivating examples [START_REF] Müller | Functional modelling and classification of longitudinal data[END_REF].

To explore the potential of adopting an operatorvalued kernel feature space approach, we are interested in the problem of functional classification in the case where there are multiple functions per observation. This study is valuable from a variety of perspectives. Our motivating example is the practical problem of sound recognition which is of great importance in surveillance and security applications [START_REF] Dufaux | Automatic sound detection and recognition for noisy environment[END_REF][START_REF] Istrate | Information extraction from sound for medical telemonitoring[END_REF][START_REF] Rabaoui | Using one-class svms and wavelets for audio surveillance[END_REF]. This problem can be tackled by classifying incoming signals representing the environmental sounds into various predefined classes. In this setting, a preprocessing step consists in applying signal-processing techniques to generate a set of features characterizing the signal to be classified. These features form a so-called feature vector which contains discrete values of different functional parameters providing information about temporal, frequential, cepstral and energy characteristics of the signal. In standard machine learning methods, the feature vector is considered to be a subset of R n by concatenating samples of the different functional features, and this has the drawback of not considering any dependencies between different values over subsequent time-points within the same functional datum. Employing these methods implies that permuting time points arbitrarily, which is equivalent to exchanging the order of the indexes in a multivariate vector, should not change the result of statistical analysis. Taking into account the inherent sequential nature of the data and using the dependencies along the time-axis should lead to higher quality results [START_REF] Lee | Functional data analysis: classification and regression[END_REF]. In our work, we use a functional data analysis approach based on modeling each sound signal by a vector of functions (in (L 2 ) p for example, where L 2 is the space of square integrable functions and p is the number of functional parameters) rather than by a vector in R n , with the hope of improving performance by: (1) considering the relationship between samples of a function variable and thus the dynamic behavior of the functional data, (2) capturing discriminative characteristics between functions contrary to the concatenation procedure.

During the last decade, various kernel-based methods have become very popular for solving classification problems. Among them, the regularized least squares classification (RLSC) is a simple regularization algorithm which achieves good performance, nearly equivalent to that of the well-known Support Vector Machines (SVMs) [START_REF] Rifkin | Regularized least squares classification[END_REF][START_REF] Lian | Nonlinear functional models for functional responses in reproducing kernel Hilbert spaces[END_REF][START_REF] Rifkin | In defense of one-vsall classification[END_REF][START_REF] Zhang | Svm vs regularized least squares classification[END_REF]. By using operator-valued kernels, we extend the RLSC algorithm to cover situations where there are multiple functions per observation. One main obstacle for this extension involves the inversion of the block operator kernel matrix (kernel matrix where the block entries are linear operators). In contrast to the situation in the multivariate case, this inversion is not always feasible in infinite dimensional Hilbert spaces. In this paper, we attempt to overcome this problem by characterizing a class of operator-valued kernels and performing an eigenvalue decomposition of this kernel matrix.

The remainder of this paper is organized as follows.

In section 2, we review concepts of operator-valued kernels and their corresponding reproducing kernel Hilbert spaces and discuss some ideas for understanding the associated feature maps. Using these kernels, we propose a functional regularized least squares classification algorithm in section 3. It is an extension of the classical RLSC to the case where there are multiple functions per observation. The proposed algorithm is experimentally evaluated on a sound recognition task in section 4. Finally, section 5 presents some conclusions and future work directions.

Operator-valued kernels and associated feature spaces

In the machine learning literature, operator-valued kernels were first introduced by Micchelli and Pontil (2005b) to deal with the problem of multi-task learning. In this context, these kernels, called multitask kernels, are matrix-valued functions and their corresponding vector-valued reproducing kernel Hilbert spaces are used to learn multiple tasks simultaneously [START_REF] Evgeniou | Learning multiple tasks with kernel methods[END_REF]. Most often, multitask kernels are constructed from scalar-valued kernels which are carried over to the vector-valued setting by a positive definite matrix. More details and some examples of multi-task kernels can be found in [START_REF] Caponnetto | Universal multi-task kernels[END_REF]. In addition, some works have studied the construction of such kernels from a functional data analysis point of view (continuous data). For example, in [START_REF] Kadri | Nonlinear functional regression: a functional RKHS approach[END_REF], the authors showed how infinite dimensional operator-valued kernels can be used to perform nonlinear functional regression in the case where covariates as well as responses are functions.

Since we are interested in the problem of functional classification, we use a similar framework to that of [START_REF] Kadri | Nonlinear functional regression: a functional RKHS approach[END_REF], but we consider the more general case where we have multiple functions as inputs to the learning module. Kernel-based learning methodology can be extended directly from the vector-valued to the functional-valued case. The principle of this extension is to replace vectors by functions and matrices by linear operators; scalar products in vector space are replaced by scalar products in function space, which is usually chosen as the space of square integrable functions L 2 on a suitable domain. In the present paper, we focus on infinite dimensional operator-valued kernels and their corresponding functional-valued RKHS.

An operator-valued kernel K : X × X -→ L(Y ) is the reproducing kernel of a Hilbert space of functions from an input space X which takes values in a Hilbert space Y . L(Y ) is the set of all bounded linear operators from Y into itself. Input data are represented by a vector of functions, so we consider the case where X ⊂ (L 2 ) p and Y ⊂ L 2 . Function-valued RKHS theory is based on the one-to-one correspondence between reproducing kernel Hilbert spaces of function-valued functions and positive operator-valued kernels. We start by recalling some basic properties of such spaces. We say that a Hilbert space F of functions f : X -→ Y has the reproducing property, if ∀x ∈ X the linear functional f -→ f (x), y Y is continuous for any x ∈ X and y ∈ Y . By the Riesz representation theorem it follows that for a given x ∈ X and for any choice of y ∈ Y , there exists an element h y x ∈ F, s.t.

∀f ∈ F h y x , f F = f (x), y Y
We can therefore define the corresponding operator-

valued kernel K : X × X -→ L(Y ) such that K(x 1 , x 2 )y 1 , y 2 Y = h y1 x1 , h y2 x2 F It follows that h y1 x1 (x 2 ), y 2 Y = h y1 x1 , h y2 x2 F = K(x 1 , x 2 )y 1 , y 2 Gy
and thus we obtain the reproducing property

K(x, .)y, f F = f (x), y Y (1)
Consequently, we obtain that K(., .) is a positive definite operator-valued kernel as defined below: (see proposition 1 in (Micchelli & Pontil, 2005b) for the proof)

Definition: We say that K(x 1 , x 2 ), satisfying K(x 1 , x 2 ) = K(x 2 , x 1 ) * (the superscript * indicates the adjoint operator), is a positive definite operatorvalued kernel if given an arbitrary finite set of points {(x i , y i )} i=1,...,n ∈ X ×Y , the corresponding block matrix K with

K ij = K(x i , x j )y i , y j Y is positive semi- definite.
Importantly, the converse is also true. Any positive operator-valued kernel K(x 1 , x 2 ) gives rise to an RKHS F K , which can be constructed by considering the space of function-valued functions f having the form f (.) = n i=1 K(x i , .)y i and taking completion with respect to the inner product given by

K(x 1 , .)y 1 , K(x 2 , .)y 2 F = K(x 1 , x 2 )y 1 , y 2 Y .
In the following, we present an example of functionvalued RKHS with functional inputs and the associated operator-valued kernel.

Example. Let X = H and Y = L 2 (Ω), where H is the Hilbert space of constants in [0, 1] and L 2 (Ω) the space of square integrable functions on Ω. We denote by M the space of L 2 (Ω)-valued functions on H whose

norm g 2 M = Ω H [g(v)(x)] 2 dvdx is finite.
Let (F; ., . F ) be the space of functions from H to L 2 (Ω) such that:

   F = {f, ∃f = df (v) dv ∈ M, f (u) = u 0 f (v)dv} f 1 , f 2 F = f 1 , f 2 M F is a RKHS with kernel K(u, v) = M ϕ(u,v)
. M ϕ is the multiplication operator associated with the function ϕ where ϕ(u, v) is equal to u if u(x) ≤ v(x) ∀x ∈ Ω and v otherwise. It is easy to check that K is Hermitian and nonnegative. Now we show that the reproducing property holds for any

f ∈ F, w ∈ L 2 (Ω) and u ∈ H f, K(u, .)w F = f , [K(u, .)w] M = Ω H [f (v)](x)[K(u, v)w] (x)dvdx = Ω u 0 [f (v)](x)w(x)dvdx = Ω [f (u)](x)w(x)dx = f (u), w L 2 (Ω)
Similar to the scalar case, operator-valued kernels provide an elegant way of dealing with nonlinear algorithms by reducing them to linear ones in some feature space F nonlinearly related to input space. A feature map associated with an operator-valued kernel K is a continuous function

Φ : X × Y -→ L(X, Y ) such that, for every x 1 , x 2 ∈ X and y 1 , y 2 ∈ Y K(x 1 , x 2 )y 1 , y 2 Y = Φ(x 1 , y 1 ), Φ(x 2 , y 2 ) L(X,Y )
where L(X, Y ) is the set of mappings from X to Y . By virtue of this property, Φ is called a feature map associated with K. Furthermore, from (1), it follows that in particular

K(x 1 , .)y 1 , K(x 2 , .)y 2 F = K(x 1 , x 2 )y 1 , y 2 Y
which means that any operator-valued kernel admits a feature map representation with a feature space F ⊂ L(X, Y ), and corresponds to a dot product in another space.

From this feature map perspective, we study the geometry of a feature space associated with an operatorvalued kernel and we compare it with the one obtained by a scalar-valued kernel. More precisely, we consider two reproducing kernel Hilbert spaces (RKHS) F and H. F is a RKHS of function-valued functions on X with values in Y . X ⊂ (L 2 ) p , Y ⊂ L 2 and let K be the reproducing operator-valued kernel of F. H is also a RKHS, but of scalar-valued functions on X with values in R, and k its reproducing real-valued kernel. The mappings Φ K and Φ k associated, respectively, with the kernels K and k are defined as follows

Φ y K : (L 2 ) p → L((L 2 ) p , L 2 ), x → K(x, .)y and Φ k : (L 2 ) p → L((L 2 ) p , R), x → k(x, .)
These feature maps can be seen as a mapping of the input data x i , which are vector of functions in (L 2 ) p , into a feature space in which the dot product can be computed using the kernel functions. This idea leads to design nonlinear methods based on linear ones in the feature space. In a supervised classification problem for example, since kernels could map input data into a higher dimensional space, kernel methods deal with this problem by finding a linear separation in the feature space between data which can not be separated linearly in the input space. We now compare the dimension of feature spaces obtained by the maps Φ K and Φ k . To do this, we adopt a functional data analysis point of view where observations are composed of sets of functions. Direct understanding of this FDA viewpoint comes from the consideration of the "atom" of a statistical analysis. In a basic course in statistics, atoms are "numbers", while in multivariate data analysis the atoms are vectors and methods for understanding populations of vectors are the focus. FDA can be viewed as the generalization of this, where the atoms are more complicated objects, such as curves, images or shapes represented by functions [START_REF] Zhao | The functional data analysis view of longitudinal data[END_REF]. Based on this, the dimension of the input space is p since x i ∈ (L 2 ) p is a vector of p functions. The feature space obtained by the map Φ k is a space of functions, so its dimension from a FDA point of view is one. The map Φ K projects the input data into a space of operators L(X, Y ). This means that using the operator-valued kernel K corresponds to mapping the functional data x i into a higher, possibly infinite, dimensional space (L 2 ) d with d → ∞. In a binary functional classification problem, we have higher probability to achieve linear separation between the classes by projecting the functional data into a higher dimensional feature space rather than into a lower one, that is why we think that it is more suitable to use operatorvalued than scalar-valued kernels in this context.

Functional regularized least squares classification

In this section, we show how to extend the regularized least squares classification algorithm (RLSC) to functional contexts using operator-valued kernels. To use these kernels for a classification problem, we consider the labels to be functions in some function space rather than real values as usual. The functional classification problem can then be framed as that of learning a function-valued function f : X -→ Y where X ⊂ (L 2 ) p and Y ⊂ L 2 . The RLSC algorithm [START_REF] Rifkin | Regularized least squares classification[END_REF] is based on solving a Tikhonov minimization problem associated with a square loss function, and then an estimate f * of f in a Hilbert space F with reproducing operator-valued kernel K : X ×X -→ L(Y ) is obtained by minimizing

f * = arg min f ∈F n i=1 y i -f (x i ) 2 Y + λ f 2 F (2)
By the representer theorem (Micchelli & Pontil, 2005a;[START_REF] Kadri | Nonlinear functional regression: a functional RKHS approach[END_REF], the solution of this problem has the following form

f * (x) = n j=1 K(x, x j )β j , β j ∈ Y (3)
Substituting (3) in (2) , we come up with the following minimization over the scalar-valued functions β i rather than the function-valued function f

min βv∈(Y ) n n i=1 β i - n j=1 K(x i , x j )β j 2 Y +λ n i,j K(x i , x j )β i , β j Y (4)
β v is the vector of functions (β i ) i=1,...,n ∈ (L 2 ) n . The problem (4) can be solved in three ways. Assuming that the observations are made on a regular grid {t 1 , . . . , t m }, one can first discretize the functions x i and y i and then solve the problem using multivariate data analysis techniques. However, this has the drawback, as well known in the FDA literature, of not considering the relationships that exist between samples. The second way consists in considering the output space Y to be a scalar valued reproducing Hilbert space. In this case, the functions β i can be approximated by a linear combination of a scalar kernel

β i = m l=1 α il k(s l , .
) and then the problem (4) becomes a minimization problem over the real values α il . Another possible way to solve the minimization ( 4) is to compute its derivative using the directional derivative and setting the result to zero to find an analytic solution of the problem. It follows that β v satisfies the system of linear operator equations

(K + λI)β v = y v ( 5
)
where

K = [K(x i , x j )] n i,j=1 is a n × n block operator matrix (K ij ∈ L(Y )
) and y v the vector of functions (y i ) n i=1 ∈ (L 2 ) n . In this work, we are interested in this third approach which extends the classical RLSC algorithm to functional data analysis domain. One main obstacle for this extension is the inversion of the block operator kernel matrix K. Block operator matrices generalize block matrices to the case where the block entries are linear operators between infinite dimensional Hilbert spaces. In contrast to the situation in the multivariate case, inverting such matrices is not always feasible in infinite dimensional spaces. To overcome this problem, we study the eigenvalue decomposition of a class of block operator kernel matrices obtained from operator-valued kernels having the following form

K(x i , x j ) = G(x i , x j )T, ∀x i , x j ∈ X (6)
where G is a scalar-valued kernel and T is an operator in L(Y ). This kernel construction is adapted from (Micchelli & Pontil, 2005a;b). Choosing T depends on the context. For multi-task kernels, T is a finite dimensional matrix which model relations between tasks. In FDA, [START_REF] Lian | Nonlinear functional models for functional responses in reproducing kernel Hilbert spaces[END_REF] suggested the use of the identity operator, while [START_REF] Kadri | Nonlinear functional regression: a functional RKHS approach[END_REF] showed that it will be more useful to choose other operators than identity that are able to take into account functional properties of the input and output spaces. They introduced a functional extension of the Gaussian kernel based on the multiplication operator. In this work, we are interested in kernels constructed from the integral operator. This seems to be a reasonable choice since functional linear model (see Eq. ( 7)) are based on this operator [START_REF] Ramsay | Functional Data Analysis[END_REF]

y(s) = α(s) + x(t)ν(s, t)dt (7) Algorithm 1 Functional RLSC Input data x i ∈ (L 2 ([0, 1])) p , size n labels y i ∈ L 2 ([0, 1]), size n Eigendecomposition of G G = G(x i , x j ) n i,j=1 ∈ R n×n eigenvalues α i ∈ R, size n eigenvectors v i ∈ R n , size n Eigendecomposition of T T ∈ L(Y ) Initialize k: number of eigenfunctions eigenvalues δ i ∈ R, size k eigenfunctions w i ∈ L 2 ([0, 1]), size k Eigendecomposition of K = G ⊗ T K = K(x i , x j ) n i,j=1 ∈ (L(Y )) n×n eigenvalues θ i ∈ R, size n × k θ = α ⊗ δ eigenfunctions z i ∈ (L 2 ([0, 1])) n , size n × k z = v ⊗ w Solution of (4) β = (K + λI) -1 y Initialize λ: regularization parameter β = n×k i=1 (θ i + λ) -1 n j=1 z i,j , y j z i
where α and ν are the functional parameters of the model. So we consider the following positive definite operator-valued kernel

(K(x i , x j )y)(t) = G(x i , x j ) Ω e -|t-s| y(s)ds (8)
where y ∈ Y and {s, t} ∈ Ω = [0, 1]. Note that a similar kernel was proposed in [START_REF] Caponnetto | Universal multi-task kernels[END_REF] for linear spaces of functions from R to Y . The n×n block operator kernel matrix K of operator-kernels having the form (6) can be expressed by a Kronecker product between the matrix G = G(x i , x j ) n i,j=1 in R n×n and the operator T ∈ L(Y )

K =    G(x 1 , x 1 )T . . . G(x 1 , x n )T . . . . . . . . . G(x n , x 1 )T . . . G(x n , x n )T    = G ⊗ T
In this case, the eigendecomposition of the matrix K can be obtained from the eigendecompositions of G and T (see Algorithm 1). Let θ i and z i be, respectively, the eigenvalues and the eigenfunctions of K, the inverse operator K -1 is given by

K -1 y v = i θ -1 i y v , z i z i , ∀y v ∈ (L 2 ) n
Now we are able to solve the system of linear operator equation ( 5) and the functions β i can be computed from eigenvalues and eigenfunctions of the matrix K, as described in Algorithm 1.

Experiments

Our experiments are based on a sound recognition task. The performance of the functional RLSC algorithm, described in section 3, is evaluated on a data set of sounds collected from commercial databases which include sounds ranging from screams to explosions, such as gun shots or glass breaking, and compared with the RLSC method [START_REF] Rifkin | Regularized least squares classification[END_REF].

Many previous works in the context of sound recognition have concentrated on classifying environmental sounds other than speech and music [START_REF] Dufaux | Automatic sound detection and recognition for noisy environment[END_REF][START_REF] Peltonen | Computational audiroty scene recognition[END_REF]. Such sounds are extremely versatile, including signals generated in domestic, business, and outdoor environments. A system that is able to recognize such sounds may be of great importance for surveillance and security applications [START_REF] Istrate | Information extraction from sound for medical telemonitoring[END_REF][START_REF] Rabaoui | Using one-class svms and wavelets for audio surveillance[END_REF]. The classification of a sound is usually performed in two steps. First, a pre-processor applies signal processing techniques to generate a set of features characterizing the signal to be classified. Then, in the feature space, a decision rule is implemented to assign a class to a pattern.

Database description

As in [START_REF] Rabaoui | Using one-class svms and wavelets for audio surveillance[END_REF], the major part of the sound samples used in the recognition experiments is taken from two sound libraries (Leonardo Software;Real World Computing Paternship, 2000). All signals in the database have a 16 bits resolution and are sampled at 44100 Hz, enabling both good time resolution and a wide frequency band, which are both necessary to cover harmonic as well as impulsive sounds. The selected sound classes are given in Table 1, and they are typical of surveillance applications. The number of items in each class is deliberately not equal.

Note that this database includes impulsive sounds and harmonic sounds such as phone rings (C6) and children voices (C7). These sounds are quite likely to be recorded by a surveillance system. Some classes sound very similar to a human listener: in particular, explosions (C4) are pretty similar to gunshots (C2). Glass breaking sounds include both bottle and window breaking situations. Phone rings are either electronic or mechanic alarms. Temporal representations and spectrograms of some sounds are depicted in Figure 1 and 2. Power spectra are extracted through the Fast Fourier Transform (FFT) every 10 ms from 25 ms frames. They are represented vertically at the corresponding frame indexes. The frequency range of interest is between 0 and 22 kHz. A lighter shade indicates a higher power value. These figures show that in the considered database we can have both: (1) many similarities between some sounds belonging to different classes, (2) diversities within the same sound class.

Results

Following [START_REF] Rifkin | In defense of one-vsall classification[END_REF], the 1-vs-all multiclass classifier is selected in these experiments. So we train N (number of classes) different binary classifiers, each one trained to distinguish the data in a single class from the examples in all remaining classes. We run the N classifiers to classify a new example. In section 3, we showed that operator-valued kernels can be used in a classification problem by considering the labels y i to be functions in some function space rather than real values. Similarly to the scalar case, a natural choice for y i would seem to be the Heaviside step function in L 2 ([0, 1]) scaled by a real number. The used operator valued-kernel is based on the integral operator as defined in (8). Eigenvalues δ i and eigenfunctions w i associated with this kernel are equal to 2 1+µ 2 i and µ i cos(µ i x) + sin(µ i x), respectively ; where µ i are solutions of the equation cot µ = 1 2 (µ -1 µ ). The adopted sound data processing scheme is the following. Let X be the set of training sounds, shared in N classes denoted C 1 , . . . , C N . Each class contains m i sounds, i = 1, . . . , N . Sound number j in class C i is denoted s i,j , (i = 1, . . . , N, j = 1, . . . , m i ). The preprocessor converts a recorded acoustic signal s i,j into a time/frequency localized representation. In multivariate methods, this representation is obtained by splitting the signal s i,j into T i,j overlapping short frames and computing a vector of features z t,i,j , t = 1, . . . , T i,j which characterize each frame. Since the pre-processor is a series of continuous time-localized features, it will be useful to take into account the relationships between feature samples along the time axis and consider dependencies between features. That is why we use a FDA-based approach in which features representing a sound are modeled by functions z i,j (t).

In this work, Mel Frequency Cepstral Coefficients (MFCCs) are used to describe the spectral shape of each signal. These coefficients are obtained using 23 channels Mel filterbank and a Hamming analysis window of length 25 ms with 50% overlap. We choose to use 13 MFCC features and the energy parameter measured along the sound signal. So, each sound is characterized by 14 functional parameters: 13 cepstral functions and one energy function.

Performance of the Functional RLSC based classifier is compared to the results obtained by the RLSC algorithm, see Table 2 and3. The performance is measured as the percentage number of sounds correctly recognized and it is given by (W r /T n ) × 100%, where W r is the number of well recognized sounds and T n is the total number of sounds to be recognized. The use of the Functional RLSC is fully justified by the results presented here, as it yields consistently lower error rates and a high classification accuracy for the major part of the sound classes.

Conclusion

This paper has put forward the idea that by viewing operator-valued kernels from a feature map perspective, we can design more general kernel methods wellsuited for complex data. Based on this, we have extended the regularized least squares classification algorithm to functional data analysis contexts where input data are real-valued functions rather than finite dimensional vectors. Through experiments on sound recognition, we have shown that the proposed approach is efficient and improves the classical RLSC method in a sound classification dataset. Further investigations will consider larger classes of operator-valued kernels. It would also be interesting to study learning methods for choosing the operator-valued kernel.

Figure 1 .

 1 Figure 1. Structural similarities between two different classes.

Figure 2 .

 2 Figure 2. Structural diversity inside the same sound class and between classes.

Table 1 .

 1 Classes of sounds and number of samples in the database used for performance evaluation.

	Classes	Number Train Test	Total Duration(s)
	Human screams	C1	40	25	65	167
	Gunshots	C2	36	19	55	97
	Glass breaking	C3	48	25	73	123
	Explosions	C4	41	21	62	180
	Door slams	C5	50	25	75	96
	Phone rings	C6	34	17	51	107
	Children voices	C7	58	29	87	140
	Machines	C8	40	20	60	184
	Total		327	181	508	18mn 14s

Table 2 .

 2 Confusion Matrix obtained when using the Functional Regularized Least Squares classification algorithm

Table 3 .

 3 Confusion Matrix obtained when using the Regularized Least Squares Classification (RLSC) algorithm

		C1	C2	C3	C4	C5	C6	C7	C8
	C1	92	4	4.76	0	5.27	11.3	6.89	0
	C2	0	52	0	14	0	2.7	0	0
	C3	0	20	76.2	0	0	0	17.24	5
	C4	0	16	0	66	0	0	0	0
	C5	4	8	0	4	84.21	0	6.8	0
	C6	4	0	0	0	10.52	86	0	0
	C7	0	0	0	8	0	0	69.07	0
	C8	0	0	19.04	8	0	0	0	95
			Total Recognition Rate = 77.56%

Acknowledgments

This work was supported by Ministry of Higher Education and Research, Nord-Pas-de-Calais Regional Council and FEDER through the 'Contrat de Projets Etat Region (CPER) 2007-2013'. H.K. is supported by Junior Researcher Contract No. 4297 from the Nord-Pas-de-Calais region. A.R. is supported by the PASCAL2 Network of Excellence, ICT-216886, ANR Project ASAP ANR-09-EMER-001 and the IN-RIA ARC MABI.


Hachem Kadri

Philippe Preux

Emmanuel Duflos
email: emmanuel.duflos@inria.fr

Stéphane Canu
email: scanu@insa-rouen.fr

Multiple functional regression with both discrete and continuous covariates

In this paper we present a nonparametric method for extending functional regression methodology to the situation where more than one functional covariate is used to predict a functional response. Borrowing the idea from Kadri et al. (2010a), the method, which support mixed discrete and continuous explanatory variables, is based on estimating a function-valued function in reproducing kernel Hilbert spaces by virtue of positive operator-valued kernels.

Introduction

The analysis of interaction effects between continuous variables in multiple regression has received a significant amount of attention from the research community. In recent years, a large part of research has been focused on functional regression where continuous data are represented by real-valued functions rather than by discrete, finite dimensional vectors. This is often the case in functional data analysis (FDA) when observed data have been measured over a densely sampled grid. We refer the reader to Ramsay andSilverman (2002, 2005) and [START_REF] Ferraty | Nonparametric Functional Data Analysis[END_REF] for more details on functional data analysis of densely sampled data or fully observed trajectories. In this context, various functional regression models [START_REF] Ramsay | Functional data analysis[END_REF] have been proposed according to the nature of explanatory (or covariate) and response variables, perhaps the most widely studied is the generalized functional linear model where covariates are functions and responses are scalars [START_REF] Cardot | Functional linear model[END_REF][START_REF] Cardot | Spline Estimators for the Functional Linear Model[END_REF][START_REF] James | Generalized linear models with functional predictors[END_REF][START_REF] Müller | Generalized functional linear models[END_REF][START_REF] Preda | Regression models for functional data by reproducing kernel Hilbert spaces methods[END_REF].

In this paper, we are interested in the case of regression models with a functional response. Two subcategories of such models have appeared in the FDA literature: covariates are scalars and responses are functions also known as "functional response model" [START_REF] Faraway | Regression analysis for a functional response[END_REF][START_REF] Chiou | Functional response models[END_REF]; both covariates and responses are functions [START_REF] Ramsay | Some tools for functional data analysis[END_REF][START_REF] He | Extending correlation and regression from multivariate to functional data. Asymptotics in statistics and probability[END_REF][START_REF] Cuevas | Linear functional regression: the case of fixed design and functional response[END_REF][START_REF] Prchal | Spline estimator for the functional linear regression with functional response[END_REF][START_REF] Antoch | Functional linear regression with functional response: application to prediction of electricity consumption[END_REF]. In this work, we pay particular attention to this latter situation which corresponds to extending multivariate linear regression model to the functional case where all the components involved in the model are functions. Unlike most of previous works which consider only one functional covariate variable, we wish to perform a regression analysis in which multiple functional covariates are used to predict a functional response. The methodology which is concerned with solving such task is referred to as a multiple functional regression.

Previous studies on multiple functional regression [START_REF] Han | Novel Perspectives On Market Valuation of Firms Via Functional Regression[END_REF][START_REF] Matsui | Regularized functional regression modeling for functional response and predictors[END_REF][START_REF] Valderrama | Forecasting pollen concentration by a two-Step functional model[END_REF] assume a linear relationship between functional covariates and responses and model this relationship via a multiple functional linear regression model which generalizes the model in [START_REF] Ramsay | Some tools for functional data analysis[END_REF] to deal with more than one covariate variable. However, extensions to nonparametric models have not been considered. Nonparametric functional regression (Ferraty andVieu, 2002,2003) is addressed mostly in the context of functional covariates and scalar responses. More recently, [START_REF] Lian | Nonlinear functional models for functional responses in reproducing kernel hilbert spaces[END_REF] and Kadri et al. (2010a) showed how function-valued reproducing kernel Hilbert spaces (RKHS) and operator-valued kernels can be used for the nonparametric estimation of the regression function when both covariates and responses are curves. Building on these works, we present in this paper a nonparametric multiple functional regression method where several functions would serve as predictors. Furthermore, we aim at extending this method to handle mixed discrete and functional explanatory variables. This should be helpful for situations where a subset of regressors are comprised of repeated observations of an outcome variable and the remaining are independent scalar or categorical variables. In [START_REF] Antoch | Functional linear regression with functional response: application to prediction of electricity consumption[END_REF] for example, the authors discuss the use of a functional linear regression model with a functional response to predict electricity consumption and mention that including the knowledge of special events such as festive days in the estimation procedure may improve the prediction.

The remainder of this paper is organized as follows. Section 2 reviews the multiple functional linear regression model and discusses its nonparametric extension. This section also describes the RKHS-based estimation procedure for the nonparametric multiple functional regression model. Section 3 concludes the paper.

Multiple functional regression

Before presenting our nonparametric multiple function regression procedure, we start this section with a brief overview of the multiple functional linear regression model [START_REF] Matsui | Regularized functional regression modeling for functional response and predictors[END_REF][START_REF] Valderrama | Forecasting pollen concentration by a two-Step functional model[END_REF]. This model extends functional linear regression with a functional response [START_REF] Ramsay | Some tools for functional data analysis[END_REF][START_REF] Ramsay | Functional data analysis[END_REF] to deal with more than one covariate and seeks to explain a functional response variable y(t) by several functional covariates x k (s). A multiple functional linear regression model is formulated as follows:

y i (t) = α(t) + p k=1 Is x ik (s)β k (s, t)ds + ǫ i (t), t ∈ I t , i = 1, . . . , n, (1)
where α(t) is the mean function, p is the number of functional covariates, n is the number of observations, β k (s, t) is the regression function for the k-th covariate and ǫ i (t) a random error function. To estimate the functional parameters of this model, one can consider the centered covariate and response variables to eliminate the functional intercept α. Then, β k (., .) are approximated by a linear combination of basis functions and the corresponding real-valued basis coefficients can be estimated by minimizing a penalized least square criterion. Good candidates for the basis functions include the Fourier basis [START_REF] Ramsay | Functional data analysis[END_REF] and the B-spline basis [START_REF] Prchal | Spline estimator for the functional linear regression with functional response[END_REF].

It is well known that parametric models suffer from the restriction that the inputoutput relationship has to be specified a priori. By allowing the data to model the relationships among variables, nonparametric models have emerged as a powerful approach for addressing this problem. In this context and from functional input-output data (x i (s), y i (t)) n i=1 ∈ (G x ) p × G y where G x : I s -→ R and G y : I t -→ R , a nonparametric multiple functional regression model can be defined as follows:

y i (t) = f (x i (s)) + ǫ i (t), s ∈ I s , t ∈ I t , i = 1, . . . , n,
where f is a linear operator which perform the mapping between two spaces of functions. In this work, we consider a slightly modified model in which covariates could be a mixture of discrete and continuous variables. More precisely, we consider the following model

y i (t) = f (x i ) + ǫ i (t), i = 1, . . . , n, (2)
where x i ∈ X is composed of two subsets x d i and x c i (s). x d i ∈ R k is a k × 1 vector of discrete dependent or independent variables and x c i (s) is a vector of p continuous functions, so each x i contains k discrete values and p functional variables.

Our main interest in this paper is to design an efficient estimation procedure of the regression parameter f of the model ( 2). An estimate f * of f ∈ F can be obtained by minimizing the following regularized empirical risk

f * = arg min f ∈F n i=1 y i -f (x i ) 2 Gy + λ f 2 F (3)
Borrowing the idea from Kadri et al. (2010a), we use function-valued reproducing kernel Hilbert spaces (RKHS) and operator-valued kernels to solve this minimization problem. Function-valued RKHS theory is the extension of the scalar-valued case to the functional response setting. In this context, Hilbert spaces of function-valued functions are constructed and basic properties of real RKHS are restated. Some examples of potential applications of these spaces can be found in [START_REF] Kadri | Function-Valued Reproducting Kernel Hilbert Spaces and Applications[END_REF] and in the area of multi-task learning (discrete outputs) see [START_REF] Evgeniou | Learning multiple tasks with kernel methods[END_REF]. Function-valued RKHS theory is based on the one-to-one correspondence between reproducing kernel Hilbert spaces of function-valued functions and positive operator-valued kernels. We start by recalling some basic properties of such Spaces. We say that a Hilbert space F of functions X -→ G y has the reproducing property, if ∀x ∈ X the evaluation functional f -→ f (x) is continuous. This continuity is equivalent to the continuity of the mapping f -→ f (x), g Gy for any x ∈ X and g ∈ G y . By the Riesz representation theorem it follows that for a given x ∈ X and for any choice of g ∈ G y , there exists an element

h g x ∈ F , s.t. ∀f ∈ F h g x , f F = f (x)
, g Gy We can therefore define the corresponding operator-valued kernel K(., .) ∈ L(G y ), where L(G y ) denote the set of bounded linear operators from G y to G y , such that

K(x, z)g 1 , g 2 Gy = h g1 x , h g2 z F It follows that h g1 x (z), g 2 Gy = h g1 x , h g2 z F = K(x, z)g 1 , g
2 Gy and thus we obtain the reproducing property K(x, .)g, f F = f (x), g Gy It is easy to see that K(x, z) is a positive kernel as defined below: Definition: We say that K(x, z), satisfying K(x, z) = K(z, x) * , is a positive operatorvalued kernel if given an arbitrary finite set of points {(x i , g i )} i=1,...,n ∈ X × G y , the corresponding block matrix K with K ij = K(x i , x j )g i , g j Gy is positive semi-definite.

Importantly, the converse is also true. Any positive operator-valued kernel K(x, z) gives rise to an RKHS F K , which can be constructed by considering the space of functionvalued functions f having the form f (.) = n i=1 K(x i , .)g i and taking completion with respect to the inner product given by K(x, .)g 1 , K(z, .)

g 2 F = K(x, z)g 1 , g 2 Gy .
The functional version of the Representer Theorem can be used to show that the solution of the minimization problem (3) is of the following form:

f * (x) = n j=1 K(x, x j )g j (4)
Substituting this form in (3) , we arrive at the following minimization over the scalarvalued functions g i rather than the function-valued function f

min g∈(Gy) n n i=1 y i - n j=1 K(x i , x j )g j 2 Gy + λ n i,j K(x i , x j )g i , g j Gy (5)
This problem can be solved by choosing a suitable operator-valued kernel. Choosing K presents two major difficulties: we need to construct a function from an adequate operator, and which takes as arguments variables composed of scalars and functions. [START_REF] Lian | Nonlinear functional models for functional responses in reproducing kernel hilbert spaces[END_REF] considered the identity operator, while in [START_REF] Kadri | Nonlinear functional regression: a functional RKHS approach[END_REF] the authors showed that it will be more useful to choose other operators than identity that are able to take into account functional properties of the input and output spaces. They also introduced a functional extension of the Gaussian kernel based on the multiplication operator. Using this operator, their approach can be seen as a nonlinear extension of the functional linear concurrent model [START_REF] Ramsay | Functional data analysis[END_REF]. Motivated by extending the functional linear regression model with functional response, we consider in this work a kernel K constructed from the integral operator and having the following form:

(K(x i , x j )g)(t) = [k x d (x d i , x d j ) + k x c (x c i , x c j )] k y (s, t)g(s)ds (6)
where k x d and k x c are scalar-valued kernels on R k and (G x ) p respectively and k y the reproducing kernel of the space G y . Choosing k x d and k y is not a problem. Among the large number of possible classical kernels k x d and k y , we chose the Gaussian kernel. However, constructing k x c is slightly more delicate. One can use the inner product in (G x ) p to construct a linear kernel. Also, extending real-valued functional kernels such as those in [START_REF] Rossi | Support vector machine for functional data classification[END_REF] to multiple functional inputs could be possible.

To solve the problem (5), we consider that G y is a real-valued RKHS and k y its reproducing kernel and then each function in this space can be approximated by a finite linear combination of kernels. So, the functions g i (.) can be approximated by m l=1 α il k y (t l , .) and solving (5) returns to finding the corresponding real variables α il . Under this framework and using matrix formulation, we find that the nm × 1 vector α satisfies the system of linear equation (K + λI)α = Y (7)

where the nm × 1 vector Y is obtained by concatenating the columns of the matrix (Y il ) i≤n, l≤m and K is the block operator kernel matrix (K ij ) 1≤i,j≤n where each K ij is a m × m matrix.

Conclusion

We study the problem of multiple functional regression where several functional explanatory variables are used to predict a functional response. Using function-valued RKHS theory, we have proposed a nonparametric estimation procedure which support mixed discrete and continuous covariates. In future, we will illustrate our approach and evaluate its performance by experiments on simulated and real data.

Acknowledgments H.K. is supported by Junior Researcher Contract No. 4297 from the the Nord-Pas de Calais region.


Mohamed Kharboutly

Alexandre Melis

Aude Bolopion

Nicolas Chaillet

Member, IEEE Michaël Gauthier

2D robotic control of a planar dielectrophoresis

come

I. INTRODUCTION

This article deals with the robotic control of a non-contact dielectrophoresis system which can be considered as an original robotic structure compared to the current industrial robot. The first industrial robot UNIMATE [START_REF] Devol | Programmed article transfer[END_REF] based on standard joints was commercialized in 1961 (see figure 1). Nowadays more than one million of robots are in use all over the world. In the 1980's the use of compliant structures in robotics [START_REF] Schutter | Compliant robot motion i. a formalism for specifying compliant motion tasks[END_REF] was started to enable high precision positioning making them, at present, the most widely used structure for microscale robots [START_REF] Bacher | Flexures for high precision robotics[END_REF], [START_REF] Richard | Concept of modular flexure-based mechanisms for ultra-high precision robot design[END_REF]. However, transmission of movement in such robots is obtained via the movements of mechanical parts which largely limits throughput due to inertial effects. In the 2000's, LightWeight Robots [START_REF] Hirzinger | Advances in robotics: The dlr experience[END_REF], [START_REF] Schreiber | The fast research interface for the kuka lightweight robot[END_REF] have been developed by KUKA [START_REF] Kuka | [END_REF] to reduce robot inertia. However, the impact of inertia is still important in the small scales (micro-nano) where the inertia of the object is highly negligible compared to the one of the robots. We propose to develop robots that use non-contact transmission of movement to manipulate micro-nano-objects. Besides eliminating the inertia of a robotic structure, this approach also eliminates friction and adhesion (between the tweezer and the component) which are highly detrimental to a robot performance and life time. These 'No Weight Robots' NWR are at the cross-road between parallel robot and current non-contact manipulation. Firstly, NWR consists of moving components by applying forces coming from several physical field sources which 1 FEMTO-ST Institute, AS2M departement. UMR CNRS 6174 -UFC/ENSMM/UTBM, 24 rue Alain Savary, 25000, Besanc ¸on, France, EU.

have a similar effect to parallel robotics [START_REF] Pierrot | Optimal design of a 4-dof parallel manipulator: From academia to industry[END_REF], [START_REF] Bellakehal | Force/position control of parallel robots using exteroceptive pose measurements[END_REF] where the platform is moved by several mechanical forces coming from several robotic legs. The use of non-contact forces, rather than mechanical forces, changes the robot design drastically. In this regard, existing robotic approaches cannot be transferred to NWR. Secondly, current noncontact manipulation has been achieved mostly by open loop for object positioning or self-assembly [START_REF] Pawashe | Multiple magnetic microrobot control using electrostatic clamping[END_REF][START_REF] Floyd | Two-dimensional contact and non-contact micro-manipulation in liquid using an untethered mobile magnetic micro-robot[END_REF][START_REF] Ivan | First experiments on magpier: A planar wireless magnetic and piezoelectric microrobot[END_REF][START_REF] Kharboutly | Modeling the trajectory of a micro particle in a dielectrophoresis device[END_REF][START_REF] Sun | Ac electrokinetic micro-and nano-particle manipulation and characterization[END_REF][START_REF] Kuzyk | Dielectrophoresis at the nanoscale[END_REF][START_REF] Arai | Synchronized laser micromanipulation of multiple targets along each trajectory by single laser[END_REF][START_REF] Onda | Multi-beam bilateral teleoperation of holographic optical tweezers[END_REF] . The only exception concerns laser trapping which has been experimented in closed-loop by Arai et al. [START_REF] Arai | Synchronized laser micromanipulation of multiple targets along each trajectory by single laser[END_REF], [START_REF] Onda | Multi-beam bilateral teleoperation of holographic optical tweezers[END_REF]. However, laser trapping induces forces around tens of picoNewtons limiting the achievable throughput. The dielectrophoresis proposed in this paper generate forces around thousand times higher [START_REF] Kharboutly | Modeling the trajectory of a micro particle in a dielectrophoresis device[END_REF], [START_REF]Predictive control of a micro bead's trajectory in a dielectrophoresis-based device[END_REF]. Providing robotic control strategies will enable active and reprogrammable trajectory control and guarantee the final position of a manipulated object. This paper introduces a numerical model of a micro-bead's behavior in a dielectrophoresis system in the next section. The 2D trajectory control based on an inverted model is described in the fourth section and experimental validation is presented in the last section. In this section, we present a 3D dielectrophoretic force simulator applied on a micro-bead, which will be used as reference system. Secondly, a 2D direct dynamic model based on the 3D dielectrophoretic force simulator, will be presented and used to establish a 2D control law presented in the next section.

A. Dielectrophoresis force similator

In order to compute the electric field and then the dielectrophoretic force applied to a micro-object in an electrode structure, a numerical simulator is needed. This numerical simulator must be able to compute the dielectrophoretic force generated by very complex geometries in a very short time. For one hand, corresponding analytic equations are very complex and hard to be established. For a second hand, the finite element modeling (FEM) solution is limited to a long computation time and specially when electric voltage changes frequently. Thus, we propose to use the hybrid numeric simulator proposed in [START_REF] Kharboutly | Modeling the trajectory of a micro particle in a dielectrophoresis device[END_REF] gathering the ability of the FEM solution to simulate complex electrodes geometry and the short computation time of the analytical equations. According to [START_REF] Hughes | Nanoelectromechanics in Engineering and Biology[END_REF], the dielectrophoretic force -→ F DEP applied to the micro-bead's center X(x, y, z) with respect to the electric field -→ E (X, U ) can be written as:

- → F DEP (X, U ) = 2πǫ m r 3 Re[K(ω)]∇( - → E 2 (X, U )), (1)
where

K(ω) = ǫ * p -ǫ * m ǫ * p + 2ǫ * m , (2)
and ǫ * p and ǫ * m are respectively the complex permittivity of the particle and the medium with:

ǫ * = ǫ + j σ ω , (3)
ǫ is the relative permittivity, σ is the conductivity and ω is the angular velocity of the electric field. Thus, if we consider a configuration of n electrodes, by applying n -1 sinusoidal electric voltages identified by there magnitudes U = [U 1 , ...U n-1 ] and there angular velocity ω, the electric field -→ E (X, U ) can be computed using the hybrid method described in [START_REF] Kharboutly | Modeling the trajectory of a micro particle in a dielectrophoresis device[END_REF]. This hybrid method consists in computing the electric field -→ E (X, U ) by integrating the surface charge density on the electrodes. In fact the electric charge density Q and the magnitudes of the applied voltages U on the electrodes are linearly related:

Q = n-1 i=1 (C i U i ), (4)
where U i is the magnitude of the applied voltage on the ith electrode and C i is the elementary inter-capacitance between the electrodes influenced by the ith electrode.

The inter-capacitance between the electrodes depends on only the geometric shape of the electrodes and the electric permittivity of the medium. The C i is simulated using FEM Once the matrix of the electric charge density Q is computed, the electric field can be calculated analytically in a point X(x, y, z) in the medium. In fact, with each value Q i,j of the computed matrix Q corresponds a x i,j , y i,j point on the electrodes (z i,j = 0 because of the electrodes are in the x, y plane). Thus, the expression of the electric field -→ E at the point X(x, y, z) is:

- → E (x, y, z) = i   j Q i,j - → r 4πǫ m - → r 3   , (5)
where r = [xx i,j , yy i,j , z], and the DEP force can be also computed analytically with respect to (1). The figure 3 resumes the DEP modeling simulator (DMS) block. The block's inputs are the geometric shape of the electrodes, the applied voltages and the micro-bead's current position. This block generates the computed x, y and z components of the dielectrophoretic force applied to the micro-bead in its center.

B. 3D direct dynamic model

The dynamic of a micro-bead in motion under dielectrophoretic force field, in a liquid medium is ruled by the following dynamic equation:

- → F DEP (X) + - → F Drag + - → P = m - → Ẍ , ( 6
)
where X is the space coordinates of the micro-bead's center X(x, y, z), Ẋ is its velocity, Ẍ its acceleration, -→ F DEP (X) is the applied dielectrophoretic force on the center of the micro-bead, -→ P is its apparent weight (sum of the weight and the buoyancy), m is its mass and -→ F Drag is the viscosity friction created on the micro-bead. Generally, in the microscale, micro-manipulation in a liquid medium with dynamic viscosity ν is characterized by a Reynolds number much smaller than 1. In this case the micro-bead's inertia impact is very small compared to the viscosity friction -→ F Drag . Thus the inertia term m -→ Ẍ can be neglected and the dynamic equation becomes:

- → F DEP (X) + - → F Drag ( Ẋ) + - → P = 0. (7)
In the micron scale the Stokes approach of the viscosity friction is valid, -→ F Drag ( Ẋ) becomes:

- → F Drag ( Ẋ) = -6πνR - → Ẋ , ( 8
)
where ν is the dynamic viscosity and R the radius of the micro-bead. The dynamic equation is thus:

- → Ẋ = - → F DEP (X) + - → P 6πνR . ( 9
)
The diagram in the figure 4 illustrates the 3D direct dynamic modeling. Having the applied electric voltages and the electrodes geometry as input, the direct modeling simulator computes the corresponding micro-bead's trajectory. In generally, the micro-bead's behavior in dielectrophoretic force field is characterized by its high dynamics and nonlinearity. This numeric simulator is experimentally validated in [START_REF] Kharboutly | Modeling the trajectory of a micro particle in a dielectrophoresis device[END_REF] where we have shown that the dynamics are very high and the time response of the micro-bead is less than 3ms.

Moreover the behavior of the micro-bead is subjected to a high nonlinearity and especially when the micro-bead approaches the electrodes.

C. 2D simplified dynamic model

In order to reduce the complexity of the computation, we will consider that the electrodes surface is planar in the x, y plane. The 3D dielectrophoretic dynamic modeling simulator is designed to run on a classic P C with high performance (typically GHz) and it is not optimized to be integrated directly into a controller card with lower calculation performance (typically M Hz). Thus, a reduction of the 3D simulator in 2D is proposed. We assume that the micro-bead will move only in a 2D horizontal plane parallel to the electrodes surface at a height equal to its radius R. The impact of this assumption is going to be discussed at the end of the paper. The 2D simulator uses a similar approach to the 3D DMS presented above. In this 2D model, a database of the elementary spacial force is created. This database links the 2D dielectrophoretic force directly to the applied voltages, which will reduce sufficiently the computation time. Using the linear relationship between the electric field -→ E and the applied voltages V , the dielectrophoretic force can be written as a second order equation with respect to the electric voltages. Using the following electric voltages vector (see in figure 5) :

V = [-u y , U ref -u x , u y , U ref + u x ]; (10)
the 2D dielectrophoretic force [F DEPx , F DEPy ] can be written as the following:

F DEPx = f 11 u 2 x + f 12 u 2 y + f 13 u x u y + f 14 U ref u x + f 15 U ref u y + f 16 U 2 ref F DEPy = f 21 u 2 x + f 22 u 2 y + f 23 u x u y + f 24 U ref u x + f 25 U ref u y + f 26 U 2 ref . ( 11
)
U ref is a reference voltage, u x and u y are the varying voltages and f i,j are spacial functions in x and y essentially dependent on the electrodes geometries. Discrete values of these functions will be computed in a x, y grid points using the 3D simulator and stored in a database and a quadratic interpolation is used to evaluate these functions in an arbitrary (x, y) point. Using this 2D numeric model, the 2D direct dynamic model becomes: The weight -→ P has been removed as its projection on the x, y plane is null. Consequently, the computation time is reduced and few arithmetic iterations are executed in a very short time, even with the interpolation procedure. Indeed, 40 CPU clock cycles are needed to compute the 2 components of the dielectrophoretic force in a grid point, and 180 CPU clock cycle in an interpolated position. In the other hand, using the 3D dielectrophoretic force computation, at least 10 4 CPU clock cycle are needed. Thus if we consider that the micro-bead's time response is 3ms and for a successful tracking at least 10 control sequence are generated, a controller card with 1M Hz clock takes 0.2ms to compute the dielectrophoretic force using the 2D model. The same controller card may takes more then 10ms when using the 3D dynamic modeling.

ẋ ẏ = 1 6πνR F DEPx F DEPy ( 12
)

III. 2D TRAJECTORY TRACKING

To control the micro-bead's trajectory in a 2D dielectrophoretic local periodic structure, an elementary control law for tracking trajectories is presented in this section where a micro-bead moves in one structure. The behavior of a micro-bead in a dielectrophoretic system is characterized by its high dynamics as presented in [START_REF] Kharboutly | Modeling the trajectory of a micro particle in a dielectrophoresis device[END_REF] and the nonlinearity with respect to the applied voltages as shown in the equation [START_REF] Floyd | Two-dimensional contact and non-contact micro-manipulation in liquid using an untethered mobile magnetic micro-robot[END_REF]. This elementary control law must takes into consideration this two problematics. Consequently a simple proportional integrator control is not sufficient especially when the micro-bead approaches the electrodes where the nonlinearity becomes very high. The analytic inversion of the 2D model ( 12) is not possible due to the strong coupling between the control variables u x and u y . One way to solve this problem is to use the Newton-Raphson numeric method which is able to find the values of the control variables, u x and u y to follow a reference trajectory.

A. 2D Inverse dynamic model

Newton-Raphson is a method for finding successively better approximations to the roots of a real-valued functions. By sampling the 2D dynamic model ( 12) and knowing the trajectory [x(t), ŷ(t)] with respect to the time we are able to compute the appropriate control variable u x (t) and u y (t) using the Newton-Raphson method as illustrated in the figure 6:

[ ̂x(t) ̂y(t)] Invers dynamics (Newton-Raphson) [ u (t) x u (t) y ] 3D simulator [ x(t) y(t) ]
Fig. 6. The Newton-Raphson method is used to find the control variables ux and uy

By sampling the dynamic equation ( 12) using a sampling period T we obtain:

xk+1 ŷk+1 = T 6πνR F DEPx (u xk , u y k ) F DEPy (u xk , u y k ) + x k y k ( 13
)
where xk+1 and ŷk+1 are the next trajectory point at the date kT . Applying the Newton-Raphson method to this model consists in finding iteratively a series of u x and u y . At the date kT we have:

u x k+1 u y k+1 = u x k u y k -J(u x k , u y k ) -1 F (u x k , u y k ) G(u x k , u y k )
where u x0 and u y0 are the last computed control variable, J is the jacobin matrix:

J(u x k , u y k ) =    ∂F (ux,uy) ∂ux ux k ,uy k ∂F (ux,uy) ∂uy ux k ,uy k ∂G(ux,uy) ∂ux ux k ,uy k ∂G(ux,uy) ∂uy ux k ,uy k    (14) and
F (u x , u y ) = F DEPx (u x , u y ) -6πνR(x k+1 -x k ) G(u x , u y ) = F DEPy (u x , u y ) -6πνR(ŷ k+1 -y k ) (15)
The iterations clasically stops when:

u x l+1 -u x l δ u and u y l+1 -u y l δ u
where δ u is an error threshold.

B. Numeric application and experimental results

Considering the electrodes geometry presented in the figure 7 submerged in a ultra pure water, where the trajectory of a micro-bead made of silica will be controlled. The table I contains the numeric values of the system physical parameters.

Firstly, we are tracking a square reference trajectory with 1s period, presented in the figure 8(a).

Applying the Newton-Raphson method, a series of u x and u y control variables are computed and presented in the figure 8(b). These computed electric voltages are transmitted to a digital-analogic convertor, then they are amplified and applied to the electrodes where the particle is placed (see figure 7). The micro-bead's position is captured by a high speed camera acquisition at 300 images per seconds.

The figure 9 shows the real trajectory of the micro-bead when applying the u x and u y series already computed. The relative error between the real trajectory and the reference is less then 8%. This results shows the ability to control the trajectory of a micro-bead in dielectrophoretic system using open loop control strategy.

IV. DISCUSSION

In order to show the limitation of our current control law a second experiment with a trajectory 10 times faster is studied (figure 10(a)). The 2D controller computes a new series of u x and u y presented in the figure 10(b).

The figure 11 shows the real trajectory of the micro-bead when applying the new u x and u y series. As this figure shows, the real trajectory does not follow the reference and the error is bigger then 50%. This is essentially due to the limitation of the 2D controller where the dynamic model is limited to 2D.

In reality the micro-bead moves in the three directions x, y and z but in the 2D dynamic mode, it supposes to move in a horizontal plane. The main difference between both experiments is the applied voltages u x and u y . In the second case, the computed voltage is greater then 50V . When applying these voltages, experiments show that the height of the micro-beads is significantly different from the plane assumed in the 2D dynamic model. Thus, in the first experiment, the real height of the micro-bead is very close to the plane of the 2D dynamic modeling consequently the error between the real trajectory and the reference is less then 8%. In the second experiment, the micro-bead's motion is 100µm higher than the plane of the 2D dynamic model, consequently the error between the real trajectory and the reference is bigger then 50%. However, the object trajectory computed using the 3D dynamic model is very close to the experimental measurement (see figure 11). This indicates that our 3D simulator is reliable and could be used to develop an new 3D control law based on a new simplified 3D dynamic model. The current method proposed in this paper has shown its relevance for low altitude trajectories, the extension of this approach in 3D will be developed in future works.

V. CONCLUSION

In order to control the trajectory of a micro-object for long distance and high speed, an elementary open-loop positioning control for a micro-bead is presented using dielectrophoresis. A 3D dielectrophoretic force simulator has been firstly presented. We have shown that the full 3D dynamic model is too complex to be introduced in a controller card. Thus a reduced 2D dynamic model based on the 3D dynamic model was developed. This model is then used to establish a 2D control law. We have shown that inverting this model cannot be done in an analythic way and we have proposed to use the Newton-Raphson numerical method, in order to compute the appropriate control variable with respect to a reference trajectory. Two experiments was presented including two reference trajectories with different speeds. Experiments have shown that the 2D controller has succeeded to track the low speed trajectory, and failed to track the high speed trajectory. This difference has been also discussed. Experiments also have confirmed that the original 3D dynamic modeling is reliable, and could be used in developing a new control law dedicated to 3D trajectories.

Fig. 1 .

 1 Fig. 1. Movement transmission used in robotics: (i) standard joints used in a majority of robots; (ii) compliant joints based on mechanical deformation used in high precision positionning systems; (iii) the third alternative: movement transmission based on non contact forces

  software. These simulations are executed in preprocessing which reduces the total time of the force computation. If we consider the planar electrodes drawn in the figure2(red lines), the number of electrodes n is equal to 4 and they are placed in the x, y plane. To compute the electric charge density Q with respect to the applied voltages U = [U 1 , U 2 , U 3 ], n -1 = 3 FEM simulations are required. The figures 2(a) and 2(b) show the elementary inter-capacitances C 1 and C 3 . The figure 2(c) shows how the electric charge density Q is analytically computed with respect to the applied voltages U = [75V, 0, 75V ] and the elementary inter-capacitances C 1 and C 3 .

Fig. 2 .

 2 Fig. 2. The electric charge density computed on the electrodes by applying the following electric voltages: U = [75V, 0, 75V ].

Fig. 3 .

 3 Fig. 3. DEP modeling simulator (DMS).

Fig. 4 .

 4 Fig. 4. A dynamic modeling and DMS are used to compute the microbead's 3D trajectory.

Fig. 5 .

 5 Fig. 5. Geometry of the electrodes and applied voltages: definition of control parameters ux and uy.

Fig. 7 .

 7 Fig. 7. Experimental electrode used to apply the dielectrophoretic motion. The square presents the reference trajectory.

(

  a) x(t) and ŷ(t) low speed trajectory. The computed ux and uy values for tracking the low speed trajectory.

Fig. 8 .

 8 Fig. 8. The control variable ux and uy computed for tracking the low speed trajectory.

Fig. 9 .

 9 Fig. 9. The real trajectory made by the micro-bead when applying the computed ux and uy for the low speed trajectory.

  x(t) and ŷ(t) high speed trajectory. The computed ux and uy values for tracking the high speed trajectory.

Fig. 10 .

 10 Fig. 10. The control variable ux and uy computed for tracking the high speed trajectory.

Fig. 11 .

 11 Fig.11. The real trajectory made by the micro-bead when applying the computed ux and uy for the high speed trajectory.

ACKNOWLEDGMENT

The authors would like to thank all the contributors in this collaborative work under the european project FAB2ASM (contract FoF-NMP-2010-260079): Efficient and Precise 3D Integration of Heterogeneous Microsystems from Fabrication to Assembly, (http://fab2asm.eu/). This article is also supported by the european project IRSES EC-NANOMAN, European Chinese exchange program project on Nanomanipulation.


Franc ¸ois

Nicolas Vauchelet
email: vauchelet@ann.jussieu.fr

Numerical simulation of a hyperbolic model for chemotaxis after blow up

Keywords: chemotaxis, non-local conservation equations, numerical analysis, finite volume scheme, measure-valued solutions 2010 AMS subject classifications: Primary: 35Q92, 35L60, 65M08; Secondary: 35D30

. In this paper, we investigate the numerical discretization of this system. A scheme based on a finite volume approach is proposed and the convergence of the numerical solution towards the unique duality solution is stated. Numerical simulations are provided that shows the behaviour of solutions after blow up.

Introduction

Chemotaxis is the phenomenom in which a population of cells rearranges its structure according to the behaviour of some chemical present in the environment. In this work, we focus on the following macroscopic model for chemotaxis:

∂ t ρ + ∂ x (a(∂ x S)ρ) = 0, (1)
a(∂ x S) = -c φ(c∂ x S), (2)
-∂ xx S + S = ρ. (3)
In this system, ρ(t, x) ≥ 0 denotes the density of cells at time t ≥ 0 and position x ∈ R, and S(t, x) is the concentration of the chemical, which is here produced by the bacteria at a rate proportional to their density, and diffuses in the system. The positive constant c corresponds to the individual velocity of cells, whereas a is a macroscopic velocity. The real-valued function φ describes the influence of the chemoattractant concentration S on the global motion of bacteria. When φ is nonincreasing, hence a nondecreasing, cells attract themselves, we are in the case of positive chemotaxis, the chemical is then called the chemoattractant, otherwise, we are in the repulsive case and speak about the chemorepellant. We complement this system with the boundary conditions ρ(t = 0, x) = ρ ini (x), lim x→±∞ ρ(t, x) = 0, lim x→±∞ S(t, x) = 0. (4)

Introducing the elementary interaction kernel K solving -∂ xx K + K = δ 0 , i.e. K = 1 2 e -|x| , the latter system reduces to the nonlocal scalar conservation law

∂ t ρ + ∂ x (a(∂ x K * ρ)ρ) = 0. ( 5
)
When a is nondecreasing, this equation is known as the aggregation equation. It is now classical that regular solutions blow up in finite time when K is not smooth (see e.g. [START_REF] Bertozzi | Blow-up in multidimensional aggregation equation with mildly singular interaction kernels[END_REF]). Thus measure solutions have to be considered, together with a suitable definition of the product a(∂ x K * ρ)ρ, when ρ is a measure. System (1)-( 3) can be obtained thanks to a hydrodynamical limit of the so-called Othmer-Dunbar-Alt model. This system has been introduced to take care of the run and tumble process: the motion of cells is due to the alternance of a swim phase in a given direction (run) and of a reorientation phase during which cells take a new direction for the next run (tumble). Denoting f (t, x, v) the distribution function at time t, position x ∈ R and velocity v ∈ {-c, c}, the kinetic system of equations writes in one space dimension (see e.g. [START_REF] Dolak | Kinetic models for chemotaxis: Hydrodynamic limits and spatiotemporal mechanisms[END_REF])

∂ t f ε + v∂ x f ε = 1 ε ((1 + φ(-v∂ x S ε ))f ε (-v) -(1 + φ(v∂ x S ε ))f ε (v)), (6)
-∂ xx S ε + S ε = ρ ε = f ε (v) + f ε (-v). ( 7
)
The real c represents the constant velocity of cells. The left-hand side describes the run phase whereas the right-hand side models the tumble phase. The operator T [S] := 1 2 (1 + φ(-v∂ x S)) is called the turning rate and corresponds to the rate of reorientation of cells. Obvioulsy we should have 0 ≤ T [S] ≤ 1. The parameter ε is a scaling factor which is assumed here to be very small (ε ≪ 1). This corresponds to the phenomenon of domminant taxis. Other scalings are possible, that lead eventually to drift-diffusion equations such as the classical Keller-Segel model.

The hydrodynamic limit consists in letting ε → 0, and it is easily checked that formally the solution (ρ ε , S ε ) of the kinetic system converges to some solution (ρ, S) of system (1)-(3). A complete proof of this result is given in [START_REF] James | Chemotaxis : from kinetic equations to aggregate dynamics[END_REF], together with a global-in-time existence and uniqueness of duality solutions for system (1)- [START_REF] Bouchut | Duality solutions for pressureless gases, monotone scalar conservation laws, and uniqueness[END_REF]. The main drawback of this method is that, since duality solutions are presently defined only in one space dimension, we are limited to x ∈ R. On the other hand, existence of measurevalued solutions for the aggregation equation ( 5) has been obtained in [START_REF] Carrillo | Global-in-time weak measure solutions and finite-time aggregation for nonlocal interaction equations[END_REF] in any space dimension. The authors make use of optimal transport technique: in this geometric approach, the solution appears as a gradient flow for the interaction energy. However, this latter tool is not so convenient to build numerical schemes and to obtain numerical simulations.

The aim of this paper is therefore to propose a numerical scheme based on a finite volume approach to get numerical simulations of solutions to the macroscopic model ( 1)-( 3). The main ingredient is a proper definition of the flux J = a(∂ x S)ρ. The outline is the following. In the next Section we recall some definitions and the existence result stated in [START_REF] James | Chemotaxis : from kinetic equations to aggregate dynamics[END_REF]. In Section 3 we present the scheme and state some of its main properties in the attractive case. Finally Section 4 provides numerical simulations, both in the attractive and the repulsive case.

Existence of duality solutions

Let us denote by M b (R) the set of bounded Radon measures and by P 1 (R) the set of nonnegative measures in M b (R) with finite first moment, that is R |x|dµ(x) < ∞. Duality solutions have been introduced in [START_REF] Bouchut | One-dimensional transport equations with discontinuous coefficients[END_REF] to solve scalar conservation laws with discontinuous coefficients. More precisely, it gives sense to measure valued solutions of the scalar conservation law

∂ t ρ(t, x) + ∂ x (b(t, x)ρ(t, x)) = 0, where b ∈ L ∞ ((0, T ) × R) satisfies the so-called one-sided Lipschitz (OSL) condition ∂ x b(t, .) ≤ β(t) for β ∈ L 1 (0, T ), in the distributional sense. ( 8
)
We refer to [START_REF] Bouchut | One-dimensional transport equations with discontinuous coefficients[END_REF] for the precise definition and general properties of these solutions. In this framework we define duality solutions of the studied system in the spirit of [START_REF] Bouchut | Duality solutions for pressureless gases, monotone scalar conservation laws, and uniqueness[END_REF]: 2) is satisfied in the weak sense :

Definition 2.1 We say that (ρ, S) ∈ C([0, T ]; M b (R)) × C([0, T ]; W 1,∞ (R)) is a duality solution to (1)-(3) if there exists a ∈ L ∞ ((0, T ) × R) and α ∈ L 1 loc (0, T ) satisfying ∂ x a ≤ α in D ′ (R), such that (i) for all 0 < t 1 < t 2 < T ∂ t ρ + ∂ x (aρ) = 0 in the sense of duality on ]t 1 , t 2 [, (ii) equation (
∀ ψ ∈ C 1 (R), ∀ t ∈ [0, T ], R (∂ x S∂ x ψ + Sψ)(t, x) dx = ψ(x) ρ(t, dx), (iii) a = a(∂ x S) a.e.
The OSL estimate suggests that the velocity field a has to be compressive, that is to satisfy the OSL estimate. Therefore we assume to be in the case of positive chemotaxis. Then cells attract themselves and a is non-decreasing; more precisely we assume

a ∈ C 1 (R), 0 ≤ a ′ ∈ L ∞ (R), ∀ x ∈ R |a(x)| ≤ c. (9)
The latter estimate means that the collective displacement of cells should not be faster than the velocity of each individual cell; it is a direct consequence of the fact that the turning rate

T [S] satisfies 0 ≤ T [S] ≤ 1.
Moreover we have a natural one-sided estimate on the potential S when ρ is nonnegative. In fact, from (3) ∂ xx S ≤ S. This estimate can be considered as an entropy estimate for the scalar conservation law (1) which is crucial for the proof of uniqueness of solutions (see Theorem 5.1 of [START_REF] James | Chemotaxis : from kinetic equations to aggregate dynamics[END_REF]). From now on, we will define by A the antiderivative of a which vanishes at 0. We have from the chain rule that, when it is defined, the product

a(∂ x S)ρ = a(∂ x S)(-∂ xx S + S) = -∂ x (A(∂ x S)) + a(∂ x S)S.
Then a natural definition of the flux is given by

J = -∂ x (A(∂ x S)) + a(∂ x S)S. (10)
We are now in position to state the existence and uniqueness result of [START_REF] James | Chemotaxis : from kinetic equations to aggregate dynamics[END_REF].

Theorem 2.2 Let us assume that ρ ini is given in P 1 (R) and that ( 9) is satisfied. Then, for all T > 0 there exists a unique duality solution (ρ, S) with 0 ≤ ρ ∈ P 1 (R) of ( 1)-( 3) which satisfies in the distributional sense :

∂ t ρ + ∂ x J = 0, ( 11
)
where J is the flux defined in (10). Moreover, there exists a universal representative, denoted a, such that a = a(∂ x S) a.e. and aρ = J, in the sense of measures.

Then, we have ρ = X # ρ ini where X is the backward flow corresponding to a(∂ x S).

Remark 1 This theorem is proved in [START_REF] James | Chemotaxis : from kinetic equations to aggregate dynamics[END_REF] with a weaker assumption than [START_REF] Perthame | PDE models for chemotactic movements: parabolic, hyperbolic and kinetic[END_REF]. However, this stronger assumption is needed here for the numerical analysis.

3 Numerical scheme

Discretization

Let us consider a uniform space discretization with step δx and denote by δt the time step, and set λ = δt/δx. Then t n = nδt and x i = x 0 + iδx. We assume that (ρ n i ) 0≤i≤N is an approximation of (ρ(t n , x i )) 0≤i≤N . We obtain an approximation ρ n+1 i of ρ(t n+1 , x i ) by using the following Lax-Friedrichs discretization of equations ( 11)-(10):

ρ n+1 i = ρ n i (1 -λc) + λ 2 c(ρ n i-1 + ρ n i+1 ) + λ 2 (J n i-1/2 -J n i+1/2 ), (12)
J n i+1/2 = - A(∂ x S n i+1 ) -A(∂ x S n i ) δx + a n i+1/2 S n i+1 + S n i 2 . ( 13
)
We recall that c denotes the constant modulus of the velocity of cells, and the velocity a(∂ x S) is discretized as

a n i+1/2 =      0 if ∂ x S n i+1 = ∂ x S n i , A(∂ x S n i+1 ) -A(∂ x S n i ) ∂ x S n i+1 -∂ x S n i otherwise, (14)
complemented with the standard centered finite difference

∂ x S n i+1 = S n i+2 -S n i 2δx
. We couple this equation with the following standard discretization of the equation for the chemoattractant (3)

- S n i+1 -2S n i + S n i-1 δx 2 + S n i = ρ n i . (15)
In order to avoid the treatment of boundary conditions, we assume that the solutions are compactly supported in the computational domain. Then from now on, we have that

ρ n 0 = S n 0 = S n 1 = J n -1/2 = 0 and ρ n N = S n N = J n N +1/2 = 0.

Numerical analysis

Before stating and proving our convergence result, we start by a Lemma which proves a CFL-like condition for the scheme : Lemma 3.1 Let us assume that (9) holds and that the condition

λ := δt δx ≤ 2 3c , (16)
is satisfied. Then the scheme defined in ( 12)-( 15) is nonnegative.

Proof. Let us assume that ρ n i ≥ 0. From ( 14)-( 15) we deduce that we can rewrite (13) as

J n i+1/2 = a n i+1/2 ρ n i+1 + ρ n i 2 . ( 17
)
Thus we can rewrite (12)

ρ n+1 i = ρ n i 1 -λc + λ 4 (a n i-1/2 -a n i+1/2 ) + λ 2 c + a n i-1/2 2 ρ n i-1 + λ 2 c - a n i+1/2 2 ρ n i+1 . (18)
Moreover, by assumption [START_REF] Perthame | PDE models for chemotactic movements: parabolic, hyperbolic and kinetic[END_REF] we have |a n i+1/2 | ≤ c for all i = 0, . . . , N . Therefore, if λ ≤ 2 3c , all the coefficients in front of ρ n i-1 , ρ n i and ρ n i+1 are nonnegative. We conclude that the scheme is nonnegative.

Let us define

ρ δ (t, x) = n∈N N i=0 ρ n i 1 [nδt,(n+1)δt)×[x i ,x i+1 ) (t, x),
and S δ , ∂ x S δ , J δ and a δ are defined in a similar way thanks to (S n i ) i , (∂ x S n i ) i , (J n i+1/2 ) i and (a n i+1/2 ) i .

Theorem 3.2 Let us assume that we are given ρ ini ∈ P 1 (R) and define ρ 0 i =

x i+1 x i ρ ini (dx) ≥ 0. Under assumption [START_REF] Perthame | PDE models for chemotactic movements: parabolic, hyperbolic and kinetic[END_REF], if ( 16) is satisfied, then the discretization (ρ δ , S δ ) converges towards the solution (ρ, S) of Theorem 2.2 as δt and δx go to 0.

Proof. We first notice that since the matrix of the linear system ( 15) is a M -matrix, we have that if (ρ n i ) i is nonnegative, then S n i ≥ 0 for all i = 0, . . . , N . Let us define M n i = δx i j=0 ρ n j and M n+1 i = δx i j=0 ρ n+1 j . Since the scheme ( 12) is conservative, we have

M n N = M 0 N . Clearly, ρ n i = (M n i -M n i-1
)/δx and from (17) we have

J n i+1/2 = a n i+1/2 (M n i+1 - M n i-1
)/(2δx). Then we deduce from (12) that

M n+1 i = (1 -λc)M n i + λ 2 c - a n i+1/2 2 M n i-1 + λ 2 c + a n i+1/2 2 M n i+1 . (19)
If λ satisfies (16), we have λc < 1 moreover, by assumption [START_REF] Perthame | PDE models for chemotactic movements: parabolic, hyperbolic and kinetic[END_REF],

|a n i+1/2 | ≤ c. Thus M n+1 i is a convex combination of M n i-1 , M n i and M n i+1 .
We deduce that provided condition ( 16) is satisfied, we have 0

≤ ρ n i = (M n i -M n i-1
)/δx and equation ( 18) implies a BV (R) estimate on (M n i ) i . Moreover, for all i = 0, . . . , N , we have

M n i ≤ M n N = M 0 N which provides a L ∞ estimate on (M n i ) i . Summing (15), we get M n i = δx i j=0 S n j -2∂ x S n i . (20)
We deduce that 0

≤ δx N j=0 S n j ≤ M n N . Then we have a L ∞ ∩ BV (R) estimate on (∂ x S n i ) i . Let us define M δ (t, x) = n∈N N i=0 M n i 1 [nδt,(n+1)δt)×[x i ,x i+1 ) (t, x).
Using standard arguments, we have a L ∞ ∩BV ((0, T )×R) estimate on M δ . It implies the convergence, up to a subsequence, of

M δ in L 1 loc (R + × R) towards a function M ∈ L ∞ ∩ BV ((0, T ) × R)
when δt and δx go to 0 and satisfy (16). By the same token we have the strong convergence in

L 1 loc (R + × R) of ∂ x S δ towards ∂ x S and by definition of ∂ x S i we have the convergence in L 1 loc (R + , W 1,1 loc (R)) of S δ towards S. Let us define, in the weak sense, ρ = ∂ x M ∈ M b (R). Obviously, noting that ρ n i = (M n i -M n i-1
)/δx, we deduce that ρ is the limit in the distributional sense of ρ δ . Passing to the limit in the equation (20) we deduce that S is a weak solution of equation ( 3).

Moreover we have

A(∂ x S n i+1 ) -A(∂ x S n i ) ∂ x S n i+1 -∂ x S n i = a(θ n i ), θ n i ∈ (∂ x S n i , ∂ x S n i+1 ).
Using assumption [START_REF] Perthame | PDE models for chemotactic movements: parabolic, hyperbolic and kinetic[END_REF], we deduce that the sequence (a δ ) δ is bounded in L ∞ , thus we can extract a subsequence converging in L ∞ -weak * towards a. From the L 1 loc convergence of (∂ x S δ ) δ , we deduce that a = a(∂ x S) a.e. Then, from (13), we have the convergence in the sense of distributions of J δ towards J = -∂ x (A(∂ x S)) + a(∂ x S)S a.e. Finally, taking the limit in the distributional sense of equation (12) we deduce that ρ is a distributional solution to (11)- (10). By uniqueness of this solution, we deduce that (ρ, S) is the unique duality solution of Theorem 2.2.

Numerical simulations

We present in this section some numerical results obtained thanks to an implementation of the scheme (12)-(15). Even though the scheme was designed to ensure convergence in the attractive case, we evidence that it gives interesting results in the repulsive case as well.

In all the following, the computational domain is assumed to be [-2.5, 2.5] and the velocity c is normalized to 1.

Attractive case

Let us consider the function a(x) = 2/π Arctan(10x), which satisfies clearly [START_REF] Perthame | PDE models for chemotactic movements: parabolic, hyperbolic and kinetic[END_REF]. In Figure 1, we plot the dynamics for the smooth initial data:

ρ ini (x) = e -10(x-1.25) 2 + 0.8e -20x 2 + e -10(x+1) 2 .
We notice that the blow-up occurs fastly. Then after a small time, solutions are formed by 3 peaks which can be considered as numerical Dirac masses. Then the Dirac masses move and collapse in finite time. This behaviour is very similar to the one observed in [START_REF] James | Chemotaxis : from kinetic equations to aggregate dynamics[END_REF] where a particle method has been implemented. Figure 1: Dynamics of the cell density for an initial data given by a sum of 3 regular bumps. We notice the fast blow up of solutions, then the obtained aggregates collapse together.

Repulsive case

It is very interesting to see the numerical results obtained by this scheme in the repulsive case, i.e. when the function a is non-increasing. In Figure 2, we present the result for a(x) = -2/π Arctan(10x) (left) and a(x) = -2/π Arctan(50x). The initial data is ρ ini (x) = e -10x 2 . In this case, the velocity x → a(∂ x S) does not satisfy the one-sided Lipschitz estimate [START_REF] Nieto | High field limit for Vlasov-Poisson-Fokker-Planck equations[END_REF]. Therefore we cannot define measure solutions in the sense of duality. However, we can prove, using the arguments in e.g. [START_REF] Nieto | High field limit for Vlasov-Poisson-Fokker-Planck equations[END_REF] that if ρ ini ∈ L 1 ∩ W 1,∞ (R), we have global in time existence of solutions in L 1 ∩ W 1,∞ (R). We observe in Figure 2 that support of the solution increases in time. It corresponds to the fact that cells repulse themselves. The particle scheme proposed in [START_REF] James | Chemotaxis : from kinetic equations to aggregate dynamics[END_REF] does not allow to get satisfying numerical results in the repulsive case. In fact in this latter scheme, we approximate the solution by a finite sum of Dirac masses. However, as we noticed above, duality solutions can not be defined in this case. Therefore we do not know how to solve the system for an initial data given by a sum of Dirac masses.

Finally, Figure 3 displays the dynamic of cells density in the repulsive case a(x) = -2/π Atan(10x) and for the initial data ρ ini (x) = e -10(x-0.7) 2 + e -10(x+0.7) 2 .

Figure 2 :

 2 Figure 2: Dynamics of the cell density in the repulsive case, i.e. for a non-increasing function a(x) = -2/π Atan(kx). Left : k = 10; Right : k = 50.

Figure 3: Dynamics of the cell density for an initial data given by a sum of 2 regular bumps in the repulsive case for a(x) = -2/π Atan(10x).


Lise Cary

Joël Casanova

Amira Mekni

Noureddine Gaaloul

Catherine Guerrot

Emmanuelle Petelet-Giraud

Effect of artificial recharge by treated wastewater on the quality of the Korba coastal aquifer (Cape Bon, Tunisia): insights from Boron isotopes

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

Effect of artificial recharge by treated wastewater on the quality of the Korba coastal aquifer (Cape Bon, Tunisia): insights from Boron isotopes.

L. Cary (1), J. Casanova (1), A. Mekni (2,3), N. Gaaloul (2), C. Guerrot (1), and E. (1) BRGM (Bureau de Recherches géologiques et minières), Orléans, France (l.cary@brgm.fr), ( 2 Mediterranean coastal aquifers are highly stressed, partly because of intensive agricultural activities in semi-arid or arid context. Water quality of groundwater is often degraded by different processes: salt water intrusion, diffuse pollution... In this context, the Korba aquifer on the eastern coast of Cape Bon (Tunisia) undergoes overexploitation. Since the 60's, reversal of hydraulic gradient and saline intrusion were studied (Ennabli, 1980;Kouzana et al. 2009aKouzana et al. , 2009b;;Kouzana et al., 2010) and the groundwater hydrogeological functioning was modeled (Kerrou et al., 2010;Panicoti et al. 2001aPanicoti et al. , 2001b;;Zghibi et al., 2011). Crucial for local agriculture, aquifer management led to the implementation of artificial recharge with treated wastewaters at the end of 2008. Since then, groundwater quality was monitored to trace the effectiveness of artificial recharge, based on boron isotopes, to better determine the different system components. Samples were taken from recharge control piezometers and in farmers' wells. Cation exchange was demonstrated to be one of the main processes controlling the groundwater quality. Sodium, potassium and boron were clearly in deficit compared to a mixing line with seawater, whereas calcium, strontium and bicarbonates were in excess. Treated wastewaters were concentrated in chloride and sodium and sometimes more than in wells, equivalent to a mixing of 6 to 14 % with sea water. The 2009 isotopic signatures in wells were clearly distributed on a mixing line between a nitrate-polluted and salinized pole especially constituted by abandoned wells and which had already reached equilibrium governed by cation exchange or sorption processes, and fresh groundwater. This mixing line explained the intermediate composition of a less salinized groundwater pole of intensively pumped wells. In 2010, in a context of low rainwater inputs, a general decrease of δ11B in the same wells may indicate mixing with deeper groundwaters and a continuous isotopic re-equilibration between the two previous end-members according to the different inputs, e. g. rainwater recharge or river flows. δ11B increased between 2010 and 2011 up to their 2009 value. Such a variation between 2010 and 2011 cannot be linked to a simple cation exchange but appear more reasonably linked to a higher rainfall recharge after a rainy winter. The mixing of treated wastewater was estimated at a maximum of 30 % with groundwaters and seemed very local. According to fluxes directions, artificial recharge could clearly not influence some wells composition. The impact of artificial recharge seemed thus very limited in terms of quality: recharge waters were salted and poorly contribute to refresh the system. Ultimately, the degradation of irrigated soils and crops diversity, currently limited by salt tolerance, will not improve. Moreover, B isotopes showed that artificial recharge seemed to influence some piezometers at vicinity but had low impact, if any, on close farmers' wells.

  ) INRGREF (National Institute of research in Rural Engineering of Water and Forestry), Tunis, Tunisia, (3) INAT (Agronomic National Institute of Tunisie), Tunis, Tunisia


S Blagojevic

Wolfram Kloppmann

L Cary
email: l.cary@brgm.fr

Z Jovanovic

R Stikic

An example of treated waste water use for soil irrigation in the SAFIR project

HAL is

The safe use of treated domestic wastewater for irrigation needs to address the risks for humans (workers, exposed via contact with irrigation water, soil, crops and food, consumers, exposed via ingestion of fresh and processed food), for animals (via ingestion of crops on soil), for the crops and agricultural productivity (via salinity and trace element uptake), for soil (via accumulation or release of pollutants) as well as for surface, groundwaters and the associated ecosystems (via runoff and infiltration, [START_REF] Kass | The impact of freshwater and wastewater irrigation on the chemistry of shallow groundwater: a case study from the Israeli Coastal Aquifer[END_REF][START_REF] Bouwer | Groundwater problems caused by irrigation with sewage effluent[END_REF]. In this context, the European FP6 SAFIR project (Safe and High Quality Food Production using Low Quality Waters and Improved Irrigation Systems and Management) investigates the geochemical quality of the root zone soil, knowing it is the main transit and storage compartment for pollutants. The type of reaction (sorption, co-precipitation. . . ) and the reactive mineral phases also determine the availability of trace elements for the plant and determine the passage towards crops and products. Reactions of the infiltrating water with the soil solid phase are important for the solute cycling, temporary fixation and remobilisation of trace pollutants. Therefore the soil water quality was directly or indirectly assessed. Direct measurements of soil water were made through porous cups.

The experiments were carried out during the growing season of 2006, 2007 and 2008 in a vegetable commercial farm, located at 10 km north of Belgrade. The soil is silty clayey, and developed on alluvial deposits. It was classified as humogley according to USDA Soil Classification. The climate of the field side is a continental type with hot and dry summers and cold and rainy winters. As in the rest of Serbia, farm suffers from water deficits during the main growing season. The initial soil quality was assessed through a sampling campaign before the onset of first year irrigation; the soil quality was then monitored throughout three years. Soil sampling focused on the fully irrigated plots because the potential impact of irrigation water quality on soil and plant quality are expected higher for fully irrigated soils compared to other irrigation strategies. Samples were taken within the soil volume of potential influence around each of the drip emitters. Potato (Solanum tuberosum) variety Liseta was used for investigation. The seeds tubers were planted in the similar period in all three seasons (middle of Spring) at the depth of 10 cm. Two irrigation methods were used in all three seasons: drip and furrow systems. Water for irrigation was supplied from a channel which is located 100m away from the experimental field. For all experiments, three sampling campaigns were foreseen for each of the three irrigation seasons: at pre-planting, at the end of irrigation, and at harvest. After three campaigns, the results show a variability of the elements concentrations in water and soil between the three years. The soil appears significantly depleted in CaO (a mean of -40 %), MgO (-20%), Na (-30%), and Sr (-10%) and Pb (-12%). On the contrary, concentrations of Mn, Ni, V and Li slightly increase (+15 to 20%) whereas SiO2, Al2O3, Fe2O3, Cu and Cr do not significantly increase (a mean of + 10%). Knowing that potatoes do need between 40 to 50 kg per ha of CaO and 15 to 30 kg per ha of MgO [START_REF] Soltner | Les grandes productions végétales[END_REF], potato absorption of Ca and Mg may be the main sink for both elements. A statistical analysis (ACP) shows precisely a Ca-Mg-Sr pole which explains more than 90 % of the second component; the first component being explained by Al2O3, SiO2, Fe2O3 and TiO2 at the same percentage.

  Antonious, G.F. and Snyder, J.C., 2007. Accumulation of heavy metals in plants and potential phytoremediation of lead by potato, Solanum tuberosum L. Journal of environmental science and health. Part A, Toxic/hazardous substances & environmental engineering, 42(6): 811-816.