Control of Thruster-Assisted, Bipedal Legged Locomotion of the Harpy Robot

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

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

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

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 $i\in \left\{1(stanceleg),2(swingleg)\right\}$ and stance leg-end positions $(p_1={[p_{tx},p_{ty},p_n]}^{\top })$. All of these variables are stacked inside the vector $q_s={[q_y,q_p,q_r,q_{iThigh},q_{iKnee},q_{iHip},q_{iThrust},p_1^{\top }]}^{\top }$ to form the states vector $x_s={[q_s^{\top },{\dot{q}}_s^{\top }]}^{\top }$.

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₁, 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 $u={[u_j^{\top },{\lambda }_i^{\top },F_1^{\top }]}^{\top }$ embodies the joint control actions $u_j={[u_{iThigh},u_{iKnee},u_{iHip},u_{iThrust}]}^{\top }$, thruster wrench λ_i and stance leg GRF F₁.

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 ${\lambda }_{net}={[F^B,M^B]}^{\top }$. This wrench vector consists of the net thruster force and moment on the body frame. This map is given by the adjoint transformation ${\lambda }_{net}={\sum }_{i}A{d}_{iB}^{\top }{\lambda }_i$, which can be expanded as

${\lambda }_{net}=\sum _{i=1}^2\left[\begin{array}{cc}\hfill R_i\hfill & \hfill 0\hfill \\ \hfill R_{i}S\left(d_i\right)\hfill & \hfill R_i\hfill \end{array}\right]{\lambda }_i(1)$

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 ${\vec{d}}_i\times F_i\left(d_i\right.$ is the distance from the torso to the thruster). This full-model will be used to validate the control design approach. Also, it will be employed to obtain other low-order models, including five-link, three-link and VLIP models. All of these models can only move in the sagittal plane of locomotion and are shown in Figure 4.

FIGURE 4

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 $\mathcal{K}(q_s,\dot{q_s})$ and potential $\mathcal{V}(q_s)$ energies of the system which leads to $\mathcal{L}(q_s,\dot{q_s})=\mathcal{K}(q_s,\dot{q_s})-\mathcal{V}(q_s)$. The resulting equations of motion are expressed in vector form as following

$D_s\left(q_s\right){\ddot{q}}_s+C_s\left(q_s,{\dot{q}}_s\right){\dot{q}}_s+G_s\left(q_s\right)=B_s\left(q_s\right)u,(2)$

where D_s (q_s) is the symmetric inertial matrix and is only dependent on q_s, the $C_s(q_s,{\dot{q}}_s){\dot{q}}_s$ matrix contains the Coriolis terms, G_s (q_s) contains gravity terms, and the control matrix B_s(q_s) maps the inputs to the generalized coordinate accelerations ${\ddot{q}}_s$. Consequently, the full-model can be written in a state-space form as

${\dot{x}}_s=\left[\begin{array}{c}\hfill {\dot{q}}_s\hfill \\ \hfill D_s^{-1}\left(-C_s-G_s+B_{s}u\right)\hfill \end{array}\right]=f_s\left(x_s\right)+g_s\left(x_s\right)u,(3)$

where the state vector is denoted by $x_s={[q_s^{\top },{\dot{q}}_s^{\top }]}^{\top }$. Other planar models are obtained similarly after considering the modeling assumptions explained above. Next, the derivation of the non-instantaneous double support model follows.

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_2={[p_{tx},p_{ty},p_n]}^{\top }$ makes contact with the ground. Then two leg-ends (p₁ and p₂) 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 $x_s^{+}$ and GRF. The widely employed assumptions of rigid ground models such as instantaneous duration and inelasticity are considered to capture the system’s behavior at the impact moment. As a result, there are no changes in the joints’ position (i.e., $q_s^{-}=q_s^{+}$). Only angular velocity spikes occur at the impact moment. To cover every aspect of thruster-assisted legged locomotion and for the sake of completeness, we will briefly explain this process.

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₂, $q_d={[q_s^{\top },p_2^{\top }]}^{\top }$. Then, it is assumed that the angular momentum is conserved about the previous stance leg-end p₁ at the moment of impact. The equations of motion are re-formulated to include the impulsive GRF, $\delta J_e^{\top }{F}_e$, at each leg-end $p_e={[p_1^{\top },p_2^{\top }]}^{\top }$

$D_d\left(q_d\right){\ddot{q}}_d+C_d\left(q_d\right){\dot{q}}_d+G_d\left(q_d\right)=B_d\left(q_d\right)\left[\begin{array}{c}\hfill u_j\hfill \\ \hfill {\lambda }_{net}\hfill \end{array}\right]+\delta J_e^{\top }{F}_e(4)$

where $J_e={[{(\partial p_1/\partial q_d)}^{\top }{(\partial p_2/\partial q_d)}^{\top }]}^{\top }$ and $F_e={[F_1^{\top },F_2^{\top }]}^{\top }$. In this equation, $F_e\in {\mathbb{R}}^6$ is a Lagrange multiplier that assumes both legs are fixed to the ground at the moment of impact. Based on this assumption, we can then formulate this as a zero velocity constraint at the foot-ends, i.e., $\dot{p_e}=J_e{\dot{q}}_d=0$. Combining this with the conservation of angular momentum, the post-impact states and GRF can be solved as follows

$\left[\begin{array}{c}\hfill {\dot{q}}_d^{+}\hfill \\ \hfill F_e\hfill \end{array}\right]={\left[\begin{array}{cc}\hfill D_d\left(q_d^{-}\right)\hfill & \hfill -J_e{\left(q_d^{-}\right)}^{\top }\hfill \\ \hfill J_e\left(q_d^{-}\right)\hfill & \hfill 0^{6\times 6}\hfill \end{array}\right]}^{-1}\left[\begin{array}{c}\hfill D_d\left(q_d^{-}\right){\dot{q}}_d^{-}\hfill \\ \hfill 0^{6\times 1}\hfill \end{array}\right](5)$

The inertial matrix D_d (q_d) is a square, symmetric and positive definite, and $J_e(q_d)D_d{(q_d)}^{-1}{J}_e{(q_d)}^{\top }$ is full rank Westervelt et al. (2018). As a result, the matrix inversion shown above is possible, which yields the impact map Δ(⋅). After the impact event, both feet stay fixed to the ground. This fixation results in a non-instantaneous DS phase, which is assumed to last for a fixed time interval. This condition is enforced by considering appropriate constraints in the system and is skipped here.

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 $L_{f}y=\frac{\partial y}{\partial x}f$ and holonomic constraints 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 $\mathrm{\Gamma }=\{{[q_s^{\top },{\dot{q}}_s^{\top }]}^{\top }\in {\mathbb{R}}^{2n}|G(q_s,{\dot{q}}_s)=0\}$ where $G={[h_s^{\top }(x_s),L_{f_s}^{\top }{h}_s(x_s)]}^{\top }\in {\mathbb{R}}^{2(n-1)}$ and $rank\{L_{f_s}{h}_s(x_s)\}=n-1$ can be defined. The rank condition guarantees that Γ, i.e., $q_s=h_s^{-1}(0)$, is 1-dimensional.

Since there is only one degree of underactuation in all of our planar models, therefore, $h_s^{-1}(0)$ takes the form of a closed curve. We can find a transformation for the configuration variable, $q_s^{\prime }=H_1{q}_s$, such that the actuated coordinates are stacked on top the underactuated coordinate. With this in mind, we consider the following parametric descriptions q₁ = r₁ (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

$h_s\left(q_s\right)=H_2{H}_1{q}_s-r\left(q_n\right)=q_j-r\left(q_n\right)(6)$

Here q_j denotes the actuated joint variables, and H₂ extracts q_j from $q_s^{\prime }$. Note that $r={[r_1,\dots ,r_{n-1}]}^{\top }$ and the matrix $H_2\in {\mathbb{R}}^{(n-1)\times n}$ can take a trivial form if each joint independently is derived with a single actuator only. The classical problem of enforcing VCs considers fixed gait parameters or, at best, in event-based methods, permits dynamics-free updates of the gait parameters at the boundaries of the gait cycles. As such, by deforming the manifold Γ at the boundaries (i.e., after impact moments), gait characteristics, including average walking speed, step length, etc., are regulated.

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 $q_s=h_s^{-1}(0)$ is deformed. With this observation, we aim to apply the thrusters in two different ways, as discussed earlier in Section 1. In a nutshell, the idea is to continuously deform $q_s=h_s^{-1}(0)$ such that the following conditions are satisfied. First, we want $h_s^{-1}(0)$ to remain continuous and form closed curves, i.e. 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 $C_{eq}(q_n,{\dot{q}}_n)=0$ and inequality constraints $0\le C_{ineq}(q_n,{\dot{q}}_n)$, to be satisfied. This problem (we refer the readers to Westervelt et al. (2003) for more details) can take the following constrained ordinary differential equation form:

$\begin{array}{cc}\hfill {\ddot{q}}_n\hfill & \hfill =a_1\left(q_n\right)+a_2\left(q_n\right){\dot{q}}_n^2+\mathrm{e}\mathrm{ff}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{s}\mathrm{o}\mathrm{f}\mathrm{d}\mathrm{e}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{m}\mathrm{i}\mathrm{n}\mathrm{g}{h}_s^{-1}\left(0\right)\hfill \\ \hfill 0\hfill & \hfill \ne \frac{\partial h_s}{\partial x_s}{\left[0^{\top },{\left[D_s^{-1}\left(q_n\right)B_s\left(q_n\right)\right]}^{\top }\right]}^{\top }\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{d}\mathrm{e}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{m}\mathrm{e}\mathrm{d}{\mathrm{h}}_{\text{s}}^{-1}(0)\hfill \\ \hfill q_n\left(t+T\right)\hfill & \hfill =q_n\left(t\right)\hfill \\ \hfill {\dot{q}}_n\left(t+T\right)\hfill & \hfill ={\dot{q}}_n\left(t\right)\hfill \\ \hfill 0\hfill & \hfill =C_{eq}\left(q_n,{\dot{q}}_n\right)\hfill \\ \hfill 0\hfill & \hfill \le C_{ineq}\left(q_n,{\dot{q}}_n\right)\hfill \end{array}(7)$

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 $|\frac{F_T}{F_N}|<\mu$ can form the equality and inequality constraints, where 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 $q_s=h_s^{-1}(0)$ is deformed. Positive invariance allows finding the control input u such that when q_s (0) is on $q_s=h_s^{-1}(0)$ the velocities ${\dot{q}}_s(t)$ remains tangent to $q_s=h_s^{-1}(0)$ yielding the solutions 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 $\mathrm{\Sigma }:{[{\dot{x}}_1,{\dot{x}}_2]}^{\top }={[x_1,x_2+u]}^{\top }$ and part of ${\mathbb{S}}^1$ set – unit circle in ${\mathbb{R}}^2$ – which is $\mathrm{\Gamma }=\left\{x_1,x_2|0<x_1<1,G(x_1,x_2)=x_1^2+x_2^2-1=0\right\}$. It is straightforward to show that $u=\frac{-x_1^2}{x_2}-x_2$ will stabilize the solutions of Σ, $x(t)={[x_1(t),x_2(t)]}^{\top }$, on Γ. In this case, the restriction dynamics f|_Γ, which is prescribed by u, takes the following form $f{|}_{\mathrm{\Gamma }}={[x_1,-x_1^2/x_2]}^{\top }$. The reason that in our definition of Γ we considered part of ${\mathbb{S}}^1$ is evident here, i.e., x₂ ≠ 0.

Now, assume that we have the means to deform Γ and that the goal is to satisfy the constraint given by x₁ ≤ 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 $\frac{\partial G}{\partial x}$ and 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 $q_s=h_s^{-1}(0)$ that naturally do not violate the transversality condition (i.e., $\frac{\partial G}{\partial x}\perp ̸g(x)$).

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 $g={[0^{(n-1)\times n},{(D_s^{-1}{B}_s)}^{\top }]}^{\top }$, two possible options are considered: 1) scaling ω_iq_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

$\left[\begin{array}{c}\hfill q_j-r\left(q_n\right)+\omega \left(t\right)\hfill \\ \hfill {\dot{q}}_j-r^{\prime }\left(q_n\right){\dot{q}}_n+\dot{\omega }\left(t\right)\hfill \end{array}\right]\iff \left[\begin{array}{c}\hfill q_j-r\left(q_n\right)+{\int }_0^T\nu \left(\tau \right)\mathrm{d}\tau \hfill \\ \hfill {\dot{q}}_j-r^{\prime }\left(q_n\right){\dot{q}}_n+\nu \left(t\right)\hfill \end{array}\right](8)$

where $r^{\prime }=\frac{\partial r(q_n)}{\partial q_n}$, ω(t) and ν(t) are used for manipulating h_s (x_s) and $L_{f_s}{h}_s(x_s)$, respectively. Notice that in our approach the primer variables are time-varying. Constant priming terms yield a discrete collection of parameterized systems, i.e., 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^{\prime }={[r^{\top }(q_n,\omega ),q_n]}^{\top }$. For a fixed q_n, it is possible to show that

$B_s^{*}{D}_s\left(q_n\right){\left[r^{\prime }\left(\omega \left(t\right),q_n\right),1\right]}^{\top }(9)$

where B* = [0^1×(n−1), 1] and Eq. 9 can vanish on a point on Γ_ω which implies g_s (x_s) and $\frac{\partial h_s(x_s)}{\partial x_s}$ can be orthogonal. While the orthogonality condition may get violated at only a finite number of points on Γ_ω, these violations are not acceptable in our proposed framework and we avoid that. Without loss of generality, consider Mth order polynomials (e.g., Bezier forms with degree M as in ref. [40]) to define

$r\left(q_n,\omega \left(t\right)\right)=diag\left(\omega \left(t\right)\right)N{Q}_n(10)$

where $diag(\omega (t))\in {\mathbb{R}}^{(n-1)\times (n-1)}$ is a matrix with ω(t) as diagonal entries, $N\in {\mathbb{R}}^{(n-1)\times (M+1)}$ is a matrix of constant coefficients and $Q_n={[1,q_n,q_n^2,\dots ,q_n^M]}^{\top }$. Then Eq. 9 is re-written in the following simplified form

$D_3\left(q_n\right)diag\left(\omega \left(t\right)\right)N{Q}_n^{\prime }+D_4\left(q_n\right)=0(11)$

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 $Q_n^{\prime }=\partial Q_n/\partial q_n={[\mathrm{0,1,2}{q}_n,\dots ,M{q}_n^{(M-1)}]}^{\top }$. Meaning, using this choice of VC priming, it is possible that for some values ω(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 $\dot{y}={\dot{q}}_j-r^{\prime }(q_n){\dot{q}}_n+\omega (t)$ yields relative degree $\left\{2,\dots ,2\right\}$ on every points on Γ_ω. It is also noticeable that this choice of manipulating y has no effects on $null\left\{\frac{\partial h_s}{\partial q_s}\right\}={\left[r^{\prime \top }(q_n),1\right]}^{\top }$ which means that at least the primer has no influence over the transversality condition as long as 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

$\left[\begin{array}{c}\hfill \dot{y}\hfill \\ \hfill \ddot{y}\hfill \end{array}\right]=\left[\begin{array}{c}\hfill {\dot{q}}_j-r^{\prime }\left(q_n\right){\dot{q}}_n\hfill \\ \hfill L_{f_s}^2{h}_s\hfill \end{array}\right]+\left[\begin{array}{c}\hfill 0\hfill \\ \hfill \left[I,-r^{\prime }\left(q_n\right)\right]D_s^{-1}\left(q_s\right)B_s\hfill \end{array}\right]u+\left[\begin{array}{c}\hfill I\hfill \\ \hfill 0\hfill \end{array}\right]\omega \left(t\right)(12)$

where $L_{f_s}^2{h}_s=-\left[I,-r^{\prime }(q_n)\right]D_s^{-1}(q_s)\left(C_s(q_s){\dot{q}}_s+G_s(q_s)\right)-r^{\prime \prime }(q_n){\dot{q}}_n^2$. As a result, the closed-loop system has to possess strong disturbance rejection properties.

The roles of this disturbance term will be beneficial for us, though, as it will adjust the equilibria of h_s (x_s) and $L_{f_s}{h}_s(x_s)$ under stabilizing controllers with adjustable (and measurable) basins of attraction to successfully satisfy state, control, and GRF constraints. To do this, we need to design an update law for the primer variable vector ω(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 $u=-{\left(L_g{L}_f{h}_s(x_s)\right)}^{-1}(L_f^2{h}_s(x_s)+K_{P}y+K_D\dot{y})$ where $K_i\in {\mathbb{R}}^{(n-1)\times (n-1)}$ are constant matrices and instead will remain focused on deforming Γ_ω 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

${\mathrm{\Gamma }}_{\omega }=\left\{x_s\in X\mid h_s\left(x_s\right)=K_P^{-1}{K}_D\omega ,L_{f_s}{h}_s\left(x_s\right)=-\omega \right\}(13)$

where the equilibrium points for h_s (x_s) and $L_{f_s}{h}_s(x_s)$ under ω are obtained by solving

${\left[\begin{array}{c}\hfill h_s\left(x_s\right)\hfill \\ \hfill L_{f_s}{h}_s\left(x_s\right)\hfill \end{array}\right]}_{\omega }=\left[\begin{array}{cc}\hfill K_P^{-1}{K}_D\hfill & \hfill K_P^{-1}\hfill \\ \hfill -I\hfill & \hfill 0\hfill \end{array}\right]\left[\begin{array}{c}\hfill \omega \hfill \\ \hfill 0\hfill \end{array}\right](14)$

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 $q_s^{\prime }=H_1{q}_s={\left[q_j^{\top },p_1^{\top },q_n\right]}^{\top }$ where $q_j\in {\mathbb{R}}^m$ are actuated joint angles, m is the number of these joints in the planar models and p₁ is stance leg contact point as defined previously. The control matrix in the Euler-Lagrange equations can take the following form

$B_s=\left[\begin{array}{ccc}\hfill 0\hfill & \hfill 0\hfill & \hfill {\left(\frac{\partial p_{T,COM}}{\partial q_j}\right)}^{\top }\hfill \\ \hfill I^{m\times m}\hfill & \hfill 0\hfill & \hfill {\left(\frac{\partial p_{T,COM}}{\partial p_1}\right)}^{\top }\hfill \\ \hfill 0\hfill & \hfill I^{2\times 2}\hfill & \hfill {\left(\frac{\partial p_{T,COM}}{\partial q_n}\right)}^{\top }\hfill \end{array}\right](15)$

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\in null\left\{\frac{\partial p_{T,COM}}{\partial q_n}\right\}$. As result, the models shown in Figure 4 will be equivalent and the thruster forces 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 $\tan (q_n)=\frac{F_x^B}{F_z^B}$, the GRF components $F_x^B$ and $F_z^B$ can be shifted to the leg-end point. This leaves us with an underactuated 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., $F^B\in null\left\{\frac{\partial p_{T,COM}}{\partial q_n}\right\}$) is that thruster-assisted legged locomotion over any surface with friction coefficients as small as near zero values becomes feasible. However, the major challenge, when walking over a very slippery surface, i.e., $\tan (q_n)=\frac{F_x^B}{F_z^B}=\frac{F_T}{F_N}<<1$, is that smaller COM swing motions (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 $q_s=h_s^{-1}(0)$ is stabilizable, for the previously defined $B_s^{*}$ the following will hold $B_s^{*}{D}_s(q_n){\left[r^{\prime \top }(\omega (t),q_n),1\right]}^{\top }=\alpha$ where α ≠ 0. As a result, the following equation can be obtained at every point on Γ_ω and it governs the dynamics of the ${[q_n,{\dot{q}}_n]}^{\top }$ curve

${\ddot{q}}_n=-{\alpha }^{-1}\left({\beta }_1{\left({\dot{q}}_n\right)}^2+{\beta }_2\right)(16)$

In this equation ${\beta }_2=B_s^{*}{G}_s(q_n)$ and β₁ is given by

${\beta }_1=B_s^{\ast }\left(D_s\left(q_n\right){\left[r^{\prime \prime \top },0\right]}^{\top }+{\mathrm{\Sigma }}_{i=1}^n\left[r^{\prime \top },1\right]Q_i\left(q_n\right){\left[r^{\prime \top },1\right]}^{\top }\right)(17)$

where Q_i (q_n) are the Christoffel Symbols. A similar algebraic relationship for the constraints ${[u_j^{\top },{\lambda }_{net}^{\top }]}^{\top }$ is obtained. Following the relationships given by Eqs. 6 and 8, constraints equations can be written in terms of y and $\dot{y}$, which is skipped here to avoid a cluttered notation.

Next, we will steer y and $\dot{y}$ using the primer variable ω(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-$\dot{y}$ space. While a similar analysis is possible in the q_n-${\dot{q}}_n$ space, the y-$\dot{y}$ space can offer some geometric tools that can help formulate the problem as a classical motion planning problem. Obviously, the result of such adjustments would be the deviation from the nominal gait trajectories, which can affect, e.g., hybrid invariance. We will use the thrusters’ actions to deal with these issues.

Since we assumed a pre-stabilized system – in fact all of the above derivations only make sense if $q_j=r(q_n)+{\int }_0^t\omega (\tau )d\tau$ and ${\dot{q}}_j=r^{\prime }(q_n){\dot{q}}_n+\omega (t)$ – it is reasonable to evaluate the constraints $c_l\le {[u_j^{\top },F_1^{\top }]}^{\top }\le c_u$ (c_l and c_u are constraint lower and upper bounds) based on the steady-state solutions, i.e., y_ω and ${\dot{y}}_{\omega }$, and ignore the transient solutions, i.e., ${\left[y^{\top },{\dot{y}}^{\top }\right]}^{\top }=e^{At}{z}_0+{\int }_0^t{e}^{A(t-\tau )}{B}_{\omega }w(\tau )d\tau$, where $z_0={\left[y_0^{\top },{\dot{y}}_0^{\top }\right]}^{\top }$ is the initial condition and

$A=\left[\begin{array}{cc}\hfill 0^{\left(n-1\right)\times \left(n-1\right)}\hfill & \hfill I^{\left(n-1\right)\times \left(n-1\right)}\hfill \\ \hfill -K_P\hfill & \hfill -K_D\hfill \end{array}\right](18)$

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 ${\dot{y}}_{\omega }$ is guaranteed – currently we assume that this holds true and later we ensure the GAS property is achieved – by the GAS property of the controller, any transient constraint violations are quickly compensated.

Other than simplifying the nonlinear constraint satisfaction problem given in Eq. 7, considering ${[y_{\omega }^{\top },{\dot{y}}_{\omega }^{\top }]}^{\top }$ has another interesting result which will be explained below. Consider the set

$Y_{\omega }=\left\{z\in {\mathbb{R}}^{2\left(n-1\right)}|\left[I^{\left(n-1\right)\times \left(n-1\right)},-K_P^{-1}{K}_D\right]z=0\right\}(19)$

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 $z_{\omega }={[y_{\omega }^{\top },{\dot{y}}_{\omega }^{\top }]}^{\top }$ in the set defined by Y_ω. Then, one can slide ${[y_{\omega }^{\top },{\dot{y}}_{\omega }^{\top }]}^{\top }$ over the hyperplane defined by 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., $\mathrm{m}\mathrm{i}\mathrm{n}\left\{‖{[y_{\omega }^{\top },{\dot{y}}_{\omega }^{\top }]}^{\top }-0‖\right\}$. The last condition implies that the modified gaits should closely match the predefined gaits. The nominal (predefined) gait parameters are obtained off-line without taking into account the constraint equations given by Eq. 7.

5.4 Primer Variable Vector, ω(t), Update Policy and Achieving Global Asymptotic Stability Property Over Y_ω

Consider the support hyperplane given by

$c_i={\left(\nabla C_i\left(z,\omega \left(t\right)\right){|}_{z_{\omega }^{\prime }}\right)}^{\top }\left(z-z_{\omega }^{\prime }\right)=0(20)$

at the point $z_{\omega }^{\prime }$ where the nonlinear constraint C_i(z, ω(t)) – i.e., the ith 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

FIGURE 5. Geometric interpretation of the level set {ζ_w|V ≤ V_max}. For a constant reference r, the update law $\dot{w}$ to the manipulated reference 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

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., $\frac{\partial G}{\partial x}\perp ̸g(x)$ at any point on Γ; 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 $z_{\omega }^{\top }=0$ and constraint error margins (shown in Figure 6).

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)‖₂, ‖ω(t)‖₂ ≠ 0, and step length, $\delta \in {\mathbb{R}}^{+}$.

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:

$\dot{\omega }\left(t\right)=k\rho \delta (21)$

where $k\in {\mathbb{R}}^{+}$ is the tuning parameter.

Consider ζ = P^−0.5z and ζ_ω = P^−0.5z_ω where P is a positive definite matrix. The maximum level set value V_max,i (ω(t)) for the quadratic Lyapunov function $V(\zeta ,{\zeta }_{\omega })=V(\zeta ,\omega (t))={\left(\zeta -{\zeta }_{\omega }\right)}^{\top }P\left(\zeta -{\zeta }_{\omega }\right)$ is obtained by doing the change of variable given above. In this way, 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 $c_i(z,\omega (t))=c_{i,1}^{\top }z+c_{i,2}$ to the point z_ω can be better imagined when 2-norm ‖.‖₂ is used, shown in Figure 6. The unit vector $e_{c_i}$ along the shortest distance between z_ω and the constraint hyperplane c_i = 0 is given by $e_c=\nabla C_i{|}_{z_{\omega }}/\Vert \nabla C_i{|}_{z_{\omega }}{\Vert }_2$ and the largest level set value is given by $V_{max,i}(\omega (t))=c_{i,2}/\sqrt{c_{i,1}^{\top }{c}_{i,1}}$.

Now consider the quadratic norm ‖ζ‖_P. Notice that using the change of variable ζ = P^−0.5z 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 $e_{c_i}$ which now is the maximum distance traveled towards the constraint hyperplane from ζ_ω in a unit ball of ‖.‖_P is given by

$e_{c_i}=-{\left({\gamma }^{\top }{P}^{-1}\gamma \right)}^{-0.5}{P}^{-1}\gamma (22)$

where $\gamma =\nabla C_i{|}_{z_{\omega }}$. To obtain the equation, we considered the dual norm of the quadratic norm ‖γ‖_p which is ‖γ‖_* = ‖P^−0.5γ‖₂ and this algebraic relationship $e_{c_i}=-P^{-1}\gamma /\Vert \gamma {\Vert }_{*}$.

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 $\dot{V}(\zeta ,\omega (t))$ is given by

$\dot{V}\left(\zeta ,\omega \left(t\right)\right)=-{\left(\zeta -{\zeta }_{\omega }\right)}^{\top }Q\left(\zeta -{\zeta }_{\omega }\right)+2{\left(\zeta -{\zeta }_{\omega }\right)}^{\top }P{\dot{\zeta }}_{\omega }(23)$

where Q is a positive definite matrix. For $\dot{\omega }(t)\ne 0$, there exists 0 < κ ≪ 1 such that

$\Vert 2{\left(\zeta -{\zeta }_{\omega }\right)}^{\top }{P}^{0.5}{A}^{-1}{B}_{\omega }\frac{\dot{\omega }\left(t\right)}{\kappa }\Vert \le {\lambda }_{\mathit{min}}\left(Q\right)\Vert \zeta -{\zeta }_{\omega }{\Vert }^2(24)$

As a result, the resulting inequality $\dot{V}(u,\omega (t))<(\kappa -1){\lambda }_{\mathit{min}}(Q)\Vert z{\Vert }^2<0$ holds because it is assumed κ ≪ 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 $x_d={[q_d^{\top },{\dot{q}}_d^{\top }]}^{\top }$.

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 kth 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 $J_d\ddot{q}+{\dot{J}}_d\dot{q}=0$ results in the following differential algebraic equation

$\left[\begin{array}{cc}\hfill D_d\left(q_d\right)\hfill & \hfill -J{\left(q_d\right)}^T\hfill \\ \hfill J_d\left(q_d\right)\hfill & \hfill 0^{\left(6\times 6\right)}\hfill \end{array}\right]\left[\begin{array}{c}\hfill {\ddot{q}}_d\hfill \\ \hfill F_e\hfill \end{array}\right]=\left[\begin{array}{c}\hfill B_d\eta -H_d\left(q_d,{\dot{q}}_d\right)\hfill \\ \hfill -{\dot{J}}_d\left(q_d,{\dot{q}}_d\right)\hfill \end{array}\right](25)$

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 ${\dot{x}}_d=f_d(x_d)+g_d(x_d)\eta$. This model is discretized at each sample time. To cast this as an NMPC problem a cost function ϕ(x_d, η) to be minimized is formulated resulting in the following optimization problem

$\begin{array}{cc}\hfill \underset{\eta \left[k\right]}{\mathbf{m}\mathbf{i}\mathbf{n}}\varphi \left(x_d,\eta \right)\hfill & \hfill =\sum _{k=1}^N{W}_r\Vert \left|\right.\left|\right.\left(x_d\left[k\right]-r_d\left[k\right]\right)\left|\right.{\left|\right.}_2^2+W_{\eta }\left|\right.\left|\right.\delta \eta \left[k\right]\left|\right.{\left|\right.}_2^2\hfill \\ \hfill \mathbf{s}\mathbf{u}\mathbf{b}\mathbf{j}.\mathbf{t}\mathbf{o}:\hfill \\ \hfill \hfill & \hfill C_{eq}\left(x_d,\eta \right)=0\hfill \\ \hfill \hfill & \hfill C_{ineq}\left(x_d,\eta \right)=0\hfill \end{array}(26)$

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 $x_s^{+}$, and the discretized dynamics ${\dot{x}}_d$.

$C_{eq}\left(x_d,\eta \right):\left\{\begin{array}{cc}{x}_d\left[0\right]\hfill & =R_s^d{x}_s^{+}\hfill \\ x_d\left[k+1\right]\hfill & =f\left(x_d\left[k\right]\right)+g\left(x_d\left[k\right]\right)\eta \left[k\right]\hfill \end{array}\right.(27)$

In Eq. 27, the matrix $R_s^d$ swaps the roles of the legs. The inequality equations include the difference between the final DS state and initial SS state, GRF constraints and the limits imposed on the state and input vectors.

$C_{ineq}\left(x_d,\eta \right):\left\{\begin{array}{cc}\hfill |x_d\left[N\right]-x_s\left[0\right]|& \le ϵ\hfill \\ \hfill \left|\right.|F_T\left[k\right]|\left|\right.& <\mu \left|\right.{|F}_N\left[k\right]|\left|\right.\hfill \\ \hfill 0& <{|F}_N\left[k\right]|\hfill \\ \hfill \left|\right.{|x}_d\left[k\right]\left|\right.|& <x_{dmax}\hfill \\ \hfill \left|\right.|\eta \left[k\right]\left|\right.|& <{\eta }_{max}\hfill \end{array}\right.(28)$

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.

FIGURE 7

FIGURE 7. An overview of Harpy’s controller architecture.

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.

TABLE 1

TABLE 1. Model parameters.

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

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

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₀) 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

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

FIGURE 11

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.

FIGURE 12

FIGURE 12. Illustrates the reference joint angle trajectories in the 3-link ROM.

FIGURE 13

FIGURE 13. Illustrates the reference joint velocity trajectories in the 3-link ROM.

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₁-${\dot{q}}_1$) corresponding to each control limits are shown in Figure 10. It can be seen that $h_s^{-1}(0)=q_s$ is deformed to satisfy the constraints which can affect hybrid invariance in the locomotion system.

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

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

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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.

Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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.

Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

Google Scholar

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

Google Scholar

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.

Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

Xin, G., Wolfslag, W., Lin, H.-C., Tiseo, C., and Mistry, M. (2020). An Optimization-Based Locomotion Controller for Quadruped Robots Leveraging Cartesian Impedance Control. Front. Robot. AI 7, 48. doi:10.3389/frobt.2020.00048

PubMed Abstract | CrossRef Full Text | Google Scholar