ORIGINAL RESEARCH article
Sec. Field Robotics
From Walking to Running: 3D Humanoid Gait Generation via MPC
- Dipartimento di Ingegneria Informatica, Automatica e Gestionale, Sapienza University of Rome, Rome, Italy
We present a real time algorithm for humanoid 3D walking and/or running based on a Model Predictive Control (MPC) approach. The objective is to generate a stable gait that replicates a footstep plan as closely as possible, that is, a sequence of candidate footstep positions and orientations with associated timings. For each footstep, the plan also specifies an associated reference height for the Center of Mass (CoM) and whether the robot should reach the footstep by walking or running. The scheme makes use of the Variable-Height Inverted Pendulum (VH-IP) as a prediction model, generating in real time both a CoM trajectory and adapted footsteps. The VH-IP model relates the position of the CoM to that of the Zero Moment Point (ZMP); to avoid falling, the ZMP must be inside a properly defined support region (a 3D extension of the 2D support polygon) whenever the robot is in contact with the ground. The nonlinearity of the VH-IP is handled by splitting the gait generation into two consecutive stages, both requiring to solve a quadratic program. Thanks to a particular triangular structure of the VH-IP dynamics, the first stage deals with the vertical dynamics using the Ground Reaction Force (GRF) as a decision variable. Using the prediction given by the first stage, the horizontal dynamics become linear time-varying. During the flight phases, the VH-IP collapses to a free-falling mass model. The proposed formulation incorporates constraints in order to maintain physically meaningful values of the GRF, keep the ZMP in the support region during contact phases, and ensure that the adapted footsteps are kinematically realizable. Most importantly, a stability constraint is enforced on the time-varying horizontal dynamics to guarantee a bounded evolution of the CoM with respect to the ZMP. Furthermore, we show how to extend the technique in order to perform running on tilted surfaces. We also describe a simple technique that receives input high-level velocity commands and generates a footstep plan in the form required by the proposed MPC scheme. The algorithm is validated via dynamic simulations on the full-scale humanoid robot HRP-4, as well as experiments on the small-sized robot OP3.
Humanoid robots have become more and more advanced in the last decade. Using legged locomotion, they are in principle able to perform very complex and dynamic motions, such as those needed to traverse nonlevel terrain and realize running gaits. However, the design of online walking and running pattern generation algorithms is challenging from a control standpoint, as there are several hard requirements such as the need to always maintain dynamic balance, whether the robot is walking or running.
When walking on flat ground, balance is usually guaranteed by keeping the Zero Moment Point1 (ZMP) within the convex hull of the contact surfaces, that is, the support polygon. To achieve real time control, a popular approach is to reduce the complexity of humanoid dynamics using the Linear Inverted Pendulum (LIP) model, obtained by assuming constant CoM height and neglecting angular momentum variations around the CoM. The linearity of the LIP model has allowed to efficiently define control schemes based on optimization, such as preview-based approaches for tracking the desired ZMP trajectory (Kajita et al., 2003), or MPC schemes that encode the balance requirement through constraints on the ZMP (Wieber, 2006).
1.1 Related Work
Variation of the CoM height is a key requirement for traversing nonlevel grounds. Removing the assumption of constant CoM height in the LIP leads to the Variable-Height Inverted Pendulum (VH-IP) model (Koolen et al., 2016), in which height variations have the effect of modifying the natural frequency of the pendulum. The varying natural frequency can be seen as an additional input of the model, which introduces a nonlinearity that entails higher complexity in the gait generation algorithm. Another possibility is to consider it as a time-varying parameter in the dynamic equation (Herdt et al., 2012), design the vertical motion offline, and generate the horizontal trajectories with a linear time-varying MPC. Other approaches (Hu et al., 2014; Brasseur et al., 2015) aim instead at bounding the effect of the difference between the VH-IP and LIP on the ZMP, in a way to guarantee that the latter is always inside the support polygon. It is also possible to find specific trajectories for which the vertical CoM motion has LIP-like dynamics (Englsberger et al., 2013; Zamparelli et al., 2018). This approach can be very effective but it constrains the CoM trajectories that can be generated.
Regardless of the employed model, gait generation schemes must necessarily deal with the intrinsic instability of humanoid dynamics. In the LIP, this instability is reflected in the presence of an unstable mode. This component is often called Divergent Component of Motion (DCM) (Takenaka et al., 2009a), or Capture Point (Pratt et al., 2006), and plays a major role in the stabilization of gait generation schemes. In the effort to extend this concept to the VH-IP, if one follows the same approach as for the LIP, it is hard to enforce stability on the DCM itself as it will remain coupled with the CoM dynamics (Hopkins et al., 2014). However, a truly decoupling change of coordinates can be found, (Caron et al., 2019a), but it is a time-varying change of coordinates with no available closed-form expression, as it is implicitly defined by solving a nonlinear differential equation.
A separate line of research has been directed at the generation of running motions, with early work dating back several decades (Raibert, 1986). To generate a running gait, the controller must be able to produce a variable CoM height, but it also needs to account for flight phases, in which contact with the ground is lost. As such, none of the aforementioned control schemes for 3D walking can be directly adopted to generate a running gait. Biomechanic studies on human running (Blickhan, 1989; Novacheck, 1998) have consolidated the use of the Spring Loaded Inverted Pendulum (SLIP) to model running dynamics, which has frequently been involved in the generation of running motions for legged robots (Wensing and Orin, 2013; Secer and Saranli, 2018; Secer and Cinar, 2020). In these works, control is realized by means of a step-to-step regulation of apex states (flight states with zero vertical velocity). Xiong and Ames, (2018) derive a spring-mass model from the robot dynamics to generate hopping motions by optimizing the second-order derivative of the leg length. Winkler et al. (2018) perform offline trajectory optimization to compute dynamic motions, involving flight phases and locomotion over tilted planes. The method makes use of a convenient polynomial parameterization to efficiently optimize phase durations, foot locations, and contact forces while enforcing dynamic consistency through Centroidal Dynamics. More recently, Li et al. (2021) propose an optimization-based CoM trajectory planning using a SLIP model during stance phases for quadrupedal hopping on two legs. SLIP dynamics are highly nonlinear and thus unsuitable for predictive control unless considerably approximated.
Although the ZMP is not defined during flight phases due to the loss of contact with the ground, balance during support phases can still be enforced by means of a ZMP criterion. ZMP-based approaches for running are proposed with VH-IP modeling of the dynamics during support phases and a free-falling mass model during flight phases (Kajita et al., 2007; Takenaka et al., 2009b; Tajima et al., 2009). The vertical motion is generated independently, the ZMP trajectory is preassigned and the CoM trajectory is obtained by numerically solving the dynamics in order to meet the next step boundary conditions. Englsberger et al. (2016) generate running gaits aimed at approximating natural human behavior, using measured data and a spatially independent polynomial encoding of the CoM trajectory. More recently, Ahn and Cho, (2021) proposed an MPC formulation for executing jumping motions where the vertical trajectory is generated offline. Also, Sugihara et al. (2021) introduced another ZMP-based MPC approach to generate 3D walking and running; constraints are enforced approximately in order to obtain closed-form solutions to speed up the execution.
We recently proposed an Intrinsically Stable MPC (IS-MPC) with an explicit stability constraint embedded in the formulation (Scianca et al., 2020). This constraint is in charge of dealing with the unstable character of humanoid dynamics, by guaranteeing that the CoM trajectory is always bounded with respect to the ZMP. Since the LIP is used as a prediction model, the DCM dynamics can be decoupled by using a simple change of coordinates, and a closed-form condition on the long-term evolution can be found. Zamparelli et al. (2018) used similar ideas and proposed an IS-MPC extension that is capable of generating 3D trajectories for the CoM. This can be realized by constraining the allowed CoM evolution to a specific set of trajectories that present LIP-like dynamics along all three axes. While this approach is very effective in many circumstances, it limits the number of trajectories that can be generated and it is thus not applicable to every situation. Furthermore, it cannot be extended to generate running motions.
In this article, we present a unified MPC algorithm that generates online both 3D walking and running gaits as well as transitions between them. Motions are generated in such a way to follow as closely as possible a footstep plan, in an environment constituted by horizontal patches at different heights, known as world of stairs.
The main contributions of this article can be summarized by the following points:
• we describe a scheme for generating a humanoid gait with a variable CoM height using the VH-IP model. The proposed architecture uses a two-stage MPC scheme and only requires solving standard linear-quadratic QP problems by means of the fact that the natural frequency present in the horizontal dynamics can be regarded as a time-varying parameter. With respect to the recent state of the art methods, our scheme has the advantage of maintaining a linear formulation without limiting the possible vertical CoM motions by constraining the pendulum frequency (Zamparelli et al., 2018; Sugihara et al., 2021);
• by properly shaping the vertical GRF using constraints in the first MPC stage, we introduce flight phases in order to allow the scheme to generate running motions;
• we propose a variation of our stability constraint (Scianca et al., 2020) to extend its applicability to the case of linear time-varying systems, and discuss how it is included in the formulation, both in the case of walking with variable height and running. Differently from Caron et al. (2019a), the proposed constraint is linear and can be efficiently enforced in a QP;
• we show how to extend the proposed scheme in order to allow running over tilted surfaces.
The proposed approach allows to generate gaits involving stair climbing and running for traversing complex terrain configurations. Furthermore, thanks to the ability to generate a CoM trajectory with variable height, it can realize a transition from walking to crouching (Kamioka et al., 2015), for example, in order to pass below an overhanging obstacle. The capabilities of the scheme will be demonstrated in dynamic simulations, as well as experiments for the case of walking. As a supplementary contribution, we will also present a simple footstep planner that, in response to high-level velocity commands, can be used to generate the required input for the proposed scheme.
The article is organized as follows: Section 2 describes the dynamics of balance on nonlevel ground. The VH-IP model is introduced in Section 3. The overall approach is sketched out in Section 4, and then the MPC is presented in detail in Section 5, while a simple footstep planner is discussed in Section 6. Simulations, as well as an experiment on the ROBOTIS OP3 humanoid robot, are shown in Section 7 and Section 8. In Section 9 we will briefly discuss an extension that allows the robot to run on tilted surfaces. Section 10 concludes the article.
2 Dynamic Balance
The ability to maintain dynamic balance is a basic requirement of any gait generation scheme. On horizontal ground, balance may be guaranteed using the ZMP criterion, which prescribes that the ZMP (the point where the horizontal component of the moment of the ground reaction forces becomes zero) should be at all times within the support polygon, that is, the convex hull of the contact surfaces.
When moving to 3D locomotion, where the contact surfaces are not necessarily coplanar, a support polygon cannot be defined. At the same time, one needs to consider that the definition of ZMP is actually satisfied along a line, and therefore this point is not necessarily located on the ground (Sardain and Bessonnet, 2004). Based on this, several authors have proposed extensions of the ZMP criterion to the 3D case; in particular, Sugihara et al. (2021) require that the ZMP belongs to a pyramid with vertex at the CoM and base defined as the convex hull of the projections of the active contact surfaces on an arbitrary plane, typically chosen as a horizontal plane below all contact surfaces.
In this article, we are going to consider as admissible ZMP region for balance the convex hull of the active contact surfaces:
where p is2 the generic 3D point and vj, j = 1, … , Nc, are the vertexes of all active contact surfaces. It is easy to verify that
FIGURE 1. Left: In a world of stairs, where all contact surfaces are parallel to the ground, the admissible ZMP region
3 System Modeling
In this article, we will use the VH-IP model (Caron et al., 2019a) to characterize the relationship between the CoM and the ZMP.
Under the assumption of zero angular momentum around the CoM, the Newton and Euler (with respect to the ZMP) equations for the humanoid lead to
where m is the robot mass, pc = (xc, yc, zc) is the CoM, g = (0, 0, − g) is the gravity acceleration, f is the ground reaction force (GRF in the following) and pz = (xz, yz, zz) is the ZMP (see Figure 2).
where fz is the vertical component of the GRF. The horizontal dynamics are instead obtained as:
represents the (variable) squared natural frequency of the pendulum.
Together, Eqs 3–5 constitute the dynamic model of our system. We consider as control inputs the horizontal ZMP components xz, yz, and the vertical GRF fz, while the profile of the vertical ZMP component zz is assumed to be uniquely determined3 once the footstep plan is specified (see Section 5.1).
Note the triangular structure of our model: the vertical dynamics (3) are independent of the horizontal dynamics (4–5). Once fz is chosen, the evolution of zc is determined, and the same is true for λ; at this point, the horizontal dynamics can be interpreted as a linear time-varying system. We will exploit this fact when designing our MPC controller.
The assumption that fz is a control input requires that at least one contact with the ground is active. During flight phases, which are an integral part of a running gait, the GRF is identically zero. Setting f = 0 in (1) gives
indicating that the motion of the CoM is entirely controlled by gravity. Note that Eq. 7 coincides with (3–5) when fz = 0.
4 Proposed Approach
For a humanoid robot moving in a world of stairs, we wish to design a control scheme for replicating as closely as possible an assigned footstep plan. In general, this plan will consist of a sequence of candidate footstep poses (3D positions and orientations) with specified timings. In addition, the plan also specifies for each step a reference value for the CoM height relative to the ZMP, and whether the step itself should be performed by walking or running. A simple algorithm to compute such a footstep plan will be discussed in Section 6.
A block scheme of the proposed approach is shown in Figure 3. The footstep plan is the input to the IS-MPC (Intrinsically Stable MPC) block, which consists of two sequential stages. The first is a quadratic program (QP-z) where fz in (3) is chosen over a certain control horizon to generate a vertical CoM trajectory which tracks the reference height as closely as possible. As a byproduct, the corresponding values of λ can be computed via (6) thanks to the triangular structure of the VH-IP model. At this point, it is possible to enter the second stage, which is another quadratic program (QP-xy) over the same control horizon, where xz, and yz in (4–5) are chosen and the footsteps are adapted to generate a horizontal CoM trajectory xc, yc resulting in a gait which is both dynamically balanced and internally stable. A qualitative representation of the trajectories produced by each stage is shown in Figure 4, with reference to a typical example.
FIGURE 4. Left: a qualitative representation of the final trajectories after QP-z is solved in an example case. Here, the reference CoM/ZMP distance is kept constant for the initial 3 steps and then lowered in order to pass below an obstacle. Right: a qualitative representation of the trajectories generated after QP-xy is solved.
The generated 3D CoM trajectory, along with an appropriate trajectory of the swing foot4, is sent to a kinematic control block where joint commands are computed using standard pseudoinversion techniques.
5 IS-MPC for 3D Walking and Running
In this section, we will describe in detail the IS-MPC scheme for 3D walking and running based on the VH-IP model of Section 3.
The scheme works in discrete time over time intervals of duration δ. A sampled variable at tk will be denoted by the k superscript, for example,
To simplify the exposition, we will assume that all footsteps in the plan have the same orientation, taken without loss of generality to be the same of the x axis. This has the effect of decoupling the motion along the sagittal and coronal plane, making it possible to present the method for the x component only, with the understanding that a similar reasoning holds for the y component, except when specified. However, it is relatively easy to extend the proposed method to the case of variable orientation. When the orientation of the footsteps is not constant, the x and y components are coupled by a rotation matrix. As a consequence, a single QP problem involving both x and y components of the ZMP must be formulated, which however remains linear-quadratic provided that the orientation is not a decision variable in the QP, as shown in (Scianca et al., 2020).
At each time tk, the IS-MPC block receives the current footstep plan, that is a sequence of footstep timings
FIGURE 5. A schematic visualizing the time discretization as well as the two horizons used in the scheme.
On the basis of the current footstep plan, IS-MPC computes a 3D CoM trajectory and an adapted set of footsteps over a control horizon5 Tc = C ⋅ δ, with Tc ≤ Tp (see Figure 5). As mentioned in the previous section, this is obtained via two sequential stages, QP-z and QP-xy. Both these programs hinge upon the definition of a suitable ZMP constraint for balance, which will therefore be introduced next.
5.1 Moving ZMP Constraint
As discussed in Section 2, the robot is dynamically balanced if the ZMP is inside the convex hull
The center of the moving constraint region is denoted by
where the αj(t)’s are piecewise-linear functions of time, ranging between 0 and 1 and such that
FIGURE 6. The profile of functions αj(t) used to express the center of the moving constraint region, for a case with F =2. In this example, α0(t) relates the center of the moving constraint to the position of the current support foot, while α1(t) and α2(t) relate it to the position of two predicted footsteps.
Note that, in our formulation, the moving constraint region is planar. Thus, the vertical coordinate zz of the ZMP must at all times coincide with zmc, which is uniquely determined based on the footstep plan. In fact, Eq. 8 entails that zmc only depends on the vertical coordinates of the footsteps, which are fixed and cannot be modified by IS-MPC.
5.2 Vertical QP
The goal of QP-z is to generate the vertical CoM trajectory. The prediction model is given by (3), with the decision variable fz being piecewise-constant over intervals of duration δ:
In order to avoid slipping while walking, the GRF must satisfy a condition of sufficient friction, commonly expressed through a friction cone (Ahn and Cho, 2021). Enforcing this directly as a constraint, however, would introduce a nonlinear coupling between the vertical GRF and the horizontal CoM accelerations. To avoid this, we impose the simplified constraint
where μ is the friction coefficient and
In correspondence of flight phases (the robot is running), constraint (9) is replaced by
in compliance with model (7).
Collect all decision variables in a vector
along with the predicted state variables
and reference values for the absolute CoM height
In the latter,
QP-z is formulated as
where αz, βz are positive weights. The first term of the cost function takes care of tracking the reference CoM height, while the second term penalizes sudden height variations. The last term is included to reduce the difference between consecutive samples of the vertical GRF, so as to generate a smoother GRF profile.
As discussed in Section 4, the vertical GRF produced in output by QP-z is not directly applied to the system. Instead, the associated vertical CoM trajectory will be sent to the kinematic controller and converted to suitable joint commands for the robot.
5.3 Horizontal QP
The goal of QP-xy is to complete the generation of the CoM trajectory with the horizontal components. Thanks to the assumption that all footsteps have the same orientation, QP-xy itself can be split into two decoupled programs QP-x and QP-y. We will present in detail the formulation of QP-x, with the understanding that QP-y is formally identical except when explicitly noted.
When the robot is walking, its CoM obeys the horizontal dynamics (4–5), where the time-varying λ should be computed via (6) using the vertical CoM trajectory generated by QP-z. In particular, we can sample-and-hold λ(t) to obtain a piecewise-time-invariant prediction model:
In time intervals where the robot is running, its CoM obeys the flight-phase dynamics (7), and thus the prediction model becomes
which, as already noticed, is equivalent to setting λk+i = 0 in (12).
QP-x will enforce five kinds of constraints, namely: ZMP constraints, ground patch constraints, footstep kinematic constraints, swing foot constraints, and finally a special stability constraint.
,x and dz,y. Following the discussion in Section 5.1, in a generic time interval we guarantee balance by confining the ZMP to a moving constraint region having the same shape:
where the coordinates (xmc, ymc) of the moving constraint region are given by (8).
5.3.2 Ground Patch Constraint
To maintain linearity of the constraints, we do not allow footstep adaptations that would move the footstep to a ground patch located at a different height from that of the candidate footsteps. For this reason, we enforce the following constraint:
In this article, we assume the full knowledge of the environment and therefore the above computation can also be performed offline. For algorithmic solutions to this problem see, for example, (Deits and Tedrake, 2014; Deits and Tedrake, 2015).
5.3.3 Kinematic Constraints
The third type of constraint ensures that footsteps are placed so as to be kinematically realizable by the robot. This means enforcing a box constraint on the step length:
where σ ∈ (0, 1),
5.3.4 Swing Foot Constraint
To prevent the leg joint speed from exceeding a feasible range during the swing foot motion, we enforce the following swing foot constraint (Herdt et al., 2010)
5.3.5 Stability Constraint
All models of humanoid robot dynamics have an underlying unstable behavior, due to the balancing nature of these systems. As a consequence, keeping the ZMP inside the support polygon is not sufficient to generate a successful gait. In fact, a perfectly acceptable ZMP trajectory could be associated with a diverging CoM trajectory, making the resulting motion infeasible for the robot in practice. The goal of this section is to derive a constraint capable of guaranteeing the internal stability of the scheme, that is, keeping the CoM trajectory bounded with respect to the ZMP.
In previous works (Scianca et al., 2016; Scianca et al., 2020), we approached the above problem by enforcing a stability condition for the LIP model, in which zc − zz is assumed to be constant. Let us consider the same model as a starting point:
where λLIP = g/(zc − zz) is constant. Using the new coordinate
the unstable dynamics is isolated and expressed as:
The instability of system (17) implies that, for a generic ZMP trajectory, xu would diverge with respect to xz, even if the ZMP is inside the support polygon. However, at time tk, the system can be initialized in such a way that the free evolution of (19) exactly cancels the component of the forced evolution that would diverge with respect to xz, allowing the CoM trajectory to remain bounded with respect to the ZMP. This special initialization depends on the future ZMP trajectory, and is identified by the following proposition.
Proposition 1. Assume that |xz (t′) − xz(t)| ≤ a + b (t′ − t), ∀t′ ≥ t, for some a, b > 0, and that
then, system (17) is internally stable, that is, xc is bounded with respect to xz:
Proof. See Supplementary Appendix S1.The first hypothesis in this proposition requires that the future ZMP trajectory is bounded by a linear function. This allows the ZMP to be discontinuous, but it does not allow arbitrarily large skips. Note that this requirement is in our case always satisfied due to the presence of the kinematic constraint. In fact, a bounding linear function for the future ZMP can easily be derived by considering the maximum allowed step size and the step duration. We omit however the details of this calculation as it is fairly tedious, and the specific expression of this bound offers no additional insight.This result is not directly applicable to the present setting, because system (12) is time-varying and therefore the coordinate change (18) does not actually isolate the unstable dynamics. To recover a favorable situation, we will assume in our prediction model that λk+i = λLIP for t ≥ tk+C, consistently with the fact that the control horizon does not extend beyond tk+C. This is equivalent to defining a modified prediction model alternative to (12) given by
Proposition 2. Assume that |xz (t′) − xz(t)| ≤ a + b (t′ − t), ∀t′ ≥ t, for some a, b > 0, and that the current state
Then, system (21) is internally stable, that is, xc is bounded with respect to xz:
Proof. System (21) is equivalent to system (17) for t ≥ tk+C. Hence, Prop. 1 states that if
The thesis follows immediately.Condition (22) is noncausal, as it depends on the entire future ZMP trajectory. To find a causal condition we proceed, similarly to (Scianca et al., 2020), by conjecturing an anticipative tail, that is, a ZMP trajectory
To convert this causal condition into a constraint, it is necessary to express the integral in the left hand side in terms of the MPC decision variables, that is the samples of the piecewise-constant ZMP. This requires explicitly stating the form of the state transition matrix Φ(tk+C, t) and the input matrix B(t). In particular, it is sufficient to specify these terms at each sampled time instant tk+i along the control horizon, recalling that their expression will vary depending on whether tk+i belongs to a support or a flight phase. The state transition matrix from tk+i to tk+C is given by:
where in support and flight phases we have, respectively,
Similarly, for the input matrix we have, respectively,
in support and flight phases.Substituting these into (24) leads to the stability constraint
5.3.6 QP-x Formulation
Collect the candidate footstep coordinates over the control horizon in vector
and the decision variables in vectors
QP-x is formulated as:
The first term of the cost function depends on all decision variables, directly on the ZMP and indirectly on the footsteps (through
An analogous problem QP-y is formulated for the y component, the only difference being that the kinematic constraint is given by the second expression in (15). Note that the ZMP constraints are obviously not enforced in time intervals belonging to flight phases; for the same reason, no ZMP decision variables are associated with such intervals. The constraints on the footstep placement are instead activated or deactivated depending on the gait type (walking/running) and the specific phase (single/double support, support/flight). Figure 8 recapitulates this alternation.
The first ZMP sample
We conclude this section with a comment regarding the feasibility of QP-x and QP-y. Even if the stability constraint impacts the feasibility of the scheme (a detailed analysis for the planar walking case has been provided in Scianca et al. (2020)), there can be situations in which other constraints may conflict. For instance, reaching the assigned ground patch could conflict with the kinematic constraints if a strong deviation from the footstep plan has been previously required; in these situations, the horizontal QP can become unfeasible. This is a limitation that can be overcome by using a motion replanner rather than just relying on the local footstep adjustment capabilities of the controller.
6 A Velocity-Driven Footstep Planner
This section describes a possible method for determining the required inputs for the proposed scheme in response to a high-level forward velocity command. We assume that the 3D map of the environment, called “ elevation map”, is known. In addition to determining footstep placements and timing, we propose a simple criterion to decide for each step whether it will be realized by walking or running, and the corresponding reference CoM height relative to the ZMP.
The basic idea is that a change in velocity should be reflected in a change of both the spacing between the footsteps and their timing (Scianca et al., 2020). The goal is, given a desired velocity v, to determine a step length L and a step duration T such that L/T = v. The method assumes that a triplet of cruise parameters
A deviation from
By setting ΔL = αΔT with α > 0, one can derive an expression for the step duration realizing v
which has to be divided into phases (single/double or support/flight). Such division can be directly chosen based on typical values for the ratio between phase durations in human walking and running. As a general rule based on human biomechanics (Novacheck, 1998), the double support duration should constitute the 20–40% of the step duration. Although this proportion should be inverted for flight phases during the run, here we will use for simplicity the same proportion in running gaits. The associated step length is easily calculated as L = vT.
This procedure outputs a larger L if a larger velocity v is commanded. It stands to reason that shorter steps should be realized by walking, whereas larger steps require running. In practice one may identify a step length threshold Lmax and command the scheme to generate a running step (mode = run) when L > Lmax, and a walking step (mode = walk) otherwise. Being this simple planner conceived for reference forward velocity tracking, the candidate’s footsteps on the y direction simply alternate with a fixed foot distance with respect to the initial CoM position.
So far, a triplet
Finally, we need to assign the reference CoM height relative to the ZMP during the step, which will be used in (11). In particular, we choose a constant reference value
In order to validate the proposed approach, we now present dynamic simulations on an HRP-4 humanoid robot in the DART dynamic environment. The actual torque limits of the platforms were not considered, as the primary objective here is to demonstrate the dynamic stability of a physically simulated robot performing dynamic gaits generated by our IS-MPC scheme. It should be noted in the fact that at the time of conducting this research, just a small percentage of the currently available humanoid robots would be able to perform dynamic running in a 3D environment, for example, Atlas by Boston Dynamics or ongoing projects such as the MIT humanoid (Chignoli et al., 2021) and Kangaroo by PAL Robotics8.
In all the simulations the kinematic controller runs at 100 Hz, while the vertical and horizontal QPs are solved using hpipm (Frison and Diehl, 2020). The simulations use the following parameters: sampling time δ = 0.01 s, control horizon Tc = 0.7 s, preview horizon Tp = 1.8 s,
Clips of all simulations are included in Supplementary Video S1.
7.1 Walking on Flat Ground With Variable CoM Height
For this simulation, the ground is flat and therefore considered as a single horizontal patch. The candidate footstep positions are equally spaced by 0.15 m along x, and displaced by ± 0.18 m on y. Every step lasts 0.7 s, divided in 0.5 s for the single support and 0.2 s for the double support. The desired vertical CoM/ZMP displacement h varies throughout the simulation. It is initially 0.7 m for the first four steps, then is raised to 0.77 m for four subsequent steps. At each of the following five steps, it is progressively decreased by 0.02 m until reaching the final value of 0.58 m.
The simulation results are reported in Figures 9, 10. The stroboscopic motion of Figure 9 also shows the actual CoM trajectory; note the tracking of the variable reference height. In the first plot of Figure 10, we display the CoM and ZMP trajectory in the (x, y) plane for the first eleven seconds of the gait. The CoM trajectory is clearly bounded with respect to the ZMP, proving the effectiveness of the proposed stability constraint. The second plot reports the vertical GRF computed by QP-z, which generates the variations of the CoM height. Note its smooth profile which is due to the last term of the cost function of QP-z. See the Supplementary Video S1.
FIGURE 10. Walking on flat ground with variable CoM height: CoM/ZMP trajectories (left) and vertical GRF (right).
7.2 Running on Flat Ground
For the running case, the spacing of the candidate footsteps along x is increased to 0.3 m, and the step duration is now 0.3 s for the single support phase and 0.15 s for the flight phase. At the start, an initial walking step is performed in preparation for running. The desired vertical CoM/ZMP displacement is kept constant at h = 0.7 m. The arm swing is enforced by a dedicated arm controller; its frequency is the same as the gait, while the amplitude is empirically assigned, similar to what was done in Kajita et al. (2002).
The results are shown in Figures 11, 12. The vertical CoM motion clearly exhibits the typical behavior observed in human biomechanics (Blickhan, 1989; Novacheck, 1998), consisting of the alternation between absorption and a generation phase: absorption is characterized by a horizontal deceleration and vertical descent of the CoM, while during the generation phase the CoM accelerates horizontally and rises in height until reaching an apex. Figure 12 shows that the ZMP trajectory is discontinuous (consistent with the fact that the ZMP is not defined during the flight phases) and the CoM trajectory is bounded. The vertical GRF has a smooth profile while presenting small jumps when starting or concluding a support phase, due to constraints (9). However, the last term in the cost function of QP-z pushes the vertical GRF at the transition towards the minimum admissible value, which in itself is quite low with respect to the GRF peaks during running.
Overall, we were able to reach a maximum forward velocity of 1.25 m/s. In our opinion, this value is due to: i) the simplified model and ii) the kinematic controller. The effect of the simplified model could be alleviated by using a stabilizer based on the ZMP feedback. As for the kinematic controller, this can be improved by properly addressing the foot impact phase.
The results are illustrated in more detail in Supplementary Video S1, in which we also included a push recovery simulation to further show the versatility of the scheme. In this simulation, HRP-4 is pushed twice by forces lasting 0.1 s of (50, 75, 0) N and (0, 75, 0) N, respectively. By adapting the footstep positions the robot can withstand the disturbances and proceed with the running gait. The same pushes would destabilize the robot if step adaptation was removed.
7.3 Walking and Running in a World of Stairs
In this simulation, the robot moves in a world of stairs. An additional overhanging obstacle was placed along the path, which requires the robot to lower its CoM in order to pass below it. The footstep plan is generated by the planner described in Section 6, where we have set α = 0.4,
Simulation results are reported in Figures 13, 14. The robot starts walking (mode = walk) with a step timing of 0.94 s (divided in 0.74 s for the single support phase and 0.2 s for the double support phase) and a step length along the sagittal direction of 0.28 m. When the reference velocity changes, the step length increases by 0.48 m, triggering a transition to running (mode = run); the corresponding durations are 0.24 s for the single support phase and 0.2 s for the flight phase. The sagittal velocity plot shows a satisfactory tracking of the reference sagittal velocity, while the CoM height follows the ground patch profile during both stair climbing and running.
FIGURE 14. Walking and running in a world of stairs: reference and actual sagittal velocity (left), reference and actual CoM height (right).
See the Supplementary Video S1.
In order to validate the proposed scheme on a physical platform, we tested the algorithm on the small humanoid OP3 by ROBOTIS. The same parameters adopted for the HRP-4 simulations were used, except for the sampling time δ = 0.025 s (used both by the MPC and the kinematic controller),
The environment is composed of two ground patches, that is, the floor on which the robot is initially standing and slightly elevated ground. Along its path, the robot will encounter an overhanging obstacle, which has to be avoided by lowering the CoM. The footstep plan is generated with the following parameters: α = 0.4,
The algorithm runs on the robot onboard computer (INTEL NUC with Intel Core i3-7100U 2 × 2.40 GHz processor and 8 GB 2133 MHz RAM) in real time. Each QP has a total of 28 decision variables (Tc/δ = 0.7/0.025), while the number of constraints varies depending on which are active at a particular time (see Figure 8): QP-z always has exactly 28 constraints, while both QP-x and QP-y oscillate around 60.
Snapshots of the experiment are reported in Figure 15. This experiment validates the proposed scheme on a small-sized humanoid, demonstrating its practical usability on a different platform with respect to the HRP-4 adopted for the dynamic simulations; see also the Supplementary Video S1.
FIGURE 15. Snapshots from the experiment. OP3 walks forward and steps onto an elevated ground patch then lower its CoM to crouch in order to pass below an obstacle. The cable attached to the robot’s back is used for power supply only and does not sustain the robot while walking.
9 Running Over Tilted Surfaces
The proposed method has been presented in the context of a world of stairs, where all contact surfaces are parallel. In this section, we show that it can be easily extended to the case of running over tilted surfaces. Although it would be possible to devise an extension for performing both walking and running motions, considering the case of running simplifies the problem. The reason is that, in a running gait, the robot is in contact with the tilted ground with a single foot at a time, and there is no double support phase. This is convenient because introducing double support over tilted surfaces would require adapting the definition of moving constraint to the case of nonparallel contact surfaces.
In Section 5.2 we considered friction by means of a constraint on the vertical GRF (9). Avoiding slipping on tilted surfaces requires higher friction because the tangential component of the GRF must also compensate for the tangential component of gravity. In order to account for this, it is sufficient to increase the minimum required vertical GRF by multiplying the right side of (9) by an appropriate coefficient which depends on the maximum slope of the ground. Since the environment is fully known, this can also be done on a per-footstep basis, by increasing the minimum GRF only when the foot is on a tilted surface, resulting in a less conservative GRF profile.
In order to maintain a formulation that is as close as possible to that presented in Section 5, we express the position of the CoM and ZMP in a frame that has the same orientation of the tilted patch (the z axis is the normal to the patch plane). In this frame, the vertical CoM dynamics assumes the expression:
and the horizontal dynamics become:
Below, we modify the gait generation scheme described in Sections 5.2 and 5.3 by writing the prediction model based on (25–27) and performing an indirect compensation of the disturbance in the stability constraint (Smaldone et al., 2019).
In particular, the prediction model becomes:
while the corresponding stability condition is given by the following result.
Proposition 3. Assume that |xz (t′) − xz(t)| ≤ a + b (t′ − t), ∀t′ ≥ t, for some a, b > 0, and that the current state
where the disturbance matrix is D = (0 1)T. Then, system (28) is internally stable, that is, xc is bounded with respect to xz:
Proof. Following the same steps of the proof of Prop. 2 for system (28) and using (Smaldone et al., 2019)
we get (29).Condition (29) can be written as a stability constraint as follows:
and used in the horizontal QP whenever required in order to step on tilted surfaces. Recall that the disturbance terms in (30) depend only on the tilting angle of the inclined planes and this information is assumed to be available from the footstep plan9.
9.1 Running Over Tilted Surfaces: Simulation
In this simulation the robot must run in an environment with two ramps (inclined by ± 8°) placed across the x path. A footstep plan is generated in response to the input velocity v = 0.35 m/s, which is switched to 1.1 m/s before the first ramp and then set back to the initial value after the second. We adopted the following parameters for the footstep planner: α = 0.4,
Results are shown in Figures 16, 17. The CoM/ZMP trajectories, representing five steps of the running gait, highlight the typical effect of the proposed form of indirect disturbance compensation, consisting of a CoM trajectory that leans against the effect of the perturbation introduced by the tilted surface.
FIGURE 17. Running over tilted surfaces. CoM/ZMP trajectories (left) and CoM height (right). The plots only refer to the portion of the simulation before and after the two tilted surfaces (represented as rectangles in the CoM/ZMP plots). The black arrow represents the equivalent disturbance caused by the tangential component of the gravity when stepping on the tilted patch.
We presented a gait generation algorithm that extends our IS-MPC to both the case of 3D walking and running gaits. The proposed scheme is articulated as a two-stage MPC, and features a new stability constraint for the VH-IP model, ensuring bounded CoM trajectories with respect to the ZMP. A simple planner which provides the required input for the scheme is provided, as well as an extension to the case of running over tilted surfaces. Dynamic simulation on HRP-4 and an experiment on the OP3 humanoid are presented in order to validate the framework.
The strengths of the presented method can be summarized as:
• the effect of the CoM height is fully accounted for in the balance condition, without the need to bound height variations or approximating the horizontal dynamics;
• the resulting scheme is completely formulated as a set of quadratic optimizations subject to linear constraints, which allows for a very efficient MPC implementation, as demonstrated by it running onboard on the OP3 humanoid robot;
• the stability constraint included in the formulation provides a guarantee of bounded CoM/ZMP trajectories;
• the generated trajectories are very natural and compatible with what is observed in human walking/running, as reported by studies in biomechanics.
Future work will be aimed at extending the framework in order to include features to achieve robustness, as well as the integration of the scheme with a more sophisticated footstep planner. Furthermore, the time-varying formulation of the stability constraint can be abstracted from the specific setting of variable height gait generation, and possibly applied to a different scenario modeled as a time-varying unstable system.
Data Availability Statement
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.
All authors contributed to the conceptual development of this work. FS and NS implemented the proposed method and wrote the initial draft. All authors revised and approved the submitted version.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/frobt.2022.876613/full#supplementary-material
1Other studies use the Centroidal Moment Pivot (CMP) point (Popovic et al., 2005) instead of the ZMP for gait generation and/or push recovery, for example, Griffin et al. (2017) and Shafiee-Ashtiani et al. (2017). The CMP purely encodes linear forces and thus does not require to assume the derivative of the centroidal angular momentum to be zero.
2Throughout the article, points are represented by their coordinate vectors.
3An alternative (Zamparelli et al., 2018) is to choose zz as a control input in such a way that λ is constant. The approach adopted in this article is less constraining in terms of the CoM trajectories.
4For compactness, we do not discuss in detail the module which provides the swing foot trajectory. It uses polynomial functions whose endpoints are adjusted in order to match the footstep positions, for example, see (Scianca et al., 2017).
5Having two different horizons for planning and control is a common occurrence, as the window over which the planner operates is generally independent of the window over which the MPC optimizes. In fact, the first can range from a few steps up to the total duration of the locomotion task, while the second is usually much shorter in order to make the problem computationally tractable online. Furthermore, as shown in Scianca et al. (2020), having Tc < Tp is necessary to achieve recursive feasibility, as the information in the preview horizon is used to mitigate the effect of the limited control horizon on the stability of the generated trajectories.
6The choice of rectangular regions simplifies the exposition, but the reasoning can be generalized to any convex shape.
8A practical method to validate the result of our scheme on an actual robotic platform is to run a set of dynamic simulations. The generated motions can be evaluated based on the joint torques necessary to realize them, which can be directly obtained by querying the simulator. If the torques turn out to be excessive, one can try to modify some of the parameters in order to find a set that allows to produce a realizable motion (our results indicate that the duration of flight phases and step length are the most critical).
9This information is in practice also used to generate a swing foot trajectory that is consistent with the slope of the terrain.
Aboudonia, A., Scianca, N., De Simone, D., Lanari, L., and Oriolo, G. (2017). “Humanoid Gait Generation for Walk-To Locomotion Using Single-Stage MPC,” in 17th IEEE-RAS Int. Conf. on Humanoid Robots, 178–183. doi:10.1109/humanoids.2017.8239554
Brasseur, C., Sherikov, A., Collette, C., Dimitrov, D., and Wieber, P.-B. (2015). “A Robust Linear MPC Approach to Online Generation of 3D Biped Walking Motion,” in 2015 IEEE-RAS Int. Conf. on Humanoid Robots, 595–601. doi:10.1109/humanoids.2015.7363423
Caron, S., Kheddar, A., and Tempier, O. (2019b). “Stair Climbing Stabilization of the HRP-4 Humanoid Robot Using Whole-Body Admittance Control,” in 2019 IEEE Int. Conf. on Robotics and Automation, 277–283. doi:10.1109/icra.2019.8794348
Chignoli, M., Kim, D., Stanger-Jones, E., and Kim, S. (2021). “The MIT Humanoid Robot: Design, Motion Planning, and Control for Acrobatic Behaviors,” in 2020 IEEE-RAS Int. Conf. on Humanoid Robots, 1–8. doi:10.1109/humanoids47582.2021.9555782
Deits, R., and Tedrake, R. (2015). “Computing Large Convex Regions of Obstacle-free Space through Semidefinite Programming,” in Algorithmic Foundations of Robotics XI (Springer), 109–124. doi:10.1007/978-3-319-16595-0_7
Englsberger, J., Kozlowski, P., Ott, C., and Albu-Schäffer, A. (2016). Biologically Inspired Deadbeat Control for Running: from Human Analysis to Humanoid Control and Back. IEEE Trans. Robot. 32, 854–867. doi:10.1109/tro.2016.2581199
Englsberger, J., Ott, C., and Albu-Schäffer, A. (2013). “Three-dimensional Bipedal Walking Control Using Divergent Component of Motion,” in 2013 IEEE/RSJ Int. Conf. on Intelligent Robots and Systems, 2600–2607. doi:10.1109/iros.2013.6696723
Griffin, R. J., Wiedebach, G., Bertrand, S., Leonessa, A., and Pratt, J. (2017). “Walking Stabilization Using Step Timing and Location Adjustment on the Humanoid Robot, Atlas,” in 2017 IEEE/RSJ Int. Conf. on Intelligent Robots and Systems, 667–673. doi:10.1109/iros.2017.8202223
Herdt, A., Diedam, H., Wieber, P.-B., Dimitrov, D., Mombaur, K., and Diehl, M. (2010). Online Walking Motion Generation with Automatic Footstep Placement. Adv. Robot. 24, 719–737. doi:10.1163/016918610x493552
Herdt, A., Perrin, N., and Wieber, P.-B. (2012). “LMPC Based Online Generation of More Efficient Walking Motions,” in 2012 IEEE-RAS Int. Conf. on Humanoid Robots, 390–395. doi:10.1109/humanoids.2012.6651549
Hopkins, M. A., Hong, D. W., and Leonessa, A. (2014). “Humanoid Locomotion on Uneven Terrain Using the Time-Varying Divergent Component of Motion,” in 2014 IEEE-RAS Int. Conf. on Humanoid Robots, 266–272. doi:10.1109/humanoids.2014.7041371
Hu, K., Ott, C., and Lee, D. (2014). “Online Human Walking Imitation in Task and Joint Space Based on Quadratic Programming,” in 2014 IEEE Int. Conf. on Robotics and Automation, 3458–3464. doi:10.1109/icra.2014.6907357
Kajita, S., Kanehiro, F., Kaneko, K., Fujiwara, K., Harada, K., Yokoi, K., et al. (2003). “Biped Walking Pattern Generation by Using Preview Control of Zero-Moment Point,” in 2003 IEEE Int. Conf. on Robotics and Automation, 1620–1626.
Kamioka, T., Watabe, T., Kanazawa, M., Kaneko, H., and Yoshiike, T. (2015). “Dynamic Gait Transition Between Bipedal and Quadrupedal Locomotion,” in 2015 IEEE/RSJ Int. Conf. on Intelligent Robots and Systems, 2195–2201. doi:10.1109/iros.2015.7353671
Koolen, T., Posa, M., and Tedrake, R. (2016). “Balance Control Using Center of Mass Height Variation: Limitations Imposed by Unilateral Contact,” in 2016 IEEE-RAS Int. Conf. on Humanoid Robots, 8–15. doi:10.1109/humanoids.2016.7803247
Li, S., Chen, H., Zhang, W., and Wensing, P. M. (2021). “Quadruped Robot Hopping on Two Legs,” in 2021 IEEE/RSJ Int. Conf. on Intelligent Robots and Systems, 7448–7455. doi:10.1109/iros51168.2021.9636120
Popovic, M. B., Goswami, A., and Herr, H. (2005). Ground Reference Points in Legged Locomotion: Definitions, Biological Trajectories and Control Implications. Int. J. robotics Res. 24, 1013–1032. doi:10.1177/0278364905058363
Scianca, N., Cognetti, M., De Simone, D., Lanari, L., and Oriolo, G. (2016). “Intrinsically Stable MPC for Humanoid Gait Generation,” in 16th IEEE-RAS Int. Conf. on Humanoid Robots, 101–108. doi:10.1109/humanoids.2016.7803336
Scianca, N., Modugno, V., Lanari, L., and Oriolo, G. (2017). “Gait Generation via Intrinsically Stable MPC for a Multi-Mass Humanoid Model,” in 2017 IEEE-RAS Int. Conf. on Humanoid Robots, 547–552. doi:10.1109/humanoids.2017.8246926
Secer, G., and Cinar, A. L. (2020). “A Momentum-Based Foot Placement Strategy for Stable Postural Control of Robotic Spring-Mass Running with Point Feet,” in 2020 IEEE/RSJ Int. Conf. on Intelligent Robots and Systems, 3620–3626. doi:10.1109/iros45743.2020.9340980
Shafiee-Ashtiani, M., Yousefi-Koma, A., and Shariat-Panahi, M. (2017). “Robust Bipedal Locomotion Control Based on Model Predictive Control and Divergent Component of Motion,” in 2017 IEEE Int. Conf. on Robotics and Automation, 3505–3510. doi:10.1109/icra.2017.7989401
Smaldone, F. M., Scianca, N., Modugno, V., Lanari, L., and Oriolo, G. (2019). “Gait Generation Using Intrinsically Stable MPC in the Presence of Persistent Disturbances,” in 19th IEEE-RAS Int. Conf. on Humanoid Robots, 651–656. doi:10.1109/humanoids43949.2019.9035068
Sugihara, T., Imanishi, K., Yamamoto, T., and Caron, S. (2021). “3D Biped Locomotion Control Including Seamless Transition Between Walking and Running via 3D ZMP Manipulation,” in 2021 IEEE Int. Conf. on Robotics and Automation, 6258–6263. doi:10.1109/icra48506.2021.9561503
Takenaka, T., Matsumoto, T., and Yoshiike, T. (2009a). “Real-time Motion Generation and Control for Biped Robot - 1st Report: Walking Gait Pattern Generation,” in 2009 Int. Conf. on Intelligent Robots and Systems, 1084–1091.
Takenaka, T., Matsumoto, T., Yoshiike, T., and Shirokura, S. (2009b). “Real-time Motion Generation and Control for Biped Robot - 2nd Report: Running Gait Pattern Generation,” in 2009 IEEE/RSJ Int. Conf. on Intelligent Robots and Systems, 1092–1099.
Wensing, P. M., and Orin, D. E. (2013). “High-speed Humanoid Running through Control with a 3D-SLIP Model,” in 2013 IEEE/RSJ Int. Conf. on Intelligent Robots and Systems, 5134–5140. doi:10.1109/iros.2013.6697099
Wieber, P. B. (2006). “Trajectory Free Linear Model Predictive Control for Stable Walking in the Presence of Strong Perturbations,” in 2006 6th IEEE-RAS Int. Conf. on Humanoid Robots, 137–142. doi:10.1109/ichr.2006.321375
Winkler, A. W., Bellicoso, C. D., Hutter, M., and Buchli, J. (2018). Gait and Trajectory Optimization for Legged Systems through Phase-Based End-Effector Parameterization. IEEE Robot. Autom. Lett. 3, 1560–1567. doi:10.1109/lra.2018.2798285
Xiong, X., and Ames, A. D. (2018). “Bipedal Hopping: Reduced-Order Model Embedding via Optimization-Based Control,” in 2018 IEEE/RSJ Int. Conf. on Intelligent Robots and Systems, 3821–3828. doi:10.1109/iros.2018.8593547
Keywords: humanoid robot, model predictive control, legged locomotion, running, variable height inverted pendulum model
Citation: Smaldone FM, Scianca N, Lanari L and Oriolo G (2022) From Walking to Running: 3D Humanoid Gait Generation via MPC. Front. Robot. AI 9:876613. doi: 10.3389/frobt.2022.876613
Received: 15 February 2022; Accepted: 11 May 2022;
Published: 16 August 2022.
Edited by:Enrico Mingo Hoffman, Pal Robotics S.L., Spain
Reviewed by:Yue Hu, University of Waterloo, Canada
Johannes Englsberger, German Aerospace Center (DLR), Germany
Francisco Javier Andrade Chavez, University of Waterloo, Canada
Copyright © 2022 Smaldone, Scianca, Lanari and Oriolo. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Nicola Scianca, firstname.lastname@example.org