- Silicon Synapse Lab., Department of Electrical and Computer Engineering, Northeastern University, Boston, MA, United States

Fast constraint satisfaction, frontal dynamics stabilization, and avoiding fallovers in dynamic, bipedal walkers can be pretty challenging. The challenges include underactuation, vulnerability to external perturbations, and high computational complexity that arise when accounting for the system full-dynamics and environmental interactions. In this work, we study the potential roles of thrusters in addressing some of these locomotion challenges in bipedal robotics. We will introduce a thruster-assisted bipedal robot called *Harpy*. We will capitalize on Harpy’s unique design to propose an optimization-free approach to satisfy gait feasibility conditions. In this thruster-assisted legged locomotion, the reference trajectories can be manipulated to fulfill constraints brought on by ground contact and those prescribed for states and inputs. Unintended changes to the trajectories, especially those optimized to produce periodic orbits, can adversely affect gait stability and hybrid invariance. We will show our approach can still guarantee stability and hybrid invariance of the gaits by employing the thrusters in Harpy. We will also show that the thrusters can be leveraged to robustify the gaits by dodging fallovers or jumping over large obstacles.

## 1 Introduction

Raibert’s hopping robots Raibert et al. (1984), and Boston Dynamic’s BigDog Raibert et al. (2008) are amongst the most successful examples of legged robots, as they can hop or trot robustly even in the presence of significant unplanned disturbances. Other than these successful examples, many bipedal and anthropomorphic robots have also been introduced. Boston Dynamics’ dynamic humanoid, ATLAS, has pushed the limits of dynamic legged locomotion with its 28 hydraulically actuated joints. This robot has showcased impressive mobility feats, including jumping over obstacles and dynamic flip-turns. Unfortunately, there are no publications laying out the details of Boston Dynamics works. Therefore, little is known about the controllers used in ATLAS.

Another successful example is Cassie. Agility Robotics developed this bipedal robot based on an earlier prototype called ATRIAS led by Oregon State University Apgar et al. (2018). The biped can negotiate unstructured environments such as ramps and stairs inside buildings robustly and efficiently. With 20 Degrees of Freedom (DOF) and ten actuators, Cassie is a hard-to-control floating base possessing 6 Degrees of Under-actuation (DOU). With a smaller number of DOU than Cassie and ATRIAS and a restricted frontal dynamics, Michigan’s MABEL possesses pogo-stick-style feet, compliant legs, and an anthropomorphic morphology. This robot has shown stable and natural running gaits similar to humans Sreenath et al. (2011). Completely blind and relied on no visual feedback, MABEL has showcased stable gaits, which involve traversing along rough terrains even when no information about the bumps’ whereabouts in its path is available Park et al. (2012).

Fully actuated systems such as Valkyrie, ASIMO, Mahru, and Yobotics-IHMC have shown impressive performance similar to the above examples. NASA’s Johnson Space Center led a team of partners from academia and industry and developed NASA’s first bipedal humanoid, Valkyrie Paine et al. (2015). Valkyrie has showcased the successful completion of sophisticated human-style tasks in the DARPA Robotic Challenge (DRC). Other examples such as Honda’s ASIMO Hirose and Ogawa (2006) and Samsung’s Mahru III Kwon et al. (2007) have demonstrated capabilities such as quasi-static walking, running, dancing, and climbing stairs inside buildings. Or, Yobotics-IHMC lower body, humanoid biped has shown recovery from severe pushes Pratt et al. (2009).

Despite all of these accomplishments, state-of-the-art bipedal robots are prone to fall-over and cannot negotiate highly rough terrains. Even humans known for their natural, efficient, and robust locomotion feats cannot recover from unpredictable situations such as severe external pushes, scuffing, or slippage on icy surfaces.

Our main goal is to enhance bipedal systems’ robustness through distributed arrays of thrusters and nonlinear control. This paper will report our recent efforts in dynamic modeling and designing closed-loop feedback for the thruster-assisted locomotion of a legged robot called Harpy, which is shown in Figure 1. Currently, Harpy’s hardware is being developed at Northeastern University (NU). This bipedal robot is equipped with a total number of eight custom-made joint actuators, a pair of coaxial thrusters fixed to Harpy’s torso, and a light body structure fabricated out of reinforced carbon fibers.

**FIGURE 1**. Illustration of *Harpy’s* CAD model. Harpy is a bipedal robot with a pair of thrusters and is being developed by the authors to study robust, efficient and agile thruster-assisted legged locomotion.

The overarching goal is to integrate aerial and legged systems’ capabilities in a single platform. Aerial robots possess fast mobility mainly because they can fly over obstacles. This property of aerial systems can be highly suitable for applications such as automated package delivery or monitoring and inspection from vantage points. While legged robots may have difficulty negotiating extremely bumpy terrains, e.g., semi-collapsed buildings in the aftermath of an earthquake, these systems maintain a superior energetic efficiency of locomotion compared to aerial drones.

Unlike aerial robots, legged systems can safely operate near humans and animals and are very suitable for many applications such as search and rescue inside buildings, digital agriculture, monitoring livestock for disease symptoms, assisting workers in construction zones, etc. In all of these applications, the fast-rotating blades in rotary-wing robots can result in severe laceration injuries. As a result, there are strict rules on aerial systems’ operation near humans.

From a control design perspective, thruster-assisted legged locomotion has not been explored previously. In this work, we are particularly interested in the closed-loop feedback design implications of thrusters, and our objective is to apply thrusters to achieve two distinguished goals. We will briefly discuss these goals below.

### 1.1 Disturbance Rejection in the Sagittal and Frontal Planes of Locomotion

One objective is to use the thrusters and robustify the gaits. Bipedal platforms are incredibly vulnerable to disturbance and fall-over. Environment and model disturbances are common contributors to these systems’ tip-over instability (rotational and postural) Goswami and Kallem (2004), which can lead to complete loss or severe damages to the system. Because of the inherent complexity of disturbance rejection, a handful of mainstream strategies have been successfully applied to solve such a vulnerability Gubina et al. (1974); Pratt and Tedrake (2006); Kajita and Tani (1991); Hill and Fahimi (2015); Mandava and Vundavilli (2018); Song et al. (2018); Arcos-Legarda et al. (2019). None of these strategies have been able to mitigate tip-over scenarios entirely.

Rotational stability can be achieved if the external forces and moments acting on the system yield zero centroidal moments, i.e., the system’s overall angular momentum is conserved Vukobratović and Stepanenko (1972). External forces and moments may arise from various sources, such as gravity, immature ground contacts, and unexpected external disturbances. To retain rotational stability the orientation of the Ground Reaction Force (GRF) should be regulated relative to the Center of Mass (COM).

Bipedal robots cannot directly control the magnitude and orientation of the GRF because of the unilaterality of the contact forces. However, they can modulate the GRF indirectly through the nonlinear couplings in these systems. The unilaterality of GRF yields a considerably limited tip-over recovering and disturbance rejection capability. For instance, when the extent of the disturbance force is moderately small, e.g., a gentle push, small body movements can be leveraged to retain the robot’s rotational balance.

Upper body movements can be applied to induce restoring moments at the hips and ankles by transferring angular moments. Traditional forward and inverse kinematics frameworks that place virtual components such as dampers and springs in strategic locations in the robot have been successfully applied to describe these robot-environment interactions Pratt et al. (2001). Still, their use is limited to minor disturbances.

Other widely applied methods based on Center of Pressure (COP) control are limiting. Only limited to quasi-steady gaits, these frameworks locate and regulate the COP of a bipedal robot within the support polygon through various strategies, e.g., applying ankle torques, to oppose the tip-over moment. That said, the effectiveness of these methods is limited. The movement of the COP position along the feet, which is proportional to the torque generated at the ankle, can end to the boundaries of the support polygon, and tip-over can occur. As a result, the support polygon’s size is often interpreted as the decision (stability) margin in the state space of the COM. The position of the COM relative to the support polygon is constantly monitored to decide whether the COP-based controller is enough to retain stability or not Pratt et al. (2006).

More aggressive disturbance rejection mechanisms such as a flywheel effect allow for additional momentum contribution required to retain a bipedal robot’s stability against more prominent external disturbances, e.g., a hard impulsive push force applied to the system Stephens (2007). However, the flywheel method’s effectiveness can be limited mainly because significant resisting torques must be applied to the joints to terminate the flywheel effect after the disturbance is attenuated, which can yield other issues such as the violation of contact forces at leg-end points. Larger disturbances can be handled by taking a few compensating steps, just like humans retaining their balance under sudden pushes. While this approach is very promising, its effectiveness is highly affected by other limitations in the step kinematics and impact dynamics (the topology of terrain).

### 1.2 Dynamic Walking and Performance Satisfaction

Another objective is to explore performance satisfaction paradigms superior to existing methods. We will briefly explain this objective here.

Dynamic bipedal walkers possess small support polygons. Small stability margin and underactuation at the contact points result in severe challenges to achieve energetically efficient and robust gaits. These two antagonistic properties are often achieved at the cost of sacrificing one another, and satisfying efficiency and robustness at the same time is still an open problem. Examples such as the robust locomotion feats of ATLAS robot achieved through the application of energy-hungry hydraulic actuators and the efficient gaits of Cornell’s passive walker with a minimum number of electric actuators are two notable examples that fall on the two opposite ends of the performance spectrum.

The problem of simultaneously providing asymptotic stability, optimizing desired performance indices, and satisfying gait constraints in legged systems has been studied extensively Westervelt and Grizzle (2007); Galloway et al. (2015); Dai and Tedrake (2016); Feng et al. (2014). However, the challenges remain, and there are severe limits with existing methods. For instance, model-based methods such as Hybrid Zero Dynamics (HZD), which have offered rigorous ways to assign attributes such as a minimized total cost of transport and robustness in an off-line fashion, are restrictive for dynamic scenarios involving joint trajectory re-planning. Or, in work by, despite increased cost-of-transport optimization accuracy without additional computational time using a variable-time-interval trajectory optimization method, obtaining optimal solutions remain restricted to simple models and slow gaits.

To better understand existing options, we categorize performance satisfaction in dynamic systems – regardless of being a legged or non-legged robot (e.g., robot manipulators) – into three broad categories. Namely: 1) trajectory optimization (TO), 2) optimization-based (nonlinear) controls, and 3) reference trajectory manipulation.

TO problems’ goal is to generate optimal trajectories that satisfies constraints on states, inputs, and GRF while ensuring that the trajectories lead to stable walking gaits. TO problems for the legged robot are challenging to solve due to their nonlinear dynamics, high degrees of freedom and the system’s hybrid nature brought on by ground impact. Previous works such as Posa et al. (2014); Hereid et al. (2016); Medeiros et al. (2020) have proposed methods to transcribe this as a Nonlinear Programming (NLP) problem through direct collocation methods where polynomial splines are used to approximate the continuous dynamics and thus reducing computational complexity without needing to account for the actual dynamics. Other works such as Hereid et al. (2015) have instead proposed utilizing multiple shooting methods to break the original problem down into smaller steps without approximations.

In both cases, however, the system’s dynamics need to be considered along with contact dynamics to generate the trajectories. These works Pardo et al. (2017); Xin et al. (2020), to dodge the need to define ground contact dynamics explicitly, have employed null space projection methods. Meaning, zero acceleration constraints are enforced on feet ends in Hereid et al. (2016). The issue remains, and these methods are extremely computationally expensive and cannot be implemented in real-time, taking a few minutes to solve the TO problem. For works that use reduced-order models such as the centroidal dynamics or utilize Zero Moment Point (ZMP) based methods as in Pratt et al. (2006); Winkler et al. (2015), experimental results of online optimization are available. These results are restricted to pseudo-static gaits rather than dynamic ones.

In the second category (i.e., optimization-based control schemes), the goal is to compute constraint-aware feedback stabilizing control loops. This goal is most commonly achieved through a predictive framework, usually by creating a linearized model over a finite-time horizon. For instance, in Bledt et al. (2018); Fahmi et al. (2019); Angelini et al. (2019), in a hierarchical framework, reduced-order models around the COM are used to generate reference acceleration for the low-level tracking controller. The downfall of these options is the need to linearize and simplify the robot’s underlying dynamics to make it feasible in real-time, and as a result, not all constraints can be considered. In Hutter et al. (2014); Xin et al. (2020); Galloway et al. (2015); Nguyen and Sreenath (2016), the desired control inputs are computed taking into account the full-dynamics, and then optimization is carried out on a simplified least square or QP problem for tractability.

A different approach, which is the focus of the third category in our list, is to remove optimization from the control strategy and instead modify the reference trajectories to obey desired constraints. This idea was popularized through Reference Governors (RG) Kapasouris et al. (1988), where an efficient online optimization method is employed. The core to this idea is that the reference trajectory that the controller must follow can be manipulated while keeping it close to the original trajectory in the event that boundaries created by constraints are to be violated. Since its inception, this idea has spawned many variations Kolmanovsky et al. (2012); Garone et al. (2017) including an optimization-free approach known as Explicit Reference Governor (ERG) Nicotra and Garone (2015). Besides the possibility of utilizing an optimization-free form, the other major advantage with RG is that it acts as an add-on scheme to an existing controller without the need for any modification on the control scheme. However, these methods have only proven useful in relatively simple nonlinear systems as their applications to high-dimensional nonlinear problems are nearly impossible.

Other emerging paradigms such as Approximate Dynamic Programming (ADP), reinforcement learning, decoupled approaches to design control for nonlinear stochastic systems, etc., can potentially remedy the challenges in the future. For instance, employed a learning method to understand the stability and robustness of stability against the external disturbances of a passive biped walker. They used a multi-objective, multi-modal particle swarm optimization algorithm to find stable initial conditions for their biped walker model. These approaches are currently far from providing any practical solutions to the problem of performance satisfaction in dynamic bipedal walkers. They are shown to be only effective on simpler practical robots, mainly those that can only demonstrate quasi-static gaits.

## 2 Contributions and the Summary of the Proposed Framework

While gait design for complex, multi-DOF legged robots based on full or reduced-order models has been addressed extensively, optimization-free gait re-design in a reactive fashion and within gait strides are often considered only for quasi-static walkers. In dynamic walkers, such reactive gait adjustments take very complex forms involving optimization problems to ensure gait feasibility conditions.

One primary reason that within-strides gait adjustment in dynamic walkers involves optimization is that these bipeds possess small support polygons that leave small to no stability margins and make gait adjustment very intricate. With this observation and to combat the challenges associated with dynamic walking, our primary objective is to apply novel thrust-vectoring-based control actions in bipedal legged locomotion. Despite their potential merits, as we discussed above, the application of these thrusters has remained unexplored and existing studies, are limited to only flight control and not thruster-assisted legged locomotion.

Specifically speaking, the contribution of this work achieved due to thrusters’ presence in our robot is four-pronged.

• The first contribution of this work is that it employs thruster actions during the gait cycles for several vital reasons. These objectives include securing frontal dynamics stability, avoiding sagittal plane fallovers, and making aggressive jumps over obstacles. Previously, to achieve these goals, bipedal-legged robots entirely relied on indirect methods such as posture control (at stance phase) and rules of conservation of angular momentum (at flight phase during running) to avoid these fallover scenarios.

• This work introduces a decoupled view towards satisfying two antagonistic properties of efficiency and robustness. Meaning it shows that one can consider part of the gait cycle for performance satisfaction and the rest of the step period for gait robustification. The application of this decoupled view towards performance satisfaction is impossible in conventional bipeds. As a result, the decoupled idea remained unexplored. A major consequence of applying the decoupled view is that it allows the application of simpler computational algorithms for joint trajectory planning (or re-planning), control, and GRF constraint satisfaction, which constitute the next contribution of this work.

• Introducing an efficient controller, which secures the feasibility of the gaits and constraints, is the third contribution of this work. This work presents a systematic approach based on reference governors, which are relied on provable Lyapunov stability properties, to minimize computational overhead. This approach permits performing motion planning in the state-space of the Zero Dynamics (ZD) of our bipedal robot, yielding a significantly lower cost of computation than widely used optimization-based methods applied on full- or reduced-order dynamical models. To comprehensively describe this method, we will extend our previous works Dangol et al. (2020); Dangol and Ramezani (2020) by conjointly employing state-space motion planning tools and celebrated bipedal robot control design frameworks such as the HZD method.

We take a simulation-based approach to validate our method. We will use the simulators developed for our Harpy platform, capitalize on the thrusters’ action to resolve and re-design Harpy’s gait parameters during its Single Support (SS) phase. The SS phase is the most extended phase in Harpy’s gait cycle. Then, we will assume well-tuned supervisory controllers within our reference governor model and focus on fine-tuning Harpy’s joints’ desired trajectories to satisfy gait feasibility constraints, including saturating controls and retaining feasible contact forces.

• Finally, this work shows that achieving hybrid invariance in the face of significant external perturbations, in a finite-time fashion, using thrust vectoring is possible. Here, hybrid invariance is violated for an important reason. Since we devise intermediary filters based on RGs, there can be a mismatch between the robot’s states at the gaits’ boundaries. Put differently, Harpy’s gait modifications and impact events can lead to severe deviations from its nominal periodic orbits. In this work, we demonstrate that owing to Harpy’s thrusters hybrid invariance will be achieved by either applying the thrusters throughout the whole gait cycle or employing them intermittently.

This work is organized as follows. In Section 3, we will briefly introduce our platform, Harpy. We will avoid reporting hardware-related matters in this work as they are beyond the scope of this publication. Harpy’s hardware will be reported in subsequent publications. In Section 4, Harpy’s planar dynamics, which consist of continuous and impact models, will be derived. The robot’s SS phase will be modeled as a Reduced-Order Model (ROM) following standard conventions and assumptions found in textbooks Westervelt and Grizzle (2007). Then, a two-point DS model followed by an impact map for a non-instantaneous phase will be obtained. Hybrid invariance will be achieved using a Nonlinear Model Predictive Control (NMPC) scheme during the DS phase. Gait design based on HZD method and satisfying performance constraints based on RG will be explained in detail. Last, our simulation results will be discussed in Section 6 as we will report the performance of the proposed approach.

## 3 Quick Overview of Harpy Platform

This section will discuss Harpy’s physical properties and elaborate on our philosophy regarding the design. Also, we will briefly discuss the weight considerations and mass allocation and their impact on modeling.

### 3.1 Physical Properties

The system, shown in Figure 1, weighs roughly 4 lbs. It is 2.4 ft tall. The legs are composed of a parallelogram mechanism. With a highly integrated and unified actuation and sensing applied in the design of Harpy, which will be explained below, minimum use of metal components, including fasteners, housings, support structures, etc., is achieved. Parts are designed for modularity, energy efficiency, and low weight. It is worth noting that high energy-to-weight ratios are critical characteristics of birds which are our design role models. Birds are capable of showcasing legged and aerial locomotion to perfection.

### 3.2 Upper Body Design

The torso, shown in Figure 2, hosts two brush-less DC actuators with harmonic drives to increase the output torque. These actuators are used to move the legs in the frontal plane of motion. Two co-axial thrusters are attached to a separate actuator and can move relative to the torso in the sagittal plane of locomotion. We will use the thrusters for the following purposes: 1) to directly regulate the ground contact forces, 2) to recover the system when incidental tip-overs occur, 3) to generate the lift force required for jumping over obstacles, 4) to stabilize the frontal dynamics and 5) to secure hybrid invariance. Additionally, the torso encapsulates power electronics, including the amplifiers for actuators and thrusters, sensors, communication nodes, and a computer for online processing and closed-loop feedback.

**FIGURE 2**. **(A)** Illustrates the current status of the ongoing Harpy hardware developments at Northeastern University. **(B,C)** Show successful mono-pedal jumping tests with and without thrusters. A gantry supports the half-finished prototype in these tests.

The robot is underactuated when the thrusters are off and can be overactuated when active. Indeed, when Harpy is unconstrained, i.e., it is not in contact with the environment, it possesses 14 DOF. Six DOF are associated with the translation and orientation of the robot. There are two DOF at the locations where the hips are attached to the torso.

### 3.3 Payload Reduction and Added Mass Challenge

The mechanical design and actuation approach applied in Harpy is different from existing platforms. Instead of the generous incorporation of metal components, it has been tried to rely on composite material, fiber-inlay additive manufacturing, and embedding methods. The design approach based on embedding all of the components within the composite structure of the robot has been applied to achieve an ultra-light robot. Gaining a minimized total weight is essential in Harpy’s design. Harpy can be looked upon as a self-manipulating system therefore, large payloads can lead to the need for stronger actuators or thrusters. These actuators are often heavier, and consequently, extra mass can be introduced in the design.

### 3.4 Spring-Loaded Inverted Pendulum Model Design Considerations

The actuators and Harpy’s structure are designed to deliver a large range of motion on all DOF. Low-inertia limbs are considered to capture Spring-Loaded Inverted Pendulum (SLIP) model characteristics. Each leg, shown in Figure 1, possesses eight DOF and after considering the mechanical couplings three DOF remain. Back-driveable, harmonic drive actuators replace extra torque or force sensors or series compliance which can help reduce the overall mass of Harpy’s legs.

Low limb inertia and low reflected actuator inertia make the robot capable of extremely fast leg-swings. Besides, the hip and leg actuators are located so that their axes of rotation intersect at the hip joint. This design consideration helps reduce the moment of inertia in the frontal plane. The inherent compliance in the legs, particularly in the 4-bar linkage, can potentially reduce mechanical bandwidth in the legs. This property can affect foot placement performance. In a recently completed platform that uses a similar leg mechanism called Husky, fast and precise foot placement was achieved, supporting the feasibility of Harpy’s leg design.

## 4 Harpy Modeling

This section will briefly outline Harpy configuration space, underactuated and actuated coordinates, full-model, and ROMs. We will break down the gait cycle into two distinct phases: 1) an SS phase during which only one leg, called the stance leg, is in constant contact with the ground and 2) a DS phase during which both legs are fixated to the ground for a very short period. The transitions, including SS to DS and DS to SS, will be marked by the swing leg touch-down and lift-off, respectively as shown in Figure 3.

**FIGURE 3**. Illustration of the SS and DS phases. Also, shows the thrusters utilized to regulate the GRF.

Assuming point-foot walking, we will design predefined periodic gaits for the system following the well-established design framework of HZD as discussed in Westervelt et al. (2003). We will avoid discussing gait parameterization, design, and ZD derivations as the procedure has been exhaustively reported in the literature.

### 4.1 Configuration Space and Control Inputs

We assume a body coordinate frame attached to the torso COM with the x- and *z*-axis pointing forwards and upwards, respectively. During the SS phase, Harpy’s full-model possesses 14 DOF. The joint angles include roll (*q*_{r}), pitch (*q*_{p}), yaw (*q*_{y}), thigh (*q*_{iThigh}), knee (*q*_{iKnee}), hip (*q*_{iHip}), thruster (*q*_{iThrust}) angles for

For our planar models, the leg labels are switched at the moment of impact. The hip *q*_{iHip}, thigh *q*_{iThigh} and *q*_{iKnee} angles represent leg abduction-adduction, swing and flexion-extension motions, respectively. The thruster angle, *q*_{iThrust}, represents the rotation around the torso *y*-axis. By considering *p*_{1}, we permit the leg-end to bounce and slide on the contact surface.

Leg joints and thruster joints are actuated for a total of eight DOF, while attitude angles and stance leg-end positions constitute the underactuated coordinates, i.e., six DOF. The control vector *λ*_{i} and stance leg GRF *F*_{1}.

Each co-axial thruster introduces two control actions, including a net thrust force and yaw moment in the thruster coordinate frame. These thruster wrench vectors are mapped (based on Harpy’s geometry) to a net wrench vector in the body coordinate frame

where *R*_{i} is a rotation about *q*_{iThrust}, and *S* (*d*_{i}) is a skew-symmetric matrix that represents the distance component of the cross product

**FIGURE 4**. Illustrates the stick diagram of the 3D and 2D models. The full 3D model is shown on the **(A)**. The equivalent five-, three-link and *Variable-Length, Inverted Pendulum (VLIP)* models are shown in **(B–D)**, respectively.

The planar five-link configuration vector embodies seven variables. Torso angle (*q*_{p}) is measured relative to the vertical line to the ground surface. The thigh and knee angles are defined similarly to the 3D model while the hip joint angles along with the yaw and roll angles about the world coordinate frame are ignored. As a result, only the thruster forces *F*_{xThrust} and *F*_{zThrust}, expressed in the torso coordinate frame, are considered part of the control input. Stance leg-end positions are parameterized with the normal (*p*_{n}) and tangent (*p*_{t}) terms to complete a 7-DOF configuration variable vector. The configuration variables for planar three-link and VLIP models are also constructed similarly as shown in Figure 4.

### 4.2 Restricting Frontal and Transversal Dynamics Using Thrusters

The thrusters are located at each side of the body and are aligned with the body COM about the sagittal axis. This design allows the use of the thruster net force and torque around the COM to stabilize the robot’s roll and yaw angles. To restrict the motions of the 3D model to the sagittal plane of locomotion, we employ thruster actions and a simple closed-loop controller. The controller stabilizes the frontal and transversal dynamics. As a result, 2D gaits designed for the sagittal plan of locomotion can be easily applied to the 3D model. For instance, frontal dynamics stabilization is achieved by employing opposite thruster forces at the left and right thrusters.

### 4.3 Single Support Model

At single support, the system retains only a single point of contact with the ground. This contact point is not fixed. A simple static friction model describes the contact force. Regarding the mass allocation in the system, the torso is modeled as a single point mass with no moment of inertia. Since the hip and thigh actuators are located near each other, a single point mass is considered for both of them. Similarly, the knee actuator is modeled as a point mass with no moment of inertia. All of the connecting rods are assumed to be mass-less.

The mass allocation in the planar models is slightly different from the 3D model. For instance, in our five-link model, a single point mass is fixated to the hip joint and is equivalent to the overall mass at the hip in the 3D model. In the three-link model, each leg mass is equivalent to the knee actuator mass, which is fixated to the lower leg in the five-link model. The hip and thigh from the 3D model are summed up and fixated at the hip joint in the three-link model. Last, in the VLIP model, the system’s total mass is considered a point mass.

The governing equations of motion are derived using the methods of Lagrange by taking the kinetic

where *D*_{s} (*q*_{s}) is the symmetric inertial matrix and is only dependent on *q*_{s}, the *G*_{s} (*q*_{s}) contains gravity terms, and the control matrix *B*_{s}(*q*_{s}) maps the inputs to the generalized coordinate accelerations

where the state vector is denoted by

### 4.4 Non-Instantaneous Double Support Model

The end of the SS phase is marked by an impulsive impact event when the swing leg-end *p*_{1} and *p*_{2}) will be fixed to the ground during the DS phase until the new swing leg lifts off. We will assume fixed DS time intervals for all gait cycles. This switching from the SS to DS phase is captured by defining the switching conditions *p*_{2n} = 0 and *p*_{2t} > 0. Switching happens when the swing leg-end height, *p*_{2n}, is zero and its tangential (horizontal) position *p*_{2t} > 0, i.e., the swing leg lands in front of the stance leg.

The map to describe leg-end impact is modeled as in the literature Hurmuzlu and Marghitu (1994) which solves for the post-impact states

In order to formulate this impact map Δ(*x*_{s}), the model given by Eq. 3 is considered. The vector *q*_{s} is augmented to include the new stance leg-end position *p*_{2}, *p*_{1} at the moment of impact. The equations of motion are re-formulated to include the impulsive GRF,

where

The inertial matrix *D*_{d} (*q*_{d}) is a square, symmetric and positive definite, and

## 5 Control

Here, we outline an overview of our approach to satisfy state, control and GRF constraints during gait cycles by *deforming* the ZD Manifolds whose restriction dynamics governes the gaits. The notion of directional derivatives *y* = *h*_{s} (*x*_{s}) = *h*_{s} (*q*_{s}) will be adopted in a similar way as mainstream publications in this field.

The holonomic constraints (*y*) are widely known as Virtual Constraints (VCs) since they are enforced by closed-loop feedback Westervelt et al. (2018). Based on these constraints, celebrated control invariant sets of the form

Since there is only one degree of underactuation in all of our planar models, therefore, *q*_{1} = *r*_{1} (*q*_{n}), *…*, *q*_{n−1} = *r*_{n−1} (*q*_{n}) where *q*_{n} in all of our planar models can be the last entry of *q*_{s}. As a result, the output function almost similarly takes the following form for all of the models

Here *q*_{j} denotes the actuated joint variables, and *H*_{2} extracts *q*_{j} from

However, these methods happen in a discrete-time fashion, and the resulting closed-loop system possesses a very small basin of attraction. Part of the reason the continuous deformation of Γ has never been considered before is that standard bipedal robots, with their small support polygon, cannot achieve stable periodic orbits when *r* (*q*_{n}(*t*)) = *r* (*q*_{n} (*t* + *T*)), where *T* is the gait period.

Second, we want Γ to remain stabilized at all times; otherwise, enforcing VCs will be impossible. Last, we want gait feasibility conditions, including the equality

where the first line governs the restriction dynamics (*a*_{i} (*q*_{n}) are nonlinear terms) and the second line is the condition for the stabilizability of Γ. Widely considered gait feasibility constraints such as |*x*_{s}| < *x*_{max}, |*u*| < *u*_{max}, 0 < *F*_{N} and *F*_{T}, *F*_{N} are the tangential and normal components of GRF and *μ* is the coefficient of friction.

In our approach, the positive invariance property plays a key role and is closely dependent on how *u* such that when *q*_{s} (0) is on *q*_{s}(*t*) remain on the manifold Γ for all *t* > 0.

If this property is guaranteed, the constraint satisfaction problem can be transformed into a motion planning problem in the state space of the internal dynamics that can be conveniently tackled using simple path-planning tools. This particularly becomes important and handy when fast and reactive gait planning is needed in dynamic walkers. We will further elaborate on this concept with a simple example.

### 5.1 Motivation Behind Valid Deformations of ${q}_{s}={h}_{s}^{-1}(0)$

To motivate our idea, we will consider *f*|_{Γ}, which is prescribed by *u*, takes the following form *x*_{2} ≠ 0.

Now, assume that we have the means to deform Γ and that the goal is to satisfy the constraint given by *x*_{1} ≤ 1/2. The deformation shown in red in Figure 4D cannot result in a stabilizable and consequently positive invariant Γ set. Notice that the control vector field in this example is *g* = [0,1]^{⊤} which is shown in Figure 4D. Please note that possible options where *g*(*x*) are not orthogonal – or *g*(*x*) is transversal to the deformed manifold (Γ) – are not unique.

With this simple example describing the challenge, our first step would be to identify the valid deformations of

### 5.2 Choices and Options on How to Manipulate ${h}_{s}^{-1}(0)$

Consider the state-space representation of the system dynamics given by Eq. 3. Here, we will explore our options and choices in order to manipulate Eq. 6. While many options can be considered, we will focus on two general frameworks, as explained below.

With this in mind that the control vector field is given by *ω*_{i}*q*_{n} and 2) shifting *q*_{n} + *ω*_{i}. We are interested in knowing how the stabilizability of *Γ*_{ω} is affected if these two options are applied to manipulate the state-dependent equilibrium point of the actuated coordinates. In equilibrium shifting scenario (case ii), notice that the position and velocity priming can be trivially equivalent as is evidently seen below

where *ω*(*t*) and *ν*(*t*) are used for manipulating *h*_{s} (*x*_{s}) and *x*_{s} = *f*_{s} (*x*_{s}) + *g*_{s} (*x*_{s})*u*_{ω}, which is not desired here. To see this, it is enough to consider the numerical parameters *ω*, which are used to parameterize *r* (*q*_{n}, *ω*), as auxiliary control input in discrete Poincare maps. These auxiliary control inputs can only provide discontinuous means of priming *Γ*_{ω} at its boundaries.

##### 5.2.1 Stabilizability Condition and First Scenario

We will show that by choosing (i), the continuous manipulation of the primer parameters can violate the transversality condition (*g*_{s} (*x*_{s})⊥̸*∂h* (*x*_{s})/*∂x*_{s}). To do this, consider *y* = *q*_{j} − *r* (*q*_{n}, *ω*). Then, consider the fact that the inertia matrix *D*_{s} (*q*_{n}) only depends on angles which are evaluated here on *Γ*_{ω}, i.e., *q*_{n}, it is possible to show that

where *B** = [0^{1×(n−1)}, 1] and Eq. 9 can vanish on a point on *Γ*_{ω} which implies *g*_{s} (*x*_{s}) and *Γ*_{ω}, these violations are not acceptable in our proposed framework and we avoid that. Without loss of generality, consider *M*th order polynomials (e.g., Bezier forms with degree *M* as in ref. [40]) to define

where *ω*(*t*) as diagonal entries,

where *D*_{i} (*q*_{n}) are the block matrices of *D*_{s} (*q*_{n}). This algebraic relationship can be easily solved for a fixed *q*_{n} and *ω*(*t*) where *ω*(*t*) the transversality condition can be violated on *Γ*_{ω}.

##### 5.2.2 Stabilizability Condition and Second Scenario

We will now consider the other option, i.e., scaling. It is straightforward to show that *Γ*_{ω}. It is also noticeable that this choice of manipulating *y* has no effects on *y* = *q*_{j} − *r* (*q*_{n}) is a valid VC. Notice that the role of the primer *ω*(*t*) in this form is comparable to the role of a disturbance term in the system given below

where

The roles of this disturbance term will be beneficial for us, though, as it will adjust the equilibria of *h*_{s} (*x*_{s}) and *ω*(*t*) such that the finite-time convergence of the solutions to *Γ*_{ω} in the closed-loop system is unaffected. In addition, while the evolution of the trajectories in the constraint-admissible subspace (this will be explained shortly) of the state space is secured, we want *q*_{j} closely track predefined *r* (*q*_{n}, *ω*(*t*)) when possible. While stronger stability properties (e.g., global exponential stability) are desirable, our major concern will be the finite-time enforcement of the changed virtual constraints.

As far as the design of *u* is concerned, any nonlinear controller (or linear controllers if the nonlinear terms are bounded and the bounds are known) can satisfy the VCs. We will limit ourselves to the following modest feedback law *Γ*_{ω} in order to satisfy our constraints. Hence, we will assume stabilizing supervisory controllers that guarantee the enforcement of the virtual constraints. However, their disturbance rejection properties have to be carefully considered.

The control law given above can generate Global Asymptotic Stability (GAS) at the equilibrium point of the system given by Eq. 12 when the primer variable vector is time-invariant, i.e., *ω*(*t*) = *ω*. Subsequently, *Γ*_{ω} takes the following form

where the equilibrium points for *h*_{s} (*x*_{s}) and *ω* are obtained by solving

Of course, realizing GAS property under *ω*(*t*), i.e., when the primer variable is time-varying, as we shall see later, will not be achieved trivially and requires a regulated rate of change in the primer variable vector. We will discuss this in the proof of GAS property of the closed-loop system later in Section 5.4.

A traditional approach to deal with a time-varying disturbance is to synthesize a family of linear controllers at each equilibrium point of the system given by Eq. 12 and then synthesize the stabilizing controller based on interpolating between the linear controllers by gain scheduling. However, nonlinear update laws will produce better results. Using them, it would be possible to directly incorporate the constraints, particularly when gait re-planning has to happen quickly.

Next, with the role of the primer variable *ω*(*t*) in our output function *y* = *h*_{s} (*q*_{s}, *ω*(*t*)) set, we will show that states, joint torques and GRF can be expressed in terms of the primer variable.

### 5.3 Constraints Derivation

Consider the configuration variable vector *m* is the number of these joints in the planar models and *p*_{1} is stance leg contact point as defined previously. The control matrix in the Euler-Lagrange equations can take the following form

where *p*_{T,COM} is the physical location of the thruster action *F*^{B}. Notice that based on how the thruster actions *F*^{B} are incorporated in the Euler-Lagrange equations (Eq. 2), the coordinate *q*_{n}, which was primarily assumed to be underactuated, will be actuated. While this is not a problem, to keep our desired normal form – i.e., *q*_{n} remains underactuated – and to allow the direct regulation of GRF, *F*^{B} is found such that *F*^{B} can be translated to the base of the VLIP model based on the law of translating forces and moments.

In other words, if it is assumed that *q*_{n} angle as it was planned and the possibility of directly regulating GRF. Note that aligning *λ*_{net} with the virtual leg in the model takes place by employing the thruster joint actuators *F*^{B} which are responsible to move the thrusters with respect to the body.

One interesting interpretations of the orthogonality condition (i.e., *q*_{n}) are allowed which can be problematic given that in our design framework all joint motions are parameterized based on *q*_{n}. In other words, very small Δ*q*_{n} means no locomotion.

Based on the adjustment we made in the control matrix *B*_{s} and that *α* ≠ 0. As a result, the following equation can be obtained at every point on *Γ*_{ω} and it governs the dynamics of the

In this equation *β*_{1} is given by

where *Q*_{i} (*q*_{n}) are the Christoffel Symbols. A similar algebraic relationship for the constraints *y* and

Next, we will steer *y* and *ω*(*t*) in Eq. 12 in order to make sure the solutions of Eq. 16 stay within the constraint-admissible space. To do this, consider *y*-*q*_{n}-*y*-

Since we assumed a pre-stabilized system – in fact all of the above derivations only make sense if *c*_{l} and *c*_{u} are constraint lower and upper bounds) based on the steady-state solutions, i.e., *y*_{ω} and

This assumption may cause the intermittent violation of the constraints, which is expected. However, because finite-time convergence to the constraint-admissible sub-spaces at the neighborhood of *y*_{ω} and

Other than simplifying the nonlinear constraint satisfaction problem given in Eq. 7, considering

which is the locus of all of the steady-state solutions of the system Eq. 12. We will show that it is possible to create positive invariant sets around any point *Y*_{ω}. Then, one can slide *Y*_{ω} using a suitable update policy to ensure the following objectives are achieved.

First objective is the realization of positive invariance property in a neighborhood around the steady-state solutions on *Y*_{ω}. Second objective is to satisfy gait feasibility constraints as outlined before. Last, we want to maintain a minimal distance from the origin, i.e.,

### 5.4 Primer Variable Vector, *ω*(*t*), Update Policy and Achieving Global Asymptotic Stability Property Over *Y*_{ω}

Consider the support hyperplane given by

at the point *C*_{i}(*z*, *ω*(*t*)) – i.e., the *i*th entry in the constraint vector *C* (*z*, *ω*(*t*)) given by Eq. 7 – and *Y*_{ω} intersect (see Figures 5, 6). Since there are multiple constraints, we consider the constraint with the shortest distance from *z*_{ω} on *Y*_{ω}. While there are many ways to simplify the geometric representation of the nonlinear constraints *C*_{i}(*z*, *ω*(*t*)) at the neighborhood of *z*_{ω} on *Y*_{ω}, the approach adopted here is very efficient with a minimum computation overhead.

**FIGURE 5**. Geometric interpretation of the level set {*ζ*_{w}|*V* ≤ *V*_{max}}. For a constant reference *r*, the update law *w*, places the new value *w** at the edge of the level set *V* = *V*_{max} while trying to be as close to the actual reference as possible without violating the constrains *C*_{i}(*z*, *ω*(*t*)), which are depicted in red.

**FIGURE 6**. **(A–C)** Illustration of constraint satisfaction based on manifold deformation; Also, it shows the concept with various norm definitions. **(D)** Illustration of the relationship between the stabilizability of the manifold Γ (unit circle is used as an example) and the transversality condition, i.e., *G*(*x*) is the algebraic equation that defines Γ and *g*(*x*) is the control vector field used to stabilize the trajectories on Γ.

Our approach is motivated by classical methods widely used in constraint optimization problems based on support hyperplanes. These methods are usually very conservative, but in our case their use is justified as the system is highly nonlinear, unpredictable. That said, the location and distance of these hyperplanes from *z*_{ω} has to be carefully defined.

We will use Lyapunov functions with quadratic forms to define level sets around *z*_{ω}. The largest invariant set confined in the constraint-admissible space will define the maximum distance between *z*_{ω} and *c*_{i} (*z*, *ω*(*t*)). These invariant sets, geometrically, can take various forms using non-quadratic energy functions. It is important to note that the geometry of the invariant set around *z*_{ω} directly affects the performance of this approach, i.e., the convergence rate to

Motivated by deepest gradient descent method which is a key mechanism of manipulating pre-compensated systems, we consider two search parameters including a step direction, *ρ* = − *ω*(*t*)/‖*ω*(*t*)‖_{2}, ‖*ω*(*t*)‖_{2} ≠ 0, and step length,

The step direction *ρ* is trivially obtained based on the observation that the primer variable *ω*(*t*) and consequently the steady-state solutions *z*_{ω} should remain possibly very close to the origin. The step length *δ* is obtained based on the difference between the values of a candidate Lyapunov functions and the largest invariant sets around *z*_{ω}. Based on these two search parameters, the following update law for the primer variable vector *ω*(*t*) is considered:

where

Consider *ζ* = *P*^{−0.5}*z* and *ζ*_{ω} = *P*^{−0.5}*z*_{ω} where *P* is a positive definite matrix. The maximum level set value *V*_{max,i} (*ω*(*t*)) for the quadratic Lyapunov function *V*_{max,i} (*ω*(*t*)) can geometrically be interpreted with unit circles around *z*_{ω} in *z*-space. The size of the largest level set *V*_{max,i} (*ω*(*t*)) that is not violating the closest constraint *z*_{ω} can be better imagined when 2-norm ‖.‖_{2} is used, shown in Figure 6. The unit vector *z*_{ω} and the constraint hyperplane *c*_{i} = 0 is given by

Now consider the quadratic norm ‖*ζ*‖_{P}. Notice that using the change of variable *ζ* = *P*^{−0.5}*z* in the constraint equation *c*_{i} (*z*, *ω*(*t*)) = 0, 2-norm and the quadratic norm become interchangeable. In this way, the maximum level set value for any candidate, quadratic Lyapunov function can be obtained. Therefore, the step length in Eq. 21 is given by *δ* = *V*_{max,i} (*ω*(*t*)) − *V* (*ζ*, *ω*(*t*)). After the change of variable, the unit vector *ζ*_{ω} in a unit ball of ‖.‖_{P} is given by

where *γ*‖_{p} which is ‖*γ*‖_{*} = ‖*P*^{−0.5}*γ*‖_{2} and this algebraic relationship

Now, we make sure under the update law given by Eq. 21 the constraint-admissible set defined at any point on *Y*_{ω} is a positive invariant set. Taking any Lyapunov function of the form given above and using Barbashin-Krasovskii-LaSalle principle, for a constant primer variable *ω* the GAS property of the closed-loop system given by Eq. 12 is readily achievable.

Take the time-varying primer variable *ω*(*t*) governed by the update law given above, then

where *Q* is a positive definite matrix. For *κ* ≪ 1 such that

As a result, the resulting inequality *κ* ≪ 1. This result shows that if the primer dynamics is very fast, i.e., 1/*κ* in Eq. 21 is a very large value, then the closed-loop system given by Eq. 23 is GAS.

Next, we briefly show how the thruster can be used to achieve hybrid invariance.

### 5.5 Thruster-Assisted Impact Invariance

One critical aspect of HZD-based methods is selecting *y* = *h*_{S} (*x*_{s}) that leads to a hybrid zero dynamics. In these methods, the assumption of underactuation at the contact points does not leave any options for control better than the deadbeat hybrid extension of these systems to achieve impact invariance. This section briefly discusses another use of the thruster actions to secure hybrid invariance, complementary to HZD-based methods.

At the transition from SS to DS phase, the two-point impact, as discussed in Section 4.4, renders all of the body joints except the torso joint fixated to their pre-impact values. Hence, large deviations in joint velocities from the reference trajectories and subsequently from the ZD manifold (*Γ*_{ω}) will be resulted.

Since the joint actuators are not able to make corrections needed to steer the states back to *Γ*_{ω} the thrusters can be leveraged in the DS phase to achieve hybrid invariance. Here, impact invariance such that the initial DS phase state *x*_{d,0} returns to the initial SS phase state *x*_{s,0} is sought where

As opposed to the SS phase, the constraints in the DS phase take a more complex form. Additionally, we have to ensure that the final state of the DS phase (*x*_{d,f}) matches the initial states at the SS phase (*x*_{s,0}).

We apply a Nonlinear Model Predictive Control (NMPC) scheme to steer the post-DS states back to the ZD manifold. This scheme is known for being costly. However, the duration of the DS phase is significantly shorter than the SS phase. Note that a reference trajectory for each DS state *r*_{d} [*k*] is generated at every *k*th sample over the duration of the DS phase. The reference trajectory can be a simple linear trajectory between the post-impact states *x*_{d,0} and the initial SS phase state *x*_{s,0}.

The continuous DS phase model as described in Eq. 4, along with the kinematic constraint

where *η* = [*u*^{⊤}, *λ*_{net}] is an augmented input vector that consists of the control action from the SS phase along with the net thrust force along the torso, which can be expressed in the state-space form given by *ϕ*(*x*_{d}, *η*) to be minimized is formulated resulting in the following optimization problem

In Eq. 26, *W*_{r} is the weighting term for the cost associated with tracking the reference trajectory *r*_{d} [*k*], *W*_{η} similarly contains the weight that penalizes rapid changes *δ* in the input *η*, and *N* denotes the total simulation time steps. The constraints of the optimization problem are denoted by *C*_{eq} (*x*_{d}, *η*), *C*_{ineq} (*x*_{d}, *η*) which represent the equality and inequality constraints, respectively. The equality terms include the initial DS state, which is obtained from the post-impact SS state

In Eq. 27, the matrix

Here *ϵ* is a relaxation term applied to the final state. With all constraints satisfied, the NMPC guides the DS states towards the initial condition of SS phase, resulting in impact invariance.

## 6 Simulation Results

In this section, we will report the result of our simulation works. An overview of the control architecture explained previously is given in Figure 7. We consider walking on flat ground when Harpy makes point contacts with the surface. Using the full-dynamics model explained above, we will apply the thruster actions and proposed framework to demonstrate the desired capabilities, including, frontal stabilization, trajectory manipulation to stay in constraint-admissible regions of the state space, hybrid invariance and high jumps over obstacles.

Since one of the motivations behind our approach is to reduce computation overhead, we will briefly demonstrate the results from integrating our ROMs (i.e., the VLIP, 3-link, etc., models explained previously) in this process. We will demonstrate that this integration can further reduce the costs.

We note that the gaits are designed based on the HZD approach. In the SS phases, the desired trajectories *h*_{d} (*x*_{s}) are parameterized as Bezier polynomials with the coefficients tuned offline. All model parameters used in the simulation closely match Harpy’s properties and are listed in Table 1.

The full-dynamics model can be constrained in its frontal plan of locomotion using the thrusters. We note that frontal dynamics stabilization due to underactuation can be very challenging. We ran the simulations using the full-dynamics model of Harpy where the legs are modeled with all DOF including the hip frontal, hip sagittal, and knee sagittal joints. Figure 8 shows the stick-diagram results from several key states during walking. The thrusters not only are used to stabilize the frontal dynamics, but also to jump over obstacles, and assist the robot to land stably.

**FIGURE 8**. Illustration of thruster-assisted locomotion over a rough terrain simulated on Harpy’s full-dynamics. Blue and red arrows show the thruster actions and GRF, respectively.

The simulation results such as the robot’s states and GRF are shown in Figure 9. In general, planar gaits are unstable in 3D systems. The simulation result demonstrate that the thruster-assisted robot is capable of tracking the target body trajectory and walk stably even when the 2D gait is used. The thrusters can assist the robot to walk with a stable and nearly zero roll and heading angles throughout the simulation and can be seen in Figure 8.

**FIGURE 9**. Illustrates the state values and forces during walking and jumping feats simulated with Harpy’s 3D full-dynamics. The reference joint trajectories (*z*_{0}) are manipulated (*z*_{ω}) such that the applied constraint equations do not violate the constraint as shown at *t* = [9.5, 11.5]s.

In Figures 9, 10, the thrusters are used to satisfy the constraints including saturating the control actions and enforcing no-slip conditions. First, a desired trajectory for the COM of the full-dynamics is considered and using the priming approach explained above it is manipulated during the walking phase to satisfy GRF feasibility constraints. The primer is disabled during the flight phase. As shown in Figure 9, the GRF constraints are violated at around *t* = [9.5, 11.5]. Figure 11 depicts a close-up view of the priming performance at *t* = [9.5, 11.5].

**FIGURE 10**. Illustrates the effects of reducing the control thresholds on the control actions and their corresponding phase portraits.

**FIGURE 11**. Illustrates the closeup view of the constraint violations during 3D walking and jumping at *t* = [9.5, 11.5]s. The reference trajectories violate the constraints and the controller manipulates the applied references to prevent the constraints from being violated.

In addition, we considered constraining the control actions. To do this, the control inputs were saturated in the low-dimensional space of the ROMs and then were mapped to the full-dynamics. The priming problem to saturate control inputs in the three-link ROM is considered. Figures 12, 13 show the desired joint angles and velocities, respectively.

In Figure 10, the results after applying the primer are shown. The first two figures from left show the effects of changing the limits on the control inputs during the SS phase. The limits are adjusted from high to low values resulting in the peak control to decrease below the desired thresholds. However, a large gap can be seen between the threshold and maximum control input. Note that the geometries of the positive invariant sets around steady-state responses in the output dynamics are not optimal in our work as explained before and shown in Figure 5. The phase portraits **(***q*_{1}**-****)** corresponding to each control limits are shown in Figure 10. It can be seen that

Hybrid invariance can be achieved using the thruster actions. Figure 14 shows the feasibility conditions during the DS phase where NMPC is applied. The friction cone condition is feasible which indicates that both feet are fixated to the ground throughout the DS phase.

**FIGURE 14**. Illustrates the friction constraints during the DS phase of leg 1 and 2. In this plot, the intermediate SS phases are omitted. Dashed blue lines are the thresholds.

Note that the static friction coefficient is assumed to be *μ* = 0.3. Also note that the normal forces spike to 60 N during the DS phase. The total weight of the biped is only 0.8 kg which suggests that the unusual behavior is only possible because of the thruster actions. Figure 15 shows the control actions for the thrusters during the DS phase.

**FIGURE 15**. Illustrates the thruster actions during the DS phase. In this plot, the intermediate SS phases are omitted as the thrusters are inactive.

## 7 Concluding Remarks

In this work, we studied the roles of thrusters in addressing the prevailing challenges in bipedal robotics, potentially paving the way to see new designs with unexplored capabilities. Some of the common challenges faced by bipedal robots include constraint satisfaction, frontal dynamics stabilization, and avoiding fallovers. Combatting these issues can be pretty challenging in these systems due to underactuation and high vulnerability to external perturbations.

In this publication, we introduced a thruster-assisted bipedal robot called Harpy. Harpy platform possesses two legs and two coaxial thrusters attached to its torso. We capitalized on Harpy’s unique design to propose an optimization-free approach to satisfy Harpy’s gait feasibility conditions, including control and contact forces. The reference trajectories were manipulated based on deforming stabilizable zero-dynamics manifolds to fulfill constraints brought on by ground contacts and those prescribed by states and inputs without violating hybrid invariance.

In standard bipedal robots, unintended changes to the restriction dynamics over the zero dynamics manifold, especially those optimized to produce periodic orbits, can adversely affect gait stability. Not only we showed that the manipulation of the system trajectories is possible in Harpy full-dynamics, but also hybrid invariance can be realized by employing the thrusters, something that is often achieved in a limited fashion by event-based regulators.

In addition, we demonstrated that the thrusters can be utilized to robustify the gaits by dodging fallovers or even jumping over large obstacles.

## Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

## Author Contributions

PD produced the simulation results with major contributions from ES; AR wrote this paper, and the paper was reviewed by PD, ES, and AR together. The mathematical framework reported in this paper is developed by AR.

## Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

## References

Angelini, F., Xin, G., Wolfslag, W. J., Tiseo, C., Mistry, M., Garabini, M., et al. (2019). “Online Optimal Impedance Planning for Legged Robots,” in 2019 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS) (Macau, China: IEEE), 6028–6035. doi:10.1109/iros40897.2019.8967696

Arcos-Legarda, J., Cortes-Romero, J., Beltran-Pulido, A., and Tovar, A. (2019). Hybrid Disturbance Rejection Control of Dynamic Bipedal Robots. *Multibody Syst. Dyn.* 46 (3), 281–306. doi:10.1007/s11044-019-09667-3

Apgar, T., Clary, P., Green, K., Fern, A., and Hurst, J. W. (2018). Fast Online Trajectory Optimization for the Bipedal Robot Cassie. *Robotics: Sci. Syst.* 101, 14. doi:10.15607/rss.2018.xiv.054

Bledt, G., Powell, M. J., Katz, B., Di Carlo, J., Wensing, P. M., and Kim, S. (2018). “Mit Cheetah 3: Design and Control of a Robust, Dynamic Quadruped Robot,” in 2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS) (Madrid, Spain: IEEE), 2245–2252. doi:10.1109/iros.2018.8593885

Dai, H., and Tedrake, R. (2016). “Planning Robust Walking Motion on Uneven Terrain via Convex Optimization,” in IEEE-RAS International Conference on Humanoid Robots (Humanoids) (Cancun, Mexico: IEEE), 579–586. doi:10.1109/humanoids.2016.7803333

Dangol, P., Ramezani, A., and Jalili, N. (2020). “Performance Satisfaction in Midget, a Thruster-Assisted Bipedal Robot,” in 2020 American Control Conference (ACC) (Denver, CO, USA: IEEE), 3217–3223. doi:10.23919/acc45564.2020.9147448

Dangol, P., and Ramezani, A. (2020). Towards Thruster-Assisted Bipedal Locomotion for Enhanced Efficiency and Robustness. *arXiv preprint arXiv:2005.00347*. doi:10.1016/j.ifacol.2020.12.2721

Fahmi, S., Mastalli, C., Focchi, M., and Semini, C. (2019). Passive Whole-Body Control for Quadruped Robots: Experimental Validation over Challenging Terrain. *IEEE Robot. Autom. Lett.* 4, 2553–2560. doi:10.1109/lra.2019.2908502

Feng, S., Whitman, E., Xinjilefu, X., and Atkeson, C. G. (2014). “Optimization Based Full Body Control for the Atlas Robot,” in IEEE-RAS International Conference on Humanoid Robots (Madrid, Spain: IEEE), 120–127. doi:10.1109/humanoids.2014.7041347

Galloway, K., Sreenath, K., Ames, A. D., and Grizzle, J. W. (2015). Torque Saturation in Bipedal Robotic Walking through Control Lyapunov Function-Based Quadratic Programs. *IEEE Access* 3, 323–332. doi:10.1109/access.2015.2419630

Garone, E., Di Cairano, S., and Kolmanovsky, I. (2017). Reference and Command Governors for Systems with Constraints: A Survey on Theory and Applications. *Automatica* 75, 306–328. doi:10.1016/j.automatica.2016.08.013

Goswami, A., and Kallem, V. (2004). Rate of Change of Angular Momentum and Balance Maintenance of Biped Robots. *IEEE Xplore* 4, 3785–3790. doi:10.1109/robot.2004.1308858

Gubina, F., Hemami, H., and McGhee, R. B. (1974). On the Dynamic Stability of Biped Locomotion. *IEEE Trans. Biomed. Eng.* 21, 102–108. doi:10.1109/tbme.1974.324294

Hereid, A., Cousineau, E. A., Hubicki, C. M., and Ames, A. D. (2016). “3d Dynamic Walking with Underactuated Humanoid Robots: A Direct Collocation Framework for Optimizing Hybrid Zero Dynamics,” in 2016 IEEE International Conference on Robotics and Automation (ICRA) (Stockholm, Sweden: IEEE), 1447–1454. doi:10.1109/icra.2016.7487279

Hereid, A., Hubicki, C. M., Cousineau, E. A., Hurst, J. W., and Ames, A. D. (2015). “Hybrid Zero Dynamics Based Multiple Shooting Optimization with Applications to Robotic Walking,” in 2015 IEEE International Conference on Robotics and Automation (ICRA) (Seattle, WA, USA: IEEE), 5734–5740. doi:10.1109/icra.2015.7140002

Hill, J., and Fahimi, F. (2015). Active Disturbance Rejection for Walking Bipedal Robots Using the Acceleration of the Upper Limbs. *Robotica.* 33 (2), 264–281. doi:10.1017/S0263574714000332

Hirose, M., and Ogawa, K. (2006). Honda Humanoid Robots Development. *Phil. Trans. R. Soc. A.* 365, 11–19. doi:10.1098/rsta.2006.1917

Hurmuzlu, Y., and Marghitu, D. B. (1994). Rigid Body Collisions of Planar Kinematic Chains with Multiple Contact Points. *Int. J. Robotics Res.* 13, 82–92. doi:10.1177/027836499401300106

Hutter, M., Sommer, H., Gehring, C., Hoepflinger, M., Bloesch, M., and Siegwart, R. (2014). Quadrupedal Locomotion Using Hierarchical Operational Space Control. *Int. J. Robotics Res.* 33, 1047–1062. doi:10.1177/0278364913519834

Kajita, S., and Tani, K. (1991). “Study of Dynamic Biped Locomotion on Rugged Terrain-Derivation and Application of the Linear Inverted Pendulum Mode,” in Proceedings. 1991 IEEE International Conference on Robotics and Automation (Sacramento, CA, USA: IEEE Computer Society), 1405–1406.

Kapasouris, P., Athans, M., and Stein, G. (1988). “Design of Feedback Control Systems for Stable Plants with Saturating Actuators,” in Proceedings of the 27th IEEE Conference on Decision and Control (Austin, TX, USA: IEEE). doi:10.1109/CDC.1988.194356

Kolmanovsky, I., Kalabić, U., and Gilbert, E. (2012). Developments in Constrained Control Using Reference Governors. *IFAC Proc. Volumes* 45, 282–290. doi:10.3182/20120823-5-nl-3013.00042

Kwon, W., Kim, H. K., Park, J. K., Roh, C. H., Lee, J., Park, J., et al. (2007). “Biped Humanoid Robot Mahru Iii,” in IEEE-RAS International Conference on Humanoid Robots (Pittsburgh, PA, USA: IEEE), 583–588. doi:10.1109/ichr.2007.4813930

Mandava, R. K., and Vundavilli, P. R. (2018). Near Optimal PID Controllers for the Biped Robot While Walking on Uneven Terrains. *Int. J. Autom. Comput.* 15 (6), 689–706. doi:10.1007/s11633-018-1121-3

Medeiros, V. S., Jelavic, E., Bjelonic, M., Siegwart, R., Meggiolaro, M. A., and Hutter, M. (2020). Trajectory Optimization for Wheeled-Legged Quadrupedal Robots Driving in Challenging Terrain. *IEEE Robot. Autom. Lett.* 5, 4172–4179. doi:10.1109/lra.2020.2990720

Nguyen, Q., and Sreenath, K. (2016). “Exponential Control Barrier Functions for Enforcing High Relative-Degree Safety-Critical Constraints,” in 2016 American Control Conference (ACC) (Boston, MA, USA: IEEE), 322–328. doi:10.1109/acc.2016.7524935

Nicotra, M. M., and Garone, E. (2015). “Explicit Reference Governor for Continuous Time Nonlinear Systems Subject to Convex Constraints,” in 2015 American Control Conference (ACC) (Boston, MA, USA: IEEE), 4561–4566. doi:10.1109/acc.2015.7172047

Paine, N., Mehling, J. S., Holley, J., Radford, N. A., Johnson, G., Fok, C.-L., et al. (2015). Actuator Control for the Nasa-Jsc Valkyrie Humanoid Robot: A Decoupled Dynamics Approach for Torque Control of Series Elastic Robots. *J. Field Robotics* 32, 378–396. doi:10.1002/rob.21556

Pardo, D., Neunert, M., Winkler, A. W., Grandia, R., and Buchli, J. (2017). “Hybrid Direct Collocation and Control in the Constraint-Consistent Subspace for Dynamic Legged Robot Locomotion,” in Robotics: Science and Systems. doi:10.15607/rss.2017.xiii.042

Park, H.-W., Sreenath, K., Ramezani, A., and Grizzle, J. W. (2012). “Switching Control Design for Accommodating Large Step-Down Disturbances in Bipedal Robot Walking,” in 2012 IEEE International Conference on Robotics and Automation (Saint Paul, MN, USA: IEEE), 45–50. doi:10.1109/icra.2012.6225056

Posa, M., Cantu, C., and Tedrake, R. (2014). A Direct Method for Trajectory Optimization of Rigid Bodies through Contact. *Int. J. Robotics Res.* 33, 69–81. doi:10.1177/0278364913506757

Pratt, J., Carff, J., Drakunov, S., and Goswami, A. (2006). “Capture point: A Step toward Humanoid Push Recovery,” in 2006 6th IEEE-RAS international conference on humanoid robots (Genova, Italy: IEEE), 200–207. doi:10.1109/ichr.2006.321385

Pratt, J., Chew, C.-M., Torres, A., Dilworth, P., and Pratt, G. (2001). Virtual Model Control: An Intuitive Approach for Bipedal Locomotion. *Int. J. Robotics Res.* 20, 129–143. doi:10.1177/02783640122067309

Pratt, J. E., Krupp, B., Ragusila, V., Rebula, J., Koolen, T., van Nieuwenhuizen, N., et al. (2009). “The Yobotics-Ihmc Lower Body Humanoid Robot,” in IEEE/RSJ International Conference on Intelligent Robots and Systems (St. Louis, MO, USA: IEEE), 410–411. doi:10.1109/iros.2009.5354430

Pratt, J. E., and Tedrake, R. (2006). “Velocity-based Stability Margins for Fast Bipedal Walking,” in Fast Motions in Biomechanics and Robotics (Berlin, Germany: Springer), 299–324.

Raibert, M., Blankespoor, K., Nelson, G., and Playter, R. (2008). Bigdog, the Rough-Terrain Quadruped Robot. *IFAC Proc. Volumes* 41, 10822–10825. doi:10.3182/20080706-5-kr-1001.01833

Raibert, M. H., Brown, H. B., and Chepponis, M. (1984). Experiments in Balance with a 3d One-Legged Hopping Machine. *Int. J. Robotics Res.* 3, 75–92. doi:10.1177/027836498400300207

Song, S., Tang, C., Wang, Z., and Yan, G. (2018). Active Disturbance Rejection Controller Design for Stable Walking of a Compass-like Biped. *Trans. Inst. Meas.* 40 (14), 4063–4077. doi:10.1177/0142331217741957

Sreenath, K., Park, H.-W., Poulakakis, I., and Grizzle, J. W. (2011). A Compliant Hybrid Zero Dynamics Controller for Stable, Efficient and Fast Bipedal Walking on mabel. *Int. J. Robotics Res.* 30, 1170–1193. doi:10.1177/0278364910379882

Stephens, B. (2007). “Humanoid Push Recovery,” in 2007 7th IEEE-RAS International Conference on Humanoid Robots (Pittsburgh, PA, USA: IEEE), 589–595. doi:10.1109/ichr.2007.4813931

Vukobratović, M., and Stepanenko, J. (1972). On the Stability of Anthropomorphic Systems. *Math. biosciences* 15, 1–37.

Westervelt, E., and Grizzle, J. (2007). *Feedback Control of Dynamic Bipedal Robot Locomotion. Control and Automation Series*. Boca Raton, FL, USA: CRC PressINC.

Westervelt, E. R., Grizzle, J. W., Chevallereau, C., Choi, J. H., and Morris, B. (2018). *Feedback Control of Dynamic Bipedal Robot Locomotion*. Boca Raton, FL, USA: CRC Press.

Westervelt, E. R., Grizzle, J. W., and Koditschek, D. E. (2003). Hybrid Zero Dynamics of Planar Biped Walkers. *IEEE Trans. Automat. Contr.* 48, 42–56. doi:10.1109/tac.2002.806653

Winkler, A. W., Mastalli, C., Havoutis, I., Focchi, M., Caldwell, D. G., and Semini, C. (2015). “Planning and Execution of Dynamic Whole-Body Locomotion for a Hydraulic Quadruped on Challenging Terrain,” in 2015 IEEE International Conference on Robotics and Automation (ICRA) (Seattle, WA, USA: IEEE), 5148–5154. doi:10.1109/icra.2015.7139916

Keywords: legged robots, dynamics, modeling, control, mechatronics, automation

Citation: Dangol P, Sihite E and Ramezani A (2021) Control of Thruster-Assisted, Bipedal Legged Locomotion of the Harpy Robot. *Front. Robot. AI* 8:770514. doi: 10.3389/frobt.2021.770514

Received: 04 September 2021; Accepted: 11 November 2021;

Published: 23 December 2021.

Edited by:

Giuseppe Carbone, University of Calabria, ItalyReviewed by:

Ravi Mandava, Maulana Azad National Institute of Technology, IndiaHassène Gritli, University of Tunis, Tunisia

Copyright © 2021 Dangol, Sihite and Ramezani. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Alireza Ramezani, a.ramezani@northeastern.edu