Peter E Raad
email: praad@smu.edu

Pavel L Komarov

Mihai G Burzo

An Integrated Experimental and Computational System for the Thermal Characterization of Complex Three-Dimensional Submicron Electronic Devices

The work presents the creation of an analysis engine and experimental system capable of fully characterizing the thermal behavior of complex three-dimensional active submicron electronic devices. First, the transient surface temperature field of pulsed devices is noninvasively measured with submicron spatial resolution. Next, the thermal conductivity of each thin-film layer composing the device is measured and a numerical model is built using these values. The temperature distribution map is then used as input for an ultrafast inverse computational solution to fully characterize the thermal behavior of complex three-dimensional devices. By bringing together measurement and computation, it becomes possible for the first time to non-invasively extract the transient three-dimensional thermal behavior of nanoscale embedded features that cannot otherwise be accessed. The power of the method is demonstrated by verifying that it can extract details of interest of specially constructed micro-resistors.

INTRODUCTION

The determination of temperature fields is crucial for the design, reliability analysis, and predictions of performance of microelectronic devices. As a result, there is an increased need for methods to determine the temperature of features at the submicron level [START_REF] Goodson | Short-time-scale Thermal Mapping of Microdevices using a Scanning Thermoreflectance Technique[END_REF][START_REF] Grauby | High Resolution Photothermal Imaging of High Frequency Using Visible Charge Couple Device Camera Associated with Multichannel Lock-in Scheme[END_REF]. Computational AN INTEGRATED EXPERIMENTAL AND COMPUTATIONAL SYSTEM … approaches can provide insight into the internal thermal behavior of such complex devices, but can be limited by the inherent necessity of modeling the heat sources, which in the case of selfheating microelectronic devices, are the result of electrical fields whose exact shapes and locations are difficult to specify with reasonable certainty. Moreover, such devices can actually experience irreversible changes in thermo-physical properties and/or geometries that cannot be otherwise predicted from theory or monitored. Experimental approaches can also be helpful in determining thermal behavior, but require either physical access or a visual path to the region of interest. Contact methods, for example, present the difficulties of having to access features of interest with an external probe, or in the case of embedded features, fabricate a measuring probe into the device, and then having to isolate and exclude the influence of the probe itself. Noncontact methods, on the other hand, can provide surface temperature profiles, but in and of themselves cannot impart information on internal behavior. In other words, these methods provide a two-dimensional perspective on what otherwise is, in the case of stacked complex devices, an intricate three-dimensional thermal behavior.

However, by combining an experimental method capable of mapping the surface temperature of a complex device with spatial and temporal accuracy together with a computational engine capable of rapidly and accurately resolving the geometric and material complexities of a full three-dimensional microelectronic device, it becomes possible to use the independent information from the experimental measurements to mitigate the lack of knowledge in the source model parameters, which directly affect the usefulness of the computational results. This article introduces a proof of concept of a coupled computational-experimental approach that uses a measured two-dimensional surface temperature mapping to help obtain a fully three-dimensional thermal characterization of an active micro-device.

METHODOLOGY

The overall approach combines computational and experimental methods previously developed by the authors. The transient two-dimensional surface temperature is measured by the use of the thermoreflectance thermography system [START_REF] Komarov | Thermal Characterization Of Pulse-Activated Microelectronic Devices By Thermoreflectance-Based Surface Temperature Scanning[END_REF], while the three-dimensional thermal behavior of multi-layered integrated circuits (with embedded features) is inferred by solving the inverse heat transfer problem with the self-adaptive ultra-fast numerical technique [START_REF] Raad | System and Method for Predicting the Behavior of a Component[END_REF][START_REF] Wilson | A Transient Self-Adaptive Technique for Modeling Thermal Problems with Large Variations in Physical Scales[END_REF]. To minimize the number of uncertainty sources that the inverse method must deal with, the thermophysical properties of the various thin layers are measured independently with the Transient AN INTEGRATED EXPERIMENTAL AND COMPUTATIONAL SYSTEM … Thermo-Reflectance technique [START_REF] Burzo | A Study of the Effect of Surface Metallization on Thermal Conductivity Measurements by the Transient Thermo-Reflectance Method[END_REF][START_REF] Komarov | Transient thermo-reflectance measurements of thermal conductivity and interface resistance of metallized natural and isotopically-pure silicon[END_REF][START_REF] Komarov | Performance Analysis of the Transient Thermo-Reflectance Method for Measuring the Thermal Conductivity of Single Layer Materials[END_REF]. In what follows are the details of each of these three methods.

Thermo-Reflectance Thermography System

The experimental temperature mapping system is based on the thermoreflectance (TR) method, where the change in the surface temperature is measured by detecting the change in the reflectivity of the sample. Since the change in reflectivity (i.e., thermo-reflectance coefficient, C th ) is of the order of 10 -4 -10 -5 K -1 for most electronic materials, the system had to be designed, built, and fine-tuned to achieve the levels of uncertainty that are required for the success of the measurements presented here. The C th is extremely sensitive to the wavelength of the laser and varies strongly from one material to another. Hence, in order to maximize the signal to noise ratio, we modulated the probing laser light and used a lock-in amplifier to extract the useful signal.

The measurement methodology requires two steps. First, the coefficient of thermal reflectance is determined for each of the surface materials to be scanned (calibration). Second, the changes in the surface reflectivity sample's surface. To overcome the inherently low signal to noise ratio, the activation voltage of the device is modulated, resulting in a modulated photodetector signal that can more easily yield the useful signal from the raw photodetector signal output. The photodetector signal, containing the change in surface reflectivity caused by the temperature variations of the DUT, is acquired with a lock-in amplifier (or an oscilloscope) and scaled according to the calibrated data. The limitations of the lock-in and oscilloscope will be discussed.

The calibration approach consists of determining the relationship between the changes in reflectance and surface temperature. The change in reflectance will be measured by a differential scheme involving two identical PDs in order to minimize the influence of fluctuations in the energy output of the probing laser. The sample temperature will be controlled by a thermoelectric (TE) element and measured with a thermocouple. The calibration must be performed for each of the materials on the surface of each device where a mapping of the temperature is carried out. The temperature field of the desired area of the activated device can be mapped by repeating the procedure at multiple physical locations, which is achieved by precisely moving the probing head with submicron resolution.

Computational Engine

The numerical engine is capable of simulating the transient thermal behavior of active multilayered devices whose dimensions vary over several orders of magnitude and where the thermophysical properties of the materials used are not isotropic. The thermal modeling engine is used to for solving the required inverse heat transfer problem. The measured experimental surface temperature data will be used as input. Since the heat transfer problem is parabolic in nature, a unique temperature field can be computed in the entire structure by solving the corresponding heat transfer problem.

The novel approach begins by solving the corresponding steady-state problem by the use of a grid nesting technique. The nesting technique defines a template that is then used to solve the transient problem in a multiple grid fashion. The advantage to solving with multiple grids in time is that the majority of the problem domain can be resolved in time at lower nest levels (with a coarser grid), while the finer grid resolution is reserved for the parts of the problem that demand the finer resolutions in space and time.

Since the physical dimensions of the various materials used in modeling high performance electronic devices vary greatly, a uniform mesh that resolves all of the details in three AN INTEGRATED EXPERIMENTAL AND COMPUTATIONAL SYSTEM … dimensions results in a prohibitively large computational grid. A common method for dealing with dimensional variation is to skew the mesh and concentrate more grid points in areas where higher resolutions are needed. The shortcoming of using a biased-mesh approach to resolve the geometry is that the problem geometry, and not the temperature gradients, will end up dictating the meshing. The meshing strategy used in the development of the novel numerical technique [START_REF] Raad | System and Method for Predicting the Behavior of a Component[END_REF] was set on ensuring that the method is (i) automatic and adaptive, (ii) independent of user expertise, and (iii) independent of materials (including air), geometry features, embedded vias, and heat source locations. The approach makes it possible to start with the full 3D geometry of a real device in all of its complexity, and then uses a physics-based automatic error predictor to focus the entire available computational power on only those regions that require further refinement in order to achieve the level of acceptable error prescribed by the analyst [START_REF] Wilson | A Transient Self-Adaptive Technique for Modeling Thermal Problems with Large Variations in Physical Scales[END_REF]. The new approach allows for the addition of the difficult problem geometry associated with top surface features of an IC over a large area, and then letting the error prediction technique decide if the regions need further refinement. The strength (and novelty) of the method is that it uses effective thermal properties that are consistent with the local grid spacing at the particular grid level in use. As a result, dealing with air, embedded vias, and ultra-thin multi-layered structures requires no special treatment For complex computationally demanding problems, such is the case of complex electronic ICs one of which being presented in Fig. 2, the approach affords a two order of magnitude increase in speed over conventional methods. The benefits are even more significant if one considers the facts that, in conventional methods, (i) significant effort must be expended by an experienced analyst on mesh generation and optimization; (ii) several mesh iterations are normally required to ensure that features of interest are well resolved; and (iii) several solutions with increasing resolution are required to ensure that grid convergence is achieved. When the most complex case is solved on the newest personal computers, the simulation time is further reduced into the minutes.

Thermal Properties Measurements System -Transient Thermoreflectance (TTR)

Since it has been proven [START_REF] Mihai | Thermal Transport Properties of Gold-Covered Thin-Film Silicon Dioxide[END_REF] that the thermal properties of thin-films varies from that of bulk materials it is necessary for the numerical simulation to input the real value of the thermal conductivity for all of the materials composing the system under study. The transient thermoreflectance measurement system is used to determine the real properties of the thin-film materials to be used in the numerical simulation.

The schematic in Fig. 3 depicts the heating and probing spots produced by a typical TTR system. The source of energy in the TTR method is provided by a pulsed laser with a short pulse duration. During each pulse, a given volume below the sample surface heats up to a temperature wavelength is 532 nm, pulse width is 6.1 ns, and maximum pulse energy is 0.5 mJ. The CW probing laser is an He-Ne with a 633 nm wavelength.

RESULTS

By the use of the integrated approach presented above we are able to solve the inverse conduction problem associated with a complex, multi-layered, deep submicron IC in order to infer the thermal behavior of the embedded features that cannot be otherwise achieved. To prove the concept, we begin by applying it to a simple device type for which the heat source shape is controllable and well defined. The inverse solution will focus on determining a single unknown, namely the depth of the heat source from the device surface. Therefore, we constructed a series of Aluminum micro-resistors, embedded in a layer of silicon dioxide at different depths from the surface. The actual depth of each device is measured independently with an ellipsometer and a stylus profiler. The geometry of the resistor is shown schematically in Fig. 4. A gold pad is deposited on top of the oxide layer to enable probing of the surface temperature (the oxide layer is transparent to the wavelength of the probing light irradiation).

The integrated experimental-computational approach will then be used to infer the depth of the heater by first scanning the surface temperature and then using the numerical simulation to solve the inverse problem for the unknown depth. Because in this case the problem involves only one unknown, the method will use a single-parameter optimization approach that minimizes the RMS error between the measured and computed surface temperature fields. The validity and power of the method will thus be demonstrated by comparing the depth values inferred from applying the coupled technique with those measured by the ellipsometer and profiler. The determination of the correct heat source depth yields the correct, complete, three-dimensional temperature field.

For existing devices, the highly resolved and accurate picture of the 3D temperature field would provide the ability to detect hot spots, diagnose performance, and assess reliability. In design and manufacturing of new devices, the new tool has the potential to provide a rapid approach for analyzing the thermal behavior of complex stacked structures, identify regions of excessive heat densities, and ultimately contribute to improved thermal designs, better device reliability, and shorter design cycle time. In the long-term, this tool can also be integrated with existing electrical characterization tools to provide the ability to design devices while concurrently comprehending electrical and thermal issues. The outcome of this work will contribute to dealing with the critical aspects facing the electronics industry and that are brought about by continued miniaturization, introduction of novel materials, and never-ceasing demand for higher performance.

  as a function of changes in temperature are measured at each point of interest with submicron spatial resolution and better than microsecond temporal resolution. The resulting reflectivity waveforms are combined to obtain a transient temperature field over the scanned area. A schematic of the proposed TTRG is shown in Fig. 1. The probing beam reflects from the heated surface back along the optical path to the sensitive area of a photodiode. The intensity of the reflected light depends on the reflectivity (temperature) of the

Fig. 1

 1 Fig. 1 Surface temperature mapping system

Fig. 2

 2 Fig. 2 Schematic of a representative complex IC problem

  photodetector that reads the instantaneous surface reflectivity. The changes in surface reflectivity are linearly proportional to the changes in surface temperature, within a wide but finite temperature range [10,11]. The specific system in the NETS Laboratory at SMU (http://engr.smu.edu/netsl) uses an Nd:YAG pulsed heating laser whose

Fig. 3

 3 Fig. 3 Schematic of the heating and probing spots on a sample for a typical TTR system

Fig. 4

 4 Fig. 4 Geometry of the scanned micro-resistor

 TIMA EDITIONS / THERMINIC 2005 ISBN 2-916187-01-4


Lorenzo Codecasa
email: codecasa@elet.polimi.it

Dario D ' Amore
email: damore@elet.polimi.it

Paolo Maffezzoni

Structure Function Representation of Multi-Directional Heat-Flows

A relation is established between the structure function of the RC transmission line modeling a one-port passive distributed thermal network and the spatial distribution of thermal properties in heat diffusion problems with generic multidirectional heat-flows.

Introduction

Distributed thermal networks are widely used for modeling heat conduction problems in components and packages. Structure functions, originally introduced by V. Szekely et al. [START_REF] Székely | Fine Structure of Heat Flow Path in Semiconductor Devices: a Measurement and Identification Method[END_REF], are means for characterizing the responses of one-port distributed thermal networks in terms of equivalent RC transmission lines. In fact, as recently precised in [START_REF] Codecasa | Canonical Forms of Thermal Impedances[END_REF], the response of a one-port passive distributed thermal network modeling a generic thermal problem coincides with the response of a proper RC transmission line defined by a structure function.

Since their introduction, structure functions have been used for inferring on the spatial distributions of thermal properties in components and packages [START_REF] Székely | Fine Structure of Heat Flow Path in Semiconductor Devices: a Measurement and Identification Method[END_REF][START_REF] Székely | A new evaluation method of thermal transient measurement results[END_REF][START_REF] Szekely | Restoration of physical structures: an approach based on the theory of RC networks[END_REF][START_REF] Rencz | Measuring Partial Thermal Resistances in a Heat-Flow Path[END_REF]. For a one-dimensional thermal problem in which power is injected at one boundary and temperature rise is measured at the same boundary, a one-port passive distributed thermal network naturally arises. In this case, the cumulative resistance and capacitance of the RC transmission line are exactly the cumulative thermal resistance and capacitance from the boundary at which power is injected. In this case, hereafter referred to as one-directional heat flow, by determining the structure function from the thermal response of the distributed thermal network, precise information on the spatial distribution of thermal properties is recovered [START_REF] Katchalov | verse Boundary Spectral Problems[END_REF].

No procedure is reported in literature for exactly determining a structure function from a thermal response. Moreover in practice, this technique is used in applications which are only approximated by one-directional heat-flows. The more inaccurate is this approximation, the more inaccurate is the recovered information on the spatial distribution of thermal properties. That is because no exact interpretation of the structure function is available. The situation is even more critical for a generic thermal problem, hereafter referred to as multi-directional heat-flow. For such a thermal problem, the structure function can still be associated to a one-port passive distributed thermal network, as precised in [START_REF] Codecasa | Canonical Forms of Thermal Impedances[END_REF]. However the relation between the structure function and the spatial distribution of thermal properties in the thermal problem is not known.

In this paper, which extends the results in [START_REF] Codecasa | Physical Interpretation and Numerical Approximation of Structure Functions of Components and Packages[END_REF], the structure function of a one-port passive distributed thermal network modelling a multidirectional heat-flow is considered. A procedure for exactly determining a structure function from a thermal response is shown. Moreover an exact relation between the structure function of the equivalent RC transmission line on one hand and the spatial distribution of thermal properties in the thermal problem on the other hand is established. As a result a physical interpretation of the structure function is obtained.

Precisely, firstly the problem of relating the structure function of an RC transmission line to the spatial distribution of material properties of the heat diffusion problem is transformed into the equivalent problem of relating the structure function of an LC transmission line to the spatial distribution of material properties of a wave propagation problem. Secondly the relation between the structure function and the response of such LC transmission line is outlined. Lastly a family of wave fronts and a family of rays are naturally associated to the wave propagation problem. The structure function is related to the thermal properties along rays between wave-fronts. Using this result, a strategy is proposed for naturally exploiting structure functions in practical applications.

The rest of this paper is organized as follows. In section 2 the definition of structure function is recalled. In section 3 the companion wave propagation problem is introduced. A procedure for determining the structure function from the wavepropagation problem is outlined in section 4 and applied in 5. Wave-fronts and rays of the wave propagation problem and the relation between structure function and spatial distribution of material properties is shown in section 6 and applied in section 7.

RC Transmission Line Representation of One-Port Passive Distributed Thermal Network

We refer to a generic component or package that extends in a bounded spatial region Ω. As is well known, the relation between the generated power density G(r, t) and the temperature rise distribution u(r, t) with respect to ambient temperature and the heat flux density q(r, t) is ruled by the First Principle of Thermodynamics and by Fourier's law as follows

∇ • q(r, t)) + c(r) ∂v ∂t (r, t) = G(r, t), (1)
q(r, t) = -k(r)∇u(r, t), (2)
in which c(r) is the volumetric heat capacity and k(r) is the thermal conductivity. These equations are completed by the conditions on the boundary ∂Ω of Ω and by the initial condition for the temperature rise u(r, t). The boundary conditions, assumed of Robin's type, are

q(r, t) • ν(r) = h(r)u(r, t), (3)
in which h(r) is the heat transfer coefficient and ν(r) is the unit vector outward normal to ∂Ω. The initial condition is assumed to be zero

u(r, 0) = 0. (4)
A one-port passive distributed thermal network is defined as in [START_REF] Codecasa | Canonical Forms of Thermal Impedances[END_REF][START_REF] Codecasa | Compact Modeling of Electron Devices for Electro-thermal Analysis[END_REF][START_REF] Codecasa | Compact Thermal Networks for Modeling Packages[END_REF], by introducing the power P (t) and the temperature rise T (t) measured at its port. The power P (t) determines the power density G(r, t) as G(r, t) = g(r)P (t) [START_REF] Codecasa | Physical Interpretation and Numerical Approximation of Structure Functions of Components and Packages[END_REF] in which g(r) is a function of support Σ. The temperature rise T (t) is

T (t) = Ω g(r)u(r, t) dr. (6)
The relation between the power P (t) and the temperature rise T (t) is represented, in the time domain, by a power impulse thermal response z RC (t) and, in the complex angular frequency domain, by a thermal impedance function Z RC (s). As proved in [START_REF] Codecasa | Canonical Forms of Thermal Impedances[END_REF], the response of a one-port passive distributed thermal network is the short circuit input response of a passive RC transmission line of total resistance R 0 = Z RC (0). Any definition of the port temperature T (t) different from Eq. ( 6) in general leads to an active distributed thermal network, which clearly is not equivalent to a passive RC transmission line as shown in Fig. 1. The non-PSfrag replacements T (t) singular RC transmission line is ruled by equations

P (t) C 1 R 2 K(R)
P (x, t) = -k(x) ∂T ∂x (x, t), ∂P ∂x (x, t) = -c(x) ∂T ∂t (x, t)
in which k(x) and c(x) are the thermal conductivity and volumetric heat capacity at coordinate x along the line, P (x, t) and T (x, t) are the power and temperature rise at coordinate x along the line. Different choices of the x coordinate correspond to different k(x) and c(x) and to different representations of the non-singular RC transmission line. Two representations are hereafter considered. In the travel time representation [START_REF] Sondhi | Determination of vocal-tract shape from impulse response at the lips[END_REF][START_REF] Gopinath | Inversion of the telegraph equation and the synthesis of nonuniform lines[END_REF], the coordinate x = τ is chosen is such a way that k(x) = c(x) = A(τ ) so that

P (τ, t) = -A(τ ) ∂T ∂τ (τ, t), ∂P ∂τ (τ, t) = -A(τ ) ∂T ∂t (τ, t)
In the structure function representation [START_REF] Protonotarios | Theory of Nonuniform RC Lines -Part I: Analytical Properties and Realizability Conditions in the Frequency Domain[END_REF][START_REF] Protonotarios | Analysis and Intrinsic Properties of the General Nonuniform Transmission Line[END_REF], the coordinate x = R is chosen in such a way that k(x) = 1 and c(x) = K(R) so that

P (R, t) = - ∂T ∂R (R, t), ∂P ∂R (R, t) = -K(R) ∂T ∂t (R, t).
Thus R is the cumulative resistance along the line and

C(R) = R 0 K(R) dR
is the cumulative capacitance. Such representations are equivalent. In fact from the travel time representation the structure function representation is recovered as follows

R = τ 0 dτ A(τ ) , (7)
K(R) = A 2 (τ ), (8)
C(R) = τ 0 A(τ ) dτ. (9)
Similarly from the structure function representation the travel time representation is recovered as follows

τ = R 0 K(R) dR, A(τ ) = K(R).

Companion Wave Propagation Problem

The heat diffusion problem of Eqs. ( 1)-( 4) and the one-port passive distributed thermal network of Eqs. ( 5), ( 6) can be formally associated to a companion wave propagation problem and to a companion one-port passive distributed network. Precisely the wave propagation problem in Ω is,

∂j ∂t (r, t) = -k(r)∇v(r, t) (10)
∇ • j(r, t)) + c(r) ∂v ∂t (r, t) = G(r, t), (11)
with boundary conditions, on ∂Ω,

∂j ∂t (r, t) • ν = h(r)v(r, t), (12)
and with initial conditions, in Ω,

v(r, 0) = 0, j(r, 0) = 0. ( 13
)
The companion one-port passive distributed network is defined by introducing the current I(t) and the voltage V (t) measured at its port, as follows. The current I(t) determines G(r, t) as

G(r, t) = g(r)I(t). (14)
The voltage V (t) is

V (t) = Ω g(r)v(r, t) dr. (15)
The relation between the current I(t) and the voltage V (t) is represented, in the time domain, by an impulse response z LC (t) and, in the complex angular frequency domain, by an impedance function Z LC (s). Since

Z LC (s) = s Z RC (s 2 ),
the response of this one-port distributed network is the short circuit input response of a passive LC transmission line obtained from the passive RC transmission line of Fig. 1 by substituting the resistors with inductors, as shown in Fig. 2. This passive LC transmission line is in general singular since it consists of the cascade connection of a PSfrag replacements V (t) lumped LC ladder and of a non-singular LC transmission line ruled by equations

I(t) C 1 R 2 K(L)
∂I ∂t (x, t) = -k(τ ) ∂V ∂x (x, t), ∂I ∂x (x, t) = -c(x) ∂V ∂t (x, t)
in which I(x, t) and V (x, t) are the current and voltage at coordinate x along the line. Different choices of the x coordinate correspond to different k(x) and c(x) and to different representations of the non-singular LC transmission line. Two representations are hereafter considered. In the travel time representation [START_REF] Sondhi | Determination of vocal-tract shape from impulse response at the lips[END_REF][START_REF] Gopinath | Inversion of the telegraph equation and the synthesis of nonuniform lines[END_REF] the coordinate

x = τ is chosen in such a way that k(x) = c(x) = A(τ ) so that ∂I ∂t (τ, t) = -A(τ ) ∂V ∂τ (τ, t), (16)
∂I ∂τ (τ, t) = -A(τ ) ∂V ∂t (τ, t) (17)
In the structure function representation [START_REF] Szekely | Restoration of physical structures: an approach based on the theory of RC networks[END_REF], the coordinate x = L is chosen in such a way that k(x) = 1 and c(x) = K(L) so that

∂I ∂t (L, t) = - ∂V ∂L (L, t), ∂I ∂L (L, t) = -K(L) ∂V ∂t (L, t).
Thus L is the cumulative inductance along the line and

C(L) = L 0 K(L) dL
is the cumulative capacitance. Such representations are equivalent. In fact from the travel time representation the structure function is recovered as follows

L = τ 0 dτ A(τ ) , K(L) = A 2 (τ ), C(L) = τ 0 A(τ ) dτ.  TIMA EDITIONS / THERMINIC 2005 ISBN 2-916187-01-4
Similarly from the structure function representation, the travel time representation is recovered as follows

τ = L 0 K(L) dL, A(τ ) = K(L).
The spatial distributions of material properties k(r) and c(r) are common to the heat diffusion problem and to the companion wave propagation problem. Similarly the structure function representation is common to the one-port passive distributed thermal network and to the companion one-port passive distributed network. Hereafter, in order to relate the structure function to the spatial distributions of thermal properties k(r) and c(r), the companion wave propagation problem and the companion one-port passive distributed network are considered.

Procedure for Determining Structure Functions

The lumped LC ladder in Fig. 2 can be determined by applying the continuous fraction expansion of the Z LC (s) impedance. Precisely for k = 1, 2, . . .

C 2k-1 = lim s→∞ 1/sZ 2k-2 (s) L 2k = lim s→∞ Z 2k-1 (s)/s
in which Z 0 (s) = Z LC (s) and

Z 2k-1 (s) = 1/(1/Z 2k-2 (s) -sC 2k-1 ), Z 2k (s) = Z 2k-1 (s) -sL 2k .
The continuous fraction expansion breaks down after p lumped elements are determined. It can be easily proved that p is the maximum order of the derivatives of z LC (t) with respect to time existing at t = 0. Equivalently p can be related to the spatial regularity of the k(r), c(r) and g(r) functions.

The determination of the non-singular LC transmission line in Fig. 2 can be obtained from Z p (s) as follows. It results

A(0) = 1/ lim s→∞ Z p (s).
The other values of A(τ ) can be determined by generalizing [START_REF] Sondhi | Determination of vocal-tract shape from impulse response at the lips[END_REF]. The impulse response z p (t) is computed. Then the following Fredholm integral equation is solved for I(τ, t)

I(τ, t) + A(0) 2 τ -τ z p (|t -r|)I(τ, r) dr = 1 in which -τ ≤ t ≤ τ.

It results in

τ 0 A(τ ) dτ = τ 0 I(τ, t) dt.
By recalling Eqs. ( 7)-( 9), also R, K(R) and C(R) can be recovered from A(τ ). This computation is stopped at τ = τ 0 such that R(τ 0 ) = R 0 . A procedure for determining structure functions is thus established. It has to be noted that if all derivatives of z LC (t) exist at t = 0, then the continuous fraction expansion does not break down. In this case only an infinite lumped LC network is determined.

Analytical Example: Part I

A slab of thickness L, area S, thermal conductivity k and volumetric heat capacity c is considered. The slab extends in the region Ω in which 0 ≤ x ≤ L, the power P (t) is uniformly generated within the Σ sub-region of Ω, 0 ≤ x ≤ l ≤ L. On the surface x = L of the boundary ∂Ω the temperature is set to the ambient temperature. On the rest of the boundary ∂Ω the thermal flux is set to zero. According to Eq. ( 6), the mean temperature rise in Σ is the T (t) temperature rise of a one-port passive distributed thermal network. By means of the procedure outlined in section 4, the structure function representation of the one-port passive distributed thermal network can be derived. The results for some values of l are hereafter reported.

From the continuous fraction division one lumped capacitor is derived having capacitance

C 1 = lSc.
By solving the Fredholm integral equation the travel time representation follows from which the structure function representation is obtained. For l = L it results in

γ = 1 (1 -3ρ) 1 3
, 0 < ρ < 1 3 .

For l = 2L/3 it results in

γ =              2 3 1 -9 8 ρ 1 3 , 0 ≤ ρ ≤ 37 72 4 9 5 3 -3ρ 1 3
, 37 72 ≤ ρ ≤ 5

For l = L/3 it results in

γ =                          2 9 8 27 -2 3 ρ 1 3
, 0 ≤ ρ ≤ 7 For l = L/6 it results in Both cumulative and differential structure functions are shown in Figs. 3,4 in which γ = C/LSc, ρ = R/(L/kS), κ = K/kcS 2 . As expected, as l/L approaches 0 the structure function approaches that of a uniform RC transmission line.

γ =                                                                1 18 1 27 -1 6 ρ 1 3 , 0 ≤ ρ ≤ 7 36 1 6 2 9 -1 2 ρ 1 3 , 7 36 ≤ ρ ≤ 10

Relation between Spatial Distributions of Thermal Properties and Structure Functions

Let I(t) be a unit impulse and let v(r, t), j(r, t) the solution of Eqs. ( 10)-( 13).

Let ω k (t) be the sub-region of Ω in each r of which v(r, r) = 0 at some r ≤ t and let ∂ω k (t) be its boundary. Similarly let ω k (t) be the subregion of Ω in each r of which j(r, r) = 0 at some r ≤ t and let ∂ω k (t) be its boundary. It results ω c (t) ⊃ ω k (t) and ∂ω k (t) ⊃ ∂ω c (t). From Eqs. (10)-( 13) it follows v(r, 0) = g(r) c(r) .

Thus ω c (0) is the support of g(r), that is the Σ region, while ω k (0) is the sub-region of Σ support

(∇ψ k (r)) 2 = k(r) k(r) (19)
such that ψ k (r) = 0 defines the surface ∂ω k (0). The family of rays orthogonal both to the family of wave fronts ∂ω c (t) and ∂ω k (t) is also considered.

It results in

Proposition 1 For a ray γ(t) going either from ∂ω c (0) to ∂ω c (t) or from ∂ω k (0

) to ∂ω k (t) it is t = γ(t) c(r) k(r) dγ.
From the theory of characteristic lines [START_REF] Kline | Electromagnetic Theory and Geometrical Optics[END_REF], it follows that, v(r, t) is affected by all and only the values of c(r) in ω c (t) and of k(r) in ω k (t). Moreover, as a consequence of the particular definition of V (t) given by Eq. ( 15), the following fundamental result can be proved Proposition 2 V (2t) is affected by all and only the values of c(r) in ω c (t) and of k(r) in ω k (t).

Let V (τ, t) and I(τ, t) be solutions of Eqs. ( 16), (17) when I(t) is a unit impulse. From section 4 it follows Proposition 3 V (2t) is affected by all and only the capacitances and inductances of the lumped LC ladder network and by the values of A(τ ) in the interval 0 ≤ τ ≤ t.

As a consequence of Propositions 1, 2 and 3, and of Eqs. ( 7)-( 9), the following main result is derived Proposition 4 For each 0 < R ≤ R 0 , the restriction of the structure function K(ρ) to 0 ≤ ρ ≤ R, is affected by all and only the values of c(r) in ω c (t) and of k(r) in ω k (t), being

t = γ(t) c(r) k(r) dγ = R 0 K(ρ) dρ (20)
in which γ(t) is any ray going either from ∂ω c (0) to ∂ω c (t) or from ∂ω k (0) to ∂ω k (t).

Eq. (20) expresses the physical interpretation of structure functions for generic multi-directional heat-flows. The natural manner for exploiting this result in applications is the following. Let C be an heat diffusion problem. Let N be a one-port passive distributed thermal network modeling C and let K(R) be its structure function. Let ∂ω c (t) and ∂ω k (t) be the wave fronts of the companion wave propagation problem of C, with their rays γ(t) going from 5: Propagation within Ω of the ∂ω c (t) wave front of v(r, t) and of the ∂ω k (t) wave front of j(r, t). ∂ω c (0) to ∂ω c (t) and from ∂ω k (0) to ∂ω k (t). Also let C be a second heat diffusion problem whose thermal properties present a spatially localized difference with respect to C. Let N be the one-port passive distributed thermal network modeling C and let K (R) be its structure function. If R is the smallest value at which K(R) and K (R) differ, either ∂ω c (t) is the smallest surface at which the volumetric heat capacity of C and C differ or ∂ω k (t) is the smallest surface at which the thermal conductivity of C and C differ, being

(t) (t) (t) (t) (t) (t) (t) (t) PSfrag replacements Ω Ω γ γ γ ∂Ω ∂Ω ω c (t) ω c (t) ω c (t) ω k (t) ω k (t) ω k (t) ω k (t) Σ Σ Σ Σ Σ Figure
t = γ(t) c(r) k(r) dγ = R 0 K(ρ) dρ = R 0 K (

Figure 1 :

 1 Figure 1: Equivalent passive RC transmission line.

Figure 2 :

 2 Figure 2: Equivalent passive LC transmission line.

Figure 3 :Figure 4 :

 34 Figure 3: Cumulative Structure Functions.

7Figure 6 :

 6 Figure 6: Differential Structure Functions of C and C .

  ρ) dρ. )\γ(t 1 ) is any ray going either from ∂ω c (t 1 ) to ∂ω c (t 2 ) or from ∂ω k (t 1 ) to ∂ω k (t 2 ).

	modified heat diffusion problems C , the quantity
	γ(t2)\γ(t1)	c(r) k(r)	dγ
	can be determined in which γ(t 2
			Moreover, by applying this procedure to different


L Penco

B Clement
email: brice.clement@inria.fr

V Modugno
email: modugno@diag.uniroma1.it

E Mingo Hoffman

G Nava

D Pucci

Nikos G Tsagarakis
email: nikos.tsagarakis@iit.it

J.-B Mouret

S Ivaldi

Robust Real-time Whole-Body Motion Retargeting from Human to Humanoid

Transferring the motion from a human operator to a humanoid robot is a crucial step to enable robots to learn from and replicate human movements. The ability to retarget in realtime whole-body motions that are challenging for the humanoid balance is critical to enable human to humanoid teleoperation. In this work, we design a retargeting framework that allows the robot to replicate the motion of the human operator, acquired by a wearable motion capture suit, while maintaining the whole-body balance. We introduce some dynamic filter in the retargeting to forbid dangerous motions that can make the robot fall. We validate our approach through several experiments on the iCub robot, which has a significantly different body structure and size from the one of the human operator.

I. INTRODUCTION

Current approaches for motion generation and control in humanoid robots essentially rely on optimizers and planners that search for the optimal sequence of joint commands (typically joint torques or velocities) according to an objective function that fulfills multiple tasks under several constraints [START_REF] Siciliano | Robotics: Modelling, Planning and Control[END_REF], [START_REF] Ivaldi | Special issue on whole-body control of contacts and dynamics for humanoid robots[END_REF]. For instance, it is common to have tracking tasks in the objective function to specify that the end-effectors should follow a desired trajectory [START_REF] Mansard | A versatile generalized inverted kinematics implementation for collaborative working humanoid robots: The stack of tasks[END_REF]. Unfortunately, designing and tuning the desired trajectories to realize complex tasks is time-consuming and often requires the expert knowledge of the controller/planner and of the humanoid kinematics/dynamics, which prevents an easy deployment of new tasks [START_REF] Modugno | Learning soft task priorities for safe control of humanoid robots with constrained stochastic optimization[END_REF].

An alternative is to follow an imitation approach: a human performs a movement and the robot attempts to reproduce it [START_REF] Akgun | Novel interaction strategies for learning from teleoperation[END_REF]. Kinesthetic teaching is now a mature approach for robotic arms and industrial manipulators equipped with torque sensing: it allows the human operator to show the robot the desired trajectories by physically manipulating the robot links. While this approach is relatively straightforward for robotic arms, it can be hardly done for humanoids to demonstrate whole-body movements (because it is not possible to physically guide all the robot links at once while ensuring the robot balance). Motion retargeting in this sense Fig. 1: Teleoperation of the iCub robot by an operator equipped with the XSens motion capture suit.

represents the extension of the kinesthetic teaching concept at a whole-body level.

However, real-time retargeting (or teleoperation) of entire whole-body movements is challenging for humanoids. First, direct mapping is not possible because of important differences in kinematics (e.g., joint limits, limb lengths) and dynamics (e.g., mass distribution, inertia); second, the robot needs to maintain its balance while imitating the human, so it has to trade-off between imitation and feasibility/safety. Third, since we do not know a priori which motion is going to be retargeted on the robot by the operator, it is not possible to tune offline the controller and/or the retargeting parameters for specific motions: rather, we need to provide a generic solution that is able to handle a variety of motions.

Our requirement is to have a retargeting framework that can handle a variety of whole-body motions shown by the human operator, and generates appropriate robust motions (not only feasible) for the humanoid robot. This work is motivated by our application in the project AnDy [START_REF] Ivaldi | Anticipatory models of human movements and dynamics: the roadmap of the andy project[END_REF], where we would like to show directly to the humanoid human-like movements that represents demonstrations of collaborative policies. In this paper we present our current developments to teleoperate the iCub robot (see Figure 1).

The paper is organized as follows. In Section II we discuss the previous works on retargeting and teleoperation of humanoids. In Sections III we present on overview of the framework and in Section IV the corresponding retargeting method. We discuss our teleoperation experiments in Section V and outline the future works in Section VI.

II. RELATED WORK

Human-to-humanoids motion retargeting is an active research field. The first approaches to the retargeting problem consisted in attaching virtual markers to the robot and applying the inverse kinematics to adjust the human reference motion to fit the robot constraints [START_REF] Pollard | Adapting human motion for the control of a humanoid robot[END_REF], [START_REF] Ott | Motion capture based human motion recognition and imitation by direct marker control[END_REF], [START_REF] Miura | Robot motion remix based on motion capture data towards human-like locomotion of humanoid robots[END_REF].

Several balance controllers have been then proposed for maintaining dynamic stability, since the mass distribution often differs between the human and the robot [START_REF] Dariush | Online transfer of human motion to humanoids[END_REF], [START_REF] Miura | Robot motion remix based on motion capture data towards human-like locomotion of humanoid robots[END_REF]. Dariush et al [START_REF] Dariush | Online transfer of human motion to humanoids[END_REF] considered only upper body movements and used task space control with a separate ZMP-controller to ensure balance. Yamane presented a control-based approach to imitate human motions with a force-controlled robot, generating online the joint trajectories under the assumptions that the feet maintain ground contact [START_REF] Yamane | Controlling humanoid robots with human motion data: Experimental validation[END_REF].

Recently, multi-task whole-body QP controllers were applied to retargeting. Otani and Bouyarmane [START_REF] Otani | Adaptive whole-body manipulation in human-to-humanoid multi-contact motion retargeting[END_REF] retargeted onto the HRP-4 robot in full-body dynamics simulation a bag retrieval, a door opening and a box lifting task. Their method relied on a-priori knowledge of the contact events and could not be formulated on-line without an off-line pre-processing of the recorded human motion sequence. Kanajar et al [START_REF] Kanajar | Climbing over large obstacles with a humanoid robot via multi-contact motion planning[END_REF] retargeted a challenging multi-contact motion (climbing over an obstacle) with the Coman robot: they formulated the tasks based on the observed contact states to generate keyframes for the contact points, then optimized according to the robot stability and mechanical constraints. A multirobot QP controller was also proposed to simulate humanrobot physical interaction by retargeting human motions on a simulated human model [START_REF] Otani | Generating Assistive Humanoid Motions for Co-Manipulation Tasks with a Multi-Robot Quadratic Program Controller[END_REF]. In these works, however, QP controllers were used to retarget motions offline and they were specifically tuned to the known demonstrated motions.

The first online retargeting works were for the majority related exclusively to upper body movements [START_REF] Brygo | Tele-manipulation with a humanoid robot under autonomous joint impedance regulation and vibrotactile balancing feedback[END_REF], [START_REF] Shin | On-line human motion transition and control for humanoid upper body manipulation[END_REF], [START_REF] Fritsche | First-person tele-operation of a humanoid robot[END_REF]. Frequently, the retargeting of the upper-body joints, which is important for manipulation, is done independently from the one of the legs, crucial for balancing and locomotion [START_REF] Dariush | Online transfer of human motion to humanoids[END_REF], [START_REF] Kim | Whole Body Motion Control Framework for Arbitrarily and Simultaneously Assigned Upper-Body Tasks and Walking Motion[END_REF]. At KIST for example, they combined their upper body teleoperation framework [START_REF] Shin | On-line human motion transition and control for humanoid upper body manipulation[END_REF] with a real-time walking pattern generator [START_REF] Hong | Real-time walking pattern generation method for humanoid robots by combining feedback and feedforward controller[END_REF] to retarget also the walking [START_REF] Kim | Whole Body Motion Control Framework for Arbitrarily and Simultaneously Assigned Upper-Body Tasks and Walking Motion[END_REF].

In the last few years real-time retargeting of general whole-body motions has been addressed. Koenemann [START_REF] Koenemann | Real-time imitation of human whole-body motions by humanoids[END_REF] has implemented a teleoperation of a physical robot (NAO) considering also changes in the contact support in realtime. His method is based on a simplified human model and an approximate posture mapping followed by a posture stabilization that allow transitions between double and single stances phases. However, this stabilization is not a robust solution for more complex humanoid robots that are more challenging to maintain stable during the motion, such as the iCub. Furthermore, the center of mass tracking formulation and the simplified model do not take into account the retargeting of torso movements -we will address this in our approach in Section IV.

Ishiguro et al [START_REF] Ishiguro | Bipedal oriented whole body master-slave system for dynamic secured locomotion with lip safety constraints[END_REF] proposed a secured system to teleoperate the robot LIP safety constraints. In their work they conducted some experiments involving dynamic upperbody and leg motions onto the legged robot CHIDORI and humanoid robot JAXON.

Off-line methods have been proposed as well to reconstruct the human motion within the physical constraints imposed by the retargeted subject kinematics and dynamics. Ayusawa and Yoshida [START_REF] Ayusawa | Motion retargeting for humanoid robots based on simultaneous morphing parameter identification and motion optimization[END_REF] proposed a simultaneous morphing parameter identification and motion optimization. In [START_REF] Borno | Robust physics-based motion retargeting with realistic body shapes[END_REF] Borno et al used instead a LQR-tree formulation to transfer the motion between 3D realistic human models and adapting it to the different body shapes.

Very efficient real-time methods that can reshape a dynamic motion demonstrated by a human and adapt the dynamics of the human to the dynamics of the robot, use an inverse dynamics control scheme with a quadratic programming optimization solver [START_REF] Ramos | Dancing humanoid robots: Systematic use of osid to compute dynamically consistent movements following a motion capture pattern[END_REF], [START_REF] Ramos | Dynamic motion capture and edition using a stack of tasks[END_REF]. However, even if they successfully handle dynamic overshoots due to large angular momentum variation, they require a suitable walking pattern generator to handle the linear momentum.

The model predictive control (MPC) instead, is able to deal with great linear momentum variations by considering the full temporal evolution of the system. The optimal trajectory is then updated real-time with sensor measurements. MPC enhances the expressivity capabilities of the motion generator and gives more compliance to the robot that can react to external unexpected contact events [START_REF] Scianca | Gait generation via intrinsically stable mpc for a multi-mass humanoid model[END_REF], [START_REF] Scianca | Intrinsically stable mpc for humanoid gait generation[END_REF]. The only drawback of using an MPC approach in a teleoperation framework is the underlying computational cost.

Our approach is based on an inverse kinematics (IK) control scheme with a QP optimization solver. In its simplicity, it is able to guarantee the stability of the teleoperated robot even when the retargeted movements challenge its balance. We use a ZMP retargeting correction that provides robustness and allows generality of the retargeted motion.

We evaluate our framework with several experiments with the iCub humanoid robot being teleoperated by a human operator (see also attached video). Differently from [START_REF] Fritsche | First-person tele-operation of a humanoid robot[END_REF] that teleoperated only the upperbody of the iCub, we control its entire body and we are able to retarget whole-body movements that are challenging for the balance (e.g., squat).

We formerly explored an alternative approach based on torque control that allows physical interaction with the teleoperated robot 1 . The framework was built upon the CoDyCo torque controller [START_REF] Romano | The codyco project achievements and beyond: Toward human aware whole-body controllers for physical human robot interaction[END_REF] that computes the joint torques that minimize the forces at the contacts [START_REF] Nori | icub whole-body control through force regulation on rigid non-coplanar contacts[END_REF]. The reason for choosing a desired set of contact forces is that we can use the contact wrench as a fictitious control input of the centroidal momentum dynamics. Instantaneous values for forces are computed so as to follow a desired trajectory of the center of mass and to reduce the system's angular momentum. The torque control allows a more flexible interaction with the robot, in particular it enables direct physical interaction between the human and the robot that is teleoperated. However, in the case of iCub, the absence of joint torque sensing 2 and saturation of force/torque sensors limit the torque control loop, to the point that velocity control seems more appropriate for quasi-static movements. Moreover, the structure of the torque controller ensures stability for complex motions at the price of a reduced mimicry and a notable delay 3 . In this paper we overcome these limits with a novel framework for online motion retargeting with a velocity-based QP controller that we describe in the next section.

III. FROM MOTION CAPTURE TO ROBOT MOTION A. Motion capture

The first step for a motion retargeting technique is to track the human pose. Recent developments in human motion capture allow now high-fidelity and high-frequency tracking data. Motion-capture is widely used nowadays in various fields including physiotherapy, surveillance, computer graphics and foremost in the cinema, using external cameras or wearable sensors. For our experiments we used the Xsens MVN system [START_REF]Xsens the leading innovator in 3d motion tracking technology[END_REF]. It is a wearable system consisting of 17 IMUs, providing a real-time estimation of the human posture.

Once the data is acquired from the motion capture system, it can be mapped to feasible corresponding values for the robot that are set as references for the multi-task QP controller (see Fig. 2).

B. QP controller based on stack-of-tasks

We opted for a velocity-based QP controller based on OpenSoT [START_REF] Hoffman | Robot control for dummies: Insights and examples using OpenSoT[END_REF], an open-source library implementing QP controllers based on the stack-of-tasks [START_REF] Hoffman | Whole-Body Compliant Control of iCub: first results with OpenSoT[END_REF]. Nowadays, QP controllers have become widespread thanks to their flexibility and various formulations have been proposed in the literature, [START_REF] Fava | Multi-contact motion retargeting from human to humanoid robot[END_REF], [START_REF] Otani | Adaptive whole-body manipulation in human-to-humanoid multi-contact motion retargeting[END_REF], [START_REF] Otani | Generating assistive humanoid motions for co-manipulation tasks with a multi-robot quadratic program controller[END_REF], which take into account joint velocity, acceleration and torque control as well as contact forces. The general form of QP-based controllers is:

min x || Ax -b || 2 W s.t. l ≤ Cx ≤ u (1)
where x are the controlled variables (e.g. joint velocities), A is the Jacobian task matrix, b is the desired task to perform, 2 Joint torques are estimated online with a model-based approach that exploits the sensor readings from the force/torque sensors distributed on the robot [START_REF] Fumagalli | Force feedback exploiting tactile and proximal force/torque sensing[END_REF]. 3 These limitations can be observed in the video of the teleoperation with the CoDyCo torque controller: https://youtu.be/-ib0n5shuxg W is a positive definite weight matrix, C is the constraint matrix and finally l (lower) and u (upper) are constraints. QP controllers permit to handle various type of constraints that, depending on the chosen formulation, may include robot dynamics, friction cones, self-collision avoidance, joint limits and many more (see for example [START_REF] Otani | Generating assistive humanoid motions for co-manipulation tasks with a multi-robot quadratic program controller[END_REF]).

For the teleoperation this kind of controllers is appropriate since they allow to take into account both Cartesian tasks (body segment positions) and Joint space tasks while satisfying all the robot constraints. Moreover, the tasks can be specified both in a hard and soft priority fashion.

IV. RETARGETING METHOD

The QP controller allows to define Cartesian tasks and a postural task together with their subtasks. Since the iCub body dimension is close to the one of a 5 years old child, there is a great difference with the size of the human operator wearing the xSens suit. For this reason, the direct retargeting of Cartesian tasks that are expressed in the global frame does not make any sense. We tried then to assign Cartesian tasks related to the relative positions of the feet with respect to the hips and of the hands with respect to the shoulders as done in [START_REF] Koenemann | Real-time imitation of human whole-body motions by humanoids[END_REF]. However, to do so we should accurately measure the ratio between the lengths of the limbs of the robot and of the human, and when starting the streaming of the data from Xsens, the 3D Xsens skeleton initial body orientation with respect the global frame should also match precisely the one of the robot. Furthermore, the MVN Xsens skeleton estimation is affected by some noise that makes the global frame drift on the ground over time. Hence for long teleoperations the retargeting of Cartesian tasks could be compromised because of the change of orientation of the global frame.

For these reasons, we opted for retargeting only the joint angles and simplifying the initialization process for the teleoperation together with the reliability. More specifically, we do not retarget all the joint angles but we separate the postural task in several subtasks: head sub (neck joints), torso sub (torso joints), lArm sub (left arm joints), rArm sub (right arm joints).

We additionally considered a Cartesian task for the ground position of the CoM of the robot (com sub ) and a Cartesian task for the height and orientation of the floating base (base sub ). The global position of the feet is also taken into consideration, in order to keep each foot in contact with the ground whenever it is a support link (lF oot and/or rF oot). The posture of the legs is retargeted indirectly through the CoM and floating base tasks. Our resulting selected stack of tasks for the teleoperation is the following: stack = (lF oot + rF oot + head sub )/ (com sub + base sub + torso sub + lArm sub + rArm sub );

where T a + T b means that the tasks T a and T b are in a SOFT (Weighted) Priority relation while T a /T b means the tasks T a and T b are in a HARD (Null-Space) Priority relation.

In single support mode, the task of the foot that is not the support one can be removed and a leg postural task can added for the lifted leg. However, in this paper we only show experiments in double support.

As constraints, here we mostly care for joint limits and joint velocity limits; the chosen formulation is based on joint velocity.

A. Joint Angles Mapping

The Xsens model of the human has 66 DoFs that correspond to 22 spherical joints. Except for the torso, the other spherical joints can be easily visually assigned to the corresponding ones of the iCub. Less intuitive is the mapping between the individual joints composing the spherical one. Through several tests on simulation we identified the right mapping reported in Fig. 3. For the torso, we used an approximate mapping. We considered the joints that are the most involved in the motion that are jL4L3, j1T12, and jT9T8, corresponding to the vertebrae going from the second lowest lumbar vertebra (L4) to the thoracic vertebra at the level of the breastbone (T8). The joints jT1C7 and jL5S1 (together with the other joints not reported in Fig. 3) have been considered in some tests as well but their contribution is negligible. Each joint value of the torso is obtained from the sum of the corresponding rotations of the joints jL4L3, j1T12, and jT9T8. Except for the yaw of the torso, this mapping is just approximate since the actual angle between the trunk and the hips is generally lower than the one given by the sum of the three joints. However, we noticed that the difference is not so relevant, that is reasonable if we consider that most of the inclination of the torso is generated by the hip movement.

After the identification of the human-iCub joints association, we retarget the variation of the joint angles from the starting posture

∆q i,H = q i,H -q start,H (2)
q i,R = q start,R + ∆q i,H (3)
where q start,R , q start,H and q i,R , q i,H are the human and robot joint position vectors at the start and at each time step respectively and ∆q i,H is the joint position deviation vector from the starting posture of the human in each time step.

B. Center of Mass Tracking

To track the center of mass we present an improvement of the method proposed in [START_REF] Otani | Adaptive whole-body manipulation in human-to-humanoid multi-contact motion retargeting[END_REF], [START_REF] Koenemann | Real-time imitation of human whole-body motions by humanoids[END_REF], where the robot tracks a normalized offset (denoted as o) between its support feet. Let us consider the 2D ground projection of the human CoM p CoM . The position p CoM with respect to an arbitrary support foot (in our case the left foot) is projected onto the line connecting the two feet. The result is then normalized to get a value in between 0 and 1.

o = (p CoM,H -p LF oot,H ) • (p RF oot,H -p LF oot,H ) || p RF oot,H -p LF oot,H || 2 (4)
So when the human is in a symmetric pose in double support mode the offset has a value around 0.5 and when the human stands on a single foot the offset is either 0 (left foot) or 1 (right foot). The robot 2D CoM ground projection is then

To retarget also changes of the human CoM that are not on the line connecting the two feet, we first measure the maximum backward and forward CoM displacement of the human and of the robot over their support polygon (with the origin lying on the feet line), i.e.

= (∆ CoM,H -(-δ CoM back,H )) (δ CoM f orw,H -(-δ CoM back,H )) (6)
from which we get the corresponding robot CoM displacement

∆ CoM,R = o (δ CoM f orw,R + δ CoM back,R ) -δ CoM back,R (7)
then we apply this displacement in the orthogonal direction of the line connecting the two feet of the robot.

C. Controlling the Floating Base

We can consider the pelvis as the floating base for the human. To control the height of the floating base of the robot we consider the deviation of the height of the human pelvis ∆base i,H from the starting value over time

∆base i,H = base i,H -base start,H (8)
In order to follow the human motion, we expect the corresponding deviation of the floating base of the robot ∆base i,R to be proportional to the one of the human

∆base i,H = α∆base i,R (9)
α = h base,R h base,H (10)
where α is the ratio between the height of the floating base of the robot and of the pelvis of the human when in N-pose 4 .

Then we can calculate the height of the robot base at each time step as:

base i,R = base start,R + ∆base i,R . (11)
We also retarget the change of orientation of the floating base in a similar way, by computing the roll, pitch and yaw from the quaternion information given by the motion capture system.

D. ZMP Retargeting Correction

During whole body teleoperation of humanoid robots, disastrous crashes may occur if the desired CoM trajectories recorded from the human do not ensure the balance of the controlled robot when retargeted.

To this scope, we propose a QP-based "preprocessor" that adjusts in real-time the desired commanded CoM to satisfy constraints that represent a condition for dynamic balance. In order to achieve a stable CoM trajectory we employ the linear inverted pendulum model (LIPM) in combination with the Zero Moment Point (ZMP) criterion.

The ZMP is represented with a point on the ground plane where the tipping moments, generated by the gravity and the inertial forces, are equal to zero. A humanoid robot keeps its balance if the ZMP is contained inside the support polygon of the robot. 4 N pose is a resting pose, where the human stands with lowered arms close to its body.

Through the LIPM model it is possible to establish a simple relation between the ZMP and the CoM dynamics:

pCoM = g h (p CoM -p ZM P ) (12)
where g represents the gravity acceleration, h is the height of the inverted pendulum and p CoM = (x CoM , y CoM ) T , p ZM P = (x ZM P , y ZM P ) T represent respectively the CoM and ZMP positions of the LIPM on the horizontal plane. By employing the equation ( 12) is possible to cast a QP problem to obtain an optimal correction of the desired CoM that satisfies the balance condition on the humanoid CoM + T g h t-1 (p CoM -p ZM P ) lb SP < p ZM P < ub SP [START_REF] Kanajar | Climbing over large obstacles with a humanoid robot via multi-contact motion planning[END_REF] where ṗdes

CoM is the desired CoM velocity, T is the sampling time, ṗt-1

CoM , h t-1 are respectively the last CoM velocity and the last CoM height measured from the robot and lb SP and ub SP are the lower and upper bound of the support polygon of the robot. Updating the height of the LIPM at each time step with the actual robot CoM z position provides a better model to estimate the ZMP position and consequently, a more accurate correction.

V. EXPERIMENTS

We set up two experiments to validate our approach. First, through dynamics simulations on the simulated robot in Gazebo, we show how the ZMP retargeting correction is essential to retarget a motion that is dynamically stable onto the robot. Then, we tested our framework by teleoperating the real robot iCub in real-time. Videos of our experiments with the real robot can be seen at https:// youtu.be/CjLQu_6ifAE and https://youtu.be/ iZVAacyvYhM.

Simulated robot -We selected three kind of motions to show the efficacy of our method: squat motion, hip roll exaggerate motion and a grasping motion involving some torso and leg movements (see Figure 7). Figure 5 shows how the ZMP associated to the retargeted motion without the correction lies outside the support polygon, making the robot fall inevitably (that is why we show this in simulation). This is due both to mechanical limitations of the robot, which cannot achieve the same CoM displacement given the retargeted joint values (that might go beyond the robot joint limits), and to the generation of some momentum different from the human one that makes the desired robot CoM trajectory unstable. Hence, the ZMP trajectory is corrected in real-time to stay inside the support polygon and the CoM trajectory is modified accordingly (see Figure 6).

Real robot -We teleoperated the robot trying to move all its links to show the effectiveness of our framework. During the teleoperation the ZMP position of the robot always lies inside the support polygon as expected (see Figure 8). The robot joint trajectories follow the retargeted values guaranteeing the mimicry of the human motion (see Figure 9). Even though the legs joint angles are not taken into consideration in the stack, the retargeting of the height and of the orientation of the floating base makes the resulting legs joint trajectories close to the human ones. This framework allows to teleoperate the robot for a potentially unlimited amount of time. Even if the human performs some motions that are very challenging for the balance, the ZMP constraints introduced by our framework make the robot maintain its balance. The video attachment at https://youtu.be/ iZVAacyvYhM shows that even if the operator is falling, the robot does not and keeps its balance.

VI. CONCLUSIONS

We proposed a framework to teleoperate the iCub robot, based on an inverse kinematics (IK) control scheme with a QP optimization solver. It allows a robust real-time retargeting of generic whole-body motions. The generality is obtained thanks to the mapping of all the main human joints: not only the arms and the legs, but also the torso and the head. The CoM as well is retargeted in a way that does not restrict the range of possible retargeted robot motions. The robustness is given by a ZMP correction approach that guarantees the stability of the retargeted trajectory in double support. In the future we will overcome the double support limitation of this paper, integrate single support and teleoperation with locomotion.

Fig. 2 :

 2 Fig. 2: Teleoperation pipeline.

Fig. 3 :

 3 Fig. 3: Mapping between the Xsens joints and the iCub joints.

Fig. 4 :

 4 Fig. 4: Determination of the normalized offset from the ground projection of the CoM p CoM,H and the feet positions of the human.

  -ṗCoM ) T R ( ṗdes CoM -ṗCoM ) s.t. ṗCoM = ṗt-1

Fig. 5 :

 5 Fig. 5: ZMP position of the simulated robot without (red) and with (blue) the stability correction during the teleoperation while performing a squat motion (up-left), a hip roll exaggerate motion (up-right), a grasping motion involving some torso and leg movements (down).

Fig. 6 :

 6 Fig. 6: CoM desired (red) and stabilized (blue) position of the simulated robot during the teleoperation while performing a squat motion (up-left), a hip roll exaggerate motion (upright), a grasping motion involving some torso and leg movements (down).

Fig. 7 :

 7 Fig. 7: Snapshots of the teleoperated simulated robot and of the Xsens skeleton while performing a squat motion, a hip roll exaggerate motion and a grasping motion during the teleoperation.

Fig. 8 :

 8 Fig. 8: iCub: CoM (left) and ZMP (right) position trajectories during the teleoperation experiment.

Fig. 9 :

 9 Fig. 9: iCub: joint trajectories (blue) compared to the retargeted values (red) during the teleoperation experiment.

Fig. 10 :

 10 Fig. 10: Snapshots from the teleoperation experiment. Remarkably, the robot does not fall when the human operator falls back (frame on the bottom-right).

  δ CoM back,H , δ CoM f orw,H and δ CoM back,R , δ CoM f orw,R respectively. Then, we retarget the human CoM displacement ∆ CoM,H (comprised within the minimum value -δ CoM back,H and maximum value δ CoM f orw,H ) onto the robot, by computing the offset o

Some experiments with the torque controller are shown in this video: https://youtu.be/-ib0n5shuxg

ACKNOWLEDGMENTS

The authors would like to thank Augusto Cesar Francisco for his precious help when wearing the motion capture suit.

*This work received funding European Commission through the project H2020 AnDy (GA no. 731540) and CogIMon (GA no 644727), and from the European Research Council (ERC) (GA no. 637972, project "ResiBots").


Marie Charbonneau

Valerio Modugno

Francesco Nori
email: fnori@google.com

Giuseppe Oriolo

Daniele Pucci

Serena Ivaldi

Learning robust task priorities of QP-based whole-body torque-controllers

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

I. INTRODUCTION

Applications involving highly redundant robots, such as humanoids, have the potential to disrupt many industries and bring countless benefits to society. Nevertheless, the design of efficient and safe controllers for humanoid platforms is highly challenging, especially for applications involving interaction with the environment. Indeed, a great number of issues can stem from the complexity of the control algorithms, task inconsistencies, or the need of guaranteeing constraints satisfaction, in order to generate physically consistent behaviors. A highly effective solution to the whole-body control problem is to decompose a complex behavior into several elementary tasks, according to a prioritization scheme [START_REF] Ivaldi | Special issue on whole-body control of contacts and dynamics for humanoid robots[END_REF], typically either a "strict" or a "soft" task hierarchy.

With strict prioritization strategies (such as a traditional stack-of-tasks), a fixed task hierarchy is assured by geometrical conditions such as a null space task projector [START_REF] Saab | Dynamic whole-body motion generation under rigid contacts and other unilateral constraints[END_REF], [START_REF] Nava | Stability analysis and design of momentum-based controllers for humanoid robots[END_REF], or by the use of quadratic programming (QP) to compute optimal control actions [START_REF] Dafarra | A Control Architecture with Online Predictive Planning for Position and Torque Controlled Walking of Humanoid Robots[END_REF]. Conversely, prioritization in a soft task hierarchy (based on a weighted combination of tasks), can be achieved by assigning each task a weight defining its relative importance [START_REF] Salini | Synthesis of complex humanoid whole-body behavior: A focus on sequencing and tasks transitions[END_REF], [START_REF] Otani | Generating assistive humanoid motions for co-manipulation tasks with a multi-robot quadratic program controller[END_REF]. Fig. 1: Different iCub models performing the same wholebody motion with several tasks. In this paper, with the purpose to eventually ease the passage from simulated to realworld robots, we optimize the task priorities for robustness, in order to allow their transfer from one robot model to a different one, without the need of re-tuning.

Usually, task priorities are designed a priori by experts, then manually tuned to adjust the task ordering, timing, transitions, etc. It is a tedious operation which, depending on the complexity of the system and the task, can take a notable portion of a researcher's time, and can be particularly hard in the case of whole-body control of floating-base platforms.

A recent line of research seeks to tackle this issue, aiming to automatically learn whole-body task priorities [START_REF] Dehio | Multiple task optimization with a mixture of controllers for motion generation[END_REF], [START_REF] Ha | Evolutionary optimization for parameterized wholebody dynamic motor skills[END_REF], [START_REF] Modugno | Learning soft task priorities for safe control of humanoid robots with constrained stochastic optimization[END_REF]. For instance in [START_REF] Dehio | Multiple task optimization with a mixture of controllers for motion generation[END_REF], task coefficients are learned for controlling a single arm, allowing the transfer of acquired knowledge to different tasks. In [START_REF] Ha | Evolutionary optimization for parameterized wholebody dynamic motor skills[END_REF], task generalization is achieved with an evolutionary algorithm employed to learn policies for a range of tasks using parametrized motor skills. Then, in [START_REF] Modugno | Learning soft task priorities for safe control of humanoid robots with constrained stochastic optimization[END_REF], task priorities guaranteeing the non violation of system constraints are learned for bimanual motions.

Learning methods often use a random exploration over numerous experiments, which could be problematic to perform on hardware. For this reason, training is preferably performed in simulation. However, inherent differences between simulated and real robots can render an optimal solution untransferrable from one to the other. This gap between reality and simulation needs to be accounted for, in order to achieve automatic parameter tuning.

Solutions to this problem have recently been addressed by trial-and-error algorithms [START_REF] Cully | Robots that can adapt like animals[END_REF], [START_REF] Spitz | Trial-and-error learning of repulsors for humanoid qp-based whole-body control[END_REF]. In [START_REF] Cully | Robots that can adapt like animals[END_REF], prior knowledge from simulations was exploited to find acceptable behaviors on the real robot, in few trials. In [START_REF] Spitz | Trial-and-error learning of repulsors for humanoid qp-based whole-body control[END_REF], a trial-and-error learning algorithm encouraged exploration of the task space, allowing adaptation to inaccurate models, also in few trials.

Learning from simulation can also simplify data gathering, but there is no guarantee that the acquired knowledge is effective in the real world. For example, optimal control and movement primitives are combined in [START_REF] Clever | Cocomopl: A novel approach for humanoid walking generation combining optimal control, movement primitives and learning and its transfer to the real robot HRP-2[END_REF] to find solutions which can be easily deployed on the real robot, but they strongly rely on the accuracy of the simulated model. Also, since QP solvers often allow for constraints relaxation, strict constraints satisfaction is not always ensured by the frameworks presented above. In [START_REF] Modugno | Learning soft task priorities for safe control of humanoid robots with constrained stochastic optimization[END_REF], the learning procedure guarantees strict constraints fulfillment, but transferring knowledge from simulation to reality did not fully achieve the desired behavior. This implies that constraints satisfaction is beneficial, but not enough to achieve transferability: solutions which are robust rather than optimal are needed, in order to achieve a better generalization.

Looking for robust solutions is central to transfer learning, which seeks to exploit previously acquired knowledge [START_REF] Pan | A Survey on Transfer Learning[END_REF] to allow a model trained on one task to be re-purposed on a second related one. This approach can benefit from domain randomization (DR) [START_REF] Tobin | Domain randomization for transferring deep neural networks from simulation to the real world[END_REF], which consists in randomizing some aspects of the simulation to enrich the range of possible environments experienced by the learner. In [START_REF] Antonova | Reinforcement learning for pivoting task[END_REF], robust policies for pivoting a tool held in the robot's gripper were learned in simulation, given random friction and control delays, such that the learned policies proved to be effective on the real robot as well.

This work proposes to apply the idea of DR to whole-body controllers. In this context, we want to ensure that balance is maintained while performing a task, even if large differences exist between the learning domain and the testing domain.

To achieve this result, we combine a DR approach with fitness functions promoting the robustness of the learned controller, and provide a method to learn robust task priorities which achieve desired goals, while allowing to facilitate the transfer of results from simulation to reality. The effectiveness of the proposed method is demonstrated by optimizing parameters in simulation for the goal of making the iCub perform steps, and showing that it is possible to overcome issues related to the transferability problem.

The paper is organized as follows. Section II introduces the modeling of a floating-base robot for whole-body control. Section III describes the task-based whole-body controller used for validating the approach, as well as the constrained optimization algorithm used for learning task priorities. Experiments led with different models of the iCub robot are then illustrated in Section IV. Results are discussed in V, leading to the conclusion of the paper.

II. BACKGROUND

A. Notation used throughout the paper

• The set of real numbers is denoted by R • The transpose operator is denoted by (•)

• The Moore-Penrose pseudoinverse is denoted by (•) † • |v| denotes the Euclidean norm of a vector v ∈ R n • Given a time function f (t) ∈ R n , its first-and secondorder time derivatives are denoted by ḟ (t) and f (t)

B. System Modelling

The robot is modelled as described in [START_REF] Nava | Stability analysis and design of momentum-based controllers for humanoid robots[END_REF], where it is assumed to be a free-floating system composed of n+1 links connected with n 1-DOF (degree of freedom) joints. The robot configuration space is characterized by q, composed of the position and orientation of a base frame B attached to a link of the robot, and the joint angles s. The robot velocity is expressed as ν = (v B , ṡ), where v B is the velocity of the base frame expressed with respect to the inertial frame.

Applying the Euler-Poincaré formalism [16, Chapter 13.5] yields the following system equation of motion:

M (q) ν + h(q, ν) = ζτ + nc k=1 J c k f k (1)
where M ∈ R n+6×n+6 is the mass mass matrix, h ∈ R n+6 is the bias vector of Coriolis and gravity terms, τ ∈ R n is a vector representing the actuation joint torques, and ζ = (0 n×6 , 1 n ) is a selector matrix, in which 1 n ∈ R n×n denotes the identity matrix and 0 n×6 ∈ R n×6 is a zero matrix. Assuming that the robot is interacting with the environment through n c distinct wrenches (as an abuse of notation, here we do not define wrench as the dual of a twist), f k ∈ R 6 denotes the k-th external wrench applied by the environment on the frame c k of the robot.

III. METHODS

The method proposed for learning robust task priorities (outlined in Fig. 2) relies on two main parts: (i) a QPbased whole-body torque-controller which tracks desired task trajectories and sends joint torque commands to the robot, and (ii) a black-box constrained stochastic optimization procedure, posing no restriction on the structure of the learning problem. It is used to optimize task priorities w as follows: at the end of a roll-out (i.e. execution of a footstep), the fitness of the obtained robot behavior is evaluated, allowing the optimization algorithm to update the task weights.

A. QP-based whole-body torque controller

The control framework used in this paper was defined from a stack of tasks solved by quadratic programming. The following subsections define the control input and objective, before describing an optimization problem allowing the computation of a control input which achieves the objective.

1) Control input:

A control input u = τ F c is composed of joint torques τ and contact forces F c , where F c ∈ R 6nc is a stacked vector of contact wrenches. Defining J c as the contact Jacobian, stacked in the same way as F c , and B = ζ J c , the system dynamics (1) can then be reformulated as

M (q) ν + h = Bu (2)
2) Control objective: The developed controller has the following objectives, from which a stack of tasks is defined:

-Stabilize the center of mass position

X CoM ∈ R 3 -Stabilize the stance foot pose X stance ∈ R 6 -Stabilize the swing foot pose X swing ∈ R 6 -Stabilize the neck orientation X neck ∈ R 3 -Track joint positions (postural task) s -Minimize joint torques τ
The velocity ẊT of a Cartesian task T can be computed from the associated Jacobian J T and the robot velocity ν with ẊT = J T ν, where T ∈ {CoM, stance, swing, neck}. Deriving this expression yields the task acceleration:

ẌT = JT ν + J T ν (3)
In view of Eq. ( 2) and ( 3), task accelerations ẌT can be formulated as a function of the control input u:

ẌT (u) = JT ν + J T M -1 (Bu -h) (4)
Stabilization of a Cartesian task may then be attempted by minimizing XT (u) = ẌT (u) -Ẍ * T , the error on the task acceleration, where Ẍ * T is a feedback term composed of reference linear accelerations obtained with a proportionalderivative (PD) feedback control policy, and reference angular accelerations computed as in [17, section 5.11.6, p.173].

As for the postural task, similarly as above, in view of Eq. ( 2), one can obtain s(u), a formulation of the joint accelerations s as a function of the control input u:

s(u) = ζ M -1 (Bu -h) (5)
Stabilization of the postural task may then be attempted by minimizing s(u) = s(u) -sref , the error on the joint accelerations, with sref obtained through inverse kinematics by solving the following optimization problem:

sref = arg min s 1 2 T w T νT 2 + w s νs 2 (6)
In the above, the error on robot acceleration is computed

with νT = ν -J † T JT ν -Ẍ * T
for Cartesian tasks, and

with νs = ν -v * B s *
for the postural task. s * is a feedback term computed using a PD control policy, and v * B is a feedback term on the base velocity used to ensure convergence of the base acceleration to zero.

3) QP formulation: In the context where tasks may have relative priorities, soft task priorities would allow to achieve a trade-off between weighted tasks. A set of task weights is defined, in order to attribute each task a priority, with w = {w CoM , w stance , w swing , w neck , w s , w τ }

The terms w CoM , w stance , w swing , w neck ∈ R refer to weights associated to the CoM, stance foot, swing foot and neck Cartesian tasks, and w s , w τ ∈ R to the weights associated to the postural task and joint torque regularization.

The control architecture described above can thus be formulated as the following optimization problem.

u * = arg min u 1 2 cost (8a) subject to Cu ≤ b (8b)
where the constraint Eq. (8b) ensures that the contact forces remain within the associated friction cones. The cost function (8a) is computed as the weighted sum of all task objectives:

cost = T w T XT (u) 2 +w s s(u) 2 +w τ |τ (u)| 2 ,

in which

w T ∈ R refers to weights associated to Cartesian tasks.

Reorganizing the terms in the cost function, one can easily verify that it has the form of a QP problem.

B. Constrained stochastic optimization

In this work, the learning problem is cast as a black-box constrained stochastic optimization. Given a fitness function φ(w) : R n P → R (with n P the number of task weights to optimize), sets of n IC inequality constraints and n EC equality constraints h, g, the idea is to find an optimal solution w • ∈ W ⊆ R n P to the problem:

w • = arg max w φ(w) (9a)
subject to

g i (w) ≤ 0 i = 1, . . . , n IC (9b) h i (w) = 0 i = 1, . . . , n EC (9c)
The Covariance Matrix Adaptation -Evolution Strategy (CMA-ES), first described in [START_REF] Hansen | Completely derandomized selfadaptation in evolution strategies[END_REF], is a standard tool for solving black-box optimization problems. In its basic implementation, CMA-ES defines a multivariate Gaussian distribution N ( w, σ 2 , Σ) with mean w, step size σ 2 and covariance Σ. However, since CMA-ES was not originally designed for handling constraints, this work employs a constrained extension of it called (1+1)CMA-ES, which uses covariance constrained adaptation (CCA) [START_REF] Arnold | A (1+1)-CMA-ES for constrained optimisation[END_REF]. This variant was benchmarked in [START_REF] Modugno | Learning soft task priorities for safe control of humanoid robots with constrained stochastic optimization[END_REF] for bimanual tasks, as well as in [START_REF] Modugno | Safe trajectory optimization for whole-body motion of humanoids[END_REF] for trajectory optimization of a whole-body movement (standing up from a chair).

A single iteration, or generation, of (1+1)-CMA-ES comprises several stages. First, a sample w 1 is drawn according to the rule w 1 = w + σDz, where D is the Cholesky factor of the covariance matrix Σ = D D and z is a standard normal distributed vector z ∼ N (0, I). The fitness function φ is evaluated for this sample, and the Gaussian distribution is updated based on the result, as described in [START_REF] Igel | A computational efficient covariance matrix update and a (1+1)-CMA for evolution strategies[END_REF]. This update also exploits the previous generations to change the covariance and step size of the search distribution. Finally, information about the successful iterations is stored in a search path s ∈ R n P . A more detailed description of the algorithm can be found in [START_REF] Modugno | Safe trajectory optimization for whole-body motion of humanoids[END_REF].

IV. EXPERIMENTS

This section exposes experiments with the iCub robot [START_REF] Metta | The icub humanoid robot: An open platform for research in embodied cognition[END_REF]. They were performed in order to validate empirically that the method described in section III is capable of generating task priorities which (i) produce robust whole-body motions, even when contacts due to physical interaction with the environment evolve in time, and (ii) can cope with imperfections in the robot model, disturbances and noise.

1) Robot setup: Experiments were conducted in simulation using the open-source robot simulator Gazebo [START_REF] Koenig | Design and use paradigms for gazebo, an open-source multi-robot simulator[END_REF]. They were performed with the iCub robot, using 23 DOF on legs, arms and torso. The robot is equipped with 6 force/torque sensors on the arms, legs and feet, allowing whole-body dynamics computations. We tested our method using two distinct models of the iCub with different inertial properties: the first one with a tethered power supply and the second one with a battery pack on the back, as shown in Fig. 1.

2) Whole-body motion: The controller described in section III-A was developed in Matlab/Simulink using WBToolbox [START_REF] Romano | A whole-body software abstraction layer for control design of free-floating mechanical systems[END_REF]. It can be used either for a simulated or a real robot.

It is applied to the task of performing steps, i.e. performing the following sequential movements for each step: move the CoM above the stance foot, move the swing foot up by 0.025m, move the swing foot back to its initial pose, and move the CoM back to its initial position. The passage from one movement to the next takes place when Cartesian task errors are smaller than a threshold and contact forces are consistent with the desired motion, or a maximum state duration has been reached. Over the course of a step, the desired stance foot pose, neck orientation and posture remain at their initial value.

Note that over all experiments, the proportional-derivative gains used to compute reference accelerations associated to Cartesian and postural tasks are kept constant.

The next subsections describe the experiments we performed, while Table I presents them in summary form. Training and testing were performed with different robot models and stepping tasks, as well as under different randomized conditions (RC) as described in Table II.

A. Training with the tethered iCub model

In order to optimize task weights, three different fitness functions were compared. The first one, φ p , favored performance on the Cartesian tasks and less deployed effort. The second one, φ r , focused on robustness, by favoring solutions with smaller excursion of the ZMP position P ZM P with respect to the center of the support polygon O SP . The third fitness function, φ pr , was a combination of the first two:

φ p = - 1 X Tmax t end t=0 T XT 2 - 0.0001 τ max t end t=0 |τ | 2 (10a) φ r = - 1 P ZM Pmax t end t=0 |P ZM P -O SP | 2 (10b) φ pr = 1 2 (φ p + φ r ) (10c)
X Tmax , τ max and P ZM Pmax are normalization factors. Inequality constraints to the learning problem were defined on joint position limits and torque limits: they act as an extra safety built on top of the QP controller.

The task used for the optimization is to perform 1 step, as described in section IV-.2. The simulation was limited to 10 seconds, allowing the robot to perform one step and shift its weight to the side in order to start a second one, making sure that the robot remained stable after foot touchdown.

Early termination of a simulation took place in cases where the robot fell or the QP in Eq. ( 8) could not be successfully solved (which may happen when the weights being tested are far from being optimal, and lead to a configuration in which the QP solver simply can not find a solution). In these cases, a penalty of -1.5 was added to the fitness in Eq. [START_REF] Cully | Robots that can adapt like animals[END_REF].

Optimized task priorities were obtained by performing 200 iterations of (1+1)-CMA-ES [START_REF] Arnold | A (1+1)-CMA-ES for constrained optimisation[END_REF] applied to the control framework, with an exploration rate of 0.1; each learning experiment was repeated 10 times. Since task priorities as used in Eq. ( 6) and ( 8) are relative to each other, w CoM was attributed a fixed value of 1. The remaining task priorities were attributed bounds as shown in the bottom of Table I. Furthermore, the weight values w 0 = [1, 1, 1, 0.1, 1e-3, 1e-4] obtained by hand were verified to allow the tethered iCub model to successfully perform the desired stepping motion, and were therefore used as a starting point.

Moreover, to achieve robustness through domain randomization and enable the controller to eventually cope with real-world data, the robot was subjected to the randomized conditions 1 to 5 (see Table II) at each learning iteration. These conditions constitute the set of conditions j under which the controller had to perform. In particular, tests performed in simulation showed that applying a force of 10N on the chest for no more than 1 second was sufficient to destabilize the robot when using unoptimal weights. Thus, using such disturbances during learning can encourage the generation of robust task priorities.

An additional set of learning experiments was performed for φ pr without using DR, to assess the contribution of DR.

Results showed that 200 iterations were sufficient to achieve strong convergence. Weights obtained with each of the fitness functions in Eq. ( 10) are shown in Table I.

B. Testing 1: with the tethered iCub model

In order to validate the robustness achieved with the optimized weights, each set of them were tested on the same iCub model used for training, but not submitted to randomized conditions. The testing task in this case was to achieve 6 consecutive steps. The success rates achieved with the optimized weights from each fitness function are shown in Table I, and typical trajectories achieved using each fitness function are shown in the left part of Fig. 3 and in Fig. 4.

C. Testing 2: with the backpacked iCub model

In order to replicate conditions similar to performing experiments on the real robot, we prepared a second model of the iCub, with a battery pack installed on its back. This II. The column # gives the number of training experiments performed given each fitness function, while only 1 set of manually tuned gains was used. In the column DR, "Yes" means that the RC defined at the top of the table were used for DR while training; "No" means that they were disabled. The "success" columns report success rates achieved when testing the sets of trained weights.

| n F ∈ Z, n F ≤ 7}
and for each wrench i: II), with the testing task to achieve 6 consecutive steps. The success rates achieved with the optimized weights from each fitness function are shown in Table I. Successful CoM and feet trajectories obtained with φ pr , using DR and not, are shown in Fig. 3. As could be expected, due to the addition of noise, ZMP trajectories showed to be highly noisy and are therefore not shown. These results show that weights obtained with φ pr allow for a higher robustness of the controller. Additionally, the rate of success achieved with DR shows to be significantly higher than without DR, demonstrating that DR did have a measurable impact on the achieved robustness.

at time t F i (s, with 1e-2 precision) t F i | t F i ∈ R + , t F i ≤ 10 duration d F (s, with 1e-2 precision) d F i | d F i ∈ R + , d F i ≤ 1 direction (γ F i , θ F i , ϕ F i ) (rad) γ F i | γ F i ∈ R, γ F i ≤ 2π θ F i | θ F i ∈ R, θ F i ≤ 2π ϕ F i | ϕ F i ∈ R, ϕ F i ≤ 2π force magnitude F F i (N) F F i | F F i ∈ R, F F i ≤ 10 torque magnitude τ F i (Nm) τ F i | τ F i ∈ R, τ F i ≤ 10
V. DISCUSSION AND CONCLUSION In summary, the proposed method generated task priorities for successful whole-body control of different robot models, in 200 learning iterations. It was demonstrated by performing training on a first model of the iCub robot, then testing on a different model submitted to diverse working conditions.

Results achieved with φ p , a fitness function favoring task performance, were likely limited by overfitting on the model used for learning and did not allow much robustness with respect to disturbances. On the other hand, optimizing for robustness with φ r allowed higher chances of success when changing conditions, by encouraging smaller ground reaction forces and the generation of angular momentum through the torso and arms. However, φ r might have neglected the fulfillment of tasks, which are used for keeping balance. Instead, using φ pr , a fitness function combining robustness and performance, allowed to achieve superior results.

With φ pr , the swing foot placement, crucial for stability at touchdown, was given high importance, while the neck orientation task a lesser one, allowing compliance to external forces, which thus facilitates recovery from external perturbations and contact switching. As for the postural task, its allotted low priority allows it to be used as regularization (just as joint torques), instead of competing with Cartesian tasks. Such a solution is interesting, as it may not have been a priori self-evident to an expert defining task priorities. Furthermore, the ranges over which sets of optimized weights were obtained show that the problem has multiple local minima. Therefore, although task priorities require proper tuning, the controller is not highly sensitive to a single precise adjustment of task weights.

Eventually, in order to further improve performance when using a different robot model or moving on to experiments on the real robot, one may be interested in performing another round of optimization on the task priorities, for example exploiting data-efficient learning methods, as done in [START_REF] Spitz | Trial-and-error learning of repulsors for humanoid qp-based whole-body control[END_REF].

Nonetheless, the robustness achieved with the proposed method has the potential to allow higher success when passing from simulation to real-world experiments. Further works shall explore different DR conditions, such as randomizing values associated to the dynamics of the system, control delay, or measurement noise amplitude and delay. The extension of the proposed method to additional parameters such as feedback gains will also be examined.

Fig. 2 :

 2 Fig.2: Overview of the proposed method. Given task priorities w, the QP-based controller computes a control input u under a set of conditions j (e.g. desired trajectories, disturbances, noise). An outer learning loop allows the optimization of the task weights.

Fig. 3 :

 3 Fig.3: Typical CoM and feet trajectories obtained for 6 strides performed with the tethered and backpacked iCub models, given initial weights w 0 or weights optimized using DR (except for "no DR") with φ p , φ r , φ pr . Each color denotes the use of a different set of optimized weights. The x, y and z axes correspond respectively to the sagittal, frontal and vertical axes.

6 .

 6 Gaussian noise on joint velocity signals On 7. Gaussian noise on F/T sensor signals On model was subjected randomized conditions 5 to 7 (see Table

Fig. 4 :

 4 Fig. 4: Typical ZMP trajectories obtained for 6 steps performed with the tethered iCub model, from left to right: given initial weights, weights optimized with φ p , φ r , φ pr . Each color (blue, red or yellow) denotes the use of a different set of optimized weights. The x, y and z axes correspond respectively to the sagittal, frontal and vertical axes.

TABLE I :

 I Summary of performed experiments and achieved results Randomized conditions RC are defined in Table

	Scenario					Training				Testing 1	Testing 2
	Robot					tethered				tethered	backpacked
	Task					1 step				6 steps	6 steps
	RC					1, 2, 3, 4, 5				-	5, 6, 7
	Fitness	#	DR w CoM	wstance	w swing	w neck	ws	wτ	success	success
	hand tuning	1	No	1	1	1	1	1e-3	1e-4	1 /1	0 /1
	φp	10 Yes	1	0.2 ± 0.3 1.1 ± 1.2 (1 ± 3)e-3	0.5 ± 0.3	(2 ± 5)e-5	5 /10	0 /10
	φr	10 Yes	1	1.6 ± 1.2 1.8 ± 1.0	0.1 ± 0.1	(4 ± 7)e-3	1e-10	7 /10	5 /10
	φpr	10 Yes	1	0.9 ± 1.3 2.4 ± 1.1	0.6 ± 1.2	1e-6	1e-10	10 /10	10 /10
	φpr	10	No	1	1.0 ± 0.3 0.1 ± 0.3	0.2 ± 0.1	1e-6	(4 ± 5)e-6	8 /10	2 /10
	weights lower bounds:	1	1e-4	1e-4	1e-10	1e-6	1e-10
	weights upper bounds:	1	10	10	10	10	0.1

TABLE II :

 II Randomized set of conditions j

Randomized Condition (RC) in j

Random value 1. Gaussian noise on F/T sensor signals On / Off 2. Appointed swing foot Left / Right 3. Direction in which swing foot is moved Front / Back 4. X CoM d moved forward by δ (m) δ | δ ∈ R + , δ ≤ 0.02 5. n F external wrenches applied on chest {n F

This work was supported by the EU H2020 program under the Marie Sklodowska-Curie SECURE grant (n.642667), as well as the european projects An.Dy (n.731540) and Comanoid (n.645097).