Learning Stable Dynamical Systems using

Contraction Theory
Caroline Blocher2 , Matteo Saveriano1 and Dongheui Lee1
1
Human-centered Assistive Robotics, Technical University of Munich, Munich, Germany
(Tel: +49-89-289-26900; E-mail: {[Link],dhlee}@[Link])
2
Department of Civil and Environmental Engineering, Imperial College, London, UK (E-mail: c.blocher16@[Link])

Abstract—This paper discusses the learning of robot point- Gaussian Mixture Regression (GMR)
to-point motions via non-linear dynamical systems and Gaus- Spurious
Demonstrations
sian Mixture Regression (GMR). The novelty of the proposed attractors
approach consists in guaranteeing the stability of a learned GMR
dynamical system via Contraction theory. A contraction analysis
is performed to derive sufficient conditions for the global stability
of a dynamical system represented by GMR. The results of
this analysis are exploited to automatically compute a control
input which stabilizes the learned system on-line. Simple and Contracting
Control Stable motion
effective solutions are proposed to generate motion trajectories (Sec. III - Sec. IV-A)
close to the demonstrated ones, without affecting the stability x +
of the overall system. The proposed approach is evaluated on a Activation
public benchmark of point-to-point motions and compared with Function
(Sec. IV-B)
state-of-the-art algorithms based on Lyapunov stability theory. Proposed Approach
Keywords—Learning contracting systems. Stable discrete
movements. Learning from demonstration. Contraction theory. Fig. 1. System overview. A dynamical system represented by GMR generates
spurious attractors (top). The dynamical systems reaches the desired target by
applying an automatically computed control input (bottom).
1. I NTRODUCTION
between demonstrations and quadratic stability constraints,
Programming by Demonstration (PbD) [1] is a useful tool to which prevents SEDS to accurately learn the motion [9].
rapidly increase robot skills. Instead of explicitly programming The accurate motion reproduction is considered in [8] and
a robot, the user teaches the robot how to perform a certain in the SEDSII approach [9]. A common idea between [8]
task by demonstrating the correct task execution [2], [3]. and [9] is to learn a Lyapunov function, which minimizes the
Learned skills are represented in a compact form which contradictions between the stability constraints and the training
reduces memory requirements [4]. data, from demonstrations. In [8], the DS is represented by an
A recent trend in PbD research [5]–[11] suggests to repre- extreme learning machine (ELM) [15]. ELM parameters are
sent robotic skills via stable dynamical systems (DS). Stable learned by solving a constrained optimization problem, where
DS are well suited to represent point-to-point motions, since stability constraints are derived from the learned Lyapunov
they are guaranteed to converge to a specified target. Moreover, function. In SEDSII the Lyapunov function is used to compute
dynamical systems are robust to changes in the initial/target a stabilizing control input. The main advantage of SEDSII
location, and can be used in cluttered environments to generate is its applicability to any regression technique, while [8] is
collision-free paths [12]–[14]. Dynamic Movement Primitives limited to ELM. Alternative approaches [10], [11] propose to
(DMP) [5] are one of the first example of DS learned by learn a diffeomorphic transformation that maps the trajectories
demonstrations. DMP superimpose a linear DS and a non- of the DS into a space where demonstrations are accurately
linear force term learned from a single demonstration. reproduced. The approach in [11] is significantly faster than
The generation of stable motions from non-linear DS, [10], but it works only for linear DS.
represented by GMR, is considered in [6], where asymptotic Contraction theory [16] is effectively applied in [17] to
convergence to the target is guaranteed for trajectories that stabilize a neural network DS. Authors perform the contraction
remain inside the demonstration area. The global convergence analysis and use the resulting constraints to learn the neural
of a DS, represented by GMR, is guaranteed by the stable network paramters. The parameter learning is formulated as a
estimator of dynamical systems (SEDS) [7]. In SEDS, GMM constrained optimization problem.
parameters are learned by solving a constrained optimiza- The contribution of this paper is to investigate the applica-
tion problem, where stability constraints are derived from a bility of Contraction theory to stable motion generation via DS
quadratic Lyapunov function. The main advantage of SEDS and GMR. To this end, we perform the contraction analysis
is that the learned DS is proven to be globally stable [7]. of a DS represented by GMR and derive sufficient stability
For complex motions, however, contradictions may occur conditions. These conditions are leveraged to automatically
 TABLE 1 (see Definition 3 in Tab. 1). Notably, contraction of a DS can
D EFINITIONS be analyzed by using the matrix measure µ(F ) of the Jacobian.
Definition
In order to show that (1) is contracting, we have to compute
" 1:#GMR parameters
mkx
"
Σkx Σkxẋ
#
the Jacobian matrix. This computation is not trivial, due to the
k
m = k
,Σ = , Ak = Σkẋx (Σkx )−1 , term ∂hk (x)/∂x. Hence, we exploit Partial Contraction theory
mkẋ Σkẋx Σkẋ
π k N (x|mk k
x ,Σx )
[18] to simplify the problem. We define an auxiliary DS:
bk = mkẋ − Ak mkx , hk (x) = PK i i i
i=1 π N (x|mx ,Σx ) K
Definition 2: Let ||A||i be an induced matrix norm of A on Rm×m . X
ẏ = hk (x)(Ak y + bk ) (2)
The corresponding matrix measure is defined by
||I+A||i −1 k=1
µ(A) = lim→0+ 
. Useful properties of µ(A) are:
1) µ(αA) = α µ(A), ∀A ∈ Rm×m , α ≥ 0 where y is the auxiliary state variable. Proving that the
2) µ(A + B) ≤ µ(A) + µ(B), ∀A, B ∈ Rm×m auxiliary system has an equilibrium point x? and that it is
Definition 3: Jacobian F (x, t) = ∂f /∂x is uniformly negative definite globally contracting allows us to conclude that the original
if ∃β > 0|∀x ∈ C ⊆ Rm , ∀t ≥ 0, 12 F (x, t)T + F (x, t) ≤ −βI < 0


x? [18]. The Jacobian
DS globally exponentially converges to P
K
of the auxiliary system is simply Ja = k=1 hk (x)Ak .
compute a stabilizing control input given the GMR parameters.
The control action is smoothly activated or deactivated using Theorem 1. The dynamical system in (1) globally exponen-
an activation function (see Fig. 1). In this way, we improve tially converges to x? if for all k = 1, . . . , K:
the reproduction accuracy without affecting the stability of the • bk = −Ak x? (3)
controlled DS. Compared to [17], we perform the contraction k
• ∃µ(·) such that µ(A ) < 0 (4)
analysis of a GMR system, and use the resulting constraints
to compute, at run time, a stabilizing control input. Our Proof. Condition (3) guarantees that the auxiliary system (2)
approach shares with SEDSII the idea of stabilizing the DS has an equilibrium at x? . Indeed, by substituting (3) into
PK point
at run time. In contrast to SEDSII, the proposed control (2), we have k=1 h (x)Ak (y−x? ) = 0 if y = x? . Condition
k
law does not require the learning of a suitable Lyapunov (4) guarantees that the auxiliary system is globally contracting.
function. The performance of the proposed approach, called Indeed, considering (4), Definition 2,P and that 0 ≤ hx ≤ 1
Contracting GMR (C-GMR), are evaluated on a public dataset K
from Definition 1, we have µ(Ja ) ≤ k=1 hk (x)µ(Ak ) < 0
and compared with Lyapunov based approaches in [7], [9]. for all y ∈ Rm . As the auxiliary system is autonomous and
globally contracting, it has a unique equilibrium at x? [16].
2. P ROPOSED A PPROACH
Partial Contraction theory allows us to conclude that the DS
2.1. Problem Definition (1) globally exponentially converges to x? .
In this section, we recap on DS represented by GMR from
[7]. Let us consider N demonstrations D = {xt,i , ẋt,i }T,N 2.3. Global Exponential Stabilizer
t=1,i=1
of a point-to-point motion, where xt,i , ẋt,i ∈ Rm are re- Stability conditions in Theorem 1 allow the usage of differ-
spectively the position and the velocity. We assume that the ent matrix measures, Lyapunov theory requires the usage of
demonstrations are drawn from a non-linear and smooth DS the Euclidean norm. We exploit this advantage of Contraction
ẋ = f (x, θ), where θ ∈ Rd is a vector of parameters theory to automatically compute a stabilizing controller.
depending on the regression technique. GMR is used in this The DS in (1) can be stabilized with the control input:
work to model f (x, θ). Therefore, the non-linear function K   XK  
f (x, θ) is parametrized by θ = {π k , mk , Σk }K
X
k=1 , i.e. by the ẋ = hk (x) Ak x + bk + hk (x) U k x − Ak x? (5)
priors π k , the means mk and the covariance matrices Σk (see k=1
| {z }
k=1
| {z }
Tab. 1). f (x, θ) can be written as [7]: System Control

K
X To prove that (5) globally exponentially
PK converges to x? , we
ẋ = f (x, θ) = hk (x)(Ak x + bk ) (1) consider the auxiliary DS ẏ = k=1 h
k
(x)(A y + bk ) +
k
PK k k k ?
PK k k ?
k=1 h (x)(U y − A x ). The term −
k=1
k=1 h (x)A x is
k k k
where A , b and h (x) are defined in Tab. 1. As shown in introduced to satisfyPcondition (3). The auxiliary DS can be
K k k k


[7], the DS in (1) can generate spurious equilibria, which limit re-written as ẏ = k=1 h (x) A + U (y − x? ), which
the applicability of GMR in point-to-point motion learning. underlines that the control matrices U k have to guarantee
µ(Ak + U k ) < 0, k = 1, .., K (condition (4)).
2.2. Contraction Analysis To compute the control matrices U k , we propose an auto-
Contraction theory [16] is a novel approach to analyze the matic procedure based on the measure µ1 (·) associated to the
stability of non-linear DS. Contraction theory is based on l1 -norm. The procedure is summarized in Algorithm 1. Given
the idea that moving along a trajectory of the contracting a square matrix A ∈ Rm×m , Algorithm 1 finds a diagonal
system, the pointwise distance to its neighboring trajectories matrix U ∈ Rm×m such that C = A+U has negative diagonal
exponentially decreases. In particular, a DS is contracting if elementsPcii < 0, i = 1, .., m and is diagonally dominant
m
the generalized Jacobian F (x, t) is uniformly negative definite |cii | > j=1 |cij |, i = 1, .., m. In other words, Algorithm
Algorithm 1 Find a matrix U such that µ1 (A + U ) < 0 Target
Require: A ∈ Rm×m , p > 1, m // dimension of state space Mean
Covariance
U ← Matrix of zeros in Rm×m
for d =P1 to m do

x2
s = i6=d |adi |
if add > 0 and add < s then
udd ← − s − p ∗ add
else if add > 0 and add > s then x1
udd ← − 2 ∗ add
Fig. 3. The state space partitioned into three regions. Dr is the demonstration
else if add < 0 and |add | < s then area, while the union Cr ∪ Br is a contracting region.
udd ← − s
else 2.4. Accurate Reproduction of Stable Motions
udd ← 0 // no need to modify add
We propose a modification of the control law in (5) which
end if improves the reproduction accuracy and preserves the conver-
end for gence properties of the controlled DS. As shown in Fig. 3,
return U the state space is divided into three regions: the contraction
region Cr , a ball of radius r centered at the target Br and the
demonstration area Dr . Given Cr , Br and Dr , we propose a
P stabilizing control input:
1 guarantees that µ1 (C) , maxj (cjj + i6=j |cij |) < 0. In !
K K
details, Algorithm 1 inspects A row by row. For each row, X   X  
ẋ = h (x) Ak x + bk + ω(x, t)
k
hk (x) U k x − bk
the algorithm
Pm computes the sum of the off-diagonal elements k=1 k=1
s = i6=d |adi | and considers four cases: i) if the d-th
| {z } | {z }
System Control
diagonal element of A is positive (add > 0) and add < c (6)
then udd = −s − p ∗ add . ii) If add > 0 and add < s The activation function ω(x, t) ∈ R in (6) is computed as:
then udd = − 2 ∗ add . iii) If add < 0 and |add | < s then 
ω̇(x, t) = −γ(ω(x, t) − c(x)) t < tmax
udd = −s. iv) In all the other cases there is no need to (7)
ω(x, t) = 1 t ≥ tmax
modify A. It is easy to verify that add + udd < 0 and that
|add + udd | > c, d = 1, .., m, i.e. µ1 (A + U ) < 0. where c(x) = 1 if x ∈ Cr ∪ Br , c(x) = 0 if x ∈ Dr and
The controller in (5) guarantees the convergence towards ω(x0 , t0 ) = c(x0 ), where x0 is a given initial state.
a unique equilibrium. The control gains are automatically Theorem 2. The dynamical system in (6), with activation
computed given the GMM parameters and the DS is stabilized function (7), globally asymptotically converges to x?
at run time. Figure 2 shows qualitative results of the proposed
approach when applied to stabilize a GMR model. The control Proof. Being interested in the asymptotic stability, we have
input affects the reproduced point-to-point motions and do not to analyze the convergence properties of (6) for t → ∞. We
guarantee an accurate reproduction (see Fig. 2). The reason is can notice that, for ω(x, t) = 1, the controlled DS in (6) and
that the controlled DS (5) converges with an exponentially (5) are the same. Recalling that (5) is globally exponentially
stable dynamics which contradicts the demonstrations. More stable, and considering that ω(x, t) = 1 for t ≥ tmax , we can
formally, ∃β > 0 such that kx(t)−x? k1 ≤ βkx(0)−x? k1 e−ct , conclude that (6) globally asymptotically converges to x? .
where k · k1 is the l1 -norm and µ1 (Ak + U k ) ≤ −c < 0, ∀k. The value of γ > 0 in (7) can be chosen considering that, in
An approach to alleviate this issue is presented in the next practice, ω(x, t) = c(x) after 5/γ s. We activate the controller
section. by setting c(x) = 1 inside the region Cr ∪ Br . Trajectories
that start in Cr ∪ Br follow an exponentially convergent path
200 unless they reach the equilibrium x? or they enter Dr . The
control input is smoothly deactivated (c(x) = 0) inside the
demonstration area Dr to allow an accurate reproduction. The
100 control action is activated if the system has not converged
x2 [mm]

within tmax seconds. This guarantees that the generated mo-

tion trajectories reach the target from any initial state. In our
0 experiments, we set tmax to twice the maximum time duration
of the demonstrations. This prevents the DS to converge to
Streamlines Target
spurious equilibria inside Dr for badly initialized GMR. As
Demonstrations Generated motion
-100 empirically shown in Sec. 3, a GMR with a proper number of
-50 0 50 100 150
x1 [mm] components in general does not generate spurious attractors
Fig. 2. A globally contracting motion generated by the controlled system
within Dr . For each x ∈ Dr , in fact, the GMR generates
(5). The control action guarantees global convergence to x? = [0, 0]T , but it a velocity close to the demonstrated one(s). If demonstrated
affects the accurate reproduction of the demonstrated data. velocities are zero only at the target, the generated velocity
will not drop to zero far from the given target. It is worth Streamlines Demonstrations Target Generated motion
noticing that eventual spurious attractors will only affect the
accuracy of the proposed approach, but not its stability.
The demonstration area Dr is constructed off-line using

x2
the approach in [6]. Each training point D = {xt,i }T,N t=1,i=1
is clustered into K clusters Dk , where K is the number of
Gaussian components.

A training point

xt,i belongs to Dk if
k k j j
N xt,i |µx , Σx > N xt,i |µx , Σx , ∀j 6= k, where N (·) is

x2
the probability that xt,i is generated from the k-th component.
Given the K clusters, a set of K scalars is computed as:
δ k = α min N xt,i |µkx , Σkx


(8)
xt,i ∈D k

x2
where 0 < α ≤ 1 is a constant. δ k is the probability of
the point in Dk located at the maximum distance
from µk .
Hence, the region Dtot = {xt,i : N xt,i |µkx , Σkx ≥ δ k }
contains all training data points. The demonstration area Dr

x2
is then Dr = Dtot \Br . The scalar α in (8) defines the
area of Dr . Bigger values of α result in tighter Dr around
the demonstrations. In all the experiments, we found that
α = 0.1 prevents trajectories starting in Dr to exit from the

x2
demonstration area. It is worth noticing that the construction
of Dr is computationally efficient. Indeed, the probabilities
x1 x1 x1 x1


N xt,i |µkx , Σkx in (8) are used to train the GMR and come
at no extra cost. Therefore, we just have to find the minimum
value for K vectors, multiply the obtained K values by a Fig. 4. Stable motions learned by C-GMR on the LASA dataset.
scalar α, and store the results in δ k , k = 1, .., K. target. Hence, the area among Di and T is a suitable measure
At run time, we firstly check if a state x belongs to Br , i.e. of the reproduction error. Since Di and T have a different
if kx − x? k ≤ r. If x is not part of Br , we have to determine number of points, we first leverage multi-dimensional dynamic
whether x belongs to Dr or
Cr . The state x ∈ Dr if there exists time warping (DTW) [19] to find the optimal non-linear match
a k such that N x|µkx , Σkx ≥ δ k , k = 1, .., K, otherwise x ∈ between the two trajectories and then compute the area among
Cr . This procedure is computationally

efficient. Indeed, the them. The area metric measures how well the learned DS
probabilities N x|µkx , Σkx , k = 1, .., K, are used to compute preserves the shape of the motion. To measure how the DS
hk (x) (see Tab. 1) and come at no extra cost.
Therefore, one preserves the kinematics qof the demonstrations, we use the
has only to check if the K values N x|µkx , Σkx ≥ δ k . In case PT,N
velocity error [8] Ve = N1T t,i=1 kẋt,i − f (xt,i )k2 .
x ∈ Cr ∪ Br or t > tmax , the controller is activated, which
only requires to sum up the K matrices Ak and U k . 3.2. Reproduction accuracy
We compare C-GMR with SEDS and SEDSII in terms
3. E XPERIMENTAL R ESULTS
of reproduction error. The sampling time for reproduction
In this section, we evaluate the performance of Contract- is set to 2ms. SEDSII requires to predefine the form of
ing GMR (C-GMR), in terms of reproduction accuracy and the control Lyapunov function (CLF) depending on the form
training time. In order to underline basic differences between of the motion. As we chose not to adapt any parameter to
C-GMR and the Lyapunov based approaches, we perform a each motion separately but to run the experiments unsuper-
comparison with SEDS [7] and SEDSII [9]. vised, we tested two different parametrization of CLF, namely
3.1. Benchmark and metrics CLF0 and CLF3. CLF0 has no asymmetric components, i.e
CLF 0 = xT P 0 x. CLF3 has three asymmetric components,
The LASA Handwritten1 dataset is used as a benchmark to P3
2
i.e. CLF 3 = CLF 0 + l=1 β l (x) xT P l (x − µl ) . As
test stable motion generation with DS. The dataset contains
detailed in [9], β l , P l and µl are learned from demonstrations
20 motions in 2D. Each motion ends in x? = [0, 0]T and is
by solving a constrained optimization problem.
demonstrated three times. Data are collected at 50 Hz.
A DS for each motion in the dataset is learned by GMR
Two metrics are used to measure the reproduction accuracy.
considering all available demonstrations. The dataset contains
One is the area between each demonstration Di ∈ D and the
motions with different complexity. Complex motions usually
trajectory T generated by C-GMR (SEDS or SEDSII) starting
require many components to be accurately represented. To
from the first point in Di . Both Di and T end at the same
maximize the accuracy, the number of components is com-
1 The source code for SEDS, SEDSII and the LASA dataset are available puted by means of Bayesian information criterion (BIC) [20].
on-line: [Link]/sourcecode/ C-GMR is used to stabilize f (x) at run time. As qualitatively
 101 102 0.6 0.22

Time [ms]
Velocity error [m/s]
0.2
0.3
100 101
Area [m2]

0.18
0
1 5 10 15 20 25 30 2 4 6 8 10
Gaussian components (#) Space dimension (#)
10-1 100
Fig. 6. Execution time of Algorithm 1, averaged over 100 executions.
10-2
SEDS SEDSII SEDSII C-GMR SEDS SEDSII SEDSII C-GMR
(CLF0) (CLF3) (CLF0) (CLF3) K, the right graph by fixing K = 10 and varying T . The
(a) (b)
execution time is always below 0.1 ms.
0.08
0.08
Fig. 5. Overall accuracy for SEDS, SEDII and C-GMR on the LASA dataset.

Time [ms]
0.07
Black error bars are 10% and 90% quantiles of the median value.
0.04
0.06
shown in Fig. 4, C-GMR accurately learns complex point-to- 0.01
0.05 x100
1 5 10 15 20 25 30 1 2 3 4 5 6 7 8 9 10
point motions while guaranteeing global convergence. Gaussian components (#) Training data (#)

GMR, SEDS and SEDSII require the numerical solution of

Fig. 7. Time required to compute Dr , averaged over 100 executions.
optimization problems, for which different initial conditions
lead to different results. In order to evaluate the typical We use the same configuration of the previous experiment
performance, we tested each approach 10 times for each to compare the training times of C-GMR, SEDS and SED-
motion. Overall accuracy is shown in Fig. 5. For a better SII. Overall training time is shown in Fig. 8. For a better
visualization, we use a logarithmic scale for the ordinate visualization, we use a logarithmic scale for the ordinate axis.
axis. As the reproduction error is not normally distributed, SEDS training step requires about 6.5 s (median value) for
we consider the median instead of the mean. To indicate the each motion, while SEDSII CLF0 and CLF3 spend 0.76 s and
maximal and minimal deviation from the typical performance, 2.7 s respectively. C-GMR has a (best) median training time
we provide the location of the 10% and the 90% quantiles of 0.58 s. It is worth noticing that C-GMR time in Fig. 8
(black error bars in Fig. 5). Results for C-GMR are obtained includes the time needed to compute the region Dr and to run
by setting the radius of Br equal to the 15% of the distance Algorithm 1. As shown in the previous experiment, this time
between the starting point and the target. is at least three orders of magnitude smaller than GMM time.
Quantitative results in Fig. 5 show that SEDS has median
error area of 1.3m2 and a median velocity error of 8.3m/s 102 GMM time
respectively. This means that most of the motions in the CLF time
Training Time [s]

dataset cannot be learned with the Lyapunov function xT x 101

used in SEDS. SEDSII CLF0 and CLF3 have almost the
same (best) median accuracy (0.596 m2 and 0.591 m2 for the
100
area error, 5.2 m/s and 4.9 m/s for Ve ), meaning that most
of the considered motions can be learned with the Lyapunov
function xT P 0 x. The proposed C-GMR has a smaller median SEDS SEDSII SEDSII C-GMR
accuracy (0.768 m2 and 6.1 m/s) compared to SEDSII, but (CLF0) (CLF3)
higher accuracy than SEDS. C-GMR loses accuracy when the
state reaches the ball Br around the target. Inside Br , in fact, Fig. 8. Overall training time for SEDS, SEDII and C-GMR on the LASA
dataset. Black error bars are 10% and 90% quantiles of the median value.
the motion follows an exponentially converging dynamics.
3.4. Generalization to different initial/target positions
3.3. Training time
Motion generation based on DS has the advantage to be
We investigate the extra training time introduced by C- robust to changes in the initial and target positions. In partic-
GMR, i.e. the time required to execute Algorithm 1 and to ular, stable DS are able to generate convergent trajectories for
compute Dr . The execution time of Algorithm 1 depends any initial/target position. The stability of C-GMR is mathe-
on the number K of Gaussians and on the dimension m of matically proven in Theorem 2. Nevertheless, it is interesting
the state space (A, U ∈ Rm×m ). The execution time of our to qualitatively show how the learned DS adapts to different
unoptimized Matlab implementation of Algorithm 1 is shown initial/target positions. As a proof of concept, we train C-GMR
in Fig. 6. The left graph is obtained by fixing m = 2 (as in on the S-shape motion in the LASA dataset. Figure 9(a) shows
the LASA dataset) and varying K, the right graph by fixing the trajectories generated from different initial positions. As
K = 10 and varying m. The execution time is below 0.6 ms. expected, all the trajectories converge towards the target. Black
The time required to construct Dr depends on K and on solid lines are trajectories that start or enter in Dr . Figure 9(a)
the number of training points T . The execution time of our shows how the system smoothly transits between the different
unoptimized Matlab implementation is shown in Fig. 7. The region. This is due to the smooth activation function in (7).
left graph is obtained by keeping T = 300 fixed and varying Black dashed lines are trajectories that start and remain in
 Generated motion Demonstrations Target Initial position overall system. The resulting algorithm, called Contracting
250 250 GMR (C-GMR), has been tested on a public dataset consisting
of complex 2D motion, showing promising results in terms of
accuracy and training time. We also performed an experimental
x2 [mm]

x2 [mm]
75 75
comparison with Lyapunov based approaches, useful to under-
stand the main differences between the different approaches.
-100 -100
Our next research will focus on extending C-GMR to other
-50 75 200 -50 75 200
x1 [mm] x1 [mm] regression techniques, and on testing C-GMR on real robots.
(a) (b)
ACKNOWLEDGEMENTS
Fig. 9. Generalization of the S-shape motion. (a) Trajectories generated from This work has been supported by the Technical University
different initial positions. (b) Trajectories converging to a different goal.
of Munich, International Graduate School of Science and
Cr ∪ Br and converge exponentially. Figure 9(b) shows DS Engineering (IGSSE).
trajectories converging to the new goal x? = [−25, 25]T .
R EFERENCES
3.5. Discussion [1] A. Billard, S. Calinon, R. Dillmann, and S. Schaal, “Robot programming
by demonstration,” in Springer Handbook of Robotics. Springer Berlin
The experimental evaluation shows that the proposed C- Heidelberg, 2008, pp. 1371–1394.
GMR is less accurate than SEDSII. This is because SEDSII [2] D. Lee and C. Ott, “Incremental kinesthetic teaching of motion prim-
is optimized along the entire motion, while C-GMR loses itives using the motion refinement tube,” Autonomous Robots, vol. 31,
no. 2, pp. 115–131, 2011.
accuracy when the state reaches the ball Br (see Fig. 3) around [3] M. Saveriano, S. An, and L. Lee, “Incremental kinesthetic teaching
the target. In general, contraction theory allows to consider a of end-effector and null-space motion primitives,” in International
coordinate transformation Θ(x). In our future research, we will Conference on Robotics and Automation, 2015, pp. 3570–3575.
[4] S. Calinon, F. Guenter, and A. Billard, “On learning, representing and
focus on learning from demonstration a suitable coordinate generalizing a task in a humanoid robot,” Transactions on Systems, Man,
transformation Θ(x) that maps the input space into a space and Cybernetics. Part B, vol. 37, no. 2, pp. 286–298, 2007.
where demonstrations can be more accurately reproduced. [5] A. Ijspeert, J. Nakanishi, P. Pastor, H. Hoffmann, and S. Schaal,
“Dynamical Movement Primitives: learning attractor models for motor
Regarding the training time, C-GMR outperforms SEDSII. behaviors,” Neural Computation, vol. 25, no. 2, pp. 328–373, 2013.
C-GMR, in fact, stabilizes the dynamical system at run time, [6] S. M. Khansari-Zadeh and A. Billard, “BM: An iterative algorithm
to learn stable non-linear dynamical systems with Gaussian Mixture
without learning steps. The control input is automatically Models,” Int. Conf. on Rob. and Aut., pp. 2381–2388, 2010.
computed given GMR parameters, and (as shown in Sec. [7] ——, “Learning stable non-linear dynamical systems with Gaussian
3.3), the time to compute the control input is practically Mixture Models,” Trans. on Rob., vol. 27, no. 5, pp. 943–957, 2011.
[8] A. Lemme, F. Reinhart, K. Neumann, and J. J. Steil, “Neural learning of
negligible. A possible application of C-GMR is incremental vector fields for encoding stable dynamical systems,” Neurocomputing,
learning, where GMR parameters are continuously updated as vol. 141, pp. 3–14, 2014.
new demonstrations of the task are provided. On the contrary, [9] S. M. Khansari-Zadeh and A. Billard, “Learning control Lyapunov
function to ensure stability of dynamical system-based robot reaching
SEDSII is preferable in applications which require high accu- motions,” Rob. And Auton. Systems, vol. 62, no. 6, pp. 752–765, 2014.
racy and for which the training time is not a limitation. [10] K. Neumann and J. J. Steil, “Learning robot motions with stable
A general assumption behind stable motion generation with dynamical systems under diffeomorphic transformations,” Robotics and
Autonomous Systems, vol. 70, pp. 1–15, 2015.
autonomous DS is that the demonstrations do not explicitly [11] P. Perrin and P. Schlehuber-Caissier, “Fast diffeomorphic matching
show spurious equilibria. In case demonstrated velocities drop to learn globally asymptotically stable nonlinear dynamical systems,”
to zero far from the given target, the desired behavior is Systems & Control Letters, vol. 96, pp. 51–59, 2016.
[12] M. Saveriano and L. Lee, “Point cloud based dynamical system modu-
probably to stop for a certain time and then to reach the goal. lation for reactive avoidance of convex and concave obstacles,” in Intl
A unique autonomous DS is not suitable to execute this task. Conf. on Intelligent Robots and Systems, 2013, pp. 5380–5387.
A solution is to use one autonomous DS for each attractor and [13] M. Saveriano, F. Hirt, and L. Lee, “Human-aware motion reshaping
using dynamical systems,” Pattern Recognition Letters, 2017.
to switch among them as soon as one equilibrium is reached. [14] M. Saveriano and L. Lee, “Distance based dynamical system modulation
for reactive avoidance of moving obstacles,” in Intl Conf. on Robotics
4. C ONCLUSION and Automation, 2014, pp. 5618–5623.
[15] Z. Huang and C. Siew, “Extreme learning machine: Theory and appli-
This work investigated the possibility of applying Con- cations,” Neurocomputing, vol. 70, no. 1-3, pp. 489–501, 2006.
traction theory to learn point-to-point motions via dynamical [16] W. Lohmiller and J. Slotine, “On Contraction analysis for nonlinear
systems and Gaussian mixture regression. We derived suffi- systems,” Automatica, vol. 34, no. 6, pp. 683–696, 1998.
[17] H. Ravichandar and A. Dani, “Learning contracting nonlinear dynamics
cient stability conditions for a dynamical system represented from human demonstrations for robot motion planning,” in Dynamics,
by GMR by applying contraction analysis. An automatic and Systems and Control Conference, 2015.
computationally efficient procedure is developed to compute [18] W. Wang and J. Slotine, “On Partial Contraction analysis for coupled
nonlinear oscillators,” Biol. Cybern., vol. 92, no. 1, pp. 38–53, 2005.
a stabilizing control input given the GMR parameters. At run [19] P. Sanguansat, “Multiple multidimensional sequence alignment using
time, the control action is smoothly activated or deactivated generalized dynamic time warping,” Transactions on Mathematics,
using an activation function. The combination of the control vol. 11, no. 8, pp. 668–678, 2012.
[20] G. Schwarz, “Estimating the dimension of a model,” The annals of
action and the activation function significantly improves the statistics, vol. 6, no. 2, pp. 461–464, 1978.
reproduction accuracy without affecting the stability of the