Provably Robust Learning-Based Approach for High-Accuracy
Tracking Control of Lagrangian Systems

arXiv:1804.01031v1 [cs.RO] 3 Apr 2018

Mohamed K. Helwa, Adam Heins, and Angela P. Schoellig
Abstract— Inverse dynamics control and feedforward linearization techniques are typically used to convert the complex
nonlinear dynamics of robots to a set of decoupled double
integrators, and then a standard, outer-loop controller can be
used to calculate the commanded acceleration for the linearized
system. However, these methods typically depend on having a
very accurate system model, which is not available in practice.
While this challenge has been addressed in the literature using
different learning approaches, most of these approaches do
not provide safety guarantees in terms of stability of the
learning-based control system. In this paper, we provide a novel,
learning-based control approach based on Gaussian processes
(GPs) that ensures both stability of the closed-loop system and
high-accuracy tracking. In our proposed approach, we use GPs
to approximate the error between the commanded acceleration
and the actual acceleration of the system, and then use the
predicted mean and variance of the GP to calculate an upper
bound on the uncertainty of the linearization. This uncertainty
bound is then used in a robust, outer-loop controller to ensure
stability of the overall system. Moreover, we show that using our
proposed strategy, the tracking error converges to a ball with a
radius that can be made arbitrarily small through appropriate
control design. Furthermore, we verify the effectiveness of
the proposed approach via simulations on (i) a 2 degree-offreedom (DOF) planar manipulator, and (ii) a 6 DOF industrial
manipulator using the Gazebo simulator.

I. I NTRODUCTION
High-accuracy tracking is an essential requirement in
advanced manufacturing, self-driving cars, medical robots,
and autonomous flying vehicles, among other application
areas. To achieve high-accuracy tracking for these complex, typically high-dimensional, nonlinear robotic systems,
a typical approach is to use inverse dynamics control [1]
or feedforward linearization techniques [2] to convert the
complex nonlinear dynamics into a set of decoupled double
integrators, and then a standard, linear, outer-loop controller,
e.g., a proportional-derivative (PD) controller, can be used to
make the decoupled linear system track the desired trajectory [1]. However, these linearization techniques depend on
having accurate system models, which are difficult to obtain
from the first principles of physics.
To address this problem, for many decades robust control techniques have been used to design the outer-loop
controllers to account for the uncertainties in the linearization [3]. However, the selection of the uncertainty bounds in
the robust controller design is challenging. On the one hand,
The authors are with the Dynamic Systems Lab (www.dynsyslab.org),
Institute
for
Aerospace
Studies,
University
of
Toronto,
Canada.
E-mail:
mohamed.helwa@robotics.utias.utoronto.ca,
adam.heins@robotics.utias.utoronto.ca,
schoellig@utias.utoronto.ca.
This research was supported by NSERC grant RGPIN-2014-04634 and
OCE/SOSCIP TalentEdge Project #27901.

Fig. 1. Block Diagram of our proposed control strategy; GP regression
models learn the uncertainty in the linearization, and then using the
mean and variance of the GP, one can calculate an upper bound on the
uncertainty to be used in a robust, outer-loop controller. The symbol q
is the actual position vector, qd is the desired position vector, aq is the
commanded acceleration vector, u is the force/torque input vector, and η is
the uncertainty vector.

selecting high bounds typically results in a conservative
behavior, and hence, a large tracking error. On the other
hand, relatively small uncertainty bounds may not represent
the true upper bounds of the uncertainties, and consequently,
stability of the overall system is not ensured. Alternatively,
several approaches have been proposed for learning the
inverse system dynamics from collected data for the practical
cases where the system models are not available or not
sufficiently accurate; see [4]–[7]. Combining a-priori model
knowledge with learning data has also been studied in [4],
[8]. However, these learning approaches typically neglect the
learning regression errors in the analysis, and they do not
provide a proof of stability of the overall, learning-based control system, which is crucial for safety-critical applications
such as medical robots. The limitations of the robust control
and the learning-based techniques show the urgent need for
novel, robust, learning-based control approaches that ensure
both stability of the overall control system and high-accuracy
tracking. This sets the stage for the research carried out in
this paper.
In this paper, we provide a novel, robust, learning-based
control technique that achieves both closed-loop stability
and high-accuracy tracking. In particular, we use Gaussian
processes (GPs) to approximate the error between the commanded acceleration to the linearized system and the actual
acceleration of the robotic system, and then use the predicted
mean and variance of the GP to calculate an upper bound on
the uncertainty of the linearization. This uncertainty bound is
then used in a robust, outer-loop controller to ensure stability
of the overall system (see Figure 1). Moreover, we show that
using our proposed strategy, the tracking error converges to a
ball with a radius that can be made arbitrarily small through
appropriate control design, and hence, our proposed approach
also achieves high-accuracy tracking. Furthermore, we verify

the effectiveness of the proposed approach via simulations on
(i) a 2 DOF planar manipulator using MATLAB Simulink,
and (ii) a 6 DOF industrial manipulator using the Gazebo
simulator.
This paper is organized as follows. Section II provides a
summary of some recent related work. Section III describes
the considered problem, and Section IV provides the proposed approach for solving the problem. Section V provides
theoretical guarantees of the effectiveness of the proposed
approach. Section VI provides the simulation results, and
Section VII concludes the paper.
Notation and Basic Definitions: For a set S, S̄ denotes
its closure and S ◦ its interior. The notation Bδ (y) denotes a
ball of radius δ centered at a point y. A matrix P is positive
definite if it is symmetric and all its eigenvalues are positive.
For a vector x, kxk denotes its Euclidean norm. We say a
function f (x) is smooth if its partial derivatives of all orders
exist and are continuous. We say that the solutions of ẋ =
f (t, x) are uniformly ultimately bounded with ultimate bound
b if there exist positive constants b, c, and for every 0 < a <
c, there exists T (a, b) ≥ 0 such that kx(t0 )k ≤ a implies
kx(t)k ≤ b, for all t ≥ (T + t0 ), where t0 is the initial time
instant. A kernel is a symmetric function k : A × A → R.
A reproducing kernel Hilbert space (RKHS) corresponding
to
a kernel k(., .) includes functions of the form f (a) =
Pm
j=1 αj k(a, aj ) with m ∈ N, αj ∈ R and representing
points aj ∈ A.
II. R ELATED W ORK
The study of safe learning dates back to the beginning
of this century [9]. In [10] and [11], Lyapunov-based reinforcement learning is used to allow a learning agent to safely
switch between pre-computed baseline controllers. Then, in
[12], risk-sensitive reinforcement learning is proposed, in
which the expected return is heuristically weighted with the
probability of reaching an error state. In several other papers,
including [13], [14] and [15], safe exploration methods are
utilized to allow the learning modules to achieve a desired
balance between ensuring safe operation and exploring new
states for improved performance. In [9], a general framework
is proposed for ensuring safety of learning-based control
strategies for uncertain robotic systems. In this framework,
robust reachability guarantees from control theory are combined with Bayesian analysis based on empirical observations. The result is a safety-preserving, supervisory controller
of the learning module that allows the system to freely
execute its learning policy almost everywhere, but imposes
control actions to ensure safety at critical states. Despite its
effectiveness for ensuring safety, the supervisory controller
in this approach has no role in reducing tracking errors.
Focusing our attention on safe, learning-based inverse
dynamics control, we refer to [16], [17]. In these papers,
Gaussian processes (GPs) are utilized to learn the errors in
the output torques of the inverse dynamics model online.
In [16], the GP learning is combined with a state-of-the-art
gradient descent method for learning feedback terms online.
The main idea behind this approach is that the offsets learnt

online by the gradient descent method would correct for fast
perturbations, while the GP is responsible for correcting slow
perturbations. This allows for exponential smoothing of the
GP hyperparameters, which increases the robustness of the
GP at the cost of having slower reactiveness. Nevertheless,
[16] does not provide a proof of the robust stability of
the closed-loop system. In [17], the variance of the GP
prediction is utilized to adapt the parameters of an outer-loop
PD controller online, and the uniform ultimate boundedness
of the tracking error is proved under some assumptions on
the structure of the PD controller (e.g., the gain matrix
was assumed to be diagonal, which imposes a decentralized
gain control scheme). The results of [17] are verified via
simulations on a 2 DOF manipulator.
Our approach differs from [17] in several aspects. First,
we do not use an adaptive PD controller in the outer loop,
but add a robustness term to the output of the outer-loop
controller. Second, while [17] uses the GP to learn the error
in the estimated torque from the nominal inverse dynamics,
in our approach, we learn the error between the commanded
and actual accelerations. This can be beneficial in two
aspects: (i) This makes our approach applicable to industrial
manipulators that have onboard controllers for calculating the
torque and allow the user to only send commanded acceleration/velocity; (ii) this makes our approach also applicable
beyond inverse dynamics control of manipulators; indeed,
our proposed approach can be applied to any Lagrangian
system for which feedforward/feedback linearization can be
used to convert the nonlinear dynamics of the system to
a set of decoupled double integrators, such as a quadrotor
under a feedforward linearization, see Section 5.3 of [18].
Third, while [17] shows uniform ultimate boundedness of
the tracking error, it does not provide discussions on the
size of the ultimate ball. In this work, we show that using
our proposed approach, the size of the ball can be made
arbitrarily small through the control design. Fourth, in our
approach, we do not impose any assumption on the structure
of the outer-loop PD controller and decentralized, outer-loop
control is not needed for our proof. Finally, we verify our
approach using simulations on a 6 DOF manipulator.
III. P ROBLEM S TATEMENT
In this paper, we consider Lagrangian systems, which
represent a wide class of mechanical systems [19]. In what
follows, we focus our attention on a class of Lagrangian
systems represented by:
M (q(t))q̈(t) + C(q(t), q̇(t))q̇(t) + g(q(t)) = u(t),

(1)

where q = (q1 , · · · , qN ) is the vector of generalized coordinates (displacements or angles), q̇ = (q̇1 , · · · , q̇N ) is
the vector of generalized velocities, u = (u1 , · · · , uN ) is
the vector of generalized forces (forces or torques), N is
the system’s degree of freedom, M , C, and g are matrices
of proper dimensions and smooth functions, and M (q) is a
positive definite matrix. Fully-actuated robotic manipulators
are an example of Lagrangian systems that can be expressed

by (1). Despite focusing our discussion on systems represented by (1), we emphasize that our results in this paper can
be easily generalized to a wider class of nonlinear Lagrangian
systems for which feedback/feedforward linearization can be
utilized to convert the dynamics of the system into a set of
decoupled double integrators, as we will discuss in detail in
Remark 5.3.
For the nonlinear system (1) with uncertain matrices M ,
C, and g, we aim to make the system positions and velocities
(q(t), q̇(t)) track a desired smooth trajectory (qd (t), q̇d (t)).
For simplicity of notation, in our discussion, we drop the
dependency on t from q, qd , their derivatives, and u. Our
goal is to design a novel, learning-based control strategy
that is easy to interpret and implement, and that satisfies
the following desired objectives:
(O1) Robustness: The overall, closed-loop control system
satisfies robust stability in the sense that the tracking
error has an upper bound under the system uncertainties.
(O2) High-Accuracy Tracking: For feasible desired trajectories, the tracking error converges to a ball around
the origin that can be made arbitrarily small through
the control design. For the ideal case, where the preassumed system parameters are correct, the tracking
error should converge asymptotically to the origin.
(O3) Adaptability: The proposed strategy should incorporate
online learning to continuously adapt to online changes
of the system parameters and disturbances.
(O4) Generalizability of the Approach: The proposed approach should be general enough to be also applicable
to industrial robots that have onboard controllers for
calculating the forces/torques and allow the user to
send only commanded acceleration/velocity. It is also
desirable that the approach is general enough to be
implementable on other Lagrangian systems beyond (1).
IV. M ETHODOLOGY
In this section, we present our proposed methodology, and
then in the next section, we show that it satisfies objectives
(O1)-(O4).
A standard approach for solving the tracking control
problem for (1) is inverse dynamics control. In particular, if
the matrices M , C, and g are all known, then the following
inverse dynamics control law
u = C(q, q̇)q̇ + g(q) + M (q)aq

(2)

converts the complex nonlinear dynamic system (1) into
q̈ = aq ,

(3)

where aq is the commanded acceleration, a new input to
the linearized system (3) to be calculated by an outer-loop
control law, e.g., a PD controller (see Figure 1). However,
the standard inverse dynamics control (2) heavily depends on
accurate knowledge of the system parameters. In practice,
the matrices M , C, and g are not perfectly known, and
consequently, one has to use estimated values of these
matrices M̂ , Ĉ, and ĝ, respectively, where M̂ , Ĉ, and ĝ

are composed of smooth functions. Hence, in practice, the
control law (2) should be replaced with
u = Ĉ(q, q̇)q̇ + ĝ(q) + M̂ (q)aq .

(4)

Now by plugging (4) into the system model (1), we get
q̈ = aq + η(q, q̇, aq ),

(5)

where η(q, q̇, aq ) = M −1 (q)(M̃ (q)aq + C̃(q, q̇)q̇ + g̃(q)),
with M̃ = M̂ − M , C̃ = Ĉ − C, and g̃ = ĝ − g.
Because of the term η, the dynamics (5) resulting from the
inverse dynamics control are still nonlinear and coupled. To
control the uncertain system (5), on the one hand, robust
control methods are typically very conservative, while on the
other hand, learning methods do not provide robust stability
guarantees.
Hence, in this paper, we combine ideas from robust
control theory with ideas from machine learning, particularly
Gaussian processes (GPs) for regression, to provide a robust,
learning-based control strategy that satisfies objectives (O1)(O4). The main idea behind our proposed approach is to use
GPs to learn the uncertainty vector η(q, q̇, aq ) in (5) online.
Following [17], we use a set of N GPs, one for learning each
element of η, to reduce the complexity of the regression.
A main advantage of GP regression is that it does not only
provide an estimated value of the mean µ, but it also provides
an expected variance σ 2 , which represents the accuracy of
the regression model based on the distance to the training
data. The punchline here is that one can use both the mean
and variance of the GP to calculate an upper bound ρ on kηk
that is guaranteed to be correct with high probability, as we
will show later in this section. One can then use this upper
bound to design a robust, outer-loop controller that ensures
robust stability of the overall system. Hence, our proposed
strategy consists of three parts:
(i) Inner-loop controller: We use the inverse dynamics
control law (4), where M̂ , Ĉ, and ĝ are estimated values of
the system matrices from an a-priori model.
(ii) GPs for learning the uncertainty: We use a set of
N GPs to learn the uncertainty vector η in (5). We start by
reviewing GP regression [15], [20]. A GP is a nonparametric
regression model that is used to approximate a nonlinear
function J(a) : A → R, where a ∈ A is the input vector.
The ability of the GP to approximate the function is based
on the assumptions that function values J(a) associated with
different values of a are random variables, and that any finite
number of these variables have a joint Gaussian distribution.
The GP predicts the value of the function, J(a∗ ), at an
arbitrary input a∗ ∈ A from a set of n observations Dn :=
ˆ j )}n , where J(a
ˆ j ), j ∈ {1, · · · , n}, are assumed
{aj , J(a
j=1
to be noisy measurements of the function’s true values. That
ˆ j ) = J(aj ) + ω 0 , where ω 0 is a zero mean Gaussian
is, J(a
noise with variance σω2 . Assuming, without loss of generality
(w.l.o.g.), a zero prior mean of the GP and conditioned on
the previous observations, the mean and variance of the GP
prediction are given by:
µn (a∗ ) = kn (a∗ )(Kn + In σω2 )−1 Ĵn ,

(6)

σn2 (a∗ ) = k(a∗ , a∗ ) − kn (a∗ )(Kn + In σω2 )−1 kTn (a∗ ), (7)
ˆ 1 ), · · · , J(a
ˆ n )]T is the vector
respectively, where Ĵn = [J(a
of observed, noisy function values. The matrix Kn ∈ Rn×n
is the covariance matrix with entries [Kn ](i,j) = k(ai , aj ),
i, j ∈ {1, · · · , n}, where k(ai , aj ) is the covariance function defining the covariance between two function values
J(ai ), J(aj ) (also called the kernel). The vector kn (a∗ ) =
[k(a∗ , a1 ), · · · , k(a∗ , an )] contains the covariances between
the new input and the observed data points, and In ∈ Rn×n
is the identity matrix. The tuning of the GP is typically done
through the selection of the kernel function and the tuning of
its hyperparameters. For information about different standard
kernel functions, please refer to [20].
We next discuss our implementation of the GPs. The GPs
run in discrete time with sampling interval Ts . The inputs to
each GP regression model are the same (qk , q̇k , aq k ), and the
output is an estimated value of an element of the η vector
at the sampling instant k. For the training data for each GP,
n observations of (q, q̇, aq ) are used as the labeled input
together with n observations of an element of the vector
q̈ − aq + ω as the labeled output, where ω ∈ RN is Gaussian
noise with zero mean and variance diag(σω2 1 , · · · , σω2 N ); see
(5). In this paper, we use the squared exponential kernel
1

k(ai , aj ) = ση2 e− 2 (ai −aj )

T

M̄ −2 (ai −aj )

,

(8)

which is parameterized by the hyperparameters: ση2 , the prior
variance, and the positive length scales l1 , · · · , l3N which are
the diagonal elements of the diagonal matrix M̄ . Hence, the
expected mean and variance of each GP can be obtained
through equations (6)-(8). Guidelines for tuning the GP
hyperparameters σω2 , ση2 , l1 , · · · , l3N can be found in [15].
As stated before, a main advantage of GP regression is that
the GP provides a variance, which represents the accuracy of
the regression model based on the distance between the new
input and the training data. One can then use the predicted
mean and variance of the GP to provide a confidence interval
around the mean that is guaranteed to be correct with high
probability. There are several comprehensive studies in the
machine learning literature on calculating these confidence
intervals. For completeness, we review one of these results,
particularly Theorem 6 of [20]. Let aaug = (q, q̇, aq ), and
ηi (aaug ) denote the i-th element of the unknown vector
function η.
Assumption 4.1: The function ηi (aaug ), i ∈ {1, · · · , N },
has a bounded RKHS norm kηi kk with respect to the covariance function k(a, a0 ) of the GP, and the noise ωi added
to the output observations, i ∈ {1, · · · , N }, is uniformly
bounded by σ̄.
The boundedness of the RKHS norm implies that the
functions are well-behaved in the sense that they are regular
with respect to the kernel [20].
Lemma 4.1 (Theorem 6 of [20]): Suppose that Assumption 4.1 holds. Let δp ∈ (0, 1). Then, Pr{∀aaug ∈
A, kµ(aaug ) − ηi (aaug )k ≤ β 1/2 σ(aaug )} ≥ 1 − δp ,
where Pr stands for the probability, A ⊂ R3N is
compact, µ, σ 2 are the GP mean and variance eval-

uated at aaug conditioned on n past observations, and
β = 2kηi k2k + 300γ ln3 ((n + 1)/δp ). The variable γ ∈
R is the maximum information gain and is given by
γ = max{aaug,1 ,··· ,aaug,n+1 }∈A 0.5 log(det(I + σ̄ −2 Kn+1 )),
where det is the matrix determinant, I ∈ R(n+1)×(n+1) is
the identity matrix, Kn+1 ∈ R(n+1)×(n+1) is the covariance
matrix given by [Kn+1 ](i,j) = k(aaug,i , aaug,j ), i, j ∈
{1, · · · , n + 1}.
The punchline here is that we know from Lemma 4.1 that
one can define for each GP a confidence interval around the
mean that is guaranteed to be correct for all points aaug ∈ A,
a compact set, with probability higher than (1 − δp ), where
2
δp is typically picked very small. Let µk,i , σk,i
, and βk,i
represent the expected mean, variance, and β parameter
in Lemma 4.1 of the i-th GP at the sampling instant k,
respectively, where i ∈ {1, · · · , N }. We select the upper
bound on the absolute value of ηi at k to be
1/2

1/2

ρk,i (µk,i , σk,i ) = max(|µk,i − βk,i σk,i |, |µk,i + βk,i σk,i |).
(9)
Then, a good estimate of the upper bound on kηk at k is
q
ρk = ρk,1 (µk,1 , σk,1 )2 + · · · + ρk,N (µk,N , σk,N )2 . (10)
(iii) Robust, outer-loop controller: We use the estimated
upper bound ρk to design a robust, outer-loop controller. In
particular, for a smooth, bounded desired trajectory qd (t), we
use the outer-loop control law
aq (t) = q̈d (t) + KP (qd (t) − q(t)) + KD (q̇d (t) − q̇(t)) + r(t),
(11)
where KP ∈ RN ×N and KD ∈ RN ×N are the proportional
and derivative matrices of the PD control law, respectively,
and r ∈ RN is an added vector to the PD control law that
will be designed to achieve robustness. Let e(t) := (q(t) −
qd (t), q̇(t)−q̇d (t)) denote the tracking error vector. From (11)
and (5), it can be shown that the tracking error dynamics are
ė(t) = Ae(t) + B(r(t) + η(q(t), qd (t), aq (t))),

(12)

where

A=

0
−KP

I
−KD




, B=

0
I


,

(13)

and I ∈ RN ×N is the identity matrix. From (12) and (13), it
is clear that the controller matrices KP and KD should be
designed to make A a Hurwitz matrix.
We now discuss how to design the robustness vector r(t).
To that end, let P be the unique positive definite matrix
satisfying AT P + P A = −Q, where Q is a positive definite
matrix. We define r(t) as follows
(
B T P e(t)
−ρ(t) kB
kB T P e(t)k > ,
T P e(t)k
r(t) =
(14)
T
−ρ(t) B P e(t)
kB T P e(t)k ≤ ,
where ρ(t) is the last received upper bound on kηk from the
GPs, i.e., we use
ρ(t) = ρk , ∀t ∈ [kTs , (k + 1)Ts ],

(15)

and  is a small positive number. It should be noted that 
is a design parameter that can be selected to ensure highaccuracy tracking, as we will discuss in the next section.
V. T HEORETICAL G UARANTEES
After discussing the proposed control strategy, we now
justify that it satisfies the desired objectives (O1)-(O4). To
that end, we require the following reasonable assumptions:
Assumption 5.1: The uncertainty vector η(q, q̇, aq ) in (5)
can be modeled by N independent GPs.
Assumption 5.2: The GPs run at a sufficiently fast sampling rate such that the calculated upper bound on kηk is
accurate between two consecutive sampling instants.
Next, we impose another assumption to ensure that the
added robustness vector r(t) will not cause the uncertainty
vector norm kη(q(t), q̇(t), aq (t))k to blow up. It is easy to
show that the uncertainty function η(q, q̇, aq ) is smooth, and
so kηk attains a maximum value on any compact set in its
input space (q, q̇, aq ). However, since from (11) and (14),
aq is a function of ρ(t), an upper bound on kηk, one still
needs to ensure the boundedness of kηk for bounded q, q̇ or
bounded tracking error e. Hence, we present the following
assumption.
Assumption 5.3: For a given, smooth, bounded desired
trajectory (qd (t), q̇d (t)), there exists ρ̄ > 0 such that kηk < ρ̄
for each e ∈ D, where D is a compact set containing
{e ∈ R2N : eT P e ≤ e(0)T P e(0)}, and e(0) is the initial
tracking error.
We now justify that Assumption 5.3 is reasonable, and
that it is automatically satisfied for small uncertainties in
the system matrices, particularly in the inertia matrix M (q).
In this discussion, we suppose that ρ(t) in (14) satisfies
ρ(t) ≤ kη(t)k + c, where c is a positive scalar, and study
whether imposing r(t) into (11) can make kη(t)k blow up.
Recall that η(q, q̇, aq ) = M −1 (q)(M̃ (q)aq +C̃(q, q̇)q̇+g̃(q)).
From (11), we have η(q, q̇, aq ) = M −1 (q)(M̃ (q)r + M̃ (q̈d +
KP (qd − q) + KD (q̇d − q̇)) + C̃(q, q̇)q̇ + g̃(q)). It is evident
that kηk ≤ kM −1 (q)M̃ (q)kkrk + kM̃ (q̈d + KP (qd − q) +
KD (q̇d − q̇)) + C̃(q, q̇)q̇ + g̃(q))k. From (14), it is easy to
verify kr(t)k ≤ kρ(t)k = ρ(t), and so kr(t)k ≤ (kη(t)k+c).
Hence, kηk ≤ kM −1 (q)M̃ (q)kkηk + kM −1 (q)M̃ (q)kc +
kM̃ (q̈d +KP (qd −q)+KD (q̇d − q̇))+ C̃(q, q̇)q̇+g̃(q))k. Now
if the uncertainty in the matrix M (q), M̃ (q), is sufficiently
small such that kM −1 (q)M̃ (q)k < 1 is satisfied, then we
have kηk ≤ 1−kM −11(q)M̃ (q)k (kM −1 (q)M̃ (q)kc + kM̃ (q̈d +
KP (qd − q) + KD (q̇d − q̇)) + C̃(q, q̇)q̇ + g̃(q))k). Since
qd , q̇d , q̈d are all bounded by assumption, if e(t) ∈ D,
a compact set, then q(t), q̇(t) are also bounded. Then, it is
straightforward to show that there exists a fixed upper bound
ρ̄ on kηk that is valid for each e ∈ D, and Assumption 5.3
is satisfied.
From Assumption 5.3 and our discussion in the previous
paragraph, we know that kη(t)k ≤ ρ̄ if e(t) ∈ D, and
consequently, it is reasonable to saturate any estimate of
ρ(t) beyond ρ̄. Hence, we suppose that the estimation of

ρ is slightly modified to be
ρ(t) = min(ρGP (t), ρ̄),

(16)

where ρGP (t) is the upper bound on the uncertainty norm,
kηk, calculated from the GPs in (9), (10), and (15). It is
straightforward to show that with the choice of ρ(t) in (16)
and for bounded smooth trajectories, the condition e(t) ∈ D
for all t ≥ 0 implies that aq (t) in (11) is always bounded,
and so aaug = (q, q̇, aq ) always lies in a compact set. To
be able to provide theoretical guarantees, we also assume
w.l.o.g. that the small positive number  in (14) is selected
sufficiently small such that
r
ρ̄
 δ1 ,
(17)
2λmin (Q)
where δ1 > 0 is such that Bδ1 (0) ⊂ {e ∈ R2N : eT P e <
e(0)T P e(0)}, and λmin (Q) > 0 is the smallest eigenvalue
of the positive definite matrix Q.
Based on Assumptions 4.1 and 5.1-5.3, we provide the
following main result.
Theorem 5.1: Consider the Lagrangian System (1) and a
smooth, bounded desired trajectory (qd (t), q̇d (t)). Suppose
that Assumptions 4.1, 5.1, 5.2, and 5.3 hold. Then, the
proposed, robust, learning-based control strategy in (4), (11),
and (14), with the uncertainty upper bound ρ calculated by
(16) and the design parameter  satisfying (17), ensures with
high probability of at least (1 − δp )N that the tracking error
e(t) is uniformly ultimately bounded with an ultimate bound
that can be made arbitrarily small through the selection of
the design parameter .
Proof: From Assumption 5.3, we know that kη(t)k ≤
ρ̄ when e(t) ∈ D, where D is a compact set containing
{e ∈ R2N : eT P e ≤ e(0)T P e(0)}. In the first part of the
proof, we assume that the upper bound ρGP (t) calculated by
(9), (10) and (15) is a correct upper bound on kη(t)k when
e(t) ∈ D. Thus, in the first part of the proof, we know that
ρ(t) calculated by (16) is a correct upper bound on kη(t)k
when e(t) ∈ D, and we use Lyapunov stability analysis to
prove that e(t) is uniformly ultimately bounded. Then, in
the second part of the proof, we use Lemma 4.1 to evaluate
the probability of satisfying the assumption that ρGP (t) is a
correct upper bound on kη(t)k when e(t) ∈ D, and hence,
the probability that the provided guarantees hold.
The first part of the proof closely follows the proof of
the effectiveness of the robust controller in Theorem 3 of
Chapter 8 of [1], and we include the main steps of the
proof here for convenience. Consider a candidate Lyapunov
function V (e) = eT P e. From (12), it can be shown that
V̇ = −eT Qe + 2wT (η + r), where w = B T P e. Then, from
(14), we need to study two cases.
For the case where kwk > , we have wT (η +r)=wT (η −
w
)≤−ρkwk + kηkkwk from the Cauchy-Schwartz inρ kwk
equality. Since {e ∈ R2N : eT P e ≤ e(0)T P e(0)} ⊂ D by
definition and from Assumption 5.3, we know that kηk ≤ ρ̄.
Also, by our assumption in this part of the proof, kηk ≤ ρGP .
Then, from (16), kηk ≤ ρ, and wT (η + r) ≤ 0. Thus, for

this case, V̇ ≤ −eT Qe, which ensures exponential decrease
of the Lyapunov function.
Next, consider the case where kwk ≤ . If w = 0, then
V̇ = −eT Qe < 0. Then, for kwk ≤  and w 6= 0, it is easy to
w
+r).
show V̇ = −eT Qe+2wT (η +r)≤−eT Qe+2wT (ρ kwk
w
T
T
From (14), we have V̇ ≤ −e Qe + 2w (ρ kwk − ρ w ).
w
It can be shown that the term 2wT (ρ kwk
− ρ w ) has a
maximum value of (ρ)/2 when kwk = /2. Thus, V̇ ≤
−eT Qe + (ρ)/2. From (16), ρ ≤ ρ̄, and consequently V̇ ≤
−eT Qe+(ρ̄)/2. If the condition eT Qe > (ρ̄)/2 is satisfied,
then V̇ < 0. Since Q is positive definite by definition, then
eT Qe ≥ λmin (Q)kek2 , where λmin (Q) > 0 is the smallest
eigenvalue of Q. Hence, if λmin (Q)kek2 > (ρ̄)/2, then
V̇ < 0.q
Thus, the Lyapunov function is strictly decreasing if
ρ̄
. Let Bδ be the ball around the origin of
kek > 2λmin
q (Q)
ρ̄
radius δ := 2λmin
(Q) , Sδ be a sufficiently small sublevel
set of the Lyapunov function V satisfying B̄δ ⊂ Sδ◦ , and Bc
be the smallest ball around the origin satisfying Sδ ⊂ B̄c .
Since the Lyapunov function V is strictly decreasing outside
B̄δ , the tracking error e(t) eventually reaches and remains
in Sδ ⊂ B̄c , and so the tracking error e(t) is uniformly
ultimately bounded, and its ultimate bound is the radius
of Bc . Note that from (17), Bδ ⊂ {e ∈ R2N : eT P e <
e(0)T P e(0)} ⊂ D, and ρ is a correct upper bound on kηk.
One can see that δ and hence the radius of Bc depend on the
choice of the design parameter . Indeed,  can be selected
sufficiently small to make Bδ and Bc arbitrarily small.
In the second part of the proof, we calculate the probability
of our assumption in the first part that ρGP (t) is a correct
upper bound on kη(t)k when e(t) ∈ D. Recall that e(t) ∈
D implies that aaug (t) is in a compact set, as discussed
immediately after (16). From Assumption 5.2, our problem
reduces to calculating the probability that ρGP is a correct
upper bound on kηk for all the sampling instants. Using the
confidence region proposed in Lemma 4.1 for calculating
the upper bound on the absolute value of each element of
η, and under Assumption 4.1, the probability that this upper
bound is correct for all samples is higher than (1 − δp ) from
Lemma 4.1. Since the N GPs are independent and the added
noise to the output observations ω is uncorrelated, then the
probability that the upper bounds on the absolute values of
all the elments of η, and hence the upper bound on kη(t)k,
are correct is higher than (1 − δp )N .
Remark 5.1: Although in practice it is difficult to estimate
the upper bound ρ̄ on kηk used in (16), one can be conservative in this choice. Then, unlike robust control techniques that
keep this conservative bound unchanged, (16) would relax
the upper bound ρ̄ when the GPs learn a less-conservative
(lower) upper bound from collected data. Having a lessconservative upper bound ρ will result in a lower tracking
error. Indeed, it can be shown that if ρ(t) ≤ ρ0 < ρ̄ for all t,
then the tracking error will converge to an ultimate ball Bc0
smaller than Bc .
Remark 5.2: While in theory  can be selected sufficiently
small to ensure arbitrarily accurate tracking as shown in the

proof of Theorem 5.1, in practice, having a smaller  will
result in a higher r(t) from (14), a higher aq (t), and a higher
control action u. Thus, achieving arbitrarily accurate tracking
may be limited by the actuation limits of the robots, which
is a typical scenario in practice. Incorporating the actuation
limits in the theoretical analysis is an interesting point for
future research.
Remark 5.3: One advantage of the proposed approach is
that it is generalizable in two different ways. First, since the
GPs in our approach learn the error between the commanded
and actual acceleration, and not the error in the force/torque,
and since the GPs do not directly alter the onboard controller,
our approach can be also applied to a wide class of industrial
manipulators that have onboard controllers computing the
joint torques and allow the user to send the commanded joint
acceleration/velocity. Second, the idea of our approach can
be easily generalized to any class of Lagrangian systems
that can be expressed, possibly after a feedforward/feedback
linearization, as a set of decoupled double integrators with
an added uncertainty vector. The GPs can be used to learn
the uncertainty vector resulting from the linearization. For
instance, as shown in Section 5.3 of [18], the nonlinear
dynamics of a quadrotor moving in the x, y plane at a
fixed height can be converted using feedforward linearization
into two decoupled double integrators. Thus, we can again
assume q̈(t) = aq (t) + η(x(t), aq (t)), where η represents the
uncertainty in the linearization, and x is the quadrotor state
vector. GPs can be used to learn η, and an upper bound on
kηk can be then used to design a robust, outer-loop controller
for calculating aq . The details are beyond the scope of this
paper.
VI. S IMULATION R ESULTS
In this section, we present several simulation results to
demonstrate the effectiveness of the proposed approach.
In particular, simulations of the proposed approach were
conducted on a 2 DOF, planar manipulator using MATLAB
Simulink as well as a 6 DOF, industrial manipulator using
the standard Gazebo simulator. It should be noted that in
the simulations, we follow a practical approach in which we
always use a confidence interval [µk − 3σk , µk + 3σk ], which
is simpler to implement and found to be effective in practice
[15]. We also set the upper bound ρ̄ in (16) to be a very
high positive number. Hence, the simulation results mainly
evaluate the effectiveness of the upper bound estimated by
the GPs.
A. 2 DOF Planar Manipulator
In the 2 DOF simulations, we use the robot dynamics
(1) for the system, where M , C, and g are as defined in
Chapter 7 of [1]. For the system parameters, a value of 1 kg
is used for each link mass and 1 kg · m2 for each link inertia.
The length of the first link is 2 m and that of the second
link is 1 m. The joints are assumed to have no mass and
are not affected by friction. Then, it is assumed that these
parameters are not perfectly known. Thus, in the inverse
dynamics controller (4), we use parameters with different

TABLE I
AVERAGE RMS T RACKING E RROR ( IN RAD) OVER 12 T RAJECTORIES
FOR D IFFERENT C ONTROLLERS ON A 2 DOF M ANIPULATOR

Uncertainty
10%
20%
30%

Nominal
0.1557
0.2807
0.3793

Fixed Robust
0.0489
0.0512
0.0534

Robust Learning
0.0065
0.0101
0.0148

levels of uncertainties, and for each uncertainty case, evaluate
our proposed control strategy. The desired trajectories are
sinusoidal trajectories with different amplitudes and frequencies. All the simulation runs are initialized at zero initial
conditions.
Since the linearization by (4) is not perfect due to uncertainties, we use 2 GPs to learn the uncertainty vector
η in (5). Each GP uses the squared exponential kernel,
and is parameterized with ση,i = 1, σω,i = 0.001, and
lj,i = 0.5, for all j ∈ {1, · · · , 6} and i ∈ {1, 2}. The
Gaussian processes use the past n = 20 observations for
prediction. For the robust controller,  is set to 0.001.
A sequence of 12 trajectories was run for 3 different
cases of model uncertainty. Each of the three cases makes
the M̂ matrix differ from the M matrix by using values
for the estimated link masses that differ from the true link
mass values. In particular, in the three uncertainty cases, the
estimated mass differs from the actual mass by 10%, 20%,
and 30% for each link, respectively.
The tracking performance was compared between three
controllers: a nominal controller with no robust control, a
robust controller with a fixed upper bound on the uncertainty
norm ρ = 1000, and our proposed, robust, learning-based
controller. The root-mean-square (RMS) error of the joint
angles was averaged over the 12 trajectories, and is presented
for each controller and uncertainty case in Table I.
It is clear that while the robust controller with a high, fixed
value for the upper bound on the uncertainty improves the
tracking performance compared to the nominal controller, it
is conservative, and thus, still causes considerable tracking
errors. The tracking errors are significantly reduced by our
proposed, robust, learning-based controller, which is able to
learn a less conservative upper bound on the uncertainty
norm. On average, our proposed controller reduces the tracking errors by 96.21% compared to the nominal controller and
by 79.84% compared to the fixed, robust controller.
A comparison of the tracking error for each controller
over a single desired trajectory is shown in Figure 2. The
brief period at the beginning of the trajectory shows the
learning period of the robust learning controller. The fixed
robust controller outperforms the robust learning controller
until the latter’s Gaussian processes have observed enough
data to make more accurate predictions.
B. 6 DOF Industrial Manipulator
The proposed approach is further tested on a 6 DOF
industrial manipulator using the Robot Operating System
(ROS) and Gazebo simulator. The model parameters are

Fig. 2. The joint angle error over time for each controller for tracking
π
π
a single trajectory qd (t) = ( 18
sin(0.4t), 18
sin(0.4t)) with 30% model
uncertainty in the link masses.
TABLE II
AVERAGE RMS T RACKING E RROR (IN RAD) OVER 6 T RAJECTORIES
ON A 6 DOF M ANIPULATOR

Uncertainty
30%
60%
90%
Average

Nominal
0.0275
0.0434
0.0539
0.0416

Robust Learning
0.0139
0.0218
0.0304
0.0221

Improvement
58.05%
49.66%
43.52%
50.41%

selected based on the UR10 robot manipulator’s open-source
Gazebo package in [21]. The saturation limits are set to
sufficiently high values to evaluate the system performance
when the desired trajectories are feasible to be perfectly
tracked under the actuation limits.
Similar to our discussion in the previous subsection, we
use in this example 6 GPs to learn the uncertainty vector η.
Each GP uses the squared exponential kernel, and is parameterized with ση,i = 8.0, σω,i = 0.001, and lj,i = 0.05,
for all j ∈ {1, · · · , 18} and i ∈ {1, · · · , 6}. Also, each GP
uses the past n = 20 observations for prediction, and runs
at a sampling rate of 125 Hz. For the robust controller, the
parameter  is set to 5 × 10−6 .
The performance of the robust learning controller is compared with that of the nominal, non-robust controller on 6
different trajectories and 3 different, large model uncertainties. Model uncertainties were introduced in the same way
as in the 2 DOF robot simulation; that is, the estimated
link masses used in the inverse dynamics controller (4)
were changed from the true values used in the simulation’s
physics environment. The three different uncertainty cases
use estimated link masses that differ from the true link
masses by 30%, 60%, and 90%, respectively. The results
of the comparison are summarized in Table II.
It is clear that our proposed controller consistently reduces
the tracking error in the presence of model uncertainties.
The performance of each controller over all six trajectories
with a fixed model uncertainty in the link masses of 60% is
presented in Table III. This table shows how our proposed
controller consistently provides a significant performance
improvement over the nominal controller across a variety
of trajectories.

TABLE III
RMS T RACKING E RROR (IN RAD) FOR 60% M ODEL U NCERTAINTY
AND DIFFERENT T RAJECTORIES

Trajectory
#1
#2
#3
#4
#5
#6
Average

Nominal
0.0389
0.0416
0.0450
0.0466
0.0444
0.0437
0.0434

Robust Learning
0.0204
0.0210
0.0211
0.0235
0.0223
0.0225
0.0218

Improvement
47.50%
49.54%
53.10%
49.49%
49.77%
48.58%
49.66%

Fig. 3. The end effector position path for Trajectory #1, #2, #3, and #4
in Table III, respectively, under 60% model uncertainty in the masses and
our proposed robust, learning-based controller. The black dots represent the
initial positions of the end effector.

The 3-dimensional paths followed by the manipulator’s
end effector under 60% model uncertainty and our proposed
controller for 4 different desired trajectories are shown in
Figure 3. As can be seen from the figure, despite the large
uncertainty, the proposed controller makes the end effector
follow the desired path very closely.
VII. C ONCLUSIONS
For an important class of Lagrangian systems, we have
provided a novel, learning-based control strategy based on
Gaussian processes (GPs) that ensures stability of the closedloop system and high-accuracy tracking of smooth trajectories. The main idea is to use GPs to estimate an upper
bound on the uncertainty of the linearized system, and then
use the uncertainty bound in a robust, outer-loop controller.
Unlike most of the existing, learning-based inverse dynamics
control techniques, we have provided a proof of the closedloop stability of the system that takes into consideration
the regression errors of the learning modules. Moreover,
we have proved that the tracking error converges to a ball
with a radius that can be made arbitrarily small through
appropriate control design. Furthermore, we have verified
the effectiveness of our approach via simulations on planar
manipulators and 6 DOF industrial manipulators using the

standard Gazebo simulator. A promising direction for future
research is to formally incorporate the saturation limits of the
system into the theoretical analysis of the proposed control
system behavior.
R EFERENCES
[1] M. W. Spong, S. Hutchinson, and M. Vidyasagar. Robot Modeling and
Control. John Wiley & Sons, Inc., 2006.
[2] V. Hagenmeyer, and E. Delaleau. Exact feedforward linearization
based on differential flatness. International Journal of Control, 76(6),
pp. 537-556, 2003.
[3] C. Abdallah, D. M. Dawson, P. Dorato, and M. Jamshidi. Survey of
robust control for rigid robots. IEEE Control Systems Magazine, 11(2),
pp. 24-30, 1991.
[4] J. Sun de la Cruz. Learning Inverse Dynamics for Robot Manipulator
Control. M.A.Sc. Thesis, University of Waterloo, 2011.
[5] R. Calandra, S. Ivaldi, M. P. Deisenroth, E. Rueckert, and J. Peters.
Learning inverse dynamics models with contacts. IEEE International
Conference on Robotics and Automation, Seattle, 2015, pp. 3186-3191.
[6] A. S. Polydoros, L. Nalpantidis, and V. Kruger. Real-time deep learning of robotic manipulator inverse dynamics. IEEE/RSJ International
Conference on Intelligent Robots and Systems, Hamburg, 2015, pp.
3442-3448.
[7] S. Zhou, M. K. Helwa, and A. P. Schoellig. Design of deep neural
networks as add-on blocks for improving impromptu trajectory tracking. IEEE Conference on Decision and Control, Melbourne, 2017, pp.
5201-5207.
[8] D. Nguyen-Tuong, J. Peters. Using model knowledge for learning
inverse dynamics. IEEE International Conference on Robotics and
Automation, Anchorage, 2010, pp. 2677-2682.
[9] J. F. Fisac, A. K. Akametalu, M. N. Zeilinger, S. Kaynama, J. Gillula,
and C. J. Tomlin. A general safety framework for learning-based
control in uncertain robotic systems. arXiv preprint arXiv:1705.01292,
2017.
[10] T. J. Perkins, and A. G. Barto. Lyapunov design for safe reinforcement
learning. The Journal of Machine Learning Research, 3, pp. 803-832,
2003.
[11] J. W. Roberts, I. R. Manchester, and R. Tedrake. Feedback controller
parameterizations for reinforcement learning. IEEE Symposium on
Adaptive Dynamic Programming and Reinforcement Learning, Paris,
2011, pp. 310-317.
[12] P. Geibel, and F. Wysotzki. Risk-sensitive reinforcement learning
applied to control under constraints. Journal of Artificial Intelligence
Research, 24, pp. 81-108, 2005.
[13] M. Turchetta, F. Berkenkamp, and A. Krause. Safe exploration in
finite markov decision processes with Gaussian processes. Advances
in Neural Information Processing Systems, Barcelona, 2016, pp. 43124320.
[14] F. Berkenkamp, R. Moriconi, A. P. Schoellig, A. Krause. Safe learning
of regions of attraction for uncertain, nonlinear systems with Gaussian
processes. IEEE Conference on Decision and Control, Las Vegas,
2016, pp. 4661-4666.
[15] F. Berkenkamp, A. P. Schoellig, and A. Krause. Safe controller optimization for quadrotors with Gaussian processes. IEEE International
Conference on Robotics and Automation, Stockholm, 2016, pp. 491496.
[16] F. Meier, D. Kappler, N. Ratliff, and S. Schaal. Towards robust online
inverse dynamics learning. IEEE/RSJ International Conference on
Intelligent Robots and Systems, Daejeon, 2016, pp. 4034-4039.
[17] T. Beckers, J. Umlauft, D. Kulić, and S. Hirche. Stable Gaussian process based tracking control of Lagrangian systems. IEEE Conference
on Decision and Control, Melbourne, 2017, pp. 5580-5585.
[18] M. K. Helwa, A. P. Schoellig. On the construction of safe controllable
regions for affine systems with applications to robotics. arXiv preprint
arXiv:1610.01243, 2016.
[19] R. M. Murray. Nonlinear control of mechanical systems: a Lagrangian
perspective. Annual Reviews in Control, 21, pp. 31-42, 1997.
[20] N. Srinivas, A. Krause, S. M. Kakade, and M. W. Seeger. Informationtheoretic regret bounds for Gaussian process optimization in the bandit
setting. IEEE Transactions on Information Theory, 58(5), pp. 32503265, 2012.
[21] universal robot package, ver. 1.1.9, link: wiki.ros.org/universal robot,
2014.

