Modesto T López-López
email: modesto@ugr.es

Pavel Kuzhir

Jaime Caballero-Hernández

Laura Rodríguez-Arco

Juan D G Duran

Georges Bossis

Yield stress in magnetorheological suspensions near the limit of maximum-packing fraction

This work deals with the magnetic field-induced static yield stress of magnetorheological (MR) suspensions with concentration near the limit of maximum-packing fraction. With this aim, homogeneous suspensions of iron microparticles with 50 vol.% concentration were prepared, and their yield stress measured as a function of the applied magnetic field. In view of the failure of existing models to predict, on the basis of realistic hypotheses, the values of the yield stress of highly concentrated MR suspensions, we developed a new model. Our model considers that field application induces body-centered tetragonal (BCT) structures.

Upon shearing, these structures deform in such a way that interparticle gaps appear between neighboring particles of the same chain, whereas the approach of particles of parallel chains ensures the mechanical stability of the whole multi-chain structure. Based on this hypothesis, and using finite element method simulations of interparticle magnetic interactions, our model is able to quantitatively predict the yield stress of highly concentrated MR suspensions.

Furthermore, estimations show that the main contribution to the field-dependent part of the

I. INTRODUCTION

The maximum packing fraction of particles in suspension is conventionally defined as the volume fraction of particles or particle aggregates in closest-packing at which the suspension viscosity approaches infinity [START_REF] Zhou | Yield stress and maximum packing fraction of concentrated suspensions[END_REF]]. For monodisperse hard spheres, the random loose maximum-packing fraction is approximately 0.63 0.64 [START_REF] Onoda | Random loose packings of uniform spheres and the dilatancy onset[END_REF]]. However, porosity and shape of the particles, as well as particle-particle interactions, among other factors, have strong influence on the maximum packing fraction [START_REF] Zhou | Yield stress and maximum packing fraction of concentrated suspensions[END_REF]]. For example, flocculated suspensions have a lower maximum-packing fraction due to the fact that the particles are assembled into porous aggregates, whereas suspensions of charged particles do have because of the strong electrostatic repulsion between the particles [START_REF] Russell | Colloidal dispersions[END_REF]].

Highly concentrated suspensions are of interest in many technical and industrial fields, for example in paints and cosmetics [START_REF] Clausen | The rheology and microstructure of concentrated non-colloidal suspensions of deformable capsules[END_REF]]. However, experimental rheology of highly concentrated suspensions is hampered by the poor reproducibility of measurements and the sensitivity to shear history. From the theoretical viewpoint, approaches based on the interpolation of the dilute-limit theory to a concentrated regime often fail because of the neglect of many-body interactions between particles [Larson (1994)], and taking these interactions into account makes the analytical work intractable [START_REF] Clausen | The rheology and microstructure of concentrated non-colloidal suspensions of deformable capsules[END_REF]]. Thus, realistic theoretical modeling remains an open and difficult issue in many cases.

In the particular case of suspensions of non-Brownian magnetic microparticles, known as magnetorheological (MR) suspensions, maximum packing-fraction is considerably diminished (with respect to the mentioned 0.63 0.64) because of the formation of particle aggregates induced by magnetic attraction, due to the remnant magnetization of the particles, and van der Waals forces. As a consequence, there are only a few rheological studies for particle volume concentrations higher than 45% [START_REF] Chin | Rheological properties and dispersion stability of magnetorheological (MR) suspensions[END_REF]; de [START_REF] De Vicente | Normal force study in concentrated carbonyl iron magnetorheological suspensions[END_REF]; Laun et al. (2008a); Laun et al. (2008b)]. Only two of these works [START_REF] Chin | Rheological properties and dispersion stability of magnetorheological (MR) suspensions[END_REF]; Laun et al. (2008b)] present data of the yield stress, whereas the others [START_REF] De Vicente | Normal force study in concentrated carbonyl iron magnetorheological suspensions[END_REF]; Laun et al. (2008a)] deal with the normal force. Besides, to the best of our knowledge, the existing theories (critically reviewed in section III.A below) underestimate the mechanical properties of MR suspensions for concentrations near the maximum packing-fraction or give unphysical results. Because of this, we have aimed to develop an appropriate theoretical model for the rheological properties of highly concentrated MR suspensions. In the present paper, we focus our attention on one of the most important MR properties: the static yield stress, which is commonly defined as the threshold stress required to fracture the suspension structure in its weakest point, and thus induce the flow of the suspension [START_REF] Barnes | An introduction to rheology[END_REF]]. In order to validate our model, we conducted experimental measurements of the yield stress of a suspension containing 50 vol.% of iron particles. We will see that results of the theoretical model agree quite well with experimental ones within the range of applied magnetic fields, 0<H 0 <25 kA/m.

II. EXPERIMENTAL

Silica-coated iron particles (Fe-CC, density 7.2 g•cm -3 ), supplied by BASF (Germany) were used as solid phase for the preparation of the MR suspension. According to the manufacturer, Fe-CC particles have median diameter, d50 = 5 µm. The choice of these particles was motivated by their silica coating, which ease the dispersion of the particles without requiring a surfactant. Mineral oil (density and viscosity at 25 ºC: 0.85 g•cm -3 and 39.58 ± 0.16 mPa•s, respectively) purchased from Sigma Aldrich (Germany) was used as carrier liquid. The MR suspension was prepared as follows. A small amount (approx. 4 g) of Fe-CC powder was added to a relatively large amount of mineral oil (we took 50 cm 3 ) and the resulting mixture was mechanically stirred until a homogeneous suspension was obtained. Then, a further small amount of Fe-CC powder was added to the suspension and, again, the resulting mixture was homogenized by mechanical stirring. This step was repeated as many times as required until the concentration of Fe-CC powder was so high that it was impossible to homogenize the mixture. At this point, we added approx. 1 cm 3 of mineral oil and the homogenization of the mixture was possible by mechanical stirring. By this protocol, we assured that the concentration of particles in the suspension was as close as possible (within reasonable limits) to the maximum-packing fraction. At the end of this process, the volume fraction of particles in suspension was approximately 0.50, as obtained by density measurements. This concentration is relatively far from the 0.63 0.64 maximum-packing fraction reported for monodisperse hard spheres, likely due to the formation of aggregates as a consequence of attractive forces between particles, as discussed in the introduction.

Rheological measurements were performed at 25 ºC using a rheometer MCR 300 (Physica-Anton Paar). The measuring geometry consisted of a homemade set of non-magnetic parallel plates of 20 mm in diameter, with rough surfaces in order to avoid wall slip.

Visualization experiments on the deformation field of the MR fluid [which will be communicated in future] confirm the absence of the wall slip in our experimental system. The applied magnetic field was generated with the help of a solenoid, placed in such a way that the measuring geometry was at middle height of the coil and with axis coincident with that of the coil. The gap thickness between the lower (stationary) plate and the upper (rotational) plate of the measuring geometry was fixed at 0.35 mm.

As mentioned above, the static yield stress of a suspension is the shear stress required to induce its flow. In our experiments, we determined it by subjecting the samples to shear rate ramps in the range 0.005 0.02 s -1 and extrapolating the obtained values of the shear stress to zero shear rate. To be precise, we used the following experimental protocol. The MR suspension was placed in the measuring system of the rheometer and, immediately afterwards, subjected to a shear rate ramp of 1 min of duration in the range 0 100 s -1 . This pre-shear stage was carried out in the absence of magnetic field. Its aim was to impose identical initial conditions in order to ensure reproducibility of the measurements. Then, a uniform magnetic field of chosen value, ranging between 0 and 25 kA/m, was applied with the help of the solenoid and the suspension was left at rest during 30 s a sufficient time to induce stable MR structures. Finally, the shear rate ramp was applied under the same magnetic field applied in the previous step. Note that the shear rate ramp consisted of 5 different values of the applied shear rate, each of them maintained during, at least, 30 min. The stress response was measured as a function of time during these periods of time. Once the measurement at a given applied field was accomplished, we repeated the pre-shear stage at zero magnetic field, and remade the same type of measurement at a higher applied field. The whole protocol was performed for three different freshly prepared MR samples. Results shown in this manuscript correspond to the average of these three different measurements. Note also, that we chose such a long duration of each measurement in order to achieve strains as large as γ ∼10 3 and ensure a steady-state flow of the MR suspension. This experimental protocol guarantees reproducible results, independent of the previous mechanical history of the suspension.

III. THEORY

A. Theoretical background

Before presenting in details our theoretical model it is worth checking, in an approximate manner, the agreement between the predictions of the existing theories and the experimental results for the static yield stress of highly concentrated suspensions. One of the first models dealing with this property assumed the formation of single chain structures, with affine displacement of each particle along the direction of the shear, as depicted in Fig. 1a [ Bossis et al. (2002); [START_REF] Ginder | Rheology of magnetorheological fluids: Models and measurements[END_REF]; [START_REF] Klingenberg | Studies on the steady-shear behavior of electrorheological suspensions[END_REF]]. This model will be hereinafter referred to as the affine model . As the particle chains are tilted by the straining motion, the horizontal component of the magnetic attractive force, F x , acting between particles first increases, then reaches a maximum and finally decreases with the strain. At the strain for which the force is maximum, the chains are considered to become unstable and to break, inducing the flow of the MR suspension. The yield stress is simply calculated as the maximal interparticle force F x multiplied by the number n s of chains per unit surface of the upper rheometer plate:

Ys x nF σ =⋅. ( 1 )
A more general model, based on fundamental thermodynamic relations, considers the shear stress in a sheared suspension as a result of the change in its internal energy:

/ U σ γ =∂ ∂ , ( 2
)
where γ is the strain and U is the suspension internal energy per unit volume (see Fig. 1b) [START_REF] Bossis | Yield stress in magnetorheological and electrorheological fluids: A comparison between microscopic and macroscopic structural models[END_REF]; Bossis et al. (2002)]. This model will be called here the thermodynamic model. The static yield stress is also calculated in this case as the shear stress corresponding to the maximum in the stress versus strain curve. This model can be used whatever the description of the structure and, here, one should distinguish the microscopic description from the macroscopic one. The difference between both descriptions stands in the fact that, in the former the loss of contact between two neighboring particles is taken into account for the estimation of the energy dependency with the strain, contrary to the latter in which only the change of the energy due to the inclination of the aggregates is considered. As a result, both approaches give similar predictions if the interactions between particles are long ranged, as it is the case of dipolar interactions between particles of low magnetic permeability. On the other hand, predictions are very different if strong short range forces exist, as it happens in the case of particles with high magnetic permeability [START_REF] Bossis | Yield stress in magnetorheological and electrorheological fluids: A comparison between microscopic and macroscopic structural models[END_REF]]. It is worth nothing that the affine model is just a particular case of the microscopic approach of the thermodynamic model, for which both force and energy-based calculations of the stress, Eqs.

(1) and ( 2) respectively, are completely equivalent. In the affine microscopic model, the interparticle forces are directly calculated either analytically [START_REF] Ginder | Rheology of magnetorheological fluids: Models and measurements[END_REF]] or numerically via multipolar approach [START_REF] Clercx | Many-body electrostatic interactions in electrorheological fluids[END_REF]; [START_REF] Klingenberg | Studies on the steady-shear behavior of electrorheological suspensions[END_REF]] or finite element method [Bossis et al. (2002); [START_REF] Ginder | Shear stresses in magnetorheological fluids: Role of magnetic saturation[END_REF]].

In the macroscopic approach of the thermodynamic model a mean field theory must be used, usually the Maxwell-Garnet one which does not require assumptions of any specific particle arrangement at the microscopic scale (one should only define the macroscopic features of the structure, i.e. columns, ellipsoids or stripes) and allows obtaining general results [START_REF] Bossis | Yield stress in magnetorheological and electrorheological fluids: A comparison between microscopic and macroscopic structural models[END_REF]]. As a consequence, the microscopic approach ensures a much better correspondence with experiments for MR suspensions composed of particles with high magnetic permeability, while the macroscopic one gives good predictions only for weakly magnetizable particles.

FIG. 1. Different model geometries used for stress calculations in MR fluids. The microscopic affine model [START_REF] Ginder | Rheology of magnetorheological fluids: Models and measurements[END_REF]] supposes affine displacement of particles with shear (a). Arrows denote displacement vectors of particles. The stress arises from the restoring interparticle force F x . If a given particle leaves its equilibrium position (dashed sphere in (a)), it will be immediately stacked on top of the closest particle, showing that the structure is mechanically unstable. The macroscopic thermodynamic model [START_REF] Bossis | Yield stress in magnetorheological and electrorheological fluids: A comparison between microscopic and macroscopic structural models[END_REF]] assumes formation of thick columns and ignores arrangement of particles within them (b). The stress comes from the restoring magnetic torque denoted by a bold arrow. A more realistic structure combines the features of both previous ones. Under strain, some particles experience affine motion, creating gaps between them, and the others are drawn into these gaps, maintaining the mechanical contacts with neighboring particles. One of the simplest structures corresponding to this picture is the BCT cluster shown in (c), which is found to be the most favorable for the energetic point of view [Tao and Sun (1992), [START_REF] Clercx | Many-body electrostatic interactions in electrorheological fluids[END_REF], [START_REF] Tao | Structural transformations of an electrorheological and magnetorheological fluid[END_REF]]. The stress response of such structure arises from both the longitudinal striction (due to the formation of gaps) and the restoring magnetic torque.

When these two models are applied to the system studied in the present work (suspension of strongly magnetizable particles at 0.50 volume fraction), it is found that the affine microscopic model gives the correct order of magnitude for the yield stress (360 Pa vs.

200 Pa experimental value-for the yield stress upon application of a field of 18.5 kA/m) note that we have used the approach of [START_REF] Ginder | Rheology of magnetorheological fluids: Models and measurements[END_REF] to estimate the interparticle forces. On the other hand, the macroscopic approach of the thermodynamic model strongly underestimates the yield stress (9 Pa vs. 200 Pa for a magnetic field of 18.5 kA/m). Thus, at a first view, we may be tempted to restrict our analysis to the affine model. However, it is clear that the rupture of chains by keeping an equal gap between neighboring particles as supposed in this model and depicted in Fig. 1a is not realistic see also comments on the caption of Figure 1. Furthermore, it is well known that in concentrated MR suspensions, thick columnar aggregates or other more isotropic structures are built upon field application, but single chains have never been observed [START_REF] Cutillas | A comparison between flow induced structures in electrorheological and magnetorheological fluids[END_REF]]. Therefore, from a realistic point of view, the thermodynamic model, which is based on general thermodynamic principles and assumes the formation of thick columns, is much more appropriate since it considers the correct size and shape of the aggregates. As already mentioned, the main reason for which the macroscopic model does not fit well to experimental results lies in the use of dipolar approach for interparticle interactions. This approach strongly underestimates the magnetic forces between particles and, as a consequence, gives too low values for the yield stress. In more details, the shear stress predicted by this model is given by the following relation:

2 0 22 () (1 ) H γ σµ µ µ γ ⊥ =- + E , ( 3 )
with µ 0 =4π•10 -7 H/m being the magnetic permeability of vacuum, H the magnetic field intensity inside the MR suspension, and µ E and µ ⊥ the diagonal components of the magnetic permeability tensor of the suspension with respect to a frame of reference with main axes in the directions parallel and perpendicular to the column aggregates. These last magnitudes will be called hereinafter longitudinal and transverse permeabilities for brevity. The yield stress is found as the maximum stress of the stress-strain dependency Eq. (3)-, which is reached at a critical strain, 1/ 3 crit γ =

. As seen in Eq. ( 3), the effects of the interparticle interactions in the suspension stress are taken into account through the values of these magnitudes. The inaccuracy of the dipolar approach can be corrected by including higher-order magnetic interactions in the suspension permeability. However as will be shown by our theory (see section III-B), the main contribution to the yield stress of suspensions of highly magnetizable particles does not come from the value of the permeability, but from the change of the energy with the formation of gaps and, therefore, from the derivative of the permeability with respect to the gap between the particles. Thus, Eq. ( 3) should include a term containing /

µ γ ∂∂ E and / µ γ ⊥ ∂∂ .
The most rigorous way to take into account the interparticle gap effect on the strength of a real structure is to conduct numerical simulations of the structure dynamics under slow shear flow. [START_REF] Bonnecaze | Yield stresses in electrorheological fluids[END_REF] simulated the structure of electrorheological fluids by the molecular dynamic method using multipolar interactions between particles, and calculated the electrostatic stress with the help of Eq. ( 2). However, they reported rather strong irregular oscillations of the stress response, as a result of the structure breakage and reformation, and the predicted values of the yield stress were too low.

In contrast to simulations, the microscopic models of [START_REF] Bossis | Yield stress in magnetorheological and electrorheological fluids: A comparison between microscopic and macroscopic structural models[END_REF], assuming more realistic multi-chain or zigzag-like clusters, provide a simpler picture of the interparticle gap effect on the suspension yield stress. The authors supposed that small interparticle gaps were formed under shear and calculated the stress via the derivative of the internal energy (Eq. ( 2)). The theoretical yield stress was found to strongly depend on the strain at which the interparticle gaps were formed. If the particles begin to separate from each other at zero strain (as it is the case of the affine deformation of multi-chain clusters), then the multipolar interactions predict a yield stress 2 0 6.5

Y µ H σ =
for the 50 vol.% MR suspension of the present work, which is an order of magnitude lower than the stress measured in our experiments. If the particles form zigzag-like clusters, they begin to separate at a non-zero critical strain that corresponds to the full extension of the clusters. This gives a several times increase of the yield stress as compared to affine models. However, at the critical strain, the stress experiences a jump of several orders of magnitude, inconsistent with experimental observations. Therefore, neither affine nor zigzag-like cluster models can be applied, in the form they are reported in literature, to highly concentrated MR suspensions.

In the next subsection we will develop a new model, which will be based on a realistic microscopic structure and will give correct values for the yield stress of MR suspensions in the concentrated regime. For this, we will take into account all the requirements for the particle structure that are drawn from the analysis performed above for previous models: (i) interparticle gaps must appear in the structure in order to obtain realistic values of the yield stress; (ii) the structure must be able to sustain a relatively large extension, in such a way that the critical strain would be high enough to give appropriate values of the yield stress; (iii) the failure of the structure must occur without any considerable jump in the stress. This implies a continuous formation of interparticle gaps during the shear strain of the clusters.

B. Theoretical model

Let us consider a simple shear deformation of a concentrated MR suspension confined between two parallel plates. An external magnetic field of intensity H 0 is applied perpendicularly to the plates. If the field is strong enough, it will induce percolating structures of magnetic particles. We will assume that these consist of body-centered tetragonal (BCT) structures, with four peripheral chains shifted vertically by a particle radius with respect to a central chain, and with neighboring particles being in contact at zero strain. This structure was found to be the most favorable from the energetic point of view [Tao and Sun (1992), [START_REF] Tao | Structural transformations of an electrorheological and magnetorheological fluid[END_REF]]. In addition, it was observed by laser diffraction experiments in electrorheological fluids electric counterparts of magnetorheological suspensions [START_REF] Chen | Laser diffraction determination of the crystalline structure of an electrorheological fluid[END_REF]]. Applied to our case of a highly concentrated suspension, this structure meets all the requirements mentioned at the end of the previous subsection, if we assume that the central chain deforms in an affine manner with the applied strain and the peripheral chains remain always in contact with the central one (see Fig. 1c). Since in our experiments we used surfaces with high rugosity, the structure is supposed to be always stuck to both plates, wall slip being absent, until its failure at the critical strain. Thus, when the upper plate is displaced a certain distance ∆x, the structure is strained by a magnitude γ =∆x/h, where h is the gap between the plates. We suppose that the structure deforms homogeneously until its failure and this assumption does not contradict to the condition of its mechanical stability: in the increasing branch of the stress-vs.-strain curve, any small perturbation of the homogeneous strain field should decay with the time. Under the strain, the particle structure turns along the vorticity axis and extends along its major axis (Fig. 1c), both effects contributing to the shear stress. On the one hand, the extension creates interparticle gaps and, consequently, restoring forces along the main axis of the structure. On the other hand, the structure rotation induces a restoring magnetic torque that tends to turn back the structure and align it with the magnetic field.

In order to calculate the stress response for a given stationary strain, γ, we must first obtain the magnetic permeability tensor of the suspension. Under homogeneous deformation, the main axes of magnetization of the suspension coincide with the axes of the BCT structures. Thus, in a reference frame linked to the particle structures, the magnetic permeability tensor is diagonal and has two dissimilar components: µ E along the major axis of the structures and µ ⊥ along their minor axis. From these, the magnetic permeability tensor with respect to the laboratory (rheometer) reference frame is obtained by rotation of the reference frame along the y-axis. The two components, zz µ and xz µ , which have importance in our calculations, are expressed as follows:

2 22 1 11 zz γ µµ µ γ γ ⊥ =+ ++ E , (4) 2 (
)

1 xz γ µµ µ γ ⊥ =- + E , (5)
µ E and µ ⊥ can be calculated as a function of the interparticle gap δ by solving Maxwell s equations by finite element method simulation. These calculations are reported in details in the Appendix, together with interpolation formulas for the dependencies of µ E and µ ⊥ with the relative gap, δ/a, with a being the particle radius. The relative gap δ/a is related to the strain by the following formula:

( )

2

/2 1 1 a δγ = +-. Thus, the magnetic permeability components can be expressed in terms of the strain γ.

Below the yield point, and under static straining conditions, the MR suspension may be considered as a magnetic, anisotropic elastic solid. For such continuum, the stress tensor is expressed, in its most general form, by the equation given by [START_REF] Landau | Electrodynamics of Continuous Media[END_REF]:

( ) , 1 2 ik ik i k k i ik T F FH B H B σδ γ ⎛⎞ ∂ =+ + + ⎜⎟ ∂ ⎝⎠ H # # , ( 6 )
where H and B are respectively the magnetic field intensity and the magnetic flux density inside the suspension, γ ik are the components of the strain tensor, δ ik those of the unit tensor, and F # is a thermodynamic function defined through the free energy of the suspension per unit volume, F, as follows:

0 0 FF F d ≡-⋅= - ⋅ ∫ H HB B H # , ( 7 )
with F 0 being the free energy per unit volume of the suspension in the absence of magnetic field. Note that [START_REF] Shkel | A continuum approach to electrorheology[END_REF] used similar expressions for electrorheological fluids in the small deformation limit for calculations of the dielectric tensor and the storage modulus. We shall exploit these formulas in a broader range of strains, which will allow us to calculate the yield stress of the suspension.

Magnetic particles used as solid phase in MR suspensions, usually show non-linear magnetization behavior, i.e. their magnetic permeability depends on the applied magnetic field strength. Similarly, the magnetic permeability of MR suspensions also presents a dependency with the applied magnetic field. However, in the relatively narrow range of magnetic field intensities used in our experiments, we may consider, within a good approximation, that the suspension permeability is field-independent. In this case, the magnetic flux density of the suspension presents a linear dependency with the magnetic field strength:

0 ii k k B H µ µ =
. By substitution of this into Eq. ( 7), the expression for the thermodynamic function F # reduces to:

2 22 2 2 00 0 00 0 22 11 11 1 22 2 1 2 1 FF H H F H H γ µµ µµ µµ µµ γ γ ⊥⊥ ⊥ =- - =- - + + EE E # , (8)
where cos HH θ = E and sin HHθ ⊥ = are the components of the internal magnetic field in the directions parallel and transverse to the BCT structures, respectively; θ is the strain angle (see Fig. 1b), which is related to the strain through the formula tan γ θ = . In the experimental case studied in this work, the internal magnetic field, H, in the thin layer of suspension confined between the two plates of the rheometer is related to the external applied field, H 0 , by the following expression: 2 00 2 ( 1)

zz HH H γ µµ µ γ ⊥ == + + E . ( 9
)
By substitution of Eqs. ( 5) and ( 8) into Eq. ( 6), we obtain the final expression for the shear stress (xz-component of the stress tensor) as a function of the applied strain and the magnetic field strength:

( ) ( ) 2 22 2 00 0 22 2 2 2 11 1 (1 ) 2 1 1 2 1 HH H µ µ γ γ γ σµ µµ µ µ µµ γ γγ γγ γ ⊥ ⊥ ⊥ ∂ ⎡⎤ ∂ =- - ⋅ + ⋅ + - ⎢⎥ + ∂+ ∂ + + ⎣⎦ E E E , ( 1 0 )
with the expression for / µ γ ∂∂ E given in the Appendix. The first term in Eq. ( 10) is the elastic contribution due to the restoring magnetic torque acting on the tilted BCT structures. This term is also connected to the variation of the suspension magnetic permeability with the rotation of the structures. The second term stands for the elastic contribution due to the restoring striction forces that tend to compress the structures extended by the shear. In other words, these forces tend to reduce the interparticle gaps induced by the strain. From the macroscopic point of view, this term is related to the variation of the magnetic permeability with the extension of the aggregates. The last term is the Maxwell stress, which arises from the deformation of the structures under the applied magnetic field. On the other hand, the first and the second terms come from the effect of deformation on the magnetic properties of the structures. Analysis shows that the first and the last terms of Eq. ( 10) are negligible as compared to the second one. This confirms the extreme importance of the interparticle gaps as a principle source of the stress of MR suspensions upon applied magnetic field.

IV. RESULTS

The strain dependency of both µ E and µ ⊥ is shown in Fig. 2 for the suspension under study in the present work. In the range of the strains, 0 < γ < 0.5, the longitudinal permeability, µ E , shows an important decrease with the strain. This decrease is connected with the appearance and enlargement of gaps between the particles of a same chain when the BCT structures are extended by the shear (see Fig. 1c). The size of these gaps increases proportionally to 2 γ , and this provokes a quite strong variation of µ E . On the other hand, the transverse permeability, µ ⊥ , appears to be independent of the strain. This can be easily understood since µ ⊥ depends mainly on the distance between particles of two opposite peripheral chains (dashed circles in the top view of Fig. 1c). The initial separation between these particles is equal to (3 1 )( 2) 1 . 5 aa -⋅ 0 . When the BCT structure is extended, there is an approach of these particles, but at small enough strains it is almost negligible with respect to the initial separation and, thus, it does not affect appreciably the value of µ ⊥ .

FIG. 3. Theoretical stress-strain curve for a suspension consisting of Fe-CC particles dispersed in mineral oil (volume concentration 50 %). The applied magnetic field is H 0 =18.5 kA/m. Fig. 3 shows the stress-strain dependency calculated with the help of Eq. ( 10) for the MR suspension under study (50 vol.% of Fe-CC in mineral oil) upon application of a field H 0 =18.5 kA/m. As expected, when the strain increases, the stress first increases, then reaches a maximum, and finally decreases gradually towards nearly zero values at γ >0.3. Note that the maximum of the stress takes place at a critical strain, γ cr ≈0.115, much lower than that predicted by the macroscopic model of [START_REF] Bossis | Yield stress in magnetorheological and electrorheological fluids: A comparison between microscopic and macroscopic structural models[END_REF]:

1/ 3 0.58 cr γ =≈ . This is not surprising because, in our model, the longitudinal component of the suspension permeability falls with the strain quite rapidly (see Fig. 2), provoking a rapid increase of the stress at low strains. On the other hand, in the macroscopic model of [START_REF] Bossis | Yield stress in magnetorheological and electrorheological fluids: A comparison between microscopic and macroscopic structural models[END_REF], the permeability µ zz decreases slowly with the strain, the structure rotation with respect to the applied field being the only cause for the decrease in this model. Nevertheless, the critical strain γ cr ≈0.115

predicted for the BCT structure by the present model is two-three times larger than the one calculated for single chains see Fig. 7 in the work by Bossis et al. (2002). Consequently, the critical strain obtained in the present work should be large enough to ensure reasonable predictions of the yield stress. As stated above, structures become unstable at strains γ >γ cr and are thus supposed to break at the summit of the stress-strain curve. So, the yield stress corresponds to the maximum of the stress-strain curve and is calculated by replacing the strain γ by the critical strain γ cr ≈0.115 in Eq. ( 10 . We refer to this quantity as yield stress increment. The yield stress at zero field is likely due to the remnant magnetization of the iron particles, as well as to non-magnetic colloidal interactions between them, mainly van der Waals forces, which may induce the formation of particle aggregates in the absence of applied field. The importance of these forces in suspensions of magnetic particles has been studied by different authors see for example [START_REF] Phulé | The role of the dispersed-phase remnant magnetization on the redispersibility of magnetorheological fluids[END_REF]].

As seen in Fig. 4, the yield stress increment increases gradually with the magnetic field intensity. The experimental dependence of the yield stress increment is approximately quadratic, as proved by the fact that the best fit to a power law gives an exponent of 1.91 ± 0.07 (dashed line in Fig. 4). The increment of the yield stress with the applied field is commonly explained by the fact that, at higher fields, stronger magnetic forces act between particles and, thus, stronger forces are required to separate them and to break the structures.

Concerning the predictions of our theoretical model, as observed in Fig. 4, the theory (without adjustable parameters) still underestimates the experimental data by about 50%. One of the reasons for such a discrepancy could be a polydispersity of the MR suspension. As stated by [START_REF] Kittipoomwong | Dynamic yield stress enhancement in bidisperse magnetorheological fluids[END_REF], small magnetic particles may form locally more compact aggregates and cause the larger particles to form more robust aggregates. Such aggregates resist better to the applied strain and give a higher yield stress compared to that of the monodisperse MR fluid at the same volume fraction of particles.

V. CONCLUSIONS

We have presented in this work reliable experimental data of the magnetic fieldinduced static yield stress of a MR suspension with concentration near the limit of maximumpacking fraction. Both the preparation of the MR suspension and the experimental protocol followed for the measurement of the static yield stress were carefully chosen to guarantee reproducibility of the results. In addition, we have developed a model for the static yield stress of highly concentrated MR suspensions upon magnetic fields. This model is based on the key hypothesis that interparticle gaps inside the field-induced particle structures must appear when they are subjected to shear strain, at the same time that their mechanical stability must be maintained. This perspective supposes an original contribution with respect to standard models based on the affine deformation of single chains, which seems unrealistic for highly concentrated MR suspensions. The key hypothesis taken as starting point in our model is easily fulfilled by considering a structural unit consisting of five chains of particles located at positions of a BCT structure. When this structure is sheared the particles of the central chain are supposed to move in an affine way, whereas the particles of the peripheral chains remain in contact with those of the central chain ensuring the mechanical stability of the structure.

Estimations of our model show that the main contribution to the stress of highly concentrated MR suspensions comes from the change in the magnetic permeability of these, as interparticle gaps are forced and enlarged by the shear stress. Therefore, the neglect of this particular aspect in most of the existing macroscopic models is the likely reason for their failure when applied to highly concentrated suspensions of strongly magnetizable particles. This last statement is supported by the quite good agreement obtained between predictions of our model and experimental results, taking into account that our model does not include any adjustable parameter. The theoretical model developed in this work will be, in the near future, extended to the steady state shear flow of MR suspensions taking into account eventual flow instabilities in highly concentrated regime. Finally, the simple interpolation formulas issued from the present model could be useful for engineering calculations.

APPENDIX: PERMEABILITY CALCULATIONS

We calculated the longitudinal and transverse magnetic permeabilities of the MR suspension by solving Maxwell s equations by means of finite element method (FEM) simulations. These simulations were performed with the help of the free software FEMM [Meeker (2009)]. For this aim, and although the BCT structures of our model are threedimensional (3D), we considered a planar problem since 3D-problems cannot be implemented with FEMM software. Nevertheless, we expect that the simplification of the real 3D-problem to a planar-one will not be a restriction for the validity of the results obtained in this appendix.

The representative planar cells used for the calculations of both components of the magnetic permeability tensor, µ E and µ ⊥ , are shown in Figs. A1 andA2, respectively. The vertical axis of symmetry of the first cell coincides with the major axis of the BCT aggregate.

The two closely spaced semi-circles represent particles belonging to the central chain and the two lateral circles represent particles of the peripheral chains. The whole BCT structure is relatively long (aspect ratio /(2 3) 40 ha ≈ ) and can be considered as an infinite stack of the unit cells. The horizontal dimension d of the first cell (Fig. A1) is chosen to be equal to the mean distance between BCT structures in their hexagonal arrangement: For the calculation of the transverse magnetic permeability, µ ⊥ , the external magnetic field must be applied perpendicularly to the BCT structure. However, for the convenience of definition of the boundary conditions, we can still reduce our problem to longitudinal magnetic fields. For this aim, a 90º rotation of the internal structure shown in Fig. A1 is enough. By doing it we obtain the unit cell used for the calculations of µ ⊥ (Fig. A2). In this /2 1 1 a δγ =+ -, we obtain the longitudinal magnetic permeability, as a function of the strain γ. Finally, the derivative of µ E with respect to the strain, which appears in Eq. ( 10) for the stress, is given by the following expression:

( )

FIG. 2 .

 2 FIG.2. Longitudinal, µ E , and transverse, µ ⊥ , components of the magnetic permeability tensor of a suspension consisting of Fe-CC particles dispersed in mineral oil (volume concentration 50 %) as function of the applied shear strain.

  FIG. 4. Experimental and theoretical dependencies of the yield stress increment,

  a is the particle radius and Φ=0.5 the volume fraction of particles in the suspension. In this figure, the applied magnetic field is oriented vertically and the intensity H of the internal magnetic field (averaged over the cell volume) is imposed on both lateral borders of the cell. Periodic boundary conditions are used for the upper and the lower borders of the cell.

FIG. A1 .

 A1 FIG. A1. Elementary cell of the MR suspension used for the calculation of the longitudinal magnetic permeability, µ E . The surface plot of the magnetic flux density and the magnetic field lines are shown. The red spots correspond to the regions of high magnetic flux density in the vicinity of the contact points between spheres.

  figure, the two closely spaced circles represent particles of the central chain and the two semi-

ACKNOWLEDGEMENTS

We would like to thank Professor Andrey Zubarev for helpful discussions. Financial support by projects Fis2009-07321 (Ministerio de Ciencia e Innovación, Spain), P08-FQM-3993 and P09-FQM-4787 (Junta de Andalucía, Spain), Biomag (PACA, France) and Factories of the Future (Grant No. 260073, DynExpert FP7) is gratefully acknowledged. MTL-L also acknowledges financial support by Universidad de Granada (Spain). LR-A acknowledges financial support by MECD (Spain) through its FPU program.

The non-linear magnetic properties of carbonyl iron particles are well described by the Fröhlich-Kennelly law [START_REF] Jiles | Introduction to Magnetism and Magnetic Materials[END_REF]]:

(1 ) 1

(1 )

where µ p is the magnetic permeability of the particles, µ i = 250 and M S = 1.36•10 6 A/m are, respectively, the initial permeability and the saturation magnetization of the particles, and H p is the magnetic field intensity inside the particles.

Once the magnetic field distribution is found for both cells (Figs. A1 andA2), the corresponding permeability component is calculated by the following formula:

where H is the intensity of the mean magnetic field inside the MR suspension, imposed on the lateral borders of the cell and related to the external field, H 0 , via Eq. ( 9). The integration in Eq. (A.2) is performed over the whole cell volume V. Recall that the simulations are performed considering planar geometries, and this means that the circles of Figs. A1 andA2 are taken, from the point of view of the FEM simulation, as infinite cylinders oriented perpendicular to the plane of the page. In order to evaluate the numerical error related to this, we calculated by FEM simulation the longitudinal magnetic permeability of an infinite chain composed of touching spheres, and that of an infinite chain composed of touching cylinders with their axes perpendicular to the applied magnetic field. Simulations showed that the magnetic permeability of the chain of spheres is 3.43 times smaller than that of the chain of cylinders. Consequently, all the results obtained for the multi-chain geometry of Figs. A1 and A2, including those of Eq. (A.2), were reduced by the factor 3.43. This is perhaps a rough approximation but its validity is supported by the fact that it allows obtaining a good


Franc ¸ois Keith

Nicolas Mansard

Sylvain Miossec

Abderrahmane Kheddar

From Discrete Mission Schedule to Continuous Implicit Trajectory using Optimal Time Warping

paper presents a generic solution to apply a mission described by a sequence of task on a robot while accounting for its physical constraints, without computing explicitly a reference trajectory. A naive solution to this problem would be to schedule the execution of the tasks sequentially, avoiding concurrency. This solution does not exploit fully the robot capabilities such as redundancy and have poor performance in terms of execution time or energy. Our contribution is to determine the time-optimal realization of the mission taking into account robotic constraints that may be as complex as collision avoidance. Our approach achieves more than a simple scheduling; its originality lies in maintaining the task approach in the formulated optimization of the task sequencing problem. This theory is exemplified through a complete experiment on the real HRP-2 robot.

Introduction

A robot is designed to perform missions in various application contexts. When the environment is well or partially structured most missions can be hierarchically decomposed into a set of tasks (i.e. generic sensory-motor functions) which has to be mapped into robot execution. Numerous works have been proposed to compute such a sequence of tasks from a given mission and a set of causal paradigms [START_REF] Dechter | Constraint Processing[END_REF][START_REF] Ghallab | Automated Plan-Theory and Practice[END_REF]. However, they generally produce a symbolic plan, where the only numerical precisions lie on the scheduled time data. Moreover, constraints have to be expressed under a symbolic expression. Its robotic application into the real world requires the time sequence to be refined, typically through an applicative path planner [START_REF] Lavalle | Planning Algorithms[END_REF], that will compute the trajectories to be followed by the robot [START_REF] Lamare | Integrating a temporal planner with a path planner for a mobile robot[END_REF]). Yet, the meaning of the symbolic plan is lost in the global trajectory. Such low-level methods lack of robustness to environment changes or uncertainties. Consequently, the remaining trajectory may have to be recomputed several times while the mission is being achieved. Moreover, it is difficult (and often specifically hard coded) to enhance the trajectory with symbolic data, that would help re-computing only part of the plan (Py and Ingrand Copyright c 2009, Association for the Advancement of Artificial Intelligence (www.aaai.org). All rights reserved.

2004

) or distort locally the trajectory after small environment changes (Quinlan and [START_REF] Khatib | Elastic bands: Connecting path planning and robot control[END_REF].

Rather than using a trajectory planner between the temporal reasoning and its real robotic execution, we propose to use a sensory-motor control approach based on task components. The task function [START_REF] Samson | Robot Control: the Task Function Approach[END_REF] is an elegant approach to produce intuitively sensorbased robot objectives. Based on the redundancy of the system, the approach can be extended to consider a hierarchical set of tasks [START_REF] Siciliano | A general framework for managing multiple tasks in highly redundant robotic systems[END_REF]. Hierarchy of tasks are becoming popular to build complex behavior for very redundant robot such as humanoids [START_REF] Mansard | Task sequencing for sensor-based control[END_REF][START_REF] Sentis | A whole-body control framework for humanoids operating in human environments[END_REF].

A task (i.e. a task function) can be directly linked to the symbols on which the task temporal network is reasoning. Mission decomposition is thus executable directly using the sensory-motor mapping of the task function. However, exclusive task sequencing on the robot produces generally jerky suboptimal movements which may look to humans as monotonous automated motions. This paper focuses on finding a solution to produce automatically an optimal plane/schedule that makes use of the redundancy by enabling task concurrency. It seems difficult to use temporal networks to produce a scheduling with task overlapping when the tasks concurrency is restricted by physical limitations of the robot (for example obstacles or balance of a biped robot), since the constraints are not in a discreet form. On the other hand, semi-infinite optimization [START_REF] Lee | Newton-type algorithms for dynamics-based robot movement optimization[END_REF]) is known to generate low level trajectories, while accounting for such constraints, but with insufficient robustness to environment uncertainties.

In this paper, we propose to rely on task for both the symbolic reasoning and control of the robot. In between, we propose to use semi-infinite optimization to refine the symbolic schedule and account for system constraints. Given a sequence of tasks to achieve a mission, our solution returns for each task the optimal times at which it is activated and inactivated and the optimal parameters for the task execution. The originality of our approach lies in keeping the task component in the formulation of this problem, which can roughly translate to optimizing tasks overlapping by manipulating tasks, i.e. the controllers as variables of the optimization problem.

Generic Task Sequencing

Task function formalism and Stack of Tasks

Defining the motion of the robot in terms of task simply consists in choosing several control laws to be applied on a subpart of the robot degrees of freedom (DOF). A task is defined by a vector e (typically, the error between a signal and its desired value, e = s-s * ). The Jacobian of the task is noted J = ∂e ∂q , where q is the robot configuration vector. In the following, we consider that the robot input control is the joint velocity q: ė = J q. Considering a reference behavior ė * to be executed in the task space, typically,

ė * = -λe, (1)
the control law to be applied on the robot whole body is given by the least-square solution:

q = J + ė * + Pz (2)
where J + is the least-square inverse of J, P = I -J + J is the null-space of J and z is any secondary criterion. P ensures a decoupling of the task with respect to z, which can be extended recursively to a set of n tasks, each new task being fulfilled without disturbing the previous ones:

qi = qi-1 + (J i P A i-1 ) + ( ėi -J i qi-1 ), i = 1 . . . n (3)
where q0 = 0 and P A i is the projector onto the null-space of the augmented Jacobian J A i = (J 1 , . . . J i ). The robot joint velocity realizing all the tasks in the stack is q = qn . A complete implementation of this approach is proposed in [START_REF] Mansard | Task sequencing for sensor-based control[END_REF] under the name Stack of Tasks (SoT). The structure enables to easily add or remove a task, or to swap the priority order between two tasks, during the control. Constraints (such as joints limit) can be locally taken into account. The continuity of the control law is preserved even at the instant of change.

Gain handling

The simple attractor presented in (1) produces a nice exponential decrease. However, it also produces a strong acceleration at the beginning of the task, while at the end of the task, e decreases slowly (slow convergence). A very classical 'trick' when regulating a task is to rather use an adaptive gain λ = λ(e(t)), for example:

λ(e) = (λ F -λ I ) exp -e β + λ I (4)
with λ I the gain at infinity, λ F the gain at regulation (λ I ≤ λ F ) and β the slope at regulation.

Sequence of tasks

A task sequence is a finite set of tasks sorted by order of realization, and eventually linked to each other. Any pair of tasks can be either independent (i.e. they can be achieved in parallel), or constrained (i.e. one may have to wait for another one to be achieved, in order to make sense or to be doable). The sequence can be formulated into a classical temporal network scheduling, starting at t 0 and ending at t End . Both values are finite and the sequence does not loop.

Besides, we may consider for the sake of clarity but without loss of generality that each task appears only once.

The position of a task in the sequence is defined by the time interval during which it is maintained in the SoT. For a given task i, this interval is noted [t I i , t F i ], defined with respect to t 0 . The regulation time t R i is defined by e i (t R i ) = ǫ i , with ǫ i the tolerance upon task realization.

A task sequence is characterized by a set of timeconstraints binding the schedules of two tasks e i and e j . They can be defined as follow1 : e i must begin or end once e j has begun, has ended or has been regulated. We use the following notation to describe the sets of pairs of tasks e i and e j that undergo these dependencies (e i is the direct predecessor of e j ) :

S I,I = {(e i , e j ) | t I i ≤ t I j } S R,I = {(e i , e j ) | t R i ≤ t I j } S R,F = {(e i , e j ) | t R i ≤ t F j } S F,I = {(e i , e j ) | t F i ≤ t I j } S F,F = {(e i , e j ) | t F i ≤ t F j } (5)

Continuous optimization of sequence of tasks

Given a set of hypothesis described using (5), a generic solution is now proposed to automatically compute the optimal task and schedule parameters to be executed by the SoT.

General problem formulation

An optimization problem is composed of a criterion to minimize, along with a set of variables and equality and inequality constraints to be satisfied. Our criterion is the total duration of the mission. The variables are: (i) the times t I i and t F j and (ii) the gains (λ I , λ F , β). The general optimization problem is written as follows:

min x t End (6a) subject to q = SoT x (q, t) (6b) seq(q) < 0 (6c) φ(q) < 0 (6d) ∀i, t F i ≤ t End (6e) with x = [t I 1 , t F 1 , λ I 1 , λ F 1 , β 1 , . . . , t I n , t F n , λ I n , λ F n , β n , t
End ] set of optimization variables of each task, t End the duration of the mission, and seq(q) and φ(q) the sequencing and robotic constraints.

Constraint definitions

Parameters x satisfy both the sequencing and the robotic time-constraints enumerated hereafter:

Tasks constraints, noted seq(q) are defined as follow:

• Time coherence for each task i, that is:

∀i, 0 ≤ t I i < t F i ≤ t End (7)
• Termination condition for each task i, that is:

∀i, s * i -s i (t F i ) < ǫ i (8)
• Task sequence conditions for any couple of tasks i and j:

∀(i, j) ∈ S I,I , t I i ≤ t I (9) ∀(i, j) ∈ S F,F , t F i ≤ t F j (10) ∀(i, j) ∈ S F,I , t F i ≤ t I j (11) ∀(i, j) ∈ S R,I , t R i ≤ t I j (12) ∀(i, j) ∈ S R,F , t R i ≤ t F j ( 13
)
• Gain consistency for each task i, namely:

∀i, λ I i ≤ λ F i ( 14
)
The constraints ( 7), ( 9), ( 10), ( 11) and ( 14) are linear. On the contrary, the constraint ( 8) is impossible to compute directly using x, and is determined from a simulation of the execution. Care has to be taken while resolving the condition described by ( 12) and ( 13). Indeed, discretizing t R to the closest simulation step will produce discontinuities which may disturb the optimization process. The constraint ( 12) and ( 13) are thus equivalently rewritten:

∀(i, j) ∈ S R,I , s * i -s i (t I j ) ≤ ǫ i (15) ∀(i, j) ∈ S R,F , s * i -s i (t F j ) ≤ ǫ i ( 16
)
Robot constraints : φ(q) Those constraints are mainly due to hardware intrinsic limitations of the robot:

• Joint limits, given by q min ≤ q ≤ q max (17) q min , q max are the lower and upper joint limits.

• Velocity limits, given by qmin ≤ q ≤ qmax (18) qmin , qmax are the minimal and maximal velocity.

• Collision between a given pair of objects i and j, given by

d ij ≥ 0 ( 19
)
d ij is the distance between objects i and j. Object designate those found in the mission environments and each link of the robotic system. Hence, both collision with the environment and self-collision of the robot have to be evaluated.

It can be shown that (6) defines a continuous optimization problem. However, q is in fact a vector of functions of time, hence constraints φ(q) are semi-infinite, i.e. taking place for all the values of the continuous variable t ∈ [t 0 , t End ]. Cares have thus to be taken when implementing them in classical solver softwares.

Implementation Optimization

At each optimization step, the solver chooses a new set of parameters x. It then computes the constraints. Constraints (7), ( 9), ( 10) and ( 14) can be evaluated directly. As stated previously, the other constraints can not be directly computed (since they do not write in an analytical explicit formulation). They are thus evaluated using a complete simulation of their execution. The chosen value of the current optimization variable vector x is transmitted by the solver to a simulation engine. The simulation returns the evaluation of the constraints and the optimization solver computes a new step vector x, until convergence.

Simulation

The simulation is basically a numerical integration of ( 6) using an explicit Euler integration method with a fixed step ∆t = 0.005sec. The times t I i and t F i are continuous variables that may not be aligned with the integration grid, and have to be explicitly added to ensure continuity with respect to the optimization variables. The simulation engine runs under the AMELIF framework [START_REF] Evrard | Framework for haptic interaction with virtual avatars[END_REF], an interactive dynamic simulator for virtual avatars which includes collision detection and task handling according to the SoT formalism. The execution for both simulation and real-robot control is performed by a generic control framework based on [START_REF] Mansard | A versatile generalized inverted kinematics implementation for collaborative humanoid robots: The stack of tasks[END_REF].

Experiment Temporal network

The experiments have been realize using a pick-in-box scenario: the robot has to open the fridge with one hand before grasping the object with the other hand. Similarly, it has to wait for the grasping to be completed before closing the fridge. Instead of explicitly defining these causal constraints, we simply consider the collision detection to ensure the causality. The optimal schedule should use the system redundancy by overlapping the tasks on each arm.

The task sequence is composed of 10 tasks (see Fig. This mission can not be split into smaller independent sequences. The constraints considered for this problem are thus sequencing and robotic constraints (joint position and velocity limits), and non-colliding with the fridge.

Results of the optimization

The optimization ran on a 3GHz desktop PC running under Windows OS, using MATLAB SQP solver. The sequence found is described on Fig. 2. The overlaps between the tasks of each arm appear clearly: the left arm starts moving before the fridge is open. It then starts the reaching motion even if the fridge is not completely open. The right arm starts to when the task is added in the SoT, its error is first regulated (this corresponds to the dark part (red or dark blue) of the block). From t R i , the error is nearly null and the task is kept in the SoT (light part (yellow or cyan) of the block) until t F i .

close the fridge before the left arm has completely left the fridge area. The whole task sequence lasts 47sec. Without these two overlaps, the optimal timing will only be 71sec. The SoT formalism allows to directly apply the optimized task sequence as a control law using the same task definition. The task sequence is executed by the robot HRP-2, a full-size humanoid robot with 30 actuated DOF. For the tasks requiring a haptic interaction (i.e. opening and closing the fridge) the force sensor of the robot is used to close the loop and compensate for position uncertainties. The robot manages to grasp the can without constraint violation, in a smooth manner thanks to the optimized gain. Snapshots of the execution are given in Fig. 3. See http://www. laas.fr/ ˜fkeith/keithfridge.avi for the video of the complete execution.

Conclusion

We devise a method which allows to automatically pass from a symbolic plan to a complete motion generator that takes into account all kind of robot constraints, such as joint limits or collision. The result optimize both the behavior and the overlapping scheduling of a sequence of tasks composing a robotic mission. The solution derives from an optimization formulation of the tasks scheduling keeping the formalism built on the top of a task-function based control. The solution is generic, and may applied to any problem of motion generation based on a schedule of tasks under continuous semi-infinite constraints, while optimizing various criteria. For the time being, our method still needs a predefined ordered sequence. As a future work the autonomy will be improved by determining automatically the ordered sequence and compute all the necessary subtasks from definitions of actions/objects associations.

  1): e 0 -open the right gripper ; e 1 -move the right arm to the fridge ; e 2 -close the right gripper ; e 3 -open the fridge ; e 4 Close the fridge ; e 5 -open the left gripper ; e 6 -move the left gripper in the fridge area ; e 7 -move the left gripper to the can ; e 8 -close the left gripper ; e 9 -lift the can ; e 10 -remove the can out of the fridge.

Figure 1 :

 1 Figure 1: Sequence of tasks for taking the can in the fridge

Figure 3 :

 3 Figure 3: HRP-2 grasping a can in the fridge.

contrary to Allen Logic, that only considers the start and end date, here is also considered the regulation time t R


François Keith
email: keith@lirmm.fr

Nicolas Mansard
email: nmansard@laas.fr

Sylvain Miossec
email: sylvain.miossec@bourges.univ-orleans.fr

Abderrahmane Kheddar
email: kheddar@lirmm.fr

Optimized Time-Warping Tasks Scheduling for Smooth Sequencing

Keywords: Scheduling algorithms, Task-based control, Optimization, Robotics

This paper presents an optimal formulation for sequencing a set of robotic tasks. Tasks description and execution are based on the task formalism. A naive solution would be to arrange the execution of the tasks sequentially (potential redundancy not exploited, not optimal), or to set manually the timing and the behavior for each task (burdensome for the user, lack of autonomy, not optimal), or to arrange the tasks by priority and pill them up at once in the stack-of-tasks (high likelihood of conflicts, wrong or cumbersome execution). Our contribution is to determine the time-optimal realization of the mission taking into account robotic constraints as complex as collision avoidance; its originality lies in keeping the task formalism in the formulated optimization of the task sequencing problem. This theory is exemplified through a simulation on a 2D robot.

INTRODUCTION

A robot is designed to perform missions in various application contexts. When the environment is well or partially structured most missions can be hierarchically decomposed. That is, missions undergo functional objective decomposition into a set of processes or operations that can be defined as templates. Each operation can be decomposed into a set of tasks (i.e. generic sensorymotor functions), and each task can be easily mapped into robot execution. The whole scheme may constitute an exploitable generic skill/behavior. Yet, various levels of decomposition can be achieved depending on the envisaged software/hardware implementation, additional environment constraints, the human-machine interface, etc. In the end, the robot is assigned with a sequence of tasks to realize a given mission. A typical example of such a mission decomposition is given in Fig. 1.

Numerous works have been proposed to obtain such a sequence of tasks from a given mission and a set of causal paradigms [START_REF] Dechter | Constraint Processing, chapter 12, Temporal Constraint Network[END_REF]; [START_REF] Ghallab | Automated Planning: Theory and Practice[END_REF]). However, they generally produce a symbolic plan, where the only numerical precisions lie on the scheduled time data. Its robotic application into the real world requires the time sequence to be refined, typically through an applicative path planner [START_REF] Lavalle | Planning Algorithms[END_REF]), that will compute the trajectories to be followed by the robot [START_REF] Lamare | Integrating a temporal planner with a path planner for a mobile robot[END_REF]). Yet, the meaning of the symbolic plan Fig. 1. Typical schedule plan for the execution of a given mission: in order to grasp the green object placed behind the planar door (the yellow line segment), the mission can be divided in six steps: (i) reach and grab the door's handle, (ii) open the door, (iii) reach the object, (iv) grasp it, (v) lift the object out and (vi) close the door. For this example, one link of the stickrobot is not allowed to move.

is lost in the global trajectory. Such low-level methods lack of robustness to environment changes or uncertainties. Consequently, the remaining trajectory may have to be recomputed several times while the mission is being achieved. Moreover, it is difficult (and then often specifically hard coded) to enhance the numerical trajectory with symbolic data, that would help re-computing only part of the plan [START_REF] Py | Dependable execution control for autonomous robots[END_REF]) or distort locally the trajectory to apprehend small environment changes [START_REF] Quinlan | Elastic bands: Connecting path planning and robot control[END_REF]; [START_REF] Lamiraux | Reactive path deformation for nonholonomic mobile robots[END_REF]).

Rather than using a trajectory planner between the temporal reasoning and its robotic real execution, we suggest using a sensory-motor control approach based on task components. The task function [START_REF] Samson | Robot Control: the Task Function Approach[END_REF]) or the operational space formulation [START_REF] Khatib | A unified approach for motion and force control of robot manipulators: The operational space formulation[END_REF]) are elegant approaches to produce intuitively robot objectives. They also allow to address the control problem directly in the sensor space, improving robustness of the action execution against environment uncertainties and variability [START_REF] Espiau | A new approach to visual servoing in robotics[END_REF]; [START_REF] Chaumette | Visual servo control, pt i: Basic approaches[END_REF]). They are trajectory free, which means that it is not necessary to explicitly compute all trajectories before the execution or during the execution, namely in response to environment changes. Moreover, since a same task space is valid for a large set of robots, a control scheme based on task formalism is certainly portable and easy to modify and to maintain. In addition, these methods include a kinematics or dynamic formulation that decouples the task space from the null-space (i.e. the joint space that let the task invariant) [START_REF] Liégeois | Automatic supervisory control of the configuration and behavior of multibody mechanisms[END_REF]; [START_REF] Hanafusa | Analysis and control of articulated robot with redundancy[END_REF]).

A secondary task can then be applied in the null-space, and, recursively, a hierarchic set of tasks (or stack-of-tasks) can be considered [START_REF] Nakamura | Taskpriority based redundancy control of robot manipulators[END_REF]; [START_REF] Siciliano | A general framework for managing multiple tasks in highly redundant robotic systems[END_REF]). Hierarchy of tasks are becoming popular to build complex behavior for very redundant robot such as humanoids [START_REF] Baerlocher | An inverse kinematic architecture enforcing an arbitrary number of strict priority levels[END_REF]; [START_REF] Bolder | Visually guided whole body interaction[END_REF]; Mansard et al. (2007); [START_REF] Sian | A switching command-based whole-body operation method for humanoid robots[END_REF]; [START_REF] Sentis | A whole-body control framework for humanoids operating in human environments[END_REF]). The formalism introduced in (Mansard and Chaumette ( 2007)) proved to be efficient in dealing with complex humanoid missions: the Stack of Tasks (SoT) is mainly a hierarchy of tasks driving the robot while ensuring locally a given set of constraints to be satisfied. We make use of this formalism (Section 2).

A task (i.e. a task function [START_REF] Samson | Robot Control: the Task Function Approach[END_REF])) can be directly linked to the symbols on which the task temporal network is reasoning (e.g. reaching an object to be grasped is a task that requires the robot arm to be available, and whose post-condition is to have the gripper on the object -it is also directly described by a sensory-motor function applicable to the SoT).

Mission decomposition is thus executable directly by a SoT, which guarantees good robustness and avoid unnecessary trajectory (re)computation. However, exclusive task sequencing on the robot produces generally saccadic movements which may look to humans as monotonous automated motions. On the other hand, it is difficult for the temporal network to produce a scheduling with task overlapping when tasks' concurrency is restricted by physical limitations of the robot (for example obstacles or dynamical constraints on a humanoid). Since the problem is not in a discreet form -i.e. the resources (mainly redundancy in motion generation) are not discreet-, task scheduling techniques can not apply directly to our problem. On the contrary, semi-infinite optimization techniques [START_REF] Lee | Newton-type algorithms for dynamics-based robot movement optimization[END_REF], [START_REF] Miossec | Development of a software for motion optimization of robots-application to the kick motion of the HRP-2 robot[END_REF]) can be used to generate low level trajectories for the overall execution, providing the problem to be formulated this way. One of the drawbacks of such an approach is that the generated trajectory does not necessarily account for the robot controller and is not robust to environment uncertainties: we propose a method that encapsulates the tasks concept to solve this issue.

The overall architecture, including the optimized task sequencing, is illustrated in Fig. 2. Given a set of elementary tasks sequence to achieve a given mission, our solution re-Fig. 2. Schematic representation of components of a robotic mission execution. turns for each task, the optimal times at which it is put and removed from the SoT and also the optimal parameters for the task execution (e.g. control gain). We additionally expect from this method a smooth tasks sequencing (i.e. smooth transitions of tasks through task overlapping). The originality of our approach lies in keeping the task component in the optimization formulation of this problem, which can roughly translate to optimizing tasks overlapping by manipulating tasks, i.e. the controllers, as variables of the optimization problem. The task formulation details are first recalled in Section 2. The generic solution is then detailed in Section 3. The theory is finally exemplified through a simulation on a 2D robot, the mission consisting in getting an object placed behind a planar door.

GENERIC TASK SEQUENCING

Task function formalism and Stack of Tasks

Defining the motion of the robot in terms of task simply consists in choosing several control laws to be applied on a subpart of the robot degrees of freedom (DOF).

A task is defined by a vector e (typically, the error between a signal s and its desired value, e = ss * ). The Jacobian of the task is noted J = ∂e ∂q , where q is the robot configuration vector. In the following, we consider that the robot input control is the joint velocity q. The equation of motion is thus reduced to the kinematics: ė = J q (1) Considering a reference behavior ė * to be applied in the task space, the control law to be applied on the robot whole body is given by the least-square solution: q = J + ė * + Pz (2) where J + is the least-square inverse (Ben-Israel and Greville ( 2003)) of J, P = I -J + J is the null-space of J and z is any secondary criterion that will be applied without disturbing the main task thanks to the projection into P. A typical requested behavior is the regulation of the error, which can be obtained though an exponential decrease by setting:

ė * = -λe (3)
As mentioned earlier, (2) enables to compose a complex behavior from a set of tasks [START_REF] Siciliano | A general framework for managing multiple tasks in highly redundant robotic systems[END_REF]; [START_REF] Baerlocher | An inverse kinematic architecture enforcing an arbitrary number of strict priority levels[END_REF]; [START_REF] Sentis | A whole-body control framework for humanoids operating in human environments[END_REF]): z can be used to fulfil a secondary task, without disturbing the main task having priority. This nice decoupling can be extended recursively to a set of n tasks, each new task being fulfilled without disturbing the previous ones:

qi = qi-1 + (J i P A i-1 ) + ( ėi -J i qi-1 ), i = 1 . . . n (
4) where q0 = 0, P A i is the projector onto the null-space of the augmented Jacobian J A i = (J 1 , . . . , J i ) and J i =

J i P A i-1
is the limited Jacobian of the task i. The robot articular velocity realizing all the tasks in the stack is q = qn . A similar recursive formulation can be obtained to compute the P A i , [START_REF] Baerlocher | An inverse kinematic architecture enforcing an arbitrary number of strict priority levels[END_REF]. The complete implementation of this approach, denoted SoT, is proposed in Mansard et al. (2007). The proposed structure enables to easily add or remove a task, or to swap the priority order between two tasks, during the control. Constraints (such as joints limit) can be taken into account. The continuity of the control law is preserved even at the instant of change.

Gain handling

The simple attractor presented in (3) produces a nice exponential decrease. However, it also produces a strong acceleration at the beginning of the task, while at the end of the task, e decreases slowly (slow convergence). A very classical 'trick' when regulating a task is to rather use an adaptive gain λ = λ(e(t)), for example:

λ(e) = (λ F -λ I ) exp -e β + λ I (5)
with λ I the gain at infinity, λ F the gain at regulation (λ I ≤ λ F ) and β the slope at regulation.

Sequence of tasks

A task sequence is a finite set of tasks sorted by order of realization, and eventually linked to each other. Any pair of tasks can be either independent (i.e. they can be achieved in parallel), or constrained (i.e. one may have to wait for another one to be achieved, in order to make sense or to be doable).

The sequence can be formulated into a classical temporal network scheduling, starting at t 0 and ending at t End . Both values are finite and the sequence does not loop. Besides, we may consider for the sake of clarity but without loss of generality that each task appears only once in the sequence.

The position of a task in the sequence is defined by the time interval during which it is maintained in the SoT. For a given task i, this interval is noted [t I i , t F i ]: the task enters in the SoT at t I i and is removed at t F i . These instants are defined with respect to the beginning of the sequence at t 0 . However, they do not indicate the achievement level of the task: t F i may apply before the task regulation. Thus we define the tolerance on the task regulation ǫ i : a task is considered as regulated when e i (t) ≤ ǫ i . The regulation time t R i is defined by e i (t R i ) = ǫ i . A task sequence is characterized by a set of timeconstraints binding the schedules of two tasks e i and e j . They can be defined as follow 1 : e i must begin or end once e j has begun, has ended or has been regulated. A graphical representation of such time-constraints is given in Fig. 3.

For example, first we have to grasp an object and maintain the force closure on it (e A ) before we can move it (e B ). The task (e B ) can only start once the task (e A ) has been realized, and must end before the task (e A ).

1 contrary to Allen Logic, that only considers the start and end points of the time interval, here is also considered the regulation time t R B begins once A has begun

t I A ≤ t I B B begins once A is realized t R A ≤ t I B B begins once A has ended t F A ≤ t I B
B ends once A has ended t F A ≤ t F B Fig. 3. Four time-dependency relations are considered.

We use the following notation to describe the set of pairs of tasks e i and e j that undergo these dependencies (e i is the direct predecessor of e j ) :

S I,I = {(e i , e j ) | t I i ≤ t I j } (6a) S R,I = {(e i , e j ) | t R i ≤ t I j } (6b) S F,I = {(e i , e j ) | t F i ≤ t I j } (6c) S F,F = {(e i , e j ) | t F i ≤ t F j } (6d)

CONTINUOUS OPTIMIZATION OF SEQUENCE OF TASKS

General problem formulation

An optimization problem is composed of a criterion to minimize, and of a set of equality and inequality constraints that must be satisfied. Our chosen criterion is to minimize the regulation duration of the mission. The variables of our problem are three for each task: (i) the time of its entry and (ii) of its removal (from the SoT), and (iii) the gains (λ I , λ F , β), which describe the task execution behavior.

The general optimization problem is written as follows: min

x t End (7a)
subject to q = SoT x (q, t) (7b) seq(q) < 0 (7c) φ(q) < 0 (7d) ∀i, t F i ≤ t End (7e) The vector x gathers the optimization variables of each task:

x = [t I 1 , t F 1 , λ I 1 , λ F 1 , β 1 , . . . , t I n , t F n , λ I n , λ F n , β n ]. t
End is the duration of the mission, seq(q) and φ(q) are respectively the sequencing and the robotic constraints.

The optimization criterion t End is computed indirectly. An equivalent explicit definition could be given by t End = max i (t F i ). However this constraint is not smooth. Giving only (7b), the problem is smooth and properly defined: at the optimal solution, t End will be equal to the maximum termination time of all tasks' t F i . Vector q is in fact a vector of functions of time, hence constraints φ(q) are semi-infinite, i.e. taking place for all the values of the continuous variable t ∈ [0, t End ].

It can be shown that (7) defines a continuous optimization problem. However, it cannot be solved directly because of the semi-infinite nature of the constraints. Therefore we expanded the semi-infinite constraint into a discreet form.

Constraints

x must satisfy both the sequencing and the robotic timeconstraints enumerated hereafter:

Tasks constraints: noted seq(q) and defined as follow:

• Time coherence for each task i, that is:

∀i, 0 ≤ t I i < t F i ≤ t End (8) • Termination condition for each task i, that is: ∀i, s * i -s i (t F i ) < ǫ i (9)
• Any task sequence condition for a given couple of tasks i and j described in (6) • The control gain, namely:

∀i,

λ I i ≤ λ F i ( 10
)
The constraints (6a), (6c), ( 6d) and ( 8) are linear. On the contrary, the constraint ( 9) is impossible to compute directly using x, and is determined from a simulation of the execution. Care has to be taken while resolving the condition described by (6b). Indeed, discretizing t R to the closest simulation step will produce discontinuities which may disturb the optimization process. A rather fastidious solution to this continuity problem would be to determine this point by interpolation. Another solution is to reformulate (6b) and evaluate the regulation of the task i when the task j begins. The constraint (6b) becomes:

∀(i, j) ∈ S R,I , s * i -s i (t I j ) ≤ ǫ i (11)
Robot constraints: φ(q)

Those constraints are mainly due to hardware intrinsic limitations of the robot:

• Joint limits, given by q min ≤ q ≤ q max (12) q min , q max are respectively the lower and upper joint limits.

• Velocity limits, given by qmin ≤ q ≤ (13) qmin , qmax are respectively the minimal and maximal velocity limits for each joint.

• Collision avoidance between a given pair of objects i and j, given by d ij ≥ 0 (14) d ij is the distance between objects i and j. Objects designate those found in the mission's environments and each link of the robotic system. Hence, both collision with the environment and self-collision of the robot have to be evaluated [START_REF] Benallegue | Fast C 1 proximity queries using support mapping of sphere-toruspatches bounding volumes[END_REF].

All of those constraints are semi-infinite: the following section presents how they have been tackled.

Technical aspects of the optimization resolution

Semi-infinite constraints: In a first approach, we tried to discretize the semi-infinite constraints on the basis of the simulation steps grid. However, since the number of the grid sample points changes in function of t End , the number of constraints is variable. Subsequently a classical optimization solver can not handle it.

Let c be the evaluation value of a given constraint: (∀t ∈ [t I , t End ], c(t) < 0). We considered associating only one value to the constraint, c V , that is computed as follows: if the constraint is always satisfied, then c V is the higher value of c(t). Otherwise, it is the sum of all the violations. The reformulation c V < 0 acts similarly to the semiinfinite constraint c(t).

Scaling:

Since the constraints are not homogeneous (times, angles, velocities, distances), they have to be normalized based on the constraint values obtained while executing the sequence corresponding to the initial set of parameters x 0 . This simple scaling improves significantly the convergence of the optimization.

Absolute versus relative timing

In this parameterization, the tasks and constraints are described with an absolute time. As it is, decreasing t I i for a task i will not have any direct effect on t F i : we have also to decrease t F i then decrease t End : it is thus necessary to propagate the reduction for all the following tasks. To avoid this, another parameterization consists in describing the SoT entry time of a given task with respect (i.e. relatively) to the previous one. We introduce a relative timing: each task is now described by two delays (instead of the absolute times t I and t F ), namely:

(1) dt I : is the delay which occurs between (i) the maximum time of entry, of end or of regulation of the preceding tasks, and (ii) the SoT entry time of the task in question.

(2) dt F : is the delay between the SoT entry and the removal times of the task in question.

These two delays fulfil the following equations:

t I j = max max (i,j)∈SI,I {t I i }, max (i,j)∈SR,I {t R i }, max (i,j)∈SF,I {t F i } +dt I (15) t F i = t I i + dt F i (16)
Subsequently, the new parameter vector is noted :

x ′ = [dt I 1 , dt F 1 , λ I 1 , λ F 1 , β 1 , . . . , dt I n , dt F n , λ I n , λ F n , β n ].
If the task sequence is only a chain of tasks realized one after the other, we directly have x ′ = f (x), with f a linear function, and

t End = i dt I i + i dt F i
Considering this new set of parameters, the formulation of the optimization problem changes. The objective remains the same: minimize t End , while the constraints become:

• Constraints on the delay:

∀i, 0 ≤ dt I i (17) ∀i, 0 < dt F i (18) • Constraint of termination ∀i, s * i -s i (t F i ) ≤ ǫ (19)
• Constraint on the task sequence

∀(i, j) ∈ S F,F , t F i ≤ t F j (20)
The constraints (6b), (6c) and ( 8) are now replaced by the constraints ( 17) and ( 18). Each task's start depends on its predecessors in the sequence. Because of the possible dependence on the regulation time, the absolute times t I i , t R i and t F i can not be pre-computed, and are computed during the simulation.

IMPLEMENTATION

Presentation

At each optimization step, the solver chooses a new set of parameters x. It then computes the constraints. Constraints ( 17) and ( 18) can be evaluated directly. As stated previously, the other constraints can not be directly computed (since they do not write in an analytical formulation). They are thus evaluated using a complete simulation of their execution. The chosen value of the current optimization variable vector x is transmitted by the optimization solver to the simulation engine. The simulation returns the evaluation of the constraints and the optimization solver computes a new step vector x, until convergence. The optimization solver is chosen to be the SQP algorithm from the MATLAB optimization toolbox; the simulation engine runs also under MATLAB.

Simulation

In Section 2, we presented the computation of the jointvelocities control law (7). The simulation is basically a numerical integration of this last equation (using an explicit Euler integration method with a fixed step ∆t = 0.005sec). The starting t I i and ending t F i times of tasks are continuous variables that are not aligned with the grid. Those instants are important since they correspond to a change in the SoT and thus a change in the control. If postponing the change of control to the next time step (like on a real system) we will not have a continuous problem (hence potentially raising the same problem described in Section 3.2). To solve this problem, the entry time t a for a given task is added as an integration point during the time step [t, t + ∆t], resulting in the separation of the time step into the two smaller time steps [t, t a ] and [t a , t + ∆t].

Initialization t End = t End max = i dt I i + i dt F i t = 0 while (t < t End ) do [t I 1 , t R 1 , t

EXPERIMENT

Temporal network

The sequence of tasks (Fig. 4) describes a robot taking out an object placed behind a door. The corresponding tasks are:

• This is a complex mission that can not be split into smaller sequences. Indeed, the sequence is centered on the door: the grasping part does not make sense if it is closed. Instead of adding an explicit timing conditions between the tasks to ensure that this will never occur, we chose to consider as constraint the collision between the left arm and the door, in order to allow task overlapping (the condition t I 3 ≤ t I 4 only postpones the entry time of the task without preventing overlapping). Similarly, we considered the collision constraint between the gripper and the handle or the object. The gripper closing tasks are used to grasp and hold objects, so these tasks have to be maintained as long as necessary.

The constraints considered for this problem are thus sequencing and robotic constraints (joint position and velocity limits).

Results of the optimization

We ran the optimization on a 3GHz desktop PC running under Windows OS. No specific effort of software optimization has been made. The sequence found is described on Fig. (5).

The overlaps between the tasks of the left and the right arm appear clearly: the left arm starts to move before the door is open. It then starts to move toward the object pose even if the door is not completely open. And finally, the right arm starts to close the door before the left arm has completely left the door area. The whole task sequence lasts 46.9sec. Without these two overlaps, the robot will (i) move to and grasp the object (e 4 ) only after the door is fully opened (e 3 ) and it will close the door (e 7 ) only after the object is completely taken out (e 6 ); then the total mission would have taken 54sec. Fig. 5. Result of tasks sequencing and behavior optimization. Each task is described by two periods: the dark one is the achievement period [t I i , t R i ], the bright one is the SoT presence period [t I i , t F i ]. The obtained execution is plotted on Fig. 6 which underlines the relation between the norm of error of the tasks and velocity of the robot joints.

CONCLUSION

We devise a method which allows to optimize both the behavior and the overlapping scheduling of a sequence of tasks composing a robotic mission. The solution derives from an optimization formulation of the tasks scheduling keeping the formalism built on the top of a task-function based control. This allows to include the robot limitations as well as collision avoidance as constraints. Our method is exemplified through a complete simulation of a complex mission, where we demonstrated an improvement in the smoothness of the generated motion. For the time being, our method still needs a predefined ordered sequence. As a future work we will increase the autonomy by determining automatically the ordered sequence and compute all the necessary subtasks from definitions of actions/objects associations. We will also focus on more complex scenario, doable by a humanoid robot, using in particular visual interaction.

Fig. 4 .

 4 Fig. 4. Sequence describing the robot taking the object in two. The function handleStackOfTasks computes the velocity of the robot induced by to the tasks execution and integrates it, altogether with any other simulated objects or processes, to obtain the new positions.

  e 1 Move the right arm to the door • e 2 Close the right gripper • e 3 Open the door • e 4 Move the left gripper to the object • e 5 Close the left gripper • e 6 Remove the object out • e 7 Close the door • e 8 Open the right gripper • e 9 Move the right arm away

Fig. 6 .

 6 Fig. 6. Simulation: decrease of the errors when the final sequence is run

⋆ This work is partially supported by grants from the ImmerSence EU CEC project, Contract No. 27141 www.immersence.info (FET-Presence) under FP6.