Stéphane Caron

Quang-Cuong Pham

Yoshihiko Nakamura

Stability of surface contacts for humanoid robots: Closed-form formulae of the Contact Wrench Cone for rectangular support areas

come   L'archive ouverte pluridisciplinaire

I. INTRODUCTION

From a mechanical point of view, establishing contact with the environment amounts to constraining a certain number of Degrees Of Freedom (DOF) of an end-effector link. For instance, a sliding planar contact imposes three equality constraints (two on the orientation of the link and one for the link-to-surface distance) while a fixed contact constraints all six DOFs of the end-effector link. A stability criterion can be seen as a set of inequality constraints describing the conditions under which these equalities are preserved.

The Zero-Moment Point (ZMP) [START_REF] Vukobratović | Zero-moment pointthirty five years of its life[END_REF] criterion, namely that this point lies in the convex hull of ground contact points, may be the best-known example of stability criterion. It is known [START_REF] Hirukawa | A universal stability criterion of the foot contact of legged robots-adios zmp[END_REF] to be necessary but not sufficient for contact stability, as it is accounts for neither sliding in the horizontal plane nor rotation around the vertical axis (yaw). Yet, humanoid robots are often subject to significant yaw moments during single support phases, resulting in undesired foot rotations. Attesting the importance of this problem, recent works have been using upper-body motions to compensate for these yaws while continuing to use ZMP [START_REF] Ugurlu | Yaw moment compensation for bipedal robots intrinsic angular momentum constraint[END_REF], [START_REF] Cisneros | Yaw moment compensation by using full body motion[END_REF].

A necessary and sufficient condition for the stability of the six contact DOFs can be obtained by considering individual contact forces distributed on the contact surface, as can be found e.g., in bipedal balance control [START_REF] Ott | Posture and balance control for biped robots based on contact force optimization[END_REF] or motion planning [START_REF] Saab | Dynamic whole-body motion generation under rigid contacts and other unilateral constraints[END_REF], [START_REF] Hauser | Fast interpolation and time-optimization with contact[END_REF]. This approach is however hampered by redundancy: the vector of contact forces has many more components (three times the number of contact points) than the DOF of the contact constraint (six). This redundancy impedes the inversion of the equations of motion, notably the use of stateof-the-art Inverse Dynamics based on QR-decomposition [START_REF] Righetti | Optimal distribution of contact forces with inverse-dynamics control[END_REF] (which do not apply to redundant force variables) 1 1 . In Time-Optimal Path Parameterization (TOPP), redundancy prompted the use of further contact approximations [START_REF] Caron | Kinodynamic motion retiming for humanoid robots[END_REF] or expensive polytope projection algorithms [START_REF] Hauser | Fast interpolation and time-optimization with contact[END_REF]. In the present paper, we argue that such workarounds can be avoided by using a suitable contact representation (for instance, both [START_REF] Caron | Kinodynamic motion retiming for humanoid robots[END_REF] and [START_REF] Hauser | Fast interpolation and time-optimization with contact[END_REF] boil down to a single matrix inversion for a biped in single-support, as we will see in Section V).

The key insight here is: distributed contact forces lie in their respective friction cones if and only if the contact wrench belongs to a certain wrench cone [START_REF] Balkcom | Computing wrench cones for planar rigid body contact tasks[END_REF], which we calculate. The contact wrench naturally solves the redundancy issue, as its dimension is minimal (six). In the context of multi-contact planning, it was advocated as a generalization of ZMP in [START_REF] Hirukawa | A universal stability criterion of the foot contact of legged robots-adios zmp[END_REF] and applied to walking on rough terrains [START_REF] Hirukawa | A pattern generator of humanoid robots walking on a rough terrain[END_REF], [START_REF] Zheng | A walking pattern generator for biped robots on uneven terrains[END_REF]. However, these works made the same "sufficient friction" assumption as ZMP, meaning that the resulting criterion does not account for sliding and yaw rotations. Besides, the contact wrench was computed from individual contact forces, which we argue yields unnecessarily involved calculations: the wrench cone can be computed explicitly from the sole geometry of the end-effectors in contact.

In the present communication, we derive the closed-form formulae of the wrench cone in the case of rectangular support areas, which is of practical importance since most humanoid robot feet can be adequately approximated by rectangles. This result simplifies dynamics computations, as we will see for humanoid motions. It also provides an analytical description of "yaw friction", from which we derive a control law to prevent undesirable yaw rotations.

The rest of the paper is organized as follows. In Section II, we recall the definitions related to contact stability. In Section III, we discuss the physics of surface contact and give a theoretical justification to the model where individual contact forces are applied at the vertices of the contact polygon. Then, in Section IV, we derive a closed-form expression of the wrench cone in the case of rectangular contact areas.

We apply the resulting solution in a humanoid experiment in Section V before concluding in Section VI.

II. BACKGROUND

A. Contact Forces and Contact Wrench

Consider a robot with n degrees of freedom making N point contacts with the environment, at points C 1 , . . . , C N in the laboratory frame. The equations of motion of the robot are:

M(q) q + h(q, q) = S τ a + N i=1 J(C i ) f i , (1)
where q, q, q are the n-dimensional vectors of DOF values, velocities and accelerations, M is the n × n joint-inertia matrix, h(q, q) the n-dimensional vector of gravity and Coriolis forces. In case the robot has n a actuated joints, τ a is the n a -dimensional vector of torques at the actuated joints and S is the n a × n joint selection matrix. Finally, for each i ∈ [1, N ], f i is a 3-dimensional vector of contact force (in the laboratory frame) and J(C i ) is the 3 × n translation Jacobian calculated at point C i . We assume that both the environment and the contacting link are rigid bodies. Thus, interactions between them can be represented by a single contact wrench w = (f , τ ), which one can compute from contact forces as:

f def = i f i , (2)
τ def = i --→ OC i × f i , (3)
where O is the origin of the link frame. The contact jacobian J wr is the 6×n matrix obtained by stacking vertically J(O), the translation Jacobian computed at O, and J rot , the rotation Jacobian of the link, both taken with respect to the absolute frame. With these definitions,

J wr w = N i=1 J(C i ) f i . (4)
(See Appendix A for calculations.) Thus, the equations of motion can be rewritten as:

M(q) q + h(q, q) = S τ a + J wr w. (5)

B. Contact Stability

Assume now that the robot is in a given state (q, q). Then, accelerations q and generalized forces (actuated torques or contact forces) are bound by a complementarity condition [START_REF] Pang | Stability characterizations of rigid body contact problems with coulomb friction[END_REF]. For the sake of explanation, let us describe first the case of a single translation coordinate x, as depicted in Figure 1. Under Coulomb's friction model, either of the following situations occurs :

• Fixed contact: ẍ = 0 and the contact force obeys

|f t c | ≤ µf n c
, where f t c and f n c are the horizontal and vertical components of f c and µ is the static friction coefficient;

f a f c g m x Fig. 1.
Simple case of a body on a horizontal surface, with one DOF, under the action of three forces: its weight (mg), the resultant contact force (fc) and its inertial force (fa = [mẍ 0 0]).

• Sliding contact: ẍ > 0 and the contact force obeys

f t c = -µ k f n c
, where µ k is the kinetic friction coefficient. The acceleration ẍ and contact force f c are in a complementary relationship: when one is equality-constrained, the other is inequality-constrained.

Similarly, when a six-DOF end-effector is in contact, its translational and rotational accelerations are in a complementary relationship with some generalized contact forces γ (here, γ = w or γ = (f 1 , . . . , f N )). The contact mode describes, for each variable in a complementary relationship, which of the acceleration or force variables is equalityor inequality-constrained. In the present work, we are interested in the fixed contact mode where the position and orientation of the end-effector are equality-constrained to a reference value, all force variables being therefore inequalityconstrained.

Definition 1 (Weak Contact Stability): A contact is weakly stable when there exists a solution ( q, τ a , γ) of the equations of motion satisfying the fixed contact mode for all contacting links. That is to say,

• for each contact (i), the relative velocity and acceleration at contact are zero:

J (i)
wr q = 0 and J (i) wr q = -J(i) wr q, • actuated torques τ a are within torque limits, • complementary forces γ satisfy their inequality constraints (friction cones or wrench cone). This formulation has been widely used in the literature. In approaches based on inverse dynamics, the first condition is first enforced kinematically, then torques and complementary forces are solved [START_REF] Hauser | Fast interpolation and time-optimization with contact[END_REF], [START_REF] Righetti | Optimal distribution of contact forces with inverse-dynamics control[END_REF]. The "weakness" in the definition above refers to the notions of strong and weak stability, as stated by [START_REF] Pang | Stability characterizations of rigid body contact problems with coulomb friction[END_REF]. Strong stability happens when all solutions to the equations of motion satisfy the fixed contact mode.

Note that choosing between contact forces and the contact wrench changes the equations of motion (Equations ( 1) and ( 5)) but the underlying contact stability is the same by Equation (4). In the rest of the paper, we will always refer to contact stability in the weak sense.

III. SURFACE CONTACT

Suppose we take contact forces f 1 , . . . , f N as complementary variables to the position and orientation of the contacting link. Let f n i and f t i denote the normal and the tangential components of the contact force f i . Coulomb friction provides the complementary inequalities:

• Unilaterality:

f n i > 0, (6)
• Non-slippage:

f t i ≤ µf n i , (7)
with µ is the static coefficient of friction. However, when the contact is mediated through a surface and not through a set of points, the reality of contact is continuous. To account for this continuity, we model the action of the environment at the surface S by two quantities: a scalar field p(x, y) corresponding to normal pressure, and a two-dimensional vector field σ(x, y) for tangential mechanical stress. Figure 2-(A) illustrates these two fields for a rectangular contact area. For convenience, we also will denote by ν def = σ(x, y) + p(x, y)n, where n is the unit vector normal to the contact surface (pointing upward). The equations of motion become

M q + h = S τ a + S J(C xy ) ν(x, y)dxdy (8)
where J(C xy ) is the 3 × n translation Jacobian calculated at the point of coordinate C xy on the surface. By analogy with its discrete counterpart ( 2)-( 3), the wrench resulting from ν(x, y) is

f def = S ν(x, y)dxdy, (9)
τ def = S ---→ OC xy × ν(x, y)dxdy, (10)
where O is the origin of the link frame. Note how the resulting torque in Equation ( 10) involves infinitesimal forces ν(x, y) but not infinitesimal torques (under bounded forces, torques vanish when application points draw infinitely closer). In other words, our model of surface contacts is a continuum of point contacts. Accordingly, under Coulomb friction, the inequality constraints for ν(x, y) are:

• Unilaterality: p(x, y) > 0, (11)
• Non-slippage: σ(x, y) ≤ µ p(x, y).

In the present literature, surface contact is often modeled using sets of contact points. In the light of a continuous model, the proposition below gives a theoretical justification for this practice.

Proposition 1: Assume that the contact surface S is a convex polygon with vertices C 1 , . . . , C N . If there exists a field (x, y) → ν(x, y) satisfying complementary inequalities (11)- [START_REF] Zheng | A walking pattern generator for biped robots on uneven terrains[END_REF], then there exists contact forces applied at C 1 , . . . , C N , summing up to the same contact wrench, and satisfying complementary inequalities ( 6)- [START_REF] Hauser | Fast interpolation and time-optimization with contact[END_REF].

Proof: consider pressure and stress fields summing up to w. By convexity, one can find barycentric coordinates α 1 (x, y), . . . , α k (x, y), i.e., positive functions such that i α i (x, y) = 1 and each point C xy ∈ S can be written

C xy = i α i (x, y)C i . Then, define for each vertex C i a force f i := S α i (x, y)ν(x, y)dxdy, C 1 C 2 C 3 C 4 O X Y (A) (B)
x y z By positivity of the α i 's, it is straightforward to check that all f n i > 0 and f t i ≤ µf n i . In addition, this expression of f i ensures that the resulting wrenches are equal, i.e., (2) = (9) and (3) = [START_REF] Balkcom | Computing wrench cones for planar rigid body contact tasks[END_REF].

If Dirac fields are authorized, then one can show that the converse implication of Proposition 1 is true. If not, whether this converse implication is true remains an open question. The bottom line of this argument is that the forces at the vertices of the convex hull completely render the dynamics of the surface contact. The wrench cone that we derive in the next section will share the same property.

IV. WRENCH CONE FOR RECTANGULAR SURFACES

While the wrench cone could be computed numerically for arbitrary contact surfaces by using the face-span manipulations presented in [START_REF] Hirai | Analysis and planning of manipulation using the theory of polyhedral convex cones[END_REF], [START_REF] Balkcom | Computing wrench cones for planar rigid body contact tasks[END_REF] 2 , the analytical formula that we derive below gives insights into the nature of that cone while grounding for more efficient downstream calculations.

Consider a rectangular area (C 1 C 2 C 3 C 4 ) as depicted in Figure 2 (B). We calculate the contact wrench w at the center O in the link frame. Let us denote by (f x i , f y i , f z i ) the three components of the contact force f i at point C i , expressed in the link frame. Under the common linear approximation of friction cones, Coulomb inequalities become

|f x i |, |f y i | ≤ µf z i (13) f z i > 0 ( 14
)
The following proposition gives the analytical formulae for the wrench cone. Proposition 2 (Contact Wrench Cone): There exists a solution (f 1 , . . . , f 4 ) satisfying inequalities ( 13)-( 14) if and only if there exists a wrench w = (f x , f y , f z , τ x , τ y , τ z ) such that:

|f x | ≤ µf z (15) |f y | ≤ µf z (16) f z > 0 (17) |τ x | ≤ Y f z (18) |τ y | ≤ Xf z (19) τ z min ≤ τ z ≤ τ z max ( 20
)
2 also known as the double description method where

τ z min def = -µ(X + Y )f z + |Y f x -µτ x | + |Xf y -µτ y |, τ z max def = +µ(X + Y )f z -|Y f x + µτ x | -|Xf y + µτ y |.
Provided that the kinematic contact equation is satisfied, this condition is necessary and sufficient for weak contact stability.

The proof of this proposition is given in Appendix B. The validity of the formula was also tested empirically with a script available at [15].

Let us now detail each line of the Contact Wrench Cone (CWC). The first two inequalities (15)-( 16) correspond to the usual Coulomb friction. Inequalities ( 17), ( 18) and ( 19) are equivalent to having the CoP lie inside of the contact area. The last inequality (20) provides a bound on the admissible yaw torque that was implicitely encoded in the contact-force model.

Observe how this relation is more complex than "no rotation occurs while τ z is small enough", as it is coupled with all other components of the contact wrench. The yaw torque is bounded by τ z min and τ z max , both of which may be either positive or negative (for instance, both become negative when the CoP nears the corner τ x = Y f z , τ y = Xf z of the support polygon). Notably, the "safest" value for τ z is not zero but: We see from this rewriting that higher tangential forces or roll-pitch torques reduce the range of admissible yaw torques. In particular, when these other constraints are saturated (e.g., when the CoP reaches a corner of the support polygon), τ z safe is the only solution that prevents the contact from breaking. Therefore, τ z = τ z safe appears as a sensible control law to prevent undesired yaw rotations.

τ z safe def = (τ z min + τ z max )/2 = sgn(-f x τ x ) min(Y |f x |, µ|τ x |) + sgn(-f y

V. EXPERIMENT

We test the CWC in the dynamics simulator OpenHRP with a model of the HRP4 robot. Note that OpenHRP has its own contact model where forces are distributed along the edges (not corners) of the contact surface.

We enforce the CWC within the Time-Optimal Path Parameterization framework (TOPP), a well-known projection of system dynamics along a pre-defined path that has been used for motion planning of humanoid robots [START_REF] Hauser | Fast interpolation and time-optimization with contact[END_REF], [START_REF] Pham | Kinodynamic planning in the configuration space via admissible velocity propagation[END_REF], [START_REF] Caron | Kinodynamic motion retiming for humanoid robots[END_REF]. We considered the case of a single contact at the left foot and designed a motion that would challenge all six contact DOFs. In single contact, the contact wrench is fully determined by the unactuated rows of the equation of motion, i.e., w = (PJ wr ) -1 P [M(q) q + h(q, q)] where P = I -S is the 6 × n projector on the free-flying coordinates. The matrix PJ wr is non-singular thanks to the use of the contact wrench.

For this experiment, we designed by hand a set of eleven key postures. The geometric path was obtained by interpolating Bezier curves between these postures and using inverse kinematics to fix the position and orientation of the support foot on the ground. To time this path into a trajectory, we used the open-source TOPP solver [START_REF] Pham | A general, fast, and robust implementation of the time-optimal path parameterization algorithm[END_REF]. The solver takes as input the path and a vector representation of the system dynamics along it, which is easy to compute once one knows the projectors mentioned above (see [START_REF] Caron | Kinodynamic motion retiming for humanoid robots[END_REF] for details). As is common with numerical TOPP solvers, we added safety margins in the CWC computations input to TOPP: the foot dimensions X and Y were scaled by 70% (versus 100% in OpenHRP) and the friction coefficient set to µ = 0.4 (versus µ = 0.8 in OpenHRP).

Figure 3 shows a timelapse of the final retimed motion. A video is also available online at [15]. Because time optimality yields bang-bang control laws, the retimed motion always saturates at least one of the contact constraints. Therefore, trying to execute it faster should result in the end-effector breaking surface contact, e.g., tipping on the edge of the foot. We observed this phenomenon in the experiment, as depicted in Figure 4. Note that, as we used HRP4's stabilizer while executing the motion, it was still possible in practice to accelerate the motion by about 5% (i.e., reducing the total duration by 5% with uniform time scaling) whithout noticing any change in contact stability.

VI. CONCLUSION

In this paper, we calculated the closed-form expression of the Contact Wrench Cone for rectangular contact surfaces. This formula has several implications. First, it is very simple, making computations easier than other models based on contact forces. Second, it describes the phenomenon of yaw friction by a double inequality that is, to the best of our knowledge, novel. From these bounds, we derived a simple control law to avoid undesired foot rotations, a recurring problem for bipeds in single contact. We also assessed the validity of the CWC condition in dynamics simulations.

The mathematical derivation of the wrench cone formula (Appendix B) can be carried out for parametric shapes other than rectangles. We are surveying this application as part of future work, as well as the combination of local wrench cones into the gravito-inertial wrench [START_REF] Hirukawa | A universal stability criterion of the foot contact of legged robots-adios zmp[END_REF].

Fig. 2 .

 2 Fig. 2. Contact in the surface plane. (A) Illustration of the stress field (blue arrows) and pressure field (red discs representing magnitude). (B) Notations used in Section IV.

  τ y ) min(X|f y |, µ|τ y |), with sgn the sign function. From Equation (20), τ z may deviate from τ z safe by at most µ(X + Y )f z -max(Y |f x |, µ|τ x |) -max(X|f y |, µ|τ y |).

Fig. 3 .

 3 Fig. 3. Snapshots of the retimed motion. Total duration is 7.3 s. Time stamps are shown below each frame. The motion stresses all components of the wrench cone. The first segment stresses the pitch by moving the CoP forward. The second stresses the rool through arm motions. Meanwhile, the yaw component is stressed by chest-pitch motions. Finally, the waist performs an elliptic motion (back and forth, up and down) throughout the whole motion, thus stressing the translation of the contact foot.

Fig. 4 .

 4 Fig. 4. Zoom on the contact forces computed by OpenHRP at the left foot. Each thin purple segment corresponds to one force. Time stamps are shown below each column. First row: retimed motion with CWC enforced. Second row: 10% acceleration of the motion, CWC is not enforced. Around 4.2 s, the contact surface degenerates to a line contact on the left edge of the foot. Third row: 15% acceleration of the motion, CWC is not enforced. The contact surface degenerates to a line before contact is lost and the humanoid falls.

Department of Mechano-Informatics, The University of Tokyo, Japan, caron@ynl.t.u-tokyo.ac.jp

School of Mechanical and Aerospace Engineering, Nanyang Technological University, Singapore

When this is not the case, the authors advocate the use of Singular Value Decomposition to compute new independent variables; however, they do not mention how to compute the inequality constraints applying to these new variables.

Acknowledgments

This work was partially supported by Tier 1 grant n • RG109/14 awarded by Singapore MOE to the second author.

APPENDIX

A. Proof of Equation 4The rotation Jacobian satisfies ω = J rot q, where ω is the rotation velocity of the link. As a consequence of this, for any point C on the link and any vector u, one has

Next, the position of any point C i of the link in the absolute frame is related to that of its origin O by

The corresponding translation Jacobian is

The first term of the expression equals the translation component of J wr w. By applying (21), we see that the second term equals J rot ( --→ OC i × f i ). Factoring J rot out of the summation yields the rotation component of J wr w.

B. Calculation of the Wrench Cone

The wrench is defined by ( 2)-(3) as:

. By unilaterality [START_REF] Hirai | Analysis and planning of manipulation using the theory of polyhedral convex cones[END_REF] we have f z > 0, so we can define:

and normalize the system by dividing each row accordingly.

From the non-slippage constraint (13), we have α x i , α y i and

Then, introduce the new variables:

We can reduce the complete system in two ways. First, using the fact that the relation M from α to γ is a linear surjection from

the computation of antecedents being straightforward). Since there is no other constraint on the α i 's than their domain and relation to γ i 's, one can obtain an equivalent system by replacing α ∈ [-1, 1] 8 by γ ∈ M[-1, 1] 8 . Then, the three equations in D i 's are:

By linear combination, we can use them to rewrite M[-1, 1] 8 as:

Finally, using the same equations, one can express all D i 's as functions of, e.g., D 4 . The inequalities D ∈ [0, 1] 4 become:

This system has solutions if and only if all lower bounds are smaller than all upper bounds. Matching all pairs of lower and upper bounds one by one (we will show an example of this technique below in a more complex situation), one can check that all these inequalities boil down to

]. The complete system is now:

One can use the first three equations to eliminate the redundant variables γ x , γ y and γ y , expressing them as functions of γ x in the inequality constraints. After simplification, the resulting system is:

There exist a solution γ x if and only if all of its lower bounds are smaller than all of its upper bounds. Let us match all pairs of lower bounds (26)-( 29) and upper bounds ( 22)-(25). One can check that:

Consequently, the complete system becomes:

And (K 1 , K 2 , C 1 , C 2 ) ∈ [-1, 1] 4 . In a more concise form, these last eight inequalities can be written

We conclude by de-normalizing all inequalities.


Stéphane Caron

Quang-Cuong Pham

Yoshihiko Nakamura

ZMP support areas for multi-contact mobility under frictional constraints

Keywords: Contact stability, Humanoid locomotion, Zerotilting Moment Point (ZMP)

We propose a method for checking and enforcing multi-contact stability based on the Zero-tilting Moment Point (ZMP). The key to our development is the generalization of ZMP support areas to take into account (a) frictional constraints and (b) multiple non-coplanar contacts. We introduce and investigate two kinds of ZMP support areas. First, we characterize and provide a fast geometric construction for the support area generated by valid contact forces, with no other constraint on the robot motion. We call this set the full support area. Next, we consider the control of humanoid robots using the Linear Pendulum Mode (LPM). We observe that the constraints stemming from the LPM induce a shrinking of the support area, even for walking on horizontal floors . We propose an algorithm to compute the new area, which we call pendular support area. We show that, in the LPM, having the ZMP in the pendular support area is a necessary and sufficient condition for contact stability. Based on these developments, we implement a whole-body controller and generate feasible multi-contact motions where an HRP-4 humanoid locomotes in challenging multi-contact scenarios.

I. INTRODUCTION

The Zero-tilting Moment Point (ZMP) is the dynamic quantity thanks to which roboticists solved the problem of walking on horizontal floors. One of its key properties is that dynamic stability, i.e. the balance of gravity and inertial forces by valid contact forces, implies that the ZMP lies in the convex hull of ground contact points, the so-called support area [START_REF] Vukobratović | On the stability of anthropomorphic systems[END_REF], [START_REF] Goswami | Postural stability of biped robots and the foot-rotation indicator (fri) point[END_REF]. The support area thus provides a necessary (non-sufficient) condition for contact stability on horizontal floors.

For locomotion, the second key property of the ZMP lies in its coupling with the position of the center of mass (COM). By keeping a constant angular momentum and constraining the COM to lie on a plane, this relation simplifies into the Linear Inverted Pendulum Mode (LIPM) [START_REF] Kajita | Study of dynamic biped locomotion on rugged terrain-derivation and application of the linear inverted pendulum mode[END_REF], [START_REF] Kajita | The 3d linear inverted pendulum mode: A simple modeling for a biped walking pattern generation[END_REF]. In the LIPM, the COM moves away from the ZMP under the linear dynamics of a point-mass at the tip of an inverted pendulum. The stabilization problem is then to control the position of the tip (COM) of the pendulum by moving its fulcrum (ZMP).

These two main merits, a geometric stability condition and linearized dynamics, are as well-known as the two main limitations of the classical ZMP: it does not account for friction, and it can only be applied when all contacts are Stéphane Caron is with the Laboratoire d'Informatique, de Robotique et de Microélectronique de Montpellier (LIRMM), CNRS / Université de Montpellier. Quang-Cuong Pham is with the School of Mechanical and Aerospace Engineering, Nanyang Technological University, Singapore. Yoshihiko Nakamura is with the Department of Mechano-Informatics, University of Tokyo, Japan.

Corresponding author: stephane.caron@normalesup.org. coplanar. The latter results from the definition of the ZMP as the point on the floor. In a general multi-contact scenario, each contact defines its own surface and there is no single "floor" plane. In a classic survey paper [START_REF] Sardain | Forces acting on a biped robot. center of pressure-zero moment point[END_REF], Sardain and Bessonnet stated the problem to address as follows:

The generalization of the ZMP concept [to the case of multiple non-coplanar contacts] would be actually complete if we could define what is the pseudosupport-polygon, a certain projection of the threedimensional (3-D) convex hull (built from the two real support areas) onto the virtual surface, inside which the pseudo-ZMP stays. In this paper, we construct the areas conjectured by Sardain and Bessonnet. Our first contribution is to characterize the ZMP support area generated by valid contact forces, with no other constraint on the robot motion. We call this set the full support area. Our analysis provides a geometric construction that allows for fast calculations. Also, contrary to the assumption that "friction limits are not violated" usually made in the literature, this area fully takes friction into account.

From a control point of view, locomoting systems usually regulate their linear and angular momenta, as examplified by the Linear Pendulum Mode where the robot maintains a constant COM height and moves with a constant angular momentum. These tasks limit the whole-body momentum of the system, which consequently shrinks the ZMP support area. The second contribution of this paper is an algorithm to compute the ZMP support area for the linear-pendulum mode, which we call pendular support area. The latter takes into account both friction and angular-momentum constraints on the robot motion.

Combining these two advances, we design a whole-body controller for humanoids locomoting on arbitrary terrains. Choosing a virtual plane above the COM, we regulate the robot dynamics around that of a linear non-inverted pendulum. We showcase the applicability of the controller by locomoting a model of the HRP-4 humanoid robot in a challenging multicontact scenario involving a large step supported by a handon-wall contact.

The paper is organized as follows. In Section II, we review previous work and introduce background notions necessary to understand the present work. In Section III, we characterize the full support area. We then provide in Section IV an algorithm to calculate the pendular support area, which we apply to multi-contact locomotion in Section V. Concluding remarks are finally provided in Section VI.

II. BACKGROUND

A. Previous work 1) Stability criteria: on horizontal floors, the ZMP of dynamically stable motions lies within the convex hull of contact points (CHCP). However, when the robot makes contact with different non-coplanar surfaces, the ZMP can no longer be defined as a point on the "ground" and the CHCP has no established connection with dynamic stability. Various attempts have been made in the literature to overcome this difficulty.

One line of research [START_REF] Sardain | Forces acting on a biped robot. center of pressure-zero moment point[END_REF], [START_REF] Mitobe | Control of legged robots during the multi support phase based on the locally defined zmp[END_REF], [START_REF] Harada | Dynamics and balance of a humanoid robot during manipulation tasks[END_REF], [START_REF] Inomata | 3dzmp-based control of a humanoid robot with reaction forces at 3-dimensional contact points[END_REF] conjectured that the CHCP (a 3D volume in general) conveys the stability condition, and consequently sought to define a new point lying within this volume. Both [START_REF] Mitobe | Control of legged robots during the multi support phase based on the locally defined zmp[END_REF] and [START_REF] Inomata | 3dzmp-based control of a humanoid robot with reaction forces at 3-dimensional contact points[END_REF] assumed that the moments at centers of pressure and ZMP are all zeros, which is not the case in general 1 and thus results in a point that may not exist even in situations where stability is possible. Harada et al. [START_REF] Harada | Dynamics and balance of a humanoid robot during manipulation tasks[END_REF] considered the CHCP as a ZMP support volume when the robot makes two feet contact with a horizontal floor and hand contacts with the environment. They detailed how to project the support volume on the floor to obtain a ZMP support area. While their construction applies to the general case with nonzero angular momentum, like all approaches based on convex hull of contact points, it assumes infinite friction coefficients. In this paper, we will see how to construct support areas that also take friction into account.

A parallel line of research kept the ZMP in a plane but relaxed the constraint that it coincides with the floor [START_REF] Kagami | A fast dynamically equilibrated walking trajectory generation method of humanoid robot[END_REF], [START_REF] Sugihara | Real-time humanoid motion generation through zmp manipulation based on inverted pendulum control[END_REF], [START_REF] Shibuya | Trajectory planning of biped robot using linear pendulum mode for double support phase[END_REF], [START_REF] Sato | Walking trajectory planning on stairs using virtual slope for biped robots[END_REF]. Kagami et al. [START_REF] Kagami | A fast dynamically equilibrated walking trajectory generation method of humanoid robot[END_REF] had the insight that the ZMP could be taken relative to any plane normal, yet still assuming that this plane should pass through all contact points, which restricted their scope to a maximum of three contacts points. Sugihara et al. [START_REF] Sugihara | Real-time humanoid motion generation through zmp manipulation based on inverted pendulum control[END_REF] introduced the notion of "Virtual Horizontal Plane" (VHP) in which contact points are projected on a virtual plane via the line connecting them to the COM, and the convex hull of these points is then taken as support area. Shibuya et al. [START_REF] Shibuya | Trajectory planning of biped robot using linear pendulum mode for double support phase[END_REF] took the idea further by considering a virtual plane above the COM, while Sato et al. [START_REF] Sato | Walking trajectory planning on stairs using virtual slope for biped robots[END_REF] applied this idea to stair climbing. However, like CHCP, VHP support areas suppose infinite friction coefficients. The approach proposed in the present paper also uses a virtual plane and the linear pendulum mode, but our support areas fully account for friction, which makes them a necessary and sufficient condition for contact stability.

Breaking away from the notion of ZMP, another line of work has focused on building criteria that keep equivalence with full contact stability [START_REF] Saida | Fsw (feasible solution of wrench) for multi-legged robots[END_REF], [START_REF] Hirukawa | A universal stability criterion of the foot contact of legged robots-adios zmp[END_REF], [START_REF] Qiu | Human motions analysis and simulation based on a general criterion of stability[END_REF], [START_REF] Escande | Planning contact points for humanoid robots[END_REF], [START_REF] Caron | Leveraging cone double description for multi-contact stability of humanoids with applications to statics and dynamics[END_REF]. The seminal work by Saida et al. [START_REF] Saida | Fsw (feasible solution of wrench) for multi-legged robots[END_REF] impulsed a shift of paradigm from the ZMP to the gravito-inertial wrench. It also proposed to orient the virtual plane orthogonally to the resultant force, which reinstates the support area as convex hull of (projected) contact points -an idea that may have been overlooked by the literature so far. Next, Hirukawa et al. [START_REF] Hirukawa | A universal stability criterion of the foot contact of legged robots-adios zmp[END_REF] constructed the first full stability criterion for the gravito-inertial wrench, yet with a high number of variables including all contact forces. Later developments [START_REF] Qiu | Human motions analysis and simulation based on a general criterion of stability[END_REF], [START_REF] Escande | Planning contact points for humanoid robots[END_REF], [START_REF] Caron | Leveraging cone double description for multi-contact stability of humanoids with applications to statics and dynamics[END_REF] reduced these redundant variables to the gravito-inertial wrench using the double-description method. Compared to traditional ZMP solutions, this approach has the benefit of providing a full stability criterion, but at the cost of the non-linear dynamics of the gravito-inertial wrench. Furthermore, the nice geometric construction of the ZMP area is replaced by a considerably less intuitive six-dimensional cone. The method that we introduce in this paper reconciles full stability, linearized dynamics and a geometric support area.

2) Control: when it comes to control, the use of the ZMP is historically tied to the LIPM, which was introduced in [START_REF] Kajita | The 3d linear inverted pendulum mode: A simple modeling for a biped walking pattern generation[END_REF], [START_REF] Sugihara | Real-time humanoid motion generation through zmp manipulation based on inverted pendulum control[END_REF] and used in a wealth of subsequent works [START_REF] Mitobe | Control of legged robots during the multi support phase based on the locally defined zmp[END_REF], [START_REF] Harada | Dynamics and balance of a humanoid robot during manipulation tasks[END_REF], [START_REF] Shibuya | Trajectory planning of biped robot using linear pendulum mode for double support phase[END_REF], [START_REF] Sato | Walking trajectory planning on stairs using virtual slope for biped robots[END_REF], [START_REF] Kajita | Biped walking pattern generation by using preview control of zero-moment point[END_REF], [START_REF] Harada | An analytical method for real-time gait planning for humanoid robots[END_REF], [START_REF] Morisawa | Biped locomotion control for uneven terrain with narrow support region[END_REF], [START_REF] Tedrake | A closed-form solution for real-time zmp gait generation and feedback stabilization[END_REF]. Kajita et al. [START_REF] Kajita | Biped walking pattern generation by using preview control of zero-moment point[END_REF] brought in the technique of model predictive control as a way to generate COM trajectories from desired ZMP positions. Harada et al. [START_REF] Harada | An analytical method for real-time gait planning for humanoid robots[END_REF] proposed an analytical alternative with polynomial solutions for the coupled COM-ZMP trajectories. Recently, Tedrake et al. [START_REF] Tedrake | A closed-form solution for real-time zmp gait generation and feedback stabilization[END_REF] exhibited a closed-form solution for the linear-quadratic regulator tracking a reference ZMP. However, all of these methods only apply to locomotion on horizontal floors.

Aiming for locomotion on rough terrains, Zhao et al. [START_REF] Zhao | A three dimensional foot placement planner for locomotion in very rough terrains[END_REF] extended the LIPM to a "Prismatic Inverted Pendulum" where the COM altitude is allowed to vary linearly, while Morisawa et al. [START_REF] Morisawa | Pattern generation of biped walking constrained on parametric surface[END_REF] provided a wider derivation where the COM belongs to a general two-dimensional manifold. Rather than the ZMP, recent papers [START_REF] Morisawa | Biped locomotion control for uneven terrain with narrow support region[END_REF], [START_REF] Englsberger | Three-dimensional bipedal walking control based on divergent component of motion[END_REF] chose to control the Capture Point to stabilize the unstable dynamics of the LIPM. In terms of support areas, though, the question is the same for the Capture Point and the ZMP and was not addressed by these developments. With the method proposed in the present paper, we realize a control system with marginally stable dynamics using the ZMP as control point and a linear pendulum mode. The benefit of our approach compared to these previous works is that we are able to derive at the same time the support area corresponding to our control variable.

As with stability criteria, solutions breaking away from control points were also explored in the whole-body control literature. The main alternative is to regulate contact forces directly, resulting in force distribution schemes where desired contact forces and torques are tracked by a whole-body controller [START_REF] Hirukawa | A universal stability criterion of the foot contact of legged robots-adios zmp[END_REF], [START_REF] Hyon | Full-body compliant humanhumanoid interaction: balancing in the presence of unknown external forces[END_REF], [START_REF] Lee | Ground reaction force control at each foot: A momentum-based humanoid balance controller for non-level and non-stationary ground[END_REF], [START_REF] Ott | Posture and balance control for biped robots based on contact force optimization[END_REF], [START_REF] Righetti | Optimal distribution of contact forces with inverse-dynamics control[END_REF]. Force objectives can express whole-body tasks, such as tracking of desired COM or angular momentum, as well as local ones, such as minimizing friction forces [START_REF] Ott | Posture and balance control for biped robots based on contact force optimization[END_REF] or end-effector torques [START_REF] Lee | Ground reaction force control at each foot: A momentum-based humanoid balance controller for non-level and non-stationary ground[END_REF]. Notably, Righetti et al. [START_REF] Righetti | Optimal distribution of contact forces with inverse-dynamics control[END_REF] characterized the class of force-distribution controllers for linear-quadratic objectives in the absence of inequality constraints. Overall, force distribution schemes yield fast computations and can cope with arbitrary contact conditions, but they lack the foresight and intuition of methods based on control points and support areas. Indeed, for locomotion, support areas provide both reachable COM locations and a stability margin (the point-to-boundary distance). Finding such indicators in the high-dimensional contact-force space is still elusive. In recent developments, [START_REF] Nagasaka | Whole-body control of a humanoid robot based on generalized inverse dynamics and multicontact stabilizer that can take acount of contact constraints[END_REF], [START_REF] Audren | Model preview control in multi-contact motion-application to a humanoid robot[END_REF] added a level of foresight to their contact-force controllers via model-predictive control, while Zheng et al. [START_REF] Zheng | Generalized distance between compact convex sets: Algorithms and applications[END_REF] constructed a metric that can be used as wrench-space stability margin. In the present paper, we show that support areas can be derived in arbitrary multi-contact configurations as well, providing both COM reachability and stability margins suitable for locomotion.

B. Newton-Euler equations

Let m and G represent the total mass and center-of-mass (COM) of the robot, respectively. We write p A the vector of absolute coordinates of a point A and denote by O the origin of the absolute frame (so that p O = 0). For a link k, define:

• m k the total mass of the link; • p G k the vector of absolute coordinates of its COM G k ; • R k its orientation matrix in the absolute frame;

• ω k its angular velocity in the link frame;

• I k its inertia matrix in the link frame. The linear momentum P and angular momentum L G of the robot, taken at the COM G, are defined by:

P def = link k m k ṗG k , (1)
L G def = link k m k --→ GG k × ṗG k + R k I k ω k . (2)
The fundamental principle of dynamics states that the rate of change of the momentum is equal to the total wrench of forces acting on the system, that is:

Ṗ LG = f g 0 + contact i f i --→ GC i × f i , (3)
where f g denotes the gravity force, C i the i th contact point and f i the contact force exerted by the environment on the robot at C i . Equation ( 3) is called the Newton-Euler equations of the system, sometimes also referred to as "dynamic balance" or the "dynamic equilibrium". It can be equivalently derived from Gauss's principle of least constraint, and corresponds to the six unactuated components in the equations of motion of the system {robot + environment} [START_REF] Wieber | Holonomy and nonholonomy in the dynamics of articulated motion[END_REF]. Define the gravito-inertial wrench, taken this time at O:

w gi O def = f gi τ gi O def = f g - Ṗ p G × (f g -Ṗ) -LG . (4)
Define similarly the contact wrench:

w c O def = f c τ c O def = contact i f i p Ci × f i (5)
The Newton-Euler equation ( 3) can be written in terms of these two wrenches as:

w gi + w c = 0. (6)
This formulation separates spatial accelerations (w gi ), that are usually measured by inertial measurement units, from interaction forces with the environment (w c ), usually measured by force sensors. Since the two wrenches are simply opposites, we will only use the contact wrench (f , τ O ) in the following calculations, dropping the superscript c to alleviate notations.

C. Contact stability

We assume that all contacts between the environment and the robot are surface contacts, i.e. contacts between a flat surface of the robot and a flat surface of the environment. Wrenches exerted on each contacting link are then fully described by applying contact forces at the vertices of the contact polygon [START_REF] Caron | Stability of surface contacts for humanoid robots: Closed-form formulae of the contact wrench cone for rectangular support areas[END_REF], which warrants the formulation by contact points that we have followed so far. Define the contact normal n i at C i as the normal to the contact surface pointing from the environment towards the contacting link. Under Coulomb's model of dry friction, contact forces f i lie inside a friction cone directed by n i :

n i × f c i × n i 2 ≤ µ i (f c i • n i ). (7)
Frictional constraints restrict the range of contact wrenches w c that the robot can generate without breaking any contact: when each contact force lies in a cone C i , the contact wrench lies in the Contact Wrench Cone (CWC) C c projected from all C i 's via the mapping [START_REF] Sardain | Forces acting on a biped robot. center of pressure-zero moment point[END_REF]. Because of the connection ( 6), the gravito-inertial wrench belongs to Gravito-inertial Wrench Cone (GIWC) C gi = -C c . This link between feasible motions and keeping contacts is embodied by the notion of contact stability:

Definition 1 (Weak contact stability [START_REF] Pang | Stability characterizations of rigid body contact problems with coulomb friction[END_REF], [START_REF] Balkcom | Computing wrench cones for planar rigid body contact tasks[END_REF]). A motion of the robot is (weak-contact) stable if and only if the contact wrench it generates belongs to the CWC.

Weak contact stability is the underlying stability criterion used in recent multi-contact developments [START_REF] Hirukawa | A universal stability criterion of the foot contact of legged robots-adios zmp[END_REF], [START_REF] Qiu | Human motions analysis and simulation based on a general criterion of stability[END_REF], [START_REF] Escande | Planning contact points for humanoid robots[END_REF], [START_REF] Caron | Leveraging cone double description for multi-contact stability of humanoids with applications to statics and dynamics[END_REF]. Static stability [START_REF] Bretl | Testing static equilibrium for legged robots[END_REF] corresponds to contact stability when the whole-body momentum is zero. Note that the term "stability" is used here in the sense defined by Pang and Trinkle [START_REF] Pang | Stability characterizations of rigid body contact problems with coulomb friction[END_REF]. It should not be confused with the (a-priori unrelated) notion of Lyapunov stability.

D. Linearized wrench cones

We linearize the quadratic constraint ( 7) by approximating friction cones by friction pyramids . Denoting by f ij the ray vectors of the latter, the constraint becomes:

f i = ray j λ ij f ij , λ ij ≥ 0. (8)
The set of ray vectors {f ij } is known as the span representation of the pyramidal cone. It can be computed directly from the contact frame and friction coefficient µ i . For example, the expression of a four-sided pyramid is

n i ± µi √ 2 t i ± µi √ 2 b i , with (t i , b i , n i )
an orthonormal contact frame, where the normalization by √ 2 is added to make the friction pyramid an inner approximation of the friction cone. Injecting [START_REF] Inomata | 3dzmp-based control of a humanoid robot with reaction forces at 3-dimensional contact points[END_REF] into Equation ( 5) yields a span representation for the contact wrench cone:

f τ O = i,j λ ij f ij p Ci × f ij λ ij ≥ 0. (9)
Let us define τ O,ij = p Ci ×u ij . After re-indexing the couples i, j into a single index i (counting the same contact point C i multiple times accordingly), we get:

f τ O = i λ i f i τ O,i , λ i ≥ 0. (10)
We have thus obtained the span representation of the CWC. A motion is weak-contact stable if and only if its contact wrench can be written as [START_REF] Sugihara | Real-time humanoid motion generation through zmp manipulation based on inverted pendulum control[END_REF] for a certain set of coefficients λ i ≥ 0.

III. FULL ZMP SUPPORT AREA

Let n be a fixed unit space vector, not necessarily aligned with the gravity vector. The Zero-tilting Moment Point (ZMP) is a point Z where the moment of the contact wrench aligns with n [START_REF] Sardain | Forces acting on a biped robot. center of pressure-zero moment point[END_REF], that is, n × τ Z = 0. Consequently,

-n × (p Z × f ) + n × τ O = 0 (11) -(n • f )p Z + (n • p Z )f + n × τ O = 0 (12)
We are interested in computing the ZMP in the plane that contains O and that is orthogonal to n, hereafter denoted by

Π(O, n). The relation Z ∈ Π(O, n) is expressed by n • p Z = 0
, so that the equation above becomes:

p Z = n × τ O n • f (13)
On horizontal floors, O is taken on the floor and n is upward vertical, so that Π(O, n) coincides with the floor plane.

However, in what follows we assume that O can be located at an arbitrary fixed position in space, while n can be an arbitrary unit vector.

A. Construction of the full support area

Equation ( 13) presents the ZMP as a two-dimensional projection of the contact wrench (f , τ O ). Since contact stability is characterized by the CWC, we define the support area of the ZMP as the image of the CWC by this projection: Definition 2 (Full support area). The full support area S of the ZMP in the plane Π(O, n) is the image of the CWC by the projection (13):

S = p Z = n × τ O n • f (f , τ O ) ∈ C c (14)
The key idea to calculate this area is to use the span representation [START_REF] Sugihara | Real-time humanoid motion generation through zmp manipulation based on inverted pendulum control[END_REF] of the CWC, which enables rewriting Equation [START_REF] Saida | Fsw (feasible solution of wrench) for multi-legged robots[END_REF] as

p Z = i λ i (n × τ O,i ) i λ i (n • f i ) , λ i ≥ 0. (15)
Next, define the points Z i by:

p Zi def = n × τ O,i n • f i (16)
and denote by p i def = (n•f i ) the virtual pressure of the contact force generator f i through the virtual plane. Then,

p Z = i λ i p i p Zi i λ i p i , λ i ≥ 0. ( 17
)
On horizontal floors, n and all contact forces f i point upwards, so that all pressures p i are positive. This makes the ZMP Z a convex combination of the Z i 's from Equation [START_REF] Caron | Leveraging cone double description for multi-contact stability of humanoids with applications to statics and dynamics[END_REF]. Furthermore, Equation ( 16) simplifies to Z i = C i , i.e. , the vertices of the support area S coincide with contact points. Our definition of the full support area therefore coincides with the conventional support area on horizontal floors.

In general, however, virtual pressures p i can be either positive or negative. 2

α i def = +λ i p i σ(I + ) for i ∈ I + , α def = σ(I + ) σ(I) . ( 18
)
β i def = -λ i p i σ(I -) for i ∈ I -, β def = σ(I -) σ(I) . (19)
Equation [START_REF] Qiu | Human motions analysis and simulation based on a general criterion of stability[END_REF] becomes

p Z = 1 α -β α i∈I + α i p Zi -β i∈I - β i p Zi . ( 20
)
Define the positive-pressure polygon as the convex hull of Z i 's for i ∈ I + :

P + def = { i∈I + α i p Zi , α i ≥ 0, i α i = 1}
, and define the negative-pressure polygon P -similarly. Equation ( 20) can be rewritten as

p Z = αp Z + -βp Z - α -β , (21)
where α ≥ 0, β ≥ 0, Z + ∈ P + and Z -∈ P -. Next, let

P + -P -= p Z + -p Z -| Z + ∈ P + , Z -∈ P - (22)
denote the Minkowski difference of P + and P -. It is a polygon as difference of two polygons, so that it admits a span representation P + -P -= CONV({r 1 , . . . , r k }) as convex hull (CONV) of a set of vertices. We can now characterize the support area as follows:

Proposition 1. If all virtual pressures p i have the same sign, then the full support area S is the convex hull of the vertices {Z i }. Otherwise, when both P + and P -are nonempty, let

P + -P -= CONV({r 1 , . . . , r k }) (23)
denote their Minkowski difference. The full support area S is then the reunion of two polygonal cones C + and C -given by

C + = P + + i R + r i , (24)
C -= P -+ i R + (-r i ). ( 25
)
In particular, when P + and P -intersect with nonempty interior, the full support area S spans the whole virtual plane Π(O, n).

Proof: (Note that this proof is quite formal and may be skipped by first-time readers.) The main idea of the proof is to combine pairs of points in P + × P -to generate rays in C + and C -. Let us first clear out simple cases. When all p i 's are positive (P -= ∅), Equation [START_REF] Caron | Leveraging cone double description for multi-contact stability of humanoids with applications to statics and dynamics[END_REF] shows that p Z is a convex combination of the vertices of P + , where the weights (λ i p i )/( j λ j p j ) can be chosen freely, so that S = P + . The case where P + = ∅ is treated identically. Next, suppose that both P + and P -are nonempty. Equation ( 21) can be reformulated as

p Z = p Z + + β α -β ----→ Z -Z + = p Z -+ α β -α ----→ Z + Z -(26)
where α (resp. β) is the weight of P + (resp. P -) in Equation [START_REF] Tedrake | A closed-form solution for real-time zmp gait generation and feedback stabilization[END_REF]. Therefore, the set of points Z defined by this equation is

S def = p Z + + β α -β ----→ Z -Z + , α ≥ β ≥ 0, Z ± ∈ P ± ∪ p Z -+ α β -α ----→ Z + Z -, β ≥ α ≥ 0, Z ± ∈ P ± ( 27
)
left foo t rig ht foo t wall contact Fig. 3. Example where the full support area S is the union of two cones C + and C -(in green). There are three contacts: two feet on a horizontal floor, and a wall contact located 50 cm forward and 90 cm above ground. The virtual plane is taken in the feet's plane. Polygons P + and P -used in the geometric construction of the support area are drawn in purple, while the rays of friction cones are depicted by gray lines. In this example, C + shows how the conventional support polygon between the two feet is extended by the wall contact, while C -is a complementary support area resulting from the ability to "push down on the wall" .

Given the orderings of α and β, we can further simplify the ratios into a single positive scalar, so that S = C + ∪ C -with

C + = p Z + + λ ----→ Z -Z + , λ ≥ 0, Z ± ∈ P ± , (28)
C -= p Z -+ λ ----→ Z + Z -, λ ≥ 0, Z ± ∈ P ± . ( 29
)
The set D = P + -P -is a convex polygon as Minkowski difference of two convex polygons. To conclude, we show that C + = P + + R + D. The inclusion ⊂ is straightforward from [START_REF] Righetti | Optimal distribution of contact forces with inverse-dynamics control[END_REF]. Now, let

p C = p Z + 0 + µ(p Z + 1 -p Z -) (30)
denote any point in

P + + R + D. Define p Z + def = 1 1 + µ p Z + 0 + µ 1 + µ p Z + 1 . (31)
One can check that p C = p Z + + µ(p Z + -p Z -), where Z + belongs to P + as convex combination of two points from this convex polygon. Thus C ∈ C + , which establishes the converse inclusion ⊃. Finally, note how, when P + ∩P -has non-empty interior, D contains a neighborhood of the origin. For any pair of points (A, Z + ) ∈ Π(O, n) × P + , this implies that there exists a scaling > 0 such that (p A -p Z + ) ∈ D. Then, as C + = P + +R + D, we have p Z + + 1 ( (p A -p Z + )) = p A ∈ C + . Since A was taken arbitrarily in the virtual plane, this establishes that S = Π(O, n), i.e. the support area spans the whole virtual plane.

The geometric construction given by Proposition 1 provides fast computations of the full support area. In particular, it is faster to use this approach than to project the CWC computed by polyhedral duality methods [START_REF] Caron | Leveraging cone double description for multi-contact stability of humanoids with applications to statics and dynamics[END_REF]. (A comparison of computation times for both approaches is reported in Table I.) co nt ac t 1 co nta ct 2 contact 3 Fig. 4. Example where the support area spans the whole virtual plane. There are three contacts in total. Contacts 1 and 2 have an upward vertical normal and correspond e.g. to two feet on horizontal floor. Contact 3 is located one meter above the others and has a downward vertical normal, corresponding e.g. to a hand pushing on the ceiling above the robot. The virtual plane is taken 50 cm above contacts 1 and 2 and 50 cm below contact 3. Gray lines represent the ray generators of friction cones. In this setting, the robot's contacts are in force closure: they can generate any resultant wrench, hence any ZMP.

B. Geometric properties

Proposition 1 gives the span representation of the full support area. By construction, this area depends only on contact locations {p Ci } and on the choice of the virtual plane Π(O, n). The latter involves the choice of a fixed point O, but we will now see that the location of this point only matters as far as it defines the virtual plane Π(O, n) of the ZMP. Proof: By definition of the support area,

S = Z ∈ Π(O, n) : -→ OZ = n×τ O n•f (32)
where (f , τ O ) ∈ C c ranges over the CWC. Now, choose a point O ∈ Π(O, n) and consider

S = Z ∈ Π(O, n) = Π(O , n) : --→ O Z = n×τ O n•f . (33)
Given a contact wrench (f , τ O ), we have

--→ O Z = n × τ O n • f = n × ( --→ O O × f ) + n × τ O n • f (34) --→ O Z = --→ O O -(n • --→ O O) f n • f + -→ OZ = --→ O Z. ( 35
)
Thus, Z = Z, and since the wrench we considered is arbitrary, we have shown that S = S . By contrast, the support area does change for displacements of O along the plane normal n. Let us analyze the impact of this remaining coordinate by relaxing the assumption that (n • p Z ) = 0. We now denote by d Z the coordinate of the virtual plane, so that d

Z def = (n • p Z ) and d G def = (n • p G ).
On a horizontal floor, d Z and d G correspond to the altitude of the ZMP and COM, respectively. Given Proposition 2, we can write the virtual plane Π(d Z , n). The definition n × τ Z = 0 of the ZMP yields:

p Z = n × τ O n • f + d Z f n • f . ( 36
)
And repeating the step from Equation ( 16),

p Zi def = n × (p Ci × f i ) n • f i + d Z f i n • f i , ( 37
)
p Zi = p Ci + (d Z -d i ) f i n • f i . ( 38
)
We see from this equation that the vertices Z i are located at the intersection between the plane Π(d Z , n) and the ray (C i , f i ) of the linearized friction cone, which establishes that:

Proposition 3. The vertices of the support area are located at the intersection between the virtual plane and the rays of the friction cones.

Combining Propositions 1 and 3, we get:

Proposition 4 (Polygonal case). When all virtual pressures p i = (n•f i ) have the same sign, the support area is the convex hull of the intersections between linearized friction cones and the virtual plane.

This result is coherent with the horizontal-floor setting where the virtual plane intersects friction cones at their apexes (i.e. at contact points) and virtual pressures are all positive from contact unilaterality. A condition similar to the positivity of virtual pressures was also observed by Saida et al. (Proposition 4.1 in [START_REF] Saida | Fsw (feasible solution of wrench) for multi-legged robots[END_REF]) for some plane components of the contact wrench.

Figure 2 illustrates the geometric construction in the polygonal case. The support area (in green) is the convex hull of black points projected from friction rays (red lines). Alternatively, each individual contact surface projects its own support polygon (blue polygon in Figure 2), and the ZMP support area is the convex hull of these individual polygons. This second construction can be useful for contact planning, where the quality of contact candidates can be assessed by the expansion that they bring to the support area.

Figure 3 shows a configuration where the support area is the union of two polygonal cones. In this setting, the robot has its two feet on a horizontal floor, and pushes its hand on a wall in front of it. The positive cone C + illustrates how the wall contact expands the conventional support area between the two feet, while the negative cone C -reveals that a second support area appears behind the two floor contacts. This new area results from the ability for the robot to generate downward forces at the wall contact.

Finally, let us determine in which cases a polygonal ZMP support area can be found, which boils down to finding a suitable vector n. Proposition 5. Either contact forces can generate any arbitrary resultant force or one can find a plane normal n such that the support area is polygonal. Proof: Let C f denote the cone positively spanned by the f i 's. Then, the condition ∀i, (n

• f i ) > 0 is equivalent to n ∈ C * f , where C * f is the dual cone of C f defined by C * f = {y : ∀x ∈ C f , y • x ≥ 0}. ( 39
)
The set of solutions C * f can be computed from C f . In particular, C * f = {0} if and only if C f = R3 , i.e. the f i 's positively span the whole space. Figure 4 provides an example where contacts are in force closure. As they can generate arbitrary contact wrenches, the support area spans the whole virtual plane. Note that being able to generate arbitrary forces is a weaker condition than force closure, which usually assumes that contact forces can generate arbitrary resultant forces and momenta.

Excluding such configurations, from Proposition 5 one can choose n such that n • f > 0 for all resultant forces f generated by valid contacts. This is desirable insofar as it eliminates potential singularities n • f = 0 where the ZMP would be undefined [START_REF] Saida | Fsw (feasible solution of wrench) for multi-legged robots[END_REF]. In the classical horizontal-floor setting, taking n as the floor normal is an example of such a solution.

C. Relation with the whole-body momentum

The full support area describes the ZMPs sustainable by valid contact forces. As such, it does not take into account kinematic and dynamic limitations of the robot itself, in particular the limited changes in whole-body momentum ( Ṗ, LG ) that the robot can generate around its COM by moving its limbs. Consider the simple example depicted in Figure 5, with a single point contact. According to Proposition 4, the full support area is given by the intersection of the friction cone with the virtual plane Π(O, n). To generate a ZMP Z, the robot needs to change its angular momentum around the COM by

LG,y = -x Z f z + (d Z -d G )f x (40)
Thus, ZMPs close to the boundary of the full support area (|x Z | → d Z sin α) tend to generate higher angular momenta, a phenomenon which is amplified by the altitude d Z of the virtual plane. In practice, legged robots cannot generate a arbitrary angular momenta, as torque limits and their bounded range of motion yield additional constraints, which in turn restrict the ZMPs that the robot can generate. One way to take these limitations into account is to put additional constraints on the whole-body momentum, as we will now see with the linear pendulum mode.

IV. PENDULAR ZMP SUPPORT AREA

From a control point of view, the strong point of the ZMP lies in its direct relationship with the acceleration of the COM. This can be seen with a suitable rewriting of the Newton-Euler equations of the system: Proposition 6. The COM acceleration, ZMP position and angular momentum are bound by the equation:

pG = g + n • ( pG -g) d G -d Z -→ ZG + n × LG m(d G -d Z ) , ( 41
)
where g is the gravity vector.

Proof: 36),

Expanding τ O = p G × f + LG in Equation (
(n • f )p Z = n × (p G × f ) + n × LG + (n • p Z )f (42) (n • f ) -→ GZ = (n • -→ GZ)f + n × LG . ( 43
)
Using Equations ( 4) and ( 6), f = m( pGg), so that:

m(n • ( pG -g)) -→ GZ = m(n • -→ GZ)( pG -g) + n × LG . (44)
Equation ( 41) is a rearrangement of this last formula.

A. Linear pendulum mode

The Linear Pendulum Mode (LPM) is a particular mode of Equation [START_REF] Pham | Kinodynamic planning in the configuration space via velocity interval propagation[END_REF], obtained by regulating the COM altitude and the angular momentum to constant values: Definition 3 (Linear Pendulum Mode (LPM)). Provided with a normal vector n, the linear-pendulum assumptions are:

3 n • pG = 0 (45) LG = 0. ( 46
)
When d Z < d G , the COM is above the virtual plane with respect to the plane normal. In such situations, Equation (41) yields the well-known linear inverted pendulum mode (LIPM) [START_REF] Kajita | The 3d linear inverted pendulum mode: A simple modeling for a biped walking pattern generation[END_REF], [START_REF] Sugihara | Real-time humanoid motion generation through zmp manipulation based on inverted pendulum control[END_REF]:

Proposition 7 (Linear Inverted Pendulum Mode (LIPM)). Assume that h def = d G -d Z > 0, and define ω LIP = g/h.
The COM acceleration in the LPM is then related to the ZMP by: pG

= g -ω 2 LIP -→ GZ (47)
The LIPM was a necessity in previous works as the ZMP was assumed to lie on the floor below the COM. But now that the virtual plane coordinate d Z can be chosen freely, we may also consider the case where d Z > d G , i.e. taking the virtual plane above the COM. In this case, Equation ( 41) yields a linear non-inverted pendulum:

pG = g + ω 2 LNP -→ GZ (48)
The behavior of the ZMP with respect to COM acceleration differs between the two modes. In the LIPM, the ZMP is a repulsor of the COM: when it is maintained at a fixed position, the COM is "pushed away" from it and diverges to infinity. When selecting a virtual plane below the center of mass (d Z < d G ), ZMP trajectories will therefore follow the COM when it accelerates, and precede it when it decelerates.

In the LNPM, the ZMP is a marginal attractor of the COM: when it is maintained at a fixed position, the COM will "orbit" around it (yet without asymptotic convergence, as there is no damping term in Equation ( 48)). In a virtual plane above the center of mass (d Z > d G ), ZMP trajectories will therefore precede the COM when it accelerates, and follow it otherwise.

These two complementary behaviors reflect the fact that, in the LPM, the resultant contact force f always points from the COM to the ZMP. In what follows, we choose to take the virtual plane above the center of mass and use the LNPM.

B. Shrinking of the support area

Previous methods [START_REF] Kagami | A fast dynamically equilibrated walking trajectory generation method of humanoid robot[END_REF], [START_REF] Sugihara | Real-time humanoid motion generation through zmp manipulation based on inverted pendulum control[END_REF], [START_REF] Kajita | Biped walking pattern generation by using preview control of zero-moment point[END_REF], [START_REF] Harada | Dynamics and balance of a humanoid robot during manipulation tasks[END_REF], [START_REF] Shibuya | Trajectory planning of biped robot using linear pendulum mode for double support phase[END_REF], [START_REF] Sato | Walking trajectory planning on stairs using virtual slope for biped robots[END_REF] applied the LIPM using the convex hull of contact points (CHCP) as ZMP support area. By Proposition 4, when all contacts are coplanar the CHCP is the full support area. However, it turns out that the constraints (45)-( 46) of the linear-pendulum mode shrink the ZMP support area.

This can be seen as follows. Suppose that the ZMP is located at a vertex C k of a given contact polygon. The resultant contact force f must then be realized as one contact force f k applied at C k , while all other contact forces f j = 0 for j = k. (From Equation [START_REF] Caron | Leveraging cone double description for multi-contact stability of humanoids with applications to statics and dynamics[END_REF], p Z = p Z k = p C k implies that only the λ i 's corresponding to C k can be strictly positive.) However, in the LPM the contact force f = pGg is also directed from Z = C k to G. Hence, the ZMP support area cannot be the CHCP as soon as the COM lies outside of the friction cone of at least one ground contact point. An example of such situation is when a biped robot stretches its legs, as depicted in Figure 6.

Methods that take the Convex Hull of Contact Points (CHCP) as ZMP support areas rely by construction on the assumption of "infinite friction" to rule out such cases, at the cost of being problematic for locomotion on low-friction floors (apart from sliding, foot yaw rotations due to insufficient friction have also been observed and studied [START_REF] Cisneros | Yaw moment compensation by using full body motion[END_REF], [START_REF] Caron | Stability of surface contacts for humanoid robots: Closed-form formulae of the contact wrench cone for rectangular support areas[END_REF]). This assumption explains why the shrinking of the ZMP support area in the LPM was not highlighted in previous works.

We call pendular support area the shrunk support area that takes into consideration the Equations ( 45)-( 46) of the LPM.

C. Computation of the pendular support area

In this section, we assume that n = e Z , the upward vertical (opposite to gravity) unit vector. We now propose an algorithm Fig. 6. Situation where the pendular support area (green stripes), i.e. the ZMP support area for the Linear Pendulum Mode, is smaller than the convex hull of ground contact points (blue polygon). The robot has its legs stretched, with its two feet one meter apart on a horizontal floor (in gray). Its COM is 50 cm above ground, and the friction at contact is µ = 0.5. The ZMP cannot be located at the corner C k of the convex hull, as it would require a resultant contact force f k (in magenta) lying outside of the friction cone C k (in red). The pendular support area in this Figure has been computed using the algorithm from Section IV-C.

to compute the pendular support area based on the doubledescription method [START_REF] Fukuda | Double description method revisited[END_REF]. It turns out that computations of this set are very similar to that of the COM static-equilibrium polygon. We explain both calculations here.

1) Static-equilibrium polygon: In static equilibrium, the equations on the contact wrench are:

f = -mg (49) n × τ O = -n × (p G × mg) (50) n • τ O = 0 (51)
By expanding the triple product in Equation (50), we can rewrite them equivalently as:

E 3 0 3×3 0 1×3 n f τ O = -mg 0 (52)
p G = (n/mg) × τ O + z G n (53)
In concise form, these two equations can be written:

Aw O = b (54) p G = Cw O + d (55)
Consider the stacked vector of contact forces

f all = [f 1 • • • f n ]
. Linearized friction cones are given by linear inequalities F i f i ≤ 0. For instance, four-sided friction pyramids can be formulated as

F i =     -1 0 -µ i +1 0 -µ i 0 -1 -µ i 0 +1 -µ i     R i . (56)
Combining all F i 's in a block diagonal matrix F yields an inequality Ff all ≤ 0. Meanwhile, Equation ( 5) provides a linear mapping w O = G O f all from contact forces to the contact wrench, where G O is the so-called grasp matrix. Then, the set of valid contact forces in static equilibrium is given by:

Ff all ≤ 0 (57) AG O f all = b (58)
This expression of a polytope by linear inequalities is known as the half-space representation. By the Minkowski-Weyl theorem [START_REF] Fukuda | Double description method revisited[END_REF], all polytopes can be equivalently written in terms of linear inequalities (half-space representation) or in terms of vertices and rays (span representation). The double description method [START_REF] Fukuda | Double description method revisited[END_REF] provides a "black-box" algorithm to convert between one representation and the other. Using this tool, one can compute the span representation of this set as:

f all = i α i v i (59)
where α i > 0, i α i = 1, and the v i 's are computed vertices.

The span representation of the stability polygon is finally given by

p G = i α i (CG O v i ) + d. (60)
In the case of static stability, the solution is always a polygon, so there is no need to consider rays in this span representation. Note that this method for computing the COM staticequilibrium polygon by double-description is different from the recursive polytope projection algorithm [START_REF] Bretl | Testing static equilibrium for legged robots[END_REF].

2) Pendular support area: the Equations ( 45)-( 46) of the LPM are written, in terms of the contact wrench:

n • f = mg (61) n • τ G = 0 (62) n × τ G = 0 (63)
where g ≈ 9.81 m.s -2 is the gravity constant. Taking the resultant moment τ O rather than τ G , one can rewrite the equations above as:

1 mg z G E 3 [n×] -(n × p G ) n f τ O = p G 0 (64)
p Z = z Z -z G mg f + p G , (65)
Equation ( 64) results from (61)-(62), while Equation ( 65) is a reformulation of (63). Contrary to the previous setting, where the COM position p G resulted from the contact wrench in static equilibrium, we now assume that p G is known (e.g. measured from the instantaneous robot state). In concise form, the two equations above can be written:

A w O = b (66) p Z = C w O + d (67)
where the matrices A , C and vectors b , d now depend on p G . From there, the computations are the same as for static stability: the set of valid contact forces f all is given by

Ff all ≤ 0 (68) A G O f all = b (69)
Using the double description method, one can compute the span representation of this set as

f all = i α i v i + j λ j r j (70)
where α i > 0, i α i = 1, λ j > 0 and the v i 's (rep. r j 's) are the computed vertices (resp. rays). The span representation of the pendular support area is finally given by

p Z = i α i (C G O v i ) + j λ j (C G O r j ) + d . (71)
Similarly to the full support area, the pendular area can be either conical or polygonal. Compared with the full support area or static-equilibrium polygon, it depends not only on the set of contacts {p Ci }, but also on the instantaneous position p G of the COM. Yet, where the full support area yields only a necessary condition for contact stability (it is a 2D projection of the 6D contact wrench, the remaining four components being unconstrained), the pendular support area gives a condition both necessary and sufficient under the LPM (the four other components of the contact wrench being determined by equality constraints).

V. EXPERIMENTS A. Trajectory generation

We now design a trajectory generator for ZMP-COM trajectories based on the model-preview control formalism introduced in previous works [START_REF] Kajita | Biped walking pattern generation by using preview control of zero-moment point[END_REF], [START_REF] Audren | Model preview control in multi-contact motion-application to a humanoid robot[END_REF]. We use the ZMP as a command and COM as the output variable. First, we define the support of the ZMP trajectory p Z (γ(t)) as a line segment:

p Z (γ(t)) = γ(t)p 1 + (1 -γ(t))p 0 , (72)
where p 0 (resp. p 1 ) denotes the initial (resp. final) ZMP position. Assuming that its initial velocity is zero, the COM follows a parallel line segment p G (η(t)):

p G (η(t)) = η(t)p 1 + (1 -η(t))p 0 , (73)
where η(t) = ω 2 (γ(t) -η(t)). Next, define the state of the control problem by x(t) = [η(t) η(t)] and its command by the linear position γ(t) of the ZMP. Discretizing the time interval into K steps of duration δt, the system's dynamics become

x k+1 = cos(ωδt) 1 ω sin(ωδt) -ω sin(ωδt) cos(ωδt)
x k (74)

+ 1 -cos(ωδt) ω sin(ωδt) γ k Let X = [x 0 • • • x K ] and γ = [γ 0 • • • γ K-1 ]
. Applying (74) repeatedly, we build the matrices Φ and Ψ such that X = Φx 0 + Ψγ. We assume that the system starts with zero COM velocity, so that x 0 = 0 and X = Ψγ. We formulate the trajectory generation problem as a Quadratic Program (QP) as follows:

Objective: min w 1 c 1 (γ) + w 2 c 2 (γ)

Constraints: 0 ≤ γ ≤ 1 x K = Ψ last γ = [ 1 0 ] γ K-1 = 1 Fig. 7.
Variations of the shape of the pendular support area with the coordinate d Z of the virtual plane. In this configuration, the humanoid is making two tilted foot contacts and a left hand contact. (There is no right hand contact.) The axis directed by the plane normal n and going through the COM is depicted by a dashed line. In the proposed method, the plane altitude can be chosen freely as long as d Z = d G , however the support area becomes small when d Z and d G are close.

The objective is the weighted sum of two terms:

c 1 (γ) = 1 K k (η k -γ k ) 2 (75) c 2 (γ) = k (γ k -γ k-1 ) 2 (76)
The first one minimizes COM accelerations while the second regularizes the ZMP trajectory. The weights of the cost function are set to w 1 = 1 and w 2 = 100. Meanwhile, the constraints ensure respectively that:

• the ZMP belongs to the line segment (γ(t) ∈ [0, 1])

• the COM ends at the destination point (η K = 1) with zero velocity ( ηK = 0), • the ZMP also ends at the destination point (γ K-1 = 1). The solution to this QP provides a ZMP trajectory p Z (t), which is then integrated into a COM trajectory p G (t) by applying (48). We also added a damping term at this stage to smooth out undesired COM oscillations. On a side note, set aside the two regularization objectives, this optimization problem falls under the framework of Time-Optimal Path Parameterization (TOPP) [START_REF] Pham | A general, fast, and robust implementation of the time-optimal path parameterization algorithm[END_REF], [START_REF] Pham | Time-optimal path parameterization for redundantly actuated robots: A numerical integration approach[END_REF]. One could then trade smoothness of COM accelerations for Admissible Velocity Propagation (AVP), allowing for the integration into a kinodynamic planner of COM trajectories [START_REF] Pham | Kinodynamic planning in the configuration space via velocity interval propagation[END_REF].

With this method, contact stability is mostly enforced by checking that the segment [p 0 , p 1 ] lies inside the pendular support area computed for the extremities p 0 and p 1 of the segment. However, the COM moves as the robot performs the motion, which affects both the position and shape of the pendular support area. This phenomenon can have two undesired outcomes:

• the area becomes empty: we observed empirically that this would only happen when the COM is far away (e.g. more than one meter) from contacts, and was not a threat in practice. • the area constraints the COM motion toward divergence: past a certain COM-to-contacts distance, the support area slants away from contacts, in a way such that it contains no ZMP that could bring the pendulum back above contacts. In the middle (green) case, the COM slid to the right while the area slid by a lesser amount, so that the vertical of the COM is now on its edge: the configuration is statically marginally stable. In the rightmost (blue) case, the pendulum slid further: the vertical of the COM is no longer in the support area, and the leftmost feasible ZMP (dashed line) will steer the COM further away to the right.

Figure 8 illustrates this phenomenon. In practice, ensuring that the ZMP is well within the support area computed at the beginning p 0 and end p 1 of the segment was enough to rule out both undesirable outcomes.

B. Long stride while leaning on a wall

We implemented the whole pipeline described so far to generate multi-contact motions for a model of the HRP-4 humanoid robot. The scenario is depicted in Figure 9. The robot has to step on inclined platforms in order to reach its goal configuration on the right. Because there is no platform for its left foot in the middle of the course, the only way for it to complete the task is to put its left hand on the elevated "wall" platform while pushing with its right foot on the opposite tilted surface. Relying on these two simultaneous contacts, the humanoid can perform a long stride which would have been impossible to achieve in single-support.

As input given to solve this scenario, we assume that a contact planner provides a sequence of stances, where a stance σ i ({p Ci }, p i ) provides both a reference position of the COM p i and a set of contact points {p Ci }. The first stage of our solution computes stance-to-stance COM trajectories. To move from σ i to stance σ i+1 , the controller considers the line segment [p i , p i+1 ]. The trajectory generator is called if this segment is included in the pendular support area S(p i , d Z ) for the initial COM position. Otherwise, d Z is increased until the segment is included in S(p i , d Z ). This condition is easy to fulfill in practice, as we observed that the region S(p G , d Z ) grows like the section by the plane Π(d Z , n) of a cone passing through G (see Figure 7). In the motion depicted in Figure 9, this process lead us to take the virtual plane one meter above the center of mass of the humanoid (HRP-4 is 1.5-meter tall). Although contact planning was done by hand in this experiment, we also used the pendular support area to select COM positions p i from contact locations {p Ci }, by enforcing that the ZMP at the vertical of p i lies well inside S(p i , d Z ). Fig. 9. Snapshots of the motion generated in the Linear Pendulum Mode with a ZMP above the robot's COM. The robot has to put its left hand on a "wall" contact (in the background) and its right foot on the opposite tilted platform in order to perform an ample swing of the left leg that is otherwise impossible. In these simulations, the virtual plane is taken one meter above the robot's COM. Green polygons in this virtual plane correspond to the pendular support areas for each snapshot. Vertical lines represent the (virtual) linear pendulum. Contact forces (arrows pointing from contacts toward the robot) are computed at each time instant to cross-validate the contact-stability of the motion.

Once a reference COM trajectory p G (t) has been computed, we generate whole-body joint-angles by differential inverse kinematics (IK) under the following constraints, by decreasing task weight:

1) tracking of contacting end-effector poses, 2) tracking of COM trajectory, 3) tracking of free end-effector poses, 4) minimum variations in angular momentum 5) preferred values for some joint-angles. We applied our own IK solver for the task, which is released in the pymanoid library. 4 Similarly to [START_REF] Lee | Ground reaction force control at each foot: A momentum-based humanoid balance controller for non-level and non-stationary ground[END_REF], this solver relies on a single-layer QP problem, using the above cost function and a set of inequality constraints q ≤ min( qmax , -K s (q -q max ))

to limit joint velocities. Gains and weights used in the simulations are reported in Table II, while other simulation parameters are given in Table III. Implementation details can be checked in the library and source code released with the paper (see the motion editor distributed in [42]). ZMP tracking is not an explicit task in the list above, and is realized as a side effect of the COM tracking task. The latter comes second in the hierarchy of the differential IK solver and is therefore not fulfilled perfectly, as illustrated in Figure 10. The same holds for the task LG = 0, which has an even lower rank in the hierarchy. Yet, ZMP deviations in the generated trajectory are small enough and one can check that the ZMP always stays well within the pendular support area (see the accompanying video). We cross-validated the contact-stability of the final motion by computing, at each time step, a set of valid contact forces f all . Force vectors are depicted in Figure 9 and displayed in the accompanying video [42]. 18.5 ± 3.7 33.9 ± 12.2 67.9 ± 11.9

C. Computation times

Table I reports average computation times for both kinds of support area during the execution of the motion from Figure 9. Simulations were run on an average laptop computer with an Intel(R) Core(TM) i7-6500U CPU @ 2.50 Ghz. The source code for these simulations is located in the full support area folder of [42].

The number of samples for each average in Table I is around 1000 for one-and two-contact configurations, and around 100 for three-contact configurations. The first two lines show that computing the full support area is faster using our geometric construction than by projecting the CWC [START_REF] Caron | Leveraging cone double description for multi-contact stability of humanoids with applications to statics and dynamics[END_REF], which is expected since the latter is a more complex 6D cone. These two criteria are always computed faster than a pendular support area, as the latter adds equality constraints to the former, while neither of them takes into account limits on the angular momentum. Also, recall that the full support area is only a necessary condition for contact stability.

The first line of the second group corresponds to the algorithm described in Section IV-C. We also implemented a variant of the projection algorithm from [START_REF] Bretl | Testing static equilibrium for legged robots[END_REF] to compute the pendular support area. 5 In practice, this variant is one to three times slower than CDD, except for triple contacts where we observed a significant slowdown of CDD, along with some freezes. This point is implementation-specific rather than a limitation of the underlying double-description algorithm. 6

VI. CONCLUSION AND FUTURE WORK

In this paper, we have generalized the notion of ZMP support areas to take into account frictional constraints and multiple non-coplanar contacts. First, we derived the geometric construction of the full support area, which contains all ZMPs generated by valid contact forces. Then, we moved to the control of humanoid robots in the Linear Pendulum Mode. We noticed that the constraints stemming from this mode shrink the support area (even on horizontal floors), and proposed an algorithm to compute the new area, which we called pendular support area. Armed with these new conceptual tools, we designed a whole-body controller for locomotion across arbitrary multi-contact stances, which we demonstrated in simulation with a model of the HRP-4 humanoid robot.

Support areas are a general property of contact wrenches and can be applied in all fields related to mobility, including locomotion, grasping or workpiece fixturing. As a mathematical object, they are 2D non-linear projections of the 6D contact wrench cone. From there, one may ask: is it the most general we can do, or can a 3D projection take into account more components of the resultant moment? For the interested reader, we provide some elements of answer in Appendix A, although the question itself remains open.

In practice, we believe that the ability to take ZMPs in virtual support areas paves the way for multiple future developments. For contact planning, we noticed how individual support polygons can be used to measure the "increase of stability" of a prospective contact. For whole-body control, our approach can be followed in the more general framework of a-priori dimensionality reduction of control problems. At one end of the spectrum, the full support area corresponds to a controller with state variables for all six components of the whole-body momentum, such as [START_REF] Nagasaka | Whole-body control of a humanoid robot based on generalized inverse dynamics and multicontact stabilizer that can take acount of contact constraints[END_REF], [START_REF] Audren | Model preview control in multi-contact motion-application to a humanoid robot[END_REF]. The task of such a controller is harder, but the area depends only on contact locations. At the other end of the spectrum, the pendular support area corresponds to fixing four momentum components, thus reducing control to a two-dimensional problem. The task of the controller is then simpler, at the cost that the area depends on the position of the COM. Between these two ends of the spectrum lies a hierarchy of intermediate three-, four-or five-dimensional problems that can be computed, in a similar fashion, by a-priori reduction of equality constraints. Exploring these possibilities is the focus of our current research.

F u l l S u p p o r t A r e a ( i n c l u d e s f r i c t i oFig. 1 .

 1 Fig.1. Overview of the construction proposed in this paper. The full ZMP support area, including pressure and frictional constraints, is computed in an arbitrary virtual plane (here, above the robot's head). For locomotion, linearized pendulum dynamics are obtained by regulation of the angular momentum. The shrinking of the support area incurred by the Linear Pendulum Mode is fully taken into account. A whole-body controller based on these developments enables multi-contact locomotion in arbitrary environments.

  Let us then partition the set of generator indices I into I + def = {i | p i > 0} and I -def = {i | p i < 0}. For any K ⊂ I, denote by σ(K) def = i∈K λ i |p i | and define

Fig. 2 .

 2 Fig.2. Geometric construction of the ZMP support area in the polygonal case. Ray generators of friction cones (red lines) are traced until they intersect the virtual plane, yielding points inside the support area (black dots). The support area (in green) is the convex hull of this set of points (Proposition 4). The blue polygon corresponds to the individual support polygon of the right contact. These polygons describe the contribution of each contact to the expansion of the support area, and can e.g. be used in contact planning to evaluate how new contact candidates would expand the area.

Proposition 2 .

 2 The support area S does not depend on the coordinates of the reference point O in the virtual plane Π(O, n).

Fig. 5 .

 5 Fig. 5. Illustration of the full support area (green) for a single point contact. Having the ZMP Z (red) at the boundary of the support area generates an angular momentum LG,y (blue) around the center-of-mass G which increases with both the altitude d Z of the virtual plane and the proximity of Z to the boundary of the support area.

Proposition 8 (

 8 Linear Non-inverted Pendulum Mode (LNPM)). Assume that h def = d Z -d G > 0, and define ω LNP = g/h . The COM acceleration in the LPM is then related to the ZMP by:

Fig. 8 .

 8 Fig.8. Three pendular support areas and their corresponding pendulums at rest configurations, i.e. with the ZMP at the vertical of the COM. The HRP-4 posture corresponds to the leftmost (red) configuration, where the vertical of the COM is centered in the area. In the middle (green) case, the COM slid to the right while the area slid by a lesser amount, so that the vertical of the COM is now on its edge: the configuration is statically marginally stable. In the rightmost (blue) case, the pendulum slid further: the vertical of the COM is no longer in the support area, and the leftmost feasible ZMP (dashed line) will steer the COM further away to the right.

Fig. 10 .

 10 Fig. 10. COM and ZMP trajectories (resp. in green and red) for the motion of Figure 9. Reference trajectories are straight dashed lines between via points. Five pendular support areas corresponding to single-and double-support configurations are plotted for reference (green polygons). The accompanying video [42] shows the evolution of support areas throughout the motion.

TABLE I COMPUTATION

 I TIMES (IN MS) FOR THE FULL SUPPORT AREA AND CWC[START_REF] Caron | Leveraging cone double description for multi-contact stability of humanoids with applications to statics and dynamics[END_REF], AS WELL AS FOR THE PENDULAR SUPPORT AREA USING[START_REF] Fukuda | Double description method revisited[END_REF] (CDD) OR[START_REF] Bretl | Testing static equilibrium for legged robots[END_REF] (B&L). BOLD INDICATES FASTEST IN EACH CATEGORY.

	Stability Criterion	One contact Two contacts	Three contacts 6

CWC (CDD) 1.1 ± 4.2 4.5 ± 1.8 9.5 ± 2.6 Full support area 0.8 ± 0.4 2.6 ± 2.1 6.4 ± 3.1 Pendular area (CDD) 6.4 ± 3.3 24.8 ± 4.7 385.3 ± 251.3 Pendular area (B&L)

TABLE II GAINS

 II AND WEIGHTS USED IN THE DIFFERENTIAL IK TRACKER (N/A: NO GAIN FOR TASKS REGULATING ACCELERATIONS)

	TASK DESCRIPTION	GAIN [Hz]	WEIGHT
	Contacting end-effector	1	100
	Free end-effector	0.03	1
	Center of mass tracking	1	5
	Angular momentum variations	N/A	0.2
	Joint-limit gain Ks	50	N/A
	Velocity smoothness	N/A	1
	Preferred joint-angles	0.05	0.1

Sardain and Bessonnet [5] pointed out that the term "zero moment point" is misleading, as the moment at these points has actually a nonzero component along the surface normal. They suggested that "zero-tilting moment point" would be a more proper name.

We assume n is chosen so that none of them is zero, which is easy to do since there is only a finite set of generators {f i }.

Strictly speaking, our definition includes both n× LG = 0 and n• LG = 0, while only the former is necessary to realize the linear pendulum mode.

See the contact stability ROS package in[42] for details.

Averages for "Pendular area (CDD)" in triple-contact are reported for samples where CDD did not freeze, amounting for around 20% of all samples.

ACKNOWLEDGMENT This research was supported by 2014-2015 NEDO International R&D and Demonstrative Project in Environmental and Medical Fields / International R&D and Demonstrative Project in Robotics Fields / "Research and Development of Robot Open Platform for Disaster Response" (PI: Yoshihiko Nakamura).

APPENDIX A PERSPECTIVES ON A THREE-DIMENSIONAL ZMP

Considering Proposition 2, the ZMP decouples the moment of a wrench from the position at which it is taken, yet at the "cost" of one dimension. A natural question is then: could a three-dimensional ZMP perform a similar decoupling for all three coordinates of the moment? Unfortunately the answer seems to be negative, at least in the following sense:

In other words, at least one coordinate of the ZMP depends on the reference point O. Let us then consider the remaining moment coordinate (n • τ O ). Applying the same calculations as those following Equation [START_REF] Saida | Fsw (feasible solution of wrench) for multi-legged robots[END_REF], one can write:

We then define a spatial point including all three coordinates, which we call the n-Moment Point (n-MP for short):

The vertices of its support volume V can be computed in the same fashion as in Section III by

The geometric construction of support areas can also be applied mutatis mutandis to V:

• when all virtual pressures n • f i have the same sign, V is the convex hull of the above vertices, • otherwise, it is the union of two polyhedral convex cones built on the Minkowski difference of positive-and negative-pressure polyhedra. A complete implementation of this construction can be found in the n moment point folder of the accompanying source code [42].

The n-MP is a three dimensional spatial point equivalent to the moment τ O , in the sense that one can be computed from the other by

In other words, M represents the screw coordinates of the contact wrench along its non-central axis directed by n, with magnitude n • f and pitch n • --→ OM . However, adding the third moment coordinates makes the shape of the support volume V depend on the choice of the reference point O. Formally: Proposition 10. There is no non-empty subspace of displacements --→ OO of the reference point O, independent from the resultant f , that leaves the n-MP invariant.

Proof: Consider a displacement --→ OO of O in the plane. From Equation (78), it results in a variation --→ OO • n×f n•f of the n-MP coordinate along n. This term needs to be zero for any displacement leaving the n-MP invariant, thus --→ OO is parallel to either n or f . The former would yield a variation of the plane coordinates of the n-MP (that is to say, of the ZMP). The latter is excluded as we are looking for an invariance that is independent from the resultant force.