## Abstract

Nonlinear dynamical systems have been used in many disciplines to model complex behaviors, including biological motor control, robotics, perception, economics, traffic prediction, and neuroscience. While often the unexpected emergent behavior of nonlinear systems is the focus of investigations, it is of equal importance to create goal-directed behavior (e.g., stable locomotion from a system of coupled oscillators under perceptual guidance). Modeling goal-directed behavior with nonlinear systems is, however, rather difficult due to the parameter sensitivity of these systems, their complex phase transitions in response to subtle parameter changes, and the difficulty of analyzing and predicting their long-term behavior; intuition and time-consuming parameter tuning play a major role. This letter presents and reviews dynamical movement primitives, a line of research for modeling attractor behaviors of autonomous nonlinear dynamical systems with the help of statistical learning techniques. The essence of our approach is to start with a simple dynamical system, such as a set of linear differential equations, and transform those into a weakly nonlinear system with prescribed attractor dynamics by means of a learnable autonomous forcing term. Both point attractors and limit cycle attractors of almost arbitrary complexity can be generated. We explain the design principle of our approach and evaluate its properties in several example applications in motor control and robotics.

## 1. Introduction

In the wake of the development of nonlinear systems theory (Guckenheimer & Holmes, 1983; Strogatz, 1994; Scott, 2005), it has become common practice in several branches of science to model natural phenomena with systems of coupled nonlinear differential equations. Such approaches are motivated by the insight that coupling effects of nonlinear systems exhibit rich abilities for forming complex coordinated patterns without the need to explicitly plan or supervise the details of such pattern formation. Among the many different forms of nonlinear systems (e.g., high-dimensional, weakly coupled, strongly coupled, chaotic, Hamiltonian, dissipative), this letter addresses low-dimensional nonlinear systems, for example, as typically used to model phenomena of motor coordination or cognitive science (Kelso, 1995; Thelen & Smith, 1994).^{1} In this domain, there are often two modeling objectives. First, a model of a baseline behavior is required, as in generating a basic pattern for bipedal locomotion or reach-and-grasp in arm movement. Such behaviors are goal oriented; the focus is less on emergent coordination phenomena and more on achieving a task objective. After this baseline model has been accomplished, the second objective is to use this model to account for more complex phenomena with the help of the coupling dynamics of nonlinear systems. For instance, a typical example is the modulation of locomotion due to resonance entrainment of the pattern generator with the dynamics of a physical body (Nakanishi et al., 2004; Hatsopoulos & Warren, 1996). Another example is the coupling between motor control and perception (Dijkstra, Schoner, Giese, & Gielen, 1994; Kelso, 1995; Swinnen et al., 2004). In order to allow investigations of such second objectives, a dynamical systems model has to be found first.

Finding an appropriate dynamical systems model for a given behavioral phenomenon is nontrivial due to the parameter sensitivity of nonlinear differential equations and their lack of analytical predictability. Thus, modeling is often left to the intuition and the trial-and-error patience of the researchers. Many impressive studies have been generated in this manner (Schoner & Kelso, 1988; Schöner, 1990; Taga, Yamaguchi, & Shimizu, 1991; Schaal & Sternad, 1998; Kelso, 1995), but the lack of a generic modeling tool is unsatisfactory.

In this letter, we propose a generic modeling approach to generate multidimensional systems of weakly nonlinear differential equations to capture an observed behavior in an attractor landscape. The essence of our methodology is to transform well-understood simple attractor systems with the help of a learnable forcing function term into a desired attractor system. Both point attractor and limit cycle attractors of almost arbitrary complexity can be achieved. Multiple degrees of freedom can be coordinated with arbitrary phase relationships. Stability of the model equations can be guaranteed. Our approach also provides a metric to compare different dynamical systems in a scale-invariant and temporally invariant way.

We evaluate our approach in the domain of motor control for robotics, where desired kinematic motor behaviors will be coded in attractor landscapes and then converted into control commands with inverse dynamics controllers. Importantly, perceptual variables can be coupled back into the dynamic equations, such that complex closed-loop motor behaviors are created out of one relatively simple set of equations. Inspired by the biological concept of motor primitives (Giszter, Mussa-Ivaldi, & Bizzi, 1993; Mussa-Ivaldi, 1999), we call our system *dynamical movement primitives*, as we see them as building blocks that can used and modulated in real time for generating complex movements.

The following sections first introduce our modeling approach (see section 1), then, examine its theoretical properties (see section 2), and finally explore our approach in the example domain of motor control in various scenarios (see section 3). Matlab code is provided as supplemental material to allow readers to explore properties of the system.^{2} Early versions of the dynamical system presented in this letter have been published elsewhere in short format (Ijspeert, Nakanishi, & Schaal, 2002b, 2003) or some review articles (Schaal, Mohajerian, & Ijspeert, 2007; Schaal, Ijspeert, & Billard, 2003). Here, we review previous work and present our system in more detail, introduce examples of spatial and temporal couplings, and discuss issues related to generalization and coordinate systems. In the end, this letter presents a comprehensive and mature account of our dynamic modeling approach with discussions of related work, which will allow readers to apply or improve research on this topic.

## 2. A Learnable Nonlinear Attractor Systems

Before developing our model equations, it will be useful to clarify the specific goals pursued with this model:

Both learnable point attractor and limit cycle attractors need to be represented. This is useful to encode both discrete (i.e., point to point) and rhythmic (periodic) trajectories.

^{3}The model should be an autonomous system, without explicit time dependence.

The model needs to be able to coordinate multidimensional dynamical systems in a stable way.

Learning the open parameters of the system should be as simple as possible, which essentially opts for a representation that is linear in the open parameters.

The system needs to be able to incorporate coupling terms, for example, as typically used in synchronization studies or phase resetting studies and as needed to implement closed-loop perception-action systems.

The system should allow real-time computation as well as arbitrary modulation of control parameters for online trajectory modulation.

Scale and temporal invariance would be desirable; for example, changing the amplitude or frequency of a periodic system should not affect a change in geometry of the attractor landscape.

### 2.1. Model Development.

^{4}which, throughout this letter, we write in first-order notation, where is a time constant and and are positive constants. If the forcing term

*f*=0, these equations represent a globally stable second-order linear system with (

*z*,

*y*)=(0,

*g*) as a unique point attractor. With appropriate values of and , the system can be made critically damped (with ) in order for

*y*to monotonically converge toward

*g*. Such a system implements a stable but trivial pattern generator with

*g*as single point attractor.

^{5}The choice of a second-order system in equation 2.1 was motivated by our interest in applying such dynamical systems to motor control problems, which are most commonly described by second-order differential equations and require position, velocity, and acceleration information for control. In this spirit, the variables would be interpreted as desired position, velocity, and acceleration for a control system, and a controller would convert these variables into motor commands, which account for nonlinearities in the dynamics (Sciavicco & Siciliano, 2000; Wolpert, 1997). Section 2.1.7 expands on the flexibilities of modeling in our approach.

Choosing the forcing function *f* to be phasic (i.e., active in a finite time window) will lead to a point attractive system, while choosing *f* to be periodic will generate an oscillator. Since the forcing term is chosen to be nonlinear in the state of the differential equations and since it transforms the simple dynamics of the unforced systems into a desired (weakly) nonlinear behavior, we call the dynamical system in equation 2.1 the *transformation system*.

#### 2.1.1. A Point Attractor with Adjustable Attractor Landscape.

*f*in equation 2.1 could hypothetically be chosen, for example, as where are fixed basis functions and

*w*are adjustable weights. Representing arbitrary nonlinear functions as such a normalized linear combination of basis functions has been a well-established methodology in machine learning (Bishop, 2006) and also has similarities with the idea of population coding in models of computational neuroscience (Dayan & Abbott, 2001). The explicit time dependence of this nonlinearity, however, creates a nonautonomous dynamical system or, in the current formulation, more precisely a linear time-variant dynamical system. However, such a system does not allow straightforward coupling with other dynamical systems and the coordination of multiple degree-of-freedom in one dynamical system (e.g., as in legged locomotion; cf. section 3.2).

_{i}*x*where is a constant. Starting from some arbitrarily chosen initial state

*x*

_{0}, such as

*x*

_{0}=1, the state

*x*converges monotonically to zero.

*x*can thus be conceived of as a phase variable, where

*x*=1 would indicate the start of the time evolution and

*x*close to zero means that the goal

*g*has essentially been achieved. For this reason, it is important that

*x*=0 is a stable fixed point of these equations. We call this equation the

*canonical system*because it models the generic behavior of our model equations, a point attractor in the given case and a limit cycle in the next section. Given that equation 2.2 is a linear differential equation, there exists a simple exponential function that relates time and the state

*x*of this equation. However, avoiding the explicit time dependency has the advantage that we have obtained an autonomous dynamical system now, which can be modified online with additional coupling terms, as discussed in section 3.2.

*N*exponential basis functions , where and

*c*are constants that determine, respectively, the width and centers of the basis functions and

_{i}*y*

_{0}is the initial state

*y*

_{0}=

*y*(

*t*=0).

Note that equation 2.3 is modulated by both *g*−*y*_{0} and *x*. The modulation by *x* means that the forcing term effectively vanishes when the goal *g* has been reached, an essential component in proving the stability of the attractor equations. The modulation of equation 2.3 by *g*−*y*_{0} will lead to useful scaling properties of our model under a change of the movement amplitude *g*−*y*_{0}, as discussed in section 2.1.4. At the moment, we assume that , that is, that the total displacement between the beginning and the end of a movement is never exactly zero. This assumption will be relaxed later but allows a simpler development of our model. Finally, equation 2.3 is a nonlinear function in *x*, which renders the complete set of differential equations of our dynamical system nonlinear (instead of being a linear time-variant system), although one could argue that this nonlinearity is benign as it vanishes at the equilibrium point.

The complete system is designed to have a unique equilibrium point at (*z*, *y*, *x*)=(0, *g*, 0). It therefore adequately serves as a basis for constructing discrete pattern generators, with *y* evolving toward the goal *g* from any initial condition. The parameters *w _{i}* can be adjusted using learning algorithms (see section 2.1.6) in order to produce complex trajectories before reaching

*g*. The canonical system

*x*(see equation 2.2) is designed such that

*x*serves as both an amplitude and a phase signal. The variable

*x*monotonically and asymptotically decays to zero. It is used to localize the basis functions (i.e., as a phase signal) but also provides an amplitude signal (or a gating term) that ensures that the nonlinearity introduced by the forcing term remains transient due to asymptotical convergence of

*x*to zero at the end of the discrete movement.

Figure 1 demonstrates an exemplary time evolution of the equations. Throughout this letter, the differential equations are integrated using Euler integration with a 0.001 s time step. To start the time evolution of the equations, the goal is set to *g*=1, and the canonical system state is initialized to *x*=1. As indicated by the reversal of movement direction in Figure 1 (top left), the internal states and the basis function representation allow generating rather complex attractor landscapes.

Figure 2 illustrates the attractor landscape that is created by a two-dimensional discrete dynamical system, which we discuss in more detail in section 2.1.5. The left column in Figure 2 shows the individual dynamical systems, which act in two orthogonal dimensions, *y*_{1} and *y*_{2}. The system starts at *y*_{1}=0, *y*_{2}=0, and the goal is *g*_{1}=1, *g*_{2}=1. As shown in the vector field plots of Figure 2, at every moment of time (represented by the phase variable *x*), there is an attractor landscape that guides the time evolution of the system until it finally ends at the goal state. These attractor properties play an important role in the development of our approach when coupling terms modulate the time evolution of the system.

#### 2.1.2. A Limit Cycle Attractor with Adjustable Attractor Landscape.

Limit cycle attractors can be modeled in a similar fashion to the point attractor system by introducing periodicity in either the canonical system or the basis functions (which corresponds to representing an oscillator in Cartesian or polar coordinates, respectively). Here we present the second option (see Ijspeert et al., 2003, for an example of the first option).

*r*) and a phase signal () to the forcing term

*f*in equation 2.1: where the exponential basis functions in equation 2.7 are now von Mises basis functions, essentially gaussian-like functions that are periodic. Note that in case of the periodic forcing term,

*g*in equation 2.1 is interpreted as an anchor point (or set point) for the oscillatory trajectory, which can be changed to accommodate any desired baseline of the oscillation. The amplitude and period of the oscillations can be modulated in real time by varying, respectively,

*r*and .

Figure 3 shows an exemplary time evolution of the rhythmic pattern generator when trained with a superposition of several sine signals of different frequencies. It should be noted how quickly the pattern generator converges to the desired trajectory after starting from zero initial conditions. The movement is started simply by setting the *r*=1 and . The phase variable can be initialized arbitrarily: we chose for our example. More informed initializations are possible if such information is available from the context of a task; for example, a drumming movement would normally start with a top-down beat, and the corresponding phase value could be chosen for initialization. The complexity of attractors is restricted only by the abilities of the function approximator used to generate the forcing term, which essentially allows almost arbitrarily complex (smooth) attractors with modern function approximators.

#### 2.1.3. Stability Properties.

*g*and the forcing term

*f*in one expression results in where

*u*is a time-variant input to the linear spring-damper system. Equation 2.8 acts as a low-pass filter on

*u*. For such linear systems, with appropriate and , for example, from critical damping as employed in our work, it is easy to prove bounded-input, bounded-output (BIBO) stability (Friedland, 1986), as the magnitude of the forcing function

*f*is bounded by virtue that all terms of the function (i.e., basis functions, weights, and other multipliers) are bounded by design. Thus, both the discrete and rhythmic system are BIBO stable. For the discrete system, given that

*f*decays to zero,

*u*converges to the steady state

*g*after a transient time, such that the system will asymptotically converge to

*g*. After the transient time, the system will exponentially converge to

*g*as only the linear spring-damper dynamics remain relevant (Slotine & Li, 1991). Thus, ensuring that our dynamical systems remain stable is a rather simple exercise of basic stability theory.

Another path to prove stability for our approach was suggested by Perk and Slotine (2006), who proved that our dynamical systems equations are two hierarchically coupled systems, each fulfilling the criterion of contraction stability (Lohmiller & Slotine, 1998). Contraction theory provides that any parallel or serial arrangement of contraction stable systems will be contraction stable too, which concludes the stability proof. This property will be useful below, where we create multiple degree-of-freedom dynamical systems, which inherit their stability proof from this contraction theory argument.

#### 2.1.4. Invariance Properties.

*w*of the forcing term, and (3) the global timescaling parameter and the goal parameter

_{i}*g*, or amplitude parameter

*r*. We assume that the first and second are kept constant for a particular behavior, but the parameters of the third can be varied as we consider them natural high-level parameters that should adjust the behavior to a particular context without changing the behavior qualitatively. Thus, formally, we wish that the attractor landscape of the dynamical systems does not change qualitatively after a change of ,

*g*, or

*r*. This topic can be addressed in the framework of topological equivalence (Jackson, 1989). Mathematically, if two dynamical systems and are topologically equivalent, then there exists an orientation preserving homeomorphism and , where the notation denotes a functional mapping. When scaling the movement amplitude

*g*−

*y*

_{0}in the discrete dynamical system (see equations 2.1–2.4) with a scalar mapping

*k*, , the following scaling law is an orientation-preserving homeomorphism between the original equations using

*g*−

*y*

_{0}and the scaled differential equations using

*k*(

*g*−

*y*

_{0}): This result can be verified by simply inserting equations 2.9 into the scaled equations using the simple trick of

*k*(

*g*−

*y*

_{0})+

*ky*

_{0}as a goal. Similarly, when scaling in the rhythmic system (see equations 2.1 and 2.5–2.7), the same scaling law (see equation 2.9) allows proving topological equivalence (for simplicity, assume that

*g*=0 in equation 2.1 for the rhythmic system, which can always be achieved with a coordinate shift). For the scaling of the time constant , topological equivalence for both the discrete and rhythmic systems can be established trivially with Figure 4 illustrates the spatial (see Figure 4a) and temporal (see Figure 4b) invariance using the example from Figure 1. One property that should be noted is the mirror-symmetric trajectory in Figure 4a when the goal is at a negative distance relative to the start state. We discuss the issue again in section 3.4.

Figure 5 provides an example of why and when invariance properties are useful. The blue (thin) line in all subfigures shows the same handwritten cursive letter a that was recorded with a digitizing tablet and learned by a two-dimensional discrete dynamical system. The letter starts at a *StartPoint*, as indicated in Figure 5a, and ends originally at the goal point *Target*_{0}. Superimposed on all subfigures in red (thick line) is the letter a generated by the same movement primitive when the goal is shifted to *Target*_{1}. For Figures 5a and 5b, the goal is shifted by just a small amount, while for Figures 5c and 5d, it is shifted significantly more. Importantly, for Figures 5b and 5d, the scaling term *g*−*y*_{0} in equation 2.3 was left out, which destroys the invariance properties as described above. For the small shift of the goal in Figures 5a and 5b, the omission of the scaling term is qualitatively not very significant: the red letter “a” in both subfigures looks like a reasonable “a.” For the large goal change in Figures 5c and 5d, however, the omission of the scaling term creates a different appearance of the letter “a,” which looks almost like a letter “*u*.” In contrast, the proper scaling in Figure 5c creates just a large letter “a,” which is otherwise identical in shape to the original one. Thus, invariance properties are particularly useful when generalization properties are required that are not just confined to a very local area of the originally learned primitive. It should also be noted that a simple change of the goal state automatically creates a complex rescaling of the entire movement, without any need for explicit rescaling computations.

In conclusion, our dynamical systems equations for point attractors and limit cycles are actually a model of a family of similar behaviors, parameterized by the timing, goal, or amplitude. Importantly, the weights *w _{i}* in the forcing term remain unscaled, that is, constant under these scalings, which will allow us to use them later as means to statistically classify a behavior.

#### 2.1.5. Extension to Multiple Degrees of Freedom.

In order to obtain a system with multiple DOFs, three approaches can be pursued. First, every DOF could simply have its own complete set of dynamic equations, that is, we copy the discrete or rhythmic systems above for every DOF. While this approach is theoretically viable, it has the disadvantage that there is no coordination between the DOFs; for a longer run of the dynamical systems, for example, they may numerically diverge from each other. Moreover, disturbance rejection will normally require that in case of a disturbance of one DOF, the coordination among DOFs is not destroyed.

A second possibility would be to create coupling terms between the different DOFs, in particular, between the canonical systems, for example, as done in Taga et al. (1991). This is especially interesting for maintaining desired phase lags between DOFs during periodic movements, such as for locomotion, and for switching between different rhythmic patterns (e.g., gaits) by changing the coupling terms. However, such modeling becomes rather complex in terms of tuning coupling parameters for tight synchronization, analyzing stability, and dealing with the complex transient behavior before the system phase-locks.

A simpler multi-DOF approach is to share one canonical system among all DOFs and maintain only a separate set of transformation systems (see equation 2.1) as well as separate forcing terms for each DOF (see Figure 6). Thus, the canonical system provides the temporal coupling between DOFs, while the transformation system achieves the desired attractor dynamics for each DOF. As mentioned above, the argument of contraction stability put forward by Perk and Slotine (2006) extends to this multi-DOF approach, which is thus guaranteed to provide stable coordination among the DOFs. We have been pursuing this solution in our evaluations below. Interestingly, the canonical system becomes a central clock in this methodology, not unlike the assumed role of central pattern generators in biology (Getting, 1988; Grillner, 1981), in particular with the notion of two-level central pattern generators (McCrea & Rybak 2008), with one part of the network providing synchronized timing signals (in our case, the canonical systems) and another part of the network providing the shape of the motor patterns (the transformation systems).

Note that for different modeling purposes, different coupling approaches can be chosen. For some applications, for instance, on a humanoid robot, it might be useful to have one canonical system per limb with multiple transformation systems for controlling the joints within a limb. By adding coupling terms between canonical systems, movements of different limbs can then be synchronized with stable and adjustable phase differences.

#### 2.1.6. Learning the Attractor Dynamics from Observed Behavior.

Our systems are constructed to be linear in the parameters *w _{i}*, which allows applying a variety of learning algorithms to fit the

*w*. In this letter, we focus on a supervised learning framework. Of course, many optimization algorithms could be used too if only information from a cost function is available.

_{i}We assume that a desired behavior is given by one or multiple desired trajectories in terms of position, velocity, and acceleration triples , where .^{6} Learning is performed in two phases: determining the high-level parameters (*g*, *y*_{0}, and for the discrete system or *g*, *r*, and for the rhythmic system) and then learning the parameters *w _{i}*.

For the discrete system, the parameter *g* is simply the position at the end of the movement, *g*=*y _{demo}*(

*t*=

*P*) and, analogously,

*y*

_{0}=

*y*(

_{demo}*t*=0). The parameter must be adjusted to the duration of the demonstration. In practice, extracting from a recorded trajectory may require some thresholding in order to detect the movement onset and end. For instance, a velocity threshold of 2% of the maximum velocity in the movement may be employed, and could be chosen as 1.05 times the duration of this thresholded trajectory piece. The factor 1.05 is to compensate for the missing tails in the beginning and the end of the trajectory due to thresholding.

For the rhythmic system, *g* is an anchor point that we set to the midposition of the demonstrated rhythmic trajectory . The parameter is set to the period of the demonstrated rhythmic movement divided by . The period must therefore be extracted beforehand using any standard signal processing (e.g., a Fourier analysis) or learning methods (Righetti, Buchli, & Ijspeert, 2006; Gams, Ijspeert, Schaal, & Lenarcic, 2009). The parameter *r*, which will allow us to modulate the amplitude of the oscillations (see the next section), is set, without loss of generality, to an arbitrary value of 1.0.

The learning of the parameters *w _{i}* is accomplished with locally weighted regression (LWR) (Schaal & Atkeson, 1998). It should be emphasized that any other function approximator could be used as well (e.g., a mixture model, a gaussian process). LWR was chosen due to its very fast one-shot learning procedure and the fact that individual kernels learn independent of each other, which will be a key component to achieve a stable parameterization that can be used for movement recognition in the evaluations (see section 3.3).

*f*are to be adjusted such that

*f*is as close as possible to

*f*.

_{target}These equations are easy to implement. It should be noted that if multiple demonstrations of a trajectory exist, even at different scales and timing, they can be averaged together in the above locally weighted regression after the *f _{target}* information for every trajectory at every time step has been obtained. This averaging is possible due to the invariance properties mentioned above. Naturally, locally weighted regression also provides confidence intervals on all the regression variables (Schaal & Atkeson, 1994, 1998), which can be used to statistically assess the quality of the regression.

The regression performed with equation 2.14 corresponds to what we will call a batch regression. Alternatively, we can also perform an incremental regression. Indeed, the minimization of the cost function *J _{i}* can be performed incrementally, while the target data points

*f*(

_{target}*t*) come in. For this, we use recursive least squares with a forgetting factor of (Schaal & Atkeson, 1998) to determine the parameters

*w*.

_{i}#### 2.1.7. Design Principle.

In developing our model equations, we made specific choices for the canonical systems, the nonlinear function approximator, and the transformation system, which is driven by the forcing term. But it is important to point out that it is the design principle of our approach that seems to be the most important, and less the particular equations that we chose for our realization. As sketched in Figure 7, there are three main ingredients in our approach. The canonical system is a simple (or well-understood) dynamical system that generates a behavioral phase variable, which is our replacement for explicit timing. Either point attractive or periodic canonical system can be used, depending on whether discrete or rhythmic behavior is to be generated or modeled. For instance, while we chose a simple first-order linear system as a canonical system for discrete movement, it is easily possible to choose a second-order system instead (see Ijspeert et al., 2003) or even nonlinear equations. Similarly, the rhythmic canonical system in our case was a phase oscillator, but a van der Pol oscillator, Duffing oscillator (Strogatz, 1994), or any other oscillatory system could be a useful alternative choice. For instance, if the transient behavior of a rhythmic system from a rest state to an oscillatory behavior would be of interest, using a canonical oscillator that has a Hopf bifurcation could be beneficial.

With the input from the canonical system, the nonlinear function approximator generates a forcing term. As mentioned before, any choice of function approximator is possible (Bishop, 2006).

Finally, the transformation system should be a dynamical system that is easily analyzed and manipulated when excited with a forcing term. We used a critically damped linear spring system, but other systems, like higher-order or lower-order dynamical systems are equally possible. This choice is partially guided by which level of derivatives the output behavior of the entire dynamical system should have.

All systems together should fulfill the principle of structural equivalence, should be autonomous, and should have easily analyzable stability properties. Obviously a large variety of model equations can be generated, tailored for different contexts. We address the coupling terms in Figure 7 in section 3.2.

#### 2.1.8. Variations.

Given the general design principle from the previous section, we briefly discuss some variations of our approach that have been useful in some applications.

*g*to a new value, equation 2.1 generates a discontinuous jump in the acceleration . This can be avoided by filtering the goal change with a simple first-order differential equation: In this formulation,

*g*

_{0}is the discontinuous goal change, while

*g*is now a continuous variable. This modification does not affect the scaling properties and stability properties of the system and is easily incorporated in the learning algorithm with LWR.

*f*cannot create a discontinuity in acceleration. For this purpose, it is useful to work with a second-order canonical system (instead of equation 2.2), whose constants are chosen to be critically damped (see Ijspeert et al., 2003, for an example). The forcing term is reformulated as Note that in contrast to equation 2.3, the multiplication on the right-hand side has

*v*instead of

*x*.

*v*starts at zero at moment onset and returns to zero at the end of the movement, such that the forcing term vanishes at both movement end and beginning. The danger of a discontinuous initial acceleration is thus avoided.

*y*

_{0}and goal

*g*coincide. As can be seen in equation 2.3, the forcing term would always be zero in this case. This problem can be circumvented by defining the forcing term as

*g*and

_{fit}*y*

_{0,fit}denote the goal and start point used to fit the weights of the nonlinear function. Thus, during fitting, this quotient is always equal to +1 or −1, such that the forcing term is guaranteed to be active. The numerically very small positive constant avoids dividing by zero. There is, however, a numerical danger that after learning, a new goal different from

*y*

_{0}could create a huge magnification of the originally learned trajectory. For instance, if during learning we had

*g*−

_{fit}*y*

_{0,fit}=0, and , a change of the goal offset to the start position of

*g*−

*y*

_{0}=0.01 would create a multiplier of 100 in equation 2.18. While theoretically fine, such a magnification may be practically inappropriate. A way out is to modify the forcing term to use the maximal amplitude of a trajectory as scaling term,

*A*=max(

*y*)−min(

*y*): In this variant, the goal

*g*becomes decoupled from the amplitude

*A*, and both variables can be set independently. While it is easy to lose the strict property of structural equivalence in this way, it may be practically more appropriate for certain applications. Hoffmann, Pastor, Park, and Schaal (2009) suggested a similar approach.

*r*

_{0}creates a smooth change in

*r*. Another useful addition could be to make the input vector to the forcing term in equation 2.6 two-dimensional: This input vector is composed of first-order Fourier terms, which exploits our knowledge that the rhythmic system must be periodic.

^{7}Simple oscillations are easier to fit with this variant, while for more complex oscillations, numerically we could not see an advantage in our experiments. For this two-dimensional formulation, the weighted regression in equation 2.14 needs to become a two-dimensional weighted regression, which requires a matrix inversion instead of the simple division. This is an insignificant change, of course.

#### 2.1.9. Summary.

Table 1 summarizes the key equations of our model. Given the previous sections, the equations should be self-explanatory. The only special component is the initialization of *h _{i}* for the discrete system, which generates equal spacing in time

*t*by using the analytical solution of equation 2.2 to compute the corresponding spacing in

*x*. In the following section, we provide a variety of evaluations of our suggested approach to learning attractor systems models.

Discrete . | . | Rhythmic . | . |
---|---|---|---|

Transformation system: | Transformation system: | ||

Equation 2.1 | Equation 2.1 | ||

Canonical system: | Canonical system: | ||

Equation 2.2 | Equation 2.5 | ||

Forcing term: | Forcing term: | ||

Equation 2.3 | Equation 2.6 | ||

Optional terms: | Optional terms: | ||

Equation 2.15 | Equation 2.15 | ||

Equation 2.20 | |||

h= equal spacing in _{i} | h= equal spacing in _{i} | ||

Discrete . | . | Rhythmic . | . |
---|---|---|---|

Transformation system: | Transformation system: | ||

Equation 2.1 | Equation 2.1 | ||

Canonical system: | Canonical system: | ||

Equation 2.2 | Equation 2.5 | ||

Forcing term: | Forcing term: | ||

Equation 2.3 | Equation 2.6 | ||

Optional terms: | Optional terms: | ||

Equation 2.15 | Equation 2.15 | ||

Equation 2.20 | |||

h= equal spacing in _{i} | h= equal spacing in _{i} | ||

Notes: The high-level design parameters of the discrete system are , the temporal scaling factor, and *g*, the goal position. The design parameters of the rhythmic system are *g*, the baseline of the oscillation; , the period divided by ; and *r*, the amplitude of oscillations. The terms *C _{t}* and

*C*are coupling terms that are application dependent and explained in section 3.2. The parameters

_{c}*w*are fitted to a demonstrated trajectory using locally weighted learning. The parameters , and

_{i}*c*are positive constants. Unless stated otherwise, the values at the bottom of the table are used in this letter.

_{i}## 3. Evaluations

In the next sections, we present and review several experimental evaluations of applying our approach to learning attractor systems in the domain of motor control, using both simulation and robotic studies. The evaluations are intended to demonstrate the properties of our methodology, but also the domain-specific choices that need to be made. We will discuss imitation learning of discrete and rhythmic movement, online modulation with the help of coupling terms, synchronization and entrainment phenomena, and movement recognition based on a motor generation framework.

### 3.1. Imitation Learning.

Imitation learning (Schaal, 1999) is simply an application of supervised learning given in section 2.1.6. Ijspeert et al. (2003) and Ijspeert, Nakanishi, and Schaal (2002a) demonsrated imitation learning with a 30 degrees-of-freedom (DOFs) humanoid robot (Atkeson et al., 2000) for performing a tennis forehand, a tennis backhand, and rhythmic drumming movements. Importantly, these tasks required the coordination and phase locking of 30 DOFs, which was easily and naturally accomplished in our approach. The imitated movement was represented in joint angles of the robot or Cartesian coordinates of an end effector. Indeed, only kinematic variables are observable in imitation learning. Thus, our dynamical systems represented kinematic movement plans—accelerations as a function of position and velocity. The robot was equipped with a controller (a feedback PD controller and a feedforward inverse dynamics controller) that could accurately track the kinematic plans (i.e., generate the torques necessary to follow a particular joint angle trajectory, given in terms of desired positions, velocities, and accelerations). Movies illustrating these results can be viewed at http://biorob.epfl.ch/dmps.

There are numerous other ways to accomplish imitation learning. In robotics, one of the most common classic methods to represent movement plans is by means of third-order or fifth-order splines, which could be equally applied to the learning from demonstration approach in this section. For instance, Wada and Kawato (2004) presented an elegant algorithm that recursively fits a demonstrated trajectory with a growing number of spline nodes until an accuracy criterion is reached. However, splines are nonautonomous representations with no attractor properties, not dynamical systems. Thus, while splines are effective for imitation learning, they would not allow online modulation properties as presented in section 3.2. Rescaling the splines in space and time for generalization is possible but requires an explicit recomputing of the spline nodes.

Another alternative to fitting a dynamical system to observed data was presented by Khansari-Zadeh and Billard (2010), who used a mixture model approach to estimate the entire attractor landscape of a movement skill from several sample trajectories. To ensure stability of the dynamical system toward an attractor point, a constraint optimization problem has to be solved in a nonconvex optimization landscape. This approach is different from ours in that it creates the attractor landscape in the state-space of the observed data, while our approach creates the attractor landscape in the phase space of the canonical system. The latter is low dimensional even for a high-DOF system. A state-space mixture model for our humanoid robot above would require a 60-dimension state space and thus would create computational and numerical problems. However, state-space models can represent much more complex attractor landscapes, with different realizations of a movement in different parts of the state-space. Our approach creates inherently stereotypical movements to the goal, no matter where the movements starts. Thus, the state-space approach to fitting a dynamical systems appears to pursue a quite different set of goals than our trajectory-based approach does.

### 3.2. Online Modulation of the Dynamical Systems.

One goal of modeling with dynamical systems is to use the ability of coupling phenomena to account for complex behavior. Imitation learning from the previous section demonstrated how to initialize dynamical systems models with learning from demonstration. In this section, we discuss different methods for how dynamical system models can be modulated online to take dynamic events from the environment into account. Those online modulations are among the most important properties offered by a dynamical systems approach, and these properties cannot easily be replicated without the attractor properties of our proposed framework.

#### 3.2.1. Spatial Coupling.

*C*to equation 2.1 to become

_{t}It should be noted that a prudent choice of the coupling term is critical and often needs to be specialized for different objectives. The design of coupling terms is a research topic by itself. A typical example from the domain of motor control is obstacle avoidance with the help of potential fields (Khatib, 1986; Koditschek, 1987; Rimon & Koditschek, 1992). Here, obstacles are modeled as repelling potential fields that are designed to automatically push a control system to circumnavigate them in an online reactive way instead of deliberative preplanning. Such reactive behavior assumes that obstacles may appear in an unforeseen and sudden way, such that preplanning is not possible or useful.

**y**=[

*y*

_{1}

*y*

_{2}

*y*

_{3}]

^{T}, with as the corresponding velocity vector. The objective of a movement is to generate a reaching movement from any start state to a goal state

**g**=[

*g*

_{1}

*g*

_{2}

*g*

_{3}]

^{T}. The discrete dynamical system is initialized with a minimum jerk movement (Flash & Hogan, 1985), which is frequently used as an approximate model of smooth human movement. On the way to the goal state, an obstacle is positioned at

**o**=[

*o*

_{1}

*o*

_{2}

*o*

_{3}]

^{T}and needs to be avoided. A suitable coupling term

**C**

_{t}=[

*C*

_{t,1}

*C*

_{t,2}

*C*

_{t,3}]

^{T}for obstacle avoidance can be formulated as The angle is interpreted as the angle between the velocity vector and the difference vector (

**o**−

**y**) between the current position and the obstacle. The vector

**r**is the vector that is perpendicular to the plane spanned by and (

**o**−

**y**), and serves to define a rotation matrix

**r**, which causes a rotation of 90 degrees about

**r**(Sciavicco & Siciliano, 2000). Intuitively, the coupling term adds a movement perpendicular to the current movement direction as a function of the distance vector to the obstacle (see Hoffmann et al., 2009, for more details). The constants are chosen to and .

Figure 8 illustrates the behavior that the obstacle-avoidance coupling term generates for various trajectories starting from different initial position around the origin **y**=[0 0 0]^{T} but ending at the same goal state **g**=[1 1 1]^{T}. Depending on the start position, the coupling term creates more or less curved movements around the obstacle at **o**=[0.5 0.5 0.5]^{T}. The behavior looks intuitively natural, which is not surprising as it was derived from human obstacle-avoidance behavior (Fajen & Warren, 2003).

A more complex example of spatial coupling is given in Figure 9. Using imitation learning, a placing behavior of a cup on a target was coded in a discrete dynamical system for a 3D end effector movement of the robot, a Sarcos Slave 7 DOF robot arm. The first row of images shows the unperturbed behaviors. In the second row, the (green) target is suddenly moved to the right while the robot has already begun moving. This modification corresponds to a change of the goal parameter *g*. The third row of images demonstrates an avoidance behavior based on equation 3.2, when the blue ball comes too close to the robot's movement. We emphasize that one single discrete dynamical system created all these different behaviors; there was no need for complex abortion of the ongoing movement or replanning. More details can be found in Pastor, Hoffmann, Asfour, and Schaal (2009).

#### 3.2.2. Temporal Coupling.

Figure 10 illustrates this approach for a robot drumming example. A rhythmic 7 DOF dynamical system was initialized with a drumming pattern that consisted of one slow drumbeat followed by two quick ones (Pongas et al., 2005). An external metronome generated the beat of a rhythm to which the slow drumbeat of the robot was to synchronize. In the beginning, the robots started from immobility (). Within two beats (the time needed to extract the frequency from the acoustic signal), perfect synchronization and phase locking is achieved with a 0.15 Hz signal—very rapid synchronization. Afterward, we had the external metronome pace increase frequency slowly to 0.5 Hz to demonstrate the continuous adaptation ability of the oscillator. Figure 10 shows the elbow DOF angular acceleration of the drumming pattern, which has the most significant contribution to the whole arm movement. As can be seen, with increasing frequency, the overall acceleration amplitude of the pattern changes but not the qualitative waveform. This property is due to the invariance properties of the dynamical systems. All other DOFs of the arm demonstrate the same behavior and are equally phase-locked to the beat of the metronome.

#### 3.2.3. Temporal and Spatial Coupling.

*y*and velocity of the 1 DOF discrete dynamical system drive the time evolution of

*y*, which can be interpreted as the position of a simple point mass-controlled by a proportional-derivative controller with gains

_{a}*k*and

_{p}*k*. We use the dynamics from Figure 1 to generate the input . At time

_{v}*t*=0.35 s, the point mass is suddenly blocked from any further motion and released again at

*t*=0.9 s. Thus, for we have . Without coupling terms, the dynamical system would just continue its time evolution, regardless of what happens to the point mass. For a practical control system, this behavior is undesirable as the desired trajectory

*y*can move away significantly from the current state

*y*, such that on release of the mass, it would create a dangerously large motor command and jump to catch up with the desired target.

_{a}*e*=

*y*−

_{a}*y*. This error is used as an additive coupling term in the transformation system, which hinders the state

*y*to evolve too far away from

*y*. Equation 3.11 affects the time constant of all differential equations of the dynamical system, that is, both the canonical and the transformation systems. This modification of the time constant slows the temporal evolution of the dynamics in case of a significant tracking error. The constants are chosen to be . It should be noted that many other coupling terms could be created to achieve similar behavior and that our realization is just a simple and intuitive design of such a coupling term.

_{a}Figure 11 illustrates the behavior due to these coupling terms in comparison to the unperturbed (dashed line) time evolution of the dynamics. The top-left plot of Figure 11 also shows with the dash-dot line the position *y _{a}*. During the holding time period, the entire dynamics comes almost to a stop and then smoothly resumes after the release of the mass roughly with the same behavior as where the system had left off; the system continues with the negative dip before moving toward the goal. Without the coupling terms,

*y*would already have evolved all the way to the goal position, and the error between

*y*and

_{a}*y*would have grown very large. These types of couplings have been used successfully with the humanoid robot (see the video at http://biorob.epfl.ch/dmps). It should also be emphasized that many different behaviors could be generated with other coupling terms, and it is really up to the modeler to decide which behavioral properties to realize with a particular realization of coupling terms.

### 3.3. Movement Recognition.

Due to the temporal and spatial invariance of our representation, an interesting aspect of our dynamical systems approach arises as trajectories that are topologically similar tend to be fit by similar parameters *w _{i}*. This property opens the possibility of using our representation for movement recognition, for example, in order to recognize the intent, the style, or the sequential decomposition of movement.

^{8}It should be noted that such recognition is about spatiotemporal patterns, not just spatial patterns, as in typical character recognition research. For instance, copying the signature of another person is usually done with the aim of generating the same spatial pattern on paper. The temporal pattern, however, is usually vastly different during the copying process. This temporal difference would manifest itself in a change of

*w*in our representation, even if the spatial copy were perfect. Thus, the following example should not be confused with standard spatial character recognition.

_{i}To illustrate this idea, we carried out a simple task of fitting trajectories performed by a human user when drawing two-dimensional single-stroke patterns. The 26 letters of the Graffiti alphabet used in hand-held computers were chosen. These characters are drawn in a single stroke and are fed as a two-dimensional trajectory (*y*_{1}(*t*), *y*_{2}(*t*)) to be fitted by a 2 DOF discrete dynamical system. Five examples of each character were presented (see Figure 12 for some examples).

Similarities between two characters *a* and *b* can be measured by computing the correlation between their parameter vectors. The correlation corresponds to , where **w**_{a} and **w**_{b} are the parameter vectors of characters *a* and *b*. These vectors are the union—i.e., —of the parameter vectors for the *y*_{1}(*t*) and *y*_{2}(*t*) trajectories, respectively, and . It should be pointed out that using the LWR function approximator (see section 2.1.6) is particularly advantageous in this application. LWR learns each coefficient of **w** locally using nearest-neighbor interpolation, unlike mixture models or gaussian process regression (Schaal & Atkeson, 1998), which are global function approximators. Thus, for instance, a variation at the end of a movement will not affect the parameter values at the beginning of a movement, which creates a very robust parameter representation, suitable for the correlation-based classification above.

Figure 13 shows the correlations of all 130 instantiations of the characters. The correlation indicates similarity in the acceleration profiles. As illustrated by high correlation values in the diagonal of the correlation matrix, the correlations between instances of the same character tend to be systematically higher than correlations between instances of different characters. These similarities in weight space can therefore serve as a basis for recognizing demonstrated movements by fitting them and comparing the fitted parameters *w _{i}* with those of previously learned characters in memory. In this example, a simple one-nearest-neighbor classifier in weight space would serve the purpose. More sophisticated classifiers like support vector machines would be possible too but are beyond the scope of this letter. Using such a classifier within the whole alphabet (5 instances of each letter), we obtained a 87% recognition rate (113 out of the 130 instances had a highest correlation with an instance of the same letter). The wrongly labeled instances are sometimes similar-looking letters (e.g., a Q being recognized as a G), but sometimes significantly different (e.g., a J being recognized as an S). The latter cases usually mean that the correlation coefficient is rather low, such that the nearest neighbor is chosen by chance rather than due to a significant correlation coefficient. Such cases could be excluded by more sophisticated classifiers that would employ, for example, confidence levels in decision making.

As a baseline comparison, we repeated this spatiotemporal character recognition experiment with dynamic time warping (DTW) (Sakoe & Chiba, 1987), a dynamic programming technique to compute a matching score between two arbitrary trajectories. DTW accomplished a 79% recognition rate, which is significantly lower than our dynamical systems representation. DTW is particularly well suited to match similar trajectories that have temporal variance. For changes in spatial scale, DTW is less suitable, while the invariance properties of our dynamical systems representation are built to be immune to such changes. To demonstrate this fact, we changed the scale of the five recorded instances per character by multiplying the position values of the trajectories for the first character with 1.0, the second with 1.25, the third with 1.5, the fourth with 1.75, and the fifth with 2.0. We then repeated the character recognition experiment. As expected, the results with the dynamical systems representation remained virtually identical as without scaling. For DTW, however, the recognition rate reduced to only 25%.

Further studies are required to evaluate the quality of recognition in larger training and test sets. What we wanted to demonstrate is the ability for recognition without any specific system tuning or sophisticated classification algorithm.

### 3.4. Generalization and Coordinate Systems.

Section 2.1.4 already alluded to the invariance properties of our dynamical systems model. The invariance properties are derived from the mathematical principle of structural equivalence and thus are theoretically sound. However, these invariance properties manifest themselves quite differently as a function of different coordinate systems.

Figure 14 illustrates some of the interesting issues that can arise in generalizing a learned movement to new goal states. We used a 2 DOF discrete dynamical system to represent a movement from the origin to a goal point, where on the way to the goal, the movement makes a loop. In Figures 14a, 14b, and 14c, the original movement, which was used to fit the dynamical system, is shown in a heavy (red) line. These original movements are in the first quadrant of the coordinate systems and differ only in the orientation of the goal relative to the start; Figure 14a is at roughly a 45 degree angle, Figure 14b is at about 10 degrees, and Figure 14c is at a zero angle.

In Figures 14a and 14b, we used the formulation from Table 1 to fit the model. Then we applied the model to generate movement to six different targets, distributed with 60 degrees difference on a circle around the origin.

The generalization to new targets in Figure 14a looks, from our human intuition, reasonable. The loop gets slightly distorted for goals close to the horizontal or vertical line. Indeed, if the goal were to lie exactly on a horizontal or vertical line, the loop would collapse. Movements in the second quadrant are mirror symmetric to similar movements in the first quadrant. The same holds for the third and fourth quadrants. Between the first and third quadrants, and the second and fourth quadrants, we observe point symmetry.

In Figure 14b, we repeated the experiment of Figure 14a, just that the original movement is closer to the horizontal line. The generalization pattern looks quite different. This effect comes from the fact that the start and goal state of *y*_{2} are rather close together, such that in equation 2.3, the term (*g*_{2}−*y*_{0,2}) causes a rather large amplification of the movement in *y*_{2} when the goal moves farther away from the origin in *y*_{2} direction. Theoretically, from a spatial invariance point of view, there is nothing wrong with this generalization pattern, just that it does not look very intuitive to us anymore and that it would risk hitting joint angle limits.

In Figure 14c, we changed the coordinate system for representing the movement to a local coordinate system. Here, one axis of the movement is always aligned with the line from the start to the goal point. The second coordinate axis is perpendicular. Thus, as can be seen in the two plots on the right of Figure 14c, *y*_{2} has a zero difference between start and goal coordinate, and only *y*_{1} has a difference. Thus, we use equation 2.19 for representing *y*_{2}. In order to generalize to new goals, the movement is generated in local coordinates and requires a subsequent coordinate transformation to rotate the trajectory to global coordinates. With this strategy, the generalization pattern in Figure 14c looks the most appealing to the human eye. Mathematically, all we did was change coordinates to represent our dynamical systems model.

In conclusion, the invariance properties in our dynamical systems model are mathematically well founded, but the choice of coordinates for representing a model can make a big difference in how generalization of the dynamics systems appears. From a practical point of view, one should first carefully investigate what properties a model requires in terms of temporal and spatial invariance and then realize these properties by choosing the most appropriate variant of the dynamical systems model and the most appropriate coordinate system for modeling.^{9}

## 4. Related Work

As mentioned in section 1, the central goal of our work was to derive learnable attractor dynamical systems models for goal-directed behavior and to explore generalization and coupling phenomena with such systems as models for biological and technical systems. Thus, we roughly classify related work according to more biological or more technical domains. Besides general nonlinear systems theory, most related work comes from the field of biological and artificial motor control.

### 4.1. Neural Control of Movement.

In the field of neural control of movement, early work on dynamical systems models was suggested by Bullock and Grossberg (1989), where the features of human point-to-point reaching movements were modeled with a point attractor system. Giszter et al. (1993) and Mussa-Ivaldi (1997, 1999) developed the idea of superposition of force fields to model wiping behaviors in the frog and also the generation of human arm movement. Kelso, Scholtz, and Schoner (1988), Schöner and Kelso (1988), Schöner (1990), Haken, Kelso, Fuchs, and Pandya (1990), Kelso (1995), and Schöner and Santos (2001) have developed a large body of work of how to model movement phenomena with dynamical systems, emphasizing dynamic phenomena like phase transitions. Similarly, Kugler and Turvey (1987), Turvey (1990), and Sternad et al. (1996) investigated dynamic phenomena of biological movement, with a focus on rhythmic movements and distinction of order parameters and control parameters.

There is also extensive work on using coupled oscillators as models of biological behaviors, often discussed under the keyword “central pattern generators.” (Ijspeert, 2008) provides an extensive review on this topic.

In general, all of these previous approaches focused their modeling on rather specific systems and required a fair amount of user intuition. While drawing inspiration from these previous projects, our work goes beyond by suggesting a general approach to modeling with attractor systems, using the same form of representations in all applications, and using machine learning to adjust the model to observed data.

### 4.2. Robotics and Control Theory.

Potential field approaches create vector fields according to which a movement system is supposed to move. This idea has the same spirit as dynamical systems approaches in motor control. Deriving robot controllers from potential fields has a long tradition in robotics (Khatib, 1986; Koditschek, 1987; Tsuji, Tanaka, Morasso, Sanguineti, & Kaneko, 2002; Okada, Tatani, & Nakamura, 2002; Li & Horowitz, 1999). Potential fields represent attractor landscapes, with the movement goal acting as a point attractor. Designing such potential fields for a given behavior is often a hard problem, with few analytically sound approaches (Koditschek, 1987). Due to this lack of analytical tractability, some people have suggested recurrent neural network (Paine & Tani, 2004; Jaeger & Haas, 2004) or evolutionary methods (Ijspeert, Hallam, & Willshaw, 1999; Ijspeert, 2001) as a design strategy for nonlinear dynamical systems controllers.

From a more theoretical side, Bühler and Koditschek (1990), Rizzi and Koditschek (1994), Burridge, Rizzi, and Koditschek (1999), and Klavins and Koditschek (2001) developed a variety of control algorithms in the context of nonlinear dynamics that could be investigated both experimentally and analytically and that demonstrated very good performance. The idea of contraction theory by Lohmiller and Slotine (1998) also offers a promising set of tools for design principles with nonlinear dynamical systems.

Our suggested approach differs from such previous work in a variety of details. In contrast to the majority of previous work, we address high-dimensional dynamical systems rather than low-dimensional systems, as we work directly with families of trajectories to cover the global space. Other approaches seek global attractor landscapes out of general formulations, but it is hard to understand and manipulate the geometry of global attractors. Our trajectory-based thinking creates simpler although less flexible attractor landscapes but scales easily to higher dimensions and enables machine learning to shape the landscapes. We also develop our work more from a dynamical systems view that focuses on arbitrary vector fields, often interpreted as velocity or acceleration fields. In contrast, more control-theoretic approaches work directly with motor command generation, at the cost that motor commands do not generalize well to new situations because they are posture dependent in a robot. Our velocity or acceleration fields generalize more easily as they are purely kinematic in nature, but for robotics, we require a controller that transforms vector fields into control commands.

Probably one of the approaches that comes closest to ours is reservoir computing (Maass, Natschläger, & Markram, 2002; Jaeger & Haas, 2004). Reservoir computing uses a large recurrent network of randomly connected neurons with specific readout neurons for learning nonlinear dynamics, in particular, for time-series prediction. The parameters of the network (e.g., the synaptic weights) are kept constant except for the weights to the readout neurons, which are updated using linear regression given a desired output pattern. The approach has been used for predicting time series of complex dynamical phenomena with impressive results. In principle, the approach can learn a large class of dynamical systems, including point attractor and limit cycle dynamics. The main differences with our approach is that the underlying dynamics is much more complex than ours (with several hundreds of state variables), that reservoir computing does not offer proof of stability of learned attractors, and that it is less easy to incorporate feedback terms for online trajectory modulation. Therefore, it is not yet clear yet how reservoir computing can be used for controlling a robot (but, for some interesting steps in that direction, see Joshi & Maass, 2005; Wyffels & Schrauwen, 2009).

As a final remark, there is a rich literature on using optimization approaches to generate goal-directed movement (Kawato, 1996; Todorov, 2004), also discussed in some areas of reinforcement learning (Peters & Schaal, 2008; Theodorou, Buchli, & Schaal, 2010). Optimization approaches look for some general cost functions as a unifying priniciple to generate movement behavior, usually by offline optimization. As optimization can be rather sensitive to initial and final conditions of the movement and other environment factors, optimization does not necessarily generate movement primitives (i.e., stereotypical behavior), and every change of initial or final condition usually requires running the optimization again. Optimization is also mostly done on time-indexed trajectories. Thus, optimization approaches to movement generation are inherently different from our suggested approach to generate movement primitives with generalization and coupling properties. Nevertheless, it is possible to use optimization approaches in the context of our dynamical systems model to optimize the open parameters of the model (Peters & Schaal, 2008; Theodorou et al., 2010), which creates an interesting bridge between dynamical systems approaches and optimization approaches that exploits the best of both worlds (Schaal et al., 2007).

### 4.3. Dynamical Movement Primitives.

Our work on learnable dynamical systems originated from the desire to model elementary motor behaviors, called motor primitives, in humans and robots as attractor systems, an approach that has a long tradition in neuroscientific modeling (Kelso, 1995; Ijspeert, 2008). Initial work addressed imitation learning for robotics (Ijspeert et al., 2002a, 2002b, 2003). Subsequent investigation examined the approach for biped locomotion (Nakanishi et al., 2004), adaptive frequency modulation (Buchli, Righetti, & Ijspeert, 2006; Gams et al., 2009), reinforcement learning (Peters & Schaal, 2008; Kober & Peters, 2009), models for biological movement generation (Schaal et al., 2007; Hoffmann et al., 2009), and generating libraries of motor skills (Pastor et al., 2009; Ude, Gams, Asfour, & Morimoto, 2010). Various other projects in robotics created work that resembles our basic approach—the generation of kinematic movement primitives as attractor systems using basis function approaches to model the nonlinear dynamics (e.g., Billard & Mataric, 2001; Schaal et al., 2003; Billard, Calinon, Dillmann, & Schaal, 2008; Kulvicius, Ning, Tamosiunaite, & Worgötter, 2012).

A interesting variant was presented by Gribovskaya, Khansari-Zadeh, and Billard (2010) and Khansari-Zadeh and Billard (2010). In this work, the authors use gaussian mixture models (GMM) to directly learn nonlinear attractor landscapes for movement primitives from demonstrated trajectories as a state-space approach; they avoid any explicit or implicit timing system like our canonical system, and they do not include stabilizing dynamics as accomplished with the spring-damper system in our transformation system equations. As an advantage, they obtain a movement primitive representation that depends on only observable states, which can be considered a direct policy learning approach, as discussed in Schaal (1999). Such a representation can learn much more complex attractor landscapes, where the attractor landscape can strongly vary throughout the state-space. As a disadvantage, the authors have to spend significant numerical effort on a nonconvex optimization problem to ensure stability of their movement primitives, they have to address learning of mixture models in potentially rather high-dimensional and ill-conditioned spaces where the number of mixture components may be hard to determine, and they require many more data throughout the state-space to represent the attractor dynamics. At this point, movement recognition and periodic motion are not addressed, and it is not entirely predictable how their movement primitives generalize in areas of the state-space that have few or no data. In essence, this approach differs in a similar way from ours as state-space-based optimal control and reinforcement learning differs from trajectory-based optimization approaches (Peters & Schaal, 2008). State-space approaches have a more powerful representation but quickly degrade in high-dimensional settings. Trajectory-based approaches work well in very high dimensions and generalize well throughout these spaces, but they have less representational power. It is really up to a particular application which properties are beneficial.

The large variety of follow-up and related approaches to our initial work on dynamical movement primitives is one of the motivations to present in this letter the theory, insights, and a refined approach to learnable dynamical systems that we hope will continue to attract even more active research in the future.

## 5. Conclusion

This letter presented a general design principle for learning and modeling with attractor dynamical systems for goal-directed behavior, which is particularly useful for modeling motor behaviors in robotics but also for modeling biological phenomena. Using nonlinear forcing terms that are added to well-understood dynamical systems models, we can create a rich variety of nonlinear dynamics models for both point attractive and limit cycle systems. The nonlinear forcing term can be represented as an autonomous coupling term that can be learned with standard machine learning techniques that are linear in the open parameters. We demonstrated the properties of our approach, highlighting theoretical and practical aspects. We illustrated our method in various examples from motor control.

To the best of our knowledge, the approach we have presented is the first realization of a generic learning system for (weakly) nonlinear dynamical systems that can guarantee basic stability and convergence properties of the learned nonlinear systems and that scales to high-dimensional attractor systems. Besides the particular realizations of nonlinear system models presented in this letter, we believe it is of greater importance to highlight the design principle that we employed. This design principle seems to be applicable for many other nonlinear dynamical systems models, as well as technical applications and computational neuroscience.

Several points of our approach require highlighting. First, the proposed system of equations is not very complicated. We primarily make use of linear spring-damper differential equations, while nonlinearities are introduced with the help of standard kernel-based function approximators. Despite this simplicity, a rather powerful methodology can be developed to create nonlinear dynamical systems models.

Second, many of the individual properties of our approach can be accomplished with alternative methods. For instance, imitation learning or coding of observed trajectories is easily accomplished with classical spline methods (Sciavicco & Siciliano, 2000; Miyamoto et al., 1996), and temporal and spatial scaling of splines could be accomplished by appropriate manipulations of the spline nodes. However, splines are time-based representations that are not easily modulated online by external variables, and splines can have quite undesirable behaviors between spline nodes, particularly when the execution is perturbed. Another example is movement recognition, where we illustrated that a simple classification based on the parameters of the function approximator allowed rather successful classification. Of course, many other pattern recognition methods could accomplish similar or better performance. But what we wanted to emphasize is that the invariance properties of our parameterization of dynamical systems allow an elegant approach to identify a movement. Thus, not the individual properties of our approach but rather the ability to address many issues in one coherent and simple framework is what should be emphasized in our work.

Third, from the viewpoint of applying our work to the representation and learning of motor behaviors, we followed the less common approach of coding kinematic movement behaviors in our dynamical systems approach. An inverse dynamics tracking controller was used to convert the kinematic state variables into motor commands. Other approaches would prefer to directly code motor commands to represent motor behaviors. There is nothing in our approach that would prevent us from interpreting the output of our dynamical systems directly as motor commands, and reinforcement learning could be used to optimize such a representation (Theodorou et al., 2010). However, it is only due to kinematic representations (Bernstein, 1967) that we could make interesting use of the generalization properties of our dynamical systems. Generalizing motor commands in a similar way would not create useful behavior as the nonlinear dynamics of physical systems would significantly alter the desired behavior (Hollerbach, 1984). As tracking controllers are rather well understood by now (Sciavicco & Siciliano, 2000), many control-theoretic papers have equally turned to kinematic representations of motor behaviors (Rizzi & Koditschek, 1994; Chevallereau, Westervelt, & Grizzle, 2005). In this vein, understanding how to couple nonlinear kinematic planning behaviors with nonlinear controllers, and even a closed loop through both system, seems to be an important topic to address. Our approach can be seen as a step toward such goals.

## Acknowledgments

This research was supported in part by the European Commission grant AMARSI FP7-ICT-248311 and by National Science Foundation grants ECS-0326095, IIS-0535282,CNS-0619937, IIS-0917318, CBET-0922784, EECS-0926052, the DARPA program on Learning Locomotion, the Okawa Foundation, and the ATR Computational Neuroscience Laboratories. H.H. was funded by a grant of the German Research Foundation (DFG) HO-3771-1.

## References

*Motor control*(pp.

## Notes

^{1}

With low-dimensional, we refer to systems with less than about 100 degrees of freedom.

^{2}

The code can be downloaded from http://www-clmc.usc.edu/Resources/Software.

^{3}

Note that we borrowed the terminology *discrete trajectories* from the motor control literature (Schaal, Sternad, Osu, & Kawato, 2004) to denote point-to-point (nonperiodic or episodic) trajectories—trajectories that are not repeating themselves, as rhythmic trajectories do. This notation should not be confused with *discrete dynamical systems*, which denotes difference equations—those that are time discretized.

^{4}

As will be discussed below, many other choices are possible.

^{6}

We assume that the data triples are provided with the same time step as the integration step for solving the differential equations. If this is not the case, the data are downsampled or upsampled as needed.

^{7}

This input vector was used in Ijspeert et al. (2003).

^{8}

In recent years, significant neuroscientific research has suggested that movement recognition could be generated with a movement-generating representation (Rizzolatti & Arbib, 1998).

^{9}

See, for instance, Pastor et al. (2009).