Fast Trajectory Optimization for Legged Robots using Vertex-based ZMP Constraints

TL;DR

This paper presents a fast, unified trajectory optimization method for legged robots that uses a vertex-based ZMP support-area representation, enabling rapid generation of diverse gaits and transitions validated on a real quadruped.

Contribution

It introduces a vertex-based support-area constraint that handles arbitrary contact orientations, enabling quick, versatile motion planning for legged robots.

Findings

Generated motions in less than a second

Successfully demonstrated on a real quadruped robot

Supported diverse gaits and transitions

Abstract

This paper combines the fast Zero-Moment-Point (ZMP) approaches that work well in practice with the broader range of capabilities of a Trajectory Optimization formulation, by optimizing over body motion, footholds and Center of Pressure simultaneously. We introduce a vertex-based representation of the support-area constraint, which can treat arbitrarily oriented point-, line-, and area-contacts uniformly. This generalization allows us to create motions such quadrupedal walking, trotting, bounding, pacing, combinations and transitions between these, limping, bipedal walking and push-recovery all with the same approach. This formulation constitutes a minimal representation of the physical laws (unilateral contact forces) and kinematic restrictions (range of motion) in legged locomotion, which allows us to generate various motion in less than a second. We demonstrate the feasibility of the generated motions on a real quadruped robot.

Figures

Tables

Table 1. TABLE I: Specs of the NLP for 16- and 4-step motions

	(16 steps, $1m$) $|$ (4 steps, $0.2m$)
	Walk	Trot	Pace	Bound
Horizon $T$ [s]	6.4 $|$ 1.6	2.4 $|$ 0.6	3.2 $|$ 0.8	3.2 $|$ 0.8
Variables [-]	646 $|$ 202	387 $|$ 162	1868 $|$ 728	1868 $|$ 728
Constraints [-]	850 $|$ 270	548 $|$ 255	2331 $|$ 939	2331 $|$ 939
$t_{k+1}-t_k$ [s]	0.1	0.05	0.02	0.02
Cost term	$J_{\lambda }$	-	-	-
Time Ipopt [s]	0.25 $|$ 0.06	0.02 $|$ 0.01	0.21 $|$ 0.12	0.17 $|$ 0.04
Time Snopt [s]	0.35 $|$ 0.04	0.04 $|$ 0.01	0.54 $|$ 0.18	0.42 $|$ 0.29

Equations

$$\left[\begin{array}{c}\dot{\mathbf{c}}\\ \ddot{\mathbf{c}}\end{array}\right]$$

$$\text{find }\mathbf{x}(t),\mathbf{u}(t)$$

$$\text{subject to}\mathbf{x}(0)-{\mathbf{x}}_0$$

$$\dot{\mathbf{x}}(t)-\mathbf{f}(\mathbf{x}(t),\mathbf{u}(t))$$

$$\mathbf{h}(\mathbf{x}(t),\mathbf{u}(t))$$

$$\mathbf{x}(T)-{\mathbf{x}}_T$$

$$\mathbf{x}(t),\mathbf{u}(t)$$

$$\mathbf{x}(t)$$

$$\mathbf{u}(t)$$

$$\mathbf{u}(t)={\mathbf{u}}_0={\mathbf{c}}_0+{\dot{\mathbf{c}}}_{0}h{g}^{-1}.$$

$${\mathbf{h}}_1(\mathbf{x}(t),\mathbf{u}(t))\le 0\iff \mathbf{u}\in \mathcal{P}({\mathbf{p}}^f,{\alpha }^f,c^f),$$

$${\mathbf{h}}_2(\mathbf{x}(t))\le 0\iff {\mathbf{p}}^f\in \mathcal{R}(\mathbf{c}).$$

$$\mathbf{w}=({\mathbf{w}}_c,{\mathbf{w}}_p,{\mathbf{w}}_u)$$

$$(\text{ref}eq:initia{l}_{c}onstraints),\text{(given initial state)}$$

$$(\text{ref}eq:co{m}_{c}ontinuity),(\text{ref}eq:co{m}_{a}cce{l}_{c}onstraint),(\text{ref}eq:kinemati{c}_{c}onstraints),(\text{ref}eq:convexit{y}_{a}ll)\text{(path-constraints)}$$

$$(\text{ref}eq:termina{l}_{c}onstraints),\text{(desired final state)}$$

$$\mathbf{w}=\mathrm{argmin}(\text{ref}eq:lambd{a}_{c}ost),\text{(robustness cost)}$$

$$\mathbf{x}(t)$$

$${\mathbf{w}}_c$$

$$\sum _{i=1}^4\left[\begin{array}{c}{T}_k\\ i\end{array}\right]{\mathbf{a}}_{k,i}{T}_k^{i-1}+\left[\begin{array}{c}{\mathbf{a}}_{k,0}\\ 0\end{array}\right]=\left[\begin{array}{c}{\mathbf{a}}_{k+1,0}\\ {\mathbf{a}}_{k+1,1}\end{array}\right].$$

$${\int }_{t_k}^{t_{k+1}}\ddot{\mathbf{c}}(t)\mathrm{d}t\approx {\int }_{t_k}^{t_{k+1}}{\mathbf{f}}_2(\mathbf{x}(t),\mathbf{u}(t))\mathrm{d}t.$$

$$\ddot{\mathbf{c}}[t]={\mathbf{f}}_2(\mathbf{x}[t],\mathbf{u}[t]),\forall t\in \left\{\begin{array}{c}{t}_k,\frac{t_{k+1}-t_k}{2},t_{k+1}\end{array}\right\}$$

$$\left[\begin{array}{c}{\mathbf{p}}^f(t)\\ {\alpha }^f(t)\end{array}\right]=\sum _{i=0}^3\left[\begin{array}{c}{\mathbf{a}}_{s,i}^f\\ b_{s,i}^f\end{array}\right](t-t_s{)}^i,$$

$$\begin{array}{rl}& {\mathbf{w}}_p=\left[\begin{array}{c}{\mathbf{w}}_p^1,\dots ,{\mathbf{w}}_p^{n_f}\end{array}\right],\\ \text{where }& {\mathbf{w}}_p^f=\left[\begin{array}{c}{\mathbf{p}}_1^f,{\alpha }_1^f,\dots ,{\mathbf{p}}_{n_s}^f,{\alpha }_{n_s}^f\end{array}\right]\end{array}$$

$$-{\mathbf{r}}^{x,y}<{\mathbf{p}}^f[t]-\mathbf{c}[t]-{\mathbf{p}}_{nom}^f<{\mathbf{r}}^{x,y}.$$

$$\begin{array}{rl}\mathbf{\lambda }(t)& ={\left[\begin{array}{c}{\mathbf{\lambda }}^1(t),\dots ,{\mathbf{\lambda }}^{n_f}(t)\end{array}\right]}^T,\\ \text{where }{\mathbf{\lambda }}^f(t)& =\left[\begin{array}{c}{\lambda }_1^f(t),\dots ,{\lambda }_{n_v}^f(t)\end{array}\right]\in {\left[0,1\right]}^{n_v}.\end{array}$$

$$\mathbf{u}(t)=\sum _{f=1}^{n_f}\sum _{v=1}^{n_v}{\lambda }_v^f(t)({\mathbf{p}}^f(t)+\mathbf{R}({\alpha }^f(t)){\mathbf{v}}_v),$$

$${\mathbf{w}}_u=\left[\begin{array}{c}{\mathbf{\lambda }}_1,\dots ,{\mathbf{\lambda }}_{n_u}\end{array}\right].$$

$${\Vert {\mathbf{\lambda }}_i\Vert }_1$$

$$0\le {\lambda }_v^f[t_i]$$

$${\lambda }_v^{f\ast }(t)=\frac{c^f(t)}{n_v(t)},$$

$$J_{\lambda }({\mathbf{w}}_u)=\sum _{i=1}^{n_u}{\Vert {\mathbf{\lambda }}_i-{\mathbf{\lambda }}_i^{\ast }\Vert }_2^2.$$

$${\mathbf{q}}_{b,ref}(t)={\left[\begin{array}{cccccc}0& 0& 0& c_x(t)& c_y(t)& h\end{array}\right]}^T.$$

$${\ddot{\mathbf{q}}}_{b,ref}={\ddot{\mathbf{q}}}_{b,ff}+K_p({\mathbf{q}}_b-{\mathbf{q}}_{b,ref})+K_d({\dot{\mathbf{q}}}_b-{\dot{\mathbf{q}}}_{b,ref}).$$

$${\ddot{\mathbf{q}}}_{j,ref}={\mathbf{J}}_j^{+}\left(\ddot{\mathbf{p}}-\dot{\mathbf{J}}\dot{\mathbf{q}}-{\mathbf{J}}_b{\ddot{\mathbf{q}}}_{b,ref}\right).$$

$$\mathbf{\tau }=(\mathbf{P}{\mathbf{S}}^T{)}^{+}\mathbf{P}(\mathbf{M}{\ddot{\mathbf{q}}}_{ref}+\mathbf{C}),$$

Full text

I Introduction

Planning and executing motions for legged systems is a complex task. A central difficulty is that legs cannot pull on the ground, e.g. the forces acting on the feet can only push upwards. Since the motion of the body is mostly generated by these constrained (=unilateral) contact forces, this motion is also restricted. When leaning forward past the tip of your tows, you will fall, since your feet cannot pull down to generate a momentum that counteracts the gravity acting on your Center of Mass (CoM). Finding motions that respect these physical laws can be tackled by various approaches described in the following.

A successful approach to tackle this problem is through full-body Trajectory Optimization (TO), in which an optimal body and endeffector motion plus the appropriate inputs are discovered to achieve a high-level goal. This was demonstrated by [1, 2, 3, 4, 5, 6, 7] resulting in an impressive range of motions for legged systems. These TO approaches have shown great performance, but are often time consuming to calculate and not straight-forward to apply on a real robot. In [8] the authors generate an wide range of quadruped gaits, transitions and jumps based on a parameterized controller and periodic motions. While the resulting motions are similar to ours, the methods are very different: While our approach is based on TO with physical constraints, [8] optimizes controller parameters based mainly on motion capture data.

Previous research has shown that to generate feasible motions to execute on legged systems, non-TO approaches also work well, although the motions cannot cover the range of the approaches above. One way is to model the robot as a Linear Inverted Pendulum (LIP) and keep the Zero-Moment-Point (ZMP) [9] inside the convex hull of the feet in stance. This approach has been successfully applied to generate motions for biped and quadruped walking [10, 11, 12, 13]. However, these hierarchical approaches use predefined footholds, usually provided by a higher-level planner beforehand that takes terrain information (height, slope) into account. Although this decoupling of foothold planning and body motion generation reduces complexity, it is unnatural, as the main intention of the footholds is to assist the body to achieve a desired motion. By providing fixed foot-trajectories that the body motion planner cannot modify, constraints such as stability or kinematic reachability become purely the responsibility of the lower-level body motion planner, artificially constraining the solution. A somewhat reverse view of the above are Capture Point (CP)[14] approaches, which have been successfully used to generate dynamic trotting and push recovery motions for quadruped robots [15, 16]. A desired body motion (usually a reference CoM velocity) is given by a high-level planner or heuristic, and a foothold/Center of Pressure (CoP) trajectory must be found that generates it.

Because of the dependency between footholds and body motion, approaches that optimize over both these quantities simultaneously, while still using a simplified dynamics model, have been developed [17, 18, 19, 20, 21, 22]. This reduces heuristics while increasing the range of achievable motions, but still keeps computation time short compared to full body TO approaches. These approaches are most closely related to the work presented in this paper.

The approaches [18, 19, 20, 21] demonstrate impressive performance on biped robots. One common difficulty in these approaches however is the nonlinearity of the CoP constraint with respect to the orientation of the feet. In [19, 20] the orientation is either fixed or solved with a separate optimizer beforehand. In [21] the nonlinearity of this constraint is accepted and the resulting nonlinear optimization problem solved. However, although the orientation of the individual feet can be optimized over in these approaches, a combined support-area with multiple feet in contact is often avoided, by not sampling the constraint during the double-support phase. For biped robots neglecting the constraint in the double-support phase is not so critical, as these take up little time during normal walking. For quadruped robots however, there are almost always two or more feet in contact at a given time, so the correct representation of the dynamic constraint in this phase is essential.

We therefore extend the capabilities of the approaches above by using a vertex-based representation of the CoP constraint, instead of hyperplanes. In [23] this idea is briefly touched, however the connection between the corners of the foot geometry and the convexity variables is not made and thereby the restriction of not sampling in the double-support phase remains. Through our proposed formulation, double- and single-stance support areas can be represented for arbitrary foot geometry, including point-feet. Additionally, it allows to represent arbitrarily oriented 1D-support lines, which wasn’t possible with the above approaches. Although not essential for biped walk on flat feet, this is a core necessity for dynamic quadruped motions (trot, pace, bound). This is a reason why ZMP-based approaches were so far only used for quadrupedal walking, where 2D-support areas were present.

The approach presented in this paper combines the LIP-based ZMP approaches that are fast and work well in practice with the broader range of capabilities of a TO formulation. A summary of the explicit contributions with respect to the papers above are:

• •

We reformulate the traditional ZMP-based legged locomotion problem [10] into a standard TO formulation with the CoP as input, clearly identifying state, dynamic model and path- and boundary-constraints, which permits easier comparison with existing methods in the TO domain. Push recovery behavior also naturally emerges from this formulation.

• •

We introduce a vertex-based representation of the CoP constraint, instead of hyperplanes. which allows us to treat arbitrarily oriented point-, line-, and area-contacts uniformly. This enables us to generate motions that are difficult for other ZMP-based approaches, such as bipedal walk with double-support phases, point-feet locomotion, various gaits as well as arbitrary combinations and transitions between these.

• •

Instead of the heuristic shrinking of support areas, we introduce a cost term for uncertainties that improve the robustness of the planned motions.

We demonstrate that the problem can be solved for multiple steps in less than a second to generate walking, trotting, bounding, pacing, combinations and transitions between these, limping, biped walking and push-recovery motions for a quadruped robot. Additionally, we verify the physically feasibility of the optimized motions through demonstration of walking and trotting on a real $80\text{\,}\mathrm{kg}$ hydraulic quadruped.

II Method

II-A Physical Model

We model the legged robot as a Linear Inverted Pendulum (LIP), with its CoM ${\mathbf{c}\!=\!(c_{x},c_{y})}$ located at a constant height $h$. The touchdown position of the pendulum with the ground (also known as ZMP or CoP) is given by ${\mathbf{u}\!=\!(u_{x},u_{y})}$ as seen in Fig. 1. The CoM acceleration $\ddot{\mathbf{c}}$ is predefined by the physics of a tipping pendulum

[TABLE]

The second-order dynamics are influenced by the CoM position $\mathbf{c}$, the CoP $\mathbf{u}$ and gravity $g$. This model can be used to describe a legged robot, since the robot can control the torques in the joints, thereby the contact forces and through these the position of the CoP. Looking only at the x-direction (left image in Fig. 1), if the robot decides to lift the hind leg, the model describing the system dynamics is a pendulum in contact with the ground at the front foot $\mathbf{p}^{RF}$, so ${\mathbf{u}\!=\!\mathbf{p}^{RF}\!=\!(p_{x},p_{y})}$. Since this pendulum is nearly upright, the CoM will barely accelerate in x; the robot is balancing on the front leg. However, lifting the front leg can be modeled as placing the pendulum at ${\mathbf{u}=\mathbf{p}^{LH}}$, which is strongly leaning and thereby must accelerate forward in x. By distributing the load between the legs, the robot can generate motions corresponding to a pendulum anchored anywhere between the contact points, e.g. $\mathbf{u}\in\mathcal{P}$ (see Fig. 1). Therefore, the CoP $\mathbf{u}$ is considered the input to the system and an abstraction of the joint torques and contact forces.

II-B Trajectory Optimization Problem

We want to obtain the inputs $\mathbf{u}(t)$ that generate a motion $\mathbf{x}(t)$ from an initial state $\mathbf{x}_{0}$ to a desired goal state $\mathbf{x}_{T}$ in time $T$ for a robot described by the system dynamics $\mathbf{f}(\mathbf{x},\mathbf{u})$, while respecting some constraints $\mathbf{h}(\mathbf{x},\mathbf{u})\leq 0$ and optimizing a performance criteria $J$. This can be formulated as a continuous-time TO problem

[TABLE]

The dynamics are modeled as those of a LIP (1), whereas the state and input for the legged system model are given by

[TABLE]

which includes the CoM position and velocity and the position and orientation of the ${n_{f}}$ feet. The input $\mathbf{u}(t)$ to move the system is the generated CoP, abstracting the usually used contact forces or joint torques.

II-C Specific Case: Capture Point

We briefly show that this general TO formulation, using the LIP model, also encompasses Capture Point methods to generate walking motions. Consider the problem of finding the position to step with a point-foot robot to recover from a push. With the initial position $\mathbf{c}_{0}$ and the initial velocity $\dot{\mathbf{c}}_{0}$ generated by the force of the push we have ${\mathbf{x}_{0}\!=\!(\mathbf{c}_{0},\dot{\mathbf{c}}_{0})}$. The robot should come to, and remain, at a stop at the end of the motion, irrespective of where and when, so we have ${\dot{\mathbf{c}}_{T\to\infty}\!=\!0}$. We parametrize the input by the constant parameter ${\mathbf{u}(t)\!=\!\mathbf{u}_{0}}$, as we only allow one step with a point-foot. We allow the CoP to be placed anywhere, e.g. no path constraints (2d) and do not have a preference as to how the robot achieves this task, e.g. ${J(\mathbf{x},\mathbf{u})\!=\!0}$.

Such a simple TO problem can be solved analytically, without resorting to a mathematical optimization solver (see Appendix -A). The point on the ground to generate and hold the CoP in order to achieve a final steady-state maintaining zero CoM velocity becomes

[TABLE]

This is the one-step Capture Point (CP), originally derived by [14] and the solution of our general TO formulation (2) for a very specific case (e.g one step/control input, zero final velocity).

II-D General Case: Legged Locomotion Formulation

Compared to the above example, our proposed formulation adds the capabilities to represent motions of multiple steps, time-varying CoP, physical restrictions as to where the CoP can be generated and preferences which of the feasible motions to choose. This TO formulation is explained on a high-level in the following, corresponding to Fig. 2, whereas more specific details of the implementation are postponed to the next section.

II-D1 Unilateral Forces

We clearly differentiate between the CoP $\mathbf{u}$ and the feet positions $\mathbf{p}^{f}$, which only coincide for a point-foot robot with one leg in contact. Generally, the footholds affect the input bounds of $\mathbf{u}$. We use $\mathbf{u}$ to control the body, but must at the same time choose appropriate footholds to respect the unilateral forces constraint. Traditional ZMP approaches fix the footholds $\mathbf{p}^{f}$ in advance, as the combination of both the CoP $\mathbf{u}$ and the footholds make this constraint nonlinear. We accept this nonlinearity and the higher numerical complexity associated with it. This gives us a much larger range of inputs $\mathbf{u}$, as we can “customize” our bounds $\mathcal{P}$ by modifying the footholds according to the desired task. Therefore, the first path constraint of our TO problem is given by

[TABLE]

where $\mathcal{P}$ represents the convex hull of the feet in contact as seen in Fig. 2 and $c^{f}\in\{0,1\}\in\mathbb{Z}$ is the indicator if foot $f$ is in contact.

We implement this convex hull constraint by weighing the vertices/corners of each foot in contact. This extends the capabilities of traditional representations by line segments/hyperplanes to also model point- and line-contacts of arbitrary orientation. We use predefined contact sequences and timings $c^{f}(t)$, to only optimize over real-valued decision variables $w\in$ and not turn the problem into a mixed-integer Nonlinear Programming Problem (NLP). Simply by adapting this contact schedule $c^{f}(t)$, the optimizer generates various gaits as well as combinations and transitions between these, for which previously separate frameworks were necessary.

II-D2 Kinematic Reachability

When modifying the footholds to enclose the CoP, we must additionally ensure that these stay inside the kinematic range $\mathcal{R}$ of the legs (Reachability). This constraint that depends on both the CoM $\mathbf{c}$ and foothold positions $\mathbf{p}^{f}$ is formulated for every leg $f$ as

[TABLE]

Allowing the modification of both these quantities simultaneously characterizes the legged locomotion problem more accurately and reduces heuristics used in hierarchical approaches.

II-D3 Robust Motions

With the above constraints the motion will comply to physics and the kinematics of the system. This feasible motion is assuming a simplified model, a perfect tracking controller and an accurate initial state. To make solutions robust to real world discrepancies where these assumptions are violated, it is best to avoid the borders of feasible solutions, where the inequality constraints are tight (${\mathbf{h}\!=\!\mathbf{0}}$). This can be achieved by artificially shrinking the solutions space by a stability margin (e.g. ${\mathbf{h}\!\leq\!\mathbf{m}}$). For legged locomotion this is often done by shrinking the support area to avoid solutions were the CoP is placed at the marginally-stable border [11].

We do not restrict the solution space, but choose the more conservative of the feasible motions through a performance criteria $J_{\lambda}$. This soft constraint expresses “avoid boundaries when possible, but permit if necessary”. The robot is allowed to be at marginally stable states, but since there are many uncertainties in our model and assumptions, it is safer to avoid them. This cost does not require a hand-tuned stability margin and the solution can still be at the boundaries when necessary. However, especially for slow motions (e.g. walking) where small inaccuracies can accumulate and cause the robot to fall, this cost term is essential to generate robust motions for real systems.

III Implementation

There exist different methods to solve Optimal Control problems (2), namely Dynamic Programming (Bellman Optimality Equation), indirect (Maximum Principle) and direct methods [24]. In direct methods the continuous time TO problem is represented by a finite number of decision variables and constraints and solved by a nonlinear programming solver. If the decision variables $\mathbf{w}$ fully describe the input $\mathbf{u}(t)$ and state $\mathbf{x}(t)$ over time, the method is further classified as a simultaneous direct method, with flavors Direct Transcription and Multiple Shooting. In our approach we chose a Direct Transcription formulation, e.g. optimizing state and controls together. This has the advantage of not requiring an ODE solver, constraints on the state can be directly formulated and the sparse structure of the Jacobian often improves convergence. The resulting discrete formulation to solve the continuous problem in (2) is given by

[TABLE]

where $\mathbf{w}_{c}$ are the parameters describing the CoM motion, $\mathbf{w}_{p}$ the feet motion (swing and stance) and $\mathbf{w}_{u}$ the position of the CoP. This section describes in detail how we parametrize the state $(\mathbf{w}_{c},\mathbf{w}_{p})$ and input $\mathbf{w}_{u}$, formulate the path-constraints and defined the cost (21).

III-A Center-of-Mass Motion

This section explains how the continuous motion of the CoM can be described by a finite number of variables to optimize over, while ensuring compliance with the LIP dynamics.

III-A1 CoM Parametrization

The CoM motion is described by a spline, strung together by $n$ quartic-polynomials as

[TABLE]

with coefficients $\mathbf{a}_{k,i}\in 2$ and $t_{k}$ describing the global time at the start of polynomial $k$.

We ensure continuity of the spline by imposing equal position and velocity at each of the $n\!-\!1$ junctions between polynomial $k$ and $k\!+\!1$, so $\mathbf{x}[t_{k+1}^{-}]=\mathbf{x}[t_{k+1}^{+}]$. Using $T_{k}=t_{k+1}\!-\!t_{k}$ we enforce

[TABLE]

III-A2 Dynamic Constraint

In order to ensure consistency between the parametrized motion and the dynamics of the system (1), the integration of our approximate solution $\ddot{\mathbf{c}}(t)$ must resemble that of the actual system dynamics, so

[TABLE]

Simpson’s rule states that if $\ddot{\mathbf{c}}(t)$ is chosen as a $2^{nd}$-order polynomial (which is why $\mathbf{c}(t)$ is chosen as $4^{th}$-order) that matches the system dynamics $\mathbf{f}_{2}$ at the beginning, the center and at the end, then (11) is bounded by an error proportional to $(t_{k+1}-t_{k})^{4}$. Therefore we add the following constraints for each polynomial

[TABLE]

(see Appendix -B for a more detailed formulation). By keeping the duration of each polynomial short (\sim$$50\text{\,}\mathrm{m}\mathrm{s}), the error of Simpson’s integration stays small and the $4^{th}$-order polynomial solution $\mathbf{c}(t)$ is close to an actual solution of the Ordinar Differential Equation (ODE) in (1).

This formulation is similar to the ”collocation” constraint [25]. Collocation implicitly enforces the constraints (12) at the boundaries through a specific parametrization of the polynomial, while the above formulation achieves this through explicit constraints in the NLP. Reversely, collocation enforces that $\frac{\partial{\mathbf{c}(t)}}{\partial{t}}=\dot{\mathbf{c}}(t)$ through the explicit constraint, while our formulation does this through parametrization in (8).

III-B Feet Motion

III-B1 Feet Parametrization

We impose a constant position $\mathbf{p}^{f}_{s}$ and orientation $\alpha^{f}_{s}$ if leg $f$ is in stance. We use a cubic polynomial in the ground plane to move the feet between two consecutive contacts

[TABLE]

where $(t\!-t_{s})$ is the elapsed time since the beginning of the swing motion. The vertical swingleg motion does not affect the NLP and is therefore not modeled. The coefficients $\mathbf{a}_{s,i}\in 2$ and $b_{s,i}\in$ are fully determined by the predefined swing duration and the position and orientation of the enclosing contacts $\begin{Bmatrix}\mathbf{p}^{f}_{s},\alpha^{f}_{s}\end{Bmatrix}$ and $\begin{Bmatrix}\mathbf{p}^{f}_{s+1},\alpha^{f}_{s+1}\end{Bmatrix}$. Therefore the continuous motion of all ${n_{f}}$ feet can be parametrized by the NLP decision variables

[TABLE]

are the parameters to fully describe the motion of a single leg $f$ taking ${n_{s}}$ steps.

III-B2 Range-of-Motion Constraint

To ensure a feasible kinematic motion, we must enforce $\mathbf{p}^{f}\in\mathcal{R}(\mathbf{c})$, which is the gray area in Fig. 3. We approximate the area reachable by each foot through a rectangle $[-\mathbf{r}^{x,\!y},\mathbf{r}^{x,\!y}]$, representing the allowed distance that a foot can move from its nominal position $\mathbf{p}^{f}_{nom}$ (center of gray area). The foothold position for each foot $f$ is therefore constrained by

[TABLE]

Contrary to hierarchical approaches, this constraint allows the optimizer to either move the body to respect kinematic limits or place the feet at different positions. A constraint on the foot orientation can be formulated equivalently.

III-C Center of Pressure Motion

To represent the continuous CoP trajectory, we parameterize it through the load carried by each endeffector. This parametrization is used to formulate a novel convexity constraint based on vertices instead of hyperplanes. Finally this section introduces a cost that keeps the CoP from marginally stable regions and improves robustness of the motion.

III-C1 CoP Parametrization

The CoP $\mathbf{u}(t)$ is not parametrized by polynomial coefficients or discrete points, but by the relative load each corner of each foot is carrying. This load is given by

[TABLE]

${n_{v}}$ represents the number of vertices/corners of foot $f$. For the square foot in Fig. 3, four lambda values represent one foot and distribute the load amongst the corners. These multipliers represent the percentage of vertical force that each foot is carrying, e.g. $|\!|\boldsymbol{\lambda}^{f}(t)|\!|_{1}=0.9$ implies that leg $f$ is carrying 90% of the weight of the robot at time $t$. Using these values, the CoP is parameterized by

[TABLE]

where $\mathbf{R}(\alpha^{f})\in\mathbb{R}^{2\times 2}$ represents the rotation matrix corresponding to the optimized rotation $\alpha^{f}$ of foot $f$ (13). $\mathbf{v}_{v}$ represents the fixed position (depending on the foot geometry) of corner $v$ of the foot expressed in the foot frame. For a point-foot robot with ${\mathbf{v}_{v}\!=\!\mathbf{0}}$, (17) simplifies to $\mathbf{u}=\sum_{f=1}^{n_{f}}\lambda^{f}\mathbf{p}^{f}$.

We represent $\boldsymbol{\lambda}(t)$ for the duration of the motion by piecewise-constant values $\boldsymbol{\lambda}_{i}=\boldsymbol{\lambda}(t_{i})$ discretized every $20\text{\,}\mathrm{m}\mathrm{s}$, resulting in $n_{u}$ nodes. Therefore the CoP $\mathbf{u}$ can be fully parameterized by $\mathbf{w}_{p}$ and the additional NLP decision variables

[TABLE]

III-C2 Unilateral Forces Constraint

We represent the essential input constraint (6), which ensures that only physically feasible forces inside the convex hull of the contacts are generated, for $i=1,\dots,n_{u}$ as

[TABLE]

where $c^{f}\,{\in}\,\{0,1\}\,{\in}\,\mathbb{Z}$ is the indicator if foot $f$ is in contact. The constraints (17) and (19a) allow $\mathbf{u}$ to be located anywhere inside the convex hull of the vertices of the current foot positions, independent of whether they are in contact. However, since only feet in contact can actually carry load, (19b) enforces that a leg that is swinging (${c^{f}\!=\!0}$) must have all the corners of its foot unloaded. These constraints together ensure that the CoP lies inside the green area shown in Fig. 3.

III-C3 Robust Walking Cost

To keep the CoP away from the edges of the support-area we could constrain $\lambda_{v}^{f}$ of each leg in stance to be greater than a threshold, causing these legs in contact to never be unloaded. This conceptually corresponds to previous approaches that heuristically shrink support areas and thereby reduce the solution-space for all situations. We propose a cost that has similar effect, but still permits the solver to use the limits of the space if necessary.

The most robust state to be in, is when the weight of the robot is equally distributed amongst all the corners in contact, so

[TABLE]

where $n_{v}(t)={n_{v}}\sum_{f=1}^{{n_{f}}}c^{f}(t)$ is the total number of vertices in contact at time $t$, predefined by the contact sequence $c(t)$. This results in the CoP to be located in the center of the support areas. The deviation of the input values from the optimal values $\lambda^{*}$ over the entire discretized trajectory (18) is then given by

[TABLE]

For a support triangle (${\lambda^{f*}_{v}\!=\!\frac{1}{3}}$) this cost tries to keep the CoP in the center and for a line (${\lambda^{f*}_{v}\!=\!\frac{1}{2}}$) in the middle. For quadruped walking motions this formulation generates a smooth transition of the CoP between diagonally opposite swing-legs, while still staying away from the edges of support-areas whenever possible.

IV Tracking the Motion

The motion optimization part of our approach is largely robot independent. The only robot specific information needed to run the framework are the robot height, the number of feet, their geometry and their kinematic range. For execution however, the optimized motion must be translated into joint torques $\boldsymbol{\tau}$ using a fully-body dynamics model. This section discusses this generation summarized by Fig. 4.

IV-A Generating full-body reference accelerations

The 6–Degrees of Freedom (DoF) base pose is reconstructed using zero desired orientation (in Euler angles x,y,z), the optimized CoM motion $\mathbf{c}$ (assuming the geometric center of the base coincides with the CoM), and the constant base height $h$ as

[TABLE]

In order to cope with uncertainties it is essential to incorporate feedback into the control loop. We do this by adding an operational space PD-controller on the base that creates desired 6D base accelerations according to

[TABLE]

The derivate of the pose, the base twist $\dot{\mathbf{q}}_{b}\in\mathbb{R}^{6}$ represents the base angular and linear velocities and $\ddot{\mathbf{q}}_{b,f\!f}$ is the optimized CoM acceleration from the NLP. This controller modifies the planned body motion if the current state deviation from the reference state.

In order to obtain the desired joint accelerations that correspond to the planned Cartesian motion of the feet we can use the relationship $\ddot{\mathbf{p}}(t)=\mathbf{J}\ddot{\mathbf{q}}+\dot{\mathbf{J}}\dot{\mathbf{q}},$ where $\dot{\mathbf{q}},\ddot{\mathbf{q}}\in 6+n$ represent the full body state (base + joints) and $\mathbf{J}=\begin{bmatrix}\mathbf{J}_{b}&\mathbf{J}_{j}\end{bmatrix}\in\mathbb{R}^{3{n_{f}}\times(6+n)}$ the Jacobian that maps full-body velocities to linear foot velocities in world frame. Rearranging this equation, and using the Moore–Penrose pseudoinverse $\mathbf{J}_{j}^{+}$, gives us the reference joint acceleration

[TABLE]

IV-B Inverse Dynamics

The inverse dynamics controller is responsible for generating required joint torques $\boldsymbol{\tau}$ to track the reference acceleration $\ddot{\mathbf{q}}_{ref}$, which is physically feasible based on the LIP model. This is done based on the rigid body dynamics model of the system, which depends on the joint torques, but also the unknown contact forces. To eliminate the contact forces from the equation, we project it into the space of joint torques by $\mathbf{P}=\mathbf{I}-\mathbf{J}_{c}^{+}\mathbf{J}_{c},$ where $\mathbf{J}^{T}_{c}$ is the contact Jacobian that maps Cartesian contact forces to joint torques [26, 27]. This allows us to solve for the required joint torques through

[TABLE]

where $\mathbf{M}$ is the joint space inertia matrix, $\mathbf{C}$ the effect of Coriolis forces on the joint torques and $\mathbf{S}$ the selection matrix which prohibits from actuating the floating base state directly. We found it beneficial to also add a low-gain PD-controller on the joint position and velocities. This can mitigate the effects of dynamic modeling errors and force tracking imperfections.

V Results

We demonstrate the performance of this approach on the hydraulically actuated quadruped robot HyQ [28]. The robot weighs approximately $80\text{\,}\mathrm{kg}$, moves at a height of about $0.6\text{\,}\mathrm{m}$ and is torque controlled. Base estimation [29] is performed on-board, fusing Inertial Measurement Unit (IMU) and joint encoder values. Torque tracking is performed at $1000\text{\,}\mathrm{Hz}$, while the reference position, velocity and torque set-points are provided at $250\text{\,}\mathrm{Hz}$. The C++ dynamics model is generated by [30].

V-A Discussion of generated motions

This section analyses the different motions generated by changing the sequence and timings of contacts $c(t)$. There is no high-level footstep planner; the footholds are chosen by the optimizer to enable the body to reach a user defined goal state $\mathbf{x}_{T}$. The results where obtained using C++ code interfaced with Interior Point Method (Ipopt [31]) or Sequential Quadratic Programming (Snopt [32]) solvers on an Intel Core i7/$2.8\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$ Quadcore laptop. The Jacobians of the constraint and the gradient of the cost function are provided to the solver analytically, which is important for performance. We initialize the decision variables $\mathbf{w}$ with the quadruped standing in default stance for a given duration. The shown motions correspond to the first columns (e.g. 16 steps) in Table I. The reader is encouraged to view the video at https://youtu.be/5WLeQMBuv30, as it very intuitively demonstrates the performance of this approach. Apart from the basic gaits, the video shows the capability of the framework to generate gradual transitions between them, bipedal walking, limping and push-recovery.

V-A1 Walk

Fig. 5(a) shows a walk of multiple steps, with the two support areas highlighted for swinging RF$\to$LH. The effect of the cost term $J_{\lambda}$ is clearly visible, as the CoP is accumulated away from the support area borders by left-right swaying of the body. Only when switching diagonally opposite legs the CoP lies briefly at the marginally stable border, but then immediately shifts to a more conservative location. Without the cost term, the CoM motion is a straight line between $\mathbf{x}_{0}$ and $\mathbf{x}_{T}$, causing the real system to fail.

V-A2 Trot

Fig. 5(b) shows a completely different pattern of support areas and CoP distribution. During trotting only line-contacts exist, so the possible places to generate the CoP is extremely restricted compared to walking. Notice how the CoP lies close to the CoM trajectory during the middle of the motion, but deviates quite large back/forward during the start/end of the motion (e.g. the robot pushing off from the right-front (green) leg in the second to last step). This is because the distance between the CoP and the CoM generates the acceleration necessary for starting and stopping, whereas in the middle the robot is moving with nearly constant velocity.

V-A3 Pace/Bound/Biped Walk

Specifying legs on the same side to be in contact, with a short four-leg transition period between them produces the motion shown in Fig. 5(c). This can also be viewed as biped walking with line-feet (e.g. skis), with the constraint enforced also during the double-stance phase. The first observation is the sideways swaying motion of the CoM. This is necessary because the support areas do not intersect (as in the trot) the CoM trajectory. Since the CoP always lies inside these left and right support areas, they will accelerate the body away from that side until the next step, which then reverses the motion. We found that the LIP model with fixed zero body orientation does not describe such a motion very well, as the inherent rotation (rolling) of the body is not taken into account. In order to also demonstrate these motions on hardware, the LIP model must be extended by the angular body motion. Specifying the front and hind legs to alternate between contact generates a bound Fig. 5(d). The lateral shifting motion of the pace is now transformed to a forward backward motion of the CoM due to support areas. In case of an omni directional robot a bounding gait can simply be considered a side-ways pace.

VI Conclusion

This paper presented a TO formulation using vertex-based support-area constraints, which enables the generation of a variety of motions for which previously separate methods were necessary. In the future, more decision variables (e.g. contact schedule, body orientation, foothold height for uneven terrain), constraints (e.g. friction cone, obstacles) and more sophisticated dynamic models can be incorporated into this formulation. Additionally, we plan to utilize the speed of the optimization for Model Predictive Control (MPC).

-A Derivation of Capture Point

Consider the differential equation describing a LIP (linear, constant coefficients, second order) in x-direction

[TABLE]

The general solution to the homogeneous part of the equation can be construct by the Ansatz $c(t)=e^{\alpha t}$ which leads to the characteristic equation $\alpha^{2}e^{\alpha t}-\frac{g}{h}e^{\alpha t}=0,$ resulting in $\alpha=\pm\sqrt{\frac{g}{h}}.$ Assuming constant input $u_{0}$ leads to the partial solution $c_{p}(t)=u_{0}$, and the space of solutions for the entire ODE is given by

[TABLE]

where $\beta_{1},\beta_{2}\in\mathbb{R}$ are the free parameters describing the motion. Imposing the initial position $c(0)=\beta_{1}+\beta_{2}+u_{0}\overset{!}{=}c_{0}$ and velocity $\dot{c}(0)=\alpha\beta_{1}-\alpha\beta_{2}\overset{!}{=}\dot{c}_{0}$ we obtain

[TABLE]

As $t\to\infty$ we require the velocity to remain at zero (pendulum at rest). Since $\lim_{t\to\infty}e^{-\alpha t}=0$, and $\alpha\neq 0$ we have

[TABLE]

which is known as the one-step Capture Point originally derived in [14].

-B Dynamic Constraint

The system dynamics constraint (11) enforced through $\ddot{\mathbf{c}}[t]=\mathbf{f}_{2}(\mathbf{x}[t],\mathbf{u}[t])$, with the local polynomial time {\underaccent{\bar}{t}\!=\!(t-t_{k})}, are formulated as

[TABLE]