T Biben

L Joly

Wetting on nanorough surfaces

We present in this letter a free-energy approach to the dynamics of a fluid near a nano-structured surface. The model accounts both for the static phase equilibrium in the vicinity of the surface (wetting angles, Cassie-Wenzel transition) and the dynamical properties like liquid slippage at the boundary. This method bridges the gap between phenomenological phase-field approaches and more macroscopic lattice-Boltzmann models.

Viscous dissipation is a major problem in micro or even nanofluidic systems [1]. Large pressure gradients are indeed necessary to produce a significant flow at small scale. This effect is well illustrated by the Poiseuille law relating the flux to the size of the duct: for a cylindrical duct of radius R, one can check that the flux decreases like R 4 when a given pressure gradient is applied. The Poiseuille law is however obtained by assuming no-slip boundary conditions (BC) for the liquid at the walls, and recent experiments have shown that this might not always be the case [START_REF] Lauga | Handbook of experimental fluid dynamics[END_REF]. In the presence of partial slip BC, one can expect the situation to be less penalizing, specially if the slip length associated to the boundary is of the order of the duct radius or even larger [START_REF] Bocquet | [END_REF]. Superhydrophobic surfaces have thus attracted a considerable interest recently as potential candidates for highly slipping surfaces: if the liquid is repelled by the surface, one can imagine that a thin layer of gas (air or vapor) at the surface could produce a lubricating effect. The situation is however more complex since natural or artificial superhydrophobic surfaces are extremely rough: roughness is indeed a key ingredient in superhydrophobicity, since it favors the trapping of air or vapor bubbles at the boundaries [4]. This effect is due to capillarity, and the question to understand the coupling between the three ingredients, hydrodynamics, capillarity and surface roughness is a true challenge. From the experimental point of view, the situation is still very controversial, with measured slip lengths varying on several orders of magnitude [START_REF] Bocquet | [END_REF]. It is thus crucial to have numerical models able to investigate the static and the dynamic properties of a fluid at a structured boundary.

Several approaches have been considered (see [START_REF] Bocquet | [END_REF] for a recent review), ranging from molecular dynamics at the nanometric scales [5] to hydrodynamic models (including the lattice-Boltzmann approach [6]) at micronic scales [7,8,9,10,11,12], and a very recent phase-field model at intermediate scales [13]. The method we shall consider in this letter is intermediate between phasefield phenomenological approaches, not accounting for the specificities of the gas phase (low viscosity and high compressibility) and the more macroscopic lattice Boltzmann approaches. It accounts explicitly for a possible liquid-vapor coexistence in the vicinity of a structured surface, and also accounts for an hydrodynamic flow with liquid slippage at the surface. This model is able to re-produce both static properties such as wetting angles or Cassie-Wenzel states, but also to predict the effective slip lengths of structured interfaces by accounting explicitly for the low viscosity gas layer that forms at the boundary. Moreover, this approach gives access to the intrusion dynamics of the fluid inside the pores of the structured surface, which is an important phenomenon when the liquid is submitted to large pressures.

The model is constructed as follows:

∂ρ ∂t + ∇ • (ρv) = -C δΩ δρ or = C ′ ∆ δΩ δρ (1)
is the transport equation for the local density field ρ(r, t), where t is time and r denotes the position; v(r, t) is the local velocity field and Ω[ρ] = F [ρ] -µ ρ(r)dr is the grand-potential functional associated to the thermodynamic equilibrium between the liquid and the vapor phase (F [ρ] is the free-energy). ρ is thus a physical order parameter here, and not simply a mathematical object to identify the phases as in usual phase-field approaches [13]. The convective transport of the density field is included in the left term, while the local thermodynamic equilibrium is contained in the right terms: the first writing corresponds to the Allen-Cahn model [14], interesting for the study of the phase diagram (non-conserved dynamics, fixed chemical potential µ), and the second writting is the Cahn-Hilliard model [15], corresponding to a more realistic conserved dynamics, that we shall use to investigate the non-equilibrium properties. C ′ is related to the molecular diffusion constant D in a bulk phase by the relation: D = C ′ f "(ρ) where f " denotes the second derivative of the bulk free energy f (ρ) at the bulk density ρ. In this work we used D = 2300µm 2 /s in the liquid phase, to match the auto-diffusion constant in pure water (D = 2270µm 2 /s). The free energy functional is

F [ρ] = dr f (ρ) + w 2 2 |∇ρ| 2 + ρV wall (r) , with f (ρ) = k B T ρ log ρb
1-ρb -1 -aρ 2 the free energy density of the bulk van der Waals theory (a is the mean field attractive energy between two atoms, and b is the close packing inverse density) . The square gradient term accounts for the natural thickness of the liquid-vapor interface and gives rise to the surface tension γ LV = w ρL ρV 2(f (ρ) -µρ + P )dρ [START_REF] Rowlinson | Molecular Theory of Capillarity[END_REF] where P is the bulk pressure, ρ L and ρ V are the coexisting liquid and vapor densities respectively. w is related to the in-terfacial thickness. V wall (r) is the interaction potential applied by the wall. Equation (1) has to be solved with the corresponding equation for the velocity field:

∂ρ m v ∂t + ∇ • (ρ m vv) = ∇ • σ -ρ∇ δF [ρ] δρ + f wall , (2)
where ρ m (r, t) is the local mass density (ρ m (r, t) = M ρ(r, t), with M the molecular mass), σ = η(ρ)(∇v + ∇v t ) is the local viscous stress tensor (with η(ρ) the shear viscosity). Although more accurate prescriptions can be considered, we use in this work the simple ansatz η(ρ) = ρ m ν, where the kinematic viscosity ν is assumed to be the same in the liquid and the vapor phase. The second term -ρ∇ δF [ρ] δρ is the thermodynamic force field applied by the density field ρ(r) to the flow. This term contains both the compressibility of the fluid, fixed by the van der Waals theory, and the capillary force when interfaces are present. The last term of equation ( 2) accounts for the interaction with the wall.

We shall now discuss the interactions between the fluid and the wall. The two fields ρ(r, t) and v(r, t) are defined in the whole parallelepipedic resolution box, and are thus defined and solved inside the walls. While the static properties of the wall are controlled by V wall (r) (expulsion of the fluid from the wall and wetting properties), the boundary conditions for the velocity field are controlled by f wall . A prerequisite for a correct description of the dynamical behavior of the liquid at superhydrophobic surfaces is to model properly the intrinsic slip properties of the liquid on a flat solid surface. Our prescription for f wall ensures partial slip BC in this case. These boundary conditions correspond to the continuity of the tangential stress σ z at the fluid/solid interface: σ z in the fluid equals -λv (z s ), the friction stress applied by the wall, where z is the normal direction, the tangential direction under consideration, λ is the friction coefficient and v (z s ) the slip velocity at the wall (z = z s ). The partial slip BC thus expresses:

v (z s ) = 1 λ σ z = η λ ∂v ∂z z=zs = b ∂v ∂z z=zs , (3)
where b = η/λ is the slip length, and η is the shear viscosity of the fluid. If b, η and the slip velocity are known, the stress at the interface is -ηv (z s )/b. The bare interaction between an atom of the wall and an atom of the fluid is the Lennard-Jones potential in our model: V LJ (r) = 4ǫ (σ/r) 12 -(σ/r) 6 . The global potential applied by the wall on a single atom of the fluid is V wall (r) = ρ wall (r ′ )V (|r ′ -r|)dr ′ , where ρ wall (r) is the density of atoms in the wall (assumed to be a constant inside the wall and zero elsewhere), and V (r) is the bare interaction potential [V LJ (r)] multiplied by the pair distribution function between an atom of the wall and an atom of the fluid. The main property of this pair distribution function is to vanish when r → 0 since two atoms cannot overlap, while V LJ (r) → +∞ in this limit. The product of these two functions thus goes to a well defined limit in r = 0. As a result, the potential energy inside the wall is large, but finite. We shall denote by V cut this energy and use the simple ansatz: V wall (r) = min(V cut , ρ wall (r ′ )V LJ (|r ′ -r|)dr ′ ). This prescription fixes the interaction to be of the Lennard-Jones type close to the boundary with the fluid, and to be V cut inside the wall. V cut fixes the value of the fluid density inside the wall, chosen to be negligible as compared to the density of the vapor phase (ρσ 3 10 -3 ). In this limit, the results become independent of V cut .

The force field f wall exerted by the wall on the fluid is chosen as a friction force:

f wall = -kǫρ(r)v(r) dr ′ ρ wall (r ′ )e -3 2 (r ′ -r) 2 σ 2
. This force is proportional to the interaction energy ǫ between the solid and the fluid [START_REF] Bocquet | [END_REF][START_REF] Barrat | [END_REF]. This prescription enabled us to reproduce correctly the evolution of the BC with the wetting properties. In particular, we obtain a transition from a no-slip BC in a wetting situation (θ ∼ 0 • ) to a partial slip BC with b ∼ 10 nm (depending on the value of k) in a non-wetting situation (θ ∼ 120 • ).

The two equations ( 1) and ( 2) can be solved numerically on a regular lattice (cubic unit cell in 3D or square in 2D) with periodic BC at the edges of the resolution box on the appropriate quantities, as we shall discuss. The wall is defined in the resolution box by the function ρ wall (r). We consider a slab geometry corresponding to a fluid confined between two walls, the bottom wall is textured with a periodic structure (crenels in this study), while the top wall can either be structured in static studies (Fig. 1-c), or planar in dynamic situations (Fig. 1-d). The density ρ(r) is periodic in the slab (x-y) direction, but also in the z direction since it goes to a constant value (nearly vanishing) inside the walls. The velocity field is periodic in the x-y direction, but not in the z direction where a linear shear v s (r) is applied. In this case, the periodic BC is taken on u(r) = v(r)v s (r) in any directions. By using a simple Euler scheme, equations ( 1) and ( 2) can be solved iteratively starting from an initial configuration for ρ(r) and v(r) . Taking advantage of the periodicity of ρ(r) and u(r), we use an implicit method in the Fourier space to improve the numerical stability, but simpler schemes can be used as well.

To illustrate the ability of the model to account for the small scale physics of the interfaces, we consider two situations. The first study corresponds to a static case, where v s (r) = 0. The wetting properties can be probed either by measuring the contact angles of a drop placed on the surface, or more simply by measuring the three surface tensions of the liquid-solid (γ SL ), vapor-solid (γ SV ) and liquid-vapor (γ LV ) interfaces. The model presented here gives a direct access to the free-energy or the grandpotential (an advantage compared to molecular simulations); the surface tensions can thus be measured directly by considering a fluid confined between two structured walls as depicted in Fig. 1-c. Provided the distance between the walls is large compared to the typical relaxation length of the density profile, the walls can be considered as independent. The relaxation length scale of the density profile is given by the interfacial thickness w, and we use a distance between the two walls of the order of 50w , which ensures an accurate determination of the surface tensions.The liquid-solid [resp. vapor-solid] surface tension is measured by starting the simulation with a liquid [resp. vapor] phase between the walls. For planar interfaces this prescription is very efficient and the two surface tensions can be measured in the full range of variation of the interaction potential V LJ , which is done by varying ǫ. From these two quantities, and the knowledge of γ LV , the contact angle with the planar surface can be calculated through the quantity Γ = γSV-γSL γLV . Γ corresponds to cos θ when |Γ| ≤ 1 and is an increasing function of ǫ (attractions favor wetting of the liquid phase). We shall consider Γ rather than ǫ in the discussions, since the contact angle is the experimentally relevant quantity. When the wall is not planar, and has a crenel shape as depicted in Fig. 1-a, the situation is more complex since various states can be observed. We used in this case two different types of initial configurations: Wenzel states corresponding to a monophasic system, or Cassie states corresponding to a diphasic situation (see Fig. 1-c). Of course, the equilibrium result should not depend on the initial state, this prescription gives us however a way to probe metastability in the system. The wetting properties of the structured wall can be measured as well by determining the effective surface tensions γ eff SV and γ eff LV for the equilibrated configurations, the variation of the corresponding wetting parameter Γ eff is plotted in Fig. 2. This figure illustrates very well the Cassie-Wenzel transition between the imbibed and super-hydrophobic states [4]. Moreover, several important features can be noted from these results: first we observe a strong metastability, with the presence of superhydrophobic metastable states, and next, we observe an asymmetry in the wetting properties when we go from the non-wetting region to the wetting one. The macroscopic model describing the transition between Cassie and Wenzel states is compared to the numerical results in Fig. 2. This model is quite accurate in the non-wetting region, which is not a surprise since we have the same geometry: the interfaces are planar at coexistence. In the wetting regime on the contrary, we observe discrepancies with the macroscopic theory. The difference comes from corner effects: the interaction potential between the fluid and the wall has particular values at the corners that generate a corner energy. These corner effects are mainly visible when a liquid phase fills the crenels (wetting situation), in gas phases the low value of the density leads to negligible contributions. If ǫ U and ǫ D are the excess free-energies in the liquid phase at the upper and the lower corners (see in Fig. 1-c the points U and D), the Wenzel theory can be modified as follows:

Γ eff = Γ(1 + 2d/L) -2 ǫ U + ǫ D Lγ LV (4)
where ǫ U and ǫ D are measured independently. This theory reproduces very well the effective contact angle (Fig. 2). Please note that ǫ U and ǫ D are not constant, they vary slightly with cos θ, which is visible in Fig. 2. We now turn to the dynamical properties of the surfaces. The dynamics of a fluid close to a textured wall can be accounted for using effective partial slip BC [START_REF] Bocquet | [END_REF], where both the slip length b and the hydrodynamic position of the wall z s are unknown. For a smooth surface, it has been shown that the BC applies at about one molecular layer inside the fluid, for both no-slip and slip cases [START_REF] Barrat | [END_REF]. However for superhydrophobic surfaces, it is not clear where the BC should apply. In this work we therefore measured both b and z s , by probing our system using both Couette and Poiseuille flows [START_REF] Barrat | [END_REF]. The geometry of the system corresponds to Fig. 1-d), where the top wall is planar and the bottom one is structured. This hybrid geometry allows for the determination of the effective boundaries for both the planar and the textured interfaces at the same time. The figure corresponds to a Couette flow, for which the top wall moves with a velocity U perpendicular to the crenels, while the bottom wall is fixed. The Poiseuille situation corresponds to fixed walls, but a pressure gradient is applied. We shall investigate the variation of both b and z s as a function of the liquid pressure. Rather than the pressure itself, we can equivalently fix the average density in the slab by using a conserved dynamics, and measure the pressure afterwards. We shall use δ , the recess length of the liquid in the crenel (taken from the top of the crenels, as depicted in Fig. 1d), as the control parameter. δ can either be negative or positive. The slip length of the planar (bare) liquid-wall interface is taken to be around 13 nm at coexistence, and is observed to vary very weakly with pressure (Fig. 3). On the contrary, the slip length of the structured wall varies in larger proportions, between 35 nm for the upper limit down to 0 for the lower one. Pressure has the effect to reduce slippage, and even to cancel it. The effective slip length crosses the bare one when δ c ≃ 2 nm, and is lower above. This result is largely independent of the crenel depth since the depth we consider in the present simulation (30 nm) is large compared to the crossover length δ c .

To conclude, the model presented here is very flexible and can account for complex BC. Although we focussed on crenelated walls, we can treat any type of geometry, 3D shapes such as posts or random surfaces, or conic shapes to model surfaces coated with nanotubes [18]. The dynamics on crenelated surfaces revealed the strong sensitivity of the effective slip length to the recess of the meniscus δ inside the crenels, cancelling slippage as soon as δ ≃ 2 nm. This critical value increases with the periodicity L, although it remains a small fraction of it. However, increasing L is not necessarily a good idea since the corresponding intrusion pressure will decrease. Although nanostructured surfaces have lower slip lengths than microstructured ones in general, the slip properties resist better to pressure. There is thus a compromise to find between slippage and pressure resistance.

Figure 1 :

 1 Figure 1: Schematic views of the systems under consideration.

Figure 2 :

 2 Figure 2: Wetting properties of a crenelated surface (see text for details).

Figure 3 :

 3 Figure 3: Evolution of the slip length b with the recess δ of the meniscus inside the crenel, for a flow perpendicular to the crenel. Inset: Evolution of the shear plane position zs as a function of δ. For the SH surface, zw is defined as the top of the crenels.


J.-C Aval
email: aval@labri.fr

X Viennot
email: viennot@labri.fr

THE PRODUCT OF TREES IN THE LODAY-RONCO ALGEBRA THROUGH CATALAN ALTERNATIVE TABLEAUX

The aim of this note is to show how the introduction of certain tableaux, called Catalan alternative tableaux, provides a very simple and elegant description of the product in the Hopf algebra of binary trees defined by Loday and Ronco. Moreover, we use this description to introduce a new associative product on the space of binary trees.

Introduction

Loday and Ronco defined in [START_REF] Loday | Hopf algebra of the planar binary trees[END_REF] an interesting Hopf algebra structure on the linear span of rooted planar binary trees. This algebra is defined as a sub-algebra of the Malvenuto-Reutenauer Hopf algebra of permutations. Let S n be the symmetric group and k be a ground field. We denote by k[S n ] the group algebra. Malvenuto and Reutenauer construct in [START_REF] Malvenuto | Duality between quasi-symmetric functions and the Solomon descent algebra[END_REF] a Hopf algebra structure on

k[S ∞ ] = n≥0 k[S n ].
It is worth to recall here that the Malvenuto-Reutenauer algebra contains the sum of Solomon descent algebras Sol ∞ = n≥0 Sol n with Sol n of dimension 2 n-1 .

In [START_REF] Loday | Hopf algebra of the planar binary trees[END_REF], Loday and Ronco define a sub-Hopf algebra of k[S ∞ ]:

k[Y ∞ ] = n≥0 k[Y n ]
where Y n is the set of planar binary trees with n internal vertices. The aim of this work is to present a very simple presentation for the product of two trees in k[Y ∞ ] through the use of Catalan alternative tableaux. These objects were introduced by X. Viennot [START_REF] Viennot | Alternative tableaux, permutations and partially asymmetric exclusion process, Isaac Newton institute[END_REF] as a special case of alternative tableaux, which are in bijection with permutation tableaux. Permutations tableaux were introduced by E. Steingrimsson and L. Williams [START_REF] Steingrimsson | Permutation tableaux and permutation patterns[END_REF], as a subclass of Γ-diagram defined by A. Postnikov [8]. This notion was used by S. Corteel and L. Williams [START_REF] Corteel | Tableaux combinatorics for the asymmetric exclusion process[END_REF][START_REF] Corteel | Tableaux combinatorics for the asymmetric exclusion process and Askey-Wilson polynomials[END_REF] in the study of the physical model named PASEP (partially asymmetric exclusion process), see for example the seminal paper by B. Derrida and al. [START_REF] Derrida | Exact solution of a 1D asymmetric exclusion model using a matrix formulation[END_REF]. These tableaux are also related to the study of total positivity for Grassmannian [START_REF] Williams | Enumeration of totally positive Grassmann cells[END_REF]. Both permutation and alternative tableaux are in bijection with permutations, see for example P. Nadeau [START_REF] Nadeau | The structure of alternative tableaux[END_REF]. The advantage of alternative tableaux is to preserve the symmetry between rows and columns. This new interpretation of the Loday-Ronco product motivates the introduction of a new associative product, that we call the # product, on the space of binary trees. This new product is studied in [START_REF] Chapoton | Some dendriform functors[END_REF] by F. Chapoton. This article is constructed as follows: in Section 2 we recall the definition of the Loday-Ronco algebra, then we introduce in Section 3 the Catalan alternative tableaux and prove that they are in bijection with binary trees; we state and prove in Section 4 the main result of this work, and we introduce the new # product in Section 5.

The Loday-Ronco Hopf algebra

We recall the definition of the Loday-Ronco product of binary trees. Since this product is inherited from the Malvenuto-Reutenauer product of permutations, we shall first recall the definition of the product in k[S ∞ ], denoted by * . We refer to [START_REF] Malvenuto | Duality between quasi-symmetric functions and the Solomon descent algebra[END_REF] for more details and only recall briefly the definition.

Let u = u 1 u 2 . . . , u k be a k-tuple of distinct integers. We define the standardization of u and denote it by Std(u) as the unique permutation σ ∈ S k that preserves the relative order of the u i 's, i.e. Remark 1. The product * that we consider is sometimes known as the product in the dual Malvenuto-Reutenauer algebra. But it is the one used in [START_REF] Loday | Hopf algebra of the planar binary trees[END_REF] to define the Loday-Ronco algebra, that we shall now describe.

σ i < σ j ⇐⇒ u i < u j .
Let Y n denote the set of binary trees with n internal vertices. We recall that the cardinality of Y n is given by the n-th Catalan number

C n = 1 n+1 2n n .
Let Ỹn denote the set of increasing binary trees, i.e. of binary trees such that each internal vertex has a distinct label in {1, . . . , n}, and such that the labels increase along the tree.

It is well known that increasing binary trees are in bijection with permutations: to obtain the permutation from the tree, you just have to read the labels from left to right.

Below is an example of a plane binary tree with 8 internal vertices, with a increasing binary tree with the same underlying tree, and with the corresponding permutation σ ∈ S n . We denote by Ψ : S n → Y n the composition of the bijection S n ≃ Ỹn with the projection Ỹn → Y n which consists in forgetting the labels. The induced linear map We also define for any tree T the set

Ψ : k[S n ] → k[Y n ] has a linear dual Ψ * : k[Y n ] → k[S n ]
Z T = {σ ∈ S n / Ψ(σ) = T } so that Ψ * (T ) = σ∈Z T σ.
The inclusion map Ψ * gives rise to a graded linear map

Ψ * : k[Y ∞ ] → k[S ∞
] and the main result in the construction of the Loday-Ronco algebra may now be stated as (Theorem 3.1 in [START_REF] Loday | Hopf algebra of the planar binary trees[END_REF]):

Theorem 2. The image of the inclusion map

Ψ * : k[Y ∞ ] → k[S ∞ ] is a sub-Hopf algebra of k[S ∞ ]. So, k[Y ∞ ]
inherits a structure of Hopf algebra.

Trees and Catalan alternative tableaux

We now present the Catalan alternative tableaux. Let us denote by N the set of nonnegative integers. A Catalan alternative tableau in given by • a path in N × N from {0} × N to N × {0} made of (0, 1) and (1, 0) steps. The length of the path is called the size of the tableau, and the cells below the path are simply called the cells of the tableau. The path defining the tableau can be called the shape of the tableau. • a set of blue and red dots in the cells of the tableau such that:

(1) there is no dot below a red dot;

(2) there is no dot on the left of a blue dot;

(3) any cell of the tableau is either below a red dot, or on the left of a blue dot. Let us give an example of a Catalan alternative tableau of size 23.

0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1
It is possible to directly check that Catalan alternative tableaux are enumerated by Catalan numbers (whence their name), but we shall use the following proposition, more adapted to our context. Proposition 1. The Catalan alternative tableaux of size n -1 are in bijection with binary trees with n internal nodes.

Proof. We refer to [START_REF] Viennot | Catalan tableaux and the asymmetric exclusion process[END_REF] (algorithm 2.2) for a formal proof and give here only the idea of the construction. In fact, in that paper, the algorithm was given in term of "Catalan permutation tableaux", the subclass of permutation tableaux corresponding to Catalan alternative tableaux, and discussed in [START_REF] Steingrimsson | Permutation tableaux and permutation patterns[END_REF].

We start with Catalan alternative tableau of size n -1 and rotate it:

The thick line represents the tree under construction. Recursively, for any "Up-Down" pattern in the thick line, we operate a shift, and two cases are to be distinguished:

• if the corresponding corner in the tableau is blue:

t s u s t u
and we erase the row of the corner in the tableau;

• if the corresponding corner in the tableau is red:

s u t s t u
and we erase the column of the corner in the tableau. On the example, we obtain: = It is not difficult to verify that this construction is a bijection.

For a permutation σ ∈ S n , its Up-Down sequence (cf. [START_REF] Viennot | Up-down sequences of permutations, paths and canopy of binary trees[END_REF]) is the vector Q(σ) = (q 1 , . . . , q n-1 ) ∈ {-1, +1} n-1 such that

q i = +1 iff σ i+1 > σ i .
It is clear that for any tree T , all the σ in Z T have the same Up-Down sequence, which we may call the Up-Down sequence of T , also called canopy of the binary tree T in [START_REF] Viennot | Up-down sequences of permutations, paths and canopy of binary trees[END_REF]. Now we may view the shape of a Catalan alternative tableau of size n -1 as a vector in {-1, +1} n-1 (horizontal steps correspond to "-1" entries and vertical steps to +1 entries).

We have the following property:

Proposition 2. The shape of the tableau associated to a tree T through the bijection described in Proposition 1 is the Up-Down sequence of T , as well as the common Up-Down sequence of any permutation σ in Z T .

Now the algorithm described above may be extended to labelled trees: we may put n labels on the shape of a Catalan alternative tableau of size n -1: If we keep the labels of the nodes when we apply the algorithm to get a tree from the tableau, we obtain a labelled tree In the previous example, we may say that the tableau was labelled with the permutation 792531648. As a consequence of the bijection, we get the following property. Proposition 3. Let σ be a permutation of size n. We may label a tableau C with σ, then apply the bijective algorithm. The labelled tree that we obtain is increasing if and only if:

• the shape of C is the Up-Down sequence of σ;

• the position of the red and blue dots in C is the only one which gives the binary tree Ψ(σ).

The product of trees through Catalan alternative tableaux

Now we come to the main result of this work. where C 1 and C 2 are the Catalan alternative tableaux associated respectively to T 1 and T 2 , and the question mark (?) represent any (valid) placement of (red and blue) dots in the rectangles.

Proof. By definition of Ψ * , we have:

(1) Ψ * (T 1 * T 2 ) = Ψ * (T 1 ) * Ψ * (T 2 ) = σ 1 ∈Z T 1 * σ 2 ∈Z T 2 σ 2 = σ∈S σ.
Let σ be an element of S. By definition of the product * in the Malvenuto-Reutenauer algebra, σ is of the form: σ = τ 1 .τ 2 (concatenation) with the letters appearing in τ 1 and τ 2 form a partition of {1, . . . , n}, and

(2) Ψ(τ 1 ) = T 1 and Ψ(τ 2 ) = T 2 .
Thus if Ψ(σ) = T ,the Up-Down sequence of T Q(T ) is either Q(T 1 )UpQ(T 2 ) or Q(T 1 )DownQ(T 2 ). Hence the form of the Catalan alternative tableau C associated to T is one of the two given in the Theorem 3. We label the shape of C with the entries of σ. The red and blue dots in C have to be placed in a position such that by applying the bijective algorithm, we obtain an increasing binary tree. But if we apply the algorithm to the part of C that carries the entries of τ 1 (respectively τ 2 ), Propositions 1 and 3 imply that the (red and blue) dots of C in the corresponding subparts of C have to be placed in the same configuration than in C 1 (respectively C 2 ). Thie implies that C has the required form.

Conversely, let T be a tableau of the form described in the Theorem 3, and σ ∈ Z T . By cutting σ in two parts u 1 and u 2 of lengths the sizes of C 1 and C 2 , we may write: σ = u 1 .u 2 with Std(u 1 ) = τ 1 and Std(u 2 ) = τ 2 .

It is again a simple application of Propositions 1 and 3 that we have: Ψ(τ 1 ) = T 1 and Ψ(τ 2 ) = T 2 , which was to be proved to complete the proof of Theorem 3.

The # product of binary trees

In light of Theorem 3, it seems natural to introduce a new product on k[Y ∞ ] as follows.

Definition 4. We define the # product of two binary trees T 1 and T 2 , associated respectively to Catalan alternative tableaux C 1 and C 2 by:

T 1 # T 2 = T
where the sum is taken over the trees T associated to Catalan alternative tableaux in the set:

C 1 ? C 2
and the question mark (?) represent any (valid) placement of (red and blue) dots in the rectangles.

It is clear that this defines an associative product on k[Y ∞ ]. It is worth to note that for T 1 ∈ Y k and T 2 ∈ Y l , then the product T 1 # T 2 is in Y k+l-1 (in this case the number of internal edges is preserved). We give below an example of this product, that should be checked by the reader.

# = + +

  For example, Std(3275) = 2143. Conversely, for σ ∈ S k a permutation and A = {a 1 , a 2 , . . . , a k } a set of k (distinct) integers, we define σ |A the k-tuple with distinct entries in A such that Std(σ |A ) = σ. With this notation we may define the product * in k[S ∞ ] as follows. Let σ ∈ S k and τ ∈ S l . We set σ * τ = A⊔B={1,2,...,k+l} σ |A .τ |B where ⊔ denotes the disjoint union, and . stands for concatenation. For example : 12 * 213 = 12 435+13 425+14 325+15 324+23 415+24 315+25 314+34 215+35 214+45 213.

Theorem 3 .

 3 Let T 1 and T 2 be two binary trees. Their product in the Loday-Ronco algebraT 1 * T 2 = Tis given by taking the sum over the trees T associated to Catalan alternative tableaux in the union

Acknowledgement. The authors sincerely thank F. Chapoton for fruitful discussions.

This research has been supported by the ANR (project MARS/06-BLAN-0193).


Abstract-IAC was initially introduced as a developmental mechanisms allowing a robot to self-organize developmental trajectories of increasing complexity without pre-programming the particular developmental stages. In this paper, we argue that IAC and other intrinsically motivated learning heuristics could be viewed as active learning algorithms that are particularly suited for learning forward models in unprepared sensorimotor spaces with large unlearnable subspaces. Then, we introduce a novel formulation of IAC, called R-IAC, and show that its performances as an intrinsically motivated active learning algorithm are far superior to IAC in a complex sensorimotor space where only a small subspace is neither unlearnable nor trivial. We also show results in which the learnt forward model is reused in a control scheme.

Index Terms-active learning, intrinsically motivated learning, developmental robotics, artificial curiosity, sensorimotor learning.

I. INTRINSICALLY MOTIVATED EXPLORATION AND LEARNING

Developmental robotics approaches are studying mechanisms that may allow a robot to continuously discover and learn new skills in unknown environments and in a lifelong time scale [START_REF] Weng | Autonomous mental development by robots and animals[END_REF], [START_REF] Lungarella | Developmental robotics: A survey[END_REF]. A main aspect is the fact that the set of these skills and their functions are at least partially unknown to the engineer who conceive the robot initially, and are also task-independent. Indeed, a desirable feature is that robots should be capable of exploring and developing various kinds of skills that they may re-use later on for tasks that they did not foresee. This is what happens in human children, and this is also why developmental robotics shall import concepts and mechanisms from human developmental psychology.

A. Learn from the Real Experimentations

Like children, the "freedom" that is given to developmental robots to learn an open set of skills also poses a very important problem: as soon as the set of motors and sensors is rich enough, the set of potential skills become extremely large and complicated. This means that on the one hand, it is impossible to try to learn all skills that may potentially be learnt because there is not enough time, and also that there are many skills or goals that the child/robot could imagine but never be actually learnable, because they are either too difficult or just not possible (for example, trying to learn to control the weather by producing gestures is hopeless). This kind of problem is not at all typical of the existing work in machine learning, where usually the "space" and the associated "skills" to be learnt and explored are well-prepared by a human engineer. For example, when learning hand-eye coordination in robots, the right input and output spaces (e.g. arm joint parameters and visual position of the hand) are typically provided as well as the fact that hand-eye coordination is an interesting skill to learn. But a developmental robot is not supposed to be provided with the right subspaces of its rich sensorimotor space and with their association with appropriate skills: it would for example have to discover that arm joint parameters and visual position of the hand are related in the context of a certain skill (which we call hand-eye coordination but which it has to conceptualize by itself) and in the middle of a complex flow of values in a richer set of sensations and actions.

B. Intrinsic motivations

Developmental robots have a sharp need for mechanisms that may drive and self-organize the exploration of new skills, as well as identify and organize useful sub-spaces in its complex sensorimotor experiences. In psychology terms, this amount to trying to answer the question "What is interesting for a curious brain?". Among the various trends of research which have approached this question, of particular interest is work on intrinsic motivation. Intrinsic motivations are mechanisms that guide curiosity-driven exploration, that were initially studied in psychology [START_REF] White | Motivation reconsidered: The concept of competence[END_REF]- [START_REF] Deci | Intrinsic Motivation and Self-Determination in Human Behavior[END_REF] and are now also being approached in neuroscience [START_REF] Schultz | Getting Formal with Dopamine and Reward[END_REF]-|8]. Machine learning researchers have proposed that such mechanism might be crucial for self-organizing developmental trajectories as well as for guiding the learning of general and reusable skills in machines and robots [START_REF] Oudeyer | Intrinsic Motivation Systems for Autonomous Mental Development[END_REF][START_REF] Barto | Intrinsically motivated learning of hierarchical collections of skills[END_REF]. Experiments have been conducted in real-world robotic setups, such as in [START_REF] Oudeyer | Intrinsic Motivation Systems for Autonomous Mental Development[END_REF] where an intrinsic motivation system was shown to allow for the progressive discovery of skills of increasing complexity, such as reaching, biting and simple vocal imitation with and AIBO robot. In these experiments, the focus was on the study of how developmental stages could self-organize into a developmental trajectory without a direct pre-specification of these stages and their number.

II. ROBUST INTELLIGENT ADAPTIVE CURIOSITY (RIAC) AS ACTIVE LEARNING

The present paper aims to propose a new version of the algorithm called Intelligent Adaptive Curiosity (IAC) presented in [START_REF] Barto | Intrinsically motivated learning of hierarchical collections of skills[END_REF], and to show that it can be used as an efficient active learning algorithm to learn forward models in a complex unprepared sensorimotor space. This algorithm, based on intrinsic motivations heuristics, implements an active and adaptive mechanism for monitoring and controlling the growth of complexity in exploration and incremental learning.

R-IAC : Robust Intrinsically Motivated Active Learning

Adrien Baranes and Pierre-Yves Oudeyer INRIA Bordeaux Sud-Ouest -FLOWERS team

In [START_REF] Oudeyer | Intrinsic Motivation Systems for Autonomous Mental Development[END_REF], it was presented focusing on its ability to generate organized developmental stages and trajectories within a cognitive modeling endeavour. Here, we rather take an engineering approach to study how IAC and a new formulation, called Robust-IAC (R-IAC), can efficiently drive the robot to learn fast and correctly a forward model.

A. Developmental Active Learning

An essential activity of epigenetic robots is to learn forward models of the world, which boils down to learning to predict the consequences of its actions in given contexts. This learning happens as the robot collects learning examples from its experiences. If the process of example collection is disconnected from the learning mechanism, this is called passive learning. In contrast, researchers in machine learning have proposed algorithms allowing the machine to choose and make experiments that maximize the expected information gain of the associated learning example [START_REF] Cohn | Active learning with statistical models[END_REF], which is called "active learning". This has been shown to dramatically decrease the number of required learning examples in order to reach a given performance in data mining experiments [START_REF] Hasenjager | Active Learning in Neural Networks[END_REF], which is essential for a robot since physical action costs time and energy. We argue that intrinsically motivated learning with IAC can be considered as an "active learning" algorithm. We will show that some of them allow very efficient learning in unprepared spaces with the typical properties of those encountered by developmental robots, outperforming standard active learning heuristics.

The typical active learning heuristics consist in focusing the exploration in zones where unpredictability or uncertainty of the current internal model are maximal, which involves the online learning of a meta-model that evaluates this unpredictability or uncertainty.

Unfortunately, it is not difficult to see that it will fail completely in unprepared robot sensorimotor spaces. Indeed, the spaces that epigenetic robots have to explore are typically composed of unlearnable subspaces, such as for example the relation between its joints values and the motion of unrelated objects that might be visually perceived. Classic active learning heuristics will push the robot to concentrate on these unlearnable zones, which is obviously undesirable.

Based on psychological theories proposing that exploration is focused on zones of optimal intermediate difficulty or novelty [START_REF] Berlyne | Conflict, Arousal, and Curiosity[END_REF], [START_REF] Csikszentmihalyi | Creativity-Flow and the Psychology of Discovery and Invention[END_REF], intrinsic motivation mechanisms have been proposed, pushing robots to focus on zones of maximal learning progress [START_REF] Oudeyer | Intrinsic Motivation Systems for Autonomous Mental Development[END_REF]. As exploration is here closely coupled with learning, this can be considered as active learning. We will now present the IAC system together with its novel formulation R-IAC. After this, we will evaluate their active learning performances in an inhomogeneous sensorimotor space with unlearnable subspaces.

B. Prediction Machine and Analysis of Error Rate

We consider a robot as a system with motor channels M and sensori channels S (M and S can be low-level such as torque motor values or touch sensor values, or higher level such as a "go forward one meter" motor command or "face detected" visual sensor"). Real valued action/motor parameters are represented as a vector 𝐌(𝐭), and sensors, as 𝐒(𝐭), at a time t. 𝐒𝐌(𝐭) represents a sensorimotor context, i.e. the concatenation of both motors and sensors vectors.

We also consider a Prediction Machine PM, as a system based on a learning algorithm (neural networks, KNN, etc.), which is able to create a forward model of a sensorimotor space based on learning examples collected through selfexperiments. Experiments are defined as series of actions, and consideration of sensations detected after actions are performed. An experiment is represented by the set (𝐒𝐌(𝐭), 𝐒(𝐭 + 𝟏)), and denotes the sensori consequence S(t+1) that is observed when actions encoded in M(t) are performed in the sensori context S(t). This set is called a "learning exemplar". After each trial, the prediction machine PM gets this data and incrementally updates the forward model that it is encoding, i.e. the robot incrementally increases its knowledge of the sensorimotor space. In this update process, PM is able to compare, for a given context 𝐒𝐌(𝒕), differences between predicted sensations 𝐒 (𝒕 + 𝟏) (estimated using the created model), and real consequences S(𝒕 + 𝟏). It is then able to produce a measure of error 𝒆(𝒕 + 𝟏), which represents the quality of the model for ensorimotor context 𝐒𝐌(𝒕).

Then, we consider a module able to analyze learning evolutions over time, called Prediction Analysis Machine PAM, Fig. 1. In a given subregion R of the sensorimotor space, this system monitors the evolution of errors in predictions made my PM by computing its derivative, i.e. the learning progress, LP=e N -e F in this particular region over a sliding time window (see Fig 1). LP is then used as a measure of interestingness used in the action selection scheme outlined below. The more a region is characterized by learning progress, the more it is interesting, and the more the system will perform experiments and collect learning examplars that fall into this region. Of course, as exploration goes on, the learnt forward model becomes better in this region and learning progress might decrease, leading to a decrease in the interestingness of this region. To precisely represent the learning behavior inside the whole sensorimotor space and differentiate its various evolutions in various subspaces/subregions, different PAM modules, each associated to a different subregion 𝑅 𝑖 of the sensorimotor space, need to be built. Therefore, the learning progress provided as the output values of each PAM (called Di the Fig. 2.) become representative of the interestingness of the associated region 𝑅 𝑖 . Initially, the whole space is considered as one single region 𝑅 0 , associated to one PAM, which will be progressively split into subregions with their own PAM as we will now describe.

C. The Split Machine

The Split Machine SpM possesses the capacity to memorize all the experimented learning exemplars (𝐒𝐌(𝐭), 𝐒(𝐭 + 𝟏)), and the corresponding errors values 𝒆(𝒕 + 𝟏). It is both responsible for identifying the region and PAM corresponding to a given SM(t), but also responsible of splitting (or creating in R-IAC where parent regions are kept in use) sub-regions from existing regions.

1) Region Implementation

We use a tree representation to store the list of regions as shown in Fig. 3. The main node represents the whole space, and leafs are subspaces. 𝐒(𝐭) and 𝐌(𝐭) are here normalized into [0;1] n . The main region (first node), called 𝑅 0 , represents the whole sensorimotor space. Each region stores all collected examplars that it covers. When a region contains more than a fixed number T split of exemplars, we split it into two ones.

When this criterion has been reached by a region, the split algorithm is executed, splitting just in one dimension at a time. An example of split execution is shown in Fig. 3, using a two dimensions input space.

3) R-IAC Split Algorithm

In R-IAC, the splitting mechanism is based on comparisons between the learning progress in the two potential child regions. The principal idea is to perform the separation which maximizes the dissimilarity of learning progress comparing the two created regions. This leads to the direct detection of areas where the learning progress is maximal, and to separate them from others (see Fig. 4). This contrasts with IAC where regions where built independently of the notion of learning progress. Reusing the notations of the previous section, in R-IAC the split of 𝜑 𝑛 into 𝜑 𝑛 +1 and 𝜑 𝑛+2 is done by choosing 𝑗 and 𝑣 𝑗 such that:

(𝐿𝑃 𝐞 𝒕 + 𝟏 | 𝐒𝐌 𝒕 , 𝐒 𝒕 + 𝟏 ∈ 𝜑 𝑛+1 -𝐿𝑃 𝐞 𝒕 + 𝟏 | 𝐒𝐌 𝒕 , 𝐒 𝒕 + 𝟏 ∈ 𝜑 𝑛 +2 ) 2 is maximal, where 𝐸 where 𝐸 is a set of errors values 𝑒 𝑖 with errors indexed by their relative order i of encounter (e.g. error e( 9) corresponds to a prediction made by the robot before another prediction which resulted in e [START_REF] Barto | Intrinsically motivated learning of hierarchical collections of skills[END_REF]: this implies that the order of examplars collected and associated prediction errors are stored in the system). 𝐸 is the cardinal of this set, and 𝐿𝑃 𝐸 is the learning progress responsible of the computation of learning progress estimations.

D. Action Selection Machine

We present here an implementation of Action Selection Machine ASM. The ASM decides of actions 𝐌 𝒕 to perform, given a sensori context 𝐒 𝒕 . (See Fig. We now present the different exploration modes used by the Action Selection Machine, in IAC and R-IAC algorithm:

1) Random Babbling Exploration Mode (1)

The random babbling mode corresponds to a totally random exploration, which does not consider previous actions and context. This mode appears in both IAC and R-IAC algorithm, with a probability 𝒑 𝟏 = 30%.

2) Learning Progress Maximization Exploration Mode (2)

The learning progress maximization mode aims to maximize the informational gain obtained after each experiment. To do this, it considers data computed by all PAM. The main idea is to consider that regions which have maximum learning progress are potentially the more interesting to explore. In the IAC algorithm, mode 2 action selection is straightforward: the leaf region which learning progress is maximal is found, and a random action within this region is chosen. In the R-IAC algorithm, we take into account the fact that many regions may have close learning progress values by taking a probabilistic approach. Furthermore, instead of focusing on the leaf regions like in IAC, R-IAC continues to monitor learning progress in node regions and select them if they have more learning progress. Let us give more details:

a) Probabilistic approach
The probabilistic mechanism is based on the consideration of a probability to choose a region proportional to the learning progress. We have, for a set of derivatives 𝐷 = 𝐷 0 , 𝐷 1 … 𝐷 𝑛 representing each region 𝑅 1 , 𝑅 2 … 𝑅 𝑛 using PAM, the probabilities 𝑷 𝒏 to choose the region 𝑅 𝑛 as :

𝑷 𝒏 = 𝐷 𝑛 -𝑚𝑎𝑥 𝐷 𝑖 𝐷 𝑖 -𝑚𝑎𝑥 𝐷 𝑖 𝒊

b) Multiresolution Monitoring of Learning Progress

In the IAC algorithm, the estimation of learning progress only happens in leaf regions. In R-IAC, learning progress is monitored in all regions created during the system's life time, which allows to track learning progress at multiple resolution in the sensorimotor space. This implies that when a new exemplar is available, R-IAC updates the evaluation of learning progress in all regions that cover this exemplar (but only if the exemplar was chosen randomly, i.e. not with mode 3 as described below).

In R-IAC mode 2, when a region has been chosen with the probabilistic approach and mutiresolution scheme, a random action is chosen within this region. Mode 2 is typically chosen with a probability 𝒑 𝟐 = 60% in both R-IAC and 𝒑 𝟐 = 70% in IAC (which means this is the dominant mode).

3) Error Maximization Exploration Mode (3)

Mode 3 combines a traditional active learning heuristics with the concept of learning progress: in mode 3, a region is first chosen with the same scheme as in R-IAC mode 2. But once this region has been chosen, an action in this region is selected such that the expected error in prediction will be maximal. This is currently implemented through a k-nearest neighbor regression of the function SM(t)  e(t+1) which allows to find the point of maximal error, to which is added random noise (to avoid to query several times exactly the same point). Like shown in Fig. 2, we store, for each region, coordinates of the generated point (called "MAX"). Mode 3 is typically chose with a probability 𝒑 𝟑 = 60% in R-IAC (and does not appear in IAC).

E. Consequences of Learning Progress Examination

The examination of learning progress allows a control of the learning complexity. Let us imagine three typical situations to illustrate precisely this phenomenon:  The system is exploring a simple area: The learning rate is high during a brief instant, and then, it is approximately null. The derivative is thus constant, and the probability to explore this kind of area is low.  The system is exploring a difficult area:The learning rate is varying very rapidly. The derivative is thus about zero, and the probability is low, like in the previous case.  The system is exploring a zone of mean difficulty: The learning rate is increasing, the derivative is thus negative, and the probability is depending on its absolute value. Observing these three examples, representing possible situations encountered, we argue that the learning progress is a guide toward areas of intermediate difficulty.

III. THE HAND-EYE-CLOUDS EXPERIMENT

We will now compare the performances of IAC and R-IAC as active learning algorithms to learn a forward model in a complex 6-dimensional sensorimotor space that includes large unlearnable zones. Both algorithms will also be compared with baseline random exploration.

In this experiment, a simulated robot has two 2-D arms with two joints controlled by motor inputs 𝒒 𝟏𝟏 , 𝒒 𝟏𝟐 , 𝒒 𝟐𝟏 , 𝒒 𝟐𝟐 . On the tip of one of the two arms is attached a square camera capable to detect the sensori position (𝒙, 𝒚) of point-blobs (relative to the square). These point-blobs can be either the tip of the other arm or clouds in the sky (see figure 5). This means that when the right arm is positioned such that the camera is over the clouds, which move randomly, the relation between motor configurations and perception is quasi-random. If on the contrary the arms are such that the camera is on top of the tip of the other arm, then there is an interesting sensorimotor relationship to learn. Formally, the system has the relation:

(𝒙, 𝒚) = 𝑬(𝒒 𝟏𝟏 , 𝒒 𝟏𝟐 , 𝒒 𝟐𝟏 , 𝒒 𝟐𝟐 )

where (𝒙, 𝒚) is computed as follows:

(1) Nothing has been detected : the camera has been placed over the white wall, (𝒙, 𝒚) = (-10, -10); (2) The hand appears inside the camera: The value of the relative position of the hand in the camera referential 𝑪 is taken. According to the camera size, the x and y values are in the interval [0; 6]; (3) The camera is looking at the window: Two random values playing the role of random clouds displacement are chosen for output. The interval of outputs corresponds to camera size. (4) The camera is looking at the window and sees both hand and cloud: we choose a random output value, like if just a cloud had been detected.

This setup can be thought to be similar to the problems encountered by infants discovering their body: they do not know initially that among the blobs moving in their field of view, some of them are part of their "self" and can be controlled, such as the hand, and some other are independent of the self and cannot be controlled (e.g. cars passing in the street or clouds in the sky).

Thus, in this sensorimotor space, the "interesting" potentially learnable subsuspace is next to a large unlearnable subspace (unlearnable), and also next to a large very simple subspace (when the camera is looking neither the clouds not the tip of the other arm). The primary challenge is thus to avoid the noisy area, and to detect others.

Results.

In these experiments, the parameters of IAC and R-IAC are Tsplit=250, the learning progress window is 50, 𝒑 𝟏 = 0.3, 𝒑 𝟐 = 0.6, 𝒑 𝟑 = 0.1. Experiments span a duration of 100000 actions. The learning algorithm that is used to learn the forward model is an incremental version of Gaussian Mixture Regression [START_REF] Calinon | On Learning, Representing and Generalizing a Task in a Humanoid Robot[END_REF].

A first study of what happens consists in monitoring the distance between the center of the eye (camera), and the hand (tip of the other arm). A small distance means that the eye is looking the hand, and a high, that it is focusing on clouds (noisy part) or on the white wall. Fig. 6 shows histograms of these distances. We firstly observe the behavior of the Random exploration algorithm. The curve shows that the system is, in majority, describing actions with a distance of 22, corresponding to the camera looking at clouds or at the white wall. Interestingly, the curve of the IAC algorithm is similar but slightly displaced towards shorter distance: this shows that IAC pushed the system to explore the "interesting" zone a little more. We finally observe that RIAC shows a large difference with both IAC and Random exploration: the system spends three times more time in a distance inferior to 8, i.e. exploring sensorimotor configurations in which the camera is looking at the other arm's tip. Thus, the difference between R-IAC and IAC is more important than the difference between IAC and Random exploration.

Then, we evaluated the quality of the forward model learnt using the three exploration algorithms. We considered this quality in two respects: 1) the capability of the model to predict the position of the hand in the camera given motor configurations for which the hand is within the field of view of the robot; 2) the capacity to use the forward model to control the arm: given a right arm configuration and a visual objective, we tested how far the forward model could be used to drive the hand to reach this visual objective.The first kind of evaluation was realized by first building a test database of 1000 random motor configurations for which the hand is within the field of view, and then using it for testing the learnt models built by each algorithm at various stages of their lifetime (the test consisted in predicting the position of the hand in the camera given joint confiurations). Results are reported on the right of figure 7. The second evaluation consisted in generating a set of (𝒙, 𝒚) 𝑪 , 𝒒 𝟐𝟏 , 𝒒 𝟐𝟐 values that are possible given the morphology of the robot, and then use the learnt forward models to try to move the left arm to reach the (𝒙, 𝒚) 𝑪 objectives corresponding to particular 𝒒 𝟐𝟏 , 𝒒 𝟐𝟐 values. The distance between the reached point and the objective point was each time measured, and results are reported in the left graph of figure 7.

Both curves on figure 7 confirm clearly the qualitative results of the previous figure: R-IAC outperforms significantly IAC, which is only slighlty better than random exploration. We have thus shown that R-IAC is much more efficient in such an example of complex inhomogeneous sensorimotor space. IV. CONCLUSION IAC was initially introduced as a developmental mechanism allowing a robot to self-organize developmental trajectories of increasing complexity without pre-programming the particular developmental stages. In this paper, we have argued that IAC and other intrinsically motivated learning heuristics could be viewed as active learning algorithms, and were based on heuristics that are more suited than traditional active learning algorithms for operation in unprepared sensorimotor spaces with large unlearnable subspaces. Then, we have introduced a novel formulation of IAC, called R-IAC, and shown that its performances as an intrinsically motivated active learning algorithm were far superior to IAC in a complex sensorimotor space where only a small subspace was interesting. We have also shown results in which the learnt forward model was reused in a control scheme.

Fig. 1 .

 1 Fig. 1. Internal mechanism of the Prediction Analysis Machine PAM associated to a given subregion R of the sensorimotor space. This module considers errors detected in prediction by the Prediction Machine PM, and returns a value representative of the learning progress in the region. Learning progress is the derivative of errors analyzed between a far and a near past in a fixed length sliding window.

2 ) 1 ) 2 ) 3 ) 2 S

 21232 IAC Split Algorithm In the IAC algorithm, the idea was to find a split such that the two sets of exemplars into the two subregions would minimize the sum of the variances of 𝐒 𝒕 + 𝟏 components of exemplars of each set, weighted by the number of exemplars of each set. The idea was to split in the middle of zones of maximal change in the function SM(t)  S(t+1). Mathematically, we consider 𝜑 𝑛 = 𝐒𝐌 𝒕 , 𝐒 𝒕 + 𝟏 𝒊 as the set of exemplars possessed by region 𝑅 𝑛 . Let us denote 𝑗 a cutting dimension and 𝑣 𝑗 , an associated cutting value. Then, the split of 𝜑 𝑛 into 𝜑 𝑛 +1 and 𝜑 𝑛+2 is done by choosing 𝑗 and 𝑣 𝑗 such that: (All the exemplars 𝐒𝐌(𝒕), 𝐒(𝒕 + 𝟏) 𝒊 of 𝜑 𝑛+1 have a 𝑗 𝑡ℎ component of their 𝐒𝐌 𝒕 smaller than 𝑣 𝑗 (All the exemplars 𝐒𝐌(𝒕), 𝐒(𝒕 + 𝟏) 𝒊 of 𝜑 𝑛+2 have a 𝑗 𝑡ℎ component of their 𝐒𝐌 𝒕 greater than 𝑣 𝑗 (The quantity : 𝜑 𝑛+1 . 𝜎 𝐒 𝒕 + 𝟏 | 𝐒𝐌 𝒕 , 𝐒 𝒕 + 𝟏 ∈ 𝜑 𝑛+1 + 𝜑 𝑛 +2 . 𝜎 𝐒 𝒕 + 𝟏 | 𝐒𝐌 𝒕 , 𝐒 𝒕 + 𝟏 ∈ 𝜑 𝑛+2 is minimal, where where S is a set of vectors, and S , its cardinal.

Fig. 3 .

 3 Fig. 3. The sensorimotor space is iteratively and recursively split into subspaces, called "regions". Each region 𝑅 𝑛 is responsible for monitoring the evolution of the error rate in the anticipation of consequences of the robot's actions, if the associated contexts are covered by this region.

Fig. 2 .

 2 Fig. 2. General architecture of IAC and R-IAC. The prediction Machine is used to create a forward model of the world, and measures the quality of its predictions (errors values). Then, a split machine cuts the sensorimotor space into different regions, whose quality of learning over time is examined by Prediction Analysis Machines. Then, an Action Selection system, is used to choose experiments to perform.

  2.). The ASM heuristics is based on a mixture of several modes, which differ between IAC and R-IAC. Both IAC and R-IAC algorithms share the same global loop in which modes are chosen probabilistically: Loop:  Action Selection Machine ASM: given S(t), execute an action 𝐌 𝒕 using the mode (𝒏) with probability 𝒑 𝒏 and based on data stored in the region tree;  Prediction Machine PM: Estimate the predicted consequence 𝑺 𝒕+𝟏 using the prediction machine PM ;  External Environment: Measure the real consequence 𝑺 𝒕+𝟏  Prediction Machine PM: Compute the error 𝒆 𝒕 + 𝟏 = 𝒂𝒃𝒔(𝑺 𝒕+𝟏 -𝑺 𝒕+𝟏 );  Update the prediction machine PM with 𝐒𝐌 𝐭 , 𝐒 𝐭 + 𝟏  Split Machine SpM: update the region tree with 𝐒𝐌 𝐭 , 𝐒 𝐭 + 𝟏 and 𝒆 𝒕 + 𝟏 ;  Prediction Analysis Machine PAM: update evaluation of learning progress in the regions that cover 𝐒𝐌 𝐭 , 𝐒 𝐭 + 𝟏 End Loop

Fig. 4 .

 4 Fig.4. Evolution of the sensorimotor regions over time. The whole space is progressively subdivided in such a way that the dissimilarity of each subregion in terms of learning progress is maximal.

Fig. 6 .

 6 Fig. 6. Histogram of distances repartitions between hand and eye, after 100000 experiments, comparing Random, IAC and R-IAC exploration methods.

Figure 7

 7 Figure 7 Left: evolution of performances in control based on the model learnt by Random exploration (blue line), IAC exploration (red line) and R-IAC exploration (black line). Right : evolution of the generalization capabilities in prediction of the learnt forward models with Random exploration (blue), IAC (red), and R-IAC (black)


Samer Ammoun
email: samer.ammoun@mines-paris.org

Fawzi Nashashibi
email: fawzi.nashashibi@inria.fr

Real time trajectory prediction for collision risk estimation between vehicles

Keywords: Trajectory prediction, risk assessment, cooperative driving, vehicle dynamic model

In this paper, we present our approach for collision risk estimation between vehicles. The vehicles are equipped with GPS receivers and communication devices. Our approach consists on using the knowledge given trough communication tool to predict the trajectories of the surrounding vehicles. Based on these trajectories, we identify the configurations of the collisions between vehicles. The risk is calculated using several indicators that are reflecting not only the possible collisions but also the dangerousness of these collisions. Our algorithm is tested on crossroads using scenarios involving real prototypes producing realistic scenarios.

I. INTRODUCTION

The reduction of the amount of accidents on the roads is still a priority in the transportation domain in general and particularly in the ITS 1 domain. In fact, less accidents means more secure, more efficient and stress less driving. Actually, vehicles are considered as the main source of danger on roads since they are generally the major actors in the accidents. For these reasons, our work is aiming to give the vehicle some intelligence in order to be able to reduce the risk of accidents. As a first step, the vehicle gets the capacity to sense its environments via its own sensors. In the second step, the vehicle will be able to estimate and to deal with the risk incurred during the navigation. In this paper, we are focusing on the vehicle-vehicle accident cases. The vehicle-static objects accidents are considered as trivial and implicit cases for this algorithm. For the case of vehiclevulnerable objects 2 , our method is still valid with a need to propose a model of displacement for the vulnerable objects.

Regarding the possible help that the new technologies may offer to the driver, we consider that the nature of the risk is very related to the degree of the driver's responsibility. In other words, it depends on the degree of automation of the vehicle. We identify three levels to deal with the risk corresponding to the degree of automation of the vehicle :

1) Driver assistance system : in this case, the driver has a full control of the vehicle. The vehicle's intelligence 1 Intelligent Transportation Systems 2 as pedestrian or bicycles consists in the elaboration of adapted recommendations to the driver in order to help him to deal with risky situations. Some active safety features are now taking place into the vehicles but the driver is still the unique master of its own vehicle. The role of the vehicle is to sense the environment and to analyze the information about the neighborhood (no matter if the information is processed and analyzed on board the vehicle or exchanged via communication). The vehicle will perform periodically a prediction of the future situation of its neighborhood and detect any future risky configuration. The complexity of the prediction is related to the fact that the prediction is computed based on the recent information about the environment. The prediction does not take into consideration the next behavior of the driver who has the total liberty to change the future trajectory of the vehicle using the available actuators (brake pedal, steering wheel, gear and accelerator pedal). Our presented work is situated at this level of risk estimation. 2) Fully automated driving with standalone planning : in this case, the driving is fully automated and no role is reserved for the driver. The role of the vehicle is to sense its environment, to perform a planning of the next trajectory and to supervise the right execution of this path. The goal of the motion planner is to reach a final position with some constraints: collision free, lane keeping, etc. Regarding that the planning is performed locally, each vehicle will continue to estimate any possible collision with the trajectories of the neighbor vehicles. 3) Cooperative planning: here the vehicles, within the same neighborhood, are generating cooperatively riskfree trajectories. The navigation is then 100 percent safe. Because some hazardous situations are still possible, the vehicles keep running a collision detection algorithm.

In this paper, we will present our approach for performing risk estimation related to collisions between non-automated vehicles. In the first section, we will present an overview of our approach. In the following sections, we will detail different parts of our algorithm. Before concluding, we will present some results of our algorithm experimented on real scenarios.

II. GLOBAL APPROACH FOR RISK ASSESSMENT

Our approach consists of reducing the risk of collision through four sequential actions : environment modeling, trajectory prediction, collision detection and risk management. these steps are illustrated on the block diagram of the figure 1. In this work, we suppose that we dispose of the current state of the vehicle environment. In fact, modeling the environment of navigation consists on proposing a representation of the surrounding environment. This model should offer not only the static geometry of the roads but also the positions and the attributes of the dynamic objects on (and next to) the roads. This representation is provided to the vehicle using (separately or combined) two different approaches regarding the source of information :

• Standalone approach : each vehicle is responsible of building its own environment map based on the output of its internal and perception sensors. The drawback of this method is that the knowledge of the vehicle is limited quantitatively by the range of its sensors and qualitatively by the precision and the richness of the delivered information. In this approach, we are relying completely on the detection and the classification algorithms to give a description of the environment (static or mobile objects, types of vehicles, pedestrians, etc.) • Cooperative approach : based on inter-vehicle communication, the vehicles are able to exchange information about its respective environment. Each vehicle is thus receiving an accurate and rich information from the surrounding vehicles. The drawback of this method is that it is depending on the number of vehicles equipped with communication devices. With this method, a nonequipped vehicle is totally invisible to other vehicle.

At least during the period of deployment of the communication devices in vehicles, it is important to combine the two methods to take advantages of both approaches.

Based on this local map and in order to estimate any possible risky configurations, we proceed by projecting the current environment state model into the future. This prediction step is considered as the most delicate step in our approach. In fact, we proceed by predicting the behavior of each vehicle accounting on the current inputs. Nonetheless, we do not attempt to predict in any case the future behavior of the driver. In fact, any change in the driver's behavior is going to potentially influence the inputs of the system and consequently the future trajectory of the vehicle.

Each environment state model (and naturally the predicted trajectories) is taking into account the errors and the imprecisions on the positions and on the attributes of the different entities in the scene. These errors are due to the errors on the sensors outputs and/or the errors of detection and classification algorithms.

Once predicted for the ego vehicle and for each neighboring vehicle, the future trajectories are utilized to compute the possible collision between each couple of trajectories. This approach considers that the risk of collision is independent between each couple of vehicles. Obviously, the problem is more complex since we should take into consideration the whole situation in order to estimate this risk. In our application, we made the choice of exploiting the collision detection not for path planning but only for collision prevention.

Based on geometric and dynamic models of the vehicle, the dangerousness of the collision is estimated using several risk indicators.

In order to avoid accidents, two actions are possible with different degrees of implication in the decision and therefore with differences from the responsibility and legacy points of view.

• The informative action : the ADAS system will inform the driver about the dangerousness of the situation and potentially propose some maneuvers to reduce the risk. This information is delivered through an appropriate HMI. In this paper, we do limit our intervention on this level. • The active behavior : in this system, the vehicle takes the initiative and operates on the actuators to realize the adequate to perform collision avoidance In the next sections, we will describe in detail each part of this algorithm.

III. TRAJECTORY PREDICTION

A trajectory is not only a list of geographic positions of the vehicle. It is a spatio-temporal representation of the displacement of the vehicle. In other terms, it is a dynamic representation of the path. The trajectory prediction is considered as a key step in our methodology. Taking into consideration that it is very difficult to predict the human behavior, we assume that the aim of our approach is to predict at the time t 0 the future displacement of the vehicle during the time slice [t 0 , t 0 + h] in case that the driver maintains the same driving profile. h is prediction period in seconds. Whatever the prediction algorithm is, the period between the incoming of new inputs on the actuators of the vehicle and the integration of these inputs into our prediction model will always be considered as a period of uncertainty for the risk estimator. For this reason, the faster is the algorithm of prediction the more reactive is the system and the shorter is this period of uncertainty. The quality of the prediction can be estimated by measuring the errors committed between the predicted and the real trajectories of the vehicle.

In the next sections, we will expose the most used methods to compute the future trajectories of vehicles.

A. Geometric approach

Based on the fact that the roads are designed with specific geometric models, the trajectory of the vehicle can be modeled with geometric features since it is constrained to follow the roads. Some research works propose to use a polynomial representation of the trajectories. First, third and fifth degree polynomials are the most used regarding the complexity of the speed and acceleration models [2]. Others are using one or more clotoids to model the trajectory taking into account the continuity and the derivability of the entire profile [START_REF] Khosla | Accurate estimation of forward path geometry using twoclothoid road model[END_REF]. The advantage of this shape-based trajectory generation is the computation facilities that assure a reduced time of estimation. It is also used to predict some trajectory before the beginning of the maneuver [START_REF]An analysis of the lane changing manoeuvre on roads : the contribution of inter-vehicle cooperation via communication[END_REF]. The main drawback is that this predictor is not based on any model of the vehicle. Moreover, the shape of the roads is not able to determine totally the behavior of the vehicle especially on some parts of the roads where there are many possibilities as on crossroads.

B. Dynamic approach

In this approach, the model tends to reproduce the dynamic behavior of the vehicle. Many models exists regarding the degree of complexity of the representation. Between the complicated models (like Ackeremann model) and the simplest one (as the bicycle model), the shapes of the equations are different regarding the inputs of the system. Different levels of inputs can be considered as the acceleration, speed, yaw angle, steering wheel rate, brake and acceleration pedal pressure [START_REF] Sayyed | Cooperative decentralized intersection collision avoidance using extended Kalman filtering[END_REF], gear state etc. In spite of the complexity of this method, the main advantage is that the future trajectory is computed based on a vehicle model. The position, speed and acceleration are direct outputs of the predictor.

We choose to use the bicycle dynamic model for its simplicity. We integrate the model using a Kalman linear filter with the positions, the speed and the acceleration as inputs.

The path prediction is performed using a linear Kalman filter. By using this recursive filter, we will be able to:

• Filter the noisy GPS position.

• Estimate the position of the vehicle between two outputs of the GPS receiver. • Estimate a horizon of displacement regarding the motion model of the vehicle. • Estimate the positioning errors: this filter determines the ellipse of uncertainty which is the ellipse where the vehicle had the maximum probability to exist. In the prediction step, the state vector X and the covariance matrix P are estimated in each iteration using the equations:

Xk/k-1 = A k-1 Xk-1/k-1 + B k-1 u k-1 ;
(1)

P X k/k-1 = A k-1 P X k-1/k-1 A t k-1 + R v k-1 (2)
In our case, X is expressed as follow

X =         x y V x V y a x a y         (3)
The prediction matrix A is :

        1 0 ∆T 0 0 0 0 1 0 ∆T 0 0 0 0 1 0 ∆T 0 0 0 0 1 0 ∆T 0 0 0 0 1 0 0 0 0 0 0 1         (4)
In our system, we don't have a control input u so no need for the transformation matrix B. P is the error covariance matrix and R v is the covariance of the process noise which is assumed to be a zero mean normal distribution.

The correction state is performed each time we received a new GPS position. The observation vector Z is equal to:

Z k =     x y V x V k     (5)
Z is linked to the state vector through the following equation:

Z k = C k X k (6)
Where

C k =     1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0     (7)
is the observation matrix.

The optimal Kalman gain is then calculated using this equation:

K k = P X k/k-1 C t k C k P X k/k-1 C t k + R n -1 (8)
And regarding the calculated gain, we update the state vector X and the error covariance matrix P :

Xk/k = Xk/k-1 + K k (Z k -C k Xk/k-1 ) (9)
P k/k = (I -K k C k ) × P k/k-1 (10)
I is the identity matrix. To be able to use this filter, we are assuming that the observation noise is a zero mean Gaussian white noise.

IV. POSITIONING UNCERTAINTY AND VEHICLE

BOUNDARIES

It is very common that a GPS based representation of the position of the vehicle is modeled using an elliptic shape due to the probabilistic uncertainty model. This shape is representing the probability of presence. Thus, the collision between vehicles is performed by calculating the intersection between the corresponding ellipses. Because the intersection calculation between ellipses is not analytically easy, we use a method commonly utilized in the 3D modeling domain. This method consists on modeling the ellipse by a series of aligned circles as shown in the figure 2. The detection of the collision between the vehicles is therefore performed by calculating the distances between the respective circles of each vehicle. The use of this representation is also useful for the determination of risk indicators as explained in next paragraph. A collision is detected if one of the circles of the vehicle is intersecting at least one of the other vehicles circles as illustrated on the figure 2. In case we detect a possible collision, the main parameter to study is the time that remains before the first impact between the vehicles. This time is important because it delivers the duration in witch the driver can act to avoid the collision or minimize its dangerousness. This parameter is called time to impact TTC. It is calculated as the duration between the current time and the instant of the first impact between the vehicle if the respective vehicles are keeping the current speed vectors. In our study, we compare this TTC to the duration of 2 seconds which is commonly recognize as the reaction time of the driver (1 sec) and of the vehicle (1 sec). However, this parameter does not give any indication about the confidence of this detection or the possible way to minimize the risk related to this accident.

For these reasons, we propose three more kinds of indicators. Those parameters should be analyzed together in order to constitute a clear view of the possible collision :

• The number of circles implicated in the collision : This parameter give us an idea about the dangerousness of the collision. In fact, a collision that involves a larger number of circles is considered more dangerous. • The duration of the predicted collision : this indicator gives an idea about the robustness of the estimation. The longer is this time the more sure is the collision.

• The configuration of the collision : this parameter corresponds to the position of the circles implicated in the detected collision. This parameter is useful to propose a safe method to minimize the risk and to customize the avoidance maneuver regarding the configuration during the predicted collision.

For the moment, the TTC is still the main risk indicator. We aims to propose a risk that combines at least the TTC, the number of circles and the time of the collision in one risk function. The configuration of the collision is only useful in the case that we propose maneuver to avoid the collision.

VI. EXPERIMENTATION

Many experiences were led with our system using the fleet of vehicles LaRA 3 . In this paper, we are presenting the tests performed on crossroads. The two vehicles LaRA1 and LaRA2 are cooperating via communication by exchanging their positions and their speed. The experience on real vehicles is conducted in the Velizy suburb of Paris. The two vehicles are equipped with similarly with a Trimble Ag332 GPS receivers and a Dlink 2100 Wifi access point. The frequency of the communication frames was harmonized with the GPS data frequency (10 Hz). In this experience, we are not using any high level protocol of communication. In fact we are using the libpcap library [5] to address directly the communication device. For more information about the communication procedure, please refer to [START_REF] Ammoun | Crossroads risk assessment using GPS and Inter-vehicle communications[END_REF] Each vehicle LaRA is receiving the positions and the speed of the other vehicle. It is performing a prediction of its own trajectory and of the trajectory of distant vehicle. The figure 3 presents the trajectories of both vehicles as predicted on LaRA2 for a temporal horizon of 10 sec. The blue (resp. the red) ellipses present the positions of LaRA1 (resp. LaRA2). The analysis of these graphs shows that during this test, the time to collision is getting smaller as the vehicles are approaching the crossroad which means that the collision was becoming more unavoidable. At the same time, the period of collision is increasing; this means that the certainty of the estimation of the collision is getting higher. The percentage of circles involved in the collision is also increasing which means that the collision will be more and more dangerous. During this test, we did not go intentionally to the real collision, one of the vehicles brakes to let the other vehicle pass. This behavior is illustrated on the graph by the reduction of the dangerousness and the probability of the collision. Figure 5 shows an example of human machine interface that can be offered to the driver. On the left, we can see the dynamic map of the environment with the position of the local vehicle LaRA1 as obtained from the local GPS and the position of the distant vehicle LaRA2 as received by communication. On the same map, we can also see the predicted trajectories of both vehicles as estimated in LaRA1. On the right side, the risk indicator (associated only to the TTC) is presented by a 10 levels graphical indicator. On top of the HMI, the video is taken from LaRA1 during the experience. This algorithm is running in real time on both vehicles. One of the main difficulties that we met was the communication signal loss in a no-line of sight configuration (typically on crossroads). This loss constrains the time to prevent the driver, reduces the change to anticipate the risk and make our algorithm more sensitive to the GPS error. Equipping the central part of the crossroads by communication relay could be the solution of this problem.

VII. CONCLUSION

In this paper, we have presented our approach for risk assessment between vehicles. Our approach is based on predicting the future trajectory of each vehicle. The collision is detected using a special representation of the vehicle form with regard to the GPS errors. The enhancement of the risk estimator will be performed by enhancing the quality of each of his sequential parts. One of the main axes of enhancements is to propose an accurate model of trajectory prediction by proposing an accurate dynamic model of the vehicle and respectively more complex equations for the Kalman filtering. In this paper, we have also presented many risk indicators and we are working currently on merging these indicators in an unique function which describes the dangerousness of the accidents. The studied scenario is representing one risky situation on crossroads. The next step will be to study more complex behavior on crossroads like left(resp. right) turn across path.

Figure 1 .

 1 Figure 1. A block diagram summarizing the risk assessment steps

Figure 2 .

 2 Figure 2. The collision detection is performed by detecting the collision between the circles around each vehicle

Figure 3 .

 3 Figure 3. The trajectory of LaRA1 (blue) and LaRA2(red) during the experience. The estimated trajectories is represented by a series of uncertainty ellipses

Figure 4 .

 4 Figure 4. From top to bottom: the time to collision, the duration of the collision and the percentage of implicated circles during the experience of roads crossing

Figure 5 .

 5 Figure 5. A prototype of HMI to present to the driver

ACKNOWLEDGMENT

  This work was conducted within the research project INTERSAFE-2. INTERSAFE-2 is part of the 7th Framework Programme, funded by the European Commission. The partners of INTERSAFE-2 thank the European Commission for supporting the work of this project.