Locomotion Planning through a Hybrid Bayesian Trajectory Optimization

TL;DR

This paper introduces a Bayesian optimization approach for locomotion planning in legged robots, enabling learning of contact schedules from high-level task descriptors to improve planning efficiency and alignment with human intuition.

Contribution

The paper presents a novel bi-level optimization framework using Gaussian processes to learn contact schedules, reducing computational complexity in motion planning.

Findings

Successfully applied to single-legged hopping scenarios

Learns contact transitions that match human intuition

Performs competitively against heuristic methods

Abstract

Locomotion planning for legged systems requires reasoning about suitable contact schedules. The contact sequence and timings constitute a hybrid dynamical system and prescribe a subset of achievable motions. State-of-the-art approaches cast motion planning as an optimal control problem. In order to decrease computational complexity, one common strategy separates footstep planning from motion optimization and plans contacts using heuristics. In this paper, we propose to learn contact schedule selection from high-level task descriptors using Bayesian optimization. A bi-level optimization is defined in which a Gaussian process model predicts the performance of trajectories generated by a motion planning nonlinear program. The agent, therefore, retains the ability to reason about suitable contact schedules, while explicit computation of the corresponding gradients is avoided. We delineate the algorithm in its general form and provide results for planning single-legged hopping. Our method is capable of learning contact schedule transitions that align with human intuition. It performs competitively against a heuristic baseline in predicting task appropriate contact schedules.

Equations

$$\underset{\mathbf{x}}{\min }$$

$$\displaystyle\hat{\mathbf{y}}\mkern-5.0mu=\mkern-5.0mu\operatorname*{argmin}\limits_{\mathbf{y}}\big{\{}c_{\text{LL}}(\mathbf{x},\mathbf{y})|\mathbf{g}_{\text{LL}}(\mathbf{x},\mathbf{y})=\mathbf{0},\,\mathbf{h}_{\text{LL}}(\mathbf{x},\mathbf{y})\leq\mathbf{0}\big{\}},$$

$${\mathbf{g}}_{\text{UL}}\left(\mathbf{x},\mathbf{y}\right)=0,$$

$${\mathbf{h}}_{\text{UL}}\left(\mathbf{x},\mathbf{y}\right)\le 0,$$

$$\underset{\mathbf{y}}{\min }$$

$${\mathbf{g}}_{\text{LL}}\left({\mathbf{x}}_0,\mathbf{y}\right)=0,$$

$${\mathbf{h}}_{\text{LL}}\left({\mathbf{x}}_0,\mathbf{y}\right)\le 0,$$

$$M\left(\mathbf{x},\mathbf{y}\right):=$$

$$+{\sigma }_3\sum _{j=1}^m{\left(\max \left(0,h_{LL,j}\left(\mathbf{x},\mathbf{y}\right)\right)\right)}^2,$$

$$M\left(\mathbf{x}\right)\sim \mathcal{G}\mathcal{P}\{\mu \left(\mathbf{x}\right),k\left(\mathbf{x},{\mathbf{x}}^{\prime }\right)\},$$

$$\mu \left({\mathbf{x}}^{\ast }\right)$$

$${\sigma }^2\left({\mathbf{x}}^{\ast }\right)$$

$${\mathbf{s}}_k=\underset{\mathbf{s}\in \mathbf{S}}{\mathrm{argmin}}\left({\mu }_{k-1}\left({\mathbf{z}}_k,\mathbf{s}\right)-{\beta }_k^{1/2}{\sigma }_{k-1}\left({\mathbf{z}}_k,\mathbf{s}\right)\right),$$

$$\mathbb{E}{\left[{\left(\frac{\Delta \mathbf{m}}{\mathbf{m}}\right)}^2\right]}_k$$

$$\mathbf{K}\left(\mathbf{x},{\mathbf{x}}^{\prime }\right)={\mathbf{K}}_Z\left(\mathbf{z},{\mathbf{z}}^{\prime }\right)\circ {\mathbf{K}}_S\left(\mathbf{s},{\mathbf{s}}^{\prime }\right),$$

$$\mathbf{q}$$

$$\mathbf{u}$$

$$\mathbf{M}\left(\mathbf{q}\right)\ddot{\mathbf{q}}+\mathbf{b}\left(\mathbf{q},\dot{\mathbf{q}}\right)+\mathbf{g}\left(\mathbf{q}\right)={\mathbf{S}}^{\top }\mathbf{u}+{\mathbf{J}}_{\text{c}}^{\top }{\mathbf{\lambda }}_{\text{c}},$$

$${\mathbf{J}}_{\text{c}}\ddot{\mathbf{q}}+{\dot{\mathbf{J}}}_{\text{c}}\dot{\mathbf{q}}=0,$$

$$\left[\begin{array}{cc}\mathbf{M}& -{\mathbf{J}}_{\text{c}}^{\top }\\ -{\mathbf{J}}_{\text{c}}& 0\end{array}\right]\left[\begin{array}{c}\ddot{\mathbf{q}}\\ {\mathbf{\lambda }}_{\text{c}}\end{array}\right]=\left[\begin{array}{c}{\mathbf{S}}^{\top }\mathbf{u}-\mathbf{b}-\mathbf{g}\\ {\dot{\mathbf{J}}}_{\text{c}}\dot{\mathbf{q}}\end{array}\right].$$

$$\begin{array}{r}{\mathbf{x}}_{\text{hm}}=\left[x_{\text{hm,1}},\dots ,x_{\text{hm},n_{\text{S}}}\right],\\ {\mathbf{z}}_{\text{hm}}=\left[h(x_{\text{hm,1}}),\dots ,h(x_{\text{hm},n_{\text{S}}})\right].\end{array}$$

$$\mathbf{Z}$$

$$\mathbf{S}$$

$${\mathbf{S}}_{\text{Baseline}}=\{\begin{array}{l}\left[3,0,0,0,0\right],\\ \left[4,3,5,0,0\right],\\ \left[5,4,6,0,0\right],\\ \left[4,3,3,3,4\right],\\ \left[5,4,3,4,6\right].\end{array}$$

Full text

I Introduction

Deployment of legged systems into real-world scenarios requires autonomous locomotion planning strategies. A challenge in legged locomotion is the necessity to reason about contact schedules during planning. We are interested in the efficient automation of contact schedule selection based on high-level descriptors of a motion task.

The literature provides a wide spectrum of approaches to locomotion planning problems. The separation into simplified models for footstep planning together with robust whole-body control for tracking has found widespread application [1, 2]. In [3], footstep plans were optimized using a mixed-integer formulation incorporating a reachability heuristic. In [4], the stochastic Covariance Matrix Adaptation Evolution Strategy (CMA-ES) [5] optimized footholds on a pre-labeled terrain costmap. Sampling-based methods employing Rapidly-exploring Random Trees (RRT) have solved complex motion planning problems while being guided by heuristics and requiring post-processing [6, 7, 8]. The strict separation of footstep planning and whole-body tracking requires a conservative approach for the planning of the prior task to ensure feasibility of the latter task.

In dynamic locomotion tasks, inconsistent contact planning can entail performance degradation or even instability. The motion planning should thus reason about the robot dynamics and the contact states concurrently. A continuous-time differential dynamic programming (DDP) based method was provided in [9], while [10, 11] proposed direct collocation-based approaches that optimize over the contact schedule while planning dynamically consistent motions. These algorithms optimize over the contact duration and can implicitly alter the number of contact phases by setting step durations to zero. However, this seems unlikely from a numeric standpoint. Furthermore, including the contact schedule in the optimization variables increases both the computational burden and the risk to not converge to a feasible solution. Similar numerical concerns arise in contact invariant optimization [12, 13], where contact interaction is implicitly defined by complementary constraints. This raises the question of whether contact schedule selection can be decoupled from trajectory optimization, without sacrificing model consistency in favor of heuristics.

Learning contact schedule selection to guide the trajectory optimization presents itself as a viable option.

Techniques combining optimal control with learning have found widespread interest in the robotics community. Optimization is efficient in finding good local solutions, while data-driven techniques offer the ability to model complex relations. In [14], controller adaptation to model inaccuracies was framed as a reinforcement learning problem. In [15] and [16], robustness of neural network policies was demonstrated for locomotion tasks. In [17], a classifier predicts contact mode sequences for a manipulation task under linearized dynamics to simplify the quadratic program solved online. The strongly nonlinear locomotion tasks considered here make the optimization prone to become locally trapped or to not converge at all. Efficient data acquisition should, therefore, avoid extensive sampling in regions likely to produce infeasible problems. This exploration-exploitation trade-off can be achieved through a combination of Bayesian Optimization (BO) with Gaussian Process (GP) models, as discussed in [18, 19] and applied to contextual bandit problems in [20]. Previously, promising results with these approaches were demonstrated for applications such as automatic locomotion controller optimization in [21, 22].

In this paper, we propose an efficient method to automate contact schedule selection in optimization-based legged locomotion. The planning problem is cast as a bi-level optimization.

At the lower level (LL), a direct collocation method solves the motion planning problem for the state and input trajectories under a fixed contact schedule. At the upper level (UL), BO is used to optimize a GP model which maps the contact scheduling policy to the trajectory optimization performance. The temporal component of gait planning is therefore decoupled from optimizing continuous state and input trajectories. Employing a learning-based method in the UL allows for sample-efficient offline training that avoids predicted infeasible actions, which can slow down computation significantly [23]. Unlike in [24], we furthermore avoid taking gradients with respect to the contact schedule and reduce the risk of getting trapped in local optima. We evaluate our approach on planning trajectories for the single-legged hopping robot shown in Fig. 1.

II Methods

Optimization-based locomotion planning attempts to compute optimal state and input trajectories, as well as the underlying contact sequence and durations. This proceeds in accordance with the system dynamics and the environment, as contact imposes constraints on the achievable motions. In general, operating on these continuous and discrete variables defines a mixed logic optimization problem [25], which is computationally expensive to solve. Here, we employ a bi-level optimization that separates gait planning (discrete) from trajectory planning (continuous) [24]. At the LL, we define a nonlinear program (NLP) to find optimal state and input trajectories under fixed gaits. At the UL, we leverage a gradient-free method capable of treating the NLP score as a black box function. The overall problem takes the form

[TABLE]

where $\mathbf{x}$ are the UL decision variables, $\mathbf{y}$ are the LL decision variables, $c_{i}$ denotes the cost, $\mathbf{g}_{\text{i}}$ are the equality constraints, and $\mathbf{h}_{\text{i}}$ are the inequality constraints, where ${i\in\{\text{LL},\text{UL}\}}$.

LL optimization

At the LL we consider the gait to be given and set ${\mathbf{x}=\mathbf{x}_{0}}$. We then determine the optimal state and input trajectories for the resulting locomotion task. To this end, we formulate an optimization problem using a direct collocation method [10]. The resulting NLP is written in its general form as

[TABLE]

where ${c_{\text{LL}}=f_{\text{O}}}$ is the LL objective function. The LL decision variables vector ${\mathbf{y}=\left[\mathbf{y}_{1}^{\top},\dots,\mathbf{y}_{N}^{\top}\right]^{\top}}$ stacks the individual node vectors ${\mathbf{y}_{i}=\left[\mathbf{q}^{\top}_{i},\dot{\mathbf{q}}^{\top}_{i},\mathbf{u}^{\top}_{i}\right]^{\top}}$ consisting of generalized positions $\mathbf{q}$, generalized velocities $\dot{\mathbf{q}}$, and input torques $\mathbf{u}$. The constraints include system dynamics, ground contact, friction, and joint limits. The NLP in (2) returns the optimized trajectory nodes $\mathbf{y}_{\text{opt}}$ together with the associated cost $f_{\text{O,opt}}$ and constraint values $\mathbf{g}_{\text{LL,opt}}$ and $\mathbf{h}_{\text{LL,opt}}$.

Merit function

In order to quantify the quality of solutions provided by the LL optimization, we define a merit function $M\left(\mathbf{x},\mathbf{y}\right)$. The merit function weighs trajectory cost and constraint violation based on

[TABLE]

where the first term penalizes the trajectory cost, while the second and the third terms penalize constraint violation of the equality and inequality constraints, respectively. In general, their ratios should ensure a separation between feasible and infeasible runs, while avoiding steep gradients in the merit function.

UL optimization

At the UL we consider the gait to be variable and select $\mathbf{x}$ such that $c_{\text{UL}}$ is minimal. Here, we set ${c_{\text{UL}}=M\left(\mathbf{x},\mathbf{y}\right)}$ to convey that an optimal gait should lead to both minimal cost and constraint violation. The UL optimization is constrained by the LL solution such that ${M\left(\mathbf{x}\right)=M\left(\mathbf{x},\mathbf{y}=\text{NLP}\left(\mathbf{x}\right)\right)}$. We therefore choose a gradient-free method from BO to solve it. To this end, we model the merit function as a GP according to

[TABLE]

where $\mu\left(\mathbf{x}\right)$ is a mean function and $k\left(\mathbf{x},\mathbf{x}^{\prime}\right)$ is a kernel function. The kernel function serves as a similarity measure and relates distance in pairing space ${\mathbf{d}_{\mathbf{x}}=\mathbf{x}-\mathbf{x}^{\prime}}$ to distance in merit space ${d_{m}=m-m^{\prime}}$. We denote the set of observed samples by $\{\mathbf{X},\mathbf{m}\}$, where $\mathbf{X}$ is the matrix of observed UL decision variables and $\mathbf{m}$ the vector of corresponding merit scores, and a new query by $\{\mathbf{x}^{*},m^{*}\}$. The GP model then uses the observations to infer predictions for new queries, while providing estimates of the associated uncertainty. It does so via the predicted mean $\mu$ and predicted standard deviation $\sigma$, defined as

[TABLE]

where $\mathbf{K}$ is the resulting kernel matrix and $\sigma_{n}$ encodes the noise-level in the observations.

Context and actions

Following the nomenclature in [20], we encode planning tasks using context $\mathbf{z}$. The context here includes the distance to the goal location together with the terrain heightmap. The available gaits are referred to as actions $\mathbf{s}$ which encode a contact sequence and the associated contact durations. We define a context-action pair as ${\mathbf{x}=\mathbf{z}\times\mathbf{s}}$. These parameterize the UL optimization and subjected to $\mathbf{z}$ we try to select $\mathbf{s}$ such that $M\left(\mathbf{x}\right)$ is minimal.

III Bayesian Optimization Algorithm

The concepts introduced in the previous section are combined into a single BO algorithm. In the following, we introduce the GP regression formulation with the associated termination criterion and further specify both the kernel function and the merit function employed. The corresponding pseudo-code is provided in Algorithm 1.

GP regression

We employ GP regression to sequentially solve the UL optimization in (1). Specifically, we use a variation of the contextual Gaussian Process upper confidence bound algorithm (CGP-UCB) presented in [20]. At iteration $k$, context $\mathbf{z}_{k}$ is sampled uniformly at random from context space $\mathbf{Z}$ and the posterior distribution is evaluated over the entire action space $\mathbf{S}$. Action $\mathbf{s}_{k}$ is then selected using upper confidence bound (UCB) sampling according to the acquisition function

[TABLE]

where $\mu_{k-1}$ and $\sigma_{k-1}$ denote the posterior mean and standard deviation of the previous iteration, and $\beta_{k}=\log\left(k\right)$ is an exploration-exploitation trade-off variable that encourages ongoing exploration. Based on the resulting $\mathbf{x}_{k}$, we solve the NLP and compute a merit score $m_{k}$ from the optimized trajectory. These are then added to the set of observed samples and the GP posterior belief is updated. This process continues until a termination criterion is satisfied.

Termination criterion

We assert convergence of the regression by monitoring a discounted version of the relative prediction error. Specifically, we use a first-order lowpass filter on the squared relative residuals (fSRR). At iteration $k$, the expectation of the fSRR is computed as

[TABLE]

where $k$ indicates the current iteration index, $m_{k}$ is the observed merit score, and $\Delta m_{k}$ is the corresponding difference between the observation and the prediction. The scalar $\rho$ defines how much emphasis is placed on the more recent predictions in computing the fSRR. Data acquisition is terminated once the fSRR falls below a pre-defined threshold provided as

[TABLE]

Kernel function

The GP model requires a kernel function $k\left(\mathbf{x},\mathbf{x}^{\prime}\right)$ as a similarity metric. We follow [20] and define one kernel over context space and another kernel over action space. Our overall kernel is then a multiplicative combination of the two, encoding that samples are only similar if both their context and their action are similar:

[TABLE]

with kernel hyperparameters $\mathbf{\sigma}_{Z}$ and $\mathbf{\sigma}_{S}$. For the underlying kernel functions, we choose ARD-Matern3/2 kernels, as these have limited smoothness and therefore have some capability of accommodating potential steps in the sample data. Such discontinuities can arise during transitions from one contact schedule to another, or when one contact schedule suddenly becomes infeasible.

Merit refinement

Failed trajectories can lead to large constraint penalties in the merit function. Such outliers should be avoided, as GP models presume some smoothness characteristics of the functions they attempt to approximate. We therefore subject the merit function $M$ in (3) to a sigmoid function to arrive at the new merit function ${M^{\prime}=\tanh(M)}$. The sigmoid function acts as a soft range limiter and ensures that the merit scores have an upper bound, while the transform is continuous and retains the merit score order.

IV Modelling

We evaluate the performance of the proposed algorithm on a simulated robot. In the following, the robot model and the environmental model are discussed.

IV-A System

The agent considered is a single-legged hopper, modelled in the sagittal plane (Fig. 1). It consists of a base (B), thigh, shank and actuated joints at the hip (H) and knee (K). The generalized coordinates $\mathbf{q}$ and joint torques $\mathbf{u}$ are

[TABLE]

where $x_{i}$ is a horizontal position, $z_{i}$ is a vertical position, $\varphi_{i}$ is a joint angle, and $u_{i}$ is the corresponding joint torque. The associated equations of motion (EoM) are written as

[TABLE]

where ${\mathbf{M}\in\mathbb{R}^{4\times 4}}$ denotes the inertia matrix, ${\mathbf{b}\in\mathbb{R}^{4}}$ groups Coriolis and centrifugal effects, ${\mathbf{g}\in\mathbb{R}^{4}}$ is the gravitational force vector, ${\mathbf{S}\in\mathbb{R}^{2\times 4}}$ is the selection matrix, ${\bm{\lambda}_{\text{c}}\in\mathbb{R}^{2}}$ refers to the ground reaction force, and ${\mathbf{J}_{\text{c}}\in\mathbb{R}^{2\times 4}}$ is the corresponding Jacobian.

The system dynamics are represented as a hybrid system by projecting the rigid body equations into the null space of the contact constraints. During flight, no external force is present and equation (13) can be solved to obtain the generalized accelerations. During stance, friction limited ground contact occurs at the point foot. The associated contact constraint is derived to yield

[TABLE]

which is used to solve (13) simultaneously for the accelerations and contact force. The stance dynamics are

[TABLE]

The resulting contact forces are furthermore constrained in our trajectory optimization (TO) to lie in the friction cone defined by the local terrain normal.

IV-B Environment

The environment is encoded via a heightmap representation. Terrain height $h$ is sampled along the horizontal with discretization $\Delta x$ on the interval ${x\in[x_{\text{min}},x_{\text{max}}]}$. For a total of $n_{\text{S}}$ samples, the resulting heightmap is

[TABLE]

A $C^{1}$- continuous representation of the terrain is recovered by leveraging shape-preserving piecewise cubic interpolation. Terrain height and gradient at intermediate locations are then computed from the interpolation polynomials.

V Experiments

Our main interest is in whether the GP model can learn reasonable contact schedule transitions for a given task. We are furthermore interested in whether the resulting predictor can perform competitively against a heuristic baseline. In the following, we introduce the general experimental setup and provide results for planning motions on both flat ground and uneven terrain.

V-A Preliminaries

The following describes the general setup used for both the flat ground scenario and the uneven terrain scenario. First, the representation of context and action is provided. Then, the algorithm used for assessing model performance is outlined.

V-A1 Context-action space

The input space of the regression is spanned by the context space $\mathbf{Z}$ and the action space $\mathbf{S}$. A context ${\mathbf{z}\in\mathbf{Z}}$ consists of a goal distance and an $n_{\text{T}}$-dimensional vector of terrain features, denoting varying heightmap values. We consider goal distances of up to {z_{\text{D}}=1.0\mathrm{m}} and terrain height magnitudes of up to {|z_{\text{T}}|=0.2\mathrm{m}}, such that

[TABLE]

An action ${\mathbf{s}\in\mathbf{S}}$ consists of an $n_{\text{p}}$-dimensional vector of alternating contact/flight phase durations. As we employ collocation nodes with constant time length, phase durations are encoded by the number of nodes. We consider up to two jumps, ${n_{\text{p}}=5}$, and phase node numbers in the range of ${n_{\text{N}}=\left[3,6\right]}$, as these encourage dynamic motions with lower jerk. We then have

[TABLE]

The UL optimization determines the desired gait by setting certain phases to zero. For example, a single stance phase is encoded by action ${\mathbf{s}_{\text{A}}=\left[n_{1},0,0,0,0\right]}$, a single jump by ${\mathbf{s}_{\text{B}}=\left[n_{1},n_{2},n_{3},0,0\right]}$, and a double jump by ${\mathbf{s}_{\text{C}}=\left[n_{1},n_{2},n_{3},n_{4},n_{5}\right]}$. Our nomenclature requires ${\mathbf{s}_{\text{D}}=\left[n_{1},0,0,n_{4},n_{5}\right]}$ to be encoded by ${\mathbf{s}_{\text{B}}}$.

V-A2 Performance metric

Performance of a trained GP model is assessed on a set of ${n_{\mathbf{z}}}$ contexts sampled uniformly at random. For each context, the model predicts an optimal action and the merit score is determined by running the optimization. The procedure is outlined in Algorithm 2.

V-B Scenario 1: flat ground

The scenario of hopping over flat ground is considered first. We refer to this flat terrain as FT with ${n_{\text{T}}=0}$, i.e., the terrain is not part of the context for this experiment.

V-B1 Learning contact transitions

We trained a GP model on flat terrain via Algorithm 1 and refer to it as GP-FT. Fig. 2 provides the fSRR termination criterion over the number of iterations. It is apparent that the error metric is steadily decreasing, which means that the predicted merit scores become more accurate over time. Discontinuities in the fSRR indicate larger updates in the posterior belief, which can happen when exploring new parts of the context-action space. Fig. 3 provides an overview of the learned contact sequence transitions over different goal distances, while Fig. 4 displays a selection of learned contact schedules. From Fig. 3 it is evident that the GP model is capable of learning natural contact sequence transitions that progress from stance over a single jump to a double jump. Referring to the contact schedules in Fig. 4, the model furthermore learned to vary contact durations within the contact sequences. The overall trend of adapting the phase durations or adding steps to the contact schedule as goal distance increases matches what human intuition would suggest.

V-B2 Performance

We compare the performance of our learned model to a heuristic baseline. The baseline runs the trajectory optimization for a fixed set of $5$ contact schedules and selects the best action based on the resulting merit scores. The individual schedules were selected before model training based on intuition. Initial and final stance phases were chosen sufficiently long to initiate and terminate dynamic hopping. Intermediate stance phases are shorter to keep motions fluid. Flight phases are short to limit energy expenditure. The action-set available to the baseline is

[TABLE]

The performance of the GP-FT model was evaluated by competing against the baseline according to Algorithm 2 over ${n_{\mathbf{z}}=100}$ runs. The results are provided in Fig. 5 which illustrates the merit scores of different goal distances for the baseline and GP-FT approaches. In general, the trained GP-FT has a better performance than the baseline method. At the phase transitions, performance seems to converge with the baseline having a slight edge at the stance to single jump and single jump to double jump transition.

V-C Scenario 2: uneven terrain

The scenario of hopping over uneven terrain is considered next. Fig. 6 shows some of the randomly generated terrains used for training of the GP. Here, we vary $3$ terrain nodes to create obstacles. Thus, in the following, we use ${n_{\text{T}}=3}$ features as terrain context to describe the heightmap.

V-C1 Learning contact transitions

We trained a GP model on randomly generated rough terrain via Algorithm 1 and refer to it as GP-RT. This proceeded in two stages. First, the model was trained on flat ground until convergence to an intermediate fSRR value. This initializes the model with unperturbed samples from the underlying system dynamics. Then, the model was trained on rough terrain until convergence to a final fSRR value. At each iteration, a terrain was randomly sampled from a Gaussian distribution.

V-C2 Performance

The performance of the GP-RT model was evaluated by competing against the GP-FT model according to Algorithm 2 over ${n_{\mathbf{z}}=1000}$ runs. The results are provided in Fig. 7. The merit scores were smoothed using a moving average filter for interpretability. The GP-RT model incurs lower merit scores than the GP-FT model on average. This can in part be attributed to the GP-RT model selecting actions that fail less frequently than those selected by the GP-FT model. Here, GP-FT generated motions fail in about $6\%$ of the evaluated runs, while GP-RT motions fail only in $3\%$ of them. Note that we restrict our analysis to a maximum of two jumps and randomly sampled context may generate truly infeasible objectives (e.g., Fig. 8 bottom with the maximum goal distance {z_{\text{D}}=1.0\mathrm{m}}).

In general, the GP-FT model acts aggressively and selects very dynamic motions, or short contact schedules, at the risk of failing. The GP-RT model acts more conservatively and adapts its actions to the observed terrain to avoid failing. For this reason, we observed the GP-RT actions to perform slightly worse than the GP-FT actions on terrains with limited roughness. However, GP-RT motions demonstrate an increased level of robustness on terrains with more severe roughness. This behavior can be attributed to the underlying probabilistic model used for modeling the context (i.e., terrain heightmap). In the GP-RT case, the optimized gait policy always assumes some level of uncertainty over the observed terrain. To avoid failure the GP-RT policy chooses gait sequences which are less prone to uncertainties such as premature TO convergence or falling into a local minimum.

An example of this robust behavior is provided in Fig. 8. Here, the GP-RT model selects a double jump schedule, while the GP-FT model selects a single jump schedule. On the smaller obstacle (top), the GP-RT performs slightly worse than the GP-FT. On the higher obstacle (bottom), the GP-RT performs significantly better than the GP-FT. This is because the lower level NLP for the latter does not converge. Learning the terrain features is therefore advantageous in selecting contact schedules competitively.

VI Conclusion

In this paper, we proposed a method that learns to select contact schedules based on high-level task descriptors. The problem is cast as a bi-level optimization, where contact planning proceeds in the UL and constrained TO in the LL. The performance of a trajectory resulting from a specific contact schedule task is quantified using a merit function. The merit function is modeled as a GP and contact schedule selection is learned via a BO approach.

It was shown that the GP model is capable of learning contact schedule transitions that align with human intuition. The trained model is capable of outperforming a heuristic baseline in predicting task appropriate contact schedules. During training, the GP model learns about both the underlying system dynamics and the interaction with the specific terrain. It hereby does not necessarily find the global optimum of the strongly non-convex problem, but learns to select contact schedules that the NLP solver performs well with.

It was demonstrated that a model trained on rough terrain outperforms a model trained on only flat terrain, highlighting that the method can incorporate terrain information into the decisions. However, we were impressed by how close the model trained on flat terrain comes to the rough terrain performance. It shows that our bi-level formulation provides a certain robustness. Even when the UL does not modify its decision to the terrain, the LL optimizes with full terrain information and manages to find reasonable solutions. The main difference is seen in scenarios with extreme terrain features, where the TO cannot converge with the flat terrain gait while the rough terrain gait succeeds in the task.

Future work will include the extension to longer gait sequences over more complex terrains. Moreover, we wish to use the method in a 3D setting as well. For both aspects, the scalability of the GPs to the larger input space has to be addressed, as exact GP regression has ${\mathcal{O}\left(n^{3}\right)}$ complexity. A promising strategy would be to employ a Deep Neural Network to imitate the GP model and thus to achieve better scalability at the expense of longer training sessions.