Numerically-Robust Inductive Proof Rules for Continuous Dynamical Systems

. We formulate numerically-robust inductive proof rules for unbounded stability and safety properties of continuous dynamical systems. These induction rules robustify standard notions of Lyapunov functions and barrier certiﬁcates so that they can tolerate small numerical errors. In this way, numerically-driven decision procedures can establish a sound and relative-complete proof system for unbounded properties of very general nonlinear systems. We demonstrate the eﬀectiveness of the proposed rules for rigorously verifying unbounded properties of various nonlinear systems, including a challenging powertrain control model.


Introduction
Infinite-time stability and safety properties of continuous dynamical systems are typically established via inductive arguments over continuous time. For instance, proving stability of a dynamical system is similar to proving termination of a program. A system is stable at the origin in the sense of Lyapunov, if one can find a Lyapunov function (essentially a ranking function) that is everywhere positive except for reaching exactly zero at the origin, and never increases over time along the direction of the system dynamics [11]. Likewise, proving unbounded safety of a dynamical system requires one to find a barrier function (or differential invariant [19]) that separates the system's initial state from the unsafe regions, and whenever the system states reach the barrier, the system dynamics always points towards the safe side of the barrier [21]. In both cases, once a candidate certificate (Lyapunov or barrier functions) is proposed, the verification problem is reduced to checking the validity of a universally-quantified first-order formula over real-valued variables. The standard approaches for the validation step use symbolic quantifier elimination [4] or Sum-of-Squares techniques [18,17,24]. However, these algorithms are either extremely expensive or numerically brittle. Most importantly, they can not handle systems with nonpolynomial nonlinearity, and thus fall short of a general framework for verifying practical systems of significant complexity.
The standard approach of checking invariance conditions in program analysis is to use Satisfiability Modulo Theories (SMT) solvers [16]. However, to check the inductive conditions for nonlinear dynamical systems, one has to solve nonlinear SMT problems over real numbers, which are highly intractable or undecidable [23]. Recent work on numerically-driven decision procedures provides a promising direction to bypass this difficulty [5,6]. They have been used for many bounded-time verification and synthesis problems for highly nonlinear systems [12]. However, the fundamental challenge with using numerically-driven methods in inductive proofs is that numerical errors make it impossible to verify the induction steps in the standard sense. Take the Lyapunov analysis of stability properties as an example. A dynamical system is stable if there exists a function that vanishes exactly at the origin and its derivatives strictly decreases over time.
Since any numerical error blurs the difference between strict and non-strict inequality, one can conclude that numerically-driven methods are not suitable for verifying these strict constraints. However, proving a system is stable within an arbitrarily tiny neighborhood around the origin is all we really need in practice. Thus, there is a discrepancy between what the standard theory requires and what is needed in practice, or what can be achieved computationally. To bridge this gap, we need to rethink about the fundamental definitions.
In this paper, we formulate new inductive proof rules for continuous dynamical systems for establishing robust notions of stability and safety. These proof rules are practically useful and computationally certifiable in a very general sense. For instance, for stability, we define the notion of ε-stability that requires the system to be stable within an ε-bounded distance from the origin, instead of exactly at the origin. When ε is small enough, ε-stable systems are practically indistinguishable from stable systems. We then define the notion of ε-Lyapunov functions that are sufficient for establishing ε-stability. We then rigorously prove that the ε-Lyapunov conditions are numerically stable and can be correctly determined by δ-complete decisions procedures for nonlinear real arithmetic [7]. In this way, we can rely on various numerically-driven SMT solvers to establish a sound and relative-complete proof systems for unbounded stability and safety properties of highly nonlinear dynamical systems. We believe these new definitions have eliminated the core difficulty for reasoning about infinite-time properties of nonlinear systems, and will pave the way for adapting a wide range of automated methods from program analysis to continuous and hybrid systems. In short, the paper makes the following contributions: -We define ε-stability and ε-Lyapunov functions in Section 3. We prove that finding ε-Lyapunov functions is sufficient for establishing ε-stability. -We define two types of robust proof rules for unbounded safety in Section 3, which we call Type 1 and Type 2 ε-barrier functions. The former relies on strict contraction, and the latter relies on reachable-set computation to guarantee bounded escape.
-We prove that δ-complete decision procedures provide a sound and relativecomplete proof system for the proposed numerically-robust induction rules, in both Sections 3 and 4. We demonstrate the effectiveness of the proposed methods on various nonlinear systems in Section 5. Section 2 covers the basic definitions and Section 6 concludes the paper.
Related Work. Several lines of work have proposed relaxed and practical notions to capture the spirit of the stability requirements. Early work from the 1960s introduced practical stability, which defined bounds on system behaviors over finite time horizons [14,26,27,2]. These methods can show whether a system leaves a safe set or enters a goal set over a finite time horizon based on Lyapunov-like functions. Stability defined in this sense is equivalent to estimating the reachable set over a finite time horizon. Thus, the shortcoming is that it may not capture the desired behavior of the system over unbounded time. Similarly, notions of boundedness and ultimate boundedness specify limits on the system behaviors [11]. Boundedness specifies whether the system remains within a given bounded region. Ultimate boundedness specifies that the system eventually returns to the given bounded region. These properties can be established based on Lyapunov-like conditions. Related notions have been generalized to switched systems [29,30]. Also, the related notion of region stability defines systems that eventually enter and remain within a specified set [20]. We present stability concepts that unify and extend the above notions. A related relaxation of the traditional notions of stability includes almost Lyapunov functions [15], which allow the strict stability conditions to be neglected in a region near the equilibrium point. The challenge of applying this technique in practice is that the size and shape of the neglected region are not specified a priori, so a constructive technique for specifying a stability region is not straightforward. Our work is related to efforts to construct and check robust barrier certificates using Lyapunov-like functions to ensure that controllers satisfy safety constraints [28]. This work provides a framework in which to specify analytic constraints on controller behaviors. By contrast, our work focuses on providing constraints that can be checked fully automatically. Our notion of ε-barrier functions is closely related to t-barrier certificates from [1], though we choose to focus on distance bounds from the barrier (ε) rather than time bounds that indicate how long it takes for behaviors to re-enter the barrier once it has left (t).

Dynamical Systems
Throughout the paper, we use the following definition of an n-dimensional autonomous dynamical system: where an open set D ⊆ R n is the state space, init ⊆ D is a set of initial states, and f : D → R n is a vector field specified by Lipschitz-continuous functions on each dimension. For notational simplicity, all variable and function symbols can represent vectors. When vectors are used in logic formulas, they represent conjunctions of the formulas for each dimension. For instance, when x = (x 1 , . . . , x n ), we write x = 0 to denote the formula x 1 = 0 ∧ · · · ∧ x n = 0. For any system defined by (1), we write its solution function as Note that F usually does not have an analytic form. However, since f is Lipschitzcontinuous, F exists and is unique. We will often use Lie derivatives to measure the change of a scalar function along the flow defined by another vector field:

First-order Language over the Reals L RF
We will make extensive use of first-order formulas over real numbers with Type 2 computable functions [25] to express and infer properties of nonlinear dynamical systems. Definition 2 introduces the syntax of these formulas.
Definition 2 (Syntax of L R F ). Let F be the class of all Type 2 computable functions over real numbers. We define: ), where f ∈ F, possibly constant; We regard ¬ϕ as an operation that is defined inductively as usual. For instance, ¬(t > 0) is defined as −t ≥ 0, and ¬(∃x i ϕ) is defined as ∀x i ¬ϕ. For any L R F terms u and v, variable x, and L R F predicate ϕ, we write ∃ [u,v] xϕ and ∀ [u,v] respectively, which applies to open intervals too. Next, Definition 3 introduces syntactic perturbation of formulas in L R F . Definition 3 (δ-Strengthening and Robust Formulas [7]). Let δ ∈ Q + be arbitrary. Let ϕ be an arbitrary L R F formula. The δ-strengthening of ϕ, denoted by ϕ +δ , is obtained from ϕ by replacing every atomic predicate of the form Definition 4 (δ-Complete Decision Procedures [7]). Let S be a class of L R F -sentences. We say a decision procedure is δ-complete over S iff for any ϕ ∈ S, the procedure correctly returns one of the following answers: true : ϕ is true.
δ-false : ϕ +δ is false. When the two cases overlap, either decision can be returned.
It follows that if ϕ is δ-robust, then a δ-complete decision procedure can correctly determine the truth value of ϕ.

Robust Proofs for Stability
We first focus on stability. We will define the notion of ε-stability, as a relaxation of the standard Lyapunov stability, and then define ε-Lyapunov functions, which are sufficient for proving ε-stability in a robust way.

Stability and Lyapunov Functions
Conventionally, ε and δ are used to best highlight the connection with ε-δ conditions for continuity. We will mostly reserve the use of ε for defining conditions that are robust under ε-bounded numerical errors. Thus, we replace ε by τ in the standard definitions to avoid confusion.
Definition 5 (Stability). We say the system in (1) is stable at the origin in the sense of Lyapunov, iff for any τ -ball neighborhood of the origin, there exists a δ-ball around the origin, such that, if the system starts within the δ-ball then it never escapes the τ -ball. We capture the definition by the following L R F -formula: Definition 6 (Lyapunov Function). Consider a dynamical system given in the form of (1), and let V : D → R be a differentiable function. We say V is a non-strict Lyapunov function for the system, iff the following predicate is true: Proposition 1. For any dynamical system defined by f , if there exists a Lyapunov function V , then the system is stable. Namely, LF(f, V ) → Stable(f ).

Epsilon-Stability
The standard definitions of stability requires a system to stabilize within arbitrarily small neighborhoods around the origin. However, very small neighborhoods are practically indistinguishable from the origin. Thus, it is practically sufficient to prove that a system is stable within some sufficiently small neighborhood. We capture this intuition by making a minor change to the standard definition, by simply putting a lower bound ε on the τ parameter in Definition 5. As a result, the system is required to exhibit the same behavior as standard stable systems outside the ε-ball, but can behave arbitrarily within the ε-ball (for instance, oscillate around the origin). The formal definition is as follows: Definition 7 (Epsilon-Stability). Let ε ∈ R + be arbitrary. We say a dynamical system in (1) is ε-stable at the origin in the sense of Lyapunov, iff it satisfies the following condition: In words, for any τ ≥ ε, there exists δ such that all trajectories that start within the δ-ball will stay within a τ -ball around the origin.
Note that the only difference with the standard definition is that τ is bounded from below by a positive ε instead of 0. The definition is depicted in Figure 1c, which shows the difference with the standard notion in Figure 1a. Since the only difference with the standard definition is the lower bound on the universally quantified τ , it is clear that ε-stability is strictly weaker than standard stability.
Thus, any system that is stable in the standard definition is also ε-stable for any ε ∈ R + . On the other hand, one can always choose small enough ε such that an ε-stable system is practically indistinguishable from stable systems in the standard definition.

Epsilon-Lyapunov Function
We now define the corresponding notion of Lyapunov function that can be used for proving ε-stability. The robustness problem in the standard definition comes from the singularity of the origin. With the relaxed notion of stability, the system may oscillate within some ε-neighborhood of the origin. With the relaxation, we now have room for constructing a few nested neighborhoods that can trap the trajectories in a way that is robust under sufficiently small perturbations. To achieve this, we make use of balls of different sizes, as shown in the following definition. We write B ε to denote open ε-balls around the origin.
Definition 8 (Epsilon-Lyapunov Functions). Let V : D → R be a differentiable scalar function defined for the system in (1), and let ε ∈ R + be an arbitrary value. We say V is an ε-Lyapunov function for the system, iff it satisfies the following conditions: 1. Outside the ε-ball, there is some positive lower bound on the value of V .
Namely, there exists α ∈ R + such that for any x ∈ D \ B ε , V (x) ≥ α. 2. Inside the ε-ball, there is a strictly smaller ε -ball in which the value of V is bounded from above, to create a gap with its values outside the ε-ball.
Formally, there exists ε ∈ (0, ε) and β ∈ (0, α) such that for all In sum, the three conditions can be expressed with the following L R F -formula: It is important to note that ε , α, β, and γ, are not fixed constants, but existentially quantified variables. Thus the condition can hold true for infinitely many values of these parameters, which is critical to robustness. The only free variable in the formula is ε, used in B ε and the bound for ε . Note also that neither of LF ε (f, V ) and the standard definition LF(f, V ) implies the other.
is seemingly more complex than the standard Lyapunov conditions in Definition 6 because of the extra existential quantification. In Theorem 3, we show that it does not add computational complexity in checking the conditions.
The key result is that the conditions for an ε-Lyapunov function are sufficient for establishing ε-stability. Theorem 1. If there exists an ε-Lyapunov function V for a dynamical system defined by f , then the system is ε-stable. Namely, LF ε (f, V ) → Stable ε (f ).
Proof. Let τ ≥ ε be arbitrary, and let α, γ ∈ R + , β ∈ (0, α), and ε ∈ (0, ε) be as specified by the definition of LF ε (f, V ). Let x 0 ∈ B ε be an arbitrary point. For any t ∈ R ≥0 , let x(t) := F (x 0 , t) be the system state as defined in (2). We use contradiction to prove for any t ∈ R + , inequality x(t) < ε ≤ τ holds. Since ε < ε and F (x 0 , .) is continuous, we know t 1 and t 2 with the following conditions exists (∂B ε and ∂B ε are boundaries of the corresponding balls): We know V (x(t 1 )) ≤ β < α ≤ V (x(t 2 )) and hence V (x(t 1 )) < V (x(t 2 )) are both true; however, this is in contradiction with the mean value theorem and the fact that Remark 2. Proof of Theorem 1 shows that once state of the system enters B ε , it never leaves B ε . However, it would be still possible for the state to leave B ε . One the other hand, since closure of B ε \ B ε is bounded, and for every x in this area, V is continuous at x and ∇ f V (x) ≤ −γ, no trajectory can be trapped in the closure of B ε \ B ε . Therefore, even though state of the system might leave B ε , it will visit inside of this ball infinitely often. Example 1. Consider the time-reversed Van der Pol system given by the following dynamics. Figure 3 shows the vector field of this system around the origin.
, and P is the 9 × 9 constant matrix given in [8], is a 6-degree polynomial that can be obtained using simulation-guided techniques from [10]. Using dReal [9] with δ := 10 −25 and the Euclidean norm, we are able to prove that z T P z is a 10 −12 -Lyapunov function. Table 1 lists the parameters used for this proof.

Automated Proofs with Delta-Decisions
We now prove that unlike the conventional conditions, the new inductive proof rules are numerically robust. It follows that δ-decision procedures provide a sound and relative-complete proof system for establishing the conditions in the following sense: -(Soundness) A δ-complete decision procedure is always correct when it confirms the existence of an ε-Lyapunov function. -(Relative Completeness) For a given ε-inductive certificate, there exists δ > 0 such that a δ -complete procedure is able to verify it, for any 0 < δ ≤ δ. To prove these properties, the key fact is that the continuity of the functions in the induction conditions ensures that there is room for numerical errors in the conditions. Consequently, the formulas allow δ-perturbations in their parameters. This is captured by Lemma 1, and the proof is given in [8].
Note that if a formula φ is δ-robust then for every δ ∈ (0, δ), φ is δ -robust as well. The soundness and relative-completeness then follow naturally.
Theorem 2 (Soundness). If a δ-complete decision procedure confirms that LF ε (f, V ) is true then V is indeed an ε-Lyapunov function, and f is ε-stable.
Theorem 3 (Relative Completeness). For any ε ∈ R + , if LF ε (f, V ) is true then there exists δ ∈ Q + such that any δ-complete decision procedure must return that LF ε (f, V ) is true.
We remark that the quantifier alternation used in Definition 8 can be eliminated without extra search steps. It confirms that we only need to run SMT solving to handle the universally quantified subformula. The reason is that the α, β, and γ parameters can be found by estimating the range of V (x) and ∇ f V (x) in the different neighborhoods. In fact, we can rewrite LF ε (f, V ) in the following way to eliminate the use of α, β, and γ: Note that in this form the universal quantification is implicit in the sup and inf operators. In this way, the formula is existentially quantified on only ε , which can then be handled by binary search. This is an efficient way of checking the conditions in practice. We also remark that without this method, the original formulation with multiple parameters can be directly solved as ∃∀-formulas as well using more expensive algorithms [13].

Robust Proofs for Safety
In this section, we define two types of ε-barrier functions that are robust to numerical perturbations. Proving unbounded safety requires the use of barrier functions. The idea is that if one can find a barrier function that separates initial conditions from the set of unsafe states, such that no trajectories can cross the barrier from the safe to the unsafe side, then the system is safe. Here we use a formulation similar to the that of Prajna [21]. The standard conditions on barrier functions include constraints on the vector field of the system at the exact boundary of the barrier set, which introduces robustness problems. We show that it is possible to avoid these problems using two different formulations, which we call Type 1 and Type 2 ε-barrier functions. Type 1 ε-barrier functions strengthen the original definition and requires strict contraction of the barrier. Instead of only asking the system to be contractive exactly on the barrier's border, we force it to be contractive when reaching any state within a small distance from the border. Type 2 ε-barrier functions allow the system to escape the barrier for a controllable distance and a limited period of time. It should then return to the interior of the safe region. Type 1 ε-barriers can be seen as a subclass of Type 2 ε-barriers. The benefit for allowing bounded escape is that the shape of the barrier no longer needs to be an invariant set, which can be particularly helpful when the shape of the system invariants cannot be determined or expressed symbolically. The downside to Type 2 ε-barriers is that checking the corresponding conditions requires integration of the dynamics, which can be expensive but can still be handled by δ-complete decision procedures. The intuition behind the two definitions is shown in Figure 2 and will be explained in detail in this section.

Safety and Barrier Functions
Before formally introducing robust safety and ε-barrier functions, we define the safety and barrier functions first. It is easy to see that the robustness problem with the barrier functions is similar to that of Lyapunov functions: if the boundary is exactly separating the safe and unsafe regions then the inductive conditions are not robust, since deviations in the variables by even a small amount from the barrier will make it impossible to complete the proof.
Definition 9 (Safety). Let B : D → R be a scalar function defined for the system in (1). We say B ≤ 0 defines a safe (or forward invariant) set for the system, iff the following formula is true: Definition 10 (Barrier Function). Let B : X → R be a differentiable scalar function defined for the system in (1). We say B is a barrier function for the system, iff the following formula is true:

Type 1: Strict Contraction
In the standard definition, the boundary of the barrier set is typically a manifold defined by equality, which is not numerically robust. To avoid this problem, we need the barrier boundary to be belt-shaped in the sense that there is a clear gap between the safe and unsafe regions. The idea is as shown in Figure 2c: we need a second and stronger barrier defined by B = −ε for some reasonable ε, so that the system is clearly separated from B = 0. The formal definition is as follows.
Definition 11 (ε-Barrier Certificates). Let ε ∈ R + be arbitrary. A differentiable scalar function B : D → R is an ε-barrier function iff the following conditions are true: Formally, the condition is defined as It should be intuitively clear from the definition that the existence of ε-barrier functions is sufficient for establishing invariants and safety properties. The new requirement is that the system stays robustly within the barrier, by the area defined by −ε ≤ B(x) ≤ 0. Proof. Assume Barrier ε (f, init, B) is true. It is easy to see Barrier(f, init, B +ε), as specified in Definition 10, is also true. Therefore, using Proposition 3, we know Safe(f, init, B + ) and hence Safe(f, init, B) are both true.
It is clear that there is room for numerically perturbing the size of the area and still obtaining a robust proof. The proof is similar to the one for Lemma 1 as shown in [8].
Example 2 (Type 1 ε-Barrier for timed-reversed Van der Pol). Consider the time-reversed Van der Pol system introduced in Example 1. We use the same example to demonstrate the effect of numerical errors in proving barrier certificates. The level sets of the Lyapunov functions in the stable region are barrier certificates; however, for the barriers that are very close to the limiting cycle, numerical sensitivity becomes a problem. In experiments, when ε = 10 −5 and δ = 10 −4 , we can verify that the level set z T P z = 90, is a Type 1 ε-barrier. Table 2 lists parameters used in this proof. Figure 3 (Left) shows the direction field for the timed-reversed Van der Pol dynamics, the border of the set z T P z ≤ 90, which we prove is a type 1 ε-barrier, and the boundary of set z T P z ≤ 110, which is clearly not a barrier, since it is outside of the limit cycle.
The conditions for ε-Lyapunov and ε-barrier functions look very similar, but there is an important difference. In the case of Lyapunov functions, we do not evaluate the Lie derivative of the balls. Thus, the balls do not define barrier sets. On the other hand, the level sets of Lyapunov functions always define barriers.
Remark 3. The ε-barrier functions can also be used as a sufficient condition for ε-stability, if a barrier can be found within the ε-ball required in ε-stability.  Remark 4. A technical requirement for proving robustness of the ε-barrier conditions is that ¬init defines a simple set that can be over-approximated, such that for every ε ∈ R + , there is δ ∈ R + such that for any point that satisfies ¬init +δ there is an ε-close point that satisfies ¬init. A sufficient condition for this restriction is that init be of the form where a i , b i ∈ Q are arbitrary constants, and ϕ is a quantifier-free formula with only strict inequalities [22].

Type 2: Bounded Escape
We now introduce the second set of conditions for establishing ε-invariant sets. This set of conditions can be used only when the ε-variations are considered. This notion is inspired by the notion of k-step invariants [3] for discrete-time systems.
The ε-margin that we allow at the boundary of the invariants allows us to exploit more techniques. Using reachable set computation, we can directly check if all states stay within the barrier set at each step. To ensure that the conditions are inductive and useful, we need to impose the following two requirements: -(Contraction) Similar to the strengthening in barrier certificates, we require that the system does not sit at the boundary: the dynamics at the boundary should be contracting. The difference with Type 1 ε-barriers is that, this condition is not imposed through the vector field on the boundary. Instead, it is a reachability condition: after some amount of time, all states should return to the interior of an appropriate set. -(Bounded Escape) Before reaching back to the invariant set, we allow the system to step outside the invariant, but only up to a bounded distance from the boundary. The intuition is depicted in Figure 2d. In the formal definition, we parameterize the conditions with the time for contraction and the maximum deviation from the invariant set, as follows.
2. There exists ε > ε such that any state in B(x) ≤ −ε will enter B(x) ≤ −ε after time T . 3. During time [0, T ], the system may step outside of B(x) ≤ −ε but there exists some ε * ∈ (0, ε] such that all states stay within B(x) ≤ −ε * . In all, we define the conditions with the following formula Theorem 6, shows that conditions in Definition 12 ensure that the system never leaves the invariant B ≤ 0. The key is the second condition: induction works because all states come back to the interior of the set defined by B ≤ −ε.
With the third condition only, we cannot perform induction because the set may keep growing.
Proof. For the purpose of contradiction, suppose starting from x 0 ∈ init, the system is unsafe. Using continuity of the barrier B and the solution function F , let t ∈ R ≥0 be a time at which B(x(t)) = 0, where x(t) is by definition F (x 0 , t). By the 1 st property in Definition 12, we know B(x 0 ) ≤ −ε < 0. Using continuity of B and F , let t ∈ [0, t) be the supremum of all times at which B(x(t )) = −ε. By the 3 rd property in Definition 12, we know t−t > T , and by the 2 nd property in Definition 12, we know B(x(t + T )) ≤ −ε < −ε. Using continuity of B and F , we know there is a time t ∈ (t + T, t) at which B(x(t )) = −ε. However, this is in contradiction with t being the supremum.
Theorem 7. For any ε ∈ R + , there exists δ ∈ Q + such that Barrier T,ε (f, init, B) is a δ-robust formula. Example 3. We use this example to show how Type 2 ε-barriers can be used to establish safety. Consider the following system.
Let init be the set {x | −0.1 ≤ x 1 ≤ 0.1, −0.1 ≤ x 2 ≤ 0.1}, and let U , the unsafe set, be the set {x | −2.0 ≤ x 1 ≤ −1.1, −2.0 ≤ x 2 ≤ −1.1}. The system is stable and safe with respect to the designated unsafe set. However, the safety cannot be shown using any invariant of the form B(x) := x 2 1 + x 2 2 − c ≤ 0, where c ∈ Q + is a constant, in the standard definition. This is because the vector field on the boundary of such sets do not satisfy the inductive conditions. Nevertheless, we can show that for c = 1, B(x) is a Type 2 ε-barrier. The dReal query verifies the conditions with ε = 0.1. Since U (x) → B(x) > and init(x) → B(x) < −ε , we know that the system cannot reach any unsafe states. Figure 3 (Right), illustrates the example. The green set at the center represents init, and the red set represents unsafe set U . The B(x) = 0 level set is not invariant, as evidenced in the figure by the forward images at t = 0.14 and t = 0.28 leaving the set; however, as the dReal query proves, the reachable set over 0 ≤ t ≤ 10 does not leave the B(x) = 1.0 level set and is completely contained in the B(x) = −0.1 level set by t = 0.4. Since U (x) → B(x) > 1.0 and init(x) → B(x) < −0.1, then the system cannot reach any state in U .

Experiments
In this section, we show examples of nonlinear systems that can be verified to be ε-stable or safe with ε-barriers.   Table 2: Results for the ε-barrier functions. Each barrier function B(x) is of the form z T P z − , where z is a vector of monomials over x. We indicate the highest degree of the monomials used in z, the size of the P , the level used for each barrier function, and the value of ε and γ used to the check ∇ f B(x) < −γ.  Table 2 contains parameters we use to verify requirements of Definition 11 for Type 1 ε-barrier functions in our examples. The ε-Lyapunov functions in these examples are of the form V (x) := z T P z, where z is a vector of products of the state variables and P is a constant matrix obtained using simulation-guided techniques from [10]. All the P matrices are given in [8].
Time-Reversed Van der Pol. The time-reversed Van der Pol system has been used as an example in the previous sections. Figure 3 (Left) shows the direction field of this system around the origin. Using dReal with δ := 10 −25 , we are able to establish a 10 −12 -Lyapunov function and a 10 −5 -barrier function.
Normalized Pendulum. A standard pendulum system has continuous dynamics containing a transcendental function, which causes difficulty for many techniques. Here, we consider a normalized pendulum system with the following dynamics, in which x 1 and x 2 represent angular position and velocity, respectively. In our experiment, using δ = 10 −50 , we can prove that function V := x T P x is ε-Lyapunov, where ε := 10 −12 .
Using δ := 0.01, we are able to prove that for any value ∈ [0.1, 10], the function B(x) := x T P x − , with x being the system state, and P a constant matrix given in [8], is a Type 1 0.01-barrier function.
Moore-Greitzer Jet Engine. Next, we consider a simplified version of the Moore-Greitzer model for a jet engine. The system has the following dynamics, in which x 1 and x 2 are states related to mass flow and pressure rise.
Using dReal with δ := 0.1, we are able to prove that for any value ∈ [1, 10], the function B(x) := z T P z − , with x being the system state, z being the vector of monomials defined in the previous section, and P a constant matrix given in [8], is a Type 1 0.1-barrier function.
Powertrain Control System. Next, we consider a closed-loop model of a powertain control (PTC) system for an automotive application. The system dynamics consist of four state variables, two associated with a plant and two for a controller. The plant models fuel and air dynamics of an internal combustion engine and the controller is designed to regulate the air-fuel (A/F) ratio within a given range of an optimal value, referred as stoichiometric value. Two states related to the plant represent the manifold pressure, p, and the ratio between actual A/F ratio and stoichiometric value, r. The two associated with the controller are the estimated manifold pressure, p est , and the internal state of the PI controller, i. The system is highly nonlinear, with the following dynamicṡ c3 + c4c2p + c5c2p 2 + c6c 2 2 p c13(c3 + c4c2pest + c5c2p 2 est + c6c 2 2 pest)(1 + i + c14(r − c16)) − r ṗest = c1 2û1

Conclusion
We formulated new inductive proof rules for stability and safety for dynamical systems. The rules are numerically robust, making them amenable to verification using automated reasoning tools such as those based on δ-decision procedures. We presented several examples demonstrating the value of the new approach, including safety verification tasks for highly nonlinear systems. The examples show that the framework can be used to prove stability and safety for examples that were out of reach for existing tools. The new framework relies on the ability to generate reasonable candidate Lyapunov functions, which are analogous to ranking functions from program analysis. Future work will include improved techniques for efficiently generating the ε-Lyapunov and ε-barrier functions and related theoretical questions.