arXiv:1705.08551v3 [stat.ML] 13 Nov 2017

Safe Model-based Reinforcement Learning with
Stability Guarantees
Felix Berkenkamp
Department of Computer Science
ETH Zurich
befelix@inf.ethz.ch

Matteo Turchetta
Department of Computer Science,
ETH Zurich
matteotu@inf.ethz.ch

Angela P. Schoellig
Institute for Aerospace Studies
University of Toronto
schoellig@utias.utoronto.ca

Andreas Krause
Department of Computer Science
ETH Zurich
krausea@ethz.ch

Abstract
Reinforcement learning is a powerful paradigm for learning optimal policies from
experimental data. However, to find optimal policies, most reinforcement learning
algorithms explore all possible actions, which may be harmful for real-world systems. As a consequence, learning algorithms are rarely applied on safety-critical
systems in the real world. In this paper, we present a learning algorithm that
explicitly considers safety, defined in terms of stability guarantees. Specifically,
we extend control-theoretic results on Lyapunov stability verification and show
how to use statistical models of the dynamics to obtain high-performance control
policies with provable stability certificates. Moreover, under additional regularity
assumptions in terms of a Gaussian process prior, we prove that one can effectively
and safely collect data in order to learn about the dynamics and thus both improve
control performance and expand the safe region of the state space. In our experiments, we show how the resulting algorithm can safely optimize a neural network
policy on a simulated inverted pendulum, without the pendulum ever falling down.

1

Introduction

While reinforcement learning (RL, [1]) algorithms have achieved impressive results in games, for
example on the Atari platform [2], they are rarely applied to real-world physical systems (e.g., robots)
outside of academia. The main reason is that RL algorithms provide optimal policies only in the
long-term, so that intermediate policies may be unsafe, break the system, or harm their environment.
This is especially true in safety-critical systems that can affect human lives. Despite this, safety in RL
has remained largely an open problem [3].
Consider, for example, a self-driving car. While it is desirable for the algorithm that drives the
car to improve over time (e.g., by adapting to driver preferences and changing environments), any
policy applied to the system has to guarantee safe driving. Thus, it is not possible to learn about the
system through random exploratory actions, which almost certainly lead to a crash. In order to avoid
this problem, the learning algorithm needs to consider its ability to safely recover from exploratory
actions. In particular, we want the car to be able to recover to a safe state, for example, driving at a
reasonable speed in the middle of the lane. This ability to recover is known as asymptotic stability
in control theory [4]. Specifically, we care about the region of attraction of the closed-loop system
under a policy. This is a subset of the state space that is forward invariant so that any state trajectory
that starts within this set stays within it for all times and converges to a goal state eventually.
31st Conference on Neural Information Processing Systems (NIPS 2017), Long Beach, CA, USA.

In this paper, we present a RL algorithm for continuous state-action spaces that provides these kind
of high-probability safety guarantees for policies. In particular, we show how, starting from an initial,
safe policy we can expand our estimate of the region of attraction by collecting data inside the safe
region and adapt the policy to both increase the region of attraction and improve control performance.
Related work Safety is an active research topic in RL and different definitions of safety exist [5, 6].
Discrete Markov decision processes (MDPs) are one class of tractable models that have been analyzed.
In risk-sensitive RL, one specifies risk-aversion in the reward [7]. For example, [8] define risk as
the probability of driving the agent to a set of known, undesirable states. Similarly, robust MDPs
maximize rewards when transition probabilities are uncertain [9, 10]. Both [11] and [12] introduce
algorithms to safely explore MDPs so that the agent never gets stuck without safe actions. All these
methods require an accurate probabilistic model of the system.
In continuous state-action spaces, model-free policy search algorithms have been successful. These
update policies without a system model by repeatedly executing the same task [13]. In this setting, [14] introduces safety guarantees in terms of constraint satisfaction that hold in expectation.
High-probability worst-case safety guarantees are available for methods based on Bayesian optimization [15] together with Gaussian process models (GP, [16]) of the cost function. The algorithms
in [17] and [18] provide high-probability safety guarantees for any parameter that is evaluated on
the real system. These methods are used in [19] to safely optimize a parametric control policy on a
quadrotor. However, resulting policies are task-specific and require the system to be reset.
In the model-based RL setting, research has focused on safety in terms of state constraints. In
[20, 21], a priori known, safe global backup policies are used, while [22] learns to switch between
several safe policies. However, it is not clear how one may find these policies in the first place.
Other approaches use model predictive control with constraints, a model-based technique where the
control actions are optimized online. For example, [23] models uncertain environmental constraints,
while [24] uses approximate uncertainty propagation of GP dynamics along trajectories. In this
setting, robust feasability and constraint satisfaction can be guaranteed for a learned model with
bounded errors using robust model predictive control [25]. The method in [26] uses reachability
analysis to construct safe regions in the state space. The theoretical guarantees depend on the solution
to a partial differential equation, which is approximated.
Theoretical guarantees for the stability exist for the more tractable stability analysis and verification
under a fixed control policy. In control, stability of a known system can be verified using a Lyapunov
function [27]. A similar approach is used by [28] for deterministic, but unknown dynamics that are
modeled as a GP, which allows for provably safe learning of regions of attraction for fixed policies.
Similar results are shown in [29] for stochastic systems that are modeled as a GP. They use Bayesian
quadrature to compute provably accurate estimates of the region of attraction. These approaches do
not update the policy.
Our contributions We introduce a novel algorithm that can safely optimize policies in continuous
state-action spaces while providing high-probability safety guarantees in terms of stability. Moreover,
we show that it is possible to exploit the regularity properties of the system in order to safely learn
about the dynamics and thus improve the policy and increase the estimated safe region of attraction
without ever leaving it. Specifically, starting from a policy that is known to stabilize the system
locally, we gather data at informative, safe points and improve the policy safely based on the improved
model of the system and prove that any exploration algorithm that gathers data at these points reaches
a natural notion of full exploration. We show how the theoretical results transfer to a practical
algorithm with safety guarantees and apply it to a simulated inverted pendulum stabilization task.

2

Background and Assumptions

We consider a deterministic, discrete-time dynamic system
xt+1 = f (xt , ut ) = h(xt , ut ) + g(xt , ut ),

(1)

with states x ∈ X ⊂ Rq and control actions u ∈ U ⊂ Rp and a discrete time index t ∈ N. The true
dynamics f : X × U → X consist of two parts: h(xt , ut ) is a known, prior model that can be
obtained from first principles, while g(xt , ut ) represents a priori unknown model errors. While the
model errors are unknown, we can obtain noisy measurements of f (x, u) by driving the system to
the state x and taking action u. We want this system to behave in a certain way, e.g., the car driving
2

on the road. To this end, we need to specify a control policy π : X → U that, given the current state,
determines the appropriate control action that drives the system to some goal state, which we set as
the origin without loss of generality [4]. We encode the performance requirements of how to drive
the system to the origin through a positive cost r(x, u) that is associated with states and actions and
has r(0, 0) = 0. The policy aims to minimize the cumulative, discounted costs for each starting state.
The goal is to safely learn about the dynamics from measurements and adapt the policy for performance, without encountering system failures. Specifically, we define the safety constraint on the
state divergence that occurs when leaving the region of attraction. This means that adapting the
policy is not allowed to decrease the region of attraction and exploratory actions to learn about the
dynamics f (·) are not allowed to drive the system outside the region of attraction. The region of
attraction is not known a priori, but is implicitly defined through the system dynamics and the choice
of policy. Thus, the policy not only defines performance as in typical RL, but also determines safety
and where we can obtain measurements.
Model assumptions In general, this kind of safe learning is impossible without further assumptions.
For example, in a discontinuous system even a slight change in the control policy can lead to drastically
different behavior. Moreover, to expand the safe set we need to generalize learned knowledge about
the dynamics to (potentially unsafe) states that we have not visited. To this end, we restrict ourselves
to the general and practically relevant class of models that are Lipschitz continuous. This is a typical
assumption in the control community [4]. Additionally, to ensure that the closed-loop system remains
Lipschitz continuous when the control policy is applied, we restrict policies to the rich class of
Lπ -Lipschitz continuous functions ΠL , which also contains certain types of neural networks [30].
Assumption 1 (continuity). The dynamics h(·) and g(·) in (1) are Lh - and Lg Lipschitz continuous
with respect to the 1-norm. The considered control policies π lie in a set ΠL of functions that
are Lπ -Lipschitz continuous with respect to the 1-norm.
To enable safe learning, we require a reliable statistical model. While we commit to GPs for the
exploration analysis, for safety any suitable, well-calibrated model is applicable.
Assumption 2 (well-calibrated model). Let µn (·) and Σn (·) denote the posterior mean and covariance matrix functions of the statistical model of the dynamics (1) conditioned on n noisy measurements.
1/2
With σn (·) = trace(Σn (·)), there exists a βn > 0 such that with probability at least (1 − δ) it holds
for all n ≥ 0, x ∈ X , and u ∈ U that kf (x, u) − µn (x, u)k1 ≤ βn σn (x, u).
This assumption ensures that we can build confidence intervals on the dynamics that, when scaled by
an appropriate constant βn , cover the true function with high probability. We introduce a specific
statistical model that fulfills both assumptions under certain regularity assumptions in Sec. 3.
Lyapunov function To satisfy the specified safety constraints for safe learning, we require a tool
to determine whether individual states and actions are safe. In control theory, this safety is defined
through the region of attraction, which can be computed for a fixed policy using Lyapunov functions [4]. Lyapunov functions are continuously differentiable functions v : X → R≥0 with v(0) = 0
and v(x) > 0 for all x ∈ X \ {0}. The key idea behind using Lyapunov functions to show stability
of the system (1) is similar to that of gradient descent on strictly quasiconvex functions: if one can
show that, given a policy π, applying the dynamics f on the state maps it to strictly smaller values
on the Lyapunov function (‘going downhill’), then the state eventually converges to the equilibrium
point at the origin (minimum). In particular, the assumptions in Theorem 1 below imply that v is
strictly quasiconvex within the region of attraction if the dynamics are Lipschitz continuous. As a
result, the one step decrease property for all states within a level set guarantees eventual convergence
to the origin.
Theorem 1 ([4]). Let v be a Lyapunov function, f Lipschitz continuous dynamics, and π a policy. If
v(f (x, π(x))) < v(x) for all x within the level set V(c) = {x ∈ X \ {0} | v(x) ≤ c}, c > 0, then
V(c) is a region of attraction, so that x0 ∈ V(c) implies xt ∈ V(c) for all t > 0 and limt→∞ xt = 0.
It is convenient to characterize the region of attraction through a level set of the Lyapunov function,
since it replaces the challenging test for convergence with a one-step decrease condition on the
Lyapunov function. For the theoretical analysis in this paper, we assume that a Lyapunov function is
given to determine the region of attraction. For ease of notation, we also assume ∂v(x)/∂x 6= 0 for
all x ∈ X \ 0, which ensures that level sets V(c) are connected if c > 0. Since Lyapunov functions
are continuously differentiable, they are Lv -Lipschitz continuous over the compact set X .
3

In general, it is not easy to find suitable Lyapunov functions. However, for physical models, like the
prior model h in (1), the energy of the system (e.g., kinetic and potential for mechanical systems) is a
good candidate Lyapunov function. Moreover, it has recently been shown that it is possible to compute
suitable Lyapunov functions [31, 32]. In our experiments, we exploit the fact that value functions in
RL are Lyapunov functions if the costs are strictly positive away from the origin. This follows directly
from the definition of the value function, where v(x) = r(x, π(x)) + v(f (x, π(x)) ≤ v(f (x, π(x))).
Thus, we can obtain Lyapunov candidates as a by-product of approximate dynamic programming.
Initial safe policy Lastly, we need to ensure that there exists a safe starting point for the learning
process. Thus, we assume that we have an initial policy π0 that renders the origin of the system in (1)
asymptotically stable within some small set of states S0x . For example, this policy may be designed
using the prior model h in (1), since most models are locally accurate but deteriorate in quality as
state magnitude increases. This policy is explicitly not safe to use throughout the state space X \ S0x .

3

Theory

In this section, we use these assumptions for safe reinforcement learning. We start by computing the
region of attraction for a fixed policy under the statistical model. Next, we optimize the policy in order
to expand the region of attraction. Lastly, we show that it is possible to safely learn about the dynamics
and, under additional assumptions about the model and the system’s reachability properties, that this
approach expands the estimated region of attraction safely. We consider an idealized algorithm that is
amenable to analysis, which we convert to a practical variant in Sec. 4. See Fig. 1 for an illustrative
run of the algorithm and examples of the sets defined below.
Region of attraction We start by computing the region of attraction for a fixed policy. This is an
extension of the method in [28] to discrete-time systems. We want to use the Lyapunov decrease condition in Theorem 1 to guarantee safety for the statistical model of the dynamics. However, the posterior
uncertainty in the statistical model of the dynamics means that one step predictions about v(f (·)) are
uncertain too. We account for this by constructing high-probability confidence intervals on v(f (x, u)):
Qn (x, u) := [v(µn−1 (x, u)) ± Lv βn σn−1 (x, u)]. From Assumption 2 together with the Lipschitz
property of v, we know that v(f (x, u)) is contained in Qn (x, u) with probability at least (1 − δ). For
our exploration analysis, we need to ensure that safe state-actions cannot become unsafe; that is, an
initial set of safe set S0 remains safe (defined later). To this end, we intersect the confidence intervals:
Cn (x, u) := Cn−1 ∩ Qn (x, u), where the set C is initialized to C0 (x, u) = (−∞, v(x) − L∆v τ )
when (x, u) ∈ S0 and C0 (x, u) = R otherwise. Note that v(f (x, u)) is contained in Cn (x, u) with
the same (1 − δ) probability as in Assumption 2. The upper and lower bounds on v(f (·)) are defined
as un (x, u) := max Cn (x, u) and ln (x, u) := min Cn (x, u).

Given these high-probability confidence intervals, the system is stable according to Theorem 1 if
v(f (x, u)) ≤ un (x) < v(x) for all x ∈ V(c). However, it is intractable to verify this condition
directly on the continuous domain without additional, restrictive assumptions about the model.
Instead, we consider a discretization of the state space Xτ ⊂ X into cells, so that kx − [x]τ k1 ≤ τ
holds for all x ∈ X . Here, [x]τ denotes the point in Xτ with the smallest l1 distance to x. Given this
discretization, we bound the decrease variation on the Lyapunov function for states in Xτ and use the
Lipschitz continuity to generalize to the continuous state space X .
Theorem 2. Under Assumptions 1 and 2 with L∆v := Lv Lf (Lπ + 1) + Lv , let Xτ be a discretization of X such that kx − [x]τ k1 ≤ τ for all x ∈ X . If, for all x ∈ V(c) ∩ Xτ with c > 0, u = π(x),
and for some n ≥ 0 it holds that un (x, u) < v(x) − L∆v τ, then v(f (x, π(x))) < v(x) holds for all
x ∈ V(c) with probability at least (1 − δ) and V(c) is a region of attraction for (1) under policy π.
The proof is given in Appendix A.1. Theorem 2 states that, given confidence intervals on the statistical
model of the dynamics, it is sufficient to check the stricter decrease condition in Theorem 2 on the
discretized domain Xτ to guarantee the requirements for the region of attraction in the continuous
domain in Theorem 1. The bound in Theorem 2 becomes tight as the discretization constant τ
and |v(f (·)) − un (·)| go to zero. Thus, the discretization constant trades off computation costs for
accuracy, while un approaches v(f
√ (·)) as we obtain more measurement data and the posterior model
uncertainty about the dynamics, βn σn decreases. The confidence intervals on v(f (x, π(x)) − v(x)
and the corresponding estimated region of attraction (red line) can be seen in the bottom half of Fig. 1.

Policy optimization So far, we have focused on estimating the region of attraction for a fixed policy.
Safety is a property of states under a fixed policy. This means that the policy directly determines
4

π30

π15

policy π0

σ (x, u)

S30

V(c15 )

u0

−L∆v τ

l0

state x
(a) Initial safe set (in red).

state x
(b) Exploration: 15 data points.

state x

∆v(x, π(x))

v(x)

action u

D30

(c) Final policy after 30 evaluations.

Figure 1: Example application of Algorithm 1. Due to input constraints, the system becomes unstable
for large states. We start from an initial, local policy π0 that has a small, safe region of attraction (red
lines) in Fig. 1(a). The algorithm selects safe, informative state-action pairs within Sn (top, white
shaded), which can be evaluated without leaving the region of attraction V(cn ) (red lines) of the
current policy πn . As we gather more data (blue crosses), the uncertainty in the model decreases
(top, background) and we use (3) to update the policy so that it lies within Dn (top, red shaded) and
fulfills the Lyapunov decrease condition. The algorithm converges to the largest safe set in Fig. 1(c).
It improves the policy without evaluating unsafe state-action pairs and thereby without system failure.
which states are safe. Specifically, to form a region of attraction all states in the discretizaton Xτ
within a level set of the Lyapunov function need to fulfill the decrease condition in Theorem 2 that
depends on the policy choice. The set of all state-action pairs that fulfill this decrease condition is
given by

Dn = (x, u) ∈ Xτ × U | un (x, u) − v(x) < −L∆v τ ,
(2)

see Fig. 1(c) (top, red shaded). In order to estimate the region of attraction based on this set, we
need to commit to a policy. Specifically, we want to pick the policy that leads to the largest possible
region of attraction according to Theorem 2. This requires that for each discrete state in Xτ the
corresponding state-action pair under the policy must be in the set Dn . Thus, we optimize the policy
according to
πn , cn =

argmax c,
π∈ΠL ,c∈R>0

such that for all x ∈ V(c) ∩ Xτ : (x, π(x)) ∈ Dn .

(3)

The region of attraction that corresponds to the optimized policy πn according to (3) is given
by V(cn ), see Fig. 1(b). It is the largest level set of the Lyapunov function for which all state-action
pairs (x, πn (x)) that correspond to discrete states within V(cn ) ∩ Xτ are contained in Dn . This
means that these state-action pairs fulfill the requirements of Theorem 2 and V(cn ) is a region of
attraction of the true system under policy πn . The following theorem is thus a direct consequence
of Theorem 2 and (3).
Theorem 3. Let Rπn be the true region of attraction of (1) under the policy πn . For any δ ∈ (0, 1),
we have with probability at least (1 − δ) that V(cn ) ⊆ Rπn for all n > 0.
Thus, when we optimize the policy subject to the constraint in (3) the estimated region of attraction is
always an inner approximation of the true region of attraction. However, solving the optimization
problem in (3) is intractable in general. We approximate the policy update step in Sec. 4.
Collecting measurements Given these stability guarantees, it is natural to ask how one might obtain
data points in order to improve the model of g(·) and thus efficiently increase the region of attraction.
This question is difficult to answer in general, since it depends on the property of the statistical model.
In particular, for general statistical models it is often not clear whether the confidence intervals
contract sufficiently quickly. In the following, we make additional assumptions about the model and
reachability within V(cn ) in order to provide exploration guarantees. These assumptions allow us to
highlight fundamental requirements for safe data acquisition and that safe exploration is possible.
5

We assume that the unknown model errors g(·) have bounded norm in a reproducing kernel Hilbert
space (RKHS, [33]) corresponding to a differentiable
kernel k, kg(·)kk ≤ Bg . These are a class of
P∞
well-behaved functions of the form g(z) = i=0 αi k(zi , z) defined through representer points zi
and weights αi that decay sufficiently fast with i. This assumption ensures
p that g satisfies the
Lipschitz property in Assumption 1, see [28]. Moreover, with βn = Bg + 4σ γn + 1 + ln(1/δ) we
can use GP models for the dynamics that fulfill Assumption 2 if the state if fully observable and the
measurement noise is σ-sub-Gaussian (e.g., bounded in [−σ, σ]), see [34]. Here γn is the information
capacity. It corresponds to the amount of mutual information that can be obtained about g from nq
measurements, a measure of the size of the function class encoded by the model. The information
capacity has a sublinear dependence on n for common kernels and upper bounds can be computed
efficiently [35]. More details about this model are given in Appendix A.2.
In order to quantify the exploration properties of our algorithm, we consider a discrete action
space Uτ ⊂ U. We define exploration as the number of state-action pairs in Xτ × Uτ that we can
safely learn about without leaving the true region of attraction. Note that despite this discretization,
the policy takes values on the continuous domain. Moreover, instead of using the confidence intervals
directly as in (3), we consider an algorithm that uses the Lipschitz constants to slowly expand the safe
set. We use this in our analysis to quantify the ability to generalize beyond the current safe set. In
practice, nearby states are sufficiently correlated under the model to enable generalization using (2).
Suppose we are given a set S0 of state-action pairs about which we can learn safely. Specifically, this
means that we have a policy such that, for any state-action pair (x, u) in S0 , if we apply action u in
state x and then apply actions according to the policy, the state converges to the origin. Such a set
can be constructed using the initial policy π0 from Sec. 2 as S0 = {(x, π0 (x)) | x ∈ S0x }. Starting
from this set, we want to update the policy to expand the region of attraction according to Theorem 2.
To this end, we use the confidence intervals on v(f (·)) for states inside S0 to determine state-action
pairs that fulfill the decrease condition. We thus redefine Dn for the exploration analysis to
[
 0
Dn =
z ∈ Xτ × Uτ | un (x, u) − v(x) + L∆v kz0 − (x, u)k1 < −L∆v τ . (4)
(x,u)∈Sn−1

This formulation is equivalent to (2), except that it uses the Lipschitz constant to generalize safety.
Given Dn , we can again find a region of attraction V(cn ) by committing to a policy according to (3).
In order to expand this region of attraction effectively we need to decrease the posterior model
uncertainty about the dynamics of the GP by collecting measurements. However, to ensure safety
as outlined in Sec. 2, we are not only restricted to states within V(cn ), but also need to ensure that
the state after taking an action is safe; that is, the dynamics map the state back into the region of
attraction V(cn ). We again use the Lipschitz constant in order to determine this set,
[ 
Sn =
z0 ∈ V(cn ) ∩ Xτ × Uτ | un (z) + Lv Lf kz − z0 k1 ≤ cn }.
(5)
z∈Sn−1

The set Sn contains state-action pairs that we can safely evaluate under the current policy πn without
leaving the region of attraction, see Fig. 1 (top, white shaded).
What remains is to define a strategy for collecting data points within Sn to effectively decrease model
uncertainty. We specifically focus on the high-level requirements for any exploration scheme without
committing to a specific method. In practice, any (model-based) exploration strategy that aims to
decrease model uncertainty by driving the system to specific states may be used. Safety can be
ensured by picking actions according to πn whenever the exploration strategy reaches the boundary
of the safe region V(cn ); that is, when un (x, u) > cn . This way, we can use πn as a backup policy
for exploration.
The high-level goal of the exploration strategy is to shrink the confidence intervals at state-action
pairs Sn in order to expand the safe region. Specifically, the exploration strategy should aim to visit
state-action pairs in Sn at which we are the most uncertain about the dynamics; that is, where the
confidence interval is the largest:
(xn , un ) = argmax un (x, u) − ln (x, u).

(6)

(x,u)∈Sn

As we keep collecting data points according to (6), we decrease the uncertainty about the dynamics
for different actions throughout the region of attraction and adapt the policy, until eventually we
6

Algorithm 1 S AFE LYAPUNOV L EARNING
1: Input: Initial safe policy π0 , dynamics model GP(µ(z), k(z, z0 ))
2: for all n = 1, . . . do
3:
Compute policy πn via SGD on (7)
4:
cn = argmaxc c, such that ∀x ∈ V(cn ) ∩ Xτ : un (x, πn (x)) − v(x) < −L∆v τ
5:
Sn = {(x, u) ∈ V(cn ) × Uτ | un (x, u) ≤ cn }
6:
Select (xn , un ) within Sn using (6) and drive system there with backup policy πn
7:
Update GP with measurements f (xn , un ) + n
have gathered enough information in order to expand it. While (6) implicitly assumes that any
state within V(cn ) can be reached by the exploration policy, it achieves the high-level goal of any
exploration algorithm that aims to reduce model uncertainty. In practice, any safe exploration scheme
is limited by unreachable parts of the state space.
We compare the active learning scheme in (6) to an oracle baseline that starts from the same initial
safe set S0 and knows v(f (x, u)) up to  accuracy within the safe set. The oracle also uses knowledge
about the Lipschitz constants and the optimal policy in ΠL at each iteration. We denote the set that this
baseline manages to determine as safe with R (S0 ) and provide a detailed definition in Appendix A.3.
Theorem 4. Assume σ-sub-Gaussian measurement noise and that the model error g(·)
in (1) has RKHSpnorm smaller than Bg .
Under the assumptions of Theorem 2,
with βn = Bg + 4σ γn + 1 + ln(1/δ), and with measurements collected according to (6), let
∗

0 )|+1)
n∗ be the smallest positive integer so that β 2 nγ ∗ ≥ Cq(|R(S
where C = 8/ log(1 + σ −2 ).
L2v 2
n∗ n
Let Rπ be the true region of attraction of (1) under a policy π. For any  > 0, and δ ∈ (0, 1), the
following holds jointly with probability at least (1 − δ) for all n > 0:

(i) V(cn ) ⊆ Rπn

(ii) f (x, u) ∈ Rπn ∀(x, u) ∈ Sn .

(iii) R (S0 ) ⊆ Sn ⊆ R0 (S0 ).

Theorem 4 states that, when selecting data points according to (6), the estimated region of attraction V(cn ) is (i) contained in the true region of attraction under the current policy and (ii) selected
data points do not cause the system to leave the region of attraction. This means that any exploration
method that considers the safety constraint (5) is able to safely learn about the system without leaving
the region of attraction. The last part of Theorem 4, (iii), states that after a finite number of data
points n∗ we achieve at least the exploration performance of the oracle baseline, while we do not
classify unsafe state-action pairs as safe. This means that the algorithm explores the largest region
of attraction possible for a given Lyapunov function with residual uncertaint about v(f (·)) smaller
than . Details of the comparison baseline are given in the appendix. In practice, this means that any
exploration method that manages to reduce the maximal uncertainty about the dynamics within Sn is
able to expand the region of attraction.
An example run of repeatedly evaluating (6) for a one-dimensional state-space is shown in Fig. 1. It
can be seen that, by only selecting data points within the current estimate of the region of attraction,
the algorithm can efficiently optimize the policy and expand the safe region over time.

4

Practical Implementation and Experiments

In the previous section, we have given strong theoretical results on safety and exploration for an
idealized algorithm that can solve (3). In this section, we provide a practical variant of the theoretical
algorithm in the previous section. In particular, while we retain safety guarantees, we sacrifice
exploration guarantees to obtain a more practical algorithm. This is summarized in Algorithm 1.
The policy optimization problem in (3) is intractable to solve and only considers safety, rather
than a performance metric. We propose to use an approximate policy update that that maximizes
approximate performance while providing stability guarantees. It proceeds by optimizing the policy
first and then computes the region of attraction V(cn ) for the new, fixed policy. This does not
impact safety, since data is still only collected inside the region of attraction. Moreover, should the
optimization fail and the region of attraction decrease, one can always revert to the previous policy,
which is guaranteed to be safe.
7

V(c50 )

angle [rad]

angular velocity [rad/s]

0.30

5

V(c0 )

0

unsafe region

−5

safely optimized policy
initial policy

0.25
0.20
0.15

π0

0.10

π50

0.05
0.00

−1.0

−0.5

0.0

0.5

1.0

0.0

angle [rad]

0.5

1.0

1.5

time [s]

(a) Estimated safe set.

(b) State trajectory (lower is better).

Figure 2: Optimization results for an inverted pendulum. Fig. 2(a) shows the initial safe set (yellow)
under the policy π0 , while the green region represents the estimated region of attraction under the
optimized neural network policy. It is contained within the true region of attraction (white). Fig. 2(b)
shows the improved performance of the safely learned policy over the policy for the prior model.
In our experiments, we use approximate dynamic programming [36] to capture the performance of the
policy. Given a policy πθ with parameters θ, we compute an estimate of the cost-to-go Jπθ (·) for the
mean dynamics µn based on the cost r(x, u) ≥ 0. At each state, Jπθ (x) is the sum of γ-discounted
rewards encountered when following the policy πθ . The goal is to adapt the parameters of the policy
for minimum cost as measured by Jπθ , while ensuring that the safety constraint on the worst-case
decrease on the Lyapunov function in Theorem 2 is not violated. A Lagrangian formulation to this
constrained optimization problem is


X
πn = argmin
r(x, πθ (x)) + γJπθ (µn−1 (x, πθ (x)) + λ un (x, πθ (x)) − v(x) + L∆v τ , (7)
πθ ∈ΠL

x∈Xτ

where the first term measures long-term cost to go and λ ≥ 0 is a Lagrange multiplier for the safety
constraint from Theorem 2. In our experiments, we use the value function as a Lyapunov function
candidate, v = J with r(·, ·) ≥ 0, and set λ = 1. In this case, (7) corresponds to an high-probability
upper bound on the cost-to-go given the uncertainty in the dynamics. This is similar to worst-case
performance formulations found in robust MDPs [9, 10], which consider worst-case value functions
given parametric uncertainty in MDP transition model. Moreover, since L∆v depends on the Lipschitz
constant of the policy, this simultaneously serves as a regularizer on the parameters θ.
To verify safety, we use the GP confidence intervals ln and un directly, as in (2). We also use
confidence to compute Sn for the active learning scheme, see Algorithm 1, Line 5. In practice, we
do not need to compute the entire set Sn to solve (3), but can use a global optimization method or
even a random sampling scheme within V(cn ) to find suitable state-actions. Moreover, measurements
for actions that are far away from the current policy are unlikely to expand V(cn ), see Fig. 1(c). As
we optimize (7) via gradient descent, the policy changes only locally. Thus, we can achieve better
data-efficiency by restricting the exploratory actions u with (x, u) ∈ Sn to be close to πn , u ∈
[πn (x) − ū, πn (x) + ū] for some constant ū.

Computing the region of attraction by verifying the stability condition on a discretized domain suffers
from the curse of dimensionality. However, it is not necessary to update policies in real time. In
particular, since any policy that is returned by the algorithm is provably safe within some level set,
any of these policies can be used safely for an arbitrary number of time steps. To scale this method to
higher-dimensional system, one would have to consider an adaptive discretization for the verification
as in [27].
Experiments A Python implementation of Algorithm 1 and the experiments based on TensorFlow [37] and GPflow [38] is available at https://github.com/befelix/safe_learning.

We verify our approach on an inverted pendulum benchmark problem. The true, continuous-time
dynamics are given by ml2 ψ̈ = gml sin(ψ) − λψ̇ + u, where ψ is the angle, m the mass, g the
gravitational constant, and u the torque applied to the pendulum. The control torque is limited, so that
the pendulum necessarily falls down beyond a certain angle. We use a GP model for the discrete-time
dynamics, where the mean dynamics are given by a linearized and discretized model of the true
8

dynamics that considers a wrong, lower mass and neglects friction. As a result, the optimal policy for
the mean dynamics does not perform well and has a small region of attraction as it underactuates the
system. We use a combination of linear and Matérn kernels in order to capture the model errors that
result from parameter and integration errors.
For the policy, we use a neural network with two hidden layers and 32 neurons with ReLU activations
each. We compute a conservative estimate of the Lipschitz constant as in [30]. We use standard
approximate dynamic programming with a quadratic, normalized cost r(x, u) = xT Qx + uT Ru,
where Q and R are positive-definite, to compute the cost-to-go Jπθ . Specifically, we use a piecewiselinear triangulation of the state-space as to approximate Jπθ , see [39]. This allows us to quickly
verify the assumptions that we made about the Lyapunov function in Sec. 2 using a graph search. In
practice, one may use other function approximators. We optimize the policy via stochastic gradient
descent on (7).
The theoretical confidence intervals for the GP model are conservative. To enable more data-efficient
learning, we fix βn = 2. This corresponds to a high-probability decrease condition per-state, rather
than jointly over the state space. Moreover, we use local Lipschitz constants of the Lyapunov function
rather than the global one. While this does not affect guarantees, it greatly speeds up exploration.
For the initial policy, we use approximate dynamic programming to compute the optimal policy for
the prior mean dynamics. This policy is unstable for large deviations from the initial state and has poor
performance, as shown in Fig. 2(b). Under this initial, suboptimal policy, the system is stable within
a small region of the state-space Fig. 2(a). Starting from this initial safe set, the algorithm proceeds to
collect safe data points and improve the policy. As the uncertainty about the dynamics decreases, the
policy improves and the estimated region of attraction increases. The region of attraction after 50 data
points is shown in Fig. 2(a). The resulting set V(cn ) is contained within the true safe region of the
optimized policy πn . At the same time, the control performance improves drastically relative to the
initial policy, as can be seen in Fig. 2(b). Overall, the approach enables safe learning about dynamic
systems, as all data points collected during learning are safely collected under the current policy.

5

Conclusion

We have shown how classical reinforcement learning can be combined with safety constraints in terms
of stability. Specifically, we showed how to safely optimize policies and give stability certificates
based on statistical models of the dynamics. Moreover, we provided theoretical safety and exploration
guarantees for an algorithm that can drive the system to desired state-action pairs during learning. We
believe that our results present an important first step towards safe reinforcement learning algorithms
that are applicable to real-world problems.
Acknowledgments
This research was supported by SNSF grant 200020_159557, the Max Planck ETH Center for
Learning Systems, NSERC grant RGPIN-2014-04634, and the Ontario Early Researcher Award.

References
[1] Richard S. Sutton and Andrew G. Barto. Reinforcement learning: an introduction. MIT press,
1998.
[2] Volodymyr Mnih, Koray Kavukcuoglu, David Silver, Andrei A. Rusu, Joel Veness, Marc G.
Bellemare, Alex Graves, Martin Riedmiller, Andreas K. Fidjeland, Georg Ostrovski, Stig Petersen, Charles Beattie, Amir Sadik, Ioannis Antonoglou, Helen King, Dharshan Kumaran, Daan
Wierstra, Shane Legg, and Demis Hassabis. Human-level control through deep reinforcement
learning. Nature, 518(7540):529–533, 2015.
[3] Dario Amodei, Chris Olah, Jacob Steinhardt, Paul Christiano, John Schulman, and Dan Mané.
Concrete problems in AI safety. arXiv:1606.06565 [cs], 2016.
[4] Hassan K. Khalil and J. W. Grizzle. Nonlinear systems, volume 3. Prentice Hall, 1996.
9

[5] Martin Pecka and Tomas Svoboda. Safe exploration techniques for reinforcement learning –
an overview. In Modelling and Simulation for Autonomous Systems, pages 357–375. Springer,
2014.
[6] Javier García and Fernando Fernández. A comprehensive survey on safe reinforcement learning.
Journal of Machine Learning Research (JMLR), 16:1437–1480, 2015.
[7] Stefano P. Coraluppi and Steven I. Marcus. Risk-sensitive and minimax control of discrete-time,
finite-state Markov decision processes. Automatica, 35(2):301–309, 1999.
[8] Peter Geibel and Fritz Wysotzki. Risk-sensitive reinforcement learning applied to control under
constraints. J. Artif. Intell. Res.(JAIR), 24:81–108, 2005.
[9] Aviv Tamar, Shie Mannor, and Huan Xu. Scaling Up Robust MDPs by Reinforcement Learning.
In Proc. of the International Conference on Machine Learning (ICML), 2014.
[10] Wolfram Wiesemann, Daniel Kuhn, and Berç Rustem. Robust Markov Decision Processes.
Mathematics of Operations Research, 38(1):153–183, 2012.
[11] Teodor Mihai Moldovan and Pieter Abbeel. Safe exploration in Markov decision processes. In
Proc. of the International Conference on Machine Learning (ICML), pages 1711–1718, 2012.
[12] Matteo Turchetta, Felix Berkenkamp, and Andreas Krause. Safe exploration in finite markov
decision processes with gaussian processes. pages 4305–4313, 2016.
[13] Jan Peters and Stefan Schaal. Policy gradient methods for robotics. In Proc. of the IEEE/RSJ
International Conference on Intelligent Robots and Systems, pages 2219–2225, 2006.
[14] Joshua Achiam, David Held, Aviv Tamar, and Pieter Abbeel. Constrained policy optimization.
In Proc. of the International Conference on Machine Learning (ICML), 2017.
[15] Jonas Mockus. Bayesian approach to global optimization, volume 37 of Mathematics and Its
Applications. Springer, Dordrecht, 1989.
[16] Carl Edward Rasmussen and Christopher K.I Williams. Gaussian processes for machine
learning. MIT Press, Cambridge MA, 2006.
[17] Jens Schreiter, Duy Nguyen-Tuong, Mona Eberts, Bastian Bischoff, Heiner Markert, and Marc
Toussaint. Safe exploration for active learning with Gaussian processes. In Machine Learning
and Knowledge Discovery in Databases, number 9286, pages 133–149. Springer International
Publishing, 2015.
[18] Yanan Sui, Alkis Gotovos, Joel W. Burdick, and Andreas Krause. Safe exploration for optimization with Gaussian processes. In Proc. of the International Conference on Machine Learning
(ICML), pages 997–1005, 2015.
[19] Felix Berkenkamp, Angela P. Schoellig, and Andreas Krause. Safe controller optimization for
quadrotors with Gaussian processes. In Proc. of the IEEE International Conference on Robotics
and Automation (ICRA), pages 493–496, 2016.
[20] J. Garcia and F. Fernandez. Safe exploration of state and action spaces in reinforcement learning.
Journal of Artificial Intelligence Research, pages 515–564, 2012.
[21] Alexander Hans, Daniel Schneegaß, Anton Maximilian Schäfer, and Steffen Udluft. Safe
exploration for reinforcement learning. In Proc. of the European Symposium on Artificial
Neural Networks (ESANN), pages 143–148, 2008.
[22] Theodore J. Perkins and Andrew G. Barto. Lyapunov design for safe reinforcement learning.
The Journal of Machine Learning Research, 3:803–832, 2003.
[23] Dorsa Sadigh and Ashish Kapoor. Safe control under uncertainty with Probabilistic Signal
Temporal Logic. In Proc. of Robotics: Science and Systems, 2016.
10

[24] Chris J. Ostafew, Angela P. Schoellig, and Timothy D. Barfoot. Robust constrained learningbased NMPC enabling reliable mobile robot path tracking. The International Journal of Robotics
Research (IJRR), 35(13):1547–1536, 2016.
[25] Anil Aswani, Humberto Gonzalez, S. Shankar Sastry, and Claire Tomlin. Provably safe and
robust learning-based model predictive control. Automatica, 49(5):1216–1226, 2013.
[26] Anayo K. Akametalu, Shahab Kaynama, Jaime F. Fisac, Melanie N. Zeilinger, Jeremy H.
Gillula, and Claire J. Tomlin. Reachability-based safe learning with Gaussian processes. In
Proc. of the IEEE Conference on Decision and Control (CDC), pages 1424–1431, 2014.
[27] Ruxandra Bobiti and Mircea Lazar. A sampling approach to finding Lyapunov functions for
nonlinear discrete-time systems. In Proc. of the European Control Conference (ECC), pages
561–566, 2016.
[28] Felix Berkenkamp, Riccardo Moriconi, Angela P. Schoellig, and Andreas Krause. Safe learning
of regions of attraction in nonlinear systems with Gaussian processes. In Proc. of the Conference
on Decision and Control (CDC), pages 4661–4666, 2016.
[29] Julia Vinogradska, Bastian Bischoff, Duy Nguyen-Tuong, Henner Schmidt, Anne Romer, and
Jan Peters. Stability of controllers for Gaussian process forward models. In Proceedings of the
International Conference on Machine Learning (ICML), pages 545–554, 2016.
[30] Christian Szegedy, Wojciech Zaremba, Ilya Sutskever, Joan Bruna, Dumitru Erhan, Ian Goodfellow, and Rob Fergus. Intriguing properties of neural networks. In Proc. of the International
Conference on Learning Representations (ICLR), 2014.
[31] Huijuan Li and Lars Grüne. Computation of local ISS Lyapunov functions for discrete-time systems via linear programming. Journal of Mathematical Analysis and Applications, 438(2):701–
719, 2016.
[32] Peter Giesl and Sigurdur Hafstein. Review on computational methods for Lyapunov functions.
Discrete and Continuous Dynamical Systems, Series B, 20(8):2291–2337, 2015.
[33] Bernhard Schölkopf. Learning with kernels: support vector machines, regularization, optimization, and beyond. Adaptive computation and machine learning. MIT Press, Cambridge, Mass,
2002.
[34] Sayak Ray Chowdhury and Aditya Gopalan. On kernelized multi-armed bandits. arXiv preprint
arXiv:1704.00445, 2017.
[35] Niranjan Srinivas, Andreas Krause, Sham M. Kakade, and Matthias Seeger. Gaussian Process
Optimization in the Bandit Setting: No Regret and Experimental Design. IEEE Transactions on
Information Theory, 58(5):3250–3265, 2012.
[36] Warren B. Powell. Approximate dynamic programming: solving the curses of dimensionality.
John Wiley & Sons, 2007.
[37] Martín Abadi, Ashish Agarwal, Paul Barham, Eugene Brevdo, Zhifeng Chen, Craig Citro,
Greg S. Corrado, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Ian Goodfellow,
Andrew Harp, Geoffrey Irving, Michael Isard, Yangqing Jia, Rafal Jozefowicz, Lukasz Kaiser,
Manjunath Kudlur, Josh Levenberg, Dan Mane, Rajat Monga, Sherry Moore, Derek Murray,
Chris Olah, Mike Schuster, Jonathon Shlens, Benoit Steiner, Ilya Sutskever, Kunal Talwar, Paul
Tucker, Vincent Vanhoucke, Vijay Vasudevan, Fernanda Viegas, Oriol Vinyals, Pete Warden,
Martin Wattenberg, Martin Wicke, Yuan Yu, and Xiaoqiang Zheng. TensorFlow: Large-Scale
Machine Learning on Heterogeneous Distributed Systems. arXiv:1603.04467 [cs], 2016.
[38] Alexander G. de G. Matthews, Mark van der Wilk, Tom Nickson, Keisuke Fujii, Alexis
Boukouvalas, Pablo León-Villagrá, Zoubin Ghahramani, and James Hensman. GPflow: a
Gaussian process library using TensorFlow. Journal of Machine Learning Research, 18(40):1–6,
2017.
11

[39] Scott Davies. Multidimensional triangulation and interpolation for reinforcement learning. In
Proc. of the Conference on Neural Information Processing Systems (NIPS), pages 1005–1011,
1996.
[40] Andreas Christmann and Ingo Steinwart. Support Vector Machines. Information Science and
Statistics. Springer, New York, NY, 2008.

A
A.1

Proofs
Stability verification

Lemma 1. Using Assumptions 1 and 2, let Xτ be a discretization of X such that kx − [x]τ k1 ≤ τ
for all x ∈ X . Then, for all x ∈ X , we have with probability at least 1 − δ that

v(µn−1 ([z]τ )) − v([x]τ ) − v(f (z)) − v(x) ≤ Lv βn σn−1 ([z]τ ) + (Lv Lf (Lπ + 1) + Lv )τ, (8)
where z = (x, π(x)) and [z]τ = ([x]τ , π([x]τ )).
Proof. Let z = (x, π(x)), [z]τ = ([x]τ , π([x]τ )), and µ = µn−1 , σ = σn−1 . Then we have that

v(µ([z]τ )) − v([x]τ ) − v(f (z)) − v(x) ,

= v(µ([z]τ )) − v([x]τ ) − v(f (z)) + v(x) ,

= v(µ([z]τ )) − v(f ([z]τ )) + v(f ([z]τ )) − v(f (z)) + v(x) − v([x]τ ) ,

≤ v(µ([z]τ )) − v(f ([z]τ )) + v(f ([z]τ )) − v(f (z)) + v(x) − v([x]τ ) ,
≤ Lv kµ([z]τ ) − f ([z]τ )k1 + Lv kf ([z]τ ) − f (z)k1 + Lv kx − [x]τ k1 ,
≤ Lv βn σ([z]τ ) + Lv Lf k[z]τ − zk1 + Lv kx − [x]τ k1 ,

where the last three inequalities follow from Assumptions 1 and 2 to last inequality follows
from Lemma 3. The result holds with probability at least 1 − δ. By definition of the discretization
and the policy class ΠL we have on each grid cell that
kz − [z]τ k1 = kx − [x]τ k1 + kπ(x) − π([x]τ )k1 ,
≤ τ + Lπ kx − [x]τ k1 ,
≤ (Lπ + 1)τ,

where the equality in the first step follows from the definition of the 1-norm. Plugging this into the
previous bound yields

v(µ([z]τ )) − v([x]τ ) − v(f (z)) − v(x) ≤ Lv βn σ([z]τ ) + (Lv Lf (1 + Lπ ) + Lv ) τ,
which completes the proof.
Lemma 2. v(f (x, u)) ∈ Qn holds for all x ∈ X , u ∈ U, and n > 0 with probability at least (1 − δ).
Proof. The proof is analogous to Lemma 1 and follows from Assumptions 1 and 2.
Corollary 1. v(f (x, u)) ∈ Cn holds for all x ∈ X , u ∈ U, and n > 0 with probability at least
(1 − δ).
Proof. Direct consequence of the fact that Lemma 2 holds jointly for all n > 0 with probability at
least 1 − δ.
Lemma 1 show that the decrease on the Lyapunov function on the discrete grid Xτ is close to that
on the continuous domain X . Given these confidence intervals, we can now establish the region of
attraction using Theorem 1:
12

Theorem 2. Under Assumptions 1 and 2 with L∆v := Lv Lf (Lπ + 1) + Lv , let Xτ be a discretization of X such that kx − [x]τ k1 ≤ τ for all x ∈ X . If, for all x ∈ V(c) ∩ Xτ with c > 0, u = π(x),
and for some n ≥ 0 it holds that un (x, u) < v(x) − L∆v τ, then v(f (x, π(x))) < v(x) holds for all
x ∈ V(c) with probability at least (1 − δ) and V(c) is a region of attraction for (1) under policy π.
Proof. Using Lemma 1 it holds that v(f (x, π(x)) − v(x) < 0 for all continuous states x ∈ V(c)
with probability at least 1 − δ, since all discrete states xτ ∈ V(c) ∩ X fulfill Theorem 2. Thus we
can use Theorem 1 to conclude that V(c) is a region of attraction for (1).
Theorem 3. Let Rπn be the true region of attraction of (1) under the policy πn . For any δ ∈ (0, 1),
we have with probability at least (1 − δ) that V(cn ) ⊆ Rπn for all n > 0.
Proof. Following the definition of Dn in (2), it is clear from the constraint in the optimization
problem (3) that for all x ∈ Dn it holds that (x, πn (x)) ∈ Dn or, equivalently that un (x, π(x)) −
v(x) < −L∆v τ , see (2). The result V(cn ) ⊆ Rπn then follows from Theorem 2.
Note that the initialization of the confidence intervals Q0 ensures that the decrease condition is always
fulfilled for the initial policy.
A.2

Gaussian process model

One particular assumption that satisfies both the Lipschitz continuity and allows us to use GPs as a
model of the dynamics is that the model errors g(x, u) live in some reproducing kernel Hilbert space
(RKHS, [40]) corresponding to a differentiable kernel k and have RKHS norm smaller than Bg [35].
In our theoretical analysis, we use this assumption to prove exploration guarantees.
A GP(µ(z), k(z, z0 )) is a distribution over well-behaved, smooth functions f : X × U → R (see Remark 1 for the vector-case, Rq ) that is parameterized by a mean function µ and a covariance function
(kernel) k, which encodes assumptions about the functions [16]. In our case, the mean is given by
the prior model h, while the kernel corresponds to the one in the RKHS. Given noisy measurements
of the dynamics, fˆ(z) = f (z) +  with z = (x, u) at locations An = {z1 , . . . , zn }, corrupted
by independent, Gaussian noise  ∼ N (0, σ 2 ) (we relax the Gaussian noise assumption in our
analysis), the posterior is a GP distribution again with mean, µn (z) = kn (z)T (Kn + σ 2 I)−1 yn ,
covariance kn (z, z0 ) = k(z, z0 ) − kn (z)T (Kn + σ 2 I)−1 kn (z0 ), and variance σn2 (z) = kn (z, z).
The vector yn = [fˆ(z1 ) − h(z1 ), . . . , fˆ(zn ) − h(zn )]T contains observed, noisy deviations from
the mean, kn (z) = [k(z, z1 ), . . . , k(z, zn )] contains the covariances between the test input z and the
data points in Dn , Kn ∈ Rn×n has entries [Kn ](i,j) = k(zi , zj ), and I is the identity matrix.

Remark 1. In the case of multiple output dimensions (q > 1), we consider a function with onedimensional output f 0 (x, u, i) : X × U × I → R, with the output dimension indexed by i ∈
I = {1, . . . , q}. This allows us to use the standard definitions of the RKHS norm and GP model.
T
In this case,
P we define the GP posterior distribution as µn (z) = [µn (z, 1), . . . , µn (z, q)] and
σn (z) = 1≤i≤q σn (z, i), where the unusual definition of the standard deviation is used in Lemma 3.
Given the previous assumptions, it follows from [28, Lemma 2] that the dynamics in (1) are Lipschitz
continuous with Lipschitz constant Lf = Lh + Lg , where Lg depends on the properties (smoothness)
of the kernel.
Moreover, we can construct high-probability confidence intervals on the dynamics in (1) that fulfill Assumption 2 using the GP model.
Lemma 3. ([35, Theorem 6]) Assume σ-sub-Gaussian noise andpthat the model error g(·)
in (1) has RKHS norm bounded by Bg . Choose βn = Bg + 4σ γn + 1 + ln(1/δ). Then,
with probability at least 1 − δ, δ ∈ (0, 1), for all n ≥ 1, x ∈ X , and u ∈ U it holds that
kf (x, u) − µn−1 (x, u)k1 ≤ βn σn−1 (x, u).
Proof. From [34, Theorem 2] it follows that |f (x, u, i) − µn−1 (x, u, i)| ≤ βn σn (x, u, i) holds with
probability at least 1 − δ for all 1 ≤ i ≤ q. Following Remark 1, we can model the multi-output
function as a single-output function over an extended parameter space. Thus the result directly
transfers by definition of the one norm and our definition of σn for multiple output dimensions
13

in Remark 1. Note that by iteration n we have obtained nq measurements in the information
capacity γn .
That is, the true dynamics are contained within the GP posterior confidence intervals with high
probability. The bound depends on the information capacity,
γn =

max

A⊂X ×U ×I : |A|=nq

(9)

I(yA ; fA ),

which is the maximum mutual information that could be gained about the dynamics f from samples.
The information capacity has a sublinear dependence on n(6= t) for many commonly used kernels
such as the linear, squared exponential, and Matérn kernels and it can be efficiently and accurately
approximated [35]. Note that we explicitly account for the q measurements that we get for each of
the q states in (9).
Remark 2. The GP model assumes Gaussian noise, while Lemma 3 considers σ-sub-Gaussian noise.
Moreover, we consider functions with bounded RKHS norm, rather than samples from a GP. Lemma 3
thus states that even though we make different assumptions than the model, the confidence intervals
are conservative enough to capture the true function with high probability.
A.3

Safe exploration

Remark 3. In the following we assume that Dn and Sn are defined as in (4) and (5).
Baseline As a baseline, we consider a class of algorithms that know about the Lipschitz continuity
properties of v, f , and π. In addition, we can learn about v(f (x, u)) up to some arbitrary statistical
accuracy  by visiting state x and obtaining a measurement for the next state after applying action u,
but face the safety restrictions defined in Sec. 2. Suppose we are given a set S of state-action pairs
about which we can learn safely. Specifically, this means that we have a policy such that, for any
state-action pair (x, u) in S, if we apply action u in state x and then apply actions according to the
policy, the state converges to the origin. Such a set can be constructed using the initial policy π0
from Sec. 2 as S0 = {(x, π0 (x)) | x ∈ S0x }.

The goal of the algorithm is to expand this set of states that we can learn about safely. Thus, we need
to estimate the region of attraction by certifying that state-action pairs achieve the −L∆v τ decrease
condition in Theorem 2 by learning about state-action pairs in S. We can then generalize the gained
knowledge to unseen states by exploiting the Lipschitz continuity,

Rdec (S) = S0 ∪ z ∈ Xτ × Uτ | ∃(x, u) ∈ S : v(f (x, u))−v(x)++L∆v kz−(x, u)k1 <−L∆v τ ,
(10)
where we use that we can learn v(f (x, u)) up to  accuracy within S. We specifically include S0 in
this set, to allow for initial policies that are safe, but does not meet the strict decrease requirements
of Theorem 2. Given that all states in Rdec (S) fulfill the requirements of Theorem 2, we can estimate
the corresponding region of attraction by committing to a control policy π ∈ ΠL and estimating the
largest safe level set of the Lyapunov function. With D = Rdec (S), the operator

Rlev (D) = V argmax c, such that ∃π ∈ ΠL : ∀x ∈ V(c) ∩ Xτ , (x, π(x)) ∈ D
(11)

encodes this operation. It optimizes over safe policies π ∈ ΠL to determine the largest level set,
such that all state-action pairs (x, π(x)) at discrete states x in the level set V(c) ∩ Xτ fulfill the
decrease condition of Theorem 2. As a result, Rlev (Rdec (S)) is an estimate of the largest region
of attraction given the -accurate knowledge about state-action pairs in S. Based on this increased
region of attraction, there are more states that we can safely learn about. Specifically, we again use
the Lipschitz constant and statistical accuracy  to determine all states that map back into the region
of attraction,

R (S) = S∪ z0 ∈ Rτlev (Rdec (S)) × Uτ | ∃z ∈ S : v(f (z))++Lv Lf kz−z0 k1 ≤ max v(x) ,
x∈Rlev (Rdec (S))

(12)

where Rτlev (D) = Rlev (D) ∩ Xτ . Thus, R (S) ⊇ S contains state-action pairs that we can visit to
learn about the system. Repeatedly applying this operator leads the largest set of state-action pairs
that any safe algorithm with the same knowledge and restricted to policies in ΠL could hope to reach.
14

Algorithm 2 Theoretical algorithm
1: Input: Initial safe policy S0 , dynamics model GP(µ(z), k(z, z0 ))
2: for all n S
= 1, . . . do 
3:
Dn = (x,u)∈Sn−1 z0 ∈ Xτ × Uτ | un (x, u) − v(x) + L∆v kz0 − (x, u)k1 < −L∆v τ ,
4:
πn , cn = argmaxπ∈ΠL ,c∈R>0 c, such that for all x ∈ V(c) ∩ Xτ : (x, π(x)) ∈ Dn

S
5:
Sn = z∈Sn−1 z0 ∈ V(cn ) ∩ Xτ × Uτ | un (z) + Lv Lf kz − z0 k1 ≤ cn }
 0
S
6:
= z∈Sn−1 z ∈ Rτlev (Dn ) × Uτ | un (z) + Lv Lf kz − z0 k1 ≤ maxx∈Rlev (Dn ) v(x)}
7:
(xn , un ) = argmax(x,u)∈Sn un (x, u) − ln (x, u)
8:
Sn = {z ∈ V(cn ) × Uτ | un (z) ≤ cn }
9:
Update GP with measurements f (xn , un ) + n
Specifically, let R0 (S) = S and Ri+1 (S) = R (Ri (S)). Then R (S) = limi→∞ Ri (S) is the set
of all state-action pars on the discrete grid that any algorithm could hope to classify as safe without
leaving this safe set. Moreover, Rlev (R (S)) is the largest corresponding region of attraction that
any algorithm can classify as safe for the given Lyapunov function.
Proofs In the following we implicitly assume that the assumptions of Lemma 3 hold and that βn is
defined as specified within Lemma 3. Moreover, for ease of notation we assume that S0x is a level set
of the Lyapunov function v(·).
Lemma 4. V(cn ) = Rlev (Dn ) and cn = maxx∈Rlev (Dn ) v(x)
Proof. Directly by definition, compare (3) and (11).
Remark 4. Lemma 4 allows us to write the proofs entirely in terms of operators, rather than having
to deal with explicit policies. In the following and in Algorithm 2 we replace V(cn ) and cn according
to Lemma 4. This moves the definitions closer to the baseline and makes for an easier comparison.
We roughly follow the proof strategy in [18], but deal with the additional complexity of having safe
sets that are defined in a more difficult way (indirectly through the policy). This is non-trivial and
the safe sets are carefully designed in order to ensure that the algorithm works for general nonlinear
systems.
We start by listing some fundamental properties of the sets that we defined below.
Lemma 5. It holds for all n ≥ 1 that
(i) ∀z ∈ Xτ × Uτ , un+1 (z) ≤ un (z)

(ii) ∀z ∈ Xτ × Uτ , ln+1 (z) ≥ ln (z)

(iii) S ⊆ R =⇒ Rlev (S) ⊆ Rlev (R)
(iv) S ⊆ R =⇒ Rdec (S) ⊆ Rdec (R)
(v) S ⊆ R =⇒ R (S) ⊆ R (R)

(vi) S ⊆ R =⇒ R (S) ⊆ R (R)

(vii) Sn ⊇ Sn−1 =⇒ Dn+1 ⊇ Dn

(viii) D1 ⊇ S0

(ix) Sn ⊇ Sn−1

(x) Dn ⊇ Dn−1

Proof. (i) and (ii) follow directly form the definition of Cn .
(iii) Let π ∈ ΠL be a policy such that for some c > 0 it holds for all x ∈ V(c) ∩ Xτ that
(x, π(x)) ∈ S. Then we have that (x, π(x)) ∈ R, since S ⊆ R. Thus it follows that with
cs = argmax c s.t. ∃π ∈ ΠL : ∀x ∈ V(c) ∩ Xτ , (x, π(x)) ∈ S
(13)
15

and
cr = argmax c s.t. ∃π ∈ ΠL : ∀x ∈ V(c) ∩ Xτ , (x, π(x)) ∈ R
we have that cr ≥ cs . This implies V(cr ) ⊇ V(cs ). The result follows.

(14)

(iv) Let z ∈ Rdec (S). Then there exists (x, u) ∈ S such that v(f (x, u)) − v(x) +  + L∆v kz −
(x, u)k1 < −L∆v τ . Since S ⊆ R we have that (x, u) ∈ R as well and thus z ∈ Rdec (R).
(v) S ⊆ R =⇒ Rlev (Rdec (S)) ⊆ Rlev (Rdec (R)) due to (iii) and (iv). Since z0 ∈ R (S), there
must exist an z ∈ S such that v(f (z)) +  + Lv Lf kz − z0 k1 ≤ maxx∈Rlev (Rdec (S)) v(x).
Since S ⊆ R it follows that z ∈ R. Moreover,
max
v(x) ≤
max
v(x)
(15)
x∈Rlev (Rdec (S))

x∈Rlev (Rdec (R))

follows from Rlev (Rdec (S)) ⊆ Rlev (Rdec (R)), so that we conclude that z0 ∈ R (R).

(vi) This follows directly by repeatedly applying the result of (v).

(vii) Let z0 ∈ Dn . Then ∃(x, u) ∈ Sn−1 : un (x, u) − v(x) + L∆v kz0 − (x, u)k1 < −L∆v τ . Since
Sn ⊇ Sn−1 it follows that (x, u) ∈ Sn as well. Moreover, we have
un+1 (x, u) − v(x) + L∆v kz0 − (x, u)k1
≤ un (x, u) − v(x) + L∆v kz0 − (x, u)k1 < −L∆v τ
since un+1 is non-increasing, see (i). Thus z0 ∈ Dn+1 .

(viii) By definition of C0 we have for all (x, u) ∈ S0 that u0 (x, u) < v(x) − L∆v τ . Now we have
that
u1 (x, u) − v(x) + L∆v k(x, u) − (x, u)k1 ,
= u1 (x, u) − v(x),
≤ u0 (x, u) − v(x), by Lemma 5 (i)
< − L∆v τ,
which implies that (x, u) ∈ D1 .

(ix) Proof by induction. We consider the base case, z ∈ S0 , which implies that z ∈ D1 by (viii).
Moreover, since S0x is a level set of the Lyapunov function v by assumption, we have that
Rlev (S0 ) = S0x . The previous statements together with (iii) imply that z ∈ Rτlev (D1 ) × Uτ ,
since D1 ⊇ S0 by (viii). Now, we have that
(i)

u1 (z) + Lv Lf kz − zk1 = u1 (z) ≤ u0 (z).
Moreover, by definition of C0 , we have for all (x, u) ∈ S0 that
u0 (x, u) < v(x) − L∆v τ < v(x).
As a consequence,
u0 (x, u) ≤ max v(x),
(x,u)∈S0

=
≤

(16)

max

v(x),

(17)

max

v(x),

(18)

x∈Rlev (S0 )
x∈Rlev (D1 )

where the last inequality follows from (iii) and (viii). Thus we have z ∈ S1 .

For the induction step, assume that for n ≥ 2 we have z0 ∈ Sn with Sn ⊇ Sn−1 . Now since
z0 ∈ Sn we must have that z0 ∈ Rτlev (Dn ) × Uτ . This implies that z0 ∈ Rτlev (Dn+1 ) × Uτ , due
to Lemma 5 (iii) and (vii) together with the induction assumption of Sn ⊇ Sn−1 . Moreover,
there must exist a z ∈ Sn−1 ⊆ Sn such that
un+1 (z) + Lv Lf kz − z0 k1 , ≤ un (z) + Lv Lf kz − z0 k1 ,
(19)
≤ max v(x),
(20)
x∈Rlev (Dn )

≤

max

x∈Rlev (Dn+1 )

v(x),

(21)

which in turn implies z ∈ Sn+1 . The last inequality follows from Lemma 5 (iii) and (vii)
together with the induction assumption that Sn ⊇ Sn−1 .
16

(x) This is a direct consequence of (vii), (viii), and (ix).

Given these set properties, we first consider what happens if the safe set Sn does not expand after
collecting data points. We use these results later to conclude that the safe set must either expand or
that the maximum level set is reached. We denote by
zn = (xn , un )
the data point the is sampled according to (6).
Lemma 6. For any n1 ≥ n0 ≥ 1, if Sn1 = Sn0 , then for any n such that n0 ≤ n < n1 , it holds that
s
C1 qβn2 γn
2βn σn (zn ) ≤
,
(22)
n − n0
where C1 = 8/ log(1 + σ −2 ).
Proof. We modify the results for q = 1 by [35] to this lemma, but use the different definition for βn
from [34]. Even though the goal of [35, Lemma 5.4] is different from ours, we can still apply their
reasoning to bound the amplitude of the confidence interval of the dynamics. In particular, in [35,
Lemma 5.4], we have rn = 2βn σn−1 (zn ) with zn = (xn , un ) according to Lemma 3. Then
(23)

2
rn2 = 4βn2 σn−1
(zn ),

= 4βn2
≤ 4βn2 q

q
X
i=1
q
X

!2
σn−1 (zn , i)
2
σn−1
(zn , i)

(24)

,
(Jensen’s ineq.),

(25)

i=1

≤ 4βn2 qσ 2 C2

q
X

2
log(1 + σ −2 σn−1
(zn , i)),

(26)

i=1

where C2 = σ −2 / log(1 + σ −2 ). The result then follows analogously to [35, Lemma 5.4] by noting
that
n
X
rj2 ≤ C1 βn2 qγn
∀n ≥ 1
(27)
j=1

according to the definition of γn in this paper and using the Cauchy-Schwarz inequality.
The previous result allows us to bound the width of the confidence intervals:
Corollary 2. For any n1 ≥ n0 ≥ 1, if Sn1 = Sn0 , then for any n such that n0 ≤ n < n1 , it holds
that
s
C1 qβn2 γn
un (zn ) − ln (zn ) ≤ Lv
,
(28)
n − n0
where C1 = 8/ log(1 + σ −2 ).

Proof. Direct consequence of Lemma 6 together with the definition of C and Q.

Corollary 3. For any n ≥ 1 with C1 as defined in Lemma 6, let Nn be the smallest integer satisfying
C L2 q
Nn
≥ 12 v and Sn+Nn = SNn , then, for any z ∈ Sn+Nn it holds that
β2
γn+Nn
n+Nn

un (z) − ln (z) ≤ .
Proof. The result trivially follows from substituting Nn in the bound in Corollary 2.
17

(29)

Lemma 7. For any n ≥ 1, if R (S0 ) \ Sn 6= ∅, then R (Sn ) \ Sn 6= ∅.
Proof. As in [18, Lemma 6]. Assume, to the contrary, that R (Sn ) \ Sn = ∅. By definition R (Sn ) ⊇ Sn , therefore R (Sn ) = Sn . Iteratively applying R to both sides, we get in the
limit R (Sn ) = Sn . But then, by Lemma 5,(vi) and (ix), we get
(30)

R (S0 ) ⊆ R (Sn ) = Sn ,
which contradicts the assumption that R (S0 ) \ Sn 6= ∅.

Lemma 8. For any n ≥ 1, if R (S0 ) \ Sn 6= ∅, then the following holds with probability at least
1 − δ:
Sn+Nn ⊃ Sn .
(31)
Proof. By Lemma 7, we have that R (Sn ) \ Sn 6= ∅. By definition, this means that there exist
z ∈ R (Sn ) \ Sn and z0 ∈ Sn such that
v(f (z0 )) +  + Lv Lf kz − z0 k1 ≤

max

x∈Rlev (Rdec (Sn ))

v(x)

(32)

Now we assume, to the contrary, that Sn+Nn = Sn (the safe set cannot decrease due to Lemma 5
(ix)). This implies that z ∈ Xτ × Uτ \ Sn+Nn and z0 ∈ Sn+Nn = Sn+Nn −1 . Due to Corollary 2, it
follows that
un+Nn (z0 ) + Lv Lf kz − z0 k1

(33)

≤ v(f (z0 )) +  + Lv Lf kz − z0 k1
≤
max
v(x)
x∈Rlev (Rdec (Sn ))

=

max

x∈Rlev (Rdec (Sn+Nn ))

v(x)

by (32)

(34)
(35)

by (iii), (iv) and (ix)

(36)

Thus, to conclude that z ∈ Sn+Nn according to (5), we need to show that
Rlev (Dn+Nn ) ⊇ Rlev (Rdec (Sn )).
To this end, we use Lemma 5 (iii) and show that
Dn+Nn ⊇ Rdec (Sn+Nn ). Consider (x, u) ∈ Rdec (Sn+Nn ), we know that there exists a (x0 , u0 ) ∈
Sn+Nn = Sn+Nn −1 such that
−L∆v τ > v(f (x0 , u0 )) − v(x0 ) +  + L∆v k(x, u) − (x0 , u0 )k1 ,
0

0

0

0

0

≥ un+Nn (x , u ) − v(x ) + L∆v k(x, u) − (x , u )k1 ,

(37)
(38)

where the second inequality follows from Corollary 2. This implies that (x, u) ∈ Dn and thus
Dn+Nn ⊇ Rdec (Sn+Nn ). This, in turn, implies that z ∈ Sn+Nn , which is a contradiction.
Lemma 9. For any n ≥ 0, the following holds with probability at least 1 − δ:
Sn ⊆ R0 (S0 ).

(39)

Proof. Proof by induction. For the base case, n = 0, we have S0 ⊆ R0 (S0 ) by definition.

For the induction step, assume that for some n ≥ 1, Sn−1 ⊆ R0 (S0 ). Let z ∈ Sn . Then, by
definition, ∃z0 ∈ Sn−1 such that
un (z0 ) + Lv Lf kz − z0 k1 ≤

max

x∈Rlev (Dn )

v(x),

(40)

v(x)

(41)

which, by Corollary 1, implies that
v(f (z0 )) + Lv Lf kz − z0 k1 ≤

max

x∈Rlev (Dn )

Now since z0 ∈ R0 (S0 ) by the induction hypothesis, in order to conclude that z ∈ R0 (S0 ) we need
to show that Rlev (Dn ) ⊆ Rlev (Rdec (R(S0 ))) .
18

Let (x, u) ∈ Dn , then there exist (x0 , z0 ) ∈ Sn−1 such that

un−1 (x0 , u0 ) − v(x0 ) + L∆v k(x, u) − (x0 , u0 )k1 < −L∆v τ,

(42)

v(f (x0 , u0 )) − v(x0 ) + L∆v k(x, u) − (x0 , u0 )k1 < −L∆v τ,

(43)

which, by Corollary 1, implies that

which means that (x, u) ∈ Rdec (R0 (S0 )) since Sn−1 ⊆ R0 (S0 ) and therefore (x0 , u0 ) ∈ R0 (S0 )
holds by the induction hypothesis. We use (iii) to conclude that Rlev (Dn ) ⊆ Rlev (Rdec (R(S0 ))),
which concludes the proof.
Lemma 10. Let n∗ be the smallest integer, such that n∗ ≥ |R0 (S0 )|Nn∗ . Then, there exists n0 ≤ n∗
such that Sn0 +Nn0 = Sn0 holds with probability 1 − δ.
Proof. By contradiction. Assume, to the contrary, that for all n ≤ n∗ , Sn ⊂ Sn+Nn . From Lemma 5
(ix) we know that Sn ⊆ Sn+Nn . Since Nn is increasing in n, we have that Nn ≤ Nn∗ . Thus, we
must have
S0 ⊂ SNn∗ ⊂ S2Nn∗ · · · ,
(44)
so that for any 0 ≤ j ≤ |R0 (S0 )|, it holds that |SjTn∗ | > j. In particular, for j = |R0 (S0 )|, we get
|SjNn∗ | > |R0 (S0 )|,

(45)

which contradicts SjNn∗ ⊆ R0 (S0 ) from Lemma 9.

Corollary 4. Let n∗ be the smallest integer such that
C1 L2v q|R0 (S0 )|
n∗
≥
,
β n∗ γ n∗
2

(46)

then there exists a n0 ≤ n∗ such that Sn0 +Nn0 = Sn0 .
Proof. A direct consequence of Lemma 10 and Corollary 3.
A.4

Safety and policy adaptation

In the following, we denote the true region of attraction of (1) under a policy π by Rπ .
Lemma 11. Rlev (Dn ) ⊆ Rπn for all n ≥ 0.
Proof. By definition, we have for all (x, u) ∈ Dn that the exists (x0 , u0 ) ∈ Sn−1 such that
−L∆v τ ≥ un (x0 , u0 ) − v(x0 ) + L∆v k(x, u) − (x0 , u0 )k1 ,

≥ v(f (x0 , u0 )) − v(x0 ) + L∆v k(x, u) − (x0 , u0 )k1 ,
≥ v(f (x, u)) − v(x),

where the first inequality follows from Corollary 1 and the second one by Lipschitz continuity,
see Lemma 1.
By definition of Rlev in (11), it follows that for all x ∈ Rlev (Dn ) ∩ Xτ we have that (x, πn (x)) ∈ Dn .
Moreover, Rlev (Dn ) is a level set of the Lyapunov function by definition. Thus the result follows
from Theorem 2.
Lemma 12. f (x, u) ∈ Rπn ∀(x, u) ∈ Sn .
Proof. This holds for S0 by definition. For n ≥ 1, by defition, we have for all z ∈ Sn there exists an
z0 ∈ Sn−1 such that
max

x∈Rlev (Dn )

v(x) ≥ un (z0 ) + Lv Lf kz − z0 k1

≥ v(f (z0 )) + Lv Lf kz − z0 k1
≥ v(f (z))

where the first inequality follows from Corollary 1 and the second one by Lipschitz continuity,
see Lemma 1. Since Rlev (Dn ) ⊆ Rπn by Lemma 11, we have that f (z) ∈ Rπn .
19

Theorem 4. Assume σ-sub-Gaussian measurement noise and that the model error g(·)
in (1) has RKHSpnorm smaller than Bg .
Under the assumptions of Theorem 2,
with βn = Bg + 4σ γn + 1 + ln(1/δ), and with measurements collected according to (6), let
∗

0 )|+1)
n∗ be the smallest positive integer so that β 2 nγ ∗ ≥ Cq(|R(S
where C = 8/ log(1 + σ −2 ).
L2v 2
n∗ n
Let Rπ be the true region of attraction of (1) under a policy π. For any  > 0, and δ ∈ (0, 1), the
following holds jointly with probability at least (1 − δ) for all n > 0:

(i) V(cn ) ⊆ Rπn

(ii) f (x, u) ∈ Rπn ∀(x, u) ∈ Sn .

(iii) R (S0 ) ⊆ Sn ⊆ R0 (S0 ).

Proof. See Lemmas 11 and 12 for (i) and (ii), respectively. Part (iii) is a direct consequence
of Corollary 4 and Lemma 9.

20

