Fast Falsification of Hybrid Systems using Probabilistically Adaptive Input

TL;DR

This paper introduces a probabilistically adaptive algorithm for rapidly finding falsifying inputs in hybrid systems, significantly improving efficiency over existing methods by dynamically adjusting search granularity.

Contribution

The paper presents a novel 'Las Vegas tree search' algorithm that adaptively refines input space discretization based on local complexity for effective falsification.

Findings

Outperforms existing falsification techniques on standard benchmarks

Efficiently finds falsifying inputs with fewer evaluations

Adapts to problem difficulty through probabilistic discretization

Abstract

We present an algorithm that quickly finds falsifying inputs for hybrid systems, i.e., inputs that steer the system towards violation of a given temporal logic requirement. Our method is based on a probabilistically directed search of an increasingly fine grained spatial and temporal discretization of the input space. A key feature is that it adapts to the difficulty of a problem at hand, specifically to the local complexity of each input segment, as needed for falsification. In experiments with standard benchmarks, our approach consistently outperforms existing techniques by a significant margin. In recognition of the way it works and to distinguish it from previous work, we describe our method as a "Las Vegas tree search".

Tables

Table 1. Table 1 . Size of 𝒜 l subscript 𝒜 𝑙 \mathcal{A}_{l} for input dimensionality  n 𝑛 n and level  l 𝑙 l given that the size of 𝐭 l subscript 𝐭 𝑙 \mathbf{t}_{l} is one.

$n$	$l=0$	$1$	$2$	$3$	$4$	$5$	$6$	$7$	$8$	$9$	$10$
$2$	4	4	9	20	44	96	208	448	960	2048	4352
$3$	8	12	30	73	174	408	944	2160	4896	11008	24576

Table 2. Table 2 . Experimental Results comparing the success rate out of of 50 trials (“succ.” , higher is better) and number of iterations averaged over successful trials (“iter.”, lower is better, where M=arithmetic mean, SD=standard deviation) of uniform random sampling, Breach/CMA-ES, and FalStar / algorithm   1 for a maximum of 300 iterations per trial. The best results for each requirement are highlighted.

				Breach:			FalStar:
	Random			CMA-ES			aLVTS
	succ.	iter.		succ.	iter.		succ.	iter.
Formula	/50	M	SD	/50	M	SD	/50	M	SD
AT1	43	106.6	83.9	50	39.7	23.6	50	8.5	6.7
AT2 $\left(i=3\right)$	50	41.0	36.7	50	13.2	9.1	50	33.4	27.5
AT2 $\left(i=4\right)$	49	67.0	60.8	6	17.8	15.9	50	23.4	22.5
AT3	19	151.1	98.1	50	145.2	63.0	50	86.3	52.1
AT4 $\left(a\right)$	36	117.3	71.8	50	97.0	47.7	50	22.8	10.6
AT4 $\left(b\right)$	2	117.7	9.2	49	46.7	58.0	50	47.6	23.5
Summary AT	199	95.3	47.9	255	42.8	29.0	300	29.2	19.4
AFC27	15	129.1	90.8	41	121.0	49.3	50	3.9	4.3

Equations

$$\phi =\neg \phi \mid \phi \vee \phi \mid \phi \wedge \phi \mid \phi {\text{U}}_I\phi \mid {\square }_I\phi \mid {◊}_I\phi \mid \mu ,$$

$$\rho (\mu ,\mathbf{y},t)$$

$$\rho (\neg \phi ,\mathbf{y},t)$$

$$\rho ({\phi }_1\vee {\phi }_2,\mathbf{y},t)$$

$$\rho ({\phi }_1\wedge {\phi }_2,\mathbf{y},t)$$

$$\rho ({\phi }_1{\text{U}}_I{\phi }_2,\mathbf{y},t)$$

$$\rho ({◊}_I\phi ,\mathbf{y},t)$$

$$\rho ({\square }_I\phi ,\mathbf{y},t)$$

$$\rho (\phi ,\mathbf{y})>0⟹\mathbf{y}\models \phi \text{ and }\rho (\phi ,\mathbf{y})<0⟹\mathbf{y}/\models \phi .$$

$$\text{Find }{\mathbf{u}}^{\ast }=\underset{\mathbf{u}\in ([0,T]\to {\mathbb{R}}^n)}{\arg \min }\rho (\phi ,\mathcal{M}(\mathbf{u}))\text{ such that }\rho (\phi ,\mathcal{M}({\mathbf{u}}^{\ast }))<0.$$

$$\text{Lower:}\underline{\rho }(\phi ,\mathbf{y})=\underset{{\mathbf{y}}^{\prime }}{\min }\rho (\phi ,\mathbf{y}{\mathbf{y}}^{\prime })\text{Upper:}\overline{\rho }(\phi ,\mathbf{y})=\underset{{\mathbf{y}}^{\prime }}{\max }\rho (\phi ,\mathbf{y}{\mathbf{y}}^{\prime })$$

$$\text{Find }{\mathbf{u}}^{\ast }\text{ such that }\overline{\rho }(\phi ,\mathcal{M}({\mathbf{u}}^{\ast }))<0.$$

$${\mathbf{p}}_l=\{(2j+1)/2^l\mid j\in {\mathbb{N}}_0\le (2^l-1)/2\}.$$

$${\mathbf{v}}_{i,l}={\underline{v}}_i+{\mathbf{p}}_l\times ({\overline{v}}_i-{\underline{v}}_i).$$

$${\mathbf{b}}_l=\{(b_1,...,b_n)\in {\mathbb{N}}_0^n\mid b_1+\cdots +b_n=l\}.$$

$${\mathcal{A}}_l=\bigcup _{(b_1,\dots ,b_n)\in {\mathbf{b}}_l}{\mathbf{t}}_l\times {\mathbf{v}}_{1,b_1}\times \cdots \times {\mathbf{v}}_{n,b_n}.$$

$${unexplored}_l(\mathbf{u})$$

$${explored}_l(\mathbf{u})$$

$$w_l=\frac{|{unexplored}_l(\mathbf{u})|+|{explored}_l(\mathbf{u})|}{2^l\cdot |{\mathcal{A}}_l|},$$

$$={\square }_{[0,30]}v<120$$

$$\text{AT2}(i)$$

$$\displaystyle=\lnot\big{(}\Box_{[10,30]}v\in[50,60]\big{)}$$

$$\text{AT4}(\overline{v},\underline{\omega })$$

$$={\square }_{[11,50]}(rise\vee fall)⟹{\square }_{[1,5]}|\mu |<0.008,$$

Full text

cyber-physical system, hybrid system, testing, falsification, stochastic optimization, temporal logic, Las Vegas tree search

††copyright: rightsretained††conference: ; ; ††journalyear: 2018††ccs: Computer systems organization Embedded and cyber-physical systems††ccs: Mathematics of computing Stochastic control and optimization††ccs: Theory of computation Random search heuristics

1. Introduction

The falsification problem we consider seeks a (time-bounded) input signal that causes a hybrid system model to violate a given temporal logic specification. A popular way to address this is to first construct a “score function” (Jegourel et al., 2013) that quantifies how much of the specification has been satisfied during the course of an execution. The falsification can then be treated as an optimization problem, which can be solved using standard algorithms. This approach, especially using “robustness” (Fainekos and Pappas, 2009) as the score function, has been successfully applied, resulting in a number of now mature tools (Donzé, 2010; Annpureddy et al., 2011).

Despite its apparent success, the commonly used robustness semantics (Fainekos and Pappas, 2009) is in general not a perfect optimization function. Greedy hill climbing may lead to local optima, hence robustness is only an heuristic score function (Jegourel et al., 2014) with respect to the falsification problem. In practice, standard optimization algorithms overcome this limitation by including stochastic exploration. The most sophisticated of these can also model the dynamics of the system (e.g., (Akazaki, 2016)), in order to estimate the most productive direction of input signal space to explore. There is, however, “no free lunch” (Wolpert and Macready, 1997), and high performance general purpose optimization algorithms are not necessarily the best choice. For example, such algorithms often optimize with respect to the entire input trace, without exploiting the time causality of the problem, i.e., the fact that a good trace (one that eventually falsifies the property) may be dependent on a good trace prefix.

The contribution of this paper is a fast randomized falsification algorithm (section 3.1) that exploits the time-causal structure of the problem and that adapts to local complexity. In common with alternative approaches, our algorithm searches a discretized space of input signals, but in our case the search space also includes multiple levels of spatial and temporal granularity (section 3.2). The additional complexity is mitigated by an efficient tree search that probabilistically balances exploration and exploitation (section 3.3).

The performance of our algorithm benefits from the heuristic idea to explore simple (coarse granularity) inputs first, then gradually switch to more complex inputs that include finer granularity. Importantly, the finer granularity tends only to be added where it is needed, thus avoiding the exponential penalty of searching the entire input space at the finer granularity. While it is always possible to construct pathological problem instances, we find that our approach is very effective on benchmarks from the literature, including the Powertrain benchmark (Jin et al., 2014), which is considered a difficult challenge for falsification. Our experimental results (section 4) demonstrate that our algorithm can consistently beat the best existing methods, in terms of speed and reliability of finding a falsifying input.

2. Preliminaries

In this work we represent a deterministic black-box system model as an input/output function $\mathcal{M}\colon([0,T]\to\mathbb{R}^{n})\to([0,T]\to\mathbb{R}^{m})$. In general, $\mathcal{M}$ comprises continuous dynamics with discontinuities. $\mathcal{M}$ takes a time-bounded, real-valued input signal $\bm{u}\colon[0,T]\to\mathbb{R}^{n}$ of length $|\bm{u}|=T$ and transforms it to a time bounded output signal $\bm{y}\colon[0,T]\to\mathbb{R}^{m}$ of the same length, but potentially different dimensionality. The dimension $n$ of the input indicates that at each moment $t\leq T$ within the time horizon $T$, the value $\bm{u}(t)\in\mathbb{R}^{n}$ of the input is an $n$-dimensional real vector (analogously for the output).

We denote by $\bm{u}_{1}\bm{u}_{2}$ the concatenation of signals $\bm{u}_{1}$ and $\bm{u}_{2}$ that have the same dimensions, such that $\bm{u}_{1}{}\bm{u}_{2}\colon[0,T_{1}+T_{2}]\to\mathbb{R}^{n}$. Concatenation of more than two signals follows naturally and is denoted $\bm{u}_{1}\bm{u}_{2}\bm{u}_{3}\cdots$. A constant input signal segment is written $(t,v)$, where $t$ is a time duration and $v\in\mathbb{R}^{n}$ is a vector of input values. A piecewise constant input signal is the concatenation of such segments.

In this work we adopt the syntax and robustness semantics of STL defined in (Donzé and Maler, 2010). The syntax of an STL formula is thus given by

[TABLE]

where the logical connectives and temporal operators have their usual Boolean interpretations and equivalences, $I$ is the interval of time over which the temporal operators range, and atomic formulas $\mu\equiv f(x_{1},\dots,x_{m})\geq 0$ are predicates over the spatial dimensions of a trace. The robustness of trace $\bm{y}$ with respect to formula $\varphi$, denoted $\rho(\varphi,\bm{y})$, is calculated inductively according to the following robustness semantics, using the equivalence $\rho(\varphi,\bm{y})\equiv\rho(\varphi,\bm{y},0)$.

[TABLE]

An important characteristic of the robustness semantics is that it is faithful to standard boolean satisfaction, such that

[TABLE]

Together, these equations justify using the robustness semantics $\rho(\varphi,\mathcal{M}(\bm{u}))$ to detect whether an input $\bm{u}$ corresponds to the violation of a requirement $\varphi$. This correspondence is exploited to find such falsifying inputs through global hill-climbing optimization:

[TABLE]

Of course, finding an adequate falsifying input $\bm{u}^{*}$ is generally hard and subject to the limitations of the specific optimization algorithm used.

Sound approximations of the lower and upper bounds of the robustness of a prefix $\bm{y}$ can sometimes be used to short-cut the search. We thus define lower and upper bounds in the following way.

[TABLE]

A lower bound $\underline{\rho}(\varphi,\mathcal{M}(\bm{u})))>0$ can be used to detect that a prefix cannot be extended to a falsifying trace (e.g., after the deadline for a harmful event has passed). An upper bound $\overline{\rho}(\varphi,\mathcal{M}(\bm{u}))<0$ similarly implies $\mathcal{M}(\bm{u}{}\bm{u}^{\prime})\not\models\varphi$ for all $\bm{u}^{\prime}$, concluding that input $\bm{u}$ is already a witness for falsification (e.g., a limit is already exceeded).

3. Approach

We wish to solve the following falsification problem efficiently:

[TABLE]

Our approach is to repeatedly construct input signals $\bm{u}=\bm{u}_{1}{}\bm{u}_{2}{}\bm{u}_{3}{}\cdots$, where $\bm{u}_{i}$ is drawn from a predetermined search space of candidate input segments, $\mathcal{A}$. The choice is probabilistic, according to a distribution $\mathcal{D}$ that determines the search strategy, i.e., which inputs are likely to be tried next given a partially explored search space. The construction of each input is done incrementally, to take advantage of the potential short cuts described at the end of section 2.

Algorithm 1 “adaptive Las Vegas Tree Search” (aLVTS) codifies the high level functionality of this probabilistic approach, described in detail in section 3.1.

The effectiveness of our algorithm in practice comes from the particular choices of $\mathcal{A}$ and $\mathcal{D}$, which let the search gradually adapt to the difficulty of the problem. The set $\mathcal{A}$ (defined in section 3.2) contains input segments of diverse granularity, which intuitively corresponds to how precise the input must be in order to find a falsifying trace. The distribution $\mathcal{D}$ (defined in section 3.3) initially assigns high probabilities to the “coarsest” input segments in $\mathcal{A}$. Coarse here means that the segments tend to be long in relation to the time horizon $T$ and large in relation to the extrema of the input space. The algorithm probabilistically balances exploration and exploitation of segments, but as the coarser segments become fully explored, and the property has not been falsified, the algorithm gradually switches to finer-grained segments.

3.1. Algorithm

Algorithm 1 searches the space of input signals constructed from piecewise constant (over time) segments, which are chosen at random according to the distribution defined by $\mathcal{D}$ in algorithm 1. This distribution is a function of the numbers of unexplored and explored edges at different levels of granularity, and thus defines the probabilities of exploration, exploitation and adaptation. The precise calculation made by $\mathcal{D}$ is described in section 3.3.

As the search proceeds, the algorithm constructs a tree whose nodes each correspond to a unique input signal prefix. The edges of the tree correspond to the constant segments that make up the input signal. The root node corresponds to time 0 and the empty input signal (algorithm 1).

To each node identified by an input signal prefix $\bm{u}$ is associated a set of unexplored edges, $\mathit{unexplored}(\bm{u})\subseteq\mathcal{A}$, that correspond to unexplored input signal segments, and a set of explored edges, $\mathit{explored}(\bm{u})\subseteq\mathcal{A}$, that remain inconclusive with respect to falsification. Initially, all edges are unexplored (algorithm 1 and algorithm 1). Once an edge has been chosen (algorithm 1), the unique signal segment associated to the edge may be appended to the signal prefix associated to the node, to form an extended input signal. If the chosen edge is unexplored, it is removed from the set of unexplored edges (algorithm 1) and the extended input signal $\bm{u}{}\bm{u}^{\prime}$ is transformed by the system into an extended output signal (algorithm 1). If the requirement is falsified by the output signal ($\bm{y}$ in algorithm 1), the algorithm immediately quits and returns the falsifying input signal (algorithm 1). If the requirement is satisfied, with no possibility of it being falsified by further extensions of the signal (algorithm 1), the algorithm quits the current signal (algorithm 1) and starts a new signal from the root node (algorithm 1). This is the case, in particular, when the length of the signal exceeds the time horizon of the formula as a consequence of the definition of $\underline{\rho}$ in eq. 4. If the requirement is neither falsified nor satisfied, the edge is added to the node’s set of explored edges (algorithm 1). Regardless of whether the chosen edge was previously explored or unexplored, if the signal remains inconclusive, the extended input signal becomes the focus of the next iterative step (algorithm 1).

3.2. Definition of $\mathcal{A}$

Given an input signal segment $\bm{u}^{\prime}$ of length $t$ time units and value $(v_{1},\dots,v_{n})\in\mathbb{R}^{n}$, let $\underline{v}_{i}$ and $\overline{v}_{i}$ denote the minimum and maximum possible values, respectively, of dimension $i\in\{1,\dots,n\}$. For each integer level $l\in\{0,\dots,l_{\mathrm{max}}\}$, we define the set of possible proportions of the interval $[\underline{v}_{i},\overline{v}_{i}]$ as

[TABLE]

The numerators of all elements are coprime with the denominator, $2^{l}$, hence $\mathbf{p}_{i}\cap\mathbf{p}_{j}=\varnothing$ for all $i\neq j$. By definition, $\mathbf{p}_{0}$ also includes [math]. Hence, $\mathbf{p}_{0}=\{0,1\}$, $\mathbf{p}_{1}=\left\{{\footnotesize\frac{1}{2}}\right\}$ and $\mathbf{p}_{2}=\left\{{\footnotesize\frac{1}{4}},{\footnotesize\frac{3}{4}}\right\}$, etc. The set of possible values of dimension $i$ at level $l$ is thus given by

[TABLE]

Rather than making the granularity of each dimension independent, we interpret the value of $l$ as a granularity “budget” that must be distributed among the (non-temporal) dimensions of the input signal. The set of possible per-dimension budget allocations for level $l$ is given by

[TABLE]

For example, with $n=2$, $\mathbf{b}_{3}=\left\{(0,3),(1,2),(2,1),(3,0)\right\}$. If we denote the set of possible time durations at level $l$ by $\mathbf{t}_{l}$, then the set of possible input segments at level $l$ is given by

[TABLE]

Note that while $\mathbf{t}_{l}$ is not required here to share the granularity budget, this remains a possibility. Our implementation (section 4.1) actually specifies $\mathbf{t}_{l}$ by defining a fixed number of control points per level, $(k_{0},\ldots,k_{l_{\mathrm{max}}})$, such that the $\mathbf{t}_{l}=\{T/k_{l}\}$ are singleton sets. The sizes of various $\mathcal{A}_{l}$ for different choices of $n$ and $l$, assuming $|\mathbf{t}_{l}|=1$, are given in table 1.

In summary, an input segment $\bm{u}^{\prime}=(t,v_{1},\ldots,v_{n})\in\mathcal{A}_{l}$ has $t\in\mathbf{t}_{l}$ and corresponding budget allocation $b_{1}+\cdots+b_{n}=l$, with the value $v_{i}$ for each dimension given by $v_{i}=\underline{v}_{i}+p_{i}(\overline{v}_{i}-\underline{v}_{i})$, where $p_{i}=(2j_{i}+1)/2^{b_{i}}$, for some $j_{i}$, defines the proportion between minimum $\underline{v}_{i}$ and maximum $\overline{v}_{i}$.

By construction, $\mathcal{A}_{i}\cap\mathcal{A}_{j}=\varnothing$, for all $i\neq j$. Hence, we define

[TABLE]

The set of all possible input signal segments is given by $\mathcal{A}=\mathcal{A}_{0}\cup\mathcal{A}_{1}\cup\cdots\cup\mathcal{A}_{l_{\mathrm{max}}}$.

Figure 1 depicts the construction of $\mathcal{A}$ for two dimensions. The majority of candidate input points is concentrated on the outer contour, corresponding to an extreme choice for one dimension and a fine-grained choice for the other dimension. While this bias appears extreme, as layers are exhausted, finer-grained choices become more likely. For example, after two points from $\mathcal{A}_{0}$ have been tried, all remaining points from both levels in the second panel would be equally probable.

3.3. Definition of $\mathcal{D}$

The distribution $\mathcal{D}(\bm{u})$ is constructed implicitly. First, a granularity level $l\in\{0,\dots,l_{\mathrm{max}}\}$ is chosen at random, with probability in proportion to the fraction of edges remaining at the level, multiplied by an exponentially decreasing scaling factor. Defining the overall weight of level $l$ as

[TABLE]

level $l$ is chosen with probability $w_{l}/\sum^{l_{\mathrm{max}}}_{i=0}w_{i}$.

Having chosen $l$, one of the following strategies is chosen uniformly at random:

1. (1)

select $\bm{u}^{\prime}\in\mathit{unexplored}_{l}(\bm{u})$, uniformly at random; 2. (2)

select $\bm{u}^{\prime}\in\mathit{explored}_{l}(\bm{u})$, uniformly at random; 3. (3)

select $\bm{u}^{\prime}\in\mathit{explored}_{l}(\bm{u})$, uniformly at random from those that minimise $\rho(\varphi,\mathcal{M}(\bm{u}\bm{u}^{\prime}))$; 4. (4)

select $\bm{u}^{\prime}\in\mathit{explored}_{l}(\bm{u})$, uniformly at random from those that minimise $\rho(\varphi,\mathcal{M}(\bm{u}\bm{u}^{\prime}\bm{u}^{*}))$, where $\bm{u}^{*}$ denotes any, arbitrary length input signal suffix that has already been explored from $\bm{u}\bm{u}^{\prime}$.

Strategy 1 can be considered pure exploration, while strategies 3–4 are three different sorts of exploitation. In the case that $\mathit{unexplored}_{l}(\bm{u})$ or $\mathit{explored}_{l}(\bm{u})$ are empty, their corresponding strategies are infeasible and a strategy is chosen uniformly at random from the feasible strategies. If for all $\bm{u}^{\prime}\in\mathit{explored}(\bm{u})$, $\mathit{explored}(\bm{u}\bm{u}^{\prime})=\varnothing$, then strategy 4 is equivalent to strategy 3, but it is not infeasible.

3.4. Design Choices

In contrast to Monte Carlo tree search (MCTS), which typically makes a good decision by reaching a goal state multiple times during playout, our algorithm expects to reach a goal state only once. In place of a statistic about reaching goal states, we therefore adopt a heuristic, namely the lowest robustness. Using this heuristic is typically much more effective than choosing entirely at random, but simply substituting the heuristic value for the statistic in the UCT formula (upper confidence bound applied to trees (Kocsis and Szepesvári, 2006), successfully used with MCTS) is plausible, but nevertheless questionable in our application. Moreover, the deterministic way that UCT balances exploration and exploitation is often not optimal in the context of falsification.

Our design choices are guided by the observation that simple, coarse-grained input signals either immediately solve typical falsification problems, or need only small modifications to do so. Since such signals may be quickly explored exhaustively, our level scaling factors are designed to make this happen with high probability.

Unexplored edges have unknown potential, while every explored edge does not yet lead to falsification, but does not exclude it. The only exception to this occurs when all the edges of a node have been explored and discarded. In this unlikely case the edge leading to such a node exists and can be traversed, but the node and its trace are immediately rejected (algorithm 1).

While not explicit in the presentation of algorithm 1, our approach is deliberately incremental in the evaluation of the system model. In particular, we can re-use partial simulations to take advantage of the fact that traces share common prefixes. Hence, for example, one can associate to every visited $\bm{u}$ the terminal state of the simulation that reached it, using this state to initialize a new simulation when subsequently exploring $\bm{u}\bm{u}^{\prime}$. This idea also works for the calculation of robustness. We note, however, that incremental simulations may be impractical. For example, suspending and re-starting Simulink can be more expensive than performing an entire simulation from the start. Under such circumstances, the inner loop of our algorithm can be replaced by a recursion that executes the model on the complete input, returning the trace from the innermost recursion. The checks in algorithms 1 and 1 can be done on the appropriate prefix of the output returned from the recursive call (which takes place in algorithm 1).

4. Evaluation

We briefly describe our implementation and the benchmarks, then present the results of an experimental comparison with random sampling and Breach.

4.1. Implementation

We have implemented algorithm 1 in the prototype tool FalStar, using the Scala programming language and interfacing to MATLAB/Simulink through a Java API. The code is publicly available on github,111https://github.com/ERATOMMSD/falstar including all files and instructions required to re-run our experiments.

The main data structure is an explicit representation of the search tree. Its nodes are labelled by inputs and the scores relevant for strategies 3 and 4 of section 3.3. Similarly, we keep track of the sets $\mathit{unexplored}_{l}$ and $\mathit{explored}_{l}$ for each level $l$.

We do not re-run the simulation for every new input segment, as implied by algorithm 1 of algorithm 1. Instead, we cache the output traces within the tree data structure, alongside their inputs. As mentioned in section 3.4, the cost of pausing and resuming a simulation in Simulink incurs a prohibitively large overhead, hence we have implemented a recursive algorithm that executes exactly one simulation up to the full time horizon $T$ per trial/iteration of the outer loop.

The implementation supports initial parameters with respect to the input signals, such as engine speed $\omega$ that should remain constant, by augmenting the first level of the search tree with an additional dimension for each of the parameters.

Currently, the possible time durations $\mathbf{t}_{l}$ for each level $l$ can be specified statically as singleton sets, using a sequence of numbers of control points $(k_{0},\ldots,k_{l_{\mathrm{max}}})$, such that $\mathbf{t}_{l}=\{T/k_{l}\}$. The particular choices depend on the requirements of the case study, as noted below. We make use of both progressive subdivisions (increasing $k_{l}$ with $l$) and constant subdivision (keeping $k_{l}$ the same for all $l$).

We use the algorithm from (Donzé et al., 2013) to compute the robustness of STL formulas efficiently, and the approach from (Dreossi et al., 2015) to compute the respective upper and lower bounds.

FalStar exposes its functionality through a simple specification and scripting language. In the spirit of the SMT-LIB file format (Barrett et al., 2017), the language is S-expression-based, to facilitate integration with other tools. Using this interface one can declare models, input signal specifications, as well as STL requirements. In addition to solving falsification problems and generating reports, FalStar can be used to simulate systems for a given input or to compute robustness values of given STL requirements for specified output traces. To aid comparisons, FalStar can also call Breach as an alternative solver, or generate a standalone MATLAB script to run Breach.

4.2. Benchmarks

The benchmarks here are based on automotive Simulink models that are commonly used in the falsification literature. We focus on those benchmarks that have at least one time-varying input, in contrast to just searching for an initial configuration.

Automatic Transmission

This benchmark was proposed in (Hoxha et al., 2014) and consists of a Simulink model and its related requirements. The model has two inputs: throttle and brake. Depending on the speed and engine load, an appropriate gear is automatically selected. Besides the current gear $g$, the model outputs the speed $v$ and engine $\omega$. We consider the following requirements:

[TABLE]

There are a few noteworthy details. The robustness score of AT2 remains $\infty$ until $g=i$ is found by unguided search, and will jump abruptly to a finite value at that time. To falsify AT3, a fairly precise speed between 50 and 60 has to be reached and maintained, which requires a fine grained control of the input, which could be a challenge for aLVTS. To falsify AT4 one must balance keeping the engine RPM $\omega\leq\underline{\omega}$ while accelerating sufficiently to reach speed $v>\overline{v}$ before the deadline at time $10$.

The input signal for the benchmarks is piecewise constant with 4 control points for random sampling/Breach, which is sufficient to falsify all requirements. We choose 6 levels, with 2,2,3,3,3,4 control points, respectively, corresponding to a time granularity of input segment durations between 15 (= $\frac{30}{2}$, coarsest) to 7.5 (= $\frac{30}{4}$, finest).

Powertrain Control

The benchmark was proposed for hybrid systems verification and falsification in (Jin et al., 2014). It comprises a detailed model of a fuel control system that maintains the air-to-fuel $\mathit{AF}$ ratio close to a reference value $\mathit{AF}_{\text{ref}}$. The model has two control algorithms: During startup and when high power is requested, it operates in feed-forward mode, whereas otherwise measurements of the actual air-to-fuel are fed back into the controller.Falsification tries to detect amplitude and duration of spikes in the air-to-fuel ratio that occur as a response either to mode switches or to changes in the throttle $\theta$ at a given engine speed $\omega$. The relation between the three quantities is nontrivial and in part determined by time delays. The input $\theta\in[0,62.1)$ (normal mode) and $\theta\in[61.2,81.2)$ (power mode) varies throughout the trace, whereas $\omega\in[900,1100]$ is initially chosen and kept constant. Requirements are expressed with respect to the normalized error $\mu=|\mathit{AF}-\mathit{AF}_{\text{ref}}|/\mathit{AF}_{\text{ref}}$, where $\mathit{AF}_{\text{ref}}$ depends dynamically on the current mode.

We consider requirement 27 from (Jin et al., 2014), denoted AFC27 in the following, which states that after a falling or rising edge of the throttle, $\mu$ should return to the reference value within $1$ time unit and stay close to it for some time. We represent this as

[TABLE]

where rising and falling edges of $\theta$ are detected by $\mathit{rise}=\theta<8.8\land\Diamond_{[0,\epsilon]}\ 40.0<\theta$ and $\mathit{fall}=40.0<\theta\land\Diamond_{[0,\epsilon]}\ \theta<8.8$ for $\epsilon=0.1$. The concrete bound of 0.008 in AFC27 and the edge detection parameters are taken from the report of the ARCH friendly competition on falsification (which featured this problem in 2017 (Dokhanchi et al., 2017) and 2018 (Dokhanchi et al., 2019)) as a balance between difficulty of the problem and ability to find falsifying traces. The interval $[11,50]$ over which the requirement is checked begins $1$ time unit after the transition from startup to normal mode at time $10$. Without this short delay, falsification would trivially discover the large spike following the mode switch.

The input signal is piecewise constant with 10 control points, again following (Dokhanchi et al., 2017, 2019), and 5 levels with 10 control points each for aLVTS.

4.3. Experimental Results

We compare the performance of aLVTS with uniform random sampling (both implemented in FalStar222https://github.com/ERATOMMSD/falstar commit 43f5ca7) and with with the state-of-the-art stochastic global optimization algorithm CMA-ES (Igel et al., 2007) implemented in the falsification tool Breach.333https://github.com/decyphir/breach release version 1.2.9 We do not make a comparison with Breach’s Nelder-Mead algorithm, since it has significantly poorer performance than CMA-ES on almost all of the benchmarks. The machine and software configuration was: CPU Intel i7-3770, 3.40GHz, 8 cores, 8Gb RAM, 64-bit Ubuntu 16.04 kernel 4.4.0, MATLAB R2018a, Scala 2.12.6, Java 1.8.

We compare two performance metrics: success rate (how many falsification trials were successful in finding a falsifying input) and the number of iterations made, which corresponds to the number of simulations required and thus indicates time. To account for the stochastic nature of the algorithms, the experiments were repeated for $50$ trials. For a meaningful comparison of the number of iterations until falsification, we tried to maximize the falsification rate for a limit of $300$ iterations per trial.

The number of iterations of the top-level loop in algorithm 1 in our implementation corresponds exactly to one complete Simulink simulation up to the time horizon (cf. section 4.1). For random sampling and CMA-ES, the number of iterations likewise corresponds to samples taken by running exactly one simulations each. Hence the comparison is fair and, as the overhead is dominated by simulation time, the numbers are roughly proportional to wall-clock times.

Table 2 summarizes our results in terms of success rate, and average number of iterations (M) and standard deviation (SD) of successful trials. The unambiguously (possibly equal) best results are highlighted in blue. Where the lowest average number of iterations was achieved without finding a falsifying input for every trial, we highlight in grey the lowest average number of iterations for $100\%$ success. We thus observe that aLVTS achieves the best performance in all but one case, AT2 $(i=3)$. Importantly, within the budget of 300 iterations per trial, aLVTS achieves a perfect success rate. CMA-ES is successful in 296 trials out of the total 350, with sub maximal success for AT2 $(i=4)$ and AFC27. In comparison, random sampling succeeds in only 214 trials, with sub maximal success in all but AT2.

The number of iterations required for falsification varies significantly between the algorithms and between the benchmarks. For the automatic transmission benchmarks, as an approximate indication of relative performance, CMA-ES requires about 50% more iterations as aLVTS (geometric mean: 42.8 vs. 29.2), and random sampling requires again twice as many as CMA-ES. Except for AT2, variance of aLVTS is lower than that of CMA-ES.

Looking at the individual benchmarks, CMA-ES is only consistently faster than aLVTS on AT2 $(i=3)$. Despite having fewer average iterations on AT2 $(i=4)$, CMA-ES cannot be considered faster than aLVTS because most of its trials fail to find a falsifying input—the numbers of iterations are not comparable. For the powertrain model (AFC27), the performance of aLVTS is more than an order of magnitude better: 3.9 iterations on average, compared to 121 for CMA-ES.

Figure 2 compares all trial runs for AFC27, ordered by the number of iterations required for falsification. Similar plots for the automatic transmission benchmarks are shown in fig. 3. The shape of each curve gives an intuition of the performance and consistency of its corresponding algorithm. In general, fewer iterations and more successful results are better, so it is immediately clear that our aLVTS performs significantly better than both random sampling and CMA-ES.

4.4. Discussion

For AT1, aLVTS quickly finds the falsifying input signal, as the required throttle of $100$ and brake of [math] are contained in level 0 and are very likely to be tried early on. In contrast, even though this is a problem that is well-suited to hill-climbing, CMA-ES has some overhead to sample its initial population, as clearly visible in fig. 3(a).

While CMA-ES deals very well with AT2 for $i=3$, it struggles to find falsifying inputs for $i=4$ (cf. figs. 3(b) and 3(c)). We attribute this to the fact that reaching gear 4 by chance occurs rarely in the exploration of CMA-ES when the robustness score is uninformative. aLVTS not only explores the spatial dimensions, but takes opportunistic jumps to later time points, which increases the probability of discovering a trace (prefix) where the gear is reached.

A priori, one would expect CMA-ES to perform well with AT3 and AT4, exploiting its continuous optimization to fine tune inputs between conflicting requirements. E.g., AT3 requires that $v$ is both above $50$ and below $60$; AT4 requires that $v$ is high while maintaining low $\omega$, which is proportional to $v$. One would similarly expect the limited discrete choices made by aLVTS to hinder its ability to find falsifying inputs. Despite these expectations, our results demonstrate that in most situations aLVTS converges to a falsifying input more consistently and with fewer iterations than CMA-ES. We speculate that this is because CMA-ES is too slow to reach the “sweet spots” in the input space, where its optimization is efficient.

For AT3, there are a few instances where where aLVTS does not quickly find a good prefix towards the corridor $v\in[50,60]$ at time 10 (the rightmost points in fig. 3(d)), which can be explained by the probabilistic nature of the search.

Regarding two instances of AT4 in figs. 3(e) and 3(f), , the graph for aLVTS is generally smoother and shallower, whereas CMA-ES shows consistent performance for only some of the trials but takes significantly more time on the worst 10 trials. We remark that the two parameter settings seem to pose opposite difficulty for the two algorithms, as CMA-ES is significantly quicker for two thirds of the trials for the second instance. It is unclear what causes this variation in performance of CMA-ES.

The plateaux apparent in some of the results using CMA-ES are difficult to explain, but suggest some kind of procedural logic or counter to decide termination. In contrast, the curves for random sampling and aLVTS are relatively smooth, reflecting their purely probabilistic natures.

5. Related Work

The idea to find falsifying inputs using robustness as an optimization function for global optimization originates from (Fainekos and Pappas, 2009) and has since been extended to the parameter synthesis problem (e.g., (Jin et al., 2015)). Two mature implementations in MATLAB are S-Taliro (Annpureddy et al., 2011) and Breach (Donzé, 2010), which have come to define the benchmark in this field. We describe their optimization algorithms below. Approaches to make the robustness semantics more informative are (Eddeland et al., 2017; Akazaki and Hasuo, 2015). These use integrals instead of min/max in the semantics of temporal operators.

Underminer (Balkan et al., 2017) is a recent falsification tool that learns the (non-) convergence of a system to direct falsification and parameter mining. It supports STL formulas, SVMs, neural nets, and Lyapunov-like functions as classifiers. Other global approaches include (Adimoolam et al., 2017), which partitions the input space into sub-regions from which falsification trials are run selectively. This method uses coverage metrics to balance exploration and exploitation. Comprehensive surveys of simulation based methods for the analysis of hybrid systems are given in (Kapinski et al., 2016; Bartocci et al., 2018).

Users of S-Taliro and Breach can select from a range of optimization algorithms (Uniform Random, Nelder-Mead, Simulated Annealing, Cross-Entropy, CMA-ES). These cover a variety of trade-offs between exploration of the search space and exploitation of known good intermediate results. As such, their performance varies with the structure of the problem at hand.

In general, global optimization needs to find a solution in an unstructured combinatorial search space, where the different parameters are assumed, a priori, to be independent. Good combinations of parameter choices have to be discovered first, before the optimization can work effectively (cf. AT2 and AT3). This incurs an exponential cost of exhausting the different combinations, which is mitigated in different ways by the various algorithms. Nelder-Mead is seeded by a constant number of samples and tries to obtain better points in the search space by making linear combinations of the previous ones, according to a fixed scheme. The price paid is not just lack of exploration, but also the danger of selecting mid points between good and bad values in some dimensions, as a consequence of the weighted average of linear combinations. To this end, Cross-Entropy methods and CMA-ES discover relations between parameters and exploit these by decorrelating the choices across independent dimensions (i.e., not trying all possible combinations), as explained nicely in the context of falsification in (Sankaranarayanan and Fainekos, 2012). Simulated Annealing, on the other hand, relies on fairly unguided random exploration, and thus requires many iterations to converge, which can be prohibitively expensive for complex models.

In comparison, our approach explicitly takes into account the temporal structure of the falsification problem. The exponential search space is pruned by focusing on the early part of the input signal, by investigating a small number of coarse choices first, and by ignoring those prefixes which are unpromising according to the robustness heuristic. As demonstrated by our experiments, the strength of optimization algorithms to discover precise solutions in the continuous input space does not seem efficient when the mass of falsifying traces is large and exploration is the key factor.

The characteristic of our approach to explore the search space incrementally is shared with rapidly-exploring random trees (RRTs). The so-called star discrepancy metric guides the search towards unexplored regions and a local planner extends the tree at an existing node with a trajectory segment that closely reaches the target point. RRTs have been used successfully in robotics (LaValle and Kuffner Jr, 2001) and also in falsification (Dreossi et al., 2015). On the other hand, the characteristic of taking opportunistic coarse jumps in time is reminiscent of stochastic local search (Deshmukhand et al., [n. d.]) and multiple-shooting (Zutshi et al., 2014).

Monte Carlo tree search (MCTS) has been applied to a model of aircraft collisions in (Lee et al., 2015); and in a falsification context more recently (Zhang et al., 2018b) to guide global optimization, building on the previous idea of time-staging (Zhang et al., 2018a). That work noted the strong similarities between falsification using MCTS and statistical model checking (SMC) using importance splitting (Jegourel et al., 2013). The robustness semantics of STL, used in (Zhang et al., 2018a; Zhang et al., 2018b) and the present approach to guide exploration, can be seen as a “heuristic score function” (Jegourel et al., 2014) in the context of importance splitting. All these approaches construct trees from traces that share common prefixes deemed good according to some heuristic. The principal difference is that importance splitting aims to construct a diverse set of randomly-generated traces that all satisfy a property (equivalently, falsify a negated property), while falsification seeks a single falsifying input. The current work can be distinguished from standard MCTS and reinforcement learning (Sutton and Barto, 2018) for similar reasons. These techniques tend to seek optimal policies that make good decisions in all situations, unnecessarily covering the entire search space.

6. Conclusion

We have presented a probabilistic algorithm that finds inputs to a hybrid system that falsify a given temporal logic requirement. On standard benchmarks our algorithm significantly and consistently outperforms existing methods. While the falsification problem is inherently hard (no theoretically best solution can exist), we believe that we have found an approach that gives good results in practice, with probabilistic certainty. We believe the reason it works well stems from a property shared by many falsification problems, namely, that there tends to be a significant mass of falsifying inputs. By design it scales with the difficulty of the problem, finding trivial solutions (extreme inputs) immediately and has a small constant overhead for “almost easy” problems. This pays off, in particular, for system models that are expensive to simulate, such as the powertrain benchmark. The approach admits incremental computation of simulations and can easily be parallelized.

As future work, we wish to extend the method with some form of optimization, e.g., by interpolation of previously seen traces, reminiscent of the linear combinations computed by Nelder-Mead. We expect that this would help to “fine-tune” results with small positive robustness into full solutions, mitigating the limitation inherent in the discretization of the input space. Likewise, extrapolation of results could be used to propagate the robustness from one level to the next, reducing the need to look at different branches with similar inputs. Finally, we plan to apply the approach to more benchmarks, specifically those that combine discrete and continuous input domains.

Acknowledgements.

The authors are supported by ERATO HASUO Metamathematics for Systems Design Project (No. Grant #JPMJER1603), Sponsor Japan Science and Technology Agency http://www.jst.go.jp; and Grants-in-Aid Grant #No. 15KT0012, Sponsor Japan Society for the Promotion of Science https://www.jsps.go.jp.