Giovanni De Magistris
email: giovannidemagistris@hotmail.it

Sylvain Miossec

Adrien Escande

Abderrahmane Kheddar

Design of Optimized Soft Soles for Humanoid Robots

Keywords: Shape Optimizations, Compliant soles, Humanoid robots, Walking

We describe a methodology to design foot soles for a humanoid robot given walking gait parameters (i.e. given centerof-mass and zero-moment-point trajectories). In order to obtain an optimized compliant sole, we devised a shape optimization framework which takes -among other inputs, an initial rough (simplified) shape of the sole and refines it through successive optimization steps under additional constraints and a cost function. The shape is optimized based on the simulation of the sole deformation during an entire walking step, taking time dependent input of the walking pattern generator into account. Our shape optimization framework is able to minimize the impact of the foot with the ground during the heel-strike phase and to limit foot rotation in case of perturbations. Indeed, low foot rotation enforces a vertical posture and secures the balance of the humanoid robot. Moreover, weight restriction (formulated as a constraint on the sole volume) is added to our optimization problem.

Introduction

Bipedal walking on flat ground is a well-documented and thoroughly studied problem in the humanoid research field. An excellent review is presented in [START_REF] Wieber | Modeling and control of legged robots[END_REF]. Recent studies show that humanoid walking on flat floors reached a good level of maturity. Indeed, humanoid robots can walk with human-in-loop [START_REF] Bussy | Human-humanoid haptic joint object transportation case study[END_REF], with human stylized gaits [START_REF] Miura | Robot motion remix based on motion capture data -Towards human-like locomotion of humanoid robots[END_REF][START_REF] Mombaur | From human to humanoid locomotion-an inverse optimal control approach[END_REF] or with emotional behavior gaits [START_REF] Felis | Modeling and identification of emotional aspects of locomotion[END_REF]. In the latter works, the precise control of the contact interaction between the humanoid robot feet and the ground plays an important role. Reducing the impact at the landing is important to perform smooth walking gaits. For example, reducing impacts and having smooth walking is critical in human-inthe-loop walking under hard visual servoing constraints [START_REF] Agravante | Collaborative human-humanoid carrying using vision and haptic sensing[END_REF] or for precise stylized or emotional walking gaits.

To reduce the impact of the foot landing, different studies illustrate the importance of the foot shape design. Adamczyk et al. [START_REF] Adamczyk | The advantages of a rolling foot in human walking[END_REF] show that in a human passive walking locomotion, the energy loss during heel contact with the ground varies with (i) the foot shape, and (ii) the first contact point between the foot and the ground [START_REF] Kwan | Optimal foot shape for a passive dynamic biped[END_REF]. They explain how a rolling foot shape considerably reduces impact-shock when the heel touches the ground. Recently, Fukushima et al. [START_REF] Fukushima | Gait Stabilization of Passive Dynamic Walking by Foot Shape Optimization[END_REF] show that foot shape, once optimized, allows passive walking on various slopes. Miyakoshi and Cheng [START_REF] Miyakoshi | Examining Human Walking Characteristics with a Telescopic Compass-like Biped Walker Model[END_REF] propose a model for biped walking with a curved surface foot (using an array of point contacts). This model improves the stability of walking and reduces impact forces during foot contact maintaining the same human-like walking characteristics. Aerts et al. [START_REF] Aerts | The mechanical properties of the human heel pad: a paradox resolved[END_REF] found by experimental measurements that the shock absorption in humans happens at the foot pads. In most people, the heel pad is 25mm thick, and compresses by 12mm during stance. To make the foot more controllable, the material is hyperelastic -it starts out being compliant, and then stiffen during stance. Furthermore, Sardain and Bessonnet [START_REF] Sardain | Gait Analysis of a Human Walker wearing Robot Feet as Shoes[END_REF] show how the sole mechanism influences the stabilization of biped robot walking.

The large majority of humanoid robots are embedded with force sensing devices at the ankles. In order to protect the robot structure and the force sensors from high impact forces during landing phases, two types of shock-absorbing mechanisms are commonly used.

The first is an internal passive mechanism. This type of shock-absorbing mechanism is present in the feet of many humanoid robots such as the Honda's Asimo [START_REF] Hirai | The development of Honda humanoid robot[END_REF] and the Kawada's HRP families [START_REF] Kaneko | IEEE International Conference on Robotics and Automation[END_REF]. The Asimo's shock-absorbing system smooths the impact force deforming rubber bushing inserted into guides. HRP robots have three rubber bushing with dampers that are attached between the foot and the leg. Unfortunately, this non-measurable compliance makes the attitude of the robot difficult to control, especially in complex maneuvers [START_REF] Vaillant | Multi-contact vertical ladder climbing with an HRP-2 humanoid[END_REF]. For this reason, a dedicated software-compensator is necessary to stabilize the walking [START_REF] Kajita | Balancing a humanoid robot using backdrive concerned torque control and direct angular momentum feedback[END_REF][START_REF] Nakaoka | Constraint-based dynamics simulator for humanoid robots with shock absorbing mechanisms[END_REF]. Bruneau et al. [START_REF] Bruneau | Dynamic walk of a bipedal robot having flexible feet[END_REF] explain the difficulty in modeling and controlling these flexibilities.

The second shock-absorbing system consists in using external passive mechanisms. Yamaguchi et. al use a multisensor foot mechanism with shock absorbing material [START_REF] Yamaguchi | Multisensor foot mechanism with shock absorbing material for dynamic biped walking adapting to unknown uneven surfaces[END_REF].

The foot sensors detect the path surface during the landing phase, which is combined with an adaptive biped walking control method and a flexible material to reduce the impact forces. Of course, combination of both is possible. In fact, all humanoid robots have a thin patterned rubber sole mounted (glued) on the part of their feet in contact with the ground. Because it is thin, large part of the landing shocks are absorbed at the ankle compliance.

In this paper, we aim at designing an optimized compliant sole added under the rigid structure of each foot to reduce the impact force at the landing phase, while accounting for the difficulty of the control after the landing phase (by increasing rotational stiffness). In addition, the soft soles must also be designed to absorb landing impact shocks and ground unevenness with a better contactsurface casting. During walking, the sole's compliance depends on the contact area variations. This external flexible sole can also offer other benefits such as better adherence properties and adapting, through deformation, to the terrain irregularities. We already mentioned this problem and reported very preliminary results in [START_REF] David | Study of an external passive shock-absorbing mechanism for walking robots[END_REF], Fig. 1. We decided to pursue this idea further. Our long-term objective is to optimize the sole together with the parameters of the walking pattern generator (WPG) to design soles that can adapt (in form and stiffness) to specific conditions and prerequisite needs for walking. In order to achieve this goal, we decided to proceed in three phases:

1. devise an approach by which we optimize the sole shape for given walking parameters and prerequisites, 2. extend the latter approach to combine tuning the walking parameters and sole shapes in a single optimization problem, 3. use the latter approach to obtain an optimized sole shape for different terrain conditions and requests.

This paper deals with the first phase of our work-plan. In Sec. 2, we explain our framework; in Sec. 3, we illustrate our optimization problem; in Sec. 4 and 5, we detail each module of the simulation framework and the methods used to define the sole geometry and the mechanical properties of the shape optimization problem; in Sec. 6 we describe the cost function of our optimization process; in Sec. 7 we present some solutions obtained using our approach.

Sole design framework

Shape optimization is very important in different engineering disciplines. When load carrying efficiency is tightly linked to the shape, the optimization of the structure is necessary. The common goal of a shape optimization problem is to find an object shape that minimizes an objective function and satisfies some constraints [START_REF] Allaire | Conception optimale de structures[END_REF][START_REF] Arora | Review of formulations for structural and mechanical system optimization[END_REF]. Many examples show the importance of shape optimization: for example the design of fillets, structures with holes [START_REF] Haftka | Structural shape optimization -A survey[END_REF], general 2-D bodies [START_REF] Espiga | Shape optimization of elastic homogeneous 2D bodies by the boundary element method[END_REF], curvilinear-stiffened panels [START_REF] Dang | Optimization and Postbuckling Analysis of Curvilinear-Stiffened Panels Under Multiple-Load Cases[END_REF], hull design [START_REF] Percival | Hydrodynamic optimization of ship hull forms[END_REF], electromagnetics [START_REF] Akcelik | Adjoint methods for electromagnetic shape optimization of the low-loss cavity for the International Linear Collider[END_REF], and automotive aerodynamics [START_REF] Muyl | Hybrid method for aerodynamic shape optimization in automotive industry[END_REF]. Particularly relevant are the studies on aerospace engineering [START_REF] Holst | Genetic Algorithms Applied to Multi-Objective Aerospace Shape Optimization[END_REF]. These shape optimization problems adopt different approaches to find optimal design of statically-loaded elastic structures. In our study, the force applied to the structure/object (i.e. the sole) depends on both time and space. The simulation of the foot/sole movement in a specific instant of time depends on the previous ones due to the sliding friction. To take into account these factors, we developed a dedicated shape optimization framework (Fig. 2) that optimizes the sole shape for a given walking pattern parameters, see Fig. 1a. Our aim is to find the sole's shape that minimizes a certain cost function computed over the simulation of the foot movement during the contact phase, while bounding the volume of the sole.

Recently, adaptive or refinement-based parametrization have been developed to improve the efficiency of shape optimization algorithms [START_REF] Kohli | Shape optimization using adaptive shape refinement[END_REF][START_REF] Desideri | Nested and selfadaptive Bezier parameterizations for shape optimization[END_REF][START_REF] Nagy | Isogeometric sizing and shape optimisation of beam structures[END_REF][START_REF] Han | An Evolutionary Geometry Parametrization for Aerodynamic Shape Optimization[END_REF]. Zingg et al. [START_REF] Zingg | Comparative Evaluation of Genetic and Gradient-Based Algorithms Applied to Aerodynamic Optimization[END_REF] show that gradient-based optimization is more efficient than gradient-free optimization in aerodynamics. This difference becomes more important with a higher number of design variables. In our optimization problem, the decision variables (i.e. the parameters) are the spline coefficients describing the sole deformation. The number of these coefficients is chosen by the user. Since we already envision extending our approach to a single optimization problem of the sole shape and the walking pattern generator at the same time, we expect the number of parameters to increase significantly. For this reason, we use a gradient-based formulation of the sole optimization problem.

The Fig. 2 illustrates our framework. The sole shape and the walking pattern are optimized using two distinct modules (Sec. 1). The Zero Moment Point (ZMP) position and force obtained from the walking pattern generator are the input to the simulation of the movement and the deformation. The simulation module will be detailed later.

Problem formulation

Our objective is to obtain a sole shape Ω, seen as a connected subset of R 3 , which lowers at best (i) the impact force during the heel-strike phase (shock-absorbing),

Shape optimizer

Initial shape

Optimized shape

Geometry parametrization

Update shape New shape Simulation of movement and deformation Walking Pattern Generator

Simulator X ini , p ini X opt p F des Z des s X act X act c(Ω)
Figure 2: Framework to optimize the sole shape for the given walking parameters. X ini , Xact and Xopt are respectively the initial, actual (during shape optimization process) and optimized sole nodes positions defining its shape. p ini and p are respectively the initial and intermediate (during shape optimization) sole shape parameters. s is the function transforming X ini in Xact. The Walking Pattern Generator (also formulated as an optimization problem) is an independent module that return the desired force F des and ZMP positions Z des as input to the simulation of the sole movement and deformation. The geometry parametrization module describes the shape of sole. The Shape optimizer finds numerically minimums of the cost function c(Ω).

and (ii) the foot rotation during the whole foot movement, while tracking a given resultant wrench of support reaction. We formulate this goal as an optimization program. We specify a pair of functions (F des (t), Z des (t)) representing the total contact forces and ZMP position over the time interval [t 0 , t f ] during which we want the sole to be in contact with the floor. We defined the cost function as follows, using the Cartesian stiffness matrix K c :

c(Ω) = th t0 σ max (K tr c (Ω, t)) dt -w tf t0 σ min (K rot c (Ω, t)) dt def = c tr -w c rot
(1) where σ max and σ min denotes respectively the largest and smallest singular value of a matrix, which implies we optimize the stiffness in the worst direction. The two terms correspond respectively to:

1. minimization of the translational part of the Cartesian stiffness matrix K c (maximize the compliance) during the heel-strike phase to reduce the impact force, and 2. maximization of the rotational part of the Cartesian stiffness matrix K c during the whole foot movement to decrease foot rotation in case of perturbations. Low foot rotation enforces the vertical posture and secures the balance of the robot during walking.

t h is the end of the heel-strike phase that we define as lasting 10% of the whole contact phase, when the vertical component of force almost reached its maximum, that is at the end of double support. Note that c and K c are parametrized by F des and Z des . We also specify a maximum volume v max (and hence the weight) for the sole, otherwise we would obtain very large and bulky soles. We thus have the following program:

min. Ω c(Ω), s.t. v(Ω) ≤ v max (2)
where v(Ω) denotes the volume of the sole.

Sole description

To make the problem (2) tractable, we need to describe the shape Ω with a finite number of parameters, and to simulate the movement and deformations undergone by a given shape. We do so by choosing a specific parametrization of the shape Ω, using Finite Element Model (FEM) to include the mechanical properties of the sole and a special algorithm to simulate the movement and deformations.

In the remaining of this paper, we divide the surface ∂Ω into two non-empty disjoint parts:

∂Ω = S D ∪ S S , with S D ∩ S S = ∅ (3)
S D is the surface of the sole attached to the rigid part of the foot (Dirichlet surface) and S S is the surface of the sole that could be in contact with the ground. We also denote by Ω I := Ω \ ∂Ω the interior volume of Ω.

Geometry parametrization

An important part of the shape optimization problem is the choice of basis in the parametrization. In this paper, we choose to describe the sole geometry by giving an initial shape that is deformed according to a function parametrized by a 2D spline. The coefficients of this spline are the optimization decision variables. The obtained shape is used to simulate the whole foot movement imposed by the walking pattern generator outputs.

Therefore, we describe Ω by (i) an initial volume Ω 0 and (ii) a deformation of this volume according to a parametrized function s. Fig. 3 gives a 2D illustration of this process.

The shape Ω 0 is discretized into finite elements and we define a point O s on S D (Fig. 1b), typically at its center.

To each node X i of the FEM with Cartesian coordinates (x i , y i , z i ) in a frame attached to O s , we associate spherical coordinates (r i , θ i , ϕ i ) with the following convention (Fig. 3):

O s θi θi θ i r θ ϕ x y z (x , y , z ) s s s (x , y , z ) i i i X i '=O s +s(θ i )(X i -O s ) O s X i X i ' s(θ i ) θ s( ) θ
r i = x 2 i + y 2 i + z 2 i θ i = atan2 (y i , -z i ) , ϕ i = atan2 x i , y 2 i + z 2 i
Then the coordinates of the node after deformation are

X i = s(p, θ i , ϕ i )X i (4)
We define s as a 2D B-spline with parameters p:

s(p, θ, ϕ) = n θ j=0 nϕ k=0 p j,k N j (θ)N k (ϕ) (5)
where the N j and N k are B-spline basis function of degree 2 scaled appropriately so that s is spatially defined on

- π 2 , π 2 2
. The variables n θ and n ϕ are user-defined; they control the precision of the spline. The spherical mapping has the singular points (the poles) for φ = ± π 2 . In proximity of the poles, close points X i have very different θ i and thus may have very different images X i through equation (4). To avoid this problem, we impose that in proximity of the pole coordinates the spline function is constant, by constraining the p j,0 (respectively n ϕ ) to have the same value.

We used a spherical parametrization, rather than a Cartesian one, since it allows a 2D shape parametrization of the surface, rather than a 3D one. The Cartesian parameterization would then have more parameters and would be redundant (tangential deformations to the surface of the body having no effect on the shape).

To avoid that the sole size is too large or too small, we bound the shape optimization parameters p j,k .

Finite Element Model

To simulate the deforming motion given by the flexible sole and obtain a desired ZMP force and position, we use FEM and contact modeling from mechanical laws (Sec. 5).

The simplest FEM is the static linear elastic model under the assumption of small deformations [START_REF] Bonnet | The finite element method in solid mechanics[END_REF]. This model is valid when the node displacements w.r.t the object size is relatively low and the sole deformation response is fast enough. In this paper, we assume that under a given force, the dynamic sole deformations are faster than the external dynamic forces and that the sole has small deformations.

In the static case for linear elasticity, a characteristic stiffness matrix is extracted from each element. These elementary matrices are then assembled in a large stiffness matrix K of the sole.

Once discretized by the finite element method, the displacements at each point of the mesh depend on external forces applied to the same points. Therefore, we can write the following linear relationship:

K l U = l F (6)
where l F ∈ R N is the vector of generalized forces (nodal forces) and l U ∈ R N is the vector of the node displacements. The superscript l before a letter denotes that the vector is expressed in the sole frame and no superscript means that the vector is expressed in the world frame.

The displacement of the Dirichlet nodes are known and for this reason, we can write the system (6) as:

  K DD K DI K DS K ID K II K IS K SD K SI K SS     0 l U I l U S   =   l F D l F I l F S   (7)
where I, D and S denotes respectively nodes of Ω I , S D and S S .

To simplify the mechanical model of sole, as in [START_REF] Bro-Nielsen | Real-time volumetric deformable models for surgery simulation using finite elements and condensation[END_REF], we add that internal nodes are not subject to external forces l F I = 0. Therefore, keeping the two last lines, the system becomes: From ( 8), we can write a condensed linear elasticity law [START_REF] Bro-Nielsen | Real-time volumetric deformable models for surgery simulation using finite elements and condensation[END_REF], to be used in first analysis of our algorithm:

K II K IS K SI K SS l U I l U S = 0 l F S (8)
l U S = K SS -K SI K -1 II K IS -1 l F S = K -1 S l F S l U I = -K -1 II K IS l U S (9)
where K S is the stiffness surface matrix.

As in [START_REF] Bro-Nielsen | Real-time volumetric deformable models for surgery simulation using finite elements and condensation[END_REF], we consider that external forces apply to the sole surface nodes only. Hence, to find the displacement of each nodes, we divide the simulation into two steps:

1. find the displacement l U S induced by the contact with the ground (Sec. 5), and then 2. compute the displacement l U I using the second equation in (9).

Contacts

We describe how we model the contact between the sole and the (flat) ground. To model the contact, we use a relationship between the two unknowns of the system (9): the contact force and the relative position between the sole node position and the ground position (Fig. 4). This relative position is called contact space displacement vector δ. For a point contact α, the normal component of δ is:

δ n,α = n (P α -Q α ) (10)
where n = 0 0 1 is the normal vector. We denote with the subscript n the normal component of the vector v n = nv and with the subscript t the tangential components of the vector v t = tv, where t = 1 0 0 0 1 0 .

We need to compute the forces F applying to the contact points P , which impose non-interpenetration between the sole and the ground. We use the two following laws:

1. Signorini's law establishes for each contact node α the complementarity between F α and δ α along each contact normal [START_REF] Signorini | Sopra alcune questioni di elastostatica[END_REF]:

0 ≤ F n,α ⊥ δ n,α ≥ 0 (11)
2. Coulomb's law characterizes the dry friction. If F α is the contact force applied to the contact node α and µ the friction coefficient, we can distinguish two complementarity conditions for stick and slip motions:

(a) F t,α < µ|F n,α | and δ t,α = 0 (stick condi- tion) (b) F t,α = -µF n,α δt,α δt,α (slip condition)
Based on the framework in Fig. 2, we want to find a possible position and orientation of the sole together with the configuration of the contact points to obtain F des and Z des given by the walking pattern generator. Our algorithm is illustrated in Fig. 5.

Frictional contact problem (algorithm 1)

Walking pattern generator Sole Position at q -1 equation ( 25)

Update sole position and orientation (equation ( 24)) equation ( 27)

(F des , Z des ) (O q-1 l , Υ q-1 ) (δ, F ) NO YES (O l , Υ) YES NO P Figure 5
: Framework for contact handling to obtain F des and Z des . q -1 is the previous time step.

Signorini's law

To find the contact space displacements and forces, we can formulate the contact problem as a linear complementarity problem (LCP). Various methods, direct or iterative are possible to solve LCP [START_REF] Murty | Linear Complementarity, Linear and Nonlinear Programming[END_REF]. For our applications, we use an iterative Gauss-Seidel solver based on [START_REF] Duriez | Realistic haptic rendering of interacting deformable objects in virtual environments[END_REF]. Using this contact resolution, LCP takes into account the coupling between the different contact points.

Given the foot orientation and position defined respectively by three angles Υ = (θ, φ, ψ) and a vector O l , the sole surface node positions in the absolute frame P free are given by the following relationship (Fig. 4):

P free = [O l • • • O l ] T + H l P free 1 • • • l P free m T def = O l + H l P free (12
) where l P free is the vector of the free node positions in the sole frame. H = diag(R, . . . , R), where R is the foot orientation matrix defined by three Euler rotation angles Υ = (θ, φ, ψ). The subscript m is the number of nodes of suface S S of the sole that can potentially be in contact with the ground. The superscript T denotes the transpose operator.

When the sole comes into contact with the ground, we can express the displacement U S as:

U S = P -P free (13)
Using the linear relationship (9) tranformed in the global frame, the position of each surface node is composed of two terms:

P = HK -1 S H T F + P free def = W F + O l + H l P free (14)
W is the Delassus operator [START_REF] Moreau | Numerical Treatment of Contact and Friction: The Contact Dynamics Method[END_REF]. By definition of δ α , linking a surface node to a reference point Q α on the ground, we can write:

δ = P -Q = W F +O l +H l P free -Q def = W F +δ free ( 15
)
where

Q = Q 1 . . . Q m .
Once δ is computed as presented next, l U S is obtained by:

l U S = H T U S = H T (δ + Q) -l P free (16)
Then, the normal component of the displacement vector [START_REF] Vaillant | Multi-contact vertical ladder climbing with an HRP-2 humanoid[END_REF] for each contact point α is:

δ n,α = nO l + nR l P free α + n m β=1 (W αβ F β ) def = δ free n,α + n m β=1 (W αβ F β ) (17)
where

W αβ def = W t W n def = W tt W tn W nt W nn is the 3 × 3 Delas-
sus' operator coupling the contacts α and β. Q α does not appear in [START_REF] Nakaoka | Constraint-based dynamics simulator for humanoid robots with shock absorbing mechanisms[END_REF] since we consider the absolute frame on the flat floor, and hence nQ α = 0. Using [START_REF] Aerts | The mechanical properties of the human heel pad: a paradox resolved[END_REF] and ( 17), the Signorini's problem can be solved using a LCP in the frictionless case.

Coulomb's law

We now explain how we add the Coulomb's law. In the 3D case, this law adds a non-linearity to the initial contact problem. To solve friction contact cases, two approaches are commonly used in the literature: (i) keeping the problem linear by using a multi-sided pyramid approximates of the friction cones [START_REF] Murty | Linear Complementarity, Linear and Nonlinear Programming[END_REF] and (ii) try solving the non-linear problem like in [START_REF] Jean | Numerical methods for three dimensional dynamical problems[END_REF].

Using the first method, we can stack multiple contacts with friction in a global LCP. However, if we choose to model the friction with pyramids of 4 facets, each contact generates six new lines in the LCP and the computational cost of the problem can become very high. In addition, the accuracy is very low with 4 facets pyramids.

The non-linear approach in [START_REF] Jean | Numerical methods for three dimensional dynamical problems[END_REF] is based on the Newton's method. For the paper to be self-contained, we explain this method for one node in contact (we also derive gradients from them later-on).

From [START_REF] Vaillant | Multi-contact vertical ladder climbing with an HRP-2 humanoid[END_REF], we find the tangential component of the contact space displacement:

δ t,α = tO l,α + tR l P free α + t m β=1 (W αβ F β ) -tQ α ≈ tO l,α + tR l P free α + t m β=1 (W αβ F β ) -tP (q-1) α def = δ free t,α + t m β=1 (W αβ F β ) ( 18
)
where t is the tangential vector, P (q-1) α is the node position in absolute frame at the previous time step and q is the current time step. In case of a node already in contact at q -1, the contact reference point is the previous node position P (q-1) α . In case of a new contact, we project on the ground the node position P (q-1) α as an approximation of the reference contact position. In fact, real contact takes place at an intermediate time between previous and current time step, and a node may have moved a little in the tangent direction. Precise contact position could be found, but at the price of a variable step size and event driven simulation, which is computationally costly when many contacts occur. Additionally, this approximation is good enough when the time-step size is sufficiently small. From ( 17) and ( 18), we obtain:

δ α = O l + R l P free α - P (q-1) t,α 0 + m β=1 W αβ F β def = δ free α + m β=1 W αβ F β (19)
Using ( 19) and Coulomb's law, we build for each contact α an application Φ ∈ R 2 × R × R 2 × R. This method to solve friction contacts is taken by [START_REF] Jean | Documentation sur LMGC: Logiciel de Mécanique Gérant le Contact[END_REF][START_REF] Duriez | Realistic haptic rendering of interacting deformable objects in virtual environments[END_REF]. In this section, we detail these equations below to compute the gradient of our algorithm in Sec. 5.4.

           Φ 1 (δ α , F α ) = δ t,α -δ free t,α -W tn F n,α -W tt F t,α Φ 2 (δ α , F α ) = δ n,α -δ free n,α -W nn F n,α -W nt F t,α Φ 3 (δ α , F α ) = F t,α - µF n,α (F t,α -q t δ t,α ) F t,α -q t δ t,α Φ 4 (δ α , F α ) = q n δ n,α (20)
The first three equations (two for Φ 1 and one for Φ 2 ) are given by [START_REF] Yamaguchi | Multisensor foot mechanism with shock absorbing material for dynamic biped walking adapting to unknown uneven surfaces[END_REF] and the last three equations (two for Φ 3 and one for Φ 4 ) by the Coulomb's law. To solve [START_REF] David | Study of an external passive shock-absorbing mechanism for walking robots[END_REF], we need to find a couple (δ α , F α ) such as Φ(δ α , F α ) = 0.

Alart and Curnier [START_REF] Alart | Time-stepping for threedimensal rigid body dynamics[END_REF] explains that a Newton method can be applied to this type of function. Detail about the convergence are given in [START_REF] Alart | Critères d'injectivité et de surjectivité pour certaines applications de R n dans lui-même; application à la mécanique du contact[END_REF]. Based on [START_REF] Jean | Documentation sur LMGC: Logiciel de Mécanique Gérant le Contact[END_REF], we choose the coefficients q n and q t as:

q n = 1 W nn , q t = λ min λ 2 max ( 21
)
where λ min and λ max are respectively the min and max eigenvalue of the matrix W tt .

To solve [START_REF] David | Study of an external passive shock-absorbing mechanism for walking robots[END_REF], we also use the Newton's method:

δ (r) α F (r) α = δ (r-1) α F (r-1) α -∂Φ δ (r-1) α , F (r-1) α -1 Φ δ (r-1) α , F (r-1)
α
where ∂Φ is the Jacobian matrix of the application Φ and it is illustrated in AppendixA.

If the first derivatives are sufficiently smooth and the initial point is not too far from one of the roots of the equations, the convergence of the Newton's Method is quadratic.

Gauss-Seidel method for contacts

In multi-contact case (m contacts), we need to take into account the coupling between the contact nodes. Using [START_REF] Yamaguchi | Multisensor foot mechanism with shock absorbing material for dynamic biped walking adapting to unknown uneven surfaces[END_REF], we have 3m equations with 3m force unknowns and 3m movement unknowns.

Rewriting [START_REF] Yamaguchi | Multisensor foot mechanism with shock absorbing material for dynamic biped walking adapting to unknown uneven surfaces[END_REF] as:

δ α -W αα F α = m β=1,β =α W αβ F β + δ free α ( 22
)
As in [START_REF] Duriez | Realistic haptic rendering of interacting deformable objects in virtual environments[END_REF], we can use an iterative Gauss-Seidel method to obtain a real-time solution of the contact problem under the Signorini's and Coulomb's laws. This method is well-adapted to sparse matrices and it converges, at least locally, when the system is over-constrained and there are multiple solutions. The convergence of this method is shown in [START_REF] Jourdan | A Gauss-Seidel like algorithm to solve frictional contact problems[END_REF].

The principle is to test the contacts one by one, and always in the same order, blocking all others to their most recently calculated value (algorithm 1). 2 , 3 are respectively the desired accuracy for Signorini's law and Coulomb's law.

Algorithm 1: F = contact(F , δ free ) Data: W , µ Parameters: 2 , 3 for 1 ≤ α ≤ m do δ α ← m β=1,β =α W αβ F β + δ free α if δ n,α < 2 then F α ← -W -1 αα δ α if F α,t > µ F α,n then F α,n ← -δ free α,n /W nn δ α ← W αα F α + δ free α do c ← [∂Φ (δ α , F α )] -1 [Φ (δ α , F α )] δ α F α ← δ α F α -c while Φ (δ α , F α ) > 3 else F α ← 0
To find the appropriate contact node forces and then the node positions for a given foot 6D position (translation and rotation), we can use the frictional contact algorithm 3 using Gauss-Seidel method [START_REF] Duriez | Realistic haptic rendering of interacting deformable objects in virtual environments[END_REF]. q is the time step and r is the total number of simulation steps; 1 is the desired accuracy for the contact problem. The initialization procedure of the frictional contact algorithm is described in algorithm 2.

Algorithm 2: (F , δ free ) = init(O l , Υ) Data: O ini l , l P free , P (q-1) , F (q-1)

compute R(Υ)

if q = 1 then for α = 1 . . . m do F α ← 0 δ free α ← O l + -tO ini l nR l P free α if q = 2 . . . r then for α = 1 . . . m do F α ← F (q-1) α δ free α ← R l P free α + O l - P (q-1) t,α 0 Algorithm 3: (F , δ) = AlgoContacts(O l , Υ) Parameters: 1 (F , δ free ) ← init(O l , Υ) //algorithm 2 do F old ← F F ← contact(F , δ free ) //algorithm 1 while F -F old F > 1
compute δ using equation 19

Desired force and ZMP position algorithm

In the previous section, we detailed the frictional contact problem where the inputs are the foot position and orientation and the outputs are the contact forces and the contact space displacements. In this section, we propose an algorithm to compute foot position and orientation to obtain the desired F des (t) and Z des (t) and null vertical torque at ZMP Γ Z des ,n = 0 mentioned in Fig. 5.

Let f be the function computing the resultant force F tot , the ZMP Z and vertical torque at desired ZMP Γ Z des ,n from (F , δ). We want to solve the equation

f (AlgoContacts(O l , Υ)) =   F des (t) Z des (t) 0   (23)
We solve it by Newton method. Therefore, we derive a Newton step by finding a differential relationship between (F tot , Z, Γ Z des ,n ) and (O l , Υ):

  F des -F tot Z des -Z -Γ Z des ,n   =   dF tot dZ dΓ Z des ,n   = J dO l dΥ (24)
where dΥ = dθ dφ dψ T ; J is the Jacobian matrix that links ZMP Force and position with foot orientation and position.

In AppendixB, we detail the terms of equation [START_REF] Espiga | Shape optimization of elastic homogeneous 2D bodies by the boundary element method[END_REF]. The method to obtain the node positions for a desired ZMP position and force is illustrated in algorithm 4. 4 is the desired accuracy for ZMP position and force. In algorithm 4, we relax the convergence of AlgoContact 3 through a user-defined parameter λ ∈ [0 1]:

F -F old F > λ F des -F tot Z des -Z -Γ Z des ,n (25)
Using this relaxed method, we reduce the simulation time by 30%. The convergence of our algorithm is achieved when the following condition is met:

F -F old F < 1 (26)
and: 27) is true compute U using equation ( 9) compute P free using equation ( 12) P ← U + P free

F des -F tot Z des -Z -Γ Z des ,n < 4 (27)
dO l dΥ ← J -1   F des -F tot Z des -Z -Γ Z des ,n   O l ← dO l + O l Υ ← dΥ + Υ while equation (

Initial conditions

To solve the problem, we have to set the condition at the first time step. We impose an initial position O ini l and yaw orientation of the foot just before the first contact that will be used for the reference contact positions Q. Initial roll and pitch orientation of the foot just before first contact is not imposed since it could be discontinuous with the roll and pitch orientation found from desired ZMP position. Initial roll and pitch orientation of the foot are then taken to be the same as those obtained from the desired ZMP position. Hence, tangential reference contact positions are:

Q t = tO ini l + tR ini (θ, φ, ψ ini ) l P free α ( 28
)
where ψ ini is the desired fixed yaw orientation of the foot at first contact with the ground. For the normal component, δ ini n,α = δ n,α . For the tangential term, equation ( 18) becomes:

δ ini t,α = tO l,α + tR l P free α + t m β=1 (W αβ F β ) -tQ α ≈ t(O l,α -O ini l,α ) + t(R -R ini ) l P free α +t m β=1 (W αβ F β ) (29)
From ( 17) and ( 29), we obtain:

δ α = O l,α + R l P free α - t(O ini l,α + R inil P free α ) 0 + m β=1 W αβ F β (30)
The Jacobian matrix for the first step is detailed in AppendixC.

Cost function

As explained in Sec. 3, our cost function is based on the Cartesian stiffness computed in the ZMP position. Here, we explain how we can find this matrix. For an iterative algorithm we need to find a differential relationship between (F tot , Γ Z ) and ( z O l , ω), where ω is angular velocity and z denotes the ZMP frame. We consider z O l rather than O l so that reference point for small rotation is Z rather than the absolute origin. Defining a small rotation dω = ω dt, we want the following matrix K c :

dF tot dΓ Z = K c z dO l dω = K tr c K tr,rot c K rot,tr c K rot c z dO l dω (31)
The detailed description with computational details of the Cartesian stiffness is in AppendixD.

Results

We first present the conditions of the sole optimization problems, followed by a presentation of the obtained soles and discussions.

Optimization description

We used smooth ZMP and force trajectories of Fig. 6 obtained with the walking pattern generator [START_REF] De Magistris | Optimized humanoid walking with soft soles[END_REF] based on a minimization of the robot energy consumption. Pattern generator was used on a single cyclic step. To obtain the whole ZMP trajectory under each foot, the ZMP repartition is computed during the double-support-phase of the walking gait. The ZMP trajectory is defined as a twice differentiable piecewise polynomial function whose parameters are obtained by solving a QP optimization problem. The force trajectory is a by-product of the computation. More details are given in [START_REF] Pajon | A humanoid walking pattern generator for sole design optimization[END_REF] and [START_REF] De Magistris | Optimized humanoid walking with soft soles[END_REF].

Table 1 summarizes the different cases used to obtain optimized sole shapes. They were chosen to study the influence of different parameters on the sole shape obtained.

Here is a detailed description of studied cases:

1. A reference case that can be applied to the robot. 2. A case giving more weight to the rotation criterion for better stability during the walk. Used to see the influence of relative weight given to impact force and foot rotation reduction. 3. A case with a different material, to see the influence of the choice of material on the sole shape and criteria. 4. A case with a different ZMP trajectory simply going from heel to toe in the middle of the foot, to see the influence of the ZMP trajectory. 5. A case with a non-modified surface attached to the rigid lower part of the robot foot S D (Dirichlet nodes with no displacement) during optimization, as opposed to the other cases where the shape of this surface of Dirichlet nodes is modified and optimized (but kept flat). In this case this surface remains constant during optimization, in order to obtain a sole that can be used on real robot without changing the existing rigid part of the foot. This case allows to see the cost of this constraint on the foot shape and criteria.

To obtain dynamic properties close to the original robot set-up, we impose that the existing robot foot weight with ankle flexibility is similar to our proposed new foot with the flexible sole. This is enforced by the constraint on the total volume in problem (2)1 .

Table 1: Different cases for the different optimized sole shapes. The Young modulus of neoprene is 1 MPa and of butyl rubber is 1.4 MPa [START_REF]Part I of the Engineering Tripos[END_REF]. The cost function is the equation ( 1); WPG denotes the ZMP force and trajectory obtained using the walking pattern generator in [START_REF] Pajon | A humanoid walking pattern generator for sole design optimization[END_REF].

Case WPG Material

Cost Dirichlet In the case of the neoprene for a maximum sole weight of 1 kg, the maximum volume is:

v max = m max ρ = 1 kg 1230 kg m =3 = 813 cm 3 ( 32
)
where ρ is the material density [START_REF]Part I of the Engineering Tripos[END_REF]. In the case of butyl rubber v max = 1098 cm 3 . Starting from the initial shape in Fig. 1b, we obtain convergence to an optimized sole shape. We then refine the optimized shape using iso2mesh's library2 and we run a second optimization with the refined mesh as initial shape, etc. During the shape optimization process, the mesh could be invalid because of some reversed tetrahedrons (with a negative volume) [START_REF] Hecht | FreeFem++ Manual[END_REF]. In this case, we stop automatically the optimization, remesh and restart. To obtain one optimized sole, our optimization takes between 1 day and 1.5 days of computation with a partially optimized code on a personal laptop. The library iso2mesh is very useful in our case because it is able to create quality volumetric mesh from isosurface patches. Fig. 8, 9, 10, 11, 12 shows the optimized soles for the different cases in Table 1 and a dry friction of 0.8.

Discussions

Notice for all optimized soles, that to minimize the impact during the heel-strike phase, the optimizer creates a flexible appendix with minimum connection with the rigid part of the foot, in a shape similar to the human heel. This "heel" part is different for the different cases in Table 1. In particular, notice that for the trajectory in Fig. 6, this part is in the back-left part of the sole, for the trajectory in Fig. 7 it is in the back-center part of the sole. One can conclude that the optimizer creates a heel where the ZMP trajectory is along the sole.

The difference of the heel part for different cost functions is also remarkable. In the case 1-Ref, this part is more rounded, shorter and with less connection with the rigid part of the foot in comparison with the optimized sole in the case 2-Rot. This interesting result shows that if we want to enforce lower foot rotation to keep the vertical posture, we need a less rounded, longer heel with rigid support. Those characteristics to minimize foot rotation are also observed at the front of the optimized sole, where only rotational part of the Cartesian stiffness is maximized. We then obtained a large front part of the optimized sole similar to front part of human feet.

To fulfill the volume constraint of the problem (2), the optimizer whittle many sole parts down. Especially, less material is present in the middle x-axis, which increases rotational stiffness along x-axis.

More precisely, from Fig. 13 and Table 2, comparing the different cases, one can observe the followings:

1. For a given sole, Cartesian stiffness in translation and rotation are rather linked, since their evolution with respect to time is similar in Fig. 13. Hence minimizing translational stiffness and maximizing rotational stiffness are antagonistic objectives. Different compromises can be found between both criteria. Since translational stiffness is minimized only at the beginning, we can observe that stiffnesses are lower at the beginning and larger later. 2. For the case 1-Ref compared to the initial sole, c tr has been divided by almost 5, while c rot is about the same. So shape optimization improves one of the two sole criteria by a large factor. 3. For the case 2-Rot compared to the case 1-Ref, the rotational stiffness is improved while the impact absorption is worsen, which is another compromise with more importance given to maximizing the rotational stiffness, as expected with increased weight on this criterion. 4. For the case 3-Rubber compared to case the 1-Ref, what is lost in translational stiffness minimization is gained in rotational stiffness maximization, which is another compromise. Criteria for this butyl-rubber sole are another compromise close to neoprene sole, so this change of material has no significant impact on the sole performance. 5. For the case 4-Straight compared to the case 1-Ref.

Both c tr and c rot are worse. Having a straight ZMP trajectory seems worse than a trajectory starting and finishing in a corner of the foot. This is caused by a poor rotational stiffness around the y-axis due to an entirely soft heel. It shows that ZMP trajectory has a strong influence on sole shape optimization. In future works, we plan to optimize the sole shape and the ZMP trajectories simultaneously. 6. For the case 5-Fixed compared to the case 1-Ref, c tr and c rot are about 40% worse. It shows that it is worth optimizing the rigid part of the foot and redesign robot feet accordingly.

Additionally, as explained in Sec. 4.1, we bound the shape optimization parameters p j,k . This bound in combination with the initial parallelepiped shape can create some sharp shapes in the optimized soles. We notice that for example, in the front part in Fig. 8, 9, 10, 11, 12.

These results can be visualized in the accompanying video available at https://youtu.be/wivyn5klAMw. This video presents the deformation simulations of the soles. One can notice that the deformation of the heel appendices at the beginning of foot contact are undergoing deformations of the order of cm, with rotations of about 45 • . The small deformation hypothesis is compromised, and results in this case are not precise. For more precise results, we plan to use a large-deformation FEM code in the future. 1.

Conclusion

In this paper, we devised a framework to optimize the sole shape of a humanoid robot for a given cost function.

To obtain faster and easier to implement algorithms, some model approximations have been done. In particular, we used a simplest FEM for static linear elastic model under the assumption of small deformations to model the soles. The obtained soles appear having some similar characteristics with the human foot. We also showed that the shape optimization improves significantly the criteria considered for the sole, and that the ZMP trajectory and the shape of the contact with rigid part of the foot have a strong influence on the performance of the sole, but not the different material tested.

The next step of this research is to extend the latter approach to optimize the walking pattern parameters and sole shape in a single optimization problem formulation. This proved to be pertinent and was investigated in simpler 2D predefined forms in [START_REF] Yamane | Effect of foot shape on locomotion of active biped robots[END_REF].

In the papers [START_REF] De Magistris | Humanoid walking with compliant soles using a deformation estimator[END_REF][46], we investigated adding a plain rectangular compliant sole to a humanoid robot. In the near future, we will build up and test the optimized soles on a real robot to validate the procedure and criteria considered. We will also use more accurate FEM, for example considering large-displacement and nonlinear material behaviors.

In the current state of our framework, the gradient of the cost function is obtained by finite differences, which impose numerous simulations of the sole movement for each iteration of the shape optimization process. We are currently implementing the analytical gradient. We expect this will speed-up the presented computation by a factor 20-100.

After those improvements, we also plan to develop a controller for stabilization with a sole prior to experiments on the HRP-4 humanoid robot.

Acknowledgments

The research leading to these results has received funding from the EU 7th Framework Program (FP7/2007-2013) under grant agreement no 611909 (KoroiBot).

AppendixA. Jacobian of the application Φ

∂Φ =       I 2×2 0 2×1 -W tt -W tn 0 1×2 1 -W nt -W nn ∂Φ 3 ∂δ t,α 0 2×1 ∂Φ 3 ∂F t,α ∂Φ 3 ∂F n,α 0 1×2 q n 0 1×2 0       (A.1)
in which:

∂Φ 3 ∂δ t,α = q t µF n,α Π (F t,α -q t δ t,α ) , ∂Φ 3 ∂F t,α = I 2×2 -µF n,α Π (F t,α -q t δ t,α ) , ∂Φ 3 ∂F n,α = - µ (F t,α -q t δ t,α ) F t,α -q t δ t,α Π(V ) = ∂ ∂V V V = 1 V 3 (v 2 ) 2 -v 1 v 2 -v 1 v 2 (v 1 ) 2
AppendixB. Relationship between ZMP forceposition and foot position-orientation From (24), we can straightforwardly obtain F tot in function of the force applied at the contact node α:

F tot = m α=1 F α (B.1)
Then for Z and Γ Z des ,n we start with equation of resultant torque at a point A given by:

Γ A = m β=1 (P β -A) × F β (B.2)
Using the definition of the ZMP and the vertical resultant torque at Z des we get:

  Γ Z,t1 = 0 Γ Z,t2 = 0 Γ Zdes,n   = m β=1   (P t2,β -Z t2 ) • F n,β -(P t1,β -Z t1 ) • F n,β (P t1,β -Z des,t1 ) • F t2,β -(P t2,β -Z des,t2 ) • F t1,β   (B.
3) The subscripts t 1 and t 2 denote respectively the first and second component of the vector.

We now have relations (B.1) and (B.3) between (F tot , Z, Γ Z des ,n ) and nodes position and force, and relation ( 19), Signorini's and Coulomb's laws between nodes position and force as well as sole position and orientation. To find the differential relationship [START_REF] Espiga | Shape optimization of elastic homogeneous 2D bodies by the boundary element method[END_REF], we have to take the differential of previously cited equations.

Considering F α , δ α , O l and (θ,φ,ψ) as variables, the differential of ( 19) is:

dδ α = dO l + m β=1 W αβ dF β + ν θ ν φ ν ψ dΥ (B.4)
where ν x for each angle x ∈ {θ, φ, ψ}:

ν x = ∂R ∂x l P free α + m β=1 ∂R ∂x C αβ R T + RC αβ ∂R ∂x T F β
Following the Coulomb's law, we can distinguish two cases for each contact node α:

1. Stick: We have three equations given by (B.4) and three equations given by δ α = 0 → dδ α = 0. 2. Slip: We have three equations given by (B. [START_REF] Mombaur | From human to humanoid locomotion-an inverse optimal control approach[END_REF] and three equations given by the differential of the Coulomb's law (A.1).

Using those considerations, we can write the following system that relates contact forces and contact space displacements to the foot position and orientation:

A 11 A 12 A 21 A 22 dF dδ = B 11 B 12 0 3s×3 0 3s×3 dO l dΥ (B.5)
where: dδ includes only the slipping contacts, m is the number of contacts, s is the number of slipping contacts, p(j) j ∈ [1, ..., s] is the j th slip contact and:

                                             A 11 = -W , A 12 = a 1 • • • a j • • • a s , a j = 0 1 3×3 • • • I p(j) 3×3 • • • 0 s 3×3 T A 21 = c T 1 • • • c T j • • • c T s T , c j =   0 1 3×3 • • •   ∂Φ 3 ∂F t,j ∂Φ 3 ∂F n,j 0 1×2 0   p(j) • • • 0 m 3×3   A 22 = diag (d j ) , d j =   ∂Φ 3 ∂δ t,j 0 2×1 0 1×2 q n   B 11 = I • • • I T B 12 = b 1 • • • b i • • • b m T , ∀i ∈ [1 m] b i = ν θ i ν φ i ν ψ i T
The differential of (B.3) gives: where for p(j) j ∈ [1 s] is the j th slip contact: Considering Γ Z des ,n rather than Γ Z,n gives a simpler differential (B.6) since Z des is constant. Since Z converges to Z des , the same solution is obtained in the end.

                                         D 11 = I 1 3×3 • • • I i 3×3 • • • I m 3×3 , ∀i ∈ [1 m] D 21 = d 1 • • • d i • • • d m , ∀i ∈ [1 m] d i =      0 0 1 F tot,n (P it1 -Z t1 ) 0 0 1 F tot,n (P it2 -Z t2 ) (-P it2 + Z t2 ) (P it1 -Z t1 ) 0      D 22 = e 1 • • • e j • • • e s , j

AppendixC. Jacobian for the first step

Considering F α , δ α , O l and (θ,φ,ψ) as variables, the differential of ( 30) is:

dδ α =

AppendixD. Details of the cost function

To obtain the equation [START_REF] Desideri | Nested and selfadaptive Bezier parameterizations for shape optimization[END_REF], we start by differentiating [START_REF] Yamaguchi | Multisensor foot mechanism with shock absorbing material for dynamic biped walking adapting to unknown uneven surfaces[END_REF] with respect to t instead of (θ, φ, ψ). To do so, we use the following derivation of a vector applied on a rotation matrix dR = dωR (D.1)

where dω denotes the skew matrix of dω. (D.1) applied to the differential of [START_REF] Yamaguchi | Multisensor foot mechanism with shock absorbing material for dynamic biped walking adapting to unknown uneven surfaces[END_REF] gives after factorization in dω: where A ij are same as in (B.5) and:

         B ω11 = I 1 • • • I i • • • I m T ∀i ∈ [1 m] B ω12 = b 1 • • • b i • • • b m T ∀i ∈ [1 m] b T i = -( RP free i ) + m β=1 W iβ F β -( W iβ F β )
Using differential of (B.1) and (B.2) with (D.3), we have the following relationship:

dF tot dΓ = G dF dδ = G 11 0 3×3s G 21 G 22 A -1 B ω dO l dω (D.4)
where for p(j) j ∈ [1 s] is the j th slip contact:

                           G 11 = I 1 3×3 • • • I i 3×3 • • • I m 3×3 , ∀i ∈ [1 m] G 21 = g 1 • • • g i • • • g m , ∀i ∈ [1 m] g i =   0 0 (P it2 -Z it2 ) 0 0 -(P it1 -Z it1 ) -(P it2 -Z it2 ) (P it1 -Z it1 ) 0   G 22 = ζ 1 • • • ζ j • • • ζ s , ζ j =   0 F p(j)
,n 0 -F p(j),n 0 0 F p(j),t2 -F p(j),t1 0   Using (D.4) and the transport law for the screw, we can find the Cartesian stiffness as:

dF tot dΓ Z = GA -1 B ω dO l dω = GA -1 B ω Ξ z

Figure 1 :

 1 Figure 1: (a): Rectangular parallelepiped soles mounted on HRP-4's feet; (b): Initial sole shape where Os is the center of the surface attached to the foot

Figure 3 :

 3 Figure 3: Geometry parametrization module in 2D case at the left; convention for the 3D spherical coordinates at the right.

Figure 4 :

 4 Figure 4: Contact between the sole and the ground.

Figure 6 :Figure 7 :

 67 Figure 6: ZMP trajectory at left and ZMP force at right for the energy-minimizing case.

Figure 13 :

 13 Figure 13: Cost functions for the five cases in Table1.

Figure 8 :Figure 9 :Figure 10 : 1 )Figure 11 :Figure 12 :

 891011112 Figure 8: Optimized sole for the case 1-Ref: ZMP trajectory illustrated in Fig.6, material neoprene and w = 1 in the cost function[START_REF] Wieber | Modeling and control of legged robots[END_REF]

∂PF

  t2,β ∂δ t2,β dδ t2,β F n,β + m β=1(P t2,β -Z t2 )dF n,β -dZ t2 m β=1 F n,β = 0 -m β=1 ∂P t1,β ∂δ t1,β dδ t1,β F n,β -m β=1 (P t1,β -Z t1 )dF n,β + dZ t1 m β=1 F n,β = 0 m β=1 ∂P t1,β ∂δ t1,β dδ t1,β F t2,β + m β=1 (P t1,β -Z des,t1 )dF t2,β m β=1 ∂P t2,β ∂δ t2,β dδ t2,β F t1,β -m β=1 (P t2,β -Z des,t2 )dF t1,β = dΓ Zdes,n(B.6)With ∂P t,β ∂δ t,β = I (obtained by differentiating (15)) it can be presented in the following matrix form:  des -F tot Z des -Z -Γ Z des ,n

  .5) with (B.7), we finally get J in (24):J = DA -1 B (B.8)

W 6 )

 6 dO l,α + m β=1 αβ dF β + ν θ modification of (B.5) for the initial step only for b i as follows:b i = ν θ i -For the first time step, the differential of (15) gives:dP β = dδ β + t ∂R ini∂θ l now have additional terms that finally gives initial step J ini = J + J in[START_REF] Espiga | Shape optimization of elastic homogeneous 2D bodies by the boundary element method[END_REF], where:

dδ α = dO l + m β=1 W

 β=1 αβ dF β -( RP free α )dω + m β=1 W αβ F β -( W αβ F β ) dω (D.2)Therefore in matrix form we have:A 11 A 12 A 21 A 22 dF dδ = B ω11 B ω12 0 3s×3 0 3s×3 dO l dω (D.3)

  n -O l,t2 + Z t2 -O l,n 0 O l,t1 -Z t1 O l,t2 -Z t2 -O l,t1 + Z t1 0   (D.6)

  ). c tr opt is minimized, so increase is worse. c rot opt is maximized, so increase is better and decrease is worse.

	Case	c tr ini	c tr opt	c tr opt comp.	c rot ini	c rot opt	c rot opt comp.	c ini	c opt
	1-Ref	3.88 • 10 6 8.06 • 10 5	/	1.79 • 10 5 1.74 • 10 5	/	3.70 • 10 6	6.33 • 10 5
	2-Rot	3.88 • 10 6 4.65 • 10 6	×5.77	1.79 • 10 5 3.50 • 10 5	×2.01	-1.75 • 10 8 -3.45 • 10 8
	3-Rubber	5.39 • 10 6 9.59 • 10 5	+19.0%	2.48 • 10 5 2.14 • 10 5	+23.0%	5.15 • 10 6	7.45 • 10 5
	4-Straight 9.14 • 10 6 1.02 • 10 6	+26.6%	3.41 • 10 5 5.84 • 10 4	/2.98	8.80 • 10 6	9.62 • 10 5
	5-Fixed	3.88 • 10 6 1.21 • 10 6	+50.1%	1.79 • 10 5 1.16 • 10 5	-33.3%	3.70 • 10 6	1.09 • 10 6

This constraint is implemented as an inequality instead of an equality; yet this is equivalent at the solution, because the constraint is active (but an inequality enhances the convergence).

http://iso2mesh.sourceforge.net/


T Levi

A Tixier-Mita

B D Ségard

H Toshiyoshi

H Fujita

T Fujii

BIOMIMETIC MICROFLUIDIC NEURONS FOR HYBRID EXPERIMENTS

Keywords: Neuromorphic engineering, Neuron, Biomimetic

This paper reports a new way to explore in the neuromorphic engineering, the biomimetic artificial neuron using microfluidic techniques. This device is able to mimic the electrical activity of one biological neuron. Usually these artificial neurons are made in Silicon but this device could replace the electronic one and solve most of the issues of biocompatibility. This PDMS device is composed of two chambers for intra and extra-cellular modelling, different PDMS channels, selective permeable membrane for positive ionic exchange, quake valves and electrodes for recording the membrane potential. We obtain an electrical membrane potential similar to the biological neuron.

INTRODUCTION

The neuromorphic engineering is the design of biomimetic artificial neural systems. Mostly of these systems are silicon-based [START_REF] Mahowald | A Silicon neuron[END_REF], [START_REF] Le Masson | Feedback inhibition controls spike transfer in hybrid thalamic circuits[END_REF]. The main goal of those systems is the design of tools for biomedical applications like neuroprosthesis [START_REF] Ambroise | Real-time biomimetic Central Pattern Generators into FPGA for hybrid experiments[END_REF] and for the understanding of the human nervous system [START_REF] Ir | Biomaterials. Electronic dura mater for long-term multimodal neural interfaces[END_REF]. Nevertheless using silicon neurons bring some bio-compatibility issues (reject, power consumption). That's why a new way of neuromorphic engineering should be explored: the artificial neural systems based on microfluidic techniques (Figure 1). This new research does not exist yet in the state of the art.

THEORY

The main goal of this artificial neuron is to simplify the interactions with biological neurons and to reduce the bio-compatibility issues keeping the same ionic currents (K+, etc.) than in biology and using biocompatible material like PDMS. The biological neuron transmits electrical signals based on ion flow through their plasma membrane. Action potentials are propagated along axons and represent the fundamental electrical signals by which information are transmitted from one place to another in the nervous system.

EXPERIMENTAL

Based on this physiological behavior and Hodgkin-Huxley formalism [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF], we propose a microfluidic structure composed of chambers representing the intra and extracellular environments, connected by channels actuated by Quake valves [START_REF] Unger | Monolithic Microfabricated Valves and Pumps by Multilayer Soft Lithography[END_REF]. The Quake valve are activated for closing the channel and then to put the membrane potential to 0V. The activation of these valves is done thank to air pressure injection (Figure 2). These channels are equipped with selective permeable membranes (Nafion) to mimic the exchange of species found in the biological neuron. Thick PDMS membrane is used to create the Quake valve membrane. Integrated gold microelectrodes are used to measure the potential difference between the intracellular and extracellular environments: the membrane potential (Figure 3).

RESULTS AND DISCUSSION

We obtain an electrical behavior similar to the biological neuron (Figure 4). The two action potentials are created by two air pulses of 100ms which close the Quake valve (one at 5.2s and one at 6.3s). The repolarizing phase has slower velocity than the rising phase like we observe in biology. The main parameter for controlling the shape of the action potential is the difference of ion concentration between the two chambers. In our case, we use different concentrations of KCl (1mmol.L-1 and 10µmol.L-1), which involve a 120mV of difference of potential (Nernst equation). We could also control the spike frequency with the programmed time sequence of the quake valve controller.

CONCLUSION

The preliminary results of this microfluidic device demonstrate that it is able to mimic the electrical activity of one biological neuron, in terms of amplitude, shape and frequency.

The next step of our work is to make hybrid experiments in the same chip: neuron culture in one part of the PDMS device and the artificial neurons in the other part. Figure 5 describes this PDMS device. In the left part, gold microelectrode arrays is used for stimulation or recording of biological cells like neurons or muscle cells. In the right part, the microfluidic neuron is used for stimulating the biological part.

Figure 1 :

 1 Figure 1: Microfluidic neuron position and future hybrid experiments between real neurons and microfluidic neurons.

Figure 2 :

 2 Figure 2: On the left picture, the Quake valve is not activated, on the right one, it is activated by air pressure.These channels are equipped with selective permeable membranes (Nafion) to mimic the exchange of species found in the biological neuron. Thick PDMS membrane is used to create the Quake valve membrane. Integrated gold microelectrodes are used to measure the potential difference between the intracellular and extracellular environments: the membrane potential (Figure3).

Figure 3 :

 3 Figure 3: Description of the different layers of the microfluidic neuron. In the bottom left PDMS channel, we use a 1mmol.L -1 KCL and in the upper right PDMS channel, a 10µmol.L -1 .

Figure 4 :

 4 Figure 4: Membrane potential of the microfluidic neuron depending of the activation of Quake valve.

Figure 5 :

 5 Figure 5: PDMS platform for hybrid experiments with gold electrode array for cell culture (left) and microfluidic neuron (right).

ACKNOWLEDGEMENTS

EUJO-LIMMS activities are also supported by the Core-to-Core Program of the Japan Society for the Promotion of Science (JSPS) as a matching fund to the 7th Framework Programme of the European Union.

We would like also to thank Vincent Raimbault from IMS Lab., University of Bordeaux, France and Ayako Araki from IIS, The University of Tokyo, Japan for their support in microfluidic knowledge.