An MPC Walking Framework With External Contact Forces
Sean Mason Nicholas Rotella Stefan Schaal Ludovic Righetti
cs.RO


A “Vector-like chiral” fourth family to explain muon anomalies
    [
    December 30, 2023
==============================================================




In this work, we present an extension to a linear Model Predictive Control (MPC) scheme that plans external contact forces for the robot when given multiple contact locations and their corresponding friction cone. To this end, we set up a two-step optimization problem. In the first optimization, we compute the Center of Mass (CoM) trajectory, foot step locations, and introduce slack variables to account for violating the imposed constraints on the Zero Moment Point (ZMP). We then use the slack variables to trigger the second optimization, in which we calculate the optimal external force that compensates for the ZMP tracking error. This optimization considers multiple contacts positions within the environment by formulating the problem as a Mixed Integer Quadratic Program (MIQP) that can be solved at a speed between 100-300 Hz. Once contact is created, the MIQP reduces to a single Quadratic Program (QP) that can be solved in real-time (< 1kHz). Simulations show that the presented walking control scheme can withstand disturbances 2-3x larger with the additional force provided by a hand contact.





§ INTRODUCTION

In bipedal locomotion, online pattern generators have seen much success because of their ability to adapt to continually changing environment and internal states. In order to generate new plans quickly, researchers typically utilize a simplified dynamic model for fast planning of center of mass (CoM) trajectories. The plans are then tracked using a whole-body controller typically utilizing either inverse kinematics or inverse dynamics. Historically, the Linear Inverted Pendulum Model (LIPM) has been the most prominent model for biped walking because it is linear and yet effectively captures the dynamics for walking on flat surfaces. The low dimensionality and linearity of the LIPM allow the dynamics of the model to be regulated through linear control methods for compact formulations and efficient computations. 

In <cit.>, Kajita et. al introduced the preview control of the LIPM, which optimizes the CoM trajectory over a time-horizon. For a linear system with no inequality constraints, this is equivalent to solving the well known finite horizon Linear Quadratic Regulator (LQR), which yields a time-varying feed-forward and feedback policy. If the problem includes inequality constraints, the solution no longer provides a feedback policy that stabilize the trajectory. (One approach used to rectify this is)  Model Predictive Control (MPC). In an MPC framework, feedback is generated by continuously resolving the constrained optimization from the measured state and applying the first feed-forward command. In <cit.>, Center of Pressure (CoP) constraints in the form of linear inequalities are introduced and the problem is solved in a receding horizon manner to produce CoM trajectories robust to perturbations. This approach was further extended to additionally optimize over the footstep locations <cit.>, allowing the robot to take recovery steps and track a reference velocity. Researchers have further developed MPC methods using the LIPM to adjust footstep timing <cit.>, plan 3D trajectories <cit.>, and control the Divergent Component of Motion (DCM) <cit.>. A nonlinear extension also allows the robot to deal with walking around obstacles <cit.>.









The problem of designing receding horizon controllers that handle rough terrain and allow the use of hand contacts remain a largely open area. 
More general models such as the centroidal momentum dynamics have gained popularity to generate patterns allowing multi-contact behaviors <cit.> but they are not convex and are computationally more demanding. While significant progress has been made with these algorithms, they have not been applied for receding horizon control of legged robots.

Simpler models still allowing hand contacts have also been explored. In <cit.>, hand contacts with the environment were introduced in an MPC problem by deriving a non-linear dynamical model of the CoM that included an external wrench. To solve the non-linear optimization problem a Newton scheme was presented in which each step of the optimization was bounded such that the intermediate solutions were always feasible. While external hand contacts were used for stability, the contact locations and timing were predetermined. Additionally, the contact considered was a grasp and rather than a push, thus neglecting friction cone constraints. In <cit.>, an external wrench is also included in the MPC scheme to design a walking pattern generator for physical collaboration. To drive cooperation with a second agent, the robot assumes the role of either a leader or follower in the task of collectively carrying an object. Contact forces are included in the dynamics and cost function such that the problem remains linear in the constraints and quadratic in the cost.

In this paper, we extend simplified walking pattern generators to exploit hand contacts when stability cannot be maintained by stepping alone. We treat the hand contact as a complement to footstep planning. To limit computational complexity and allow receding horizon control, we propose a two-part optimization that automatically calculates how to use hand contacts when the CoM and footstep adjustment alone are not sufficient. The contact locations are chosen among multiple contact points, considering both the time to reach the contact and friction cone constraints. Our decomposition maintains the computational simplicity of preview controllers based on the LIPM dynamics in the first optimization. The second part then selects a hand contact location and contact force to stabilize the robot when necessary by solving a low-dimensional mixed integer quadratic program (MIQP). The result is a controller that reactively decides when and how to use the hands of the robot to further stabilize the walking controller. Push recovery simulation experiments demonstrate that our approach can significantly improve the stability of a walking robot under strong perturbations.



§ NEWTON-EULER DYNAMICS


We begin by deriving the LIPM dynamics with the addition of an external contact force f_c = [f_c^x ,f_c^y,f_c^z], corresponding to a hand contact located at point p = [p^x ,p^y,p^z]. Similar to <cit.>, we separate foot contact forces, f_i, from all other forces in the Newton-Euler Equations (<ref>), (<ref>).


    m(c̈ +g) = f_c + ∑ f_i


    c × m(c̈ +g) + L̇ = p × f_c + ∑ s_i × f_i
,
where c is the CoM, L is the angular momentum taken about the CoM, s_i is the location of the co-planar foot contacts, and the coordinate frame has the z axis aligned with the gravity vector. Dividing the Euler equation by the foot contact forces in the z direction, f_i^z, of the Newton equation and rearranging yields:


    ∑ s_i × f_i/∑ f_i^z  = c × m(c̈ +g) - (p × f_c )+L̇/m(c̈^z +g) -f_c^z.










Under the LIPM assumptions, the angular momentum and the CoM height remain constant (L̇ = 0, c̈^z = 0), with the ground plan at zero (s^z_i = 0). We can then write the equations for the Zero Moment Point (ZMP) which we will refer to as Z_hand (i.e. ZMP with external hand contact).


    Z_hand^x,y = ∑ s^x,y_i f^z_i/∑ f_i^z  =  - m c^z c̈^x,y - mgc^x,y - p^z f^x,y_c + p^x,y f^z_c/mg -f_c^z

Notice that if there is no contact force (i.e. f_c = 0) this term reduces to the familiar ZMP equations for the LIPM, Z_lipm, <cit.> :


    Z^x,y_lipm = ∑ s^x,y_i f^z_i/∑ f_i^z  = c^x,y -c^z c̈^x,y/g.




§ SHIFTING THE ZMP SUPPORT POLYGON

To express the difference between the ZMP with and without a hand contact, we introduce the variable Δ Z such that:


    Δ Z = Z_hand - Z_lipm

or equivalently

    Z_hand = Z_lipm + Δ Z.


By expressing the ZMP this way, one can see how the traditional ZMP from the LIPM shifts with the addition of an external force. As mentioned in <cit.>, we can additionally think about the implication of this force on the ZMP bounds. That is, we can rewrite


    Z≤ Z_lipm + Δ Z ≤Z

as

    Z  - Δ Z ≤ Z_lipm≤Z -  Δ Z.




Eq. (<ref>) implies that it is possible to shift the effective support polygon that bounds Z_lipm through the addition of an external force. Solving for Δ Z yields the nonlinear expression


    Δ Z^x,y = ( p^z/mg -f^z_c) f^x,y_c + ( -p^x,y  -c^z c̈^x,y/g + c^x,y/mg -f^z_c) f^z_c.




§ OPTIMIZATION FORMULATION



Due to constraints on footstep locations (e.g. kinematic limits and environmental obstructions) and CoM motion (e.g. acceleration limits), it is not always feasible to keep the ZMP within the specified bounds. We remedy these situations by shifting the ZMP bounds by specifying an external force with the hand. One approach is to embed Eq. (<ref>) into the optimization resulting in nonlinear dynamics. It also introduces the variables f_ci, i=1, , N, and p for a total of (3N+3) new decision variables, where N is the number of time-steps in the MPC preview horizon. This problem grows even further if one would like to consider points on different surfaces with different constraints (i.e. friction cone and boundary constraints). Because of the increased dimensionality, complexity, and thus the time needed to solve this larger nonlinear optimization in an online manner, we propose a more light-wight approach that remains linear by breaking the problem into two separate optimizations. In the first optimization, we decide both when an external force is necessary and what value Δ Z^x,y is desired to shift the support polygon. In the second optimization, we decide both the location of the hand contact p and the force f_c. While separating the optimization into two stages will result in a sub-optimal solution, we would argue that in the context of a quick recovery motion, feasibility is more important than optimality. The result is a walking controller that hand contacts in a reactive manner when the control over the CoM and footstep locations alone is deemed insufficient. An outline of the algorithm logic is shown in Alg. 1.





 §.§ Stage 1: MPC of CoM Trajectory and Footstep Locations



In the first stage of the optimization, we formulate an MPC scheme to minimize CoM jerk, ZMP tracking error, and CoM velocity tracking error. The scheme is based on the well known approach from <cit.> so we omit the details for brevity. The only difference in our approach is that we write the ZMP constraints with the addition of a trajectory of slack variables, S_k, k = 1  N. We can write this constraint in the same form as Eq. (<ref>), omitting the subscript k, as:


    Z  - S ≤ Z_lipm≤Z - S,


where Z and Z correspond to the upper and lower bounds of a conservative support polygon (see Fig. <ref>). The amount we scale the support polygon boundaries is directly linked to how early the hand contact will be triggered to help. When the ZMP leaves the conservative safety boundaries, ||S|| > 0, we treat this as an indication that additional support is necessary to maintain balance, triggering the use of a hand contact. In order to account for the fact that it will take some time for the hand to reach the contact point location, we constrain the first N_D time steps of the slack variable S to be zero, where N_D is representative of the time-steps it will take for the hand to make contact. Furthermore, we add a weighted quadratic term to minimize S in the cost function to favor solutions that only use footstep and CoM adjust for balance.



 §.§ Stage 2: Hand Contact


Once the first optimization signifies the need for a hand contact, ||S||>0,  a second optimization calculates the necessary contact force and location to shift the ZMP boundaries by Δ Z such that the ZMP is centered within the support polygon (see Fig. <ref>). While the hand has not yet created a contact, we formulate a MIQP that includes all the reachable contact point locations, along with the corresponding friction cones, to determine which, if any, is the optimal contact location to enact the shift specified by the Δ Z trajectory. Once a contact has been made, we remove all other contacts from the optimization and solve for the contact force as single QP that can be solved in real-time. We begin by formulating the problem for when the hand is already in contact.



 §.§ Determining the Contact Force


In Eq. (<ref>), we can see that with the values of Δ Z^x,y determined, we have two equations and three unknowns. Geometrically, the solutions for the 3D force vector form a line can be rewritten in parametric form as:


    f^x,y_c =  p^x,y  +c^z c̈^x,y/g - c^x,y - Δ Z^x,y/p^z f^z_c + mgΔ Z^x,y/p^z.
   
 
 Note that f_c^x,y is a linear function of f_c^z and thus the only decision variable is f_c^z. We include linearized friction cone constraints, by first changing the frame of our variables to be aligned with the surface normal of the contact:
 
 
    [ f^b; f^t; f^n ] = R [ f^x; f^y; f^z ]

 
 and write the following constraints
 
    f^b ≤ u f^n        ,       f^b ≥ -μ f^n 
    
             f^t ≤ u f^n        ,       f^t ≥ -μ f^n
    
             f^n ≤ f^n_max       ,       f^n ≥ 0       ,


where μ is the static coefficient of friction that defines the linearized fiction cone. The cost function to optimize is then:

    f^zminimize       (f^b)^2 + (f^t)^2 + κ(f^n)^2 
       subject to      Friction Cone Constraints (<ref>).


Reducing the gain κ favors solutions that are closer to the center of the friction cone and thus more robust to slipping. This optimization formulation specified by Eq. (<ref>,<ref>,<ref>) will be referred to as QP1. Because the formulation of QP1 shifts the support polygon by exactly Δ Z, the optimization does not consider that the ZMP does not have to be centered in the support polygon. To address this, we include the addition of the new slack variable S_z = [S_z^x, S_z^y] that is bounded as follows:
 
    S_z≤ S_z ≤S_z.


Rewriting Eq. (<ref>) as a function of S_z and f_z results in:




  
    f^x,y_c  =  ( p^x,y  +c^z c̈^x,y/g - c^x,y - Δ Z^x,y/p^z) f^z_c 
     +( 1/p^z) S_z^x,y f^z_c  + ( mg/p^z) S_z^x,y + (mg Δ Z^x,y/p^z).

Unfortunately,  by introducing the new decision variable S_z we have added the bilinear term S_z f^z_c. In order to keep the problem linear and solvable with a QP, we introduce a convex relaxation of the bilinear term by means of a McCormick envelope <cit.>. In doing so, we create the new variables W^x,y = S_z^x,y f^z_c, more compactly written as W, and add the following new constraints to the system:

    S_z  f^z_c        +       S_z  f^z_c       -       S_z f^z_c      ≤       W 
    S_z  f^z_c        +       S_z  f^z_c       -       S_z f^z_c      ≤       W
    S_z  f^z_c        +       S_z  f^z_c       -       S_z f^z_c      ≥       W
    
     S_z  f^z_c       +       S_z  f^z_c        -       S_z f^z_c      ≥       W.



This new formulation will be referred to as QP2. 





 §.§ Determining the Contact Location


In order to determine the contact location, we reason about reachable contact locations and their corresponding friction cone in the form of a MIQP. To define whether a contact has been chosen, we introduce the binary variables b_i ∈0,1,    i=1, , N_c. We then optimize over the time steps k̂ corresponding to the non-zero elements in the Δ Z trajectory.


    f^z_ik̂, b_iminimize       (f^b)_ik̂^2 + (f^t)_ik̂^2 + κ(f^n)_ik̂^2 
       subject to



    ∑_1^N_c b_i = 1



    b_i ⟹{   0,        f^z_ik̂ = 0
       1,         (Friction Cone Constraints)_i
      .
 



§ EXPERIMENTS





 §.§ Optimization Benchmarks

The experiments were conducted in this section using the optimization library Gurobi <cit.>. The average solve while searching among contacts, shown in Table <ref>, was found to vary linearly with the length of the Δ Z trajectory  and the number of contacts points. Once in contact, the optimization problem collapses down to a single QP which can be solved in 0.11 ms.

 

  
  Δ Z Time-steps      4|c||Number of Contact Points (N_c)   

     5     10     20     40   

 16      5.2 ms      9.7 ms     18.7 ms     35.9 ms  

 8       3.2 ms     5.2 ms     9.8 ms     19.8 ms 

 4       1.8 ms     2.9 ms     5.4 ms     10.2 ms 

 
table Solving speed for multiple contact points considering different horizon lengths.







 §.§ QP1 vs QP2

As shown in Fig. <ref>, by solving for contact forces using the formulation of QP2 as opposed to QP1 we expand the feasible solution space from a line to a volume. To compare how the solutions of QP1 and QP2, we performed a parameter sweep over μ and Δ Z and evaluated 5,000 random friction cones orientations and CoM accelerations for the following ranges:  μ =[.1,1] and Δ Z = [-.05,.05], S_z = [-.02, .02]. In total, there were 2,000,000 test cases evaluated, resulting in the following statistics:

 

  
     Success Rate     Average Force (N)

 Force QP1     36.1%     28.0 

 Force QP2     73.6%     11.7





On average, formulation of QP2 favorably found optimal solutions for approximately twice the number test cases at approximately half the force. Because of the convex relaxation of the bilninear term used in QP2, we additionally need to check that the solutions satisfy the original problem constraints. Over the values tested, the maximum force violation was 0.037 N, which is negligible given the accuracy of force tracking capabilities of the robot seen in experiments. 



 §.§ Maximum Push Experiments


We evaluated the proposed control scheme in simulation through a series of push recovery tests on a 31 DoF torque controlled humanoid. We are particularly interested in the relationship between the max impulse the robot can recover from with and without an external force through hand contact with the environment. Additionally, we are interested in how the number of recovery steps before hand contact affects this maximum impulse value. To test this, the robot was pushed towards a wall while walking in place from three separate distances that allow for zero, one, and two recovery steps before hand contact. Additionally, the robot footstep locations were constrained to remain at a minimum distance from the wall (0.35m) so that the upper-body would not make contact unless the arm was purposefully extended (see Fig. <ref>). In order to track the CoM and footstep trajectory, we use an optimization-based whole-body inverse dynamics controller formulated as a QP, similar to <cit.>, to resolve for joint torques at 1kHz. To track the desired forces, we implemented a simple force feedback policy that was added to the resulting joint torques from the whole-body QP solver (τ_qp): 

    τ = τ_qp + J_c^T (f_c + P (f_c - f̂_̂ĉ))


, where J_c is the Jacobian from the CoM to the hand in contact, f̂_̂ĉ is the measured force at the end-effector in global coordinates, and P is a force-feedback gain. In the future, we plan to include these desired forces and friction cone constraints directly in the whole-body QP. The following parameters were used in these experiments: MPC horizon =  1.6 s,  MPC Δ t = 0.1 s, robot mass = 62.5 kg,  CoM height = 0.78m, single support duration = 1.0 s, double support duration = 0.1 s, μ = 1, and the support polygon boundary used in Sec. <ref> scaled the foot dimensions by a factor of 0.1.

 The results of the push experiments are shown in Fig. <ref>. On average, using an external hand contact increased the maximum disturbance impulse by 233%, 156%, and 13% for the experiments where the robot was allowed zero, one, and two recovery steps respectively. While addition of a hand contact greatly increased the walking controllers robustness to external impulses for zero and one recovery steps, in the 2-step recovery situation, the robot had a very high acceleration by the time the wall became close enough to make contact. This resulted in high desired contact forces that were difficult to track. The experiments confirm that because of the short optimization horizon and the reactive nature of this formulation, that approach presented is most effective in close proximity of contacts or when there is limited control authority through stepping. 



Additionally, it is important to note that perfect force tracking was not achieved. Fig. <ref> shows the tracking of the external hand force in one of the experiments. This is important because it is very improbable that on a real humanoid perfect force tracking while making and breaking contacts with the environment will be achieved. The robustness of the approach presented comes from constantly replanning the contact force from the measured robot state. 

Further demonstrations of the behavior of this control scheme is shown in the accompanying video [<http://bit.ly/masonIcra2018>]. The video shows the response to pushes with and without hand contact, choosing between surfaces with different friction properties, and utilizing contacts on surfaces that are rotated at different angles.





§ CONCLUSION

In this work, we presented an extension to an MPC walking control scheme that reactively plans a desired contact location and force without any predetermined timing. In doing so, the formulation considers friction properties of the contact surfaces and the time it takes for the robot to move its hand to the contact location. To this end, we introduced a two part optimization scheme in which the first optimization plans the CoM and footsteps trajectories while the second optimization reasons about the hand contacts. The second optimization is initiated when the first optimization violates some conservative ZMP constraints. This separation into two optimization problems allowed us to reach real-time computation while maintaining low computational complexity. As shown in our experiments, there is a great benefit of using hand contact together with a walking controller when within a one recovery step proximity. Because this control scheme separates walking from hand contacts, there is the possibility to pair the second stage with different LIPM based walking pattern generators and other methods to trigger the need for hand contact. Potential directions for future research would be to use the DCM of the LIPM or incorporate an external force estimator <cit.> to trigger the need for hand contact. 

IEEEtran



