François Digne
email: digne@u-picardie.fr

AND Ivan Marin
email: marin@math.jussieu.fr

Jean Michel
email: jmichel@math.jussieu.fr

¸ois Digne

The center of pure complex braid groups

Keywords:

Introduction

Let W denote a (finite) complex reflection group of rank n ≥ 1, that is a finite subgroup of GL n (C) generated by pseudo-reflections, that is elements of GL n (C) which fix an hyperplane. We assume W to be irreducible, in order to simplify statements. To such a group is associated an hyperplane arrangement A, made of the collection of the fixed hyperplanes associated to the pseudo-reflections in W . Let X = C n \ A denote the hyperplane complement. The groups P = π 1 (X) and B = π 1 (X/W ) are known as the pure braid group and braid group associated to W . There is a short exact sequence 1 → P → B → W → 1. In the case W is a real reflection group, the group B is an Artin-Tits group, with prototype the usual braid group B associated to W = S n . An extension of many results on Artin-Tits groups to the more general setup introduced here has been proposed in [START_REF] Broué | Complex reflection groups, braid groups, Hecke algebras[END_REF]. Several of them are still conjectural, such as the determination of the center of these groups. The goal of this note is to clarify the status of this question. Recall that in [START_REF] Broué | Complex reflection groups, braid groups, Hecke algebras[END_REF] are introduced (infinite order) elements β ∈ Z(B), π ∈ Z(P ), with β |Z(W )| = π such that the image of β in W generates Z(W ).

Our purpose is to summarize what can be stated on this topic, as follows:

Theorem 1.1. The center Z(B) is infinite cyclic, and generated by β.

Theorem 1.2. The center Z(P ) is infinite cyclic, and generated by π.

Theorem 1.3. There exists a short exact sequence

1 → Z(P ) → Z(B) → Z(W ) → 1
These three results were conjectured in [START_REF] Broué | Complex reflection groups, braid groups, Hecke algebras[END_REF], and proved there for the infinite series G(de, e, n), as well as for the (easy) case of groups of rank 2. Actually, for a given group W , theorem 1.2 implies by general arguments theorem 1.3 and theorem 1.1, as every element in Z(W ) can be lifted to an element in Z(B). The remaining cases were the exceptional groups G 24 , G 25 , G 26 , G 27 , G 29 , G 31 , G 32 , G 33 , G 34 in Shephard-Todd notation. As noted in [START_REF] Broué | Complex reflection groups, braid groups, Hecke algebras[END_REF], theorem 1.1 for the so-called "Shephard" groups G 25 , G 26 , G 32 is true because they have the same B as some Coxeter group (however, contrary to what is claimed in [START_REF] Broué | Complex reflection groups, braid groups, Hecke algebras[END_REF], this does not prove the other two results !). Note also that in [2, proposition 2.23] the bottom-right square in the diagram is not commutative and the bottom sequence is not exact. Theorem 1.1 was proved by D. Bessis in [START_REF] Bessis | Finite complex reflection arrangements are K(π, 1)[END_REF] for the well-generated groups G 24 , G 25 , G 26 , G 27 , G 29 , G 32 , G 33 , G 34 . The remaining case of G 31 can be obtained by a short argument due to D. Bessis (personal communication) that we reproduce here for the convenience of the reader (see section 3).

Once theorem 1.1 is known, theorem 1.3 reduces to theorem 1.2, and more precisely to the statement Z(P ) ⊂ Z(B). Actually we prove the following stronger theorem in the subsequent sections.

Theorem 1.4. If U is a finite index subgroup of B, then Z(U ) ⊂ Z(B).

Garside theory

A Garside monoid is a cancellative monoid which is generated by its atoms (that is elements which have no proper divisors), such that any two elements have a least common right-multiple and a greatest common left-divisor and such that there exist an element ∆ (a "Garside element") whose left-and right-divisors are the same and generate the monoid. Here a left-divides b, which will be denoted by a b, means that there exists c with b = ac and similarly for right-divisibility. We refer to [START_REF] Dehornoy | Groupes de Garside[END_REF] or [START_REF] Digne | Garside and locally Garside categories[END_REF] for the basic notions on Garside theory. We will use in particular the following property: if M is a Garside monoid, let α(x) denote the left-gcd of an element x ∈ M with ∆, then α(xy) = α(xα(y)) for any y ∈ M .

In this section, we consider a Garside monoid M which satisfies in addition the following properties:

(i) There is an additive length function on M such that the atoms have length 1. This implies that the atoms are precisely the elements of length 1 and that the only element of length 0 is 1. (ii) For every couple of atoms s = t, their right-lcm ∆ s,t is balanced (meaning that its left-and right-divisors are the same). (iii) For any couple of atoms s = t and any positive integer n, the left-gcd of ∆ s,t and s n is equal to s. Note that many of the monoids which have been studied for complex braid groups and for Artin-Tits groups satisfy the above conditions. In particular we have Proposition 2.1. Let M be one of the following monoids:

(M1) The classical monoid of positive elements in the Artin-Tits groups associated to finite Coxeter groups (see [START_REF] Dehornoy | Groupes de Garside[END_REF]). (M2) The "dual" monoid of [START_REF] Bessis | Finite complex reflection arrangements are K(π, 1)[END_REF]. (M3) The "parachute" monoid of [START_REF] Corran | A new Garside structure for the braid groups of type (e, e, r)[END_REF]. (M4) The "dual" monoids for Artin-Tits groups of type à and C of [START_REF] Digne | Présentations duales pour les groupes de tresses de type affine Ã[END_REF] and [START_REF] Digne | A Garside presentation for Artin-Tits groups of type Cn[END_REF]. (M5) The monoids f (h, m) (for h, m ≥ 1) presented by generators x 1 , . . . , x m and relations

x 1 x 2 . . . x m x 1 . . .

h terms = x 2 x 3 . . . x m x 1 . . . h terms = . . .
Then M is a Garside monoid and satisfies (i), (ii) and (iii).

Proof. In case (M1) since in the classical monoid we have ∆ s,t = sts . . . for some e ≥ 2 and these are the only decompositions of ∆ st , we get that ∆ st is balanced. Moreover the only divisors of s n are smaller powers of s, of which only s divides ∆ s,t .

In cases (M2) and (M4) all divisors of ∆ are balanced, in particular the right-lcm of two atoms. Moreover in any decomposition of ∆ into a product of atoms, all atoms are different, hence the same property holds for divisors of ∆. If s is an atom we have thus s 2 ∆, so that α(s 2 ) = s 2 . Since α(s 2 ) has length at most 2 and is different from s 2 it has to be of length 1, hence equal to s and by induction we get α(s n ) = α(sα(s n-1 )) = α(s 2 ) = s. This implies that the left-gcd of s n with ∆ s,t , which divides α(s n ), is equal to s.

In the monoid M for G(e, e, r) introduced in [START_REF] Corran | A new Garside structure for the braid groups of type (e, e, r)[END_REF], any pair of dictinct atoms can be embedded in a submonoid M ′ which is of type (M1) or (M2) (see [4, section 6.3]). This embedding maps right-lcms in M ′ to right-lcms in M and left-divisibility in M ′ is the restriction of left-divisibility in M (see [4, lemmas 5.1 and 5.3]).

This implies that for properties (ii) and (iii) we are either in the first or second situation, depending on the choice of the couple of atoms .

The monoids (M5) have been investigated by M. Picantin in his thesis [START_REF] Picantin | Petits Groupes Gaussiens[END_REF]. Property (i) is clear, as the presentation is homogeneous. It is readily checked that the left-gcd of two distinct atoms is the Garside element ∆ = x 1 x 2 . . . x m x 1 . . . (h terms), which proves (ii). Finally, (iii) follows from the fact that every divisor of ∆ but itself has a unique decomposition as a product of atoms.

Complex braid groups however may be the group of fractions of Garside monoids which do not fulfill our conditions. An example, pointed to us by M. Picantin, is given by the monoid presented by generators σ 1 , . . . , σ n-1 , τ 1 , . . . , τ n and relations

σ i σ i+1 σ i = σ i+1 σ i σ i+1 , σ i σ j = σ j σ i for |j -i| ≥ 2, σ i τ i τ i+1 = τ i τ i+1 σ i , σ i τ i = τ i+1 σ i , σ i τ j = τ j σ i for j ∈ {i, i + 1}.
It is a Garside monoid for the Artin group of type B n , for which the right lcm σ i τ i = τ i+1 σ i of σ i and τ i+1 is clearly not balanced, hence condition (ii) is not satisfied. This presentation dates back to [START_REF] Chow | On the algebraic braid group[END_REF], and the properties of the associated monoid are studied in [START_REF] Picantin | A (not so) wild Garside structure for the hyperoctaedral braid groups[END_REF].

As pointed out to us by E. Godelle, there exist Garside monoids satisfying (i) and (ii) but not (iii), such as M =< a, b | a 2 = b 2 >, which provides a counterexample to the next proposition. In [14, proposition 4.4.3] the next proposition is proved for a classical Artin-Tits monoid (see also [15, lemma 2.2]). Proposition 2.2. Let M be a Garside monoid satisfying (i), (ii) and (iii) above. Let r be an atom of M and b, z ∈ M be such that r j b = bz for some j ≥ 1; then z = t j for some atom t with rb = bt.

Proof of the proposition.

Note that by properties (i) and (iii), given two atoms r = t, the element δ r,t defined by ∆ r,t = rδ r,t satisfies δ r,t = 1 and r δ r,t .

We first prove the following lemma.

Lemma 2.3. Let r be an atom in M and b ∈ M such that r b and let j ≥ 1 be such that r j b is divisible by some atom different from r; then there exists an atom s = r such that δ r,s b.

Proof. The proof is by induction on j ≥ 1. Let u = r be an atom with u r j b. Then ∆ r,u r j b; writing ∆ r,u = rδ r,u , we get δ r,u r j-1 b. If j = 1 the conclusion holds with s = u, and this covers the case j = 1.

In the case j > 1, since δ r,u = 1 and r δ r,u , we have v δ r,u for some atom v = r hence v r j-1 b and we can apply the induction assumption to j -1 < j and get the result.

We also need the following lemma. As ∆ r,s is balanced, it conjugates r to some atom u, by lemma 2.4, hence, cancelling a power of r, we get rδ r,s = δ r,s u, whence δ r,s u j b ′ = δ r,s b ′ z. Cancelling δ r,s and applying the induction assumption we get z = t j for some atom

t with ub ′ = b ′ t, so rb = rδ r,s b ′ = δ r,s ub ′ = δ r,s b ′ t = bt.
We first prove a general corollary for a Garside group, that is for the group of fractions G of a Garside monoid M . It is a property of Garside monoids that G is generated as a monoid by M and the powers of ∆.

Corollary 2.5. Let G be the group of fractions of a Garside monoid M satisfying properties (i), (ii) and (iii) then if U is a subgroup of G containing a nontrivial power of each atom, the center of U is contained in the center of G.

Proof. If x ∈ Z(U ), then for any atom s, there exists a positive integer N such that x commutes with s N . Let i be large enough so that x∆ i ∈ M . Since conjugation by ∆ permutes the atoms by lemma 2.4, we have s N x∆ i = xs N ∆ i = x∆ i u N for some atom u such that s∆ i = ∆ i u. Proposition 2.2 applied with b = x∆ i , r = s, j = N implies that sx∆ i = x∆ i t, for some atom t, hence we have t N = u N , which gives t = u by property (iii), since the left-gcd of t N = u N with ∆ t,u has to be equal to t and to u. Hence sx∆ i = x∆ i u and sx = xs.

A particular case of the above corollary is the following, using the fact that for monoids of type (M4) in proposition 2.1, the squares of the atoms lie in the pure group.

Corollary 2.6. The center of pure Artin-Tits groups of type à or C is trivial.

Another application of proposition 2.2 gives theorem 1.4 for all finite complex reflection groups except G 31 .

Corollary 2.7. Assume W is a complex reflection group different from G 31 and let U be a finite index subgroup of B. Then Z(U ) ⊂ Z(B).

Proof. It is well-known that all possible braid groups B can be obtained from a 2-reflection group, so we can restrict to this case. Corollary 2.5 gives the result for all complex braid group for which a monoid satisfying (i), (ii) and (iii) is known. This covers all of them except for a complex reflection group of type G 31 and the infinite series G(de, e, r) for d > 1 and e > 1. Indeed, one can use

• the classical monoid for Coxeter groups and Shephard groups (groups which have same X/W as a Coxeter group), as well as for G 13 , which has the same braid group as the Coxeter group I 2 (6) • the parachute monoid for the G(e, e, r)

• the dual monoid for the G(e, e, r), G(d, 1, r) and the exceptional groups of rank at least 3 which are not G 31 • the monoids f (4, 3) and f (5, 3) for G 12 and G 22 . Since the only 2-reflection exceptional groups are G 12 , G 13 and G 22 this indeed covers everything but G 31 and the G(de, e, r), for d > 1, e > 1.

But the braid group associated with G(de, e, r) is a subgroup of finite index in the braid group associated to G(de, 1, r) (and thus of the classical braid group on r strands), whence the result in this case.

Springer theory and G 31

We let B n denote the braid group associated to G n , and P n = Ker(B n → G n ) the corresponding pure braid group. In particular B 37 denotes the Artin-Tits group of type E 8 . By Springer theory (see [START_REF] Springer | Regular elements of finite reflection groups[END_REF]), G 31 appears as the centralizer of a regular element c of order 4 in G 37 , as a consequence of [1, thm. 12. For the convenience of the reader we reproduce it here. We let π denote the positive generator of Z(P 37 ). We have c4 = π. On the other hand, it can be checked that G 37 has a regular element of order 24, and as a consequence of [START_REF] Broué | Sur certains éléments réguliers des groupes de Weyl et les variétés de Deligne-Lusztig associées[END_REF] there exists d ∈ B 37 with d24 = π whose image in G 37 is a regular element d of order 24. Moreover, by [1, 12.5 (ii)], d6 is conjugated to c, so up to conjugating d we can assume d6 = c, and in particular d ∈ B 31 . So Z(B 31 ) lies inside the centralizer of d, which is by another application of [1, 12.5 (iii)] the braid group of the centralizer in G 37 of d. Since 24 divides exactly one reflection degree of G 37 , Springer's theory says that d generates its centralizer, as its centralizer is a reflection group whose single reflection degree is 24. This implies that the centralizer of d is the cyclic group generated by d, thus any x ∈ Z(B 31 ) is a power da of d. Since its image in

G 37 should lie in Z(G 31 ) =< d 6 >, Z(B 31 ) =< d6 > is a cyclic group, isomorphic
to Z as it is infinite, for instance because β 31 ∈ Z(B 31 ) has infinite order. But 4 , hence d6 = β 31 and Z(B 31 ) =< β 31 >≃ Z. This concludes the proof of theorem 1.1.

β 31 ∈ Z(B 31 ) =< d6 > satisfies β 4 31 = π 31 = π 37 = ( d6 )
In order to prove theorem 1.4, we will need an explicit description of this embedding B 31 ֒→ B 37 . A presentation for B 37 and B 31 is given by the Coxeter-like where the circle means x 1 x 2 x 3 = x 2 x 3 x 1 = x 3 x 1 x 2 . This latter presentation was conjectured in [START_REF] Broué | Complex reflection groups, braid groups, Hecke algebras[END_REF] and proved in [START_REF] Bessis | Finite complex reflection arrangements are K(π, 1)[END_REF].

We choose for regular element c = (s 4 s 2 s 3 s 1 s 4 s 3 s 5 s 6 s 7 s 8 ) 6 , and we use the algorithms of [START_REF] Franco | Computation of Centralizers in Braid groups and Garside groups[END_REF] included in the GAP3 package CHEVIE (see [23]) to find generators for B 31 , considered as the centralizer of c in B 37 . From this we get a description of the embedding B 31 ֒→ B 37 as follows:

x 1 → (s 2 s 3 s 1 s 5 ) -1 s 1 s 4 (s 2 s 3 s 1 s 5 )
x 2 → (s 4 s 2 s 3 s 5 s 6 s 5 s 7 ) -1 s 2 s 5 (s 4 s 2 s 3 s 5 s 6 s 5 s 7 )

x 3 → (s 5 s 6 s 7 ) -1 s 1 s 4 (s 5 s 6 s 7 )

x 4 → (s 2 s 5 s 6 ) -1 s 4 s 6 (s 2 s 5 s 6 )

x 5 → (s 3 s 1 s 5 s 6 ) -1 s 4 s 8 (s 3 s 1 s 5 s 6 ) One can check that β 31 = (x 4 x 1 x 2 x 3 x 5 ) 6 → c.

Krammer representations and G 31

We use the generalized Krammer representation R : B 37 ֒→ GL 120 (Q[q, q -1 , t, t -1 ]) for the Artin-Tits group of type E 8 defined in [START_REF] Digne | On the linearity of Artin braid groups[END_REF] and [START_REF] Cohen | Linearity of Artin groups of finite type[END_REF]. We use the definitions of [START_REF] Digne | On the linearity of Artin braid groups[END_REF] (note however the erratum given in [START_REF] Marin | Sur les représentations de Krammer génériques[END_REF]). One has β 37 = β 2 31 and R(β 37 ) = q 30 t. The standard generators s i of B 37 are mapped to semisimple matrices with eigenvalues q 2 t (once), -q (28 times) and 1 (91 times). We embed Q[q ±1 , t ±1 ] into Q[q ±1 , u ±1 ] under t → u 2 . Then R(β 31 ) can be diagonalized, has two 60-dimensional eigenspaces, corresponding to the eigenvalues ±q 15 u. An explicit base change, that we choose in GL 120 (Q[q ±1 , u ±1 ) so that it specializes to the identity when q → 1, u → 1, provides another faithful representation R : B 37 → GL 120 ([q ±1 , u ±1 ]) with R(B 31 ) ⊂ GL 60 (Q[q ±1 , u ±1 ]) × GL 60 (Q[q ±1 , u ±1 ]). Let U be a finite index subgroup of B 31 , and let N ≥ 1 such that x 2N i = (x 2 i ) N ∈ U for all i. Note that, since the images of the x i 's in G 31 have order 2, we have x 2 i ∈ P 31 . We prove that the centralizer of R(U ) in GL 60 ×GL 60 is equal to the centralizer of R(B 31 ). This is the case as soon as R j (KU ) = Mat 60 (K), where j ∈ {1, 2}, where K is an arbitrary extension of Q(q, u), R = R 1 × R 2 is the obvious decomposition, and KU is the group algebra of U . Actually, it is clearly enough to prove this for a specialization of R j to given values of t, u. More precisely, if A is a unital ring with field of fractions K and if we have a morphism Q[q ±1 , u ±1 ] → A, letting R ′ j : B 31 → GL 60 (A) denote the induced representation we have that R ′ j (KU ) = Mat 60 (K) implies that R j (Q(q, u)U ) = Mat 60 (Q(q, u)).

Let j ∈ {1, 2}. We use the morphism

Q[q ±1 , u ±1 ] → Q[[h]
] which maps q → e h , u → e 7h (the choice of 7 being rather random), and denote R ′ j : P 31 → Mat 60 (Q[[h]]) the induced representation. We prove that the unital algebra generated over Q((h)) by the R ′ j (x 2N i ) full Mat 60 (Q((h))). As x 2 i ∈ P 31 ⊂ P 37 we have

R ′ j (x 2 i ) ≡ 1 modulo h (see e.g. [16]) hence R ′ j (x 2 i ) N ≡ 1 modulo h. We define y i ∈ Mat 60 (Q) by h -1 (R ′ j (x 2 i ) -1) ≡ y i modulo h and thus h -1 (R ′ j (x 2 i ) N -1)
≡ N y i . By Nakayama's lemma it is now sufficient to prove that the y i generate Mat 60 (Q). It turns out that the y i belong to Mat 60 (Z). By another application of Nakayama's lemma it is sufficient to check that the reduction mod p of the y i generates Mat 60 (F p ) for some prime p. For a given p, the determination by computer of the dimension of the subalgebra of Mat 60 (F p ) generated by elements y 1 , . . . , y 5 is easy : starting from the line F 0 = F p Id ⊂ Mat 60 (F p ) we compute the vector space F r+1 = F r + F r y 1 + . . . F r y 5 until dim F r = dim F r+1 . For p = 37 we get 3600 = dim Mat 60 (F p ), thus concluding the argument.

By faithfulness of R this proves that Z(U ) commutes with all B 31 , thus Z(U ) ⊂ Z(B 31 ), and this concludes the proof of theorem 1.4.

Miscellaneous consequences

We notice that the representations R 1 and R 2 are deduced one from the other under the field automorphism of Q(q, u) defined by q → q, u → -u. Since R is faithful, this implies that the representations R j are faithful representations of B 31 . We note that the R j have dimension 60, which is the number of reflections in G 31 , and thus the dimension of the representation involved in conjecture 1 of [START_REF] Marin | Krammer representations for complex braid groups[END_REF].

This enables us to prove this conjecture for G 31 , that is the following theorem:

Theorem 5.1. B 31 can be embedded into GL 60 (K) as a Zariski-dense subgroup, for K a field of characteristic 0, and P 31 is residually torsion-free nilpotent.

Proof. Since P 31 embeds in P 37 , the statement about residual torsion-free nilpotence is a consequence of the corresponding statement for P 37 , proved in [START_REF] Marin | Sur les représentations de Krammer génériques[END_REF], [START_REF] Marin | On the residual nilpotence of pure Artin group[END_REF]. Embedding Q(q, u) into Q((h)) through q → e h , u → e αh , for α a transcendent number, we can assume R j (B 31 ) ⊂ GL 60 (Q[α]((h))).

In addition we know that R j (x 2 i ) = exp(ha i ) for some a i ∈ gl 60 (Q[α][[h]]) with a i ≡ ỹi mod h for some ỹi ∈ gl 60 (Q[α]) which specialize to y i under α → 7. A computer calculation similar to the above shows that the Lie algebra generated by the y i is gl 60 (Q), and thus that the Lie algebra generated by the a i over C((h)) is gl 60 (C((h))). By [START_REF] Marin | L'algèbre de Lie des transpositions[END_REF] lemme 21 this proves that the Zariski closure of R j (B 31 ) contains GL 60 (C((h))), which proves the theorem.

We remark that this conjecture, if true for all braid groups, would immediately imply theorem 1.4. Note however that the representations R j have same dimension as, but are not isomorphic to the representation constructed in [START_REF] Marin | Krammer representations for complex braid groups[END_REF]. This can be seen from the eigenvalues of the generators : since R(x 1 ) is conjugated to R(s 1 s 4 ) and since s 1 s 4 = s 4 s 1 , the 3 eigenvalues of the R(s i ) provide at most 9 eigenvalues for R(x 1 ), and by checking each of them we find that R(x 1 ) has for eigenvalues q 2 t (twice), q 2 (6 times), -q (44 times) and 1 (68 times). As a consequence, the R j (x 1 ) both have 4 eigenvalues, to be compared with the 3 eigenvalues the generators have in the construction of [START_REF] Marin | Krammer representations for complex braid groups[END_REF].

e=

  tst . . .

  e

Lemma 2 . 4 .

 24 Let b ∈ M be a balanced element. Then, for any atom r b there exists an atom t b with rb = bt, and r → t provides a permutation of the atoms dividing b. Proof. For r b we write b = rb ′ . Since b is balanced we have b ′ b hence b = b ′ t for some t ∈ M (which obviously divides b). Thus rb = rb ′ t = bt and r → t is well-defined and clearly injective from the set of divisors of b to itself. Because of the length function it maps atoms to atoms, and it is surjective by the obvious reverse construction. We show now the proposition by induction on the length of b, the case of length 0 being trivial. If r b then, by simplifying b by r, we get the result by induction. We thus can assume b = 1 and r b. By lemma 2.3 there is an atom s = r such that δ r,s b. Writing b = δ r,s b ′ we get r j δ r,s b ′ = δ r,s b ′ z.

  5 (iii)], B 31 can be identified with the centralizer of a lift c ∈ B 37 of c, in such a way that the natural diagram B 31 / / B 37 G 31 / / G 37 commutes. The proof of theorem 1.1 for G 31 was communicated to us by D. Bessis.


Oscar E Ramos

Layale Saab
email: lsaab@laas.fr

Sovannara Hak
email: shak@laas.fr

Nicolas Mansard
email: nmansard@laas.fr

Dynamic Motion Capture and Edition using a Stack of Tasks

This paper presents a complete methodology to quickly reshape a dynamic motion demonstrated by a human and to adapt the dynamics of the human to the dynamics of the robot. The method uses an inverse dynamics control scheme with a quadratic programming optimization solver. The motion data recorded using a motion capture system is introduced into the control scheme as a reference posture task to be followed by the joints trajectory respecting the dynamic limitations as well as the contact constraints. The motion is further modified using arbitrary tasks to let the robot imitate the original motion more closely or to make voluntary changes for aesthetic reasons. The results show the method applied to the humanoid robot HRP-2 imitating a human "pop dance".

I. INTRODUCTION

Human motion imitation by a humanoid robot is a challenging task involving coordination, control and stabilization of the robot. Coordination problems are inherent to kinematically redundant robots, and humanoids are known to be highly redundant systems. Control difficulties arise from the complex humanoid tree-like structure as well as its unstable nature due to its vertical position. The robot must not only reproduce some task, but it should not fall while performing the task, keeping its dynamic balance. These constraints make imitation of captured motion a more complicated problem in robotics than in computer animation.

A classical approach for motion imitation in computer graphics is inverse kinematics [START_REF] Gleicher | Retargetting motion to new characters[END_REF], [START_REF] Grochow | Style-based inverse kinematics[END_REF] since the characters usually do not possess physical dynamic characteristics. The main tools are related to motion editing and retargetting captured motions [START_REF] Yamane | Natural motion animation through constraining and deconstraining at will[END_REF]. Even though sometimes dynamics is taken into account [START_REF] Popovic | Editing dynamic properties of captured human motion[END_REF], its main objective is to give more realism to the animation rather than to make the motion feasible. The motion generated with methods for computer animation is not guaranteed to be suitable for a direct reproduction on a robot.

In robotics, capture and analysis of motion has become an active research area during these recent years to generate motion by imitation [START_REF] Suleiman | On human motion imitation by humanoid robot[END_REF]. This imitation paradigm that lets the robot acquire the way of performing a task by observing human demonstrations has been called learning from observation (LFO) [START_REF] Nakaoka | Learning from observation paradigm: Leg task models for enabling a biped humanoid robot to imitate human dances[END_REF] and was initially introduced for robotic manipulators with the name assembly plan from observation (APO) [START_REF] Ikeuchi | Toward an assembly plan from observation. Task recognition with polyhedral objects[END_REF]. The starting point is usually the motion acquired from a human expert using a motion capture system. Chinese Kungfu is an example of this type of imitation [START_REF] Zhao | Kinematics mapping and similarity evaluation of humanoid motion based on human motion capture[END_REF]. Also, Japanese traditional dancing has been replayed on a humanoid robot imitating the motion of a professional dancer [START_REF] Nakaoka | Generating whole body motions for a biped humanoid robot from captured human dances[END_REF], [START_REF] Nakaoka | Leg motion primitives for a dancing humanoid robot[END_REF], Fig. 1: Inverse Dynamics Cascade Scheme [START_REF] Nakaoka | Learning from observation paradigm: Leg task models for enabling a biped humanoid robot to imitate human dances[END_REF]. Moreover, modifications to the upper body motion have been realized for the Japanese dancing [START_REF] Shiratori | Temporal scaling of upper body motion for sound feedback system of a dancing humanoid robot[END_REF]. Another motion sequence resembling a dance was presented in [START_REF] Pollard | Adapting human motion for the control of a humanoid robot[END_REF] adapting the captured joint trajectories to generate feasible motion in a robot, and robot modeling is recreated in [START_REF] Miura | Robot motion remix based on motion capture data towards human-like locomotion of humanoid robots[END_REF] where walking and turning motion are created based on the motion obtained from a professional model.

Optimization is a classical solution for reshaping the captured motion before imitation. For example, upper body motion is produced solving a separate optimization problem for each frame of the motion in [START_REF] Safonova | Optimizing human motion for the control of a humanoid robot[END_REF]. A constraint that uses trajectory optimization and filtering to preserve the main characteristics of the original motion and to respect at the same time the physical limitations of the humanoid robot is proposed in [START_REF] Ruchanurucks | Humanoid robot motion generation with sequential physical constraints[END_REF]. A similar approach is proposed in [START_REF] Suleiman | On human motion imitation by humanoid robot[END_REF] but the optimization problem is reformulated using Lie algebra for the dynamic equations of motion. In this case, the ZMP is controlled using the cart table model [START_REF] Kajita | Biped walking pattern generation by using preview control of zero-moment point[END_REF].

In this paper, we propose to use a generic hierarchical optimization solver to solve simultaneously the dynamic reshaping and the motion edition. The considered solver consists of the inverse-dynamics control cascade proposed initially in [START_REF] Saab | Generation of Dynamic Motion for Anthropomorphic Systems under Prioritized Equality and Inequality Constraints[END_REF] and extended in [START_REF] Saab | Generic Dynamic Motion Generation with Multiple Unilateral Constraints[END_REF]. The flexibility of the scheme allows the addition of arbitrary tasks on the operational space that modify the trajectory of the joints, to generate a more similar motion, or to change some part of it, as Figure 1 shows. This paper is organized as follows. Section II presents the dynamic model considering the contact constraints. In Section III we recall the inverse dynamics solver applied for the generation of motion by resolving a set of tasks. The tasks used for motion imitation are introduced in Section IV. Section V presents the method used to reshape the observed human motion to the robot kinematics using an optimization algorithm. Finally, the results obtained in simulation for the humanoid robot HRP-2 imitating a dance motion are presented in Section VI.

II. DYNAMIC MODEL

A. Dynamic model with contact conditions

In this work, we consider a humanoid robot in contact with its environment. The full representation of the humanoid robot is given by the generalized coordinates q = (x b , q a ) with x b as the position and orientation of the robot with respect to the world, and q a as the vector including the n actuated joints. The arbitrary center of the body in contact (e.g. the body center of mass or the center of the joint attached to it) will be represented by x c and the frame attached to it by {C}. By abuse, x c will denotes the 6D position (translation and orientation) of the frame {C}. The contact forces are then expressed in terms of the generalized 6D contact force φ c = (f c , τ c ) written at {C}. Considering a single body in contact, the dynamic model is:

Aq + b + J T c φ c = S T τ (1)
where J c is the contact Jacobian of x c expressed at the same point as φ c , A is the whole-body inertia matrix, b is the vector including Coriolis, centrifugal and gravity forces, τ is the actuated torque vector, and S is a matrix that selects the actuated joints. This can be generalized for more than two bodies in contact as shown later.

An additional non-holonomic restriction to be satisfied when two rigid bodies are in contact is that no relative motion should occur between them [START_REF] Yamane | Dynamics filter-concept and implementation of online motion generator for human figures[END_REF]: ẋc = 0 and ẍc = 0 1 . At the acceleration level, this is ẍc = J c q + Jc q = 0, which implies

J c q = -Jc q (2)
B. Force constraint for rigid planar contacts 1) Unilateral model: In the following, the only contact model considered is the rigid multi-planar contact. For each contact, at least three non aligned points must lie on the planar surface. Consider only one body in rigid contact and let the points belonging to the contact surface be p i , i = 1, ..., n p with n p ≥ 3. At each vertex p i , the force is represented by

f i = (f ix , f iy , f iz ). The torque expressed at {C} is c τ i = c p i × f i = c pi f i ,
where c pi is the skew-matrix obtained from c p i , and c p i is the point p i with respect to the reference frame {C}. The total force at the center of contact is then

f c = np i=1 f i and the total torque is τ c = np i=1 c pi f i : φ c = f c τ c = I I ... I c p1 c p2 ... c p np    f 1 . . . f np    . (3)
1 We consider rigid contact (ẍ = 0, contact of class 6) which implies a minimum of three contact points. For other classes of contact (e.g. punctual contact) or non-rigid contact (i.e. ẍ is obtained as a non-zero function of the exerted forces), the same structure can be used, but it is not considered here.

A sufficient condition to guarantee a rigid contact between a body of the robot and the planar surface is that the reaction force at the interface must be directed towards the robot; equivalently, the normal components f iz of each contact force must be positive [START_REF] Saab | Generation of Dynamic Motion for Anthropomorphic Systems under Prioritized Equality and Inequality Constraints[END_REF]:

f ⊥ = (f 1z , f 2z , . . . , f npz ) ≥ 0 (4)
2) Constraint reduction: The contact can be represented by [START_REF] Yamane | Natural motion animation through constraining and deconstraining at will[END_REF] and inequality [START_REF] Popovic | Editing dynamic properties of captured human motion[END_REF]. Using the planar hypothesis, it is possible to reduce the size of these conditions. Consider a reference frame {O} as the projection of {C} on the contact surface so that the Z axis of {O} is normal to the surface, and let o c = (0, 0, o c z ) be the origin of frame {C} expressed with respect to {O}. The total force and torque with respect to the new frame are given by

φ o = f o τ o = I 0 o ĉ I f c τ c . (5)
Each contact point in the reference frame {O} is o p i = ( o p ix , o p iy , 0), with o p iz = 0 as they also lie in the contact surface. Using (5) and the relation

τ o = np i=1 o pi f i , it is straightforward to verify that:         f cx f cy f cz -o c z f cy + τ cx o c z f cx + τ cy τ cz         =         np i f ix np i f iy np i f iz np i o p iy f iz - np i o p ix f iz - np i o p iy f ix + i o p ix f iy         (6)
Since the rigid contact condition on the force only constraints the normal component f iz , the first, second and sixth rows of (6) impose no restriction to the dynamic evolution and can be arbitrarily chosen. With these simplifications, the rigidcontact constraint including ( 4) and ( 6) is given by [START_REF] Saab | Generic Dynamic Motion Generation with Multiple Unilateral Constraints[END_REF]:

Q c1 Q c2 φ c f ⊥ = 0 and f ⊥ ≥ 0 (7)
with

Q c1 =   0 0 -1 0 0 0 0 o c z 0 -1 0 0 -o c z 0 0 0 -1 0   (8)
Q c2 =   1 1 . . . 1 o p 1y o p 2y . . . o p np y -o p 1x -o p 2x . . . -o p np x   (9)
Multiple contact: The model presented above is valid for one body in contact with one plane. However, it is straightforward to extend it to multiple bodies in contact with multiple planes, by considering one couple of conditions ( 2) and ( 7) and one couple of variables Φ c and f ⊥ for each contact. In that case, the complete dynamic equation of the robot is:

Aq + b + J T c φ c = S T τ (10)
The equations ( 2) and ( 7), considered one for each body in contact, complete the dynamic model.

ZMP-based stability condition: Considering frame {O} in the interface between the supporting feet and the ground, and assuming that the robot is in contact with a single plane, the ZMP is computed as the average of the contact points o p i weighted by the normal component of the contact forces f iz at each point:

p zmp = np i o p ix f iz np i f iz , np i o p iy f iz np i f iz , 0 (11)
The unilateral condition guarantees that f ⊥ ≥ 0, and therefore, the ZMP will exist inside the support polygon, ensuring the dynamic stability condition. Consequently, [START_REF] Ikeuchi | Toward an assembly plan from observation. Task recognition with polyhedral objects[END_REF] ensures the fulfillment of the classical ZMP condition.

III. DYNAMIC INVERSION USING A STACK OF TASKS

The dynamic model presented above allows a dynamic inversion from a reference joint trajectory, that is, it finds the torques to exert on the robot joints given a reference trajectory q * . However, it is often much straightforward to express the robot objectives in dedicated spaces, as proposed by the task function approach [START_REF] Samson | Robot control: the task function approach[END_REF]. A classical way to compose a movement from a set of tasks is to introduce a strict hierarchy between tasks to ensure a safe behavior in case of conflict. This is a well known solution in kinematics [START_REF] Siciliano | A general framework for managing multiple tasks in highly redundant robotic systems[END_REF] and in dynamics [START_REF] Sentis | Control of free-floating humanoid robots through task prioritization[END_REF]. In [START_REF] Escande | Fast resolution of hierarchized inverse kinematics with inequality constraints[END_REF] a hierarchical solver including inequalities was introduced for kinematics, and it was extended in [START_REF] Saab | Generic Dynamic Motion Generation with Multiple Unilateral Constraints[END_REF] for dynamics. In the following, the principles of this solver are quickly recalled, and the classical tasks to be used are briefly introduced.

A. Generic QP cascade scheme Consider a linear system with n linear inequalities (or equalities)

A 1 x ≤ b 1 , A 2 x ≤ b 2 , . . . , A n x ≤ b n .
The system is solved for x by the hierarchical QP solver using the following procedure. The first stage is:

min x,w1 w 1 (12) A 1 x -b 1 ≤ w 1
where w 1 is a slack variable. The second stage minimizes the slack variable w 2 :

min x,w2 w 2 (13) A 2 x -b 2 ≤ w 2 A 1 x -b 1 ≤ w * 1
where w * 1 is the slack variable determined in the first stage [START_REF] Pollard | Adapting human motion for the control of a humanoid robot[END_REF]. In case of conflict between the constraints, the slack variable w 2 allows to violate A 2 x ≤ b 2 , giving a greater priority to A 1 x ≤ b 1 . The subsequent stages are formulated in a similar way, and the n-th stage can be written as: where

min x,wn w n (14) A n x -b n ≤ w n A a x -b a ≤ w a *
A a = (A 1 , • • • , A n-1 ), b a = (b 1 , • • • , b n-1 ) and w a * = (w * 1 , • • • , w * n-1
). The augmented slack variable w a * contains the slack variables obtained in the previous n -1 stages. This hierarchical QP is denoted as ( 12)

≺ (13) ≺ • • • ≺ (14).

B. Operational tasks

For any observable s(q) of the robot configuration, a desired position and orientation s * can be specified without loss of generality by the task e = s -s * . A usual choice for the task rate of change is ë = -λ p e -λ v ė with λ p > 0 and λ v = 2 λ p to make the error decrease exponentially. At the acceleration level, the relation between the task space and the joint space is:

J q = ë -J q (15)
where J = ∂e ∂q is the task Jacobian.

C. QP for the dynamic inversion

With the notation presented, the generic hierarchical QP proposed in [START_REF] Saab | Generic Dynamic Motion Generation with Multiple Unilateral Constraints[END_REF] for the resolution of the inverse dynamics using the task space is written as (1) ≺ (2) ≺ (7) ≺ (15), and is shown in Figure 2. To deal with all these constraints, the optimization vector determined by the QP solver is x = (q, τ, φ c , f ⊥ ). The operational space dynamic tasks correspond to the accelerations of the controlled task e, which is ë, and Figure 2 shows the case of N tasks e 1 , ..., e N . These tasks can be added or removed at any time and with any priority. Additionally, the contact constraints can also be added or removed at any time making the whole cascade very flexible when handling whole-body motion.

IV. TASK SET FOR IMITATION

The solver detailed previously can be used for various kinds of tasks. Examples can be found in [START_REF] Saab | Generic Dynamic Motion Generation with Multiple Unilateral Constraints[END_REF]. In the frame of imitation and edition, we have mainly used four types of tasks:

• The first task regulates the posture (actuated part of the robot configuration) to follow the motion of the demonstrator.

• The other tasks were used to edit the motion and include: the control of the position and orientation of an operational point, the control of the trajectory of a particular joint, and the control of the position on the ground plane of one foot sliding. These tasks are described in detail in this section.

A. Posture Task

The desired trajectory for the joint angles is obtained using the motion capture system (as it will be explained in Section V). This acquired motion can be directly replayed on a human-like character used in the field of computer graphics. However, it is not straightforward to regenerate the data by a robot as it does not consider the dynamic model or the constraints associated with it. In this case, the inverse dynamics control scheme can be used as it considers the contact constraints and the dynamic stability conditions. Particularly, a task that follows the joints evolution as closely as possible while satisfying the prior dynamic conditions is implemented and will be called the posture task. Let q jk = (q j , . . . , q k ) ∈ R k-j+1 with k, j ∈ N, 1 ≤ j, k ≤ n and k ≥ j represent the vector containing the angular values that need to be controlled, and q * jk the corresponding desired configuration. A posture task is expressed as e = q jk -q * jk . At the velocity level, this will lead to ė = J jk q with the Jacobian J jk ∈ R (k-j+1)×(6+n) selecting only the desired joints: [START_REF] Kajita | Biped walking pattern generation by using preview control of zero-moment point[END_REF] where

J jk = 0 (k-j+1)×(6+j-1) I k-j+1 0 (k-j+1)×(n-k)
I k-j+1 ∈ R (k-j+1)×(k-j+1
) is the identity matrix and 0 k1×k2 ∈ R k1×k2 are zero matrices. In the special case of controlling all the actuated joints, the identity matrix will have its largest size, lying in R n×n . Finally, the reference acceleration is given as a tracking of the joint trajectory obtained by the motion capture system, using a PD controller:

q * jk = -λ p (q jk -q * jk ) -λ v ( qjk -q * jk ) (17)
with λ p and λ v defined as before.

B. Arbitrary Edition of the Generated Motion

The posture task reproduces the desired motion at the joint level satisfying the dynamic constraints. However, the PD acts as a low-pass filter, generating non desired movements at some points or erasing some delicate or very dynamic movements, typically due to fast oscillatory motions. Nevertheless, the structure of the stack of tasks can be used to overcome this problem by adding operational tasks that enhance or even modify the original motion. Moreover, the priorities of the tasks can be modified to better achieve the desired motion.

An operational task controls directly the position and/or orientation of different operational points. In our case, we consider the head, chest, waist, feet and hands of the robot as the controlled end-effectors. Each operational point x i = (p ix , p iy , p iz , ϕ ix , ϕ iy , ϕ iz ) ∈ R 6 is expressed in Cartesian coordinates to specify its position and the roll, pitch, yaw angles for its orientation. Despite this 6D representation, we are free to control either both position and orientation, or only some axis of interest. For this sake, a diagonal selection matrix S x is defined as:

S x = diag(s 1 , s 2 , s 3 , s 4 , s 5 , s 6 ) (18)
where s j is a binary, 1 or 0, the former one to control that particular element of position or orientation, and the latter one to leave it unconstrained. Thus, the differential relation is expressed as:

S x ėi = S x J q (19)
and the reference x * i to be followed only considers the controlled elements of position and/or orientation.

We use this operational task to edit the motion in two cases:

• Specification of target points: A new desired target for a chosen operational point can be specified without defining the desired trajectory to reach it. This point can be determined using forward kinematics applied on a certain operational point to compute its position from the kinematic optimization, which will be described in Section V, or it can be arbitrarily set. • Specification of a trajectory: Let the trajectory for the operational point x be called x o (t). The new trajectory x n (t) that will be set as the desired trajectory for that operational point will be x n (t) = x o (t) + x m (t), where x m (t) is the trajectory modification that can be done on any of the six degrees of freedom of x. This trajectory modification x m (t) can be time varying or constant, according to the requirements. It is important to point out that the added operational task must have higher priority than the joints posture task it would interfere with. For instance, if an operational task is added to the hand, the priority of the arm posture task must be reduced. Alternatively, the task can be removed, but it is preferred to be kept, as it will serve as a "guide" for the new trajectory. If these priority considerations are not taken into account, the desired motion would be blinded by the solutions satisfying the posture task with higher priority, then, the desired effect would not be achieved. Other tasks (typically, the gaze) could have been considered, even if not meaningful in our case of pop dancing.

C. Foot Sliding

At some points, the observed motion revealed slight changes in the contact which becomes not exactly rigid. Thus, to let a foot slide, the rigid contact constraints stating ẋc = 0 and ẍc = 0 can be relaxed constraining the motion to the horizontal XY plane. For the contact i, this restriction at the velocity level can be formulated as v z = 0, ω x = 0, ω y = 0 allowing the other velocity elements to take arbitrary values. This guarantees no translation on the vertical axis Z or rotation about it. To account for these constraints, the dynamic model would need to be reformulated. Alternatively, another simpler and more flexible approach is to remove the contact constraint and add a task that restricts the motion in Z direction, imposing no restrictions to the X or Y axis, which will limit the motion to the XY plane. Considering x * as the desired task, this particular case of the 6D task ( 19) is x * = (k x , k y , H z , 0, 0, k rz ) and S x = diag(0, 0, 1, 1, 1, 0). where S x is the selection matrix defined in [START_REF] Saab | Generic Dynamic Motion Generation with Multiple Unilateral Constraints[END_REF], k x , k y , k rz ∈ R are arbitrary values and H z is the known height of the foot operational point with respect to the world frame. When the foot is in contact with the ground, this constraint restricts the height of the foot to remain constant, and the rotation about the X and Y axis to be null. This task does not exactly correspond to a sliding, since sliding can accept forces orthogonal to the motion. The proposed solution is thus more restrictive than necessary. However, the visual effect would be the same, and it will keep the motion feasible and dynamically consistent.

V. MULTIBODY MOTION ACQUISITION

Our motion capture system is a 10-camera with marker tracking system [START_REF]Motion analysis[END_REF]. It provides the 6D position, translation and orientation, of a set of unconstrained bodies in space, typically the limbs of the demonstrator, at a frequency of 200 Hz and with a precision of 2 mm. For this set of 6D position, we need to recompute the joint position of the demonstrator, knowing its geometric model. This is achieved by solving a non linear optimization problem where the geometric model of the demonstrator is known, and the joint position is the optimization variable.

A. Optimization problem formulation

The geometric model gives the position and orientation of each joint q i with respect to the robot world reference frame {W }, as the transformation matrix W T qi (q) which is a function of the joint vector q. Then, the problem of finding the suitable joint values q * (t) for the robot at time t can be reduced to minimizing the difference between these transformation matrices as:

q * (t) = arg min q n i W T * qi (t) W T qi (q) (20) s.t. q i ≤ q i ≤ q i ( 21
)
where W T * qi (t) is the 6D position of the body i given by the motion capture system, q i , q i are the minimum and maximum angular values of a joint i, respectively, and the condition to satisfy the joint limits for every joint i ∈ {1, ..., n} at every optimization process has been added. The symbol represents a distance operator in the matrix group. We used a weighted norm of the translation-axis-angle vectorial writing of the transformation matrices given as follows. The homogeneous matrices can be decomposed in rotation and translation parts:

W T * qi (t) = R * qi (t) p * qi (t) 0 1 (22)
W T qi (q) = R qi (q) p qi (q) 0 1 . (23)
For the part corresponding to the position, the difference to minimize is directly represented by the L -2 norm as p * qi (t) -p qi (q) 2 . For the orientation part, the difference is measured by the rotation angle θ i , which is obtained using the axis/angle representation of the product R di = R qi (q)R * qi (t) -1 ∈ SO(3). If the rotation matrices are close enough, the angle of R di about the arbitrary axis is very small. Representing the product of the rotation matrices as

R di =   n x o x a x n y o y a y n z o z a z   (24)
the angle will be

θ i = atan2 (s a , c a ) (25)
where

s a = (n y -o x ) 2 + (n z -a x ) 2 + (o z -a y ) 2 2 (26)
c a = n x + o y + a z -1 2 . ( 27
)
The angle θ i is always between 0 and π. Finally, the operator in ( 20) is given by:

W T * qi (t) W T qi (q) = w pi p * qi (t) -p qi (q) 2 +w oi θ i (28)
where w pi is the weight corresponding to the position part of the joint i and w oi = 1 -w pi is the weight for the orientation part. The weights are included to provide more flexibility to the model and are experimentally determined to give more priority to one of the parts. The W T qi (q) should be computed from the geometric model of the demonstrator. In our case, we have used directly the model of the robot to compute W T qi (q). The optimization method has proved to be robust enough to handle the approximation, and the dynamic solver is sufficient to recover all the resulting noise and inaccuracy.

B. Calibration of the body frame

In [START_REF] Samson | Robot control: the task function approach[END_REF] it is supposed that the same frames are used to express both the motion capture observation and the geometric model. This is of course not the case in practice, since the motion-capture system uses its own arbitrary frames to express the multi-body motion. The observation of the motion capture has to be then reformulated in the proper frame. To this end, the person whose motion wants to be recovered starts with a position that is well known for the robot. Using this known configurations, the transformation matrices between each node m i and the corresponding joint q i in the robot is obtained from a classical "calibration" step and represented by mi T qi , which remains constant for all the process as long as the markers do not have relative motion with respect to the body they are attached to. This matrix considers the differences in orientation between the markers and the frames defined in the kinematic model of the robot for each joint, as well as the differences in the segment lengths of the robot and the dancer. Another matrix that is obtained during the calibration process is the one relating the origin of the motion capture system {W m } and the robot reference frame {W } represented by W T Wm . Let the position and orientation of each node m i with respect to the motion capture reference frame be represented by the homogeneous transformation Wm T mi (t) which varies in time according to the motion to be imitated. With these transformation matrices, the desired configuration of the robot joint i is:

W T * qi (t) = W T Wm Wm T mi (t) mi T qi (29)
VI. RESULTS To validate the proposed method, the motion of a dancer was acquired using the motion capture system, which provides the position and orientation of the frames associated with each node. The motion was retargeted to a modified-HRP2 model (ie we added one degree of freedom to both chest and neck joints, obtaining a kinematic structure closer to the one of HRP4, to obtain a nicer final motion), and edited to corrected the retargeting error and introduced new undemonstrated features.

The kinematic optimization was applied to obtain the corresponding joint trajectories, which were then validated dynamically using the posture task in the inverse dynamics solver. The joints trajectory is first obtained from the motion capture, and it is, of course, not stable nor dynamically consistent as Figure 3b shows. A first dynamic motion is then obtained using only the joint trajectory as references. This motion is displayed in Figure 3c. The motion is stable; however, the geometric and dynamic retargetting have lost some data and produced some errors compared to the initial demonstrated trajectory. Three editions were thus applied:

• The knee oscillations (smoothed by the PD) were enhanced.

• The right hand motion was corrected.

• An additional motion of the left foot consisting in sliding on the ground, which was not present in the initial demonstration, was introduced. The sequence showing the time moments when these modifications were added is shown in Figure 4. The following subsections provide more details of these three modifications.

A. Knee Oscillation

The knees constitute a particular case as the dancer permanently moved them but at the dynamic level this motion was strongly weaken. To correct this problem we referred to the task approach where the motion of the knee joint was analyzed. Figure 5a shows the joint evolution obtained after the kinematic optimization for the right knee. Between iterations 2000 and 2800, and between 6200 and 7300, the motion of the joint is oscillatory and those are the moments corresponding directly to the observed motion at the dancer's right knee. A scalogram using the Gaussian wavelet was constructed and is shown in Figure 5b, where it is observed (in red circle) that there are salient frequencies at those points. It was determined that the scale a corresponding to the maximum values at the desired positions is 36. The frequency and the scales are related by f = fsfw a , where f s is the sampling frequency and f w is the center frequency of the wavelet. For the Gaussian derivative of order 4, f w = 0.5, and considering that the sampling frequency used during the acquisition is 200 Hz, the resulting frequency is 2.7 Hz. Then, a task on the knee was added at that frequency to resemble more the motion. Figure 6a shows the evolution of the joint at the right knee (joint number 4) in the joint space as a function of the iteration time, and Figure 6b shows the evolution in the X axis of the operational space. Note that there is a difference in the iteration time of Figure 5 and Figure 6 as the former one corresponds to the kinematic optimization of the motion sampled at 200Hz and the latter corresponds to the dynamic level that runs at 1000Hz (the iteration relation is 5 to 1). The red line shows the evolution of the joint when there is only the posture task in the leg, and the blue line shows the evolution when the operational task is added to the knee, but both lines are obtained using the dynamic control scheme. Even though the joint is not directly controlled, it is observed in Figure 6a that with the operational task addition, the joint at the knee presents an oscillation with higher amplitude, whereas with only the task posture the oscillation is weak. The results in the X axis of the operational space show a clear consistent oscillation with similar amplitude where the knee task was added Both the joint space and the operational space show the effect of the task at the knee.

B. Right hand motion

The fast up and down motion of the right arm was also smoothed as a consequence of the PD used. This was especially noted when the arm could not reach the upper positions that the dancer performed. Then, an operational task to raise more the arm was introduced. The result is shown in Figure 7. The trajectory of the right arm in the Z axis using only the posture task is shown in red, whereas the trajectory with the operational correction is shown in blue. The corrected trajectory, improved the upper positions of the right hand. The task in the right hand was also used to avoid the auto collision of the hand with the head, as can be noticed of Fig. 3, fourth thumbnails.

C. Foot sliding introduction

We artificially introduced a sliding movement of the right foot, to prove that extra features can be added as desired. To introduce the sliding effect on a foot, the ZMP of the motion obtained using only the dynamic posture was analyzed. Figure 8 shows in blue the trajectory of the ZMP in the Y axis as a function of the time. The wide red lines show the limits of the support polygon in the Y axis, and the dashed red lines show the boundaries of the foot (inside the support polygon). Between the iteration times 46940 and 48310, the ZMP lies completely in the area corresponding to the right foot. Then, it was at this time that the sliding task was introduced for the left foot, guaranteeing the dynamic stability of the robot, as the sliding foot cannot be considered anymore part of the support polygon.

VII. CONCLUSION

A method for the imitation of whole-body motion for humanoid robots has been presented. The contribution of this work is to propose a complete methodology to quickly reshape a dynamic motion demonstrated by a human expert, adapt the dynamics of the human body to the own dynamics of the robot and modify or edit as desired the initial motion to introduce extra features that were not demonstrated. It allows to build complex dynamic behaviors, based on a composition of tasks and constraints that are used as basic bricks for motion generation. The method was satisfactorily applied to the imitation of dancing motion, but it can be, in general, applied to the imitation of any type of motion. The obtained motion is dynamically consistent, and could be directly applied on the real humanoid robot.

Fig. 2 :

 2 Fig. 2: Dyanmic QP Cascade Scheme.

Fig. 5 :

 5 Fig. 5: Scalogram of the right knee joint evolution

Fig. 3 :

 3 Fig. 3: Results for the robot imitating the dance performed by a human.

Fig. 4 :

 4 Fig. 4: Task sequence

Fig. 6 :Fig. 7 :

 67 Fig. 6: Temporal evolution of the right knee

Fig. 8 :

 8 Fig. 8: Evolution of the ZMP with only the posture task


Eric Falcon

Claude Laroche

Observation of depth-induced properties in wave turbulence on the surface of a fluid

Keywords: 47, 35, -i -Hydrodynamic waves 47, 27, -i -Turbulent flows 92, 10, Hm -Ocean waves

We report the observation of changes in the wave turbulence properties of gravitycapillary surface waves due to a finite depth effect. When the fluid depth is decreased, a hump is observed on the wave spectrum in the capillary regime at a scale that depends on the depth. The possible origin of this hump is discussed. In the gravity regime, the wave spectrum still shows a power law but with an exponent that strongly depends on the depth. A change in the scaling of the gravity spectrum with the mean injected power is also reported. Finally, the probability density function of the wave amplitude rescaled by its rms value is found to be independent of the fluid depth and to be well described by a Tayfun distribution.

Introduction. -Wave turbulence is ubiquitous in nature. It ranges from surface or internal waves in oceanography, Alfvén waves in solar winds, Rossby waves in geophysics, elastic or spin waves in solids, waves in optics, and quantum waves in Bose condensates (for reviews see [1,2]). Wave turbulence theory (a statistical theory describing an ensemble of weakly nonlinear interacting waves) predicts analytical solutions of the kinetic equations of weak turbulence at the equilibrium or in a stationary out-of equilibrium regime in various systems involving wave dynamics [START_REF] Nazarenko | Wave Turbulence[END_REF]. Surprisingly, well-controlled laboratory experiments on wave turbulence were scarce [START_REF] Wright | [END_REF] up to last years where new observations was reported such as intermittency [5], fluctuations of the energy flux [6], finite size effect of the system [7,8] and the full space-time power spectrum of wave amplitudes [9]. Several questions still are open, notably about the validity domain of the theory in experiments (horizontal finite size effects, role of strongly nonlinear coherent structures), and the possible existence of weak turbulence solutions for nondispersive (ie, for 2D acoustic waves) or weakly dispersive systems [10,11]. Indeed, the lack of dispersion could lead to cumulative nonlinear effects leading to shock wave formation. It is thus of primary interest to know the evolution of the wave energy spectrum when the dispersion relation of a system is (a) Corresponding author: eric.falcon@univ-paris-diderot.fr changed from a dispersive to a nondispersive regime.

A simple way to experimentally do this is to change the depth of a fluid on which surface waves propagate. Indeed, gravity waves are known to become nondispersive in a shallow water limit. In this limit, weak turbulence predicts power spectra of wave amplitude much less steep than in the deep regime for both gravity [12,13] and capillary [14] wave turbulence. The prediction for the gravity regime has been tested in few laboratory experiments using a large-scale wave flume with a sloping bottom [15,16]. At a more applied level, in situ observations exist in oceanography when ocean waves propagate from deep water into shallow coastal areas [16]. Notably, when ocean surface waves approach the shore, bottom friction and depth-induced wave breaking are no more negligible, and near-resonant triad interactions could play a dominant role instead of the 4-wave interaction process [17,[START_REF] Massel | Ocean Surface Waves[END_REF].

In this letter, we study gravity-capillary wave turbulence on the surface of a fluid within a constant depth tank. When the depth is decreased from a deep to a thin fluid layer, a hump is observed on the power spectrum in the capillary regime at a scale that depends on the depth. In the gravity regime, the wave spectrum still shows a power-law but with an exponent that depends on the depth. The scaling of the power spectrum with the mean injected power is measured for different fluid depths as well as the probability distribution of the wave ampli-tudes. This latter is found to be independent of the fluid depth and to be well described by a Tayfun distribution.

Experimental set-up. -The experimental setup, described in [7], consists of a rectangular plastic vessel, 20 × 20 cm 2 , filled with mercury (unless otherwise stated) up to a depth h. h is varied from 3 mm up to 22 mm with a 0.1 mm accuracy. Mercury is used because of its low kinematic viscosity ν = 10 -7 m 2 /s. Surface waves are generated on the fluid by the horizontal motion of a rectangular (L × H cm 2 ) plunging plastic wave maker driven by an electromagnetic shaker. We take L = 13 cm and H = 3.5 cm, the wave maker being at 0.2 mm from the bottom vessel. The wave maker is driven with low-frequency random vibrations (typically from 1 to 5 Hz). The amplitude of the surface waves η(t) at a given location is measured by a capacitive wire gauge (plunging perpendicularly to the fluid at rest) [7]. η(t) is analogically low-pass filtered at 1 kHz and is recorded during 300 s by using an acquisition card with a 2 kHz sampling rate. The power injected I(t) into the fluid by the wave maker is also measured using a force sensor and a velocity sensor both fixed on the wave maker [6]. The system is driven by the rms value of the velocity fluctuations σ V of the wave maker [7]. The mean power injected by the wave maker I is measured as a function of the depth h for the same value of the forcing amplitude σ V . As shown in the inset of Fig. 1 (for both water or mercury as the working fluid), one finds that I ∼ ρh for h > h * = 2.5 mm. ρ is the density of the fluid (the mercury density being 13.6 times the water one). h * is the depth below which the dewetting of the fluid on the vessel bottom occurs that forms islands of fluid surrounded by domains without fluid. One has also I ∼ ρhσ 2

V (not shown here). Thus, typical experiments are performed either keeping constant the forcing amplitude σ V when changing the fluid depth, or adjusting σ V to have a constant mean injected power I when the depth is decreased. We have performed experiments with two different fluids (mercury or water). Experimental results with mercury are shown below, but qualitative similar results have been observed with water.

Dispersion relation of linear waves. -For an arbitrary fluid depth h, the dispersion relation of inviscid linear surface waves on a fluid reads

ω 2 = (gk + γ ρ k 3 ) tanh(kh) (1)
where ω is the wave pulsation, k its wavenumber, g = 9.81 m/s 2 the acceleration of the gravity, γ = 0.4 N/m the surface tension, and ρ = 13.6 × 10 3 kg/m 3 the density of the mercury. Gravity waves are dominant in Eq. ( 1) for wavelengths λ ≫ 2πl c ≃ 1 cm [capillary length l c ≡ γ/(ρg) ≃ 1.7 mm] whereas capillary waves dominates when λ ≪ 2πl c . The deep-water limit corresponds to kh ≫ 1 that is for λ ≪ 2πh, whereas the shallow-water limit corresponds to λ ≫ 2πh. Typically, for our range of frequencies 4 ≤ f ≤ 200 Hz, Eq. ( 1) leads to 0.1 ≤ λ/(2πh) ≤ 2.3 for Wave turbulence power spectrum. -The power spectra of wave amplitude S(f ) are shown in Fig. 1 for different depths h at a constant forcing amplitude σ V . For large depths (e.g. h = 22 mm), it displays similar results than those found previously in the deep regime: two power laws corresponding to the gravity and capillary wave turbulence regimes separated by a crossover frequency f c related to l c [7]. When h is decreased, important changes are observed in the spectrum shape (see Fig. 1). First, a hump appears on the power spectrum at a frequency denoted f p (see Fig. 1). f p is found to decrease when h is decreased. Second, the power spectrum still shows a frequency-power law in the gravity regime (frequencies below f c ) but with a strongly decreasing exponent. Third, in the capillary regime (between f c and f p ), a power-law can be also fitted (except for h = 4 mm) in a narrower inertial range due to the presence of the hump (see Fig. 1). Finally, at small enough depth, a third power-law can be fitted on the high frequency part of the spectrum (between f p and a frequency given by the intercept of the spectrum and the signal to noise level -see Fig. 1). This power-law vanishes when the deep regime is reached. Similar observations are found either when performing experiments at constant I forcing instead of constant σ V forcing, or when using an optical non-invasive method (laser vibrometer) instead of the capacitive wire gauge. Moreover, these observations do not significantly depend on the probe location on the basin surface except close to the boundaries.

Let us now find the scaling of the spectrum amplitude S(f ) with I for a fixed depth. In deep regime (h = 20 mm), S(f ) is found to scale as I 0.8±0.1 for both gravity and capillary regimes on almost one decade in I . This result is consistent with ones found in previous measurements [7,[START_REF] Xia | [END_REF] whereas weak turbulence predicts a scaling in I 1/3 for 4-wave interactions (gravity in deep regime) and in I 1/2 for 3-wave interactions (capillary regime) [START_REF] Nazarenko | Wave Turbulence[END_REF]. For small depths (h = 5 mm), S(f ) is found to scale as I 0.5±0.1 for gravity waves, whereas no scaling is found for the capillary regime due to the presence of the hump. The observed change of scaling with the depth thus suggests that a change of the wave-interaction process occurs for gravity waves.

Figure 2 shows the evolution of the frequency exponents of power spectra S(f ) as a function of the depth for the gravity and capillary regimes for different forcing parameters. When h is decreased, the absolute values of exponents continuously decrease both for the capillary and gravity regimes. The decrease of the exponents when h decreases is qualitatively consistent with the expected one from the weak turbulence predictions in both depth limits. Indeed, weak turbulence predicts in both depth limits: S grav (f ) ∼ f -4 (deep) [20] and ∼ f -4/3 (shallow) [12,13] for the gravity regime (see dashed-red [light grey] lines in Fig. 2), and S cap (f ) ∼ f -17/6 (deep) [21] low) [14] for the capillary regime (see dashed-black lines in Fig. 2). The power-law modification of the gravity regime could be related to a change of wave interaction process with h (see above). For the capillary regime, it is related to the presence of the hump (see below). Finally, as shown in the inset of Fig. 2, when h is increases, the exponent of the third power-law ending the cascade is found to strongly decreases in absolute value, and then saturates to roughly the same exponent as the capillary one.

Crossover and hump frequencies. -Figure 3 shows the evolutions of the crossover frequency f c between gravity and capillary regimes and of the hump frequency f p as a function of h for different forcing parameters. f p is measured on the power spectrum as the frequency corresponding to the end of the power-law capillary cascade and the beginning of the third slope (see Fig. 1). This cutoff f p is found to strongly increase with h at small depths, and then seems to saturate at larger depths (h ≥ 10 mm). Contrarily to f c (h) (see below), the evolution of f p (h) is not given by the depth dependence of the dispersion relation of Eq. ( 1) at any fixed wavenumber. Moreover, the evolution of f p (h) is not due to the fact that less and less energy is injected into the waves when h is decreased ( I ∼ h -see above). Indeed, experiments performed with a constant I forcing whatever h gives similar spectra as the ones in Fig. 1 and similar f p (h) evolution (see ♦-symbols in Fig. 3). This evolution is not also due to the viscous damping rate ∼ 2νk 2 (ω) [START_REF] Landau | Mécanique des Fluides[END_REF] or to the bottom friction one ∼ νω(k)/h [START_REF] Landau | Mécanique des Fluides[END_REF] since both are decreasing functions of h. A possible origin of f p is the following. At the cut-off frequency f p of the capillary spectrum, the nonlinear energy transfer is changed by viscous dissipa- tion, thus it can be assumed that the typical non linear time scale of wave interaction, τ nl , is of the order of the typical viscous time scale τ v : τ nl (f p ) ∼ τ v (f p ) [24]. Thus, if the non linear time scale depends on the fluid depth, f p estimated from this above condition should also depend on the depth. A change of the wave interaction process during the transition from a deep to a thin fluid layer may modify the nonlinear time scale that thus should depend on the depth. However, this qualitative interpretation does not explain how the energy flux flowing across the scales accumulates at the hump scale as an analog of a bottleneck effect [START_REF] Ryzhenkova | [END_REF]. We cannot currently measure τ nl (and its possible depth dependence) with a temporal measurement of the wave amplitude at a given location, but it could be it with a spatio-temporal measurement [9]. This thus would deserve further studies.

Fig. 3 also shows the evolution of the gravity-capillary crossover frequency f c (h). For small h, f c is found to slightly increase with h, and then to saturate in the deep regime to a value known to be larger than the expected one

f c = 1 2π √ 2gk c ≃ 17 Hz with k c ≡ ρg/γ [7].
Here, the deep limit value is f c ≃ 28 Hz that corresponds to a wave number kc = 1.56k c . By using the relation dispersion of Eq. ( 1), f c then is found to be well fitted in same value of the forcing amplitude. The PDFs are found to be asymmetrical due to strong steepness of the waves as usual in laboratory experiments [7,25] or in oceanography [START_REF] Massel | Ocean Surface Waves[END_REF]26]. No dependence of this asymmetry with the depth is observed in our experiment. Indeed, first, the skewness S is roughly constant whatever the depth: S = 0.55 ± 0.1. Second, when these PDFs are normalized to their rms amplitude values, σ η , they all roughly collapse on a single non-Gaussian distribution as shown in Fig. 5. This means that the distribution depends only on σ η . The bottom inset of Fig. 5 shows that the rms wave amplitude σ η increases linearly with the depth h. These rescaled distributions are found to be well fitted by a Tayfun distribution (the secondorder correction to the Gaussian distribution) that reads p

[η] = ∞ 0 exp {[-x 2 -(1 -c) 2 ]
/(2s 2 )}/(πsc)dx where c ≡ 1 + 2sη + x 2 , η ≡ η/σ η and s the mean steepness of the waves [27]. Note that we have checked that the Tayfun distribution fit well our experimental normalized PDF performing for much larger depth (i.e., 18 ≤ h ≤ 140 mm) that is for depth much larger that the wave amplitude η. This also suggests that this typical distribution does not depend on the possible change of the wave interaction process during the transition from a deep to a thin fluid layer.

Conclusion. -We have performed experiments of gravity-capillary wave turbulence on the surface of a fluid and have focused on the effects of finite depth. When the fluid depth is decreased from a deep to a thin fluid layer, a hump is observed on the wave spectrum at a depthdependent scale located in the capillary regime. The presence of this hump strongly reduces the inertial range of the capillary wave turbulence. The mechanism of formation of this hump and its dependence on the depth are open problems. In the gravity regime, the power spectrum is still found to be a power-law but with an exponent that continuously decreases with decreasing depth in a qualitatively agreement with the weak turbulence predictions in both deep and shallow water limits. The scaling of the spectrum with the mean injected power is also measured and suggests that a change of the wave-interaction process occurs for gravity waves during the transition from a deep to a thin fluid layer. Finally, the probability distributions of the wave amplitude rescaled by its rms value are measured for different fluid depths. The distributions are found to be independent of the depth and to be well described by a Tayfun distribution whatever the fluid depth.

  Fig. 1: (color online) Power spectra of wave amplitude S(f ) for different fluid depths from h = 22, 11, 8, 6 and 4 mm (solid lines from top to bottom). Dashed lines are spectrum slopes reported on Fig. 2. fc is the crossover frequency between gravity and capillary regime and fp the hump frequency (see text). Forcing bandwidth and amplitude: 1 -4 Hz and σV = 0.8 (arb. unit). Curves have been vertically shifted for clarity by a factor 1, 1.5, 2.5, 6 and 12 (from bottom to top). Inset: Mean injected power rescaled by the fluid density I /ρ versus h for water ( ) and mercury (•). Dashed line is a linear fit. Forcing parameters: 1 -6 Hz and σV = 0.8 (arb. unit).

Fig. 2 :

 2 Fig.2:(color online) Exponents of power spectra as a function of the depth for gravity [red (light grey) symbols] and capillary (black symbols) regimes. Forcing parameters: 1 -4 Hz ( ), 1 -6 Hz (•) both at σV = 0.8 (arb. unit), or at constant I = 0.02 (arb. unit) (♦). Dashed lines are theoretical exponents for both regimes in the deep or shallow water limits (see text). Inset: Third exponent of the spectrum as a function of h fitted between fp and the noise level frequency.

  and f -2 (shal-

Fig. 3 :

 3 Fig. 3: (color online) Hump frequency fp (full symbols) and crossover frequency fc (open symbols) as a function of the fluid depth h for different forcings: 1 -4 Hz ( ), 1 -6 Hz (•) both at σV = 0.8 (arb. unit), or at constant I = 0.02 (arb. unit) (♦). Dashed line: linear fit of slope 9 Hz/mm. Dot-dashed line is given by the depth dependence of the dispersion relation at a fixed wavenumber kc = 1.56 ρg/γ (see text).

Fig. 4 :

 4 Fig. 4: (color online) Probability density functions of wave amplitude, η, for different values of the fluid depth, from h = 3, 4, 5, 6, 8, 10, 12, 15 to 18 mm (see arrow). Forcing parameters: 1 -6 Hz and σV = 0.8 (arb. unit).

Fig. 3 Fig. 5 :

 35 Fig. 3 by fc (h) = (g kc + γ ρ k3 c ) tanh[ kc (h -h * )]/(2π). h * ≃ 2.5 mm is the depth for which the fluid dewetting occurs on the container bottom and is of the order of the capillary length as expected for a capillary meniscus related phenomenon. Consequently, the evolution of the crossover frequency f c (h) is given by the depth dependence of the dispersion relation at a fixed wavenumber kc .Distribution of wave amplitudes. -The probability density functions (PDFs) of the wave amplitude are shown in Fig.4for different values of h for the

We thank M. Berhanu for fruitful discussions. This work has been supported by ANR Turbonde BLAN07- 3-197846.