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.

## I. PREDICTABILITY AND STABILITY

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 colleagues^{15,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 Milton^{23} 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.

## II. EXISTING METHODS FOR STABILITY ANALYSIS OF HUMAN MOVEMENT

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.

### A. Lyapunov exponents

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

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 $ \u03f5 0 $ to the initial condition. The time evolution of the perturbation $\u03f5(t)$ is governed by the equation

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

where $i=1,\u2026,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}

### B. Poincaré maps and Floquet multipliers

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}

### C. Variability analysis as proxy for stability

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 period^{49,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 Marin^{50} 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.

### D. Model-based local linear stability analysis

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.

## III. CONTRACTION THEORY

Contraction theory^{31} 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.

### A. Overview

Consider a nonlinear dynamical system

where $f$ is an $n\xd71$ nonlinear vector function and $x$ is the $n\xd71$ 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, $\delta x$, between the trajectories at a fixed time, as illustrated in Fig. 1.

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

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

This means that the virtual displacement $ \delta x $ converges exponentially to zero if the Jacobian $ \u2202 f \u2202 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

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

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

where $K$ is the generalized Jacobian

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 $\Theta (x,t)$ such that $ \Theta T \Theta >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.

### B. Favorable properties of contraction analysis

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.

## IV. THE EXPERIMENTAL TASK: TRANSPORTING A CUP OF COFFEE

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.

### A. A simplified model of the task dynamics

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.

The dynamics of the simplified model are

where $x(t)$ denotes the horizontal position of the cart, $\varphi (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 \xa8 $. 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$.

### B. Virtual implementation of the task

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\xd72.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 \u02d9 , x \xa8 ,\varphi , \varphi \u02d9 , \varphi \xa8 $ 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.

## V. CONTRACTION ANALYSIS

### A. Contraction regions

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 \u02d9 , \varphi , \varphi \u02d9 ) T = ( x 1 , x 2 , x 3 ) T $, and the resulting state space equations are

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

where

and

where

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 $\Theta (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 used^{31}

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

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.

## VI. RESULTS

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.

### A. Descriptive analysis of human data

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.

. | Subject 1 . | Subject 2 . | Subject 3 . | Subject 4 . | |
---|---|---|---|---|---|

Baseline | Early | $1.92 s \xb10.45$ | $1.77 s \xb10.45$ | $3.05 s \xb10.42$ | $1.88 s \xb10.47$ |

Late | $1.58 s \xb10.08$ | $1.30 s \xb10.29$ | $2.13 s \xb10.36$ | $1.43 s \xb10.08$ | |

Assistive | Early | $2.09 s \xb10.45$ | $1.56 s \xb10.53$ | $2.18 s \xb10.45$ | $1.61 s \xb10.46$ |

Late | $1.68 s \xb10.20$ | $1.35 s \xb10.20$ | $1.98 s \xb10.36$ | $1.55 s \xb10.33$ | |

Resistive | Early | $2.35 s \xb10.12$ | $1.32 s \xb10.27$ | $3.09 s \xb10.80$ | $1.71 s \xb10.50$ |

Late | $2.50 s \xb10.23$ | $1.33 s \xb10.22$ | $2.82 s \xb10.51$ | $1.29 s \xb10.14$ |

. | Subject 1 . | Subject 2 . | Subject 3 . | Subject 4 . | |
---|---|---|---|---|---|

Baseline | Early | $1.92 s \xb10.45$ | $1.77 s \xb10.45$ | $3.05 s \xb10.42$ | $1.88 s \xb10.47$ |

Late | $1.58 s \xb10.08$ | $1.30 s \xb10.29$ | $2.13 s \xb10.36$ | $1.43 s \xb10.08$ | |

Assistive | Early | $2.09 s \xb10.45$ | $1.56 s \xb10.53$ | $2.18 s \xb10.45$ | $1.61 s \xb10.46$ |

Late | $1.68 s \xb10.20$ | $1.35 s \xb10.20$ | $1.98 s \xb10.36$ | $1.55 s \xb10.33$ | |

Resistive | Early | $2.35 s \xb10.12$ | $1.32 s \xb10.27$ | $3.09 s \xb10.80$ | $1.71 s \xb10.50$ |

Late | $2.50 s \xb10.23$ | $1.33 s \xb10.22$ | $2.82 s \xb10.51$ | $1.29 s \xb10.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.

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

### B. Evaluating human solutions against 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 \u2212 $ 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 \u2212 $ and $ P + $ need not be fixed in state space, which is spanned by $\varphi $, $ \varphi \u02d9 $, and $ x \u02d9 $. 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.

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

## VII. DISCUSSION

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 programming^{70} and sums-of-squares programming^{71} 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, Aminzare^{74} 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.

## VIII. CONCLUSIONS

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.

## Acknowledgments

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.

## References

*et al.*,“