Data-Driven Approximate Abstraction for Black-Box Piecewise Affine
Systems

arXiv:1801.09289v2 [] 30 Jan 2018

Gang Chen and Zhaodan Kong
Abstract— How to effectively and reliably guarantee the
correct functioning of safety-critical cyber-physical systems in
uncertain conditions is a challenging problem. This paper
presents a data-driven algorithm to derive approximate abstractions for piecewise affine systems with unknown dynamics.
It advocates a significant shift from the current paradigm of
abstraction, which starts from a model with known dynamics.
Given a black-box system with unknown dynamics and a linear
temporal logic specification, the proposed algorithm is able to
obtain an abstraction of the system with an arbitrarily small
error and a bounded probability. The algorithm consists of
three components, system identification, system abstraction,
and active sampling. The effectiveness of the algorithm is
demonstrated by a case study with a soft robot.

I. I NTRODUCTION
The proliferation of cyber-physical systems (CPSs) brings
how to effectively and reliably guarantee their correct behaviors to the forefront of problems we as control engineers
need to address. One natural choice to attain correct functioning is to consider formal methods techniques, such as
model checking [1], [2], which have been successfully used
in the formal verification and synthesis of digital circuits
and software codes [3]. In recent years, we have seen
many efforts of extending formal methods to engineering
applications, e.g., automobiles [4], [5] and robotics [6], [7].
One crucial component of formal methods is a precise and
potentially concise mathematical model of the system under
investigation. However, in reality we rarely have full knowledge of complex CPSs during their design and even testing
phase. Thus how to attain formal guarantee for systems with
partially or fully unknown dynamics becomes a problem of
practical significance.
In this paper, we aim to address this problem in the
context of abstraction [1]. Given a system model T (which
can potentially have infinitely many states) and a formal
specification φ written, for instance, in linear temporal logic
(LTL), an abstract model of T is a simpler model T 0 ;
checking whether the simpler model T 0 satisfies φ suffices to
decide whether T satisfies φ [8], [9], [10], [11]. For systems
that can be described by discrete state models, abstraction
can be achieved by using the concepts of simulation and bisimulation [1], [2]. In control community, recently there have
been many successful efforts pertaining to the abstraction of
systems of more realistic dynamics, such as those that are
piecewise affine [12], [13], [14], [15]. All these studies, as
far as we know, assume models of known dynamics, which
Authors are with the Department of Mechanical and Aerospace Engineering, University of California, Davis. Z. Kong is the corresponding author
(email: zdkong@ucdavis.edu).

significantly impedes the application of abstraction in the
analysis and design of systems with inherent uncertainties,
e.g., those needing to interact with a variety of human users
and be deployed in a variety of environments.
One principle way of mitigating uncertainties is to utilize
machine learning techniques. Actually, the integration of formal methods and machine learning has shown great potential
in the formal specification, design, verification, and validation of CPS [16], [17], [18], [19], [20], [21]. In this paper, we
will focus on how to combine machine learning techniques,
particularly system identification and active learning, and
formal methods techniques to generate approximate abstractions for systems with black-box (unknown), piecewise affine
(PWA) dynamics. PWA models partition the state space into
a finite number of polyhedral regions and consider affine
dynamics in each region [22]. It has been shown that PWA
models can approximate nonlinear dynamics with arbitrary
accuracy [12]. Moreover, there exist efficient techniques for
the identification of PWA systems, e.g., optimization-based
methods and clustering-based methods [23].
The major contribution of this paper is that it addresses
many theoretical and algorithmic issues pertaining to the
integration of existing approximate abstraction techniques
[24] and system identification techniques [23]. Given a
system with unknown PWA dynamics, the paper shows that,
by following the algorithm prescribed in the paper, it is
possible to extract an abstract model with an arbitrarily small
error and a bounded probability (under certain mild assumptions). Even though the paper focuses on PWA systems, the
preliminary results obtained in it can potentially pave the
way for future developments for systems with more complex
dynamics.
The remainder of the paper is organized as follows.
In Section II, we provide preliminaries and notation used
throughout the paper. In Section III, we formally introduce
the abstraction problem that will be solved in the paper.
Section IV presents our data-driven approximate abstraction
algorithm, together with proofs demonstrating the effectiveness of our algorithm. Section V uses a soft robot system
as an example to showcase our proposed algorithm. We
conclude with final remarks in Section VI.
II. P RELIMINARIES AND N OTATION
A N dimensional polytope X is defined as the convex
hull of at least N + 1 affinely independent vectors in RN . A
complete partition of X is a set of open polytopes Xi , i ∈ I
(I is a finite index set) in RN such that Xi1 ∩ Xi2 = ∅
for all i1 , i2 ∈ I, i1 6= i2 and cl(X ) = ∪i∈I cl(Xi ), where

cl(Xi ) denotes the closure set of Xi . According to the Hrepresentation, each Xi , i ∈ I can be represented as Xi =
{x ∈ RN : Hi x ≺ Ki }, where ≺ denotes componentwise
inequality.
A piecewise affine (PWA) system [22] can be written as
follows:
xk+1

=
f (xk ) + e

 A1 x + b1 if x ∈ X1
..
f (x) =
.


As x + bs if x ∈ Xs

(1)

where xk is the state of the system at step k; f : X → RN
is a PWA map; e ∈ N (0, σe2 ) is an independently, identically
distributed and zero mean Gaussian noise with standard
deviation σe ; s is the number of modes; Ai , bi are the
parameters of the i-th mode (Ai , i = 1, · · · , s is assumed
to be nonsingular in this paper); and all s modes together
constitute a complete partition of X .
A transition system is a tuple T = (Q, δ, O, o), where Q
is the state space; δ : Q → 2Q (2Q is the powerset of Q)
is a transition map assigning a state q ∈ Q to its next state
q 0 ∈ Q; O is the set of observations; and o : Q → O is
an observation map assigning each q ∈ Q an observation
o(q) ∈ O [1]. We denote a region of the state space as
P ⊂ Q. The embedding transition system of a PWA system
S described by Eqn. (1) is a tuple Te = (Qe , δe , Oe , oe ),
where Qe = ∪i∈I Xi ; δe : x → x0 iff there exists i ∈ I such
that the transition from x to x0 satisfies Eqn. (1); Oe = I;
and oe (x) = i iff x ∈ Xi [24].
For a transition system T , including embedding transition
systems of PWA systems, the successor of a region P ⊂ Q
is define as the set of states that can be reached from the
states in P in one step, i.e., P ost(P ) = {q ∈ Q | ∃p ∈
P with p → q}. The predecessor of a region P ⊂ Q can be
defined similarly as P re(P ) = {q ∈ Q | ∃q ∈ P with q →
p}. A state q ∈ Q is called reachable if there exists a finite
execution ending at q. We denote all the reachable states of
T by Reach(T ). Given an LTL formula φ over O and a
system T , if all the traces originating from a region P ⊂ Q
satisfy φ, then we denote the situation as T (P ) |= φ. Let
XTφ = {X ⊂ Q, T (X) |= φ} denote the largest region of T
from which all the traces satisfy φ [1].
The reachability metric over two transition systems T1 and
T2 is defined as [25]:
d(T1 , T2 ) = h(Reach(T1 ), Reach(T2 )),
where h is the Hausdorff distance. Given two transition
systems T1 and T2 with the same observation set O and
a reachability metric d defined over them, a relation Sσ ⊆
Q1 × Q2 is called a σ−approximate simulation relation of
T1 by T2 [25] if for all (q1 , q2 ) ∈ Sσ :
• d(o(q1 ), o(q2 )) ≤ σ,
0
0
• ∀q1 = P ost(q1 ), there exists q2 = P ost(q2 ), such that
0
0
(q1 , q2 ) ∈ Sσ .
Moreover, T1 is said to be σ−approximately simulated by
T2 , denoted T1 ≺σ T2 .

III. P ROBLEM S TATEMENT
Formally, in this paper, we wish to solve the following
problem:
Problem 1: Given a PWA system S with unknown dynamics, an LTL specification φ, and a bound σ > 0, find a
finite transition system T̂ such that p(T̂ ≺σ T ) > 1 − δ,
where p(.) stands for probability, T is the true abstract
transition system of S, and δ is bounded.
Remark 1: By unknown dynamics, we mean that the
following system characteristics are unknown: (i) the number
of modes, s, (ii) the parameters related to the dynamics of
each mode, {Ai , bi , i = 1, · · · , s}, (iii) the parameters related
to the partitions (regions) of the state space, {Hi , Ki , i =
1, · · · , s}, and (iv) the standard deviation of the Gaussian
noise, σe . But we assume that our algorithm, which will
be presented in the next section, can use the system as a
black-box simulator to generate samples. This is a reasonable
assumption since during system design and testing phases,
engineers can always get access to a full-scale system model,
a scaled system model, or a computer simulation to generate
samples [26], [21].
Remark 2: Notice that the requirement p(T̂ ≺σ T ) >
1 − δ is inspired by the concept of probably approximately
correct (PAC) models in machine learning [27]. It simply
says that the probability that the transition system T̂ (obtained by using our algorithm) is a σ-approximate simulation
by the PWA system S is higher than 1 − δ. In other words,
given a PWA system with unknown dynamics, we intend to
find out its approximate abstract transition system T̂ with a
high enough confidence.
IV. DATA -D RIVEN A BSTRACTION A LGORITHM
Black-Box PWA
System S

Sample Pairs (xi,yi)

System Identification

YES
Refine?

NO

LTL Formula φ

System Abstraction

NO

Active Sampling

YES
Terminated?

Abstraction T

Fig. 1: Architecture of our data-driven abstraction algorithm.
Fig. 1 illustrates the basic architecture of our data-driven
abstraction algorithm to solve Problem 1. The inputs of the
algorithm are a black-box PWA system S with unknown
dynamics and an LTL specification φ; the output of the
system is a transition system T . The algorithm can be
roughly divided into three components: system identification,
system abstraction, and active sampling. The goal of the

system identification component is to derive an estimated
PWA model Ŝ based on the data sampled from the black-box
system S (serving as a simulator); the goal of the system
abstraction component is to derive a transition system T
given the identified PWA model Ŝ and the specification
φ; one important procedure of the system identification
component is refinement, which refines the estimated model
Ŝ, until no significant improvement can be achieved, based
on the currently available data and the current abstraction
T ; finally, if no satisfactory abstraction T can be found
after the refinement, the active sampling component will be
implemented to draw new data points with the help of the
black-box simulator S.
In the following sub-sections, we will present each of the
three components. Proofs regarding the effectiveness of our
algorithm will be provided at the end of the section.

Algorithm 1: System Identification Initialization
Input : A bound σ̂, a set of samples
D := {(yk , xk )}, k ∈ {1, · · · , K}, a ratio
0 < r < 1, and a large number J
Output: Âi , b̂i , Di , Ĥi , K̂i , i = 1, · · · , s
1: Set l = 0 and I0 = {1, · · · , K};
2: for |Il+1 | ≥ rK do
3:
Randomly generate a set of parameters
{Aj , bj }, j ∈ {1, · · · , J};
4:
Construct sets
Σj = {k yk − (Aj xk + bj ) k≤ σ̂, k ∈ Il } for each
j ∈ {1, · · · , J};
5:
Set l = l + 1 and Σmax = argmaxJj=1 |Σj |;
6:
With Σmax , estimate Âl and b̂l by solving
|Σmax |

(Âl , b̂l ) = argmin
A,b

A. System Identification
Given a black-box PWA system S, or subsequently a set
of K samples D := {yk , xk }, k = 1, · · · , K, the system
identification component identifies a PWA model, specified
by the number of modes s as well as the mode parameters
Âi , b̂i , Ĥi and K̂i with i = 1, · · · , s. The problem itself is
a well-studied problem. Specifically, we need to find (i) a
minimum positive integer, s, (ii) a set of parameter matrices,
{Âi }si=1 and {Ĥi }si=1 , and (iii) a set of parameter vector
{b̂i }si=1 and {K̂i }si=1 (notice that {Ĥi }si=1 and {b̂i }si=1
together constitute a complete partition {Xi }si of the PWA
system’s state space S), such that the estimated parameters
are the solution of the following minimization problem:
(Âi , b̂i , Ĥi , K̂i ) =

K
1 X
c(yk − fˆ(xk )) (2)
(xk ,yk )∈(D∩Xi ) K

argmin

k=1

where fˆ(.) is specified by Âi , b̂i , Ĥi and K̂i , i = 1, · · · , s
and c is a given penalty function, which is chosen to be
c(·) = || · ||2 in this paper. Notice that solving the identification problem involves the simultaneous solving of two subproblems, data classification and parameter estimation. Once
the data points have been classified into clusters {Di }si=1
such that (yk , xk ) ∈ Di , i.e., (yk , xk ) is attributed to the ith mode, mode parameters Âi , b̂i , Ĥi and K̂i can be easily
estimated by solving Eqn. (2).
Our system identification component is modified from the
method proposed in [22]. It consists of two main procedures:
initialization and refinement. One major difference between
our method and the one in [22] is that we utilize the current
abstract transition system to guide the refinement.
1) Initialization: The pseudo code of the initialization
procedure is shown in Alg. 1. The steps are rather selfexplanatory. Here we would like to provide a few simple
comments for clarification. We need to randomly generate
a matrix Âj and a vector b̂j (Line 3) in each loop, which
is quite inefficient; thus the termination condition |Il+1 | ≥
rK (Line 2) can be set loosely, i.e., with a rather large
r. Such a practice is reasonable, given the fact that the
initialization procedure is only meant to generate some good

X

c(yk − (Axk + b));

k=1

Set Dl = {k ∈ Il :k yk − (Âl xk + b̂l ) k≤ σ̂} and
Il+1 = Il \ Dl ;
8: Find boundaries specified by {Ĥ}li=1 and {K̂}li=1
between sets {Di }li=1 ; set s = l.
7:

enough partitions (modes), which will be further refined in
the refinement procedure. As for the boundaries specified
by {Ĥ}li=1 and {K̂}li=1 between sets {Di }li=1 (Line 8),
standard support vector machines (SVM) regression methods
[28] can be used to compute them.
2) Refinement: Given the random nature of the way
Âi , b̂i , i = 1, · · · , s are generated in Alg. 1, it is quite
unlikely that we are able to identify all the correct modes
with the initialization procedure. Potentially there are two
main issues:
• Undecidable data points: these are the data points that
belong to more than one mode, i.e., they satisfy k yk −
(Âi xk + b̂i ) k≤ σ̂ for more that one i = 1, · · · , s;
• Unfeasible data points: these are the data points that
don’t belong to any mode, i.e., there is no i = 1, · · · , s
such that k yk − (Âi xk + b̂i ) k≤ σ̂.
We refine the identified system model by eliminating these
two types of data points as follows:
For undecidable points, we can reassign them based on
their maximum likelihood with respect to all modes, i.e., we
assign each undecidable point (xk , yk ) to the optimal mode
according to the following rule:
CR(k) = argmax CC(xk , yk )i /CC(xk , yk ),

(3)

i=1,··· ,s

where CC(xk , yk ) is the total number of data points that
are within D and inside the hyper ball region centered at
(xk , yk ) and with a predefined radius ρ, and CC(xk , yk )i is
the number of data points that are in CC(xk , yk ) and belong
to mode i.
For unfeasible points, given the current abstract transition
system T , we discard those points that meet the following

simply using the definition of embedding transition system.
Third, an observation map oe is created by partitioning the
state space of the system Te into uniform grids. Fourth, given
the system Te , the observation map oe , and the LTL formula
φ (or its corresponding Buchi automaton Bφ ), an initial
transition system T0 is constructed by following standard
abstraction procedures, such as those prescribed in [1]. Fifth,
a product automaton P is constructed as P = T0 × Bφ ,
which concerns both the initial transition system T0 and
the specification φ. The product automaton is a tuple P =
(Sp , Sp0 , δp , Fp ), where Sp is the set of states, Sp0 is the
set of initial states, δp is the transition map, and Fp is the
acceptance condition. Finally, refinement is conducted by
solving a deterministic Rabin game.

Algorithm 2: System Identification Refinement
Input : Current abstract transition system T ,
parameters Âi , b̂i , Di , Ĥi , K̂i , i = 1, · · · , s
obtained in Alg. 1, and a bound σ̂
Output: Âi , b̂i , Ĥi , K̂i , i = 1, · · · , s
Set β ≥ 0, µ ≥ 0, κ ≥ 0, 0 < θ < 1, l = 1;
while Not terminated do
1: Compute (i∗ , j ∗ ) = argmin1≤i<j≤s βi,j with
βi,j =k Âi − Âj k;
2: if βi∗ ,j ∗ ≤ θ t β then
3:
Merge modes i∗ and j ∗ ; Set s = s − 1;
4:
Recompute Di∗ = {k ∈ {1, · · · , K} :
k yk − (Âi∗ xk + b̂i∗ ) k≤ σ̂};
5: Use CR() and rule Dis() to reassign points;
/ ∗ CR() and Dis() are defined in the text as
Eqn. (3) and Eqn. (4), respectively∗/
6: Compute i∗ = argmini=1,··· ,s |Di |/|D|;
7: if |Di |/|D| ≤ θ t µ then
8:
Discard mode i∗ ; let s = s − 1;
9:
Go to Step 4;
s
10: Store {Âi }si=1 as {Âold
i }i=1 ;
11: Update Âi , b̂i , Ĥi , K̂i with new {Di }si=1 ;
12: if k Âi − Âold
k≤ κ then
i
13:
Terminated.
14: else
15:
l = l + 1;
end

Algorithm 3: System Abstraction Refinement
Input : Current abstract transition system T , current
product automaton P, an initial state q of T ,
and a ratio 0 < η < 1
Output: Refined abstract transition system T̂ and
refined product automaton P̂
Initialize T̂ = T , P̂ = P, and Su = ∅;
while |Su | ≥ η|Q| do
1:
2:

3:
4:
5:
6:

condition:
Dis(k)
= {(xk , yk ) ∈ D|d(yk , o−1 (P ost(o(xk )))) ≥ σ̂},

(4)

where d(., .) is the Hausdorff distance, o(xk ) maps a continuous state xk to a discrete state of the transition system T (it
would be helpful for the readers to review the definitions in
Section II), P ost(o(xk )) maps the region corresponding to
o(xk ) to its successor region, and finally o−1 (P ost(o(xk ))))
finds the set of continuous, PWA states corresponding to the
region P ost(o(xk )). The remaining unfeasible data points
are reassigned according to the rule CR(·) (Eqn. (3)). The
pseudo code of the whole refinement procedure is shown in
Alg. 2.
B. System Abstraction
Given an estimated PWA model Ŝ (or fˆ), parameterized
by Âi , b̂i , Ĥi and K̂i with i = 1, · · · , s, and an LTL
formula φ, the goal of the system abstraction component is
to generate a good enough abstraction T . We roughly follow
the approximate abstraction procedures described in [24] to
design and implement the system abstraction component.
Here we are just going to provide a rough outline of the
abstraction algorithm. Interested readers can refer to [24] for
more details. First, a deterministic Buchi automaton Bφ is
constructed from the formula φ. Second, the corresponding
embedding transition system Te is constructed for Ŝ by

7:

8:
9:
10:
11:

Compute S> and S⊥ for P̂;
Set Su := Pp \ (S> ∪ S⊥ );
/∗See text for the definitions of S> , S⊥ ,
and Su ∗ /
for all (q, g) ∈ Su do
if q ∈ Su then
Set q := q;
for all (∃qr ∈ q, q 0 ∈ P ost(qr ) and
(qr ∩ P re(q 0 )) 6= ∅) do
Construct states q1 , q2 such that
q1 := qr ∩ P re(q 0 ),
q2 := qr \ P re(q 0 );
q := (q \ qr ) ∪ {q1 , q2 };
q := q;
Update δ and o for T̂ ;
Update P̂ and T̂ .

end
The pseudo code of the abstraction refinement procedure is
shown in Alg. 3. Given the current abstract transition system
T and the current product automaton P, a state q ∈ Q of T
falls into one of the following three categories:
• S> : the set of states from which all traces are accepted
by P,
• S⊥ : the set of states from which no trace is accepted
by P,
• Su : the set of states from which some but not all traces
are accepted by P.
The goal of the refinement is thus to eliminate Su . This leads
to the termination condition of the refinement procedure as
|Su | < r|Q|, i.e., the refinement will be terminated once the
volume of Su , |Su |, is smaller than a fraction r (specified

by the user) of the volume of Q (the state space of T ).
C. Active Sampling
The system identification component and the system abstraction component described in the last two sub-sections
are based on a fixed data set D := {(yk , xk )}K
k=1 . It is
quite obvious that the quality of the system identification
and the system abstraction depends on the quality of the data
set. To improve the quality of the system identification (in
other words, to decrease the number of data points needed
for the system identification), we use an active learning
algorithm developed by our group [29] to sample high quality
(or “informative”) data points for the system identification
component after the initial round (see Fig. 1).
The strategy to find the next data point to sample for round
t + 1 has two steps. In the first step, the best candidate for
each mode is identified as follows:
1/2

xi := argmax(Ψi,t (x) + λt %i,t (x))

(5)

x∈X̂i

where Ψi,t (x) is the Gaussian process regression mean of the
prediction error defined over the data points in Di , %i,t (x)
is the Gaussian process regression variance of the prediction
error defined over the data points in Di , and λt is a regularization factor. In the second step, the active learning algorithm chooses mode i∗ = argmini=1,··· ,s (maxx∈X̂i Ψi,t (x))
and the corresponding best candidate xi∗ to sample.

second condition in the definition of approximate relation is
satisfied. Therefore, T̂ ≺ε T .
Lemma 2: Given an estimation fˆ(.) (or Ŝ) of the PWA
system f (.) (or S), if k fˆ(x) − y k≤ η1 holds with
probability > 1 − α1 , where α1 ∈ (0, 1), η1 > 0, then
the system abstraction component described in Section IVB can derive an abstract transition system T̂ such that it
is (η1 + ε)−approximately simulated by the real abstract
transition system with a probability greater than 1 − α1 .
Proof: Here for the sake of clarify, the notations in
this proof that are slightly different from the ones used
in the other parts of the paper. Let’s use f and fˆ to
denote the real and estimated dynamics of the PWA system,
respectively. Moreover, let’s call their true abstract transition
systems as Treal and T̂real , respectively. Finally, let’s call the
approximate abstract transition systems obtained by using
the system abstraction component described in Section IVB as Tapp and T̂app , respectively. According to Lemma
1, we have Tapp ≺ε Treal and T̂app ≺ε T̂true . Then set
k eB k=k fˆ(x) − y k≤ η1 , such that y = fˆ(x) + y − fˆ(x) =
fˆ(x) + eB , which can be seen a PWA system with bounded
noise, following [24], we have T̂real ≺η1 Treal . Since the
relationship ≺. is transitive, we have T̂app ≺η1 +ε Treal . The
conclusion regarding probability follows easily.
Lemma 3: Provided with Assumption 1, k fˆ(x) − y k≤
η2 holds with
≥ 1 − α2 , where α2 = 1 −
R η2 probability
2
√ 1
exp(−υ
/(2(σ
e + C)))dυ.
−η2
2π(σe +C)

D. Theoretical Results Regarding the Effectiveness of Our
Algorithm
We add the following assumption regarding the performance of the PWA system identification.
Assumption 1: Assume the prediction error of the system
identification component described in Section IV-A can be
characterized by a zero mean Gaussian with a bounded
variance, i.e., f (x) − fˆ(x) ∼ N (0, σp (x)2 ) and σp (x) ≤ C,
where f (x) is the real PWA dynamics and fˆ(x) is the
estimated PWA dynamics.
Remark 3: The system identification of PWA system is
still an open problem and has been proven to be NP-hard
[30], [31]. In [22], the authors were able to demonstrate that,
for a fixed data set, the error is bounded, i.e., |f (x)− fˆ(x)| <
σ for any σ > 0. Thus, we believe our assumption here, even
though unproven, is still reasonable.
We have the following three lemmas regarding the integration of identification and abstraction (without active sampling
in the loop).
Lemma 1: Given a PWA system S with known dynamics,
for any bound ε > 0, Alg.3 can derive an abstract transition
system T̂ that is ε−approximately simulated by the real
abstract transition system T of S, i.e., T̂ ≺ε T .
Proof:
Consider a relationship V (q1 , q2 )
=
{(q1 , q2 )|d(q1 , q2 ) ≤ ε}, and set r|Q| ≤ ε. As (q1 , q2 ) ∈ V
implies d(q1 , q2 ) ≤ ε, the first condition in the definition of
approximate relation is satisfied. Then for all d(q1 , q2 ) < ε,
the conclusion in [24] shows that the result of Alg.3 can
guarantee that d(P ost(q1 ), P ost(q2 )) ≤ |Su | ≤ ε. Then the

Proof: Set eT := y− fˆ(x) = f (x)+e− fˆ(x), where e is
the true Gaussian noise, e ∼ N (0, σe2 ). Since f (x) − fˆ(x) ∼
N (0, σp (x)2 ), we have eT ∼ N (0, σp2 +σe2 ). The probability
that k fˆ(x) − y k≤ η2 can then be computed. When σp = C,
the probability reach
value, such that α2 =
R η2 the minimum
2
1− √ 1
exp(−υ
/(2(σ
+ C)))dυ.
e
−η2
2π(σe +C)

We have the following two lemmas regarding the performance of our active sampling component.
Lemma 4: For any ρ ∈ (0, 1), if λt = 2B +
300γt log3 (t/ρ), where γt is the maximum information gain
defined in [29], then
1/2

k Ψi,t (x) k≤ λt %i,t (x)

(6)

holds with probability > 1 − ρ.
Proof: The lemma can be proved by following similar
steps as the proof of Lemma 2 in [29].
Lemma 5: For any σ > 0, if λt = 2B + 300γt log3 (t/ρ),
where t is the number of points in Di , then there exists
t ≤ T,
1/2
λt %i,t (x) ≤ σ, ∀x ∈ Di
(7)
holds with finite number of sample time T .
Proof: According to Theorem 5 in [32], for any
bounded Xi ⊂ X , the information gain for exponential kernel is γT = O((logT )d+1 ), where d is
the
of X . Lemma 7.1 in [32] shows that
PT dimension
2α
−2 2
min{σ
%
i,t−1 (x), α} ≤ log(1+α) γT , ∀α > 0. As the
t=1
information gain γT is convergent, such that the variance
%i,t (x) will be convergent to zero. As λt is bounded,

1/2

λt %t (x) will be convergent to zero, then the lemma has
been proved.
Finally, we can prove that the algorithm presented in this
paper can solve Problem 1 with the relatively moderate
Assumption 1.
Theorem 1: Given a PWA system S with unknown dynamics, for any σ > ε > 0, the algorithm described in
this paper (with the assumption that the standard deviation
of the prediction error for the system identification component is bounded by C) can obtain an approximate abstract
transition system T̂ that is σ-approximately simulated by
the true abstract transition system T of S with a bounded
probability Rthat is greater than 1 − δ, where δ = 1 −
σ−ε
√ 1
exp(−υ 2 /(2(σe + C)))dυ.
−σ+ε
2π(σe +C)

Proof: Set δ = α2 and η2 = σ − ε. According to
Lemma 3, Lemma 4, and Lemma 5, the system identification
component together with the active sampling component can
achieve an estimation k fˆ(x) − y k≤ η2 with probability
> 1 − α2 . Then with this estimation result, according to
Lemma 1 and Lemma 2, the system abstraction component
can generate an abstract transition system T̂ that is (η2 +
ε)−approximately simulated by the real abstract transition
system T with probability > 1 − δ.
V. C ASE S TUDY
In this section, we will use a soft robot driven by series
pneumatic artificial muscles as an example to demonstrate
our algorithm.
A. Model and LTL Specification

Fig. 2: (a) Geometrical model of the soft robot driven by
series pneumatic artificial muscles (sPAM) used in the case
study. (b) Simplified planar robot with x = (x(1), x(2))T as
the coordinate vector of the moving platform.
In this case study, we will focus on a particular type of
soft robots, those that are driven by series pneumatic artificial
muscles (sPAM) [33]. Fig. 2 shows an example of such a
robot. It has two polyethylene tubing sPAMs, each of which
is controlled by a corresponding actuator with pressurized
air. Even though the dynamics of such a robot is nonlinear,
in [34], the authors have shown that its closed-loop behavior
can be approximated by piecewise affine dynamics. Here let’s

assume that the robot under investigation has dynamics as
follows:
xk+1 = f (xk ) + e
(8)
and
 

1
0


x if 0 ≤ xk (1) ≤ 0.3

0 0.98 


f (x) =
0.83 0.12
0.01


x
+
if xk (1) ≥ 0.3

0.12 0.81
0.03
where x = (x(1), x(2))T is the coordinate of the moving
platform and e is a Gaussian noise with a zero mean and a
standard deviation of 0.1. Please keep in mind that the model
is unknown to our algorithm but the model itself can be used
by our algorithm as a simulator to generate samples.
The specification that needs to be verified is written as
an LTL formula φ := (π1 ∧ ♦π2 ), where “π1 := x(1) is
below 0.3", “π2 := x(2) is above 0.6", and  is the temporal
operator “Always”, ♦ is the temporal operator “Eventually”.
Put together, φ specifies that “it should always be true that
x(1) is below 0.3 and eventually x(2) is above 0.6”.
B. Implementation Results
The algorithm proposed in this paper is implemented as a
Matlab tool. The tool takes an LTL formula and a black-box
PWA system as inputs and it outputs an abstract transition
system. The active learning used in the tool is the Gaussian
Process Adaptive Confidence Bound (GP-ACB) algorithm
proposed and implemented by our group in [29]. Moreover,
a Gaussian kernel function is used in GP-ACB. In order to
demonstrate the effectiveness of our algorithm, here we will
provide two sets of implementation results, one regarding the
system identification component and the other one regarding
the entire algorithm.
1) Performance of the System Identification Component:
Here we compare the performance of our system identification component with an existing state-of-the-art black-box
system identification tool called the Hybrid Identification
Tool (HIT) [35]. One thing we would like to point out here
is that the evaluation of our system identification component
is conducted in conjunction with other components, e.g.,
abstraction and active sampling, in the loop. It should be
expected that, at worst, our component should have the
same performance as the HIT. However, it may also be
expected that factors such as active-learning based sampling
and abstraction-guided refinement can potentially improve
the identification performance, which turns out to be the case,
at least for this particular case study.
We use two metrics to quantify system identification
errors: parameter estimation error and region estimation
error. Given the parameters of a real PWA model (unknown
to the investigated algorithms), the parameter estimation
error is the sum of the Euclidean distances between the real
parameters and the estimated parameters. The region estimation error is defined the sum of the following Hausdorff
distance:
e(X̂i , Xi ) = max{ sup inf d(x, y), sup inf d(x, y)},
x∈Xi y∈X̂i

y∈X̂i

x∈Xi

(a)

(b)

Fig. 3: Comparison results of the system identification component used in this paper (green) and the HIT algorithm
(red): (a) Average parameter estimation error with respect
to the number of samples; (b) Average region (partition)
estimation error with respect to the number of samples.

where e(X̂i , Xi ) the region estimation error related to the
ith mode with Xi as the real region or partition and X̂i
as the estimated one, d(x, y) is the Euclidean distance
between x and y, and sup and inf stand for supremum and
infimum, respectively. To calculate the Hausdorff distance,
we randomly generate 100 samples inside the state space X̂i
and the state space Xi .
The comparison results based on 5 trials are shown in Fig.
3. It shows that, averagely speaking, the system identification
component proposed in this paper has a comparable or
better performance than HIT. Particularly, Fig. 3b shows
that our algorithm has a faster convergent rate regarding
the region estimation error. This is probably due to the fact
that, generally speaking, active learning, which is used in our
algorithm, out-performs its randomly sampling counterpart,
which is used in HIT.
2) Performance of the Entire Algorithm: As the true
abstract transition system of the system, described by Eqn.
(8), is unknown, here we use the abstract transition system

obtained by using the algorithm proposed in [24] as a
benchmark. We will call this abstract transition system as
T ∗ . In [24], the authors have shown that even though their
algorithm cannot attain the true abstract transition system, it
still can get an abstract transition system that is arbitrarily
close to the real one. Of course, we should point out that, in
order to get this abstract transition system T ∗ , the algorithm
in [24] should have access to the system model, i.e., T ∗
is generated with a known model. In parallel, we run our
algorithm to extract an abstract transition system T̂ without
access to the dynamics of the model, i.e., T̂ is generated
with a black-box system with unknown dynamics. Then in
order to demonstrate the effectiveness of our algorithm, we
need to show that T̂ should approach T ∗ . This turns out to
be the case, at least for this particular case study.
Based on the system dynamics, Eqn. (8), and the LTL
specification φ, we implement the abstraction algorithm
proposed in [24] and obtain an abstract transition system
T ∗ , which will be used as a benchmark. Then the algorithm
proposed in this paper is applied to the same system, but with
unknown dynamics, and the same LTL specification. The
output of the algorithm is another abstract transition system
T̂ . Here we use a metric that is inspired by the concept of
approximate simulation to quantify the difference between
the two transition systems. To be more specific, we set the
abstraction error to σ if T ∗ is σ−approximately simulated by
T̂ . Moreover, the metric σ is normalized to σ̄ by the volume
of the state space, i.e., σ̄ := σ/|X |.
TABLE I: Comparison result with respect to different number
of samples and a fixed number of refinement steps, which is
set to 20.

σ̄

Number of samples
20
40
60
0.046
0.042
0.035

TABLE II: Comparison result with respect to different number of refinement steps and a fixed number of samples in
active sampling component, which is set to 10.

σ̄

Refinement Steps
5
10
20
0.077
0.068
0.050

The comparison results are shown in Table I and Table II.
Table I shows the comparison results with respect to different
number of samples and a fixed number of refinement steps in
the system abstraction component. Table II shows the comparison results with respect to different number of refinement
steps and a fixed number of samples in the active sampling
component. The tables show the average normalized errors σ̄
based on 5 trials. The results show that the larger the number
of refinement steps, the smaller the abstraction error; and the
larger the number of samples, the smaller the abstraction
error. The results also show that, even with a black-box
system, our algorithm can attain an approximate abstract
transition system T̂ that is close to the benchmark abstract
transition system T ∗ .

Fig. 4: State space partitions of the abstract transition system
obtained by using our algorithm for the case study. The
number of refinement steps and the number of samples in
the active sampling component are both set at 20.

VI. C ONCLUSIONS
In this paper, we proposed a data-driven approximate
abstraction algorithm for piecewise affine systems with unknown dynamics. We demonstrated both theoretically and
empirically that given a black-box PWA system and an LTL
specification, we were able to derive an abstract transition
system that is approximately simulated by the true abstract
transition system. We demonstrated the effectiveness of our
proposed algorithm with a soft robot as a case study.
ACKNOWLEDGMENTS
This work was partially supported by the Hyundai Motor
Company.
R EFERENCES
[1] C. Baier, J.-P. Katoen, and K. G. Larsen, Principles of model checking.
MIT press, 2008.
[2] R. Alur, Principles of cyber-physical systems. MIT Press, 2015.
[3] D. A. Peled, Software reliability methods.
Springer Science &
Business Media, 2013.
[4] M. Kamali, L. A. Dennis, O. McAree, M. Fisher, and S. M. Veres,
“Formal verification of autonomous vehicle platooning,” Science of
Computer Programming, 2017.
[5] H. Abbas, B. Hoxha, G. Fainekos, and K. Ueda, “Robustness-guided
temporal logic testing and verification for stochastic cyber-physical
systems,” in Cyber Technology in Automation, Control, and Intelligent
Systems (CYBER), 2014 IEEE 4th Annual International Conference
on. IEEE, 2014, pp. 1–6.
[6] Y. K. Lopes, S. M. Trenkwalder, A. B. Leal, T. J. Dodd, and
R. Groß, “Supervisory control theory applied to swarm robotics,”
Swarm Intelligence, vol. 10, no. 1, pp. 65–97, 2016.
[7] A. Ulusoy, T. Wongpiromsarn, and C. Belta, “Incremental controller
synthesis in probabilistic environments with temporal logic constraints,” The International Journal of Robotics Research, vol. 33,
no. 8, pp. 1130–1144, 2014.
[8] A. Girard, G. Pola, and P. Tabuada, “Approximately bisimilar symbolic
models for incrementally stable switched systems,” IEEE Transactions
on Automatic Control, vol. 55, no. 1, pp. 116–126, 2010.
[9] E. Haghverdi, P. Tabuada, and G. J. Pappas, “Bisimulation relations
for dynamical, control, and hybrid systems,” Theoretical Computer
Science, vol. 342, no. 2-3, pp. 229–261, 2005.

[10] P. Prabhakar, G. Dullerud, and M. Viswanathan, “Stability preserving
simulations and bisimulations for hybrid systems,” IEEE Transactions
on Automatic Control, vol. 60, no. 12, pp. 3210–3225, 2015.
[11] G. Pola, P. Pepe, and M. D. Di Benedetto, “Symbolic models for timevarying time-delay systems via alternating approximate bisimulation,”
International Journal of Robust and Nonlinear Control, vol. 25, no. 14,
pp. 2328–2347, 2015.
[12] B. Yordanov, J. Tumova, I. Cerna, J. Barnat, and C. Belta, “Temporal
logic control of discrete-time piecewise affine systems,” IEEE Transactions on Automatic Control, vol. 57, no. 6, pp. 1491–1504, 2012.
[13] E. A. Gol, X. Ding, M. Lazar, and C. Belta, “Finite bisimulations for
switched linear systems,” IEEE Transactions on Automatic Control,
vol. 59, no. 12, pp. 3122–3134, 2014.
[14] J. A. DeCastro, V. Raman, and H. Kress-Gazit, “Dynamics-driven
adaptive abstraction for reactive high-level mission and motion planning,” in Robotics and Automation (ICRA), 2015 IEEE International
Conference on. IEEE, 2015, pp. 369–376.
[15] G. Pola and M. D. Di Benedetto, “Symbolic models and control of
discrete-time piecewise affine systems: An approximate simulation
approach,” IEEE Transactions on Automatic Control, vol. 59, no. 1,
pp. 175–180, 2014.
[16] T. Dreossi, A. Donzé, and S. A. Seshia, “Compositional falsification
of cyber-physical systems with machine learning components,” arXiv
preprint arXiv:1703.00978, 2017.
[17] K. Urayama, M. C. Fu, and S. I. Marcus, “Simulation-based work load
and job release control for semiconductor manufacturing,” in Decision
and Control (CDC), 2015 IEEE 54th Annual Conference on. IEEE,
2015, pp. 7329–7334.
[18] S. A. Seshia, S. Hu, W. Li, and Q. Zhu, “Design automation of cyberphysical systems: Challenges, advances, and opportunities,” IEEE
Transactions on Computer-Aided Design of Integrated Circuits and
Systems, 2016.
[19] S. Haesaert, A. Abate, and P. M. Van den Hof, “Data-driven and
model-based verification: A bayesian identification approach,” in Decision and Control (CDC), 2015 IEEE 54th Annual Conference on.
IEEE, 2015, pp. 6830–6835.
[20] A. Kozarev, J. Quindlen, J. How, and U. Topcu, “Case studies in
data-driven verification of dynamical systems,” in Proceedings of the
19th International Conference on Hybrid Systems: Computation and
Control. ACM, 2016, pp. 81–86.
[21] Z. Kong, A. Jones, and C. Belta, “Temporal logics for learning and
detection of anomalous behavior,” IEEE Transactions on Automatic
Control, vol. 62, no. 3, pp. 1210–1222, 2017.
[22] A. Bemporad, A. Garulli, S. Paoletti, and A. Vicino, “A bounded-error
approach to piecewise affine system identification,” IEEE Transactions
on Automatic Control, vol. 50, no. 10, pp. 1567–1580, 2005.
[23] A. Garulli, S. Paoletti, and A. Vicino, “A survey on switched and
piecewise affine system identification,” IFAC Proceedings Volumes,
vol. 45, no. 16, pp. 344–355, 2012.
[24] B. Yordanov, J. Tumova, I. Cerna, J. Barnat, and C. Belta, “Formal
analysis of piecewise affine systems through formula-guided refinement,” Automatica, vol. 49, no. 1, pp. 261–266, 2013.
[25] A. Girard and G. J. Pappas, “Approximation metrics for discrete
and continuous systems,” IEEE Transactions on Automatic Control,
vol. 52, no. 5, pp. 782–798, 2007.
[26] X. Jin, A. Donzé, J. V. Deshmukh, and S. A. Seshia, “Mining
requirements from closed-loop control models,” IEEE Transactions on
Computer-Aided Design of Integrated Circuits and Systems, vol. 34,
no. 11, pp. 1704–1717, 2015.
[27] L. Valiant, Probably Approximately Correct: Nature’s Algorithms for
Learning and Prospering in a Complex World. Basic Books (AZ),
2013.
[28] C. Cortes and V. Vapnik, “Support-vector networks,” Machine learning, vol. 20, no. 3, pp. 273–297, 1995.
[29] G. Chen, Z. Sabato, and Z. Kong, “Active learning based requirement
mining for cyber-physical systems,” in Decision and Control (CDC),
2016 IEEE 55th Conference on. IEEE, 2016, pp. 4586–4593.
[30] L. Ljung, “Perspectives on system identification,” Annual Reviews in
Control, vol. 34, no. 1, pp. 1–12, 2010.
[31] F. Lauer, “On the complexity of piecewise affine system identification,”
Automatica, vol. 62, pp. 148–153, 2015.
[32] N. Srinivas, A. Krause, S. M. Kakade, and M. Seeger, “Gaussian
process optimization in the bandit setting: No regret and experimental
design,” arXiv preprint arXiv:0912.3995, 2009.

[33] F. Daerden and D. Lefeber, “The concept and design of pleated
pneumatic artificial muscles,” International Journal of Fluid Power,
vol. 2, no. 3, pp. 41–50, 2001.
[34] G. Andrikopoulos, G. Nikolakopoulos, I. Arvanitakis, and S. Manesis,
“Piecewise affine modeling and constrained optimal control for a pneumatic artificial muscle,” IEEE Transactions on Industrial Electronics,
vol. 61, no. 2, pp. 904–916, 2014.
[35] G. Ferrari-Trecate, “Hybrid identification Toolbox (HIT),” 2005.

