Active Sampling-Based Binary Verification of
Dynamical Systems

arXiv:1706.04268v2 [] 16 Jan 2018

John F. Quindlen1∗

Ufuk Topcu2†

Girish Chowdhary3‡

Jonathan P. How4§

Nonlinear, adaptive, or otherwise complex control techniques are increasingly relied
upon to ensure the safety of systems operating in uncertain environments. However, the
nonlinearity of the resulting closed-loop system complicates verification that the system
does in fact satisfy those requirements at all possible operating conditions. While analytical proof-based techniques and finite abstractions can be used to provably verify the
closed-loop system’s response at different operating conditions, they often produce conservative approximations due to restrictive assumptions and are difficult to construct in
many applications. In contrast, popular statistical verification techniques relax the restrictions and instead rely upon simulations to construct statistical or probabilistic guarantees.
This work presents a data-driven statistical verification procedure that instead constructs
statistical learning models from simulated training data to separate the set of possible
perturbations into “safe” and “unsafe” subsets. Binary evaluations of closed-loop system
requirement satisfaction at various realizations of the uncertainties are obtained through
temporal logic robustness metrics, which are then used to construct predictive models of
requirement satisfaction over the full set of possible uncertainties. As the accuracy of
these predictive statistical models is inherently coupled to the quality of the training data,
an active learning algorithm selects additional sample points in order to maximize the
expected change in the data-driven model and thus, indirectly, minimize the prediction error. Various case studies demonstrate the closed-loop verification procedure and highlight
improvements in prediction error over both existing analytical and statistical verification
techniques.

I.

Introduction

Verification is crucially important for nonlinear and adaptive control systems as they are expected to safely
operate in complex, uncertain environments. In many applications, the control system must ensure the overall
system adheres to a wide range of certification criteria. These criteria include not only closed-loop stability,
but other spatial-temporal performance metrics such as aircraft handling and maneuver requirements.1, 2
While the entire purpose of the controller is to force the system to satisfy the necessary requirements,
the inherent nonlinearity of an adaptive or nonlinear control approach complicates analysis of the closedloop system’s response. Given slightly different operating conditions, the nonlinearities in the closed-loop
dynamics may cause the system to produce drastically different trajectories. In order to ensure successful
operation of the system at different possible operating conditions, verification procedures must demonstrate
that the closed-loop dynamical system complies with those certification criteria.
Various verification procedures exist, but they typically fall within two general categories: analytical
proof- or abstraction-based methods and statistical techniques. If closed-form differential or difference equations of the model are known, then it may be possible to construct analytical certificates that provably verify
the stability or boundedness of the closed-loop dynamical system.3–5 These certificates are based upon analytical functions such as Lyapunov, barrier, or storage functions and separate the bounded set of all possible
perturbations into two sets: those from which the system is proven to successfully meet the requirements
∗1 Graduate

Student, Department of Aeronautics and Astronautics, Massachusetts Institute of Technology
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, Aerospace Controls Laboratory (Director), MIT
†2 Assistant

and those from which the system cannot be proven to meet the requirements. These analytical certificates
are extremely powerful tools for verifying the safety of a control system, but are difficult to implement on
many systems.
One problem with these analytical certificates is that an appropriate certificate is not always obvious nor
easy to construct. Even with recent simulation-guided certificate construction techniques,6–8 it is often difficult to construct a certificate that accurately approximates the true (but unknown) set of all safe operating
conditions. The main limitation is that these certificates are dependent upon the analytical function used
during construction. When the function that best captures the full set of safe perturbations is unknown in
advance, then the reliance upon a less-suitable one will result in suboptimal estimates of the true set. In
some applications, such as adaptive controllers, only relatively simple functions are currently available for
verification, even though results9 suggest more complex functions are needed to less-conservatively capture
the true set of safe conditions.
Statistical verification procedures10–14 do not suffer the same limitations as analytical methods since
they utilize simulations to directly construct statistical or probabilistic estimates for requirement satisfaction. While statistical predictions are weaker than analytical guarantees, statistical techniques do not require
the same restrictive assumptions and apply to a wider class of systems. In particular, this work presents
data-driven statistical verification techniques15 as an alternative to analytical barrier certificates. Datadriven verification applies machine learning techniques to construct statistical certificates of requirement
satisfaction. At its core, this approach interprets verification as a binary data-classification problem. Statistical models train on simulations of the closed-loop system and classify potential operating conditions as
either “safe,” meaning the resulting trajectory will ultimately satisfy the requirement, or “unsafe,” meaning
it will not. In this manner, data-driven verification possesses the same applicability of analytical certificates without the reliance upon a particular analytical function. The removal of this restriction also enables
data-driven verification to apply to a much larger class of systems where analytical Lyapunov-like functions
are not applicable. In fact, statistical data-driven methods can be applied to virtually any system with a
suitable simulation model or black-box oracle that is able to provide binary safe/unsafe labels for training
queries.
While data-driven verification methods are able to relax the conservativeness of analytical certificates
and apply to a wider class of problems, the accuracy of their predictions is fundamentally limited by the
quality of the simulation data. Predictions in regions without adequate coverage by training data will have
limited accuracy. The most straightforward solution to improve the quality of the predictions is simply to
saturate the space of possible conditions with a large, structured training dataset using design of experiments
techniques.16, 17 While this would help reduce the number of misclassification errors, it neglects the fact
simulations can be expensive to obtain. Especially with a high-fidelity model, each simulation of the closedloop system may be computationally expensive to obtain, and thus it is impractical and inefficient to blindly
saturate the space with a large training dataset. Instead, it is best to only select the most informative
training points. The problem is that the ideal training set is unknown until after a large set of simulation
data has already been obtained.
This work incorporates active learning18–22 into data-driven verification in order to improve the informativeness of the training data without relying upon an exceedingly large dataset. The expected model change
metric19, 21 identifies unobserved perturbations which would induce the greatest expected change upon the
current statistical classification model. The resulting sampling algorithms describe two closed-loop verification procedures to maximize the expected informativeness of the training dataset for a designated number
of simulations. These new closed-loop verification procedures directly addresses the limitations of earlier
data-driven verification work,15 which relied upon passive, open-loop generation of large training datasets.
The paper is structured as follows. As background, the next two sections provide an overview of the
certification problem and discuss the formulation and limitations of analytical verification techniques. These
limitations motivate the data-driven approach towards verification in Section IV. Section V presents the
closed-loop sampling procedures that iteratively update the data-driven verification model and improves its
accuracy. Subsequent case studies of nonlinear and adaptive control systems demonstrate the improvement
in accuracy and sample efficiency.

II.

Problem Formulation

Consider the deterministic nonlinear system
ẋ(t) = f (x(t), u(t), θ)

(1)

subject to uncertain operating conditions θ ∈ Rp , where x(t) ∈ Rn is the state vector, X is the set of all
feasible states, u(t) ∈ Rm is the control input vector, and Xu is the set of all feasible control inputs. In this
problem, the goal is to verify whether a particular control policy is able to ensure the closed-loop system
meets all the necessary performance requirements; therefore, the controller which generates inputs u(t) is
assumed to be known and given a priori.
The parametric uncertainties θ include all the different operating conditions of the system. These uncertainties are treated as time-invariant parameters and may arise from a variety of sources. For example,
the perturbation vector θ may include uncertainties about the initial state x(0) or system variabilities like
vehicle mass and inertia. Regardless of their source, the parametric uncertainties θ are considered “known
unknowns” and are assumed to fall within known bounded set Θ.
Assumption 1 The set of all possible parametric uncertainties θ ∈ Θ is assumed to be a known,
compact, and uncountable set Θ ⊂ Rp .
For instance, aircraft mass and C.G. location will not be perfectly known during flight, but the range of
values in which they will lie can be computed in advance. Thus, verification will perform simulations at
various conditions to check whether the aircraft will meet the necessary requirements for all possible mass
and loading conditions.
The trajectory of the system is given by Φ(x(t)|x0 , θ), which defines the time evolution of the state
vector x(t) over the time interval t ∈ [0, Tf inal ]. The nominal initial state vector x0 is assumed to be known
and fixed while the perturbation vector θ is a particular instantiation from the full set Θ. As previously
mentioned, the actual initial state x(0) may vary, but x(0) can be modeled as the sum of the fixed, nominal
x0 and varying θ terms. Given the deterministic system in (1), the closed-loop system’s trajectory response
is completely determined by these two terms x0 and θ. Since x0 is fixed, the satisfaction of the performance
requirements will only vary due to changes in θ.
A.

Requirements Modeling

For simulation-based verification, the performance requirements are assumed to be given ahead of time by
relevant certification authorities such as the FAA. These certification requirements may include relatively
straightforward considerations such as stability of the state response, i.e. x(t)→∞ as t→∞, and physical
bounds on the allowable states such as never-exceed speed of an aircraft. The criteria might also include
more complex spatial-temporal requirements like the avoidance of certain states x(t) at particular times. In
order to capture the wide range of possible criteria of interest, it is assumed the requirements can be written
with temporal logic specifications. In short, temporal logic simply provides a mathematical framework for
defining requirements and determining whether trajectory Φ(x(t)|x0 , θ) satisfies those requirements. The
following material will briefly discuss temporal logic and highlight its application to the simulation-based
binary verification problem.
A requirement is specified by a temporal logic formula ϕ. This formula consists of predicate(s) ζ, which are
functions of the state and/or control input, as well as boolean and temporal operations on those predicates.
Satisfaction of the temporal logic predicate ζ is signified by inequality ζ > 0. For instance, if a requirement
states that the state x1 (t) must remain below x1 = 2, then the predicate at time ti is ζ[ti ] = 2 − x1 (ti )
and the satisfaction of the requirement corresponds to ζ[ti ] > 0. Temporal operations on the predicates are
used to encapsulate time dependencies in the requirements. The three common operators are [t1 ,t2 ] , ♦[t1 ,t2 ] ,
and U[t1 ,t2 ] , which express that a predicate must hold “for all” time within interval [t1 , t2 ], “at some point”
between t1 and t2 , and “until” another formula is true within [t1 , t2 ]. Meanwhile, the boolean operators ¬,
∧, and ∨ can be used alongside temporal operators to express negation, conjunction, and disjunction. One
of the most important points about temporal logic is that these boolean operators can be used to construct
more complex formula from simpler ones. For example, formula ϕ3 = [t1 ,t2 ] ϕ1 ∧ ♦[t2 ,t3 ] ϕ2 states that ϕ1
must hold for all times between t1 and t2 and ϕ2 must occur at some point between t2 and t3 in order for
the formula ϕ3 to be satisfied.

Once a trajectory Φ has been generated, a comparison of the state response to the specifications in ϕ
determines whether or not the trajectory satisfies the requirement. Tuple (Φ, t) |= ϕ signifies that trajectory
Φ satisfies requirement formula ϕ at time t, while (Φ, t) |= ¬ϕ signifies failure to meet the criteria. This
paper uses metric temporal logic (MTL)23 to produce binary {−1, 1} measurements of satisfaction from
these boolean results. If the trajectory satisfies the requirement, i.e. (Φ, t) |= ϕ, then measurement y = 1,
else y = −1 when the trajectory fails to satisfy the requirement in ϕ. The measurements are rewritten
as y(θ) ∈ {−1, 1} to emphasize the trajectory and the corresponding satisfaction of the requirement is an
explicit function of θ and only changes with θ. While this paper assumes binary measurements, parallel
work24 addresses the case when continuous measurements of satisfaction are available. These continuous
measurements, obtained using signal temporal logic (STL), provide additional information which can be used
to approach the problem with a different data-driven verification method. While this additional information
allows for improved quantification of prediction confidence over the methods that will be shown in Section
IV, STL is a more restrictive assumption and is not available in all scenarios. Additionally, active samplingbased verification with binary/discrete measurements can be applied to a much larger class of problems,
such as black-box or even human-based labeling of requirement satisfaction.
B.

Problem Description

The goal of the verification procedure is to predict whether a trajectory initialized with the queried θ
vector will satisfy the requirements or not. This binary result naturally leads to the following two sets: set
Θsaf e ⊂ Θ, which contains all θ corresponding to y(θ) = 1, and its complement, unsafe set Θf ail ⊂ Θ
that contains all θ corresponding to y(θ) = −1. The verification problem can thus be viewed as a binary
classification problem, where the goal is to identify whether queried θ vectors belongs to Θsaf e or Θf ail .
Problem 1 The objective is to provide a classifier that separates Θ into estimates of Θsaf e and Θf ail .
By definition, Θsaf e ∪ Θf ail = Θ and Θsaf e ∩ Θf ail = ∅. It is assumed Θsaf e 6= ∅ and Θf ail 6= ∅.
While Assumption 1 states the true set Θ is uncountable, queries can only be made at discrete locations
and therefore a simulation-based analysis of Θ will only ever be able to observe a countable subset of Θ.
With that in mind, it is assumed there is a particular countable approximation of Θ that should be used for
statistical verification.
Assumption 2 There exists a sufficiently-fine discretization of the uncountable set Θ, called Θd , which
is acceptable for verification purposes.
During the data-driven verification process, simulations will be performed at locations θ ∈ Θd . While
the resolution of a sufficiently-fine discretization will change according to the application and possibly the
certification authority (such as the FAA), it is assumed the size |Θd | is so large it will be infeasible to perform
simulations at every θ ∈ Θd . To capture this effect, there is assumed to be a cap on the computational budget
allocated to data-driven verification.
Assumption 3 The simulation-based data-driven verification procedure is constrained by a computational budget and this budget manifests as a bound on the number of simulations, Ntotal , allocated to the
verification procedure.
If the simulation model used to generate trajectories Φ is of high fidelity and/or complexity, then the
computational overhead required to compute each trajectory will be non-negligible. Thus, it is unlikely a
large number of simulations can be performed within some desired time or budget. Even when the simulations
are relatively cheap, there will still be some feasible upper bound and only a finite number of measurements
can be obtained. Due to the variety of possible causes and constraints, Ntotal is merely intended as a
simple constraint to capture the underlying limit on the number of simulations allocated to the verification
technique, irrespective of the root cause.

III.

Prior Work: Analytical Barrier Certificates

The following section will overview prior work9 in analytical barrier certificates that are the most relevant
analytical verification techniques to the problems of interest. This discussion will serve to not only highlight

the strengths of these techniques, but also their limitations. Two of the examples in Section VI will compare
analytical certificates against the results from the paper’s data-driven verification procedures.
A.

Approach

Analytical barrier certificate techniques are a common solution for verification of nonlinear systems. In
particular, these approaches are applied to systems with stability or boundedness requirements, where the
state x(t) must remain outside of a set of failure states Xf ail for the entire trajectory length. While stability
and boundedness covers a wide range of possible problems, it is important to note that this does not include
some of the more complex spatial-temporal requirements discussed in Section A. Thus, these analytical
techniques do not apply to all possible systems of interest. While this does limit their application, analytical
barrier certificates are extremely useful and popular techniques whenever appropriate.
At their core, Lyapunov function-based barrier certificates4, 6, 7 rely upon a continuously-differentiable
scalar function V (x) : Rn → R called a Lyapunov function to analyze the performance of the system and
bound the trajectory response. These Lyapunov functions construct invariant sublevel sets that ensure the
trajectory never enters Xf ail . Given a scalar η > 0 and a Lyapunov function V (x), an η-sublevel set ΩV,η
defines the set of all x ∈ X for which V (x) is bounded below η: ΩV,η := {x ∈ Rn | V (x) ≤ η}. Furthermore,
if V (x(0)) ≤ η and V (x(t)) ≤ η ∀t > 0, then this sublevel set is called an invariant sublevel set. If it can
be proven that V (x) > η ∀x ∈ Xf ail , then any trajectory which starts within the invariant sublevel set ΩV,η
will never leave and the trajectory is guaranteed to satisfy the requirement.
This invariant sublevel set can then be used to form a barrier certificate to estimate Θsaf e . Given a known
V (x) and β that define an invariant sublevel set ΩV,β , this invariant sublevel identifies a set of θ values,
b V,β , which are proven to be safe, i.e. Θ
b V,β ⊂ Θsaf e . The barrier certificate essentially separates
called Θ
b V,β and the remaining θ values which cannot be proven safe with ΩV,β .
Θ into two classes: those θ ∈ Θ
b V,β are classified as elements of Θf ail . An example of
As they cannot be proven safe, the remaining θ 6∈ Θ
Lyapunov function-based barrier certificates applied to the well-studied unstable Van der Pol oscillator25 is
shown in Figure 1. In this figure, a 6th order Lyapunov function defines an invariant set (solid black line)
which accurately approximates the true separation boundary between Θsaf e and Θf ail (red dotted line). θ
b V,β are proven to be stable (the performance requirement).
vectors on or within this Θ
B.

Limitations

While the strength of a theoretically-proven certificate is obvious, these techniques suffer from a variety of
limitations which restrict their applicability and help motivate the work in data-driven statistical verification.
As mentioned earlier, not all verification problems can be addressed with analytical certificates. Even if they
do apply, a number of issues may arise which restrict the utility of analytical certificates.
Lyapunov function-based certificates are conservative by design and can result in vastly overly-conservative
b V,β . The underlying source of the conpredictions with poor construction of the invariant sublevel set θ ∈ Θ
b V,β are not necessarily elements of Θf ail ,
servativeness with barrier certificates is that all the remaining θ 6∈ Θ
but because nothing can be proven about these θ, they are conservatively classified as elements of Θf ail .
One of the ways to reduce the conservativeness of the barrier certificate is to provide the largest bound β for
which ΩV,β is still proven safe. In most cases, this optimal bound is unknown a priori, but simulation-guided
techniques6–8 can be used to identify the optimal or near-optimal values. These procedures take the Lyapunov function V (x) and generate simulation traces to identify invariant sublevel sets in order to maximize
β. This simulation-based method was used to construct the certificates found in Figure 1.
Another source of the conservativeness is the choice of Lyapunov function itself. If the order or structure
b V,β . This is
of V (x) is poorly chosen, then even the corresponding optimal β will result in a conservative Θ
seen in the 2nd and 4th order function-based certificates in Figure 1. These two certificates were constructed
with the same simulation-guided method, but were restricted to simpler forms of the Lyapunov function.
Even the simulation-guided β-maximization technique is handicapped by the incorrect order of the Lyapunov
function and cannot accurately predict the true separation boundary. These results highlight the importance
of proper selection of the Lyapunov function as the certificate is highly dependent upon this term.
The problem is that in many applications, the “correct” V (x) is not immediately obvious or even available. For instance, if the correct order of a polynomial Lyapunov function is not known, then it may
be desirable to start with lower-order functions as the computational complexity is lower and then run

Figure 1: 2nd, 4th, and 6th order Lyapunov function-based barrier certificates25 for the region-of-attraction of an
unstable Van der Pol oscillator. In this problem, the perturbations θ are simply the initial conditions of the second
order dynamics, i.e. θ = [x1 (0), x2 (0)]T . The estimated separation boundaries predicted by the barrier certificates
are shown as the black lines, while the true boundary is given by the dotted red line. The lower order functionbased certificates each produce a set of provably-safe θ values, but these sets are conservative and fail to accurately
approximate the true boundary.

simulation-guided techniques.7 Likewise, higher complexity Lyapunov functions might not be available even
with the current state-of-the-art. For example, model reference adaptive control (MRAC) procedures mostly
rely upon quadratic Lyapunov functions to prove stability and convergence.26, 27 Recent techniques using
these quadratic Lyapunov functions alongside simulation-guided barrier certificate construction9 were able
to obtain a barrier certificate, but this certificate was not able to determine the true Θsaf e , as seen later
in Figure 8. During these situations, it is desirable to apply a supplementary approach not limited by the
Lyapunov function V (x) to explore Θ and provide feedback as to whether a different Lyapunov function
should be used.
Although the theoretical guarantees provided by analytical certificates are extremely powerful, the issues
discussed in this subsection highlight the limitations of even state-of-the-art certificate construction techniques. The new work presented in the following sections is intended to supplement the existing simulationguided certificate construction methods or replace them when they aren’t applicable.

IV.

Data-Driven Verification

In place of barrier certificates, data-driven verification will generate statistical classifiers to provide estimates of Θsaf e and Θf ail . The key insight for data-driven verification is that the verification problem can
be re-envisioned as a data classification problem. In particular, since there are only two labels (“safe” or
“unsafe”) corresponding to satisfaction of the requirements, this is a binary classification problem.15 The
resulting statistical classifier must identify whether a particular θ vector belongs to Θsaf e or Θf ail given a
model constructed from a finite training set of observed simulation trajectories. The following section provides background and describes the initial construction of the statistical certificate during the data-driven
verification process.

A.

Initial Training Dataset

Before anything can be done, the statistical classification model first requires an initial training dataset of
simulation trajectories and their binary measurements of requirement satisfaction y(θ) ∈ {−1, 1}. Various
techniques can be used to generate the initial training points. The most straightforward approach is to
passively select training locations at which the simulations will be performed. Passive sampling refers to
design of experiments (DOE)16, 17 or Monte Carlo sampling techniques28 which do not require any prior
information and do not explicitly consider the current statistical model when selecting upcoming training
locations. This black-box, model-agnostic approach makes them appropriate for any type of problem, particularly when there are no available analytical certificates. The problem with passive sampling procedures
is that they may “waste” samples in uninformative or well-modeled regions of Θd since samples are taken in
all areas of Θd with equal weight or probability.
When data-driven verification augments an existing analytical barrier certificate, the data-driven procedure can exploit this prior information to improve the expected informativeness of the passive training
dataset. If an analytical certificate was already obtained using simulation-guided construction techniques,25
then the simulation data produced from this method can be used as the statistical classifier’s training dataset.
The “SimLFG” sampling algorithm25 used by these analytical techniques randomly selects points along potential sublevel sets of V (x) as it iteratively searches for an appropriate β. The end result is a distribution
of sample points, most of which are clustered in close proximity to the boundary of the returned barrier
certificate. In cases where an appropriate V (x) and β are given, but were not constructed using those
simulation-guided techniques, the same SimLFG algorithm can still be used by itself to generate the training
dataset. The benefit of generating the samples in this manner is that training samples are concentrated
mostly on or outside the current theoretically-proven boundary of the analytical certificate. This helps avoid
the placement of samples within a region of θ values already known to satisfy the requirements and thus
“wasting” them in uninformative areas. However, if this information is not available, then the samples must
be generated using design of experiments or Monte Carlo sampling.
B.

Support Vector Machines As Statistical Certificates

The data-driven procedure will construct a statistical model to predict the satisfaction of the requirements
at queries of θ from the initial training dataset, labeled L, consisting of training locations ΘL and corresponding measurements y. Given the binary nature of the problem (y(θ) ∈ {−1, 1}), there are a number
of appropriate machine learning/data-mining techniques such as relevance vector machines,29 kernel logistic
regression,30 decision trees,31 and random forests,30 but the most common and heavily-utilized binary classification approach is support vector machines (SVMs).32–34 Support vector machines form sets of hyperplanes
that separate Θ into estimates of Θsaf e and Θf ail when applied to this verification problem. Additionally,
nonlinear support vector machines can efficiently handle arbitrary datasets without modification, even when
the training data L is not linearly separable. In these cases, the “kernel trick”33 projects the data into a
higher-dimensional representation where the dataset is linearly separable. Initial work in data-driven verification15 has already demonstrated the ability of SVM-based classifiers to accurately separate Θ in a variety
of dynamical systems with linearly- and nonlinearly-separable datasets.
Since it will be unknown in advance if the data is linearly separable, nonlinear SVMs should be used as
the default as they will work in both cases. The nonlinear SVM used for binary classification is given by
H(θ) ∈ R,
Nsv
X
H(θ) =
αj yj κ(θ j , θ) + b,
(2)
j=1

where the sign of H(θ) returns a predicted label yb(θ) for a given query point θ,

yb(θ) = sign H(θ) .

(3)

During the SVM training process, the procedure selects Nsv support vectors from the training set L and
forms the separating hyperplanes used to construct H(θ). These Nsv ≤ |L| support vectors are specifically
chosen from the training set to produce the optimal classifier given the current observations. Each of these
support vectors θ j has an associated non-zero weight αj , while non-support vector elements of L can be
viewed as having weights α = 0. The SVM classifier places a nonlinear kernel function κ(θ j , θ) at each
support vector location to project the dataset from Rp to a higher-dimensional space. Different types of

kernel functions are possible, but this work uses the common isotropic radial basis function (RBF) kernel.
This isotropic RBF kernel function is given by
 −||θ − θ||2 
j
,
κ(θ j , θ) = exp
γ2

(4)

where θ j is the location of the indicated support vector and scalar hyperparameter γ is the kernel width.
In order to compute the optimal set of support vectors from L and their associated weights α, the
support vector machine training algorithm solves the following optimization program. Equation 5 gives the
Lagrangian dual of the primal quadratic program for soft-margin nonlinear SVMs,32
maximize
α

s. t.

Nsv
X
j=1
Nsv
X

N

αj −

N

sv X
sv
1X
αi αj yi yj κ(θ j , θ i )
2 j=1 i=1

(5)

αj yj = 0

(6)

0 ≤ αj ≤ Cα .

(7)

j=1

Rather than compute the solution using the primal form, it is generally easier to solve the problem in this
dual form. Before the training process can begin, two hyperparameters must be set by the programmer: the
α box constraint Cα from (7) and the kernel width γ from (4). The choice of hyperparameters will affect the
selection and number of support vectors. In this paper, the kernel width is kept fixed to γ = 1, although it
can also be estimated at part of a hyperparameter optimization step within the training process. Equation
2 also includes an optional bias term b, but this bias is set to b = 0 hereafter. The end result of this training
procedure is the statistical classifier shown in (2). The classifier and the resulting predictions yb(θ) for all
b saf e and Θ
b f ail .
θΘd will then define two sets, Θ
C.

Prediction Error

As with analytical certificates, data-driven certificates are susceptible to misclassifications when labeling
unseen locations in Θ. Like before in Section III, some points that were labeled as “unsafe” by the classifier,
b f ail , may actually belong to Θsaf e ; however, data-driven classifiers also introduce the dangerous
θ ∈ Θ
possibility that unobserved unsafe locations θ ∈ Θf ail will be incorrectly labeled as “safe” and included in
b saf e . While the possibility of these errors cannot be completely eliminated, a number of steps can be taken
Θ
during training to estimate their likelihood and minimize their impact.
First, both types of misclassification errors have unequal consequences. In general, it is much more
acceptable to accidentally label θ ∈ Θsaf e as “unsafe” than the reverse. The scalar constraint term Cα from
the training process can be recast as a matrix to reflect the unequal cost,
"
#
CF N
0
Cα =
.
(8)
0
CF P
Higher CF P values will place higher training penalties on false positives, that is incorrect labeling of unsafe
points as “safe.” Meanwhile, lower CF N will loosen the penalties on the false negatives, when safe points are
incorrectly labeled as “unsafe.” In this manner, the optimization program will adjust the support vectors
and their weights to avoid false positives, potentially at the cost of more false negatives.
Additionally, while the Cα constraint will help control the number of misclassification errors, it is impossible to completely eliminate the possibility they exist. Therefore, it is advantageous to estimate the likelihood
of misclassification errors, particularly in unobserved spaces of Θ. The primary method to estimate the likelihood of misclassification erros is to compute the numerical generalization error on an independent validation
dataset. This independent validation dataset contains observed datapoints and their corresponding measurements where the accuracy of the classifier’s predictions can be tested. The most common approaches are
leave-one-out and k-fold cross-validation approaches which segment the training data L into a set of training
points actually used to construct H(θ) and an independent testing set used to validate the performance of
H(θ)34 on “held-out” data.

These two validation approaches can be used to estimate the likelihood of misclassification error, but
they are limited in some aspects. First, the accuracy of the cross-validation error with respect to the true
generalization error over Θ is limited by the quality of L. If the observed training dataset does not adequately
cover a region in Θ, then the independent testing set will fail to accurately incorporate its effect upon the
generalization error. Likewise, the creation of independent, held-out datasets from the training set, reduces
the amount of already-limited training data that can be used to construct the classifier. Simply put, all the
data held out from L to create the independent validation set, could have been used to produce a better
classifier.
While their presence is noted, these limitations are an unfortunate effect of the reliance upon binary
measurements. As mentioned in Section II-A, parallel work24 has developed a similar data-driven approach
when continuous-valued measurements of requirement satisfaction are available. In that work, prediction
confidence is computed online without having to segment L due to the extra layer of information provided
by continuous measurements. While this addresses most of the misclassification issues discussed above,
continuous-valued signal temporal logic (STL) measurements are not always available. For systems with
binary MTL measurements, STL-like prediction confidence can be estimated without a separate validation
set using Platt scaling.35 This method fits a logistic regression model to the SVM classifier and produces
probabilistic estimates of prediction confidence much like systems with continuous STL measurements.24
Unfortunately, as with kernel logistic regression, Platt scaling assumes linearly separable data, which limits
its accuracy when applied to nonlinearly-separable datasets. Regardless of their respective drawbacks, many
methods exist for the computation of prediction error in the data-driven verification process.

V.

Closed-Loop Verification

Section IV described a method for constructing the statistical certificate given an initial set of training data. As with all statistical learning approaches, data-driven verification is fundamentally limited by
the training dataset used to construct the certificate. If the data is limited to only a small region of Θ,
then the statistical certificate’s accuracy will degrade past this small region. Previous work in data-driven
verification15 successfully sidestepped this issue by generating very large training datasets of randomlyselected points that implicitly covered Θ with adequate resolution. These large training sets pose a problem
if simulation trajectories are computationally expensive to obtain or the number of simulations budgeted
for the verification procedure is small. In an ideal scenario, simulations would only be performed at the
most-informative θ locations distributed along the true boundary (or boundaries) that separates sets Θsaf e
and Θf ail . This training set would minimize the predictive error without relying upon a large number
of simulation samples; however, such a training dataset would require the true boundary to be known in
advance.
The work in this section exploits active learning to cluster samples near the estimated boundary and
approximate the ideal training set without any prior knowledge of the true boundary. The main contribution
of this work is a novel verification process that iteratively selects potentially-informative training points in
order to maximize the accuracy of the predictions. Due to the iterative nature of the data-driven verification
process, we call the approach closed-loop verification. The following subsections will present methods to
identify informative locations for future simulations as well as detail two closed-loop verification procedures.
A.

Selection Criteria

Active learning describes a closed-loop sampling procedure that iteratively trains a statistical model and
selects new sample points to improve the current objective.18–22 Rather than passively select training locations, either in a grid or random fashion, active learning uses the current model information to discriminate
between potential sample locations and select new training locations which are expected to best improve
the current model. Depending on the objective, the “best” possible next sample points might be different
and there are a number of different sample selection criteria corresponding to various objectives. A comprehensive overview of some of these different objectives and selection criteria is found in Settles’ book.18 At
its core, active learning is a closed-loop process that exploits the most recent information, regardless of the
particular selection criteria. This allows the sampling procedure to minimize the expected prediction error
of the classifier given a designed number of allowable simulations, and thus avoiding the reliance upon large
datasets used in the earlier open-loop data-driven verification procedures.15

1.

Expected Model Change

Due to the desire to minimize the rate of misclassification errors, this work uses the expected model change
metric21 to rank the prospective sample locations and identify which location is expected to most improve
the current data-driven verification model. The objective is to select the sample point that, if measured,
would induce the largest expected change upon the current model once the model is retrained with this new
information. For SVMs, it was shown21 that the point which maximizes expected model change can be found
through the gradient of the Lagrangian dual objective function in (5). If a new sample location, labeled θ +1 ,
and its hypothetical measurement y(θ +1 ) are added to the current SVM, the resulting dual objective would
be given by
Nsv +1 NX
NX
sv +1
sv +1
1 X
αi αj yi yj κ(θ j , θ i ) .
(9)
D(α) =
αj −
2 j=1 i=1
j=1
With the weighting term α+1 initialized to zero, the gradient of (9) with respect to α+1 is then
Nsv
X
∂D(α)
= 1 − y(θ +1 )
αj yj κ(θ j , θ +1 ) .
∂α+1
j=1

(10)

Assuming the bias b from (2) is set to 0, this gradient is simply
∂D(α)
= 1 − y(θ +1 )H(θ +1 ) ,
∂α+1

(11)

where H(θ +1 ) is the output of the current SVM model at location θ +1 . Since α+1 must be non-negative,
the model should only be updated if (11) is positive, meaning y(θ +1 )H(θ +1 ) < 1. Note that in order for
y(θ +1 )H(θ +1 ) < 0, there must be disagreement between the sign of current prediction H(θ +1 ) and the
actual measurement y(θ +1 ). If this actual y(θ + 1) was indeed known, then the optimal sample θ which
maximizes the model change would induce the largest gradient in (11),


θ = argmax 1 − y(θ +1 )H(θ +1 ) .
(12)
θ +1 ∈Θ

The problem with the selection metric in (12) is that y(θ +1 ) is not known until after a simulation has already
been performed at that location. In place of the actual model change, the expected model change is used as
the selection metric. This expected model change is given by


(13)
θ = argmax 1 − yb(θ +1 )H(θ +1 ) ,
θ +1 ∈Θ

with yb(θ +1 ) = E[y(θ +1 )] as the expected measurement at θ +1 . Shown earlier in (3), the expected measurement yb(θ +1 ) is actually the sign of the output of the current classifier H(θ +1 ), forcing yb(θ +1 )H(θ +1 ) > 0 to
always hold true. Therefore, the “best” sample point according to the expected model change (13) is simply
the point whose classifier output H(θ) is closest to zero, i.e. |H(θ)| ≈ 0, since yb(θ +1 )H(θ +1 ) will never be
negative. This premise forms the basis of the active sampling procedure used for closed-loop verification.
B.

Sequential Sampling

Algorithm 1 details the sequential closed-loop verification procedure. The procedure selects samples based
upon the aforementioned expected model change metric, rewritten as
θ = argmin |H(θ)|.

(14)

θ∈U

The procedure assumes samples are taken from a set of unobserved samples locations, set U, such that
ΘL ∩ U = ∅. Given the existence of the sufficiently-fine Θd and an initial training set L taken from Θd , the
unobserved set U is simply the remainder U = Θd \ ΘL . Additionally, it is assumed the discretization Θd is
so fine that the final size of U is still significantly larger than Ntotal , i.e. |U|  Ntotal .
The closed-loop process is as follows. Given the initial training set L and classification model H(θ), the
procedure in Algorithm 1 selects the remaining T = Ntotal − |L| samples from U. Once the next training
location θ has been selected (line 3), a simulation is performed using that θ vector and the binary MTL
measurement is obtained (line 4). The new measurement is added to the training set L (line 5) and the
model is retrained with this new dataset (line 6). This process repeats until the size of L reaches Ntotal .

Algorithm 1 Sequential active learning using the expected model change selection metric
1:
2:
3:
4:
5:
6:
7:

C.

Input: training set L, available sample locations U, initial classifier H(θ), # of additional samples T
for i=1:T do
Select θ from U according to (14)
Run simulation at θ, obtain measurement y(θ)
Add {θ, y(θ)} to training set L, remove θ from U
Retrain H(θ) with new L
end for
Batch Sampling

The procedure described in Algorithm 1 is a sequential process - it iteratively updates the model one point
at a time. While this process will correctly guide the selection of datapoints as intended, it ignores two
important considerations. First, sequential sampling procedures do nothing to address the nontrivial cost of
retraining the classifier. As the number of additional training points in L grows, so does the computational
effort required to recompute classifier model. If the model is retrained after each the arrival of every single
measurement, then the computational cost associated just with retraining in high. Additionally, sequential
sampling algorithms do not exploit the parallelism inherent in many simulation environments. In particular,
multiple simulations can often be performed in parallel due to multiple processor cores or computers. Batch
active learning processes21 address both these considerations by selecting multiple θ values before retraining
H(θ). If the number of remaining samples T is broken into TB batches of M points, the computational cost
of retraining decreases from O(N 2 T ) to O(N 2 TB /M ).
While the selection of a batch of points between retraining steps will reduce the cost associated with
retraining the SVM, it introduces another concern. If the same criteria from (14) is blindly implemented
without care, the fine resolution of points in U can result in the selection of multiple neighboring points
in close proximity as they will have similar values according to (14). In order to prevent the selection of
redundant points, a diversity measure must be incorporated into the selection criteria. This paper uses the
angle between the hyperplanes φ(θ) induced by the samples as the diversity metric. The angle associated
with a kernel κ is shown in (15). From a diversity perspective, the goal is to maximize the angle between
samples.
|κ(θ i , θ j )|
| cos(∠(φ(θ i ), φ(θ j )))| = p
(15)
κ(θ i , θ i )κ(θ j , θ j )
The batch selection criteria metric modifies the sequential sample selection criteria in (14) with the diversity
measure from (15). This new selection metric becomes

θ = argmin λ|H(θ)| + (1 − λ)max | cos(∠(φ(θ), φ(θ j )))| ,
(16)
θ∈U

θ j ∈S

which is a weighted combination of the expected model change and diversity in the batch. Here, the batch
sampling metric selects the M points in the batch sequentially and set S ∈ Θ refers to the points previously
selected for the current batch. The scalar term λ ∈ [0, 1) controls the weighting between the metrics (14)
and (15). As λ grows, the diversity metric has less effect upon θ and the selection criteria converges towards
(14). This work uses λ = 0.7 as it has demonstrated good empirical results. This balance of utility and
diversity ensures that each subsequent sample chosen by (16) is adequately different from all the preceding
points in S.
The complete batch active learning algorithm is given in Algorithm 2. On line 1, the procedure starts
with the initial training set L and classifier H(θ) along with the set of available sample points U. The
algorithm the iteratively selects the “best” data point according to (16) (line 5) and adds this θ to set S
(line 6). Once S is full, the simulations are performed at these θ locations (line 9). On line 10, set S
and the corresponding labels are added to the training set L and this new information is used to retrain
H(θ). The closed-loop verification process repeats lines 3-11 until a designated number of iterations has
been reached. As will be illustrated in Section VI, Algorithm 2 successfully distributes the sample locations
across regions of high utility (large expected model change) and implicitly improves the prediction error over
passive sampling procedures.
Although Algorithm 2 requires more intermediate steps between retraining the SVM, the net computation
cost is lower than in Algorithm 1. As mentioned in the preceding paragraphs, the training step in Algorithm

Algorithm 2 Batch active learning using the expected model change selection metric
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:

Input: unlabeled set U, training set L, batch size M , classifier H(θ), empty set S
for each iteration do
Reset: d = 1
while d ≤ M do
Select θ from U according to (16)
Add θ to S, remove θ from U
d=d+1
end while
Run simulations ∀θ ∈ S, label the results
Add set S to training set L, reinitialize empty S
Retrain H(θ) with new L
end for

1 (line 6) and Algorithm 2 (line 11) requires O(N 2 ) operations for each iteration. In Algorithm 1, the
prediction step (line 3) requires O(N |U|) operations to weight the prospective sample locations. Note that
this does not include the cost of actually performing simulations, as that will vary from example to example.
For comparison, the batch sample selection step in line 5 of Algorithm 2 requires an additional O(M |U|)
operations to calculate the diversity measure for all M points in the batch. Although this requires additional
operations, the total computational cost for the batch procedure is lower for an equal number of samples.
Figure 2 demonstrates the computational complexity for the two procedures when applied to Example
VI-B. Algorithm 2 produces a lower complexity than the sequential procedure and this improvement in
computational complexity increases with larger batch sizes M .

Computational Complexity

10

8

10 8
Algorithm 1
Algorithm 2 (M=2)
Algorithm 2 (M=5)
Algorithm 2 (M=10)

6

4

2

0
50

100

150

200

250

# Samples
Figure 2: Computational complexity of the sequential (Algorithm 1) and batch (Algorithm 2) closed-loop verification
procedures when applied to Example VI-B. By reducing the number of retraining steps for the sample number of
total samples, Algorithm 2 lowers the computational cost required for closed-loop verification. The exact reduction
varies according the particular example, but a larger batch size M will lower the complexity.

VI.

Results

The closed-loop data-driven verification procedure is demonstrated on a number of case studies. The
results for each example are summarized, but the first two case studies also have a link to accompanying videos
and supplementary material. These videos best illustrate the iterative nature of the approach. Additionally,
the first example highlights that data-driven verification can be applied to a wider array of problems, not
just adaptive control problems. The third example also demonstrates data-driven analysis on a verification

(a) Initial L and the barrier certificate produced with a (b) Starting SVM model trained using the initial training set
quadratic Lyapunov function

Figure 3: Initial training set L and data-driven SVM certificate compared against the analytical certificate produced
by the 2nd order Lyapunov function.

problem which cannot be addressed with current analytical barrier certificate techniques.
A.

Unstable Van der Pol Oscillator

The first example is not an adaptive control system, but rather the well-studied Van der Pol oscillator problem
shown previously in Figure 1. This example not only demonstrates the use of data-driven verification
to supplement existing analytical barrier certificates, but also highlights its application to a wider class
of systems than adaptive systems. The nonlinear dynamics in (17) have an unstable limit cycle and an
asymptotically stable equilibrium point at the origin,
ẋ1 = −x2
ẋ2 = x1 +

(17)
(x22

− 1)x2 .

The verification goal is to identify initial conditions θ = [x1 (0) x2 (0)]T that will cause the system to either
converge to the origin or diverge. The region of θ values for which the trajectory converges, Θsaf e , is also
known as the region-of-attraction (ROA).
This problem has been well-studied in analytical verification work and suitable V (x) and β terms are
readily available.7, 25 In order to demonstrate the ability of the data-driven procedure to rapidly improve upon
an analytical certificate produced by a poorly-chosen Lyapunov function, a quadratic Lyapunov function is
used as the starting condition. It was shown in Figure 1 that a 6th order Lyapunov function is required to
accurately approximate the true boundary of the ROA; therefore, the analytical certificate produced by a
quadratic Lyapunov function will fail to approximate the true Θsaf e . The discretized set Θd consists of 16,000
points covering θ1 : [−3, 3] and θ2 : [−3, 3]. Even with the maximum level set possible with the quadratic
V (x), the analytical certificate incorrectly labels nearly 12% of Θ as unsafe because these θ locations fall
outside the level set defined by the quadratic function. This is shown in Figure 3(a).
This example considers the batch closed-loop verification procedure from Algorithm 2. The initial SVM
model is formed from 50 training samples obtained from the SimLFG sampling algorithm used to produce
the simulation-guided analytical barrier certificate25 seen in Figure 3(a). This initial SVM model is pictured
in Figure 3(b). Once this initial model has been constructed, 20 iterations of Algorithm 2 are performed in
batches of M = 10 points. The selection of the first batch of 10 points is illustrated in Figure 4. The influence
of the diversity measure is clearly visible at the conclusion of the process in Figure 4(d). The samples are
spread roughly evenly across the areas of high expected model change in order to prevent redundancy in the
batch.
Even after only 5 iterations (100 training points in total), the statistical classifier converged to a close
approximation of the true shape of the ROA, seen in Figure 5. The active learning procedure is also

(a) Selection of the 1st sample in the batch of 10 points

(b) Selection of the 2nd sample

(c) Selection of the 3rd sample

(d) Completed batch of 10 samples

Figure 4: Selection of the first batch of 10 points according to the expected model change criteria in (16).

compared against a passive procedure that randomly selects batches of 10 samples (without replacement)
from U. Figure 6 displays the mean and 1-σ distribution of the true misclassification error according to 50
repetitions of both of these procedures. While both procedures already outperform the 12% misclassification
error of the analytical certificate at the very start before obtaining additional samples, active learning
ultimately has both a lower mean classification error and 1-σ distribution than passive data-driven verification
in subsequent iterations. A video of the active learning procedure can be viewed onlilne here: https:
//youtu.be/aXTi99edLoA.
Although the results from Figure 6 quantify the improvement active learning offers over passive approaches, this true misclassification error is unavailable during the actual process as it requires the full set
Θd to be known in advance. Instead, the misclassification error can be approximated using the k-fold cross
validation and independent validation set methods described in Section IV-C. The estimation errors produced
by these two methods are shown in Figure 7. While k-fold cross validation is more desirable from a computation standpoint as it operates directly on the current training set, the results in Figure 7(a) indicate that it
is not suitable for use with active learning and closed-loop verification. The estimated misclassification error
increases as the number of samples grows because the closed-loop verification algorithm clusters samples
along the approximate boundary. This clustering will result in the misclassification of a higher proportion
of samples than if the data is randomly spread across Θd . In order to obtain an accurate approximation of
the misclassification error, an independent validation set must be obtained in parallel to L.

Figure 5: Statistical model of Θsaf e /Θf ail for the Van der Pol example after only 5 iterations of Algorithm 2.

True Error (%)

20

Barrier Certificate
Passive DOE
Closed-loop w/ EMC

15
10
5
0
50

100

150

200

250

# Samples
Figure 6: Comparison of the true misclassification error using both active and passive, randomized sampling for the
Van der Pol example. Note that both procedures outperform the analytical certificate, which has a misclassification
error of nearly 12%.

15
Passive DOE
Closed-loop w/ EMC

10

30

20
Passive DOE
Closed-loop w/ EMC

10

0

Estimation Error (%)

Estimation Error (%)

40

5
0
-5
-10

-10

-15
50

100

150

200

250

50

100

# Samples

150

200

250

# Samples

(a) K-fold cross-validation

(b) Independent validation set

Figure 7: Comparison of the estimation error produced by k-fold cross-validation and a randomized, independent
validation set.

B.

Concurrent Learning Adaptive Control System

The second case study is a model reference adaptive control system. In this problem, a concurrent learning
model reference adaptive controller27 (CL-MRAC) is used to control an uncertain second order system.
ẋ1 = x2 ,

(18)

ẋ2 = (−0.2 + θ1 )x1 + (−0.2 + θ2 )x2 + u.
Unlike the previous example, the perturbations [θ1 , θ2 ] are uncertain parameters in the system. The adaptive
controller estimates these parameters while simultaneously controlling the states to track a desired reference
trajectory produced by the linear system in (19)
ẋm1 = xm2
ẋm2 =

−ωn2 xm1

(19)
− 2ζn ωn xm2 +

ωn2 zcmd (t)

,

with ζn = 0.5 and ωn = 1. This reference model is excited with reference commands zcmd = 1 between 0
and 2 seconds, zcmd = 1.5 between 10 and 12 seconds, and zcmd = −1.5 between 20 and 22 seconds. At all
other times in the 40 second trajectory length (Tf = 40), the reference command defaults to 0. Note that
this reference signal is not persistently exciting.
The control input u(t) consists of three components: the reference input urm , the feedback input upd ,
and the adaptive input uad .
u(t) = urm (t) + upd (t) − uad (t)
(20)
In the absence of uncertainties, the reference input and a linear feedback controller are suitable for ensuring
closed-loop tracking.
urm = −ωn2 xm1 − 2ζn ωn xm2 + ωn2 zcmd (t)

(21)

upd = Kp e1 (t) + Kd e2 (t)

(22)

The tracking error e(t) = xm (t) − x(t) is the difference between the desired reference trajectory and the
actual state trajectory. In this example, the feedback controller gains are set to Kp = 1.5 and Kd = 1.3.
Due to the presence of the uncertainties, the third uad (t) component is necessary to ensure the closed-loop
system correctly tracks the desired reference trajectory.
"
#
h
i
x1 (t)
b
b
uad (t) = θ1 (t) θ2 (t)
(23)
x2 (t)

The terms θbT = [θb1 (t) θb2 (t)] are the estimated parameters that are updated online according to the
concurrent learning adaptive law27
pX
max
˙
θb = −Γx(t)e(t)T P B − Γc
xk xTk θe

(24)

k=1

with B T = [0 1], Γ = 2, Γc = 0.2, and pmax = 20. The symmetric positive-definite matrix P is obtained
from the Lyapunov equation AT P + P A = −I where A is the nominal open-loop plant (with θ = 0). This
e = θ(t)
b − θ will
adaptive law ensures that both the tracking error e(t) and the parameter estimation error θ(t)
asymptotically converge to 0 without requiring zcmd (t) to be persistently exciting. The cornerstone of this
convergence is the history stack in the second half of (24) and its periodic update according to the singular
27
value maximizing
The stored datapoints xk are specifically chosen from the trajectory to
Ppmax algorithm.
T
guarantee k=1 xk xk > 0 and are updated only if the new datapoint will improve the rate of convergence.
While this adaptive control law is used to improve the closed-loop response of the system, it transforms
the second-order linear plant into a nonlinear system. The history stack in particular greatly complicates
analysis of the closed-loop response due its periodic, but non-uniform updates of the saved datapoints.
Although the CL-MRAC controller ensures asymptotic stability of the closed-loop system, safety of the
system in this example is measured by whether the state x1 (t) remains within unit error of the reference
xm1 (t). Described in metric temporal logic format, this performance requirement states
ϕbound = [0,Tf ] (1 − |e1 [t]| ≥ 0) .

(25)

Even though the system is guaranteed to the stable with respect to the tracking error e(t), this stability does
not directly guarantee e1 (t) will remain within the specified bounds. The goal of the verification process
is to identify which θ vectors will not cause the trajectory to exceed those bounds and which will. Similar
to before, an analytical barrier certificate was previously obtained in prior work9 and is based upon the
quadratic Lyapunov function commonly used for adaptive systems. Unlike the first example, there are no
alternative Lyapunov functions available; a general higher-order Lyapunov function for CL-MRAC systems
has not been found.
Twenty iterations of the closed-loop data-driven procedure are performed in batches of 10 points. In
this problem, the unlabeled set U consists of 36,400 points covering θ1 : [−8, 8] and θ2 : [−10, 10]. The
initial training set L consists of 50 randomly-selected points and is used to form the initial SVM, pictured
in Figure 8. The conservativeness of the 2nd order Lyapunov function-based barrier certificate is readily
apparent in that figure. Figure 9 overlays the selection criteria from (16) over Θ and displays the 10 points
that make up the first batch. As before, the selected sample locations are spread across regions of high
expected model change. Figure 10 shows the data-driven classifier after 20 iterations. While the nonconvex
shape of the Θsaf e /Θf ail boundary is significantly more challenging than the shape of the Van der Pol
oscillator’s limit cycle, the statistical model still converges toward the true boundary and minimizes the
prediction error. The mean and 1-σ distribution of the true misclassification error for 50 repetitions are
shown in Figure 11 and reinforce this observation. A video of the active learning procedure can be viewed
online here:
https://youtu.be/ID_fZ2SKpIk.
Although active learning is shown to reduce the true misclassification error over analytical certificates,
there does exist the possibility of unsafe misclassification errors - θ locations that actually fail to meet the
requirements accidentally labeled as “safe.” In the training process described in Section IV, the Cα weighting
term from (7) can be adjusted so that these unsafe errors are penalized higher than the reverse. This will
introduce conservativeness into the construction of the classifier, seen in Figure 12. While this does increase
the true misclassification error from 2.58% at Cα = 1 to 9.50% at Cα = 5, it reduces the percentage of
unsafe misclassification errors from 1.30% to 0. This also holds true when the results are compared on an
independent validation set of equal size of L - estimated total error increases from 4.44% to 9.60%, but unsafe
error decreases from 2.80% to 0. The choice of Cα is left to the designer, but Cα 6= 1 should only be used
after completion of the active sampling procedure to avoid handicapping the search for the true boundary.
C.

Adaptive Control System with Control Saturation

The third example is a variant of the same CL-MRAC system from the preceding example with new performance requirements and more complex constraints. The open-loop system plant is the same as in (18);

Figure 8: The initial training set L and SVM for the system in (18).

Figure 9: Selection criteria for the first batch of 10 samples using Algorithm 2.

Figure 10: Statistical certificate after 20 iterations of Algorithm 2.

True Error (%)

15

Passive DOE
Closed-loop w/ EMC

10
5
0
50

100

150

200

250

# Samples
Figure 11: Comparison of the true misclassification error for the system in (18) using both active and passive, randomized sampling. Both of these outperform the analytical certificate, which misclassifies 35% of Θsaf e as “unsafe”.

Figure 12: Comparison of the data-driven classifier given different penalties on unsafe misclassification errors. The
terms Cα = 1, 2, 5 refer to the ratio of training penalty placed on mislabeling unsafe locations as “safe” against
mislabeling safe locations as “unsafe”.

however, the control input u(t) is saturated with limits −umax ≤ u(t) ≤ umax . Likewise, the perturbations
θ are different in this new problem. In addition to the same parametric uncertainties [θ1 , θ2 ], there is also
uncertainty in the initial position x1 (0). Various levels of possible control saturation limits umax are also
considered during the verification process as a 4th dimension of θ, although umax would generally be known
at run-time.
The same CL-MRAC control law from (20) and (24) is used, but the presence of the control saturation
can lead to instability in the adaptation of θb if left unaddressed. In order to counter the adverse effects
of control saturation on the adaptation, pseudo-control hedging (PCH)36 is used to augment the adaptive
system. Pseudo-control hedging creates a hedge input νh from the difference between the desired control
input before saturation udes (t) and the control saturation limit umax . This desired control input before
saturation udes (t) is produced by the same equation (20) while the actual control input after saturation is
u(t).
udes (t) = urm (t) + upd (t) − uad (t)
(26)

νh =




umax − udes (t)

if udes (t) > umax

0
otherwise



−umax − udes (t) if udes (t) < −umax


if udes (t) > umax

umax

u(t) =

udes (t)



−umax

otherwise

(27)

(28)

if udes (t) < −umax

The issue that arises due to saturation is the controller can no longer perfectly track the reference trajectory
xm (t) when the control input u(t) is saturated. Instead, the pseudo-control hedge νh modifies the reference
model in order to prevent the saturation from negatively affecting the tracking error e(t) = xm (t) − x(t). In
this example, the PCH-modified reference model becomes
ẋm1 = xm2

(29)

ẋm2 = −ωn2 xm1 − 2ζn ωn xm2 + ωn2 zcmd (t) − νh .
In addition to the changes to the plant and control architecture, the closed-loop system is also expected
to satisfy a new set of temporal requirements. More specifically, the system’s trajectory must satisfy three
conditions in order for it to be labeled as “safe”. First, the trajectory must fall within the interval x1 (t) ∈
[0.7, 1.3] at some point between t = 2 and t = 3 seconds.
ϕ1 = ♦[2,3] (x1 [t] − 0.7 ≥ 0) ∧ ♦[2,3] (1.3 − x1 [t] ≥ 0)

(30)

Next, the trajectory must satisfy a similar requirement between t = 12 and t = 13 seconds, where the
trajectory should reach x1 (t) ∈ [1.1, 1.7] at some point in the time interval.
ϕ2 = ♦[12,13] (x1 [t] − 1.1 ≥ 0) ∧ ♦[2,3] (1.7 − x1 [t] ≥ 0)

(31)

Finally, the trajectory should reach some point in the interval x1 (t) ∈ [−1.6, −1.2] at t = 22.5 seconds. This
requirement is modeled as
ϕ3 = [22.4,22.6] (x1 [t] + 1.6 ≥ 0) ∧ [22.4,22.6] (−1.2 − x1 [t] ≥ 0) .

(32)

The trajectory is only labeled as “safe” if it is able to satisfy all three of the requirements.
ϕ = ϕ1 ∧ ϕ2 ∧ ϕ3

(33)

Although the nominal reference model was designed to meet all three temporal requirements, PCH adjustments to the reference model will change the resulting reference trajectory and can accidentally lead both
the true and reference trajectories to fail to meet one or all of the temporal requirements. Therefore, it
is critically important to examine whether the actual trajectory x(t) satisfies those requirements and not
simply focus on the tracking error e(t) as in the second example problem.
Just as in the previous two case studies, active data-driven verification learns a statistical certificate of
the satisfaction of performance requirements subject to a bound on the number of simulations. Unlike the
previous two examples, an analytical barrier certificate cannot be computed due to the complexities of the
closed-loop system and performance requirements. This highlights that data-driven verification can be applied to a wider class of possible certification requirements than is possible with the existing barrier certificate
techniques. The procedure constructs a 4-dimensional discretized grid Θd of 1.285 million possible sample
locations, with each dimension corresponding to θ1 : [−5, 5], θ2 : [−5, 5], x1 (0) : [−1, 1] and umax : [3, 8].
From this grid, the active sampling procedure obtains 50 initial training samples and performs 45 iterations
in batches of 10 points. At the conclusion of the process, the statistical certificate is constructed using only
500 training samples and predicts performance satisfaction over the remaining set of unobserved values. The
mean and 1-σ prediction error averaged over 100 randomly-initialized repetitions is shown in Figure 13.
Just as with the previous case studies, active learning is able to outperform passive data-driven verification
with both a lower mean and standard deviation from that mean as the number of iterations grows. At the
conclusion of the process, the active learning procedures produces an average true misclassification error of
5.17% and estimated validation error of 5.15% compared against 9.16%/9.11% for the passive procedure.
Even for this large 4D system, active data-driven verification is shown to be a suitable tool for examining
the performance of a complex adaptive system.

VII.

Conclusion

This work has presented a data-driven approach for verification of uncertain dynamical systems. The goal
of the procedure is to estimate the boundary that separates safe perturbations from unsafe ones. The core of
the data-driven verification procedure is a support vector machine (SVM) model to classify all the points in Θ
and an active sampling procedure to iteratively improve upon it. Section IV described the construction of the
statistical certificate and the importance of the training data. Meanwhile, closed-loop algorithms for active
selection of future training points were presented in Section V. These closed-loop sampling processes select
sample points that will maximize the expected change in the model in order to minimize the misclassification
error. The utility of the closed-loop data-driven verification procedures was successfully demonstrated on
multiple case studies of nonlinear and adaptive control systems. When compared against passive procedures

Passive DOE
Closed-loop w/ EMC

True Error (%)

20
15
10
5
0
100

200

300

400

500

# Samples
Figure 13: Comparison of the prediction error for the constrained adaptive control system example using both active
and passive, randomized sampling procedures. Note that error bars are only shown for even-numbered iterations for
ease of viewing. Unlike the previous two examples, an analytical barrier certificate was not found and cannot be
compared against the results.

that rely upon random sampling, active learning results in both lower a mean and distribution of the
prediction errors.
Data-driven verification is a supplementary approach for analytical verification of dynamical systems.
In many cases, it is either not obvious or possible to obtain the analytical certificate that identifies the
full boundary, resulting in overly-conservative approximations. In these cases, statistical classifiers can
be used a direct replacement of these barrier certificates or provide feedback for subsequent corrections.
Additionally, data-driven verification can be applied to a larger class of problems than possible with analytical
techniques. In some applications, an equivalent analytical certificate does not exist because there are no
suitable analytical functions to provably bound the response.
While the lack of a governing analytical function allows data-driven methods to apply to a wider range of
possible systems, it is important to recognize that it comes at a cost. Unlike analytically-verified certificates,
data-driven certificates can accidentally mislabel unsafe perturbations as “safe”; therefore, safe predictions
from a data-driven certificate are not as strong and may be incorrect. The desire to minimize the likelihood of
these misclassifications was one of the primary motivations for active learning. Ongoing research has focused
on additional sampling methods and retraining steps to combine with the active learning procedure in this
work. Ultimately, this paper presents a foundation for more advanced data-driven verification methods to
certify the safety of complex dynamical systems subject to various sources of uncertainty.

References
1 Department

of Defense MIL-HDBK-1797, “Flying Qualities of Piloted Aircraft,” .
Aviation Administration FAR Part 25, “Airworthiness Standards: Transport Categority Airplanes,” .
3 Moore, J. and Tedrake, R., “Control Synthesis and Verification for a Perching UAV using LQR-Trees,” IEEE Conference
on Decision and Control, 2012.
4 Prajna, S., Optimization-Based Methods for Nonlinear and Hybrid Systems Verification, Ph.D. thesis, California Institute
of Technology, 2005.
5 Kapinski, J., Deshmukh, J., Jin, X., Ito, H., and Butts, K., “Simulation-Guided Approaches for Verification of Automotive
Powertrain Control Systems,” American Control Conference, 2015.
6 Kapinski, J., Deshmukh, J., Sankaranarayanan, S., and Arechiga, N., “Simulation-guided Lyapunov Analysis for Hybrid
Dynamical Systems,” Hybrid Systems: Computation and Control, 2014.
7 Topcu, U., Quantitative Local Analysis of Nonlinear Systems, Ph.D. thesis, University of California, Berkeley, 2008.
8 Reist, P., Preiswerk, P. V., and Tedrake, R., “Feedback-Motion-Planning with Simulation-Based LQR-Trees,” International Journal of Robotics Research, Vol. 35, 2016, pp. 1393–1416.
9 Quindlen, J. F., Topcu, U., Chowdhary, G., and How, J. P., “Region-of-Convergence Estimation for Learning-Based
Adaptive Controllers,” American Control Conference, 2016.
2 Federal

10 Clarke, E. M. and Zuliani, P., “Statistical Model Checking for Cyber-Physical Systems,” International Symposium for
Automated Technology for Verification and Analysis, 2011.
11 Zhang, Y., Sankaranarayanan, S., and Somenzi, F., “Statistically Sound Verification and Optimization for Complex
Systems,” Automated Technology for Verification and Analysis, Sydney, Australia, November 2014, pp. 411–427.
12 Zuliani, P., Platzer, A., and Clarke, E. M., “Bayesian statistical model checking with application to simulink/stateflow
verification,” Hybrid Systems: Computation and Control, 2010.
13 Hoxha, B., Bach, H., Abbas, H., Dokhanchi, A., Kobayashi, Y., and Fainekos, G., “Towards Formal Specification
Visualization for Testing and Monitoring of Cyber-Physical Systems,” International Workshop on Design and Implementation
of Formal Tools and Systems, 2014.
14 Dang, T. and Shalev, N., “Test Coverage Estimation Using Threshold Accepting,” Automated Technology for Verification
and Analysis, Sydney, Australia, November 2014, pp. 115–128.
15 Kozarev, A., Quindlen, J. F., How, J., and Topcu, U., “Case Studies in Data-Driven Verification of Dynamical Systems,”
Hybrid Systems: Computation and Control, 2016.
16 Montgomery, D. C., Design and Analysis of Experiments, Wiley and Sons, 8th ed., 2013.
17 Lin, C. D. and Tang, B., Latin Hypercubes and Space-Filling Designs, CRC Press, 2015, pp. 593–625.
18 Settles, B., Active Learning, Morgan and Claypool, 2012.
19 Brinker, K., “Incorporating Diversity in Active Learning with Support Vector Machines,” International Conference on
Machine Learning, 2003.
20 Ali, A., Caruana, R., and Kapoor, A., “Active Learning with Model Selection,” AAAI Conference on Artificatial Intelligence, 2014.
21 Kremer, J., Pedersen, K. S., and Igel, C., “Active Learning with Support Vector Machines,” Data Mining and Knowledge
Discovery, Vol. 4, No. 4, July 2014, pp. 313–326.
22 Settles, B. and Craven, M., “An Analysis of Active Learning Strategies for Sequence Labeling Tasks,” Conference on
Empirical Methods in Natural Language Processing, 2008.
23 Maler, O. and Nickovic, D., Monitoring Temporal Properties of Continuous Signals, Springer Berlin Heidelberg, Berlin,
Heidelberg, 2004, pp. 152–166.
24 Quindlen, J. F., Topcu, U., Chowdhary, G., and How, J. P., “Active Sampling for Closed-loop Statistical Verification of
Uncertain Nonlinear Systems,” 2017, Preprint.
25 Topcu, U., Packard, A., and Seiler, P., “Local stability analysis using ssimulation and sum-of-squares programming,”
Automatica, Vol. 44, No. 10, 2008, pp. 2669–2675.
26 Lavretsky, E. and Wise, K. A., Robust and Adaptive Control, Springer, 2013.
27 Chowdhary, G. V., Concurrent Learning for Convergence in Adaptive Control Without Persistency of Excitation, Ph.D.
thesis, Georgia Institute of Technology, December 2010.
28 Rubinstein, R. Y. and Kroese, D. P., The Cross-Entropy Method: A Unified Approach to Combinatorial Optimization,
Monte-Carlo Simulation, and Machine Learning, Springer, 2004.
29 Tipping, M. E., “Sparse Bayesian Learning and the Relevance Vector Machine,” Journal of Machine Learning Research,
Vol. 1, June 2001, pp. 211–244.
30 Bishop, C. M., Pattern Recognition and Machine Learning (Information Science and Statistics), Springer, 1st ed., 2007.
31 Bombara, G., Vasile, C.-I., Penedo, F., Yasuoka, H., and Belta, C., “A Decision Tree Approach to Data Classification
using Signal Temporal Logic,” International Conference on Hybrid Systems: Computation and Control, 2016.
32 Scholkopf, B. and Smola, A. J., Learning with Kernels: Support Vector Machines, Regularization, Optimization, and
Beyond, MIT Press, 2002.
33 Vapnik, V. N., Statistical Learning Theory, Wiley-Interscience, 1998.
34 Anguita, D., Ghio, A., Ridella, S., and Sterpi, D., “K-Fold Cross Validation for Error Rate Estimate in Support Vector
Machines,” International Conference on Data Mining (DMIN), 2009.
35 Platt, J., “Probabilistic Outputs for Support Vector Machines and Comparisons to Regularized Likelihood Methods,”
Advances in Large Margin Classifiers, 1999, pp. 61–74.
36 Kannan, S. K. and Johnson, E. N., “Model Reference Adaptive Control with a Constrained Linear Reference Model,”
IEEE Conference on Decision and Control, IEEE, 2010, pp. 48–53.

