Exciting families of passive gaits in an elastic quadruped via natural motion manifold control

Abstract

Motivated by the need for efficiency and robustness in repetitive robotic tasks such as locomotion, this study introduces the concept of Natural Motion Manifolds (NMMs) and presents a control method to stabilize and excite motions based on these structures. By considering the intersection of a Poincaré section with a surface comprising a continuum of autonomous evolutions, the proposed controller extends the linearized Poincaré map control from a single orbit to a family of orbits. This allows us to derive simple controllers to excite intrinsic nonlinear resonances and exploit the natural dynamics when varying the energy target (or the running velocity). We validated the method through simulations and experiments on a serial elastic quadruped. Relying on natural dynamics and minimal motor commands, we could implement a bounding gait at desired velocities without needing dynamic compensations. The experiments provide a thorough validation of the feasibility and the benefits of controlled, predictable, and purposeful oscillatory behavior via explicit excitation of a quadruped’s natural dynamics.

1. Introduction

Numerous studies on biological systems have shown that repetitive and rhythmic motion, such as walking and running, is largely guided by the intrinsic mechanical response of the physical structure (Alexander, 2002; Ahlborn and Blake, 2002). Especially during fast locomotion, a carefully tuned interaction between the pull of gravity, the limbs’ inertia, and the legs’ visco-elastic response ¹ is advantageous to minimize energy expenditure, enhance stability, and reduce control complexity (Daley, 2018; Geyer et al., 2006; Ludeke and Iwasaki, 2020; Poulakakis et al., 2006). Similarly, taking full advantage of natural dynamics in robotics promises to yield increased performance, less expensive hardware, and simpler control principles (Kashiri et al., 2018). Since the pioneering works of Raibert (Raibert et al., 1986) and McGeer (Mcgeer, 1990), substantial strides have been made in developing robots that include passive compliant elements (Badri-Spröwitz et al., 2022; Calisti et al., 2017; Hutter et al., 2013; Zhang et al., 2020). Both theoretical and experimental evidence confirm that leveraging elasticity by storing and releasing energy results in increased efficiency and robustness (Alexander, 1990; Chen et al., 2019; Hurst and Rizzi, 2008; Vu et al., 2015).

However, practical exploitation of the elasticity and natural motion presents significant control challenges. Section 2 discusses research directed to addressing it. As of today, significant limitations still persist. With the exception of specific works related to Hybrid Zero Dynamics (HZD) (Da et al., 2016; Da and Grizzle, 2019), existing control strategies primarily focus on the control of a single gait trajectory (Cao et al., 2015). This strongly limits leveraging the natural dynamics for speed transitions, disturbance recovery, and the methodological development of effective and efficient controllers to fully harness nonlinear intrinsic resonances. On the other hand, methods based on the excitations of resonance limit the analysis to linear approximations. However, locomotion entails large-amplitude oscillations in multi-body systems, which, coupled with elasticity and segmented designs, generate substantial nonlinear behavior. In this context, sufficiently soft springs are necessary to store enough energy for fully passive locomotion; however, this form of compliant actuation is limited by low bandwidth, which poses challenges for active corrections during highly dynamic motions. Consequently, most techniques are unable to robustly handle the limited control authority inherent in serial soft elastic actuators (necessary for efficiency and versatility (Hutter et al., 2013)). Furthermore, the experimental validation of compliant robot designs and controllers exploiting natural dynamics to realize dynamic gaits remains limited.

We propose a new methodology combining geometry, nonlinear dynamics, and closed-loop control that simultaneously addresses all these open challenges. In doing so, this work aligns itself with a series of activities that we have been pursuing in the past years. Using control to exploit intrinsic oscillations in elastic systems, including robots, was explored in Stratmann et al. (2017) and Lakatos (2018) with a bang-bang controller. This simple and robust method showed promising results but relies on heuristics and empirical tuning. A systematic approach using nonlinear modal theory (Albu-Schäffer and Della Santina, 2020) achieved forward hopping in an elastic monopod (Calzolari et al., 2022) and efficient pronking gaits in a quadruped (Calzolari et al., 2023). However, the more we applied nonlinear modal theory, the more we realized that some constraints inherent in that framework were over-constraining the solution space. To harness natural dynamics effectively, a broader range of oscillation families is needed.

This work goes the extra step of describing, stabilizing, and exciting a whole family of periodic gaits exploiting natural dynamics. We introduce Natural Motion Manifolds, a further generalization of Eigenmanifolds. They are one-parameter families of conservative natural oscillations in general hybrid Lagrangian systems. At the same time, these geometrical objects generalize and provide extra structure to works looking into continuation algorithms (Raff et al., 2022a). This, in turn, allows us to borrow ideas that we developed for the Eigenmanifold case (Della Santina and Albu-Schaeffer, 2021) and combine them with a novel extension of the Poincaré linearization method. The result is a robust and effective control methodology that we can use to excite hyper-efficient ² natural evolutions in legged systems and other robots.

As a result, we can, for the first time, experimentally demonstrate the excitation of entire families of quasi-passive locomotion gaits on a full-fledged quadruped robot equipped with soft serial elastic actuators (see Figure 1).

Figure 1.

The figure illustrates the key concepts presented in this article. Simplified models of robots equipped with soft springs are used to identify Natural Motion Manifolds (NMM), which are geometrical structures collecting a family of open-loop periodic solutions of the nonlinear system. In the case of locomotion, these orbits correspond to passive gaits. The intersection of the NMM with a Poincaré section is used to construct a coordinate system decoupling transversal and energy dynamics, and to derive simple controllers to stabilize and excite the intrinsic nonlinear oscillations: higher energy levels are realized by “climbing up” the NMM. The experiments on a serial elastic quadruped achieved a bounding gait at two desired velocities (and transition between them) by completely relying on the natural dynamics and minimal motor commands.

To summarize, this work contributes to the state of the art in elastic locomotion with

(1) The definition of the Natural Motion Manifold (NMM), a geometric concept related to families of conservative periodic orbits.

(2) The empirical characterization of NMM for conservative and non-conservative mechanical systems.

(3) A generalization of the linearized Poincaré map control approach that can be used to stabilize orbits on an NMM.

(4) A method that leverages contributions (1) and (3) to identify, stabilize, and excite passive gaits in a series elastic quadruped. It requires no empirical tuning and allows continuous adjustment of energy and speed.

(5) The first demonstration of experimental excitation of families of locomotion gaits of the natural dynamics of the quadruped.

The remainder of the paper is organized as follows. Section 2 discusses the state of the art in locomotion via natural dynamics excitation. Section 3 introduces NMMs and extended Poincaré sections, and Section 4 a general control strategy to excite evolutions contained in NMMs—that is, contribution (1) to (3). Additionally, we provide a discussion on the analogies and differences with HZD-related works (Da and Grizzle, 2019) concerning the use of multiple periodic solutions to enhance locomotion control. We then narrow our focus to locomotion by introducing contribution (4) in Section 5. Finally, Section 6 reports the extensive experimental results in efficient locomotion of an elastic quadruped, contribution (5).

2. Related work

Locomotion is an inherently difficult task, requiring, in general, careful planning, powerful actuators, and many sensors. This task is even more challenging when the additional goal of exploiting the robot’s natural dynamics is added. Due to the extensive nature of the existing literature, we will focus on a selection of papers that present the most relevant concepts for our study, with particular emphasis on works that exploit compliance to enhance locomotion performance in legged systems.

2.1. Identification and control of passive gaits

A framework for identifying passive and actuated gaits was presented in Remy et al. (2011). This provided the foundations for works that subsequently demonstrated the theoretical existence and connectedness of natural gaits (Gan et al., 2018; Raff et al., 2022a), showed that selecting gaits based on speed leads to increased efficiency (Xi et al., 2016), and developed algorithms to find gaits and optimally deal with energy losses (Raff et al., 2022b).

Controllers that focus on the stabilization of passive gaits include methods based on foot placement and energy controllers (Ahmadi and Buehler, 2006), constant touchdown angle with stance-retraction (Poulakakis et al., 2005), time delay orbit controllers (Owaki et al., 2013) in combination with energy-preserving strategies (Hyon and Emura, 2004), applications of chaos controllers and Poincaré maps (Huang et al., 2017), and combinations with learning (Wu et al., 2023). Alternatively, there exist methods exploring the combination of passive dynamics with model-based control, such as Renjewski et al. (2015) and Beck et al. (2022). Additionally, Cao et al. (2015) and Cao and Poulakakis (2016) suggested switching between fixed-points to robustly transition to different gaits and speeds for simulated models. Different from these works, the approach proposed in this article explicitly uses the natural dynamics to design the controllers and excite nonlinear oscillations, sustaining and varying the speed of a gait in hardware.

2.2. Methods based on emergent behavior

Several approaches aim to harness natural dynamics by leveraging the emergent behavior induced by an external excitation. For example, by adding elasticity in the feet (Spröwitz et al., 2013) and in the spine (Wang et al., 2017), central pattern generators can be used to self-stabilize locomotion in quadrupeds (Ruppert and Badri-Spröwitz, 2022). Alternatively, sinusoidal open-loop excitations can be used (Nie et al., 2016). Efficient forward hopping of a curved beam is achieved in Reis and Iida (2014) using linear free vibrations. Bio-inspired reflexes (Masuda and Ishikawa, 2018) and negative hysteresis (Wanders et al., 2015) are also used to excite resonance. Vibrations are also used in crawling soft robots (Calisti et al., 2017). Most of the aforementioned methods are based on linear analyses to compute the natural frequencies, a perspective that remains actively supported (Kashiri et al., 2017; Ludeke and Iwasaki, 2020). Research on adaptive natural oscillators (Buchli et al., 2006; Khoramshahi et al., 2017) demonstrates the feasibility of exciting periodic orbits by tuning the frequency of sinusoidal oscillators to the resonant frequency of machines. However, convergence issues may arise when handling hybrid systems (Khoramshahi et al., 2017), and, importantly, nonlinearities significantly influence oscillations at large amplitudes. As a result, sine waves may be inadequate to excite locomotion machines that feature nonlinear behavior (Lakatos et al., 2013). Instead, a targeted impulse that accounts for the nonlinear orbit can be more successful, especially for line-shaped oscillations that cross an equilibrium (Lakatos et al., 2015). This research prompted the design and control of the compliant quadruped BERT (Lakatos et al., 2018), which uses a simple bang-bang controller to excite nonlinear oscillations by leveraging the direction of this lower-dimensional subspace. Yet, this latter approach focuses on orbits with line-shaped topology, which is quite restrictive compared to the oscillations considered here.

2.3. Control-oriented methods

A distinctly different class of approaches includes algorithms rooted in control-related principles (such as stability, dimensionality reduction, and input minimization) that progressively integrated elasticity into their frameworks for overall improvements. A prominent example is compliant HZD (Poulakakis, 2008), which exploits physical springs by aligning virtual constraints with the elastic response so as to encode templates such as SLIP as the zero dynamics of higher-dimensional segmented legged robots (Poulakakis and Grizzle, 2009). Some relevant examples can be found in Sreenath et al. (2013), Hubicki et al. (2018), Ma et al. (2017), and Reher and Ames (2021), where compliant elements in the legs are leveraged to achieve improved stability and higher speeds. Typically, virtual constraints are realized via feedback compensations that may interfere with the natural system response. This is mitigated by relaxing the constraints along certain subspaces to obtain a natural passive response (Sadeghian et al., 2017; Sreenath et al., 2011), or via robot design (Hubicki et al., 2018). However, enforcing reduced models to increase tractability still requires partial compensations. Additionally, in Rezazadeh et al. (2015), a fixed natural frequency calculated for the leg at rest length is used for step planning. Within the same approach, oscillations are forced via sinusoidal excitation, designed using linear arguments (Rezazadeh and Hurst, 2020).

A different idea is utilized within the Divergent Component of Motion (DCM) planning framework (Englsberger et al., 2015), where the timing between steps can be tuned to match the (linearized) intrinsic natural frequencies for compliant quadrupedal walking (Pollayil et al., 2022).

Model-based optimizations for planning and control can leverage intrinsic dynamics with energy economy as the main objective. Efficient locomotion for a soft, elastic biped is shown in Werner et al. (2017), and joint stiffness and trajectory tuning are optimized in Gasparri et al. (2018). Performance gains through a combination of optimization and iterative learning are explored in Ding et al. (2024) for a parallel elastic quadruped.

2.4. Current limitations

Broadly speaking, control-oriented solutions offer structured approaches, but the controllers they generate can sometimes be practically ineffective, particularly when a desired dynamic behavior is imposed on a compliant system (Keppler et al., 2018). These methods also face limitations inherent to their specific frameworks. For example, DCM requires an accurate, high-bandwidth torque controller, which is challenging to achieve with soft series elastic actuators (SEAs). Works related to HZD often rely on simulation-based searches (as other optimization methods) to generate feasible controls, which can suffer from issues like local minima, scaling in higher dimensions, and constraints that reduce versatility. Addressing these challenges is a focus of recent HZD research (Da and Grizzle, 2019; Ma et al., 2017). Another key consideration is that most existing control strategies focus on stabilizing a single orbit (Cao et al., 2015), which limits the ability to fully exploit the natural nonlinear dynamics for speed transitions and disturbance recovery. Notable exceptions can be found in specific HZD-based studies (Da et al., 2016; Da and Grizzle, 2019; Nguyen et al., 2017), where controlled limit cycles are interpolated based on locomotion velocity. This allows for adjusting the gait based on speed, resulting in more natural-looking gaits and improved stability. This feature is common to both of these works and the present study, albeit with notable differences discussed in Section 4.7. Overall, control methods adapted to natural compliance balance tractability with performance. However, they rely on partial model compensations that can conflict with the natural dynamics, and the oscillatory nature of the entire mechanical system is ultimately not explicitly addressed.

Emergent behavior approaches typically yield simpler strategies and are effective in generating oscillations in hardware. However, analyzing convergence to a desired motion is challenging, especially as the degrees of freedom (DoF) increase. Sinusoidal excitation based on linear analysis may be unsuitable to sustain large nonlinear oscillations (Lakatos et al., 2013).

In summary, a comprehensive framework is still needed—one that leverages natural dynamics without limiting robot design while ensuring predictability, stability, and versatility. Achieving this is challenging, particularly in bridging theoretical and practical aspects for general nonlinear cases.

3. Families of energy-conservative oscillations

In this work, we are interested in families of periodic orbits of mechanical systems that conserve energy (i.e., the natural oscillations). These trajectories encode efficient oscillations of the natural dynamics and offer a geometrical structure that can be exploited for planning and control (Albu-Schäffer and Della Santina, 2020; Della Santina and Albu-Schaeffer, 2021; Della Santina et al., 2021). We’d like to highlight that until Section 5, our discussion remains pertinent to a quite general class of mechanical systems, demonstrating applicability beyond locomotion, including potential use in other contexts, such as robot arms.

A periodic orbit is a trajectory that repeats itself exactly after a time interval T > 0. Consider a system within the class of hybrid Lagrangian systems (Burden et al., 2015; Colombo and Eyrea Irazú, 2020; Raff and Remy, 2024). Given its dynamical flow $\phi (\mathit{x},t):\mathcal{X}\times \mathbb{R}\to \mathcal{X}$ with state $\mathit{x}\in \mathcal{X}$ and time t, a periodic orbit $\mathcal{O}(\mathit{x},\phi )=\{\phi (\mathit{x},t)|t\in \mathbb{R}\}$ satisfies

$\phi (\mathit{x},t)=\phi (\mathit{x},t+T).$

We say that a periodic orbit $\mathcal{O}(\mathit{x},\phi )$ is conservative if there is a quantity E( x ) that is non-constant in an open set and it is conserved along the whole orbit, that is,

$\forall \mathit{x}\in \mathcal{O}(\mathit{x},\phi ),\frac{d}{dt}E(\mathit{x})=0.$

In this work, we focus on the energy as the conserved quantity. A family of periodic orbits can be defined as a set of orbits ${\mathcal{S}}_{\mathcal{O}}(\mathit{x},\phi ;E)$ , where each orbit has period T(E), and the family encompasses all possible orbits within the parameter space defined by E. Such a family defines a two-dimensional sub-manifold in $\mathcal{X}$ that we name “Natural Motion Manifolds.”

Definition 1

A two-dimensional invariant sub-manifold $\mathcal{M}\subset \mathcal{X}$ is a Natural Motion Manifold (NMM) if all evolutions on it are periodic and conservative.

Here, “natural motion” indicates that the orbits are autonomous, open-loop solutions of the nonlinear system. Prime examples of systems featuring such mathematical structures are smooth conservative systems, that is, systems which satisfy

$\dot{\mathit{x}}=\mathit{f}(\mathit{x}),\frac{\partial E(\mathit{x})}{\partial \mathit{x}}\mathit{f}(\mathit{x})=0.$

(1)

It is well-known that this class of systems cannot have isolated closed orbits, and they are locally part of one-parameter families of orbits when they possess a single first integral (Albu-Schäffer and Della Santina, 2020; Sepulchre and MacKay, 1997). This result was recently extended to hybrid dynamical systems (Colombo and Eyrea Irazú, 2020; Raff and Remy, 2024). In the present paper, we extend this concept by considering hybrid Lagrangian systems that may include energy losses yet exhibit a subspace of solutions in which the energy is conserved, allowing for the possibility of families of periodic solutions. Figure 2 presents examples of systems that feature such families of energy-conservative orbits.

Starting from an equilibrium or one periodic orbit, the family and its bifurcations can be evaluated using, for example, numerical continuation-based algorithms within a desired energy range (see Raff et al. (2022a)). In general, a system can feature a multitude of natural motion manifolds, which may be connected or not. Additionally, from the viewpoint of bifurcation theory (Strogatz, 2018), these manifolds may exist only within certain energy ranges (e.g., asymptotically reach a homoclinic orbit). An example of these manifolds can be found in Gan et al. (2018), where all kinds of gaits of a simple conservative biped model are connected. Another example of a family of energy-conservative periodic orbits is given by Eigenmanifolds (i.e., nonlinear modes) (Albu-Schäffer and Della Santina, 2020). Indeed, Eigenmanifolds can be seen as a sub-case of the more general NMM that is obtained when a stable equilibrium (from which the manifold spawns) is included in $\mathcal{M}$ , and they are always families of breaking trajectories. A breaking trajectory is line-shaped, with the kinetic energy and thus velocity being zero at the ends of the line. This strongly limits their application in the planning and control of locomotion gaits. ³ Neither of these two requirements is necessary for defining an NMM in Def. 1.

While the local existence of neighboring closed orbits is well-studied for conservative systems (Muñoz-Almaraz et al., 2003; Sepulchre and MacKay, 1997), the same cannot be said for non-globally conservative ones. In that case, energy fluctuations typically destroy periodic orbits or, in special cases, isolate an orbit and generate a limit cycle (Strogatz, 2018). However, we emphasize that, in principle, a system that is not globally energy-conservative can still possess families of energy-conservative orbits. This is why the orbits must be conservative in the definition of an NMM, but the system that generates them is not necessarily globally conservative. In fact, we discovered such a case for the model of the sagittal quadruped robot presented in Figure 2(c) and described in Section 5, which revealed to feature many of these families.

To narrow our focus on natural motion manifolds that admit an energy parametrization, we need to first summarize the concepts of Poincaré sections, Poincaré maps, and fixed points.

Figure 2.

Examples of systems that feature families of conservative periodic orbits: (a) a pendulum with radial spring, (b) a double pendulum (exhibiting bifurcation), and (c) a sagittal elastic quadruped with massless legs. All systems are affected by gravity.

3.1. Background: Poincaré section and maps

Given a state space of dimension n, the Poincaré section is an n − 1 dimensional sub-manifold typically used to study the properties of periodic orbits of nonlinear systems. When choosing a section for a particular orbit, the section can be freely defined along the orbit and transversal to the flow. The section is not necessarily orthogonal to the flow, but cannot be tangent to it.

Definition 2

The Poincaré forward map (or Poincaré map) x → P ( x ) takes a point x on the section Σ and returns the solution of the system at its first return to Σ. A discrete system can be defined by iterating the map

${\mathit{x}}_{\mathrm{k}+1}=\mathit{P}\left({\mathit{x}}_{\mathrm{k}}\right).$

(2)

A point x * is called a fixed-point of the system if it maps onto itself, that is, x * = P ( x *).

By analyzing the fixed points, one can gain insights into the system’s behavior (e.g., local stability) and can use it to improve control design (Dadashzadeh et al., 2014). Under additional technical conditions (Burden et al., 2015), Poincaré maps are well-defined for hybrid dynamical systems. When the system has inputs, P can be augmented to account for their influence between section hits. To this end, consider systems with an input (or a parameter) vector u . The map x → P ( x , u ), where P ( x , u ) is the point of first return to Σ of the solution with constant inputs from point x , is called the controlled Poincaré map (Fradkov and Evans, 2005), and the associated discrete system is

${\mathit{x}}_{\mathrm{k}+1}=\mathit{P}({\mathit{x}}_{\mathrm{k}},{\mathit{u}}_{\mathrm{k}}).$

(3)

3.2. Fixed-points generators

Let us assume the existence of a Poincaré section Σ that transversely intersects an NMM. This intersection, schematically presented in Figure 3, identifies a curve γ, which collects the fixed points ⁴ x *. If the manifold is smooth, γ is also smooth. If there exists a bijective mapping between the parameter domain and the points on the curve, γ can be parametrized by E

${\mathit{x}}^{*}=\mathit{P}({\mathit{x}}^{*})=\gamma (E({\mathit{x}}^{*})),$

(4)

that is, each orbit of the family is uniquely identified by the fixed points on the section using one parameter. This excludes manifolds with disjoint orbits sharing the same value E from our analysis.

Figure 3.

The intersection between a Poincaré section Σ and a manifold comprising energy-conservative autonomous evolutions identifies a smooth curve γ parametrized by energy E. This curve is utilized for control design to stabilize and excite efficient mechanical oscillations (such as locomotion gaits) at different energies.

Using the above concepts, we can define curves named “Fixed-points Generators” that completely identify an NMM.

Definition 3

Given an NMM $\mathcal{M}\subset \mathcal{X}$ and a Poincaré section Σ defined in the embedding space $\mathcal{X}$ , then we call the curve $\gamma :\mathbb{R}\to \mathrm{\Sigma }\cap \mathcal{X}$ parametrizing the intersection a Fixed-points Generator of $\mathcal{M}$ . We call E the parametrizing variable in γ(E).

In the following, we will take the energy E of the curve passing through the parametrized points as the parametrizing variable, and we center our efforts on mechanical systems with NMMs that admit these types of fixed-point generators. We assume that a desired NMM has been identified and focus on controlling these objects. The main idea is that natural motion manifolds underpin nonlinear resonance dynamics, which we can exploit to achieve efficient motion if properly excited.

4. Natural motion manifold control

We can expect that the periodic orbits described by a natural motion manifold are not attractive in general, and some form of stabilization is required to realize them in practice. This is certainly the case for purely conservative systems ⁵ (Remy, 2011). Here, we seek a control method for stabilizing a system on one of its NMMs and exploiting its structure to reach a desired energy level. A strategy to accomplish this goal, but limited to the stabilization of Eigenmanifolds, was recently proposed in Della Santina and Albu-Schaeffer (2021) using a continuous controller. This approach can be regarded as efficient for the excitation of oscillations in mechanical systems since the control effort on the manifold converges to zero. However, as mentioned in Section 3, Eigenmanfolds are limited to a certain orbit topology, which, in contrast to NMMs, limits their application in locomotion. Inspired by these methods, we propose to export some control concepts of nonlinear modal theory to the control of more general oscillations (such as locomotion gaits) and merge them with approaches for controlling chaos (Fradkov and Evans, 2005). To that end, we propose a control architecture to stabilize NMMs using the associated extended controlled Poincaré map, and consisting of the following two tasks.

(1) The stabilization task locally maintains the system close to the natural motion manifold by means of linearized Poincaré maps of the transversal dynamics. Stabilization to the manifold is essential to suppress the onset of irregular or multi-periodic oscillations. ⁶ This is different from stabilizing the fixed-point of one single orbit (OGY chaos controller (Fradkov and Evans, 2005; Huang et al., 2017)), or locally stabilizing the orbit using optimal time-varying controllers such as Beck and Sakamoto (2023), or using projection operators on the orbit in Saetre and Shiriaev (2020).

(2) The excitation task injects energy to sustain the oscillation at a desired energy level by exploiting the tangential dynamics on the manifold.

The two controllers are designed to avoid interfering with each other, with their control actions converging to zero upon reaching the target orbit on the NMM (hyper-efficiency). This scheme can be seen as an extension of the linearized Poincaré map controller from a point to a curve, combined with an open-loop excitation of unstable periodic orbits. This approach increases the linear controller’s attraction region and further reduces the stabilization control space by 1 dimension. For the derivation of the controllers, a definition of transverse coordinates for the Poincaré map adapted to a curve is introduced, which facilitates the design of both the stabilization and the excitation controllers.

4.1. Transverse coordinates system

We intend to use the fixed-point generator γ to derive the stabilization and excitation controllers. To this end, we require new coordinates to construct a control error variable more easily. A natural choice is to use the energy parametrization, which uniquely defines a point on γ and avoids utilizing non-physically consistent norms to define a distance to the curve. Thus an error δ x _k for the k‐th hit at the Poincaré section can be defined as

$\delta {\mathit{x}}_{\mathrm{k}}={\mathit{x}}_{\mathrm{k}}-\gamma (E({\mathit{x}}_{\mathrm{k}})).$

(5)

The transverse coordinates are defined by introducing a local diffeomorphism in a neighborhood of the curve γ, as shown in Figure 4. A point $\mathit{x}\in {\mathbb{R}}^{n_{\mathrm{x}}}$ expressed in the new coordinates $\mathit{z}=(\epsilon ,\mathit{y})\in {\mathbb{R}}^{n_{\mathrm{x}}}$ is

$\mathit{x}(\epsilon ,\mathit{y})=\gamma (\epsilon )+{\mathrm{\Gamma }}_{\perp }^{\text{e}}(\epsilon )\mathit{y},$

(6)

where $\epsilon \in \mathbb{R}$ parameterizes γ and its value is the energy ɛ = E( x ), ɛ ∈ [E_min, E_max], while $\mathit{y}\in {\mathbb{R}}^{n_{\text{y}}}$ are the transversal coordinates (and we have momentarily dropped the iteration subscript k). The matrix ${\mathrm{\Gamma }}_{\perp }^{\text{e}}(\epsilon )$ is a basis for the transversal space, and it is defined to be orthogonal to the variation of energy with respect to the state, that is,

$\frac{\partial E(\mathit{x})}{\partial \mathit{x}}{\mathrm{\Gamma }}_{\perp }^{\text{e}}(E(\mathit{x}))=\mathbf{0}.$

(7)

Here, ${\mathrm{\Gamma }}_{\perp }^{\text{e}}(\epsilon )$ is a valid basis for the transverse coordinate system because, by construction, $\frac{\partial E(\mathit{x})}{\partial \mathit{x}}$ and ${\gamma }^{\prime }=\frac{\partial \gamma (\epsilon )}{\partial \epsilon }$ cannot ever be orthogonal to each other (since γ is parametrized using E). The associated map defined by the Jacobian is

$\delta \mathit{x}=\left[\begin{array}{cc}\hfill {\gamma }^{\prime }(\epsilon )+\frac{\partial {\mathrm{\Gamma }}_{\perp }^{\text{e}}(\epsilon )}{\partial \epsilon }\mathit{y},\hfill & \hfill {\mathrm{\Gamma }}_{\perp }^{\text{e}}(\epsilon )\hfill \end{array}\right]\delta \mathit{z}$

(8)

where the Jacobian is always full rank along γ. Additionally, we would like to remark that transverse coordinates are sufficient for stabilization, and strict orthogonality of the coordinate system is not required in general (Manchester, 2011; Saetre and Shiriaev, 2020). The basis ${\mathrm{\Gamma }}_{\perp }^{\text{e}}(\epsilon )$ is used instead of an orthogonal component to γ′, because it allows to locally decouple ɛ from y along γ, since the evolution along y (locally to γ) takes place on a constant energy sub-manifold $\{\mathit{x}\in \mathcal{X}:E(\mathit{x})=\epsilon \}$ . From (6), the inverse transformation can be derived

$\left[\begin{array}{c}\hfill \epsilon \hfill \\ \hfill \mathit{y}\hfill \end{array}\right]=\left[\begin{array}{c}\hfill E(\mathit{x})\hfill \\ \hfill {\mathrm{\Gamma }}_{\perp }^{\text{e}}{(E(\mathit{x}))}^{\top }(\mathit{x}-\gamma (E(\mathit{x})))\hfill \end{array}\right].$

(9)

Using this transformation, we can now write the Poincaré map dynamics (3) in the new coordinates z _k = (ɛ_k, y _k). For a compact derivation, we use (3) and (6) to define the map

$g({\mathit{z}}_{\mathrm{k}},{\mathit{u}}_{\mathrm{k}})={\mathit{x}}_{\mathrm{k}+\mathrm{1}}=\mathit{P}\left(\gamma ({\epsilon }_{\mathrm{k}})+{\mathrm{\Gamma }}_{\perp }^{\text{e}}({\epsilon }_{\mathrm{k}}){\mathit{y}}_{\mathrm{k}},{\mathit{u}}_{\mathrm{k}}\right).$

(10)

Then, using (9), the update equation for ɛ is

${\epsilon }_{\mathrm{k}+\mathrm{1}}=E\left(g({\mathit{z}}_{\mathrm{k}},{\mathit{u}}_{\mathrm{k}})\right),$

(11)

and the transversal error dynamics y become

$\begin{array}{cc}\hfill {\mathit{y}}_{\mathrm{k}+\mathrm{1}}=& {\mathrm{\Gamma }}_{\perp }^{\text{e}}{(E({\mathit{x}}_{\mathrm{k}+\mathrm{1}}))}^{\top }({\mathit{x}}_{\mathrm{k}+\mathrm{1}}-\gamma (E({\mathit{x}}_{\mathrm{k}+\mathrm{1}}))\hfill \\ \hfill {\mathit{y}}_{\mathrm{k}+\mathrm{1}}=& {\mathrm{\Gamma }}_{\perp }^{\text{e}}{(E(g({\mathit{z}}_{\mathrm{k}},{\mathit{u}}_{\mathrm{k}})))}^{\top }\left(g({\mathit{z}}_{\mathrm{k}},{\mathit{u}}_{\mathrm{k}})-\gamma (E(g({\mathit{z}}_{\mathrm{k}},{\mathit{u}}_{\mathrm{k}})))\right)\hfill \end{array}$

(12)

Figure 4.

Transverse coordinates system defined on the curve γ that collects the fixed points.

4.2. Transverse linearization

The Poincaré map dynamics is linearized to derive the local controllers in the new coordinates. Let us write the forward map function (12) of for the transverse coordinate y _k concisely as

$\mathit{y}=f_{\mathrm{y}}(\epsilon ,\mathit{y},\mathit{u})={\mathrm{\Gamma }}_{\perp }^{\text{e}}{\left(E\left(g\right)\right)}^{\top }(g-\gamma \left(E\left(g\right)\right).$

The first order Taylor expansion of f_y at y = 0, u = 0 at the energy level ɛ = ɛ_k along γ yields

${\mathit{y}}_{\mathrm{k}+1}=f_{\mathrm{y}}({\epsilon }_{\mathrm{k}},\mathbf{0},\mathbf{0})+\frac{\partial f_{\mathrm{y}}}{\partial \mathit{y}}{\mathit{y}}_{\mathrm{k}}+\frac{\partial f_{\mathrm{y}}}{\partial \mathit{u}}{\mathit{u}}_{\mathrm{k}}+O({\mathit{y}}^2,{\mathit{u}}^2),$

${\mathit{y}}_{\mathrm{k}+1}\approx {\mathit{A}}_{\mathrm{y}}\left({\epsilon }_{\mathrm{k}}\right){\mathit{y}}_{\mathrm{k}}+{\mathit{B}}_{\mathrm{y}}\left({\epsilon }_{\mathrm{k}}\right){\mathit{u}}_{\mathrm{k}}\mathbf{,}$

(13)

where we leveraged the fact that, for a conservative orbit, ɛ is constant, and for all orbits on γ, f_y(ɛ_k, 0, 0) = 0.

Thanks to the choice of transverse coordinates, one can show that the dynamics of ɛ is linearly decoupled from y , it remains constant on γ, and it only depends on the inputs. The linearization of the ɛ dynamics from (11) along γ, gives the approximation

${\epsilon }_{\mathrm{k}+\mathrm{1}}\approx {\epsilon }_{\mathrm{k}}+\frac{\partial E(g)}{\partial \mathit{y}}\delta {\mathit{y}}_{\mathrm{k}}+\frac{\partial E(g)}{\partial \mathit{u}}{\mathit{u}}_{\mathrm{k}}.$

(14)

Using the coordinate basis ${\mathrm{\Gamma }}_{\perp }^{\text{e}}$ , it follows that

$\begin{array}{cc}\hfill \frac{\partial E(g)}{\partial \mathit{y}}& =\frac{\partial E}{\partial \mathit{x}}\frac{\partial P}{\partial \mathit{x}}\frac{\partial (\gamma (\epsilon )+{\mathrm{\Gamma }}_{\perp }^{\text{e}}(\epsilon )\mathit{y})}{\partial \mathit{y}}\hfill \\ \hfill & =\frac{\partial E}{\partial \mathit{x}}\frac{\partial P}{\partial \mathit{x}}{\mathrm{\Gamma }}_{\perp }^{\text{e}}(\epsilon )=\frac{\partial E}{\partial \mathit{x}}{\mathrm{\Gamma }}_{\perp }^{\text{e}}(\epsilon )=\mathbf{0},\hfill \end{array}$

(15)

where we have used the property that for a differentiable function f and a fixed point x it holds

$f(\mathit{x})=f(\mathit{P}(\mathit{x}))\Rightarrow \frac{\partial f}{\partial \mathit{x}}=\frac{\partial f}{\partial \mathit{x}}\frac{\partial \mathit{P}(\mathit{x})}{\partial \mathit{x}}.$

(16)

Therefore, the linearized update equation for the energy can be written as

${\epsilon }_{\mathrm{k}+\mathrm{1}}={\epsilon }_{\mathrm{k}}+{\mathit{B}}_{\epsilon }({\epsilon }_{\mathrm{k}}){\mathit{u}}_{\mathrm{k}}.$

(17)

4.3. Input realization

To allow linearly decoupled control of ɛ and y , we introduce the following assumption.

Assumption 1

The input vector u _k has dimension $\ge n_x$ , and the matrix $\mathit{D}={[{\mathit{B}}_{\epsilon }^{\top },{\mathit{B}}_{\mathrm{y}}^{\top }]}^{\top }$ has full row rank.

This can be achieved by defining the input forces to the system $\dot{\mathit{x}}(t)=\mathit{f}(\mathit{x}(t),\mathit{u}(t))$ between each section hit via a realization

$\mathit{u}(t)={\psi }_{\mathrm{\Sigma }}(\mathit{x}(t),{\mathit{u}}_{\mathrm{k}}),$

(18)

where u _k represents a new control input, which is held constant over the kth cycle. As an example, u _k can encode the inputs for each flow phase of a hybrid dynamical system. In general, ψ_Σ has to be designed on the basis of the specific characteristics of the NMM to fulfill Assumption 1.

The linearized dynamics (11) and (12) are now utilized to obtain control gains. The input u _k is separated between a stabilization part u _s,k and an energy injection part u _j,k

${\mathit{u}}_{\mathrm{k}}({\epsilon }_{\mathrm{k}},{\mathit{y}}_{\mathrm{k}})={\mathit{u}}_{\text{s,k}}({\epsilon }_{\mathrm{k}},{\mathit{y}}_{\mathrm{k}})+{\mathit{u}}_{\text{j,k}}({\epsilon }_{\mathrm{k}}),$

(19)

where we have made explicit the dependence on the states when this vector is updated upon a section hit.

4.4. Stabilization

This task aims to stabilize the transversal y dynamics (12) without interfering with the energy dynamics (11). Therefore, we introduce a new input ${\mathit{u}}_{s,k}^{\prime }$ which acts in the nullspace of B _ɛ (ɛ_k) using the projection

${\mathit{u}}_{\text{s,k}}={\mathit{N}}_{\epsilon }(\epsilon ){\mathit{u}}_{s,k}^{\prime },$

(20)

where N _ɛ is a nullspace basis of B _ɛ (ɛ_k) such that

$\frac{\partial E}{\partial {\mathit{u}}_{\text{s,k}}^{\prime }}={\mathit{B}}_{\epsilon }{\mathit{N}}_{\epsilon }=\mathbf{0},$

(21)

and the pair ( A _y(ɛ), B _y(ɛ) N _ɛ (ɛ)) is controllable. Consequently, the local energy variation (17) is decoupled from ${\mathit{u}}_{\text{s,k}}^{\prime }$ , and for u _j,k = 0 it results

${\epsilon }_{\mathrm{k}+\mathrm{1}}={\epsilon }_{\mathrm{k}}.$

(22)

Thus, we can regard ɛ as the operating point for the stabilization controller ${\mathit{u}}_{\text{s,k}}^{\prime }$ . Using the system matrix A _y(ɛ) and input matrix B _y(ɛ) N _ɛ (ɛ), a control gain K _y(ɛ) is designed to stabilize the decoupled transversal dynamics

${\mathit{u}}_{\text{s,k}}^{\prime }({\epsilon }_{\mathrm{k}},{\mathit{y}}_{\mathrm{k}})=-{\mathit{K}}_{\mathrm{y}}({\epsilon }_{\mathrm{k}}){\mathit{y}}_{\mathrm{k}}.$

(23)

4.4.1. Stability for the conservative case

Since any input (23) does not affect the energy, we can analyze the local stability of each fixed point at γ(ɛ) independently. Therefore, if the condition on magnitude of the largest eigenvalue $\overline{\lambda }(\cdot )$

$\left|\overline{\lambda }\left({\mathit{A}}_{\mathrm{y}}(\epsilon )-{\mathit{B}}_{\mathrm{y}}(\epsilon ){\mathit{K}}_{\mathrm{y}}(\epsilon )\right)\right|<1,\forall \epsilon \in [E_{\text{min}},E_{\text{max}}]$

(24)

is satisfied, that is, A _y(ɛ) − B _y(ɛ) K _y(ɛ) is Schur stable, the transversal error is locally asymptotically stabilized within the studied energy range. Asymptotic orbital stability of the original system on $\mathcal{M}$ follows (Fradkov and Evans, 2005).

4.5. Energy controller

The energy controller u _j,k is designed to excite the oscillations at a desired energy without pushing the system away from the manifold, that is, without interfering with the transverse dynamics. Under the assumption that the stabilizer (23) keeps the system close to the manifold, we can then “climb up” the manifold to reach a desired energy target.

We propose to use the following proportional-integral control with scalar gains β and ρ

${\mathit{u}}_{\text{j,k}}=\frac{{\mathit{N}}_{\text{y}}\left({\epsilon }_{\mathrm{k}}\right)}{{||{\mathit{B}}_{\epsilon }\left({\epsilon }_{\mathrm{k}}\right){\mathit{N}}_{\text{y}}\left({\epsilon }_{\mathrm{k}}\right)||}_2}\left(\beta (E_{\text{d}}-{\epsilon }_{\mathrm{k}})+\rho \sum _{\mathrm{j}=0}^{\mathrm{k}}{E}_{\text{d}}-{\epsilon }_{\mathrm{j}}\right),$

(25)

where ${‖\cdot ‖}_2$ is the matrix 2-norm, E_d. is the desired energy level, and N _y is a nullspace basis of B _y, chosen such that B _ɛ (ɛ_k) N _y(ɛ_k) > 0, that is, we locally impose an energy variation without influencing the transversal y dynamics. In the energy-conservative case, the proportional part is sufficient to converge to a desired energy target. Practically, however, the integral control part is necessary for compensation when energy losses occur along the nominal orbit (e.g., due to friction). This controller works reasonably well under the assumption that the losses are small and do not quickly push the system away from the manifold.

Remark 1

If the nullspace of either B _y(ɛ) or B _ɛ (ɛ) contains more than one vector, the respective nullspace bases N _y(ɛ) and N _ɛ (ɛ) should be chosen consistently to yield smooth continuous profiles with respect to energy variations.

Remark 2

Under Assumption 1, one choice for the decoupling matrices can be

${\mathit{N}}_{\epsilon }={\mathit{D}}_W^{+}\left[\begin{array}{c}\hfill \mathbf{0}\hfill \\ \hfill {\mathit{I}}_{n_{\text{y}}}\hfill \end{array}\right],{\mathit{N}}_{\mathrm{y}}={\mathit{D}}_W^{+}\left[\begin{array}{c}\hfill 1\hfill \\ \hfill {\mathbf{0}}_{n_{\text{y}}\times 1}\hfill \end{array}\right],$

where ${(\cdot )}_W^{+}$ is the generalized pseudoinverse weighted with some matrix W.

4.5.1. Stability of the energy varying case

By combining the stabilizer and injection, and canceling out the zero terms, the linearized dynamics reduces to

$\begin{array}{cc}\hfill {\epsilon }_{\mathrm{k}+\mathrm{1}}& ={\epsilon }_{\mathrm{k}}+\alpha (E_{\text{d}}-{\epsilon }_{\mathrm{k}})+\beta \sum _{j=0}^k{E}_{\text{d}}-{\epsilon }_{\text{j}},\hfill \\ \hfill {\mathit{y}}_{\mathrm{k}+\mathrm{1}}& =\left({\mathit{A}}_{\mathrm{y}}({\epsilon }_{\mathrm{k}})-{\mathit{B}}_{\mathrm{y}}({\epsilon }_{\mathrm{k}}){\mathit{K}}_{\mathrm{y}}({\epsilon }_{\mathrm{k}})\right){\mathit{y}}_{\mathrm{k}}.\hfill \end{array}$

(26)

Since ɛ is decoupled from y , its convergence to ɛ = E_d can be inferred. A conditional stability argument can be used to conclude the asymptotic stability of the transversal dynamics y by recalling (24).

To prove local exponential stability of (26), the y dynamics may be viewed as a linear parameter-varying system, where the controller is gain-scheduled with the energy as the operation point. Therefore, the related techniques can be applied to quantitatively characterize its stability properties (Shamma, 2012).

4.6. Example: Pendulum with radial spring

As a first, simple example, we present the application of the proposed approach for the stabilization and excitation of a pendulum with a radial spring along one of its natural motion manifolds. This system is schematically presented in Figure 5, and it has state $\mathit{X}={[\theta ,r,\dot{\theta },\dot{r}]}^{\top }$ and parameters m = 1 kg, k = 50 Nm/rad and radial rest length r₀ = 0.5 m. It has two inputs: a torque τ acting in the radial direction θ and a linear force f acting along r. Starting from the lower equilibrium configuration, a natural motion manifold is obtained via energy-based continuation (Figure 7) up to an energy level of 12 J (and the energy at the equilibrium is zero). In this case, the chosen manifold corresponds to an Eigenmanifold (Albu-Schäffer and Della Santina, 2020), that is, a nonlinear mode that spawns from an eigenvector of the equilibrium. We choose a Poincaré section $\mathrm{\Sigma }=\{(\theta ,r):\theta =0,\dot{\theta }>0\}$ . Then, a fixed point on the section is defined by $\mathit{x}=[r,\dot{\theta },\dot{r}]$ . Both the rotational and radial dynamics are assumed to be actuated. Based on the oscillation trajectories and the manifold geometry (Figures 5 and 7), $\dot{\theta }$ is chosen as switching variable for defining a convenient new input vector $\mathit{u}\in {\mathbb{R}}^4={[u_1,u_2,u_3,u_4]}^{\top }$ , which maps to the original system’s inputs with the following realization:

$\left[\begin{array}{c}\hfill \tau \hfill \\ \hfill f\hfill \end{array}\right]=\left[\begin{array}{c}\hfill s(u_1,u_2,\dot{\theta })\hfill \\ \hfill s(u_3,u_4,\dot{\theta })\hfill \end{array}\right].$

Here, the function s(u_i, u_j, x) = u_i − (u_i − u_j)σ(cx), with c > 0 and σ(⋅) being the sigmoid function, is used to smooth the transition between u_i and u_j when the argument x crosses zero. In this way, the controlled Poincaré map x _k+1 = P ( x _k, u _k) has sufficiently large input space to allow independent control of the energy and the transversal dynamics. We proceed to derive the controllers described in Section 4.4 and Section 4.5 for the entire energy range. Since the fixed point is three-dimensional, the new transversal coordinates y are two-dimensional. The stabilizing feedback gain K _y(ɛ) is designed via a discrete infinite horizon linear quadratic regulator (DLQR) with Q = I ₂ and R = 10⁻² I ₄. The energy regulator has gains β = 0.15 and ρ = 0. The matrix elements of the resulting gain $\mathit{K}(\epsilon )={\mathit{K}}_{\mathrm{y}}(\epsilon ){\mathrm{\Gamma }}_{\perp }^{\text{e}}{(\epsilon )}^{\top }$ , obtained by combining (23) and (9), as well as the energy injection profiles, are shown in Figure 6. To validate the controllers, the system is simulated with the initial condition X ₀ = [0, 0.65, 2, 0.25] which lies outside of the manifold, and the desired energy level is set to 8 J. Figure 7 reports the resulting trajectories and input profiles. The transversal error is quickly stabilized, the target energy level is reached while traveling on the NMM, and the inputs converge to zero.

Figure 5.

Left Plot: Coordinates and input forces for a pendulum with radial spring. Right Plot: Natural motion manifold of the radial mass-spring system projected on the xz plane. The color-coded lines represent regular periodic oscillations for a range of energy levels. In this case, the manifold coincides with an Eigenmanifold.

Figure 6.

Elements of the stabilization gain K (ɛ) and of the energy injection gain L (ɛ) (i.e., the normalization factor in (25)) for the pendulum with radial spring as a function of energy. The subscripts (⋅) _τ and (⋅) _f denote the elements of the matrix row corresponding to corrections to torque (via u₁, u₂) and linear force (via u₃, u₄), respectively.

Figure 7.

System trajectories (color-coded) of the radial mass-spring system converging to a natural oscillation at a target energy level on the identified natural motion manifold (gray surface). The red line represents the intersection of the manifold with the Poincaré section, that is, the curve γ. The color of the simulated trajectories for validation transitions from blue to yellow to represent time. After a brief transient, the trajectory is stabilized on the manifold. When the orbit reaches the target energy, the control effort converges to zero.

4.7. Relations to previous studies

As mentioned in Section 2, the idea to utilize a multitude of periodic solutions to improve stability is shared with the line of works related to Da and Grizzle (2019), and to some degree also with Cao and Poulakakis (2016), though there are several distinctions. In Da et al. (2016), the velocity-based interpolation of controllers is driven by the intuition that this enhances stability, whereas here a continuum of controllers emerges from the fact that NMMs are smooth collections of infinite solutions with different energy levels. Because of this, interpolation or switching between fixed points is not necessary, and actually it is not possible: each point in γ is assumed to represent a passive periodic orbit. Discretization is used solely for practical implementations. Another key difference is that Da and Grizzle (2019) seeks stabilization of an orbit via open-loop trajectories (not necessarily periodic or passive), while in this work, a family of passive orbits is utilized as a starting point to build controllers. Additionally, natural families are identified via numerical continuation, in contrast to optimization in Da et al. (2016), Da and Grizzle (2019), and Cao and Poulakakis (2016). Finally, the energy controller is developed using the manifold, and it is decoupled from the stabilizer; these concepts are not found in the previously mentioned works.

5. Application to locomotion

We present a methodology to make use of the NMM control on a family of serial-elastic quadrupeds featuring lightweight, ⁷ non-redundant legs. ⁸ This consists of the following three steps: (i) the derivation of a control-oriented reduced order model suitable for exploring the natural dynamics, (ii) the identification and parametrization of the natural motion manifold, and (iii) the design of the stabilization and energy controllers presented in Section 4. We validate the approach through simulations and experiments.

5.1. Modeling

In this section, we introduce the model for a family of quadruped robots and provide a reduced-order representation of their dynamics.

5.1.1. Full-order model

Consider a generic quadrupedal robot constrained to the sagittal plane as shown in Figure 8. The robot consists of a floating base main body, or trunk, with legs attached to it and equipped with serial elastic actuators (SEA) at each joint. Let $\mathit{\zeta }={[{{\mathit{x}}_{\mathrm{b}}}^{\top },{\mathit{q}}^{\top },{\mathit{\theta }}^{\top }]}^{\top }\in {\mathbb{R}}^{3+2n_{\mathrm{j}}}$ , where ${\mathit{x}}_{\mathrm{b}}={\left[\begin{array}{ccc}\hfill x,\hfill & \hfill z,\hfill & \hfill \varphi \hfill \end{array}\right]}^{\top }\in {\mathbb{R}}^3$ are the trunk’s CoM position and pitch angle, $\mathit{q}\in {\mathbb{R}}^{n_{\text{j}}}$ denotes the joint angles of the legs, while $\mathit{\theta }\in {\mathbb{R}}^{n_{\text{j}}}$ is the motor position of each actuator. Using the Lagrangian formalism, the full-order sagittal dynamics is described by (Pollayil et al., 2022)

$\begin{array}{ll}\hfill \left[\begin{array}{cc}\hfill {\mathit{M}}_{\text{b}}\hfill & \hfill {\mathit{M}}_{\text{bq}}\hfill \\ \hfill {\mathit{M}}_{\text{bq}}^{\top }\hfill & \hfill {\mathit{M}}_{\text{q}}\hfill \end{array}\right]\ddot{\mathit{\zeta }}+\mathit{c}(\mathit{\zeta },\dot{\mathit{\zeta }})+\mathit{g}(\mathit{\zeta })& =\left[\begin{array}{c}\hfill \mathbf{0}\hfill \\ \hfill {\mathit{\tau }}_{\mathrm{k}}\hfill \end{array}\right]+{\mathit{J}}_{\mathrm{c}}{(\mathit{\zeta })}^{\top }\mathit{\lambda }\hfill \end{array}$

(27)

$\begin{array}{ll}\hfill {\mathit{B}}_{\text{m}}\ddot{\mathit{\theta }}+{\mathit{\tau }}_{\mathrm{k}}& ={\mathit{\tau }}_{\text{m}},\hfill \end{array}$

(28)

where M _b = diag(m, m, I_b) is the inertia matrix of the trunk, ${\mathit{M}}_{\text{q}}\in {\mathbb{R}}^{n_{\mathrm{j}}\times n_{\mathrm{j}}}$ the inertia of the legs, ${\mathit{M}}_{\text{bq}}\in {\mathbb{R}}^{3\times n_{\mathrm{j}}}$ the inertia couplings, $\mathit{c}\in {\mathbb{R}}^{n_{\mathrm{j}}}$ collects the Coriolis and centrifugal forces, $\mathit{g}\in {\mathbb{R}}^{nj}$ the gravity vector, J _c is the contact Jacobian, and λ are the contact forces. The torque generated by a deflection of the springs is τ _k = K(θ − q), with $\mathit{K}\in {\mathbb{R}}^{n_{\mathrm{j}}\times n_{\mathrm{j}}}$ being the spring stiffness matrix. Additionally, ${\mathit{B}}_{\text{m}}\in {\mathbb{R}}^{n_{\mathrm{j}}\times n_{\mathrm{j}}}$ is the motor inertia matrix, and ${\mathit{\tau }}_m\in {\mathbb{R}}^{n_{\mathrm{j}}}$ represents the motor actuation forces. Since each joint is equipped with its own SEA, K and B are diagonal matrices. The potential energy in the springs is

$U_{\mathrm{k}}(\mathit{\theta },\mathit{q})=\frac{1}{2}{(\mathit{\theta }-\mathit{q})}^{\top }\mathit{K}(\mathit{\theta }-\mathit{q}).$

(29)

We assume that the feet are the only contact points, and no slipping or sliding can occur. When a foot strikes the ground, there is an instantaneous change in velocity

${\dot{\mathit{\zeta }}}^{+}=({\mathit{I}}_{{\text{n}}_{\mathrm{j}}}-{\mathit{M}}^{-1}{\mathit{J}}_{\mathrm{c}}^{\top }{({\mathit{J}}_{\mathrm{c}}{\mathit{M}}^{-1}{\mathit{J}}_{\mathrm{c}}^{\top })}^{-1}){\dot{\mathit{\zeta }}}^{-},$

(30)

where I _j is the j × j identity matrix, M is the inertia matrix of (27), and (⋅)⁻ and (⋅)⁺ indicate the velocities at the instants before and after the impact, respectively. If a leg is on the ground, the contact forces λ are calculated through the unilateral constraints and can be explicitly expressed by

$\mathit{\lambda }=-{({\mathit{J}}_{\mathrm{c}}{\mathit{M}}^{-1}{\mathit{J}}_{\mathrm{c}}^{\top })}^{-1}\left({\mathit{J}}_{\mathrm{c}}{\mathit{M}}^{-1}\left(\mathit{c}+\mathit{g}+\left[\begin{array}{c}\hfill \mathbf{0}\hfill \\ \hfill -{\mathit{\tau }}_{\mathrm{k}}\hfill \end{array}\right]\right)+{\dot{\mathit{J}}}_{\text{c}}\dot{\mathit{\zeta }}\right).$

(31)

The joint coordinates are arranged into $\mathit{q}={[{\mathit{q}}_{\text{st}}^{\top },{\mathit{q}}_{\text{fl}}^{\top }]}^{\top }$ , where q _st denotes the joints of the legs in contact with the ground, while q _fl are the joint angles of the legs not in contact. Following the same logic, the motor coordinates are organized as $\mathit{\theta }={[{\mathit{\theta }}_{\text{st}}^{\top },{\mathit{\theta }}_{\text{fl}}^{\top }]}^{\top }$ .

Figure 8.

Schematic representation of a serially elastic actuated robotic quadruped. Springs of stiffness k are collocated between a motor ( θ _i coordinate) and its joint (q_i coordinate). Springs play a crucial role in enhancing robustness and reducing mechanical stress during uncertain ground contacts and impacts. In this work, we are especially interested in them as potential energy storage units and how this capability can be exploited to improve locomotion efficiency.

5.1.2. Control-oriented reduced-order model

For the rest of the article, we focus on a simplified model that captures the key dynamics of interest. This model is defined under the following conditions:

C1. The legs have negligible mass and inertia, that is, M _q → 0, M _bq → 0.

C2. There are three possible locomotion phases: double support (one hind foot and one front foot in contact), single support (one foot in contact), and flight (no ground contact).

C3. When in contact with the ground, the inverse kinematics of the trunk coordinates is solvable in closed-form, that is, $\exists f_{\text{ik}}:{\mathbb{R}}^3\to {\mathbb{R}}^{n_{\text{j}}},suchthat{\mathit{q}}_{\text{st}}=f_{\text{ik}}({\mathit{x}}_{\mathrm{b}},{\mathit{x}}_{\mathrm{f}})$ .

C4. The motors are controlled in position, with target θ ^ref, and the effect of τ _k is negligible on the motor side dynamics, for example, due to high gear ratio and high enough controller gains.

C1 implies that there is no loss of kinetic energy at impacts, thus (30) simplifies to ${\dot{\mathit{\zeta }}}^{+}={\dot{\mathit{\zeta }}}^{-}$ . Moreover, $\mathit{c}(\mathit{\zeta },\dot{\mathit{\zeta }})\to 0$ . Due to C4 the motors are decoupled from the joint dynamics, and with proper tuning of the motor controllers (e.g., PID tuning), we further approximate the motor dynamics with a linear first order system which reaches its target position with settling time T _θ . Substituting (31) in the trunk dynamics of (27), and considering C1, C3, and C4, leads to the simplified sagittal 3 DoF dynamics with the decoupled motor dynamics

$\begin{array}{ll}{\mathit{M}}_{\text{b}}{\ddot{\mathit{x}}}_{\mathrm{b}}+\mathit{g}\left(z\right)& ={\frac{\partial f_{\text{ik}}\left({\mathit{x}}_{\mathrm{b}}\right)}{\partial {\mathit{x}}_{\mathrm{b}}}}^{\top }{\mathit{K}}_{\mathrm{st}}\left({\mathit{\theta }}_{\text{st}}-{\mathit{q}}_{\text{st}}\right)\end{array}$

(32)

$\begin{array}{ll}\hfill \dot{\mathit{\theta }}& =\frac{4}{T_{\mathit{\theta }}}\left({\mathit{\theta }}^{\text{ref}}-\mathit{\theta }\right).\hfill \end{array}$

(33)

This reduced-order model is mathematically equivalent to a free-floating trunk affected by the active springs’ projected forces, as shown in Figure 8.

5.1.3. Hybrid dynamics

When a foot hits or leaves the ground, the dynamics equations of the reduced-order model suddenly change. Therefore, this model is hybrid in nature, consisting of the three continuous-time dynamic models defined in C2. The phases and the transition guards are reported in Figure 9. The ground reaction forces at the hind and front feet are $F^{\text{f}}\in {\mathbb{R}}^4$ . The vertical components are denoted by $F_{\text{z}}^{\text{f}}$ . Since the legs are massless, these forces can be directly computed via the springs’ forces. The touchdown of the i-th foot is defined by the condition $z_f^i=0$ , where $z_f^i$ is the distance between the foot and the ground. Take-off occurs when the vertical ground reaction becomes positive, that is, $F_{z,i}^f>0$ . A vector ${\mathit{x}}_{\mathrm{f}}=({\mathit{x}}_{\mathrm{f}}^h,{\mathit{x}}_{\mathrm{f}}^f)\in {\mathbb{R}}^4$ registers the position of feet when they strike the ground and is updated upon each foot impact and detachment. Using x _f, one can reconstruct q via inverse kinematic, that is, q = h ( x _b, x _f) (C3). The Jacobian of the inverse kinematic function follows:

$\dot{\mathit{q}}=\mathit{H}({\mathit{x}}_{\mathrm{b}},{\mathit{x}}_{\mathrm{f}}){\dot{\mathit{x}}}_{\mathrm{b}},\mathit{H}({\mathit{x}}_{\mathrm{b}})=\frac{\partial \mathit{h}({\mathit{x}}_{\mathrm{b}},{\mathit{x}}_{\mathrm{f}})}{\partial {\mathit{x}}_b}.$

(34)

Figure 9.

The three possible flow states of the hybrid system modeling the dynamics of an elastic quadruped: double support, single support, and full flight. The transitions and guard conditions are also reported.

5.1.3.1. Double support phase

During double support, both feet are in contact. The trunk dynamics can then be expressed in trunk coordinates using (34)

${\mathit{M}}_{\text{b}}{\ddot{\mathit{x}}}_{\mathrm{b}}+\mathit{g}(z)=\mathit{H}{({\mathit{x}}_{\mathrm{b}},{\mathit{x}}_{\mathrm{f}})}^{\top }\mathit{K}\left(\mathit{\theta }-\mathit{h}({\mathit{x}}_{\mathrm{b}},{\mathit{x}}_{\mathrm{f}})\right).$

(35)

5.1.3.2. Single support phase

When only one virtual foot is in contact, the robot becomes an open serial kinematic chain. In this case, it is still possible to reconstruct the joint angles q _st of the leg in contact with the ground via an inverse kinematics function ${\mathit{q}}_{\mathrm{s}\mathrm{t}}={\mathit{h}}_{\mathrm{s}\mathrm{t}}({\mathit{x}}_{\mathrm{b}},{\mathit{x}}_{\mathrm{f}}^{\mathrm{s}\mathrm{s}})$ , where ${\mathit{x}}_{\mathrm{f}}^{\mathrm{s}\mathrm{s}}$ is the position of the foot in contact with the ground. The trunk’s dynamics becomes

${\mathit{M}}_{\text{b}}{\ddot{\mathit{x}}}_{\mathrm{b}}+\mathit{g}(z)={\mathit{H}}_{\mathrm{s}\mathrm{t}}{({\mathit{x}}_{\mathrm{b}},{\mathit{x}}_{\mathrm{f}}^{\mathrm{s}\mathrm{s}})}^{\top }\mathit{K}\left({\mathit{\theta }}_{\mathrm{s}\mathrm{t}}-{\mathit{h}}_{\mathrm{s}\mathrm{t}}({\mathit{x}}_{\mathrm{b}},{\mathit{x}}_{\mathrm{f}}^{\mathrm{s}\mathrm{s}})\right),$

(36)

where H _st is the Jacobian of h _st.

5.1.3.3. Flight phase

When in flight, the trunk follows a ballistic trajectory and conserves its momentum

${\mathit{M}}_{\text{b}}{\ddot{\mathit{x}}}_{\mathrm{b}}+\mathit{g}(z)=\mathit{0},\mathit{q}=\mathit{\theta }={\mathit{q}}_{\mathrm{a}\mathrm{t}\mathrm{k}}.$

(37)

Since the legs have no mass, we assume that the joint angles can be controlled during the flight phase to the desired attack configuration q _atk without affecting the trunk dynamics. Therefore, q _atk is considered as an input for this phase. Moreover, for the sake of clarity and intuition, q _atk is expressed using polar coordinates α = ( α _atk, r _atk), which comprises the legs attack angles ${\mathit{\alpha }}_{\mathrm{a}\mathrm{t}\mathrm{k}}\in {\mathbb{R}}^2$ as well as the attack radial lengths ${\mathit{r}}_{\mathrm{a}\mathrm{t}\mathrm{k}}\in {\mathbb{R}}^2$ . Finally, since we assume that the settling time of the motors is significantly shorter than the flight phase, we can ignore the motor dynamics during this phase and assume that the motors instantly reach their targets. As a consequence, U_k is zero during the flight phase.

5.1.4. Loss of global energy conservativeness

The total energy of the reduced-order model is

$E=\frac{1}{2}{\dot{\mathit{x}}}_{\mathrm{b}}^{\top }\mathit{M}{\dot{\mathit{x}}}_{\mathrm{b}}+U_{\mathrm{k}}(\mathit{\theta },{\mathit{x}}_{\mathrm{b}},{\mathit{x}}_{\mathrm{f}})+U_{\text{g}}({\mathit{x}}_{\mathrm{b}}).$

(38)

An important observation is that energy losses can occur in this model despite the legs having no mass. This can happen when a foot detaches from the ground without unloading all the springs of the corresponding leg (i.e., U_k ≠ 0 at take-off). Consider the elastic potential energy stored in the i-th leg ${U_{\mathrm{k}}}_{\text{i}}=1/2{\sum }_{j\in n_{\mathrm{i}}}{k}_{\mathrm{i}}{({\theta }_{\mathrm{j}}-q_{\mathrm{j}})}^2$ . We denote this energy with ${U_{\mathrm{k}}}_{\mathrm{i}}^{-}$ for the instant before take-off and with ${U_{\mathrm{k}}}_{\mathrm{i}}^{+}$ for the instant after. Since the legs of the reduced model are massless, at take-off (TO) it holds ${U_{\mathrm{k}}}_{\mathrm{i}}^{+}=0$ , and the total energy variation is

$E_{\mathrm{TO}}^{+}=E_{\mathrm{TO}}^{-}-U_{\mathrm{k}}^{-}.$

(39)

To avoid energy loss, ${U_{\mathrm{k}}}_{\mathrm{i}}^{-}$ must also be zero, implying that the leg should take off with unloaded springs.

The loss of energy conservativeness of the reduced-order model (32) is not present in the complete model (27) and (28). In theory, for the full model, the energy stored in its elastic elements is transferred to the leg at the moment of take-off, and it could be used to swing the leg forward. However, from a more practical viewpoint, the energy left in the springs gets quickly dissipated in high-frequency vibrations of the leg due to viscous friction in the joints because the inertia of the leg is very small. Since we want to actively avoid this scenario, we implicitly model this energy-loss feature when looking for energy-conservative gaits to minimize both these losses and undesired leg vibrations.

Other works such as Gan et al. (2018) search for passive gaits of purely conservative models, which typically necessitate mathematical tricks to model the leg dynamics and avoid losses (such as infinitely small foot masses and hip springs). Here, we still seek solutions in the spirit of natural dynamics exploitation, but we constructed a reduced model that, while still being abstract, is closer to the physical mechanism. In particular, by removing the dynamics of legs in the air, the search for gaits does not rely on oscillations that can only be very hard to realize in practice. The price to pay for these non-conservative models is that the manner in which the energy can be lost limits the set of energy-conservative solutions. Since U_k is a quadratic, positive definite function, it follows that ∂U_k/∂x|_{x = γ(E)} = 0. Therefore, the energy loss is linearly independent of any infinitesimal perturbation of the fixed point. In this sense, we can treat the system as “locally linearly conservative.”

5.1.5. Definition of the locomotion control inputs

For the planning and control of gaits, we introduce a flight phase input $\mathit{\alpha }\in {\mathbb{R}}^4$ and a stance phase input $\mathit{\eta }\in {\mathbb{R}}^4$ . Figure 10 schematically presents these high-level inputs over one stride. The flight input α = ( α _atk, r _atk) is introduced in Section 5.1.3 to define the desired attack configuration q _atk = q _atk( α ). The stance input η is a motor position reference offset which shifts θ ^ref by η in (33). The vector ${[\mathrm{\Delta }{\alpha }_{\text{h}},\mathrm{\Delta }{\alpha }_{\text{f}},\mathrm{\Delta }{r}_{\text{h}},\mathrm{\Delta }{r}_{\text{f}}]}^{\top }$ denotes the polar coordinates representation of η .

Figure 10.

Locomotion phases showing inputs for each phase: the legs are regulated to the attack configuration α during the flight phase, while a motor offset η is added during double support. These control inputs are held until the robot detects a full flight phase again.

5.2. Energy-conservative gaits (or passive gaits)

In the case of locomotion, the natural dynamics originate from the interaction of masses, springs, and gravity. To exploit the natural dynamics, we focus our attention on the periodic solutions of (27) in the energy-conservative case that does not require control effort, that is, we seek gaits where the motor positions θ are held constant during the whole stride. In this way, the natural elastic response is exploited, and the injection or removal of energy via the actuators is avoided. First, by using the state of the trunk and the attack configuration (the nominal motor positions) as optimization variables, an initial gait is found. Then, an energy-based continuation method is used to find neighboring conservative gaits and construct a family of orbits within an energy range. This family identifies a natural motion manifold of the system. The manifold underpins the nonlinear resonance associated with this family of gaits, which we aim to exploit to excite and sustain a gait at different speeds.

5.2.1. Definition of the Poincaré section

Similarly to other works (Gan et al., 2018; Remy, 2011), we define a convenient section Σ during the flight phase as

$\mathrm{\Sigma }=\{({\mathit{x}}_{\mathrm{b}},{\dot{\mathit{x}}}_{\mathrm{b}})|\dot{z}=0\},$

(40)

that is, the upper apex. Strictly speaking, locomotion gaits are not closed orbits of (27) or (32), since the forward position x is not periodic (Raff et al., 2022a). However, this is not an issue since, for our interests, we can simply avoid imposing periodicity along this dimension. As a result, only the following state variables are required to define the fixed point of a gait (27)

${\mathit{x}}_{\mathrm{k}}={\left[\begin{array}{cccc}\hfill z,\hfill & \hfill \varphi ,\hfill & \hfill \dot{x},\hfill & \hfill \dot{\varphi }\hfill \end{array}\right]}^{\top }.$

(41)

By means of the Poincaré maps, we introduce the definition of a locomotion gait.

Definition 4

A gait is a periodic orbit identified by a fixed point of P .

For locomotion, we are interested in gaits that fully exploit the natural dynamics, that is, open-loop unforced autonomous solutions that eliminate the need for unnecessary compensations of the dynamics. Mathematically, using the input definition in Section 5.1.5, we are looking for fixed points of the kind

${\mathit{x}}^{*}=\mathit{P}({\mathit{x}}^{*},{\mathit{\alpha }}^{*},{\mathit{\eta }}^{*}),$

(42)

where η * = 0 to avoid injecting or removing energy. The orbits corresponding to these fixed points are designated passive gaits. For the rest of this subsection, we drop the dependency on η .

5.2.2. Generation of a family of passive gaits

We now look for passive gaits of the quadruped model by finding the initial condition x * and attack configuration α * that satisfy periodicity

${\mathit{x}}^{*}=\mathit{P}({\mathit{x}}^{*},{\mathit{\alpha }}^{*}),$

(43)

Using the kinematic and dynamics parameters reported in the experimental validation Section 6.2.2, we fix α = α * and evaluate the natural motion manifold via energy-based numerical continuation, that is, we find a family of fixed points parametrized by energy x *(E) such that

${\mathit{x}}^{*}(E)=\mathit{P}({\mathit{x}}^{*}(E),{\mathit{\alpha }}^{*}).$

The results of the continuation are reported in Figure 11, which depicts the family of fixed points parametrized by energy. Additionally, Figures 12 and 13 show state space projections of the NMM and snapshots of a passive gait, respectively. Details regarding the algorithm to find this manifold are presented in Appendix B.

Figure 11.

Results of the energy continuation algorithm. Upper plot: The curves represent the fixed points belonging to a family of passive gaits parametrized by energy. Lower plot: Average forward speed and period T of the oscillations.

Figure 12.

State space projections of the natural motion manifold of the elastic quadruped found by energy continuation. The manifold is depicted with color-coded orbits of increasing energy. Each orbit represents a passive gait.

Figure 13.

Example of a passive gait belonging to a natural motion manifold, with snapshots of the trunk showing the motion from full flight (dotted black line) to double-support (solid black line) and vice versa. Single support phases are depicted with solid red lines.

Remark 3

Different from the setups proposed in Gan et al. (2018) and Raff et al. (2022a), the attack configuration is part of the optimization as a parameter (as opposed to the initial state of the legs). When looking for passive gaits, a typical approach consists of optimizing the initial condition (the state) (Gan et al., 2018), and relying completely on the passive dynamics to find a gait. However, practical experimentation with lightweight legs (see, e.g., Calzolari et al. (2022)), which is a fundamental design feature for efficient quadrupedal robots, has shown that relying on the passive swing dynamics of the leg during its aerial phase is not a robust approach, for example, any model mismatch and dissipation will heavily affect the swinging motion. As a consequence, the stability of the gait is very fragile. A more reliable approach consists of choosing a fixed attack configuration, which is held during the flight phase until touchdown. This allows us to look for more robust solutions and, at the same time, removes the necessity of modeling leg dynamics.

5.3. Controller design

Passive gaits studied via conservative systems are generally unstable (Iqbal et al., 2014; Remy, 2011), with limited stability within specific forward velocity ranges as observed in an analysis of an elastic quadruped (Poulakakis et al., 2006). Additionally, energy is typically lost during the motion due to various sources of dissipation, such as impacts and friction in the robot joints and between the feet and the ground, which further complicates the control of these solutions. Therefore, there is a need for both a stabilization and an energy re-injection mechanism to physically realize the passive gaits and exploit the natural resonance behavior.

5.3.1. Control derivation

We apply the control method presented in Section 4, and we utilize the flight phase for stabilizing the gait, while the stance phase serves for energy injection. The reason behind this choice is explained in the following. For analysis and experiments, we focus our attention on a particular class of quadrupedal robots equipped with soft springs in series with the actuators. A consequence of this setup is that the torque control bandwidth is very limited compared to the duration of the stance phases. This inherently hinders the reliable stabilization of orbits during stance phases via feedback control. Consequently, it is more straightforward to rely on the regulation of the attack configuration during the flight phase (which can be achieved more precisely), while the stance phase focuses on re-injecting energy.

Given the inputs α and η , and following Section 4.2, the linearized energy-transversal dynamics along γ at α = α *, y = 0 and η = 0 results in

$\begin{array}{cccccc}\hfill {\epsilon }_{\mathrm{k}+\mathrm{1}}& ={\epsilon }_{\mathrm{k}}\hfill & \hfill +& {\mathit{B}}_{\epsilon }^{\alpha }({\epsilon }_{\mathrm{k}})\mathrm{\Delta }{\mathit{\alpha }}_{\mathrm{k}}\hfill & \hfill +& {\mathit{B}}_{\epsilon }^{\eta }({\epsilon }_{\mathrm{k}}){\mathit{\eta }}_{\mathrm{k}}\hfill \\ \hfill {\mathit{y}}_{\mathrm{k}+\mathrm{1}}& ={\mathit{A}}_{\mathrm{y}}({\epsilon }_{\mathrm{k}}){\mathit{y}}_{\mathrm{k}}\hfill & \hfill +& {\mathit{B}}_{\text{y}}^{\alpha }({\epsilon }_{\mathrm{k}})\mathrm{\Delta }{\mathit{\alpha }}_{\mathrm{k}}\hfill & \hfill +& {\mathit{B}}_{\text{y}}^{\eta }({\epsilon }_{\mathrm{k}}){\mathit{\eta }}_{\mathrm{k}},\hfill \end{array}$

(44)

where Δ α _k = α _k − α *. Notice that the motor settling time T _θ to converge to the new desired position with offset η may take a substantial portion of the time of the stance phase, and therefore, the value of T _θ influences the overall dynamics. Hence, ${\mathit{B}}_{\epsilon }^{\eta }$ and ${\mathit{B}}_{\text{y}}^{\eta }$ will vary for different values of T _θ . During the flight phase, this is not relevant because the motors are not under load, the convergence to the desired attack configuration is quick, and the trunk dynamics are not influenced.

5.3.2. Stabilizer

Stabilization of the gait manifold is achieved via corrections of the attack configuration during the flight phase, that is, the controller utilizes Δ α as input, while the motor offset η is not used. Since the energy loss is quadratic (see Section 5.1.4), a variation of the attack configuration cannot affect the energy in the first-order approximation, that is,

${\mathit{B}}_{\epsilon }^{\alpha }=\frac{\partial E}{\partial \mathit{\alpha }}=\mathit{0}.$

Therefore, by using Δ α _k as control input for the stabilization, that is, we define u _s,k = Δ α _k, we don’t need to find the matrix N _ɛ introduced in (20) to decouple the energy dynamics. Given the system matrix A _y(ɛ_k) and input ${\mathit{B}}_{\text{y}}^{\alpha }({\epsilon }_{\mathrm{k}})$ , a DLQR is used to design the feedback

$\mathrm{\Delta }{\mathit{\alpha }}_{\mathrm{k}}=-{\mathit{K}}_{\mathrm{y}}({\epsilon }_{\mathrm{k}}){\mathit{y}}_{\mathrm{k}},$

(45)

using the cost function

$J=\sum _{\mathrm{k}=0}^{\infty }{\mathit{y}}_{\mathrm{k}}^{\top }{\mathit{Q}}_{\text{y}}{\mathit{y}}_{\mathrm{k}}^{\top }+\mathrm{\Delta }{\mathit{\alpha }}_{\mathrm{k}}^{\top }{\mathit{R}}_{\mathrm{\alpha }}\mathrm{\Delta }{\mathit{\alpha }}_{\mathrm{k}},{\epsilon }_{\mathrm{k}}=\epsilon ,$

where the state and input weighting matrices are

${\mathit{Q}}_{\text{y}}={\mathit{I}}_3,{\mathit{R}}_{\mathrm{\alpha }}=\text{diag}(1{\mathrm{e}}^1,1{\mathrm{e}}^1,1{\mathrm{e}}^3,1{\mathrm{e}}^3),$

(46)

such that the input is minimized and the radial offset is penalized in favor of relatively bigger corrections in the attack angle.

5.3.3. Energy controller

During the double support phase, the system is fully actuated, and we can reliably inject energy by commanding the motors to quickly insert potential energy into the springs. The offset η could be triggered at any instant during the stance phase or upon certain conditions. Additionally, it is not strictly required that all legs trigger the injection at the same time. The input matrices ${\mathit{B}}_{\epsilon }^{\eta }$ and ${\mathit{B}}_{\text{y}}^{\eta }$ will map these different implementations accordingly. We have, therefore, a design space for how η injects energy. In this work, we propose to trigger each leg asynchronously whenever its radial spring compression reaches a maximum. This choice is motivated by the fact that at the maximum compression, the injection of energy into the joints is more efficient, as explained in the following. For a serially elastic actuated robot joint with dynamics $m\ddot{q}=k(\theta -q)$ , being θ the motor coordinate, q the joint coordinate, k the spring stiffness, and m the joint inertia, the time derivative of the total mechanical energy in the joint $E=1/2m{\dot{q}}^2+1/2k{(\theta -q)}^2$ is

$\dot{E}=\dot{q}m\ddot{q}+(\theta -q)k(\dot{\theta }-\dot{q})=k(\theta -q)\dot{\theta }.$

(47)

If we consider an instantaneous motor position change Δθ, we can approximate the energy variation as

$\mathrm{\Delta }E=k(\theta -q)\mathrm{\Delta }\theta .$

(48)

With a properly oriented Δθ, the energy increases. Given a fixed value for Δθ, this increase is maximal when also τ_k = k(θ − q) is maximal, that is, when the spring is under maximum compression. ⁹ We apply this reasoning to the offset η _k: at the maximum leg compression, a given ΔE can be achieved with minimal η _k. Therefore, being ${\tau }_{\mathrm{i}}^{\text{r}}$ the radial spring force between the i-th shoulder and its foot, we trigger η _k when ${\tau }_{\mathrm{i}}^{\text{r}}>0$ and ${\dot{\tau }}_{\mathrm{i}}^{\text{r}}$ approaches zero, that is, the leg is at maximum radial compression. Practically, the command must be triggered when ${\dot{\tau }}_{\mathrm{i}}^{\text{r}}>-c_{\tau }$ , where c _τ is a positive constant, to account for the settling time of the motor (we want the motor to settle when the leg approximately hits the maximum compression). An approximate value for this constant can be estimated by looking at the trajectories and oscillation periods of the passive gaits, and considering the time constant of the motor.

Finally, the energy controller is constructed by treating η _k as the only input available to inject energy (i.e., u _j,k = η _k), and following (25) and Remark 2

${\mathit{\eta }}_{\mathrm{k}}={\left[\begin{array}{c}\hfill {\mathit{B}}_{\epsilon }^{\eta }({\epsilon }_{\mathrm{k}})\hfill \\ \hfill {\mathit{B}}_{\text{y}}^{\eta }({\epsilon }_{\mathrm{k}})\hfill \end{array}\right]}^{-1}\left[\begin{array}{c}\hfill 1\hfill \\ \hfill {\mathbf{0}}_{3\times 1}\hfill \end{array}\right]\left(\beta (E_{\text{d}}-{\epsilon }_{\mathrm{k}})+\rho \sum _{j=0}^k{E}_{\text{d}}-{\epsilon }_{\mathrm{j}}\right).$

(49)

Remark 4

Globally, modifying the attack configuration during flight (or the motor offset during stance) can lead to energy loss (Section 5.1.4). By avoiding this condition in both planning and control, that is, by staying close to the manifold, we not only avoid wasting energy, but we also avoid in-flight leg vibrations, which can deteriorate performance.

5.3.4. Practical implementation of the controllers

For the implementation of (45) and (49), the gain K _y(ɛ) and the matrices ${\mathit{B}}_{\text{y}}^{\eta }(\epsilon )$ and ${\mathit{B}}_{\text{y}}^{\alpha }(\epsilon )$ are evaluated at discrete points on γ, and then the matrices are interpolated, in a gain scheduling fashion to operate at the point ɛ. Additionally, since (45) could result in large values for Δ α _k, the input is scaled as follows:

$\mathrm{\Delta }{\mathit{\alpha }}_{\text{k,sat}}({\epsilon }_{\mathrm{k}},{\mathit{y}}_{\mathrm{k}})=c(\mathrm{\Delta }{\mathit{\alpha }}_{\mathrm{k}}({\epsilon }_{\mathrm{k}},{\mathit{y}}_{\mathrm{k}}))\cdot \mathrm{\Delta }{\mathit{\alpha }}_{\mathrm{k}}({\epsilon }_{\mathrm{k}},{\mathit{y}}_{\mathrm{k}}),$

(50)

where $c(\cdot ):{\mathbb{R}}^4\to \mathbb{R}$ is such that each element of Δ α _k remains within certain pre-determined bounds. This scaling is performed in an effort to satisfy the input bounds while preserving the input vector direction.

Figure 14 reports the elements of the stabilization gain $\mathit{K}(\epsilon )={\mathit{K}}_{\mathrm{y}}(\epsilon ){\mathrm{\Gamma }}_{\perp }^{\text{e}}{(\epsilon )}^{\top }$ , derived by combining (45) and (9), as a function of the energy. Figure 15 presents the designed polar and radial offsets used to inject 1 J into the system.

Figure 14.

Elements of the control gain matrix $\mathit{K}(\epsilon )={\mathit{K}}_{\mathrm{y}}(\epsilon ){\mathrm{\Gamma }}_{\perp }^{\text{e}}{(\epsilon )}^{\top }$ , parametrized by energy, that implements the controller defined in (45). K_α,i denotes the elements of the matrix row corresponding to the attack angle correction of the i-th leg, while K_r,i corresponds to the leg length correction.

Figure 15.

Motor offset η (here reported in polar and radial components) required during the stance phase to increase the energy of the oscillation by 1 J. The offset is found through (44) and (49). The direction of this vector is a function of the energy: this is required to allow energy injection along the curve γ without leaving it.

5.4. Simulations

We perform simulations of the reduced order model (32) under the controllers (45) and (49). We present three different cases to validate and illustrate the features of the controllers.

5.4.1. Stabilization

The system is initialized outside of the gait manifold on an initial state that is not a fixed point. The task of controller (45) is to stabilize the system to γ. Figure 16 presents the simulation results. In about three steps, the DLQR stabilizes the transversal coordinates to zero. The robot converges to a gait on the manifold and remains on it.

Figure 16.

Evolution of the system towards the gait manifold from an initial condition outside of it. The DLQR stabilizes the transversal coordinates in a few steps, and the system convergences to γ.

5.4.2. Reaching an energy target

In this instance, the robot is initialized on the manifold at E = 5.1 J, and we provide a desired energy target E_d = 6.5 J to be reached. Here, (45) and (25) work together to climb up the manifold and reliably reach the gait at a higher energy (and thus forward velocity). For the energy injection, we set gains β = 1, ρ = 0, and we saturate the amount of injected energy per stride to 0.1 J (to showcase the energy-climb up of the NMM). Figure 17 presents the results of the simulation, consisting of snapshots of the robot, the evolution along the manifold, and the evolution of the energy and the transversal coordinates. The targeted energy level is attained with minimal perturbation of the transversal dynamics.

Figure 17.

Projections of the gait manifold and evolution of the energy and transversal coordinates for the convergence to a target energy level.

5.4.3. Case with friction losses

Finally, we present the case in which we include joint friction losses. Friction forces are modeled as viscous friction $d_v\dot{\mathit{q}}$ with d _v = 1 × 10⁻² N s rad⁻¹. The system is initialized outside of the manifold, and the target energy is E_d = 6.5 J. Figure 18 presents the results of the simulation. Differently from the two previous test cases, energy losses affect the trajectories corresponding to the nominal gaits. As such, friction losses become a perturbation that, as expected, prevents the transversal dynamics from being stabilized to zero, and the integral action in the energy injection controller (25) is required (β = 1, ρ = 0.2). The maximum injected energy per stride is limited to 0.5 J. Despite the disturbance, the system is stabilized close to the manifold, and the energy level is reached.

Figure 18.

Stabilization and excitation of the gait manifold to reach an energy target under joint friction energy losses. In this case, an integral action is required to converge and sustain the gait at the desired energy. Due to the effect of friction (which acts as a disturbance), the transversal dynamics stabilize in a neighborhood of the manifold.

6. Experimental validation

6.1. BERT: A low-budget soft serial elastic quadruped

The quadruped robot BERT (Seidel et al., 2020) used for the experiments is presented in Figure 1. The robot features 2 DoF per leg, that is, there is no out-of-plane shoulder joint, therefore, the robot is mainly limited to sagittal motions. Each joint is equipped with a serial elastic actuator, which features a steel spring. The actuators are simple brushed DC servomotors that are controlled in position mode via an internal PID. Each SEA features two encoders: one for the motor output shaft θ_i and one for its joint position q_i. The leg kinematics, SEA designs, and experimental setup are presented in Figure 19. Additionally, the robot is equipped with an IMU, which provides the absolute orientation of the trunk, rotational velocity, and linear accelerations of the trunk’s COM. The kinematic and dynamic parameters are presented in Table 1.

Figure 19.

Upper Left: Side view of the BERT robot with its kinematics superimposed. Upper Right: Leg prototype implementing serial elastic actuators. A belt drive system decouples knee and hip joint motion (Lakatos et al., 2018). Bottom: Experimental setup consisting of a treadmill and a marker-based tracking camera system.

Table 1.

Robot parameters.

m b  (kg)	I b  (kg m2)	k (Nm rad−1)	h z  (cm)	a (cm)	l (cm)
3.0	0.08	3.0	2.0	16.5	8.0

6.2. Experimental setup

Experiments are carried out on a treadmill, and data is recorded using a tracking camera system using markers on the robot. The data tracking runs at 60 Hz, and it is used to regulate the treadmill’s speed to match the robot’s forward velocity, but it is not used for feedback in the robot controller. A steering controller is implemented to prevent the robot from falling off the treadmill due to its lack of shoulders. This controller induces a slight roll rotation of the robot’s trunk ( $\le 3$ degrees). This rotation around the x-axis corresponds to an asymmetry between the left and right legs ( $\le 5$ degrees), effectively guiding the robot towards the center of the treadmill while maintaining a forward orientation. The two hind legs and the two front legs are controlled together, respectively, with the same input (as if they were together an individual leg of the sagittal model).

6.2.1. Trunk state estimation

Using the springs as torque sensors (by measuring the deflection θ − q ), the radial force acting along the line connecting shoulder and foot is estimated. When this force is higher than a small threshold (0.25 N), the foot is considered in contact with the ground. A simple state machine handles the logic to determine the current locomotion phase. The joint position readings (and joint velocities) are used in combination with the IMU measurements for state estimation. In particular, the state estimation assumes a flat ground and utilizes joint measurements during double support to accurately estimate the state of the trunk. The IMU measurements are used to reconstruct the state as soon as the system exits the double support phase.

6.2.2. Identification of a natural motion manifold

Using the parameters reported in Table 1, we apply the methodology presented in Section 5 to model the robot, obtain an NMM, and design the proposed controllers. After validating the controllers in simulation, we encountered a kinematic limitation during preliminary hardware experiments. This limitation, caused by hard joint end-stops, prevented the realization of high-energy gaits. To meet the leg kinematic limits, the gaits are studied for a system with a 25% higher inertia with respect to CAD model estimations ¹⁰ . Figure 12 presents the resulting family of autonomous oscillations for the energy range between 4.8 J to 7 J.

6.2.3 Tuning of the controllers

We use the same DLQR gains and feed-forward injection profiles as presented in Figures 14 and 15, respectively. The saturation bounds for the output of the stabilizer (presented in (50)) are set to ±0.05 rad for the attack angles, and ±2 mm for the attack lengths. For the energy controller, the gains are tuned to β = 2.5 and ρ = 0.25.

6.2.4. Energy and apex estimation

Estimations of the energy and the upper apex of the COM trajectory are necessary for stabilization and excitation of the gait. However, they are challenging to reconstruct accurately and robustly. In our case, the most reliable way to estimate the energy was to evaluate it during double support (where we can rely on joint measurements) and hold the estimated value constant during the other phases. To robustly estimate the upper apex (the Poincaré section hit, used by the stabilization controller) we use the upper apex of the previous stride, which is retrieved by detecting when $\dot{z}$ crosses zero. Additionally, these estimations are filtered with a first-order low-pass filter with a time constant τ = 0.5 s to obtain smooth average values over the last strides.

6.3. Experimental results

We present experiments for two cases.

6.3.1. Controlling to energy targets

In this experiment, the robot starts still on the treadmill (with stand-still energy due to gravity ∼ 4.8 J), and the desired energy E_d. is set to an initial value E₁ = 5.5 J, which roughly corresponds to an average forward speed ¹¹ ∼0.14 m/s. After a short transient, the robot stabilizes and runs at E₁. Then, the desired energy is increased to E₂ = 5.9 J (average ∼0.20 m/s), and the system converges to it in a few seconds, gaining speed. Figure 20 reports the time evolution of the measured mechanical energy against the target energy, while Figure 21 presents the tracked COM height and trunk orientation, together with the reconstructed forward velocity. Snapshots of the robot performing the gait at E₂ are presented in Figure 22. The projections of the recorded evolution of the trunk state during the experiment are depicted in Figure 23. Additionally, the measured motor positions against the angles of the corresponding joints for one leg are reported in Figure 24 for a time interval after the robot stably converged to running at the energy E₂. Finally, the recorded outputs of the stabilization and injector controllers are presented in Figures 25 and 26, respectively.

Figure 20.

Evolution of the measured mechanical energy against the target energy. The target energies are color-coded. An estimation for the energy (Raw, solid gray) is calculated using the joint and IMU measurements. The value of the energy used for the controllers (Filtered, solid black) is obtained by low-pass filtering the estimated value, which is updated at the end of each stride.

Figure 21.

Recorded trajectories of the trunk height and orientation, and of the reconstructed and averaged forward velocity. The measurements demonstrate highly regular and robust motion displayed by the robot when performing the identified gait at two reference energy levels.

Figure 22.

Side and front-view snapshots of the elastic quadruped robot bounding at the desired energy 5.9 J on a treadmill. The side-view shots are reported with respect to a frame solid with the treadmill to highlight the forward motion. The COM trajectory is depicted in solid red for the stance phase and dashed for the flight phase. The resulting motion is highly dynamic and closely resembles the identified passive gait of the reduced-order model.

Figure 23.

State space projection of the recorded trunk trajectories. The color of the trajectory indicates the energy reference in that time interval (red for E₁ and blue for E₂, as in Figure 20). The experiments show that a stable and regular orbit is achieved for both energies. These feedback-induced limit cycles do not correspond precisely to the orbits from the manifold evaluated using the ideal reduced model, but they still resemble its shape. This is expected since friction, impacts, and other unmodeled effects affect the motion. The relevant result is that a small excitation designed via the NMM generates regular and predictable resonant motion.

Figure 24.

Comparison of motor angle against the corresponding joint angle, showing high regularity and large spring deflections. A small motion of the motors is sufficient to excite the system’s natural dynamics and cause large joint and trunk oscillations.

Figure 25.

Stabilization control errors and output to adjust the attack configuration. Similarly to the results of the simulations with losses reported in Figure 18, in the real scenario, the DLQR stabilizer cannot drive the error to zero due to losses and other unmodeled dynamics. However, the stabilization controller makes the manifold attractive, so that the system settles onto an orbit in a neighborhood of the natural motion manifold, like in the simulations.

Figure 26.

Radial and polar offsets from the energy regulator. The injection is higher during the transition to higher energy, and it converges to steady lower values when the desired energy level is reached.

6.3.2. Response under perturbations

A series of heavy flat objects of approximately 1 cm height and 10 cm in diameter are placed on the treadmill while the robot is running at E₂. When the robot runs over the obstacles, it is initially pushed away from the current orbit, but it subsequently returns to the original cycle after a transient. The results of the experiment are summarized in Figures 27 and 28, which show the robot recovering from a perturbation over a short time interval. See Extension 1 for the video recordings of the experiments.

Figure 27.

Experimental recovery from perturbations due to impacts with solid obstacles. The feet of the robot strike the obstacles at t ≈ 26.5 s, disrupting the desired oscillation, and driving the robot to almost jump in place. During the stabilization, the DLQR performs corrections of the attack configuration to stabilize the gait, while the injector manages the energy variations. After a short transient, the cycle is again regular at the target energy.

Figure 28.

Snapshots of the robot recovering from perturbations: the compliant response of the leg passively absorbs the impact, while the robust excitation controller restores regular oscillations at the desired speed.

6.4. Discussion

The experiments show that the robot quickly converges to highly regular periodic trajectories, largely taking advantage of the natural dynamical response of the soft elastic elements for the execution of the gait. The regularity of the oscillations of the trunk can be appreciated in Figure 21 and of the joint and motor coordinates in Figure 24. The latter plot also shows that comparatively small motor oscillatory commands (generated from the injector, Figure 26) lead to larger oscillations in the joints and, thus, of the trunk.

For the i‐th joint, the actuator work can be evaluated as

$L_{\mathrm{i}}=\int {\dot{\theta }}_{\mathrm{i}}{\tau }_{\mathrm{i}}dt,$

(51)

where τ_i is the torque acting on the joint (i.e., the spring torque). By relying on the excitation of natural dynamics, the proposed approach implicitly avoids negative actuator work, which is generally linked to inefficiencies in locomotion (Hutter et al., 2013; Remy, 2011). This is because the task of the controller is to reinject energy to sustain an oscillation (positive actuator work) rather than compensating for undesired dynamics and forcing the system to evolve along specific trajectories. To demonstrate this feature, Figure 29 presents the mechanical power exchanges and total work performed by the actuators. During the double support phase, the controller dissipates almost no energy. At each stride, the system naturally dissipates all the injected energy (positive work) via friction and impacts. The time intervals where there is negative work can be attributed mainly to the internal PID of the motors, which briefly acts as a damper at touchdowns and when the robot is in the air, since it must regulate the leg to the attack configuration.

Figure 29.

Total mechanical power of the front (FL) and hind legs (HL), demonstrating predominantly positive work performed by the actuators. The respective stance and flight phases of the legs are also reported. For the i-th joint, the power is calculated as ${\dot{\theta }}_{\mathrm{i}}{\tau }_{\mathrm{i}}$ , where τ_i is the torque acting on the joint (equivalent to the spring torque). The absence of negative work implies that the actuators are not actively removing energy, an essential aspect for achieving energy-efficient motion.

To illustrate the energetic benefits of our approach, we compare the Cost of Transport (CoT) of our method against standard rigid actuation across three scenarios: the nominal case, simulations with dissipation, and experimental results. For each scenario, these comparisons are made while the robot executes the same bounding gait at the sustained energy level E₂. We evaluate the CoT for our strategy using the positive mechanical actuator work (Remy, 2011)

${\text{CoT}}_{\text{mech}}=\frac{\int \max ({\mathit{\tau }}^{\top }\dot{\mathit{\theta }},0)dt}{mgd},$

(52)

where d is the distance traveled during the time interval of interest. The rigid case—where oscillations cannot be exploited—is considered as a benchmark

${\text{CoT}}_{\text{mech}}^{\text{rigid}}=\frac{\int \max ({\mathit{\tau }}^{\top }\dot{\mathit{q}},0)dt}{mgd}.$

(53)

Table 2 compares the costs and also presents the peak mechanical powers for the experiment. The comparison highlights the difference between the hyper-efficient nominal case corresponding to zero cost, achieved through the purely passive oscillations in NMMs, and the cost associated with the realization of the gait with non-compliant robots. Simulations with joint friction losses (from Section 5.4) show that there is almost a 90% improvement in efficiency with respect to the rigid case. Although the improvements observed in the experiments are not as high, the advantage of leveraging nonlinear oscillations with the proposed method, relative to the benchmark, is still evident. The reduction of more than 50% in peak power allows for a reduction in the actuators’ size and weight. The main obstacles to achieving the theoretical efficiency levels appear to be inefficiencies due to the robot’s design and limited motor performance. Additionally, we remark the fact that the passive gaits in the NMM are not explicitly optimized for minimal mechanical CoT in the traditional sense. This suggests that optimization of both the robot and the NMM could lead to efficient locomotion.

Table 2.

Comparison of Mechanical Costs of Transport at E₂.

CoT	Nominal	Simulation	Experiment (peak power)
Proposed (52)	0	0.17	0.77 (22 W)
Benchmark (53)	1.13	1.33	1.37 (48 W)

The validity of the natural motion manifold analysis is further confirmed by comparing the oscillation periods extracted from the experiments with the theoretical periods obtained from the passive gaits of the reduced model. This comparison is presented in Figure 30, which displays closely matched oscillation periods for both of the energy levels investigated in the experiments. This outcome is relevant and not trivial because the nonlinear nature of the oscillations results in non-constant, energy-dependent periods. This demonstrates again the ability of the simple reduced model to capture relevant features of more complex dynamics.

Figure 30.

Oscillation periods of the reduced model’s gait matching the experimental measurements. The nonlinear characteristic of the oscillation leads to a non-constant frequency, which varies depending on the energy level. Despite its simplicity, the reduced-order model can capture and accurately model this feature.

It was observed experimentally that the direction and the magnitude of the motor commands generated by the injector controller had the most effect on driving the desired oscillations, while the stabilizer mainly countered the rise of irregular or multi-periodic behavior, but it is unable alone to correct an inconsistent excitation. The resulting gaits are qualitatively very similar to the passive gaits identified for the reduced model: the foot pattern is the same, and the periods of oscillations match. However, there are some quantitative differences, such as a higher rotational speed of the trunk and a slightly lower forward speed achieved at lower energies. These discrepancies were expected since a part of the unmodeled dynamics (such as friction and impacts between the feet and the ground, which are difficult to model accurately) considerably influence how close the real gaits can come to the passive ones. This can also be seen in the simulations reported in Figure 18, which feature dissipation.

7. Conclusion

Driven by efficiency and robustness arguments for the realization of repetitive tasks, we have introduced the notion of Natural Motion Manifolds and a control method capable of locally stabilizing and exciting these geometrical structures. The proposed controller was derived via the extension of a stabilization method (based on the linearized controlled Poincaré map) from one orbit to a family of orbits. In particular, the controller exploits the intersection of a Poincaré section with a surface comprising autonomous evolutions, thus obtaining a curve used to construct a coordinate system that decouples the transversal and energy dynamics. This drives the system on passive solutions characterized by zero control effort. The method was validated both in simulation and experimentally on a serial elastic quadruped. A bounding gait at a desired velocity was realized without performing dynamic compensations, not even for gravity torques: the natural dynamics were simply excited via small motor commands. The experiments show that the robot takes advantage of the natural dynamical response of the soft elastic elements to execute the gait. To our knowledge, this is one of the first experimental results to demonstrate that an excitation command designed via the natural dynamics generates controlled, predictable, and useful resonant behavior. Moreover, the efficiency of the energy injection principle was demonstrated on hardware. These outcomes provide further insights into the usefulness of energy-conservative oscillations (e.g., passive gaits) for the planning and control of robots and the role of elastic elements in energy regeneration.

In the future, we plan to exploit the bifurcations to switch to different NMMs, thus increasing the flexibility of this idea. In the case of locomotion, high-energy oscillations associated with very dynamic and fast locomotion gaits could be efficiently excited by switching to an appropriate manifold at the right moment. Moreover, since the experimental platform featured only servo motors, which lacked torque control, the control principle should also be tested on different systems, and their efficiency should be investigated. Finally, the optimization of NMMs, such as a co-evolution of hardware and controllers to minimize the cost of transport, should be explored to achieve future enhancements to the already impressive capabilities demonstrated by quadrupeds (Hoeller et al., 2024).

We believe that, for the next technological step, robots must take full advantage of their natural dynamics under both the mechanical and control design aspects. This design philosophy does not limit itself to locomotion, but it should be applied to any robotic system in the pursuit of increased performance and lower costs. This, however, requires the inversion of a standard paradigm in robotics: instead of treating the robot as a generic motion realization machine, the intrinsic capabilities of a robot should be identified first and then fully exploited to realize a task.

Supplemental Material