800IEEE TRANSACTIONS ON ROBOTICS, VOL. 26, NO. 5, OCTOBER 2010

Task-Specific Generalization of Discreteand Periodic Dynamic Movement Primitives
Aleš Ude, Andrej Gams, Tamim Asfour, and Jun Morimoto

Abstract—Acquisition of new sensorimotor knowledge by imitation is a promising paradigm for robot learning. To be effective,
action learning should not be limited to direct replication of movements obtained during training but must also enable the generationof actions in situations a robot has never encountered before. This
paper describes a methodology that enables the generalization of
the available sensorimotor knowledge. New actions are synthesized
by the application of statistical methods, where the goal and othercharacteristics of an action are utilized as queries to create a suitable control policy, taking into account the current state of theworld. Nonlinear dynamic systems are employed as a motor representation. The proposed approach enables the generation of a wide
range of policies without requiring an expert to modify the underlying representations to account for different task-specific featuresand perceptual feedback. The paper also demonstrates that the
proposed methodology can be integrated with an active vision system of a humanoid robot. 3-D vision data are used to provide query
points for statistical generalization. While 3-D vision on humanoid
robots with complex oculomotor systems is often difficult due tothe modeling uncertainties, we show that these uncertainties canbe accounted for by the proposed approach.
Index Terms—Active vision on humanoid robots, humanoid
robots, imitation learning, learning and adaptive systems.

I. INTRODUCTIONLEARNING of behaviors that can be applied to solve a giventask, regardless of the current configuration of the external

Manuscript received September 4, 2009; revised February 22, 2010 and July
1, 2010; accepted August 4, 2010. Date of current version September 27, 2010.
This paper was recommended for publication by Associate Editor A. Ijspeert
and Editor J.-P. Laumond upon evaluation of the reviewers’ comments. This
work was supported in part by the European Union Cognitive Systems Project
PACO-PLUS under Grant FP6-2004-IST-4-027657. The first author (A. Ude)
was also supported by the National Institute of Information and Communications Technology within the JAPAN TRUST International Research Cooperation
Program.
A. Ude is with the Department of Automatics, Biocybernetics, and Robotics,Jožef Stefan Institute, 1000 Ljubljana, Slovenia, and also with the Computational Neuroscience Laboratories, Advanced Telecommunications ResearchInstitute International, Kyoto 619-0288, Japan (e-mail: ales.ude@ijs.si).
A. Gams was a Visiting Researcher at Computational Neuroscience Laboratories, Advanced Telecommunications Research Institute International, Kyoto
619-0288, Japan, in summer 2009. He is now with the Department of Automatics, Biocybernetics, and Robotics, Jožef Stefan Institute, 1000 Ljubljana,
Slovenia (e-mail: andrej.gams@ijs.si).
T. Asfour is with the Institute for Anthropomatics, Karlsruhe Institute ofTechnology, 76131 Karlsruhe, Germany (e-mail: asfour@kit.edu).
J. Morimoto is with the Computational Neuroscience Laboratories, Advanced
Telecommunications Research Institute International, Kyoto 619-0288, Japan(e-mail: xmorimo@atr.jp).
This paper has supplementary downloadable material available at
http://ieeexplore.ieee.org, provided by the authors. This includes seven QuickTime format movie clips, which show the results on real robots. The completematerial is 45.2 MB in size. Contact ales.ude@ijs.si for further questions about
this work.
Digital Object Identifier 10.1109/TRO.2010.2065430

world, is a difficult problem because the search space that needs
to be explored is potentially huge [1]. The size of the search
space depends both on the number of degrees of freedom of therobot and on the objects involved in the action. Furthermore,
external objects can affect the search space indirectly. To overcome problems arising from high-dimensional and continuous
perception-action spaces, it is necessary to guide the search process. One of the most successful paradigms that can be used forthis purpose is imitation or robot programming by demonstration[1], [2].
Robot programming by demonstration requires the acquisition of example trajectories, which can be captured in various
ways. Motion capture techniques based on optical- or magnetictracking devices were applied successfully to replicate complex
movements on humanoid robots such as dancing, which would
be difficult to program manually [3]–[6]. Computer vision techniques that aim to capture human motion without special markers are slowly getting matured [7], but the application of specialized tracking devices is still preferred in most of the recent works
in order to avoid the pitfalls of computer vision approaches. Alternatively, a robot can be physically guided through the desired
trajectory, which is recorded proprioceptively. This method requires that the robot is back-drivable [8], [9] or can compensate
for the influences of external forces [10]–[12]. The advantage ofkinesthetic guiding is that the movements are recorded directlyon the learning robot and do not need to be first transferred froma system with different kinematics and dynamics.
A simple reproduction of the trained movements is not useful for problems that involve the manipulation of objects because in such tasks the observed movements need to be adapted
to the current state of the 3-D world. Already, early researchon learning from demonstration has stressed the importance oftask segmentation and the extraction of meaningful action units[13]–[15].
The main interest of this paper is in the generalization of theavailable action knowledge regardless of how this knowledge
was initially acquired. A methodology that can adapt single
trajectories obtained by imitation was proposed by Miyamoto
et al. [16]. They made use of a spline-based representation forthe desired trajectory, which they referred to as via points. By
monitoring its performance, the robot was able to continuously
adapt the via points until it could play a fairly difficult Japanesegame of kendama or execute tennis serves. This type of approach has later been studied in the frame of modern reinforcement learning theories [9]. One of the main issues is the
speed of convergence, which is important for online applications of reinforcement learning and is an active topic of research[17], [18].

1552-3098/$26.00 © 2010 IEEE

UDE et al.: TASK-SPECIFIC GENERALIZATION OF DISCRETE AND PERIODIC DYNAMIC MOVEMENT PRIMITIVES

Numerous representations have been proposed to encode
robot trajectories for learning. Spline-based representations likethe one utilized in [16] have been employed in robotics in many
contexts [19], but their explicit dependency on time can becumbersome [20]. Hidden Markov models (HMMs) are anotherpopular methodology to encode and generalize the observed
trajectories [21]–[23]. It has been shown that HMMs can beused effectively for motion and action recognition [23] andto determine which control variables should be imitated andhow [22]. Kulić et al. [24], [25] extended these works by showing how to use HMMs to automatically cluster large databasesof movements into the constituent motor primitives. Gaussian
mixture models [26] provide another motor representation withsome advantages over HMMs during reproduction of the learned
movements. Calinon et al. [26] also demonstrated how movements can be generalized by taking into account different analytically specified constraints between hands and objects. Yamane
et al. [27] proposed a graph-theoretical approach that organizes
the captured data in a large motion database. In their work, theproblem of motion generation becomes a search problem.
A fundamentally different approach to motion representation
based on nonlinear dynamic systems as policy primitives was
proposed in [28] and [29]. The resulting control policies weretermed as dynamic movement primitives (DMPs). DMPs are
based on systems of second-order differential equations, which
encode the properties of the desired motion. Ijspeert et al. [28],[29] developed equations for periodic and discrete movements
and demonstrated the learning of tasks such as tennis strokes
and drumming. One of the most important advantages of DMPs
is the ability to take into account perturbations and to include
feedback terms. Feedback terms can be added to change thetiming [20] and/or avoid some areas of the workspace. Anotherapproach based on the idea of dynamic systems of Ijspeert et al.
is described in [8], where the final trajectory is generated as alinear combination of attractor dynamics and Gaussian mixture
models that were learned from example trajectories.
A. Contribution of This Paper and Related Work
The main purpose of this paper is to propose and experimentally evaluate a methodology for generalization of example trajectories to new situations that were not observed during
training, using a representation suitable for robot control. To thisend, every example trajectory is associated with parameters that
describe the characteristics of the task, typically its goal, andserve as query points into an example database. For instance, in
our experiments, we studied the problem of ball throwing toward a specified target. In these experiments, we first acquireda set of example robot trajectories that resulted in ball throws
toward different targets. The task of generalization is to synthesize throws toward any given target within the training space.
The parameters specifying the goal of the task are, in this case,the target positions.
Our work was inspired by motor-tape theories, in which
example movement trajectories are stored directly in memory [30], [31]. Initial implementations of this theory were basedon lookup tables [32] that cannot generalize. More modern ver-801

sions interpret the task of generalization as a problem of fitting
a multivariate function to the stored training data, where the fitted function must generalize to new situations [30]. Regression
methods based on local weighting of training data at execution
time can be applied for this purpose. They have proven to be
sufficiently efficient for online control [33].
We utilize nonlinear dynamic systems as a basic motor representation. Other researchers have shown that by changing theunderlying differential equations, DMPs can be modified in several ways to account for various perturbations that might occur
during the execution of the task. For example, it has been suggested to add terms that enable obstacle avoidance to the basic
DMP equations [34], [35]. Such modifications are very useful,but they are necessarily designed for one particular issue only.
For each new problem, an expert must redesign the underlying
dynamic system, which is unsatisfactory for robots that need tosolve new problems every day. The approach proposed in thispaper enables the generalization of DMPs to new situations using the available training movements and the goal of the task,
which can normally be specified in a natural way. Since the generalized trajectories are encoded as standard DMPs, we are still
able to apply analytical modifications of the underlying dynamic
system at execution time in order to account for unforeseen perturbations, which are not part of the training data. In addition,
Section III-C demonstrates that the proposed approach is suitable for online integration with a humanoid robot’s active vision.
It is also possible to generalize the example movements to
new situations by specifying constraints and weighting themwith respect to their importance for the task [22]. For example,
if the initial positions on the training trajectories are found to
be less important, the robot can modify its motion to reach thefinal position from a different initial position at the expense of
fidelity of reproduction. Systems like this require an analyticalmodel to describe the constraints. The approach proposed in thispaper requires no physical models of the task to be specified;generalization is done in a statistical way.
Another methodology for taking into account external variables was proposed in [36], where—instead of synthesizing the
optimal DMP parameters at execution time as in our approach—the basic DMP equations were modified by coupling them withthe perceptual input. This coupling provides additional parameters that can be used to learn a more globally valid representation
as opposed to the local one generated by our system. As shownin [36], the external variables have a similar role as query points
into the example database in our approach. It is, however, moredifficult to ensure that such a system learns all interdependencies
between the perceptual input and the DMP parameters because,in general, it is unknown how many additional degrees of freedom are needed. This problem is avoided by using local models.
Our work is also related to the theories in which generalization
takes place through a combination of motor primitives [37].
An architecture that linearly combines internal models, which
can be regarded conceptually as motor primitives, based on the
output of a responsibility estimator has been proposed in [38]and tested in simulation [39]. Work under progress has beenreported in [18], where motor primitives are linearly combined
based on the output of a gating network. Such an approach takes

802

IEEE TRANSACTIONS ON ROBOTICS, VOL. 26, NO. 5, OCTOBER 2010

generalization to a higher abstraction level away from the training data, which reduces the dimensionality of the problem but
makes accurate generalization more difficult.
II. ACTION GENERALIZATION USING DYNAMIC MOVEMENTPRIMITIVES AND LOCAL WEIGHTING
Let us assume that we have a set of example trajectories
together with the parameters characterizing the task
Z = {ydk (tk ,j ), ẏdk (tk ,j ), ÿdk (tk ,j ); qk |k = 1, . . . , M
j = 1, . . . , Tk } (1)where ydk (tj ), ẏdk (tj ), and ÿdk (tj ) are the measured positions,velocities, and accelerations on trajectory k, M is the numberof examples, and Tk the number of sampling points on eachtrajectory. qk ∈ Rn are the parameters describing the task ina given example situation and are usually related to the goalof an action. In our system, they are used as query points into
a database of example trajectories. Such data can be obtained
either by kinesthetic guiding or from human demonstrations.
The trajectories can be specified either in joint or task space. In
a DMP formulation (see Appendix A for a brief review of the
DMP representation), every degree of freedom is described byits own dynamic system but with a common phase to synchronize them. Indexing of the degrees of freedom is omitted from
(1) for clarity. The issue is how to generate a DMP specifying
a movement for every new query point q, which, in general,
will not be one of the example queries qk . As explained inAppendix A, DMPs are specified by parameters w, τ (orΩ = 1/τ in case of periodic movements), and g. Thus, we needto learn a function
G (Z) : q −→ [wT , τ, g]T .

(2)

Examples Z specified in (1) are the data used for learning. While
the problem of learning of g and τ has also been studied in [18],here, the complete learning problem is treated.
In general, the functional relationship between q and[wT , τ, g]T given a set of examples Z is unknown. In mostcases, it is difficult to find a global model that provides good approximation for the function G(Z). We therefore avoid global
model identification and rather apply regression techniques togeneralize the movements. Due to significantly different sizesof datasets involved in the calculation of parameters w on the
one hand, and g and τ on the other hand, we propose to employ
different methods to estimate them. More specifically, we applied locally weighted regression (LWR), which is a regression
method that fits local models to nearby data [40], for the estimation of w. LWR has a lower computational complexity thanmany other nonparametric regression methods including Gaussian process regression (GPR) [41]. On the other hand, due toits high accuracy, we utilized GPR [42] to estimate g and τ .
Another advantage of GPR is that there is no need to determine
the place and number of basis functions. However, its computational cost can be prohibitive especially in the case of online
learning, where the dataset gradually increases in size. This was
not a problem when learning g and τ in our experiments.

Our case studies include tasks such as reaching, ball throwing,
and drumming. In the case of reaching movements, the goalof an action (or query point) was simply the final reachingdestination in Cartesian space, whereas the trajectories and,
consequently, the attractor points of the DMPs were defined injoint space. Thus, the query and attractor points were connected
through the robot kinematics. In other cases, the queries wereless directly associated with the parameters of the DMP. For
example, in the case of ball throwing, the goal was characterized
by the position of the basket into which the ball should bethrown. This position is not directly encoded in the parametersof the discrete DMP, as given by (25) and (26), shown later.
As example of periodic movements, we studied drumming on acymbal and a drum, where the height at which the cymbal was
mounted varied. The robot adapted its drumming movements
to the varying height of the cymbal. In this case, the heightof the cymbal was used as query point, which is again notdirectly encoded in the periodic DMP equations (28) and (29),shown later. In summary, the query points normally originate
from an intuitive characterization of the task, but the functional
relationship between them and the DMP parameters may not be
straightforward.
Note that G(Z) becomes a function only by constraining thesolution trajectory to be as similar as possible to the example
trajectories. For example, there are many different ways of how
to throw a ball into a basket at a certain location. The relationship between the basket positions (query points) and DMPparameters as given in (2) becomes a function by requesting thatthe generalized throwing movements are similar to the example
throws. The similarity criterion is embedded into the regression
process in our practical implementation.
In the following, we provide a methodology to generate new
DMPs for situations that are not part of the example database. To
fully specify a DMP that defines a movement in a new situation,
we need to estimate the shape parameters w, goal g, and, respectively, the time constant τ in case of discrete movements andfrequency Ω in case of periodic movements. Section II-A first
explains the generation of DMPs using only one training trajectory. We continue the description of the estimation of parameters
w for discrete (see Section II-B) and periodic movements (seeSection II-C). Section II-D deals with the estimation of other
DMP parameters that vary from example to example (g and τ orΩ). The rest of the DMP parameters (αz , βz , and αx ) are fixed
and are determined so that the convergence of the underlying
dynamic system is guaranteed.
A. Reproduction From Single Demonstration
With the DMP representation of Appendix A, the trajectoryof any smooth movement can be approximated by adapting theparameters wi of (23) and (27). The system of two first-orderlinear equations (25) and (26) in the case of discrete movementsand (28) and (29) in the case of periodic movements can be
rewritten into one second-order equation by replacing z with τ ẏin (25) and (28), respectively
τ 2 ÿ + αz τ ẏ − αz βz (g − y) = f(3)

UDE et al.: TASK-SPECIFIC GENERALIZATION OF DISCRETE AND PERIODIC DYNAMIC MOVEMENT PRIMITIVES

with f defined as in (23) and (27). All of the above-mentioned
equations are shown later in Appendix A. The formula τ =1/Ω is used in the case of periodic movements. Note that time
constant τ must be the same for all degrees of freedom. In our
experiments with discrete movements, we used τ = tT , wheretT is the duration of the training movement. On the other hand,the attractor points g vary across the degrees of freedom. They
can be extracted directly from the data as g = yd (tT ) in thecase of discrete movements and g = j yd (tj )/T in the caseof periodic movements. Writing
F (tj ) = τ 2 ÿd (tj ) + αz τ ẏd (tj ) − αz βz (g − yd (tj ))F (t1 )
w1
f =  ... , w =  ... (4)wN
F (tT )
we obtain the following system of linear equations:
Xw = f(5)

which needs to be solved to estimate the remaining parametersof a DMP that describe the desired motion. In the case of discrete
movements, we have
Ψ1 (x1 )ΨN (x1 )
x
x
.
.
.


1
1
N
N


i=1 Ψi (x1 )i=1 Ψi (x1 )




...
...
...
X=
 (6)




Ψ1 (xT )ΨN (xT )xT . . . NxT
N
i=1 Ψi (xT )i=1 Ψi (xT )and in the case of periodic movementsΓ1 (φ1 )ΓN (φ1 ) 
.
.
.


N
N


i=1 Γi (φ1 )i=1 Γi (φ1 ) 



...
...
...
X = r.

Γ1 (φT )ΓN (φT ) . . . NN
i=1 Ψi (φT )i=1 Γi (φT )(7)Here, xi and φi are obtained by respective integration of (24)and (30), shown later. The parameters w can be calculated fromthe solution of the above systems of linear equations in a leastsquares sense.
While it is common to estimate discrete movements using single demonstrations and the aforesaid batch approach, periodic
movements are more often estimated recursively by repeating
the movement a number of times. Note that the integration of(28)–(30), shown later, requires that the frequency of movementΩ is known. Unlike time duration, this frequency is not directly
observable and must be estimated. It has been suggested to replace (30), shown later, with the system of adaptive frequency
oscillators [43], [44] (see also Appendix B) to automaticallydetermine the frequency during demonstration. The waveform
estimation can then be realized by solving the equation system(5) using recursive least squares with a forgetting factor
Pj −1 xj xTj Pj −1
1(8)
Pj −1 −Pj =λλ + xTj Pj −1 xj
wj = wj −1 + (fj − xTj wj −1 )Pj xj

(9)

803

where P0 = I, w0 = 0, fj = F (tj ), xj is the M -dimensional
column vector associated with the corresponding row of thesystem matrix X, 0 < λ ≤ 1 is the forgetting factor, and thefinal weights are given as w = wT .
It is important to note that with this system, the frequency
estimation should be done simultaneously with the recursive
estimation of the form parameters w; at each time tj , first thecurrent Ω is estimated by integrating (31)–(33), shown later. The
estimated Ω is then used to calculate the target values F (tj ) of(4), where again, τ = 1/Ω. Finally, the new estimate wj is calculated using the recursion (8) and (9). Thus, training can beorganized as a coaching process; the coach demonstrates the
movement, the robot simultaneously executes the movement
using the currently estimated parameters, and the coach stops
the demonstration once he/she is satisfied with the robot’s performance. This approach is similar to the idea proposed in [45],where the transfer of human motor skills to the robot is supported by a training system that keeps the human instructor andthe robot in a real-time control loop. Note that (8) includes theforgetting factor λ; thus, later training data have more influenceon the estimated parameters than earlier data.
In the aforesaid equations, αx , αz , and βz are constant. They
are set so that the convergence of the underlying dynamic systemis ensured [20]. Unlike in some other works where the parameters wi are estimated independently of each other [28], [29],we apply a full linear system (5) to estimate w. In this way, we
can approximate trajectories more accurately because we cantake into account the interplay between the neighboring basis
functions Ψi of (23) and (27), shown later, which benefits generalization. Note that the separate estimation of {wi } has its
advantages, especially in the presence of noise when overfitting
can become a problem or when {wi } are used for classification [29]. However, here the example trajectories are trained
movements suitable for execution on a robot. They are therefore
relatively noise-free, which alleviates the danger of overfitting.
It was therefore sufficient in our experiments to determine theparameters ci and hi by setting the distribution pattern1 andincreasing N until the desired reconstruction accuracy on all example trajectories was achieved. More advanced methods [46]
could be applied to determine N, ci , and hi if overfitting becamea problem.
B. Generalization of Discrete Movements
In one-shot learning of Section II-A, all data are relevant for
the estimation of shape parameters w. Such global optimization
approaches are problematic for generalization from multiple examples because, in general, no global task models are available.
Control policies that solve the task in situations, which are very
different from the current one, do not carry much information
about the optimal policy in the current situation. Taking theexample of ball throwing, it is reasonable to assume that the
trajectories associated with target positions close to the current

1 For a given N , c

i = 1, . . . , N .

i −1i = exp(−α x N −1 ), h i = (c2

i + 1 −c i )2, h N = h N −1 ,

804IEEE TRANSACTIONS ON ROBOTICS, VOL. 26, NO. 5, OCTOBER 2010

target are more relevant than example trajectories associatedwith more distant targets.
The synthesis of shape parameters w is based on the sampled
trajectory points included in the training dataset (1). Since thenumber of sampling points is usually very large, it is beneficial
to apply a method with relatively low computational complexity,
such as, for example, locally weighted regression. In an LWR
setting and for a given query point q, the optimal parameters
can be calculated directly from the available data by minimizing
the objective function
L(Ξ(qk , w), fk )K(d(q, qk )).

C(w; q) =(10)kBased on (5), local models are characterized byL(Ξ(qk , w), fk ) = Xk w − fk 2(11)where Ξ(qk , w) = Xk w. Xk needs to be calculated as in (6) for
discrete or in (7) for periodic movements. fk can be computed
as in (4) but using the generalized τ and g (see Section II-D).
We need to minimize the objective functionMXk w − fk 2 K(d(q, qk ))

(12)

k =1

with respect to w. Here, K is the kernel function and d is themetrics in the space of query points q.
There are many possibilities to select the weighting kernel
K [40]. We chose the tricube kernel
K(d) =(1 − |d|3 )3 ,if |d| < 10,otherwise(13)

because this kernel has finite extent and a continuous first andsecond derivative, which means that the first two derivatives ofthe prediction (as a function of query points) are also continuous.
The finite support of K reduces the computational complexity ofthe optimization problem (12) because the example trajectories
for which K vanishes do not influence generalization. In thisway, we reduce the size of the system matrix associated with theobjective function (12). As discussed in [40] and [47], the choiceof weighting function is rarely critical for the performance of
locally weighted regression. The selected kernel performed wellin our experiments.
K and distance d in the space of query points determine howmuch influence each of the example movements has on the final
estimate of the control policy. The influence of each example
movement should diminish with the distance of the query point
q from the data point qk . A standard weighted Euclidean distance can be used when query points are given in Euclideanspaced(q, qk ) = D(q − qk ), D = diag(1/ai ), ai > 0.

(14)

Ideally, D should be determined so that the error in the execution of the task is as small as possible. To reduce the resulting
problem to a 1-D optimization problem while still accounting
for the possibly varying spacing across dimensions of the query
point space, we definedai = c ∗ max

min {|qj,i − qk ,i }.

j =1,...,M k =1,...,M(15)

In this way, D is fully specified by a parameter c ∈ R, which
we calculated by minimizing the validation set error [40] (seealso Section III-D). Such an approach is sufficient if the data areapproximately equally distributed along the coordinates axes,
as it was in all our experiments. The validation set-error method
evaluates the difference between the predicted output and theobserved value in the validation set, where the data in the validation set are not used for training. The minimization of thevalidation-set error gave better results on regularly spaced training data than leave-one-out cross-validation method because it
could better take into account the in-between query points. In
all of our experiments, the data were approximately regularly
spaced and we obtained the optimal value c ≈ 2.2. With such c,
the number of example trajectories with nonzero K(d(q, qk ))was about 4n .
The computational complexity
 of solving the least-squares
system (12) is O(N 2 T ), T ≤ M
k =1 Tk and thus increases linearly with the number of data points considered by LWR andquadratically with the number of radial basis functions used in(23) and (27), shown later, respectively. Due to our choice of
weighting kernel K, we normally have K(d(q, qk )) = 0 formany k. Moreover, by cutting the support of basis functions(23) and (27), shown later, once their value falls below a certain
threshold, matrices Xk become sparse as well. The quadratic
dependence on the number of basis functions is not a problem
because this number is generally much lower than the numberof data points. In most of our experiments, there were around
10 000–50 000 data points and 25–50 basis functions for DMPs.
These facts make computational complexity sufficiently low toresolve the least-squares problem (12) online using standard
methods from sparse-matrix algebra.
C. Generalization of Periodic Movements
In Section II-A, we describe a method that estimates the parameters of a DMP to reproduce a demonstrated periodic trajectory. The simultaneous and recursive estimation of frequency Ωand form parameters w allows the robot to reproduce the demonstrated motion immediately, which is important for coaching. To
generalize the learned movements to new situations, we store
the estimated frequencies and the sampled movements from thelast few movement periods, which produced a good movementon the robot as judged by the teacher. In our experiments, weused the last five periods.
The aforementioned training process makes the following data available for generalization purposes: Trajectory data points within the last few periods of motion
j = 1, . . . , Tk }{ydk (tk ,j ), ẏdk (tk ,j ), ÿdk (tk ,j )|k = 1, . . . , M,
and the associated frequencies Ωk . Each trajectory is also
associated with query parameters qk ∈ Rn . Given a newquery point q, we could now directly the minimize objectivefunction (12) to calculate the new parameters w. However,since the data were acquired by recursive least squares with
forgetting factor λ, it is necessary to generalize recursively as
well. In this way, we ensure that the generalization process
discounts earlier training data in the same way as the coachingprocess. To achieve this, we need to parse all of the sampled

UDE et al.: TASK-SPECIFIC GENERALIZATION OF DISCRETE AND PERIODIC DYNAMIC MOVEMENT PRIMITIVES

Fig. 1.

805

Procedure for the generalization of periodic movements.

trajectory points in a proper order. For this purpose, it isnecessary to maintain a separate phase information for every
example trajectory to account for different frequencies during
the generalization process. The resulting procedure shown in
Fig. 1 ensures that the data from all of the trajectories are parsed
and discounted across example trajectories in the same way as
during coaching.

Fig. 2. Training and generalization of goal-directed actions (for one dimensionof the space in which trajectories are defined).

lations can be done off-line using, for example, the Cholesky
decomposition. Note that by writingz = [K(Q, Q) + σn2 I]−1 y

D. Estimation of Attractor Points, Timings, and Frequencies
Unlike αx , αz , βz , and N , which are kept constant acrossthe example trajectories during generalization, time constant τin the case of discrete movements and frequency Ω in the case of
periodic patterns as well as the attractor points gi , i = 1, . . . , D,
where D is the dimension of the space in which the motion is
specified, change from example to example. We extract theseparameters from the training data (see Section II-A), whichenables us to estimate the function that transforms query points
q into τ , Ω, and gi , i = 1, . . . , D, directly.
Among various paradigms that could be used for this purpose,
Gaussian process regression (GPR) has proven to be especiallyeffective. It is a Bayesian regression method that provides a
predictive distribution. GPR exhibits good generalization performance and the predictive distribution can be used to measure the uncertainty of the estimated function. Complexity ofthe model can be adaptively changed according to the newly
acquired sampled data. There is no need to explicitly place additional basis functions and parameters. Simultaneously, GPR
can avoid the danger of overfitting. It has been shown that thistechnique outperforms other regression methods on some relevant robotic problems such as estimating inverse dynamics of aseven degrees of freedom robotic arm [42]. Since the numberof training data points is much lower when estimating τ , Ω, andgi than when estimating the form parameters w (now it is thesame as the number of training trajectories, before we had toconsider all of the sampled points on the example trajectories),
the computational complexity is not a problem and we selected
GPR for the estimation of these parameters. A brief review of
GPR is given in Appendix C.
With GPR, new estimates are calculated using (36) from the
Appendix. The most computationally expensive part of this formula is the calculation of [K(Q, Q) + σn2 I]−1 , but since thismatrix depends only on the training data, the necessary calcu-

∗

∗

(16)

, Ω̄ , and ḡi∗ , associated(36) and the new parameters ȳ , i.e., τ̄
with the query Q∗ = q∗ can be written as∗Mȳ ∗ =k(q∗ , qk )zk .

(17)

k =1

Thus, similarly as in LWR technique of (12), the data areweighted based on the distance between training query points
and the current query point; hence, nearby training points influence the result more. To generate a new movement, the robot
is given a desired query point. Using GPR, gi and τ (or Ω) for
this query are calculated. The optimal parameters w are thenestimated using the algorithms of Section II-B or C. The sketchof the complete training and generalization procedure is given
in Fig. 2.
III. EXPERIMENTS
We conducted a number of experiments to demonstrate theusefulness of the proposed approach. The following tests weredone for discrete movements: 1) a reaching study in simulationand on a real robot to evaluate the accuracy of reproduction
when learning discrete movements, 2) experiments with a fullsize humanoid robot that connect the proposed approach withactive vision and also demonstrate grasping, and 3) a simulated
and a real-world ball-throwing study that confirm that also more
dynamic aspects of the task can be generalized. To validate
the generalization of periodic movements, we conducted 1) asimulated periodic-pattern-movement study and 2) a real-world
drumming experiment.
Two humanoid robots were used to validate the proposed
approach: a small humanoid robot HOAP-3 built by Fujitsu
Automation and a full-size humanoid robot CB-i built in collaboration between Sarcos, ATR, and JST [48]. Both robots enable
the acquisition of trajectories via kinesthetic guiding. CB-i is,806IEEE TRANSACTIONS ON ROBOTICS, VOL. 26, NO. 5, OCTOBER 2010

TABLE IIERRORS IN REACHING MOVEMENTS (IN DEGREES AND CENTIMETERS,RESPECTIVELY) GENERATED BY A SINGLE DMP, WHICH WAS TRAINEDTO REPRODUCE ONE OF THE EXAMPLE TRAJECTORIES

Fig. 3. Forty-five example trajectories in joint space calculated via inverse
kinematics from Cartesian space minimum-jerk trajectories (left) and the final
reaching points in Cartesian space (right). The Cartesian space positions areused as query points. The outer rectangle shows the full training area; the innerrectangle shows the inner training area where training data are available fromall sides of the query points.
TABLE IERRORS IN REACHING MOVEMENTS (IN DEGREES AND CENTIMETERS,RESPECTIVELY) SYNTHESIZED BY LWR AND GPR

probably due to the fact that we do not need to learn a full
inverse kinematics of the robot. Hence, the inverse kinematics
learning is limited to a small portion of the workspace coveredby query points, and the learning procedure is less likely to
encounter the potentially problematic areas.
The errors in Tables I and II were computed by comparing generalized joint trajectories ỹ(tj ) (calculated by integrating (24)–(26), shown later, with the analytically computed minimum-jerk trajectories y(tj ), both specified in the
robot joint space. The average (18), shown below, and maximum error (19), shown below, on the trajectories were
estimated asTerroraverage =however, a much more capable system that enabled us to study
the learning of control policies integrated with a full-fledged
active humanoid vision. In our throwing experiments, we useda seven degree-of-freedom Mitsubishi PA-10 robotic arm.
A. Simulated Reaching StudyIn the first computational study, we examined how wellCartesian minimum-jerk trajectories can be generalized. Minimum jerk trajectories are often used in robotics because they
resemble human-reaching trajectories [49]. For training, we generated 45 Cartesian minimum-jerk trajectories, which were converted into joint space of a planar 2R robot (see Fig. 3, left).
The final end-effector positions on the trajectories were used
as query points. The training query points were distributed uniformly with spacing of 0.1 m in a rectangular area with corners at
locations (0.2, −0.5) and (0.6, 0.3) m (Fig. 3, right). Joint velocities and accelerations were computed analytically.
As described in Section II, we applied LWR to synthesize
the form parameters w of the DMP, while GPR was utilized tolearn the function from query points (desired goal position in
Cartesian space) to attractor points (the final joint configuration)and to time duration. Thus, here the relationship between query
points and attractor points is given by the standard inverse kinematics. While problems with nonconvexity of the movement
space2 may arise in less-constraint situations, we did not observe such problems in our experiments. This observation is
2 These problems can be addressed by biasing [50] or re-weighting [51] thedata.

1ỹ(tj ) − y(tj )
T j =1errorm ax = max ỹ(tj ) − y(tj ).
j =1,...,T

(18)
(19)

Query points, at which the generalized and analytical trajectories were computed, were sampled across the complete training
space. Results in Table I demonstrate that minimum-jerk movements can be generalized with high precision. The generalized
trajectories accurately approximate the spatial course of movement and the final configuration. Since it can be expected thaterrors are larger on the boundary of query points used for training, we estimated the errors both within the full rectangular area
enclosed by all query points of Fig. 3 and in a more limited
inner training area enclosed by query points situated at least one
query point away from the boundary points. As expected, theerrors were significantly smaller when testing was limited to theinternal training area.
Table II shows that representation with only one DMP is toorough for a precise movement reproduction. While the final
position could be reached accurately due to the properties of
DMPs, the trajectory-reproduction accuracy (18) is worse by anorder of magnitude compared with the precision of the proposed
approach. The columns describing the Cartesian-space error inboth tables show that GPR is successful at estimating the inverse
kinematics of the robot in a limited subset of the workspace. The
error is larger in Table II than in Table I because the execution
was stopped at time τ . The DMP obtained by generalization
arrives at the goal within the specified time, whereas the DMP
that was trained to reproduce one of the example trajectories
needs more time to reach a modified attractor point.

UDE et al.: TASK-SPECIFIC GENERALIZATION OF DISCRETE AND PERIODIC DYNAMIC MOVEMENT PRIMITIVES

807

Fig. 4. Image sequence showing the teaching of reaching trajectories to
HOAP-3 humanoid robot with kinesthetic guiding. Teaching of reaching movements to CB-i was done in a similar way.
Fig. 6. Evaluation of a real-world reaching experiment. Both graphs show theresults of generalization when query points smoothly change along a dimensionof the query-point space. The generalized joint trajectories (bold) resemble
the trajectories of the nearby demonstrations (thin). The left graph shows thechanging of the query point in y direction and the right graph in the z direction.
Note that the y axis is reversed to preserve similarity to Fig. 5. Query points
associated with training trajectories are presented with full circles and the query
points for generalized movements with empty circles.

points in both the y and z directions, while x remains roughly
the same, just as in the demonstrations. Thus, we can confirm the simulation results from Section III-A in a real-world
experiment.
C. Integration With Active Vision and Grasping

Fig. 5. Query points and example trajectories for the reaching experiment on
HOAP-3. The goals, i.e., final destinations, of the demonstrated movements aremarked with circles. The robot’s arm with the shoulder (RGH) and elbow (REB)
joints is also shown.

B. Real-World Reaching
We also tested the proposed approach for the generalization ofreaching movements on a humanoid robot HOAP-3. By meansof kinesthetic guiding (see Fig. 4), we recorded a set of 25
reaching movements, which were used to build a library of
example movements. The example trajectories in task space
and the associated query points are presented in Fig. 5. All
movements roughly originated from the same starting position(P0 ) and ended on a grid roughly in the coronal plane of therobot, ∼0.15 m in front of it.
The robot’s forward kinematics was used to obtain the final end-effector’s positions in Cartesian space on the trajectory.
Gaussian process regression (see Section II-D) was applied toestimate the function from these positions, which served as
query points, to the corresponding joint angles at the final configuration on the trajectory. Since the robot is redundant withrespect to the task, GPR also learns to resolve redundancies
within the demonstrated, limited subset of the workspace using
the same redundancy resolution strategy as demonstrated during
training. Four joints (three in the shoulder and one in the elbow)
were used in this experiment.
Fig. 6 shows the generalized Cartesian trajectories projected
onto the yz plane. The generalized trajectories, which are presented in bold, resemble the example reaching trajectories associated with the nearest demonstrations, which are presented bythin solid lines. The generalization procedure continuously transitions between example trajectories when changing the query

In our next experiment, we focused on the integration of
the proposed approach with active vision and grasping, which
demonstrates high performance that can be achieved by the developed system (see Figs. 7 and 8). The role of active vision isto provide parameters characterizing the task, i.e., query points.
In the case of reaching and grasping, the task is characterized by
the position of an object to be grasped. The use of active vision
is essential if the robot is to find and grasp objects in a naturalway. In this experiment, we used a full-size humanoid robot
CB-i. Like in the case of HOAP-3, the training data were obtained by guiding the robot through 25 example trajectories.
Besides reaching toward the desired positions, the acquired trajectories also avoid the table. Unlike in the case of HOAP-3, thepositions of the target object were acquired by the robot’s own
visual system.
CB-i’s oculomotor system has seven degrees of freedom(three in the neck and two in each eye), which ensures flexibility when the robot needs to find and direct its view toward
new objects. To compute 3-D data in body coordinates by stereo
vision, a robot with such a visual system must continuously update the position and orientation of the cameras, which dependon the current joint configuration of the robot, in the robot-bodycoordinate frame. The appropriate estimation and calibration
procedures are described in [52]. This work showed that evenwith carefully designed calibration procedures, it is difficult toestimate all necessary coordinate-frame transformations withhigh accuracy. As a consequence, we cannot rely on a very
accurate 3-D vision on a full-size humanoid robot with a highdegree-of-freedom oculomotor system.
The core part that integrates active vision with the developed generalization approach is the estimation of the function that maps the desired Cartesian positions (query points) to
the final joint positions (DMP-attractor points) on the training

808IEEE TRANSACTIONS ON ROBOTICS, VOL. 26, NO. 5, OCTOBER 2010

Fig. 7. Image sequence showing the robot reaching over the edge of the table toward a target. The robot tracks the target with its own eyes and determines thetarget’s location, which is used to generate a new attractor point g for the DMP using GPR. This process runs in real time. The initially open-loop movement thatavoids the table thus gradually transforms into a closed-loop movement that follows the object. Note the different head and eyes postures in the images.

Fig. 8.

Image sequence illustrating the grasping process. The head and eyes actively follow the object. The robot applied a power grasp to pick up the object.

trajectories. This function is estimated by Gaussian process regression as outlined in Section II-D.
As shown in Fig. 7, the robot needs to avoid the table while
reaching toward the object. Based on the training trajectories
that avoid the table, the system generates new trajectories that
also avoid the table, but at the same time end at the desired
Cartesian positions. Without generalizing the shape parameters w, the robot’s arm would collide with the table. The first
three images of Fig. 7 and video discrete-reach-open.mov,which is attached to this paper and can be downloaded fromhttp://ieeexplore.ieee.org, show the successful open-loop reaching movements that avoid the table. Once the robot’s arm reaches
over the table, we start updating the attractor point g in real timeusing the results of active vision and GPR. No additional programming was needed to create the closed-loop behavior, justactive vision was allowed to update the attractor points. Theoretically, it is possible to update the shape parameters w as
well, but depending on the amount of data, this calculation may
require a few hundred milliseconds and is therefore unsuitablefor feedback control. Once the time exceeds τ , the influenceof parameters w starts diminishing and the robot generates its
motion based solely on the attractor dynamics and GPR, which
again results in collision-free movements (because the training
examples are collision-free). The last three images of Fig. 7
and the attached video discrete-reach-closed.mov show thatthe robot can successfully follow a moving object, thus demonstrating the power of the DMP representation when reacting to
the external feedback and the accuracy of GPR when filteringvisual information.
The resulting robot-hand path while following the object in aclosed-loop is depicted in Fig. 9. The mean value of the robothand position while following the object that moves on thetable (as calculated by forward kinematics) and the mean valueof the object position (as estimated by active vision) differ by
about [1.6, 4.2, 7.6] cm. Although the real modeling errors are

nonlinear, the error would be roughly similar if we just estimatedobject positions and converted them into joint angles by meansof inverse kinematics. Such an error would be too large forreaching and grasping. In our system, GPR successfully corrects
at least part of the modeling errors and the resulting attractor
points are accurate enough for these tasks. This is due to thefact that already during training, the system directly relates theobject positions as estimated by vision to the final robot-joint
configurations. The vision errors can therefore be taken into
account by GPR, the system only needs to be repeatable.
The realized grasping behavior is shown in Fig. 8 and in
the attached video discrete-grasping.mov. Active vision detects when the object stops moving and initiates the open-loop
reaching. The attractor points for the initial reaching trajectory
are generated by suitably displacing the estimated query points(based on the desired approach direction). After the initial reaching movement has finished, the vision system starts supplying
the current object position and the second-order attractor dynamics automatically generates a closed-loop approach motion
for grasping. The system is accurate enough to pick an object
from the hand of a person. Note that without active eye andhead joints, the robot would not be able to follow the object,
especially when it comes close to the body, and the behavior
could not be generated.

D. Simulated Ball Throwing
To demonstrate the performance of the approach for the synthesis of more dynamic tasks, we considered the task of throwing
a ball into a basket using the same 2R robot for simulation as
in reaching. Throwing depends not only on the positional partof the movement but on velocities as well. The trajectory of
the ball after the release is fully specified by the position andUDE et al.: TASK-SPECIFIC GENERALIZATION OF DISCRETE AND PERIODIC DYNAMIC MOVEMENT PRIMITIVES

809

TABLE IIIERRORS IN THE SYNTHESIZED BALL THROWS (DISTANCE FROM THE TARGETIN CENTIMETERS)

Fig. 9. The upper image shows the 3-D path of the robot’s hand (calculated
by forward kinematics) while reaching and tracking an object (red). The object
positions estimated by active vision are shown in green, while the blue curveshows the estimated object positions with corrected mean value so that it coincides with the mean of the hand motion. The lower image shows the samedata projected onto a plane, which is roughly parallel to the table. The proposed
approach is successful at correcting systematic vision errors.

velocity at the release timex = x0 + v0 t cos(α), y = y0 + v0 t sin(α) −

gt2
2

(20)

where (x0 , y0 ) is the release point, v0 is the linear velocity ofthe ball at the release time, and α is the initial angle of thethrow path. We considered the problem where the target basket
is placed in xy-plane. The understanding of the physics of thetask allows us to compare movement generalization results with
an ideal system. Local models were also studied in [53] to refine
the throwing performance but not to generalize the throws overthe complete training space.
The target positions, i e., the positions where the ball is supposed to land, were used as query points. The training-query
points were uniformly distributed within a rectangular area withcorners at (1.2, 0.1) and (5.2, 2.1) m, with different spacingsas shown in Table III. Based on (20), example movements withproper position and velocity at release times were analyticallygenerated for the training targets. To avoid accurately modelingthe physical grip and the release of the ball from the hand—which is rather complex to model but would not contributemuch to our analysis—we also attached the release time toevery training trajectory (besides the usual duration and attractor points). This was not necessary to do on a real robot (seeSection III-E). For evaluation, new query points were generated

Fig. 10. Simulation results. (Top left) Joint velocities originating from 45
training trajectories. Only the relevant parts of the throwing trajectories (upto the release time) are displayed. (Bottom left) Generalized velocities for 32
in-between query points. They end at the release times estimated by GPR. (Topright) Training target positions (stars) and the in-between query points (circles).
(Bottom right) The generalized velocities (thick black lines) for query pointq = [2.95, 1.35] and the velocities of nearby training trajectories that wereused for generalization. Our approach preserves the trends.

on a finer grid than for training (2 cm), and the corresponding
DMPs were calculated using the proposed generalization approach. The release times were also estimated by GPR. Duringexecution, the ball was detached from the hand at the estimatedrelease time. Using the position and velocity at release time, we
could calculate where the ball would land using (20).
Fig. 10 depicts the velocities of the most coarsely sampledtraining movements and their generalization. It can be seenthat the generalized velocities and release times are similar to
the training velocities and release times. This is confirmed in810IEEE TRANSACTIONS ON ROBOTICS, VOL. 26, NO. 5, OCTOBER 2010

Fig. 12. Accuracy of real throwing. The stars are the training targets while thedots show the input query points for generalization. For evaluation purposes,
the robot executed three throws at each target, and we calculated the average
throwing error from the desired target (query point).
Fig. 11. Validation-set error as a function of the bandwidth parameter c definedin (15). The curve with stars shows the average difference (18) between thesynthesized trajectories and validation trajectories, and the curve with pluses
shows the error of the synthesized ball throws. They have similar minima.

Table III, which shows that we can generate throws to any target
within the training space with high accuracy. Even when the
training targets were distributed with the spacing of 0.5 m, we
could achieve the average accuracy of 2.0 cm for the inner training area where a sufficient number of example movements areavailable. The inner training area was defined similarly as inFig. 3, with the width of the border area being equal to the spacing between the example query points. The error dropped to0.13 cm for training targets distributed with the spacing of 12.5
cm. Thus, we can say that our approach produces throws with respectable accuracy, even when the training examples are sparse,
and that the accuracy improves with more densely distributed
training data. These results also show that the accuracy is higher
when a more advanced integration technique such as Adams–
Bashforth–Moulton method [54] is used to integrate the system(24)–(26), shown later. Errors were significantly larger when
we applied the simpler Euler’s method. While the generalized
DMPs approximate the desired trajectories and their dynamics
with high precision, care must be taken when integrating thegeneralized DMPs to reproduce this accuracy, which is important in tasks such as ball throwing.
We also tested the automatic calculation of the bandwidthparameter c of (15) for the locally weighted regression using
the validation-set-error method. As can be seen in Fig. 11, thecriterion has two local minima. The accuracy of the results
was slightly better with the larger c, but similar accuracy andfaster computational times could be achieved with the smaller
value of c. Note that the minima are very similar if we evaluate
the reproduction accuracy of the synthesized motion or if wemeasure the accuracy of throwing. This fact is important because
the former does not require new experiments with the robot,
whereas the latter can only be calculated if the robot performs
new throws. This makes the second approach impractical for
real-world training.
E. Real-World Ball Throwing
To test real throwing, we recorded 20 throwing trajectories,
which were manually trained for different targets, and measured where the ball landed for each of these trajectories. We

Fig. 13. (Top left) Library of example Fig. 8 trajectories in Cartesian space.
(Top right) Frequencies of example movements in x and y direction. (Bottom
left) Trajectories of joint 1 of the 2R planar robot. (Bottom right) Trajectories
of joint 2 of the 2R planar robot.

Fig. 14. Mean error of movement generalization over the entire scope of thelibrary of example movements estimated with the step of 1 cm.

used three of the total seven degrees of freedom of PA-10 robot
to realize robotic throwing. A coaching process similar to theone described in [55] was applied to accumulate the training
data. Note that instead of using the actually executed movements as training data for generalization, it was better to use thecommanded trajectories. In this way, any discrepancies between
the commanded and actually executed trajectories can be taken
into account by the generalization process. In addition, the commanded trajectories are smoother. Unlike in simulations, here,
we did not need to consider the release times because the release times are implicitly encoded in the arm trajectories. In our
experiments, the ball was held only loosely by the gripper. It

UDE et al.: TASK-SPECIFIC GENERALIZATION OF DISCRETE AND PERIODIC DYNAMIC MOVEMENT PRIMITIVES

811

Fig. 15. Example periodic trajectories in Cartesian space. The blue curves show the training input, and the red curves show the estimated periodic trajectories.
To transform the demonstrated screen trajectories into Cartesian space, the values have to be multiplied by 0.2 m. All final estimated training trajectories are in thesame plot in the top-right image. The frequencies are in the bottom-right plot.

was detached from the hand due to the dynamics of motion.
The repeatability of throws was about 2–3 cm; hence, this isthe best accuracy achievable by learning. Fig. 12 shows that we
can obtain accuracy of about 2–10 cm for the throws within the
inner training area, which is quite good when we consider the
repeatability of the robot and that the spacing between training
throws was about 0.5 m. The accuracy needed to get the ball into
the basket was within few centimeters. As expected, the robot
performed worse for the targets at the edge of the training area,
where less training data are available, and the precision is lower(see Fig. 12 and the attached video discrete-throwing.mov, inwhich the robot misses the last three throws into a basket placedat the edge of the training area).
F. Simulated Periodic Pattern Study
Similarly to our analysis, in the case of discrete movements,
we conducted a simulation experiment to examine how well we
can generalize periodic Cartesian trajectories of a robot endeffector. Just like in the simulated reaching study, we used aplanar 2R robot for simulation.
As an example periodic motion, we chose a Cartesian figure8 trajectory with a varying amplitude and frequency. Several
reasons speak for this, namely we can analytically generate asignal with an arbitrary amplitude for comparison and the motion frequencies of the two separate dimensions in Cartesian
space are different. While the primary aim of the generalization
algorithm is to produce DMPs for movements that cannot be
attained by simple modulation, for the purpose of evaluation,
this section uses an example that can be easily attained by modulating the amplitude-control parameter of the periodic DMP.
The next section gives an example where generalization could
not be achieved by simple modulation.
A set of example figure-8s in Cartesian space and the trajectories in joint space is depicted in Fig. 13. Cartesian-space
trajectories are distributed evenly from the smallest trajectoryat x = 0.3 + 0.2 cos(Ωt) m and y = 0.4 + 0.1 sin(2Ωt) m to

the largest at y = 0.4 + 0.55 sin(2Ωt) m with a step of 5 cm (xdoes not change). Ten sampled joint-space movements are used
to form the example library.
Using the amplitude in the y dimension as a query-point parameter, we generalized the training-joint trajectories, i.e., wecalculated the parameters of a periodic DMP as described in Section II-C and D. Different frequencies in each dimension of theexample demonstration movements (see Fig. 13) show that we
can successfully estimate the periodic movements onto a single
period. Once the new movement is generated, frequencies caneasily be modulated with the periodic DMP frequency-control
parameter, just like it is done in a standard DMP approach(see [44] for details).
Fig. 14 shows the mean error of movement generalization overthe entire scope of the example library. As before, the distanceof the generalized trajectory from the analytically computed
trajectory increases at the edge of the database as there, example
movements are available only on one side of the query point.
The accuracy of generalization is slightly decreasing, which is
the result of the increasing amplitude.
G. Real-World Drumming Experiment
In the final experiment, we applied the proposed approachto realize drumming on humanoid robot CB-i. Of the total 39
degrees of freedom, we used all seven of the right arm forthe execution of the movement. To demonstrate generalization, we trained the robot to perform one-handed drummingon a drum and a cymbal, where the cymbal was mounted
at different heights. Such drum placements—although withmore drums—are common for drummers. Fig. 15 and the attached video periodic-examples.mov show a set of example
periodic trajectories in Cartesian space, which were recordedwith a sampling rate of 100 Hz. The trajectories were demonstrated with a mouse on a screen, scaled to Cartesian space,and mapped onto the robot’s joint angles via inverse kinematics in real time. The trainer was modifying his motion

812

Fig. 16.

IEEE TRANSACTIONS ON ROBOTICS, VOL. 26, NO. 5, OCTOBER 2010

Image sequence showing the drumming at two different heights. The cymbal is higher in the top row.

points above and below the current query point. These trajectories were used for generalization. Due to the limited support of
the weighting kernel from (13), the generalized trajectories are
generated only from the four plotted example trajectories. As the
trajectories were generated by human demonstrations, they arequite different among themselves. Nevertheless, the algorithm
was able to generate trajectories similar to the example trajectories. If the demonstrations were more uniform, the generalized
trajectories would be as well.

IV. FINAL DISCUSSION AND CONCLUSION

Fig. 17. The result of generalization for four randomly chosen query points.
The figure shows the generalized trajectories in blue. The training trajectories
associated with query points just above (red) and below (green) the current query
point are also shown. Among the ten example trajectories, these trajectories
influenced the generalization process most.

based on the real-time visual feedback coming from the robot.
The drum was placed at −0.162 m. For training, we mounted
the cymbal at ten different heights (−0.177, −0.149, −0.091,−0.068, −0.049 −0.020, 0.017, 0.054, 0.079, and 0.117 m).
Thus, the maximum height difference was about 30 cm.
Fig. 16 and the attached videos periodic-generalized1.mov
and periodic-generalized2.mov show the ability of the system to generalize the trained periodic trajectories. The resulting
movements cannot be attained trough simple modulation of aperiodic DMP because standard modulation techniques do notallow the robot to modify the movement on only one side. The
proposed approach can automatically generate the appropriate
movements from the collected data. The query point for generalization is the height difference between the drum and thecymbal. While the height difference could be estimated by vision, we simplified this experiment by measuring the differencewith a ruler and providing it to the algorithm.
Fig. 17 presents the results of generalization for four randomly selected query points. The figure shows the generalized
periodic trajectories and the trajectories associated with query

Our experiments show that the proposed approach is able togeneralize the training data to new situations within the sampledtraining space. The computational time required for generalization compared to the computational time for standard one-shotlearning of DMPs increases only linearly with a number of training trajectories that are taken into account by local weighting.
The relatively low computational complexity of the method enabled us to compute new control policies directly from the sampled data online. In this way, we avoid the pitfalls associated withthe projection of the training data into lower dimensional parameter spaces (also called latent spaces), which can lead to oversmoothing and, therefore, losing important details of the task.
In contrast to reinforcement learning (RL), where a robot actively explores the solution space to adapt its movements to newsituations, we focused on generalization to new situations fromthe available data. The two methodologies are complementary;
RL techniques can be used to provide initial samples for our
trajectory libraries. On the other hand, our method can be used
to provide a good initial approximation for RL (thus speeding up
the convergence) if higher accuracy than what can be achievedby the proposed approach is required.
Our method is appropriate if the example trajectories
smoothly transition as a function of query points, as defined in(2). Otherwise, nearby data do not provide information about the
movement associated with the new query point q. A more complex relationship might require a larger number of training trajectories. It is, however, intuitive that this relationship is smoothin real-world problems because it is unlikely that a totally

UDE et al.: TASK-SPECIFIC GENERALIZATION OF DISCRETE AND PERIODIC DYNAMIC MOVEMENT PRIMITIVES

different strategy would be used to solve a task in different
but similar situations. Instead, the movement trajectories areusually adjusted when solving the task at nearby query points.
Our experiments show that even though the available training
sets were not large in real-world test, we could successfully
generalize and achieve high performance in several tasks thatare relevant for robotics.
The proposed generalization process synthesizes DMPs in astandard form, and thus, we can still exploit all of the advantages of standard DMPs at execution time, like, for example,
modulation of the time evolution of the phase by modifying(24) [20], avoiding the joint limits by modifying (22) [44], etc.
We have also demonstrated in Section III-C that the developed
system can be used in an active feedback loop by modifying theattractor point of a DMP online, with GPR correcting some of
the modeling errors causing uncertainties in the 3-D vision data.
There exist tasks in which example movements do not transition smoothly as a function of query points. Consider, forexample, reaching movements that need to avoid an obstacle
before arriving at the final destination. If there are two sets ofexample movements, each avoiding the obstacle from a different side, then example movements that avoid the obstacle fromdifferent sides should not be used for generalization simultaneously. The proposed approach could still be used, but it would
need to be supplemented by a suitable clustering procedure that
determines sets of trajectories suitable for generalization. Issues
like clustering of example movements and the identification of
motor primitives are important for cognitive robots that should
gradually increase their competence through seamless learning.
While we did not consider such issues in this paper, our approach provides a methodology to generalize the appropriately
selected training trajectories in a natural way, thus providing animportant building block for a cognitive robot.
APPENDIX A
CONTROL POLICIES AS DYNAMIC SYSTEMS
Here, we briefly explain the theoretical fundamentals of the
motor representation developed by Ijspeert et al. [28], [29].
These authors proposed to describe a control policy by a setof nonlinear differential equations with a well-defined attractor
dynamics. Here the most current formulation as outlined in[20] is used. For a single degree of freedom denoted by y,which can either be one of the internal joint angles or one of
the external task-space coordinates, the following system of
linear differential equations with constant coefficients has been
proposed as a basis for motion specification:
τ ż = αz (βz (g − y) − z)(21)τ ẏ = z.

(22)

Provided that the parameters αz , βz , and τ > 0 are selected
appropriately, e.g., αz = 4βz , this system has a unique attractor
point at y = g, z = 0.
Differential equations (21)–(22) ensure that y converges tog and can therefore be used to realize discrete point-to-point
movements. To increase a rather limited set of trajectories that813

can be encoded by (21) and (22) and thus enable the approximation of general point-to-point movements, (21) needs to bemodified. In the case of discrete movements, one can add alinear combination of radial-basis functions to (21) [20]3N
i=1 wi Ψi (x)x, Ψi (x) = exp −hi (x − ci )2
f (x) = 
N
Ψ
(x)
i
i=1
(23)
where ci are the centers of radial basis function distributed along
the trajectory, and hi > 0. A phase variable x is used in (23)instead of time to make the dependency of f on time more
implicit. Its dynamics can be defined byτ ẋ = −αx x(24)with initial value x(0) = 1. A solution to (24) is given by
exp (−αx t/τ ), thus x tends to 0 as time increases. As shownin [20], the appealing property of using the phase variable x
instead of explicit time is that by appropriately modifying (24),
we can, for example, stop the evolution of time to account for
perturbations during trajectory execution. This results in thefollowing system of differential equations:τ ż = αz (βz (g − y) − z) + f (x)(25)τ ẏ = z(26)which can be used to approximate discrete movements of various
shapes. Since x tends to zero, the influence of the nonlinear termf (x) decreases with time and system (25) and (26) converges
to [0, g]T , just like the system (21) and (22). The control policy
specified by variable y defines what is called DMP.
In the case of periodic movements, the following linear combination of periodic functions can be used to change the dynamics of the basic second-order system [20]:N
i=1 wi Γi (φ)
r, Γi (φ) = exp (hi (cos (φ − ci ) − 1))
f (φ) = Ni=1 Γi (φ)(27)
where r is the amplitude of the oscillator, and hi > 0. Writingτ = 1/Ω, (25) and (26) are replaced byż = Ω (αz (βz (g − y) − z) + f (φ))(28)

ẏ = Ωz.

(29)

The phase variable φ has been introduced in this case to avoid
the explicit dependency on time. The phase is assumed to movewith constant speedφ̇ = Ω(30)

where Ω is the frequency of oscillation.
DMPs have been designed to provide a representation thatenables accurate encoding of the desired trajectories and at thesame time permit the modulation of different properties of theencoded trajectories. In this context, the shape parameters wi
3f

N

defined in [20] is scaled by g − y 0 , i.e., f (x) = [

N
i= 1

w i Ψ i (x)/

Ψ i (x)]x(g − y 0 ), y 0 = y(0). Thus, when the attractor point gi= 1changes, the encoded movement gets scaled. We omit this scaling factor because we are not interested in automatic scaling, which is achieved differently
in our approach. If g is kept constant, the scaling factor has no effect.

814IEEE TRANSACTIONS ON ROBOTICS, VOL. 26, NO. 5, OCTOBER 2010

are determined so that the robot can accurately follow the desired trajectory by integrating (24)–(26) or (28)–(30). The otherparameters are used for modulation. For example, by changing
g, we can adjust the final destination of a discrete movement,
while τ can be adjusted to modulate its velocity. Similarly, Ω
can be changed to modulate the frequency of oscillation in caseof periodic movements.
APPENDIX B
ADAPTIVE FREQUENCY OSCILLATORSThe frequency of oscillation is not directly observable in thedata. To estimate the frequency, Righetti et al. [43] suggested toreplace the constant-speed assumption (30) by a systemφ̇i = Ωi − Ke (t) sin (φi )(31)Ω̇i = −Ke (t) sin (φi )(32)α̇i = η cos (φi ) e (t)

(33)

L

where e (t) = yd (t) − ŷ(t), and ŷ(t) = i=1 αi cos (φi ). Notethat if e (t) = 0, the system (31)–(33) becomes equivalent to(30). It has been shown that by integrating this system, thefrequencies Ωi contained in the observed motion trajectory canbe estimated. The most significant frequency is selected as thebase or fundamental frequency Ω for the DMP (see [44] for the
integration of adaptive frequency oscillators with DMPs).
APPENDIX CGAUSSIAN PROCESS REGRESSION
A Gaussian process is defined asg(q) ∼ GP (m(q), k(q, q ))(34)

where m(q) = E(g(q)) is the mean function, and k(q, q ) =E((g(q) − m(q))(g(q ) − m(q ))) is the covariance functionof the process. Let us assume that we have a set of noisy
observations {(qk , yk )|k = 1, . . . , M }, yk = g(qk ) + ,  ∼N (0, σn2 ). Subtracting the mean from the training data, we canfurther assume that m(q) = 0. If we are given a set of query
points g(q∗ ), then the joint distribution of all outputs is given as
 

 
K(Q, Q) + σn2 I K(Q, Q∗ )y
∼ N 0,(35)
K(Q∗ , Q∗ )y∗K(Q∗ , Q)
where Q, Q∗ , y, y∗ , respectively, combine all inputs and outputs, and K(·, ·) are the associated joint covariance matrices
calculated according to (34). It can be shown [42] that the expected value ȳ∗ is given byȳ∗ = E(y∗ |Q, y, Q∗ ) = K(Q∗ , Q)[K(Q, Q) + σn2 I]−1 y(36)with the following estimate for the covariance of the prediction:cov(g∗ ) = K(Q∗ , Q∗ )− K(Q∗ , Q)[K(Q, Q) + σn2 I]−1 K(Q, Q∗ ).
One commonly used covariance function isn
1 (qi − qi )22
exp −k(q, q ) = σf
2
li2
i=1

(37)

which results in a Bayesian regression model with an infinite
number of basis functions. Here, n denotes the dimension of thequery-point space. See [42] for more details.
REFERENCES[1] S. Schaal, “Is imitation learning the route to humanoid robots?” Trends
Cogn. Sci., vol. 3, no. 6, pp. 233–242, 1999.
[2] R. Dillmann, “Teaching and learning of robot tasks via observation ofhuman performance,” Robot. Auton. Syst., vol. 47, no. 2–3, pp. 109–116,
2004.
[3] M. Riley, A. Ude, and C. G. Atkeson, “Methods for motion generation andinteraction with a humanoid robot: Case studies of dancing and catching,”
in Proc. 2000 Workshop Interactive Robot. Entertainment, Pittsburgh, PA,
2000, pp. 35–42.
[4] N. S. Pollard, J. K. Hodgins, M. Riley, and C. G. Atkeson, “Adapting
human motion for the control of a humanoid robot,” in Proc. IEEE Int.
Conf. Robot. Autom., Washington, DC, 2002, pp. 1390–1397.
[5] A. Ude, C. G. Atkeson, and M. Riley, “Programming full-body movements
for humanoid robots by observation,” Robot. Auton. Syst., vol. 47, no. 2–3,
pp. 93–108, 2004.
[6] M. Ruchanurucks, S. Nakaoka, S. Kudoh, and K. Ikeuchi, “Humanoid
robot motion generation with sequential physical constraints,” in Proc.
IEEE Int. Conf. Robot. Autom., Orlando, FL, 2006, pp. 2649–2654.
[7] T. B. Moeslund, A. Hilton, and V. Krüger, “A survey of advances invision-based human motion capture and analysis,” Comput. Vis. Image
Understanding, vol. 104, no. 2, pp. 90–126, 2006.
[8] M. Hersch, F. Guenter, S. Calinon, and A. Billard, “Dynamical systemmodulation for robot learning via kinesthetic demonstrations,” IEEE
Trans. Robot., vol. 24, no. 6, pp. 1463–1467, Dec. 2008.
[9] J. Peters and S. Schaal, “Reinforcement learning of motor skills with
policy gradients,” Neural Netw., vol. 21, pp. 682–697, 2008.
[10] M. Mistry, P. Mohajerian, and S. Schaal, “An exoskeleton robot for humanarm movement study,” in Proc. IEEE/RSJ Int. Conf. Intell. Robots Syst.,
Edmonton, Canada, 2005, pp. 4071–4076.
[11] A. Albu-Schäffer, C. Ott, and G. Hirzinger, “A unified passivity based
control framework for position, torque and impedance control of flexible
joint robots,” Int. J. Robot. Res., vol. 26, no. 1, pp. 23–39, 2007.
[12] S.-H. Hyon, J. G. Hale, and G. Cheng, “Full-body compliant human–
humanoid interaction: Balancing in the presence of unknown external
forces,” IEEE Trans. Robot., vol. 23, no. 5, pp. 884–898, Oct. 2007.
[13] Y. Kuniyoshi, M. Inaba, and H. Inoue, “Learning by watching: Extracting
reusable task knowledge from visual observation of human performance,”
IEEE Trans. Robot. Automat., vol. 10, no. 6, pp. 799–822, Dec. 1994.
[14] S.-B. Kang and K. Ikeuchi, “Toward automatic robot instruction fromperception – Temporal segmentation of tasks from human hand motion,”
IEEE Trans. Robot. Automat., vol. 11, no. 5, pp. 670–681, Oct. 1995.
[15] M. Kaiser and R. Dillmann, “Building elementary robot skills from humandemonstration,” in Proc. IEEE Int. Conf. Robot. Autom., Minneapolis,
MN, 1996, pp. 2700–2705.
[16] H. Miyamoto, S. Schaal, F. Gandolfo, H. Gomi, Y. Koike, R. Osu,
E. Nakano, Y. Wada, and M. Kawato, “A kendama learning robot based on
bi-directional theory,” Neural Netw., vol. 9, no. 8, pp. 1281–1302, 1996.
[17] F. Guenter, M. Hersch, S. Calinon, and A. Billard, “Reinforcement learning for imitating constrained reaching movements,” Adv. Robot., vol. 21,
pp. 1521–1544, 2007.
[18] J. Peters, J. Kober, K. Muelling, D. Nguyen-Tuong, and O. Kroemer,“Towards motor skill learning for robotics,” presented at the Int. Symp.
Robot. Res., Lucerne, Switzerland, 2009.
[19] L. Sciavicco and B. Siciliano, Modeling and Control of Robot Manipulators. London, U.K.: Springer-Verlag, 2000.
[20] S. Schaal, P. Mohajerian, and A. Ijspeert, “Dynamics systems vs. optimal
control—A unifying view,” Progr. Brain Res., vol. 165, no. 6, pp. 425–
445, 2007.
[21] T. Asfour, F. Gyarfas, P. Azad, and R. Dilmann, “Imitation learning of dualarm manipulation tasks in humanoid robots,” Int. J. Humanoid Robot.,vol. 5, no. 2, pp. 183–202, 2008.
[22] A. Billard, S. Calinon, and F. Guenter, “Discriminative and adaptive imitation in uni-manual and bi-manual tasks,” Robot. Auton. Syst., vol. 54,
pp. 370–384, 2006.
[23] T. Inamura, I. Toshima, H. Tanie, and Y. Nakamura, “Embodied symbol
emergence based on mimesis theory,” Int. J. Robot. Res., vol. 23, no. 4–5,
pp. 363–377, 2004.
[24] D. Kulić, W. Takano, and Y. Nakamura, “Incremental learning, clusteringand hierarchy formation of whole body motion patterns using adaptive

UDE et al.: TASK-SPECIFIC GENERALIZATION OF DISCRETE AND PERIODIC DYNAMIC MOVEMENT PRIMITIVES

hidden Markov chains,” Int. J. Robot. Res., vol. 27, no. 7, pp. 761–784,
2008.
[25] D. Kulić, W. Takano, and Y. Nakamura, “Online segmentation and clustering from continuous observation of whole body motions,” IEEE Trans.
Robot., vol. 25, no. 5, pp. 1158–1166, Oct. 2009.
[26] S. Calinon, F. Guenter, and A. Billard, “On learning, representing, andgeneralizing a task in a humanoid robot,” IEEE Trans. Syst. Man Cybern.
B, Cybern., vol. 32, no. 2, pp. 286–298, Apr. 2007.
[27] K. Yamane, Y. Yamaguchi, and Y. Nakamura, “Human motion databasewith a binary tree and node transition graphs,” presented at the Robot.,
Sci. Syst. Conf., Seattle, Washington, 2009.
[28] A. J. Ijspeert, J. Nakanishi, and S. Schaal, “Movement imitation withnonlinear dynamical systems in humanoid robots,” in Proc. IEEE Int.
Conf. Robot. Autom., Washington, DC, 2002, pp. 1398–1403.
[29] A. J. Ijspeert, J. Nakanishi, and S. Schaal, “Learning rhythmic movements
by demonstration using nonlinear oscillators,” in Proc. IEEE/RSJ Int.
Conf. Intell. Robots Syst., Lausanne, Switzerland, 2002, pp. 958–963.
[30] T. Poggio and E. Bizzi, “Generalization in vision and motor control,”Nature, vol. 431, pp. 768–774, 2004.
[31] C. G. Atkeson, J. Hale, F. Pollick, M. Riley, S. Kotosaka, S. Schaal,
T. Shibata, G. Tevatia, A. Ude, S. Vijayakumar, and M. Kawato, “Usinghumanoid robots to study human behavior,” IEEE Intell. Syst., vol. 15,
no. 4, pp. 46–56, Jul./Aug. 2000.
[32] M. Raibert, “A model for sensorimotor control and learning,” Biol. Cybern., vol. 29, pp. 29–36, 1978.
[33] C. G. Atkeson, A. W. Moore, and S. Schaal, “Locally weighted learningfor control,” Artif. Intell. Rev., vol. 11, pp. 75–113, 1997.
[34] P. Pastor, H. Hoffmann, T. Asfour, and S. Schaal, “Learning and generalization of motor skills by learning from demonstration,” in Proc. IEEE
Int. Conf. Robot. Autom., Kobe, Japan, 2009, pp. 763–769.
[35] D.-H. Park, H. Hoffmann, P. Pastor, and S. Schaal, “Movement reproduction and obstacle avoidance with dynamic movement primitives andpotential fields,” in Proc. IEEE-RAS Int. Conf. Humanoid Robots, Daejeon, Korea, 2008, pp. 91–98.
[36] J. Kober, B. Mohler, and J. Peters, “Learning perceptual coupling for
motor primitives,” in Proc. IEEE/RSJ Int. Conf. Intell. Robots Syst., Nice,
France, 2008, pp. 834–839.
[37] F. A. Mussa-Ivaldi and E. Bizzi, “Motor learning through the combinationof primitives,” Phil. Tran. Royal Soc. London B, vol. 355, pp. 1755–1769,
2000.
[38] D. Wolpert and M. Kawato, “Multiple paired forward and inverse models
for motor control,” Neural Netw., vol. 17, pp. 1317–1329, 1998.
[39] M. Haruno, D. M. Wolpert, and M. Kawato, “MOSAIC model for sensorimotor learning and control,” Neural Comput., vol. 13, pp. 2201–2220,
2001.
[40] C. G. Atkeson, A. W. Moore, and S. Schaal, “Locally weighted learning,”
Artif. Intell. Rev., vol. 11, pp. 11–73, 1997.
[41] D. Nguyen-Tuong, M. Seeger, and J. Peters, “Model learning with local
Gaussian process regression,” Adv. Robot., vol. 23, pp. 2015–2034, 2009.
[42] C. E. Rasmussen and C. Williams, Gaussian Processes for Machine Learning. Cambridge, MA: MIT Press, 2006.
[43] L. Righetti, J. Buchli, and A. J. Ijspeert, “Dynamic Hebbian learning in
adaptive frequency oscillators,” Physica D, vol. 216, pp. 269–281, 2006.
[44] A. Gams, A. J. Ijspeert, S. Schaal, and J. Lenarčič, “On-line learning andmodulation of periodic movements with nonlinear dynamical systems,”
Auton. Robots, vol. 27, no. 1, pp. 3–23, 2009.
[45] E. Oztop, J. Babič, J. G. Hale, G. Cheng, and M. Kawato, “From biologically realistic imitation to robot teaching via human motor learning,”
in Neural Information Processing: 14th International Conference (Lecture Notes in Computer Science Series). Berlin/Heidelberg, Germany:
Springer-Verlag, 2008, pp. 214–222.
[46] S. Schaal and C. G. Atkeson, “Constructive incremental learning fromonly local information,” Neural Comput., vol. 10, no. 8, pp. 2047–2084,
1998.
[47] D. W. Scott, Multivariate Density Estimation: Theory, Practice, and Visualization. New York: Wiley, 1992.
[48] G. Cheng, S.-H. Hyon, J. Morimoto, A. Ude, J. G. Hale, G. Colvin,
W. Scroggin, and S. C. Jacobsen, “CB: A humanoid research platform for
exploring neuroscience,” Adv. Robot., vol. 21, no. 10, pp. 1097–1114,
2007.
[49] T. Flash and N. Hogan, “The coordination of arm movements: An experimentally confirmed mathematical model,” J. Neurosci., vol. 5, no. 7,
pp. 1688–1703, 1985.
[50] A. D’Souza, S. Vijayakumar, and S. Schaal, “Learning inverse kinematics,” in Proc. IEEE/RSJ Int. Conf. Intell. Robots Syst., Maui, HI, 2001,
pp. 298–303.

815[51] J. Peters and S. Schaal, “Learning to control in operational space,” Int. J.
Robot. Res., vol. 27, no. 2, pp. 197–212, 2008.
[52] A. Ude and E. Oztop, “Active 3-D vision on a humanoid head,” presented
at the 14th Int. Conf. Adv. Robot., Munich, Germany, 2009.
[53] E. W. Aboaf, C. G. Atkeson, and D. J. Reinkensmeyer, “Task-level robot
learning,” in Proc. IEEE Int. Conf. Robot. Autom., Philadelphia, PA, 1988,
pp. 1309–1310.
[54] J. D. Hoffman, Numerical Methods for Engineers and Scientists, 2nd ed.
Boca Raton, FL: CRC, 2001.
[55] M. Riley, A. Ude, C. G. Atkeson, and G. Cheng, “Coaching: An approachto efficiently and intuitively create humanoid robot behaviors,” in Proc.
IEEE-RAS Int. Conf. Humanoid Robots, Genoa, Italy, 2006, pp. 567–574.
Aleš Ude received the Diploma degree in appliedmathematics from the University of Ljubljana, Ljubljana, Slovenia, in 1990, and the Ph.D. degree fromthe Faculty of Informatics, University of Karlsruhe,
Karlsruhe, Germany, in 1995.
He is currently a Senior Research Associatewith the Department of Automatics, Biocybernetics, and Robotics, Jožef Stefan Institute, Ljubljana.
He is also with Computational Neuroscience Laboratories, Advanced Telecommunications ResearchInstitute International, Kyoto, Japan. His research interests include imitation learning, perception of human activity, humanoid robot
vision, and humanoid cognition.
Dr. Ude is a recipient of the Science and Technology Agency fellowship for
postdoctoral studies with the Exploratory Research for Advanced Technology(ERATO) Kawato Dynamic Brain Project, Japan.

Andrej Gams received the Ph.D. degree inrobotics from the University of Ljubljana, Ljubljana,
Slovenia, in 2009.
He is currently a Postdoctoral Assistant with
the Department of Automatics, Biocybernetics, andRobotics, Jožef Stefan Institute, Ljubljana. He performed a part of his doctoral research with the Biologically Inspired Robotics Group, Ecole Polytechnique Federale de Lausanne, Lausanne, Switzerland.
His research interests include human arm motion andthe imitation of rhythmic tasks.

Tamim Asfour received the Diploma degree in electrical engineering in 1994 and the Ph.D. degree incomputer science in 2003, both from the Universityof Karlsruhe, Karlsruhe, Germany.
He is currently a Senior Researcher with the
Karlsruhe Institute of Technology, where he is alsoa Leader of the Humanoid Research Group, Institutefor Anthropomatics. His research interest includes
humanoid robotics.
Dr. Asfour received the Forschungszentrum Informatik prize for his outstanding Ph.D. dissertation
on sensorimotor control in humanoid robotics and the development of the humanoid robot ARMAR in 2003.

Jun Morimoto received the Ph.D. degree in information science from Nara Institute of Science andTechnology, Nara, Japan, in 2001.
From 2001 to 2002, he was a Postdoctoral Fellow
with the Robotics Institute, Carnegie Mellon University, Pittsburgh, PA. Since 2002, he has been with
the Advanced Telecommunications Research Institute International, Kyoto, Japan, where he was a Researcher with the Computational Brain Project, International Cooperative Research Project, Japan Scienceand Technology Agency, from 2004 to 2009, and iscurrently the Head of the Department of Brain Robot Interface, Computational
Neuroscience Laboratories.