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

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.

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.

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.

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.

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.

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.

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 ∈ 1 ( s t a n c e l e g ) , 2 ( s w i n g l e g ) and stance leg-end positions ( p 1 = [ p t x , p t y , p n ] ⊤ ) . All of these variables are stacked inside the vector q s = [ q y , q p , q r , q i T h i g h , q i K n e e , q i H i p , q i T h r u s t , p 1 ⊤ ] ⊤ to form the states vector x s = [ q s ⊤ , q ̇ s ⊤ ] ⊤ .

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 ⊤ , λ i ⊤ , F 1 ⊤ ] ⊤ embodies the joint control actions u j = [ u i T h i g h , u i K n e e , u i H i p , u i T h r u s t ] ⊤ , 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 λ n e t = [ F B , M B ] ⊤ . This wrench vector consists of the net thruster force and moment on the body frame. This map is given by the adjoint transformation λ n e t = ∑ i A d i B ⊤ λ i , which can be expanded as λ n e t = ∑ i = 1 2 R i 0 R i S d i R i λ 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 d ⃗ i × F i d i 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.

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.

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.

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 K ( q s , q s ̇ ) and potential V ( q s ) energies of the system which leads to L ( q s , q s ̇ ) = K ( q s , q s ̇ ) − V ( q s ) . The resulting equations of motion are expressed in vector form as following D s q s q ̈ s + C s q s , q ̇ s q ̇ s + G s q s = B s q s u , (2) where D _s (q _s ) is the symmetric inertial matrix and is only dependent on q _s , the C s ( q s , q ̇ s ) 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 q ̈ s . Consequently, the full-model can be written in a state-space form as x ̇ s = q ̇ s D s − 1 − C s − G s + B s u = f s x s + g s x s u , (3) where the state vector is denoted by x s = [ q s ⊤ , q ̇ s ⊤ ] ⊤ . Other planar models are obtained similarly after considering the modeling assumptions explained above. Next, the derivation of the non-instantaneous double support model follows.

The end of the SS phase is marked by an impulsive impact event when the swing leg-end p 2 = [ p t x , p t y , p n ] ⊤ 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 ⊤ , p 2 ⊤ ] ⊤ . 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, δ J e ⊤ F e , at each leg-end p e = [ p 1 ⊤ , p 2 ⊤ ] ⊤ D d q d q ̈ d + C d q d q ̇ d + G d q d = B d q d u j λ n e t + δ J e ⊤ F e (4) where J e = [ ( ∂ p 1 / ∂ q d ) ⊤ ( ∂ p 2 / ∂ q d ) ⊤ ] ⊤ and F e = [ F 1 ⊤ , F 2 ⊤ ] ⊤ . In this equation, F e ∈ 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., p e ̇ = J e q ̇ d = 0 . Combining this with the conservation of angular momentum, the post-impact states and GRF can be solved as follows q ̇ d + F e = D d q d − − J e q d − ⊤ J e q d − 0 6 × 6 − 1 D d q d − q ̇ d − 0 6 × 1 (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 ) ⊤ 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.

To motivate our idea, we will consider Σ : [ x ̇ 1 , x ̇ 2 ] ⊤ = [ x 1 , x 2 + u ] ⊤ and part of S 1 set – unit circle in R 2 – which is Γ = x 1 , x 2 | 0 < x 1 < 1 , G ( x 1 , x 2 ) = x 1 2 + x 2 2 − 1 = 0 . It is straightforward to show that u = − x 1 2 x 2 − x 2 will stabilize the solutions of Σ, x ( t ) = [ x 1 ( t ) , x 2 ( t ) ] ⊤ , on Γ. In this case, the restriction dynamics f|_Γ, which is prescribed by u, takes the following form f | Γ = [ x 1 , − x 1 2 / x 2 ] ⊤ . The reason that in our definition of Γ we considered part of 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 ∂ G ∂ 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., ∂ G ∂ x ⊥̸ g ( x ) ).

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 ) × n , ( D s − 1 B s ) ⊤ ] ⊤ , two possible options are considered: 1) scaling ω _i q _n and 2) shifting q _n + ω _i . We are interested in knowing how the stabilizability of Γ _ω is affected if these two options are applied to manipulate the state-dependent equilibrium point of the actuated coordinates. In equilibrium shifting scenario (case ii), notice that the position and velocity priming can be trivially equivalent as is evidently seen below q j − r q n + ω t q ̇ j − r ′ q n q ̇ n + ω ̇ t ⇔ q j − r q n + ∫ 0 T ν τ d τ q ̇ j − r ′ q n q ̇ n + ν t (8) where r ′ = ∂ r ( q n ) ∂ 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.

Consider the configuration variable vector q s ′ = H 1 q s = q j ⊤ , p 1 ⊤ , q n ⊤ where q j ∈ 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 = 0 0 ∂ p T , C O M ∂ q j ⊤ I m × m 0 ∂ p T , C O M ∂ p 1 ⊤ 0 I 2 × 2 ∂ p T , C O M ∂ q n ⊤ (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 ∈ n u l l ∂ p T , C O M ∂ q n . 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 ) = 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 ∈ n u l l ∂ p T , C O M ∂ q n ) 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 ) = F x B F z B = 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 ) r ′ ⊤ ( ω ( t ) , q n ) , 1 ⊤ = α where α ≠ 0. As a result, the following equation can be obtained at every point on Γ _ω and it governs the dynamics of the [ q n , q ̇ n ] ⊤ curve q ̈ n = − α − 1 β 1 q ̇ n 2 + β 2 (16)

In this equation β 2 = B s * G s ( q n ) and β ₁ is given by β 1 = B s ∗ D s q n r ′ ′ ⊤ , 0 ⊤ + Σ i = 1 n r ′ ⊤ , 1 Q i q n r ′ ⊤ , 1 ⊤ (17) where Q _i (q _n ) are the Christoffel Symbols. A similar algebraic relationship for the constraints [ u j ⊤ , λ n e t ⊤ ] ⊤ is obtained. Following the relationships given by Eqs. 6 and 8, constraints equations can be written in terms of y and y ̇ , which is skipped here to avoid a cluttered notation.

Next, we will steer y and 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- y ̇ space. While a similar analysis is possible in the q _n - q ̇ n space, the y- 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 ) + ∫ 0 t ω ( τ ) d τ and q ̇ j = r ′ ( q n ) q ̇ n + ω ( t ) – it is reasonable to evaluate the constraints c l ≤ [ u j ⊤ , F 1 ⊤ ] ⊤ ≤ c u (c _l and c _u are constraint lower and upper bounds) based on the steady-state solutions, i.e., y _ω and y ̇ ω , and ignore the transient solutions, i.e., y ⊤ , y ̇ ⊤ ⊤ = e A t z 0 + ∫ 0 t e A ( t − τ ) B ω w ( τ ) d τ , where z 0 = y 0 ⊤ , y ̇ 0 ⊤ ⊤ is the initial condition and A = 0 n − 1 × n − 1 I n − 1 × n − 1 − K P − K D (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 y ̇ ω 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 ω ⊤ , y ̇ ω ⊤ ] ⊤ has another interesting result which will be explained below. Consider the set Y ω = z ∈ R 2 n − 1 | I n − 1 × n − 1 , − K P − 1 K D z = 0 (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 ω = [ y ω ⊤ , y ̇ ω ⊤ ] ⊤ in the set defined by Y _ω . Then, one can slide [ y ω ⊤ , y ̇ ω ⊤ ] ⊤ 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., m i n [ y ω ⊤ , y ̇ ω ⊤ ] ⊤ − 0 . 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.

Consider the support hyperplane given by c i = ∇ C i z , ω t | z ω ′ ⊤ z − z ω ′ = 0 (20) at the point z ω ′ 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.

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 ω ⊤ = 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, δ ∈ 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: ω ̇ t = k ρ δ (21) where k ∈ R + is the tuning parameter.

Consider ζ = P ^−0.5 z and ζ _ω = P ^−0.5 z _ω where P is a positive definite matrix. The maximum level set value V _max,i (ω(t)) for the quadratic Lyapunov function V ( ζ , ζ ω ) = V ( ζ , ω ( t ) ) = ζ − ζ ω ⊤ P ζ − ζ ω 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 , ω ( t ) ) = c i , 1 ⊤ 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 = ∇ C i | z ω / ‖ ∇ C i | z ω ‖ 2 and the largest level set value is given by V m a x , i ( ω ( t ) ) = c i , 2 / c i , 1 ⊤ c i , 1 .

Now consider the quadratic norm ‖ζ‖ _P . Notice that using the change of variable ζ = P ^−0.5 z in the constraint equation c _i (z, ω(t)) = 0, 2-norm and the quadratic norm become interchangeable. In this way, the maximum level set value for any candidate, quadratic Lyapunov function can be obtained. Therefore, the step length in Eq. 21 is given by δ = V _max,i (ω(t)) − V (ζ, ω(t)). After the change of variable, the unit vector 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 = − γ ⊤ P − 1 γ − 0.5 P − 1 γ (22) where γ = ∇ C i | z ω . 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 γ / ‖ γ ‖ * .

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 V ̇ ( ζ , ω ( t ) ) is given by V ̇ ζ , ω t = − ζ − ζ ω ⊤ Q ζ − ζ ω + 2 ζ − ζ ω ⊤ P ζ ̇ ω (23) where Q is a positive definite matrix. For ω ̇ ( t ) ≠ 0 , there exists 0 < κ ≪ 1 such that ‖ 2 ζ − ζ ω ⊤ P 0.5 A − 1 B ω ω ̇ t κ ‖ ≤ λ min Q ‖ ζ − ζ ω ‖ 2 (24)

As a result, the resulting inequality V ̇ ( u , ω ( t ) ) < ( κ − 1 ) λ min ( Q ) ‖ z ‖ 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.

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 ⊤ , q ̇ d ⊤ ] ⊤ .

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 q ̈ + J ̇ d q ̇ = 0 results in the following differential algebraic equation D d q d − J q d T J d q d 0 6 × 6 q ̈ d F e = B d η − H d q d , q ̇ d − J ̇ d q d , q ̇ d (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 x ̇ d = f d ( x d ) + g d ( x d ) η . 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 m i n η k ϕ x d , η = ∑ k = 1 N W r ‖ x d k − r d k 2 2 + W η δ η k 2 2 s u b j . t o : C e q x d , η = 0 C i n e q x d , η = 0 (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 x ̇ d . C e q x d , η : x d 0 = R s d x s + x d k + 1 = f x d k + g x d k η k (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 i n e q x d , η : | x d N − x s 0 | ≤ ϵ | F T k | < μ | F N k | 0 < | F N k | | x d k | < x d m a x | η k | < η m a x (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.