RoLoMa: robust loco-manipulation for quadruped robots with arms

Deployment of robotic systems in the real world requires a certain level of robustness in order to deal with uncertainty factors, such as mismatches in the dynamics model, noise in sensor readings, and communication delays. Some approaches tackle these issues reactively at the control stage. However, regardless of the controller, online motion execution can only be as robust as the system capabilities allow at any given state. This is why it is important to have good motion plans to begin with, where robustness is considered proactively. To this end, we propose a metric (derived from first principles) for representing robustness against external disturbances. We then use this metric within our trajectory optimization framework for solving complex loco-manipulation tasks. Through our experiments, we show that trajectories generated using our approach can resist a greater range of forces originating from any possible direction. By using our method, we can compute trajectories that solve tasks as effectively as before, with the added benefit of being able to counteract stronger disturbances in worst-case scenarios.


I. INTRODUCTION
In this paper, we tackle the problem of robust locomanipulation for quadruped robots equipped with robot arms, such as the one shown in Fig. 1.Here, the challenge is not only to generate whole-body trajectories for solving complex tasks requiring simultaneous locomotion and manipulation (commonly referred to as loco-manipulation), but also to optimize the robustness of such trajectories against unknown external disturbances.This is an important problem for two reasons: first, loco-manipulation allows us to extend the workspace of an otherwise-fixed-base-manipulator through the mobility of a mobile base, such as a legged robot; and second, increasing the robustness of the overall motion against disturbances leads to legged systems that are more reliable, and that can therefore be deployed in the real world with greater confidence.
Enabling loco-manipulation is a very challenging problem because it involves repeatedly breaking and making contacts between the feet and the environment in order to move around, while maintaining balance and avoiding kinematic/actuation limits.Moreover, robust motion planning is a complex subject on its own, as it requires the derivation of good metrics that are able to quantify how robust trajectories are.Consequently, Fig. 1: Snapshots of our robot solving real-world tasks in an industrial setting: (i) turning a hand wheel, (ii) pulling a lever, (iii) opening a gate whilst standing on a ramp, and (iv) lifting a bucket by pulling a rope.The robot and the objects being manipulated have been highlighted for clarity.The overlaid yellow arrows indicate motion.Video footage: https://youtu.be/3qXNHVCagL8.combining loco-manipulation and robustness is not trivial, as it brings together the challenges from both problems.
Most of the previous research done on loco-manipulation, [1]- [3], has tackled the problem by splitting the arm from the base, planning the manipulation separately from the locomotion, and then considering arm movement as a disturbance that the base should compensate for.Furthermore, most of the existing research [3]- [6] has focused on reactive robustness at the control stage, rather than taking into account robustness proactively during planning.
The key components of our approach are (i) a trajectory optimization framework which is able to solve complex realworld tasks and which takes into account the full system dynamics, and (ii) a robustness metric derived from first principles.This allows us to calculate the largest force magnitude that the robot can counteract from any given direction, while considering ground-feet contact stability and the actuation limits of the system.Our results show that, given a contact sequence, our framework is able to plan whole-body trajecto-ries that are significantly more robust to external disturbances compared to other approaches.

A. Statement of Contributions
This paper is a direct follow-up of our previous work.In [7], the framework we proposed was able to optimize wholebody trajectories for standing balance behaviors only, i.e., not actual loco-manipulation.In contrast, this paper proposes an improved formulation which is able to optimize trajectories that involve contact changes-albeit the sequence and timings of those contacts must be prescribed.We are also able to better enforce the nonlinear dynamics of rigid-body systems, as we have incorporated our findings from [8].Furthermore, we show that our formulation can handle cases where contact positions are not enforced explicitly, which actually allows it to further maximize robustness through the adoption of more suitable feet contact positions.Finally, the significant amount of systems integration work that we have done in comparison to [7] allowed us to deploy the robot in a realistic scenario mimicking an industrial offshore platform, where we showed the robot operating uninterrupted and repeatedly solving the complex sequence of tasks highlighted in Fig. 1.This has resulted in a robust loco-manipulation system which is capable of online motion planning for deployment in realistic scenarios.

II. RELATED WORK
We now summarize previous research related to motion planning for quadruped robots equipped with arms, as well as existing research on the topic of robustness.

A. Planning and Control for Quadrupeds with Arms
Murphy et al. [1] were one of the first to investigate the use of a legged robot base to improve the capabilities of a robotic arm.They used trajectory optimization (TO) with a simplified dynamics model to generate open-loop behaviors for Boston Dynamics' robot BigDog.As a result of the coordinated motion between the robot arm and the robot base, they were able to increase the performance of lifting and throwing tasks.In their hardware experiments, they showed the robot dynamically tossing cinder blocks as heavy as 16.5 kg and as far as 4.2 m.However, they only considered standing balance behaviors (which do not change support contacts/move the feet).In contrast, we consider the full dynamics model of the robot during TO, and we also consider behaviors where the robot's feet can make and break contact with the environment.
A few years later, Zimmermann et al. [2] equipped Boston Dynamics' flagship quadruped robot Spot with a Kinova arm in order to perform dynamic grasping maneuvers.Direct control of Spot's actuated joints is not possible because of restricted access to its low-level controller.As Spot's wholebody controller for locomotion compensates for the wrench induced by the manipulator, the tracking and executed motion can differ from planned motion, resulting in poor performance.To achieve precise loco-manipulation, Zimmermann et al. treated Spot's overall behavior as a black box, and built a simplified model of the combined platform from experimental data.While their approach successfully grasped the target most of the times, it failed when the estimated position of the ball was inaccurate or when the robot started executing the planned trajectories from slight offset poses.There were also cases when the robot failed to grasp the target due to Spot's complex internal behavior not being fully captured by their simplified model.For example, when the disturbance induced by the arm to the base was sufficiently large, it could cause a delay that brought the individual trajectory components out of sync.In contrast, our approach does not suffer from this drawback because the robot we use in our experiments grants us full control over its joints, and because we optimize trajectories for the robot considering the dynamics of its whole body.
Ma et al. [3] combined manipulation using model predictive control (MPC) with a locomotion policy obtained from reinforcement learning (RL).First, they modeled the wrenches (arising from the motion of the arm) applied to the base of the robot as external disturbances that can be predicted.Then, they trained the base control policy to counteract those disturbances while (i) trying to keep a horizontal base orientation and (ii) tracking velocity commands from the MPC controller of the arm.In other words, their base policy uses wrench predictions from the arm's motion to compensate for the disturbances applied to the base of the robot.
All of the previous work mentioned thus far [1]- [3] have one thing in common: they all see the arm as a disturbance to be compensated for.In contrast, Bellicoso et al. [9] approached the problem differently, and took into account the dynamics of the whole system.The authors used a whole-body controller based on inverse dynamics, re-planned locomotion continuously in a receding-horizon fashion, and explicitly provided end-effector forces for the controller to track.They equipped ANYbotics' quadruped robot ANYmal B with a Kinova arma combination which results in a fully torque-controlled mobile manipulator, and one which users have full control over.Bellicoso et al. demonstrated that the resulting system is able to perform dynamic locomotion while executing manipulation tasks, such as opening doors, delivering payloads, and humanrobot collaboration.
In [1], the manipulation task is planned offline and separately from the locomotion planner; and in [9], the task for opening the door uses a controller that tracks gripper forces which need to be explicitly specified.Sleiman et al. [6] tackled both these weaknesses when they proposed a unified MPC framework for whole-body loco-manipulation.Their approach augments the dynamics of the object being manipulated to the centroidal dynamics and full kinematics of the robot.This allows the solver to exploit the base-limb coupling and, e.g., to use the arm as a balancing "tail".Despite using an MPC approach, their planner is not adaptive with respect to the dynamic properties of the objects being manipulated.
In our previous work [7], we proposed a motion planning framework for legged robots equipped with manipulators.In that work-and in this paper-we used the same robot arm and quadruped shown in [9], with a difference only in the number of fingers on the gripper.Our motion planning approach was similar to [6] in the sense that we formulated the planning problem for the whole body of the robot in a unified manner, i.e., for the quadrupedal base and the robot manipulator simultaneously.There are a couple of differences in the problem formulation between our previous work [7] and [6], [9]; e.g., we use a full model of the robot's articulated rigid-body dynamics 1 instead of a simplified version, which allows us to plan trajectories more faithfully to the real hardware.However, the most important difference and our main contribution is that we focus on planning motions that are not only physically feasible, but that also maximize robustness against unknown external disturbances.This aspect is something that none of the previous work mentioned [1]- [3], [6], [9] have considered.

B. Motion Robustness Against Disturbances
Del Prete et al. [4] proposed a solution to improve the robustness to joint-torque tracking errors at the control stage.The authors modeled deterministic and stochastic uncertainties in joint torques within their control framework optimization.In our case, we maximize the upper-bound force magnitude the system can withstand from any possible direction, and we do this during the planning stage.Xin et al. [5] proposed a hierarchical controller in which external forces are estimated directly.Their goal was to minimize actuator torques while enforcing constraints for the contact forces.However, in contrast to our work, they did not explicitly enforce actuator limits; and since their main focus is on control, they do not have a planner for computing elaborate whole-body behaviors.
The robot experiments shown in [3] and in [6] demonstrate some capability of resistance against external forces; but this robustness is reactive, in the sense that it is either the learned locomotion policy or the MPC that compensate for the disturbances online in a reactive fashion.In contrast, research on proactive robustness [10]- [12] attempts to take uncertainty into account at the planning stage, i.e., ahead (or just in time) of online execution.A proactive approach allows planning frameworks to increase the system's ability of counteracting disturbances by exploiting kinematic redundancy.Some examples of human kinematics pre-shaping during dynamic tasks are: leaning body backwards while pulling on a rope during a game of tug of war, leaning body forwards and stretching arms while pushing on a broken down car, and a bent-knees stance in preparation for a skateboard trick.Similarly, in robotics, we can increase robustness proactively with offline planning, by e.g.compiling a library of optimal kinematic stances and optimal trajectories, which can then be looked up during online execution.
In [10], Caron et al. proposed a feasible region (a polytope), called Gravito-Inertial Wrench Cone (GIWC), which can be used as a general stability criterion.This representation is very efficient for testing the robust static equilibrium of a legged robot, but it neglects the system's actuation limits.Subsequently, Orsolino et al. [11] proposed to extend the properties of the GIWC by incorporating the torque limits of the system.They demonstrated how the resulting polytopes can be used to e.g.optimize the center of mass (CoM) trajectory in the xyplane for the base-transfer motion of quadrupeds.Orsolino et al. formulated a reduced version of the problem, but even then the technique used to compute polytopes was prohibitively expensive.As a workaround, they computed the polytope only once for the first point of the trajectory, and used that as an approximation for the rest of the motion.We have followed this line of research in our previous work [12], where we proposed a force polytope representation, called residual force polytope, which considers not only the torque limits but also the dynamics of the system during trajectory execution.The polytope is computed from the forces and torques remaining after accounting for Coriolis, centrifugal, and gravity terms, as well as from the nominal feed-forward torques of the motion.
The polytope calculations in [11] and [12] require significant computation time and, in general, deriving explicit descriptions of a projected polytope is NP-hard [13].Zhen et al. [14] formulated a computationally-tractable approach for finding maximally-sized convex bodies inscribed in projected polytopes.Later, Wolfslag et al. [15] adapted that approach for computing the robustness of static robot configurations.Their work was an improvement over exact computations; however, it was still too complex for being considered in trajectory optimization.In our previous work [7], we adapted the technique from [14] to reformulate the problem of computing the smallest unrejectable force (SUF), which allowed us to formulate bilevel trajectory optimization problems for maximizing the robustness of the generated trajectories.In this paper, our work extends those ideas further to motion with contact changes.

A. Robot Model Formulation
We formulate the model of a legged robot in the same way as in our previous work [7], i.e., as a free-floating base B to which the limbs are attached to.For example, the robot we used for our experiments (seen in Fig. 1), has four legs and one arm attached to its base; each leg has three motors and the arm has six. 2 We describe the motion of the system with respect to (w.r.t.) a fixed inertial frame I.We represent the position of the free-floating base w.r.t. the inertial frame, and expressed in the inertial frame, as I r IB ∈ R 3 .As for the orientation of the base, we represent it using modified Rodrigues parameters (MRP) [16], [17] as ψ IB ∈ R 3 .The joint angles describing the configuration of the 6-degrees of freedom (DoF) arm and the four 3-DoF legs are stacked in a vector q j ∈ R nj , where n j = 18.Finally, we write the generalized coordinates vector q and the generalized velocities vector v as where the twist ν B = [ I v B B ω IB ] ⊤ ∈ R 6 encodes the linear and angular velocities of the base B w.r.t. the inertial frame expressed in the I and B frames, and n v = 6 + n j .
2 Please note that our formulation is not tied to the specific robot shown in Fig. 1.In fact, it is general enough such that it can be applied to any legged robot, biped or quadruped, with or without arms.
The equations of motion of a floating-base rigid-body system that interacts with the environment are written as where M (q) ∈ R nv×nv is the mass matrix, and h(q, v) ∈ R nv is the vector of Coriolis, centrifugal, and gravity terms.On the right-hand side of the equation, τ ∈ R nτ is the vector of joint torques commanded to the system, and the selection matrix S = [0 nτ ×(nv−nτ ) I nτ ×nτ ] selects which DoF are actuated.We consider that all limb joints are actuated, thus n τ = n j .The vector λ ∈ R ns denotes the forces and torques 3 experienced at the contact points, with n s being the total dimensionality of all contact wrenches.The support Jacobian J s ∈ R ns×nv maps the contact wrenches λ to joint-space torques, and it is obtained by stacking the Jacobians which relate generalized velocities to limb end-effector motion as , with n c being the number of limbs in contact.Finally, f represents any external force applied to the end-effector.This force may be the result of a push or of some unpredicted disturbance.Under nominal circumstances, this force is zero, i.e., f = 0.The Jacobian J e ∈ R 3×nv is used to map a linear force f applied at the end-effector to joint-space torques.

B. Problem Discretization
In order to plan motions for complex robot systems, we use an approach called direct transcription, which is a powerful technique for TO.
We start by converting the original motion planning problem (which is continuous in time) into a numerical optimization problem that is discrete in time.We divide the trajectory into N equally spaced segments, t I = t 1 < • • • < t M = t F , where t I and t F are the start and final instants, respectively.This division results in M = N + 1 discrete mesh points, for each of which we explicitly discretize the states of the system, as well as the control inputs.Let x k ≡ x(t k ) and u k ≡ u(t k ) be the values of the state and control variables at the k-th mesh point.We treat x k ≜ {q k , v k } and u k ≜ {τ k , λ k } as a set of nonlinear programming (NLP) variables, and formulate the basis of our trajectory optimization problem as where ξ is the vector of decision variables, is the state transition function incorporating the nonlinear dynamics of the system, and X and U are sets of feasible states and control inputs enforced by a set of equality and inequality constraints.The decision variables vector ξ results from aggregating the generalized coordinates q 1:M , generalized velocities v 1:M , joint torques τ 1:N , and contact forces 3 In the equations of motion above, λ represents the contact wrenches, i.e., forces and torques, acting on the robot.However, and similarly to [18], [19], we chose to represent only the force-component as decision variables in our formulation.This is because the quadruped used in our experiments has pointlike feet and we can assume the torque-component is negligible.

C. System Constraints
After having discretized the states and controls of the system over time as decision variables, we need to define a set of rules that restrict the motion represented by those variables.We do this by specifying a set of mathematical equalities and inequalities, so that the solver "knows" how to compute trajectories that are not only physically feasible but that also complete the tasks we want the robot to solve.

1) Domain of decision variables:
The most straightforward constraints we need to write are the lower and upper bounds of each decision variable in ξ.We constrain the joint positions, velocities, and torques to be within their corresponding lower and upper bounds.
2) Initial and final velocities: We enforce the initial and final velocities of every joint to be zero, i.e., v 1 = v M = 0. Note, however, that this is not a strict requirement of our framework but is chosen to ensure static start and end configurations.
3) End-effector poses: We enforce end-effector poses with where f fk (•) is the forward kinematics function, q k are the joint coordinates at the k-th mesh point, i refers to the i-th endeffector of the robot, and p i ∈ SE(3) is the desired pose.We use these constraints for defining the position and orientation of the robot's hand at specific mesh points, as well as to define the point contacts for the robot's feet during stance phases. 4e pre-specify the contact sequence for the feet, which can be computed e.g. using contact planners such as [20].4) Contact forces: For mesh points where the robot is not in contact with the environment (according to the prespecified contact sequences), we enforce the contact forces at the respective contact points to be zero, i.e., λ k = 0.
5) Friction cone constraints: Similarly to [10], we model friction at the contact points using an inner linear approximation with a four-sided friction pyramid.Consider the set of points {C i } where the robot is in contact with its environment.Let n i and µ i be the unit normal and the friction coefficient of the support region at each contact, respectively.A point contact remains fixed as long as its contact force f c i lies inside the linearized friction cone directed by n i : where (t i , b i ) form the basis of the tangential contact plane such that

Required inputs
Footstep locations and timings.

Required inputs
A whole-body trajectory.

Required inputs
Same as in Problem 1.

Constraints
Rigid-body dynamics equations related to the SUF magnitude.

Constraints
Aggregation of constraints from Problems 1 and 2.

Purpose
Compute feasible trajectories with full-body dynamics.

Purpose
Analyze the robustness of already-existing trajectories.

Purpose
Compute trajectories that are robust to external disturbances.

Objective
Minimization of actuator torques and contact forces, or none at all.

Objective
Maximization of SUF magnitude.

Objective
Same as in Problem 2.
Fig. 2: Block summaries of NLP Problems 1, 2, and 3. Inside each block, the summary states the purpose for using that formulation, the required inputs, the decision variables and constraints involved, and the objective function employed.The empty lines with a '-' under 'Decision variables' emphasize that NLP Problem 3 is a combination of Problems 1 and 2.

6) System dynamics:
We enforce the equations of motion (Equation 2) using inverse dynamics.This decision comes from our recent research on the benefits of inverse dynamics vs. forward dynamics for eliminating dynamics defects in direct transcription formulations.In brief, we have shown that problems formulated using inverse dynamics are faster, more robust to coarser problem discretization, and converge in fewer iterations-see [8] for more details.
The inverse dynamics problem computes the joint torques and forces required to meet desired joint accelerations at a given state, i.e., where f id (•) is the function that solves the inverse dynamics problem, and the desired joint accelerations can be calculated implicitly with v * k = (v k+1 − v k )/h.We compute q * k+1 from v k+1 , and integrate it (using semi-implicit Euler integration) to compute the next generalized coordinates q * k+1 .Finally, we define the dynamics defect constraints as D. Robustness Against Disturbances 1) NLP Problem 1: Thus far, we have modeled the robot and its full body dynamics, discretized the motion planning problem, and defined a set of rules in the form of mathematical constraints.At this stage, we have all the "ingredients" required for planning feasible trajectories that can be executed on the robot.Henceforth, we will refer to this version of the formulation as NLP Problem 1-a summary for this version of the formulation is shown on the left block in Fig. 2.
Next, we are going to build upon our previous work [7] to present two different problem formulations, NLP Problem 2 and NLP Problem 3, which can be used to analyze the robustness of known trajectories and to maximize the robustness of trajectories being computed, respectively.
2) NLP Problem 2: When we compute a robot trajectory using NLP Problem 1, we may be interested in understanding how robust those trajectories are against forces applied to e.g. the end-effector.Thus, one way of determining the robustness of that trajectory is by studying the set of forces that the endeffector is able to resist, both in terms of force magnitude and direction.The metric proposed in our previous work [7], the smallest unrejectable force (SUF), represents the smallest force magnitude (applied from any possible direction) that the robot is not able to counteract.In other words, it gives the magnitude of the largest force that the robot can counteract in a worst-case scenario.Next, we explain how to compute it.
Given a discretized robot trajectory (e.g., the output of NLP Problem 1), we can compute the SUF magnitude throughout that motion by re-formulating the NLP problem.The decision variables for such a problem are where, for each k-th mesh point, ρ k is the magnitude of the SUF and K λk is the instantaneous gain matrix mapping a force expressed in end-effector space to ground-feet contact space.Each and every ρ k is bound to R + , i.e., K λk have no explicit bounds; but they are constrained by with a λ and b λ pertaining to the alternative form of writing the friction cone constraints, i.e., A λ λ ≤ b λ . 5Akin to the equations of motion, the relationship between the SUF and the robot capabilities at every mesh point is given by S ⊤ K τ + J ⊤ s K λ + J ⊤ e ρ = 0.And this can be rewritten in a way that highlights the inherent structure of the constraint: where K τ is an instantaneous gain matrix mapping a force expressed in the end-effector frame to joint-torque space.We enforce this relationship by splitting it into two parts.For the top part of the equation, concerning the floating base, we write the following nonlinear equality: As for the bottom part, concerning the limbs of the robot, we write the following nonlinear inequality: Once we have defined the above constraints and decision variables, we can maximize the following objective function with any off-the-shelf nonlinear solver: In summary, for a given (constant) trajectory, the outcome of this nonlinear optimization problem will be the magnitude of the SUF over time (i.e., the ρ k for every mesh point of the discretized trajectory) and K λk .K τ are also an output, since they can be computed as a function of K λ and ρ without performing any inversion-as hinted by Equation 17.The outputs K λ and K τ can be used to understand and explain the specific constraint that determines the upper bound of the SUF (i.e., friction cone or torque, on which foot or motor), although that is something we do not investigate in this paper.
3) NLP Problem 3: NLP Problems 1 and 2 allow us to compute a feasible whole-body trajectory and to calculate the robustness of said trajectories, respectively.By combining NLP Problem 1 and NLP Problem 2, we obtain a single (albeit more complex) problem formulation which is able to compute whole-body trajectories that are not only feasible but also more robust against external disturbances.We call this single formulation NLP Problem 3, and it is summarized on the right block in Fig. 2. The decision variables of NLP Problem 3 are (21) 5 We refer readers to [7] for a full explanation and derivation of the terms.
The constraints are the combined constraints of NLP Problems 1 and 2. The objective function is the same as NLP Problem 2, i.e., the maximization of every mesh point's SUF added up:

E. Contact Switching
In contrast to our previous work [7], the NLP formulations in Fig. 2 consider the making and breaking of contacts between the feet of the robot and its environment.This is one of the main contributions of this paper, and we will now explain how we have enabled this.
As we have previously explained in 'System Constraints' (subsection III-C), during problem discretization, we handle feet that are in contact with the ground (stance phase) differently from feet that are moving through free space (swing phase).In short, for each mesh point of the trajectory, and for each foot of the robot: if that foot is in stance phase, we enforce it to remain still and we also enforce the contact force to lie within friction cone boundaries; otherwise, we enforce only a zero contact force constraint.Since this is done during the problem transcription process, it requires knowing the timings for contact switches a priori, as well as feet positions 6for each stance phase.The upside is that the problem complexity does not increase much, especially compared with other approaches that consider contact switching by formulating complementarity constraints (e.g.Posa et al. [21]).
The approach explained in the last paragraph is enough for enabling contact switching in NLP Problem 1.However, additional changes are needed for Problems 2 and 3, since those problem formulations involve extra constraints and decision variables regarding the maximization of the SUF.In essence, the decision variables K λk pertaining to the feet in swing phase must be set to zero, since no contact forces exist in that context.Fortunately, we need not worry about K τ k , since it is defined as a function of K λk (see Equation 17).Nonetheless, the most challenging aspect of these modifications is not in the writing of the constraints, but rather in the writing of the functions that evaluate the Jacobian of those constraints (together with their sparsity structure).To this end, we use modern automatic differentiation capabilities from Julia [22], but the process of passing the results of the Jacobian and sparsity structure evaluations to specialized NLP solvers remains tricky and requires particular attention by the programmer to the indexing of the decision variables and constraints.We formulate large but sparse NLP problems through our framework and therefore specifying the sparsity pattern is very important to attain shorter computation times, as the Jacobian of the constraints contain mostly zeroes (especially the dynamics Jacobian, which is block diagonal).

F. Computation Time
It is not trivial to formulate or solve the optimization problems described in this paper.Special care needs to be taken when transcribing the constraints and objective functions, otherwise computation becomes prohibitively expensiveespecially when considering full body dynamics models.Here, we show that our approach is computationally feasible as an offline planning tool suitable for industrial applications.
The following evaluations were carried out in a singlethreaded process on a MacBook Air (2020) with an Apple M1 chip and 8 GB of unified memory.Our framework has been implemented in Julia [22], using the rigid-body dynamics library RigidBodyDynamics.jl (RBD.jl)[23], and the Artelys Knitro [24] optimization library.We used the interior-point method of Waltz et al. [25] to solve the NLP problems.Automatic sparsity detection was used to provide sparsity patterns to the solver, and automatic differentiation (forward mode) for evaluating the Jacobians of the constraints-but we allowed Knitro to compute the Hessians via the limitedmemory quasi-Newton BFGS method.In Table I, we present the elapsed time (in seconds) for solving the three NLP problems summarized in Fig. 2, for an industrial task where the quadruped robot must rotate a wheel/valve using the arm mounted on its torso. 7The average solving time was determined by solving each problem multiple times with different initial conditions.These conditions involved random translations and rotations of the robot's base relative to the manipulated object.
Based on the table, it is evident that tackling these problems involves a high computational cost, making our approach impractical for online planning.However, the computation times are still quite reasonable to consider using our framework offline, enabling us to create a repository of behaviors that can be efficiently accessed during online execution.In relation to the size of the problems, employing a 100 Hz discretization over a trajectory duration of 4 s yields a total of 401 knot points.Table II shows the number of decision variables, constraints, and nonzeros in the Jacobian for that amount of knot points, offering some insight into the complexity of problems of such size.

IV. EXPERIMENTS AND RESULTS
We now present the experiments we carried out for evaluating our work and their respective results.This section is organized as follows: A. Describes the system integration work required for combining the quadruped robot and the robot arm, as well as our planning framework with existing controllers; B. Presents a repeatability test where we commanded the robot to turn an industrial hand wheel multiple times from different starting positions and orientations; C. Compares the robustness of two distinct trajectories for turning a hand wheel and for pulling a lever; D. Explains how the SUF can be used as a tool for analyzing existing robot trajectories; E. Demonstrates the capabilities of our framework to plan robust whole-body motions involving making and breaking of contacts; and finally, F. Tests the robustness of a loco-manipulation trajectory for lifting a bucket with incremental weights until failure.Subsections A, B, C, and F are concrete evaluations of our method on real robot hardware.Whereas the focus of subsections D and E is to demonstrate new features.

A. System Integration
In order to demonstrate the capabilities of our planning framework on the real robot, we integrated it with the existing software stacks of the ANYmal quadruped and Kinova arm.
For the quadruped, we used one of ANYbotics' software releases which comes with multiple locomotion controllers working out-of-the-box.Human operators can control the robot remotely via a joystick to tell the robot where to walk and which gait to use, or they can pre-specify a mission as a set of waypoints for the robot to walk through.Moreover, they provide an interface for specifying a custom payload attached to the robot, but this is assumed to be a static payload, such as a thermal camera, or an imaging sensor.
For controlling the robot as a whole (quadruped + arm), we settled on two operation modes: teleoperated and autonomous.When the robot is in teleoperated mode, we do not move the Kinova arm; it remains still in a "parked" configuration (stowed on top of the quadruped).Doing this allows us to consider the arm a static payload, which we can specify using the interface provided by ANYbotics.In turn, this allows us to take advantage of every existing capability provided by ANYbotics' stack.Finally, whenever we wish to operate the arm, we switch to autonomous mode.In this mode, we use a custom controller 8 for the quadruped base and Kinova's velocity controller for the arm.This is the mode we use to execute the whole-body trajectories generated with our planning framework.
In summary, we start the robot in teleoperated mode by default.In this mode, we can walk the robot on flat ground and over ramps, but we cannot move the arm.We use the teleoperated mode to walk the robot around a facility and towards points of interest.Then, once the robot reaches said points of interest, it awaits an instruction (e.g, "turn the wheel", or "pull the lever").As soon as the robot is given an instruction, it switches into autonomous mode, our planning framework is triggered and computes a whole-body trajectory for the robot.This trajectory is then passed on to the arm and quadruped controllers for synchronous execution.Once the robot finishes executing the task autonomously, it switches back into teleoperated mode and the system returns control to the human operator, who can walk the robot remotely towards the next task of the mission.
A demonstration of the entire system is available here: https://youtu.be/3qXNHVCagL8.In the video, a human teleoperates the robot to walk around a mock-up scaffolding of an offshore platform (such as an oil rig).The operator approaches different points of interest, and commands the robot to autonomously turn a hand wheel, pull a lever, push a gate whilst standing on a ramp, and pull a rope to lift a 1.1 kg bucket.These tasks are the ones shown in Fig. 1.Although we have focused on tasks relevant for industrial inspection, our framework allows us to formulate virtually any locomanipulation task through the definition of a set of constraints for the robot.

B. Repeatability Test
Our goal for this experiment was to ensure the following characteristics of our system: • the robot is able to operate continuously for extended periods of time without falling; • the operator is able to send walking commands to the robot during teleoperated mode, but not during autonomous mode; • the Vicon motion capture system calculates the relative transform between the free-floating base of the robot and the manipulation target reliably; • the planning framework computes a feasible whole-body trajectory for completing the manipulation task; • the controller is able to reach and grasp targets, and track reference trajectories accurately.To verify the points above, we carried out a repeatability test for turning an industrial hand wheel.During this test, we operated the robot continuously for 30 min without a safety harness.During this time, a human operator walked the robot to different points near the manipulation target and triggered the "turn wheel" behaviour.The test was completed successfully and our system passed all of the points listed above.Fig. 3 shows a few snapshots of the robot grasping the wheel during the test.Video footage is also available here: https://youtu.be/Ok8Pcwn_I0w.

C. SUF Optimization for Turning a Wheel and Pulling a Lever
In this experiment, our goal is to compare the resulting SUF of two different objective functions: the first trajectory (baseline) is optimized considering a cost function that minimizes torques and contact forces, while the second trajectory (proposed ) maximizes the SUF explicitly.We run this experiment for two tasks: (i) turn a hand wheel clockwise by a full revolution, and (ii) pull down a lever from its resting position.Fig. 4 shows an instance of the results of this experiment.
1) Turning the hand wheel: In Fig. 4's left column, we can see that the SUF magnitude of the proposed trajectory (orange line) is within 73 N to 76 N throughout the entire motion; whereas the SUF of the baseline trajectory (blue line) lies within 56 N to 66 N-it starts and ends at ∼56 N, increasing slightly in the middle of the trajectory where it peaks at ∼66 N. The SUF mean-percentage-increase of the proposed approach versus the baseline is approximately 24 %.Video: https://youtu.be/1M32AHuuDhI.
2) Pulling the lever: In Fig. 4's right column, the SUF magnitude of the proposed trajectory (orange line) remains within 68 N to 75 N throughout the entire motion; on the other hand, the SUF magnitude of the baseline trajectory (blue line) remains within 50 N to 69 N. The SUF meanpercentage-increase of the proposed approach versus the baseline is approximately 18 %.Video: https://youtu.be/6A9eSdfcj7A.
We repeated this experiment multiple times for different robot orientations relative to the wheel and the lever.The results were identical to those shown in Fig. 4.Both the baseline and proposed trajectories successfully completed the task.However, trajectories planned using the proposed approach outperformed the baseline version in both tasks in terms of their robustness against external disturbances.

D. The SUF as a Tool for Analyzing Trajectories
Earlier in this manuscript, when we presented Fig. 2, we explained how the SUF allows us to analyze the robustness of existing trajectories, and also how it allows us to compare multiple trajectories with each other.In this experiment, we want to further develop the idea of using the SUF as an analysis tool.Our goal is to show that, thanks to our problem formulation, we can analyze individual trajectories to understand e.g.how important of a role each of the legs play during a motion.
In Fig. 5a, we show a robot trajectory planned with our framework for picking up something from the ground; and in Fig. 5b, the solid blue line represents the SUF over time for that motion.An important question we may ask is: how important is each leg for a successful execution of the motion?We can ask a similar question from a different perspective: if there is a hardware fault on any of the quadruped legs, how much of the initial motion's robustness remains?Next, we explain how we can answer these questions.Thanks to the way we formulated the SUF constraints, we can set individual terms in K λk to zero in order to "simulate" what would happen if the torques of the motors of a specific leg of the robot were at their torque limits, or what would happen if the contact force at a specific foot was on the boundary of the friction cone-or both.We solved Fig. 2's NLP Problem 2 four times (one for each leg) with the goal of analyzing the SUF at the end-effector while considering the motor torques to be at their limits and the ground-feet contact forces to be at the boundaries of their respective friction cones.The dashed lines plotted in Fig. 5b show the results.
We can see that the magnitude of the SUF decreases, regardless of which leg is affected.We can also see that the dashed lines are slightly different from each other, which is expected since the trajectory is not perfectly symmetric and each leg has a slightly different load from every other leg.Moreover, we can actually identify which leg plays the most important role in this motion by looking at the dashed line with the lowest values (the one for which the SUF magnitude decreased the most)-it was the hind right (HR) leg.

E. SUF Optimization with Making and Breaking Contacts
In this experiment, we wanted to verify that our planning framework is able to optimize trajectories involving switching contacts.We also evaluated whether the proposed objective function is able to compute a trajectory that is more robust to external disturbances than the baseline cost.
We defined a task in which the robot starts with all four feet on pre-specified positions and then raises its right hind foot off the ground for a short time.We also constrained the gripper to remain still at a certain location for the duration of the whole motion.We then solved the task in two different ways: with NLP Problem 1 (baseline) and with NLP Problem 3 (optimized).The plot in Fig. 6 shows the magnitude of the SUF over time for the resulting trajectories.Because NLP Problem 1 only outputs a trajectory, we passed it to NLP Problem 2 afterwards in order to compute the SUF valueskeep in mind that this did not change the baseline trajectory.The planner was able to compute a feasible trajectory for solving the task using both approaches.However, as we can see in Fig. 6, the optimized trajectory is more robust than the baseline.Moreover, the plot shows that the breaking and making of contacts directly influences the magnitude of the SUF.For both trajectories, the SUF magnitude decreases when the right hind foot breaks contact with the floor; and then increases again as the foot reestablishes contact.Video: https://youtu.be/H6-g8NLGyYE.

F. Robustness Test with Incremental Weights
In this last (but central) experiment, we plan a robust locomanipulation trajectory for pulling a rope attached to a bucket and execute it on the real robot.Our goal is to show that the extended version of our framework can handle problems that require simultaneous locomotion and manipulation while maximizing robustness against disturbances.Planning loco-manipulation trajectories where robustness is considered proactively was not possible before, and this is really what we have been trying to tackle with our work in this paper.
The experimental setup (shown in Fig. 7) consisted of a rope threaded through a pulley.On one end, the rope was attached to a bucket on the ground; and on the other end, the rope was attached to a handle on a platform.The task for the robot was to grasp the handle from the top of the scaffolding and then pull the rope to lift the bucket.Importantly, in order to lift the bucket high enough, the robot had to take a few steps back while simultaneously pulling the rope with its arm.The position of the handle w.r.t. the robot was calculated with the Vicon motion capture system (similarly to the other tasks), and the footsteps were pre-specified using a static crawl gait.
To test the robustness of the trajectory planned with our framework, we asked the robot to lift the bucket in multiple trials, wherein the weight of the bucket was incremented 1 kg in each trial (and the bucket was empty on the first trial).The bucket itself weighs approximately 1.1 kg.
The robot completed 4 successful trials, in which it lifted 1.1, 2.1, 3.1, and 4.1 kg.On the 5th trial (5.1 kg) the robot was able to lift the bucket momentarily, but then the handle slipped off the gripper and we counted this as a failure.A video of the experiment is available: https://youtu.be/puy2S90_3CM, and Fig. 7 shows some snapshots of the last successful trial (i.e., trial 4), where the robot lifted the bucket containing three 1 kg weight plates.Finally, Fig. 8 shows the magnitude of the SUF over time for this task, where we can see the dips corresponding to the intervals when the robot was taking steps.
This experiment showed that our robot was capable of reliably executing a loco-manipulation trajectory planned with our framework.Moreover, it showed that the robot is indeed capable of dealing with external disturbances while executing a trajectory.It is worth remembering that for this experiment (as well as for the other experiments shown in this paper) we did not model the payload; in other words, the system does not know about the bucket or how much it weighs a priori (during planning).Instead, the weighted bucket is acting as a disturbance, and the robot's ability to lift the bucket is therefore a direct outcome of our robustness metric and the feedback terms used in our controller.

V. CONTACT LOCATION OPTIMIZATION
Throughout this work, we pre-specified the position for each foot of the robot during the loco-manipulation planning stage.For example, when the human operator commands the robot to turn the hand wheel, our framework calculates the feet positions of the robot at that instant via forward kinematics, and then the mathematical constraints of our optimization problem ensure the robot's feet remain on those positions while the wheel is being turned.But this raises an important question: when we optimize a trajectory that maximizes robustness, if we constrain the feet to specific positions, will that not limit how robust the resulting trajectory is?In other words, perhaps the resulting trajectory could be more robust if another set of feet locations had been chosen.
In light of this, we set out to investigate an approach for choosing more adequate contact locations within our motion planning framework.We now consider the additional problem of optimizing the continuous location of the feet contacts.Our goal was to understand whether our planner is capable of adapting feet positions such that the resultant whole-body trajectory can be more robust.
In order to enable optimization of feet locations, we extended our NLP formulation.For that, we expand the vector of decision variables of the NLP problem to include the xycoordinates of each foot for the beginning of each individual stance phase.We added the xy-coordinates of each foot to our problem as decision variables, assuming that the robot would stand on flat ground (i.e., we assume the z-coordinate for each foot is zero).We also modified Equation 8to consider those decision variables, since previously their right handside (i.e., p i ) were a pre-specified constant position for each foot.The extended NLP formulation is more complex than the one previously used throughout this work: it has additional decision variables (some of them coupled 9 ), and more complex constraints for enforcing the contact positions; however, it provides more flexibility to the solver, which should now be able to compute feasible trajectories that further maximize robustness by adapting foot contact locations.
To test this more-flexible NLP formulation, we defined a motion planning task in which the robot must maintain its configuration while reaching a target point in task-space with its end-effector.Then, we solved three slightly different versions of the optimization problem: Nominal -In this version of the problem, we fixed the configuration of the robot base to its default joint positions (defined by the manufacturer).We did not constrain the configuration of the Kinova arm.Therefore, the solution to this version of the problem consists of the robot mov- 9 Feet positions are represented explicitly by the new xy-coordinates, but also implicitly by the forward kinematics of the robot's configuration q.
ing only its arm to reach for the target and then holding that configuration.No objective function is provided, so this is a feasibility problem.Baseline -In this version, the solver is able to change the configuration of the whole-body of the robot.This is similar to the baseline formulation used in previous sections, but the feet locations are now decision variables.The torques and contact forces over time are minimized.Proposed -In this version, the solver is also able to change the configuration of the whole-body of the robot.This is similar to the proposed formulation used in previous sections, but the foot locations are now decision variables.The magnitude of the SUF over time is maximized.
Next, we show the results of these slightly different problem versions.The robot configurations of each version are shown in Fig. 9.The feet positions of each solution are shown in the plot of Fig. 10.After solving each version of the problem, we then computed the magnitude of the SUF at the gripper.The computed SUF spheres are also shown in Fig. 9, and labeled with their respective magnitude.
The first thing we can observe from these results is that the solver did take the liberty of optimizing the feet locations.This verifies that the solver is able to deal with our more-complex NLP formulation, albeit taking longer than the previous formulation.Secondly, the plot of the feet locations clearly shows that the trajectory optimized with the baseline approach has a smaller support polygon than the trajectory optimized with the proposed approach.This is relevant because we know that the support polygon is a good representation for the region where the CoM projection can lie for achieving static balance.However, the support polygon is only an approximation since it does not take into account robot capabilities (torques and friction at the contacts); e.g., it is not guaranteed that the robot has enough actuation power to maintain a pose whose CoM projection lies on one of the corners of the support polygon.This leads to our next observation: while the baseline approach keeps the CoM position close to the nominal version, the CoM of the trajectory optimized with the proposed approach lies further from the center of the support polygon.Our metric converges to this configuration because it takes into account the robot capabilities, and is able to find a more stable and robust pose, even though its CoM projection does not lie at the center of the support polygon.Finally, we can see that the worst-case disturbance scenario that the robot can resist is better for the proposed approach, and worse for the nominal scenario.As labeled in Fig. 9, the magnitude of the SUF for the nominal, baseline, and proposed approaches were 58 N, 57 N, and 76 N, respectively.The above results show that our formulation is able to optimize feet locations, and that our metric is able to guide the solver to solutions that have increased capabilities of resisting external forces-at least in theory, that is.However, when we deploy our method on the robot, there are always practical subtleties that may affect how the robot behaves, such as signal delay or the type of controller being used.With that in mind, we carried out an experiment to assess the actual force-rejection capabilities of the real robot.Next, we explain the experiment setup and then analyze the results.
In this experiment, we executed the trajectory optimized with each approach on the real robot and, for each case, we disturbed the robot by pulling its end-effector from three different angles.A summary of the events that occurred for each angle/trajectory is given below: value measured by the force gauge.For each row, i.e., for each trajectory, the label highlighted in red corresponds to the magnitude of the SUF found experimentally for the three distinct angles used for the experiment.The way to interpret these results is to take the smallest SUF found experimentally for each trajectory (highlighted in red and signifying the worst case perturbation), and then compare those magnitudes to find the trajectory with largest worst case SUF (indicating the highest robustness).
As highlighted in Fig. 11, the proposed trajectory exhibited the highest SUF (60 N).In percentage terms, the SUF of the proposed trajectory was approximately 33 % better than both nominal and baseline trajectories, which is a significant improvement.Moreover, the SUF magnitudes found experimentally were 21 % to 22 % smaller than the SUF predicted originally (the ones shown in Fig. 9).We understand this is due to model mismatch and unaccounted factors when we deploy the trajectory on the real robot.But the fact that the percentage decrease is similar for each trajectory is reassuring, as it tells us that the unaccounted factors of the real robot affected the SUF equally, and the proposed trajectory performed better than the nominal and baseline trajectories as we had predicted in relative terms.

VI. DISCUSSION AND FUTURE WORK
In this paper, we presented a framework for planning wholebody loco-manipulation trajectories robust to external disturbances.We integrated our framework with existing software stacks to enable easy switching between teleoperated and autonomous modes.We demonstrated the capabilities of that integration by having a human operator remotely control the robot via a joystick in a mock-up rig of an industrial site, and by having the robot autonomously plan complex and rich whole-body motions for real-world tasks within that setting, such as turning a hand wheel, pulling a lever, opening a gate whilst on a ramp, and lifting a heavy bucket by pulling a rope.We also carried out a wide range of experiments to test the reliability of the full system, analyze the SUF of existing trajectories, optimize the robustness of trajectories (including those involving making and breaking of contacts).Finally, we carried out an initial investigation on the possibility of using our framework for the purpose of optimizing feet positions.
Our experiments showed that the resulting system was reliable and versatile.The robot was able to navigate the industrial scaffolding and successfully execute all the challenges we prepared.Moreover, we showed that the SUF is a valuable metric, both for analyzing existing trajectories and for maximizing the robustness of new trajectories.Finally, our preliminary research on the optimization of feet locations showed promising results.We now conclude with a list of interesting avenues for future work.

A. Taking Into Account Force-Feedback During Execution
In this work, we have disregarded the dynamics model of objects being manipulated by the robot.We maximized robustness against external disturbances at the end-effector in order to demonstrate that we do not have to necessarily model the object dynamics (which may not be available at planning time), as the controller is able to track the reference trajectory and compensate for the object dynamics through feedback terms at the controller level.However, in extreme cases, controller feedback terms are not enough to execute the desired motion appropriately.
An immediate extension of our work would be to take into account force-feedback during execution of the task.We could plan the motion just like we have done throughout this paper before starting the actual manipulation task.However, once the robot starts executing the task, the object being manipulated will exert a set of forces that in turn apply torque at the robot's joints.An interesting direction for future work would be to take into account those torques at the joint level to estimate the force being applied to the robot as a consequence of the robot-object interaction, and then replanning the motion taking into account that force estimation.This would mean that we would still be able to maximize the SUF, but we would now have a much better model of the task taking place.

B. Decreasing Computational Time for Real-Time Execution
For many real-world tasks, such as the ones we showed with our robot, it is sufficient to take the current state of the world, plan a motion in "one-shot" for completing the task, and then executing that trajectory.However, there are scenarios where this approach might fail, especially if the state of the world changes as the robot executes the taskmeaning that the planned trajectories become invalid.Realtime control schemes (such as MPC) are usually employed for dealing with those scenarios.Currently, the computation times of our framework for planning robust trajectories are too slow to be compatible with the budget available in those real-time control schemes.
There are a few possible approaches for tackling this challenge, if we want to use our framework for such scenarios.In the first approach, we can try to decrease the computational time e.g. by making the robot model simpler, making the problem smaller by considering a more coarse discretization, or by decreasing the complexity of the problem constraints.Another approach would be to change the NLP solver used for solving the trajectory optimization problems.E.g., we could attempt to formulate our robustness optimization problem using a DDP-based approach, such as [26], instead of using a direct transcription approach that relies on off-theshelf solvers.However, the disadvantage of this approach is that, while off-the-shelf commercial solvers are able to deal with a wide range of constraints, enforcing general constrains with DDP-based approaches is pretty much an active research topic of its own.Finally, a third possible approach would be to speed-up the computation of the robustness metric itself.Perhaps there is a way the metric can be learnt, represented by a surrogate model, or stored in a look-up table for quick retrievals.

C. Whole-Body Robustness Analysis
In this work and in our previous work [7], we focused on improving the robustness to external disturbances applied at the end-effector.This was because we were interested in applications where the robot had to interact with objects without modeling the dynamics of those objects.However, the mathematical derivation of the SUF applies to any point on any rigid body of the robot mechanism.In other words, instead of computing the SUF at the end-effector, we could also compute the SUF at the robot's base, at the knees, at the elbows, and so on.An interesting path for future work would be the development of a software for analyzing robot trajectories where a robot and a trajectory are given as inputs, and then as an output we could see the SUF spheres at multiple points along the robot's kinematic chain.This kind of visualization would allow us to better understand and have a very clear visual representation of bottlenecks and weaknesses in trajectories.With such a tool, we would be able to look at a trajectory, and easily indicate that at a certain time, the robot is susceptible to failure even with small disturbances applied to specific parts of its body.

D. Robust Footstep Optimization
Another interesting direction for future research is the possibility of finding contact locations by taking into account the robustness metric.We presented preliminary work on this topic in the previous section.

APPENDIX
In Table III, we list all the videos supplementing our work, in the same order as they appear in this manuscript.A playlist is also available here: shorturl.at/oFJU0.

Fig. 3 :
Fig. 3: Snapshots of the robot grasping the hand wheel during the repeatability test.Notice how the position and orientation of the robot base relative to the wheel are different on every snapshot.Video: https://youtu.be/Ok8Pcwn_I0w.

Fig. 4 :
Fig.4: From top to bottom: initial configuration of the baseline trajectories (top), initial configuration of the proposed trajectories (middle), and plot comparing the SUF of the baseline and proposed trajectories (bottom).Left and right columns concern the tasks of turning an industrial hand wheel and pulling down a lever, respectively.

Fig. 5 :
Fig.5: The image on the top shows a motion plan for picking up an object from the floor.The plot on the bottom shows, for each leg, how the SUF decreases if we were to consider the current torques to be at actuator limits and current contact force at the friction cone boundary.

Fig. 6 :
Fig. 6: Plot of the SUF magnitude over time for two trajectories solving the same task specification.Blue shows the baseline version (solved with NLP Problem 1), and orange shows the optimised version (solved with NLP Problem 3).

Fig. 7 :Fig. 8 :
Fig.7: Snapshots of the robot lifting a bucket with weight plates, weighing 4.1 kg in total.This task requires the quadruped base to take steps at the same time as the arm moves (loco-manipulation).Despite the dynamics of the heavy bucket not being modeled, the robot is still able to complete the task, thanks to the robustness of the planned motion and to the feedback gains of the controller.

Fig. 9 :Fig. 10 :
Fig. 9: Robot configurations resultant from the extended NLP formulation which allows the solver to optimize feet locations.The radius of the sphere corresponds to the magnitude of the SUF.

Fig. 11 :
Fig.11: Snapshots of the experiment in which we disturbed the robot by pulling its end-effector.Pictures on the same row correspond to the same trajectory, whereas pictures under the same column share the same disturbance angle.The snapshots were taken from the video footage collected during the experiment, at the time the force gauge indicated the greatest force magnitude-which is shown in the label on the bottom right corner of each picture.Labels highlighted in red are the worst-case scenario for each row (out of the three angles shown).The five red circles on each snapshot represent the initial positions of the four feet and of the end-effector.Video: https://youtu.be/tUXQUqLneTE.

TABLE I :
Average computation time (in seconds) for solving the valve task at different discretization rates (50 Hz and 100 Hz).

TABLE II :
Characteristics of the NLP Problems for a 4 s-long trajectory discretized at 100 Hz (valve task).