Active Sampling for Closed-loop Statistical Verification of
Uncertain Nonlinear Systems

arXiv:1705.01471v2 [] 1 Oct 2017

John F. Quindlen1

Ufuk Topcu2

Abstract— Increasingly demanding performance requirements for dynamical systems motivates the adoption of nonlinear and adaptive control techniques. One challenge is the
nonlinearity of the resulting closed-loop system complicates
verification that the system does satisfy the requirements at all
possible operating conditions. This paper presents a data-driven
procedure for efficient simulation-based, statistical verification
without the reliance upon exhaustive simulations. In contrast
to previous work, this approach introduces a method for online
estimation of prediction accuracy without the use of external
validation sets. This work also develops a novel active sampling
algorithm that iteratively selects additional training points in
order to maximize the accuracy of the predictions while still
limited to a sample budget. Three case studies demonstrate the
utility of the new approach and the results show up to a 50%
improvement over state-of-the-art techniques.

I. I NTRODUCTION
As the demand for higher performance and/or efficiency
grows, advanced nonlinear control techniques will be relied
upon to achieve the increasingly more complex requirements
associated with such demands. While methods like model
reference adaptive control (MRAC) [1] and reinforcement
learning (RL) [2] have demonstrated large improvements
in performance and efficiency, a key challenge with complex control techniques is certifying that the closed-loop
system can actually meet the requirements at all possible
uncertainties. The nonlinear (and possibly adaptive) nature
of the controller partly contributes to this difficulty since
the controller nonlinearities can result in drastically different
trajectories given only slightly different operating conditions.
Various verification techniques have been developed to
address this problem. If closed-form differential equations of
the closed-loop system are available, then it may be possible
to construct analytical certificates [3], [4] to provably verify
the closed-loop system satisfies the necessary requirements
under certain modeling assumptions. While analytical proofs
are extremely useful, they are difficult to implement on many
systems. Simulation-guided analytical methods [5]–[7] relax
this difficulty, but are still restricted by certain assumptions
and may be overly conservative due to the reliance upon
specific analytical functions for proof construction [8].
1 Ph.D.

Candidate, Department of Aeronautics and Astronautics, Massachusetts Institute of Technology (MIT).
2 Assistant Professor, Department of Aerospace Engineering and Engineering Mechanics, University of Texas at Austin
3 Assistant Professor, Departments of Agricultural & Biological Engineering
and Aerospace Engineering, University of Illinois Urbana-Champaign
4 Richard C. Maclaurin Professor of Aeronautics and Astronautics, Laboratory for Information and Decision Systems (LIDS), MIT

Girish Chowdhary3

Jonathan P. How4

Statistical verification techniques [9], [10] construct statistical certificates directly from simulations of the closedloop system. While these certificates do not suffer from
the same limitations as proof-based techniques and apply
to a wider range of systems, they are less absolute. If the
observed simulation data fails to adequately cover the entire
space of possible uncertainties, then the accuracy of the
certificate’s predictions in those unobserved regions will be
limited. However, it is prohibitively expensive to simply
saturate the set of all possible uncertainties with a large
number of simulations, particularly when the system requires
higher-fidelity verification models.
This work couples efficient data-driven statistical verification [11] with active sampling [12] to carefully select simulations in order to minimize utilization of the model while
maximizing prediction accuracy. In particular, Section III
introduces a Gaussian process (GP) verification framework to
train GP-based prediction models on a small set of trajectory
robustness measurements and estimate the satisfaction of the
requirements at unseen perturbations. In order to quantify
the model’s confidence in these predictions, the framework
also includes a simple method for online computation of
prediction confidence using cumulative distributions. Section IV exploits this validation metric to form a closedloop verification procedure that iteratively selects informative
training locations to improve prediction accuracy. As existing
active learning selection metrics [12]–[14] are not ideally
suited to this binary verification problem, we develop a new
tailor-made selection metric based on binary classification
entropy. Results in Section V demonstrate the new statistical
verification framework’s improvement in prediction error
when applied to three nonlinear systems.
II. P ROBLEM F ORMULATION
Consider the deterministic nonlinear system
ẋ(t) = f (x(t), u(t), θ)
p

(1)

subject to parametric uncertainties θ ∈ R , where x(t) ∈ Rn
is the state vector and u(t) ∈ Rm is the control input vector.
The open-loop plant in (1) is said to be deterministic; given
the same {x, u, θ}, the output will be the same every time.
Additionally, since the objective is to verify the closed-loop
system, the control inputs u(t) are assumed to be generated
by a known, deterministic control policy u(t) = g(x(t)).
The resulting closed-loop dynamics are a function of the
parametric uncertainties θ. These may include uncertainties
about initial state x(0) or uncertain parameters such as
mass/inertia properties of a physical system. Regardless of

their source, parameters θ are treated as uncertain initial
conditions that affect the state dynamics. Although the exact
values of the parameters may be unknown at run-time, they
are assumed to fall within a known, bounded set Θ.
Assumption 1 The set of all possible perturbations θ ∈ Θ
is a known, compact, uncountable set Θ ∈ Rp .
In practice, Assumption 1 is not overly restrictive as most
physical systems will have known feasible bounds on the
operating conditions that can be used as Θ.
The trajectory of the closed-loop system is given by
Φ(x(t)|x0 , θ) and defines the evolution of state vector x(t)
over the time interval t ∈ [0, Tf ]. Here, the nominal initial
state vector x0 is assumed to be fixed and known. If the
true initial condition x(0) is uncertain (x(0) 6= x0 ), then
x(0) can be modeled as the combination of known x0 and
corresponding elements of the unknown perturbations, i.e.
x(0) = x0 + θ. The resulting trajectory Φ(x(t)|x0 , θ) will
ultimately determine whether the performance requirements
are satisfied at a particular θ condition.
A. Problem Description
The closed-loop system trajectory is expected to satisfy
certain pre-specified performance criteria. These criteria may
be supplied by a wide variety of sources and can include
relatively straightforward concepts like stability or boundedness as well as more complex spatial-temporal requirements
such as those in Section V-B. This work assumes continuous
(non-binary) measurements indicate the robustness of the
trajectory to the requirements.
Assumption 2 Scalar variable y ∈ R measures the robustness of trajectory Φ(x(t)|x0 , θ) . The sign of y indicates
satisfaction of the requirements (y > 0 signifies the trajectory
satisfied the requirements, while y ≤ 0 indicates failure).
As fixed x0 and condition θ completely define the evolution
of each trajectory, we write the measurements as y(θ) to
emphasize robustness is an explicit function of θ.
Signal temporal logic (STL) [15] commonly, but not exclusively, provides these continuous measurements of robustness
to the requirements. Signal temporal logic is a mathematical
language for expressing the requirements as functions of
logical predicates and boolean and/or temporal operators.
In comparison to other temporal logic frameworks, STL
uniquely provides a continuous-valued robustness degree
ρϕ ∈ R to quantify the robustness of a trajectory with respect
to requirement ϕ. The availability of this ρϕ or a similar
measurement to quantify a trajectory’s robustness is central
to the approach in this paper. Similar work without the
availability of non-binary measurements has been presented
in [11], [16]. While it applies to a slightly larger class of
systems, the approach suffers from many limitations that will
be discussed later in this paper.
Even though continuous y(θ) measures the trajectory’s
robustness, the satisfaction of the requirements for deterministic systems is purely binary: at given condition θ, the

corresponding trajectory will either satisfy the requirements
or it will not. Due to the binary aspect of the problem, the
feasible set Θ can be segmented into two unique sets.
Definition 1 The region of satisfaction Θsat contains all
θ ∈ Θ for which the resulting trajectory satisfies the
performance requirements,
n
o
Θsat = θ ∈ Θ : y(θ) > 0 .
(2)
Definition 2 The region of failure Θf ail contains all θ ∈
Θ for which the resulting trajectory fails to satisfy the
requirements,
n
o
Θf ail = θ ∈ Θ : y(θ) ≤ 0 .
(3)
It is assumed Θsat , Θf ail 6= ∅. By construction, Θsat ∪
Θf ail = Θ and Θsat ∩ Θf ail = ∅. While the conditions
for membership in sets Θsat , Θf ail are known, the sets
themselves are unknown in advance; it is not clear whether
an arbitrary θ belongs to Θsat or Θf ail without running a
simulation. The verification goal is to compute an estimated
b sat .
region of satisfaction, Θ
Problem 1 Given a deterministic closed-loop system and
measurements of requirement satisfaction, compute an estib sat , with Θ
b f ail = Θ \ Θ
b sat .
mated region of satisfaction Θ
III. S TATISTICAL DATA -D RIVEN V ERIFICATION
Problem 1 can be viewed as a binary classification problem: predict whether θ belongs to Θsat or Θf ail . Previous
b sat and Θ
b f ail through support
work [11], [16] constructed Θ
vector machines (SVM) classification models. While the
results demonstrated the ability of the SVM-based approach
to product accurate estimates, there was no way to efficiently
quantify confidence in the predictions online without relying
upon external validation datasets. This mainly evolves out
of the fact binary classification models only utilize binary
evaluations rather than continuous y(θ). SVM extensions
such as Platt scaling [17] produce approximate measures
of confidence without validation sets, but these confidence
estimates rely upon artificial approximations of y(θ) output
from the SVM model rather than the true y(θ) values
themselves.
Although Assumption 2 indicates the sign of y(θ) can
always convert continuous y(θ) into binary evaluations, it
discards valuable information. In comparison to binary measurements, y(θ) also quantifies “just how close” a satisfactory trajectory came to failure or “how bad” an unsatisfactory
one performed. The main contribution of this section is a
new verification framework that trains a regression model in
order to avoid losing the information stored in the robustness degree. Even though regression models replace SVMs,
the problem is still fundamentally binary classification. To
emphasize this, we label the new approach regression-based
binary verification. With the extra information, regressionbased binary verification also simultaneously provides explicit measures of prediction confidence online without any
additional validation dataset.

A. Gaussian Process Regression Model
The noise-free GP regression model is constructed from a
training dataset L of initial observations. This training dataset
consists of N pairs of θ values and their measurements
y(θ). The set of N observed perturbation conditions is D =
{θ 1 , θ 2 , . . . θ N } while the measurements are grouped in
vector y = [y(θ 1 ), y(θ 2 ), . . . y(θ N )]T . The training dataset
is then L = {D, y}. At its core, a GP defines a distribution
over possible functions that predicts the value of y(θ) in
unobserved regions of Θ. A more in-depth discussion of GPs
is found in [18]. Assuming a zero-mean prior, the end result
of the training process is a Gaussian posterior predictive
distribution at arbitrary location θ ∗ with mean µ(θ ∗ ) and
covariance Σ(θ ∗ ),
µ(θ ∗ ) = KT∗ K−1 y
Σ(θ ∗ ) = κ(θ ∗ , θ ∗ ) − KT∗ K−1 K∗ .

(4)

Term K∗ is a N × 1 vector of κ(θ ∗ , θ i ) ∀i = 1 : N while
K is a N × N matrix of κ(θ i , θ j ) ∀i, j = 1 : N . Although
there are many possible options for kernel function κ, this
paper uses the squared exponential kernel with automatic
relevance determination (ARD) and kernel hyperparameters
ψ. In comparison to the isotropic squared exponential kernel,
the ARD kernel enables the hyperparameters to vary and thus
deemphasize elements of θ with minimal impact upon y(θ)
or emphasize those with high sensitivity.
The model presented in (4) is the standard GP model.
Although this is by far the most widely-used formulation, it
has been shown to have difficulty with higher-dimensional
systems. More complicated extensions [19], [20] decompose
the GP into the sum of additive models or find sparse
approximations of the full GP. For simplicity, we focus on
verification using standard GPs, but these more complex
extensions could be readily incorporated into our approach.
B. Importance of Hyperparameter Optimization
The choice in hyperparameters ψ greatly affects the GP
regression model and ultimately the posterior predictive disb sat , Θ
b f ail . Interestingly enough, the
tributions that define Θ
majority of works related to GP-based verification [14], [21]–
[23] make little-to-no mention of the choice of hyperparameters or assume the optimal choice of hyperparameters that
perfectly replicates the underlying function for y(θ) is known
and fixed. For Problem 1, this assumption cannot be made
as little is known about the distribution of the robustness
degree over θ until after simulations are performed. Instead,
the hyperparameters must be chosen with only the current
available information, training dataset L.
The likelihood distribution of the hyperparameters with respect to the training set L is given by P(ψ|L). The main issue
is that there is no closed-form solution to compute P(ψ|L)
and sampling-based approximations or Markov Chain Monte
Carlo (MCMC) methods are intractable. Instead, maximum
likelihood estimation (MLE) approximates the hyperparameters with local optimums. These methods efficiently compute
a local maximum using steepest descent methods or gradientbased methods [18].

C. Binary Predictions with Prediction Confidence
The posterior mean µ(θ) provides the basis for the binary
predictions. The mean represents the expected value of the
b sat ,Θ
b f ail
robustness degree, µ(θ) = yb(θ), and defines Θ
using the equivalent of Definitions 1 and 2 with yb(θ). Unlike
earlier work with SVM-based predictions [11], [16], the GPbased model also quantifies the prediction confidence without requiring an external validation dataset. The Gaussian
cumulative distribution function (CDF) associated with the
GP output determines the probability of satisfaction at each
θ location,
1 1  µ(θ) 
, (5)
P(y(θ) > 0|L, ψ) = P+ (θ) = + erf p
2 2
2Σ(θ)
and indicates the likelihood of misclassification error there.
The choice of hyperparameters will also affect the CDF
describing the probability of requirement satisfaction. If the
true hyperparameters, labeled ψ ∗ , are known, then the actual
probability P(y(θ) > 0|L) = P(y(θ) > 0|L, ψ ∗ ). However,
in the more likely case when ψ ∗ is not known, the probability
of satisfaction is the marginal likelihood over all possible ψ,
Z
P(y(θ) > 0|L) =

P(y(θ) > 0|L, ψ)P(ψ|L)dψ .

(6)

Since this integral requires P(ψ|L), the computation of the
total probability (6) is intractable as well and the MLE
approximation can be used in its place.
IV. C LOSED -L OOP S TATISTICAL V ERIFICATION
Data-driven verification must contend with the dueling
objectives of saturating Θ for better accuracy and minimizing |L| due to computational costs. This work assumes
the primary source of computational cost is the simulation
model. In many applications, the model fidelity required
for accuracy verification is quite high. For instance, aircraft
simulations typically require a full nonlinear flight simulation
model to capture the interactions of flight dynamics, actuator
saturation, and control and their effects on performance.
Even with lower-fidelity models, verification may be part
of a larger process such as robust nonlinear control design,
which restricts the number of simulations that can be feasibly
allocated to each candidate controller. To model these practical limits, we assume the computational budget restricts the
number of simulations allocated to the verification procedure.
Assumption 3 The computational budget manifests as a
cap on the number of simulations, Ntotal .
As the space Θ is uncountable, the approach assumes there
exists an extremely fine discretization Θd which approximates Θ. Even with this finite discretization, |Θd |  Ntotal
and it is impossible to saturate Θd with training locations
D. The remaining, unobserved sample locations are labeled
U = Θd \ D.

1) Active Sampling: With L limited to Ntotal samples,
it is crucial to select the most informative training set. As
expected, the difficulty lies in the fact that the ideal training
set is unknown apriori. Rather than passively select training
datapoints, either randomly or in a structured grid, active
learning selects informative training samples in an iterative
manner. Active learning [12] describes a wide variety of
different procedures, all of which attempt to identify the
“best” samples to obtain in order to improve a statistical
model. Depending on procedure’s objective, the definition of
the best sample will change, even for the same exact model.
One of the most general approaches is the expected model
change (EMC) metric [13] used by our prior work in closedloop verification with binary measurements [16]. Though
mostly applied to SVM models, EMC can be extended to
GP regression models. The “best” sample θ is the point
most likely to induce the largest expect change in the model,
meaning θ = argmin |µ(θ)|. Conversely, the most common
GP-based approach is variance reduction [14]. While various
approximations or derivations exist [22], [23], the end goal is
to reduces the posterior variance of the GP model. Since it is
comparatively expensive to calculate the change in posterior
variance, the most common approximation is to select θ =
argmax Σ(θ) [14]. Regardless of their efficiency or accuracy,
none of these approaches directly address the fundamental
verification objective: to predict whether y(θ) > 0 or not.
Closed-loop statistical verification introduces a new selection metric to identify informative sample locations: binary
classification entropy. Unlike variance methods, which may
also use the term “entropy” [14], binary classification entropy
exploits the probability P(y(θ) > 0),

H(θ|L, ψ) = − P+ (θ) log2 P+ (θ)

(7)
+ (1 − P+ (θ)) log2 (1 − P+ (θ)) ,
rather than the entropy of Σ(θ).
Ideally, the best possible sample location θ would be
to minimize the total posterior entropy. Unfortunately, it
is impractical to compute the expected posterior change in
entropy for the same reason as variance reduction methods
[14]. Instead, closed-loop statistical verification selects the
sample location with the largest reduction in local entropy
!
θ = argmax

H(θ 0 |L, ψ) − H(θ 0 |L+ , ψ)

θ0

= argmax H(θ 0 |L, ψ).

(8)

θ0

Once a sample is taken, the entropy at that location is 0;
therefore, the point with the largest magnitude of classification entropy would have the largest reduction. Note that the
binary classification entropy is strictly non-negative.
2) Batch Sampling: While a sequential sampling procedure based upon (8) will correctly guide the selection of
sample points as intended, it requires a large amount of
computational effort to continuously recompute the GP after

each iteration. Batching sampling presents one of the most
practical methods for reducing the retraining cost - select
multiple samples between retraining steps. Batch active
sampling methods [12], [13], [22], [23] lower the training
cost and exploit any inherent parallel computing capabilities.
These approaches select M datapoints at once and perform
their corresponding simulations in parallel. The challenge is
to ensure adequate diversity in the chosen batch set S in
order to avoid the selection of redundant trajectories.
The second major contribution of this paper is a batch active sampling extension of (8) for closed-loop statistical verification. This approach utilizes importance weighting, essentially modified importance sampling [24], to efficiently select
samples using only the current entropy. Unlike variancebased procedures, this approach freezes the information and
selects the entire batch of M samples without any intermediate retraining. The current entropy H(θ|L, ψ) is used to
construct a probability distribution from which samples of
θ can be chosen, where regions with a large magnitude of
entropy will have a higher probability of selection,
PH (θ) =

1
H(θ|L, ψ)
ZH

where ZH =

|θ d |
X
i=1

H(θ i |L, ψ).

(9)
Although pure random sampling based upon PH (θ) generates samples that are clustered in regions with large probability, it does not ensure adequate diversity in batch S.
Particularly when batch size M is low, it is desirable to
spread samples out across regions with similarly-high probability. The problem is that the randomly-selected samples
may inadvertently cluster into one area rather than spread
out, thus lacking diversity. In order to address this issue,
importance-weighting can be augmented with random matrix
theory methods to encourage diversity in the samples while
still steering samples towards regions with high probability/entropy.
The central tool to encourage diversity is Determinantal
Point Processes (DPPs) [25]. DPPs are probabilistic models for efficient sampling that penalize correlations in the
samples and therefore can be used to discourage similarities
between samples. Algorithm 1 presents the batch active
sampling procedure using DPPs. Each iteration’s DPP forms
a θ-correlation matrix L from a large number (MT ≈ 1000)
of θ locations randomly generated according to PH (θ).
Since MT must be so large to approximate PH (θ), but the
batch size M is generally small, Algorithm 1 implements
a variant of DPPs called k-DPPs [26] to correctly obtain
M  MT samples in a sequential manner. Once a sample has been selected, the remaining samples are weighted
towards points with corresponding eigenvectors of L that
are orthogonal to that recently-selected sample. The finished
batch S disperses the selected θ locations among regions of
high entropy/probability with less redundancy. It is important
to note that while the k-DPP requires MT random samples
from PH (θ), it does not actually perform simulations at all
those locations and only introduces a small additional cost
over baseline importance sampling.

degree ρϕ (θ) given by

Algorithm 1 Batch active sampling using k-DPP.
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:

Input: training set L, available sample locations U, T
batches, batch size M
Initialize: GP regression model, entropy H(θ|L, ψ)
for i = 1 : T do
Initialize: empty set S = ∅
Transform H(θ|L, ψ) into PH (θ)
Randomly generate MT locations, construct k-DPP
Obtain M samples according to k-DPP, add to S
Run simulations ∀θ ∈ S, obtain measurements yS
Add {S, yS } to L, remove S from U
Retrain regression model GP with new L
end for
b sat , Θ
b f ail
Return: predicted Θ

V. E XAMPLES
This section demonstrates Algorithm 1 on three example
problems of interest. Two of the examples consider adaptive
control systems. The third examines a lateral-directional
autopilot and demonstrates closed-loop verification on a very
complex system with nonlinear open-loop dynamics, control
saturation, sensor models, and switching logic.
A. Example 1: Model Reference Adaptive Control
The first example is a concurrent learning model reference
adaptive control (CL-MRAC) system [27]. This problem
has been examined before with other approaches, namely
proof-based techniques that produced overly-conservative
analytical certificates [8] and data-driven certificates with binary measurements [11], [16]. The following subsection will
demonstrate that entropy-based active sampling procedures
can produce accurate data-driven certificates of satisfaction
without relying upon large training datasets.
In this problem, the open-loop plant dynamics are subject
to two uncertain parameters, θ = [θ1 , θ2 ]T . The goal of
the system is to track a linear reference model with a
desirable second-order response (state xm (t)) over a 40
second trajectory. The tracking error, the distance between
the actual and reference trajectories, is given by e(t) =
xm (t)−x(t). The adaptive controller estimates the uncertain
parameters online and augments the nominal control input
with an additional component. As a result of the adaptive law
and the uncertainties, the closed-loop dynamics are nonlinear
and difficult to analyze.
While the CL-MRAC adaptive law ensures the asymptotic
convergence of the trajectory, it is difficult to prove boundedness of the transient errors. In particular, the performance
requirement considered in this example is that the tracking
error component e1 (t) remains bounded within the unit ball
at all times,
ϕbound = [0,40] (1 − |e1 [t]| ≥ 0).

ρϕ (θ) = 0min ρϕbound [t0 ](θ).
t ∈[0,40]

(11)

The goal of the active sampling procedures is to measure
y(θ) = ρϕ (θ) at various training locations and predict
P(y(θ) > 0) at unobserved θ ∈ Θd . For this example, Θd
covers the space between θ1 : [−10, 10] and θ2 : [−10, 10]
with a total of 40,401 possible sample locations. It is unlikely
Ntotal = 40, 401 in most computationally-constrained scenarios, so the samples must be carefully selected to improve
the predictions.
Figure 1 compares the verification performance of Algorithm 1 against batch versions of the existing EMC and
variance-based methods as well as passive, random sampling.
The procedures all start with the same 100 initial training sets
of 50 randomly-selected samples and operate in batches of
size M = 5 for a total of Ntotal = 350 points. Entropy-based
sampling outperforms all the competing techniques in average misclassification error. At the conclusion of the process
after |L| = 350 training points, Algorithm 1 demonstrates a
24.33% average improvement in misclassification error over
the closest competitor, variance-based sampling, as well as a
34.6% improvement over the EMC selection metric used by
closed-loop verification with binary evaluations [16]. Since
the four strategies all start from the same initial training
set and GP model in each of the 100 runs, the results
also directly compare their performance in each of those
test cases. Figure 2 plots the ratio of the 100 runs where
the misclassification rate of Algorithm 1 either matches or
exceeds that of the indicated competing strategy. The new
active sampling method matches or outperforms the topcompeting variance-based method in 88% of the runs.
Figure 3 illustrates the importance of hyperparameter
optimization when nothing is known about the correct choice
of hyperparameters. If the active sampling algorithms naively
fix the hyperparameters to their initial values and never
update them, then the prediction errors may actually increase
with additional simulation data. This mainly serves to highlight the importance of the choice of hyperparameters and
potential issues with sampling approaches [22], [23] that
assume fixed hyperparameters.
B. Example 2: Adaptive Control with Complex Temporal
Requirements
The second example verifies the same CL-MRAC system
from Section V-A, but with a different set of performance
requirements. In this example, the adaptive system must
satisfy three different signal temporal logic specifications,
i.e. ϕ = ϕ1 ∧ ϕ2 ∧ ϕ3 . These three specifications are
ϕ1 = ♦[2,3] (x1 [t] − 0.7 ≥ 0) ∧
♦[2,3] (1.3 − x1 [t] ≥ 0),
ϕ2 = ♦[12,13] (x1 [t] − 1.1 ≥ 0) ∧

(10)

♦[2,3] (1.7 − x1 [t] ≥ 0),

The performance measurement y(θ) is the STL robustness
and

ϕ3 = [22.4,22.6] (x1 [t] + 1.6 ≥ 0) ∧
[22.4,22.6] (−1.2 − x1 [t] ≥ 0) .

(12)

Misclassification Error (%)

Misclassification Error (%)

10
Open-loop: Random
Closed-loop: EMC
Closed-loop: Variance
Closed-loop: Entropy

8
6
4
2
50

100

150

200

250

300

20
Open-loop: Random
Closed-loop: EMC
Closed-loop: Variance
Closed-loop: Entropy

15

10

5

350

200

400

# Samples

100

60
50

Open-loop: Random
Closed-loop: EMC
Closed-loop: Variance

100

150

200

250

300

1000

100

80

Open-loop: Random
Closed-loop: EMC
Closed-loop: Variance

60
200

350

400

600

800

1000

# Samples

# Samples

Fig. 2: (Example 1) Percentage of the 100 runs where the entropybased procedure directly matches or outperforms the competing
approach given the same initial training dataset and GP model.
Misclassification Error (%)

800

Fig. 4: (Example 2) Total misclassification error at batch size M =
10 over 100 randomly-initialized runs. Standard deviation intervals
correspond to 0.5σ bounds.
Outperforms (%)

Outperforms (%)

Fig. 1: (Example 1) Total misclassification error at batch size M =
5 over 100 randomly-initialized runs. Standard deviation intervals
correspond to 0.5σ bounds for easier viewing.

80

600

# Samples

Fig. 5: (Example 2) Percentage of the 100 runs where the entropybased procedure directly matches or outperforms the indicated
approach given the same initial training dataset and GP model.

15
EMC: Hyp Opt
Variance: Hyp Opt
Entropy: Hyp Opt
EMC: Static Hyp
Variance: Static Hyp
Entropy: Static Hyp

10

5

0
50

100

150

200

250

300

350

Samples

Fig. 3: (Example 1) Average misclassification error of closed-loop
verification procedures with hyperparameter optimization versus the
same algorithms with static (constant) hyperparameters.

The total robustness measurement for each trajectory is then
y(θ) = min {ρϕ1 (θ), ρϕ2 (θ), ρϕ3 (θ)}.

(13)

This example considers the same two uncertain parameters
(θ1 , θ2 ) from Section V-A, but adds a third parameter θ3 for
uncertain initial state x1 (0). The new sampling set Θd covers
the space of θ1 : [−5, 5], θ2 : [−5, 5], and θ3 : [−1, 1] with a
grid of 214,221 possible sample locations. Just as with the
previous example, the four verification procedures start with
the same initial training dataset for each of the 100 runs.
In this example, these initial sets consist of 100 random θ
locations and the algorithms operate in batches of M = 10
points until they reach a total of Ntotal = 1, 050 samples.
Figure 4 displays the total misclassification error for all
four approaches while Figure 5 depicts the ratio of tests
where Algorithm 1 matches or outperforms the other three
strategies. The results reaffirm the conclusions from the
previous example: Algorithm 1 ultimately outperforms the

other approaches in at least 95% of tests, with a 33% and
35% improvement in mean misclassification error over the
EMC- and variance-based procedures.
Figure 6 highlights the use of the prediction confidence
(5) as an online validation tool to identify likely misclassification errors without external validation. Each θ queried
by the GP prediction model will return not just the binary
b sat or Θ
b f ail ), but also a confidence level
prediction (θ ∈ Θ
in that prediction. Not surprisingly, the GP model is more
likely to misclassify points with low confidence. Figure 6
displays the misclassification error rate for points with high
(≥ 95%) confidence. In comparison to Figure 4, the rate of
misclassification error is significantly lower, demonstrating
the GP prediction model correctly identified regions where
misclassifications are likely to occur. These results also show
that the prediction error is consistent across the various
open- and closed-loop verification procedures. This last fact
confirms the utility of GP-based, data-driven verification
and its prediction confidence value, regardless of the exact
sampling strategy.
C. Example 3: Lateral-Direction Autopilot
The third example examines an aircraft autopilot for
controlling lateral-directional flight modes. In particular, the
“heading-hold” autopilot turns the aircraft to and holds a
desired reference heading. The simulation model of the
autopilot and DeHavilland Beaver airframe is provided by
the “Aerospace Blockset” toolbox in Matlab/Simulink [28].
This simulation model includes numerous nonlinearities such
as the full nonlinear aircraft dynamics, nonlinear controllers,
and actuator models with position and rate saturations. The
provided “heading hold” autopilot then must satisfy several

Misclassification Error (%)

Misclassification Error (%)

6
Open-loop: Random
Closed-loop: EMC
Closed-loop: Variance
Closed-loop: Entropy

5
4
3
2
1
200

400

600

800

1000

12

Open-loop: Random
Closed-loop: EMC
Closed-loop: Variance
Closed-loop: Entropy

10
8
6
4
2
100

200

# Samples

ϕheight = [0,Tf ] (35 − |x[t] − x[0]| ≥ 0),

(14)

where Tf = 50 seconds is the final simulation time and x[t]
is the aircraft altitude (in feet) at time t. The performance
metric is similar to the last example with
t ∈[0,Tf ]

0

[t ](θ).

(15)

The verification procedures test the satisfaction of ϕheight
against different conditions of the initial Euler angles and
longitudinal inertia, θ =[roll(0), pitch(0), heading(0), Iyy ]T .
While satisfaction of other autopilot requirements can be
explored, such as heading angle overshoot or steady-state
tracking error, the altitude-hold requirement dominated the
other requirements during an initial trade-space exploration.
The space of allowable perturbations Θ spans roll(0):
◦
◦
◦
◦
◦
◦
[−60 , 60 ], pitch(0): [4 , 19 ], heading(0): [75 , 145 ], and
inertia Iyy : [5430, 8430] (kg · m2 ). The desired reference
◦
heading was kept constant at 112 . Set Θd discretizes Θ
into a normalized 4D grid with 937,692 possible locations.
This example uses the same general procedure from the
last two examples. During each of the 100 runs, the four
procedures start with the same randomly-generated initial
training set of 100 samples. The closed-loop procedures
then select 400 additional training samples in batches of
M = 5 points. As seen in Figure 7, the entropy-based
procedure outperforms the competing approaches to even
greater extent than preceding examples. The entropy-based
procedure produces below 3% average prediction error after
500 iterations, a 37% improvement over the closest variancebased method and up to 50% over the other two. Similarly,
the results in Figure 8 show that entropy-based, closedloop verification will ultimately produce a lower prediction
error than the competing algorithms 100% of the time.
These results further highlight the consistently desirable

Outperforms (%)

performance requirements.
The requirement chosen for this example is the “altitudehold” requirement of the heading autopilot. In addition to
turning the aircraft to the desired heading angle, the autopilot
must also ensure the aircraft’s altitude remains within 35 feet
of the initial altitude when the reference command was given
[28]. This requirement can be written as

y(θ) = 0 min ρ

400

500

Fig. 7: (Example 3) Total misclassification error at batch size M =
5 over 100 randomly-initialized runs. Standard deviation intervals
correspond to 0.5σ bounds.
100

80

60
100

Open-loop: Random
Closed-loop: EMC
Closed-loop: Variance

200

300

400

500

# Samples

Fig. 8: (Example 3) Percentage of the 100 runs where the entropybased procedure directly matches or outperforms the indicated
approach given the same initial training dataset and GP model.
Misclassification Error (%)

Fig. 6: (Example 2) Misclassification error within the 95% prediction confidence levels. Unlike the total misclassification error,
this only counts prediction errors for points with high confidence.
Standard deviation intervals correspond to 0.5σ bounds.

ϕheight

300

# Samples

3.5
Open-loop: Random
Closed-loop: EMC
Closed-loop: Variance
Closed-loop: Entropy

3
2.5
2
1.5
1
0.5
100

200

300

400

500

# Samples

Fig. 9: (Example 3) Misclassification error within the 95% prediction confidence levels. Unlike the total misclassification error,
this only counts prediction errors for points with high confidence.
Standard deviation intervals correspond to 0.5σ bounds.

performance of Algorithm 1 for deterministic verification of
complex nonlinear systems.
Just like Figure 6, Figure 9 demonstrates the ability
of prediction confidence (5) to correctly identify points
with low accuracy. Once these points are removed, the
misclassification error rate drops considerably, regardless of
the sampling strategy. Where as Figures 7 and 8 highlight
the desirable impact of closed-loop statistical verification,
Figure 9 reaffirms the utility of the online validation aspect
of the new GP-based statistical verification framework.
VI. C ONCLUSION
This paper presented new data-driven verification methods for simulation-based certificates of complex nonlinear
systems. In particular, we introduced a GP-based verification
framework to predict the satisfaction of requirements over the
full space of possible operating conditions. Additionally, this

new approach exploits the availability of continuous measurements to quantify prediction confidence without external
validation sets. In many applications, the simulations themselves can be computationally expensive to obtain; therefore,
it is advantageous to minimize the number of simulations
required to obtain accurate predictions. The second contribution of the paper is closed-loop verification using binary
classification entropy for active selection of future training
simulations. Using this strategy, the examples in Section V
demonstrated up to a 50% improvement in prediction error
over existing approaches for the same number of samples.
Ultimately, data-driven verification procedures are intended for use within higher-level problems such as robust,
nonlinear planning or controller optimization. In those those
problems, verification is performed on each candidate control
policy in order to estimate their robustness and the generation
of candidate control policies is often an iterative process.
Thus, the process would typically cycle through a large number of candidate control policies and it is infeasible to test
all possible uncertainties for every candidate policy. Closedloop, data-driven verification aligns with those problems by
providing the controls engineer with the best approximation
of the robustness while restricted to a budget on the number
of allowable simulations during each verification step of the
iterative process. Upcoming work [29] has further developed
the closed-loop verification procedures and extended them to
stochastic systems.
ACKNOWLEDGMENTS
This work is supported by the Office of Naval Research
and Air Force Office of Scientific Research under grants
ONR MURI N00014-11-1-0688 and AFOSR FA9550-15-10146 .
R EFERENCES
[1] E. Lavretsky and K. A. Wise, Robust and Adaptive Control, M. Grimble and M. Johnson, Eds. Springer, 2013.
[2] R. S. Sutton and A. G. Barto, Reinforcement Learning: An Introduction. The MIT Press, 1998.
[3] J. Moore and R. Tedrake, “Control synthesis and verification for a
perching UAV using LQR-trees,” in IEEE Conference on Decision
and Control, 2012.
[4] S. Prajna, “Optimization-based methods for nonlinear and hybrid
systems verification,” Ph.D. dissertation, California Institute of Technology, 2005.
[5] J. Kapinski, J. Deshmukh, S. Sankaranarayanan, and N. Arechiga,
“Simulation-guided Lyapunov analysis for hybrid dynamical systems,”
in Hybrid Systems: Computation and Control, 2014.
[6] U. Topcu, “Quantitative local analysis of nonlinear systems,” Ph.D.
dissertation, University of California, Berkeley, 2008.
[7] P. Reist, P. V. Preiswerk, and R. Tedrake, “Feedback-motion-planning
with simulation-based LQR-trees,” 2015, under review.
[8] J. F. Quindlen, U. Topcu, G. Chowdhary, and J. P. How, “Regionof-convergence estimation for learning-based adaptive controllers,” in
American Control Conference, 2016.

[9] E. M. Clarke and P. Zuliani, “Statistical model checking for cyberphysical systems,” in International Symposium for Automated Technology for Verification and Analysis, 2011.
[10] Y. Zhang, S. Sankaranarayanan, and F. Somenzi, “Statistically sound
verification and optimization for complex systems,” in Automated
Technology for Verification and Analysis, Sydney, Australia, November
2014, pp. 411–427.
[11] A. Kozarev, J. F. Quindlen, J. How, and U. Topcu, “Case studies in
data-driven verification of dynamical systems,” in Hybrid Systems:
Computation and Control, 2016.
[12] B. Settles, Active Learning, R. J. Brachman, W. W. Cohen, and T. G.
Dietterich, Eds. Morgan and Claypool, 2012.
[13] J. Kremer, K. S. Pedersen, and C. Igel, “Active learning with support
vector machines,” Data Mining and Knowledge Discovery, vol. 4,
no. 4, pp. 313–326, July 2014.
[14] Y. Zhang, T. N. Hoang, K. H. Low, and M. Kankanhalli, “Nearoptimal active learning of multi-output gaussian processes,” in AAAI
Conference on Artificial Intelligence, 2016.
[15] O. Maler and D. Nickovic, Monitoring Temporal Properties of Continuous Signals. Springer Berlin Heidelberg, 2004, pp. 152–166.
[16] J. F. Quindlen, U. Topcu, G. Chowdhary, and J. P. How, “Active
sampling-based binary verification of dynamical systems,” in AIAA
Guidance, Navigation, and Control Conference, 2018. [Online].
Available: http://arxiv.org/abs/1706.04268
[17] J. Platt, “Probabilistic outputs for support vector machines and comparisons to regularized likelihood methods,” in Advances in Large Margin
Classifiers, 1999, pp. 61–74.
[18] C. E. Rasmussen and C. K. I. Williams, Gaussian Processes for
Machine Learning, T. Dietterich, Ed. MIT Press, 2006.
[19] K. Kandasamy, J. Schneider, and B. Poczos, “High dimensional
bayesian optimization and bandits via additive models,” in International Conference on Machine Learning, 2015.
[20] T. N. Hoang, Q. M. Hoang, and K. H. Low, “A unifying framework
of anytime sparse gaussian process regression models with stochastic
variational inference for big data,” in International Conference on
Machine Learning, 2015.
[21] G. Chen, Z. Sabato, and Z. Kong, “Active learning based requirement
mining for cyber-physical systems,” in IEEE Conference on Decision
and Control, 2016.
[22] A. Gotovos, N. Casati, and G. Hitz, “Active learning for level set estimation,” in International Joint Conference on Artificial Intelligence,
2013, pp. 1344–1350.
[23] T. Desautels, A. Krause, and J. W. Burdick, “Parallelizing explorationexploitation tradeoffs in gaussian process bandit optimization,” Journal
of Machine Learning Research, vol. 15, pp. 4053–4103, 2014.
[24] R. Y. Rubinstein and D. P. Kroese, Simulation and the Monte Carlo
Method. Wiley, 2008.
[25] A. Kulesza and B. Taskar, “Determinantal point processes,” Foundations and Trends in Machine Learning, vol. 5, no. 2-3, pp. 123–286,
2012.
[26] ——, “k-DPPs: Fixed-size determinantal point processes,” in International Conference on Machine Learning, 2011.
[27] G. Chowdhary and E. Johnson, “Concurrent learning for convergence
in adaptive control without persistency of excitation,” in IEEE Conference on Decision and Control. IEEE, 2010.
[28] C. Elliott, G. Tallant, and P. Stanfill, “An example set of cyber-physical
V&V challenges for S5,” in Air Force Research Laboratory Safe and
Secure Systems and Software Symposium (S5) Conference, Dayton,
OH, July 2016.
[29] J. F. Quindlen, U. Topcu, G. Chowdhary, and J. P. How. (2017)
Closed-loop statistical verification of stochastic nonlinear systems
subject to parametric uncertainties. Preprint. [Online]. Available:
https://arxiv.org/abs/1709.06645

