In general, the walking and running gaits of legged robots are naturally unstable and challenging to control. Hybrid systems such as these are difficult to work with due to the discontinuities that occur in state and dynamics at hybrid events, such as toe touchdown. These discontinuities violate assumptions of standard controllers designed for purely continuous systems, and work is ongoing to adapt these controllers for hybrid systems [1, 20]. One strategy for hybrid control is to cancel out the effects of hybrid events by working with an invariant subsystem [39, 7, 9]. We propose instead that the effects of hybrid events are valuable due to rich control properties that can be used to stabilize trajectories of a hybrid system.
Several works have examined the utility of controlling hybrid event conditions to improve system stability without any closed-loop continuous-domain control [33, 32, 2]. For example,  found that for the paddle juggler system, paddle acceleration at impact uniquely determines the local stability properties of a periodic trajectory, Fig. 1. Other works [9, 11]
generated open-loop swing leg trajectories that produced deadbeat hopping of a SLIP-like system. Each of these works found that controlling a hybrid system only at the moment of a hybrid event is sufficient to provide stabilization. So far, however, these results have only been produced for each specific problem structure and a clear connection between these works has yet to be established.
In this work, we propose the concept of hybrid event shaping (HES), which describes how hybrid event parameters can be chosen to affect the stability properties of a periodic orbit. We also propose methods to produce values of these hybrid event parameters to optimize a stability measure of a trajectory. This approach is tested on both existing examples from [33, 32] as well as on a new bipedal robot controller.
Ii-a Hybrid system
Hybrid systems are a class of dynamical systems which exhibit both continuous and discrete dynamics [36, 23]. Following the adaptation of  in , we define a hybrid dynamical system for continuity class as a tuple where:
is the set of discrete modes.
is the set of discrete transitions forming a directed graph structure over .
is the collection of domains where is a manifold with corners .
is a collection of time-varying vector fields with control inputswhere is the space of allowable control inputs.
is a collection of guards where for each is defined as a sublevel set of a function, i.e. .
is a map called the reset map that restricts as for each .
An execution of a hybrid system  begins at an initial state . With a particular input , the system follows the dynamics on . If the system reaches the guard surface , the reset map is applied and the system continues in domain with the corresponding dynamics defined by . An execution of a hybrid system is defined over a hybrid time domain which is a disjoint union of intervals . The flow describes how the system evolves from some initial time and state for some length of execution time .
These hybrid systems may exhibit complex behavior including sliding , branching , and Zeno phenomena . To preserve the validity of local linearization theory, which is core to this work, we assume that transition surfaces are isolated with transverse intersections [23, 16] and no Zeno executions occur.
Ii-B Saltation matrix
Hybrid systems can be divided into continuous domains and discrete hybrid events. For each of these components, the linearized variational equations  describe how perturbations of state at the beginning of the phase evolve to the end of the phase. In continuous domains, variational equations can be derived  and discretized into a mapping equivalent to the linearized discrete dynamics matrix in . The variational equations of hybrid events are described by the saltation matrix, which maps perturbations at an infinitesimal time step before impact to an infinitesimal time step afterward . This saltation matrix, , describes the transition between modes and that occurs at time with and with some input . Following the formulation from [6, 26], the saltation matrix is,
The saltation matrix approximates the first-order change in perturbations in state before the hybrid event to perturbations afterward  such that:
where h.o.t. represents higher order terms.
Ii-C Periodic Stability Analysis
A dynamical system has a periodic trajectory (orbit) with period if for some initial condition , there exists a solution where for all and . Perturbations around can be mapped to perturbations after period by a linearized mapping known as the monodromy matrix, :
such that to the first order, .
Following a formulation similar to , the monodromy matrix can be computed by sequentially composing the linearized variational equations in each continuous domain () and the saltation matrices () at each hybrid event. For a hybrid periodic orbit with domains, the monodromy matrix can be formulated as:
The monodromy matrix determines local asymptotic orbital stability (which we refer to simply as stability). For nonautonomous systems, stability is determined by the maximum magnitude of the eigenvalues,. We refer to this as the stability measure, , where a trajectory is stable when . Autonomous systems always have an eigenvalue that is equal to 1 since for non-time varying dynamics, perturbations along the flow of the orbit will by definition map back to themselves after period . Assuming non-convergence in this direction is allowable, for autonomous systems is based on the remaining eigenvalues.
Iii-a Hybrid event shaping
Hybrid events can greatly affect the stability properties of an orbit because of the discontinuous changes that are made to perturbations. The saltation matrix allows for an explicit understanding of how to perform “hybrid event shaping” (HES), i.e. choosing hybrid event parameters (such as timing, state, input, and higher order “shape parameters”) to improve the stability of a periodic trajectory. The key insight is that hybrid event shaping introduces a generalizable method to stabilize hybrid systems that is independent from continuous-domain feedback, but that can work in concert with it.
In general, the open-loop continuous variational equations of a hybrid system are functions of initial and final time, initial state, and system dynamics. However, it is challenging to alter any of these parameters because changes will propagate through the rest of the trajectory and periodicity may be violated, though we present a trajectory optimization method below that can handle this. The saltation matrix is a function of nominal event time, state and dynamics, but additionally may be a function of higher order shape parameters that do not influence the dynamics of the system. These parameters arise from the derivatives of the guards ( and ) and reset maps ( and ) but are not present in the guard, reset map, or vector field definitions themselves. Therefore, shape parameters have absolutely no effect on the nominal trajectory and can be chosen completely freely.
One example of a shape parameter is the angular velocity of a massless leg of a spring-loaded inverted pendulum. Since a massless leg does not induce any torque in the air or forces at touchdown, only the position of the leg at touchdown affects the trajectory of the body. However, leg velocity appears in the saltation matrix and has a significant effect on orbital stability .
For more complex models of robots, there may not be any physical shape parameters that can be tuned. For example, legged robots with non-massless legs can not vary leg velocity at impact without also changing the periodic orbit. These cases can be handled by adjusting the trajectory at the same time as applying HES, as we show in Sec. III-E, or by adding additional virtual hybrid events.
Iii-B Virtual hybrid events
Certain control systems naturally have discontinuities in control inputs, such as bang-bang control, sliding mode control, or systems with actuators that have discretized inputs (i.e. on-off). These discontinuities in control cause instantaneous changes in the dynamics of the system, resulting in a virtual (as opposed to physical) hybrid event. Virtual hybrid events act no differently than physical hybrid events and induce saltation matrices with shape parameters to be tuned for stability. Even for systems where discontinuous control inputs are not necessary, the addition of virtual saltation matrices and shape parameters allows for a greater authority to improve stability.
Iii-C Stability measure derivative
Our goal is to determine the optimal choice of hybrid event parameters that minimizes the stability measure of a trajectory. Since directly computing eigenvalues in closed-form is not generally feasible, one solution is to use numerical methods to perform optimization . However, this strategy becomes untenable for high dimensional problems. With the saltation matrix available in closed form (1), derivatives of the stability measure can be directly computed, allowing for the use of more efficient optimization methods to select hybrid event parameters.
Assuming that the monodromy matrix depends continuously on each shape parameter , the eigenvalues of are always continuous with respect to . For a diagonalizable , the derivative of the eigenvalues with respect to can be computed in closed form . Assume that matrix has simple (non-repeating) eigenvalues, , and let and
denote the left and right eigenvectors associated with. Then the derivative is:
For matrices with eigenvalues that repeat, the derivatives of the repeated eigenvalues can be calculated similarly with a matrix of associated eigenvectors .
can be found using the derivative product rule, which simplifies if each shape parameter only appears in one saltation matrix. We make this assumption here to improve computational efficiency, but it is not required generically. Without loss of generality, take , so that:
Substituting (6) into (5) allows us to compute the derivative of the stability measure with respect to each of the shape parameters. Eq. (5) is not valid for non-diagonalizable monodromy matrices. However, the guaranteed continuity of the stability measure allows for finite-difference methods to be used in any non-diagonalizable cases.
The derivative computation from (6) can be adapted for changes in and as well. Without loss of generality, consider again . For hybrid event time , the derivative can be computed in closed-form. Additionally, the derivatives and are non-zero and can be computed through numerical methods. The product rule expansion of consists of additional terms compared to (6) but otherwise can be computed similarly. for hybrid event state can be computed this same way.
Iii-D Shape parameter stability optimization
Optimization techniques  are able to select optimal hybrid event parameters that minimize the stability measure. Two optimization methods are presented here: the first optimizing the shape parameters without affecting the dynamics of the nominal orbit, and the second optimizing both the hybrid events and periodic orbit simultaneously.
Shape parameters are powerful because they do not appear in the dynamics of the system and have no effect on the nominal trajectory. This means that for a given periodic trajectory, the shape parameters can be chosen freely. We use an optimization framework to choose these shape parameters with the goal of optimizing the stability measure of a trajectory,
The ability to compute derivatives of the stability measure allows for this optimization to be more computationally efficient. The examples below show how this optimization method is able to reproduce swing leg retraction in a one-legged hopper system by determining optimal inputs of shape parameters to minimize the stability measure.
Iii-E Trajectory optimization with hybrid event shaping
Some systems do not have shape parameter terms in their saltation matrices or do not have enough to sufficiently improve stability. In these cases, we can change the trajectory of the system itself so that the timing, state, and input parameters of the continuous variational matrices and saltation matrices improve stability properties. However, it must be ensured that the dynamics, periodicity, and other constraints of the system are obeyed.
Trajectory optimization methods are a class of algorithms that aim to minimize a cost function while satisfying a set of constraints . For dynamical systems, these costs are generally functions of states and inputs, with constraints imposed on system dynamics and any other physical limits . For specific problems, other aspects of the system may be added into the cost or constraint functions such as design parameters or minimizing time [12, 25]. Here we propose including the stability measure in the cost function to search for optimally stable trajectories. Eq. (8) gives the simplest form of this trajectory optimization problem, with periodicity and dynamics constraints being enforced, where dynamics constraints obey continuous dynamics in each domain and reset mappings at each hybrid event .
Iv Examples and Results
Here we demonstrate how HES can improve the stability of periodic trajectories for a variety of hybrid systems without any use of continuous-domain feedback control. While continuous-domain feedback could be implemented into any of these systems and should be in practice, these examples emphasize the stabilization capabilities of HES alone.
The first two examples describe how previously discovered results, paddle juggling  and swing leg retraction , fit into an HES framework. The final example demonstrates how HES can be used even without any shape parameters and how virtual hybrid events can help stabilize a more complicated biped walking system.
Iv-a Paddle juggler
The paddle juggler system , bouncing a ball with a paddle, is known to be stabilized by impacting the ball with a paddle acceleration in a range of negative values, (10), Fig. 1. The system state consists of the ball’s vertical position and velocity such that . This periodic hybrid system can be defined with two hybrid domains (descent and ascent) connected by two guards (impact and apex). The domain represents the ball’s aerial descent phase where and represents the ball’s aerial ascent phase where . The guard set is defined when the ball impacts the paddle, where the paddle follows a twice differentiable trajectory . The reset map is defined by a partially elastic impact law, , with a coefficient of restitution . The continuous dynamics in both domains follow unactuated ballistic motion: , where is the acceleration due to gravity.
Using these definitions, the saltation matrix (1) between domains and is constructed:
Observe that appears in the saltation matrix even though it does not appear anywhere in the definition of the guards, reset maps, or vector fields of the system, making it a shape parameter that can be chosen independently from the periodic orbit.
The guard set is defined when the ball reaches the apex of its ballistic motion. Its reset map is identity and there is no change in dynamics. Therefore, is identity and has no effect on the variations of the system.
The continuous variational matrices of the system can be written exactly in closed form: where is the total time spent in the air and also the period of the system. The periodicity of the system means the ball spends an equal time, , ascending as descending.
The monodromy matrix, is constructed by composing together the continuous variational matrices and saltation matrices such that .
For a given periodic bouncing trajectory, the monodromy matrix is almost fully defined, with the exception of the shape parameter, in . Given the 2-dimensional state space of this problem, the eigenvalues for any given value can be solved for explicitly. We can then solve exactly for where , giving:
This is confirmed in , where the simple dynamics of the system allowed the return map to be computed explicitly without the saltation matrix. However, that computation is not generally possible for more complex dynamics.
Iv-B 2D hopper
The spring-loaded inverted pendulum (SLIP) is a popular model for dynamic legged locomotion [24, 14, 31]. This simple hopper model is effective at capturing dynamic properties of animal and robot locomotion  and has been used as a test bed for different types of hybrid controllers .
Iv-B1 2D hopper hybrid model
Consider a point mass body with a massless leg consisting of a spring, damper, and linear actuator all in parallel, Fig. 2. This system has two domains (flight and stance) connected by two guards (touchdown and liftoff). The actuator is activated while in the air to preload the spring, but immediately releases at touchdown and provides no forces during stance. For a periodic trajectory to occur, the actuator must preload the same amount of energy that is dissipated by the damper during stance. The only control authority that exists is of the leg angle in the air, which only affects the dynamics of the body at touchdown.
During flight (), the state of the hopper is represented by the horizontal velocity, vertical position, and vertical velocity of the body: . Horizontal position is not included because it is not periodic. During stance, the body position and is defined with the toe at the origin. Horizontal position is added back into the state of the hopper such that . In flight, the dynamics of the body follow ballistic motion, while in stance there are also forces applied by the spring and damper.
The touchdown guard is defined by the preload length of the leg and angle of the leg such that . The liftoff guard is crossed when the force exerted by the spring-damper, , becomes zero: . There is no change in physical state of the system at the hybrid events and the reset maps only characterize the change in coordinates between domains.
Given a set of model parameters, a trajectory from an initial condition depends only on and . is held fixed, but is modulated from its nominal position at time in two ways. A proportional feedback term with gain K is added to stabilize the forward velocity of the system around a nominal and the angular velocity of the massless leg is also free to be chosen.
and are shape parameters that are used to stabilize this system.
Iv-B2 2D hopper HES results
For a chosen initial apex height of with a forward velocity of , and were chosen to produce a nominal orbit, though the following results generalize for any choice of feasible values.
With fixed touchdown parameters , the system is highly unstable. The shape parameters and can be optimized following (7) to improve the stability of this orbit. Doing so results in optimal shape parameters that stabilize the trajectory, Table I. Setting is a slow retraction rate, but there exists an interval of values that give equivalently minimal stability measures for a fixed value.
|Shape parameters||K||Stability Measure|
|Optimal (zero seed)||0.129||-0.015||0.948|
|Optimal (alternate seed)||0.129||-0.589||0.948|
The results were confirmed in simulation by initializing 20 perturbed points in a radius ball around the nominal initial condition. Each of these trials was simulated for steps and the error (2-norm of the difference in perturbed state and nominal state ) at apex was recorded at each step, Fig. 3. The zero shape parameter trajectories are not shown in the figure as every trial diverged within just 5 steps.
Iv-B3 2D Hopper Discussion
The feedback term of (11) is based on the Raibert stepping controller , which was utilized to great success for stabilizing early running robots. Other works have found that this simple controller is effective on more complex models .
Another stabilizing property of legged locomotion that has been studied is swing-leg retraction [32, 18, 5]. It was noted in  that a 2D SLIP was able to run stably if it impacted the ground within a range negative angular leg velocities .
The results of a negative and positive agree with qualitative expectations from  and . While a formal equivalency can not be proven, this is significant because the HES shape parameter optimization has no a priori knowledge nor any applied constraints that would bias its results. HES is able to synthesize two independently generated controllers and produce shape parameter values that stabilize an orbit. This evidence supports the potential for HES to explore other stabilizing shape parameters that are not as well studied.
Iv-C Walking biped trajectory optimization
For a legged system with non-massless legs, the leg velocity shape parameter disappears as it is no longer independent of the trajectory. Without shape parameters, an HES trajectory optimization can choose timing, state, and input parameters along with virtual hybrid events to discover stable orbits.
Iv-C1 Walking biped hybrid model
In this example, we consider a fully-actuated compass walker  with knees, Fig. 4. This biped model consists of two legs connected by an actuated hip joint. Each leg is separated into two sections, the upper leg (thigh) and lower leg (shank), which are connected by an actuated knee joint that has a hard stop when the thigh and shank are aligned. The ankles are also actuated.
We restrict the gaits to be left-right symmetric and exclusively consist of single stance phases. The stance leg is locked such that its shank and thigh are aligned with each other until liftoff. There are 3 points of actuation at the hip, swing knee, and stance ankle. The state space is defined by three angles relative to vertical: stance leg, swing thigh, and swing shank, denoted .
This system has two domains. is the unlocked knee domain where the swing leg thigh and shank are able to swing freely while we enforce that . is the locked knee domain where the thigh and shank are constrained to be aligned with each other (). In this domain, there are only two actuators because the swing knee can not exert torque. The dynamics of this model are described in .
The kneestrike guard set, between the unlocked and locked knee domains, is and the touchdown or toestrike guard set is . The reset maps at kneestrike and toestrike are also computed in .
We add discrete changes in the inputs that induce virtual hybrid events to analyze their utility in stabilizing walking trajectories. Specifically, we choose to include 5 virtual hybrid events in and 2 more virtual hybrid events in , where the values of inputs between virtual hybrid events are decision variables for the optimization. The virtual guard functions are chosen such that for the first 5 virtual hybrid events and for the last 2 virtual hybrid events for some offset that is also chosen by the optimization.
A direct collocation method was used with the cost consisting of the stability measure and a regularization on the input. Dynamics and periodicity constraints were included along with a ground penetration constraint. The initial conditions of the system, given as the state after touchdown, were allowed to vary within a bounded range.
Iv-C2 Walking biped HES results
In this experiment, three trajectories were compared to examine how HES can generate stable trajectories and the effect that virtual hybrid events have in further improving stability. A trajectory without HES was produced as a control, with its objective to minimize energy expended by using just an input regularization term in the cost. This minimum energy (ME) trajectory is comparable to how conventional robot locomotion trajectories are generated. Two HES trajectories were generated, one with virtual hybrid events (HES w/ VHE) and one without (HES w/o VHE).
The ME trajectory is highly unstable, while the both HES trajectories are stable with the trade off of a higher input cost. Specifically, HES w/ VHE has the lowest stability measure and highest energy cost, whereas HES w/o VHE was in between for both stability and cost, Table II.
|Trajectory||Stability Measure||Energy Cost|
|HES w/o VHE||0.4715||1.337|
|HES w/ VHE||0.3266||3.450|
The stability properties of the generated trajectories were confirmed through simulation. 50 trials of each trajectory were initialized with perturbations in position and velocity between . Over a sequence of 10 steps, the state error at each step was tracked for each trial. Fig. 5 shows that after 10 steps, the HES trajectories have nearly converged back to the nominal trajectory whereas the ME trajectories diverge quickly. The HES w/ VHE trajectory converges to a smaller error after 10 steps compared to the HES w/o VHE trajectory, which supports the findings of the stability measure.
While the idea of hybrid event control is not novel, hybrid event shaping provides a generalized method to analyze the stability of hybrid orbits and select hybrid event parameters to optimize stability. HES unifies results of previous simple hybrid event controllers while also being compatible with trajectory optimization techniques to produce stable trajectories for complex systems. HES also improves on previous stability optimization methods by computing the derivative of the stability measure to increase computational efficiency. Compared to previous work, HES does not rely on prior understanding to identify and choose stabilizing hybrid event parameters.
Throughout this paper, there was no use of continuous-domain feedback that is commonly utilized in hybrid control. We believe that hybrid event shaping is a piece of the puzzle that can be used in conjunction with continuous-domain feedback to improve the success rate of robots performing dynamic behaviors in real-world settings. In the future, we would like to synthesize HES methods with continuous-domain feedback controllers to produce even more stable closed-loop trajectories.
-  (2018) A model predictive control framework for hybrid dynamical systems. IFAC-PapersOnLine 51 (20), pp. 128–133. External Links: Cited by: §I.
-  (2010) Control of underactuated planar hexapedal pronking through a dynamically embedded slip monopod. In IEEE International Conference on Robotics and Automation, pp. 4721–4727. Cited by: §I.
-  (2018) Stability analysis using monodromy matrix for impacting systems. IEICE Transactions of Fundamentals of Electronics, Communications, and Computer Science E101-A (6). Cited by: §II-C.
Practical methods for optimal control and estimation using nonlinear programming, second edition. Second edition, SIAM, . External Links: Cited by: §III-E.
-  (2010-06) Swing leg control in human running. Bioinspiration & Biomimetics 5. External Links: Cited by: §IV-B3.
-  (2018) On contraction analysis for hybrid systems. arXiv preprint arXiv:1811.03956. Cited by: §II-B, §II-B.
-  (2011-09) Dimension reduction near periodic orbits of hybrid systems. In IEEE Conference on Decision and Control, pp. . External Links: Cited by: §I.
-  (2007) Passive dynamic walking with knees: a point foot model. Master’s Thesis, Massachusetts Institute of Technology. Cited by: §IV-C1, §IV-C1, §IV-C1.
-  (2014) Deadbeat control with (almost) no sensing in a hybrid model of legged locomotion. In International Conference on Advanced Mechatronic Systems, Vol. , pp. 475–480. Cited by: §I, §I.
-  (1999-12) Templates and anchors: neuromechanical hypotheses of legged locomotion on land. Journal of Experimental Biology 202 (23), pp. 3325–3332. External Links: Cited by: §IV-B.
-  (2020-05) Planning for the unexpected: explicitly optimizing motions for ground uncertainty in running. In IEEE International Conference on Robotics and Automation, External Links: Cited by: §I.
-  (1987) Direct trajectory optimization using nonlinear programming and collocation. Journal of Guidance, Control, and Dynamics 10 (4), pp. 338–342. External Links: Cited by: §III-E.
-  (2012) Differential equations, dynamical systems, and an introduction to chaos. Academic Press. Cited by: §II-B.
-  (2006) The dynamics of legged locomotion: models, analyses, and challenges. SIAM Review 48 (2), pp. 207–304. External Links: Cited by: §IV-B.
-  (2014-07) Dynamics at a switching intersection: hierarchy, isonomy, and multiple sliding. SIAM Journal on Applied Dynamical Systems 13, pp. 1082–1105. External Links: Cited by: §II-A.
-  (2016) A hybrid systems model for simple manipulation and self-manipulation systems. The International Journal of Robotics Research 35 (11), pp. 1354–1392. External Links: Cited by: §II-A, §II-A, §II-A.
-  (2012) On manifolds with corners. In Advances in Geometric Analysis, Advanced Lectures in Mathematics, Vol. 21, pp. 225–258. Cited by: item 3.
-  (2011) The optimal swing-leg retraction rate for running. In IEEE International Conference on Robotics and Automation, Vol. , pp. 4000–4006. External Links: Cited by: §III-C, §IV-B3.
-  (2017) An introduction to trajectory optimization: how to do your own direct collocation. SIAM Review 59 (4), pp. 849–904. External Links: Cited by: §III-E, §III-E.
-  (2021-12) iLQR for piecewise-smooth hybrid dynamical systems. In IEEE Conference on Decision and Control, Note: To appear Cited by: §I, §II-A.
-  (1964) On eigenvalues of matrices dependent on a parameter. Numerische Mathematik. Cited by: §III-C.
-  (2004) Dynamics and bifurcations of non-snooth mechanical systems. Springer-Verlag Berlin Heidelberg. Cited by: §II-B, §II-C, §II-C.
-  (2003) Dynamical properties of hybrid automata. IEEE Transactions on Automatic Control 48 (1), pp. 2–17. External Links: Cited by: §II-A, §II-A.
-  (1990) The mechanics of running: how does stiffness couple with speed?. Journal of Biomechanics 23, pp. 65–78. External Links: Cited by: §IV-B.
-  (2005-01) Open-loop stable running. Robotica 23, pp. 21–33. External Links: Cited by: §III-E.
-  (1995) Calculation of lyapunov exponents for dynamic systems with discontinuities. Chaos, Solitons & Fractals 5 (9), pp. 1671–1681. External Links: Cited by: §II-B.
-  (2006) Numerical optimization. Springer. Cited by: §III-D.
-  (2009) The spring loaded inverted pendulum as the hybrid zero dynamics of an asymmetric hopper. IEEE Transactions on Automatic Control 54 (8), pp. 1779–1793. External Links: Cited by: §IV-B3.
-  (1984) Experiments in balance with a 3D one-legged hopping machine. The International Journal of Robotics Research 3 (2), pp. 75–92. External Links: Cited by: §IV-B3, §IV-B3.
-  (2007-10) A simple model for clock-actuated legged locomotion. Regular and Chaotic Dynamics 12 (5). External Links: Cited by: §IV-B.
-  (2002) A movement criterion for running. Journal of Biomechanics 35 (5), pp. 649–655. External Links: Cited by: §IV-B.
-  (2003) Swing-leg retraction: a simple control model for stable running. Journal of Experimental Biology 206 (15), pp. 2547–2555. External Links: Cited by: §I, §I, §III-A, §IV-B3, §IV-B3, §IV.
-  (2001-11) Bouncing a ball: tuning into dynamic stability. Journal of Experimental Psychology: Human Perception and Performance 27, pp. 1163–84. External Links: Cited by: Fig. 1, §I, §I, §IV-A, §IV-A, §IV.
-  (1999-11) User’s guide for DIRCOL (version 2.1): a direct collacation method for the numerical solution of optimal control problems. In Technische Universitat Darmstadt, Cited by: §III-E.
-  (1997) A brief introduction to numerical analysis. Birkhäuser Basel. Cited by: §III-C.
-  (1998) Complementarity modeling of hybrid systems. IEEE Transactions on Automatic Control 43 (4). External Links: Cited by: §II-A.
-  (1992-12) Direct and indirect methods for trajectory optimization. Annals of Operations Research 37, pp. 357–373. External Links: Cited by: §III-E, §III-E.
-  (2001) On monodromy matrix computation. Computer Methods in Applied Mechanics and Engineering 190 (18). External Links: Cited by: §II-C.
-  (2021) Impact invariant control with applications to bipedal locomotion. arXiv preprint arXiv:2103.06907. External Links: Cited by: §I.
-  (2000) Dynamical systems revisited: hybrid systems with zeno executions. In Hybrid Systems: Computation and Control, External Links: Cited by: §II-A.