I Introduction
Simplistic conservative models of legged locomotion, in which no energy is lost during a stride, are a powerful tool for both the analysis of human and animal gaits in nature and the design and control of legged robots
[16, 9, 13, 4]. With just a few degrees of freedom and a small number of physical parameters, these models can accurately predict the preferred locomotion patterns of humans [14] and provide useful templates for energyefficient robot motions [6].Despite the benefits of such models, the field is still lacking a unified approach that systematically takes advantage of the conservative nature of these models to identify and characterize the different types of periodic motions available. This becomes even more important given that the same model can exhibit multiple modes of locomotion (e.g., walking, hopping, and running). To the best of our knowledge, past works have only developed results for specific conservative models and gait type [7, 20, 8, 17]
and not a class of energetically conservative systems with hybrid dynamics and multiple modes of locomotion. The goal of this paper is to create a mathematical framework rooted in the theory of hybrid dynamical systems and nonlinear dynamics to model, classify, and create periodic motions for energetically conservative models (ECMs) of legged systems.
To this end, we generalize the methodology introduced in [7] and carefully embed it into a mathematical framework for general ECMs of legged systems. We prove that families of gaits exist for such systems and highlight the role of energy in providing a coordinatefree parameterization for these families. In order to make the approach practical, we present algorithms for identifying families of gaits through the use of numerical continuation methods and introduce a number of details for their implementation. Among others, these details include projecting the state space to the subspace of periodic motions, establishing the necessary condition for the Delassus’ matrix to preserve energy across impacts, introducing the use of additional (holonomic) constraints to avoid singular dynamics, embedding the conservative system in a oneparameter family of dissipative systems and transitioning from an eventdriven formulation to a timebased formulation.
This paper can be considered to be a direct extension of [7] which showed that a simple model exhibits all common bipedal gaits and that these form continuous families of gaits in the biped’s space of trajectories. These periodic motions all emerged from a onedimensional (1D) family of hoppinginplace gaits. Other gaits, such as walking and running, were connected to these through a series of bifurcations. Furthermore, our work builds upon the oneparameter families of periodic orbits in smooth ECMs as they are the main subject in [24] and [1]. While [24] provides conditions for the existence of this family, [1] revisits concepts of socalled Nonlinear Normal Modes (NNMs) that aim to find analytic expressions of invariant lowerdimensional submanifolds. Herein, NNMs are explicitly parameterized representations of 1D manifolds that emanate from exploiting the system’s state dependencies inflicted by the conservation of energy.
Ii Theory
Iia Dynamics of Legged Systems
In our work, we consider rigid body systems subject to unilateral contact without sliding, which are commonly used to model legged robotic systems. The state of such a system is given by the vector
, where is the number of its degrees of freedom and is the tangent bundle of the configuration space . The dynamics are given by the differentialalgebraic equation:(1a)  
(1b) 
with mass matrix , elastic forces , gravitational, centrifugal, and coriolis forces , and nonconservative torques . Since we exclude sliding contacts, the active unilateral contacts are reflected in equation (1b) as holonomic bilateral constraints , with constraint Jacobian and associated constraint forces . Hence, the set of constraints (1b) depends on which contacts are active in the current contact configuration (or mode). This mode is denoted by , where is the finite set of possible contact configurations. We assume that the constraints are independent and that, consequently, is full rank and the forces can be uniquely solved for (Theorem 5.1 [5]).
Transitions into different modes are triggered as individual contacts open and close. All possible transitions are described by the set , forming a directed graph structure over (Definition 2.2 [12]). A transition is triggered by event functions , with and , in the set . For legged systems, closing contacts (touchdown) are determined kinematically when a foot that is not constrained by equation (1b) goes below the ground. Opening contacts (liftoff) are triggered kinetically when a normal contact force component in becomes negative.
Whenever triggers an event, the projection into the new mode is defined by the discrete map with and . Since this transition is instantaneous in time, and are the states right before and after an event, respectively. For mechanical systems, this discrete map only alters velocities and the projection depends exclusively on the new contact configuration:
(2) 
As it is a common choice in legged systems [10], we only consider plastic collisions in equation (2). For such collisions, the postcollision velocities must satisfy: . For this special case, the discrete map is given by: . In the field of nonlinear mechanics, the matrix is known as the Delassus’ matrix of contact configuration [5]. It describes the inertial coupling in the active constraint space^{1}^{1}1
is called the constrained contact inertia tensor in
[12]..The resulting equations of motion are those of a hybrid dynamical system and can be written in the form:
(3) 
where follows from the differentialalgebraic equation (1).
IiB Conservation of Energy
In this paper, we only consider energetically conservative models (ECMs) in which the total energy (given by the sum of kinetic energy and potential energy ) is constant.
Lemma II.1.
The hybrid dynamics are energetically conservative if and only if

all forces in the continuous dynamics of equation (1a) are conservative forces and thus, , and

for all discrete maps it holds . This implies , since the discrete dynamics do not change the value of ; i.e., .
IiC Periodic Solutions in Energetically Conservative Hybrid Dynamical Systems
Before addressing properties of solutions to for conservative systems, we establish the following notation.
Definition (Hybrid Periodic Flow).
The flow describes a solution of starting from an initial condition such that and . A solution is periodic, if there exists a period time , such that
(4) 
Assumption (Admissible Flow).
We make Assumption A1 due to the ambiguity of the flow’s time derivative at a transition between modes, Assumptions A2 and A3 to exclude equilibria, grazing and Zeno behaviors, and Assumption A4, such that the contact sequence is fixed under local changes of .
Altogether, Assumptions A1A4 yield a welldefined flow which is unique as the semigroup property holds for all admissible flows. Specifically, under Assumption A4, the fundamental solution matrix
(5) 
Definition (Monodromy Matrix).
The local linearization of a periodic solution is called the monodromy matrix.
The monodromy matrix is an important tool to study the stability and local existence of admissible periodic flows (Chapter 7.1.1 [15]). For ECMs, it holds that:
(6)  
(7) 
for in mode . Equation (6) is the well known freedom of phase in autonomous systems, as any disturbance along the flow will remain on the same periodic motion in (Theorem 2 [24]). Furthermore, for conservative systems, total energy is flowinvariant: . For admissible periodic flows, this yields the property in equation (7) (Chapter 2.4. [24]).
Lemma II.2.
Outside of an equilibrium (A2) and must be nonzero for a mechanical system. Since holds for all , these two vectors are perpendicular:
(8) 
IiD Connected Components of Energetically Conservative Gaits
The purpose of this work is to show connections between different periodic motions that we will refer to as different gaits. To eliminate the freedomofphase that is inherent to any autonomous system, we introduce an anchor constraint to further specify the solution that constitutes a specific gait:
Definition (Gait).
A gait is a periodic solution that also fulfills the anchor constraint , where is a smooth function for which the transversality condition holds.
Remark II.1.
One can always find an anchor constraint which fulfills the transversality condition on periodic solutions that are not equilibria.
Proposition (Family of Gaits).
In the vicinity of a energetically conservative gait there exist neighboring gaits.
Proof.
Due to the periodicity, it must hold:
(9) 
where we abuse the notation of the period to indicate its general dependency on . Using the implicit function theorem, we get:
(10) 
To explore neighboring gaits, we perturb the initial state of the periodic solution (4) by an infinitesimal :
(11) 
A firstorder approximation of equation (11) yields
(12) 
where
is the identity matrix.
Remark II.2.
Remark II.3.
This Proprosition is an extension of Theorem 4 in [24] that proves that for smooth conservative dynamics, orbits are dense in the state space .
Remark II.4.
What was shown here for energy can be extended to other flow invariant functions^{2}^{2}2These socalled first integrals are considered in [19] for smooth systems.
. For example, in some mechanical systems, linear or angular momentum may be conserved. The existence of such invariants can then lead to additional left eigenvectors as in equation (
7) and hence in the kernel of in equation (12).We propose to parameterize the resulting families of connected gaits by energy level . While other parameterizations are possible (e.g., using a state variable, such as speed [7]), gives a more general coordinatefree parameterization for ECMs, since gaits are inherently constrained to an equipotential surface (Lemma II.2). This parameterization is reflected in:
(13) 
with its derivative
The set of all solutions (with admissible flow) to equation (13) for all possible energy levels constitutes the gait space .
Definition (Regular Point).
We call a solution of an implicit function with a regular point if has maximum rank.
While has full rank, there exists a set of regular points that form a locally defined 1D submanifold . Since each point of represents a periodic motion, is also called a generator for a twodimensional invariant set of solutions in the state space [1].
Definition (Generators).

A set is pathconnected if for any two points , there exists a continuous function such that and .

A set is called a generator if it is pathconnected and all points are regular.
Generators can border to a point which does not belong to an admissible solution of (e.g., solutions with grazing or with a change in contact sequence). We refer to these points as inadmissible points (IP) (Fig. 1). Alternatively, they can border to a point for which becomes rank deficient. These singularities either constitute turning points (TP) (i.e., extremal values for the parameter ) or bifurcations (BP) in which the periodic solutions of equation (13) are no longer distinct. Both types of singularities connect different generators to form a connected component of the gait space.
Definition (Connected Components).
A set is a connected component of if is pathconnected and is maximal with respect to inclusion (Definition 2.3 [23]).
How such connected components of the gait space can be efficiently computed, will be discussed in the following.
Iii Implementation
Iiia Constructing Conservative Models
To implement an energetically conservative model of legged locomotion, the properties L1 and L2 in Lemma IIB must be fulfilled. L1 can be easily satisfied by implementing ideal constraints and omitting additional joint torques in equations (1). To satisfy L2 at touchdown events, we have to account for the changes in velocity, yielding:
(14) 
Using the projection in equation (2), we can write this as
(15) 
In the general case, energy conservation would only be possible if the inverse Delassus’ matrix were zero. Loosely speaking, this is because inertia and masses involved in the projection need to vanish to conserve energy. This is problematic, as this requirement leads to singularities in the systems mass matrix .
Instead, we consider vanishing masses and inertias only as a limiting case. That is, with some abuse of notation, we define a parameterized mass matrix , with parameter such that the Delassus’ matrix reads as . This parameterization must yield
(16) 
Considering equation (15), the mechanical system is only energetically conservative in the limit of . As pointed out in chapter 2.3. of [12], massless appendages of a robot possibly yield an inconsistent relationship between accelerations and net forces in equation (1a). Hence, any rank deficiency of the mass matrix has to be corrected by constraints (, ) to ensure unique, finite dimensional dynamics. With this, it is possible to cancel out appearing singularities in the inverse mass matrix by introducing a parametric scaling with in , and such that equation (1a) can be stated as
(17) 
The resulting conservative vector field, defined by equation (17), is smooth and complete in the analytic limit of . This implies the existence and uniqueness of a flow in each hybrid mode through any initial condition.^{3}^{3}3In other words: while can become singular in the limit of , the products , , and remain finite.
IiiB Numerical Exploration
The goal of our implementation is to solve the implicit function (13) in a systematic fashion to obtain the connected component . Our primary tool for the computation of generators are numerical continuation methods [2].
The issue with numerically solving equation (13) is that it has constraints but only decision variables in and . In theory, this is no problem, as the equations in (13) are not independent due to the energetically conservative nature of the dynamics [22], as was shown above. In practice, however, this can cause issues, as fluctuations in energy can be introduced during numerical integration. When this is the case, equation (13) may not be solvable with only decision variables. To tackle this issue, we use the approach reported in [24] and add a parameter to our continuous dynamics in (3):
(18) 
With the new representation (18), the conservative system is embedded in a oneparameter family of dissipative dynamics . Analytically, a periodic orbit only exists for a vanishing perturbation (Lemma 1 [24]). Hence, solutions of with are periodic solutions of the underlying conservative system. In the numerical computation of gaits, however, we might obtain solutions with a small to compensate for small energy losses caused by numerical damping in the integration schemes.
Gaits of legged systems, are not necessarily periodic in all states. In particular, the horizontal position is aperiodic to allow for forward motion. Hence, to relax the periodicity constraint (4), we split the state into a periodic part and a nonperiodic part by introducing the constant orthonormal selection matrix .
Due to Assumption A4, the contact sequence for cannot change within a connected component . It is therefore predefined and we can explicitly assign a time duration to each mode . This allows us to move away from an eventdriven evaluation of . In this approach, the event constraints become explicit components of the root function , rather than being implicitly stated in the set . This change greatly facilitates the computation of the derivatives in . Hence, a periodic solution for a given can be obtained numerically by solving the rootfinding problem :
(19) 
where is the total number of modes in a given sequence such that , and the initial states of each mode are defined recursively:
With , we refer to the Jacobian of as . In equation (19), the period is no longer an explicit decision variable but results from .
In addition to the implicit equation (19), we define an extended root function that also includes as a free variable:
(20)  
(21) 
If is a regular point of , then characterizes a locally defined 1D solution manifold. The function is well suited for a pseudoarclength continuation which is utilized to compute generators. This approach employs a predictorcorrector (PC) method with a variable step size (Chapter 6.1 [2]), which takes small iterative steps in the tangent space of to locally trace the solution curve of regular points. This tangent space is equivalent to the kernel of at a regular point of equation (20), with the tangent vector :
(22) 
As the curve can be locally pursued in two directions, defines positive orientation [2].
In this process, the crossing of simple (codimensionone^{4}^{4}4A simple or codimensionone bifurcation point is defined by a loss of rank in , i.e., .) bifurcations are detected by a flip in direction of the tangent vector (i.e., ) [2]. The detection of turning points (TP) follows from a change in sign of (i.e., ), in which remains a regular point of equation (20).
In Algorithm 1, the curve is traversed in both directions until a special point is detected. Special points are the result of a PCstep that has crossed a BP, TP, or IP. Herein, nonsuccessful PCsteps (e.g., divergence in Newton’s method) are also considered inadmissible (IP). The algorithm returns the new generator and its associated TPs and BPs. The curve has at most 2 limiting special points. As mentioned previously, TPs and BPs are singular points that connect to different generators . Algorithm 2 constructs a subset of the space of connected components. It utilizes a breadthfirstsearch to explore different generators given the location of connected TPs and BPs. Locations of regular points in the neighborhood of simple bifurcations can be found with the bifurcation equation (Chapter 8.3 [2]). As indicated above, it is essential to have a problem specific starting point that solves equation (19) and is regular.
Iv Example: OneLegged Hopper
Iva Model Description
In this section, we highlight the application of our method to a SLIPlike onelegged hopper introduced in [7] with passive swing leg dynamics that are created by a torsional hip spring (Fig. 2). Here, however, it is derived in a more formal manner including a rigorous treatment of the previously unsolved issue of the spring leg dynamics during flight. This motion, which becomes singular for vanishing footmasses, was simply ignored in [7] and is treated here by the inclusion of additional holonomic constraints.
The model consists of a torso with mass which is constrained to purely linear motions as defined in [7]. Thus, the torso’s configuration is given by the hip position (). The leg is connected to the hip via a rotational joint (with joint angle ) that includes a torsional spring (with stiffness and no damping). We model the legs as massless linear springs with leg length , natural spring length , spring stiffness , no damping, and a point mass at the foot. The total mass of the model is . We use generalized coordinates (i.e., ) to represent the configuration of the robot.
The model has two hybrid modes: stance and flight . The corresponding constraint forces in these modes are and . These forces satisfy the constraints
(23)  
(24) 
during flight and stance, respectively. The constraint (23) fixes the leg length to during flight, whereas equation (24) implements the assumption of no sliding during stance (with a horizontal contact point position ). For the continuous dynamics in equations (1), we have
(25) 
where , , describe the hip and leg spring forces, respectively. Note, the tangential and normal contact forces , are only active during stance. Similar, only holds during flight to constraint the leg length to its natural length . This leads to impulsive forces and thus discontinuous changes in whenever the foot leaves the ground with nonzero velocity. The touchdown event is defined kinematically, while the liftoff event is triggered when changes sign from positive to negative. We restrict all motions to the mode sequence , which is started at apex transit during flight.
Eventually, we would like to bring the foot mass to zero to avoid kinetic energy losses during touchdown, similar to the method used in [7, 20, 8, 12]. To fulfill C2 and satisfy equation (16), we redefine the foot mass by . Note, for , the condition in equation (16), is equally satisfied for the stance and flight transition:
(26)  
(27) 
Further, to maintain finite continuous dynamics (1a) in the limit , we redefine the constraint forces as:
(28) 
with , introducing new auxiliary forces , . The core idea here is to separate the constraint forces into two components, where the first balances the elastic forces which are expressed by the known values of . The second component balances the inertial forces and is computed when solving the differentialalgebraic equations (1). This second component is further scaled with to yield finite values for and , even in the limit . Equivalently to [7], we prescribed a leg swing frequency by the relation
(29) 
This implies that remains a finite constant value when the foot mass is brought to zero and thus . With the modifications in equations (28), (29) and taking the limit , we arrive at the same finite dimensional dynamics reported in [7]. To allow for horizontal displacement in equation (19), the matrix selects the initial state . The remaining periodic states are selected by its orthogonal complement . In this energetically conservative model, all state and parameter values are normalized with respect to , and .To allow a comparison with [7], we set the leg stiffness to (which is equivalent to hopping with 2 legs of stiffness ) and the swing frequency to .
IvB Results
Using this model, Algorithm 2 was initialized with a vertical hopping motion at energy level (that is, with initial apex height of ). Here, the motion in and simply follows a parabolic trajectory during flight and a linear oscillation during stance. There is no movement in and . This hopping motion constitutes a regular point that solves equation (19). This initial point is connected to a locally defined 1D manifold (Fig. 3) of hopping in place motions. Towards lower energies, hopping height is reduced and this generator is bounded by a point that corresponds to a vanishing time in flight at energy level . Periodic solutions of with even lower energy do exist, yet they correspond to an oscillating inplace motion. Since there is no liftoff in this motion, going beyond this point leads to a change in contact sequence. This is an inadmissible point since it breaks Assumption A4. However, there exists a locally defined manifold with this different contact sequence ^{5}^{5}5Of interest is its connectedness to an equilibrium point (EQ) (Fig. 3). The contact sequence
admits solutions in the linear eigenspace of a 1D oscillator. These linear modes exist in the range
, where is the energy at EQ.. It can be independently computed by Algorithm 1, however, it is not in the connected component of generators with contact sequence .Carrying on with , we traverse the generator towards higher energies. is bounded by a simple bifurcation point at energy level . At this point, we find three nearby generators for which the last tangent direction of points into the new generator . consists of purely vertical hopping motions at higher energies than in . The remaining generators, and , consist of forward () and backward () hopping motions, respectively. The computation of leads to another simple bifurcation point at . We find three new connected generators –. The vertical motions in are similar to gaits in and . The generators and correspond to forward and backwards hopping motions, respectively.
Figure 4 illustrates three gaits from , and at energy level . The motions in and are qualitatively different in the leg’s angular velocity at touchdown. In , the foot touches down with while the gaits in possess a longer flight duration in which the foot touches down in a returning motion with (so called speed matching). This holds equivalently for backward hopping in and .
We stopped the exploration of – with regular points at .
It is possible to encounter more special points (BP, TP, IP) in the numerical continuation at higher energy levels.
It took approximately a minute on a laptop with an i58265U CPU @1.60GHz and 4GB RAM to generate the data^{6}^{6}6The code to generate this data can be found at https://github.com/
raffmax/ConnectingGaitsinEnergeticallyConservativeLeggedSystems presented in Figures 3 and 4.
V Discussion & Conclusion
In this paper, we introduced a formal framework and a generalized methodology for the computation of connected gaits in energetically conservative legged systems. This work extends and clarifies the methodology introduced in [7] to apply not only to the gaits of legged models but to a broader class of ECMs. In terms of theory, our work extends the results in [24] to hybrid dynamical systems and clarifies the connected structure of the gait space of energetically conservative legged systems.
Our contributions further relate the study of passive gaits to established and emerging concepts in the field of nonlinear dynamics. Similar to the generators in [1], we (locally) define 1D manifolds in which there is a unique relation between motion and energy. However, our definition of these generators is different in that these 1D manifolds do not include equilibria (A2) and they are defined for hybrid dynamical systems. As a consequence, the direct connection to linear oscillations, that occur in the linearized system at equilibrium and that is a characteristic of today’s NNMs definitions, is lost. This loss is caused by two required assumptions. The first is the transversality condition of the anchor constraint that is violated in an equilibrium. We introduced it here to impose a Poincaré section, yet it can potentially be lifted, as it is done in [23]. The second is the fixed contact sequence (A4) that prohibits the connection of an equilibrium at standstill to a forward gait.
As shown in Fig. 3, the linear modes of the 1D oscillator correspond to bouncing in place. In future work, it may be possible to formally link them to the hopping gaits characterized in this paper. To make this possible, we need to relax the requirement that the contact sequence is fixed (A4). This assumption constitutes the primary limitation of our work. It is necessary, as the core results in this paper follow from the monodromy matrix . For a fixed contact sequence, changes differentiably in neighboring periodic solutions and so does the associated tangent space. This is no longer true when Assumptions A3 (i.e. the transition event is grazed by the flow) or A4 do not hold. Without A3, a Saltation matrix [15] may become discontinuous [11], which directly propagates to discontinuities in and the associated tangent space. For legged systems, continuity in the Saltation matrix can be ensured under certain conditions [21]. Turning these conditions into systematic modeling guidelines, or finding ways to connect gaits despite these discontinuities are avenues for future work.
While the focus of this paper is on energetically conservative systems, real robot systems are not energetically conservative and sources of energy loss (heat, impacts, batteries, vibrations) cannot be completely eliminated. The benefit of our approach is in utilizing the explanatory power of ECMs. While these systems do not exist in the real world, these simple models often form the core model dynamics for trajectory generation, motion planning, and control algorithms in the field. Mapping trajectories from ECMs to more realistic models with energy loss would be an interesting extension of our work, as their passivity makes them ideal candidates for the use as templates to develop energetically economical motions for legged robotic systems.
Beyond this very practical significance, the identified passive motions are a key characteristic of a given ECM. Their study, not only in simple models of legged systems, will thus allow us to better understand the fundamental nature of gait for both, robotics and biology.
References
 [1] (2020) A review on nonlinear modes in conservative mechanical systems. Annual Reviews in Control. Cited by: §I, §IID, §V.
 [2] (2012) Numerical continuation methods: an introduction. Vol. 13, Springer Science & Business Media. Cited by: §IIIB, §IIIB, §IIIB, §IIIB.
 [3] (2001) Numerical continuation, and computation of normal forms. In In Handbook of dynamical systems III: Towards applications, Cited by: §IIIB.
 [4] (2014) Lowbandwidth reflexbased control for lower power walking: 65 km on a single battery charge. The International Journal of Robotics Research 33 (10), pp. 1305–1321. Cited by: §I.
 [5] (2016) Nonsmooth mechanics: models, dynamics and control. Springer. Cited by: §IIA, §IIA.
 [6] (200502) Efficient bipedal robots based on passivedynamic walkers. Science 307 (5712), pp. 1082–1085. Cited by: §I.
 [7] (2018) All common bipedal gaits emerge from a single passive model. Journal of The Royal Society Interface 15 (146), pp. 20180455. Cited by: §I, §I, §I, §IID, §IVA, §IVA, §IVA, §V.
 [8] (1998) The simplest walking model: stability, complexity, and scaling. Journal of biomechanical engineering 120 (2), pp. 281–288. Cited by: §I, §IVA.
 [9] (2018) Gait based on the springloaded inverted pendulum. In Humanoid Robotics: A Reference, A. Goswami and P. Vadakkepat (Eds.), pp. 1–25. Cited by: §I.
 [10] (2001) Asymptotically stable walking for biped robots: analysis via systems with impulse effects. IEEE Transactions on automatic control 46 (1), pp. 51–64. Cited by: item 1, §IIA.
 [11] (1998) The stability of periodic solutions of discontinuous systems that intersect several surfaces of discontinuity. Journal of Applied Mathematics and Mechanics 62 (5), pp. 677–685. Cited by: §IIC, §V.
 [12] (2016) A hybrid systems model for simple manipulation and selfmanipulation systems. The International Journal of Robotics Research 35 (11), pp. 1354–1392. Cited by: item 4, §IIA, §IIIA, Remark III.1, §IVA, footnote 1.
 [13] (2012) Capturabilitybased analysis and control of legged locomotion, part 1: theory and application to three simple gait models. The International Journal of Robotics Research 31 (9), pp. 1094–1113. Cited by: §I.
 [14] (2001) A simple model of bipedal walking predicts the preferred speed–step length relationship. Journal of Biomechanical Engineering 123 (3), pp. 264. Cited by: §I.
 [15] (2013) Dynamics and bifurcations of nonsmooth mechanical systems. Vol. 18, Springer Science & Business Media. Cited by: §IIC, Remark II.2, §V.
 [16] (1987) Groucho running. Journal of Applied Physiology 62 (6), pp. 2326–2337. Cited by: §I.
 [17] (2015) Stable running with asymmetric legs: a bifurcation approach. International Journal of Bifurcation and Chaos 25 (11), pp. 1550152. Cited by: §I.
 [18] (1995) Calculation of lyapunov exponents for dynamic systems with discontinuities. Chaos, Solitons & Fractals 5 (9), pp. 1671–1681. Cited by: §IIC.
 [19] (2003) Continuation of periodic orbits in conservative and hamiltonian systems. Physica D: Nonlinear Phenomena 181 (12), pp. 1–38. Cited by: footnote 2.
 [20] (2009) The relative roles of dynamics and control in bipedal locomotion. Ph.D. Thesis, University of Michigan. Cited by: §I, §IVA.
 [21] (2017) Decoupled limbs yield differentiable trajectory outcomes through intermittent contact in locomotion and manipulation. In 2017 IEEE International Conference on Robotics and Automation (ICRA), pp. 2261–2266. Cited by: §V.
 [22] (2009) Continuation of periodic solutions of dissipative and conservative systems: application to elastic pendulum. Mathematical Problems in Engineering 2009. Cited by: §IIIB.
 [23] (2021) A topological approach to gait generation for biped robots. IEEE Transactions on Robotics (), pp. 1–20. External Links: Document Cited by: §V, Definition.
 [24] (1997) Localized oscillations in conservative or dissipative networks of weakly coupled autonomous oscillators. Nonlinearity 10 (3), pp. 679. Cited by: §I, §IIC, Remark II.3, §IIIB, §V.