Bounded Verification with OntheFly Discrepancy Computation^{†}^{†}thanks: We gratefully acknowledge the feedback from anonymous referees on a previous draft of this technical report. The results presented here came about from work supported and funded by the National Science Foundation (grant: CAR 1054247 and NSF CSR 1016791) and the Air Force Office of Scientific Research (AFOSR YIP FA95501210336).
Abstract
Simulationbased verification algorithms can provide formal safety guarantees for nonlinear and hybrid systems. The previous algorithms rely on user provided model annotations called discrepancy function, which are crucial for computing reachtubes from simulations. In this paper, we eliminate this requirement by presenting an algorithm for computing piecewise exponential discrepancy functions. The algorithm relies on computing local convergence or divergence rates of trajectories along a simulation using a coarse overapproximation of the reach set and bounding the maximal eigenvalue of the Jacobian over this overapproximation. The resulting discrepancy function preserves the soundness and the relative completeness of the verification algorithm. We also provide a coordinate transformation method to improve the local estimates for the convergence or divergence rates in practical examples. We extend the method to get the inputtostate discrepancy of nonlinear dynamical systems which can be used for compositional analysis. Our experiments show that the approach is effective in terms of running time for several benchmark problems, scales reasonably to larger dimensional systems, and compares favorably with respect to available tools for nonlinear models.
1 Introduction
Verifying and falsifying nonlinear, switched, and hybrid system models using numerical simulations have been studied in detail [10, 17, 4, 14, 9]. The bounded time safety verification problem for a given model is parameterized by a time bound, a set of initial states, and a set of unsafe states and it requires one to decide if there exists a behavior of the model that reaches any unsafe set from any initial state. The simulationbased procedure for this problem first generates a set of numerical approximations of the behaviors from a finite sampling of the initial states. Next, by bloating these simulations by an appropriately large factor it computes an overapproximation of the reachable states from the initial set. If this overapproximation proves safety or produces a counterexample, then the algorithm decides, otherwise, it draws more samples of initial states and repeats the earlier steps to compute more precise overapproximation. With postprocessing of the reachtube overapproximations this basic procedure can be utilized to verify termporal precedence [11] and richer classes of properties [7].
In order to make this type of procedure sound, the bloating factor should be chosen to be large. Specifically, it should be large enough to make each bloated simulation an over approximation of the reachable states of the system not only from the sampled initial state, but also from a large enough neighborhood of that state so that the union of these neighborhoods cover the entire set of initial states. On the other hand, to make the procedure complete, or at least relatively complete modulo the precision of the machine, it should be possible to make the error due to bloating arbitrarily small for any point in time. These two opposing requirements are captured in the definition of a discrepancy function of [10]: For an dimensional dynamical system, it is any function , such that (a) it gives an upperbound on the distance between any two trajectories and of the system , and (b) it vanishes as approaches . Simply using the Lipschitz constant of the dynamic function gives one such bound, but it grows exponentially with time even for some incrementally stable models [2].
In [10], it is observed that the notion of a contraction metric [19] gives a much tighter bound and it provided heuristics for finding them for some classes of polynomial systems. Sensitivity analysis approach gives strict error bounds for linear systems [9], but for nonlinear models the bounds are less clear. We present a more detailed overview of related work in Section 6. This paper fills this gap by providing a subroutine that computes a local version of the discrepancy function which turns out to be adequate and effective for sound and relatively complete simulationbased verification. This subroutine, , itself uses a Lipschitz constant and the Jacobian of the dynamic function (the right hand side of the differential equation) and simulations of the system. The Lispchitz constant is used to construct a coarse, onestep overapproximation of the reach set of the system along a simulation. Then it computes an upper bound on the maximum eigenvalue of the symmetric part of the Jacobian over this over approximation, using a theorem from matrix perturbation theory. This gives an exponential bound on the distance between two trajectories, but roughly, the exponent is the best it can be as it is close to the maximum eigenvalue of the linear approximation of the system in the neighborhood.
We propose two practical extensions of this approach. First, we show that a linear coordinate transformation can bring about exponential improvements in the estimated distance. Secondly, we propose a technique for computing inputtostate discrepancy functions for analyzing composed systems and systems with bounded nondeterministic inputs. We report the results from a number of experiments performed with a prototype implementation of this approach applied to safety verification.
2 Background
2.1 Notations
For a vector , is the norm of and denotes its component. For , For a set , . Let represents the Minkowski sum of and . Therefore, . For sets , is their convex hull. The diameter of a compact set is .
A continuous function is smooth if all its higher derivatives and partial derivatives exist and are also continuous. It has a Lipschitz constant if for every , . A function is a class function if it is continuous, strictly increasing, and .
We denote the transpose of a matrix by . The conjugated transpose of is the matrix obtained by replacing each entry in with its complex conjugate.
Given a differentiable vectorvalued function , the Jacobian of is the matrixvalued function of all the firstorder partial derivatives of . Let be the scalar components of . The Jacobian of is: . The symmetric part of the Jacobian of matrix is defined as .
For an matrix , represents the norm of : . If , , then we say is negativesemidefinite, and write . We write if .
2.2 Safety Verification Problem
Consider an dimensional autonomous dynamical system:
(1) 
where is a Lipschitz continuous function. A solution or a trajectory of the system is a function such that for any initial point and at any time , satisfies the differential equation (1).
The boundedtime safety verification problem is parameterized by: (a) an dimensional dynamical system, that is, the function defining the right hand side of its differential equation, (b) a compact set of initial states, (c) an open set of unsafe states, and (d) a time bound . A state in is reachable from within a time interval if there exists an initial state and a time such that . The set of all reachable states in the interval is denoted by . If then we write when set is clear from the context. Given a boundedtime safety verification problem, we would like to design algorithms for deciding if any reachable state is safe, that is, if . If there exists some such that , we say the system is robustly safe. A sequence of papers [10, 11, 9] presented algorithms for solving this problem for a broad class of nonlinear dynamical, switched, and hybrid systems. In the remainder of this section, we present an overview of this approach. (Figure 1).
2.3 Simulations, Reachtubes and Annotations
The algorithm uses simulation oracles that give sampled numerical simulations of the system from individual initial states.
Definition 2.1.
A simulation of the system described in Equation (1) is a sequence of timestamped sets , satisfying:

Each is a compact set in with .

The last time and for each , , where the parameter is called the sampling period.

For each , the trajectory from at is in , i.e., , and for any , the solution .
Simulation engines generate a sequence of states and error bounds using numerical integration. Libraries like CAPD [5] and VNODELP [20] compute such simulations for a wide range of nonlinear dynamical system models and the ’s are represented by some data structure like hyperrectangles.
Closely related to simulations are reachtubes. For a set of states , a reachtube of (1) is a sequence of timestamped sets satisfying:

Each is a compact set of states.

The last time and for each , .

For any , and any time , the solution .
A reachtube is analogous to a simulation from a set of states, but they are much harder to compute. In fact, an algorithm for computing exact reachtubes readily solves the safety verification problem.
The algorithms in [10, 16] require the user to decorate the model with annotations called discrepancy functions for computing reachtubes.
Definition 2.2.
A continuous function is a discrepancy function of the system in Equation (1) if

for any pair of states , and any time ,
(2) 
for any , as , ,
If the function meets the two conditions for any pair of states in a compact set then it is called a local discrepancy function.
The annotation gives an upper bound on the distance between two neighboring trajectories as a function of their initial states and time. Unlike incremental stability conditions [2], the second condition on does not require the trajectories to converge as time goes to infinity, but only as the initial states converge. Obviously, if the function has a Lipschitz constant , then meets the above criteria. In [10, 16] other heuristics have been proposed for finding discrepancy functions. As will be clear from the following discussion, the quality of the discrepancy function strongly influences the performance of the simulationbased verification algorithm. [10, 16, 17] need user provided discrepancy and simulation engines to give verification of bounded time safety and temporal precedence properties. In this paper, we will present approaches for computing local discrepancy functions that unburdens the user from finding these annotations.
2.4 Verification Algorithm
The simulationbased verification algorithm is shown in Figure 1. It takes as input some finite description of the parameters of a safety verification problem, namely, the function , the initial set , the unsafe set , and the time bound . It has two main data stuctures: The first, returned by function , is a collection of triples such that the union of all the balls around the ’s completely cover the initial set . The second data structure incrementally gets the boundedtime reachtube from .
Initially, has a singleton cover such that , , and is a small constant for simulation precision.
In the whileloop, this verification algorithm iteratively refines the cover of and for each in , computes overapproximations of the reachtube from . The higherlevel structure of the algorithm is familiar: if the reachtube from proves to be safe, i.e., disjoint from , then the corresponding triple is removed from (Line 10). If part of the reachtube from overlaps with , then the system is declared to be unsafe (Line 12). Otherwise, a finer cover of is created, and the corresponding triples with finer parameters are added to .
Here we discuss the reachtubes computed from discrepancy and simulations. For each in , a simulation , which is a sequence of , is generated. Note that contains the trajectory from , . Then we bloat each by some factor (Line 8) such that the resulting sequence contains the reachtube from . It is shown that this bloated simulation is guaranteed to be an overapproximation of and the union of these bloated simulations is an overapproximation of . Therefore, the algorithm is sound. Furthermore, the second property of ensures that the reach set overapproximations become tighter and tighter as we make smaller and smaller. Finally it will return “SAFE” for robustly safe reachtubes or find a counter example and return “UNSAFE”. For user defined discrepancy function, the factor is obtained by maximizing over and .
Indeed this is the approach taken in the algorithm presented in [10]. In this paper, we will analyze in detail the subroutine which computes a local version of discrepancy function automatically.
The following results from [10] state two key properties of the algorithm. Although in [10] was defined globally, it is easy to check that the local version still satisfies them.
Theorem 2.3.
The Algorithm in Fig.1 is sound, that is, if it returns “SAFE” then the system is safe; when it returns “UNSAFE” there exists at least one execution from that is unsafe. The Algorithm is relatively complete, that is, if the system is robustly safe, the algorithm will terminate and return “SAFE”. If any executions from is unsafe, it will terminate and return “UNSAFE”.
3 Local discrepancy function
In this section, we present the analysis of algorithm. This algorithm computes a special type of local discrepancy in terms of timevarying exponential functions that bound from above the distance between two trajectories starting from a compact neighborhood. Roughly speaking, it computes the rate of trajectory convergence or divergence for an interval of time instances.
Definition 3.1.
Consider a compact set and a sequence of time points .
For , a piecewise exponential discrepancy function
is defined as
where are real constants.
From the definition, we can immediately get that , , where is the largest time point in the sequence before .
3.1 ComputeLDF Algorithm
Figure 2 shows the pseudocode for used in Line 7 of the verification algorithm. takes as input a parameter , an error bound for simulation , the Lipschitz constant , the Jacobian matrix of function , and a simulation . It computes a piecewise exponential local discrepancy function (LDF) for the compact set and for the time points . and returns it as an array of exponential coefficients .
The algorithm starts with the initial set and with . In each iteration of the forloop it computes exponent corresponding to the time interval . In the iteration, is updated so that is an overapproximation of the reachable states from at (Lemma 3.8). In Lines 8 and 9, a set is computed by bloating the convex hull by a factor of . The set will later be proved to be a (coarse) overapproximation of the reachtube from over the time interval (Lemma 3.2). In Lines 10–13 an upper bound on the maximum eigenvalue of the symmetric part of the Jacobian over the set , is computed as (Lemma 3.6). Then is updated as for the next iteration.
3.2 Analysis of ComputeLDF
In this section, we will prove that returns a piecewise exponential LDF of the system in Equation (1), for the compact neighborhood , and the sequence of the time points in the simulation . We establish some lemmas to prove the main theorem. First, we show in Lemma 3.2 that in the iteration of the loop, the computed is an overapproximation of the set of states that can be reached by the system from over the time interval .
Lemma 3.2.
In iteration of the loop of , .
Proof.
Let denote the actual trajectory from , where is the initial state of . By Definition 2.1 for , it is known that and .
For a fixed iteration number , consider state from Definition 2.1. We know that for any , . Now consider another state . Since is the Lipschitz constant of , using Gronwall’s inequality we have that . Since , . Therefore, . Because is arbitrarily selected from , the lemma is proved.
Next, using the generalized mean value theorem (Lemma 3.3), we get that in the iteration, the computed in Line 13 is the exponential divergence (if positive) or convergence (negative) rate of the distance between any two trajectories starting from over time .
Lemma 3.3.
For any continuously differentiable vectorvalued function , and ,
(3) 
where the integral is componentwise.
Next, we will use a wellknown theorem that gives bounds on eigenvalues of perturbed symmetric matrices, the proof of which uses the CourantFischer minimax theorem. The complete proofs of Lemma 3.3 and Theorem 3.4 can be found in the appendix.
Theorem 3.4.
If and are symmetric matrices, then
where is the largest eigenvalue of a matrix.
Corollary 3.5.
If and are symmetric matrices, then
(4) 
Since is symmetric, .
From Theorem 3.4, we have .
If is a matrixvalued function: maps a state to a matrix , and every component of is continuous over some compact closed set , then we can over by each term , over . Let be denoted by , then we know
.
Using Corollary 3.5, we next show in Lemma 3.6 that calculated in Line 13 bounds the eigenvalues of symmetric part of Jacobian matrix over .
Lemma 3.6.
In the iteration, for .
Proof.
Let be the set computed in Line 9 and be the Jacobian evaluated at the center of . Consider any point . We define the perturbation matrix . Since and are symmetric matrices, Corollary 3.5 implies that . The term computed in Line 12 is the upperbound on . Therefore, . In Line 13 set equals to . Thus, , which immediately indicates that .
By Lemma 3.3 and Lemma 3.6, we can prove as in Lemma 3.7 that calculated in Line 13 is the exponential rate of divergence or convergence of two trajectories starting from over the interval .
Lemma 3.7.
In the iteration, for any two states at time , and any time , .
Proof.
Let us fix the iteration and two states . From Lemma 3.2 it’s can be seen that for any , . Define . For a fixed time , from Lemma 3.3 we have
(5)  
Since is the Minkowski sum of two convex sets and , it is also convex. Recall that , and for any , .
Differentiating , we have
Integrating both sides over to any , we have
Up to this point all the lemmas were statements about a single iteration of the for loop, next we show that in iteration of the loop, used in Lemma 3.2 and 3.7 is the reach set from at time .
Lemma 3.8.
For , , and , where is after Line 14 is executed in the iteration, and
Proof.
In this proof, let denote the trajectory from . From the Definition 2.1 for , we know that and . Let denote after Line 9 is executed in the iteration. The lemma is proved by induction on . Note that the initial set is , and before the forloop, is set as .
When , we already have .
contains the reachtube of the system. Line 8 of the algorithm in Figure 1 is computed in this way. Now we are ready to prove the main theorem.
Theorem 3.9.
The items in array computed by
are the coefficients of a
local piecewise exponential discrepancy function (Definition 3.1).