Previous research on movement control suggested that humans exploit stability to reduce vulnerability to internal noise and external perturbations. For interactions with complex objects, predictive control based on an internal model of body and environment is needed to preempt perturbations and instabilities due to delays. We hypothesize that stability can serve as means to render the complex dynamics of the body and the task more predictable and thereby simplify control. However, the assessment of stability in complex interactions with nonlinear and underactuated objects is challenging, as for existent stability analyses the system needs to be close to a (known) attractor. After reviewing existing methods for stability analysis of human movement, we argue that contraction theory provides a suitable approach to quantify stability or convergence in complex transient behaviors. To test its usefulness, we examined the task of carrying a cup of coffee, an object with internal degrees of freedom. A simplified model of the task, a cart with a suspended pendulum, was implemented in a virtual environment to study human control strategies. The experimental task was to transport this cart-and-pendulum on a horizontal line from rest to a target position as fast as possible. Each block of trials presented a visible perturbation, which either could be in the direction of motion or opposite to it. To test the hypothesis that humans exploit stability to overcome perturbations, the dynamic model of the free, unforced system was analyzed using contraction theory. A contraction metric was obtained by numerically solving a partial differential equation, and the contraction regions with respect to that metric were computed. Experimental results showed that subjects indeed moved through the contraction regions of the free, unforced system. This strategy attenuated the perturbations, obviated error corrections, and made the dynamics more predictable. The advantages and shortcomings of contraction analysis are discussed in the context of other stability analyses.

The evaluation of stability in complex and transient behaviors in human movement presents challenges as extant stability analyses require a known attractor state. We propose a variant of stability analysis, contraction analysis, that is suitable for a manipulation task with complex continuous dynamics. This differential method of stability assessment does not require knowledge of the attractor. We demonstrate this method on the control of a cart-and-pendulum system, a simple nonlinear underactuated object, as proxy for the human task of transporting a cup of coffee. Results show that humans indeed exploit contraction regions to mitigate perturbations and enhance predictability of the interaction with the complex object. Using stability to enhance predictability alleviates the need for an exact internal model. A focus on stability, including analysis of contraction, can advance the field of human movement science by allowing scientists to address questions on complex movement control.

Physical interaction with objects is ubiquitous in activities of everyday life, and humans show remarkable dexterity in handling objects with different degrees of complexity. Relatively simple examples that involve the manipulation of a free rigid object include grasping, lifting a book, and turning a key in a keyhole. In each case, the task involves moving a rigid object against a kinematic constraint. Physical interaction becomes particularly intriguing when the objects themselves have internal degrees of freedom that add complex dynamics to the interaction. A daunting example is cracking a whip: the flexible whip has infinitely many degrees of freedom that create challenging dynamics that the hand has to control.1 A more mundane example is leading a cup of coffee to one’s mouth to drink: the transporting hand applies a force not only to the cup but also indirectly to the liquid that acts back onto the hand and requires sensitive adjustments to avoid spilling the coffee.2–5 In the face of external perturbations or internal noise, further uncertainty and risk are added. Nevertheless, humans are strikingly adept at interacting with a large variety of such objects, seemingly even without much practice as we continuously encounter novel objects.

To date, most studies on object manipulation have focused either on grip forces needed for transporting solid objects or multi-digit grasping of a static object.6–9 It is frequently assumed that humans use inverse dynamics and acquire an internal model of their limbs, the object, and the environment.10–14 For example, Gawthrop and colleagues15,16 simulated controlling an inverted pendulum as a model for human balancing and showed that a non-predictive controller was unable to explain how humans performed the task. Alternative controllers that included predictor feedback based on an internal model were able to provide a better explanation. Predicting the dynamics of the interaction is, therefore, important since it allows one for planning the appropriate actions accordingly.

However, for physically interactive tasks with objects that exhibit nonlinear underactuated dynamics, a control strategy that fully relies on an internal model appears unlikely as exact predictions are challenging. Transporting a cup of coffee creates complex fluid dynamics that act back onto the hand and the applied force needs continuous regulation. Compensating for inadequate prediction via feedback control is insufficient as relatively long feedback delays in the neuromotor system can easily lead to instabilities. Therefore, we hypothesize that instead of learning exact models, humans learn to make their interactions more predictable such that approximate models are sufficient.

Predictability is a large umbrella term that requires operationalization to afford quantitative assessment. The computational literature developed several metrics to quantify predictability and complexity of time series, for example, mutual information and Kolmogorov entropy.17 In this study, we argue that stability of a complex system is one important avenue for the human control system to attain and enhance predictability because a stable system can reject perturbations and stay on its orbit.

Previous work has shown that humans tune into stability when the task affords such a strategy. For the task of rhythmically bouncing a ball with a racket, Sternad and colleagues demonstrated that the task afforded stability.18–21 Using model-derived criteria for stability, experimental results showed that humans adopted a strategy that made the ball trajectories stable. The benefit of this strategy is that extensive error corrections are obviated as the system returns to the stable orbit after (sufficiently small) perturbations, without necessitating control. In the face of external perturbations and human variability, a stable and robust strategy makes the system more predictable.

Numerous studies have analyzed the (in-)stability and control of an inverted pendulum as it presents a simple model for postural control as well as a challenging skill when balancing a pole on one’s fingertip. Milton et al.22 and Insperger and Milton23 used local linear stability analysis to explore the robustness of different controllers for stick balancing in the face of delay and sensory uncertainty. However, this analysis method has to assume that the system is at a stable attractor state. This may be an incorrect assumption because the authors also showed that the human controller did not stabilize the upright position, i.e., it did not make it a fixed-point attractor.22 Asai et al.24 examined postural control modeled as an inverted pendulum. They showed that intermittent control of the inverted pendulum resulted in more robustness as it created larger regions of stability. Research on human locomotion has shown that walking has stability properties consistent with our observations of how little attention is required and how easily humans compensate for small uneven support surfaces.25 Several studies used Lyapunov exponents and Floquet multipliers to demonstrate degrees of stability in healthy and clinical populations.26–28 In sum, there is a considerable number of studies on stability in human behavior, but they have been confined to only very few models that assumed that the behavior was at steady state, either a fixed-point or a closed orbit. Given the richness and complexity of human behavior, it would also be useful to assess stability during a transient point-to-point movement.

Four previous studies on transporting a cup of coffee, modeled as a cart with a suspended pendulum, provided empirical results that supported our hypothesis that humans seek to make the system dynamics predictable.3,4,29,30 Using the simplified model of a cup with a ball rolling inside, implemented in a virtual environment, subjects moved the cup without losing the ball, i.e., spilling the coffee. Predictability was operationalized by a measure of the safety or energy margin between the ball and the rim of the cup.3,4 Subjects increased the safety margin which made the performance less risky and, implicitly, more predictable. Two follow-up studies on the same task examined continuous rhythmic performance, where the nonlinear system exhibited significantly more complex behavior.29,30 Subjects moved the cup and ball at a given frequency but could choose their movement amplitude. Predictability of the object dynamics was evaluated by calculating the measure of mutual information between the applied force and the object dynamics. Results showed that subjects increased the mutual information, reflecting the enhanced predictability of the object dynamics. The same results were obtained in a complementary experiment, where subjects could choose their frequency for a given amplitude.30 

Despite these forays into assessing predictability and stability, there is a need for measures of stability that can be applied to transient behaviors and more complex dynamics. We propose that contraction analysis, a differential form of stability analysis developed for nonlinear dynamical systems, is an appropriate tool.31 Contraction analysis quantifies the convergence or divergence of trajectories towards each other. Due to the relative evaluation of trajectories, this method does not require specific knowledge of the attractors of the system or assume stability at the onset. In fact, it does not even require the trajectory to be close to an attractor. It thereby differs from the strong assumptions of the Lyapunov analysis or Floquet multipliers. Contraction regions are computed based on a mathematical model of the task and are, therefore, independent of the typical restrictions of human data, such as noise and limited sampling rate.

This study examined stability or convergence in human control of a complex object such as a cup of coffee. Despite the dynamic complexity, some prediction of the complex dynamics is needed, but a precise internal model of the sloshing coffee appears challenging. The hypothesis in this study is that humans develop control strategies that establish stability or convergence to attain more predictability in the interaction with the complex object. Specifically, we examined human interaction with a simplified model of a two-dimensional cup with a rolling ball in the presence of external perturbations.

This task poses a control challenge absent in free movements due to the nonlinear underactuated dynamics of the rolling ball that introduces bidirectional forces.32 The simplified model of the task was rendered in a virtual environment interfaced by a robotic manipulandum.33 The mathematical model of a cart and pendulum was analyzed using contraction theory. Human performance was evaluated against the calculated regions of convergence and divergence. Results show that humans indeed sought the contraction regions. We interpret these results to indicate that humans make the interactions more predictable in order to obviate error corrections based on a complex internal model.

Stability in human movements has been analyzed with several mathematical tools that have advantages but cannot be applied to a complex interactive task such as the one in focus. To better appreciate the assumptions and limitations of stability analysis and the opportunities of contraction analysis for such complex and transient control tasks, a brief summary of extant approaches to assess stability is provided.

The origins of the Lyapunov exponents method can be traced back to Lyapunov’s linearization method, sometimes referred to as Lyapunov’s first method.34 These exponents are computed in a local linear stability analysis of trajectories neighboring an attractor and they provide a measure of the exponential rate of convergence/divergence of these trajectories to the attractor.

Consider a n -dimensional autonomous dynamical system

x ˙ = f ( x ) ,
(1)

where x is a n -dimensional state vector. Let x ( t ) be a trajectory of the system, starting at an initial condition x 0 on the attractor. Then, one can consider a neighboring trajectory y ( t ) by applying a small perturbation ϵ 0 to the initial condition. The time evolution of the perturbation ϵ ( t ) is governed by the equation

d ϵ d t = J ( x 0 ) ϵ ,
(2)

where J ( x 0 ) is the n × n Jacobian matrix evaluated at x 0 . The spectrum of Lyapunov exponents of the system is then defined as

λ i = lim t 1 t ln ϵ ( t ) ϵ 0 ,
(3)

where i = 1 , , n . Since the state space is n -dimensional, n initial perturbations can be applied, hence resulting in n distinct Lyapunov exponents. If at least one of the exponents is positive, then that is an indicator of local instability, which is a necessary condition for a chaotic attractor.

As can be seen from (3), computing a Lyapunov exponent requires an infinitely long dataset. However, data collected from human experiments are always restricted in length, and frequently, one cannot observe and collect data for all the state variables. Therefore, a number of algorithms have been proposed to estimate the Lyapunov spectrum or the largest Lyapunov exponent from a finite time-series.35,36 Since Lyapunov exponents express the exponential rates at which orbits in the vicinity of an attractor converge (or diverge),37 these algorithms first reconstruct the attractor from a finite time-series using time-delayed copies of the measured state variable.38 This renders the estimation of the Lyapunov exponent sensitive to the chosen delay, the dimension in which the time-series is embedded, and the ever-present noise in the collected data. In fact, it has been shown that these algorithms can yield incorrect estimates of negative exponents.39 Finally, for methods using the largest Lyapunov exponent, it was shown that a negative largest Lyapunov exponent does not necessarily imply stability.40 This is commonly known as the Perron effect.41 

Poincaré maps and Floquet Theory have been used to assess the stability of periodic orbits, as typically encountered in locomotor movements. Hurmuzlu and Basdogan used the Poincaré maps and Floquet theory to assess the stability of human locomotion,27 while Nomura et al.42 used these tools to analyze how the stability of bipedal gaits depends on the trajectories of the joint angles. The Floquet analysis is similar to the Lyapunov exponents method in that it calculates exponents using a local linearized map around the orbit to quantify its stability. However, Floquet theory originally dealt with deterministic linear systems. As human movement is inherently nonlinear and stochastic, the validity of quantifying stability of human walking using Floquet multipliers is questionable. For example, it has been shown that Floquet analysis of a simple nonlinear walking model with stochastic noise over-estimates orbital stability.43,44 In addition, human walking is not exactly periodic,45 and it has even been suggested that it might not be a stable limit cycle.46 

Many studies on human movement have used variability as an inverse estimate of stability as it is well known that experts have less variability, i.e., assumed to be more stable. For example, in locomotion research, the variance of the step cycle has been interpreted as a measure of gait stability. Greater variability in stride-to-stride step length,47 stride-to-stride step width,48 and stride period49,50 has been associated with less stability and higher risk for falls. However, low variability does not necessarily imply stability, and in fact, the relation between the two is complex and task-dependent. A study by Brach et al.51 found that both too little and too much variability in step width were associated with falling. Dingwell and Marin50 also showed that slower walking speeds were associated with increased stability, even though stride-to-stride variability was larger. In redundant tasks where there is a manifold of solutions for achieving the task, variability only needs to be channeled into task-irrelevant dimensions.52–54 Furthermore, seemingly random variability may actually be the expression of a weak chaotic attractor.55 

Understanding coordination and control not only needs evaluation of the outcome but also needs measures of the process or the execution. For example, in targeted throwing, the errors to the target are not the only reflection of expertise. Variability of the process and the mapping of the execution variables into the lower-dimensional result space can serve as a window into the stability of the behavior.56 Several lines of study have examined variability in the redundant subspace or null space by mapping higher-dimensional execution into lower-dimensional result space. For example, the uncontrolled manifold analysis (UCM) approach evaluates the anisotropy of the variability in null space using the Jacobian and its eigenvalues as a measure of stability.57,58 However, it is important to point out that stability in the context of UCM expresses the anisotropy of the variability and is not a measure of robustness to perturbations. It is concerned with the degrees of freedom that are primarily controlled by the nervous system when performing a task. Moreover, the variability metric within this framework is sensitive to the choice of coordinates for the execution variables.59 It is also important to note that the distributional properties of variability are task-dependent and do not provide an absolute metric of stability.

For cases where a continuous dynamical model is known, one can analyze the stability properties of the equilibrium points by evaluating the eigenvalues of the Jacobian of the linearized system around the equilibrium points. This is Lyapunov’s first method. If the Jacobian is strictly stable, i.e., all eigenvalues are in the left-half of the complex plane, then the equilibrium point is asymptotically stable. If the Jacobian is unstable, i.e., at least one of the eigenvalues is in the right half of the complex plane, then the equilibrium point is unstable. Finally, if the Jacobian is marginally stable, i.e., all the eigenvalues are in the left half, but at least one of them lies on the imaginary axis, then nothing can be concluded about the equilibrium point. For discrete models, the eigenvalues must lie inside the unit circle to indicate stability.

Sternad et al.18,19,60 used a local linear stability analysis to investigate whether humans utilized principles of stability when they rhythmically bounced a ball. Stability analysis of a bouncing ball model identified conditions that indicated stable behavior. These criteria were used to evaluate human performance. Results revealed that humans indeed conformed to these conditions and tuned into the stable behavior of the task.

However, the model and the movement task was explicitly simple and low-dimensional. It is difficult to analyze more complex tasks with continuous or hybrid dynamics and larger state spaces, as local linear stability analysis of every equilibrium point would be daunting. Moreover, it is possible for model-based local linear stability analysis to be non-conclusive, whenever the Jacobian is found to be marginally stable.

To summarize, extant stability analyses require and assume the behavior to be a steady-state closed orbit, a fixed-point attractor, a limit cycle, etc. In this paper, we propose a method for assessing stability of transient point-to-point movements.

Contraction theory31 has fewer restrictions than previously described approaches, as it adopts a differential view on stability of a nonlinear system by examining the convergence or divergence of neighboring trajectories to each other. It can, therefore, be applied to a larger variety of nonlinear dynamical systems. It is important to note that contraction is a strong form of stability, since it demands an exponential rate of convergence.

Consider a nonlinear dynamical system

x ˙ = f ( x , t ) ,
(4)

where f is an n × 1 nonlinear vector function and x is the n × 1 state vector. Let us look at two neighboring trajectories, x ( t ) and y ( t ) , which are two solutions of (4) with initial conditions x 0 and y 0 , respectively. The virtual displacement is an infinitesimal displacement, δ x , between the trajectories at a fixed time, as illustrated in Fig. 1.

FIG. 1.

Virtual displacement between two neighboring trajectories.

FIG. 1.

Virtual displacement between two neighboring trajectories.

Close modal

From (4), an exact differential relation in the virtual displacement can be derived as

δ x ˙ = f x ( x , t ) δ x .
(5)

The squared distance between the trajectories as time evolves is governed by

d d t ( δ x T δ x ) = 2 δ x T f x δ x .
(6)

This means that the virtual displacement δ x converges exponentially to zero if the Jacobian f x is uniformly negative definite. Regions in the state space in which this negative definiteness property is satisfied are defined as contraction regions. For a uniformly negative definite Jacobian, the following holds

β > 0 , x , t 0 , 1 2 f x + f T x β I < 0.
(7)

This means that for the Jacobian to be uniformly negative definite, all the eigenvalues of its symmetric part must be uniformly negative definite. Note, however, that this is only a sufficient condition for exponential convergence of the virtual displacement.

A necessary and sufficient condition for exponential convergence can be formulated by using a more general definition of differential length. Let us take a continuously differentiable coordinate transformation of the form

δ z = Θ ( x , t ) δ x ,
(8)

where Θ ( x , t ) is a square matrix that satisfies Θ T Θ > 0 . Θ is referred to as the contraction metric. In this new coordinate frame, a contraction region is one that satisfies

1 2 K + K T < 0 ,
(9)

where K is the generalized Jacobian

K = Θ ˙ + Θ f x Θ 1 ,
(10)

meaning that all eigenvalues of the symmetric part of the generalized Jacobian K must be uniformly negative definite.

To summarize, for a region of the state space to be contracting, there must exist a metric Θ ( x , t ) such that Θ T Θ > 0 and 1 2 ( K + K T ) < 0 over that region. Any trajectory entering a ball of constant radius around another trajectory within a contraction region remains in that ball and converges exponentially to that trajectory.

It is important to mention that finding a suitable metric is not trivial and can be the most challenging part of contraction analysis. This study will use a method based on solving a partial differential equation in the system dynamics to arrive at a suitable metric.

To highlight the desirable attributes of contraction analysis, it is compared against the existing methods of stability analysis as applied to human movement control. As shown in the following, it is more versatile and suitable for movements that are not at steady state, such as those involving dynamically complex physical interactions.

1. Computations are precise

Computing contraction regions is solely based on calculating the eigenvalues of the Jacobian of the dynamical model. One can, therefore, establish with exact precision and certainty whether or not a region of the state space is contracting. Even if the contraction metric is found numerically, the contraction regions with respect to that metric can be computed exactly. In contrast, the Lyapunov exponent can only be estimated using numerical methods and algorithms since, mathematically speaking, Lyapunov exponents are only precise for infinite time series. The performance of these methods depends on the delay and the embedding dimension. Numerical stability and sensitivity are, therefore, a big challenge for an accurate estimation of the true Lyapunov exponent.

2. Stability estimates do not require a reference

Stability is typically referenced with respect to a nominal behavior; for instance, “stability of an equilibrium/fixed point,” “stability of a periodic orbit,” “rate of convergence to an attractor,” etc. But what if the observed behavior is not in the vicinity of an attractor? What if there is no prior knowledge of the attractor? This is the likely scenario for dexterous behavior with dynamically complex physical interactions. It is, therefore, useful to adopt a differential view of stability by looking at whether or not nearby trajectories converge to each other, irrespective of what behavior they converge to and irrespective of whether they are close to an attractor. This is what contraction analysis does. “Reconstructing an attractor” is not a requirement for the analysis.

3. Analysis does not rely on data that can be noisy

The algorithms used to calculate Lyapunov exponents are sensitive to noise in the data and can, therefore, result in wrong approximations of the Lyapunov exponent. In contrast, computing a contraction region is based on a model and does not require any data input. The only limitation is that the human behavior then needs to be sufficiently close to that model.

4. Analysis does not require linearization

All existing methods for stability analysis require linearization. For movements involving physical interaction and perturbing forces, linearization is often not appropriate. Contraction analysis is based on the exact differential form (5) and does not require any approximations or linearizations, as it deals with the full model and its nonlinearities. It is, therefore, more versatile and suitable for tasks involving interaction forces and perturbations. It is important to note that while (5) resembles a linearization, it is in fact an exact mathematical definition of the rate of change of virtual displacements between neighboring trajectories at a fixed time.

5. Analysis can be applied to time-varying systems

For model-based local linear stability analysis, existing theorems are only applicable for autonomous dynamical systems. As for non-autonomous systems, one can only deduce stability of the nonlinear system if its linearization is uniformly asymptotically stable. Otherwise, no conclusions can be drawn about the original nonlinear system. Contraction analysis does not distinguish between autonomous and non-autonomous systems, since convergence or divergence of trajectories in a particular region of state space is dictated by the eigenvalues of the Jacobian matrix, regardless of whether or not it is a function of time.

The action of carrying a cup of coffee is an everyday example of a physical interaction with a dynamically complex object. Moving the cup can cause sloshing of the coffee, which in turn, exerts forces on the cup and the arm. Given the complex internal dynamics and interaction forces, the question is how humans control their interactions to ensure safe transport. We hypothesize that humans create stable trajectories so that perturbations have relatively little effect and require little error correction.61 To examine this hypothesis, we tested human performance in a virtual experiment that implemented a simplified model of a cup of coffee. Human kinematic and kinetic data were collected when interacting with this model system.

In the experiment, human subjects transported a cup with a ball rolling inside and moved it from an initial position to a final position. Subjects were instructed to transport the cup and ball as fast as possible without losing the ball, i.e., spilling the coffee. To create an additional challenge, a visible small hurdle on the line exerted a force either against the direction of movement (resistive) or in the direction of movement (assistive) of the cup. The magnitude and direction of the perturbation were kept invariant across each block of trials so that subjects could learn how best to prepare for the perturbation. The hypothesis was that humans exploited contraction regions to accommodate for these perturbations. We expected that subjects would develop specific and distinct strategies for the two types of perturbation.

For the experiment, the realistic 3-dimensional cup with the fluid dynamics of the sloshing coffee was reduced to a 2-dimensional arc and a ball rolling inside the arc. This system was confined to move on a horizontal line. Under the assumption that the ball does not roll and only slides without friction along the cup, the dynamics of the system become mathematically equivalent to the dynamics of the well-known cart-pendulum system. Even though this model vastly simplifies the full dynamics, the essential challenges of the task were retained: underactuation and nonlinearity, where the latter could induce chaotic behavior. Figure 2(a) illustrates the real and simplified task; Fig. 2(b) shows the dynamic model of the simplified task.

FIG. 2.

(a) The real task and the conceptual model. (b) The dynamical model. (c) The virtual setup. The subject interacted with the object via a robotic manipulandum and felt the perturbation when the cup passed over the visible bump, although the cup stayed on the horizontal line. (d) The display that the subject saw on the screen. The distance between the start and target box was 0.4  m.

FIG. 2.

(a) The real task and the conceptual model. (b) The dynamical model. (c) The virtual setup. The subject interacted with the object via a robotic manipulandum and felt the perturbation when the cup passed over the visible bump, although the cup stayed on the horizontal line. (d) The display that the subject saw on the screen. The distance between the start and target box was 0.4  m.

Close modal

The dynamics of the simplified model are

( m + M ) x ¨ ( t ) = l m ϕ ˙ ( t ) 2 sin [ ϕ ( t ) ] ϕ ¨ ( t ) cos [ ϕ ( t ) ] + F ( t ) b x ˙ ( t ) ,
(11)
l ϕ ¨ ( t ) = g sin [ ϕ ( t ) ] G x ¨ ( t ) cos [ ϕ ( t ) ] ,
(12)

where x ( t ) denotes the horizontal position of the cart, ϕ ( t ) is the pendulum angle measured counter-clockwise from the vertical, m is the mass of the pendulum, M is the mass of the cart, l is the length of the massless pendulum rod, and g is the gravitational acceleration. The force applied by the human subject is F ( t ) . To ensure the existence of contraction regions, the model should include some form of energy dissipation from the system. This was achieved by adding damping in the cup movement direction x , where the damping coefficient is denoted by b . Physically, this damping may arise from the human arm impedance or friction on the surface. To make the task more challenging and to eliminate the trivial case where the entire state space is a contraction region, the simulated dynamics in the virtual environment incorporated a gain G for the cup acceleration x ¨ . This gain made the ball movement more sensitive to the forces applied to the cup. For the simulations and the experiment in the virtual environment, the model parameters were fixed: M = 3.5  kg, m = 0.3  kg, l = 0.35  m, b = 20  N s/m, and G = 5 .

This simplified model of the task dynamics was rendered in a virtual environment using a robotic manipulandum interfacing with the virtual object that provided subjects with haptic feedback. Figure 2(c) illustrates the subject interacting with the virtual object. A projection screen displayed the cup (which corresponded to the cart) and the ball (which corresponded to the pendulum bob), as depicted in Fig. 2(d).

Subjects were seated about 2  m in front of a large backprojection screen ( 2.4 × 2.4  m). They physically interacted with the object via a 3 -degree-of-freedom robotic manipulandum (HapticMaster, Motekforce, NL).33 By applying a force to the handle of the robotic arm, participants controlled the horizontal x position of the virtual cup. The robotic arm was restricted to move only in the horizontal direction along the subject’s frontal plane to ensure a unidirectional motion of the cup. The robotic arm provided haptic feedback, allowing participants to sense the object’s inertia, the force of the ball on the cup, and the perturbations. The force applied by the participants to the manipulandum F and the kinematics of the cup and the ball x , x ˙ , x ¨ , ϕ , ϕ ˙ , ϕ ¨ were recorded at 120  Hz.

Participants were instructed to move the cup from the start box on the left to the target box on the right as fast as possible, while ensuring that the ball did not exceed the rim of the cup and “escape.” A perturbation of magnitude 40  N, either in the direction of motion of the cup (assistive) or against it (resistive) was applied at 60 % of the travel distance. The duration of this perturbation was 20  ms. The position of the perturbation was visually displayed as a bump [Fig. 2(d)]. The subject knew the magnitude and direction of the perturbation which give them the opportunity to learn how to navigate this perturbation. For simplicity, the virtual cup remained on the horizontal line and moved through the bump. All subjects gave informed written consent before the experiment. The protocol was approved by the Institutional Review Board of Northeastern University.

The experiment consisted of 4 blocks. Block 1 comprised 60 trials without any perturbation to allow subjects to familiarize themselves with the task. Blocks 2 and 4 comprised 60 trials each and involved a series of assistive or resistive perturbations. The order of resistive and assistive perturbations was randomized. Block 3 presented 10 unperturbed trials to separate the conditions of the two experimental blocks. At the beginning of each trial, the cup was centered in the start box and the ball rested at its equilibrium position at the bottom of the cup.

The first step in computing the contraction regions of the models (11) and (12) was to re-write the equations in state space form. The state vector was taken to be X = ( x ˙ , ϕ , ϕ ˙ ) T = ( x 1 , x 2 , x 3 ) T , and the resulting state space equations are

X ˙ = x ¨ ϕ ˙ ϕ ¨ = x ˙ 1 x ˙ 2 x ˙ 3 = f ( X ) = b x 1 + F + g m sin ( x 2 ) cos ( x 2 ) + l m x 3 2 sin ( x 2 ) G m cos 2 ( x 2 ) + m + M x 3 G cos ( x 2 ) F b x 1 + g ( m + M ) sin ( x 2 ) + G l m x 3 2 sin ( x 2 ) cos ( x 2 ) G l m cos 2 ( x 2 ) l ( m + M ) .
(13)

The Jacobian of the free unforced system ( F = 0 ) is found to be

J = 1 0.08 cos 2 ( x 2 ) 0.2 α 0.1 sin ( x 2 ) x 3 cos 2 ( x 2 ) 2.5 0 0 1 1 0.01 sec ( x 2 ) 0.005 cos ( x 2 ) β 2 sin ( 2 x 2 ) x 3 cos ( 2 x 2 ) 4.07 ,
(14)

where

α = cos ( x 2 ) [ 0.04 cos ( 2 x 2 ) + 0.07 ] x 3 2 + 4 cos ( 2 x 2 ) + 13.3 sin ( 2 x 2 ) x 1 1 2.5 cos 2 ( x 2 ) 2
(15)

and

β = γ 91.1 cos ( x 2 ) 17.8 cos ( 3 x 2 ) + [ 530.2 sin ( x 2 ) 47.6 sin ( 3 x 2 ) ] x 1 2.5 cos 2 ( x 2 ) 2 ,
(16)

where

γ = [ 2.03 cos ( 2 x 2 ) 5.6 × 10 17 cos ( 4 x 2 ) + 0.5 ] x 3 2 .
(17)

A contraction region is defined where the symmetric part of the Jacobian is uniformly negative definite. For this Jacobian, no region in the state space was found such that the eigenvalues of its symmetric part J s y m were uniformly negative definite. However, this did not rule out the existence of contraction regions, since negativity of J s y m is only a sufficient condition.

The next step was to find a contraction metric Θ ( X , t ) for which some regions of the state space were contracting. To find a suitable metric that revealed the contraction regions, the following partial differential equation was used31 

Θ X f + Θ J = Θ ,
(18)

where f was given in (13) and J was the Jacobian given by (14). This partial differential equation was solved numerically to obtain the contraction metric, which then enabled the computation of the generalized Jacobian K from (10). To deduce the contraction regions, the negativity condition (9) was tested for points in the state space within the range

0.2 x ˙ 0.7 ; 1.5 ϕ 1.5 ; 6 ϕ ˙ 6 .

These boundaries were used since human subjects’ movements were confined to this range. The points in the state space that satisfied these condition were therefore contained in a contraction region.

We hypothesized that humans employed control strategies that exploited contraction regions to diminish the effects of perturbations. To test this hypothesis, the experimental trajectories collected from four subjects were evaluated with respect to the contraction regions.

Table I compares the average times required to complete the task for 10 early and 10 late trials of each block, along with one standard deviation. As to be expected, subjects decreased the average time and their standard deviations as they became more proficient in the late block. With very few exceptions ( 1 out of 4 subjects during resistive perturbations), this improvement was seen in all subjects and all conditions.

Table I.

Average task completion times and standard deviations.

Subject 1 Subject 2 Subject 3 Subject 4
Baseline  Early  1.92 s ± 0.45   1.77 s ± 0.45   3.05 s ± 0.42   1.88 s ± 0.47  
  Late  1.58 s ± 0.08   1.30 s ± 0.29   2.13 s ± 0.36   1.43 s ± 0.08  
Assistive  Early  2.09 s ± 0.45   1.56 s ± 0.53   2.18 s ± 0.45   1.61 s ± 0.46  
  Late  1.68 s ± 0.20   1.35 s ± 0.20   1.98 s ± 0.36   1.55 s ± 0.33  
Resistive  Early  2.35 s ± 0.12   1.32 s ± 0.27   3.09 s ± 0.80   1.71 s ± 0.50  
  Late  2.50 s ± 0.23   1.33 s ± 0.22   2.82 s ± 0.51   1.29 s ± 0.14  
Subject 1 Subject 2 Subject 3 Subject 4
Baseline  Early  1.92 s ± 0.45   1.77 s ± 0.45   3.05 s ± 0.42   1.88 s ± 0.47  
  Late  1.58 s ± 0.08   1.30 s ± 0.29   2.13 s ± 0.36   1.43 s ± 0.08  
Assistive  Early  2.09 s ± 0.45   1.56 s ± 0.53   2.18 s ± 0.45   1.61 s ± 0.46  
  Late  1.68 s ± 0.20   1.35 s ± 0.20   1.98 s ± 0.36   1.55 s ± 0.33  
Resistive  Early  2.35 s ± 0.12   1.32 s ± 0.27   3.09 s ± 0.80   1.71 s ± 0.50  
  Late  2.50 s ± 0.23   1.33 s ± 0.22   2.82 s ± 0.51   1.29 s ± 0.14  

The kinematic data from one representative subject is presented in Fig. 3. The left column displays the cup velocity over the cup position for the three blocks; the second column presents the ball’s angular velocity over the cup position for the same three blocks. The solid lines are the average of the first 10 trials and the last 10 trials; the shaded band represents one standard deviation. Comparing the early with late trials in the baseline and the two perturbation blocks makes apparent that both cup and ball velocities exhibited a decrease in variability. In the two perturbation blocks, the variability of the cup and ball velocities was reduced at the moment preceding the perturbation in the later trials. The two perturbation types showed visible and different effects on the velocity of both cup and ball. As expected, the assistive perturbation increased the velocity transiently, while the resistive perturbation decreased the velocity transiently. The profiles of the cup trajectories also showed a change in later trials.

FIG. 3.

Kinematic data for one of the four subjects. The left column displays cup velocity versus cup position for the three different blocks of the experiment. The right column plots ball angular velocity versus cup position. The vertical line indicates the position of the perturbation. The solid line represents the average of the first 10 and last 10 trials of the respective block; the shaded band shows one standard deviation. The different average trajectories and the reduced variability in the late trials indicate learning. The markedly different trajectories in the resistive and assistive perturbations show that subjects developed different strategies for the two perturbations.

FIG. 3.

Kinematic data for one of the four subjects. The left column displays cup velocity versus cup position for the three different blocks of the experiment. The right column plots ball angular velocity versus cup position. The vertical line indicates the position of the perturbation. The solid line represents the average of the first 10 and last 10 trials of the respective block; the shaded band shows one standard deviation. The different average trajectories and the reduced variability in the late trials indicate learning. The markedly different trajectories in the resistive and assistive perturbations show that subjects developed different strategies for the two perturbations.

Close modal

To evaluate whether subjects learnt to exploit the contraction regions, the data were analyzed against the calculated contraction regions.

1. Assistive perturbations

To test the hypothesis that subjects exploited the contraction regions in state space, the trajectories of cup and ball were plotted in 3-dimensional state space. The calculated regions of convergence are displayed as yellow-shaded volumes in Fig. 4. The point P denotes the instant just before the perturbation while P + denotes the instant just after the perturbation. Note that even though the x -position of the perturbation remained fixed in the work space for all trials, the positions of P and P + need not be fixed in state space, which is spanned by ϕ , ϕ ˙ , and x ˙ . Figure 4(a) illustrates one of the early trials for one subject, while Fig. 4(b) displays one late trial, both with an assistive perturbation. In the early trial, the trajectory did not pass through any contraction region. However, with practice, the subject learnt to enter the contraction region just before the onset of the perturbation. This caused the perturbation to occur within the contraction region, thereby mitigating the destabilizing effect of the perturbation. The same strategy was also observed in the three other subjects.

FIG. 4.

Human trajectories when negotiating the assisting perturbation: (a) Early trial: the subject did not exploit any contraction region. (b) Late trial: the subject encountered the perturbation entirely inside a contraction region. The black arrowhead indicates the starting point of the trajectory. P denotes the instant just before the perturbation, while P + denotes the instant after the perturbation. The exact ellipse indicates pendulum oscillations for zero cup velocity at the target box.

FIG. 4.

Human trajectories when negotiating the assisting perturbation: (a) Early trial: the subject did not exploit any contraction region. (b) Late trial: the subject encountered the perturbation entirely inside a contraction region. The black arrowhead indicates the starting point of the trajectory. P denotes the instant just before the perturbation, while P + denotes the instant after the perturbation. The exact ellipse indicates pendulum oscillations for zero cup velocity at the target box.

Close modal

2. Resistive perturbations

When the perturbations were resistive, subjects used a different strategy [Fig. 5(a)]. In the early trials, subjects did not approach any contraction region, similar to the assistive perturbations. For the trial displayed in Fig. 5(a), the subject failed to keep the ball in the cup, hence the trajectory ended shortly after P + . However, as the subject learnt to navigate the perturbations, they chose the perturbation onset such that the subsequent segment of their trajectory passed through a contraction region [Fig. 5(b)]. This attenuated the transient effects of the perturbation. This pattern was again consistent in the three other subjects. These observations supported the hypothesis that humans sought convergent regions to stabilize their trajectories against perturbations.

FIG. 5.

Human trajectories for the condition with resisting perturbation: (a) Early trial: in this trial, the subject dropped the ball from the cup, hence the trajectory ends just after P + . (b) Late trial: the subsequent trajectory after the perturbation entered a contraction region. The black arrowhead indicates the starting point of the trajectory. P denotes the instant just before the perturbation, while P + denotes the instant after the perturbation. The exact ellipse indicates that the cup has reached the target box at rest, while the pendulum was still oscillating.

FIG. 5.

Human trajectories for the condition with resisting perturbation: (a) Early trial: in this trial, the subject dropped the ball from the cup, hence the trajectory ends just after P + . (b) Late trial: the subsequent trajectory after the perturbation entered a contraction region. The black arrowhead indicates the starting point of the trajectory. P denotes the instant just before the perturbation, while P + denotes the instant after the perturbation. The exact ellipse indicates that the cup has reached the target box at rest, while the pendulum was still oscillating.

Close modal

Augmenting human capabilities with tools is fundamental to our daily existence and is also a central research objective in developing robotic manipulators and augmenting devices.62,63 Yet to date, there is relatively little understanding of how humans achieve their exquisite skill in handling the vast variety of objects that extend their own bodies’ capabilities. To shed light on human control, more analysis techniques are needed to reveal the intricate interaction with such complex objects and their unique challenges. It is generally assumed that human control relies on internal models of their own body, the objects, and their interaction with the environment.10–14 However, these internal models are complex and prediction of the resulting behavior can become hard, if not impossible.1,29,61,64 Given the significant delays and ubiquitous noise in the nervous system, even small errors may amplify and severely disrupt and destabilize the interaction. We, therefore, reasoned that humans learn to interact with objects by discovering and exploiting stability properties of the system.18,61

However, stability analysis itself becomes challenging (for the scientist) when the system is far from an attractor, which makes most current analysis methods inapplicable. Following a brief review of existent stability analyses as applied to human movement, we suggested to use contraction analysis, which does not require the system to be close to an equilibrium state. Using the well-known control task of transporting a cart and a suspended pendulum, mimicking a cup of sloshing coffee, we demonstrated how contraction analysis can serve as a useful method to analyze human behavior. We submit that compared to existing methods for stability assessment, specifically in human motor control research, contraction analysis is a valuable technique to study complex interactions.

Based on a mathematical model of the task dynamics, contraction analysis identifies regions in state space where contraction is known to occur. This characterization then serves as a basis to evaluate human trajectories in state space and test whether humans are sensitive to these properties of the task dynamics. The current study revealed that humans indeed developed specific strategies to stabilize their performance. When facing assistive perturbations, the cup trajectories entered a contraction region and stayed inside throughout the duration of the perturbation. Hence, it can be inferred that the convergent dynamics assisted in absorbing the perturbation. In contrast, for resistive perturbations subjects accessed a contraction region only after the perturbation. These observations indicated that, depending on the perturbation type, subjects developed distinctive strategies that exploited the dynamic properties of the task. It is important to note that it is not clear whether these “prediction-enhancing” control strategies are also the most stable or the most energy-efficient ones or whether humans tradeoff these different criteria. However, two earlier studies designed experimental conditions that made these criteria mutually exclusive. The results showed that subjects did not minimize force but rather maximized predictability measured by mutual information.29,30

This study used a greatly simplified model system as proxy for the challenge of transporting an object with internal dynamics. However, it is important to point out that contraction analysis is not limited to such low-dimensional, deterministic, and accurate models. High-dimensional models can likewise be analyzed within this framework. Contraction analysis can also be applied to stochastic dynamical systems.65 Moreover, even if the formulation of an accurate model may become difficult and parameter values become uncertain, contraction analysis is still applicable. In contrast to the Lyapunov-based methods, contraction analysis does not require tracking of the location of the steady-state equilibria to study the robustness to parameter perturbations and uncertainties. For example, Aylward et al.66 used sums-of-squares programming methods to find bounds on the maximum allowable uncertainty on the parameters such that a given system remained contracting with respect to the same contraction metric. For the present task, an interesting next step could be to perform a robustness analysis to examine how the contraction regions change with parameter uncertainty and test whether or not the trajectories still passed through the “new” contraction regions.

It is necessary to point out that contraction represents a strong form of stability as it implies an exponential rate of convergence.67 While Lyapunov exponents and Floquet multipliers also quantify exponential stability, the model-based local linear stability analysis only identifies asymptotic stability. It may be argued that humans need not exhibit exponential convergence and this requirement may be too conservative. On the other hand, if such a behavior is identified, this can be interpreted as a strong result.

While contraction theory provides numerous advantages over existing methods for stability assessment of human movement, it also suffers from some limitations. One such limitation is that the analysis requires a mathematical model, whether precise or including parameter uncertainties, to proceed with the analysis. This might not be possible for certain tasks. The present study implemented the simplified task model in a virtual environment which had the advantage that it exactly defines the environment and exactly measured the human input. If the task were examined in a real environment, additional effects such as rolling friction would necessarily arise. Further, a full 3D version of the task has more degrees of freedom with interactions that may become hard to adequately model. To maintain theoretical stringency, experimental control and reduced complexity are, therefore, important.

Another shortcoming of contraction analysis is that it requires the specification of a metric to compute the contraction regions with respect to that metric. This is not a trivial task and is similar in challenge to constructing a Lyapunov function, where there is no systematic method for finding it.68 This study was able to solve a partial differential equation (18) to arrive at a contraction metric. However, the existence of solutions to this partial differential equation is not guaranteed since it depends on the dynamical system under consideration. Recently, inspired by the progress in sums-of-squares optimization to compute Lyapunov functions,69 semi-definite programming70 and sums-of-squares programming71 have been used to search for contraction metrics. While beyond the scope of this paper, it may be interesting to compare the metrics that these different approaches may identify.

The choice of metric is also important for control design in robotics. Recent work has shown that the existence of a specific metric, known as the “control contraction metric” (CCM), for a system implies the existence of a universally stabilizing controller.72,73 In a separate study on sychronization in biological systems, Aminzare74 derived metrics induced by non- L 2 norms, that proved relevant for contraction analysis of a cellular biological process. The CCM and non- L 2 induced metrics could serve as an initial guess for a suitable metric.

One final limitation of contraction analysis is that the choice of the metric can affect the conclusions made about the subjects’ behavior. Different metrics may reveal different contraction regions.67 One might initially find a metric with contraction regions that did not intersect the human trajectories, hence conclude that contraction was not exploited. However, because contraction is both a necessary and sufficient condition for exponential convergence, once a region has been established as contracting with respect to some metric, other metrics might expand that contraction region but cannot reduce it. The fact that the metric can be freely chosen could also open an avenue to identify what metric and what coordinates humans use to represent movement.59,75 Turning the hypothesis around and assuming that humans exploit contraction regions, several metrics may be determined and contrasted. The one that humans chose might reflect their intrinsic reference metric.

We proposed contraction theory as a potential framework for assessing stability of human movement in complex tasks. The model task was a physical interaction with a dynamically complex object, specifically moving a cup containing a rolling ball through a perturbation. This task is not steady in nature but rather transient. Contraction analysis is suitable for movements that begin and terminate far from an equilibrium or periodic steady state. The method, therefore, provides a new promising tool to examine human movement.

The analysis gave first evidence that humans indeed exploited contraction regions to accommodate perturbations during dynamically complex physical interactions. By traversing contraction regions, subjects mitigated the effects of perturbations and avoided the potentially chaotic evolution of the nonlinear dynamical system. This strategy overcame control delays and noise and rendered task performance stable and more predictable. More predictable dynamics may be more feasible for developing an internal model of the task.

D. Sternad was supported by the National Institutes of Health Grant Nos. R01-HD-087089, R01-HD-081346, and R21-DC-013095 and the National Science Foundation Grant Nos. NSF-NRI 1637854 and NSF-EAGER-1548514. N. Hogan was supported by the NIH Grant No. R01-HD-087089, the National Science Foundation Grant Nos. NSF-NRI 1637814 and NSF-EAGER-1548501, and the Eric P. and Evelyn E. Newman Foundation. J. Ebert was supported by a Department of Energy Computational Science Graduate Fellowship (DOE CSGF). No conflicts of interest, financial or otherwise, are declared by the authors.

1.
N.
Hogan
and
D.
Sternad
, “
Dynamic primitives of motor behavior
,”
Biol. Cybern.
106
,
1
13
(
2012
).
2.
H. C.
Mayer
and
R.
Krechetnikov
, “
Walking with coffee: Why does it spill?
,”
Phys. Rev. E
85
(
4
),
046117
(
2012
).
3.
C. J.
Hasson
,
T.
Shen
, and
D.
Sternad
, “
Energy margins in dynamic object manipulation
,”
J. Neurophysiol.
108
(
5
),
1349
1365
(
2012
).
4.
C. J.
Hasson
and
D.
Sternad
, “
Safety margins in older adults increase with improved control of a dynamic object
,”
Front. Aging Neurosci.
6
,
158
(
2014
).
5.
J.
Han
, “
A study on the coffee spilling phenomena in the low impulse regime
,”
Achiev. Life Sci.
10
(
1
),
87
101
(
2016
).
6.
J. R.
Flanagan
,
J.
Tresilian
, and
A. M.
Wing
, “
Coupling of grip force and load force during arm movements with grasped objects
,”
Neurosci. Lett.
152
(
1
),
53
56
(
1993
).
7.
J. R.
Flanagan
and
A. M.
Wing
, “
The role of internal models in motion planning and control: Evidence from grip force adjustments during movements of hand-held loads
,”
J. Neurosci.
17
(
4
),
1519
1528
(
1997
).
8.
M.
Santello
,
M.
Flanders
, and
J. F.
Soechting
, “
Postural hand synergies for tool use
,”
J. Neurosci.
18
(
23
),
10105
10115
(
1998
).
9.
Q.
Fu
and
M.
Santello
, “
Coordination between digit forces and positions: Interactions between anticipatory and feedback control
,”
J. Neurophysiol.
111
(
7
),
1519
1528
(
2014
).
10.
M.
Kawato
, “
Internal models for motor control and trajectory planning
,”
Curr. Opin. Neurobiol.
9
(
6
),
718
727
(
1999
)
11.
B.
Mehta
and
S.
Schaal
, “
Forward models in visuomotor control
,”
J. Neurophysiol.
88
(
2
),
942
953
(
2002
)
12.
C.
Takahashi
,
R. A.
Scheidt
, and
D.
Reinkensmeyer
, “
Impedance control and internal model formation when reaching in a randomly varying dynamical environment
,”
J. Neurophysiol.
86
(
2
),
1047
1051
(
2001
).
13.
J. R.
Flanagan
,
E.
Nakano
,
H.
Imamizu
,
R.
Osu
,
T.
Yoshioka
, and
M.
Kawato
, “
Composition and decomposition of internal models in motor learning under altered kinematic and dynamic environments
,”
J. Neurosci.
19
(
20
),
RC34
RC34
(
1999
).
14.
J. R.
Flanagan
,
P.
Vetter
,
R. S.
Johansson
, and
D. M.
Wolpert
, “
Prediction precedes control in motor learning
,”
Curr. Biol.
13
(
2
),
146
150
(
2003
).
15.
P.
Gawthrop
,
I.
Loram
, and
M.
Lakie
, “
Predictive feedback in human simulated pendulum balancing
,”
Biol. Cybern.
101
(
2
),
131
146
(
2009
).
16.
P.
Gawthrop
,
I.
Loram
,
H.
Gollee
, and
M.
Lakie
, “
Intermittent control models of human standing: Similarities and differences
,”
Biol. Cybern.
108
(
2
),
159
168
(
2014
).
17.
T. M.
Cover
and
J. A.
Thomas
,
Elements of Information Theory
(
John Wiley & Sons
,
2012
).
18.
S.
Schaal
,
C. G.
Atkeson
, and
D.
Sternad
, “
One-handed juggling: A dynamical approach to a rhythmic movement task
,”
J. Mot. Behav.
28
(
2
),
165
183
(
1996
).
19.
D.
Sternad
,
M.
Duarte
,
H.
Katsumata
, and
S.
Schaal
, “
Dynamics of a bouncing ball in human performance
,”
Phys. Rev. E
63
(
1
),
011902
(
2000
).
20.
D.
Sternad
,
M.
Duarte
,
H.
Katsumata
, and
S.
Schaal
, “
Bouncing a ball: Tuning into dynamic stability
,”
J. Exp. Psychol. Hum. Percept. Perform.
27
(
5
),
1163
(
2001
).
21.
R.
Ronsse
and
D.
Sternad
, “
Bouncing between model and data: Stability, passivity, and optimality in hybrid dynamics
,”
J. Mot. Behav.
42
(
6
),
389
399
(
2010
).
22.
J.
Milton
,
J. L.
Cabrera
,
T.
Ohira
,
S.
Tajima
,
Y.
Tonosaki
,
C. W.
Eurich
, and
S. A.
Campbell
, “
The time-delayed inverted pendulum: Implications for human balance control
,”
Chaos Interdiscip. J. Nonlinear Sci.
19
(
2
),
026110
(
2009
).
23.
T.
Insperger
and
J.
Milton
, “
Sensory uncertainty and stick balancing at the fingertip
,”
Biol. Cybern.
108
(
1
),
85
101
(
2014
).
24.
Y.
Asai
,
Y.
Tasaka
,
K.
Nomura
,
T.
Nomura
,
M.
Casadio
, and
P.
Morasso
, “
A model of postural control in quiet standing: Robust compensation of delay-induced instability using intermittent activation of feedback control
,”
PLoS One
4
(
7
),
e6169
(
2009
).
25.
S.
Bruijn
,
O.
Meijer
,
P.
Beek
, and
J.
Van Dieën
, “
Assessing the stability of human locomotion: A review of current measures
,”
J. R. Soc. Interface
10
(
83
),
20120999
(
2013
).
26.
J.
Dingwell
,
J.
Cusumano
,
P.
Cavanagh
, and
D.
Sternad
, “
Local dynamic stability versus kinematic variability of continuous overground and treadmill walking
,”
J. Biomech. Eng.
123
(
1
),
27
32
(
2001
).
27.
Y.
Hurmuzlu
and
C.
Basdogan
, “
On the measurement of dynamic stability of human locomotion
,”
J. Biomech. Eng.
116
(
1
),
30
36
(
1994
).
28.
Y.
Hurmuzlu
,
C.
Basdogan
, and
D.
Stoianovici
, “
Kinematics and dynamic stability of the locomotion of post-polio patients
,”
J. Biomech. Eng.
118
,
405
411
(
1996
).
29.
B.
Nasseroleslami
,
C. J.
Hasson
, and
D.
Sternad
, “
Rhythmic manipulation of objects with complex dynamics: Predictability over chaos
,”
PLoS Comput. Biol.
10
(
10
),
e1003900
(
2014
).
30.
P.
Maurice
,
N.
Hogan
, and
D.
Sternad
, “
Predictability, force and (anti-) resonance in complex object control
,”
J. Neurophysiol.
120
,
765
780
(
2018
).
31.
W.
Lohmiller
and
J.-J. E.
Slotine
, “
On contraction analysis for non-linear systems
,”
Automatica
34
(
6
),
683
696
(
1998
).
32.
N.
Hogan
, “
Impedance control: An approach to manipulation
,”
J. Dyn. Syst. Meas. Control
107
(
1
),
1
7
(
1985
).
33.
R.
Van der Linde
and
P.
Lammertse
, “
Hapticmaster—A generic force controlled robot for human interaction
,”
Ind. Rob. Int. J.
30
(
6
),
515
524
(
2003
).
34.
A. M.
Lyapunov
, “
The general problem of the stability of motion
,”
Int. J. Control
55
(
3
),
531
534
(
1992
).
35.
P.
Grassberger
and
I.
Procaccia
, “
Characterization of strange attractors
,”
Phys. Rev. Lett.
50
(
5
),
346
(
1983
).
36.
M. T.
Rosenstein
,
J. J.
Collins
, and
C. J.
De Luca
, “
A practical method for calculating largest Lyapunov exponents from small sets
,”
Phys. D Nonlinear Phenom.
65
(
1-2
),
117
134
(
1993
).
37.
J. B.
Dingwell
, “Lyapunov exponents,” in
Wiley Encyclopedia of Biomedical Engineering
(Wiley Online Library, 2006).
38.
F.
Takens
et al.,“
Detecting strange attractors in turbulence
,”
Lect. Notes Math.
898
(
1
),
366
381
(
1981
).
39.
C.
Yang
and
Q.
Wu
, “
On stability analysis via Lyapunov exponents calculated from a time series using nonlinear mapping—a case study
,”
Nonlinear Dyn.
59
(
1
),
239
257
(
2010
).
40.
N. V.
Kuznetsov
and
G.
Leonov
, “
On stability by the first approximation for discrete systems
,” in
Proceedings of the IEEE International Conference on Physics and Control
(
2005
), pp.
596
599
.
41.
O.
Perron
, “
Die Stabilitätsfrage bei Differentialgleichungen
,”
Math. Z.
32
(
1
),
703
728
(
1930
).
42.
T.
Nomura
,
K.
Kawa
,
Y.
Suzuki
,
M.
Nakanishi
, and
T.
Yamasaki
, “
Dynamic stability and phase resetting during biped gait
,”
Chaos Interdiscip. J. Nonlinear Sci.
19
(
2
),
026103
(
2009
).
43.
J.
Ahn
and
N.
Hogan
, “Is estimation of Floquet multipliers of human walking valid?,” in
40th Annual Northeast Bioengineering Conference (NEBEC)
(2014), pp. 1–2.
44.
J.
Ahn
and
N.
Hogan
, “
Improved assessment of orbital stability of rhythmic motion with noise
,”
PLoS One
10
(
3
),
e0119596
(
2015
).
45.
J. B.
Dingwell
and
J. P.
Cusumano
, “
Nonlinear time series analysis of normal and pathological human walking
,”
Chaos: Interdiscip. J. Nonlinear Sci.
10
(
4
),
848
863
(
2000
).
46.
T.
Nomura
,
Y.
Suzuki
,
C.
Fu
,
N.
Yoshikawa
,
K.
Kiyono
,
M.
Casadio
, and
P.
Morasso
, “Stability and flexibility during human motor control,” in
Advances in Cognitive Neurodynamics (V)
(Springer, 2016), pp. 67–73.
47.
B. E.
Maki
, “
Gait changes in older adults: Predictors of falls or indicators of fear?
,”
J. Am. Geriatr. Soc.
45
(
3
),
313
320
(
1997
).
48.
J. C.
Dean
,
N. B.
Alexander
, and
A. D.
Kuo
, “
The effect of lateral stabilization on walking in young and old adults
,”
IEEE Trans. Biomed. Eng.
54
(
11
),
1919
1926
(
2007
).
49.
J. M.
Hausdorff
, “
Gait dynamics, fractals and falls: Finding meaning in the stride-to-stride fluctuations of human walking
,”
Hum. Mov. Sci.
26
(
4
),
555
589
(
2007
).
50.
J. B.
Dingwell
and
L. C.
Marin
, “
Kinematic variability and local dynamic stability of upper body motions when walking at different speeds
,”
J. Biomech.
39
(
3
),
444
452
(
2006
).
51.
J. S.
Brach
,
J. E.
Berlin
,
J. M.
Van Swearingen
,
A. B.
Newman
, and
S. A.
Studenski
, “
Too much or too little step width variability is associated with a fall history in older persons who walk at or near normal gait speed
,”
J. Neuroeng. Rehabil.
2
(
1
),
21
(
2005
).
52.
H.
Müller
and
D.
Sternad
, “
Motor learning: Changes in the structure of variability in a redundant task
,” in
Progress in Motor Control
(Springer,
2009
), pp.
439
456
.
53.
R.
Cohen
and
D.
Sternad
, “
Variability in motor learning: Relocating, channeling and reducing noise
,”
Exp. Brain Res.
193
(
1
),
69
83
(
2009
).
54.
D.
Sternad
, “
It’s not (only) the mean that matters: Variability, noise and exploration in skill learning
,”
Curr. Opin. Behav. Sci.
20
,
183
195
(
2018
).
55.
A.
Raftery
,
J.
Cusumano
, and
D.
Sternad
, “
Chaotic frequency scaling in a coupled oscillator model for free rhythmic actions
,”
Neural Comput.
20
(
1
),
205
226
(
2008
).
56.
J. B.
Dingwell
,
J.
John
, and
J. P.
Cusumano
, “
Do humans optimally exploit redundancy to control step variability in walking?
,”
PLoS Comput. Biol.
6
(
7
),
e1000856
(
2010
).
57.
J. P.
Scholz
and
G.
Schöner
, “
The uncontrolled manifold concept: Identifying control variables for a functional task
,”
Exp. Brain Res.
126
(
3
),
289
306
(
1999
).
58.
S. S.
Ambike
and
M.
Tillman
, “
Cue-induced changes in the stability of finger force-production tasks revealed by the uncontrolled-manifold analysis
,”
J. Neurophysiol.
119
,
21
32
(
2018
).
59.
D.
Sternad
,
S.-W.
Park
,
H.
Müller
, and
N.
Hogan
, “
Coordinate dependence of variability analysis
,”
PLoS Comput. Biol.
6
(
4
),
e1000751
(
2010
).
60.
T. M.
Dijkstra
,
H.
Katsumata
,
A.
De Rugy
, and
D.
Sternad
, “
The dialogue between data and model: Passive stability and relaxation behavior in a ball bouncing task
,”
Nonlinear Stud.
11
(
3
),
319
344
(
2004
).
61.
D.
Sternad
, “
Human control of interactions with objects–variability, stability and predictability
,” in
Geometric and Numerical Foundations of Movements
(Springer,
2017
), pp.
301
335
.
62.
M. T.
Mason
, “
Toward robotic manipulation
,”
Annu. Rev. Control Rob. Auton. Syst.
1
,
1
28
(
2018
).
63.
F.
Cordella
,
A. L.
Ciancio
,
R.
Sacchetti
,
A.
Davalli
,
A. G.
Cutti
,
E.
Guglielmelli
, and
L.
Zollo
, “
Literature review on needs of upper limb prosthesis users
,”
Front. Neurosci.
10
,
209
(
2016
).
64.
N.
Hogan
and
D.
Sternad
, “
Dynamic primitives in the control of locomotion
,”
Front. Comput. Neurosci.
7
,
71
(
2013
).
65.
Q.-C.
Pham
,
N.
Tabareau
, and
J.-J.
Slotine
, “
A contraction theory approach to stochastic incremental stability
,”
IEEE Trans. Automat. Contr.
54
(
4
),
816
820
(
2009
).
66.
E. M.
Aylward
,
P. A.
Parrilo
, and
J.-J. E.
Slotine
, “
Stability and robustness analysis of nonlinear systems via contraction metrics and sos programming
,”
Automatica
44
(
8
),
2163
2170
(
2008
).
67.
Z.
Aminzare
and
E. D.
Sontag
, “
Contraction methods for nonlinear systems: A brief introduction and some open problems
,” in
IEEE 53rd Annual Conference on Decision and Control (CDC)
(
2014
), pp.
3835
3847
.
68.
H. K.
Khalil
,
Nonlinear Systems
(
Prentice-Hall
,
New Jersey
, 1996), Vol. 2(5).
69.
P. A.
Parrilo
, “
Structured semidefinite programs and semialgebraic geometry methods in robustness and optimization
,” Ph.D. thesis (
California Institute of Technology
,
2000
).
70.
W.
Lohmiller
and
J.-J. E.
Slotine
, “
Contraction analysis of non-linear distributed systems
,”
Int. J. Control
78
(
9
),
678
688
(
2005
).
71.
E.
Aylward
,
P. A.
Parrilo
, and
J.-J.
Slotine
, “
Algorithmic search for contraction metrics via SOS programming
,” in
American Control Conference (ACC)
(
2006
).
72.
I. R.
Manchester
and
J.-J. E.
Slotine
, “
Control contraction metrics: Convex and intrinsic criteria for nonlinear feedback design
,”
IEEE. Trans. Automat. Control
62
(
6
),
3046
3053
(
2017
).
73.
S.
Singh
,
A.
Majumdar
,
J.-J.
Slotine
, and
M.
Pavone
, “
Robust online motion planning via contraction theory and convex optimization
,” in
IEEE International Conference on Robotics and Automation (ICRA)
(
2017
), pp.
5883
5890
.
74.
Z.
Aminzare
, Ph.D. thesis,
Rutgers, The State University of New Jersey
,
New Brunswick
,
2015
.
75.
M. O.
Abe
and
D.
Sternad
, “
Directionality in distribution and temporal structure of variability in skill acquisition
,”
Front. Hum. Neurosci.
7
,
225
(
2013
).