Bounded Verification with On-the-Fly 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 FA9550-12-1-0336).
Simulation-based 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 piece-wise exponential discrepancy functions. The algorithm relies on computing local convergence or divergence rates of trajectories along a simulation using a coarse over-approximation of the reach set and bounding the maximal eigenvalue of the Jacobian over this over-approximation. 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 input-to-state 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.
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 simulation-based 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 over-approximation of the reachable states from the initial set. If this over-approximation proves safety or produces a counter-example, then the algorithm decides, otherwise, it draws more samples of initial states and repeats the earlier steps to compute more precise over-approximation. With post-processing of the reachtube over-approximations this basic procedure can be utilized to verify termporal precedence  and richer classes of properties .
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 : For an -dimensional dynamical system, it is any function , such that (a) it gives an upper-bound 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 .
In , it is observed that the notion of a contraction metric  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 , 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 simulation-based 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, one-step over-approximation 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 input-to-state 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.
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 vector-valued function , the Jacobian of is the matrix-valued function of all the first-order 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 negative-semidefinite, and write . We write if .
2.2 Safety Verification Problem
Consider an -dimensional autonomous dynamical system:
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 bounded-time 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 bounded-time 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.
A -simulation of the system described in Equation (1) is a sequence of time-stamped 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  and VNODE-LP  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 time-stamped 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.
A continuous function is a discrepancy function of the system in Equation (1) if
for any pair of states , and any time ,
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 , 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 simulation-based 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 simulation-based 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 bounded-time reachtube from .
Initially, has a singleton cover such that , , and is a small constant for simulation precision.
In the while-loop, this verification algorithm iteratively refines the cover of and for each in , computes over-approximations of the reachtube from . The higher-level 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 over-approximation of and the union of these bloated simulations is an over-approximation of . Therefore, the algorithm is sound. Furthermore, the second property of ensures that the reach set over-approximations 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 . In this paper, we will analyze in detail the subroutine which computes a local version of discrepancy function automatically.
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 time-varying 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.
Consider a compact set and a sequence of time points .
For , a piece-wise 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 piece-wise 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 for-loop it computes exponent corresponding to the time interval . In the iteration, is updated so that is an over-approximation 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) over-approximation 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 piece-wise 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 over-approximation of the set of states that can be reached by the system from over the time interval .
In iteration of the loop of , .
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 .
For any continuously differentiable vector-valued function , and ,
where the integral is component-wise.
Next, we will use a well-known theorem that gives bounds on eigenvalues of perturbed symmetric matrices, the proof of which uses the Courant-Fischer minimax theorem. The complete proofs of Lemma 3.3 and Theorem 3.4 can be found in the appendix.
If and are symmetric matrices, then
where is the largest eigenvalue of a matrix.
If and are symmetric matrices, then
Since is symmetric, .
From Theorem 3.4, we have .
If is a matrix-valued 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 .
In the iteration, for .
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 .
In the iteration, for any two states at time , and any time , .
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
For , , and , where is after Line 14 is executed in the iteration, and
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 for-loop, is set as .
When , we already have .
And at Line 14, . Since could be positive or negative,
and at time , is at most distance to , so .
The items in array computed by
are the coefficients of a -local piece-wise exponential discrepancy function (Definition 3.1).