Verifying Controllers Against Adversarial Examples with Bayesian
Optimization

arXiv:1802.08678v2 [cs.SY] 26 Feb 2018

Shromona Ghosh1 , Felix Berkenkamp2 , Gireeja Ranade3 , Shaz Qadeer3 , Ashish Kapoor3
Abstract— Recent successes in reinforcement learning have
lead to the development of complex controllers for realworld robots. As these robots are deployed in safety-critical
applications and interact with humans, it becomes critical to
ensure safety in order to avoid causing harm. A first step in
this direction is to test the controllers in simulation. To be
able to do this, we need to capture what we mean by safety
and then efficiently search the space of all behaviors to see
if they are safe. In this paper, we present an active-testing
framework based on Bayesian Optimization. We specify safety
constraints using logic and exploit structure in the problem in
order to test the system for adversarial counter examples that
violate the safety specifications. These specifications are defined
as complex boolean combinations of smooth functions on the
trajectories and, unlike reward functions in reinforcement
learning, are expressive and impose hard constraints on the
system. In our framework, we exploit regularity assumptions on
individual functions in form of a Gaussian Process (GP) prior.
We combine these into a coherent optimization framework using
problem structure. The resulting algorithm is able to provably
verify complex safety specifications or alternatively find counter
examples. Experimental results show that the proposed method
is able to find adversarial examples quickly.

I. I NTRODUCTION
In recent years, research in control theory and robotics
has focused on developing efficient controllers for robots
that operate in the real world. Controller synthesis techniques
such as reinforcement learning, optimal control, and model
predictive control have been used to synthesize complex
policies. However, if there is a large amount of uncertainty
about the real world environment that the system interacts
with, the robustness of the synthesized controller becomes
critical. This is particularly true in safety-critical systems,
where the actions of an autonomous agent may affect human
lives. This motivates us to provably verify the properties of
controllers in simulation before deployment in the real world.
In this paper, we present an active machine learning
framework that is able to verify black-box systems against,
or alternatively find, adversarial counter examples to a given
set of safety specifications. We test the controller safety
under uncertainty that arises from stochastic environments
and errors in modeling. In essence, we actively search for
The primary part of the research for this paper was done at Microsoft
Research, Redmond while Shromona Ghosh and Felix Berkenkamp were
interns.
1 Electrical Engineering and Computer Science Department, University
of California, Berkeley, shromona.ghosh@berkeley.edu
2 Department of Computer Science, ETH Zürich, Switzerland

befelix@inf.ethz.ch
3

Microsoft

Research,

Redmond

akapoor}@microsoft.com

{giranade,

qadeer,

Fig. 1. Framework for active testing. We query the simulation of the
system with environment parameters w to obtain trajectories ξ. We test these
for safety violations in a specification monitor. The Bayesian optimization
framework actively queries new parameters w that are promising candidates
to find counter examples that violates the safety specification.

adversarial environments under which the controller could
have to operate that lead failure modes in simulation.
Historically, designing robust controllers has been considered in control theory [1], [2]. A common issue with
these techniques is that, although they consider uncertainty,
they rely on simple linear models of the underlying system.
This means that resulting controllers are often either overly
conservative or violate safety constraints if they fail to
capture nonlinear effects.
For nonlinear models with complex dynamics, reinforcement learning has been successful for synthesizing high
fidelity controllers. Recently, algorithms based on reinforcement learning that can handle uncertainty have been
proposed [3]–[5], where the performance is measured in
expectation. A fundamental issue with learned controllers is
that it is difficult to provide formal guarantees for safety
in the presence of uncertainty. For example, a controller
for an autonomous vehicle must consider human driver
behaviors, pedestrian behaviors, traffic lights, uncertainty
due to sensors, etc. Without formally verifying that these
controllers are indeed safe, deploying them on the road could
lead to loss of property or human lives.
Formal safety certificates, i.e., mathematical proofs for
safety, have been considered in the formal methods community, where safety requirements are referred to as a
specification. There, the goal is to verify that the behaviors
of a particular model satisfies a specification ( [6], [7]).
Synthesizing controllers which satisfy a high level temporal
specification have been studied in the context of motion planning [8] and for cyber-physical systems [9]. However, these
techniques rely on simple model dynamics. For nonlinear
systems, reachability algorithms based on level set methods
have been used to approximate backward reachable sets for
safety verification [10], [11]. However, these methods suffer
from two major drawbacks: (1) the curse of dimensionality

Accepted final version. To appear in 2018 IEEE International Conference on Robotics and Automation.
c 2018 IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any current or future media, including
reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution to servers or lists, or reuse of
any copyrighted component of this work in other works.

of the state space, which limits them to low-dimensional
systems; and (2) a priori knowledge of the system dynamics.
A dual, and often simpler, problem is falsification, which
tests the system within a set of environment conditions for
adversarial examples. Adversarial examples have recently
been considered for neural networks [12]–[15], where the
input is typically perturbed locally in order to find counterexamples. In [16], the authors compute adversarial perturbations for a trained neural network policy for a subset
of white box and black-box systems. However, these local
perturbations are often not meaningful for dynamic systems.
Recently, [17], [18] have focused on testing of closedloop safety critical systems with neural networks by finding
“meaningful” perturbations.
Testing black-box systems in simulators is a well studied
problem in the formal methods community [19]–[21]. The
heart of research in black-box testing focuses on developing
smarter search techniques which efficiently samples the
uncertainty space. Indeed, in recent years, several sequential
search algorithms based on heuristics such as Simulated Annealing [21], Tabu search [22], and CMA-ES [23] have been
suggested. Although these algorithms sample the uncertainty
space efficiently, they do not utilize any of the information
gathered during previous simulations.
One active method that has been used recently for testing
black-box systems is Bayesian Optimization (BO) [24], an
optimization method that aims to find the global optimum of
an a priori unknown function based on noisy evaluations.
Typically, BO algorithms are based on Gaussian Process
(GP [25]) models of the underlying function and certain algorithms provably converge close to the global optimum [26].
It has been used in robotics to, for example, safely optimize
controller parameters of a quadrotor [27]. In the testing
setting, BO has been used to actively find counter examples
by treating the search problem as a minimization problem
in [28] over adversarial control signals. However, the authors
do not consider the structure of the problem and thereby
violate the smoothness assumptions made by the GP model.
As a result, their methods are slow to converge or may fail
to find counterexamples.
In this paper, we provide a formal framework that uses
BO to actively test and verify closed-loop black-box systems
in simulation. We model the relation between environments
and safety specification using GPs and use BO to predict
the environment scenarios most likely to cause failures
in our controllers. Unlike previous approaches, we exploit
structure in the problem in order to provide a formal way
to reason across multiple safety constraints in order to
find counterexample. Hence, our approach is able to find
counterexamples more quickly than previous approaches.
Our main contributions are:
• An active learning framework for testing and verifying
robotic controllers in simulation. Our framework can
find adversarial examples for a synthesized controller
independent of its structure or how it was synthesized.
• A common GP framework to model logical safety
specifications along with theoretical analysis on when

a system is verified.
II. P ROBLEM S TATEMENT
We address the problem of testing complex black-box
closed-loop robotic systems in simulation. We assume that
we have access to a simulation of the robot that includes the
control strategy, i.e., the closed-loop system. The simulator is
parameterized by a set of parameters w ∈ W, which model
all sources of uncertainty. For example, they can represent
environment effects such as weather, non-deterministic components such as other agents interacting with the simulator,
or uncertain parameters of the physical system, e.g., friction.
The goal is to test whether the system remains safe for
all possible sources of uncertainty in W. We specify these
safety constraints on finite-length trajectories of the system
that can be obtained by simulating the robot for a given set of
environment parameters w ∈ W. Safety constraints on these
trajectories are specified using logic. We explain this in detail
in Sec. III-A, but the result is a specification ϕ that can, in
general, be written as a requirement ϕ(w) > 0, ∀w ∈ W.
For example, ϕ can encode state or input constraints that
have to be satisfied over time.
We want to test whether there exists an adversarial
example w ∈ W for which the specification is violated,
i.e., ϕ(w) < 0. Typically, adversarial examples are found
by randomly sampling the environment and simulating the
behaviors. However, this approach does not provide any guarantees and does not allow us to conclude that no adversarial
example exist if none are found in our samples. Moreover,
since high-fidelity simulations can often be very expensive,
we want to minimize the number of simulations that we have
to carry out in order to find a counterexample.
We propose an active learning framework for testing,
where we utilize the results from previous simulation runs
to make more informed decisions about which environment
to simulate next. In particular, we pose the search problem
for a counterexample as an optimization problem,
argmin ϕ(w),

(1)

w∈W

where we want to minimize the number of queries w
until a counterexample is found or we can verify that
no counterexample exists. The main challenge is that the
functional dependence ϕ(·) between parameters in W and
the specification is unknown a priori, since we treat the
simulator as a black-box. Solving this problem is difficult
in general, but we can exploit regularity properties of ϕ(w).
In particular, in the following we use GP to model the
specification and use the model to pick parameters that are
likely to be counterexamples.
III. BACKGROUND
In this section, we introduce an overview of formal
safety specifications and Gaussian processes, which we use
in Sec. IV to verify the closed-loop black-box system.

A. Safety Specification
In the formal methods community, complex safety requirements are expressed using automatons [29] and temporal
logic [30], [31]. These allow us to specify complex constraints, which can also have temporal dependence.
Example 1. A safety constraint for a quadcopter might be
that the quadcopter cannot fly at an altitude h greater than
3 m when the battery level b is below 30%.”
In logic, we can express this as “b < 0.3 implies(→)
h < 3”, which in words says if the battery level is less than
30% the quadcopter is flying at a height less than 3 m.
Importantly, these kind specifications make no assumptions about the underlying system themselves. They just
state requirements that must hold for all simulations in W.
Formally, a logic specification is a function that tests properties of a particular trajectory. However, we will continue to
write ϕ(w) to denote the specification that tests trajectories
generated by the simulator with parameters w.
A specification ϕ consists of multiple individual constraints, called predicates, which form the basic building
blocks of the logic. These predicates can be combined using
a syntax or grammar of logical operations:
ϕ := µ | ¬µ | ϕ ∧ ψ | ϕ ∨ ψ.

(2)

where µ : Ξ → R is a predicate, and is assumed to be
a smooth and continuous function of a trajectory ξ ∈ Ξ.
The constraint µ > 0 forms the basic building block of
the overall system specification ϕ. We say a predicate is
satisfied if µ(ξ) is greater than 0 or falsified otherwise. The
operations ¬, ∧, ∨ represent negation, conjunction(and) and
disjunction(or) , respectively. These basic operations can
be combined to define complex boolean formula such as
implication, →, and if-and-only-if, ↔ using the rules
ϕ → ψ := ¬ϕ∨ψ, and ϕ ↔ ψ := (¬ϕ∧¬ψ)∨(ϕ∧ψ). (3)
Since µ is a real valued function, we can convert these
boolean logic statements into an equivalent equation with
continuous output, which defines the quantitative semantics,
µ(ξ) := µ(ξ),
¬µ(ξ) := −µ(ξ),

(ϕ ∧ ψ)(ξ) := min(ϕ(ξ), ψ(ξ)),

(ϕ ∨ ψ)(ξ) := max(ϕ(ξ), ψ(ξ)).

Example 2. Let us look at the specification in Example 1,
ϕ := (b < 0.3) → (h < 3). Applying the re-write rule (3),
this can be written as ¬(b < 0.3) ∨ (h < 3). Applying the
quantitative semantics (4), we get ϕ = max(b > 0.3, h < 3),
which consists of two predicates, µ1 = b − 0.3 and
µ2 = 3 − h. Intuitively, this means ϕ > 0, i.e., the specification is satisfied, if the battery is greater than 30% or if
the quadcopter flies at an altitude less than 3m .
B. Gaussian Process
For general black-box systems, the dependence of the
specification ϕ(·) on the parameters w ∈ W is unknown
a priori. We use a GP to approximate each predicate µ(·) in
the domain W. We detail the modeling of ϕ(·) in Sec. IV.
The following introduction about GPs is based on [25].
GPs are non-parametric regression method from machine
learning, where the goal is to find an approximation of the
nonlinear function µ : W → R from an environment w ∈ W
to the function value µ. This is done by considering the
function values µ(w) to be random variables, such that any
finite number of them have a joint Gaussian distribution.
The Bayesian, non-parametric regression is based on
a prior mean function and the kernel function k(w, w0 ),
which defines the covariance between the function values
µ(w), µ(w0 ) at two points w, w0 ∈ W. We set the prior
mean to zero, since we do not have any knowledge about the
system. The choice of kernel function is problem-dependent
and encodes assumptions about the unknown function.
We can obtain the posterior distribution of a function
value µ(w) at an arbitrary state w ∈ W by conditioning
the GP distribution of µ on a set of n past measurements, yn = (µ̂(w1 ), . . . , µ̂(wn )) at environment scenarios Wn = {w1 , . . . , wn }, where µ̂(w) = µ(w) + ω and
ω ∼ N (0, σ 2 ) is Gaussian noise. The posterior over µ(w)
is a GP distribution again, with mean mn (w), covariance
kn (w, w), and variance σn (w):
mn (w) = kn (w)(Kn + In σ 2 )−1 yn ,
kn (w, w0 ) = k(w, w0 ) − kn (w)(Kn + In σ 2 )−1 kTn (w0 ),
σn2 (w) = kn (w, w0 ),

(4)

This allows us to confirm that a logic statement ϕ holds true
for all trajectories generated by simulators W, by confirming
that the function ϕ(w) takes positive values for all w ∈ W.
In the quantitative semantics (4), the satisfaction of a
requirement is no longer a yes or no answer, but can be
quantified by a real number. The nature of this quantification
is similar to that of a reward function, where lower values
indicate a larger safety violation. This allows us to introduce
a ranking among failures: ϕ(w1 ) < ϕ(w2 ) implies w1 is
a more ”dangerous” failure case than w2 . To guarantees
safety, we have to take a pessimistic outlook, and denote
ϕ(w) ≤ 0 as a violation and ϕ(w) > 0 as satisfaction of the
specification ϕ.

(5)

where the vector kn (w) = [k(w, w1 ), . . . , k(w, wn )] contains the covariances between the new environment, w,
and the environment scenarios in Wn , the kernel matrix Kn ∈ Rn×n has entries [Kn ](i, j) = k(wi , wj ), with
i, j ∈ {1, . . . , n}, and In ∈ Rn×n is the identity matrix.
C. Bayesian Optimization (BO)
In the following we use BO in order to find the minimum
of the unknown function ϕ, which we construct using the
GP models on µ in Sec. IV. BO uses a GP model to query
parameters that are informative about the minimum of the
function. In particular, the GP-LCB algorithm from [26]
uses the GP prediction and associated uncertainty in (5)
to trade off exploration and exploitation by, at iteration n,

selecting an environment according to
wn = argmin mn−1 (w) − βn1/2 σn−1 (w),

max

(6)

w∈W

where βn determines the confidence interval. We provide an
appropriate choice for βn in Theorem 1.
At each iteration, (6) selects parameters for which the
lower confidence bound of the GP is minimal. Repeatedly
evaluating the true function ϕ at samples given by (6) improves the GP model and decreases uncertainty at candidate
locations for the minimum, such that the global minimum is
found eventually [26].

min

max

𝜇*

𝜇+

𝜇()

𝜇('

Fig. 2. Equivalent parse tree T for ϕ in (7) to the function (8). We replace
the predicates µi with their corresponding pessimistic GP predictions to
obtain a lower bound on ϕ(w).

IV. ACTIVE T ESTING FOR C OUNTEREXAMPLES
In this section, we show how to model specifications ϕ
in (1) using GPs without violating smoothness assumptions
and use this to find adversarial counterexamples.
In order to use BO to optimize (1), we need to construct
reliable confidence intervals on ϕ. However, if we were to
model ϕ as a GP with commonly-used kernels, it would
need it to be a smooth function of w. Even though the predicates, µ, are typically smooth functions of the trajectories,
and hence smooth in w, conjunction and disjunction (min
and max) in (4) are non-smooth operators that render ϕ to
become non-smooth as well. Instead, we exploit the structure
of the specification ϕ and decompose ϕ into a parse tree,
where the leaf nodes are the predicates.

predicates, using GPs and then consider the non-smoothness
through the parse tree.
For each predicate µi in the parse tree T of ϕ, we construct
1/2 i
a lower confidence bound li = min−1 (w) − βn σn−1
(w),
where mi , σ i are the mean and standard deviation of the GP
corresponding to µi . From this, we can construct a lowerconfidence interval on ϕ as T (l1 (w), . . . , lq (w)), where we
replace the ith leaf node µi of the parse tree with the
pessimistic prediction li of the corresponding GP. Similar
to (6), the corresponding acquisition function for BO uses
this lower bound to select the next evaluation point,

Definition 1 (Parse Tree T ). Given a specification formula
ϕ, the corresponding parse tree, T , has leaf nodes that
correspond to function predicates, while other nodes are max
(disjunctions) and min (conjunctions).

Intuitively, the next environment selected to simulate is
the one that minimizes the worst-case predictions on ϕ.
Effectively, we propagate the confidence intervals associated
with the GP for each predicates through the parse tree T in
order to obtain predictions about ϕ directly. Note, that (9)
does not return an environment sample that minimizes the
satisfaction of all the predicates, it only minimizes the lower
bound on ϕ.
Algorithm 1 describes our active testing procedure. The
algorithm proceeds by first computing the parse tree T from
the specification, ϕ. At each iteration n of BO, we select
new environment parameters wn according to (9). We then
simulate the system with parameters wn and evaluate each
predicate µi on the simulated trajectories. Lastly, we update
each GP with the corresponding measurement of µi . The
algorithm either returns a counterexample that minimizes (1);
or when T (l1 (w), . . . , lq (w)) is greater then zero, and we
can conclude that the system has been verified.

A parse tree is an equivalent graphical representation of ϕ.
For example, consider the specification
ϕ := (µ1 ∨µ2 ) → (µ3 ∨µ4 ) = (¬µ1 ∧¬µ2 )∨(µ3 ∨µ4 ), (7)
where the second equality follows from De-Morgan’s law.
We can obtain an equivalent function ϕ(w) with (4),
ϕ(w) = max min(−µ1 (w), −µ2 (w)),

max(µ3 (w), µ4 (w)) .

(8)

The parse tree, T , for ϕ in (8) is shown in Fig. 2. We can
use the parse tree to decompose any complex specification
into min and max functions of the individual predicates; that
is, ϕ(w) = T (µ1 (w), . . . , µq (w)).
We now model each predicate µi (w) in the parse tree T of
ϕ with a GP and combine them with the parse tree to obtain
confidence intervals on the overall specification ϕ(w) for
BO. GP-LCB as expressed in (6) can be used to search for
the minimum for a single GP. A key insight to extending (6)
across multiple GPs, is that the minimum of (1) is, with high
probability, lower bounded by the lower-confidence interval
of one of the GPs used to model the predicates of ϕ. This is
because, the max and min operators do not change the value
of the predicates, but only make a choice between them. As
a consequence, we can model the smooth parts of ϕ, i.e., the

wn = argmin T (l1 (w), . . . , lq (w)).

(9)

w∈W

A. Theoretical Results
We can transfer theoretical convergence results for GPLCB [26] to the setting of Algorithm 1. To do this, we
need to make structural assumptions about the predicates.
In particular, we assume that they have bounded norm in
the Reproducing Kernel Hilbert Space (RKHS, [32]) that
corresponds to the GP’s kernel. P
These are well-behaved
functions of the form µi (w) =
j=0 αj ki (w, wj ) with
representer points wj and weights α that decay sufficiently
quickly. We leverage theoretical results from [33] and [27]

Algorithm 1 Active Testing with Bayesian Optimization
1:
2:
3:
4:
5:
6:
7:

procedure ACTIVE T ESTING(ϕ, W, β, GPs)
Build parse tree T based on specification ϕ
for n = 0, . . . do
. Until budget or convergence
1/2
li (w) = µi (w) − βn σi (w), i = 1, . . . , q
wn = argminw∈W T (l1 (w), . . . , lq (w))
Update each GP model of the predicates with
measurements (wn , µi (wn )).
return mini ϕ(wi ), the worst result.

that allow us to build reliable confidence intervals using the
GP models from Sec. III-B. We have the following result.
Theorem 1. Assume that each predicate µi has RKHS
norm bounded by Bi and that the measurement noise is
σ-sub-Gaussian. Select δ ∈q(0, 1), wn according to (9), and
P
P
1/2
i
; µi ).
let
βn = i Bi + 4σ 1 + ln(1/δ) + i I(yn−1
If T (l1 (wn ), . . . , lq (wn )) > 0, then with probability at least
1 − δ we have that minw∈W ϕ(w) > 0 and the system has
been verified against all environments in W.
i
i
Here I(yn−1
; µi ) is the mutual information between yn−1
,
the n − 1 noisy measurements of µi , and the GP prior
of µi . This function was shown to be sublinear in n for
many commonly-used kernels in [26], see the appendix
for more details. Theorem 1 states that we can verify the
system against adversarial examples with high probability,
by checking whether the worst-case lower-confidence bound
is greater than zero. We provide additional theoretical results
about the existence of a finite n such that the system can be
verified up to  accuracy in the appendix.

V. E VALUATION
In this section, we evaluate our method on several
challenging test cases. A Python implementation of our
framework and the following experiments can be found at
https://github.com/shromonag/adversarial testing.git
In order to use Algorithm 1, we have to solve the optimization problem (9). In practice, different optimization
techniques have been proposed to find the global minimum
of the function. One popular algorithm is DIRECT [34], a
gradient-free optimization method. An alternative is to use
gradient-based methods together with random-restarts. Particularly, we sample a large number of potential environment
scenarios at random from W, and run seperate optimization
routines to minimize (9) from these.
Another challenge is that the dimensionality of the optimization problem can often be very large. However, methods
that allow for more efficient computation do exist. These
methods reduce the effective size of the input space and
thereby make the optimization problem more tractable. One
possibility is to use random embedding to reduce the input
dimension as done in Random Embedding Bayesian Optimization (REMBO [35]). We can then model the GP in
this smaller input dimension and carry out BO in the lower
dimension input space.

A. Modeling smooth functions vs non-smooth function
In the following, we show the effectiveness of modeling
smooth functions by GPs and considering the non-smooth
operations in the BO search as opposed to modeling the nonsmooth function by a single GP.
Consider the following, illustrative optimization problem,
w∗ = argmin max(sin(w) + 0.65, cos(w) + 0.65)

(10)

w∈(0,10)

We consider two modeling scenarios, one where we
model max(sin(w), cos(w)) as a single GP, and another
where we model sin(w) by one GP and cos(w) by another. We initialize the GP models for sin(w), cos(w) and
max(sin(w), cos(w)) with 5 samples chosen in random. We
then use BO to find w∗ . We were able to model smooth
functions like sin(w) and cos(w) with GPs, even with fewer
samples. At each iteration of BO, we computed the next
sample by solving for the w ∈ (0, 10) which minimized
the maximum across the two GPs. This quickly stabilizes to
the true w∗ (Fig. 3c). When we model max(sin(w), cos(w))
using a GP, in Fig. 3b, the initial 5 samples were not able to
model it well. In fact, the original function in orange is not
contained within the uncertainty bounds of the GP. Hence,
in each iteration of BO, where we chose w ∈ (0, 10) which
minimized this function, we were never able to converge w∗ .
It is not surprising to see that, given these models, BO does
not always converge when we model non-smooth functions
such as in (10).
To support our claim, we repeat this experiment 15 times
with different initial samples. In each experiment we run BO
for 50 iterations. When modeling sin(w) and cos(w) as separate GPs, BO stabilized to w∗ in about 5 iterations in all 15
experiments. However, when modeling max(sin(w), cos(w))
as a single GP, it takes over 35 iterations to converge and in
5 out of the 15 cases, it did not converge to w∗ . We show
these two different behaviors in Fig. 4.
B. Collision Avoidance with High Dimensional Uncertainty
Consider an autonomous car that travels on a straight road
with a obstacle at xobs . We require that the car can come to
a stop before colliding with an obstacle. The car has two
states; location, x, and velocity, v; and one control input
acceleration; a. The dynamics of the car is given by,
ẋ = v,

v̇ = a.

(11)

Our safety specification for collision avoidance is given by,
ϕ = min(xobs − x(t)), i.e., the minimum distance between
the position of the car and the obstacle over a horizon of
length 100. We assume that the car does not know where
the obstacle is a priori, but receives locations of the obstacle
through a sensor at each time instant, xs (t). The controller
is a simple linear
K, such that at time
 state feedback control,
T
t, a(t) = K · x(t) − xs (t), v(t) .
We assume that the car initially starts at location xinit =
0, with a velocity vinit = 3 m/s. Let the obstacle be at
xobs = 5, which is not known by the car. Instead, it receives
sensor readings for the location of the obstacle such that

1.5

f(w)

1.0

2.5

2.0

2.0

1.5

1.5

1.0

1.0

0.5

0.5

0.5
0.0
−0.5

0

2

4

6

8

10

0.0

0.0

−0.5

−0.5

0

2

4

w

(a) True functions.

6

8

10

0

2

4

6

8

10

w

w

(b) Modeling the non-smooth ϕ(w) directly.

(c) Modeling and combining smooth predicates µ.

Fig. 3. The dashed orange line in Fig. 3a represents the true, non-smooth optimization function in (10) while the green and blue line represent sin(w)
and cos(w) respectively. Modeling this function directly as a GP leads to model errors Fig. 3b, where the 95% confidence interval of the GP (blue
shaded) with mean estimate (in blue line) does not capture the true function ϕ(w) in orange. In fact, the minimum (red star) is not contained within the
shaded region, causing the optimization to diverge. BO converges to the green dot, where ϕ(w) > 0 which is not a counterexample. Instead, modeling
the two predicates individually and combining them with the parse tree, leads to the model in Fig. 3c. Here, the true function is completely captured in
the confidence interval. As a consequence, BO converges to the global minimum (the red star and green dot converge).
10

Samples returned by BO

Samples returned by BO

10
8
6
4
2

8

6

4

0
0

10

20

30

40

50

0

10

20

BO iterations

(a) Modeling as separate GPs take around 5 iterations to stabilize to
w∗ (in blue), while modeling as a single GP takes around 45 iteration
to stabilize to w∗ (in orange)

30

40

50

BO iterations

(b) Modeling as separate GPs take around 5 iterations to stabilize to
w∗ (in blue), while modeling as a single GP does not stabilize to
w∗ (in orange)

Fig. 4. The orange and blue lines in Fig. 4a and Fig. 4b show the evolution of samples returned over the BO iterations when (10) is modeled as a single
GP and multiple GPs respectively for two different initialization. We see the that when modeling as a single GP, it takes longer to stabilize to w∗ and in
some cases (Fig. 4b) does not stabilize to w∗ .

200

Number of counterexamples

xs = [4.5, 5.5]. If ϕ is negative, then x(t) > xobs for
some t which signifies collision. Moreover, we constrain the
acceleration to lie in a ∈ [−3, 3].
The domain of our uncertainty is W = [4.5, 5.5]100 , i.e.,
the sensor readings xs over the horizon H = 100. We compare across three experimental setups, first, we model the GP
in the original space of W i.e., with 100 inputs; second, we
model the GP in a lower dimension input space as described
in the preamble of this section; and third, we randomly
sample inputs and test them. We run BO for 250 iterations on
the GPs, and consider 250 random samples for the random
testing. We repeat this experiment 10 times and show our
results in Fig. 5. The green and blue bar in Fig. 5 show
the average number of counterexamples returned running BO
on the GP defined over the original input space and in the
low dimension input space. In general, active testing in the
high-dimensional input space gives the best results, which
deteriorates with an increase in compression of the input
space. Random testing, shown in red performs the worst.
This is not surprising as, (1) 250 samples is not sufficient to
cover an input space of 100 dimensions uniformly; and (2)
the samples are all independent of each other. Moreover, in
the uncompressed input case, the specification evaluated at
the worst counterexample, ϕ(w∗ ), has a mean and standard
deviation of −0.0138 and 0.004 as compared to −0.0067
and 0.0011 for random sampling.

175
150
125
100
75
50
25
0
random

15

20

25

30 uncompressed

Fig. 5. The red, blue and green bars shows the average number of
counterexamples found using random sampling; applying BO on the reduced
input space and original input space respectively for the example in Sec. VB. The black lines show the standard deviation across the experiments.

C. OpenAI Gym Environments
We interfaced our tool with environments from OpenAI
gym [36] to test controllers from Open AI baselines [37].
For brevity, we refer the details of the environments to [38].
In both case studies, we introduce uncertainty around the
parameters the controller has been trained for. The rationale
behind this is that the parameters in a simulator are an estimate of the true values. This ensures that counterexamples
found, can indeed occur in the real system.
1) Reacher: In the reacher environment, we have a 2D
robot trying to reach a target. For this environment we

Number of counterexamples

120
100
80
60
40
20
0
1

2

3

4

5

6

7

8

9

10

Experiment number

Fig. 6. The green, blue and red bars show the number of counter examples
generated when modeling µ1 , µ2 as separate GPs; modeling ϕ as a single
GP and random testing respectively for the reacher example (Sec. V-C.1).
Our modeling paradigm, finds more counterexamples compared to the other
two methods.

have six sources of uncertainty: two for the goal position,
(xgoal , ygoal ) ∈ [−0.2, 0.2]2 , two for state perturbations
(δx , δy ) ∈ [−0.1, 0.1]2 and two for velocity perturbations
(δvx , δvy ) ∈ [−0.005, 0.005]2 . The state of the reacher is
tuple with the current location, x = (x, y), velocity v =
(vx , vy ), and rotation, θ. A trajectory of the system, ξ, is a
sequence of states over time, i.e., ξ = (x(t), v(t), θ(t)), t =
0, 1, 2, . . . . Our uncertainty space is, W = [−0.2, 0.2]2 ×
[−0.1, 0.1]2 ×[−0.005, 0.005]2 . Given an instance of w ∈ W,
the trajectory, ξ, of the system is uniquely defined.
We trained a controller using the Proximal Policy Optimization (PPO) [39] implementation available at Open AI
baselines. We determine a trajectory to be safe if either the
reacher reaches the goal, or if it does not rotate unnecessarily.
This can be captured as ϕ = µ1 ∨ µ2 , where, µ1 (w)
is the minimum distance between the trajectory and the
goal position, and µ2 is total rotation accumulated over the
trajectory; and its continuous variant, ϕ = max(µ1 , µ2 ).
Using our modeling approach, we model this using two
GPs, one for µ1 and another for µ2 . We compare this to
modeling ϕ as a single GP and random sampling. We run 200
BO iterations and consider 200 random samples for random
testing. We repeat this experiment 10 times. In Fig. 6, we
plot the number of counterexamples found by each of the
three methods over 10 runs of the experiment. Modeling
the predicates by separate GPs and applying BO across
them (shown in green) consistently performs better than
applying BO on a single GP modeling ϕ (shown in blue)
and random testing (shown in red). We see the that random
testing performs very poorly, in some cases (experiment runs
4, 8, 10) finds no counterexamples.
By modeling the predicates separately, the specification
evaluated at the worst counterexample, ϕ(w∗ ), has a mean
and standard deviation of −0.1283 and 0.0006 as compared
to −0.1212 and 0.0042 when considering a single GP. This
suggests, that using our modeling paradigm BO converges
(since the standard deviation is small) to a more falsifying
counterexample (since the mean is smaller).
2) Mountain Car Environment: The mountain car environment in OpenAI gym, is a car on a one-dimensional
track, positioned between two mountains. The goal is to
drive the car up the mountain on the right. The envi-

ronment comes with one source of uncertainty, the initial state xinit ∈ [−0.6, −0.4]. We introduced four other
sources of uncertainty, for the initial velocity, vinit ∈
[−0.025, 0.025]; goal location, xgoal ∈ [0.4, 0.6]; maximum
speed, vmax ∈ [0.55, 0.75] and maximum power magnitude,
pmax ∈ [0.0005, 0.0025]. The state of the mountain car
is a tuple with the current location, x, and velocity, v. A
trajectory of the system, ξ, is a sequence of states over time,
i.e., ξ = (x(t), v(t)), t = 0, 1, 2, . . . . Our uncertainty space
is given by, W = [−0.6, −0.4]×[−0.025, 0.025]×[0.4, 0.6]×
[0.55, 0.75] × [0.0005, 0.0025]. Given an instance of w ∈ W,
the trajectory, ξ, of the system is uniquely defined.
We trained two controllers one using PPO and another
using an actor critic method (DDPG) for continuous Deep
Q-learning [40]. We determine a trajectory to be safe, if it
reaches the goal quickly or if does not deviate too much
from its initial location and always maintains its velocity
in some bound. Our safety specification can be written as
ϕ = µ1 ∨ (µ2 ∧ µ3 ), where, µ1 (w) is time taken to reach the
goal, µ2 is the deviation from the initial location and µ3 is the
deviation from the velocity bound; and its continuous variant
of ϕ = max(µ1 , min(µ2 , µ3 )). We model ϕ, by modeling
each predicate, µ, by a GP. We compare this to modeling ϕ
with a single GP and random sampling. We run 200 BO
iterations for the GPs and consider 200 random samples
for random testing. We repeat this experiment 10 times.
We show our results in Fig. 7, where we plot the number
of counterexamples found by each of the three methods
over 10 runs of the experiment for each controller. Fig. 7
demonstrates the strength of our approach. The number of
counterexamples found by our method (in green bar) is much
higher compared to random sampling (in red) and modeling
ϕ as a single GP (in blue). In Fig. 7a the blue bars are smaller
than even the ones in red, suggesting random sampling
performs better than applying BO on the GP modeling ϕ.
The is because the GP is not able to model ϕ, and is so far
away from the true model, that the sample returned by the
BO is worse than if were to sample randomly.
This is further highlighted by the value of the specification
at worst counterexample, ϕ(w∗ ). The mean and standard
deviation for ϕ(w∗ ) over the 10 experiment runs is −0.5435
and 0.028 for our method, −0.3902 and 0.0621 when ϕ
is modeled as a single GP; and −0.04379 and 0.0596 for
random sampling. A similar but less drastic result holds in
the case of the controller trained with DDPG.
VI. C ONCLUSION
We presented an active testing framework that uses
Bayesian Optimization to test and verify closed-loop robotic
systems in simulation. Our framework handles complex logic
specifications and models them efficiently using Gaussian
Processes in order to find adversarial examples faster. We
showed the effectiveness of our framework on controllers
designed on OpenAI gym environments. As future work, we
would like to extend this framework to test more complex
robotic systems and find regions in the environment parameter space where the closed-loop control is expected to fail.

100

Number of counterexamples

Number of counterexamples

120
100
80
60
40
20
0

80
60
40
20
0

1

2

3

4

5

6

7

8

9

10

Experiment number

(a) Controller trained with PPO

1

2

3

4

5

6

7

8

9

10

Experiment number

(b) Controller trained with DDPG

Fig. 7. The green, blue and red bars show the number of counter examples generated when modeling µ1 , µ2 as separate GPs, modeling ϕ as a GP and
random testing respectively for the mountain car example (Sec. V-C.2). While our modeling paradigm, finds orders of magnitude more counterexample
compared to the other two methods, we notice that modeling ϕ as a single GP performs much worse than random sampling for the controller trained with
PPO Fig. 7a and comparable for the controller trained with DDPG Fig. 7b.

ACKNOWLEDGMENTS
Research reported in this paper was sponsored by the
Army Research Laboratory and was accomplished under
Cooperative Agreement Number W911NF-17-2-0196 1 and
was accomplished under Cooperative Agreement Number
W911NF-17-2-0196; and in part by Toyota under the iCyPhy
center.
R EFERENCES
[1] S. Sastry and M. Bodson, Adaptive Control: Stability, Convergence,
and Robustness. Prentice-Hall, Inc., 1989.
[2] R. F. Stengel, Stochastic Optimal Control: Theory and Application.
John Wiley & Sons, Inc., 1986.
[3] G. Kahn, A. Villaflor, V. Pong, P. Abbeel, and S. Levine, “Uncertaintyaware reinforcement learning for collision avoidance,” CoRR, vol.
abs/1702.01182, 2017.
[4] Y. Niv, D. Joel, I. Meilijson, and E. Ruppin, “Evolution of reinforcement learning in uncertain environments: A simple explanation for
complex foraging behaviors,” 2002.
[5] P. Poupart and N. Vlassis, “Model-based bayesian reinforcement
learning in partially observable domains,” in Proc Int. Symp. on
Artificial Intelligence and Mathematics,, 2008, pp. 1–2.
[6] E. M. Clarke, Jr., O. Grumberg, and D. A. Peled, Model Checking.
MIT Press, 1999.
[7] I. Mitchell and C. J. Tomlin, “Level set methods for computation
in hybrid systems,” in International Workshop on Hybrid Systems:
Computation and Control. Springer, 2000, pp. 310–323.
[8] A. Bhatia, L. E. Kavraki, and M. Y. Vardi, “Sampling-based motion
planning with temporal goals,” in Robotics and Automation (ICRA),
2010 IEEE International Conference on. IEEE, 2010, pp. 2689–2696.
[9] V. Raman, A. Donzé, M. Maasoumy, R. M. Murray, A. SangiovanniVincentelli, and S. A. Seshia, “Model predictive control with signal
temporal logic specifications,” in Decision and Control (CDC), 2014
IEEE 53rd Annual Conference on. IEEE, 2014, pp. 81–87.
[10] I. Mitchell, A. Bayen, and C. J. Tomlin, “Computing reachable sets
for continuous dynamic games using level set methods,” Submitted
January, 2004.
[11] F. Berkenkamp, M. Turchetta, A. P. Schoellig, and A. Krause, “Safe
model-based reinforcement learning with stability guarantees,” in
Proc. of Neural Information Processing Systems (NIPS), 2017.
[12] I. J. Goodfellow, J. Shlens, and C. Szegedy, “Explaining and harnessing adversarial examples,” arXiv preprint arXiv:1412.6572, 2014.
1 The views and conclusions contained in this document are those of the
authors and should not be interpreted as representing the ocial policies,
either expressed or implied, of the Army Research Laboratory or the U.S.
Government. The U.S. Government is authorized to reproduce and distribute
reprints for Government purposes notwithstanding any copyright notation
here on.

[13] N. Papernot, P. McDaniel, I. Goodfellow, S. Jha, Z. B. Celik, and
A. Swami, “Practical black-box attacks against deep learning systems
using adversarial examples,” arXiv preprint arXiv:1602.02697, 2016.
[14] V. Behzadan and A. Munir, “Vulnerability of deep reinforcement
learning to policy induction attacks,” CoRR, vol. abs/1701.04143,
2017.
[15] N. Carlini and D. Wagner, “Towards evaluating the robustness of
neural networks,” in Security and Privacy (SP). IEEE, 2017, pp.
39–57.
[16] S. H. Huang, N. Papernot, I. J. Goodfellow, Y. Duan, and
P. Abbeel, “Adversarial attacks on neural network policies,” CoRR,
vol. abs/1702.02284, 2017.
[17] T. Dreossi, A. Donzé, and S. A. Seshia, “Compositional falsification of
cyber-physical systems with machine learning components,” in NASA
Formal Methods Symposium. Springer, 2017, pp. 357–372.
[18] K. Pei, Y. Cao, J. Yang, and S. Jana, “Deepxplore: Automated
whitebox testing of deep learning systems,” in Proceedings of the 26th
Symposium on Operating Systems Principles, ser. SOSP ’17, 2017.
[19] A. Donzé, “Breach, a toolbox for verification and parameter synthesis
of hybrid systems,” in CAV. Springer, 2010, pp. 167–170.
[20] P. S. Duggirala, S. Mitra, M. Viswanathan, and M. Potok, “C2E2: A
verification tool for stateflow models,” in Tools and Algorithms for the
Construction and Analysis of Systems, C. Baier and C. Tinelli, Eds.,
2015.
[21] Y. Annpureddy, C. Liu, G. E. Fainekos, and S. Sankaranarayanan, “Staliro: A tool for temporal logic falsification for hybrid systems.” in
TACAS, vol. 6605. Springer, 2011, pp. 254–257.
[22] J. Deshmukh, X. Jin, J. Kapinski, and O. Maler, “Stochastic local
search for falsification of hybrid systems,” in ATVA. Springer, 2015.
[23] N. Hansen, “The CMA evolution strategy: A tutorial,” arXiv preprint
arXiv:1604.00772, 2016.
[24] J. Mockus, Bayesian approach to global optimization: theory and
applications. Springer Science & Business Media, 2012, vol. 37.
[25] C. Rasmussen and C. Williams, Gaussian Processes for Machine
Learning. MIT Press.
[26] N. Srinivas, A. Krause, S. M. Kakade, and M. Seeger, “Gaussian
process optimization in the bandit setting: no regret and experimental
design,” IEEE Transactions on Information Theory, vol. 58, 2012.
[27] F. Berkenkamp, A. Krause, and Angela P. Schoellig, “Bayesian optimization with safety constraints: safe and automatic parameter tuning
in robotics.” arXiv, 2016.
[28] J. V. Deshmukh, M. Horvat, X. Jin, R. Majumdar, and V. S.
Prabhu, “Testing cyber-physical systems through bayesian optimization,” Transactions on Embedded Computing Systems, 2017.
[29] R. Alur and D. L. Dill, “A theory of timed automata,” Theoretical
computer science, vol. 126, no. 2, pp. 183–235, 1994.
[30] A. Pnueli, “The temporal logic of programs,” in Proceedings of the
18th Annual Symposium on Foundations of Computer Science, ser.
SFCS ’77, 1977, pp. 46–57.
[31] O. Maler and D. Nickovic, “Monitoring temporal properties of continuous signals,” in FORMATS/FTRTFT. Springer, 2004.
[32] I. Steinwart and A. Christmann, Support vector machines. Springer
Science & Business Media, 2008.

[33] S. R. Chowdhury and A. Gopalan, “On kernelized multi-armed bandits,” in ICML, 2017, pp. 844–853.
[34] D. E. Finkel, “Direct optimization algorithm user guide,” 2003.
[35] Z. Wang, M. Zoghi, F. Hutter, D. Matheson, N. De Freitas et al.,
“Bayesian optimization in high dimensions via random embeddings.”
in IJCAI, 2013, pp. 1778–1784.
[36] G. Brockman, V. Cheung, L. Pettersson, J. Schneider, J. Schulman,
J. Tang, and W. Zaremba, “OpenAI Gym,” 2016.
[37] “Open AI Baselines,” https://github.com/openai/baselines.
[38] “Open AI Gym Environments,” https://gym.openai.com/envs.
[39] J. Schulman, F. Wolski, P. Dhariwal, A. Radford, and
O. Klimov, “Proximal policy optimization algorithms,” CoRR,
vol. abs/1707.06347, 2017.
[40] T. P. Lillicrap, J. J. Hunt, A. Pritzel, N. Heess, T. Erez, Y. Tassa,
D. Silver, and D. Wierstra, “Continuous control with deep reinforcement learning,” CoRR, vol. abs/1509.02971, 2015.

A PPENDIX
In this section, we prove the convergence of our algorithm
under specified regularity assumptions on the underlying
predicates. Consider the specification
ϕ(w) = T (µ1 (w), . . . , µq (w)),

(12)

where q represents the number of predicates. Let the domain
of the predicate indices be represented by, I = {1, . . . , q}.
The convergence proofs for classical Bayesian optimization
in [26], [33] proceed by building reliable confidence intervals
for the underlying function and then showing, that these confidence intervals concentrate quickly enough at the location
of the optimum under the proposed evaluation strategy. For
ease of exposition, we assume that measurements of each
predicate µi are corrupted by the same measurement noise.
To leverage these proofs, we need to account for the
fact that our GP model is composed of several individual
predicates and that we obtain one measurement for each
predicates at every iteration of the algorithm.
We start by defining a composite function
f : W × I → R, which returns the function values for
the individual predicates indexed by i.
f (w, i) = µi (w)

(13)

has bounded RKHS norm Bi . In this case, the function,
f (w, i) has RKHSPnorm with respect to the kernel in (14)
q
bounded by B = i Bi .
Lemma 1. Assume that f has RKHS norm bounded by
B and the measurementspare corrupted by σ-sub-Gaussian
1/2
noise. If βn·q = B + 4σ I(yq·(n−1) ; f ) + 1 + ln 1/δ, then
the following holds for all environment scenarios, w ∈ W,
predicate indices, i ∈ I, and iterations n ≥ 1 jointly with
probability at least 1 − δ,
1/2

i
(w, i)
|f (w, i) − miq·(n−1) (w, i)| ≤ βq·n σq·(n−1)

(15)

Proof. This follows directly from [27], which extends the
results from [33] and Lemma 5.1 from [26] to the case of
multiple measurements.
The scaling factor for the confidence intervals, βn·q ,
depends on the mutual information I(yq·(n−1) ; f ) between
the GP model of f and the q measurements of the individual
predicates that we have obtained for each time step so far.
It can easily be computed as
I(yq·(n−1) ; f ) = log(1 +
=

q
n−1
XX

1
Kq·(n−1) ),
σ2
2
log(1 + σj·q
(wj , i)/σ 2 ),

(16)

j=1 i=1

where Kq·(n−1) is the kernel matrix of the single GP over the
extended parameter space and the inner sum in the second
equation indicates the fact that we obtain q measurements at
every iteration.
Based on these individual confidence intervals on µ, we
can construct confidence intervals on ϕ. In particular, let
1/2

li (w) = mq·(n−1) (w, i) − βq·n σq·(n−1) (w, i)
1/2

ui (w) = mq·(n−1) (w, i) + βq·n σq·(n−1) (w, i)

(17)

The function f (·, ·) is a single output function, which can
be modeled with a single GP with a scalar output over the
extended input space, W × I. For example, if we assume
that the predicates are independent of each other, the kernel
function for f would look like,
(
ki (w, w0 ) if i = i0
k((w, i), (w0 , i0 )) =
,
(14)
0 otherwise

be the lower and upper confidence intervals on each predicate. From this, we construct reliable confidence intervals
on ϕ(w) as follows:

where ki is the kernel function corresponding to the GP
for the i−th predicate, µi . It is straightforward to include
correlations between functions in this formulation too.
This reformulation allows us to build reliable confidence
intervals on the underlying predicates, given regularity assumptions. In particular, we make the assumption that the
function f has bounded norm in the Reproducing Kernel
Hilbert Space (RKHS, [32]) corresponding to the same
kernel k(·, ·) that is used for the GP on f .

T (l1 (w), . . . , lq (w)) ≤ ϕ(w) ≤ T (u1 (w), . . . , uq (w))
(18)

Remark 1. Note, that this model is more general then the
case where we assume that each predicate, µi , individually

Lemma 2. Under the assumptions of Lemma 1. Let T be
the parse tree corresponding to ϕ. Then the following holds
for all environment scenarios, w ∈ W and iterations n ≥ 1
jointly with probability at least 1 − δ,

Proof. This is a direct consequence of Lemma 1 and the
properties of the min and max operators.
We are now able to prove the main theorem as a direct
consequence of Lemma 2.
Theorem 1. Assume that each predicate µi has RKHS
norm bounded by Bi and that the measurement noise is
σ-sub-Gaussian. Select δ ∈ (0, 1), wn according to (9), and

q
P
P
1/2
i
; µi ).
let
βn = i Bi + 4σ 1 + ln(1/δ) + i I(yn−1
If T (l1 (wn ), . . . , lq (wn )) > 0, then with probability at least
1 − δ we have that minw∈W ϕ(w) > 0 and the system has
been verified against all environments in W.

The max operator is analogous. Thus, since the parse tree T
is composed only of min and max nodes, the regret is
bounded by the maximum error over all predicates. The result
follows.

Proof. For independent variables the mutual information
decomposes additively and following Remark 1
this is a direct consequence of Lemma 2,
since T (l1 (w), . . . , lq (w)) ≤ ϕ(w) holds for all w ∈ W
with probability at least 1 − δ.

Lemma 4. Pick δ ∈ (0, 1) and βq·n as shown in Lemma 1,
then the following holds with probability at least 1 − δ,

A. Convergence proof
In the following, we prove a stronger result about convergence of our algorithm.
The key quantity in the behavior of the algorithm is the
mutual information (16). Importantly, it was shown in [27]
that it can be upper bounded by the worst-case mutual
information, the information capacity, which in turn was
shown to be sublinear by [26]. In particular, let fW denote the
noisy measurements obtained when evaluating the function f
at points in W. The mutual information obtained by the
algorithm can be bounded according to
I(fWn ×I ; f ) ≤
≤

max
W̄⊂W,|W̄|≤n

max

I(fW̄×I ; f );

D⊂W×I,|D|≤n·q

I(fD ; f );

(19)

= γq·n ,
where γn is the worst-case mutual information that we can
obtain from n measurements,
γn =

max

D⊂W×I,|D|=n

I(fD ; f ).

(20)

This quantity was shown to be sublinear in n for many
commonly-used kernels in [26].
A key quantity to show convergence of the algorithm is
the instantaneous regret,
rn = min ϕ(w) − ϕ(wn ),
w∈W

(21)

the difference between the unknown true minimizer of ϕ and
the environment parameters wn that Algorithm 1 selects at
iteration n. If the instantaneous regret is equal to zero, the
algorithm has converged.
In the
Pnfollowing, we will show that the cumulative regret,
Rn = i=1 ri is sublinear in n, which implies convergence
of Algorithm 1.
We start by bounding the regret in terms of the confidence
intervals on µi .
Lemma 3. Fix n ≥ 1, if |f (w, i) − mq·(n−1) (w, i)| ≤
1/2
βq·n σq·(n−1) (w, i) for all w, i ∈ W × I, then the regret
1/2
is bounded by rn ≤ 2βq·n maxi σq·(n−1) (w, i).
Proof. The proof is analogous to [26, Lemma 5.2]. The
maximum standard deviation follows from the properties of
the max and min operators in the parse tree T . In particular,
let a1 , b1 , a2 , b2 ∈ R with a1 − b1 < a2 − b2 . Then for
all c1 ∈ [−b1 , b1 ] and c2 ∈ [−b2 , b2 ] we have that
a1 − b1 ≤ min(a1 + c1 , a2 + c2 ) ≤ a1 + b1 .

(22)

n
X
i=1

rn2 ≤ βq·n C1 qI(fWn ×I ; f ) ≤ βq·n C1 γq·n

(23)

where rn is the regret between the true minimizing environment scenario, w∗ and the current sample, wn ; and
C1 = 8/ log 1 + σ −2
Proof. The first inequality follows similar to [26, Lemma
5.4] and the proofs in [27]. In particular, as in [27],
2
2
rn2 ≤ 4βq·n
(wn , i)
max σq·(n−1)
i∈I

The second inequality follows from (19).
Lemma 5. Under the assumptions of Lemma 2, let δ ∈ (0, 1)
and choose wn according to (9). Then, the cumulative regret
RN over N iterations of Algorithm 1 is bounded with high
probability,
n
o
p
Pr Rn ≤ C1 βN N γq·N ∀N ≥ 1 ≥ 1 − δ
(24)
where C1 =

8
log 1+σ −2 .

PN
Proof. Since, RN =
from Cauchy-Schwartz ini=1 ri , P
N
2
equality we have, RN ≤ N i=1 ri2 . The rest follows
from Lemma 4.
We introduce some notation, let
ŵn = argminw∈{w1 ,...,wn } ϕ(w)

(25)

be the minimizing environment scenario sampled by BO in
n iterations and let
w∗ = argmin ϕ(w)

(26)

w∈W

be the unknown, optimal parameter.
Corollary 1. For any δ ∈ (0, 1) and  ∈ R+ , there exits a
n∗ ,
n∗
C1
= 2
(27)
βq·n∗ γq·n∗

such that ∀n ≥ n∗ , ϕ(w∗ ) − ϕ(ŵn ) ≤  holds with
probability at least 1 − δ.

Proof.
Pn The∗ cumulative reward over n iterations, Rn =
i=1 ϕ(w ) − ϕ(wi ) where wi is the i-th BO sample.
Defining ŵn as in (25) we have,
Rn =

n
X
i=1

≥

n
X
i=1

ϕ(w∗ ) − ϕ(wi )
ϕ(w∗ ) − ϕ(ŵn )

= n(ϕ(ŵn ) − ϕ(w∗ ))

(28)

Combining this result with Lemma 5, we have with probability greater than 1 − δ that

Rn
n
r
C1 βq·n γq·n
≤
n
∗
To find, n , we bound the RHS by ,
r
C1 βq·n∗ γq·n∗
n∗
C1
≤

⇒
≥ 2
∗
n
βq·n∗ γq·n∗

ϕ(w∗ ) − ϕ(ŵn ) ≤

(29)

(30)

For n > n∗ , the minimum ϕ(ŵn ) ≤ ϕ(ŵn∗ ) =⇒ ϕ(ŵn )−
ϕ(w∗ ) ≤ .

We are now ready to prove our main convergence theorem.

Theorem 2. Under the assumptions of Lemma 2, choose
δ ∈ (0, 1),  ∈ R+ and define n∗ using Corollary 1. If
n ≥ n∗ and ϕ(ŵn ) > , then, with probability greater than
1 − δ, the following statements hold jointly
∗
• ϕ(w ) > 0
• The closed loop satisfies ϕ, i.e., the control can safely
control the system in all environment scenarios, W
• The system has been verified against all environments,
W
Proof. This holds from Lemma 5 and Corollary 1.
From Corollary 1, we have ∀n ≥ n∗ , Pr(ϕ(w∗ ) − ϕ(ŵn ) ≤
) > 1 − δ. If ∃n ≥ n∗ , such that ϕ(ŵn ) > , then we
have Pr(ϕ(w∗ ) > 0)1 − δ, i.e., the minimum value ϕ can
achieve on the closed loop system is greater than 0. ϕ is
hence, satisfied by our system in all w ∈ W.

