Feedback MPC for Torque-Controlled Legged Robots

TL;DR

This paper introduces a feedback MPC approach using DDP with barrier functions and frequency-dependent costs to enable stable torque control in legged robots at practical update rates.

Contribution

It proposes a novel DDP-based MPC method with barrier functions and frequency-dependent costs for torque-controlled legged robots, addressing update rate limitations.

Findings

Achieves stable locomotion in simulation and hardware

Effectively handles inequality constraints with barrier functions

Reduces sensitivity to high-frequency errors

Abstract

The computational power of mobile robots is currently insufficient to achieve torque level whole-body Model Predictive Control (MPC) at the update rates required for complex dynamic systems such as legged robots. This problem is commonly circumvented by using a fast tracking controller to compensate for model errors between updates. In this work, we show that the feedback policy from a Differential Dynamic Programming (DDP) based MPC algorithm is a viable alternative to bridge the gap between the low MPC update rate and the actuation command rate. We propose to augment the DDP approach with a relaxed barrier function to address inequality constraints arising from the friction cone. A frequency-dependent cost function is used to reduce the sensitivity to high-frequency model errors and actuator bandwidth limits. We demonstrate that our approach can find stable locomotion policies for the torque-controlled quadruped, ANYmal, both in simulation and on hardware.

Equations

$$\displaystyle\mbox{\boldmath$\pi$}(\mathbf{x},t)=\mathbf{u}^{*}(t)+\mathbf{K}(t)\left(\mathbf{x}-\mathbf{x}^{*}(t)\right),$$

$$\underset{u(\cdot )}{\min }\Phi (\mathbf{x}(T))+{\int }_0^{T}L(\mathbf{x}(t),\mathbf{u}(t),t)dt,$$

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

$$\mathbf{x}(0)={\mathbf{x}}_0$$

$${\mathbf{g}}_1(\mathbf{x},\mathbf{u},t)=0$$

$${\mathbf{g}}_2(\mathbf{x},t)=0$$

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

$$\hat{L}(\mathbf{x},\mathbf{u},t)=L(\mathbf{x},\mathbf{u},t)+\mu \sum _{i=1}^{N_{in}}B(h_i(\mathbf{x},\mathbf{u},t)).$$

$$B(h)=\{\begin{array}{ll}-\ln (h),& h\ge \delta ,\\ \beta (h;\delta ),& h<\delta .\end{array}$$

$$\beta (h;\delta )=\frac{1}{2}\left({\left(\frac{h-2\delta }{\delta }\right)}^2-1\right)-\ln (\delta ).$$

$$\delta \dot{\mathbf{x}}\dot{}=\mathbf{A}\delta \mathbf{x}+\mathbf{B}\delta \mathbf{u},$$

$$\mathbf{C}\delta \mathbf{x}+\mathbf{D}\delta \mathbf{u}+\mathbf{e}=0,$$

$$\Phi (\mathbf{x}(T))$$

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

$$\frac{1}{2}\delta {\mathbf{x}}^T{\mathbf{Q}}_L\delta \mathbf{x}+\frac{1}{2}\delta {\mathbf{u}}^T\mathbf{R}\delta \mathbf{u}+\delta {\mathbf{u}}^T\mathbf{P}\delta \mathbf{x},$$

$$\mathbf{K}(t)=\left(\mathbf{I}-{\mathbf{D}}^{†}\mathbf{D}\right){\mathbf{R}}^{-1}\left({\mathbf{P}}^T+{\mathbf{B}}^T\mathbf{S}\right)+{\mathbf{D}}^{†}\mathbf{C},$$

$${\hat{\nu }}_i(\omega )$$

$$r_i(\omega )$$

$$\left[\begin{array}{c}\mathbf{f}(\mathbf{x},{\mathbf{C}}_s{\mathbf{x}}_s+{\mathbf{D}}_s\mathbf{\nu })\\ {\mathbf{A}}_s{\mathbf{x}}_s+{\mathbf{B}}_s\mathbf{\nu }\end{array}\right],\begin{array}{c}{\mathbf{g}}_1(\mathbf{x},{\mathbf{C}}_s{\mathbf{x}}_s+{\mathbf{D}}_s\mathbf{\nu },t)=0\\ \mathbf{h}(\mathbf{x},{\mathbf{C}}_s{\mathbf{x}}_s+{\mathbf{D}}_s\mathbf{\nu },t)\ge 0\end{array}.$$

$$\mathbf{\nu }(\mathbf{x},t)={\mathbf{\nu }}^{\ast }(t)+\left[\begin{array}{cc}{\mathbf{K}}_{\mathbf{\nu },\mathbf{x}}(t)& {\mathbf{K}}_{\mathbf{\nu },{\mathbf{x}}_s}(t)\end{array}\right]\left[\begin{array}{c}\mathbf{x}-{\mathbf{x}}^{\ast }(t)\\ {\mathbf{x}}_s-{\mathbf{x}}_s^{\ast }(t)\end{array}\right].$$

$$\left[\begin{array}{c}\mathbf{u}(\tilde{\mathbf{x}},t)\\ \mathbf{\nu }(\tilde{\mathbf{x}},t)\end{array}\right]=\left[\begin{array}{c}{\mathbf{C}}_s{\mathbf{x}}_s^{\ast }(t)+{\mathbf{D}}_s{\mathbf{\nu }}^{\ast }(t)\\ {\mathbf{\nu }}^{\ast }(t)\end{array}\right]+$$

$$\left[\begin{array}{cc}{\mathbf{D}}_s{\mathbf{K}}_{\mathbf{\nu },\mathbf{x}}(t)& {\mathbf{D}}_s{\mathbf{K}}_{\mathbf{\nu },{\mathbf{x}}_s}(t)+{\mathbf{C}}_s\\ {\mathbf{K}}_{\mathbf{\nu },\mathbf{x}}(t)& {\mathbf{K}}_{\mathbf{\nu },{\mathbf{x}}_s}(t)\end{array}\right]\left[\begin{array}{c}\mathbf{x}-{\mathbf{x}}^{\ast }(t)\\ {\mathbf{x}}_s-{\mathbf{x}}_s^{\ast }(t)\end{array}\right].$$

$$\displaystyle\left\{\begin{array}[]{ll}\dot{\mbox{\boldmath$\theta$}}=\mathbf{T}(\mbox{\boldmath$\theta$})\mbox{\boldmath$\omega$}\\ \dot{\mathbf{p}}=_{W}\!\mathbf{R}_{B}(\mbox{\boldmath$\theta$})\,\mathbf{v}\\ \dot{\mbox{\boldmath$\omega$}}=\mathbf{I}^{-1}\left(-\mbox{\boldmath$\omega$}\times\mathbf{I}\mbox{\boldmath$\omega$}+\sum_{i=1}^{4}{{\mathbf{r}_{EE}}_{j}(\mathbf{q})\times{\mbox{\boldmath$\lambda$}_{EE}}_{j}}\right)\\ \dot{\mathbf{v}}=\mathbf{g}(\mbox{\boldmath$\theta$})+\frac{1}{m}\sum_{i=1}^{4}{\mbox{\boldmath$\lambda$}_{EE}}_{j}\\ \dot{\mathbf{q}}=\mathbf{u}_{J},\end{array}\right.$$

$$\displaystyle\left\{\begin{array}[]{ll}{\mathbf{v}_{EE}}_{j}=\mathbf{0},&\text{if $i$ is a stance leg},\\ {\mathbf{v}_{EE}}_{j}\cdot\hat{\mathbf{n}}=c(t),\quad{\mbox{\boldmath$\lambda$}_{EE}}_{j}=\mathbf{0},&\text{if $i$ is a swing leg},\end{array}\right.$$

$${\mbox{\boldmath$\lambda$}_{EE}}_{j}\in\mathcal{C}(\hat{\mathbf{n}},\mu_{c}),$$

$$h_{cone}={\mu }_c{F}_z-\sqrt{F_x^2+F_y^2}\ge 0.$$

$$h_{cone,ϵ}={\mu }_c{F}_z-\sqrt{F_x^2+F_y^2+{ϵ}^2}\ge 0,$$

$$h_{cone,ϵ}\ge 0⟹h_{cone}\ge 0.$$

Full text

I Introduction

Model Predictive Control (MPC) has gained broad interest in the robotics community as a tool for motion control of complex and dynamic systems. The ability to deal with nonlinearities and constraints has popularized the technique for many robotic applications, such as quadrotor control [1], autonomous racing [2], and legged locomotion [3, 4, 5].

MPC strategies typically optimize an open-loop control sequence for a given cost function over a fixed time horizon. The control sequence is then executed until a new control update is calculated based on the current state estimate. While this strategy assumes that the model is exact and that there are no external disturbances, the repeated optimization provides a feedback mechanism that can correct for modeling errors provided that the control loop can be executed at a sufficiently high rate. However, for high dimensional systems such as legged robots and due to the computational restrictions of mobile platforms, the achievable update rate of the MPC loop is insufficient to effectively deal with model uncertainty and external disturbances.

As a remedy, a separately designed, light-weight motion tracker is often used in practice [6]. The motion tracker runs at a higher rate than the MPC loop and provides feedback correction to the control sequence that was designed by the MPC. For complex systems such as legged robots it is a challenging task to design a controller that tracks arbitrary motions while satisfying the many constraints arising from the locomotion task. The fundamental problem is that such motion trackers do not look ahead in the horizon and therefore cannot anticipate changes in contact configuration. As an alternative, projected (time-varying) Linear Quadratic Regulators (LQR) have been proposed as a framework to automatically design feedback controllers around a given reference trajectory [7, 8]. However, the stabilizing feedback policy is always designed in a secondary stage and often with a different objective function than the one used for computing the optimal trajectories, which leads to inconsistency between the feedback policy and the MPC trajectories.

In this work, we propose a feedback MPC approach for motion control of a legged system and show that the optimized feedback policy can directly be deployed on hardware. We achieve stable locomotion under a very low update rate (15 Hz), and the optimized feedback policy removes the need for a separate motion controller. Furthermore, the modification of the control inputs is consistent with the MPC and thus produces a continuous signal across update instances.

To be able to directly apply the feedback strategy on a legged system, the optimized policy needs to respect all the constraints of the locomotion task such as friction and unilateral constraints of contact forces. To achieve this, we propose to extend the SLQ (Sequential Linear Quadratic) algorithm [9]. We use SLQ in a real-time iteration MPC scheme [10] where the algorithm optimizes a constrained feedback policy, \mbox{\boldmath\pi}(\mathbf{x},t):\mathbb{X}\times\mathbb{R}^{+}\!\to\!\mathbb{U},

[TABLE]

where ${\mathbf{u}^{*}(t)\in\mathbb{U}}$ and ${\mathbf{x}^{*}(t)\in\mathbb{X}}$ are locally optimal input and state trajectories. $\mathbf{K}(t)$ is a time-varying LQR gain matrix which maps the state deviation from $\mathbf{x}^{*}$ to an admissible control correction. We extend the algorithm to problems with inequality constraints using a barrier function method to accurately handle the constraints arising from the friction cone. We further use the frequency-aware MPC approach introduced in [11] to render the resulting feedback policy robust to the bandwidth limitations imposed by real actuators.

We perform experiments in simulation and on a real legged system (Fig. 1) and demonstrate that our approach is able to find robust and stable locomotion polices at MPC update rates as low as 15 Hz, which facilitates onboard execution on mobile platforms with limited computational power.

I-A Related Work

Methods to incorporate robustness explicitly into the MPC methodology have been previously studied in the literature [12]. Min-max MPC [13], for example, optimizes an open-loop control sequence for the worst-case disturbance inside a predefined set. While this formulation appears attractive, it can be overly conservative due to its inability to include the notion of feedback that is inherently present in the receding-horizon implementation of the control [14].

Min-max Feedback MPC was proposed to address this shortcoming by planning over a state-dependent control policy instead of an open-loop feedforward sequence [15]. Unfortunately optimizing the feedback policy for all possible disturbance realizations does not yet scale to the problem dimensions encountered in legged robotics. However, even without considering disturbances, optimizing over the feedback policy has an additional advantage. When the update rate of the MPC loop is low, the feedback policy can provide local correction to the deviation of the real platform from the optimal trajectories. Exploiting this additional aspect of feedback MPC has not yet been fully explored in robotic applications that are subject to path constraints.

The feedback policy that minimizes a cost function for a given dynamical system and path constraints can be computed using the Hamilton-Jacobi-Bellman (HJB) equation [16]. While directly solving this equation for high dimensional systems is prohibitively complex, a variant of the dynamic programming approach known as Differential Dynamic Programming (DDP) [17] has proven to be a powerful tool in many practical applications. The SLQ method that we use in this work is a DDP-based approach which uses a Gauss-Newton approximation. Consequently, it only considers the linearized dynamics instead of a second-order approximation.

Although using the LQR gains derived from a DDP-based approach directly for motion tracking generates promising results in simulation, it dramatically fails on real hardware. This phenomenon has been reported before in other real-world applications of LQR on torque-controlled robots [8, 18]. Focchi et al. [19] have shown that instability can occur if the limitations of the low-level torque controller are neglected in the high-level control design. They have argued that the bandwidth of the low-level controller inversely relates to the achievable impedance of the high-level controller. To this end, to apply the SLQ feedback policy on hardware, we need to encode these bandwidth limitations in our optimization problem. In this work, we use the frequency-aware MPC approach introduced in [11]. This MPC formulation penalizes control actions in the frequency domain and automatically finds a trade-off between the bandwidth limitation of actuators and the stiffness of the high-level feedback policy.

I-B Contributions

We propose a whole-body MPC approach for legged robots, where the actuation commands are computed directly based on the MPC feedback policy. Specifically we present the following contributions which we empirically validate on the ANYmal platform (Fig. 1) in simulation and on real hardware:

• •

We propose to apply feedback MPC for whole-body control of a legged system. To the best of our knowledge, this is the first time that such a control scheme is applied on hardware for motion control of legged robots.

• •

The SLQ algorithm is extended to include inequality constraints through a barrier function method, which allows us to formulate friction cone constraints.

• •

We show that our feedback MPC algorithm directly designs constraint-satisfactory LQR gains without additional computational cost.

• •

A frequency domain design approach is used to incorporate actuation bandwidth limits in the MPC formulation to avoid rendering stiff gains. Thus, the feedback gains can be directly applied to the robot.

• •

We show that the feedback MPC algorithm is capable of bridging the gap between low update-rate MPC and high rate execution of torque commands using only an onboard computer with moderate computational power.

II Method

II-A Problem Definition

Consider the following nonlinear optimal control problem with cost functional

[TABLE]

where $\mathbf{x}(t)$ is the state and $\mathbf{u}(t)$ is the input at time $t$. $L(\cdot)$ is a time-varying running cost, and $\Phi(\cdot)$ is the cost at the terminal state $\mathbf{x}(T)$. Our goal is to find an input trajectory $\mathbf{u}(\cdot)$ that minimizes this cost subject to the following system dynamics, initial condition, and general equality and inequality constraints:

[TABLE]

The feedback policy which minimizes this problem can be calculated using a DDP-based method. A variant of this method for continuous-time systems known as SLQ is introduced in [9], where it solves the above optimization problem in the absence of the inequality constraints in equation (7). This method computes a time-varying, state-affine control policy based on a quadratic approximation of the optimal value function in an iterative process. The SLQ approach uses a Lagrangian method to enforce the state-input equality constraints in (5). The pure state constraints in (6) are handled by a penalty method.

In SLQ, the simulation (forward pass) and the optimization (backward pass) iterations alternate. Once the backward pass is completed, a forward pass computes a new trajectory based on the improved feedback policy. The local, Linear Quadratic (LQ), approximation of the nonlinear optimal control problem is constructed after each forward pass. The LQ model permits an efficient solution of the approximate problem by solving the Riccati differential equation. The feedback policy is then updated with an appropriate linesearch procedure in the direction of the LQ problem’s solution.

We follow the same SLQ approach and extend the method with inequality constraints through a relaxed barrier function approach.

II-B Relaxed Barrier Functions

Using a barrier function is a well know technique to absorb inequality constraint into the cost function. For each constraint in a given set of $N_{in}$ inequality constraints, a barrier term $B(h)$ is added to the cost

[TABLE]

A widely used barrier function is the logarithmic barrier used in interior-point methods. The optimal solution is approached by letting $\mu\rightarrow 0$ over successive iterations. However, a downside of the log-barrier is that it is only defined over the feasible space, and evaluates to infinity outside. Due to the rollout mechanism in the SLQ approach, one cannot ensure that successive iterations remain inside the feasible region at all time. Furthermore, the Hessian of the log-barrier goes to infinity as one approaches the constraint boundary, which results in an ill-conditioned LQ approximation.

The relaxed barrier functions previously proposed for MPC problems addresses both these issues [20] and is therefore particularly suitable for the SLQ approach. This barrier function is defined as a log-barrier function on the interior of the feasible space, and switched to a different function at a distance $\delta$ from the constraint boundary.

[TABLE]

We use the quadratic extension proposed in [21]:

[TABLE]

The relaxed barrier function, which is continuous and twice differentiable, is plotted as a function of the constraint value in Fig. 2. The quadratic extension puts an upper bound to the curvature of the barrier function, which prevents ill-conditioning of the LQ approximation. Note that by letting $\delta\rightarrow 0$, the standard logarithmic barrier is retrieved. Furthermore, it has been shown that the optimal solution can be obtained for a nonzero value of $\delta$ [22], when the gradient of the penalty term is larger than the Lagrange multiplier of the associated constraint. Optimization with the relaxed barrier function can thus be interpreted as an augmented Lagrangian approach when $h<\delta$ and as a log-barrier method for $h\geq\delta$.

II-C LQ Approximation

With the inequality constraints embedded in the cost function, we obtain the following linearization of the system dynamics in (3) and state-inputs constraints in (5) for a given state trajectory $\mathbf{x}_{k-1}(t)$ and input trajectory $\mathbf{u}_{k-1}(t)$:

[TABLE]

where $\delta\mathbf{x}=\mathbf{x}(t)-\mathbf{x}_{k-1}(t)$, and $\delta\mathbf{u}=\mathbf{u}(t)-\mathbf{u}_{k-1}(t)$ are deviations from the previous iteration, around which the LQ approximation is made. Note that the time-dependency of the matrices was dropped to shorten the notation. The quadratic approximation of the cost in (8) is given by

[TABLE]

which requires access to the second-order approximation of the barrier term and inequality constraints.

With the optimal control problem reduced to an equality constrained LQ approximation, the constrained Riccati backward pass in [9] yields the quadratic value function $V(\mathbf{x},t)=\frac{1}{2}\mathbf{x}^{T}\mathbf{S}(t)\mathbf{x}+\mathbf{s}^{T}(t)\mathbf{x}+s(t)$. This value function induces the optimal feedback policy in equation (1) with feedback gains computed as

[TABLE]

where $\mathbf{I}$ is the identity matrix and $\mathbf{D}^{\dagger}$ is the right pseudo-inverse of the full row rank matrix $\mathbf{D}$. Notice how the feedback gains ensure that the equality constraints are satisfied by projecting the first term to the nullspace of the constraints, and by adding the term $\mathbf{D}^{\dagger}\mathbf{C}$ to satisfy the constraint when the state deviates from the plan.

II-D Frequency Shaping

As discussed in Section I-A, it has been proven difficult to use feedback gains from an LQR design on a torque-controlled robot. We propose to use the frequency-dependent cost function introduced in our previous work [11], which was used to render the feedforward solution robust to high frequency disturbances. In this work, we show that it has a similar effect on the feedback structure. We briefly summarize how the problem is adapted and refer to [11] for further details.

A frequency-dependent cost on the inputs can be introduced by evaluating the cost function on auxiliary inputs $\boldsymbol{\nu}$. The auxiliary inputs $\boldsymbol{\nu}$ are defined by frequency-dependent shaping functions $r_{i}(\omega)$ applied to the system inputs such that

[TABLE]

where $i$ denotes elements associated to individual inputs, $\omega$ is the signal frequency in $\text{\,}\mathrm{rad}\text{\,}{\mathrm{s}}^{-1}$, and $\hat{\nu}_{i}(\omega)$ and $\hat{u}_{i}(\omega)$ are the Fourier transform of the auxiliary input and system input respectively. Following our previous work, we use high pass filters to achieve increased costs at higher input frequencies:

[TABLE]

The transfer function $\mathbf{s}(\omega)=\mathbf{r}^{-1}(\omega)$, with state space realization ${(\mathbf{A}_{s},\mathbf{B}_{s},\mathbf{C}_{s},\mathbf{D}_{s})}$, is constructed such that ${\hat{u}_{i}(\omega)=s_{i}(\omega)\hat{\nu}_{i}(\omega)}$. The original system is augmented with an additional filter state, $\mathbf{x}_{s}$, such that $\tilde{\mathbf{x}}=[\mathbf{x}^{T},\mathbf{x}_{s}^{T}]^{T}$, and optimization is performed w.r.t the auxiliary inputs $\boldsymbol{\nu}$. The augmented system dynamics and state-input constraints are defined as

[TABLE]

The feedback policy obtained from this augmented system is of the form

[TABLE]

After optimization, the original input is retrieved by substituting this policy into the output function of the filter, $\mathbf{u}=\mathbf{C}_{s}\mathbf{x}_{s}+\mathbf{D}_{s}\boldsymbol{\nu}$, resulting in the complete feedback policy

[TABLE]

III Implementation

We apply our approach to the kinodynamic model of a quadruped robot, which describes the dynamics of a single free-floating body along with the kinematics for each leg. The Equations of Motion (EoM) are given by

[TABLE]

where ${}_{W}\!\mathbf{R}_{B}$ and $\mathbf{T}$ are the rotation matrix of the base with respect the global frame and the transformation matrix from angular velocities in the base frame to the Euler angles derivatives in the global frame. $\mathbf{g}$ is the gravitational acceleration in body frame, $\mathbf{I}$ and $m$ are the moment of inertia about the CoM and the total mass respectively. The inertia is assumed to be constant and taken at the default configuration of the robot. ${\mathbf{r}_{EE}}_{j}$ is the position of the foot $j$ with respect to CoM. $\theta$ is the orientation of the base in Euler angles, $\mathbf{p}$ is the position of the CoM in world frame, $\omega$ is the angular rate, and $\mathbf{v}$ is the linear velocity of the CoM. $\mathbf{q}$ is the vector of twelve joint positions. The inputs of the model are the joint velocity commands $\mathbf{u}_{J}$ and end-effector contact forces {\mbox{\boldmath\lambda}_{EE}}_{j} in the body frame.

III-A Equality Constraints

The equality constraints depend on the mode of each leg at a certain point in time. We assume that the mode sequence is a predefined function of time. The resulting mode-depended constraints are

[TABLE]

where ${\mathbf{v}_{EE}}_{j}$ is the end-effector velocity in world frame. These constraints ensure that a stance leg remains on the ground and a swing leg follows the predefined curve $c(t)$ in the direction of the local surface normal $\hat{\mathbf{n}}$ to avoid foot scuffing. Furthermore, the constraints enforce zero contact force at swing legs.

III-B Inequality Constraints

Our proposed relaxed barrier method allows to model the friction cone without the commonly used polytope approximation. The cone constraint for each end-effector,

[TABLE]

is defined by the surface normal and friction coefficient ${\mu_{c}=0.7}$. After projecting the contact forces to the local frame of the surface, a canonical second-order cone constraint is found in terms of local contact forces ${\mathbf{F}=[F_{x},F_{y},F_{z}]}$. An effective cone constraint used in conjunction with barrier methods [23] is given by

[TABLE]

However, the gradient of this constraint is not defined at $\mathbf{F}=0$, which causes numerical issues close to the origin. While for interior point methods this problem can be solved by using the squared constraint, $\mu_{c}^{2}F_{z}^{2}-F_{x}^{2}+F_{y}^{2}\geq 0$, this strategy does not work well together with the relaxed barrier function due to the saddle point it introduces at the origin. Since the relaxed barrier function allows infeasible iterates, the solutions can cross the origin and end up in the negative reflection of the cone, which became a feasible region through the squaring operation. We, therefore, use the perturbed cone

[TABLE]

which is differentiable at the origin, remains infeasible for any negative $F_{z}$, and is a conservative lower bound for the original cone (21). It therefore holds that

[TABLE]

In Fig. 3 the level sets of this constraint are compared to the original cone. It can be seen that the constraint is convex and the zero crossing of $h_{cone,\epsilon}$ is strictly inside the feasible region.

III-C Torque Computation

The control inputs $\mathbf{u}$ consist of contact forces and joint velocities. These commands have to be translated to torques. When using only the feedforward trajectories, desired accelerations and contact forces are extracted and tracked by a hierarchical inverse dynamics controller [24]. When using the feedback policy, we forward simulate the system under the feedback policy for a short time and extract desired accelerations from this rollout. The inverse dynamics is then only used to convert the desired accelerations into torques, without adding additional feedback. This inverse dynamics controller is evaluated at $400\text{\,}\mathrm{Hz}$, while the SLQ-MPC algorithm runs asynchronously on a second onboard Intel i7-4600U@2.1GHz dual core processor.

Finally, each individual motor has a local, embedded, control loop. For the stance legs a torque controller is used, and for the swing legs the motors take the commanded torque as a feedforward term and close the loop over the desired position and joint velocities.

IV Results

We first show the qualitative differences between a feedback policy and a feedforward policy when operating under low update rates and disturbances. Then, the influence of the relaxed barrier function cone constraint on the planned feedback gains is shown. We also examine the structure in the obtained feedback matrices and compare the difference between those obtained with frequency shaping and those without. Finally, we show that the proposed elements, when taken together, lead to a method that can be successfully executed on the onboard hardware of a torque-controlled robot.

We use a diagonal cost on the state and control inputs for all experiments. When frequency shaping is used, we set $\alpha_{i}=0.01$, $\beta_{i}=0.2$ for the contact force inputs and $\alpha_{i}=0.01$, $\beta_{i}=0.1$ for the joint velocity inputs in (17).

IV-A Feedback MPC

We first investigate the effect of low update frequencies and the use of feedback policies from the SLQ-MPC in simulation. We introduce model errors to show the performance of the different strategies. The mass of the control model is increased by $10\text{\,}\mathrm{\char 37\relax}$ with respect to the simulation model. Each MPC controller is brought to a stable trot gait and commanded to move $1\text{\,}\mathrm{m}$ forward.

First, feedforward MPC is used with an update rate equal to the control frequency of $400\text{\,}\mathrm{Hz}$. Every control loop thus has access to the optimal solution from the current state. The resulting desired linear accelerations are shown in the top of Fig. 4. Discontinuities only arise around a contact switch; the desired accelerations are continuous otherwise. In the middle plot, the desired accelerations are shown when the update rate is restricted to an update frequency of $20\text{\,}\mathrm{Hz}$. The updates are clearly visible because every time the feedforward trajectory is updated, the accumulated deviation from the feedforward plan is reset and a new open loop trajectory is tracked. In the bottom plot, we show the performance when the SLQ feedback policy is used. The discontinuities at updates are significantly reduced and smooth trajectories comparable to MPC with high update rate are retrieved.

These experiments show that using the feedback from the MPC can recover some of the performance lost due to lower update rates. By using a policy that is consistent with future MPC updates, the discontinuities at the updates are reduced up to the validity of the linear quadratic approximation.

IV-B Feedback Gains Near Inequality Constraints

In the following experiment we prescribe a task where we require to lift the left front leg and simultaneously set the desired body location towards the front left, outside of the support polygon. This task requires the algorithm to coordinate the step with the body movement. The experiment is performed without the frequency shaped cost to allow us to focus on the effect of using the inequality constraints. In the left side of Fig. 5 we show the optimized solution without inequality constraints. Without the cone constraint, the optimal strategy is to produce negative contact forces in the right hind leg, such that the desired body position can be reached as early as possible. Furthermore, we plot the maximum gain in the row associated with each vertical contact force, which shows that the zero force feedback terms for the swing leg are compliant with the corresponding equality constraint.

The resulting solution when adding inequality constraints is shown on the right of Fig. 5. Here, we used fixed barrier parameters $\mu=0.5$, $\delta=0.1$, under which all constraints are strictly satisfied. With the inequality constraints, the contact forces remain positive on the right hind leg. Additionally, where before the feedback gains are about equal for the three remaining stance legs, now, the feedback gains for the right hind leg are significantly reduced. When the contact force approaches zero, the feedback gains go to zero as well, which ensures that the inequality constraint is not violated after applying the feedback policy at a disturbed state. This interaction between the inequality constraints and the feedback gains is a result of the barrier function. As the constraint boundary is approached, the Hessian of the barrier function w.r.t. the contact forces increases. This increases the input costs $\mathbf{R}$ in (15), and thus reduces the feedback gains.

This gradual decrease in feedback gains cannot be obtained with a clamping strategy, where the gains are unaffected, or with an active set method, where the gains instantaneously decrease to zero when the constraint becomes active [25].

As a secondary effect, we see that the feedback gains on the left front leg increase before lifting the leg. Because the right hind leg is forced to have low feedback gains in the upcoming phase, the body position before lifting the foot is of high importance, which is reflected in the gains.

In addition, we note that the distance to the constraint boundary can be regulated by choosing a different barrier scaling $\mu$. This means that using a finite $\mu$, in contrast to decreasing it to zero as done in an interior point method, can be used to trade some optimality for a larger stability margin.

IV-C Feedback Structure

We visualize the feedback matrix for ANYmal in a full stance configuration. The gains obtained without using the frequency-dependent cost function are shown in Fig. 6a. The state acting on each column of the matrix is shown above the figure, and the control input affected by each row is shown on the left. The color intensity shows the magnitude of each entry in the feedback matrix, with zero shown as white and the highest gain shown in black.

In the joint velocity part of the feedback matrix, one can see how the equality constraints that require zero velocity at the end-effectors are reflected: The joint velocity commands are highly dependent on the linear and angular velocity of the base to achieve this constraint. This empirically verifies that (15) indeed produces feedback matrices that are consistent with the constraints.

The feedback matrix in Fig 6a can be compared to the feedback matrix obtained when using the frequency shaped cost function, shown in Fig 6b. We split the feedback matrix in Fig 6b into four parts, with the vertical split between system and filter state, and horizontal split between system and auxiliary input, corresponding to the partitioning in (19). First, the left side is inspected. Here we recognize a feedback pattern similar to that in Fig 6a, obtained without the frequency-dependent cost function. Indeed since the shaping functions in (17) have unit DC-gain, the feedback matrix $K_{\boldsymbol{\nu},\mathbf{x}}$ is expected to be approximately equal to the matrix obtained without frequency shaping. Furthermore, we can recognize that the direct gains from state to contact forces are an order of magnitude smaller than before.

For the frequency shaped feedback policy, the main feedback flows from system state to auxiliary inputs. The auxiliary inputs then drive the filter states, $\mathbf{x}_{s}$, which in turn provide the adaptation of the system contact force inputs.

For the joint velocities, however, the feedback from filter states to joint velocities is zero, as expected. The fact that end-effector velocities are constrained to be zero is still reflected in the feedback matrix in exactly the same way as before. Feedback terms therefore appear in the bottom right partition of Fig 6b, to satisfy the equality constraint, $\mathbf{u}=\mathbf{C}_{s}\mathbf{x}_{s}+\mathbf{D}_{s}\boldsymbol{\nu}$ , for the rows associated with joint velocities.

This analysis shows that using the frequency-dependent cost function introduces smoothness and reduced direct gains where possible, but at the same time still respects hard equality constraints on the original system inputs.

IV-D Hardware Experiments: Disturbance Rejection

As seen in the accompanying video111A video of the experiments is available at https://youtu.be/KrTrLGDA6FQ, when LQR gains obtained from the SLQ-MPC are used on hardware, the system becomes unstable even in full stance phase. As demonstrated in Section IV-C, the frequency shaped formulation reduces the direct gains between state and input, which enables successful deployment on hardware. All hardware experiments are therefore performed with both frequency shaping and inequality constraints active.

We first perform a simple experiment to verify the qualitative difference observed in simulation between using only a feedforward policy with a conventional tracking controller or when using a feedback policy. The robot is put in a standing configuration with the desired position set to the initial position with zero velocity. Afterward, we place a $5.7\text{\,}\mathrm{kg}$ mass ($\approx 15\%$ of the total mass) on top of the robot to induce a constant disturbance. Fig. 7 shows the resulting desired and measured torque trajectories in the left front knee after the system reaches an equilibrium.

In the top plot, the points at which the feedforward policy is updated are clearly visible. At each update, the state reference is reset to the measured reference, which effectively nullifies the feedback of the tracking controller. Since the feedforward control signal does not account for the additional disturbance, the system deviates from the desired trajectory and builds up feedback in the tracking controller until the next update arrives. In the bottom plot, where the feedback policy updates arrive at the same rate as the feedforward case, there are no discontinuous jumps in desired torque. This verifies our earlier observation in simulation.

IV-E Hardware Experiments: Dynamic Walking

Finally, we demonstrate that the proposed method achieves stable walking with all computations running on the onboard computers. We use a gait known as dynamic walk. The gait pattern is shown in Fig. 8. It consists of a mixture of underactuated and overactuated contact configurations when two and three feet are on the ground, respectively. The proposed friction cone constraint ensures that the trajectory does not require negative contact forces and thus successfully navigates the intriguing pattern of support polygons.

A receding horizon of $1.0\text{\,}\mathrm{s}$ was used, for which the MPC reaches an update frequency of approximately $15\text{\,}\mathrm{Hz}$. The resulting desired and measured torque and joint velocity trajectories are shown in Fig. 9 for the left front leg. The desired and measured signals are close to each other at all time, showing that the applied feedback policy respects the bandwidth limits of the actuators. More dynamic gaits such as pace and trot motions are feasible as well. In the linked video a continuous transition between the two gaits is shown. The feedback MPC is able to skilfully coordinate leg and body motions in these unstable and underactuated situations.

V Conclusion

In this work, we proposed to use feedback MPC as an effective way to handle the slow update rate associated with the computational restrictions of mobile platforms. The sensitivity of CoM control to uncertainties and sampling period has been recently analyzed in [26], which provides a theoretical basis for our observation that stable walking is possible with low update rates.

We proposed a relaxed barrier function method to extend the SLQ algorithm to optimization problems with inequality constraints. In particular, the friction cone is implemented through a perturbed second-order cone constraint. This formulation adds a convex penalty to the cost function and avoids numerical ill-conditioning at the origin of the cone.

A frequency-aware MPC approach was used to systematically include the bandwidth limit of the actuators in the feedback policy design. This was a key factor to achieve closed-loop stability on hardware without any detuning of the low-level actuator controllers as suggested in [8]. The frequency-aware approach effectively allows to set high gains in the low-frequency spectrum and to attenuate gains in high frequency. It thus increases the robustness of the feedback policy in the presence of high-frequency disturbances. We showed that the feedback policy is consistent with the constraints of the locomotion task. We empirically confirmed that the MPC policy reduces the feedback gains near the boundaries of the friction cone to respect the inequality constraints. We also demonstrated that the optimized policy sets zero gains on the contact force of the swing legs and encodes the zero end-effector velocity constraint for stance legs to satisfy state-input equality constraints.