Guilhem Barruol
email: barruol@upf.pf

Michel Granet

Jeune E ¤quipe

A Tertiary asthenospheric £ow beneath the southern French Massif Central indicated by upper mantle seismic anisotropy and related to the west Mediterranean extension

Keywords: S-waves, wave splitting, anisotropy, Central Massif, asthenosphere, Tertiary, extension

Upper mantle flow beneath the French Massif Central is investigated using teleseismic shear wave splitting induced by seismic anisotropy. About 25 three-component stations (short period, intermediate and broadband) were installed during the period 1998^1999 in the southern Massif Central, from the Clermont Ferrand volcanic area to the Mediterranean Sea. Teleseismic shear waves (SKS, SKKS and PKS) were used to determine the splitting parameters: the fast polarization direction and the delay time. Delay times ranging between 0.7 and 1.5 s have been observed at most of the sites. The azimuths of the fast split shear waves trend homogeneously NW^SE in the southern Massif Central suggesting a homogeneous mantle flow beneath this area. The observed NW^SE direction differs from the N100 ‡E Pyrenean anisotropy further south. It does not appear to be correlated to Hercynian structures nor to the present-day motion of the plate but is well correlated to the Tertiary extension direction. We propose that the opening of the western Mediterranean induced by the rotation of the Corsica^Sardinia lithospheric block and the roll-back to the SE of the Tethys slab may have generated a large asthenospheric mantle flow beneath the southern Massif Central and a deflection of the up going plume centered beneath the northern Massif Central toward the SE.

Introduction

The Massif Central (MC) in France is a portion of the Hercynian belt that has been uplifted dur-ing Cenozoic times. It represents one of the largest volcanic provinces of western Europe whose geological evolution has been related with the Alpine orogeny [START_REF] Merle | The formation of the West European rift: a new model as exempli¢ed by the Massif Central area[END_REF] and the Cenozoic extensional regime. The MC has been a¡ected, as the whole of western Europe, by the Oligocene extension, which induced a thinning of the crust beneath the Limagne graben [START_REF] Perrier | Structure profonde du Massif Central franc °ais[END_REF]. The recent volcanism led to a shallowing of the lithosphereâ sthenosphere transition to depths of about 50 km [START_REF] Sobolev | Upper mantle temperatures and lithosphere-asthenosphere system beneath the French Massif Central constrained by seismic, gravity, petrologic and thermal observations[END_REF].

Although the volcanic activity in the MC began during Paleocene times on its periphery (Bourgogne, Causse, Bas Languedoc), the main volcanic activity occurred during Miocene and Pliocene times in the central part (Cantal, Mont Dore). The maximum geographical spreading of this volcanism occurred during the lower Pliocene (between 7 and 4 Ma) [START_REF] Maury | Le volcanisme tertiaire et quaternaire en France[END_REF]. Some of the volcanoes were still active during Quaternary times (Vivarais, Languedoc and Cha| "ne des Puys, for instance) [START_REF] Brousse | Le volcanisme en France[END_REF]. From late Miocene to present, an age progression has been described along three branches originating from the Cantal area where the oldest (between 14 and 10 Ma) volcanism is observed, whereas the most recent volcanism (younger than 1 Ma) is found in the Cha| "ne des Puys to the north, in the Languedoc area to the south, and in the Arde 'che area to the southeast.

The anomalously high topography, the negative Bouguer anomaly [START_REF] Lucazeau | Evolution ge ¤othermique et ge ¤odynamique du Massif Central franc °ais depuis l'Oligoce 'ne[END_REF], the high heat £ow [START_REF] Lucazeau | Interpretation of heat £ow data in the French Massif Central[END_REF], the surface wave attenuation [START_REF] Souriau | Attenuation of Rayleigh wave across the volcanic area of the Massif Central, France[END_REF] and the low P-wave velocity anomaly [START_REF] Granet | Massif Central (France): new constraints on the geodynamical evolution from teleseismic tomography[END_REF][START_REF] Judenherc | Two-dimensional anisotropic tomography of lithosphere beneath France using regional arrival times[END_REF] are all compatible with the presence of hot material at upper mantle depths beneath the MC. Numerical models [START_REF] Sobolev | Upper mantle temperatures and lithosphere-asthenosphere system beneath the French Massif Central constrained by seismic, gravity, petrologic and thermal observations[END_REF][START_REF] Larive | Mode ¤lisation du manteau supe ¤rieur sous le Massif Central franc °ais: inversion couple ¤e des donne ¤es de tomographie sismique et des anomalies de pesanteur[END_REF] suggest that the mantle temperature beneath the northern MC is 100^200 ‡C higher than the surrounding mantle. Although the texture analyses of mantle xenoliths [START_REF] Nicolas | Peridotites xenoliths in Massif Central basalts: textural and geophysical evi-dence for asthenospheric diapirism[END_REF] suggested the presence of several small scale asthenospheric diapirs intruding the lithosphere, the size and extension of the mantle intrusion(s) are still debated.

In 1991^1992, seismic stations were deployed in the northern MC [START_REF] Granet | Massif Central (France): new constraints on the geodynamical evolution from teleseismic tomography[END_REF]. The resulting upper mantle seismic tomography evidenced a large, low velocity body beneath the Cantal and Velay areas, visible down to about 250 km [START_REF] Granet | Massif Central (France): new constraints on the geodynamical evolution from teleseismic tomography[END_REF], interpreted as a possible signature of the hot mantle plume from which the MC volcanism originated [START_REF] Granet | Imaging a mantle plume beneath the French Massif Central[END_REF]. Interestingly, this tomographic image suggested that the low velocity anomaly could have some vertical and southward extent, but the design of this 1991^1992 experiment did not allow the characterization of the extension of this low velocity anomaly.

In order to enlarge this tomographic image at the scale of the whole MC, we conducted a temporary deployment of a seismic network in 19981 999. We installed 25 three-component seismic stations on the southern part of the MC. Our aim was to extend the 3-D velocity model not only deeper but also southward to the Mediterranean Sea and to make the connection with the Pyrenean tomographic model [START_REF] Souriau | A tomographic study of the lithosphere beneath the Pyrenees, from local and teleseismic data[END_REF]. This deployment should help to constrain the deep structures of the mantle beneath the southern MC, to map the lateral variations of the lithospheric thickness, and to better constrain the lateral and vertical extent of the mantle plume. Mapping the upper mantle £ow from seismic anisotropy and teleseismic shear wave splitting was another aim of this seismic study that we present in this paper.

For the last decade, seismic anisotropy has been used to study upper mantle tectonics [START_REF] Silver | Seismic anisotropy beneath the continents: probing the depths of geology[END_REF][START_REF] Savage | Seismic anisotropy and mantle deformation: what have we learned from shear wave splitting?[END_REF]. Shear wave splitting is directly induced by seismic anisotropy: a polarized shear wave crossing an anisotropic medium is split into two perpendicularly polarized waves that propagate at di¡erent velocities. Two anisotropy parameters may be retrieved from three-component seismic records: the di¡erence in arrival time (Nt) between the two split shear waves, which depends on the thickness and on the intrinsic anisotropy of the medium, and the orientation of the split waves' polarization planes (the azimuth P for the fast wave), which is related to the orientation of the structure.

Seismic anisotropy in the upper mantle is broadly assumed to result from intrinsic elastic anisotropy of rock-forming minerals [START_REF] Mainprice | The seismic anisotropy of the Earth's mantle: from single crystal to polycrystal[END_REF] and from their preferred orientations developed in response to tectonic £ow. Olivine, which represents the main upper mantle constituent and which is the most anisotropic upper mantle mineral phase, controls upper mantle anisotropy [START_REF] Nicolas | Formation of anisotropy in upper mantle peridotites: a review[END_REF][START_REF] Ismail | An olivine fabric database: an overview of upper mantle fabrics and seismic anisotropy[END_REF]. Since shear wave splitting is a direct result of anisotropy, and hence of rock deformation, it is possible by measuring this splitting to investigate upper mantle structure in relation to plate tectonics. We present in this paper observations of splitting of teleseismic shear waves recorded on the southern £ank of the MC. We also discuss the possible upper mantle fabric in light of the various tectonic events that could have played a role: the motion of the Eurasian plate, the Hercynian orogeny, the hotspot emplacement, and the Tertiary extension related to the opening of the western Mediterranean.

Data and results

Teleseismic events were recorded by three-component seismic stations mainly from the French LITHOSCOPE network. The main station deployment occurred between October 1998 and August 1999. Table 1 lists the station locations presented in Fig. 1. We began the experiment In order to avoid noise contamination of the splitting measurements, we kept only events characterized by signal to noise ratio of the SKS phase higher than 2. From the whole data set of teleseismic events located at a distance greater than 85 ‡ and of magnitude greater than 5.8, about 60 events were suitable for shear wave splitting measurements (Background Data Set1 , Table 2). Of these events, less than 20 were recorded during the main experiment (January^August 1999).

The event origins and locations (Background Data Set 1 , Table 2) are taken from the US Geological Survey (USGS) Preliminary Determination of Epicenters, and the phase arrivals were computed using the IASP91 Earth reference model [START_REF] Kennett | Seismic traveltime tables[END_REF].

Most individual splitting measurements (Background Data Set 1 , Table 3) were performed on earthquakes occurring at distances in the range 85^120 ‡. We measured the shear wave splitting using the Silver and Chan [START_REF] Silver | Shear wave splitting and subcontinental mantle deformation[END_REF] algorithm. This method determines the anisotropy parameters, P and Nt, that best minimizes energy on the transverse component of the seismogram for a selected time window. For each, we report the split phase on which we performed the measurement, the splitting parameters (P, Nt) with their 1c uncertainty, determined from the 95% con¢dence interval in the (P^Nt) domain, and the backazimuth of the event. We ascribe a quality factor (good, fair, or poor) to the measurements, depending on the signal to noise ratio of the initial phase, the rectilinear polarization of the particle motion in the horizontal plane after anisotropy correction, and the waveform correlation between the fast and slow split shear waves. We present in Fig. 2 three examples of splitting measurements recorded by short (CLR), intermediate (GUZ) and broadband (PDS) sensors. We generally used the SKS phase, but for some events, we used the whole SKS+SKKS-wave train. For a few events occurring at epicentral distances between 130 and 140 ‡, we successfully measured the splitting of PKS phases. We also deduced direct S-wave splitting for two events: for S-wave from event 99098 occurring at an epicentral distance of 80 ‡, we considered this event deep enough (562 km) to assume that no splitting occurred on the source side of the seismic path. The result obtained from this S phase is fully consistent with the anisotropy parameters deduced from SKS splitting measurements. We also used event 99028 occurring at a depth of 66 km and at an epicentral distance of about 83 ‡. The measurements were made on the whole S+SKS-wave train, the two phases arriving almost simultaneously at this distance. Although the S and SKS rays do not sample exactly the same mantle (they are 24 km apart at 150 km depth), the quality of the measurements at the southern stations and their similarity to the SKS splitting parameters led us to consider that this particular splitting measurement should primarily re£ect the anisotropy beneath the station C Fig. 2. Examples of splitting measurements at three stations equipped with di¡erent sensors: SKKS phase for event 98271 at CLR, with a short period (1 s) sensor, SKS phase for event 99024 at GUZ with an intermediate period sensor (5 s), SKS phase for event 99315 at PDS with a broadband sensor. For each station we show two upper traces: the initial radial (a) and transverse (b) components (energy on the transverse component). Two lower traces (c and d): the radial and transverse traces corrected for anisotropy (there is no longer energy on the transverse component). The shaded area gives the time window on which the splitting measurement is done and the dashed lines represent the predicted phase arrival times (IASP91 model). The four diagrams on the right plot the fast and slow split shear waves (continuous and dashed line, respectively) raw (e) and corrected (f) for the best-calculated delay time. Particle motions in the horizontal plane are shown below, also uncorrected (g) and corrected (h) for the anisotropy: the elliptical particle motion becomes rectilinear when the anisotropy is corrected.

and not a source side anisotropy. This S+SKS splitting measurement was the only one done at three stations from the northern group (St08, St09 and St11, reported in Table 1) that were running for a short time. We considered that the quality of these data is good enough to keep measurement and to have some con¢dence in the results, but they are only weakly representative of the region and we shall not discuss the results together with the southern stations.

The individual results of shear wave splitting are plotted in Fig. 3 in polar diagrams. The azimuths of the fast split shear wave together with the delay time is plotted on the left panel of such pairs of ¢gures, whereas the backazimuths for which no splitting were detected are plotted on the right. These are called 'null' measurements and indicate either that there is no anisotropy beneath the station or that the initial polarization direction of the SKS-wave is parallel to the fast or slow direction in the anisotropic layer. Black lines represent the best-constrained results whereas the gray lines represent measurements of lower quality. Good splitting measurements were observed at almost all sites except four. For three of them (ALB, LUC and RAB), this is likely related to a Fig. 3. For each station are presented the splitting (left) and the backazimuth for which no splitting is observed (right). For the splitting measurements, the trend of each segment represents the azimuth P of the fast split shear wave and its length is proportional to the delay time Nt (up to 3.0 s). Black lines correspond to well-constrained results, dark gray lines to fair and dashed lines to poorly constrained results (Background Data Set 1 , Table A3).

C Fig. 4. Map of the SKS splitting measurements at the surface (a), and projected along the rays at 100 km (b) and at 200 km depths (c). As in Fig. 3, the black measurements are the best constrained, and the gray measurements represent measurements of lower quality. lack of data, since very few nulls are observed during the same period. On the other hand, for the permanent station CFF at Clermont Ferrand, the absence of splitting observation could be related to the absence of anisotropy along the vertical direction beneath the station, since numerous nulls of good quality are observed. From systematic analyses of more than 10 years of data at the geoscope station SSB further east, Barruol and Ho¡mann [START_REF] Barruol | Seismic anisotropy beneath the geoscope stations from SKS splitting[END_REF] observed only two positive splitting results (arriving at the same backazimuth) and a large number of nulls, suggesting a mostly isotropic structure beneath the station. This apparent isotropy along the vertical direction may result either from vertically oriented olivine aaxis, or from large scale heterogeneities, or from the presence of two anisotropic layers with perpendicularly oriented fast axes or simply from absence of deformation of the medium [START_REF] Barruol | Seismic anisotropy beneath the geoscope stations from SKS splitting[END_REF]. As shown below, measurements of olivine and orthopyroxene crystal lattice preferred orientation of peridotite nodules sampled in the northern MC [START_REF] Ho¡mann | Upper mantle seismic anisotropy in the vicinity of mantle plumes[END_REF] do not favor this last hypothesis.

SKS splitting has a poor vertical resolution because the splitting of a teleseismic shear wave observed at the Earth's surface may be acquired anywhere along the ray between the station and the core^mantle boundary. It is, however, broadly accepted from seismological and petrophysical arguments that most of the anisotropy lies within the uppermost 400 km of the Earth. There are some qualitative arguments that can help to vertically locate the anisotropy within the uppermost mantle. At a given site, for instance, homogeneous anisotropy parameters deduced from events with di¡erent backazimuths strongly favor super-¢cial and homogeneous mantle deformation. We observe such homogeneity at the central stations and this could suggest a shallow (i.e., above 200 km) anisotropy. The lithosphere is relatively thin beneath the northern MC ( 6 60 km), and probably less than 80 km thick beneath the southern MC, so part of the recorded anisotropy may therefore lie within the asthenospheric upper mantle. Alternatively, the scattering of the anisotropy parameters observed at few southwestern stations (particularly at RAB and SAI) suggests either some super¢cial (i.e., 6 100 km deep) heteroge-neous structure beneath the station (lateral variations, for instance) or the presence of several layers of anisotropy that may induce some backazimuthal variations in the anisotropy parameters [START_REF] Silver | The interpretation of shearwave splitting parameters in the presence of two anisotropic layers[END_REF].

In order to test the extent of the upper mantle deformation, i.e., the large scale stability of the anisotropy parameters, and whether scattering at a given station re£ects lateral structure variation or more complex structures such as the presence of several anisotropic layers, we present in Fig. 4 anisotropy maps on which each individual splitting measurement is projected at a given depth along the incoming ray. This representation is both helpful, since it plots the splitting parameters at the place where it could have been partially acquired, but also imperfect because it plots a ¢nal, integrated splitting at a depth where splitting was perhaps not fully completed. The individual splitting measurements (Background Data Set 1 , Table 3) are projected at the surface (Fig. 4a), at 100 (Fig. 4b) and 200 km depth (Fig. 4c). These maps ¢rst show that the lateral mantle sampling is homogeneous despite the relatively sparse network density (about 50 km between stations) and the sparse backazimuthal coverage. Events mainly arrived from the NE and the SW. Second, except for some lower quality events, the homogeneous NW^SE trend in the central part of the network is particularly striking in the 100 and 200 km depth maps and suggests a uniform and pervasive structure beneath this area. This 'channel' structure appears to be at least 400 km long (and in fact much longer when incorporating the individual measurements made at the northern stations St08, St09 and St11 and the reprocessed data at the northern MC stations [START_REF] Babuska | Seismic anisotropy of the French Massif Central and predisposition of Cenozoic rifting and volcanism by Variscan suture hidden in the mantle lithosphere[END_REF]) and about 200 km wide. Third, the scattering observed at SAI and RAB appears from these projections to be compatible with variations in the lateral structure at mantle depth. At SAI, for instance, the three measurements with roughly N^S trending P (events 99095, 99130 and 99323) sample the same upper mantle, whereas the two NE^SW trending P (events 99024 and 99345) sample a mantle area slightly o¡set to the NE. At RAB, the absence of good measurements weakens the conclusions, but the two N^S P measurements (events 99093 and 98281) sample mantle anisotropy SW from the station, whereas the N110 ‡E trending P (event 99098) integrates anisotropy NE from the station. Although the small number of splitting measurements at these sites does not allow to test the presence of several anisotropic layers, the results from these southwestern stations are compatible with lateral heterogeneities in the upper mantle.

Despite the small data set related to the short interval of the network deployment, the overall coherency of the results at most stations allows us to calculate the mean splitting parameters at each site. We used the averaging method presented by Silver and Chan [START_REF] Silver | Shear wave splitting and subcontinental mantle deformation[END_REF] that weights each individual measurement by its 95% con¢dence interval (reported in Background Data Set 1 , Table 3), so that poorly constrained results have a smaller contribution to the mean result than the best-constrained measurements. To check the consistency of the results, we determined the mean splitting parameters (P, Nt) using two sets of data: (1) the whole set of non null measurements, and (2) the good non null measurements only (Table 1). A ¢rst group of four stations (CFF, ALB, LUC, RAB), devoid of good splitting measurements, does not allow further discussion of the di¡erence between the mean anisotropy determined from the two sets of data. For a second group of stations (i.e., PDS, CLR, CNG, DIO, ENS, GUZ, BOPU, BLRC, MIR, ARBF), the mean azimuth P (ranging between N094 ‡ at ARBF and N138 ‡ at CLR) obtained from the two di¡erent data sets are consistent within þ 10 ‡. Such a coherency is primarily observed in a NW^SE trending corridor. However, mean delay times at some sites may display variations larger than 0.2 s, depending on the quality of the data taken into account. This occurs particularly at stations where very few good measurements are observed (e.g., CNG, MIR, ENS). The third group of stations shows larger di¡erences ( s 10 ‡ in P and s 0.1 s in Nt) in both anisotropy parameters deduced from the two data sets. This re£ects a stronger scattering in the splitting observations. These stations are located on the southwestern (e.g., SAI, COU, RAB) and northeastern part of the network (e.g., RMZ, LUC).

In the central area of our network, the mean azimuth P of the fast split shear wave is homogeneously aligned along a NW^SE trending corridor (Fig. 5). More precisely, P trends NW^SE in the northwestern part of the network and rotates slightly to almost E^W in the southeastern part of the network. Along this corridor, the mean Nt deduced from the good splitting observations (Table 1) are not characterized by such smooth variations as the azimuths' P: the northern stations (PDS, CNG, CLR and BLRC) display small Nt (in the range 0.61^0.95 s), whereas southern stations (GUZ, DIO and ARBF) display much higher Nt (in the range 1.11^1.36 s). MMJ, in the Rho " ne valley, has much lower Nt (0.64 s) than surrounding stations. The southernmost station BOPU (P = N116 ‡E, Nt = 1.43) is characterized by similar parameters to those observed in the Pyrenees [START_REF] Barruol | Lithospheric anisotropy beneath the Pyrenees from shear wave splitting[END_REF]: a high Nt and P trending roughly parallel to the Pyrenean belt. Fig. 5. Map of the mean splitting parameters calculated from the good measurements at each station from the southern MC and the Pyrenees [START_REF] Barruol | Lithospheric anisotropy beneath the Pyrenees from shear wave splitting[END_REF]. There is no result at some sites (RAB, ALB, LUC), since they are devoid of good SKS splitting measurements. Also plotted are the good quality splitting measurements in the northern MC [START_REF] Babuska | Seismic anisotropy of the French Massif Central and predisposition of Cenozoic rifting and volcanism by Variscan suture hidden in the mantle lithosphere[END_REF]. The stars indicate the absence of anisotropy observed at SSB [START_REF] Barruol | Seismic anisotropy beneath the geoscope stations from SKS splitting[END_REF] and CFF.

Discussion: the origin of the anisotropy

In order to discuss the nature of the mantle deformation from which the observed splitting originated, we present in Fig. 5 our measurements together with anisotropy measurements from two other studies performed in the northern MC [START_REF] Babuska | Seismic anisotropy of the French Massif Central and predisposition of Cenozoic rifting and volcanism by Variscan suture hidden in the mantle lithosphere[END_REF] and in the Pyrenees [START_REF] Barruol | Lithospheric anisotropy beneath the Pyrenees from shear wave splitting[END_REF].

From data recorded by a few three-component stations that were running during the 1991^1992 seismic experiment, some anisotropy measurements were performed by Granet et al. [START_REF] Granet | Anisotropic measurements in the Rhine graben area and the French Massif Central[END_REF] and more recently reprocessed by Babuska et al. [START_REF] Babuska | Seismic anisotropy of the French Massif Central and predisposition of Cenozoic rifting and volcanism by Variscan suture hidden in the mantle lithosphere[END_REF]. Although the results are deduced from a single event, one of the major ¢ndings of this study is a strong variation in the anisotropy parameters obtained at stations located on both sides of the N^S trending 'Sillon Houiller', which was an active fault during the Variscan orogeny. Stations on the western side of the fault are characterized by P trending N100 ‡E, whereas stations lying on its eastern side are characterized by P trending roughly N^S. These observations led the authors to conclude that the anisotropy beneath the northern MC is most likely located in the uppermost mantle and that the 'Sillon Houiller' could represent a major lithospheric boundary.

Further south in the Pyrenees, seismic anisotropy measurements [START_REF] Barruol | Lithospheric anisotropy beneath the Pyrenees from shear wave splitting[END_REF] displayed a homogeneous pattern characterized by relatively strong delay times (generally above 1.0 and often above 1.3 s) and P trending parallel to the belt (i.e., N100 ‡E). The small scale delay time variations observed in the central part of the belt [START_REF] Barruol | Anisotropy beneath the Pyrenees range from teleseismic shear wave splitting[END_REF], but also the correlation of P with the Hercynian and Pyrenean tectonic structures [START_REF] Barruol | Lithospheric anisotropy beneath the Pyrenees from shear wave splitting[END_REF], favor an anisotropy located primarily within the lithosphere. Geological and geophysical arguments also suggest that the resulting mantle deformation could re£ect Pyrenean structures in the internal part of the belt but Hercynian structures beneath its external parts.

Our new measurements from the southern MC appear to be statistically di¡erent from those obtained in the Pyrenees: the azimuths' P are oriented more NW^SE rather than N100 ‡E in the Pyrenees and the Nts are statistically lower than 1.0 s, particularly for the northern stations. We now draw inferences about the source regions of anisotropy and the processes that generate mantle deformation fabrics through a comparison of anisotropy signatures expected from these processes and observed splitting parameters.

Plate motion-related asthenospheric £ow

At the scale of southern France, seismic anisotropy is characterized by distinct patterns (Pyrenees, northern MC, southern MC), suggesting either di¡erent sources of anisotropy and/or different tectonic processes or histories. Small scale variations in the anisotropy parameters are found in each pattern, suggesting a near surface origin of anisotropy. Although it is not possible to totally reject such a hypothesis, our observations are hardly compatible with an anisotropy primarily induced by an asthenospheric £ow related to the present-day plate motion. The hypothesis explains neither the general trend of P, which is expected to trend NE^SW for Eurasia [START_REF] Gripp | Current plate velocities relative to the hotspots incorporating the Nuvel-1 global plate motion model[END_REF], nor the short scale variation in anisotropy parameters across lithospheric boundaries (across the North Pyrenean Fault or the 'Sillon Houiller'), nor the strong variations between the large Nt observed in the Pyrenees and the smaller Nt in the northern MC.

Hercynian, lithospheric fabric

The Hercynian belt outcropping in the southern MC is made of metamorphic terranes intruded by granitic bodies. This area is composed of crustal scale, south verging nappes, described by Matte [START_REF] Matte | Tectonics and plate tectonics model for the Variscan belt of Europe[END_REF] as the southern external domains of the belt. Although the surface geology at the scale of the southern MC displays strong heterogeneities both from composition (gneiss, schists, granites, T) and from the trend of the various units, the azimuths of the nappes' displacement trend roughly from N180 ‡ to N200 ‡E [START_REF] Matte | Tectonics and plate tectonics model for the Variscan belt of Europe[END_REF]. One could argue that this structural pattern outcropping in the Hercynian crustal basement could be related to a uniform and pervasive upper mantle fabric. We wish to show, however, that there are strong arguments against this assumption.

First, the Hercynian basement is hidden in many places by Mesozoic or Cenozoic sedimenta-ry covers. More than half of our stations were not directly installed on Hercynian basement but, instead, on more recent sedimentary covers (Fig. 1), particularly beneath the southeastern part of our network. Analyzing geometrical relationships between Hercynian structures and the observed anisotropy in this part of our network is therefore impossible.

Second, there is no consensus on the expected upper mantle £ow in such nappe tectonic context. Silver [START_REF] Silver | Seismic anisotropy beneath the continents: probing the depths of geology[END_REF][START_REF] Silver | Mantle structural geology from seismic anisotropy[END_REF] proposed that during a lithospheric collision, olivine a-axes could orient normal to the maximum compressive direction. The resulting fast anisotropy axis should be oriented normal to the convergence direction, i.e., roughly NWŜ E to E^W [START_REF] Matte | Tectonics and plate tectonics model for the Variscan belt of Europe[END_REF], which could be consistent with our observed P pattern. On the other hand, one could consider that olivine a-axes should orient parallel to the shear direction and hence parallel to the convergence direction. In this context, olivine a-axes should be therefore oriented roughly N^S to NE, i.e., without any relationship to our observations. Comparing our observations to tectonic predictions is hence highly speculative.

Third, in nappe tectonic contexts, one can expect large o¡sets between an outcropping crustal deformation and, if existing, the related mantle deformation at depth. This o¡set depends on the crustal thickness and on the dip of the crustal deformation zone. Although the exact geometry of the MC sutures in the lower crust and upper mantle remains speculative, a 30 ‡ angle as presented by Matte [START_REF] Matte | Tectonics and plate tectonics model for the Variscan belt of Europe[END_REF] on MC crustal scale crosssection would mean that the outcropping Hercynian deformation at a given site could be related to a mantle deformation located 50^100 km further north (assuming crustal and lithospheric thicknesses of 30 and 60 km). This prevents any relationship between outcropping geological structures and seismic anisotropy.

Fourth, the Tertiary thermal events (the plume emplacement and the rifting of the Gulf of Lion) may have mechanically and thermally thinned the preexisting lithosphere and, therefore, part of the Hercynian fabric may have been erased.

In summary, the discussion of the upper mantle Hercynian fabric beneath the southern MC is highly speculative. The lack of evidence favoring a clear Hercynian signature in the recorded anisotropy does not mean that there is no Hercynian anisotropy, but that predicting an anisotropy pattern at upper mantle depth from the outcropping Hercynian geology is hazardous.

Plume passive spreading asthenospheric £ow

Upper mantle seismic velocity anomalies associated with a hot spot are rather well imaged [START_REF] Nataf | Seismic imaging of mantle plumes[END_REF]. Nevertheless, the associated upper mantle £ow patterns and related anisotropy patterns are poorly known. Splitting measurements performed at oceanic stations located on hotspot-related volcanoes yield at most sites an absence of detectable anisotropy [START_REF] Barruol | Seismic anisotropy beneath the geoscope stations from SKS splitting[END_REF][START_REF] Ho¡mann | Upper mantle seismic anisotropy in the vicinity of mantle plumes[END_REF] or, at other sites, fast split shear waves polarized parallel to the absolute plate motion [START_REF] Barruol | Seismic anisotropy beneath the geoscope stations from SKS splitting[END_REF][START_REF] Russo | Shear wave splitting and upper mantle deformation in French Polynesia: evidence for small-scale heterogeneity related to the society hotspot[END_REF]. In all these studies, the restriction to a single station did not allow the geometry of the mantle £ow pattern, relative to the plume interaction with the asthenosphere^lithosphere system, to be tested. The seismic network deployed in Iceland could have allowed this, but SKS splitting measurements [START_REF] Bjarnason | Initial results from the ICEMELT experiment: body wave delay times and shear-wave splitting across Iceland[END_REF] display a complex pattern that could result from the interaction between the Iceland plume and the mid-Atlantic ridge. More recently, ocean bottom seismometers were deployed around Hawaii to test such plume e¡ect on the upper mantle £ow [START_REF] Walker | Shearwave splitting to test mantle deformation models around Hawaii[END_REF]. Despite a small number of splitting measurements, the authors proposed that the anisotropy pattern is consistent with a superposition of a radial pattern, the expected plume e¡ect [START_REF] Ru « Mpker | Calculating splitting parameters for plume-type anisotropic structures of the upper mantle[END_REF], and a linear pattern induced by the plate motion, generating a 'parabolic asthenospheric £ow'.

In the case of the MC, and assuming a £at lithosphere^asthenosphere boundary, the expected £ow pattern related to the only plume should be primarily radial since the plate velocity is extremely low, of the order of few mm/yr [START_REF] Gripp | Current plate velocities relative to the hotspots incorporating the Nuvel-1 global plate motion model[END_REF]. However, this assumption is likely incorrect : the presence of various lithospheric blocks in the Hercynian belt [START_REF] Lenoir | Contrasting lithospheric mantle domains beneath the Massif Central (France) revealed by geochemistry of peridotite xenoliths[END_REF] but also the presence of the neighboring Alpine and Pyrenean roots show that the bottom of the lithosphere is not £at at the regional scale and therefore that the plume-related asthenospheric £ow could be strongly controlled by such topography. The absence of seismic ani-sotropy at some northern MC stations (CFF, for instance), together with the small delay times observed in the northern MC, could indicate either steeply dipping mantle lineations [START_REF] Ru « Mpker | Calculating splitting parameters for plume-type anisotropic structures of the upper mantle[END_REF] or upper mantle structures perturbated at the vertical of the plume upwelling. Our NW^SE trending anisotropy in the southern MC could be also consistent with a radial £ow pattern around the plume center, but we explain below that the plume upwelling has likely been de£ected toward the SE by the Tertiary extensive episode that a¡ected the Gulf of Lion.

Tertiary, extension-related, asthenospheric £ow

At the end of Eocene times (Priabonian) and in the early Oligocene (Rupelian), western Europe su¡ered an extensive episode that generated the so-called European Cenozoic rift system (ECRIS), trending N^S and marked by the Rhine graben, the Bresse and Limagne graben, the Ales basin in Languedoc and the Valencia trough in NE Spain. The extension direction of this Oligocene episode is about E^W in France. Its origin remains debated but could be in£uenced by a lithosphere buckling generated by the Alpine tectonic [START_REF] Michon | The evolution of the Massif Central rift: spatio-temporal distribution of the volcanism[END_REF].

In Languedoc and SE France, the extension continued during the late Oligocene and the whole Miocene. From sedimentological analyses, Se ¤ranne et al. [START_REF] Se ¤ranne | The Gulf of Lion continental margin (NW Mediterranean) revisited by IBS: an overview[END_REF] showed that this extensional event around the Gulf of Lion is tectonically distinct from the ECRIS event. The corresponding extension directions are orientated NW^SE in the southern MC [START_REF] Se ¤ranne | The Gulf of Lion continental margin (NW Mediterranean) revisited by IBS: an overview[END_REF]. This regional extension is interpreted as back-arc spreading episodes, induced by the roll-back of the Tethys slab subducting at that time toward the northwest beneath the Cor-sica^Sardinia lithospheric block. We wish to show that this slab retreat may have induced a NW^SE regional asthenospheric £ow that may be preserved since that time and that could dominate the present-day observed anisotropy. We present a summarized chronology of the events occurring in this area taken from Se ¤ranne et al. [START_REF] Se ¤ranne | The Gulf of Lion continental margin (NW Mediterranean) revisited by IBS: an overview[END_REF] and Faccenna et al. [START_REF] Faccenna | Episodic back-arc extension during restricted mantle convection in the Central Mediterranean[END_REF].

From 80 to 30 Ma, the general tectonic regime of the western Mediterranean is dominated by the northward motion of Africa inducing the closure of the Tethys ocean through the subduction of its oceanic part beneath the Eurasian active margin. At Oligocene times, the collision of the Adriatic promontory with Eurasia induced important changes in the geodynamic evolution of the western Mediterranean that can be described as a succession of two rifting events followed by the creation of two oceanic basins. (i) From 30 to 22 Ma occurred the rifting of the Gulf of Lion, i.e., the extension of the lithosphere located between the Languedoc and the Corsica^Sardinia block in its pre-rotation situation. This NW^SE extension is attributed to the early stage of the slab roll-back. (ii) It was followed between 22 and 17 Ma by the drifting of the Corsica^Sardinia lithospheric block. This rotation induced by the sinking of the slab in the upper mantle is accompanied by the creation of new oceanic lithosphere in the Liguro^Provenc °al basin. (iii) From 17 to 10 Ma occurred the rifting of the Corsica^Sardinia lithospheric block. (iv) This rifting is followed from 10 Ma to present by the escape of the Calabrian block, and the creation of the oceanic Tyrrhenian basin induced by the second episode of slab sinking. (v) In the present-day situation, the subduction of the African oceanic lithosphere is located beneath Calabria. The slab itself is well imaged in the tomographic model, down to 600 km depth [START_REF] Cimini | P-wave deep velocity structure of the southern Tyrrhenian subduction zone from nonlinear teleseismic traveltime tomography[END_REF]. It appears to lie horizontally in the transition zone beneath Sardinia [START_REF] Lucente | Tomographic constraints of the geodynamics evolution of the Italian region[END_REF], suggesting its inability to fall into the lower mantle.

In this general context, we propose that the episode of broad lithospheric extension that created the Gulf of Lion, and particularly, the following stage of slab sinking southeast from the Corsica^Sardinia block and the related opening of the Liguro^Provenc °al oceanic basin may have induced a pressure gradient within the asthenosphere. This may have moved to the SE the low viscosity plume that was upwelling beneath the northern MC. The plume could have been de-£ected toward the SE and therefore, the upper mantle beneath the southern MC may have been a¡ected by a regional scale asthenospheric £ow. Furthermore, the lithosphere beneath the southern MC may have su¡ered a mechanical and thermal thinning. To illustrate this, we present in Fig. 6 some schematic NW^SE cross-section between the MC and the present-day Calabria, at the upper mantle scale.

The West European rifting ends 30 Myr ago (Fig. 6a). At that time the rifting of the Gulf of Lion begins, induced by the ¢rst sinking of the Tethys slab beneath the Corsica^Sardinia block. This rifting is followed between 22 and 17 Ma (Fig. 6b) by the rotation of this lithospheric block and the creation of the Liguro^Provenc °al oceanic lithosphere. This rotation ends at 17 Ma and is followed by the break-up of the Sardinia block (Fig. 6c). A regional mantle £ow is needed during this period in order to ¢ll the space left behind the sinking slab. We propose that the slab retreat to the SE may have pulled the mantle from the NW. The presence beneath the MC of a hot, low viscous mantle material could explain why such a £ow arriving from the northwest would have been mechanically easier than pulling mantle laterally (the slab was likely at that time a long and continuous structure) or from below (implying likely a crossing of the transition zone). Furthermore, this NW^SE asthenospheric channeled £ow could also have been partially controlled by the neighboring Alpine and Pyrenean lithospheric roots. The absence of more recent tectonic or thermal events in the southern MC may explain why this asthenospheric £ow could be preserved in the upper mantle. The beginning of the maximum intensity in the MC volcanism is roughly synchronous with the end of the Corsica rotation. This is perhaps not fortuitous ; the southeastward asthenospheric £ow beneath the southern MC may have de£ected the ascending plume and, therefore, as the suction ended, the mantle plume may have freely continued its upward motion and generated the peak in the volcanic activity of the northern MC.

The present-day situation (Fig. 6d) shows the Tethys slab lying in the transition zone after more than 900 km of retreat, as imaged by seismic tomography [START_REF] Cimini | P-wave deep velocity structure of the southern Tyrrhenian subduction zone from nonlinear teleseismic traveltime tomography[END_REF][START_REF] Lucente | Tomographic constraints of the geodynamics evolution of the Italian region[END_REF]. The Calabria is separated from C Fig. 6. Schematic upper mantle cross-section adapted from Se ¤ranne [START_REF] Se ¤ranne | The Gulf of Lion continental margin (NW Mediterranean) revisited by IBS: an overview[END_REF], trending NW^SE from the northern MC to the present-day Calabria location. See text for interpretation. The question mark at the base of the plume indicates the absence of argument for making this plume originate in the transition zone or arriving from the lower mantle. We propose that the upper mantle seismic anisotropy observed in the southern MC is located in the sublithospheric mantle and was generated by an asthenospheric mantle £ow induced by the retreat and sinking of the Tethys slab beneath the Corsica^Sardinia block. its original Corsica^Sardinia lithospheric block by the Tyrrhenian oceanic domain. The opening of this oceanic basin occurred after 10 Ma and is related to a second episode of slab sinking. An asthenospheric mantle £ow may have been induced by this event but could likely be of local extent since this sinking slab was narrow [START_REF] Carminati | The role of slab detachment processes in the opening of the western-central Mediterranean basins: some geological and geophysical evidence[END_REF][START_REF] Gvirtzman | Residual topography, lithospheric structure and sunken slabs in the central Mediterranean[END_REF] and no longer part of a laterally continuous structure. The mantle £ow induced by the present-day slab roll-back likely has no direct in£uence on the asthenosphere beneath the southern MC. The southward migration of the MC volcanism occurring these last 2 Myr could therefore re£ect a passive plume spreading. The mantle deformation related to this ¢nal £ow could superimpose on the previous one and could not be distinguished.

The lithospheric thickness is a key point in the discussion of the origin of the anisotropy beneath the southern MC. In the northern MC, the large negative Bouguer anomaly, the low seismic velocities observed by seismic tomography [START_REF] Granet | Massif Central (France): new constraints on the geodynamical evolution from teleseismic tomography[END_REF], and the high heat £ows [START_REF] Lucazeau | Interpretation of heat £ow data in the French Massif Central[END_REF] all support a thin ( 6 60 km) lithosphere. In the south, despite the lack of quantitative data on the lithospheric thickness beneath our network, the Tertiary extension may have affected the lithosphere around the Gulf of Lion, ¢rst mechanically during the rifting episode that began 30 Myr ago and second, thermally during the slab roll-back episode and the asthenospheric £ow related to the plume de£ection. The base of the Hercynian lithosphere could have been thermally eroded during this Cenozoic episode. The N^S elongated shape of the negative Bouguer anomaly is compatible with this hypothesis. The presence of two anisotropic layers, characterized by a lower asthenospheric layer and an upper Hercynian lithospheric layer cannot be ruled out but is not clearly demanded by our observations.

The thickness of the anisotropic layer that can explain our observed delay times depends on the intrinsic amplitude of anisotropy and therefore on the £ow strength. Upper mantle nodules brought up at the surface by the volcanic activity provide a way to investigate their deformation fabrics, and therefore their seismic properties. The average S-wave anisotropy deduced from more than a hundred upper mantle nodules from various geo-dynamic environments [START_REF] Ismail | An olivine fabric database: an overview of upper mantle fabrics and seismic anisotropy[END_REF] is about 5.0% around the Z structural direction (i.e., normal to the foliation plane, corresponding to the vertical direction in our case if we assume our hypothetic asthenospheric £ow is dominated by horizontal foliations), requiring an anisotropic layer about 100 km thick to explain our 1.0 s delay times. This thickness can be considered as a lower bound, since these petrophysical measurements do not take into account variations in the structural orientation nor variations in the rock composition that can lower the overall seismic anisotropy. A regional analysis of upper mantle nodules [START_REF] Ho¡mann | Upper mantle seismic anisotropy in the vicinity of mantle plumes[END_REF] gives interesting insights. Seismic anisotropy calculated from olivine and orthopyroxene crystallographic preferred orientations on 10 peridotite nodules sampled in the MC volcanic rocks shows S-wave anisotropies ranging between 1.4 and 5.3% [START_REF] Ho¡mann | Upper mantle seismic anisotropy in the vicinity of mantle plumes[END_REF]. Interestingly, low anisotropies ( 6 3%) are found in samples from the northern MC, extracted above the highest temperature anomalies [START_REF] Sobolev | Upper mantle temperatures and lithosphere-asthenosphere system beneath the French Massif Central constrained by seismic, gravity, petrologic and thermal observations[END_REF] deduced from the seismic tomography [START_REF] Granet | Massif Central (France): new constraints on the geodynamical evolution from teleseismic tomography[END_REF]. This could suggest that the lithospheric fabric could have been partially erased by the temperature increase and perhaps by the £uid circulation related to the plume emplacement [START_REF] Lenoir | Contrasting lithospheric mantle domains beneath the Massif Central (France) revealed by geochemistry of peridotite xenoliths[END_REF]. Such low anisotropy could explain the absence of splitting at some northern stations like CFF. On the periphery of this anomaly, and particularly in the southern MC, upper mantle nodules S-wave anisotropy is in the range 3^5%, suggesting that the anisotropic layer could be 1001 50 km thick. The larger delay times observed at our southernmost stations are also compatible with a stronger £ow related to the close vicinity of the sinking slab, and the smaller Nt at the northern stations could represent more di¡use mantle £ow.

Conclusions

We have mapped upper mantle seismic anisotropy beneath the southern MC. Despite a rather small data set related to a short array deployment, we observe a consistent anisotropy pattern throughout this region, with contrasting characteristics compared to the neighboring Pyrenean pattern. Our observations suggest the presence of a NW^SE trending upper mantle £ow in the central part of the network, i.e., between the north MC main volcanic area and the Mediterranean shore. The homogeneity of the results in this area does not favor a complex structure, such as, for instance, several anisotropic layers, but more a simple 'channel' structure. We propose that this channel could be the trace of a recent asthenospheric £ow arriving from the NW and that could have ¢lled the space left behind the sinking Tethys slab that induced the Corsica^Sardinia rotation between 22 and 17 Ma. We also propose that this asthenospheric £ow could have de£ected the plume ascending beneath the northern MC and could have thermally eroded the bottom of the lithosphere, erasing part of the Hercynian lithospheric preexisting fabric. The absence of subsequent tectonic and/or thermal event (except a possible recent passive plume spreading) may explain why such an asthenospheric structure may still be visible in the upper mantle. In such an interpretation, the slight P rotation from NW^SE in the northern MC to almost E^W in SE France is not explained but could reveal an e¡ect of the Alpine collision. The complex anisotropy signature in the southwestern part of our network may be explained by the presence of lateral anisotropy variations, since this area could represent the southwestern boundary of this channel.

Fig. 1 .

 1 Fig. 1. Map of the station locations. The continuous line represents the contour of the Hercynian basement and the dashed line the 'Sillon Houiller'. The Cenozoic volcanic ¢elds are in black. The black squares (station MC) indicate stations with splitting measurements [25,27]. White squares represent short period (1 s) three-component stations from the GFZ Potsdam. Due to the short time of installation, anisotropy measurements were only possible at stations St08, St09 and St11. CFF and ARBF are Re-NaSS broadband stations. SSB is a geoscope station. The black circles represent short, then intermediate period (1, then 5 s, respectively) three-component stations, whereas diamonds represent broadband stations.

Table 1

 1 Station location and mean anisotropy parameters

	Station	Latitude	Longitude	P	cP	Nt	cNt	Number of	Quality
								measurements
		( ‡N)	( ‡E)	( ‡)	( ‡)	(s)	(s)
	ALB	43.895	2.456	88	16	1.45	0.02		g+f+p
				^^^^0	g
	CLR	43.882	3.323	138	15	0.89	0.10		g+f+p
				130	14	0.96	0.09		g
	CNG	44.400	3.226	137	5	1.03	0.19		g+f+p
				126	5	0.61	0.20		g
	COU	43.479	2.723	141	13	0.95	0.25		g+f+p
				96	8	0.70	0.20		g
	DIO	43.931	4.299	116	5	1.34	0.07		g+f+p
				116	6	1.36	0.04		g
	ENS	43.310	3.110	111	5	0.96	0.18		g+f+p
				121	12	0.70	0.21		g
	GUZ	43.741	3.935	106	5	1.16	0.16		g+f+p
				99	6	1.11	0.01		g
	LUC	44.656	3.890	74	24	1.54	0.16		g+f+p
				^^^^0	g
	MIR	44.144	2.172	129	9	1.22	0.09		g+f+p
				124	5	1.03	0.26		g
	MMJ	43.706	4.664	110	4	0.70	0.04		g+f+p
				95	12	0.64	0.22		g
	PDS	44.322	2.724	131	3	1.02	0.08		g+f+p
				134	4	0.95	0.03		g
	RAB	43.830	1.716	82	40	0.78	0.21		g+f+p
				^^^^0	g
	RMZ	44.400	4.496	121	9	0.72	0.16		g+f+p
				101	7	0.60	0.14		g
	SAI	43.362	2.172	16	14	1.26	0.23		g+f+p
				52	15	0.79	0.27		g
	BOPU	42.870	2.876	112	5	1.58	0.11		g+f+p
				116	6	1.43	0.10		g
	BLRC	44.286	3.767	133	6	0.78	0.14		g+f+p
				125	7	0.61	0.11		g
	CFF	45.762	3.102	142	6	1.79	0.70		g+f+p
				^^^^0	g
	St08	45.359	2.801	122	10	0.90	0.20		g+f+p
				122	10	0.90	0.20		g
	St09	44.860	3.833	145	4	1.15	0.12		g+f+p
				145	4	1.15	0.12		g
	St11	44.781	2.126	116	16	1.10	0.40		g+f+p
				116	16	1.10	0.40		g
	ARBF	43.492	5.333	94	7	1.29	0.15		g
				95	7	1.31	0.12		g+f+p
	Station location and mean splitting parameters calculated with all the results (g+f+p) or with the only good measurements (g).
	The number and quality of individual splitting measurements from which the average values are calculated is tabulated.

for three months with short period (1 s) sensors recording at 31.25 samples per second. We replaced these sensors in January 1999 by intermediate period sensors (5 s Lennartz) that re-

corded at the same rate up to August 1999. Two broadband stations from the French 'Re ¤seau Large Bande Mobile' equipped with Streickeisen STS-2 seismometers were installed in May 1999 and remained in the ¢eld at PDS and SAI up to

	January 2000. Two other broadband stations
	from Leeds University (Mars88 equipped with
	Guralp CMG3) ran from October 1998 to August
	2000 at BLRC and BOPU. Continuous recording
	allowed us to select and extract teleseismic events
	at convenient distances and magnitudes suitable
	for SKS, SKKS and PKS splitting measurements.

http://www.elsevier.com/locate.epsl

Acknowledgements

This project is a contribution to the program 'Inte ¤rieur de la Terre' supported by the CNRS/ INSU. The 'Lithoscope' and 'Re ¤seau Large Bande Mobile' programs provided seismic stations. Thanks to M. Kendall and D. Francis for the Leeds University broadband stations. We thank also the GFZ Potsdam for providing 10 mobile stations. Many thanks to all the people who helped in the ¢eld, particularly J. Dorel, U. Achauer, A. Vauchez, A. Tommasi, B. Ildefonse, P. Nicol, J. Bascou, and M. Duclos. We are grateful to P. Silver for his splitting analysis program. Constructive reviews provided by V. Babuska and S. Sobolev but also discussions with J.M. Dautria and M. Se ¤ranne were extremely helpful. Thanks to the Raynal, the Rimbaud and the Favier families for their welcome. We also thank the mayors of several villages for their help in providing site installations, in particular, at Mirandol, Alban, Courniou, Saissac, le Caylar, and Luc. We thank the 'Parc National des Ce ¤vennes' for their help at La Roche, and the French 'Ministe 're de la Culture' for allowing us to install stations on historical sites at Ense ¤rune and Montmajour. [AC]


Nicolas Perrin

Darwin Lau

Vincent Padois

Effective Generation of Dynamically Balanced Locomotion with Multiple Non-coplanar Contacts

Studies of computationally and analytically convenient approximations of rigid body dynamics have brought valuable insight into the field of humanoid robotics. Additionally, they facilitate the design of effective walking pattern generators. Going further than the classical Zero Moment Point-based methods, this paper presents two simple and novel approaches to solve for 3D locomotion with multiple non-coplanar contacts. Both formulations use model predictive control to generate dynamically balanced trajectories with no restrictions on the center of mass height trajectory. The first formulation treats the balance criterion as an objective function, and solves the control problem using a sequence of alternating convex quadratic programs. The second formulation considers the criterion as constraints, and solves a succession of convex quadratically constrained quadratic programs.

Introduction

Locomotion is a challenging task for humanoid robots, even within completely known environments. Although the dynamics of multibody systems are well understood, direct approaches to resolve the motion of all degrees of freedom for locomotion have remained computationally intractable due to the large number of degrees of freedom and nonlinear behaviour. To generate dynamically balanced locomotion trajectories efficiently, current methods are typically based on simplified models of the robot dynamics.

One of the most widespread model is the Inverted Pendulum Model (IPM) [START_REF] Kajita | The 3D linear inverted pendulum mode: A simple modeling for a biped walking pattern generation[END_REF], where the mass of the robot is assumed to be concentrated at its center of mass (CoM). Although more accurate models are sometimes required [START_REF] Lee | Reaction Mass Pendulum (RMP): An explicit model for centroidal angular momentum of humanoid robots[END_REF], the simple IPM is suitable for many situations where the rotational inertial effects of the robot N. Perrin, D. Lau & V. Padois Institut des Systèmes Intelligents et de Robotique (ISIR), CNRS UMR 7222 & Université Pierre et Marie Curie, Paris, France, e-mail: {perrin, lau, padois}@isir.upmc.fr 1 arms, legs and torso are negligible or can compensate for each other. An example of such a scenario includes walking at a moderately fast pace.

One of the main properties of multi-contact locomotion is that the robot's CoM acceleration depends only on the contact forces with the environment. Hence, an important objective for the generation of locomotion is to design a CoM trajectory that is dynamically balanced at any point in time. Dynamically balanced refers to the existence of contact forces that can produce the desired CoM acceleration while respecting the contact constraints, such as being within the friction cones.

When all the contacts are coplanar, a point of particular interest, sometimes called the Zero Moment Point (ZMP) [START_REF] Vukobratović | Zero-moment point-thirty five years of its life[END_REF], can be defined. Although this name has led to confusion about its physical nature [START_REF] Goswami | Postural stability of biped robots and the foot-rotation indicator (FRI) point[END_REF], the term ZMP will be used in the remainder of the work. A system is considered dynamically balanced if the ZMP lies within the support region, i.e. the 2D convex hull of the contact points. In the IPM, the relationship between the CoM and ZMP dynamics is defined by nonlinear differential equations. However, if the vertical displacements of the CoM are set in advance, the relationship can be decoupled and expressed as linear differential equations. This substantially simplifies the problem into one that can be more efficiently and easily implemented for locomotion trajectory generation and gait control.

The trajectory generation problem has been studied through analytical approaches [START_REF] Harada | An analytical method on real-time gait planning for a humanoid robot[END_REF], and more recently, using convex optimization with constraints on the ZMP in discrete time [START_REF] Herdt | LMPC based online generation of more efficient walking motions[END_REF][START_REF] Wieber | Trajectory free linear model predictive control for stable walking in the presence of strong perturbations[END_REF][START_REF] Wieber | Viability and predictive control for safe locomotion[END_REF][START_REF] Herdt | Walking without thinking about it[END_REF]. The efficiency of convex optimization solvers allows them to be used within model predictive control (MPC) schemes and with the potential of real-time implementation for reactive walking. However, the ZMP+MPC approach suffers from two main drawbacks: 1) the trajectory of the CoM height (zdirection) must be known or fixed in advance, and 2) the contact points must always be coplanar, making this approach unsuitable for walking in complex unstructured environments or if the arms are to also be used. Due to these two restrictions, the ZMP+MPC approach can be regarded as a 2D CoM planner that operates in the horizontal xy-plane only.

Several studies have looked into extending the concept of ZMP into 3D conditions, such as GZMP [START_REF] Harada | Dynamics and balance of a humanoid robot during manipulation tasks[END_REF] and 3DZMP [START_REF] Inomata | 3DZMP-based control of a humanoid robot with reaction forces at 3-dimensional contact points[END_REF], to handle non-coplanar contact points. These criteria have been used in control algorithms to maintain the dynamical balance of a humanoid robot interacting with its environment. However, no locomotion trajectory generation algorithm has been designed based on these notions. In [START_REF] Hirukawa | A universal stability criterion of the foot contact of legged robots -adios ZMP[END_REF], a more general criterion than the ZMP was proposed to evaluate the balance of contacts during the motion of a legged robot. This criterion was used within a preview control algorithm with additional restrictions, for example, the vertical motion of the CoM must be approximately constant. Furthermore, a preliminary phase is needed to plan the inertia and gravity wrenches appropriately, which is a difficult problem.

In this paper, two simple and novel MPC approaches to solve for 3D locomotion with multiple non-coplanar contacts are presented. The 3D condition for dynamically balanced gait allows for non-coplanar multiple contacts and no restrictions on the CoM height trajectory. Using the proposed 3D dynamically balanced criterion, the first MPC formulation treats the criterion as an objective function, where the resulting non-convex MPC problem is solved using a sequence of alternating (con-vex) quadratic programs (AQP). The second formulation considers the criterion as non-convex constraints to the problem, and is solved through a succession of convex QCQPs. The results for simple locomotion scenarios show the promise in using the proposed 3D condition to generate the CoM trajectory within the control framework of robots with multiple contacts.

Although generalizations, such as allowing for multiple contacts, non-coplanar contacts and not predetermining height trajectory, inherently increase the problem complexity, care has been taken to maintain balance between the computational efficiency and limitations of the model. Indeed, our approach is more general than traditional simplified methods based on the IPM or ZMP, but thanks to our main assumption, namely, that the CoM lies inside all the contact friction cones, the number of variables that have to be taken into account is much smaller than with direct methods, such as [START_REF] Posa | A direct method for trajectory optimization of rigid bodies through contact[END_REF]. The feature of the proposed MPC schemes is that they can be effectively and efficiently solved through a succession of convex optimization problems, when the problem formulation is of particular forms, such as bilinear problems or non-convex QCQPs.

The remainder of the paper is organised as follows: Section 2 formulates the 3D condition for dynamically balance. The two MPC formulations using the 3D condition as an objective or as constraints are formulated in Sections 3 and 4, respectively. Section 5 presents a brief discussion on the two approaches. Finally, Section 6 concludes the paper and presents areas of future work.

Conditions for Dynamically Balanced Locomotion in 3D

Expressing wrenches with respect to the CoM of the robot x, the Newton-Euler equations of motion for a multibody robot system in a fixed world frame can be written as

W gravity x + W contact x = L M ẍ , (1)
where W gravity x is the gravity wrench in 6D vector notation [START_REF] Featherstone | Rigid Body Dynamics Algorithms[END_REF], W contact x the sum of the contact wrenches, L is the angular momentum of the whole robot with respect to its CoM, and M the total mass of the robot. The gravity wrench is equal to the vector 0 T Mg T T , where g is the gravity vector, and

W contact x = N ∑ j (c j -x) × f j f j , (2)
where the c j and f j are the location and force of contact j, respectively, and N is the total number of contacts. Substituting (2) into (1) results in

∑ j (c j -x) × f j f j = L M(ẍ -g) . (3)
The system can be regarded as dynamically balanced if (3) is satisfied. The goal of the problem of CoM trajectory generation is to compute ẍ(t) such that the dynamic wrench, the right hand side term of (3), can be compensated by the sum of contact force wrenches. It is usually assumed that the contact forces belong to a cone (friction cone) spanning from the contact point.

Generally, it is not trivial to compute the set, or even a reasonable subset, of the dynamic wrenches that can be compensated by a given set of contact cones [START_REF] Trinkle | On dynamic multi-rigid-body contact problems with coulomb friction[END_REF]. However, by assuming that all contact cones contain the CoM x, as shown in Fig. 1, it is possible to choose the force of contact j to be

f j = α j (x -c j ) , α j ≥ 0 . ( 4
)
This assumption is equivalent to enforcing the constraint xc j ∈ F j , where F j is the friction cone for the contact point j. Given the contact point c j , normal vector to the contact surface n j and the coefficient of friction µ j , this constraint can be expressed as a second order cone constraint or approximated by a set of linear constraints.

Fig. 1 All the contact cones contain the direction towards the position of the robot CoM x. Choosing the force of contact j equal to α j (c j -x), α j ≥ 0, yields: M(gẍ) = ∑ j α j (c j -x), α j ≥ 0.

By this selection, substituting (4) into the criterion for dynamic balance (3) results in

0 ∑ α j (x -c j ) = L M(ẍ -g) . (5)
From (5), the system is dynamically balanced if:

1. The angular momentum L of the whole robot with respect to its CoM is negligible. 2. M(gẍ) can be expressed by a positive linear combination of the vectors c j -

x , j = 1, . . . , N.

It is worth noting that the proposed condition (5) holds in 3D without the restriction that contacts are coplanar. Furthermore, although the assumption xc j ∈ F j limits range of scenarios in which the simplified conditions for dynamic balance (5) can be used, it is important to note that the restrictions are less limited and can be applied in a larger range of scenarios compared with Zero Moment Point (ZMP) conditions.

3 Model Predictive Control with Dynamically Balancing Objective

Problem Formulation: Non-Convex Optimization Problem

The aim of this problem is to determine the CoM trajectory for the locomotion of a multi-limbed robot for a given sequence of contacts over a finite time horizon.

Although the number of contacts between the robot and the environment can theoretically be infinite, without loss of generality it can be supposed that only a finite number of points N of the robot can be in contact with the environment. For example, the contact points for a rectangular foot can be represented by the 4 vertices of the foot sole. The positions of these time-varying contact points can be denoted by c 1 (t), c 2 (t), . . . , c N (t). A binary variable κ j (t) ∈ {0, 1} can be defined to denote whether c j (t) is in contact (κ j (t) = 1) or not in contact (κ j (t) = 0) with the environment at time t. In this work, the contact trajectory c j (t) and κ j (t) are predetermined.

In this proposed model predictive control (MPC) scheme, the decision variables are the CoM positions x(t) and the contact force multipliers α j (t) ∀ j ∈ {1, . . . , N}, over the finite time horizon T , where t ∈ [t 0 ,t 0 + T ]. Discretizing the horizon at a time step of δt results in K discrete time instances, where T = Kδt and t i = t 0 + i • δt. Given the contact trajectory information c j (t) ∀ j ∈ {1, . . . , N} and κ j (t) ∀ j ∈ {1, . . . , N}, the objective is to satisfy the translational components of (5) by minimising

J 1 = K ∑ i=1 M(g -ẍ(t i )) - N ∑ j κ j (t i )α j (t i )(c j (t i ) -x(t i )) 2 .
The constraints for this problem are:

• Initial/current CoM position x(t 0 ) = x 0 • Initial/current CoM velocity ẋ(t 0 ) = ẋ0 • Final CoM position at the end of the horizon x(t K ) = x f • Max acceleration ẍ(t i ) ∞ ≤ a max , i = 0, . . . , K -1
Note that any constraints on the CoM velocity ẋ and acceleration ẍ can be expressed with respect to x by the linear relationships

ẋ(t i ) = x(t i ) -x(t i-1 ) δt , i = 1, . . . , K (6)
ẍ(t i ) = ẋ(t i+1 ) -ẋ(t i ) δt , i = 0, . . . , K -1 . (7)
As a result, all constraints are either linear equality or inequality constraints with respect to x. Thus, the MPC optimization problem to determine the CoM trajectory over the finite time horizon can be expressed as min

(x(t i )) i ,(α j (t i )) i, j J 1 (x(t i )) i , (α j (t i )) i, j s.t. ẍ(t i ) ∞ ≤ a max , i = 0, . . . , K -1 x(t 0 ) = x 0 ẋ(t 0 ) = ẋ0 x(t K ) = x f . (8)
Remark: as it is presented here, to implement the approach as a model predictive controller would require measuring the CoM velocity (ẋ 0 ), which is not always easy. If the contact forces are accessible, one could estimate their total wrench and use it to obtain an approximation of M(gẍ), which could potentially be used to make the approach more robust.

An algorithm exploiting the objective function structure

On a computational level, the objective function J 1 is non-convex and it is therefore difficult to obtain the global optimum. Finding a local minimum could be relatively quick and provide good results, but in the next section we propose another approach that exploits the form of the objective function. The alternating convex optimization approach is inspired by [START_REF] Majumdar | Control design along trajectories with sums of squares programming[END_REF], where successful results for a different application are obtained by using an alternating sequence of convex QPs (AQPs) instead of trying to solve head-on an optimization problem with a bilinear objective function. The structure of function J 1 can be expressed as

J 1 = ∑ i ε i (x(t i ), (α j (t i )) j ) 2 , (9)
with

ε i = M(g -ẍ(t i )) - N ∑ j κ j (t i )α j (t i )(c j (t i ) -x(t i )).
The functions ε i are bilinear in the variables x(t i ) and α j (t i ). Therefore, if x(t i ) are fixed, the optimization problem (8) becomes a convex QP that can be efficiently solved. Similarly, if the α j (t i ) variables are fixed, we obtain a convex QP in the x(t i ) variables. This property can be used to alternatively optimize the x(t i ) and α j (t i ) variables, using the solution of each step as the fixed variables of the next step. This ensures that the objective J 1 decreases at each iteration of the algorithm until it converges.

For the first step, an initial guess was made for the CoM trajectory that is consistent with the constraints from [START_REF] Hirukawa | A universal stability criterion of the foot contact of legged robots -adios ZMP[END_REF], and the α j (t i ) variables are optimized. The most direct motion from x 0 to x f is often a natural candidate for this initial guess, so in practice this method does not require any tuning at all.

Results

The RobOptim [START_REF] Moulard | RobOptim: an optimization framework for robotics[END_REF] framework was used to illustrate the use of the proposed MPC formulation [START_REF] Hirukawa | A universal stability criterion of the foot contact of legged robots -adios ZMP[END_REF] and alternating quadratic programs that exploit the bilinear form (9) on CoM generation for locomotion. RobOptim provides a convenient interface to try various optimization tools. One interesting aspect of the formulation of the proposed MPC is that the matrices describing the constraints are very sparse, with O(m) non-zero elements where m is the total number of variables. Thus, we chose the IPOPT optimizer [START_REF] Wächter | On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming[END_REF] that can exploit matrix sparsity. Three different scenarios are simulated and presented:

1. 3 step walk with two foot supports (coplanar contacts); 2. 3 step walk with two foot supports and one hand support (non-coplanar contacts); 3. Jump-step with flight phase. The 4.5 s trajectory for scenario 1 can be defined by the following sequence as shown in Fig. 2: • 0 ≤ t < 0.6 s: Robot initially begins in double support • 0.6 ≤ t ≤ 1.2 s: Step 1 moving the left foot • 1.2 < t < 1.9 s: Double support phase • 1.9 ≤ t ≤ 2.5 s: Step 2 moving the right foot • 2.5 < t < 3.2 s: Double support phase • 3.2 ≤ t ≤ 3.8 s: Step 3 moving the left foot • 3.8 < t ≤ 4.5 s: Double support phase With a time horizon of 4.5 s and δt = 0.1 s, the MPC scheme consists of 45 discrete steps. Without loss of generality, it should be noted that the MPC in this setup solves for the entire trajectory, where in practice the time horizon would be shorter to achieve better computational performance. The dimensions of the feet 10.45 cm × 28.3 cm and the total mass 36.66 kg correspond to the Romeo robot [START_REF]Aldebaran: The Romeo project[END_REF]. At t = 0 s and t = 4.5 s, the CoM is set to be in the middle between the feet at a height of 0.6 m, x(0) = [0.14 0 0.6] T and x(4.5) = [0.85 0 0.6] T , respectively. The initial guess for the CoM trajectory follows a straight line from its initial to its final position. Fig. 2 shows the CoM motion generated after two steps of the algorithm (i.e. one optimization of the α j (t i ) variables and then one optimization of the CoM trajectory).

Fig. 2 The resulting trajectory from the MPC scheme (8) for scenario 1. In blue, the horizontal CoM trajectory produced by our algorithm in 95 ms. In red, the projected ZMP trajectory [START_REF] Kajita | The 3D linear inverted pendulum mode: A simple modeling for a biped walking pattern generation[END_REF].

Fig. 2 shows the resulting CoM trajectory (in blue) and the projection (in red)

x zmp = λ zmp (x + M(g -ẍ)), ( 10
)
where λ zmp > 0 is a scaling factor such that z zmp = 0. The projected point ( 10) is equivalent to the ZMP point for the scenario 1 in which all contact points are always coplanar. It can be observed that the projected ZMP lies within the foot during single support phases. In this scenario, the results are similar to that of the ZMP+MPC approach, since all contacts are coplanar on the ground plane. For this scenario, the optimization problem [START_REF] Hirukawa | A universal stability criterion of the foot contact of legged robots -adios ZMP[END_REF] with horizon of K steps consists of 11K decision variables, comprised of 3 for the CoM position and 8 for the two contact supports (four contact points per support) at each time step, and 9 + K constraints from [START_REF] Hirukawa | A universal stability criterion of the foot contact of legged robots -adios ZMP[END_REF]. The trajectory was computed in 95 ms (on a 2.40GHz Intel(R) Core(TM) i7-4700MQ CPU) after solving two alternations of the AQP, and hence 4 convex QPs. Executing more alternations of the algorithm showed that the cost J 1 quickly reached small values, for example, more than 1000 times smaller than the initial value after 4 optimizations. Furthermore, it was observed that the algorithm had almost converged after the 2 first optimizations.

In the second scenario, an additional contact corresponding to the right hand of the robot on a wall is considered. The contact point is at a height of 0.6 m, and activated only during the second step. This example demonstrates the ability of the proposed 3D condition and MPC algorithm to handle non-coplanar contacts. The resulting trajectory in Fig. 3 shows how the proposed algorithm can handle this situation and generated a different CoM trajectory. It could be observed from the results that the additional hand contact point enabled the robot to avoid the sway motion to the left.

Finally, the third scenario illustrates the ability of the algorithm to determine the CoM height without restrictions unlike the classical ZMP+MPC approach. The steps were replaced by one jump defined by the following sequence of actions: at t = 0.6 s the left foot leaves the ground, and at t = 2.0 s the right foot leaves the ground. After a flight phase of 0.4 s, the left foot lands at t = 2.4 s, and then the right foot lands at t = 3.8 s. Remark: in this scenario the a max bound must be at least equal to the norm of the gravitational acceleration. Fig. 4 shows the CoM trajectory produced after 2 optimizations. This example clearly shows the benefits in relaxing the fixed height Fig. 3 The resulting trajectory from the MPC scheme (8) for scenario 2 with hand support. In blue, the horizontal CoM trajectory produced after 2 steps of the algorithm (in 96 ms). In purple, the trajectory produced after 6 steps of the algorithm (in 412 ms). These trajectories are almost exactly the same, which shows that the convergence is fast. trajectory from the ZMP+MPC approach, such that the CoM trajectory generation is able to produce a jump motion in the z-direction. With the same aim of determining the CoM trajectory as Section 3, the problem in this section will be formulated by considering the dynamically balanced condition (5) as constraints over the time horizon. Furthermore, the objective to minimise is the tracking error for a reference CoM trajectory x r (t). The reference CoM trajectory is pre-planned based on the desired motion. In this MPC problem, the decision variables are the CoM jerk vectors ... x (t i ), ∀i ∈ {1, . . . , K}. The convex objective of the problem is to minimise

J 2 = K ∑ i=1 Q 1 x(t i ) -x r (t i ) 2 + Q 2 ... x (t i ) 2 ,
where Q 1 and Q 2 are the relative weights between the tracking and jerk terms, respectively. Note that the relationship between the CoM jerk ...

x and the CoM position x, velocity ẋ and acceleration ẍ is linear [START_REF] Wieber | Trajectory free linear model predictive control for stable walking in the presence of strong perturbations[END_REF] in the form

  x(t i+1 ) ẋ(t i+1 ) ẍ(t i+1 )   = A   x(t i ) ẋ(t i ) ẍ(t i )   + B ... x (t i ) .
As a result, the objective and all constraints can be expressed with respect to the CoM jerk. As presented in Section 2, the geometrical meaning of the constraint

∑ α j (c j -x) = M(g -ẍ) , α j ≥ 0 (11)
is that the vector M(gẍ) must lie within the positive cone of the vectors (c jx) ∀ j ∈ {1, . . . , N}. Such constraints will be derived in the following for the case of one and two supports in contact. As shown in Fig. 5(a), when the system is in single support it will be assumed that there are four contact points. The position of the vertices for the rectangular contact surface relative to the center of the contact c j can be represented by a 1 , a 2 , a 3 and a 4 . For the single support, the vector M(gẍ) is inside the cone produced by the vectors (c j + a 1 -x, c j + a 2 -x, c j + a 3 -x, c j + a 4 -x) = (w 1 , w 2 , w 3 , w 4 ) if the following triple products are positive:

w 1 × w 2 • (g -ẍ) ≥ 0, w 2 × w 3 • (g -ẍ) ≥ 0 w 3 × w 4 • (g -ẍ) ≥ 0, w 4 × w 1 • (g -ẍ) ≥ 0 . ( 12
)
For the case of double support contacts, the expression of the convex cone would be non-trivial if each support was assumed to consist of 4 contact points. This is because the convex cone would be dependent on the locations of each support. Furthermore, this would increase the number of surfaces of the convex cone for the vector M(gẍ) to be checked, hence increasing the number of constraints. However, by assuming that each contact support only has two contact points, at the top a t and bottom a b , then the resulting convex cone shown in Fig. 5(b) only consists of four surfaces and is a strict subset of the 4 contact point convex hull regardless of the support location. For the double support, the vector M(gẍ) is inside the cone produced by the vectors

(c L + a t -x, c L + a b -x, c R + a b -x, c R + a t -x) = (w 1 , w 2 , w 3 , w 4 )
as with the conditions in [START_REF] Majumdar | Control design along trajectories with sums of squares programming[END_REF]. In a similar manner, constraint equations for extra scenarios with additional contacts could be considered. For both mathematical and practical considerations, it will be assumed that contacts do not overlap, such that the contact points in Figure 5(b) can form a convex cone.

It can be observed that the constraints ( 12) for both single and double support are non-convex quadratic inequality constraints if the contact trajectories are known. As a result, the MPC optimization problem to determine the CoM trajectory over the finite time horizon can be expressed as min

( ... x (t i )) i K ∑ i=1 Q 1 x(t i ) -x r (t i ) 2 + Q 2 ... x (t i ) 2 s.t. ( 12
) for κ j = 1, κ a = 0 ∀a = j (12) for κ L = 1, κ R = 1 . (13)
The optimization problem ( 13) is a non-convex quadratically constrained quadratic program (QCQP) with a convex objective function.

Feasible Point Pursuit Successive Convex Approximation

On a computational level, the objective function J 2 has non-convex constraints and therefore it is non-trivial to obtain a feasible solution, let alone the global optimum. The feasible-point-pursuit successive convex approximation (FPP-SCA) [START_REF] Mehanna | Feasible point pursuit and successive approximation of non-convex QCQPs[END_REF] is an effective approach that solves the non-convex QCQP as a succession of QP convex approximations. Considering the non-convex QCQP of the standard form with the decision variable

u ∈ R n min u u T A 0 u + q T 0 u s.t. u T A k u + q T k u ≤ b k , k = 1, . . . , L , (14)
where A 0 ∈ R n×n is a positive semidefinite matrix, if any of A k ∈ R n×n from the L constraints are not positive semidefinite, then ( 14) is non-convex and the problem is NP-hard in general. By performing a linear restriction about any point z i , the FPP-SCA approach aims to solve the following problem min u,s

u T A 0 u + q T 0 u + λ ∑ k s 2 k s.t. u T A + k u + (2z T i A - k + q T k )u ≤ b k + z T i A - k z i + s k , s k ≥ 0, k = 1, . . . , L (15)
where A + k and A - k are the positive semidefinite and negative semidefinite matrices from the decomposition

A k = A + k + A - k .
The terms s k are slack variables that represent the violation of the k-th quadratic constraint and λ 1 is a constant that gives priority to minimise the constraint violation. The vector z i ∈ R n is the initial guess vector at iteration i. In the FPP-SCA approach, the convex QCQP ( 15) is repeated, where the variable z i+1 is set as the optimal solution u * i at iteration i. By using the FPP-SCA approach to solve the MPC problem [START_REF] Mehanna | Feasible point pursuit and successive approximation of non-convex QCQPs[END_REF] with horizon of K steps, at each iteration of FPP-SCA the convex QCQP problem consists of 7K decision variables, comprised of 3 for the CoM jerk and 4 constraint slack variables at each time step, and 4K quadratic constraints [START_REF] Majumdar | Control design along trajectories with sums of squares programming[END_REF] at each time step. The initial guess z 0 can be any arbitrary vector, and in the simulations the zero vector was chosen. By solving a succession of (15), it was shown in [START_REF] Mehanna | Feasible point pursuit and successive approximation of non-convex QCQPs[END_REF] and in the results of Section 4.3 that the algorithm converges within a few successions.

Results

To demonstrate the MPC formulation [START_REF] Mehanna | Feasible point pursuit and successive approximation of non-convex QCQPs[END_REF], the CoM jerk trajectories were determined using the FPP-SCA for the three following scenarios: scenario 1 from Section 3.3, scenario 2 from Section 3.3, and walking up 2 steps of a staircase with the same hand support as in scenario 2 from Section 3.3. Scenarios 1 and 2 allow for a direct comparison between the trajectories resulting from the two different MPC schemes. The number of successive QCQPs solved was set as 5, however it was observed that convergence typically happened in less than 3 successions. The reference trajectory for both scenarios was set to be simply the forward linear motion of travel at a height of 0.6 m. The values for objective function weights in ( 13) and [START_REF] Posa | A direct method for trajectory optimization of rigid bodies through contact[END_REF] were set as Q 1 = 1, Q 2 = 10 -4 and λ = 10 4 .

The trajectory for scenario 1 is shown in Fig. 6. From the projected ZMP points using [START_REF] Kajita | The 3D linear inverted pendulum mode: A simple modeling for a biped walking pattern generation[END_REF] it is clear that the produced CoM trajectory satisfies the 3D dynamically balanced criterion. To maintain good tracking performance of the reference Step 1

SS phase SS phase

SS phase

Step 3

Step 2

Fig. 6 The resulting trajectory from the MPC scheme (13) for scenario 1. In blue, the horizontal CoM trajectory and in red, the projected ZMP trajectory [START_REF] Kajita | The 3D linear inverted pendulum mode: A simple modeling for a biped walking pattern generation[END_REF] since all contacts are coplanar. Fig. 7 On the left: the resulting trajectory from the MPC scheme (13) for scenario 2 with hand support. In blue, the horizontal CoM trajectory and in red, the projected ZMP trajectory [START_REF] Kajita | The 3D linear inverted pendulum mode: A simple modeling for a biped walking pattern generation[END_REF]. It can be observed that the results are similar to that from Fig. 3. On the right: the resulting z-direction trajectory from the MPC scheme (13) for scenarios 1 and 2. In black, the CoM trajectory and in red, the reference height trajectory. The results show that the trajectory successfully lowers from the initial height z 0 = 0.75 m to the reference height 0.6 m.

trajectory for y = 0, it is expected that the projected ZMP would be as close to the boundary of the single stance (SS) support region as possible. It can be observed that this MPC scheme generated very similar CoM trajectory results from Fig. 2. As with Fig. 2, the natural artifact of left-right swaying motion can be observed in Fig. 6 due to the existence of single support instances. As such, a hand contact was included (scenario 2) during step 2 of the motion. From the resulting trajectories shown in Fig. 7 (on the left), the MPC scheme generated a very similar, in fact near identical, behaviour to that in Fig. 3. The extra hand contact allowed the CoM to better track the reference trajectory as there is no single support phase sway required during step 2. To show the tracking of the height trajectory, the initial height of the CoM was set to be 0.75 m. Fig. 7 (on the right) shows that the MPC CoM planner was able to quickly converge to the desired reference height of 0.6 m. As such, it can be claimed that for this scenario, the FFP-SCA formulation was able to achieve both feasibility of the dynamically balancing constraints and good tracking performance.

Finally, scenario 3 demonstrates the ability for the proposed MPC formulation to solve more complex behaviours, such as walking up stairs while using the hand to hold onto the staircase rails. The steps of the stairs were set as 15 cm high and the reference height trajectory for the CoM was set to be 0.6 m above the surface of stair steps. Fig. 8 (on the left) shows the x and y direction (top-down view) CoM trajectories determined by the MPC. It can be observed that the CoM results for the x and y directions are very similar to that of walking on a flat ground as shown in Fig. 7 (on the left). However, the resulting z trajectory shown in Fig. 8 (on the right) shows significant differences with scenario 2. In addition to satisfying the dynamically balanced constraints for the locomotion on the non-coplanar contacts, the desired height trajectory represented by the red line in Fig. 8 (on the right) was able to be tracked. This example shows the robustness of the formulation and the QCQP solver to various different scenarios.

Discussion

To demonstrate the advantages and use of the 3D model for dynamically balanced locomotion in Section 2, two different example ways to use the criterion for MPC generation of CoM trajectory were presented. Both approaches have a fundamentally common point: the generalisations to allow for 3D non-coplanar multiple contacts naturally result in non-convex problems. However, regardless of the use of the criterion as an objective or constraint, the nature of the criterion is that it is in a bilinear form that can be converted into a convex quadratic function by restricting some variables. Both MPC approaches take advantage of this and then solve a succession of convex optimization problems. This is important in the proposed MPC schemes to ensure that it is still possible to implement them on a robot in real-time. Comparing between the two MPC schemes, it can be observed that using the dynamically balancing criterion in the objective function results in a problem with more decision variables than treating it as a constraint (11K vs. 7K, where K is the horizon length). However, the number of constraints is significantly less in return (9 + K vs. 4K), and the resulting problem is only a QP rather than a QCQP. As a result, the approach from Section 3 is expected to be more computationally efficient than the one from Section 4. But the constraint MPC approach provides a stricter notion of feasibility to dynamically balanced locomotion, and is less concerned with optimality. Finally, it is also worth noting that compared to the traditional ZMP+MPC approach, several restrictions have been removed, such as coplanar contacts and predetermined height trajectory. The interesting point is that if any of these are relaxed, the problem complexity is identical to that if all are relaxed. As such, the 3D formulation proposed relaxes many conditions from the ZMP+MPC approach while still maintaining a balance with the computational cost of the resulting method. As with the ZMP+MPC approach, the focus is typically more on generating dynamically balanced motion rather than optimal gait behaviour. Hence, the development of methods such as AQP and FPP-SCA provides the opportunity to generate feasible motion for more general locomotion scenarios in real-time control.

Conclusion and Future Work

We proposed a novel model for dynamically balanced trajectory generation, more general than the classical IPM+ZMP approach, but simple enough to enable fast computations of CoM trajectories through an iterative resolution of convex QPs or convex QCQPs. This claim is supported by the low number of decision variables and constraint equations shown in the problem analysis. The generalizations gained from the proposed model and MPC approach include the ability to allow for multiple non-coplanar contacts and not having to predefine the CoM height trajectory. The results of the two proposed MPC approaches support the belief that the proposed 3D model of dynamically balanced locomotion is a good candidate for real-time model predictive control for multi-contact locomotion.

In future work, we will focus on performing experiments on a real humanoid robot, and address the following points:

1. Both the AQP and FPP-SCA approaches are observed to work well in practice and converge quickly. However there is no mathematical guarantee on the optimality of the solution, hence better understanding and analysis of such methods on the particular structure of the proposed MPC formulations should be more precisely studied. 2. In the optimisation of the CoM trajectories, the contact locations c i could also be optimized without losing the convex QP structure. This allows the potential to not only compute the CoM trajectory, but also optimised contact locations.

Acknowledgements

The research presented in this paper was partially funded by the ROMEO2 project.

Fig. 4

 4 Fig. 4 On the left: the horizontal CoM trajectory generated (in 176 ms) for the jump scenario. On the right: the evolution of the CoM coordinates during the flight phase.

Fig. 5

 5 Fig.5The constraints required to achieve dynamic balance for: (a) when one support is in contact with the environment κ L = 1, κ R = 0 and (b) when two supports are in contact κ L = 1, κ R = 1.

Fig. 8 7 .

 87 Fig.8On the left: the resulting trajectory from the MPC scheme (13) for scenario 3 with hand support walking up a staircase. The horizontal CoM trajectory is shown in blue. It can be observed that the results are similar to that from Fig.7. On the right: the resulting z-direction trajectory from the MPC scheme (13) for scenario 3. In black, the CoM trajectory and in red, the reference height trajectory. The blue lines show the left and right foot trajectories indicating that the robot is walking up the staircase. The results show that the trajectory successfully tracks the desired height above the staircase platform.