# An analytical model of rolling contact and its application to the modeling of bipedal locomotion

Justin Carpentier, Andrea del Prete, Nicolas Mansard, Jean-Paul Laumond

► To cite this version:

Justin Carpentier, Andrea del Prete, Nicolas Mansard, Jean-Paul Laumond. An analytical model of rolling contact and its application to the modeling of bipedal locomotion. IMA Conference on Mathematics of Robotics, Sep 2015, Oxford, United Kingdom. hal-01182733

HAL Id: hal-01182733

https://hal.science/hal-01182733

Submitted on 5 Aug 2015

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# An analytical model of rolling contact and its application to the modeling of bipedal locomotion

By J. Carpentier, A. Del Prete, N. Mansard and J. Laumond

LAAS CNRS, Univ. de Toulouse, UPS, 7 avenue du Colonel Roche, F-31400 Toulouse, France

## Abstract

We propose an original formulation to describe the motion of a rolling object in contact with a flat surface, and propose a complete application of this model to generate optimal walk of a humanoid system. We derive an analytic formulation of the contact equations which does not require any numerical approximation. We then show how the model performs when applied to the simulation of bipedal locomotion: by replacing polyhedral models of the feet by ellipsoidal ones, the numerical algorithms run faster and provide higher quality results. Furthermore, as any convex shape may be approximated locally by an ellipsoid, the scope of the model goes beyond locomotion simulation.

## 1. Introduction

Human body is a complex machinery, both highly redundant and soft. It is difficult, if not impossible, to establish a complete computational model accounting for all physical interactions of the body immersed and acting in the real world. Such a model would have to gather sparse knowledge from life sciences (e.g., biomechanics and neuroscience) and from engineering (e.g., computer animation and robotics). As a consequence of this difficulty, it is legitimate to resort to approximated models. As an example, most studies exploring the contact relationship between the body and a surface refer to finite element methods. Such methods model the body parts by means of polyhedral surfaces. However polyhedral approximation may constitute a disputable premise. In this paper, we show how more accurate contact models may improve both the understanding of contact physics, and the performance of numerical simulations. The context of the study is human locomotion. By replacing polyhedral models of the feet with ellipsoidal ones, it is possible to derive an analytical formulation of the foot-floor relationship. As benefits, the numerical simulations based on this formulation perform faster and give feasible results when compared to existing methods.

## 2. The Rolling-Without-Slipping contact model

In this section, we derive the equations of the contact between two smooth objects living in the real Euclidian space $E^3$. These equations are generic and can be applied in any simulation context, as in robotics for the modeling of foot or finger contact. Up to a certain level, they can be derived for any locally convex surfaces, which we do first in Section 2.1. However, some terms depend on the local variation of the surface, which can not be written as an algebraic form in the general case. We then specialize the derivation for ellipsoids and obtain the full constraint in closed form in Section 2.2. Directions to extend the developments to similar shapes (i.e. which are locally equivalent to ellipsoid) are given in Section 2.3.

The rolling-without-slipping problem of two shapes has previously been addressed in
## An analytical model of rolling contact and its application to bipedal locomotion modeling2

robotics, mainly in the context of dexterous manipulation as in Cole et al. (1989), Li et al. (1989) or more recently in Kry and Pai (2003). In general, authors use minimal coordinates parameterization to work on manifolds, suffering from potential singularities when the nonholonomic constraint is derived, due mainly to the closure of the chosen minimal coordinate set. A complete atlas must then be chosen in order to avoid those singularities by changing the map. Such strategy is tedious to implement in particular when regularity is needed in dynamic equations, like in sensitivity analysis.

### 2.1. The rolling without slipping condition

The rolling without slipping (RWS) of two smooth manifolds $M_1$ and $M_2$ one on each other is a classical topic in mechanics and mathematics, see Chitour et al. (2014). From a kinematic view point, the RWS condition can be stated as: the relative velocity of the two contact points between the two manifolds is zero. This statement leads to a nonholonomic constraint, i.e. a constraint linking the configuration and the velocity of both manifolds.

#### 2.1.1. General notations

We assume the two shapes to be regular manifolds of dimension 2 composed of one sheet, embedded in the real Euclidian space $E^3$, considered similar to its underlying real vector space $\mathbb{R}^3$. Each shape is characterized by an implicit equation $\Phi_i : \mathbb{R}^3 \to \mathbb{R}$ expressed in a local coordinate frame (do Carmo (1976)), i.e. $M_i \triangleq \{x \in \mathbb{R}^3, \Phi_i(x) = 0\}$, $i = 1, 2$. Given an arbitrarily world frame $\mathcal{T}_O$ centered in $O \in \mathbb{R}^3$, $T_i \in SE(3)$ the configuration of a frame centered in $C_i \in \mathbb{R}^3$ and rigidly attached to $M_i$. $M_i$ is equivalently defined by $M_i = \{P \in \mathbb{R}^3, \Phi_i(T_i^{-1}P) = 0\}$. In order to get at most one contact point, both shapes are assumed to be convex, one of them being strictly convex (Gilbert and Johnson (1985)).

#### 2.1.2. Identification of the contact points

The location of the contact points between two surfaces depends directly on their shapes $\Phi_i$ and their configurations $T_i$ in the ambient space. This search of contact points location can be settled as the following minimization problem:

$$ \min_{\mathbf{x}_1, \mathbf{x}_2 \in \mathbb{R}^3} \frac{1}{2} \|\mathbf{x}_1 - \mathbf{x}_2\|_2^2 \quad \text{subject to} \quad T_i^{-1}\mathbf{x}_i \in M_i, i = 1, 2 \qquad (2.1) $$

where $\|\cdot\|_2$ corresponds to the classical Euclidian 2-norm, and $\mathbf{x}_i$ is the coordinates in the global frame of a point belonging to $E^3$, expressed in the global frame. In general, this problem has no closed solutions and it is necessary to solve it numerically, like in Chakraborty et al. (2008). We will see in Section 2.2 that there are some practical cases where we can find a closed-form solution.

#### 2.1.3. The rolling-without-slipping condition

With the previous notation, it is now possible to derive the complete formulation of the RWS constraint. First, we write down the velocity of the contact point $P_i$ relative to the world frame:

$$ v(P_i) \triangleq \frac{d\overrightarrow{OP_i}}{dt} = v(C_i) + \omega_i \times \overrightarrow{C_i P_i}, \qquad (2.2) $$

with the linear $v(C_i)$ and angular $\omega_i$ velocities of the frame $T_i$ wrt the world. These two vectors can be stacked into $\dot{v}(C_i) \triangleq \begin{bmatrix} v(C_i) \\ \omega_i \end{bmatrix}$, the twist or spacial velocity (Featherstone (2007)) associated to the configuration $T_i$. From now, and for the rest of the developments, we note the time derivate operator as $d_t \triangleq \frac{d}{dt}$.
Then, as we mentioned earlier, the nonholonomic constraint corresponds to a zero relative velocity between the two contact points. This can be written as:

$$d_t \overrightarrow{P_2 P_1} = d_t \overrightarrow{OP_1} - d_t \overrightarrow{OP_2} = \mathbf{0} = v(P_1) - v(P_2) \triangleq g(T_1, \hat{v}_1, T_2, \hat{v}_2), \quad (2.3)$$

where we introduce g, a function of both configurations and spatial velocities.

### 2.1.4. Derivation of the nonholonomic constraint

Eq. (2.3) is sufficient from a kinematic view point, but we need here to rewrite it in terms of the body acceleration in order to add the RWS constraint into the dynamical equations of the multi-body walking system (introduced in Section 3). This reduction is obtained by differentiating (2.3):

$$d_t g(T_1, t_1, T_2, t_2) = d_t (v(P_1) - v(P_2)) = \mathbf{0} \quad (2.4)$$

Developing the time derivate of $v(P_i)$, we obtain:

$$d_t v(P_i) = a(C_i) + \dot{\omega}_i \times \overrightarrow{C_i P_i} + \omega_i \times d_t \overrightarrow{C_i P_i}, \quad (2.5)$$

which corresponds to the acceleration of a fixed material point attached to a moving frame. The main challenge here is to compute the variation of the contact point onto the manifold $\frac{dC_i P_i}{dt}$. It does not seems possible to find a simple expression of this term except in certain cases, as the example developed in the next subsection.

### 2.2. Practical case : the ellipsoid rolling without slipping on the plane

The rolling without slipping of an ellipsoid on a plane (Fig. 1) is a practical case in which we can find an explicit formulation for (2.5). To simplify the following development, we consider that the plane $M_2$ is fixed in the global frame. The implicit equation $\Phi_1$ describing the ellipsoid $M_1$ can be defined as:

$$\Phi_1(\mathbf{x}) = \mathbf{x}^\top A_1 \mathbf{x} - 1 \quad (2.6)$$

with $A_1$ a symmetric, positive-definite matrix of dimension 3 whose axes are expressed in the local coordinate frame of the ellipsoid.

<figure><img src="image_2.png" /><figcaption>Figure 1: An ellipsoid on a plane, with $P_1$ and $P_2$ the two contact points coinciding at the same position in the world frame.</figcaption></figure>

#### 2.2.1. Identification of contact points

The normal to the ellipsoid surface is the direction of the gradient of $\Phi_1$ (do Carmo (1976)), while for the plane, it corresponds to the constant normal vector $n_2$. Following Chakraborty et al. (2008), the contact points are defined by opposite contact normal vectors, i.e. $\mathbf{x} \in M_1$ is in contact iff there exists $\lambda \in \mathbb{R}_+$ such that:

$$\nabla_x \Phi_1(\mathbf{x}) = -\lambda \, {}^1T_2 n_2 \quad (2.7)$$

where ${}^1T_2 \triangleq T_1^{-1} \circ T_2$ relates to the relative configuration between the two shapes, with $\circ$ corresponding to the composition law of $SE(3)$ group. Introducing ${}^1n_2 \triangleq {}^1T_2 n_2$ as the expression of the normal to the plane $n_2$ expressed in the ellipsoid local frame, (2.7) leads to:

$$2A_1 \, {}^1x = -\lambda \, {}^1n_2 \Leftrightarrow \, {}^1x = -\frac{\lambda}{2} A_1^{-1} \, {}^1n_2 \quad (2.8)$$
An analytical model of rolling contact and its application to bipedal locomotion modeling4

Introducing this equation in (2.6) and settings it equal to zero, we obtain:

$$ \lambda = 2\sqrt{\mathbf{n}_2^\top \mathbf{A}_1^{-1} \mathbf{n}_2} \qquad (2.9) $$

and the ellipsoid contact point is finally:

$$ \mathbf{x} = \frac{-\mathbf{A}_1^{-11}\mathbf{n}_2}{\sqrt{\mathbf{n}_2^\top \mathbf{A}_1^{-1} \mathbf{n}_2}} \qquad (2.10) $$

2.2.2. Variation of the contact point onto the ellipsoid

The variation of the contact point is obtained by differentiating (2.10). The plane is not moving and the variation of the components of $\mathbf{A}_1^{-1}$ are only due to the rotation of the frame attached to the ellipsoid, which leads to the operator written (improperly) as $\dot{\mathbf{A}}_1^{-1}$. acting on vectors of $\mathbb{R}^3$:

$$ \dot{\mathbf{A}}_1^{-1} \triangleq d_t \mathbf{A}_1^{-1} = ( [\omega_1]_\times \mathbf{A}_1^{-1} - \mathbf{A}_1^{-1} [\omega_1]_\times ) \qquad (2.11) $$

Using this operator, the variation of contact point on the ellipsoid is finally given by:

$$ d_t \mathbf{x} = \frac{\mathbf{n}_2^\top \mathbf{x} \dot{\mathbf{A}}_1^{-11} \mathbf{n}_2 - \frac{1}{2} \mathbf{n}_2^\top \dot{\mathbf{A}}_1^{-11} \mathbf{n}_2 \mathbf{x}}{\sqrt{\mathbf{n}_2^\top \mathbf{A}_1^{-11} \mathbf{n}_2}} \qquad (2.12) $$

2.3. Extension to more complex cases

We can interpret equation (2.12) as the variation of the contact point on the ellipsoid according to its own motion. When the plane also moves, (2.12) must be augmented with a second term corresponding to the relative velocity. Similarly, the same equations can be easily extended to the case of two ellipsoids in contact. The main challenge remains the computation of the contact points between the two ellipsoids. Indeed, and to our best knowledge, there is no analytical solution to this problem. The expression is then composed of the variation due to both normals $\mathbf{n}_1$ and $\mathbf{n}_2$. Finally, the equation remains the same for the shapes that locally look like an ellipsoid, typically a super-ellipsoid Barr (1981). In that case, there are no analytical solutions to obtain a contact point, this problem needs to be solved numerically. Once the contact point is known, the contact constraint (2.12) remains true, with $\Phi_1$ the ellipsoid as the second order approximation of the surface at the contact point.

3. Embedding the contact model into an articulated rigid body system

We are now interested in incorporating our contact surfaces in a free-floating base articulated system $k_i$. This is done by linking the classical and spatial acceleration of the frame attached to the shape. We define the generalized configuration of the system by $q_{k_i} \in Q_{k_i} \triangleq SE(3) \times \mathbb{R}^{n_{k_i}}$, which we represent by a vector space, selecting any appropriate representation for the base orientation. We note improperly its velocity and acceleration vectors $\dot{q}_{k_i}$, $\ddot{q}_{k_i} \in \mathbb{R}^{6+n_{k_i}}$ respectively, while they actually belong to the tangent space of the configuration $T_{q_{k_i}} Q_{k_i}$.
## 3.1. Embedding the constraint into a rigid-body articulated system

### 3.1.1. Multibody spatial algebra

Based on the notation used by Featherstone (2007), the time derivative of spatial velocity $\dot{v}(C_i)$, named the spatial acceleration, is denoted as:

$$ \hat{\mathbf{a}}_{C_i} = \begin{bmatrix} \dot{\omega}_i \\ \dot{v}_{C_i} \end{bmatrix}, \quad \text{with the notation } \mathbf{v}_{C_i} = v(C_i) \qquad (3.1) $$

The spatial acceleration is interesting because, as the derivative of the spatial velocity, it is easily linked with the joint acceleration:

$$ \hat{\mathbf{v}}_{P_i} = J_{T_{P_i}} \dot{\mathbf{q}}_{k_i} \quad \text{and:} \quad \hat{\mathbf{a}}_{P_i} = J_{T_{P_i}} \ddot{\mathbf{q}}_{k_i} + \dot{J}_{T_{P_i}} \dot{\mathbf{q}}_{k_i} \qquad (3.2) $$

with $J_{T_{P_i}} \triangleq \frac{\partial T_{P_i}}{\partial q_{k_i}}$ the jacobian of the frame $T_{P_i}$ at the configuration $q_{k_i}$. Contrary to the velocity, the linear component of $\hat{\mathbf{a}}_{C_i}$ does not correspond to the acceleration of the point $C_i$, $a(C_i)$, but we have:

$$ a(C_i) = \dot{v}_{C_i} + \omega_i \times v(C_i) \qquad (3.3) $$

Knowing that, we are able to rewrite (2.5) as:

$$ d_t v(P_i) = \mathbf{a}_{P_i}^t + \omega_i \times v(C_i) + \omega_i \times d_t \overrightarrow{C_i P_i} \qquad (3.4) $$

### 3.1.2. Linking the constraint to the generalized coordinates and their derivatives

The RWS constraint (2.4) is rewritten in terms of $\mathbf{q}$ and its derivatives by injecting (3.2) into (3.4):

$$ J^t \ddot{\mathbf{q}} + \gamma(\mathbf{q}, \dot{\mathbf{q}}) = \mathbf{0} \qquad (3.5) $$

with:

$$ \mathbf{q} = \begin{bmatrix} \mathbf{q}_{k1} \\ \mathbf{q}_{k2} \end{bmatrix}, \dot{\mathbf{q}} = \begin{bmatrix} \dot{\mathbf{q}}_{k1} \\ \dot{\mathbf{q}}_{k2} \end{bmatrix}, \ddot{\mathbf{q}} = \begin{bmatrix} \ddot{\mathbf{q}}_{k1} \\ \ddot{\mathbf{q}}_{k2} \end{bmatrix}, J^t = \begin{bmatrix} J_{T_{P1}}^t & -J_{T_{P2}}^t \end{bmatrix} \qquad (3.6) $$

$$ \text{and } \gamma(\mathbf{q}, \dot{\mathbf{q}}) = \dot{J}^t \dot{\mathbf{q}} + \omega_1 \times (v(C_1) + d_t \overrightarrow{C_1 P_1}) - \omega_2 \times (v(C_2) + d_t \overrightarrow{C_2 P_2}) \qquad (3.7) $$

## 3.2. The dynamical walker model

We consider an avatar that is a tree-like rigid-body articulated system with anthropomorphic proportions following Dumas et al. (2007). and 29 degrees of freedom (23 are actuated, the others are related to the free-floating base). In order to approximate the rolling of the feet during the locomotion, we equipped the avatar with ellipsoidal feet.

### 3.2.1. Contact phases

Locomotion consists in a succession of several contact phases, each of them corresponding to a specific interaction between the articulated system and its environment. For example, bipedal walking is a cyclic process alternating phases of single support when only one foot touches the ground, and phases of double support when both are in contact with the ground.

Thereafter, the adopted approach consists in splitting the locomotion process into its natural phases and fixing the dynamics on each of these phases as well as the terminal contact points. The main idea is to have a strict control of the dynamics of the avatar in order to guide the solver in its search process.
## An analytical model of rolling contact and its application to bipedal locomotion modeling6

### 3.2.2. Dynamical equations

Taking the same notations as previously, the dynamical equation of the articulated system can be written as:

$$M(q)\ddot{q} + b(q, \dot{q}) = S^{\top}\tau + \sum_k J_k^{\top}(q)\lambda_k \quad \text{with} \quad J_k\ddot{q} + \gamma_k(q, \dot{q}) = 0, \quad (3.8)$$

where $M(q)$ is the mass matrix, $b(q, \dot{q})$ corresponds to the Coriolis, centrifugal and gravitational effects, $S$ a selection matrix, subscript $k$ corresponds to the index of each rolling-without-slipping constraint, the total number being fixed by the considered phase. Subsequently, we skip writing $k$ subscript. Let now the state and the control of the dynamical system be $x = (q, \dot{q})$ and $u = \tau$, and the transfer function of the system is:

$$\ddot{q} = f(x, u) \triangleq M(x)^{-1} (S^{\top}u - b(x) + J(x)^{\top}\lambda(x, u)), \quad (3.9)$$

where $\lambda(x, u)$ is the contact forces obtained when applying the joint torques $u$ at the state $x$. They are obtained as the Lagrange multipliers associated with the following quadratic program (QP):

$$\min_{\dot{q}} \|\ddot{q} - M^{-1}(S^{\top}u - b)\|_{M}^{2} \quad \text{w.r.t.} \quad J\ddot{q} + \gamma = 0, \quad (3.10)$$

where $\|x\|_M \triangleq \sqrt{x^\top M x}$.

## 4. Optimal Control with rolling contacts

This section highlights the exactness and the interest of the proposed models introduced in Sections 2 and 3, by integrating it in a numerical optimal-control framework used to generate efficient walking patterns. We first introduce the numerical optimal-control program used in the implementation. A first case study is presented, studying the optimal control of a ball rolling on a flat ground. Finally, we investigate the emergence of walking patterns in various scenarios.

### 4.1. The numerical optimal control framework

The generic optimal control problem (OCP) with periodization is written as (Schultz and Mombaur (2010)):

$$\min_{\underline{x} \in \mathcal{X}, \underline{u} \in \mathcal{U}} \quad \int_{0}^{T} l(t, x(t), u(t))dt + l_{T}(T, x(T)) \\ \dot{x} = f_{\Psi(t)}(t, x(t), u(t)) \quad (4.1)$$

$$x(0) = x(T) \quad (4.2)$$

$$g_{\Psi(t)}(x(t_{\Psi(t)}^{0})) = \mathbf{0} \quad (4.3)$$

where $\underline{x}$ and $\underline{u}$ are functions of time acting respectively on the state $\mathcal{X}$ and control $\mathcal{U}$ spaces, $\Psi(t) : \mathbb{R} \mapsto \mathbb{N}$ corresponds to the index of the phase at time $t$, (4.2) is the periodicity constraint and $l$ and $l_T$ being the path and terminal cost respectively. The last equation (4.3) imposes additional constraints at the beginning of each phase ($t_{\Psi(t)}^0$ denoting the starting time of the phase $\Psi(t)$). We then transform this infinite dimensional problem into a finite one, by choosing a finite set for the control variable $\underline{u}$, like the set of piecewise constant functions. Many discretization scheme can be chosen in order to solve the OCP. In our case, we chose a multiple shooting approach through the solver MUSCOD described in Leinweber et al. (2003). Although this approach is slightly difficult to implement in practice, it stabilizes the integration scheme and provides an accurate estimation of both the dynamics and the cost function.
<figure><img src="image_3.png" /><figcaption>Figure 2: Curves showing the torques for the hip, knee and ankle joints of the right leg and snapshots of the three walking movements: flat ground, tilt of +10 deg and tilt of -10 deg.</figcaption></figure>

## 4.2. Generation of bipedal locomotion

We are now interested in studying the emergence of walking patterns by the use of ellipsoidal feet and the RWS constraint. To do so, we first fixed the path cost function $l(t, x(t), u(t)) = \|u(t)\|_P^2$, with $P$ a weighting matrix of dimension $n$ (corresponding to the degrees of freedom), whose coefficients have been suited according to standard human torques during locomotion. In order to generate normal walking pattern, only two phases are needed: the single and double supports phases. We therefore imposed these two phases, each of them having its own dynamics (i.e. one or two rolling feet on the ground). It is worth to mention here that the duration of each phase was a free-variable of the global OCP. Finally, we used (4.3) to impose the points on the ellipsoid which need to be in contact with the ground at the beginning of each phase. The experimental scenarios correspond to three different slopes: 10, 0 and -10 degrees for the tilts and each scenario was computed in less than 20 min on a modern computer.

The results are depicted Fig. 2. It is a composition of the generated patterns corresponding to one of the aforementioned scenarios, with plots related to the torques of the hip, knee and ankle joints of the right leg during the complete gait cycle. Several remarks can be raised. Firstly, the relative durations of each phase on the complete cycle are not the same: heel strike event appears around 40% of the cycle for the flat case against 55% for the two others. Secondly, the torque consumption is larger in the down-hill scenario than in the up-hill one, which may surprise. Indeed, down-hill requires slowing down the system, and therefore requires more torques.

# 5. Conclusion and Future Work

The main technical contribution of this paper is the derivation of a complete model of the rolling contact, that was yet only partially derived and only for specific parametrization, with the drawback of introducing artificial singularities in the dynamics. We applied this
An analytical model of rolling contact and its application to bipedal locomotion modeling8

model for the optimal-control of a ball rolling on plane and in the context of bipedal
locomotion modeling.

Two main applications are expected from this work. First, efficient walker robots can not be achieved with flat feet. The next easy step is to consider toe joints. Beyond one-axis toes, more complex foot shapes, based on an anthropomorphic structure or built from a different principle, are the key to a really efficient walk. Our method is an important milestone to understand how to control such feet. The second application is the study of the human walk, and in particular the estimation of the forces exerted by the legs during a natural walk. Estimation can be equivalently rewritten as an optimal control problem optimizing the likelihood of the measures. However, the accuracy of the estimation directly depends on the relevance of the modeled human body.

Acknowledgment

This work is supported by the European Research Council (ERC) through the Actanthrope
project (ERC-ADG 340050), the European project KOROIBOT (FP7-ICT-2013-10/611909)
and the French National Research Agency (ANR) project ENTRACTE.

REFERENCES

A. BARR. Superquadrics and angle-preserving transformations. *IEEE Computer Graphics and Applications*, 1981.

N. CHAKRABORTY, J. PENG, S. AKELLA, AND J. E. MITCHELL. Proximity queries between convex objects: An interior point approach for implicit surfaces. *IEEE Transactions on Robotics*, 2008.

Y. CHITOUR, M. MOLINA, AND P. KOKKONEN. The rolling problem: overview and challenges. In *Geometric Control Theory and Sub-Riemannian Geometry*. Springer, 2014.

A. B. COLE, J. E. HAUSER, AND S. S. SASTRY. Kinematics and control of multifingered hands with rolling contact. *IEEE Transactions on Automatic Control*, 1989.

M. DO CARMO. *Differential Geometry of Curves and Surfaces*. Prentice-Hall, 1976.

R. DUMAS, L. CHEZE, AND J.-P. VERRIEST. Adjustments to mcconville et al. and young et al. body segment inertial parameters. *Journal of biomechanics*, 40(3):543-553, 2007.

R. FEATHERSTONE. *Robot Dynamics Algorithms*. Springer, 2007.

E. GILBERT AND D. JOHNSON. Distance functions and their application to robot path planning in the presence of obstacles. *IEEE Journal of Robotics and Automation*, 1985.

P. G. KRY AND D. K. PAI. Continuous contact simulation for smooth surfaces. *ACM Transactions on Graphics (TOG)*, 2003.

D. LEINWEBER, I. BAUER, H. BOCK, AND J. SCHLOEDER. An efficient multiple shooting based reduced {SQP} strategy for large-scale dynamic process optimization. Part 1: theoretical aspects. *Computers and Chemical Engineering*, 2003.

Z. LI, J. CANNY, AND S. SASTRY. On motion planning for dexterous manipulation. i. the problem formulation. In *IEEE International Conference on Robotics and Automation*, 1989.

G. SCHULTZ AND K. MOMBAUR. Modeling and optimal control of human-like running. *Mechanics, IEEE/ASME Transactions on*, 15(5):783-792, 2010.