Safe Learning of Quadrotor Dynamics Using Barrier Certificates *

arXiv:1710.05472v1 [cs.LG] 16 Oct 2017

Li Wang, Evangelos A. Theodorou, and Magnus Egerstedt†

Abstract— To effectively control complex dynamical systems,
accurate nonlinear models are typically needed. However, these
models are not always known. In this paper, we present a datadriven approach based on Gaussian processes that learns models of quadrotors operating in partially unknown environments.
What makes this challenging is that if the learning process
is not carefully controlled, the system will go unstable, i.e.,
the quadcopter will crash. To this end, barrier certificates are
employed for safe learning. The barrier certificates establish
a non-conservative forward invariant safe region, in which
high probability safety guarantees are provided based on the
statistics of the Gaussian Process. A learning controller is
designed to efficiently explore those uncertain states and expand
the barrier certified safe region based on an adaptive sampling
scheme. In addition, a recursive Gaussian Process prediction
method is developed to learn the complex quadrotor dynamics
in real-time. Simulation results are provided to demonstrate the
effectiveness of the proposed approach.

I. I NTRODUCTION
Safety is crucial to many physical control dynamical
systems, such as autonomous vehicles, industrial robots,
chemical reactors, and air-traffic control systems [11], [2],
[5]. If the system reaches certain unsafe states or even fails,
both the operator and the controlled plant might be put in
serious danger. The existence of model inaccuracies and
unknown disturbances create an even greater challenge to
the design of safe controllers for these systems.
Tools such as robust control and adaptive control methods
have been developed in classic control theory to ensure the
safety and stability of the system, see [9], [4] and the references therein. Meanwhile, machine learning based control
approaches are becoming increasingly popular as a way to
deal with inaccurate models [8], [17], due to their abilities
to infer unknown models from data and actively improve
the performance of the controller with the learned model. In
contrast to classic control methods, learning based control
approaches require only limited expert knowledge and fewer
assumptions about the system [22]. However, there always
exists an inherent trade-off between safety and performance
in these methods [3]. Data-driven learning approaches rarely
provides safety guarantees, which limits their applicability
to real-world safety critical control dynamical systems [22].
The objective of this paper is to construct high probability
safety guarantees for Gaussian Process (GP) based learning
approaches using barrier certificates.
*The work by the first and third authors was sponsored by Grant No.
1544332 from the U.S. National Science Foundation.
† Li Wang and Magnus Egerstedt are with the School of Electrical and Computer Engineering, Evangelos A. Theodorou is with
the School of Aerospace Engineering, Georgia Institute of Technology, Atlanta, GA 30332, USA. Email: {liwang, magnus,

evangelos.theodorou}@gatech.edu

In order to promote the application of learning based
control methods in safety-critical systems, a number of safe
learning approaches have been proposed in the literature.
Among these methods, the use of learning Control Lyapunov
Functions (CLF) is shown to be a promising approach. A
learning from demonstration method was developed in [13]
to search for a CLF from several demonstrations, and the
learned CLF was used to stabilize the system. But the learned
controller did not consider actuator limits and other safe
operation constraints. [21] introduced a verifier to explicitly
validate the learned CLF. However, when the model of the
system is inaccurate, the verifier needs to check an infinite
number of inequalities throughout the state space, which is
computationally difficult [12]. [6] seeks to learn CLF and
maximize the safe operation region for the system with
GP model. High probability safety guarantees are provided
based on Lyapunov stability and GP statistics. In addition, a
reachability-based safe learning approach was presented in
[1] to reduce the conservativeness of reachability analysis by
learning the disturbance from data.
In contrast to the aforementioned methods, this paper
interprets the safe operation region as general invariant sets
established with barrier certificates, which permits a much
richer set of safe control options, rather than Lyapunov sublevel sets. The barrier certificates formally define a forward
invariant safe region, where all system trajectories starting
in this region remains in this region for all time [18], [26],
[2]. With the barrier certificates, the safety of the system
can be certified without explicitly computing the forward
reachable set [19]. Barrier certificates were successfully
applied to many safety critical dynamical systems, such as
adaptive cruise control [2], bipedal walking [15], quadrotor
control [25], and swarm robotics [23], [24]. In this paper, we
construct a safe operation region with barrier certificates, and
gradually expand the certified safe region as the uncertainty
of the system reduces. The unknown dynamics of the system
is represented with a GP model, which provides both the
mean and variance of the prediction. Using the statistics of
GP model, a high probability safety guarantee of the system
with inaccurate model is provided.
The search for maximum volume barrier certificates involves the validation of an infinite number of inequality
constraints, which is computationally expensive. Inspired
by the discrete sampling technique used in [5], we design
an adaptive sampling algorithm to significantly reduce the
computation intensity, i.e., the more certain regions in the
state space are sampled less without loss of safety guarantees.
In addition, a recursive learning strategy based on GP is
designed to learn the complex 3D nonlinear quadrotor dy-

namics online. The learned dynamical model of the quadrotor
is then incorporated into a differential flatness based flight
controller to improve the trajectory tracking performance.
The main contributions of this paper are threefold. First,
a safe learning strategy is developed based on barrier certificates, which admits a rich set of learning control options.
Second, an adaptive sampling algorithm is proposed to
significantly reduce the computation intensity of the learning
process. Third, an recursive learning strategy based on GP
is presented to learn the complex 3D nonlinear quadrotor
dynamics online.
The rest of this paper is organized as follows. The preliminaries of barrier certificates and GP are briefly revisited
in Section II. A safe learning strategy based on barrier
certificates is presented in Section III. Section IV contains
a real-time learning algorithm for 3D quadrotor dynamics
based on GP. Simulation results are provided in Section V,
and the paper is ended by conclusions in Section VI.
II. P RELIMINARIES OF BARRIER C ERTIFICATES AND
G AUSSIAN P ROCESS
Preliminary results regarding the two fundamental tools,
i.e., barrier certificates and Gaussian Process, used to formulate the safe learning strategy are presented in this section.

certificates provide a more permissive notion of safety. As a
result, barrier certificates based safe learning controllers have
more freedom to efficiently explore those unknown states.
This fact can be illustrated with the following example.
Example 1: Consider an autonomous dynamical system

  
x2 + 0.8x22
ẋ1
=
,
(3)
ẋ2
−x1 − x2 + x21 x2
the safe region of this system is estimated with both the
Lyapunov sublevel set and barrier certificates.
Since (3) is a polynominal system, the safe sets can
be computed directly with Sum-of-Squares programs using YALMIP [14] and SMRSOFT [7] solvers. Both the
Lyapunov function and barrier certificates are limited to
second order polynomials. The safe region estimated with
the optimal polynomial Lyapunov function is
A1 = {x | V ∗ (x) ≤ 1},
where V ∗ (x) = 1.343x21 + 0.5155x1x2 + 1.152x22.
The safe region estimated with barrier certificates is
A2 = {x | h∗ (x) ≥ 0},
where h∗ (x) = 1 − 0.4254x1 − 0.3248x2 − 0.7549x22 −
0.8616x21 − 0.2846x1x2 .

A. Barrier Certificates and Set Invariance
Consider a control affine dynamical system

where x ∈ X ⊆ Rn and u ∈ U ⊆ Rm are the state and control
of the system, f : Rn → Rn and g : Rn → Rm are Lipschitz
continuous. Let the safe set of the system be encoded as the
superlevel set of a smooth function h : Rn → R,
C = {x ∈ Rn | h(x) ≥ 0}.

(2)

The function h(x) is termed a Control Barrier Function
(CBF), if there exists an extended class-κ function (κ (0) = 0
and strictly increasing) such that


∂h
∂h
f (x) + g(x)u + κ (h(x)) ≥ 0,
sup
∂x
∂x
u∈U

A2 (h∗ (x) = 0)
A1 (V ∗ (x) = 1)

2

(1)

1

x2

ẋ = f (x) + g(x)u,

0
-1
-2
-2

-1

0

1

2

x1

Fig. 1: Estimates of safe regions for system (3). The regions
enclosed by the dashed red ellipse and solid green ellipse are
estimated safe regions with optimal polynomial Lyapunov
function V ∗ (x) and barrier certificates h∗ (x), respectively.

for all x ∈ E with C ⊆ E .
Given a CBF, the barrier certified safe control space S(x)
is defined as


∂h
∂h
From Fig. 1, it can be observed that the barrier certified
f (x) + g(x)u + κ (h(x)) ≥ 0 , x ∈ E .
S(x) = u ∈ U |
∂x
∂x
safe region A2 is much larger than the Lyapunov based safe
With barrier certificates, the invariance property of C is region A1 . Consequently, safe learning controller based on
established with the following theorem,
barrier certificates are allowed to explore more states of the
Theorem [26]: Given a set C ⊂ Rn defined by (2) and a CBF system. In this paper, we will leverage the non-conservative
h defined on E , with C ⊆ E ⊂ Rn , any Lipschitz continuous safety guarantee of barrier certificates to allow a much richer
controller u : E → R such that u ∈ S(x) for the system (1) set of safe learning control options.
renders the set C forward invariant.
This type of barrier certificates expands the certified safe B. Gaussian Processes
A GP is a nonparametric regression method that can capcontrol space significantly by allowing h(x) to decrease
within C as opposed to strictly increasing [26], [2]. Com- ture complex unknown functions [20]. With a GP, every point
pared with Lyapunov sublevel set based safe region, barrier in the state space is associated with a normally distributed

random variable, which allows us to derive high probability
statements about the system.
Adding some unknown dynamics d(x) to the original class
of control-affine systems (1), we now consider a system with
partially unknown dynamics in this paper, i.e.,
ẋ = f (x) + g(x)u + d(x),

(4)

where x ∈ X ⊆ Rn and u ∈ U ⊆ Rm are the state and control
of the system. Although the proposed method applies to
general dynamical systems, here we restrict our attention
to the class of systems that can be addressed with existing
computation tools. It is also assumed that d(x) is Lipschitz
continuous. This assumption is necessary, because we want
to generalize the learned dynamics to states that are not
explored before.
Since the unmodeled dynamics d(x) is n dimensional, each
dimension is approximated with a GP model G P(0, k(x, x′ ))
with a prior mean of zero and a covariance function of
k(x, x′ ), where k(x, x′ ) is the kernel function to measure the
similarity between any two states x, x′ ∈ X . In order to make
GP inferences on the unknown dynamics, we need to get
ˆ
measurements of d(x). This measurement d(x)
is obtained
indirectly by subtracting the inaccurate model prediction
[ f (x) + g(x)u] from the noisy measurement of the system
dynamics [ẋ + N (0, σn2 )]. Since any finite number of data
points form a multivariate normal distribution, we can obtain
the posterior distribution of d(x∗ ) at any query state x∗ ∈ X
by conditioning on the past measurements [20].
Given a collection of w measurements yw =
ˆ 1 ), d(x
ˆ 2 ), ..., d(x
ˆ w )]T , the mean m(x∗ ) and variance
[d(x
2
σ (x∗ ) of d(x∗ ) at the query state x∗ are
m(x∗ ) =
σ 2 (x∗ ) =

k∗T (K + σn2 I)−1 yw ,
k(x∗ , x∗ ) − k∗T (K + σn2 I)−1 k∗ ,

(5)
(6)

where ⌊K⌋(i, j) = k(xi , x j ) is the kernel matrix, and k∗ =
[k(x1 , x∗ ), k(x2 , x∗ ), ..., k(xw , x∗ )]T .
With the learned system dynamics based on GP, a high
probability confidence interval of the unmodeled dynamics
d(x) can be established as
D(x) = {d | m(x) − kδ σ (x) ≤ d ≤ m(x) + kδ σ (x)},

(7)

where kδ is a design parameter to get (1 − δ ) confidence,
δ ∈ (0, 1). For instance, 95.5% and 99.7% confidence are
achieved at kδ = 2 and kδ = 3, respectively.
III. S AFE L EARNING

WITH

BARRIER C ERTIFICATES

In order to ensure that the learning based controller never
enters the unsafe region, we will learn barrier certificates for
the system and use the learned certificates to regulate the
controller. As discussed in Section II, the barrier certificates
certify a safe region that is forward invariant. We can first
start with an conservative barrier certificate with certified
safe region C0 (x), then gradually expand this certified safe
region with the collected data until it stops growing. This
incremental learning process is visualized in Fig. 2.

Fig. 2: Incremental learning of the barrier certificates. The
green region C0 and the yellow regions Cn are the initial and
final barrier certified safe regions, respectively. The barrier
certified safe region gradually grows as more and more data
points are sampled in the state space.

More concretely, the goal of the learning process is to
maximize the volume of the barrier certified safe region C
by adjusting h(x), i.e.,

s.t.

max vol(C )
h(x)


∂h
( f (x) + g(x)u + d) + γ h(x) ≥ 0,
max min
u∈U d∈D (x) ∂ x
∀x ∈ C .

Since u and d are independent from each other, we can
rewrite this optimization problem into
max vol(C )
h(x)




∂ hk
∂h
g(x)u + min
d
s.t. max
(8)
∂x
u∈U
d∈D (x) ∂ x
∂h
+
f (x) + γ h(x) ≥ 0, ∀x ∈ C
∂x
Using the high confidence interval D(x) in (7), the barrier
certificates constraint can be considered as
max vol(C )
h(x)


∂h
∂h
∂h
s.t. max
g(x)u + m(x) − kδ
σ (x)
(9)
u∈U
∂x
∂x
∂x
∂h
f (x) + γ h(x) ≥ 0, ∀x ∈ C .
+
∂x
When more data points are collected about the system
dynamics, the uncertainty σ (x) will gradually decrease. As
a result, more states will satisfy the barrier certificates
constraint. The goal of the exploration task is to actively
collect data to reduce σ (x) and maximize the volume of C .
It should be pointed out that the barrier certified region
maximization problem (9) is a non-convex, infinite dimensional optimization problem, which is intractable to solve
in practice. We will make two simplifications to make it
solvable, namely by employing adaptive sampling of the state
space and parameterization of the shape of C .
A. Adaptive Sampling of the State Space
Due to the Lipschitz continuity of the system dynamics,
the safety of the system in X can be evaluated by only

sampling a finite number of points in X . Inspired by [5],
we will show that we can adaptively sample the state space
without losing safety guarantees. Similar to Lemma 4 in [5],
it can be shown that h(x) and ḣ(x) are Lipschitz continuous
in x with Lipschitz constants Lh and Lḣ , respectively.
Let Xτ ⊂ X be a discretization of the state space X .
The closest point in Xτ to x ∈ X is denoted as [x]τ , where
kx − [x]τ k ≤ τ2 .
Lemma 3.1: If the following condition holds for all x ∈
Xτ ,


∂h
∂h
∂h
max
g(x)u + m(x) − kδ
σ (x)
u∈U
∂x
∂x
∂x
∂h
+ f (x) + γ h(x) ≥ (Lḣ + γ Lh )τ ,
(10)
∂x
then the safety barrier constraint


∂h
( f (x) + g(x)u + d) + γ h(x) ≥ 0
(11)
max min
u∈U d∈D (x) ∂ x
is satisfied for all x ∈ X with probability (1 − δ ), δ ∈ (0, 1).
Proof: With the definition of the high confidence
interval D(x), (10) can be rewritten as


∂h
( f (x) + g(x)u + d) + γ h(x) ≥ (Lḣ + γ Lh )τ ,
max min
u∈U d∈D (x) ∂ x
with a probability of (1 − δ ), for all x ∈ Xτ . This is
equivalent to
ḣ(x) + γ h(x) ≥ (Lḣ + γ Lh )τ ,
for all x ∈ Xτ .
Because of the Lipschitz continuity of h(x) and ḣ(x), we
have for any x ∈ X ,
ḣ(x) + γ h(x) ≥
≥

(ḣ([x]τ ) − Lḣ τ ) + γ (h([x]τ ) − Lh τ )
0.

This means that the safety barrier constraint is satisfied for
any x ∈ X , if (10) holds for all x ∈ Xτ .
With the discretization of the state space, we only need
to sample a finite number of points to validate the barrier
certificates. However, the number of required sampling points
is still very large. The following adaptive sampling strategy
further reduces the number of sampling points required.
Proposition 3.2: If the following condition is satisfied at
x∈X,


∂h
∂h
∂h
g(x)u + m(x) − kδ
σ (x)
max
∂x
∂x
∂x
u∈U
∂h
(12)
+ f (x) + γ h(x) ≥ (Lḣ + γ Lh )kτ τ ,
∂x
with kτ ≥ 0, then the safety barrier constraint (11) is satisfied
for all y ∈ X such that kx − yk ≤ kτ τ .
Proof: The proof is similar to lemma 3.1.
Leveraging the Lipschitz continuity of the barrier certificates,
we can adaptively sample the state space without losing
safety guarantees. Sparse sampling is performed at places
with large safety margin, while dense sampling is only
required at places with small safety margin.

B. Parameterization of the Barrier Certificates
Because maximizing the volume of C is a non-convex
problem in general, we can parameterize the barrier certificate hµ (x) with µ to simplify the optimization problem. For
example, hµ (x) can be formulated as 1 − Z(x)T µ Z(x), where
Z(x) is the vector of monomials, and µ is a positive semidefinite matrix. Then maximizing vol(C ) is equivalent to
minimize the trace of µ . Further simplification can be made
to fix the shape of C (by optimizing only with the known
dynamics) and enlarge the level set of barrier certificates.
With the shape parameterization and adaptive sampling
technique, the barrier certificate maximization problem (9)
can be written as
max vol(C )
µ


∂ hµ
∂ hµ
∂ hµ
σ (x)
g(x)u +
m(x) − kδ
s.t. max
∂x
∂x
∂x
u∈U
∂ hµ
+
f (x) + γ hµ (x) ≥ (Lḣ + γ Lh )τ , ∀x ∈ C ∩ Xτ .
∂x
(13)
In order to increase the learning efficiency during the
exploration phase, the most uncertain state in C is sampled,
xnext = argmax σ (x).

(14)

x∈C ∩Xτ

It is assumed that a nominal exploration controller û can
always be designed to drive the system from the current state
x to xnext , i.e., û = GoTo(x, xnext ). Then the safety barrier
certificates are enforced through a QP-based controller to
“rectify” the nominal control such that the system is always
safe,
u∗ = argmin J(u) = ku − ûk2
u∈U

∂h
∂h
∂h
g(x)u + m(x) − kδ
σ (x)
(15)
∂x
∂x
∂x
∂h
f (x) + γ h(x) ≥ 0.
+
∂x
Therefore, the actual exploration controller u∗ tries to stay
as close as possible to the desired controller û, while always
honoring the safety requirements. The exploration phase ends
when the safe region C does not grow any more. The learned
maximum barrier certificates can be further used to regulate
other control tasks the system want to achieve.
s.t.

C. Overview of the Safe Learning Algorithm
An overview of the barrier certificates based safe learning
algorithm is provided in Algorithm 1. At the beginning, a
conservative barrier certified safe region C0 is provided. The
most uncertain state xnext is computed based on the current
GP model. Then, the QP based controller (15) is used to
ensure that the system is driven to xnext without ever leaving
Cn . After updating the GP model with the sampled data
at xnext , the barrier certificate optimization problem (13) is
solved. The adaptive sampling technique (12) is adopted here
to reduce the number of states to be sampled. This process
is repeated until the safe region Cn stops growing.

Algorithm 1 Barrier Certificates based Safe Learning
Input: Initial safe set C0 ⊆ X , GP model G P(0, k(x, x′ )),
discretization Xτ , tolerance ε
Output: Final safe set Cn
Initialization : n = 0, x = x0
1: repeat
2:
n = n+1
3:
Find xnext with (14)
4:
Design nominal controller û = GoTo(x, xnext )
5:
Rectify û with (15) and drive to xnext
6:
Sample xnext , update GP
7:
Expand vol(Cn ) with (13) and adaptive sampling (12)
8: until vol(Cn )-vol(Cn−1 )≤ ε
9: return Cn

IV. O NLINE L EARNING

OF

Q UADROTOR DYNAMICS

The safe learning approach developed in Section III relies
on a learning controller that drives the system to explore
interested states. The challenge of designing this learning
controller is that the 3D quadrotor system considered in this
paper is highly nonlinear and unstable. In this section, we
will present a recursive learning controller based on GP to
learn the complex quadrotor dynamics online.
A. Differential Flatness of 3D Quadrotor Dynamics
The quadrotor is a well-modelled dynamical system with
forces and torques generated by four propellers and gravity
[27]. The relevant coordinate frames and Euler angles (roll
φ , pitch θ , and yaw ψ ) are illustrated in Fig. 3. The world,
body, and intermediate frames (after yaw angle rotation) are
denoted by the subscripts w, b, and c, respectively.

Fig. 3: Quadrotor coordinate frames.
The Euler angles are defined with the ZY X convention.
Hence, the rotation matrix from the body frame to the world
frame can be written as


cθ cψ sφ sθ cψ − cφ sψ cφ sθ cψ + sφ sψ
R = cθ sψ sφ sθ sψ + cφ cψ cφ sθ sψ − sφ cψ  ,
−sθ
sφ cθ
cφ cθ
where sθ and cθ stand for sin θ and cos θ , respectively.

Here, we adopt the quadrotor model used in [10] to
describe the nonlinear quadrotor dynamics,

r̈
= gzw + m1 Rzw fz ,





 φ̇ 
1 sφ t θ
cφ t θ
(16)
 



cφ
−sφ  ω ,
 θ̇  = 0



ψ̇
0 sφ scθ cφ scθ

where zw = [0 0 1]T , and r = [x, y, z]T , m, and g are the
position of the center of mass, the mass, and the gravitational
acceleration of the quadrotor, respectively. t θ and scθ are
short for tan θ and sec θ . The control inputs of the quadrotor
are the body rotational rates (ω = [ωx , ωy , ωz ]T ) and the thrust
( f ). It is assumed that the body rotational rates of quadrotor
are directly controllable through the fast response onboard
controller, due to the small rotational inertia and high torque
features of quadrotors [10].
Similar to [27], the dynamics in (16) is differentially
flat with the flat output chosen as η = [rT , ψ T ]T . The full
state q = [rT , ṙT , θ , φ , ψ ]T and control u = [ f , ω T ]T can... be
represented as an algebraic function of [η T , η̇ T , η̈ T , η T ].
With the differential flatness property, quadrotor trajectory
planning can be simplified as smooth parametric curves.
Given a desired trajectory ηd (t) ∈ C3 that is three times
T ] can
differentiable, the feed forward control uFF = [ fFF , ωFF
be derived by inverting the dynamics in (16),

fF F = −mkr̈d − gzw k,


 



1
0
−sθd
φ̇d
 


ω
=
0 cφd sφd cθd   θ̇d 
FF



0 −sφd cφd cθd
ψ̇d
q
where θd = atan2(βa , βb ), φd = atan2(βc , βa2 + βb2 ), βa =
−ẍd cos ψd − ÿd sin ψd , βb = −z̈d + g, and βc = −ẍd sin ψd +
ÿd cos ψd .
Differential flatness only gives the feed forward control
uFF . In addition, the unknown model error and tracking error
need to be handled by a feedback control uFB . The actual
control applied to the quadrotor is u = uFF + uFB , where

fFB = K p < Rzw , rd − r > +Kd < Rzw , ṙd − ṙ >,










yd − y
φd − φ
φ̇d − φ̇







ωFB = K p  θd − θ  + Kd  θ̇d − θ̇  + K̄ p x − xd 



0
ψd − ψ
ψ̇d − ψ̇
Note that with an inaccurate model, a high-gain feedback
controller is needed to counteract both the model error and
disturbances. As a better model is learned over time, only
a low-gain feedback controller is needed with an improved
tracking performance [16].
B. Learning based Control Using Gaussian Process
The previous section deals with precise quadrotor models.
But it is often difficult to acquire accurate parameters for
quadrotor systems. In addition, the model (16) neglects the
uncertain effects of damping, drag force, and wind disturbances. Here, we will use GP models to learn the unmodeled

dynamics. The unmodeled dynamics can be captured with six
GPs along each dimension in the state space, i.e.,



′ ))
G
P
(0,
k(q,
q

1






r̈
= gzw + m1 Rzw fz + G P 2 (0, k(q, q′ )) ,




G P 3 (0, k(q, q′ ))
  



′ ))

φ̇
φ
t
θ
c
φ
t
θ
G
P
(0,
k(q,
q
1
s
4



  





cφ
−sφ  ω + G P 5 (0, k(q, q′ )) ,
 θ̇  = 0



G P 6 (0, k(q, q′ ))
ψ̇
0 sφ scθ cφ scθ

where the input to the GPs is q = [rT , ṙT , θ , φ , ψ ]T , and the
observations for the GPs are s = [r̈T , φ̇ , θ̇ , ψ̇ ]T , respectively.
At a new query point q∗ , the mean mi (q∗ ) and variance
σi2 (q∗ ) of the unknown dynamics can be inferred with (5).
Based on the learned dynamics, a differential flatness based
feed forward controller can be derived as,

fF F = −mkr̈d − [m1(q), m2 (q), m3 (q)]T − gzw k,







1
0
−sθd
φ̇d − m4(q)




ωFF = 0 cφd sφd cθd   θ̇d − m5(q)  ,


0 −sφd cφd cθd
ψ̇d − m6 (q)
q
where θd = atan2(β̄a , β̄b ), φd = atan2(β̄c , β̄a2 + β̄b2 ),
β̄a = −(ẍd − m1 (q)) cos ψd − (ÿd − m2 (q)) sin ψd ,
β̄b = −(z̈d − m3 (q)) + g, and β̄c = −(ẍd − m1 (q)) sin ψd +
(ÿd − m2 (q)) cos ψd .
C. Recursive Online GP Learning
One issue with the GP regression is that the time complexity of GP inference is O(N 3 ), where N is the number of
data points. The majority of the time is used to compute the
inverse of the kernel matrix K. While various approximation
methods can be used to reduce the GP inference time, it is
still challenging to perform online GP inference for complex
dynamically systems like quadrotor. Here, we propose a
recursive online GP Learning method to compute the exact
GP inference.
As the quadrotor moves forward, we will actively add
multiple relevant data points into the kernel matrix at each
time step. At the same time, the data points that contribute the
least to the inference are deleted. The recursive data addition
and deletion operations are described as following.
1) Adding Multiple New Data to the Kernel Matrix: Let
the kernel matrix at the ith time step be Ki , we can save the
matrix inverse result from the previous step as Li = (Ki +
σn2 I)−1 . Denote the number of new data to be added as M.
With the new data yi+1 and kernal vector ki+1 , we have
 −1
−1
Li
ki+1
Li+1 =
kT
ci+1 + σn2 I
 i+1
T Lk
−1 T
Li + Li ki+1 (ci+1 + σn2 I − ki+1
i i+1 ) ki+1 Li
=
2
T
−1
T
−(ci+1 + σn I − ki+1Li ki+1 ) ki+1 Li

T Lk
−1
Li ki+1 (ci+1 + σn2 I − ki+1
i i+1 )
.
T Lk
(ci+1 + σn2 I − ki+1
i i+1 )
Notice that inversion operation only needs to be performed
on a M × M matrix rather than a large N × N matrix.

2) Deleting Multiple Old Data from the Kernel Matrix:
After deleting M data points from the old Kernel matrix
inversion Li = (Ki + σn2 I)−1 , the new inverse of the kernel
matrix becomes L̄i = (K̄i + σn2 I)−1 .
First, the data to be deleted is permuted to the bottom
of the kernel matrix with a permutation matrix Pπ , where
π : N → N is a permutation of N elements. The permuted
kernel matrix is Kip = Pπ Ki PπT , which can be written into a
block matrix form,


K̄i Ei
P
Ki = T
,
Ei Fi
where Ei , Fi are the known parts to be deleted. Similarly,
LPi

=
=

Pπ Li PπT
−1
 −1
L̄i
Ei
.
EiT Fi + σn2 I

Since LPi is known, it can be written into block matrix form
with the same block dimensions with (IV-C.2),


Ai Bi
P
Li = T
.
Bi Ci
With the block matrix inversion rule, L̄i can be recovered as
L̄i = Ai − BiCi−1 BTi ,
which means to perform the deletion operation, the only
matrix inverse required is Ci−1 ∈ RM×M .
With the recursive data addition and deletion method, the
GP inference can be obtained efficiently online.
V. S IMULATION R ESULTS
The GP based learning algorithm is validated on a simulated quadrotor model in two examples, i.e., online learning
of quadrotor dynamics and learning safety barrier certificates.
In the simulation, the actual weight of the quadrotor is 1.4
times the weight used in the computation. In addition, an
unknown constant wind of 0.1g is applied in the environment as illustrated in Fig. 4. Since the standard fixed pitch
quadrotor cannot generate reverse thrust, the thrust control
is limited to fz ∈ [−1.8mg, 0]. This simulation setup is very
challenging, because the learning based quadrotor controller
needs to deal with very inaccurate model and limited thrust.
A. Online Learning of Quadrotor Dynamics
In the first example, the quadrotor is commanded to track
a nominal trajectory (illustrated in Fig. 4) using a differential
flatness based controller with the given inaccurate model. A
PD controller is wrapped around to stabilize the quadrotor.
During the simulation, the quadrotor is intentionally pushed
to unknown regions that has not been explored before. This
will help us evaluate the scalability of the algorithm.
The desired trajectory of the quadrotor is given as
η̂ = [r̂(t)T , ψ̂ (t)] ∈ C3 , while the actual trajectory is η =
[r(t)T , ψ (t)]. In practice, the actual trajectory might deviate
significantly from the desired trajectory when the model is

Time per iteration (ms)

-3
-2.5
-2

Z

-1.5
-1
-0.5
0

20
15
10
5
0
0

10

3

0.5

20

30

40

Time (s)

2
-3

Fig. 6: Recursive GP inference time per iteration.

1

-2

0

-1
0

-1
1

-2

2

Y

-3

3

X

Fig. 4: A simulated quadrotor flies in an unknown wind field
with an inaccurate model.
very inaccurate. To track the desired trajectory, the nominal
trajectory is designed with a pole placement controller,
... ...
r i = r̂ i − K · [(ri − r̂i ), (ṙi − r̂˙i ), (r̈i − r̂¨i )]T .

0.4

0.4

0.3

0.3

Tracking error (m)

Tracking error (m)

In the simulation, the sample size of the recursive GP
model is fixed at 300 data points. At each time step, the
most irrelevant data point is thrown away, and the most
relevant data point is added to the GP model. The data
relevance is decided by the kernel function k(q, q∗ ), where
q = [rT , ṙT , θ , φ , ψ ]T . It can observed that the tracking error
of the learning based controller is significantly smaller than
the tracking error without GP inference, as shown in Fig. 5.

0.2
eX
eY
eZ

0.1
0

0.2

eX
eY
eZ

0.1
0
-0.1

-0.1
0

10

20

30

40

Time (s)

(a) Tracking error without GP

0

10

20

30

40

Time (s)

(b) Tracking error with GP

Fig. 5: Tracking error of the differential flatness based flight
controller with and without GP inference.
With the recursive learning strategy, it is demonstrated in
Fig. 6 that the GP inference time is always kept below 20ms.
Thus, the recursive GP inference method is very suitable for
online learning of quadrotor dynamics.
By pushing the quadrotor to unexplored regions, we can
found that learning with q′ = [ṙT , θ , φ , ψ ]T yields much
better scalability than learning with q = [rT , ṙT , θ , φ , ψ ]T .
The reason might be the position r is not as important as
other features in the current simulation setup.
B. Learning Safety Barrier Certificates
In this example, the motion of the quadrotor is constrained
within an ellipsoid safe region, i.e.,
y2
(z + 0.8)2
x2
+
+
≤ 1.
0.16 0.16
0.36

The quadrotor is controlled to fly back and forth on a
vertical path inside the ellipsoid. The goal is to learn how
aggressively the quadrotor can fly in the z direction with an
inaccurate model and limited thrust.
The barrier certificates are parameterized as
(z + 0.8)2
− µ ż2
0.36
x2
y2
ẋ2
ẏ2
−
−
−
−
≥ 0,
0.16 0.16 0.25 0.25
where µ is the barrier parameter to regulate how fast the
quadrotor can fly in the z direction. Small values of µ
correspond to large admissible speed ż, which means more
aggressive flight behavior. Thus, the objective of the learning
process is to minimize µ with the collected data.
To reduce the number of required sample points, the
adaptive sampling strategy developed in Section III-A was
adopted. An illustrative example of the adaptive sampling
strategy is given in Fig. 7. It can be observed that places
closer to the boundary of the safe region (z = −1.2 and
z = −0.2) are sampled much denser than the place closer to
the center of the safe region (z = 0). Furthermore, downward
speed (ż > 0) is sampled much denser than the upward speed
(ż < 0). This might be caused by the lack of reverse thrust
to counter the unmodeled dynamics.
A conservative barrier certificate (µ = 6.3) is provided at
the beginning of the learning process. Then, the quadrotor
gradually explores the safe region C0 and expands it to Cn
(µ = 0.6), as illustrated in Fig. 8. The nominal exploration
controller is always regulated by the barrier certificates using
the QP-based controller in (15). During the learning process,
the quadrotor never leaves the barrier certified safe region.
h µ (r) =

1−

VI. C ONCLUSIONS
A safe learning algorithm based on barrier certificates was
developed in this paper. The learning controller is regulated
by the barrier certificates, such that the system never enters
the unsafe region. The unmodel dynamics of the system
was approximated with a Gaussian Process, from which a
high probability safety guarantee for the dynamical system
was derived. The barrier certified safe region is gradually
expanded as the uncertainty of the system dynamics is
reduced with more data. This safe learning technique was
applied on a quadrotor system with 3D nonlinear dynamics.
The computation time of this learning method is reduced significantly with an adaptive sampling strategy and a recursive

C1 (µ = 1.23)
C2 (µ = 0.92)

1

1

0.5

0.5

0

0

-0.5

-0.5

-1

-1

-1.5

-1.5

-1.5

-1

-0.5

0

z
Fig. 7: Adaptive sampling of the state space. The region
enclosed by the solid green ellipse C1 is the current safe
region, while the region enclosed by the dashed red ellipse
C2 is the optimized next safe region. The green cross markers
and red asterisk markers are the data points already sampled
and to be sampled, respectively. The red circles centered
at those sample points are the confident safe regions based
on (12). All the unexplored region between C1 and C2 are
covered by the circular confident safe region.

GP inference method. Simulation results demonstrated the
effectiveness of the proposed method.
R EFERENCES
[1] A. K. Akametalu, J. F. Fisac, J. H. Gillula, S. Kaynama, M. N.
Zeilinger, and C. J. Tomlin. Reachability-based safe learning with
gaussian processes. In Decision and Control (CDC), 2014 IEEE 53rd
Annual Conference on, pages 1424–1431. IEEE, 2014.
[2] A. D. Ames, X. Xu, J. W. Grizzle, and P. Tabuada. Control barrier
function based quadratic programs for safety critical systems. IEEE
Transactions on Automatic Control, 62(8):3861–3876, 2017.
[3] A. Aswani, H. Gonzalez, S. S. Sastry, and C. Tomlin. Provably
safe and robust learning-based model predictive control. Automatica,
49(5):1216–1226, 2013.
[4] A. Bemporad and M. Morari. Robust model predictive control: A
survey. Robustness in identification and control, pages 207–226, 1999.
[5] F. Berkenkamp, R. Moriconi, A. P. Schoellig, and A. Krause. Safe
learning of regions of attraction for uncertain, nonlinear systems with
gaussian processes. In Decision and Control (CDC), 2016 IEEE 55th
Conference on, pages 4661–4666. IEEE, 2016.
[6] F. Berkenkamp, M. Turchetta, A. P. Schoellig, and A. Krause. Safe
model-based reinforcement learning with stability guarantees. 2017.
[7] G. Chesi. Domain of attraction: analysis and control via SOS
programming, volume 415. Springer Science & Business Media, 2011.
[8] M. P. Deisenroth, D. Fox, and C. E. Rasmussen. Gaussian processes
for data-efficient learning in robotics and control. IEEE Transactions
on Pattern Analysis and Machine Intelligence, 37(2):408–423, 2015.
[9] H. Fukushima, T.-H. Kim, and T. Sugie. Adaptive model predictive
control for a class of constrained linear systems based on the comparison model. Automatica, 43(2):301–308, 2007.
[10] M. Hehn and R. DAndrea. Real-time trajectory generation for
quadrocopters. IEEE Transactions on Robotics, 31(4):877–892, 2015.
[11] N. Hovakimyan, C. Cao, E. Kharisov, E. Xargay, and I. M. Gregory.
L1 adaptive controlfor safety-critical systems. IEEE Control Systems,
31(5):54–104, 2011.
[12] Y. Ito, K. Fujimoto, and Y. Tadokoro. Second-order bounds of gaussian
kernel-based functions and its application to nonlinear optimal control
with stability. arXiv preprint arXiv:1707.06240, 2017.

C0 (µ = 6.3)
Cn (µ = 0.6)

1.5

ż

ż

1.5

-1.5

-1

-0.5

0

z
Fig. 8: Initial and final barrier certificates. The regions enclosed by the solid green ellipse (C0 ) and dashed red ellipse
(Cn ) are the initial and final barrier certified safe regions,
respectively. The green cross markers and red asterisk makers
are the sampled data points.

[13] S. M. Khansari-Zadeh and A. Billard. Learning control lyapunov
function to ensure stability of dynamical system-based robot reaching
motions. Robotics and Autonomous Systems, 62(6):752–765, 2014.
[14] J. Lofberg. Yalmip: A toolbox for modeling and optimization in
matlab. In Computer Aided Control Systems Design, 2004 IEEE
International Symposium on, pages 284–289. IEEE, 2005.
[15] Q. Nguyen, A. Hereid, J. W. Grizzle, A. D. Ames, and K. Sreenath. 3d
dynamic walking on stepping stones with control barrier functions. In
Decision and Control (CDC), 2016 IEEE 55th Conference on, pages
827–834. IEEE, 2016.
[16] D. Nguyen-Tuong, J. R. Peters, and M. Seeger. Local gaussian process
regression for real time online model learning. In Advances in Neural
Information Processing Systems, pages 1193–1200, 2009.
[17] Y. Pan and E. Theodorou. Probabilistic differential dynamic programming. In Advances in Neural Information Processing Systems, pages
1907–1915, 2014.
[18] S. Prajna, A. Jadbabaie, and G. J. Pappas. Stochastic safety verification
using barrier certificates. In Decision and Control, 2004. CDC. 43rd
IEEE Conference on, volume 1, pages 929–934. IEEE, 2004.
[19] S. Prajna, A. Jadbabaie, and G. J. Pappas. A Framework for Worstcase and Stochastic Safety Verification Using Barrier Certificates.
Automatic Control, IEEE Transactions on, 52(8):1415–1428, 2007.
[20] C. E. Rasmussen. Gaussian processes for machine learning. 2006.
[21] H. Ravanbakhsh and S. Sankaranarayanan. Learning lyapunov (potential) functions from counterexamples and demonstrations. In Robotics:
Science and Systems (RSS), 2017.
[22] J. Vinogradska, B. Bischoff, D. Nguyen-Tuong, H. Schmidt, A. Romer,
and J. Peters. Stability of controllers for gaussian process forward
models. In Proceedings of the 33rd International Conference on
International Conference on Machine Learning, pages 545–554, 2016.
[23] L. Wang, A. D. Ames, and M. Egerstedt. Multi-objective compositions
for collision-free connectivity maintenance in teams of mobile robots.
In Decisions and Control Conference (CDC), pages 2659–2664, 2016.
[24] L. Wang, A. D. Ames, and M. Egerstedt. Safe certificate-based
maneuvers for teams of quadrotors using differential flatness. In IEEE
International Conference on Robotics and Automation (ICRA), pages
3293–3298, 2017.
[25] G. Wu and K. Sreenath. Safety-critical control of a planar quadrotor. In
2016 American Control Conference (ACC), pages 2252–2258. IEEE,
2016.
[26] X. Xu, P. Tabuada, J. W. Grizzle, and A. D. Ames. Robustness of
control barrier functions for safety critical contro. In Analysis and
Design of Hybrid Systems, 2015 IFAC Conference on. IEEE, Oct 2015.
[27] D. Zhou and M. Schwager. Vector field following for quadrotors
using differential flatness. In 2014 IEEE International Conference on
Robotics and Automation (ICRA), pages 6567–6572. IEEE, 2014.

