Safe Learning of Regions of Attraction for Uncertain, Nonlinear
Systems with Gaussian Processes

arXiv:1603.04915v3 [] 16 Aug 2017

Felix Berkenkamp, Riccardo Moriconi, Angela P. Schoellig, and Andreas Krause
Abstract— Control theory can provide useful insights into
the properties of controlled, dynamic systems. One important
property of nonlinear systems is the region of attraction (ROA),
a safe subset of the state space in which a given controller
renders an equilibrium point asymptotically stable. The ROA
is typically estimated based on a model of the system. However,
since models are only an approximation of the real world, the
resulting estimated safe region can contain states outside the
ROA of the real system. This is not acceptable in safety-critical
applications. In this paper, we consider an approach that learns
the ROA from experiments on a real system, without ever
leaving the true ROA and, thus, without risking safety-critical
failures. Based on regularity assumptions on the model errors
in terms of a Gaussian process prior, we use an underlying
Lyapunov function in order to determine a region in which an
equilibrium point is asymptotically stable with high probability.
Moreover, we provide an algorithm to actively and safely
explore the state space in order to expand the ROA estimate.
We demonstrate the effectiveness of this method in simulation.

I. I NTRODUCTION
Estimating the region of attraction (ROA) of an equilibrium point is an important problem when analyzing nonlinear
systems. It specifies the region within which a given control
law renders an equilibrium point asymptotically stable. Thus,
the system can be operated safely within the ROA. This is
of practical importance, since controllers for safety-critical
systems are required to guarantee safety over a certain
domain of operation, before they can be implemented on
the real system. Typically, ROAs are estimated based on a
model of the system. However, due to model errors, ROAs of
the real system can be drastically different, raising questions
about the viability of such purely model-based methods.
In this paper, we address the problem of safely estimating
the ROA from experiments on a real system, while always
remaining in the real ROA with high probability. As a result,
the equilibrium point of the closed-loop system is attractive
for any state visited throughout the learning process. This
avoids failures and makes the method applicable to safetycritical domains. To expand the estimated safe region over
time, we actively select states that improve our estimate of
the underlying uncertain model.
Felix Berkenkamp and Andreas Krause are with the Learning & Adaptive
Systems Group (www.las.ethz.ch), Department of Computer Science, ETH
Zurich, Switzerland. Email: {befelix, krausea}@ethz.ch
Riccardo Moriconi is with the Department of Mechanical Engineering,
ETH Zurich, Switzerland. Email: rmoricon@ethz.ch
Angela P. Schoellig is with the Dynamic Systems Lab (www.dynsyslab.
org), Institute for Aerospace Studies, University of Toronto, Canada. Email:
schoellig@utias.utoronto.ca
This research was supported by SNSF grant 200020 159557, NSERC
grant RGPIN-2014-04634, and the Connaught New Researcher Award.

Related work: In the literature, ROAs are computed based
on system models [1]. One typical approach is to use the level
sets of a Lyapunov function in order to quantify the ROA [2].
Given a Lyapunov function, the problem of finding the ROA
reduces to a line search for the maximum safe level set of
the Lyapunov function. For polynomial systems, Lyapunov
functions can be found efficiently by solving a system
of linear matrix inequalities (LMIs) [3]. A relaxation to
Lyapunov-like functions for ROA computation is given in [4]
and a review of numerical methods to compute Lyapunov
functions can be found in [5]. The approach in [6] considered
quadratic systems and the problem of testing whether a given
polytope lies inside the ROA. Beyond Lyapunov functions,
the ROA can be found via constraint solving [7] or by
optimizing over state trajectories [8]. The ROA can also be
approximated using sampling [9].
An area that is of particular relevance to this paper is
the computation of ROAs for systems with uncertainties.
In [10], quadratic Lyapunov functions for uncertain, linear
systems were used, while [11] considered polynomial Lyapunov functions. A more general approach for systems with a
polynomial dependence on uncertainties within a polytope is
shown in [12]. Similar ideas are used in robust control [13].
All the previous methods for estimating ROAs are based
on a model with fixed uncertainties. The method proposed
in this paper uses learning in order to reduce the model
uncertainty over time. Learning the dynamics of nonlinear
systems from data has been considered before. In particular,
Gaussian processes (GPs [14]) have been of interest in this
area, since they provide uncertainty information about the
estimated dynamic model. This can be used to derive highprobability safety guarantees. For example, [15] combines
GP regression with robust control theory to provide stability
and performance guarantees for a linear controller applied to
a nonlinear system, while [16] uses reachability analysis to
compute robust invariant sets.
In our method, we use ideas from Bayesian optimization [17], set of optimization methods that aim to find the
global optimum of an a priori unknown function based on
noisy evaluations. Bayesian optimization algorithms that are
based on GP models provably converge close to the global
optimum [18]. In this setting, [19], [20] consider safety
in terms of high-probability constraint satisfaction, with a
generalization to multiple constraints in [21].
Our contribution: In this paper, we combine ideas from GP
learning, safe Bayesian optimization, and ROA computation
based on Lyapunov functions for uncertain systems. Using
previous results on model-based ROA computation as a

Appeared in Proc. of the IEEE Conference on Decision and Control, 2016, pp. 4661 – 4666, doi: 10.1109/CDC.2016.7798979 .
c 2016 IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any current or future media, including
reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution to servers or lists, or reuse of
any copyrighted component of this work in other works.

starting point, we compute ROA estimates for the real system
by approximating the model uncertainties with a GP on a
discrete domain. We use ideas from safe Bayesian optimization in order to actively learn about the real dynamics from
experiments executed in the estimated ROA. As we learn
more about the dynamics of the real system, the uncertainty
in the estimate decreases and the ROA expands, until we
reach a desired accuracy. As a result, we are able to learn an
accurate estimate of the ROA of a general, nonlinear system
through experiments, without leaving the ROA and, thus,
without violating safety requirements.
II. P ROBLEM S TATEMENT
We consider a nonlinear, continuous-time system,
ẋ(t) = f (x(t), u(t)) + g(x(t), u(t)),
{z
} |
{z
}
|
a priori model

(1)

unknown model

where x(t) ∈ X ⊆ R is the state at time t within a connected set X and u(t) ∈ U ⊆ Rp is the control input. The
system dynamics consist of a known model f (x, u) and an
unknown model g(x, u). The latter accounts for unknown
dynamics and model uncertainties. We assume a control
policy u = π(x) is given, which has been designed for the
prior model, f (·). The resulting, closed-loop dynamics are
denoted by fπ (x) := f (x, π(x)) and gπ (x) := g(x, π(x)).
The goal is to estimate the ROA of (1) under the control
policy, π(x), based on experiments on the real system. We
want to actively learn the a priori unknown dynamics, gπ ,
without leaving the ROA. This is, of course, impossible
without further assumptions about the model in (1).
We assume that the unknown model, gπ , has low ‘complexity’, as measured under the norm of a reproducing kernel
Hilbert space (RKHS, [22]). An RKHS Hk (X ) is a complete
subspace of L2 (X ), the space of square integrable
P functions,
that includes functions of the form gπ (x) = i αi k(x, xi )
with αi ∈ R and representer points xi ∈ X , where k
(and the subscript k) refer to a symmetric, positivedefinite kernel function k(·, ·). The RKHS has an inner product h·, ·ik , which obeys the reproducing property: hgπ (·), k(x, ·)ik = gπ (x) for all gπ ∈ Hk (X ). The induced RKHS norm kgπ k2k = hgπ , gπ ik measures smoothness
of gπ with respect to the kernel k(·, ·), [18]. For universal
kernels [22], members of Hk (X ) can uniformly approximate
any continuous function on any compact subset of X . We
make assumptions about the complexity of gπ .
q

Assumption 1. The function gπ has bounded RKHS norm
with respect to a continuously differentiable, bounded kernel k(x, x0 ); that is, kgπ kk ≤ Bg .
Moreover, we make the following, standard assumption on
the a priori model, fπ (x):

Assumption 2. The function fπ is Lipschitz continuous with
Lipschitz constant Lf and bounded in X with kfπ k∞ ≤ Bf .
In the following sections, we use a Lyapunov function in
order to compute the ROA. In practice, one may compute
good Lyapunov functions during the experimentation phase

using, for example, the methods in [5]. Here, we consider a
fixed, given Lyapunov function. We assume the following:
Assumption 3. The origin is an equilibrium point of (1)
with fπ (0) = gπ (0) = 0.
Assumption 4. A fixed, two-times continuously differentiable
Lyapunov function V (x) is given. Moreover, there exists
a constant c0 > 0 such that ∂V∂x(x) (fπ (x) + gπ (x)) < 0 for
all x ∈ S0 = V(c0 ), where V(c) = {x ∈ X | V (x) ≤ c}.

Assumption 3 implies that the origin is an equilibrium
point of (1) with the unknown model gπ (x). Without this
assumption, there is no hope to prove that the origin is
attractive, but showing boundedness would be possible. Assumption 4 implies that we have chosen the policy, π(x), and
the Lyapunov function, V (x), such that the origin of (1) is
locally asymptotically stable. This is, for example, fulfilled
if the origin of the a priori dynamics ẋ = fπ (x) is locally
π (x)
is zero at the origin; that
asymptotically stable and ∂g∂x
is, the prior dynamics dominate close to the origin.
Lastly, we require measurements of gπ (x) in order to
learn about the unknown dynamics. That is, we must be
able to measure x and ẋ (fπ (·) is known). In practice, we
may consider a discrete-time approximation of ẋ, which
only requires measurements of the state, x, rather than its
derivative, ẋ. We discuss practical aspects in Sec. VII.
Assumption 5. We have access to measurements,
ĝπ (x) = ẋ − fπ (x) + ω, which are corrupted by zeromean, independent and bounded noise, kωk ≤ σ.
For ease of notation, we consider a one-dimensional
system of the form of (1) in the following. We explain how
to generalize the results to multiple dimensions in Sec. VI.

III. G AUSSIAN P ROCESSES (GP S )
The function gπ (x) in (1) is unknown a priori. Given Assumption 1, we use a GP to approximate this unknown
function over its domain X , see [18] and Lemma 1 below.
The following introduction to GPs is based on [14].
GPs are a nonparametric regression method from machine
learning, where the goal is to find an approximation of a
nonlinear map gπ : X → R from a state x to the function
value gπ (x). This is accomplished by considering the function values gπ (x) to be random variables, so that any finite
number of them have a joint Gaussian distribution [14].
The Bayesian, nonparametric regression is based on a prior
mean function and a covariance function, k(x, x0 ), which
defines the covariance of any two function values, gπ (x)
and gπ (x0 ), x, x0 ∈ X . The latter is also known as the
kernel. In this work, the mean is zero, since any prior
knowledge about the dynamics is captured by fπ (x); that
is, gπ (x) ∼ GP(0, k(x, x0 )). In general, the choice of kernel
function is problem-dependent and encodes assumptions
about the unknown function. A discussion of kernel choices
for dynamic systems can be found in [23]. In the following,
we use the same kernel as in Assumption 1.
We can obtain the posterior distribution of a function
value gπ (x) at an arbitrary state x ∈ X by conditioning the

GP distribution of gπ on a set of n past measurements,
yn = (ĝπ (x1 ), . . . , ĝπ (xn )) at states Dn = {x1 , . . . , xn },
where ĝπ (x) = ẋ − fπ (x) + ω with ω ∼ N (0, σ 2 ). The
posterior over gπ (x) is a GP distribution again, with
mean µn (x), covariance kn (x, x0 ), and variance σn (x):
µn (x) = kn (x)(Kn + In σ 2 )−1 yn ,
0

0

kn (x, x ) = k(x, x ) − kn (x)(Kn +

In σ 2 )−1 knT (x0 ),

(2)
(3)

(4)

where the vector kn (x) = k(x, x1 ), . . . , k(x, xn ) contains
the covariances between the new input, x, and the states
in Dn , Kn ∈ Rn×n is the positive-definite kernel matrix with
[Kn ](i,j) = k(xi , xj ), i, j ∈ {1, . . . , n}, and In ∈ Rn×n is
the identity matrix.
Assumption 1 allows us to model gπ (x) as a GP. In
particular, we have the following result:
σn2 (x)

= kn (x, x),

Lemma 1. Supposed that kgπ k2k ≤ Bg and that
the zero-mean noise ω is uniformly bounded by σ.
Choose βn = 2Bg + 300γn log3 (n/δ), where γn is defined
in the paragraph below. Then, for all n ≥ 1 and x ∈ X , it
holds with probability at least 1 − δ, δ ∈ (0, 1), that
|gπ (x) − µn−1 (x)| ≤ βn1/2 σn−1 (x).

(5)

Proof. See [18, Theorem 6].
Lemma 1 allows us to make high-probability statements
about the function values of gπ (x), even though the GP
model makes different assumptions than we do in Assumption 1 (e.g., about the noise, ω). The scalar βn depends on the
information capacity, γn = maxx1 ,...,xn I(gπ , yn ), which is
the maximal mutual information that can be obtained about
the GP prior from n noisy samples yn at states Dn . In [18],
it was shown that γn has a sublinear dependence on n for
many commonly used kernels and that it can be efficiently
approximated up to a constant. As a result, even though βn is
increasing with n, we are able to learn about the true values
of gπ (x) over time by making appropriate choices for data
samples in Dn , see Sec. V.

Since the closed-loop dynamics in (1) are Lipschitz continuous based on Lemma 2 and Assumption 2, existence and
uniqueness of the solutions of (1) are guaranteed.
The goal of this section is to quantify the ROA based
on the Lyapunov function in Assumption 4. From Lyapunov
stability theory we have the following [2]:
Lemma 3. The origin of the dynamics in (1) is asymptotically stable within a level set, V(c) = {x ∈ X | V (x) ≤ c}
with c ∈ R>0 , if, for all x ∈ V(c),
V̇ (x) =


∂V (x)
fπ (x) + gπ (x) < 0.
∂x

(6)

The goal of this section is to use the GP model of gπ (x)
from Sec. III in order to find the largest c such that (6) holds
true within V(c) with high probability. The existence of such
a value c0 > 0 is ensured by Assumption 4.
In order to quantify the ROA using (6), we need to
evaluate V̇ (x), which depends on the GP model of gπ (x).
From Sec. III we know that the GP posterior distribution
over gπ is Gaussian. As a consequence, since V̇ (x) in (6) is
affine in gπ (x), V̇ (x) itself is a GP. In particular, we have
that V̇ (x) is normally distributed with mean µn,V̇ (x) and
standard deviation σn,V̇ (x) given by:

∂V (x)
µn (x) + fπ (x) ,
∂x
∂V (x)
σn,V̇ (x) =
σn (x),
∂x

µn,V̇ (x) =

(7)
(8)

where µn (x) and σn (x) are the GP predictions of the
unknown dynamics, gπ (x), from (2) and (4). We can use (7)
and (8) to directly make predictions about V̇ (x) in (6) based
on measurements of gπ (x).
Since (7) and (8) provide a probabilistic estimate of V̇ (x),
we cannot expect to make deterministic statements about
stability. Instead, we consider the confidence intervals of
the GP model of V̇ (x). We denote the lower and upper
confidence intervals after (n − 1) measurements of gπ (x) by
ln (x) := µV̇ ,n−1 (x) − βn1/2 σV̇ ,n−1 ,

(9)

IV. LYAPUNOV S TABILITY

un (x) := µV̇ ,n−1 (x) + βn1/2 σV̇ ,n−1 ,

In this section, we show how the assumptions in Sec. II
can be used in order to compute the ROA for the nonlinear
system in (1) based on a GP model of gπ from Sec. III. For
now, we base our analysis on a general GP model of gπ (x)
from Sec. III with n arbitrary measurements of gπ . We
actively select states for new measurements in Sec. V.
We start by observing the following:

respectively. The confidence intervals are parameterized by
the scalar βn . In the following, we assume that βn is chosen
according to Lemma 1, which allows us to say that V̇ (x)
takes values within the interval [ln (x), un (x)] with high
probability (at least 1 − δ).
Based on these confidence intervals, we can see from (6)
that V̇ (x) < 0 holds within V(c) for some c > 0 with high
probability, if un (x) < 0 for all x ∈ V(c). Unfortunately, it
is impossible to evaluate the upper bound (10) everywhere in
a continuous domain. Nevertheless, it is possible to evaluate
predictions of V̇ (x) from (7) and (8) at a finite number of
points. In the following, we exploit the continuity properties
of V̇ (x) in order to derive stability properties within a
continuous domain based on a finite number of predictions.
In particular, we use the following property of V̇ :

Lemma 2. The function gπ is Lipschitz continuous with
Lipschitz constant Lg and bounded by kgπ k∞ ≤ Bg kkk∞ .

Proof. Boundedness by [22, Lemma 4.23]. From Assumption 1, |gπ (x) − gπ (x0 )|2 ≤ kgπ k2k dk (x, x0 ), where
dk (x, x0 ) = k(x, x) + k(x0 , x0 ) − 2k(x, x0 ) is the kernel
metric [22, Lemma 4.28]. Since k(·, ·) is continuously differentiable and bounded, |gπ (x) − gπ (x0 )|2 ≤ L2g |x − x0 |2 ,
∂k
where L2g = 2Bg kkk∞ k ∂x
k∞ .

(10)

Lemma 4. The function V̇ (x) is Lipschitz continuous with
Lipschitz constant L.
Proof. Based on the Lipschitz continuity of (1) from Assumption 2 and Lemma 2, we expand |V̇ (x) − V̇ (x0 )| to
∂V (x)
∂V (x0 )
(fπ (x) + gπ (x)) −
(fπ (x0 ) + gπ (x0 ))
∂x
∂x
∂V (x) ∂V (x0 )
≤ fπ (x) + gπ (x)
−
∂x
∂x0
0
∂V (x )
+
fπ (x) + gπ (x) − fπ (x0 ) − gπ (x0 ) ,
∂x0
≤ (Bf + Bg kkk∞ )L∂V |x − x0 | + LV (Lg + Lf )|x − x0 |,
:= L|x − x0 |,

2

V (x)
where LV = k ∂V∂x(x) k∞ and L∂V = k ∂ ∂x
k∞ are the Lip2
schitz constants of V and its first derivative. These are
guaranteed to exist by Assumption 4, since a continuous
function obtains a maximum over a bounded domain.

The continuity of V̇ (x) allows us to evaluate predictions
of V̇ (x) from (7) and (8) only at a finite number of points,
without loosing guarantees.
Lemma 5. Let Xτ ⊂ X be a discretization of X with
|x − [x]τ | ≤ τ /2 for all x ∈ X . Here, [x]τ denotes the
closest point in Xτ to x ∈ X . Choosing βn according
to Lemma 1, the following holds with probability at
least (1 − δ) for all x ∈ X and all n ≥ 1:
V̇ (x) − µV̇ , n−1 ([x]τ ) ≤ βn1/2 σV̇ , n−1 ([x]τ ) + Lτ. (11)
Proof. The result follows from the Lipschitz continuity
of V̇ (x) in Lemma 4 and is similar to [18, Lemma 5.7].
Lemma 5 provides high-probability bounds on V̇ on the
continuous domain X using the GP confidence intervals (7)
and (8) on the discrete set Xτ . It achieves this by using
the Lipschitz property of V̇ to generalize from Xτ to X .
The discretization accuracy τ trades off the accuracy for
the reduced computational cost of computing the confidence
intervals for all states in Xτ . Combining the previous results
lets us argue about the stability:
Lemma 6. With a discretization of X , Xτ , according
to Lemma 5 and with βn according to Lemma 1, the origin
of (1) is asymptotically stable within V(c) for some c > 0
with probability at least (1 − δ) if, for all x ∈ V(c) ∩ Xτ ,
µV̇ , n−1 (x) + βn1/2 σV̇ , n−1 (x) = un (x) < −Lτ.

(12)

Proof. This is a consequence of Lemmas 1, 3 and 5.
Given the previous results, after (n − 1) evaluations of gπ ,
it suffices to maximize c such that the condition in Lemma 6
holds for all discretized states in V(c):
Theorem 1. Under the assumptions of Lemma 6, pick
cn = argmax c,
c∈R>0

subject to (12) for all x ∈ V(c) ∩ Xτ .

(13)

Algorithm 1: Safe ROA exploration
Inputs: Domain X and discretization with τ , Xτ
GP prior k(x, x0 )
Initial safe set S0 ⊆ X
1 for n = 1, . . . do
cn ← argmax c, subject to

2

3
4
5

c>0

un (x) < −Lτ, for all x ∈ V(c) ∩ Xτ
Sn ← S0 ∪ V(cn )
xn ← argmaxx∈Sn σn−1 (x)
Update GP with measurement of ĝπ (xn )

Then, Sn = S0 ∪ V(cn ) is contained within the ROA of (1)
for all n ≥ 1 with probability at least (1 − δ).

Proof. The set S0 from Assumption 4 is contained in the
ROA deterministically, while V(cn ) is contained with probability at least (1 − δ) by Lemma 6. The result follows.

Theorem 1 provides a way to compute an estimate of
the ROA via an efficient binary search to solve the optimization problem (13). The discretization in Lemma 6
allows us to consider high-probability confidence intervals
on Xτ only, while we generalize using the Lipschitz constant
from Lemma 4. A finer discretization (smaller value of τ )
decreases the conservativeness of the ROA estimate. This is
different from the number of experiments, n, which decreases
the uncertainty in the model estimate as we obtain more
data. In particular, as more information about gπ becomes
available from measurements, the ROA increases beyond the
initial, deterministic safe set in Assumption 4.
V. ACTIVE LEARNING
In the previous section, we provided an estimate of the
ROA based on a GP model of the unknown dynamics
of gπ (x) in (1), which was informed by arbitrarily selected
measurements/experiments. In this section, we actively expand the ROA estimate by selecting new, safe states within
the current ROA estimate at which to obtain measurements.
As we update the model with new data, we can compute
the current estimate of the ROA, Sn , using Theorem 1.
In order to expand the ROA from the initial estimate S0 ,
we need to obtain measurements at states that are relevant
for learning gπ . Here, we assume that we have access to a
control method that drives us to desired states without leaving
the ROA. Any method that applies the policy π(x) on the
boundary of Sn fulfills this safety requirement [16]. Based
on ideas from safe Bayesian optimization [19], we use the
uncertainty estimate of the GP. We select
xn = argmax σn−1 (x)

(14)

x∈Sn

as the next state to evaluate according to (4) within a current
estimate of the safe set, Sn . At xn , the GP model is the
most uncertain about the unknown dynamics, gπ (x). The
idea behind this selection criterion is that we want to learn
about the most uncertain state in order to increase the ROA
estimate. By decreasing the uncertainty about the unknown

fπ (x) + gπ (x)

True dynamics
Estimated dynamics
True V̇ (x)

S20

Estimated V̇ (x)

V̇ (x)

V (x)

−Lτ
S4
S0
State x
a: Initial safe set, S0 .

State x

State x

b: Exploration within safe level set, Sn .

0

c: Maximum level set found.

Fig. 1. One-dimensional example of Algorithm 1. Initially, in Fig. 1a, the estimate of the dynamics is uncertain (top row, blue shade is the 2σ confidence
interval). As a result, the ROA (red part of the Lyapunov function, bottom plot) consists only of the initial safe set, S0 ; larger level sets of the Lyapunov
function are unsafe (blue part). The algorithm actively selects new states at which to obtain measurements of the dynamics, which causes the safe set to
increase in Fig. 1b. Eventually, this procedure leads to the largest safe level set in Fig. 1c. Only states inside the ROA are evaluated during the learning.

dynamics, gπ (x), over time, the safe region expands eventually. We summarize the method in Algorithm 1 and an
example is shown in Fig. 1.
It follows from Theorem 1 that every state chosen by (14)
lies inside the ROA of (1) with high probability. For a
variant of Algorithm 1 it is possible to prove an even
stronger result: that the maximum, safe level set can be
found to some accuracy (with probability at least (1 − δ))
after a finite number of experiments. The corresponding
variant of Algorithm 1 builds up the ROA estimate using the
Lipschitz constant. More details are given in the appendix
in Theorem 2.
VI. E XTENSION TO M ULTIPLE D IMENSIONS
So far, we have considered a one-dimensional system
in (1). The preceding analysis directly transfers to multiple
dimensions, at the expense of more cumbersome notation.
The only non-trivial parts of the extension are GP models
with vector predictions and the choice of βn in Lemma 1.
The main observation is that vector-valued functions can be
modeled as a single GP over an extended state space, X × I,
where integer elements in I index the output dimension
of gπ (x). At each iteration, we obtain q measurements
of gπ (x); one for each dimension of the state. As a re0
sult, βn0 = 2Bg + 300γnq
log3 (nq/δ) increases at a faster
0
rate of n = nq compared to βn in Lemma 1. The information capacity, γ 0 , is measured relative to the combined
function over X × I. Refer to [21] for more details.
VII. D ISCUSSION
Algorithm 1 allows for the safe exploration of the ROA
without leaving the ROA for a general, uncertain, nonlinear
system. While the discretization in Lemma 5 may be conservative and computing the predictions of the GP model of V̇
at the discretization points is computationally intensive, this

needs to be compared to the generality of the statements we
have made. In particular, the only assumption we have made
about the unknown dynamics, gπ (x) in (1), is Assumption 1,
which is very general.
In the analysis, we have assumed that a method exists,
which drives the system to any state xn selected by (14),
without leaving the ROA. In practice, this must be further
restricted by reachability properties, see [16]. These can be
included as an additional constraint; that is, if R is the safely
reachable set, we select xn from Sn ∩ R in (14).
Algorithm 1 can be made more data-efficient by only
evaluating states close to the boundary of the level set,
without loosing guarantees. See [19] and [21] for details.
A discrete-time variant of Algorithm 1 that determines
safety with V (fπ,d (x) + gπ,d (x)) − V (x) < 0 only requires
measurements of the states, x, rather than derivatives, ẋ,
in Assumption 5. While mapping a GP model of gπ (x)
through a nonlinear Lyapunov function renders this method
not analytically tractable, methods based on approximate
uncertainty propagation may work well in practice.
VIII. E XPERIMENTS
In this section. we evaluate Algorithm 1 in a simulated
experiment. We only provide a high-level overview of the
experiment. For details refer to the documentation and source
code at http://github.com/befelix/lyapunov-learning.
We consider an inverted pendulum with angle θ,
mass m = 0.15 kg, length l = 0.5 m, and friction coefficient µ = 0.05 Nms/rad. The dynamics are given by
mgl sin θ(t) − µθ̇(t) + u(t)
,
(15)
ml2
where u(t) is the torque applied to the pendulum. The torque
is limited so that the real system cannot recover from states
with |θ| > 30 deg. The state is x = (θ, θ̇). We assume that
θ̈(t) =

Angular velocity θ̇ [rad/s]

R EFERENCES

4
2
0
Prior safe set
True safe set
Estimated safe set

−2
−4

−0.3

−0.2

−0.1

0.0

0.1

0.2

0.3

Angle θ [rad]

Fig. 2. Prior, true, and estimated ROA level sets after 100 data points
(blues crosses). The prior model estimates a safe region (purple) that is
larger than the true ROA (green), thus it includes unsafe states. In contrast,
Algorithm 1 provides a conservative estimate (yellow), since it considers
states with V̇ (x) ≤ −Lτ , rather than V̇ (x) ≤ 0. The level set could be
increased by discretizing the space with a smaller value of τ .

we only know a linear approximation of the dynamics in (15)
for the upright equilibrium point, where additionally friction
is neglected and the mass is 0.05 kg lighter. We use a Linear
Quadratic Regulator based on this linear approximation in
order to control the origin, x = 0, and use the corresponding
quadratic Lyapunov function to determine the ROA of (15).
We model the error in (15) as a GP, see Sec. III. In
particular, we use the product of a linear and a Matérn kernel
for the GP model, which encodes nonlinear functions that are
two-times differentiable and take values that are bounded by
a linear function from above and below with high probability.
Details about kernel choice for dynamic systems can be
found in [23] and a one-dimensional sample function of the
kernel used here is in the upper plot of Fig. 1a. In practice,
the more assumptions are made about the model error, the
faster the learning algorithm will converge [23].
In the analysis, we assumed gπ to have bounded RKHS
norm in Assumption 1. Here, we model gπ as a sample
1/2
function of the GP and set βn = 2 for all iterations. We
use τ = 0.002 and high-probability, local Lipschitz constants
encoded by the kernel with Lemma 4. The initial safe set
is S0 = {(θ, θ̇) ∈ R2 | |θ| ≤ 5 deg, |θ̇| ≤ 10 deg /s}.
The results can be seen in Fig. 2. While the prior model
with the wrong mass and friction parameters estimates a safe
set that is too large (includes unsafe states), Algorithm 1
provides a conservative estimate. As we gain more data
about the dynamics and if we discretize with smaller values
of τ , the estimate improves and, in the limit, converges to
the true level set. Overall, Algorithm 1 provides a powerful
tool to learn the ROA of general nonlinear systems from
experiments, without leaving the safe region encoded by the
Lyapunov function.
IX. C ONCLUSION
We introduced a learning algorithm that, based on an
initial, approximate model and a corresponding Lyapunov
function, is able to learn about the real ROA through experiments on the real system, without leaving the true ROA with
high probability. While some of the assumptions in Sec. II
may be restrictive for practical application, the results in this
paper should be understood as a theoretical foundation for
learning algorithms that learn without risking instability.

[1] A. Zecevic and D. D. Siljak, Control of Complex Systems, ser.
Communications and Control Engineering. Springer, 2010.
[2] H. K. Khalil and J. W. Grizzle, Nonlinear Systems. Prentice Hall,
1996, vol. 3.
[3] B. Tibken, “Estimation of the domain of attraction for polynomial
systems via LMIs,” in Proc. of the IEEE Conference on Decision and
Control (CDC), vol. 4, 2000, pp. 3860–3864.
[4] S. Ratschan and Z. She, “Providing a basin of attraction to a target
region of polynomial systems by computation of Lyapunov-like functions,” SIAM Journal on Control and Optimization, vol. 48, no. 7, pp.
4377–4394, 2010.
[5] P. Giesl and S. Hafstein, “Review on computational methods for
Lyapunov functions,” Discrete and Continuous Dynamical Systems,
Series B, vol. 20, no. 8, pp. 2291–2337, 2015.
[6] F. Amato, C. Cosentino, and A. Merola, “On the region of attraction of
nonlinear quadratic systems,” Automatica, vol. 43, no. 12, pp. 2119–
2123, 2007.
[7] H. Burchardt and S. Ratschan, “Estimating the region of attraction
of ordinary differential equations by quantified constraint solving,” in
Proc. of the WSEAS International Conference on Dynamical Systems
and Control, 2007, pp. 241–246.
[8] D. Henrion and M. Korda, “Convex computation of the region of
attraction of polynomial control systems,” IEEE Transactions on
Automatic Control, vol. 59, no. 2, pp. 297–312, 2014.
[9] R. Bobiti and M. Lazar, “A sampling approach to constructing Lyapunov functions for nonlinear continuous-time systems,” in Proc. of
the Conference on Decision and Control (CDC), 2016.
[10] A. Trofino, “Robust stability and domain of attraction of uncertain
nonlinear systems,” in Proc. of the American Control Conference
(ACC), vol. 5, 2000, pp. 3707–3711.
[11] D. Coutinho, A. Trofino, and M. Fu, “Guaranteed cost control of
uncertain nonlinear systems via polynomial Lyapunov functions,”
IEEE Transactions on Automatic Control, vol. 47, no. 9, pp. 1575–
1580, 2002.
[12] U. Topcu, A. K. Packard, P. Seiler, and G. Balas, “Robust regionof-attraction estimation,” IEEE Transactions on Automatic Control,
vol. 55, no. 1, pp. 137–142, 2010.
[13] K. Zhou and J. C. Doyle, Essentials of Robust Control. Prentice Hall,
1998, vol. 104.
[14] C. E. Rasmussen and C. K. Williams, Gaussian Processes for Machine
Learning. MIT Press, 2006.
[15] F. Berkenkamp and A. P. Schoellig, “Safe and robust learning control
with Gaussian processes,” in Proc. of the European Control Conference (ECC), 2015, pp. 2501–2506.
[16] A. K. Akametalu, S. Kaynama, J. F. Fisac, M. N. Zeilinger, J. H.
Gillula, and C. J. Tomlin, “Reachability-based safe learning with
Gaussian processes,” in Proc. of the IEEE Conference on Decision
and Control (CDC), 2014, pp. 1424–1431.
[17] J. Mockus, Bayesian Approach to Global Optimization, ser. Mathematics and Its Applications, M. Hazewinkel, Ed. Springer, 1989,
vol. 37.
[18] N. Srinivas, A. Krause, S. M. Kakade, and M. Seeger, “Gaussian
process optimization in the bandit setting: No regret and experimental
design,” in Proc. of the International Conference on Machine Learning
(ICML), 2010, pp. 1015–1022.
[19] Y. Sui, A. Gotovos, J. W. Burdick, and A. Krause, “Safe exploration
for optimization with Gaussian processes,” in Proc. of the International
Conference on Machine Learning (ICML), 2015, pp. 997–1005.
[20] F. Berkenkamp, A. P. Schoellig, and A. Krause, “Safe controller
optimization for quadrotors with Gaussian processes,” in Proc. of the
IEEE International Conference on Robotics and Automation (ICRA),
2016, pp. 493–496, (to appear), arXiv:1509.01066 [cs.RO].
[21] F. Berkenkamp, A. Krause, and Angela P. Schoellig, “Bayesian optimization with safety constraints: Safe and automatic parameter tuning
in robotics.” arXiv, 2016, arXiv:1602.04450 [cs.RO].
[22] A. Christmann and I. Steinwart, Support Vector Machines, ser. Information Science and Statistics. Springer, 2008.
[23] F. Berkenkamp, A. Krause, and A. P. Schoellig, “Kernel selection
for reliable and data-efficient learning of Gaussian process dyanmic
models.” arXiv, 2016.

A PPENDIX
A. Full Exploration Proof
We consider a variant of Algorithm 1 that determines
the safe set, Sn , starting from S0 in Assumption 4, using
the Lipschitz constant directly. This allows us to prove
exploration guarantees.
First, let us consider the best ROA approximation that
we could hope to achieve. We have a probabilistic model
of gπ (x), which means we cannot expect to explore the
full ROA. Instead, we consider learning the true ROA up
to some error, . An algorithm with knowledge up to this
error, could use the Lipschitz constant of V̇ from Lemma 4
in order to learn which discrete states close to S0 fulfill the
requirement of V̇ (x) < −Lτ from Lemma 6. We define a
general operator on a set that determines these states,
R,V̇ (S) = S ∪ {x ∈ Xτ | ∃x0 ∈ Rl (S) :

(16)
0

V̇ (x) +  + Lkx − x k < −Lτ },

where Rl (S) is an operator that selects the maximum level
set of the Lyapunov function contained in S,


Rl (S) = V argmax c, subject to V(c) ∩ Xτ ⊆ S .
c>0

(17)
Thus, (16) adds a state x ∈ Xτ to S if there exists a state x0
in the current level set estimate of the ROA that can be used
in order to determine that V̇ (x) < −Lτ via the Lipschitz
constant. If we apply this operator iteratively, in the limit
we reach the maximum ROA that can be determined using
i−1
i
knowledge up to . That is, with R,
(S) = R,V̇ (R,
(S))
V̇
V̇
1
and R,V̇ (S) = R,V̇ (S), the maximum safe set that could
be determined is
i
(S).
Rl (R,V̇ (S0 )), where R,V̇ (S) = lim R,
V̇
i→∞

(18)

With the baseline defined, in the following we
provide an algorithm that achieves this baseline using the results from Sec. IV. Instead of defining
the ROA directly in terms of the GP confidence interval of V̇ , Qn,V̇ (x) = [ln (x), un (x)] as in Theorem 1, we consider the intersection of these intervals,
Cn (x) = Qn,V̇ (x) ∩ Cn−1 (x). We initialize these intervals so
that points in S0 are safe; that is, C0 (x) = [−∞, 0) if x ∈ S0
and [−∞, ∞] otherwise. As a consequence of Lemma 1,
we have that V̇ (x) ∈ Cn (x) for all n ≥ 1 with probability
at least (1 − δ). Besides ensuring that points in S0 are
always considered safe, this definition guarantees that the
estimated safe set will not decrease over time. We define ln,c (x) = minx∈X Cn (x) and un,c (x) = maxx∈X Cn (x).
Based on these confidence intervals, after (n − 1) measurements we can quantify the following data points as
having a value V̇ (x) that is sufficiently negative,
Sn,V̇ = S0 ∪ {x ∈ X | ∃x0 ∈ Rl (Sn−1,V̇ ) :
0

(19)
0

un,c (x ) + Lkx − x k < −Lτ },

and estimate the ROA based on this set as
Sn = Rl (Sn,V̇ ).

(20)

Intuitively this selection criterion is similar to the one of the
baseline in (16), except that (19) uses the GP confidence
intervals instead of perfect model knowledge. This allows
us to prove that we explore the maximum level set up to 
accuracy:
Theorem 2. Under the assumptions of Theorem 1, choose βn
as in Lemma 1 and let n∗ be the smallest positive integer so
that
CL2∂V (|R0,V̇ (S0 )| + 1)
n∗
≥
,
(21)
β n∗ γ n∗
2
where C = 8/ log(1 + σ −2 ). For any  > 0, and δ ∈ (0, 1),
under the selection criterion (14) within Sn from (20), the
following holds jointly with probability at least (1 − δ):
(i) Sn is contained in the ROA of (1) for all n ≥ 1,
(ii) Rl (R,V̇ (S0 )) ⊆ Sn∗ ⊆ Rl (R0,V̇ (S0 )).

Proof. Statement (i) follows from Theorem 1. For (ii), we
have from [19, Lemma 5] that if the safe set does not
increase, then the maximum uncertainty of gπ and V̇ within
the safe set is bounded by /k∂V /∂xk∞ = /L∂V and ,
respectively, after a finite number of iterations. At that point,
either the safe set increases similarly to the baseline, or we
have
explored the full safe set [19, Lemma 7]. Applying this

|R0,V̇ (S0 )| + 1 times provides the result [19, Cor. 4].
That is, after a finite number of evaluations, n∗ , we explore
at least as much as the baseline up to accuracy , but not more
than we could determine as safe if we knew the function
perfectly; that is, the baseline with  = 0.

